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

1"""Functions to construct sparse matrices 

2""" 

3 

4__docformat__ = "restructuredtext en" 

5 

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

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

8 

9import numbers 

10from functools import partial 

11import numpy as np 

12 

13from scipy._lib._util import check_random_state, rng_integers 

14from ._sputils import upcast, get_index_dtype, isscalarlike 

15 

16from ._sparsetools import csr_hstack 

17from ._csr import csr_matrix 

18from ._csc import csc_matrix 

19from ._bsr import bsr_matrix 

20from ._coo import coo_matrix 

21from ._dia import dia_matrix 

22 

23from ._base import issparse 

24 

25 

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

27 """ 

28 Return a sparse matrix from diagonals. 

29 

30 Parameters 

31 ---------- 

32 data : array_like 

33 Matrix diagonals stored row-wise 

34 diags : sequence of int or an int 

35 Diagonals to set: 

36 

37 * k = 0 the main diagonal 

38 * k > 0 the kth upper diagonal 

39 * k < 0 the kth lower diagonal 

40 m, n : int, tuple, optional 

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

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

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

44 format : str, optional 

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

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

47 

48 See Also 

49 -------- 

50 diags : more convenient form of this function 

51 dia_matrix : the sparse DIAgonal format. 

52 

53 Examples 

54 -------- 

55 >>> import numpy as np 

56 >>> from scipy.sparse import spdiags 

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

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

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

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

61 [1, 2, 0, 4], 

62 [0, 2, 3, 0], 

63 [0, 0, 3, 4]]) 

64 

65 """ 

66 if m is None and n is None: 

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

68 elif n is None: 

69 m, n = m 

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

71 

72 

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

74 """ 

75 Construct a sparse matrix from diagonals. 

76 

77 Parameters 

78 ---------- 

79 diagonals : sequence of array_like 

80 Sequence of arrays containing the matrix diagonals, 

81 corresponding to `offsets`. 

82 offsets : sequence of int or an int, optional 

83 Diagonals to set: 

84 - k = 0 the main diagonal (default) 

85 - k > 0 the kth upper diagonal 

86 - k < 0 the kth lower diagonal 

87 shape : tuple of int, optional 

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

89 to contain the diagonals is returned. 

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

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

92 appropriate sparse matrix format is returned. This choice is 

93 subject to change. 

94 dtype : dtype, optional 

95 Data type of the matrix. 

96 

97 See Also 

98 -------- 

99 spdiags : construct matrix from diagonals 

100 

101 Notes 

102 ----- 

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

104 off-diagonals. 

105 

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

107 

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

109 + ... 

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

111 

112 Repeated diagonal offsets are disallowed. 

113 

114 .. versionadded:: 0.11 

115 

116 Examples 

117 -------- 

118 >>> from scipy.sparse import diags 

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

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

121 array([[1, 0, 1, 0], 

122 [1, 2, 0, 2], 

123 [0, 2, 3, 0], 

124 [0, 0, 3, 4]]) 

125 

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

127 specified): 

128 

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

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

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

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

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

134 

135 

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

137 works as well: 

138 

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

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

141 [ 0., 0., 2., 0.], 

142 [ 0., 0., 0., 3.], 

143 [ 0., 0., 0., 0.]]) 

144 """ 

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

146 if isscalarlike(offsets): 

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

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

149 diagonals = [np.atleast_1d(diagonals)] 

150 else: 

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

152 else: 

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

154 

155 offsets = np.atleast_1d(offsets) 

156 

157 # Basic check 

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

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

160 

161 # Determine shape, if omitted 

162 if shape is None: 

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

164 shape = (m, m) 

165 

166 # Determine data type, if omitted 

167 if dtype is None: 

168 dtype = np.common_type(*diagonals) 

169 

170 # Construct data array 

171 m, n = shape 

172 

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

174 for offset in offsets]) 

175 M = max(0, M) 

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

177 

178 K = min(m, n) 

179 

180 for j, diagonal in enumerate(diagonals): 

181 offset = offsets[j] 

182 k = max(0, offset) 

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

184 if length < 0: 

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

186 try: 

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

188 except ValueError as e: 

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

