Coverage for /usr/lib/python3/dist-packages/scipy/sparse/_coo.py: 16%

299 statements  

« prev     ^ index     » next       coverage.py v7.4.4, created at 2025-06-14 15:25 +0200

1""" A sparse matrix in COOrdinate or 'triplet' format""" 

2 

3__docformat__ = "restructuredtext en" 

4 

5__all__ = ['coo_array', 'coo_matrix', 'isspmatrix_coo'] 

6 

7from warnings import warn 

8 

9import numpy as np 

10 

11from ._matrix import spmatrix, _array_doc_to_matrix 

12from ._sparsetools import coo_tocsr, coo_todense, coo_matvec 

13from ._base import issparse, SparseEfficiencyWarning, _spbase, sparray 

14from ._data import _data_matrix, _minmax_mixin 

15from ._sputils import (upcast, upcast_char, to_native, isshape, getdtype, 

16 getdata, downcast_intp_index, 

17 check_shape, check_reshape_kwargs) 

18 

19import operator 

20 

21 

22class _coo_base(_data_matrix, _minmax_mixin): 

23 """ 

24 A sparse matrix in COOrdinate format. 

25 

26 Also known as the 'ijv' or 'triplet' format. 

27 

28 This can be instantiated in several ways: 

29 coo_array(D) 

30 with a dense matrix D 

31 

32 coo_array(S) 

33 with another sparse matrix S (equivalent to S.tocoo()) 

34 

35 coo_array((M, N), [dtype]) 

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

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

38 

39 coo_array((data, (i, j)), [shape=(M, N)]) 

40 to construct from three arrays: 

41 1. data[:] the entries of the matrix, in any order 

42 2. i[:] the row indices of the matrix entries 

43 3. j[:] the column indices of the matrix entries 

44 

45 Where ``A[i[k], j[k]] = data[k]``. When shape is not 

46 specified, it is inferred from the index arrays 

47 

48 Attributes 

49 ---------- 

50 dtype : dtype 

51 Data type of the matrix 

52 shape : 2-tuple 

53 Shape of the matrix 

54 ndim : int 

55 Number of dimensions (this is always 2) 

56 nnz 

57 Number of stored values, including explicit zeros 

58 data 

59 COO format data array of the matrix 

60 row 

61 COO format row index array of the matrix 

62 col 

63 COO format column index array of the matrix 

64 

65 Notes 

66 ----- 

67 

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

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

70 

71 Advantages of the COO format 

72 - facilitates fast conversion among sparse formats 

73 - permits duplicate entries (see example) 

74 - very fast conversion to and from CSR/CSC formats 

75 

76 Disadvantages of the COO format 

77 - does not directly support: 

78 + arithmetic operations 

79 + slicing 

80 

81 Intended Usage 

82 - COO is a fast format for constructing sparse matrices 

83 - Once a matrix has been constructed, convert to CSR or 

84 CSC format for fast arithmetic and matrix vector operations 

85 - By default when converting to CSR or CSC format, duplicate (i,j) 

86 entries will be summed together. This facilitates efficient 

87 construction of finite element matrices and the like. (see example) 

88 

89 Canonical format 

90 - Entries and indices sorted by row, then column. 

91 - There are no duplicate entries (i.e. duplicate (i,j) locations) 

92 - Arrays MAY have explicit zeros. 

93 

94 Examples 

95 -------- 

96 

97 >>> # Constructing an empty matrix 

98 >>> import numpy as np 

99 >>> from scipy.sparse import coo_array 

100 >>> coo_array((3, 4), dtype=np.int8).toarray() 

101 array([[0, 0, 0, 0], 

102 [0, 0, 0, 0], 

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

104 

105 >>> # Constructing a matrix using ijv format 

106 >>> row = np.array([0, 3, 1, 0]) 

107 >>> col = np.array([0, 3, 1, 2]) 

108 >>> data = np.array([4, 5, 7, 9]) 

109 >>> coo_array((data, (row, col)), shape=(4, 4)).toarray() 

110 array([[4, 0, 9, 0], 

111 [0, 7, 0, 0], 

112 [0, 0, 0, 0], 

113 [0, 0, 0, 5]]) 

114 

115 >>> # Constructing a matrix with duplicate indices 

116 >>> row = np.array([0, 0, 1, 3, 1, 0, 0]) 

117 >>> col = np.array([0, 2, 1, 3, 1, 0, 0]) 

118 >>> data = np.array([1, 1, 1, 1, 1, 1, 1]) 

119 >>> coo = coo_array((data, (row, col)), shape=(4, 4)) 

120 >>> # Duplicate indices are maintained until implicitly or explicitly summed 

121 >>> np.max(coo.data) 

122 1 

123 >>> coo.toarray() 

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

125 [0, 2, 0, 0], 

126 [0, 0, 0, 0], 

127 [0, 0, 0, 1]]) 

128 

129 """ 

