Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.9/dist-packages/scipy/sparse/_construct.py: 10%

1"""Functions to construct sparse matrices and arrays 

2""" 

3 

4__docformat__ = "restructuredtext en" 

5 

6__all__ = ['spdiags', 'eye', 'identity', 'kron', 'kronsum', 

7 'hstack', 'vstack', 'bmat', 'rand', 'random', 'diags', 'block_diag', 

8 'diags_array', 'block'] 

9 

10import numbers 

11from functools import partial 

12import numpy as np 

13 

14from scipy._lib._util import check_random_state, rng_integers 

15from ._sputils import upcast, get_index_dtype, isscalarlike 

16 

17from ._sparsetools import csr_hstack 

18from ._bsr import bsr_matrix 

19from ._coo import coo_matrix, coo_array 

20from ._csc import csc_matrix, csc_array 

21from ._csr import csr_matrix, csr_array 

22from ._dia import dia_matrix, dia_array 

23 

24from ._base import issparse, sparray 

25 

26 

27def spdiags(data, diags, m=None, n=None, format=None): 

28 """ 

29 Return a sparse matrix from diagonals. 

30 

31 Parameters 

32 ---------- 

33 data : array_like 

34 Matrix diagonals stored row-wise 

35 diags : sequence of int or an int 

36 Diagonals to set: 

37 

38 * k = 0 the main diagonal 

39 * k > 0 the kth upper diagonal 

40 * k < 0 the kth lower diagonal 

41 m, n : int, tuple, optional 

42 Shape of the result. If `n` is None and `m` is a given tuple, 

43 the shape is this tuple. If omitted, the matrix is square and 

44 its shape is len(data[0]). 

45 format : str, optional 

46 Format of the result. By default (format=None) an appropriate sparse 

47 matrix format is returned. This choice is subject to change. 

48 

49 See Also 

50 -------- 

51 diags : more convenient form of this function 

52 dia_matrix : the sparse DIAgonal format. 

53 

54 Examples 

55 -------- 

56 >>> import numpy as np 

57 >>> from scipy.sparse import spdiags 

58 >>> data = np.array([[1, 2, 3, 4], [1, 2, 3, 4], [1, 2, 3, 4]]) 

59 >>> diags = np.array([0, -1, 2]) 

60 >>> spdiags(data, diags, 4, 4).toarray() 

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

62 [1, 2, 0, 4], 

63 [0, 2, 3, 0], 

64 [0, 0, 3, 4]]) 

65 

66 """ 

67 if m is None and n is None: 

68 m = n = len(data[0]) 

69 elif n is None: 

70 m, n = m 

71 return dia_matrix((data, diags), shape=(m, n)).asformat(format) 

72 

73 

74def diags_array(diagonals, /, *, offsets=0, shape=None, format=None, dtype=None): 

75 """ 

76 Construct a sparse array from diagonals. 

77 

78 Parameters 

79 ---------- 

80 diagonals : sequence of array_like 

81 Sequence of arrays containing the array diagonals, 

82 corresponding to `offsets`. 

83 offsets : sequence of int or an int, optional 

84 Diagonals to set: 

85 - k = 0 the main diagonal (default) 

86 - k > 0 the kth upper diagonal 

87 - k < 0 the kth lower diagonal 

88 shape : tuple of int, optional 

89 Shape of the result. If omitted, a square array large enough 

90 to contain the diagonals is returned. 

91 format : {"dia", "csr", "csc", "lil", ...}, optional 

92 Matrix format of the result. By default (format=None) an 

93 appropriate sparse array format is returned. This choice is 

94 subject to change. 

95 dtype : dtype, optional 

96 Data type of the array. 

97 

98 Notes 

99 ----- 

100 The result from `diags_array` is the sparse equivalent of:: 

101 

102 np.diag(diagonals[0], offsets[0]) 

103 + ... 

104 + np.diag(diagonals[k], offsets[k]) 

105 

106 Repeated diagonal offsets are disallowed. 

107 

108 .. versionadded:: 1.11 

109 

110 Examples 

111 -------- 

112 >>> from scipy.sparse import diags_array 

113 >>> diagonals = [[1, 2, 3, 4], [1, 2, 3], [1, 2]] 

114 >>> diags_array(diagonals, offsets=[0, -1, 2]).toarray() 

115 array([[1, 0, 1, 0], 

116 [1, 2, 0, 2], 

117 [0, 2, 3, 0], 

118 [0, 0, 3, 4]]) 

119 

120 Broadcasting of scalars is supported (but shape needs to be 

121 specified): 

122 

123 >>> diags_array([1, -2, 1], offsets=[-1, 0, 1], shape=(4, 4)).toarray() 

124 array([[-2., 1., 0., 0.], 

125 [ 1., -2., 1., 0.], 

126 [ 0., 1., -2., 1.], 

127 [ 0., 0., 1., -2.]]) 

128 

129 

130 If only one diagonal is wanted (as in `numpy.diag`), the following 

131 works as well: 

132 

133 >>> diags_array([1, 2, 3], offsets=1).toarray() 

134 array([[ 0., 1., 0., 0.], 

135 [ 0., 0., 2., 0.], 

136 [ 0., 0., 0., 3.], 

137 [ 0., 0., 0., 0.]]) 

138 """ 

