construct.py 25 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842
  1. """Functions to construct sparse matrices
  2. """
  3. from __future__ import division, print_function, absolute_import
  4. __docformat__ = "restructuredtext en"
  5. __all__ = ['spdiags', 'eye', 'identity', 'kron', 'kronsum',
  6. 'hstack', 'vstack', 'bmat', 'rand', 'random', 'diags', 'block_diag']
  7. import numpy as np
  8. from scipy._lib._numpy_compat import get_randint
  9. from scipy._lib.six import xrange
  10. from .sputils import upcast, get_index_dtype, isscalarlike
  11. from .csr import csr_matrix
  12. from .csc import csc_matrix
  13. from .bsr import bsr_matrix
  14. from .coo import coo_matrix
  15. from .dia import dia_matrix
  16. from .base import issparse
  17. def spdiags(data, diags, m, n, format=None):
  18. """
  19. Return a sparse matrix from diagonals.
  20. Parameters
  21. ----------
  22. data : array_like
  23. matrix diagonals stored row-wise
  24. diags : diagonals to set
  25. - k = 0 the main diagonal
  26. - k > 0 the k-th upper diagonal
  27. - k < 0 the k-th lower diagonal
  28. m, n : int
  29. shape of the result
  30. format : str, optional
  31. Format of the result. By default (format=None) an appropriate sparse
  32. matrix format is returned. This choice is subject to change.
  33. See Also
  34. --------
  35. diags : more convenient form of this function
  36. dia_matrix : the sparse DIAgonal format.
  37. Examples
  38. --------
  39. >>> from scipy.sparse import spdiags
  40. >>> data = np.array([[1, 2, 3, 4], [1, 2, 3, 4], [1, 2, 3, 4]])
  41. >>> diags = np.array([0, -1, 2])
  42. >>> spdiags(data, diags, 4, 4).toarray()
  43. array([[1, 0, 3, 0],
  44. [1, 2, 0, 4],
  45. [0, 2, 3, 0],
  46. [0, 0, 3, 4]])
  47. """
  48. return dia_matrix((data, diags), shape=(m,n)).asformat(format)
  49. def diags(diagonals, offsets=0, shape=None, format=None, dtype=None):
  50. """
  51. Construct a sparse matrix from diagonals.
  52. Parameters
  53. ----------
  54. diagonals : sequence of array_like
  55. Sequence of arrays containing the matrix diagonals,
  56. corresponding to `offsets`.
  57. offsets : sequence of int or an int, optional
  58. Diagonals to set:
  59. - k = 0 the main diagonal (default)
  60. - k > 0 the k-th upper diagonal
  61. - k < 0 the k-th lower diagonal
  62. shape : tuple of int, optional
  63. Shape of the result. If omitted, a square matrix large enough
  64. to contain the diagonals is returned.
  65. format : {"dia", "csr", "csc", "lil", ...}, optional
  66. Matrix format of the result. By default (format=None) an
  67. appropriate sparse matrix format is returned. This choice is
  68. subject to change.
  69. dtype : dtype, optional
  70. Data type of the matrix.
  71. See Also
  72. --------
  73. spdiags : construct matrix from diagonals
  74. Notes
  75. -----
  76. This function differs from `spdiags` in the way it handles
  77. off-diagonals.
  78. The result from `diags` is the sparse equivalent of::
  79. np.diag(diagonals[0], offsets[0])
  80. + ...
  81. + np.diag(diagonals[k], offsets[k])
  82. Repeated diagonal offsets are disallowed.
  83. .. versionadded:: 0.11
  84. Examples
  85. --------
  86. >>> from scipy.sparse import diags
  87. >>> diagonals = [[1, 2, 3, 4], [1, 2, 3], [1, 2]]
  88. >>> diags(diagonals, [0, -1, 2]).toarray()
  89. array([[1, 0, 1, 0],
  90. [1, 2, 0, 2],
  91. [0, 2, 3, 0],
  92. [0, 0, 3, 4]])
  93. Broadcasting of scalars is supported (but shape needs to be
  94. specified):
  95. >>> diags([1, -2, 1], [-1, 0, 1], shape=(4, 4)).toarray()
  96. array([[-2., 1., 0., 0.],
  97. [ 1., -2., 1., 0.],
  98. [ 0., 1., -2., 1.],
  99. [ 0., 0., 1., -2.]])
  100. If only one diagonal is wanted (as in `numpy.diag`), the following
  101. works as well:
  102. >>> diags([1, 2, 3], 1).toarray()
  103. array([[ 0., 1., 0., 0.],
  104. [ 0., 0., 2., 0.],
  105. [ 0., 0., 0., 3.],
  106. [ 0., 0., 0., 0.]])
  107. """
  108. # if offsets is not a sequence, assume that there's only one diagonal
  109. if isscalarlike(offsets):
  110. # now check that there's actually only one diagonal
  111. if len(diagonals) == 0 or isscalarlike(diagonals[0]):
  112. diagonals = [np.atleast_1d(diagonals)]
  113. else:
  114. raise ValueError("Different number of diagonals and offsets.")
  115. else:
  116. diagonals = list(map(np.atleast_1d, diagonals))
  117. offsets = np.atleast_1d(offsets)
  118. # Basic check
  119. if len(diagonals) != len(offsets):
  120. raise ValueError("Different number of diagonals and offsets.")
  121. # Determine shape, if omitted
  122. if shape is None:
  123. m = len(diagonals[0]) + abs(int(offsets[0]))
  124. shape = (m, m)
  125. # Determine data type, if omitted
  126. if dtype is None:
  127. dtype = np.common_type(*diagonals)
  128. # Construct data array
  129. m, n = shape
  130. M = max([min(m + offset, n - offset) + max(0, offset)
  131. for offset in offsets])
  132. M = max(0, M)
  133. data_arr = np.zeros((len(offsets), M), dtype=dtype)
  134. K = min(m, n)
  135. for j, diagonal in enumerate(diagonals):
  136. offset = offsets[j]
  137. k = max(0, offset)
  138. length = min(m + offset, n - offset, K)
  139. if length < 0:
  140. raise ValueError("Offset %d (index %d) out of bounds" % (offset, j))
  141. try:
  142. data_arr[j, k:k+length] = diagonal[...,:length]
  143. except ValueError:
  144. if len(diagonal) != length and len(diagonal) != 1:
  145. raise ValueError(
  146. "Diagonal length (index %d: %d at offset %d) does not "
  147. "agree with matrix size (%d, %d)." % (
  148. j, len(diagonal), offset, m, n))
  149. raise
  150. return dia_matrix((data_arr, offsets), shape=(m, n)).asformat(format)
  151. def identity(n, dtype='d', format=None):
  152. """Identity matrix in sparse format
  153. Returns an identity matrix with shape (n,n) using a given
  154. sparse format and dtype.
  155. Parameters
  156. ----------
  157. n : int
  158. Shape of the identity matrix.
  159. dtype : dtype, optional
  160. Data type of the matrix
  161. format : str, optional
  162. Sparse format of the result, e.g. format="csr", etc.
  163. Examples
  164. --------
  165. >>> from scipy.sparse import identity
  166. >>> identity(3).toarray()
  167. array([[ 1., 0., 0.],
  168. [ 0., 1., 0.],
  169. [ 0., 0., 1.]])
  170. >>> identity(3, dtype='int8', format='dia')
  171. <3x3 sparse matrix of type '<class 'numpy.int8'>'
  172. with 3 stored elements (1 diagonals) in DIAgonal format>
  173. """
  174. return eye(n, n, dtype=dtype, format=format)
  175. def eye(m, n=None, k=0, dtype=float, format=None):
  176. """Sparse matrix with ones on diagonal
  177. Returns a sparse (m x n) matrix where the k-th diagonal
  178. is all ones and everything else is zeros.
  179. Parameters
  180. ----------
  181. m : int
  182. Number of rows in the matrix.
  183. n : int, optional
  184. Number of columns. Default: `m`.
  185. k : int, optional
  186. Diagonal to place ones on. Default: 0 (main diagonal).
  187. dtype : dtype, optional
  188. Data type of the matrix.
  189. format : str, optional
  190. Sparse format of the result, e.g. format="csr", etc.
  191. Examples
  192. --------
  193. >>> from scipy import sparse
  194. >>> sparse.eye(3).toarray()
  195. array([[ 1., 0., 0.],
  196. [ 0., 1., 0.],
  197. [ 0., 0., 1.]])
  198. >>> sparse.eye(3, dtype=np.int8)
  199. <3x3 sparse matrix of type '<class 'numpy.int8'>'
  200. with 3 stored elements (1 diagonals) in DIAgonal format>
  201. """
  202. if n is None:
  203. n = m
  204. m,n = int(m),int(n)
  205. if m == n and k == 0:
  206. # fast branch for special formats
  207. if format in ['csr', 'csc']:
  208. idx_dtype = get_index_dtype(maxval=n)
  209. indptr = np.arange(n+1, dtype=idx_dtype)
  210. indices = np.arange(n, dtype=idx_dtype)
  211. data = np.ones(n, dtype=dtype)
  212. cls = {'csr': csr_matrix, 'csc': csc_matrix}[format]
  213. return cls((data,indices,indptr),(n,n))
  214. elif format == 'coo':
  215. idx_dtype = get_index_dtype(maxval=n)
  216. row = np.arange(n, dtype=idx_dtype)
  217. col = np.arange(n, dtype=idx_dtype)
  218. data = np.ones(n, dtype=dtype)
  219. return coo_matrix((data,(row,col)),(n,n))
  220. diags = np.ones((1, max(0, min(m + k, n))), dtype=dtype)
  221. return spdiags(diags, k, m, n).asformat(format)
  222. def kron(A, B, format=None):
  223. """kronecker product of sparse matrices A and B
  224. Parameters
  225. ----------
  226. A : sparse or dense matrix
  227. first matrix of the product
  228. B : sparse or dense matrix
  229. second matrix of the product
  230. format : str, optional
  231. format of the result (e.g. "csr")
  232. Returns
  233. -------
  234. kronecker product in a sparse matrix format
  235. Examples
  236. --------
  237. >>> from scipy import sparse
  238. >>> A = sparse.csr_matrix(np.array([[0, 2], [5, 0]]))
  239. >>> B = sparse.csr_matrix(np.array([[1, 2], [3, 4]]))
  240. >>> sparse.kron(A, B).toarray()
  241. array([[ 0, 0, 2, 4],
  242. [ 0, 0, 6, 8],
  243. [ 5, 10, 0, 0],
  244. [15, 20, 0, 0]])
  245. >>> sparse.kron(A, [[1, 2], [3, 4]]).toarray()
  246. array([[ 0, 0, 2, 4],
  247. [ 0, 0, 6, 8],
  248. [ 5, 10, 0, 0],
  249. [15, 20, 0, 0]])
  250. """
  251. B = coo_matrix(B)
  252. if (format is None or format == "bsr") and 2*B.nnz >= B.shape[0] * B.shape[1]:
  253. # B is fairly dense, use BSR
  254. A = csr_matrix(A,copy=True)
  255. output_shape = (A.shape[0]*B.shape[0], A.shape[1]*B.shape[1])
  256. if A.nnz == 0 or B.nnz == 0:
  257. # kronecker product is the zero matrix
  258. return coo_matrix(output_shape)
  259. B = B.toarray()
  260. data = A.data.repeat(B.size).reshape(-1,B.shape[0],B.shape[1])
  261. data = data * B
  262. return bsr_matrix((data,A.indices,A.indptr), shape=output_shape)
  263. else:
  264. # use COO
  265. A = coo_matrix(A)
  266. output_shape = (A.shape[0]*B.shape[0], A.shape[1]*B.shape[1])
  267. if A.nnz == 0 or B.nnz == 0:
  268. # kronecker product is the zero matrix
  269. return coo_matrix(output_shape)
  270. # expand entries of a into blocks
  271. row = A.row.repeat(B.nnz)
  272. col = A.col.repeat(B.nnz)
  273. data = A.data.repeat(B.nnz)
  274. row *= B.shape[0]
  275. col *= B.shape[1]
  276. # increment block indices
  277. row,col = row.reshape(-1,B.nnz),col.reshape(-1,B.nnz)
  278. row += B.row
  279. col += B.col
  280. row,col = row.reshape(-1),col.reshape(-1)
  281. # compute block entries
  282. data = data.reshape(-1,B.nnz) * B.data
  283. data = data.reshape(-1)
  284. return coo_matrix((data,(row,col)), shape=output_shape).asformat(format)
  285. def kronsum(A, B, format=None):
  286. """kronecker sum of sparse matrices A and B
  287. Kronecker sum of two sparse matrices is a sum of two Kronecker
  288. products kron(I_n,A) + kron(B,I_m) where A has shape (m,m)
  289. and B has shape (n,n) and I_m and I_n are identity matrices
  290. of shape (m,m) and (n,n) respectively.
  291. Parameters
  292. ----------
  293. A
  294. square matrix
  295. B
  296. square matrix
  297. format : str
  298. format of the result (e.g. "csr")
  299. Returns
  300. -------
  301. kronecker sum in a sparse matrix format
  302. Examples
  303. --------
  304. """
  305. A = coo_matrix(A)
  306. B = coo_matrix(B)
  307. if A.shape[0] != A.shape[1]:
  308. raise ValueError('A is not square')
  309. if B.shape[0] != B.shape[1]:
  310. raise ValueError('B is not square')
  311. dtype = upcast(A.dtype, B.dtype)
  312. L = kron(eye(B.shape[0],dtype=dtype), A, format=format)
  313. R = kron(B, eye(A.shape[0],dtype=dtype), format=format)
  314. return (L+R).asformat(format) # since L + R is not always same format
  315. def _compressed_sparse_stack(blocks, axis):
  316. """
  317. Stacking fast path for CSR/CSC matrices
  318. (i) vstack for CSR, (ii) hstack for CSC.
  319. """
  320. other_axis = 1 if axis == 0 else 0
  321. data = np.concatenate([b.data for b in blocks])
  322. constant_dim = blocks[0].shape[other_axis]
  323. idx_dtype = get_index_dtype(arrays=[b.indptr for b in blocks],
  324. maxval=max(data.size, constant_dim))
  325. indices = np.empty(data.size, dtype=idx_dtype)
  326. indptr = np.empty(sum(b.shape[axis] for b in blocks) + 1, dtype=idx_dtype)
  327. last_indptr = idx_dtype(0)
  328. sum_dim = 0
  329. sum_indices = 0
  330. for b in blocks:
  331. if b.shape[other_axis] != constant_dim:
  332. raise ValueError('incompatible dimensions for axis %d' % other_axis)
  333. indices[sum_indices:sum_indices+b.indices.size] = b.indices
  334. sum_indices += b.indices.size
  335. idxs = slice(sum_dim, sum_dim + b.shape[axis])
  336. indptr[idxs] = b.indptr[:-1]
  337. indptr[idxs] += last_indptr
  338. sum_dim += b.shape[axis]
  339. last_indptr += b.indptr[-1]
  340. indptr[-1] = last_indptr
  341. if axis == 0:
  342. return csr_matrix((data, indices, indptr),
  343. shape=(sum_dim, constant_dim))
  344. else:
  345. return csc_matrix((data, indices, indptr),
  346. shape=(constant_dim, sum_dim))
  347. def hstack(blocks, format=None, dtype=None):
  348. """
  349. Stack sparse matrices horizontally (column wise)
  350. Parameters
  351. ----------
  352. blocks
  353. sequence of sparse matrices with compatible shapes
  354. format : str
  355. sparse format of the result (e.g. "csr")
  356. by default an appropriate sparse matrix format is returned.
  357. This choice is subject to change.
  358. dtype : dtype, optional
  359. The data-type of the output matrix. If not given, the dtype is
  360. determined from that of `blocks`.
  361. See Also
  362. --------
  363. vstack : stack sparse matrices vertically (row wise)
  364. Examples
  365. --------
  366. >>> from scipy.sparse import coo_matrix, hstack
  367. >>> A = coo_matrix([[1, 2], [3, 4]])
  368. >>> B = coo_matrix([[5], [6]])
  369. >>> hstack([A,B]).toarray()
  370. array([[1, 2, 5],
  371. [3, 4, 6]])
  372. """
  373. return bmat([blocks], format=format, dtype=dtype)
  374. def vstack(blocks, format=None, dtype=None):
  375. """
  376. Stack sparse matrices vertically (row wise)
  377. Parameters
  378. ----------
  379. blocks
  380. sequence of sparse matrices with compatible shapes
  381. format : str, optional
  382. sparse format of the result (e.g. "csr")
  383. by default an appropriate sparse matrix format is returned.
  384. This choice is subject to change.
  385. dtype : dtype, optional
  386. The data-type of the output matrix. If not given, the dtype is
  387. determined from that of `blocks`.
  388. See Also
  389. --------
  390. hstack : stack sparse matrices horizontally (column wise)
  391. Examples
  392. --------
  393. >>> from scipy.sparse import coo_matrix, vstack
  394. >>> A = coo_matrix([[1, 2], [3, 4]])
  395. >>> B = coo_matrix([[5, 6]])
  396. >>> vstack([A, B]).toarray()
  397. array([[1, 2],
  398. [3, 4],
  399. [5, 6]])
  400. """
  401. return bmat([[b] for b in blocks], format=format, dtype=dtype)
  402. def bmat(blocks, format=None, dtype=None):
  403. """
  404. Build a sparse matrix from sparse sub-blocks
  405. Parameters
  406. ----------
  407. blocks : array_like
  408. Grid of sparse matrices with compatible shapes.
  409. An entry of None implies an all-zero matrix.
  410. format : {'bsr', 'coo', 'csc', 'csr', 'dia', 'dok', 'lil'}, optional
  411. The sparse format of the result (e.g. "csr"). By default an
  412. appropriate sparse matrix format is returned.
  413. This choice is subject to change.
  414. dtype : dtype, optional
  415. The data-type of the output matrix. If not given, the dtype is
  416. determined from that of `blocks`.
  417. Returns
  418. -------
  419. bmat : sparse matrix
  420. See Also
  421. --------
  422. block_diag, diags
  423. Examples
  424. --------
  425. >>> from scipy.sparse import coo_matrix, bmat
  426. >>> A = coo_matrix([[1, 2], [3, 4]])
  427. >>> B = coo_matrix([[5], [6]])
  428. >>> C = coo_matrix([[7]])
  429. >>> bmat([[A, B], [None, C]]).toarray()
  430. array([[1, 2, 5],
  431. [3, 4, 6],
  432. [0, 0, 7]])
  433. >>> bmat([[A, None], [None, C]]).toarray()
  434. array([[1, 2, 0],
  435. [3, 4, 0],
  436. [0, 0, 7]])
  437. """
  438. blocks = np.asarray(blocks, dtype='object')
  439. if blocks.ndim != 2:
  440. raise ValueError('blocks must be 2-D')
  441. M,N = blocks.shape
  442. # check for fast path cases
  443. if (N == 1 and format in (None, 'csr') and all(isinstance(b, csr_matrix)
  444. for b in blocks.flat)):
  445. A = _compressed_sparse_stack(blocks[:,0], 0)
  446. if dtype is not None:
  447. A = A.astype(dtype)
  448. return A
  449. elif (M == 1 and format in (None, 'csc')
  450. and all(isinstance(b, csc_matrix) for b in blocks.flat)):
  451. A = _compressed_sparse_stack(blocks[0,:], 1)
  452. if dtype is not None:
  453. A = A.astype(dtype)
  454. return A
  455. block_mask = np.zeros(blocks.shape, dtype=bool)
  456. brow_lengths = np.zeros(M, dtype=np.int64)
  457. bcol_lengths = np.zeros(N, dtype=np.int64)
  458. # convert everything to COO format
  459. for i in range(M):
  460. for j in range(N):
  461. if blocks[i,j] is not None:
  462. A = coo_matrix(blocks[i,j])
  463. blocks[i,j] = A
  464. block_mask[i,j] = True
  465. if brow_lengths[i] == 0:
  466. brow_lengths[i] = A.shape[0]
  467. elif brow_lengths[i] != A.shape[0]:
  468. msg = ('blocks[{i},:] has incompatible row dimensions. '
  469. 'Got blocks[{i},{j}].shape[0] == {got}, '
  470. 'expected {exp}.'.format(i=i, j=j,
  471. exp=brow_lengths[i],
  472. got=A.shape[0]))
  473. raise ValueError(msg)
  474. if bcol_lengths[j] == 0:
  475. bcol_lengths[j] = A.shape[1]
  476. elif bcol_lengths[j] != A.shape[1]:
  477. msg = ('blocks[:,{j}] has incompatible row dimensions. '
  478. 'Got blocks[{i},{j}].shape[1] == {got}, '
  479. 'expected {exp}.'.format(i=i, j=j,
  480. exp=bcol_lengths[j],
  481. got=A.shape[1]))
  482. raise ValueError(msg)
  483. nnz = sum(block.nnz for block in blocks[block_mask])
  484. if dtype is None:
  485. all_dtypes = [blk.dtype for blk in blocks[block_mask]]
  486. dtype = upcast(*all_dtypes) if all_dtypes else None
  487. row_offsets = np.append(0, np.cumsum(brow_lengths))
  488. col_offsets = np.append(0, np.cumsum(bcol_lengths))
  489. shape = (row_offsets[-1], col_offsets[-1])
  490. data = np.empty(nnz, dtype=dtype)
  491. idx_dtype = get_index_dtype(maxval=max(shape))
  492. row = np.empty(nnz, dtype=idx_dtype)
  493. col = np.empty(nnz, dtype=idx_dtype)
  494. nnz = 0
  495. ii, jj = np.nonzero(block_mask)
  496. for i, j in zip(ii, jj):
  497. B = blocks[i, j]
  498. idx = slice(nnz, nnz + B.nnz)
  499. data[idx] = B.data
  500. row[idx] = B.row + row_offsets[i]
  501. col[idx] = B.col + col_offsets[j]
  502. nnz += B.nnz
  503. return coo_matrix((data, (row, col)), shape=shape).asformat(format)
  504. def block_diag(mats, format=None, dtype=None):
  505. """
  506. Build a block diagonal sparse matrix from provided matrices.
  507. Parameters
  508. ----------
  509. mats : sequence of matrices
  510. Input matrices.
  511. format : str, optional
  512. The sparse format of the result (e.g. "csr"). If not given, the matrix
  513. is returned in "coo" format.
  514. dtype : dtype specifier, optional
  515. The data-type of the output matrix. If not given, the dtype is
  516. determined from that of `blocks`.
  517. Returns
  518. -------
  519. res : sparse matrix
  520. Notes
  521. -----
  522. .. versionadded:: 0.11.0
  523. See Also
  524. --------
  525. bmat, diags
  526. Examples
  527. --------
  528. >>> from scipy.sparse import coo_matrix, block_diag
  529. >>> A = coo_matrix([[1, 2], [3, 4]])
  530. >>> B = coo_matrix([[5], [6]])
  531. >>> C = coo_matrix([[7]])
  532. >>> block_diag((A, B, C)).toarray()
  533. array([[1, 2, 0, 0],
  534. [3, 4, 0, 0],
  535. [0, 0, 5, 0],
  536. [0, 0, 6, 0],
  537. [0, 0, 0, 7]])
  538. """
  539. nmat = len(mats)
  540. rows = []
  541. for ia, a in enumerate(mats):
  542. row = [None]*nmat
  543. if issparse(a):
  544. row[ia] = a
  545. else:
  546. row[ia] = coo_matrix(a)
  547. rows.append(row)
  548. return bmat(rows, format=format, dtype=dtype)
  549. def random(m, n, density=0.01, format='coo', dtype=None,
  550. random_state=None, data_rvs=None):
  551. """Generate a sparse matrix of the given shape and density with randomly
  552. distributed values.
  553. Parameters
  554. ----------
  555. m, n : int
  556. shape of the matrix
  557. density : real, optional
  558. density of the generated matrix: density equal to one means a full
  559. matrix, density of 0 means a matrix with no non-zero items.
  560. format : str, optional
  561. sparse matrix format.
  562. dtype : dtype, optional
  563. type of the returned matrix values.
  564. random_state : {numpy.random.RandomState, int}, optional
  565. Random number generator or random seed. If not given, the singleton
  566. numpy.random will be used. This random state will be used
  567. for sampling the sparsity structure, but not necessarily for sampling
  568. the values of the structurally nonzero entries of the matrix.
  569. data_rvs : callable, optional
  570. Samples a requested number of random values.
  571. This function should take a single argument specifying the length
  572. of the ndarray that it will return. The structurally nonzero entries
  573. of the sparse random matrix will be taken from the array sampled
  574. by this function. By default, uniform [0, 1) random values will be
  575. sampled using the same random state as is used for sampling
  576. the sparsity structure.
  577. Returns
  578. -------
  579. res : sparse matrix
  580. Notes
  581. -----
  582. Only float types are supported for now.
  583. Examples
  584. --------
  585. >>> from scipy.sparse import random
  586. >>> from scipy import stats
  587. >>> class CustomRandomState(np.random.RandomState):
  588. ... def randint(self, k):
  589. ... i = np.random.randint(k)
  590. ... return i - i % 2
  591. >>> np.random.seed(12345)
  592. >>> rs = CustomRandomState()
  593. >>> rvs = stats.poisson(25, loc=10).rvs
  594. >>> S = random(3, 4, density=0.25, random_state=rs, data_rvs=rvs)
  595. >>> S.A
  596. array([[ 36., 0., 33., 0.], # random
  597. [ 0., 0., 0., 0.],
  598. [ 0., 0., 36., 0.]])
  599. >>> from scipy.sparse import random
  600. >>> from scipy.stats import rv_continuous
  601. >>> class CustomDistribution(rv_continuous):
  602. ... def _rvs(self, *args, **kwargs):
  603. ... return self._random_state.randn(*self._size)
  604. >>> X = CustomDistribution(seed=2906)
  605. >>> Y = X() # get a frozen version of the distribution
  606. >>> S = random(3, 4, density=0.25, random_state=2906, data_rvs=Y.rvs)
  607. >>> S.A
  608. array([[ 0. , 0. , 0. , 0. ],
  609. [ 0.13569738, 1.9467163 , -0.81205367, 0. ],
  610. [ 0. , 0. , 0. , 0. ]])
  611. """
  612. if density < 0 or density > 1:
  613. raise ValueError("density expected to be 0 <= density <= 1")
  614. dtype = np.dtype(dtype)
  615. mn = m * n
  616. tp = np.intc
  617. if mn > np.iinfo(tp).max:
  618. tp = np.int64
  619. if mn > np.iinfo(tp).max:
  620. msg = """\
  621. Trying to generate a random sparse matrix such as the product of dimensions is
  622. greater than %d - this is not supported on this machine
  623. """
  624. raise ValueError(msg % np.iinfo(tp).max)
  625. # Number of non zero values
  626. k = int(density * m * n)
  627. if random_state is None:
  628. random_state = np.random
  629. elif isinstance(random_state, (int, np.integer)):
  630. random_state = np.random.RandomState(random_state)
  631. if data_rvs is None:
  632. if np.issubdtype(dtype, np.integer):
  633. randint = get_randint(random_state)
  634. def data_rvs(n):
  635. return randint(np.iinfo(dtype).min, np.iinfo(dtype).max,
  636. n, dtype=dtype)
  637. elif np.issubdtype(dtype, np.complexfloating):
  638. def data_rvs(n):
  639. return random_state.rand(n) + random_state.rand(n) * 1j
  640. else:
  641. data_rvs = random_state.rand
  642. ind = random_state.choice(mn, size=k, replace=False)
  643. j = np.floor(ind * 1. / m).astype(tp, copy=False)
  644. i = (ind - j * m).astype(tp, copy=False)
  645. vals = data_rvs(k).astype(dtype, copy=False)
  646. return coo_matrix((vals, (i, j)), shape=(m, n)).asformat(format,
  647. copy=False)
  648. def rand(m, n, density=0.01, format="coo", dtype=None, random_state=None):
  649. """Generate a sparse matrix of the given shape and density with uniformly
  650. distributed values.
  651. Parameters
  652. ----------
  653. m, n : int
  654. shape of the matrix
  655. density : real, optional
  656. density of the generated matrix: density equal to one means a full
  657. matrix, density of 0 means a matrix with no non-zero items.
  658. format : str, optional
  659. sparse matrix format.
  660. dtype : dtype, optional
  661. type of the returned matrix values.
  662. random_state : {numpy.random.RandomState, int}, optional
  663. Random number generator or random seed. If not given, the singleton
  664. numpy.random will be used.
  665. Returns
  666. -------
  667. res : sparse matrix
  668. Notes
  669. -----
  670. Only float types are supported for now.
  671. See Also
  672. --------
  673. scipy.sparse.random : Similar function that allows a user-specified random
  674. data source.
  675. Examples
  676. --------
  677. >>> from scipy.sparse import rand
  678. >>> matrix = rand(3, 4, density=0.25, format="csr", random_state=42)
  679. >>> matrix
  680. <3x4 sparse matrix of type '<class 'numpy.float64'>'
  681. with 3 stored elements in Compressed Sparse Row format>
  682. >>> matrix.todense()
  683. matrix([[0.05641158, 0. , 0. , 0.65088847],
  684. [0. , 0. , 0. , 0.14286682],
  685. [0. , 0. , 0. , 0. ]])
  686. """
  687. return random(m, n, density, format, dtype, random_state)