130 _format = 'coo' 

131 

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

133 _data_matrix.__init__(self) 

134 

135 if isinstance(arg1, tuple): 

136 if isshape(arg1): 

137 M, N = arg1 

138 self._shape = check_shape((M, N)) 

139 idx_dtype = self._get_index_dtype(maxval=max(M, N)) 

140 data_dtype = getdtype(dtype, default=float) 

141 self.row = np.array([], dtype=idx_dtype) 

142 self.col = np.array([], dtype=idx_dtype) 

143 self.data = np.array([], dtype=data_dtype) 

144 self.has_canonical_format = True 

145 else: 

146 try: 

147 obj, (row, col) = arg1 

148 except (TypeError, ValueError) as e: 

149 raise TypeError('invalid input format') from e 

150 

151 if shape is None: 

152 if len(row) == 0 or len(col) == 0: 

153 raise ValueError('cannot infer dimensions from zero ' 

154 'sized index arrays') 

155 M = operator.index(np.max(row)) + 1 

156 N = operator.index(np.max(col)) + 1 

157 self._shape = check_shape((M, N)) 

158 else: 

159 # Use 2 steps to ensure shape has length 2. 

160 M, N = shape 

161 self._shape = check_shape((M, N)) 

162 

163 idx_dtype = self._get_index_dtype((row, col), maxval=max(self.shape), check_contents=True) 

164 self.row = np.array(row, copy=copy, dtype=idx_dtype) 

165 self.col = np.array(col, copy=copy, dtype=idx_dtype) 

166 self.data = getdata(obj, copy=copy, dtype=dtype) 

167 self.has_canonical_format = False 

168 else: 

169 if issparse(arg1): 

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

171 self.row = arg1.row.copy() 

172 self.col = arg1.col.copy() 

173 self.data = arg1.data.copy() 

174 self._shape = check_shape(arg1.shape) 

175 else: 

176 coo = arg1.tocoo() 

177 self.row = coo.row 

178 self.col = coo.col 

179 self.data = coo.data 

180 self._shape = check_shape(coo.shape) 

181 self.has_canonical_format = False 

182 else: 

183 #dense argument 

184 M = np.atleast_2d(np.asarray(arg1)) 

185 

186 if M.ndim != 2: 

187 raise TypeError('expected dimension <= 2 array or matrix') 

188 

189 self._shape = check_shape(M.shape) 

190 if shape is not None: 

191 if check_shape(shape) != self._shape: 

192 raise ValueError('inconsistent shapes: %s != %s' % 

193 (shape, self._shape)) 

194 index_dtype = self._get_index_dtype(maxval=max(self._shape)) 

195 row, col = M.nonzero() 

196 self.row = row.astype(index_dtype, copy=False) 

197 self.col = col.astype(index_dtype, copy=False) 

198 self.data = M[self.row, self.col] 

199 self.has_canonical_format = True 

200 

201 if dtype is not None: 

202 self.data = self.data.astype(dtype, copy=False) 

203 

204 self._check() 

205 

206 def reshape(self, *args, **kwargs): 

207 shape = check_shape(args, self.shape) 

208 order, copy = check_reshape_kwargs(kwargs) 

209 

210 # Return early if reshape is not required 

211 if shape == self.shape: 

212 if copy: 

213 return self.copy() 

214 else: 

215 return self 

216 