139 # if offsets is not a sequence, assume that there's only one diagonal 

140 if isscalarlike(offsets): 

141 # now check that there's actually only one diagonal 

142 if len(diagonals) == 0 or isscalarlike(diagonals[0]): 

143 diagonals = [np.atleast_1d(diagonals)] 

144 else: 

145 raise ValueError("Different number of diagonals and offsets.") 

146 else: 

147 diagonals = list(map(np.atleast_1d, diagonals)) 

148 

149 offsets = np.atleast_1d(offsets) 

150 

151 # Basic check 

152 if len(diagonals) != len(offsets): 

153 raise ValueError("Different number of diagonals and offsets.") 

154 

155 # Determine shape, if omitted 

156 if shape is None: 

157 m = len(diagonals[0]) + abs(int(offsets[0])) 

158 shape = (m, m) 

159 

160 # Determine data type, if omitted 

161 if dtype is None: 

162 dtype = np.common_type(*diagonals) 

163 

164 # Construct data array 

165 m, n = shape 

166 

167 M = max([min(m + offset, n - offset) + max(0, offset) 

168 for offset in offsets]) 

169 M = max(0, M) 

170 data_arr = np.zeros((len(offsets), M), dtype=dtype) 

171 

172 K = min(m, n) 

173 

174 for j, diagonal in enumerate(diagonals): 

175 offset = offsets[j] 

176 k = max(0, offset) 

177 length = min(m + offset, n - offset, K) 

178 if length < 0: 

179 raise ValueError("Offset %d (index %d) out of bounds" % (offset, j)) 

180 try: 

181 data_arr[j, k:k+length] = diagonal[...,:length] 

182 except ValueError as e: 

183 if len(diagonal) != length and len(diagonal) != 1: 

184 raise ValueError( 

185 "Diagonal length (index %d: %d at offset %d) does not " 

186 "agree with array size (%d, %d)." % ( 

187 j, len(diagonal), offset, m, n)) from e 

188 raise 

189 

190 return dia_array((data_arr, offsets), shape=(m, n)).asformat(format) 

191 

192 

193def diags(diagonals, offsets=0, shape=None, format=None, dtype=None): 

194 """ 

195 Construct a sparse matrix from diagonals. 

196 

197 Parameters 

198 ---------- 

199 diagonals : sequence of array_like 

200 Sequence of arrays containing the matrix diagonals, 

201 corresponding to `offsets`. 

202 offsets : sequence of int or an int, optional 

203 Diagonals to set: 

204 - k = 0 the main diagonal (default) 

205 - k > 0 the kth upper diagonal 

206 - k < 0 the kth lower diagonal 

207 shape : tuple of int, optional 

208 Shape of the result. If omitted, a square matrix large enough 

209 to contain the diagonals is returned. 

210 format : {"dia", "csr", "csc", "lil", ...}, optional 

211 Matrix format of the result. By default (format=None) an 

212 appropriate sparse matrix format is returned. This choice is 

213 subject to change. 

214 dtype : dtype, optional 

215 Data type of the matrix. 

216 

217 See Also 

218 -------- 

219 spdiags : construct matrix from diagonals 

220 

221 Notes 

222 ----- 

223 This function differs from `spdiags` in the way it handles 

224 off-diagonals. 

225 

226 The result from `diags` is the sparse equivalent of:: 

227 

228 np.diag(diagonals[0], offsets[0]) 

229 + ... 

230 + np.diag(diagonals[k], offsets[k]) 

231 

232 Repeated diagonal offsets are disallowed. 

233 

234 .. versionadded:: 0.11 

235 

236 Examples 

237 -------- 

238 >>> from scipy.sparse import diags 

239 >>> diagonals = [[1, 2, 3, 4], [1, 2, 3], [1, 2]] 

240 >>> diags(diagonals, [0, -1, 2]).toarray() 

241 array([[1, 0, 1, 0], 

242 [1, 2, 0, 2], 

243 [0, 2, 3, 0], 

244 [0, 0, 3, 4]]) 

245 

246 Broadcasting of scalars is supported (but shape needs to be 

247 specified): 

248 

249 >>> diags([1, -2, 1], [-1, 0, 1], shape=(4, 4)).toarray() 

250 array([[-2., 1., 0., 0.], 

251 [ 1., -2., 1., 0.], 

252 [ 0., 1., -2., 1.], 

253 [ 0., 0., 1., -2.]]) 

254 

255 

256 If only one diagonal is wanted (as in `numpy.diag`), the following 

257 works as well: 

258 

259 >>> diags([1, 2, 3], 1).toarray() 

260 array([[ 0., 1., 0., 0.], 

261 [ 0., 0., 2., 0.], 

262 [ 0., 0., 0., 3.], 

263 [ 0., 0., 0., 0.]]) 

264 """ 

