Coverage for /usr/lib/python3/dist-packages/scipy/sparse/_dia.py: 17%

1"""Sparse DIAgonal format""" 

2 

3__docformat__ = "restructuredtext en" 

4 

5__all__ = ['dia_array', 'dia_matrix', 'isspmatrix_dia'] 

6 

7import numpy as np 

8 

9from ._matrix import spmatrix, _array_doc_to_matrix 

10from ._base import issparse, _formats, _spbase, sparray 

11from ._data import _data_matrix 

12from ._sputils import (isshape, upcast_char, getdtype, get_sum_dtype, validateaxis, check_shape) 

13from ._sparsetools import dia_matvec 

14 

15 

16class _dia_base(_data_matrix): 

17 """Sparse matrix with DIAgonal storage 

18 

19 This can be instantiated in several ways: 

20 dia_array(D) 

21 with a dense matrix 

22 

23 dia_array(S) 

24 with another sparse matrix S (equivalent to S.todia()) 

25 

26 dia_array((M, N), [dtype]) 

27 to construct an empty matrix with shape (M, N), 

28 dtype is optional, defaulting to dtype='d'. 

29 

30 dia_array((data, offsets), shape=(M, N)) 

31 where the ``data[k,:]`` stores the diagonal entries for 

32 diagonal ``offsets[k]`` (See example below) 

33 

34 Attributes 

35 ---------- 

36 dtype : dtype 

37 Data type of the matrix 

38 shape : 2-tuple 

39 Shape of the matrix 

40 ndim : int 

41 Number of dimensions (this is always 2) 

42 nnz 

43 Number of stored values, including explicit zeros 

44 data 

45 DIA format data array of the matrix 

46 offsets 

47 DIA format offset array of the matrix 

48 

49 Notes 

50 ----- 

51 

52 Sparse matrices can be used in arithmetic operations: they support 

53 addition, subtraction, multiplication, division, and matrix power. 

54 

55 Examples 

56 -------- 

57 

58 >>> import numpy as np 

59 >>> from scipy.sparse import dia_array 

60 >>> dia_array((3, 4), dtype=np.int8).toarray() 

61 array([[0, 0, 0, 0], 

62 [0, 0, 0, 0], 

63 [0, 0, 0, 0]], dtype=int8) 

64 

65 >>> data = np.array([[1, 2, 3, 4]]).repeat(3, axis=0) 

66 >>> offsets = np.array([0, -1, 2]) 

67 >>> dia_array((data, offsets), shape=(4, 4)).toarray() 

68 array([[1, 0, 3, 0], 

69 [1, 2, 0, 4], 

70 [0, 2, 3, 0], 

71 [0, 0, 3, 4]]) 

72 

73 >>> from scipy.sparse import dia_array 

74 >>> n = 10 

75 >>> ex = np.ones(n) 

76 >>> data = np.array([ex, 2 * ex, ex]) 

77 >>> offsets = np.array([-1, 0, 1]) 

78 >>> dia_array((data, offsets), shape=(n, n)).toarray() 

79 array([[2., 1., 0., ..., 0., 0., 0.], 

80 [1., 2., 1., ..., 0., 0., 0.], 

81 [0., 1., 2., ..., 0., 0., 0.], 

82 ..., 

83 [0., 0., 0., ..., 2., 1., 0.], 

84 [0., 0., 0., ..., 1., 2., 1.], 

85 [0., 0., 0., ..., 0., 1., 2.]]) 

86 """ 

87 _format = 'dia' 

88 

89 def __init__(self, arg1, shape=None, dtype=None, copy=False): 

90 _data_matrix.__init__(self) 

91 

92 if issparse(arg1): 

93 if arg1.format == "dia": 

94 if copy: 

95 arg1 = arg1.copy() 

96 self.data = arg1.data 

97 self.offsets = arg1.offsets 

98 self._shape = check_shape(arg1.shape) 

99 else: 

100 if arg1.format == self.format and copy: 

101 A = arg1.copy() 

102 else: 

103 A = arg1.todia() 

104 self.data = A.data 

105 self.offsets = A.offsets 

106 self._shape = check_shape(A.shape) 

107 elif isinstance(arg1, tuple): 

108 if isshape(arg1): 

109 # It's a tuple of matrix dimensions (M, N) 

110 # create empty matrix 

111 self._shape = check_shape(arg1) 

112 self.data = np.zeros((0,0), getdtype(dtype, default=float)) 

113 idx_dtype = self._get_index_dtype(maxval=max(self.shape)) 

114 self.offsets = np.zeros((0), dtype=idx_dtype) 

115 else: 

116 try: 

117 # Try interpreting it as (data, offsets) 

118 data, offsets = arg1 

119 except Exception as e: 