217 nrows, ncols = self.shape 

218 

219 if order == 'C': 

220 # Upcast to avoid overflows: the coo_array constructor 

221 # below will downcast the results to a smaller dtype, if 

222 # possible. 

223 dtype = self._get_index_dtype(maxval=(ncols * max(0, nrows - 1) + max(0, ncols - 1))) 

224 

225 flat_indices = np.multiply(ncols, self.row, dtype=dtype) + self.col 

226 new_row, new_col = divmod(flat_indices, shape[1]) 

227 elif order == 'F': 

228 dtype = self._get_index_dtype(maxval=(nrows * max(0, ncols - 1) + max(0, nrows - 1))) 

229 

230 flat_indices = np.multiply(nrows, self.col, dtype=dtype) + self.row 

231 new_col, new_row = divmod(flat_indices, shape[0]) 

232 else: 

233 raise ValueError("'order' must be 'C' or 'F'") 

234 

235 # Handle copy here rather than passing on to the constructor so that no 

236 # copy will be made of new_row and new_col regardless 

237 if copy: 

238 new_data = self.data.copy() 

239 else: 

240 new_data = self.data 

241 

242 return self.__class__((new_data, (new_row, new_col)), 

243 shape=shape, copy=False) 

244 

245 reshape.__doc__ = _spbase.reshape.__doc__ 

246 

247 def _getnnz(self, axis=None): 

248 if axis is None: 

249 nnz = len(self.data) 

250 if nnz != len(self.row) or nnz != len(self.col): 

251 raise ValueError('row, column, and data array must all be the ' 

252 'same length') 

253 

254 if self.data.ndim != 1 or self.row.ndim != 1 or \ 

255 self.col.ndim != 1: 

256 raise ValueError('row, column, and data arrays must be 1-D') 

257 

258 return int(nnz) 

259 

260 if axis < 0: 

261 axis += 2 

262 if axis == 0: 

263 return np.bincount(downcast_intp_index(self.col), 

264 minlength=self.shape[1]) 

265 elif axis == 1: 

266 return np.bincount(downcast_intp_index(self.row), 

267 minlength=self.shape[0]) 

268 else: 

269 raise ValueError('axis out of bounds') 

270 

271 _getnnz.__doc__ = _spbase._getnnz.__doc__ 

272 

273 def _check(self): 

274 """ Checks data structure for consistency """ 

275 

276 # index arrays should have integer data types 

277 if self.row.dtype.kind != 'i': 

278 warn("row index array has non-integer dtype (%s) " 

279 % self.row.dtype.name) 

280 if self.col.dtype.kind != 'i': 

281 warn("col index array has non-integer dtype (%s) " 

282 % self.col.dtype.name) 

283 

284 idx_dtype = self._get_index_dtype((self.row, self.col), maxval=max(self.shape)) 

285 self.row = np.asarray(self.row, dtype=idx_dtype) 

286 self.col = np.asarray(self.col, dtype=idx_dtype) 

287 self.data = to_native(self.data) 

288 

289 if self.nnz > 0: 

290 if self.row.max() >= self.shape[0]: 

291 raise ValueError('row index exceeds matrix dimensions') 

292 if self.col.max() >= self.shape[1]: 

293 raise ValueError('column index exceeds matrix dimensions') 

294 if self.row.min() < 0: 

295 raise ValueError('negative row index found') 

296 if self.col.min() < 0: 

297 raise ValueError('negative column index found') 

298 

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

300 if axes is not None: 

301 raise ValueError("Sparse matrices do not support " 

302 "an 'axes' parameter because swapping " 

303 "dimensions is the only logical permutation.") 

304 

305 M, N = self.shape 

306 return self.__class__((self.data, (self.col, self.row)), 

307 shape=(N, M), copy=copy) 

308 

309 transpose.__doc__ = _spbase.transpose.__doc__ 

310 

311 def resize(self, *shape): 

312 shape = check_shape(shape) 

313 new_M, new_N = shape 

314 M, N = self.shape 

315 

316 if new_M < M or new_N < N: 

317 mask = np.logical_and(self.row < new_M, self.col < new_N) 

318 if not mask.all(): 