190 raise ValueError( 

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

192 "agree with matrix size (%d, %d)." % ( 

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

194 raise 

195 

196 return dia_matrix((data_arr, offsets), shape=(m, n)).asformat(format) 

197 

198 

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

200 """Identity matrix in sparse format 

201 

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

203 sparse format and dtype. 

204 

205 Parameters 

206 ---------- 

207 n : int 

208 Shape of the identity matrix. 

209 dtype : dtype, optional 

210 Data type of the matrix 

211 format : str, optional 

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

213 

214 Examples 

215 -------- 

216 >>> from scipy.sparse import identity 

217 >>> identity(3).toarray() 

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

219 [ 0., 1., 0.], 

220 [ 0., 0., 1.]]) 

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

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

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

224 

225 """ 

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

227 

228 

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

230 """Sparse matrix with ones on diagonal 

231 

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

233 is all ones and everything else is zeros. 

234 

235 Parameters 

236 ---------- 

237 m : int 

238 Number of rows in the matrix. 

239 n : int, optional 

240 Number of columns. Default: `m`. 

241 k : int, optional 

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

243 dtype : dtype, optional 

244 Data type of the matrix. 

245 format : str, optional 

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

247 

248 Examples 

249 -------- 

250 >>> import numpy as np 

251 >>> from scipy import sparse 

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

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

254 [ 0., 1., 0.], 

255 [ 0., 0., 1.]]) 

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

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

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

259 

260 """ 

261 if n is None: 

262 n = m 

263 m,n = int(m),int(n) 

264 

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

266 # fast branch for special formats 

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

268 idx_dtype = get_index_dtype(maxval=n) 

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

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

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

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

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

274 elif format == 'coo': 

275 idx_dtype = get_index_dtype(maxval=n) 

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

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

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

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

280 

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

282 return spdiags(diags, k, m, n).asformat(format) 

283 

284 

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

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

287 

288 Parameters 

289 ---------- 

290 A : sparse or dense matrix 

291 first matrix of the product 

292 B : sparse or dense matrix 

293 second matrix of the product 

294 format : str, optional 

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

296 

297 Returns 

298 ------- 

299 kronecker product in a sparse matrix format 

300 

301 

302 Examples 

303 -------- 

304 >>> import numpy as np 

305 >>> from scipy import sparse 

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

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

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

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

310 [ 0, 0, 6, 8], 

311 [ 5, 10, 0, 0], 

312 [15, 20, 0, 0]]) 

313 

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

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

316 [ 0, 0, 6, 8], 

317 [ 5, 10, 0, 0], 

318 [15, 20, 0, 0]]) 

319 

320 """ 

321 B = coo_matrix(B) 

322 

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

324 # B is fairly dense, use BSR 

325 A = csr_matrix(A,copy=True) 

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

327 

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

329 # kronecker product is the zero matrix 

330 return coo_matrix(output_shape).asformat(format) 

331 

332 B = B.toarray() 

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

334 data = data * B 

335 

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

337 else: 

338 # use COO 

339 A = coo_matrix(A) 

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

341 

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

343 # kronecker product is the zero matrix 

344 return coo_matrix(output_shape).asformat(format) 

345 

346 # expand entries of a into blocks 

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

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

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

350 

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

352 row = row.astype(np.int64) 

353 col = col.astype(np.int64) 

354 

355 row *= B.shape[0] 

356 col *= B.shape[1] 

357 

358 # increment block indices 

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

360 row += B.row 

361 col += B.col 

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

363 

364 # compute block entries 

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

366 data = data.reshape(-1) 

367 

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

369 

370 

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

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

373 

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

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

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

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

378 

379 Parameters 

380 ---------- 

381 A 

382 square matrix 

383 B 

384 square matrix 

385 format : str 

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

387 

388 Returns 

389 ------- 

390 kronecker sum in a sparse matrix format 

391 

392 Examples 

393 -------- 

394 

395 

396 """ 

397 A = coo_matrix(A) 

398 B = coo_matrix(B) 

399 

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

401 raise ValueError('A is not square') 

402 

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

404 raise ValueError('B is not square') 

405 

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

407 

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

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

410 

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

412 

413 

414def _compressed_sparse_stack(blocks, axis): 