265 A = diags_array(diagonals, offsets=offsets, shape=shape, dtype=dtype) 

266 return dia_matrix(A).asformat(format) 

267 

268 

269def identity(n, dtype='d', format=None): 

270 """Identity matrix in sparse format 

271 

272 Returns an identity matrix with shape (n,n) using a given 

273 sparse format and dtype. 

274 

275 Parameters 

276 ---------- 

277 n : int 

278 Shape of the identity matrix. 

279 dtype : dtype, optional 

280 Data type of the matrix 

281 format : str, optional 

282 Sparse format of the result, e.g., format="csr", etc. 

283 

284 Examples 

285 -------- 

286 >>> from scipy.sparse import identity 

287 >>> identity(3).toarray() 

288 array([[ 1., 0., 0.], 

289 [ 0., 1., 0.], 

290 [ 0., 0., 1.]]) 

291 >>> identity(3, dtype='int8', format='dia') 

292 <3x3 sparse matrix of type '<class 'numpy.int8'>' 

293 with 3 stored elements (1 diagonals) in DIAgonal format> 

294 

295 """ 

296 return eye(n, n, dtype=dtype, format=format) 

297 

298 

299def eye(m, n=None, k=0, dtype=float, format=None): 

300 """Sparse matrix with ones on diagonal 

301 

302 Returns a sparse (m x n) matrix where the kth diagonal 

303 is all ones and everything else is zeros. 

304 

305 Parameters 

306 ---------- 

307 m : int 

308 Number of rows in the matrix. 

309 n : int, optional 

310 Number of columns. Default: `m`. 

311 k : int, optional 

312 Diagonal to place ones on. Default: 0 (main diagonal). 

313 dtype : dtype, optional 

314 Data type of the matrix. 

315 format : str, optional 

316 Sparse format of the result, e.g., format="csr", etc. 

317 

318 Examples 

319 -------- 

320 >>> import numpy as np 

321 >>> from scipy import sparse 

322 >>> sparse.eye(3).toarray() 

323 array([[ 1., 0., 0.], 

324 [ 0., 1., 0.], 

325 [ 0., 0., 1.]]) 

326 >>> sparse.eye(3, dtype=np.int8) 

327 <3x3 sparse matrix of type '<class 'numpy.int8'>' 

328 with 3 stored elements (1 diagonals) in DIAgonal format> 

329 

330 """ 

331 if n is None: 

332 n = m 

333 m,n = int(m),int(n) 

334 

335 if m == n and k == 0: 

336 # fast branch for special formats 

337 if format in ['csr', 'csc']: 

338 idx_dtype = get_index_dtype(maxval=n) 

339 indptr = np.arange(n+1, dtype=idx_dtype) 

340 indices = np.arange(n, dtype=idx_dtype) 

341 data = np.ones(n, dtype=dtype) 

342 cls = {'csr': csr_matrix, 'csc': csc_matrix}[format] 

343 return cls((data,indices,indptr),(n,n)) 

344 elif format == 'coo': 

345 idx_dtype = get_index_dtype(maxval=n) 