120 raise ValueError('unrecognized form for dia_array constructor') from e 

121 else: 

122 if shape is None: 

123 raise ValueError('expected a shape argument') 

124 self.data = np.atleast_2d(np.array(arg1[0], dtype=dtype, copy=copy)) 

125 self.offsets = np.atleast_1d(np.array(arg1[1], 

126 dtype=self._get_index_dtype(maxval=max(shape)), 

127 copy=copy)) 

128 self._shape = check_shape(shape) 

129 else: 

130 #must be dense, convert to COO first, then to DIA 

131 try: 

132 arg1 = np.asarray(arg1) 

133 except Exception as e: 

134 raise ValueError("unrecognized form for" 

135 " %s_matrix constructor" % self.format) from e 

136 A = self._coo_container(arg1, dtype=dtype, shape=shape).todia() 

137 self.data = A.data 

138 self.offsets = A.offsets 

139 self._shape = check_shape(A.shape) 

140 

141 if dtype is not None: 

142 self.data = self.data.astype(dtype) 

143 

144 #check format 

145 if self.offsets.ndim != 1: 

146 raise ValueError('offsets array must have rank 1') 

147 

148 if self.data.ndim != 2: 

149 raise ValueError('data array must have rank 2') 

150 

151 if self.data.shape[0] != len(self.offsets): 

152 raise ValueError('number of diagonals (%d) ' 

153 'does not match the number of offsets (%d)' 

154 % (self.data.shape[0], len(self.offsets))) 

155 

156 if len(np.unique(self.offsets)) != len(self.offsets): 

157 raise ValueError('offset array contains duplicate values') 

158 

159 def __repr__(self): 

160 format = _formats[self.format][1] 

161 return "<%dx%d sparse matrix of type '%s'\n" \ 

162 "\twith %d stored elements (%d diagonals) in %s format>" % \ 

163 (self.shape + (self.dtype.type, self.nnz, self.data.shape[0], 

164 format)) 

165 

166 def _data_mask(self): 

167 """Returns a mask of the same shape as self.data, where 

168 mask[i,j] is True when data[i,j] corresponds to a stored element.""" 

169 num_rows, num_cols = self.shape 

170 offset_inds = np.arange(self.data.shape[1]) 

171 row = offset_inds - self.offsets[:,None] 

172 mask = (row >= 0) 

173 mask &= (row < num_rows) 

174 mask &= (offset_inds < num_cols) 

175 return mask 

176 

177 def count_nonzero(self): 

178 mask = self._data_mask() 

179 return np.count_nonzero(self.data[mask]) 

180 

181 def _getnnz(self, axis=None): 

182 if axis is not None: 

183 raise NotImplementedError("_getnnz over an axis is not implemented " 

184 "for DIA format") 

185 M,N = self.shape 

186 nnz = 0 

187 for k in self.offsets: 

188 if k > 0: 

189 nnz += min(M,N-k) 

190 else: 

191 nnz += min(M+k,N) 

192 return int(nnz) 

193 

194 _getnnz.__doc__ = _spbase._getnnz.__doc__ 

195 count_nonzero.__doc__ = _spbase.count_nonzero.__doc__ 

196 

197 def sum(self, axis=None, dtype=None, out=None): 

198 validateaxis(axis) 

199 

200 if axis is not None and axis < 0: 

201 axis += 2 

202 

203 res_dtype = get_sum_dtype(self.dtype) 

204 num_rows, num_cols = self.shape 

205 ret = None 

206 

207 if axis == 0: 

208 mask = self._data_mask() 

209 x = (self.data * mask).sum(axis=0) 

210 if x.shape[0] == num_cols: 

211 res = x 

212 else: 

213 res = np.zeros(num_cols, dtype=x.dtype) 

214 res[:x.shape[0]] = x 

215 ret = self._ascontainer(res, dtype=res_dtype) 

216 

217 else: 

218 row_sums = np.zeros((num_rows, 1), dtype=res_dtype) 

219 one = np.ones(num_cols, dtype=res_dtype) 

220 dia_matvec(num_rows, num_cols, len(self.offsets), 

221 self.data.shape[1], self.offsets, self.data, one, row_sums) 

222 

223 row_sums = self._ascontainer(row_sums) 

224 

225 if axis is None: 

226 return row_sums.sum(dtype=dtype, out=out) 

227 

228 ret = self._ascontainer(row_sums.sum(axis=axis)) 

229 

230 if out is not None and out.shape != ret.shape: 

231 raise ValueError("dimensions do not match") 

232 

233 return ret.sum(axis=(), dtype=dtype, out=out) 

234 

235 sum.__doc__ = _spbase.sum.__doc__ 

236 

237 def _add_sparse(self, other): 

238 