415 """ 

416 Stacking fast path for CSR/CSC matrices 

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

418 """ 

419 other_axis = 1 if axis == 0 else 0 

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

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

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

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

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

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

426 last_indptr = idx_dtype(0) 

427 sum_dim = 0 

428 sum_indices = 0 

429 for b in blocks: 

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

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

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

433 sum_indices += b.indices.size 

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

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

436 indptr[idxs] += last_indptr 

437 sum_dim += b.shape[axis] 

438 last_indptr += b.indptr[-1] 

439 indptr[-1] = last_indptr 

440 if axis == 0: 

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

442 shape=(sum_dim, constant_dim)) 

443 else: 

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

445 shape=(constant_dim, sum_dim)) 

446 

447 

448def _stack_along_minor_axis(blocks, axis): 

449 """ 

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

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

452 """ 

453 n_blocks = len(blocks) 

454 if n_blocks == 0: 

455 raise ValueError('Missing block matrices') 

456 

457 if n_blocks == 1: 

458 return blocks[0] 

459 

460 # check for incompatible dimensions 

461 other_axis = 1 if axis == 0 else 0 

462 other_axis_dims = set(b.shape[other_axis] for b in blocks) 

463 if len(other_axis_dims) > 1: 

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

465 f'{other_axis_dims}') 

466 constant_dim, = other_axis_dims 

467 

468 # Do the stacking 

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

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

471 

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

473 # concatenation for np.int32: 

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

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

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

477 # abundance of caution. 

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

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

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

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

482 if data_cat.size > 0: 

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

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

485 .astype(idx_dtype)) 

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

487 indices = np.empty_like(indices_cat) 

488 data = np.empty_like(data_cat) 

489 csr_hstack(n_blocks, constant_dim, stack_dim_cat, 

490 indptr_cat, indices_cat, data_cat, 

491 indptr, indices, data) 

492 else: 

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

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

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

496 

497 if axis == 0: 

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

499 shape=(sum_dim, constant_dim)) 

500 else: 

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

502 shape=(constant_dim, sum_dim)) 

503 

504 

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

506 """ 

507 Stack sparse matrices horizontally (column wise) 

508 

509 Parameters 

510 ---------- 

511 blocks 

512 sequence of sparse matrices with compatible shapes 

513 format : str 

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

515 by default an appropriate sparse matrix format is returned. 

516 This choice is subject to change. 

517 dtype : dtype, optional 

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

519 determined from that of `blocks`. 

520 

521 See Also 

522 -------- 

523 vstack : stack sparse matrices vertically (row wise) 

524 

525 Examples 

526 -------- 

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

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

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

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

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

532 [3, 4, 6]]) 

533 

534 """ 

535 return bmat([blocks], format=format, dtype=dtype) 

536 

537 

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

539 """ 

540 Stack sparse matrices vertically (row wise) 

541 

542 Parameters 

543 ---------- 

544 blocks 

545 sequence of sparse matrices with compatible shapes 

546 format : str, optional 

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

548 by default an appropriate sparse matrix format is returned. 

549 This choice is subject to change. 

550 dtype : dtype, optional 

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

552 determined from that of `blocks`. 

553 

554 See Also 

555 -------- 

556 hstack : stack sparse matrices horizontally (column wise) 

557 

558 Examples 

559 -------- 

560 >>> from scipy.sparse import coo_matrix, vstack 

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

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

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

564 array([[1, 2], 

565 [3, 4], 

566 [5, 6]]) 

567 

568 """ 

569 return bmat([[b] for b in blocks], format=format, dtype=dtype) 

570 

571 

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

573 """ 

574 Build a sparse matrix from sparse sub-blocks 

575 

576 Parameters 

577 ---------- 

578 blocks : array_like 

579 Grid of sparse matrices with compatible shapes. 

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

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

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

583 appropriate sparse matrix format is returned. 

584 This choice is subject to change. 

585 dtype : dtype, optional 

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

587 determined from that of `blocks`. 

588 

589 Returns 

590 ------- 

591 bmat : sparse matrix 

592 

593 See Also 

594 -------- 

595 block_diag, diags 

596 

597 Examples 

598 -------- 

599 >>> from scipy.sparse import coo_matrix, bmat 

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

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

602 >>> C = coo_matrix([[7]]) 

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

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

605 [3, 4, 6], 

606 [0, 0, 7]]) 

607 

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

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

610 [3, 4, 0], 

611 [0, 0, 7]]) 

612 

613 """ 