346 row = np.arange(n, dtype=idx_dtype) 

347 col = np.arange(n, dtype=idx_dtype) 

348 data = np.ones(n, dtype=dtype) 

349 return coo_matrix((data, (row, col)), (n, n)) 

350 

351 data = np.ones((1, max(0, min(m + k, n))), dtype=dtype) 

352 return diags(data, offsets=[k], shape=(m, n), dtype=dtype).asformat(format) 

353 

354 

355def kron(A, B, format=None): 

356 """kronecker product of sparse matrices A and B 

357 

358 Parameters 

359 ---------- 

360 A : sparse or dense matrix 

361 first matrix of the product 

362 B : sparse or dense matrix 

363 second matrix of the product 

364 format : str, optional 

365 format of the result (e.g. "csr") 

366 

367 Returns 

368 ------- 

369 kronecker product in a sparse matrix format 

370 

371 

372 Examples 

373 -------- 

374 >>> import numpy as np 

375 >>> from scipy import sparse 

376 >>> A = sparse.csr_matrix(np.array([[0, 2], [5, 0]])) 

377 >>> B = sparse.csr_matrix(np.array([[1, 2], [3, 4]])) 

378 >>> sparse.kron(A, B).toarray() 

379 array([[ 0, 0, 2, 4], 

380 [ 0, 0, 6, 8], 

381 [ 5, 10, 0, 0], 

382 [15, 20, 0, 0]]) 

383 

384 >>> sparse.kron(A, [[1, 2], [3, 4]]).toarray() 

385 array([[ 0, 0, 2, 4], 

386 [ 0, 0, 6, 8], 

387 [ 5, 10, 0, 0], 

388 [15, 20, 0, 0]]) 

389 

390 """ 

391 B = coo_matrix(B) 

392 

393 if (format is None or format == "bsr") and 2*B.nnz >= B.shape[0] * B.shape[1]: 

394 # B is fairly dense, use BSR 

395 A = csr_matrix(A,copy=True) 

396 output_shape = (A.shape[0]*B.shape[0], A.shape[1]*B.shape[1]) 

397 

398 if A.nnz == 0 or B.nnz == 0: 

399 # kronecker product is the zero matrix 

400 return coo_matrix(output_shape).asformat(format) 

401 

402 B = B.toarray() 

403 data = A.data.repeat(B.size).reshape(-1,B.shape[0],B.shape[1]) 

404 data = data * B 

405 

406 return bsr_matrix((data,A.indices,A.indptr), shape=output_shape) 

407 else: 

408 # use COO 

409 A = coo_matrix(A) 

410 output_shape = (A.shape[0]*B.shape[0], A.shape[1]*B.shape[1]) 

411 

412 if A.nnz == 0 or B.nnz == 0: 

413 # kronecker product is the zero matrix 

414 return coo_matrix(output_shape).asformat(format) 

415 

416 # expand entries of a into blocks 

417 row = A.row.repeat(B.nnz) 

418 col = A.col.repeat(B.nnz) 

419 data = A.data.repeat(B.nnz) 

420 

421 if max(A.shape[0]*B.shape[0], A.shape[1]*B.shape[1]) > np.iinfo('int32').max: 

422 row = row.astype(np.int64) 

423 col = col.astype(np.int64) 

424 

425 row *= B.shape[0] 

426 col *= B.shape[1] 

427 

428 # increment block indices 

429 row,col = row.reshape(-1,B.nnz),col.reshape(-1,B.nnz) 

430 row += B.row 

431 col += B.col 

432 row,col = row.reshape(-1),col.reshape(-1) 

433 

434 # compute block entries 

435 data = data.reshape(-1,B.nnz) * B.data 

436 data = data.reshape(-1) 

437 

438 return coo_matrix((data,(row,col)), shape=output_shape).asformat(format) 

439 

440 

441def kronsum(A, B, format=None): 

442 """kronecker sum of sparse matrices A and B 

443 

444 Kronecker sum of two sparse matrices is a sum of two Kronecker 

445 products kron(I_n,A) + kron(B,I_m) where A has shape (m,m) 

446 and B has shape (n,n) and I_m and I_n are identity matrices 

447 of shape (m,m) and (n,n), respectively. 

448 

449 Parameters 

450 ---------- 

451 A 

452 square matrix 

453 B 

454 square matrix 

455 format : str 

456 format of the result (e.g. "csr") 

457 

458 Returns 

459 ------- 

460 kronecker sum in a sparse matrix format 

461 

462 Examples 

463 -------- 

464 

465 

466 """ 