239 # Check if other is also of type dia_array 

240 if not isinstance(other, type(self)): 

241 # If other is not of type dia_array, default to 

242 # converting to csr_matrix, as is done in the _add_sparse 

243 # method of parent class _spbase 

244 return self.tocsr()._add_sparse(other) 

245 

246 # The task is to compute m = self + other 

247 # Start by making a copy of self, of the datatype 

248 # that should result from adding self and other 

249 dtype = np.promote_types(self.dtype, other.dtype) 

250 m = self.astype(dtype, copy=True) 

251 

252 # Then, add all the stored diagonals of other. 

253 for d in other.offsets: 

254 # Check if the diagonal has already been added. 

255 if d in m.offsets: 

256 # If the diagonal is already there, we need to take 

257 # the sum of the existing and the new 

258 m.setdiag(m.diagonal(d) + other.diagonal(d), d) 

259 else: 

260 m.setdiag(other.diagonal(d), d) 

261 return m 

262 

263 def _mul_vector(self, other): 

264 x = other 

265 

266 y = np.zeros(self.shape[0], dtype=upcast_char(self.dtype.char, 

267 x.dtype.char)) 

268 

269 L = self.data.shape[1] 

270 

271 M,N = self.shape 

272 

273 dia_matvec(M,N, len(self.offsets), L, self.offsets, self.data, x.ravel(), y.ravel()) 

274 

275 return y 

276 

277 def _mul_multimatrix(self, other): 

278 return np.hstack([self._mul_vector(col).reshape(-1,1) for col in other.T]) 

279 

280 def _setdiag(self, values, k=0): 

281 M, N = self.shape 

282 

283 if values.ndim == 0: 

284 # broadcast 

285 values_n = np.inf 

286 else: 

287 values_n = len(values) 

288 

289 if k < 0: 

290 n = min(M + k, N, values_n) 

291 min_index = 0 

292 max_index = n 

293 else: 

294 n = min(M, N - k, values_n) 

295 min_index = k 

296 max_index = k + n 

297 

298 if values.ndim != 0: 

299 # allow also longer sequences 

300 values = values[:n] 

301 

302 data_rows, data_cols = self.data.shape 

303 if k in self.offsets: 

304 if max_index > data_cols: 

305 data = np.zeros((data_rows, max_index), dtype=self.data.dtype) 

306 data[:, :data_cols] = self.data 

307 self.data = data 

308 self.data[self.offsets == k, min_index:max_index] = values 

309 else: 

310 self.offsets = np.append(self.offsets, self.offsets.dtype.type(k)) 

311 m = max(max_index, data_cols) 

312 data = np.zeros((data_rows + 1, m), dtype=self.data.dtype) 

313 data[:-1, :data_cols] = self.data 

314 data[-1, min_index:max_index] = values 

315 self.data = data 

316 

317 def todia(self, copy=False): 

318 if copy: 

319 return self.copy() 

320 else: 

321 return self 

322 

323 todia.__doc__ = _spbase.todia.__doc__ 

324 

325 def transpose(self, axes=None, copy=False): 

326 if axes is not None: 

327 raise ValueError("Sparse matrices do not support " 

328 "an 'axes' parameter because swapping " 

329 "dimensions is the only logical permutation.") 

330 

331 num_rows, num_cols = self.shape 

332 max_dim = max(self.shape) 

333 

334 # flip diagonal offsets 

335 offsets = -self.offsets 

336 

337 # re-align the data matrix 

338 r = np.arange(len(offsets), dtype=np.intc)[:, None] 

339 c = np.arange(num_rows, dtype=np.intc) - (offsets % max_dim)[:, None] 

340 pad_amount = max(0, max_dim-self.data.shape[1]) 

341 data = np.hstack((self.data, np.zeros((self.data.shape[0], pad_amount), 

342 dtype=self.data.dtype))) 

343 data = data[r, c] 

344 return self._dia_container((data, offsets), shape=( 

345 num_cols, num_rows), copy=copy) 

346 

347 transpose.__doc__ = _spbase.transpose.__doc__ 

348 

349 def diagonal(self, k=0): 

350 rows, cols = self.shape 

351 if k <= -rows or k >= cols: 

352 return np.empty(0, dtype=self.data.dtype) 

353 idx, = np.nonzero(self.offsets == k) 

354 first_col = max(0, k) 

355 last_col = min(rows + k, cols) 

356 result_size = last_col - first_col 

357 if idx.size == 0: 

358 return np.zeros(result_size, dtype=self.data.dtype) 

359 result = self.data[idx[0], first_col:last_col] 

360 padding = result_size - len(result) 

361 if padding > 0: 

362 result = np.pad(result, (0, padding), mode='constant') 

363 return result 

364 