614 

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

616 

617 if blocks.ndim != 2: 

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

619 

620 M,N = blocks.shape 

621 

622 # check for fast path cases 

623 if (format in (None, 'csr') and all(isinstance(b, csr_matrix) 

624 for b in blocks.flat)): 

625 if N > 1: 

626 # stack along columns (axis 1): 

627 blocks = [[_stack_along_minor_axis(blocks[b, :], 1)] 

628 for b in range(M)] # must have shape: (M, 1) 

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

630 

631 # stack along rows (axis 0): 

632 A = _compressed_sparse_stack(blocks[:, 0], 0) 

633 if dtype is not None: 

634 A = A.astype(dtype) 

635 return A 

636 elif (format in (None, 'csc') and all(isinstance(b, csc_matrix) 

637 for b in blocks.flat)): 

638 if M > 1: 

639 # stack along rows (axis 0): 

640 blocks = [[_stack_along_minor_axis(blocks[:, b], 0) 

641 for b in range(N)]] # must have shape: (1, N) 

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

643 

644 # stack along columns (axis 1): 

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

646 if dtype is not None: 

647 A = A.astype(dtype) 

648 return A 

649 

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

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

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

653 

654 # convert everything to COO format 

655 for i in range(M): 

656 for j in range(N): 

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

658 A = coo_matrix(blocks[i,j]) 

659 blocks[i,j] = A 

660 block_mask[i,j] = True 

661 

662 if brow_lengths[i] == 0: 

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

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

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

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

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

668 raise ValueError(msg) 

669 

670 if bcol_lengths[j] == 0: 

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

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

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

674 f'dimensions. ' 

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

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

677 raise ValueError(msg) 

678 

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

680 if dtype is None: 

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

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

683 

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

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

686 

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

688 

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

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

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

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

693 

694 nnz = 0 

695 ii, jj = np.nonzero(block_mask) 

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

697 B = blocks[i, j] 

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

699 data[idx] = B.data 

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

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

702 nnz += B.nnz 

703 

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

705 

706 

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

708 """ 

709 Build a block diagonal sparse matrix from provided matrices. 

710 

711 Parameters 

712 ---------- 

713 mats : sequence of matrices 

714 Input matrices. 

715 format : str, optional 

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

717 is returned in "coo" format. 

718 dtype : dtype specifier, optional 

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

720 determined from that of `blocks`. 

721 

722 Returns 

723 ------- 

724 res : sparse matrix 

725 

726 Notes 

727 ----- 

728 

729 .. versionadded:: 0.11.0 

730 

731 See Also 

732 -------- 

733 bmat, diags 

734 

735 Examples 

736 -------- 

737 >>> from scipy.sparse import coo_matrix, block_diag 

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

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

740 >>> C = coo_matrix([[7]]) 

741 >>> block_diag((A, B, C)).toarray() 

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

743 [3, 4, 0, 0], 

744 [0, 0, 5, 0], 

745 [0, 0, 6, 0], 

746 [0, 0, 0, 7]]) 

747 

748 """ 

749 row = [] 

750 col = [] 

751 data = [] 

752 r_idx = 0 

753 c_idx = 0 

754 for a in mats: 

755 if isinstance(a, (list, numbers.Number)): 

756 a = coo_matrix(a) 

757 nrows, ncols = a.shape 

758 if issparse(a): 

759 a = a.tocoo() 

760 row.append(a.row + r_idx) 

761 col.append(a.col + c_idx) 

762 data.append(a.data) 

763 else: 

764 a_row, a_col = np.divmod(np.arange(nrows*ncols), ncols) 

765 row.append(a_row + r_idx) 

766 col.append(a_col + c_idx) 

767 data.append(a.ravel()) 

768 r_idx += nrows 

769 c_idx += ncols 

770 row = np.concatenate(row) 

771 col = np.concatenate(col) 

772 data = np.concatenate(data) 

773 return coo_matrix((data, (row, col)), 

774 shape=(r_idx, c_idx), 

775 dtype=dtype).asformat(format) 