467 A = coo_matrix(A) 

468 B = coo_matrix(B) 

469 

470 if A.shape[0] != A.shape[1]: 

471 raise ValueError('A is not square') 

472 

473 if B.shape[0] != B.shape[1]: 

474 raise ValueError('B is not square') 

475 

476 dtype = upcast(A.dtype, B.dtype) 

477 

478 L = kron(eye(B.shape[0],dtype=dtype), A, format=format) 

479 R = kron(B, eye(A.shape[0],dtype=dtype), format=format) 

480 

481 return (L+R).asformat(format) # since L + R is not always same format 

482 

483 

484def _compressed_sparse_stack(blocks, axis, return_spmatrix): 

485 """ 

486 Stacking fast path for CSR/CSC matrices or arrays 

487 (i) vstack for CSR, (ii) hstack for CSC. 

488 """ 

489 other_axis = 1 if axis == 0 else 0 

490 data = np.concatenate([b.data for b in blocks]) 

491 constant_dim = blocks[0].shape[other_axis] 

492 idx_dtype = get_index_dtype(arrays=[b.indptr for b in blocks], 

493 maxval=max(data.size, constant_dim)) 

494 indices = np.empty(data.size, dtype=idx_dtype) 

495 indptr = np.empty(sum(b.shape[axis] for b in blocks) + 1, dtype=idx_dtype) 

496 last_indptr = idx_dtype(0) 

497 sum_dim = 0 

498 sum_indices = 0 

499 for b in blocks: 

500 if b.shape[other_axis] != constant_dim: 

501 raise ValueError(f'incompatible dimensions for axis {other_axis}') 

502 indices[sum_indices:sum_indices+b.indices.size] = b.indices 

503 sum_indices += b.indices.size 

504 idxs = slice(sum_dim, sum_dim + b.shape[axis]) 

505 indptr[idxs] = b.indptr[:-1] 

506 indptr[idxs] += last_indptr 

507 sum_dim += b.shape[axis] 

508 last_indptr += b.indptr[-1] 

509 indptr[-1] = last_indptr 

510 # TODO remove this if-structure when sparse matrices removed 

511 if return_spmatrix: 

512 if axis == 0: 

513 return csr_matrix((data, indices, indptr), 

514 shape=(sum_dim, constant_dim)) 

515 else: 

516 return csc_matrix((data, indices, indptr), 

517 shape=(constant_dim, sum_dim)) 

518 

519 if axis == 0: 

520 return csr_array((data, indices, indptr), 

521 shape=(sum_dim, constant_dim)) 

522 else: 

523 return csc_array((data, indices, indptr), 

524 shape=(constant_dim, sum_dim)) 

525 

526 

527def _stack_along_minor_axis(blocks, axis): 

528 """ 

529 Stacking fast path for CSR/CSC matrices along the minor axis 

530 (i) hstack for CSR, (ii) vstack for CSC. 

531 """ 

532 n_blocks = len(blocks) 

533 if n_blocks == 0: 

534 raise ValueError('Missing block matrices') 

535 

536 if n_blocks == 1: 

537 return blocks[0] 

538 

539 # check for incompatible dimensions 

540 other_axis = 1 if axis == 0 else 0 

541 other_axis_dims = {b.shape[other_axis] for b in blocks} 

542 if len(other_axis_dims) > 1: 

543 raise ValueError(f'Mismatching dimensions along axis {other_axis}: ' 

544 f'{other_axis_dims}') 

545 constant_dim, = other_axis_dims 

546 

547 # Do the stacking 

548 indptr_list = [b.indptr for b in blocks] 

549 data_cat = np.concatenate([b.data for b in blocks]) 

550 

551 # Need to check if any indices/indptr, would be too large post- 

552 # concatenation for np.int32: 

553 # - The max value of indices is the output array's stacking-axis length - 1 

554 # - The max value in indptr is the number of non-zero entries. This is 

555 # exceedingly unlikely to require int64, but is checked out of an 

556 # abundance of caution. 

557 sum_dim = sum(b.shape[axis] for b in blocks) 

558 nnz = sum(len(b.indices) for b in blocks) 

559 idx_dtype = get_index_dtype(maxval=max(sum_dim - 1, nnz)) 