319 self.row = self.row[mask] 

320 self.col = self.col[mask] 

321 self.data = self.data[mask] 

322 

323 self._shape = shape 

324 

325 resize.__doc__ = _spbase.resize.__doc__ 

326 

327 def toarray(self, order=None, out=None): 

328 """See the docstring for `_spbase.toarray`.""" 

329 B = self._process_toarray_args(order, out) 

330 fortran = int(B.flags.f_contiguous) 

331 if not fortran and not B.flags.c_contiguous: 

332 raise ValueError("Output array must be C or F contiguous") 

333 M,N = self.shape 

334 coo_todense(M, N, self.nnz, self.row, self.col, self.data, 

335 B.ravel('A'), fortran) 

336 return B 

337 

338 def tocsc(self, copy=False): 

339 """Convert this matrix to Compressed Sparse Column format 

340 

341 Duplicate entries will be summed together. 

342 

343 Examples 

344 -------- 

345 >>> from numpy import array 

346 >>> from scipy.sparse import coo_array 

347 >>> row = array([0, 0, 1, 3, 1, 0, 0]) 

348 >>> col = array([0, 2, 1, 3, 1, 0, 0]) 

349 >>> data = array([1, 1, 1, 1, 1, 1, 1]) 

350 >>> A = coo_array((data, (row, col)), shape=(4, 4)).tocsc() 

351 >>> A.toarray() 

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

353 [0, 2, 0, 0], 

354 [0, 0, 0, 0], 

355 [0, 0, 0, 1]]) 

356 

357 """ 

358 if self.nnz == 0: 

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

360 else: 

361 M,N = self.shape 

362 idx_dtype = self._get_index_dtype( 

363 (self.col, self.row), maxval=max(self.nnz, M) 

364 ) 

365 row = self.row.astype(idx_dtype, copy=False) 

366 col = self.col.astype(idx_dtype, copy=False) 

367 

368 indptr = np.empty(N + 1, dtype=idx_dtype) 

369 indices = np.empty_like(row, dtype=idx_dtype) 

370 data = np.empty_like(self.data, dtype=upcast(self.dtype)) 

371 

372 coo_tocsr(N, M, self.nnz, col, row, self.data, 

373 indptr, indices, data) 

374 

375 x = self._csc_container((data, indices, indptr), shape=self.shape) 

376 if not self.has_canonical_format: 

377 x.sum_duplicates() 

378 return x 

379 

380 def tocsr(self, copy=False): 

381 """Convert this matrix to Compressed Sparse Row format 

382 

383 Duplicate entries will be summed together. 

384 

385 Examples 

386 -------- 

387 >>> from numpy import array 

388 >>> from scipy.sparse import coo_array 

389 >>> row = array([0, 0, 1, 3, 1, 0, 0]) 

390 >>> col = array([0, 2, 1, 3, 1, 0, 0]) 

391 >>> data = array([1, 1, 1, 1, 1, 1, 1]) 

392 >>> A = coo_array((data, (row, col)), shape=(4, 4)).tocsr() 

393 >>> A.toarray() 

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

395 [0, 2, 0, 0], 

396 [0, 0, 0, 0], 

397 [0, 0, 0, 1]]) 

398 

399 """ 

400 if self.nnz == 0: 

401 return self._csr_container(self.shape, dtype=self.dtype) 

402 else: 

403 M,N = self.shape 

404 idx_dtype = self._get_index_dtype( 

405 (self.row, self.col), maxval=max(self.nnz, N) 

406 ) 

407 row = self.row.astype(idx_dtype, copy=False) 

408 col = self.col.astype(idx_dtype, copy=False) 

409 

410 indptr = np.empty(M + 1, dtype=idx_dtype) 

411 indices = np.empty_like(col, dtype=idx_dtype) 

412 data = np.empty_like(self.data, dtype=upcast(self.dtype)) 

413 

414 coo_tocsr(M, N, self.nnz, row, col, self.data, 

415 indptr, indices, data) 

416 

417 x = self._csr_container((data, indices, indptr), shape=self.shape) 

418 if not self.has_canonical_format: 

