Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/scipy/sparse/_dia.py: 15%

1"""Sparse DIAgonal format""" 

2 

3__docformat__ = "restructuredtext en" 

4 

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

6 

7import numpy as np 

8 

9from ._base import isspmatrix, _formats, spmatrix 

10from ._data import _data_matrix 

11from ._sputils import (isshape, upcast_char, getdtype, get_index_dtype, 

12 get_sum_dtype, validateaxis, check_shape) 

13from ._sparsetools import dia_matvec 

14 

15 

16class dia_matrix(_data_matrix): 

17 """Sparse matrix with DIAgonal storage 

18 

19 This can be instantiated in several ways: 

20 dia_matrix(D) 

21 with a dense matrix 

22 

23 dia_matrix(S) 

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

25 

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

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

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

29 

30 dia_matrix((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_matrix 

60 >>> dia_matrix((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_matrix((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_matrix 

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_matrix((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 isspmatrix_dia(arg1): 

93 if copy: 

94 arg1 = arg1.copy() 

95 self.data = arg1.data 

96 self.offsets = arg1.offsets 

97 self._shape = check_shape(arg1.shape) 

98 elif isspmatrix(arg1): 

99 if isspmatrix_dia(arg1) and copy: 

100 A = arg1.copy() 

101 else: 

102 A = arg1.todia() 

103 self.data = A.data 

104 self.offsets = A.offsets 

105 self._shape = check_shape(A.shape) 

106 elif isinstance(arg1, tuple): 

107 if isshape(arg1): 

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

109 # create empty matrix 

110 self._shape = check_shape(arg1) 

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

112 idx_dtype = get_index_dtype(maxval=max(self.shape)) 

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

114 else: 

115 try: 

116 # Try interpreting it as (data, offsets) 

117 data, offsets = arg1 

118 except Exception as e: 

119 raise ValueError('unrecognized form for dia_matrix constructor') from e 

120 else: 

121 if shape is None: 

122 raise ValueError('expected a shape argument') 

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

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

125 dtype=get_index_dtype(maxval=max(shape)), 

126 copy=copy)) 

127 self._shape = check_shape(shape) 

128 else: 

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

130 try: 

131 arg1 = np.asarray(arg1) 

132 except Exception as e: 

133 raise ValueError("unrecognized form for" 

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

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

136 self.data = A.data 

137 self.offsets = A.offsets 

138 self._shape = check_shape(A.shape) 

139 

140 if dtype is not None: 

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

142 

143 #check format 

144 if self.offsets.ndim != 1: 

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

146 

147 if self.data.ndim != 2: 

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

149 

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

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

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

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

154 

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

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

157 

158 def __repr__(self): 

159 format = _formats[self.getformat()][1] 

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

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

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

163 format)) 

164 

165 def _data_mask(self): 

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

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

168 num_rows, num_cols = self.shape 

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

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

171 mask = (row >= 0) 

172 mask &= (row < num_rows) 

173 mask &= (offset_inds < num_cols) 

174 return mask 

175 

176 def count_nonzero(self): 

177 mask = self._data_mask() 

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

179 

180 def getnnz(self, axis=None): 

181 if axis is not None: 

182 raise NotImplementedError("getnnz over an axis is not implemented " 

183 "for DIA format") 

184 M,N = self.shape 

185 nnz = 0 

186 for k in self.offsets: 

187 if k > 0: 

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

189 else: 

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

191 return int(nnz) 

192 

193 getnnz.__doc__ = spmatrix.getnnz.__doc__ 

194 count_nonzero.__doc__ = spmatrix.count_nonzero.__doc__ 

195 

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

197 validateaxis(axis) 

198 

199 if axis is not None and axis < 0: 

200 axis += 2 

201 

202 res_dtype = get_sum_dtype(self.dtype) 

203 num_rows, num_cols = self.shape 

204 ret = None 

205 

206 if axis == 0: 

207 mask = self._data_mask() 

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

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

210 res = x 

211 else: 

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

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

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

215 

216 else: 

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

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

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

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

221 

222 row_sums = self._ascontainer(row_sums) 

223 

224 if axis is None: 

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

226 

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

228 

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

230 raise ValueError("dimensions do not match") 

231 

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

233 

234 sum.__doc__ = spmatrix.sum.__doc__ 

235 

236 def _add_sparse(self, other): 

237 

238 # Check if other is also of type dia_matrix 

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

240 # If other is not of type dia_matrix, default to 

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

242 # method of parent class spmatrix 

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

244 

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

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

247 # that should result from adding self and other 

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

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

250 

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

252 for d in other.offsets: 

253 # Check if the diagonal has already been added. 

254 if d in m.offsets: 

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

256 # the sum of the existing and the new 

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

258 else: 

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

260 return m 

261 

262 def _mul_vector(self, other): 

263 x = other 

264 

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

266 x.dtype.char)) 

267 

268 L = self.data.shape[1] 

269 

270 M,N = self.shape 

271 

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

273 

274 return y 

275 

276 def _mul_multimatrix(self, other): 

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

278 

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

280 M, N = self.shape 

281 

282 if values.ndim == 0: 

283 # broadcast 

284 values_n = np.inf 

285 else: 

286 values_n = len(values) 

287 

288 if k < 0: 

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

290 min_index = 0 

291 max_index = n 

292 else: 

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

294 min_index = k 

295 max_index = k + n 

296 

297 if values.ndim != 0: 

298 # allow also longer sequences 

299 values = values[:n] 

300 

301 data_rows, data_cols = self.data.shape 

302 if k in self.offsets: 

303 if max_index > data_cols: 

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

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

306 self.data = data 

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

308 else: 

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

310 m = max(max_index, data_cols) 

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

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

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

314 self.data = data 

315 

316 def todia(self, copy=False): 

317 if copy: 

318 return self.copy() 

319 else: 

320 return self 

321 

322 todia.__doc__ = spmatrix.todia.__doc__ 

323 

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

325 if axes is not None: 

326 raise ValueError(("Sparse matrices do not support " 

327 "an 'axes' parameter because swapping " 

328 "dimensions is the only logical permutation.")) 

329 

330 num_rows, num_cols = self.shape 

331 max_dim = max(self.shape) 

332 

333 # flip diagonal offsets 

334 offsets = -self.offsets 

335 

336 # re-align the data matrix 

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

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

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

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

341 dtype=self.data.dtype))) 

342 data = data[r, c] 

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

344 num_cols, num_rows), copy=copy) 

345 

346 transpose.__doc__ = spmatrix.transpose.__doc__ 

347 

348 def diagonal(self, k=0): 

349 rows, cols = self.shape 

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

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

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

353 first_col = max(0, k) 

354 last_col = min(rows + k, cols) 

355 result_size = last_col - first_col 

356 if idx.size == 0: 

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

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

359 padding = result_size - len(result) 

360 if padding > 0: 

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

362 return result 

363 

364 diagonal.__doc__ = spmatrix.diagonal.__doc__ 

365 

366 def tocsc(self, copy=False): 

367 if self.nnz == 0: 

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

369 

370 num_rows, num_cols = self.shape 

371 num_offsets, offset_len = self.data.shape 

372 offset_inds = np.arange(offset_len) 

373 

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

375 mask = (row >= 0) 

376 mask &= (row < num_rows) 

377 mask &= (offset_inds < num_cols) 

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

379 

380 idx_dtype = get_index_dtype(maxval=max(self.shape)) 

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

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

383 if offset_len < num_cols: 

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

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

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

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

388 dtype=self.dtype) 

389 

390 tocsc.__doc__ = spmatrix.tocsc.__doc__ 

391 

392 def tocoo(self, copy=False): 

393 num_rows, num_cols = self.shape 

394 num_offsets, offset_len = self.data.shape 

395 offset_inds = np.arange(offset_len) 

396 

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

398 mask = (row >= 0) 

399 mask &= (row < num_rows) 

400 mask &= (offset_inds < num_cols) 

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

402 row = row[mask] 

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

404 data = self.data[mask] 

405 

406 A = self._coo_container( 

407 (data, (row, col)), shape=self.shape, dtype=self.dtype 

408 ) 

409 A.has_canonical_format = True 

410 return A 

411 

412 tocoo.__doc__ = spmatrix.tocoo.__doc__ 

413 

414 # needed by _data_matrix 

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

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

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

418 """ 

419 if copy: 

420 return self._dia_container( 

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

422 ) 

423 else: 

424 return self._dia_container( 

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

426 ) 

427 

428 def resize(self, *shape): 

429 shape = check_shape(shape) 

430 M, N = shape 

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

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

433 

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

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

436 # explicitly clear values that were previously hidden 

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

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

439 self.data[mask] = 0 

440 

441 self._shape = shape 

442 

443 resize.__doc__ = spmatrix.resize.__doc__ 

444 

445 

446def isspmatrix_dia(x): 

447 """Is x of dia_matrix type? 

448 

449 Parameters 

450 ---------- 

451 x 

452 object to check for being a dia matrix 

453 

454 Returns 

455 ------- 

456 bool 

457 True if x is a dia matrix, False otherwise 

458 

459 Examples 

460 -------- 

461 >>> from scipy.sparse import dia_matrix, isspmatrix_dia 

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

463 True 

464 

465 >>> from scipy.sparse import dia_matrix, csr_matrix, isspmatrix_dia 

466 >>> isspmatrix_dia(csr_matrix([[5]])) 

467 False 

468 """ 

469 from ._arrays import dia_array 

470 return isinstance(x, dia_matrix) or isinstance(x, dia_array)