560 stack_dim_cat = np.array([b.shape[axis] for b in blocks], dtype=idx_dtype) 

561 if data_cat.size > 0: 

562 indptr_cat = np.concatenate(indptr_list).astype(idx_dtype) 

563 indices_cat = (np.concatenate([b.indices for b in blocks]) 

564 .astype(idx_dtype)) 

565 indptr = np.empty(constant_dim + 1, dtype=idx_dtype) 

566 indices = np.empty_like(indices_cat) 

567 data = np.empty_like(data_cat) 

568 csr_hstack(n_blocks, constant_dim, stack_dim_cat, 

569 indptr_cat, indices_cat, data_cat, 

570 indptr, indices, data) 

571 else: 

572 indptr = np.zeros(constant_dim + 1, dtype=idx_dtype) 

573 indices = np.empty(0, dtype=idx_dtype) 

574 data = np.empty(0, dtype=data_cat.dtype) 

575 

576 if axis == 0: 

577 return blocks[0]._csc_container((data, indices, indptr), 

578 shape=(sum_dim, constant_dim)) 

579 else: 

580 return blocks[0]._csr_container((data, indices, indptr), 

581 shape=(constant_dim, sum_dim)) 

582 

583 

584def hstack(blocks, format=None, dtype=None): 

585 """ 

586 Stack sparse matrices horizontally (column wise) 

587 

588 Parameters 

589 ---------- 

590 blocks 

591 sequence of sparse matrices with compatible shapes 

592 format : str 

593 sparse format of the result (e.g., "csr") 

594 by default an appropriate sparse matrix format is returned. 

595 This choice is subject to change. 

596 dtype : dtype, optional 

597 The data-type of the output matrix. If not given, the dtype is 

598 determined from that of `blocks`. 

599 

600 Returns 

601 ------- 

602 new_array : sparse matrix or array 

603 If any block in blocks is a sparse array, return a sparse array. 

604 Otherwise return a sparse matrix. 

605 

606 If you want a sparse array built from blocks that are not sparse 

607 arrays, use `block(hstack(blocks))` or convert one block 

608 e.g. `blocks[0] = csr_array(blocks[0])`. 

609 

610 See Also 

611 -------- 

612 vstack : stack sparse matrices vertically (row wise) 

613 

614 Examples 

615 -------- 

616 >>> from scipy.sparse import coo_matrix, hstack 

617 >>> A = coo_matrix([[1, 2], [3, 4]]) 

618 >>> B = coo_matrix([[5], [6]]) 

619 >>> hstack([A,B]).toarray() 

620 array([[1, 2, 5], 

621 [3, 4, 6]]) 

622 

623 """ 

624 blocks = np.asarray(blocks, dtype='object') 

625 if any(isinstance(b, sparray) for b in blocks.flat): 

626 return _block([blocks], format, dtype) 

627 else: 

628 return _block([blocks], format, dtype, return_spmatrix=True) 

629 

630 

631def vstack(blocks, format=None, dtype=None): 

632 """ 

633 Stack sparse arrays vertically (row wise) 

634 

635 Parameters 

636 ---------- 

637 blocks 

638 sequence of sparse arrays with compatible shapes 

639 format : str, optional 

640 sparse format of the result (e.g., "csr") 

641 by default an appropriate sparse array format is returned. 

642 This choice is subject to change. 

643 dtype : dtype, optional 

644 The data-type of the output array. If not given, the dtype is 

645 determined from that of `blocks`. 

646 

647 Returns 

648 ------- 

649 new_array : sparse matrix or array 

650 If any block in blocks is a sparse array, return a sparse array. 

651 Otherwise return a sparse matrix. 

652 

653 If you want a sparse array built from blocks that are not sparse 

654 arrays, use `block(vstack(blocks))` or convert one block 

655 e.g. `blocks[0] = csr_array(blocks[0])`. 

656 

657 See Also 

658 -------- 

659 hstack : stack sparse matrices horizontally (column wise) 

660 

661 Examples 

662 -------- 

663 >>> from scipy.sparse import coo_array, vstack 

664 >>> A = coo_array([[1, 2], [3, 4]]) 

665 >>> B = coo_array([[5, 6]]) 

666 >>> vstack([A, B]).toarray() 

667 array([[1, 2], 

668 [3, 4], 

669 [5, 6]]) 

670 

671 """ 

672 blocks = np.asarray(blocks, dtype='object') 