419 x.sum_duplicates() 

420 return x 

421 

422 def tocoo(self, copy=False): 

423 if copy: 

424 return self.copy() 

425 else: 

426 return self 

427 

428 tocoo.__doc__ = _spbase.tocoo.__doc__ 

429 

430 def todia(self, copy=False): 

431 self.sum_duplicates() 

432 ks = self.col - self.row # the diagonal for each nonzero 

433 diags, diag_idx = np.unique(ks, return_inverse=True) 

434 

435 if len(diags) > 100: 

436 # probably undesired, should todia() have a maxdiags parameter? 

437 warn("Constructing a DIA matrix with %d diagonals " 

438 "is inefficient" % len(diags), SparseEfficiencyWarning) 

439 

440 #initialize and fill in data array 

441 if self.data.size == 0: 

442 data = np.zeros((0, 0), dtype=self.dtype) 

443 else: 

444 data = np.zeros((len(diags), self.col.max()+1), dtype=self.dtype) 

445 data[diag_idx, self.col] = self.data 

446 

447 return self._dia_container((data, diags), shape=self.shape) 

448 

449 todia.__doc__ = _spbase.todia.__doc__ 

450 

451 def todok(self, copy=False): 

452 self.sum_duplicates() 

453 dok = self._dok_container((self.shape), dtype=self.dtype) 

454 dok._update(zip(zip(self.row,self.col),self.data)) 

455 

456 return dok 

457 

458 todok.__doc__ = _spbase.todok.__doc__ 

459 

460 def diagonal(self, k=0): 

461 rows, cols = self.shape 

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

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

464 diag = np.zeros(min(rows + min(k, 0), cols - max(k, 0)), 

465 dtype=self.dtype) 

466 diag_mask = (self.row + k) == self.col 

467 

468 if self.has_canonical_format: 

469 row = self.row[diag_mask] 

470 data = self.data[diag_mask] 

471 else: 

472 row, _, data = self._sum_duplicates(self.row[diag_mask], 

473 self.col[diag_mask], 

474 self.data[diag_mask]) 

475 diag[row + min(k, 0)] = data 

476 

477 return diag 

478 

479 diagonal.__doc__ = _data_matrix.diagonal.__doc__ 

480 

481 def _setdiag(self, values, k): 

482 M, N = self.shape 

483 if values.ndim and not len(values): 

484 return 

485 idx_dtype = self.row.dtype 

486 

487 # Determine which triples to keep and where to put the new ones. 

488 full_keep = self.col - self.row != k 

489 if k < 0: 

490 max_index = min(M+k, N) 

491 if values.ndim: 

492 max_index = min(max_index, len(values)) 

493 keep = np.logical_or(full_keep, self.col >= max_index) 

494 new_row = np.arange(-k, -k + max_index, dtype=idx_dtype) 

495 new_col = np.arange(max_index, dtype=idx_dtype) 

496 else: 

497 max_index = min(M, N-k) 

498 if values.ndim: 

499 max_index = min(max_index, len(values)) 

500 keep = np.logical_or(full_keep, self.row >= max_index) 

501 new_row = np.arange(max_index, dtype=idx_dtype) 

502 new_col = np.arange(k, k + max_index, dtype=idx_dtype) 

503 

504 # Define the array of data consisting of the entries to be added. 

505 if values.ndim: 

506 new_data = values[:max_index] 

507 else: 

508 new_data = np.empty(max_index, dtype=self.dtype) 

509 new_data[:] = values 

510 

511 # Update the internal structure. 

512 self.row = np.concatenate((self.row[keep], new_row)) 

513 self.col = np.concatenate((self.col[keep], new_col)) 

514 self.data = np.concatenate((self.data[keep], new_data)) 

515 self.has_canonical_format = False 

516 

517 # needed by _data_matrix 

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

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

520 but with different data. By default the index arrays 

521 (i.e. .row and .col) are copied. 

522 """ 

523 if copy: 

524 return self.__class__((data, (self.row.copy(), self.col.copy())), 

525 shape=self.shape, dtype=data.dtype) 

526 else: 

527 return self.__class__((data, (self.row, self.col)), 

528 shape=self.shape, dtype=data.dtype) 

529 

530 def sum_duplicates(self): 

531 """Eliminate duplicate matrix entries by adding them together 