365 diagonal.__doc__ = _spbase.diagonal.__doc__ 

366 

367 def tocsc(self, copy=False): 

368 if self.nnz == 0: 

369 return self._csc_container(self.shape, dtype=self.dtype) 

370 

371 num_rows, num_cols = self.shape 

372 num_offsets, offset_len = self.data.shape 

373 offset_inds = np.arange(offset_len) 

374 

375 row = offset_inds - self.offsets[:,None] 

376 mask = (row >= 0) 

377 mask &= (row < num_rows) 

378 mask &= (offset_inds < num_cols) 

379 mask &= (self.data != 0) 

380 

381 idx_dtype = self._get_index_dtype(maxval=max(self.shape)) 

382 indptr = np.zeros(num_cols + 1, dtype=idx_dtype) 

383 indptr[1:offset_len+1] = np.cumsum(mask.sum(axis=0)[:num_cols]) 

384 if offset_len < num_cols: 

385 indptr[offset_len+1:] = indptr[offset_len] 

386 indices = row.T[mask.T].astype(idx_dtype, copy=False) 

387 data = self.data.T[mask.T] 

388 return self._csc_container((data, indices, indptr), shape=self.shape, 

389 dtype=self.dtype) 

390 

391 tocsc.__doc__ = _spbase.tocsc.__doc__ 

392 

393 def tocoo(self, copy=False): 

394 num_rows, num_cols = self.shape 

395 num_offsets, offset_len = self.data.shape 

396 offset_inds = np.arange(offset_len) 

397 

398 row = offset_inds - self.offsets[:,None] 

399 mask = (row >= 0) 

400 mask &= (row < num_rows) 

401 mask &= (offset_inds < num_cols) 

402 mask &= (self.data != 0) 

403 row = row[mask] 

404 col = np.tile(offset_inds, num_offsets)[mask.ravel()] 

405 idx_dtype = self._get_index_dtype( 

406 arrays=(self.offsets,), maxval=max(self.shape) 

407 ) 

408 row = row.astype(idx_dtype, copy=False) 

409 col = col.astype(idx_dtype, copy=False) 

410 data = self.data[mask] 

411 # Note: this cannot set has_canonical_format=True, because despite the 

412 # lack of duplicates, we do not generate sorted indices. 

413 return self._coo_container( 

414 (data, (row, col)), shape=self.shape, dtype=self.dtype, copy=False 

415 ) 

416 

417 tocoo.__doc__ = _spbase.tocoo.__doc__ 

418 

419 # needed by _data_matrix 

420 def _with_data(self, data, copy=True): 

421 """Returns a matrix with the same sparsity structure as self, 

422 but with different data. By default the structure arrays are copied. 

423 """ 

424 if copy: 

425 return self._dia_container( 

426 (data, self.offsets.copy()), shape=self.shape 

427 ) 

428 else: 

429 return self._dia_container( 

430 (data, self.offsets), shape=self.shape 

431 ) 

432 

433 def resize(self, *shape): 

434 shape = check_shape(shape) 

435 M, N = shape 

436 # we do not need to handle the case of expanding N 

437 self.data = self.data[:, :N] 

438 

439 if (M > self.shape[0] and 

440 np.any(self.offsets + self.shape[0] < self.data.shape[1])): 

441 # explicitly clear values that were previously hidden 

442 mask = (self.offsets[:, None] + self.shape[0] <= 

443 np.arange(self.data.shape[1])) 

444 self.data[mask] = 0 

445 

446 self._shape = shape 

447 

448 resize.__doc__ = _spbase.resize.__doc__ 

449 

450 

451def isspmatrix_dia(x): 

452 """Is `x` of dia_matrix type? 

453 

454 Parameters 

455 ---------- 

456 x 

457 object to check for being a dia matrix 

458 

459 Returns 

460 ------- 

461 bool 

462 True if `x` is a dia matrix, False otherwise 

463 

464 Examples 

465 -------- 

466 >>> from scipy.sparse import dia_array, dia_matrix, coo_matrix, isspmatrix_dia 

467 >>> isspmatrix_dia(dia_matrix([[5]])) 

468 True 

469 >>> isspmatrix_dia(dia_array([[5]])) 

470 False 

471 >>> isspmatrix_dia(coo_matrix([[5]])) 

472 False 

473 """ 

474 return isinstance(x, dia_matrix) 

475 

476 

477# This namespace class separates array from matrix with isinstance 

478class dia_array(_dia_base, sparray): 

479 pass 

480 

481dia_array.__doc__ = _dia_base.__doc__ 

482 

483class dia_matrix(spmatrix, _dia_base): 

484 pass 

485 

486dia_matrix.__doc__ = _array_doc_to_matrix(_dia_base.__doc__)