673 if any(isinstance(b, sparray) for b in blocks.flat): 

674 return _block([[b] for b in blocks], format, dtype) 

675 else: 

676 return _block([[b] for b in blocks], format, dtype, return_spmatrix=True) 

677 

678 

679def bmat(blocks, format=None, dtype=None): 

680 """ 

681 Build a sparse array or matrix from sparse sub-blocks 

682 

683 Note: `block` is preferred over `bmat`. They are the same function 

684 except that `bmat` can return a deprecated sparse matrix. 

685 `bmat` returns a coo_matrix if none of the inputs are a sparse array. 

686 

687 Parameters 

688 ---------- 

689 blocks : array_like 

690 Grid of sparse matrices with compatible shapes. 

691 An entry of None implies an all-zero matrix. 

692 format : {'bsr', 'coo', 'csc', 'csr', 'dia', 'dok', 'lil'}, optional 

693 The sparse format of the result (e.g. "csr"). By default an 

694 appropriate sparse matrix format is returned. 

695 This choice is subject to change. 

696 dtype : dtype, optional 

697 The data-type of the output matrix. If not given, the dtype is 

698 determined from that of `blocks`. 

699 

700 Returns 

701 ------- 

702 bmat : sparse matrix or array 

703 If any block in blocks is a sparse array, return a sparse array. 

704 Otherwise return a sparse matrix. 

705 

706 If you want a sparse array built from blocks that are not sparse 

707 arrays, use `block()`. 

708 

709 See Also 

710 -------- 

711 block 

712 

713 Examples 

714 -------- 

715 >>> from scipy.sparse import coo_array, bmat 

716 >>> A = coo_array([[1, 2], [3, 4]]) 

717 >>> B = coo_array([[5], [6]]) 

718 >>> C = coo_array([[7]]) 

719 >>> bmat([[A, B], [None, C]]).toarray() 

720 array([[1, 2, 5], 

721 [3, 4, 6], 

722 [0, 0, 7]]) 

723 

724 >>> bmat([[A, None], [None, C]]).toarray() 

725 array([[1, 2, 0], 

726 [3, 4, 0], 

727 [0, 0, 7]]) 

728 

729 """ 

730 blocks = np.asarray(blocks, dtype='object') 

731 if any(isinstance(b, sparray) for b in blocks.flat): 

732 return _block(blocks, format, dtype) 

733 else: 

734 return _block(blocks, format, dtype, return_spmatrix=True) 

735 

736 

737def block(blocks, *, format=None, dtype=None): 

738 """ 

739 Build a sparse array from sparse sub-blocks 

740 

741 Parameters 

742 ---------- 

743 blocks : array_like 

744 Grid of sparse arrays with compatible shapes. 

745 An entry of None implies an all-zero array. 

746 format : {'bsr', 'coo', 'csc', 'csr', 'dia', 'dok', 'lil'}, optional 

747 The sparse format of the result (e.g. "csr"). By default an 

748 appropriate sparse array format is returned. 

749 This choice is subject to change. 

750 dtype : dtype, optional 

751 The data-type of the output array. If not given, the dtype is 

752 determined from that of `blocks`. 

753 

754 Returns 

755 ------- 

756 block : sparse array 

757 

758 See Also 

759 -------- 

760 block_diag : specify blocks along the main diagonals 

761 diags : specify (possibly offset) diagonals 

762 

763 Examples 

764 -------- 

765 >>> from scipy.sparse import coo_array, block 

766 >>> A = coo_array([[1, 2], [3, 4]]) 

767 >>> B = coo_array([[5], [6]]) 

768 >>> C = coo_array([[7]]) 

769 >>> block([[A, B], [None, C]]).toarray() 

770 array([[1, 2, 5], 

771 [3, 4, 6], 

772 [0, 0, 7]]) 

773 

774 >>> block([[A, None], [None, C]]).toarray() 

775 array([[1, 2, 0], 

776 [3, 4, 0], 

777 [0, 0, 7]]) 

778 

779 """ 

780 return _block(blocks, format, dtype) 

781 

782 

783def _block(blocks, format, dtype, return_spmatrix=False): 

784 blocks = np.asarray(blocks, dtype='object') 

785 

786 if blocks.ndim != 2: 

787 raise ValueError('blocks must be 2-D') 

788 

789 M,N = blocks.shape 

790 

791 # check for fast path cases 