532 

533 This is an *in place* operation 

534 """ 

535 if self.has_canonical_format: 

536 return 

537 summed = self._sum_duplicates(self.row, self.col, self.data) 

538 self.row, self.col, self.data = summed 

539 self.has_canonical_format = True 

540 

541 def _sum_duplicates(self, row, col, data): 

542 # Assumes (data, row, col) not in canonical format. 

543 if len(data) == 0: 

544 return row, col, data 

545 # Sort indices w.r.t. rows, then cols. This corresponds to C-order, 

546 # which we rely on for argmin/argmax to return the first index in the 

547 # same way that numpy does (in the case of ties). 

548 order = np.lexsort((col, row)) 

549 row = row[order] 

550 col = col[order] 

551 data = data[order] 

552 unique_mask = ((row[1:] != row[:-1]) | 

553 (col[1:] != col[:-1])) 

554 unique_mask = np.append(True, unique_mask) 

555 row = row[unique_mask] 

556 col = col[unique_mask] 

557 unique_inds, = np.nonzero(unique_mask) 

558 data = np.add.reduceat(data, unique_inds, dtype=self.dtype) 

559 return row, col, data 

560 

561 def eliminate_zeros(self): 

562 """Remove zero entries from the matrix 

563 

564 This is an *in place* operation 

565 """ 

566 mask = self.data != 0 

567 self.data = self.data[mask] 

568 self.row = self.row[mask] 

569 self.col = self.col[mask] 

570 

571 ####################### 

572 # Arithmetic handlers # 

573 ####################### 

574 

575 def _add_dense(self, other): 

576 if other.shape != self.shape: 

577 raise ValueError('Incompatible shapes ({} and {})' 

578 .format(self.shape, other.shape)) 

579 dtype = upcast_char(self.dtype.char, other.dtype.char) 

580 result = np.array(other, dtype=dtype, copy=True) 

581 fortran = int(result.flags.f_contiguous) 

582 M, N = self.shape 

583 coo_todense(M, N, self.nnz, self.row, self.col, self.data, 

584 result.ravel('A'), fortran) 

585 return self._container(result, copy=False) 

586 

587 def _mul_vector(self, other): 

588 #output array 

589 result = np.zeros(self.shape[0], dtype=upcast_char(self.dtype.char, 

590 other.dtype.char)) 

591 coo_matvec(self.nnz, self.row, self.col, self.data, other, result) 

592 return result 

593 

594 def _mul_multivector(self, other): 

595 result = np.zeros((other.shape[1], self.shape[0]), 

596 dtype=upcast_char(self.dtype.char, other.dtype.char)) 

597 for i, col in enumerate(other.T): 

598 coo_matvec(self.nnz, self.row, self.col, self.data, col, result[i]) 

599 return result.T.view(type=type(other)) 

600 

601 

602def isspmatrix_coo(x): 

603 """Is `x` of coo_matrix type? 

604 

605 Parameters 

606 ---------- 

607 x 

608 object to check for being a coo matrix 

609 

610 Returns 

611 ------- 

612 bool 

613 True if `x` is a coo matrix, False otherwise 

614 

615 Examples 

616 -------- 

617 >>> from scipy.sparse import coo_array, coo_matrix, csr_matrix, isspmatrix_coo 

618 >>> isspmatrix_coo(coo_matrix([[5]])) 

619 True 

620 >>> isspmatrix_coo(coo_array([[5]])) 

621 False 

622 >>> isspmatrix_coo(csr_matrix([[5]])) 

623 False 

624 """ 

625 return isinstance(x, coo_matrix) 

626 

627 

628# This namespace class separates array from matrix with isinstance 

629class coo_array(_coo_base, sparray): 

630 pass 

631 

632coo_array.__doc__ = _coo_base.__doc__ 

633 

634class coo_matrix(spmatrix, _coo_base): 

635 pass 

636 

637coo_matrix.__doc__ = _array_doc_to_matrix(_coo_base.__doc__)