776 

777 

778def random(m, n, density=0.01, format='coo', dtype=None, 

779 random_state=None, data_rvs=None): 

780 """Generate a sparse matrix of the given shape and density with randomly 

781 distributed values. 

782 

783 Parameters 

784 ---------- 

785 m, n : int 

786 shape of the matrix 

787 density : real, optional 

788 density of the generated matrix: density equal to one means a full 

789 matrix, density of 0 means a matrix with no non-zero items. 

790 format : str, optional 

791 sparse matrix format. 

792 dtype : dtype, optional 

793 type of the returned matrix values. 

794 random_state : {None, int, `numpy.random.Generator`, 

795 `numpy.random.RandomState`}, optional 

796 

797 If `seed` is None (or `np.random`), the `numpy.random.RandomState` 

798 singleton is used. 

799 If `seed` is an int, a new ``RandomState`` instance is used, 

800 seeded with `seed`. 

801 If `seed` is already a ``Generator`` or ``RandomState`` instance then 

802 that instance is used. 

803 This random state will be used 

804 for sampling the sparsity structure, but not necessarily for sampling 

805 the values of the structurally nonzero entries of the matrix. 

806 data_rvs : callable, optional 

807 Samples a requested number of random values. 

808 This function should take a single argument specifying the length 

809 of the ndarray that it will return. The structurally nonzero entries 

810 of the sparse random matrix will be taken from the array sampled 

811 by this function. By default, uniform [0, 1) random values will be 

812 sampled using the same random state as is used for sampling 

813 the sparsity structure. 

814 

815 Returns 

816 ------- 

817 res : sparse matrix 

818 

819 Notes 

820 ----- 

821 Only float types are supported for now. 

822 

823 Examples 

824 -------- 

825 >>> from scipy.sparse import random 

826 >>> from scipy import stats 

827 >>> from numpy.random import default_rng 

828 >>> rng = default_rng() 

829 >>> rvs = stats.poisson(25, loc=10).rvs 

830 >>> S = random(3, 4, density=0.25, random_state=rng, data_rvs=rvs) 

831 >>> S.A 

832 array([[ 36., 0., 33., 0.], # random 

833 [ 0., 0., 0., 0.], 

834 [ 0., 0., 36., 0.]]) 

835 

836 >>> from scipy.sparse import random 

837 >>> from scipy.stats import rv_continuous 

838 >>> class CustomDistribution(rv_continuous): 

839 ... def _rvs(self, size=None, random_state=None): 

840 ... return random_state.standard_normal(size) 

841 >>> X = CustomDistribution(seed=rng) 

842 >>> Y = X() # get a frozen version of the distribution 

843 >>> S = random(3, 4, density=0.25, random_state=rng, data_rvs=Y.rvs) 

844 >>> S.A 

845 array([[ 0. , 0. , 0. , 0. ], # random 

846 [ 0.13569738, 1.9467163 , -0.81205367, 0. ], 

847 [ 0. , 0. , 0. , 0. ]]) 

848 

849 """ 

850 if density < 0 or density > 1: 

851 raise ValueError("density expected to be 0 <= density <= 1") 

852 dtype = np.dtype(dtype) 

853 

854 mn = m * n 

855 

856 tp = np.intc 

857 if mn > np.iinfo(tp).max: 

858 tp = np.int64 

859 

860 if mn > np.iinfo(tp).max: 

861 msg = """\ 

862Trying to generate a random sparse matrix such as the product of dimensions is 

863greater than %d - this is not supported on this machine 

864""" 

865 raise ValueError(msg % np.iinfo(tp).max) 

866 

867 # Number of non zero values 

868 k = int(round(density * m * n)) 

869 

870 random_state = check_random_state(random_state) 

871 

872 if data_rvs is None: 

873 if np.issubdtype(dtype, np.integer): 

874 def data_rvs(n): 

875 return rng_integers(random_state, 

876 np.iinfo(dtype).min, 

877 np.iinfo(dtype).max, 

878 n, 

879 dtype=dtype) 

880 elif np.issubdtype(dtype, np.complexfloating): 

881 def data_rvs(n): 

882 return (random_state.uniform(size=n) +