792 if (format in (None, 'csr') and 

793 all(issparse(b) and b.format == 'csr' for b in blocks.flat) 

794 ): 

795 if N > 1: 

796 # stack along columns (axis 1): must have shape (M, 1) 

797 blocks = [[_stack_along_minor_axis(blocks[b, :], 1)] for b in range(M)] 

798 blocks = np.asarray(blocks, dtype='object') 

799 

800 # stack along rows (axis 0): 

801 A = _compressed_sparse_stack(blocks[:, 0], 0, return_spmatrix) 

802 if dtype is not None: 

803 A = A.astype(dtype) 

804 return A 

805 elif (format in (None, 'csc') and 

806 all(issparse(b) and b.format == 'csc' for b in blocks.flat) 

807 ): 

808 if M > 1: 

809 # stack along rows (axis 0): must have shape (1, N) 

810 blocks = [[_stack_along_minor_axis(blocks[:, b], 0) for b in range(N)]] 

811 blocks = np.asarray(blocks, dtype='object') 

812 

813 # stack along columns (axis 1): 

814 A = _compressed_sparse_stack(blocks[0, :], 1, return_spmatrix) 

815 if dtype is not None: 

816 A = A.astype(dtype) 

817 return A 

818 

819 block_mask = np.zeros(blocks.shape, dtype=bool) 

820 brow_lengths = np.zeros(M, dtype=np.int64) 

821 bcol_lengths = np.zeros(N, dtype=np.int64) 

822 

823 # convert everything to COO format 

824 for i in range(M): 

825 for j in range(N): 

826 if blocks[i,j] is not None: 

827 A = coo_array(blocks[i,j]) 

828 blocks[i,j] = A 

829 block_mask[i,j] = True 

830 

831 if brow_lengths[i] == 0: 

832 brow_lengths[i] = A.shape[0] 

833 elif brow_lengths[i] != A.shape[0]: 

834 msg = (f'blocks[{i},:] has incompatible row dimensions. ' 

835 f'Got blocks[{i},{j}].shape[0] == {A.shape[0]}, ' 

836 f'expected {brow_lengths[i]}.') 

837 raise ValueError(msg) 

838 

839 if bcol_lengths[j] == 0: 

840 bcol_lengths[j] = A.shape[1] 

841 elif bcol_lengths[j] != A.shape[1]: 

842 msg = (f'blocks[:,{j}] has incompatible column ' 

843 f'dimensions. ' 

844 f'Got blocks[{i},{j}].shape[1] == {A.shape[1]}, ' 

845 f'expected {bcol_lengths[j]}.') 

846 raise ValueError(msg) 

847 

848 nnz = sum(block.nnz for block in blocks[block_mask]) 

849 if dtype is None: 

850 all_dtypes = [blk.dtype for blk in blocks[block_mask]] 

851 dtype = upcast(*all_dtypes) if all_dtypes else None 

852 

853 row_offsets = np.append(0, np.cumsum(brow_lengths)) 

854 col_offsets = np.append(0, np.cumsum(bcol_lengths)) 

855 

856 shape = (row_offsets[-1], col_offsets[-1]) 

857 

858 data = np.empty(nnz, dtype=dtype) 

859 idx_dtype = get_index_dtype(maxval=max(shape)) 

860 row = np.empty(nnz, dtype=idx_dtype) 

861 col = np.empty(nnz, dtype=idx_dtype) 

862 

863 nnz = 0 

864 ii, jj = np.nonzero(block_mask) 

865 for i, j in zip(ii, jj): 

866 B = blocks[i, j] 

867 idx = slice(nnz, nnz + B.nnz) 

868 data[idx] = B.data 

869 np.add(B.row, row_offsets[i], out=row[idx], dtype=idx_dtype) 

870 np.add(B.col, col_offsets[j], out=col[idx], dtype=idx_dtype) 

871 nnz += B.nnz 

872 

873 if return_spmatrix: 

874 return coo_matrix((data, (row, col)), shape=shape).asformat(format) 

875 return coo_array((data, (row, col)), shape=shape).asformat(format) 

876 

877 

878def block_diag(mats, format=None, dtype=None): 

879 """ 

880 Build a block diagonal sparse matrix or array from provided matrices. 

881 

882 Parameters 

883 ---------- 

884 mats : sequence of matrices or arrays 

885 Input matrices or arrays. 

886 format : str, optional 

887 The sparse format of the result (e.g., "csr"). If not given, the result