Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/scipy/sparse/_bsr.py: 14%
310 statements
« prev ^ index » next coverage.py v7.3.2, created at 2023-12-12 06:31 +0000
1"""Compressed Block Sparse Row matrix format"""
3__docformat__ = "restructuredtext en"
5__all__ = ['bsr_matrix', 'isspmatrix_bsr']
7from warnings import warn
9import numpy as np
11from ._data import _data_matrix, _minmax_mixin
12from ._compressed import _cs_matrix
13from ._base import isspmatrix, _formats, spmatrix
14from ._sputils import (isshape, getdtype, getdata, to_native, upcast,
15 get_index_dtype, check_shape)
16from . import _sparsetools
17from ._sparsetools import (bsr_matvec, bsr_matvecs, csr_matmat_maxnnz,
18 bsr_matmat, bsr_transpose, bsr_sort_indices,
19 bsr_tocsr)
22class bsr_matrix(_cs_matrix, _minmax_mixin):
23 """Block Sparse Row matrix
25 This can be instantiated in several ways:
26 bsr_matrix(D, [blocksize=(R,C)])
27 where D is a dense matrix or 2-D ndarray.
29 bsr_matrix(S, [blocksize=(R,C)])
30 with another sparse matrix S (equivalent to S.tobsr())
32 bsr_matrix((M, N), [blocksize=(R,C), dtype])
33 to construct an empty matrix with shape (M, N)
34 dtype is optional, defaulting to dtype='d'.
36 bsr_matrix((data, ij), [blocksize=(R,C), shape=(M, N)])
37 where ``data`` and ``ij`` satisfy ``a[ij[0, k], ij[1, k]] = data[k]``
39 bsr_matrix((data, indices, indptr), [shape=(M, N)])
40 is the standard BSR representation where the block column
41 indices for row i are stored in ``indices[indptr[i]:indptr[i+1]]``
42 and their corresponding block values are stored in
43 ``data[ indptr[i]: indptr[i+1] ]``. If the shape parameter is not
44 supplied, the matrix dimensions are inferred from the index arrays.
46 Attributes
47 ----------
48 dtype : dtype
49 Data type of the matrix
50 shape : 2-tuple
51 Shape of the matrix
52 ndim : int
53 Number of dimensions (this is always 2)
54 nnz
55 Number of stored values, including explicit zeros
56 data
57 Data array of the matrix
58 indices
59 BSR format index array
60 indptr
61 BSR format index pointer array
62 blocksize
63 Block size of the matrix
64 has_sorted_indices
65 Whether indices are sorted
67 Notes
68 -----
69 Sparse matrices can be used in arithmetic operations: they support
70 addition, subtraction, multiplication, division, and matrix power.
72 **Summary of BSR format**
74 The Block Compressed Row (BSR) format is very similar to the Compressed
75 Sparse Row (CSR) format. BSR is appropriate for sparse matrices with dense
76 sub matrices like the last example below. Block matrices often arise in
77 vector-valued finite element discretizations. In such cases, BSR is
78 considerably more efficient than CSR and CSC for many sparse arithmetic
79 operations.
81 **Blocksize**
83 The blocksize (R,C) must evenly divide the shape of the matrix (M,N).
84 That is, R and C must satisfy the relationship ``M % R = 0`` and
85 ``N % C = 0``.
87 If no blocksize is specified, a simple heuristic is applied to determine
88 an appropriate blocksize.
90 Examples
91 --------
92 >>> from scipy.sparse import bsr_matrix
93 >>> import numpy as np
94 >>> bsr_matrix((3, 4), dtype=np.int8).toarray()
95 array([[0, 0, 0, 0],
96 [0, 0, 0, 0],
97 [0, 0, 0, 0]], dtype=int8)
99 >>> row = np.array([0, 0, 1, 2, 2, 2])
100 >>> col = np.array([0, 2, 2, 0, 1, 2])
101 >>> data = np.array([1, 2, 3 ,4, 5, 6])
102 >>> bsr_matrix((data, (row, col)), shape=(3, 3)).toarray()
103 array([[1, 0, 2],
104 [0, 0, 3],
105 [4, 5, 6]])
107 >>> indptr = np.array([0, 2, 3, 6])
108 >>> indices = np.array([0, 2, 2, 0, 1, 2])
109 >>> data = np.array([1, 2, 3, 4, 5, 6]).repeat(4).reshape(6, 2, 2)
110 >>> bsr_matrix((data,indices,indptr), shape=(6, 6)).toarray()
111 array([[1, 1, 0, 0, 2, 2],
112 [1, 1, 0, 0, 2, 2],
113 [0, 0, 0, 0, 3, 3],
114 [0, 0, 0, 0, 3, 3],
115 [4, 4, 5, 5, 6, 6],
116 [4, 4, 5, 5, 6, 6]])
118 """
119 format = 'bsr'
121 def __init__(self, arg1, shape=None, dtype=None, copy=False, blocksize=None):
122 _data_matrix.__init__(self)
124 if isspmatrix(arg1):
125 if isspmatrix_bsr(arg1) and copy:
126 arg1 = arg1.copy()
127 else:
128 arg1 = arg1.tobsr(blocksize=blocksize)
129 self._set_self(arg1)
131 elif isinstance(arg1,tuple):
132 if isshape(arg1):
133 # it's a tuple of matrix dimensions (M,N)
134 self._shape = check_shape(arg1)
135 M,N = self.shape
136 # process blocksize
137 if blocksize is None:
138 blocksize = (1,1)
139 else:
140 if not isshape(blocksize):
141 raise ValueError('invalid blocksize=%s' % blocksize)
142 blocksize = tuple(blocksize)
143 self.data = np.zeros((0,) + blocksize, getdtype(dtype, default=float))
145 R,C = blocksize
146 if (M % R) != 0 or (N % C) != 0:
147 raise ValueError('shape must be multiple of blocksize')
149 # Select index dtype large enough to pass array and
150 # scalar parameters to sparsetools
151 idx_dtype = get_index_dtype(maxval=max(M//R, N//C, R, C))
152 self.indices = np.zeros(0, dtype=idx_dtype)
153 self.indptr = np.zeros(M//R + 1, dtype=idx_dtype)
155 elif len(arg1) == 2:
156 # (data,(row,col)) format
157 self._set_self(
158 self._coo_container(arg1, dtype=dtype, shape=shape).tobsr(
159 blocksize=blocksize
160 )
161 )
163 elif len(arg1) == 3:
164 # (data,indices,indptr) format
165 (data, indices, indptr) = arg1
167 # Select index dtype large enough to pass array and
168 # scalar parameters to sparsetools
169 maxval = 1
170 if shape is not None:
171 maxval = max(shape)
172 if blocksize is not None:
173 maxval = max(maxval, max(blocksize))
174 idx_dtype = get_index_dtype((indices, indptr), maxval=maxval,
175 check_contents=True)
176 self.indices = np.array(indices, copy=copy, dtype=idx_dtype)
177 self.indptr = np.array(indptr, copy=copy, dtype=idx_dtype)
178 self.data = getdata(data, copy=copy, dtype=dtype)
179 if self.data.ndim != 3:
180 raise ValueError(
181 'BSR data must be 3-dimensional, got shape=%s' % (
182 self.data.shape,))
183 if blocksize is not None:
184 if not isshape(blocksize):
185 raise ValueError('invalid blocksize=%s' % (blocksize,))
186 if tuple(blocksize) != self.data.shape[1:]:
187 raise ValueError('mismatching blocksize=%s vs %s' % (
188 blocksize, self.data.shape[1:]))
189 else:
190 raise ValueError('unrecognized bsr_matrix constructor usage')
191 else:
192 # must be dense
193 try:
194 arg1 = np.asarray(arg1)
195 except Exception as e:
196 raise ValueError("unrecognized form for"
197 " %s_matrix constructor" % self.format) from e
198 arg1 = self._coo_container(
199 arg1, dtype=dtype
200 ).tobsr(blocksize=blocksize)
201 self._set_self(arg1)
203 if shape is not None:
204 self._shape = check_shape(shape)
205 else:
206 if self.shape is None:
207 # shape not already set, try to infer dimensions
208 try:
209 M = len(self.indptr) - 1
210 N = self.indices.max() + 1
211 except Exception as e:
212 raise ValueError('unable to infer matrix dimensions') from e
213 else:
214 R,C = self.blocksize
215 self._shape = check_shape((M*R,N*C))
217 if self.shape is None:
218 if shape is None:
219 # TODO infer shape here
220 raise ValueError('need to infer shape')
221 else:
222 self._shape = check_shape(shape)
224 if dtype is not None:
225 self.data = self.data.astype(dtype, copy=False)
227 self.check_format(full_check=False)
229 def check_format(self, full_check=True):
230 """check whether the matrix format is valid
232 *Parameters*:
233 full_check:
234 True - rigorous check, O(N) operations : default
235 False - basic check, O(1) operations
237 """
238 M,N = self.shape
239 R,C = self.blocksize
241 # index arrays should have integer data types
242 if self.indptr.dtype.kind != 'i':
243 warn("indptr array has non-integer dtype (%s)"
244 % self.indptr.dtype.name)
245 if self.indices.dtype.kind != 'i':
246 warn("indices array has non-integer dtype (%s)"
247 % self.indices.dtype.name)
249 idx_dtype = get_index_dtype((self.indices, self.indptr))
250 self.indptr = np.asarray(self.indptr, dtype=idx_dtype)
251 self.indices = np.asarray(self.indices, dtype=idx_dtype)
252 self.data = to_native(self.data)
254 # check array shapes
255 if self.indices.ndim != 1 or self.indptr.ndim != 1:
256 raise ValueError("indices, and indptr should be 1-D")
257 if self.data.ndim != 3:
258 raise ValueError("data should be 3-D")
260 # check index pointer
261 if (len(self.indptr) != M//R + 1):
262 raise ValueError("index pointer size (%d) should be (%d)" %
263 (len(self.indptr), M//R + 1))
264 if (self.indptr[0] != 0):
265 raise ValueError("index pointer should start with 0")
267 # check index and data arrays
268 if (len(self.indices) != len(self.data)):
269 raise ValueError("indices and data should have the same size")
270 if (self.indptr[-1] > len(self.indices)):
271 raise ValueError("Last value of index pointer should be less than "
272 "the size of index and data arrays")
274 self.prune()
276 if full_check:
277 # check format validity (more expensive)
278 if self.nnz > 0:
279 if self.indices.max() >= N//C:
280 raise ValueError("column index values must be < %d (now max %d)" % (N//C, self.indices.max()))
281 if self.indices.min() < 0:
282 raise ValueError("column index values must be >= 0")
283 if np.diff(self.indptr).min() < 0:
284 raise ValueError("index pointer values must form a "
285 "non-decreasing sequence")
287 # if not self.has_sorted_indices():
288 # warn('Indices were not in sorted order. Sorting indices.')
289 # self.sort_indices(check_first=False)
291 def _get_blocksize(self):
292 return self.data.shape[1:]
293 blocksize = property(fget=_get_blocksize)
295 def getnnz(self, axis=None):
296 if axis is not None:
297 raise NotImplementedError("getnnz over an axis is not implemented "
298 "for BSR format")
299 R,C = self.blocksize
300 return int(self.indptr[-1] * R * C)
302 getnnz.__doc__ = spmatrix.getnnz.__doc__
304 def __repr__(self):
305 format = _formats[self.getformat()][1]
306 return ("<%dx%d sparse matrix of type '%s'\n"
307 "\twith %d stored elements (blocksize = %dx%d) in %s format>" %
308 (self.shape + (self.dtype.type, self.nnz) + self.blocksize +
309 (format,)))
311 def diagonal(self, k=0):
312 rows, cols = self.shape
313 if k <= -rows or k >= cols:
314 return np.empty(0, dtype=self.data.dtype)
315 R, C = self.blocksize
316 y = np.zeros(min(rows + min(k, 0), cols - max(k, 0)),
317 dtype=upcast(self.dtype))
318 _sparsetools.bsr_diagonal(k, rows // R, cols // C, R, C,
319 self.indptr, self.indices,
320 np.ravel(self.data), y)
321 return y
323 diagonal.__doc__ = spmatrix.diagonal.__doc__
325 ##########################
326 # NotImplemented methods #
327 ##########################
329 def __getitem__(self,key):
330 raise NotImplementedError
332 def __setitem__(self,key,val):
333 raise NotImplementedError
335 ######################
336 # Arithmetic methods #
337 ######################
339 def _add_dense(self, other):
340 return self.tocoo(copy=False)._add_dense(other)
342 def _mul_vector(self, other):
343 M,N = self.shape
344 R,C = self.blocksize
346 result = np.zeros(self.shape[0], dtype=upcast(self.dtype, other.dtype))
348 bsr_matvec(M//R, N//C, R, C,
349 self.indptr, self.indices, self.data.ravel(),
350 other, result)
352 return result
354 def _mul_multivector(self,other):
355 R,C = self.blocksize
356 M,N = self.shape
357 n_vecs = other.shape[1] # number of column vectors
359 result = np.zeros((M,n_vecs), dtype=upcast(self.dtype,other.dtype))
361 bsr_matvecs(M//R, N//C, n_vecs, R, C,
362 self.indptr, self.indices, self.data.ravel(),
363 other.ravel(), result.ravel())
365 return result
367 def _mul_sparse_matrix(self, other):
368 M, K1 = self.shape
369 K2, N = other.shape
371 R,n = self.blocksize
373 # convert to this format
374 if isspmatrix_bsr(other):
375 C = other.blocksize[1]
376 else:
377 C = 1
379 from ._csr import isspmatrix_csr
381 if isspmatrix_csr(other) and n == 1:
382 other = other.tobsr(blocksize=(n,C), copy=False) # lightweight conversion
383 else:
384 other = other.tobsr(blocksize=(n,C))
386 idx_dtype = get_index_dtype((self.indptr, self.indices,
387 other.indptr, other.indices))
389 bnnz = csr_matmat_maxnnz(M//R, N//C,
390 self.indptr.astype(idx_dtype),
391 self.indices.astype(idx_dtype),
392 other.indptr.astype(idx_dtype),
393 other.indices.astype(idx_dtype))
395 idx_dtype = get_index_dtype((self.indptr, self.indices,
396 other.indptr, other.indices),
397 maxval=bnnz)
398 indptr = np.empty(self.indptr.shape, dtype=idx_dtype)
399 indices = np.empty(bnnz, dtype=idx_dtype)
400 data = np.empty(R*C*bnnz, dtype=upcast(self.dtype,other.dtype))
402 bsr_matmat(bnnz, M//R, N//C, R, C, n,
403 self.indptr.astype(idx_dtype),
404 self.indices.astype(idx_dtype),
405 np.ravel(self.data),
406 other.indptr.astype(idx_dtype),
407 other.indices.astype(idx_dtype),
408 np.ravel(other.data),
409 indptr,
410 indices,
411 data)
413 data = data.reshape(-1,R,C)
415 # TODO eliminate zeros
417 return self._bsr_container(
418 (data, indices, indptr), shape=(M, N), blocksize=(R, C)
419 )
421 ######################
422 # Conversion methods #
423 ######################
425 def tobsr(self, blocksize=None, copy=False):
426 """Convert this matrix into Block Sparse Row Format.
428 With copy=False, the data/indices may be shared between this
429 matrix and the resultant bsr_matrix.
431 If blocksize=(R, C) is provided, it will be used for determining
432 block size of the bsr_matrix.
433 """
434 if blocksize not in [None, self.blocksize]:
435 return self.tocsr().tobsr(blocksize=blocksize)
436 if copy:
437 return self.copy()
438 else:
439 return self
441 def tocsr(self, copy=False):
442 M, N = self.shape
443 R, C = self.blocksize
444 nnz = self.nnz
445 idx_dtype = get_index_dtype((self.indptr, self.indices),
446 maxval=max(nnz, N))
447 indptr = np.empty(M + 1, dtype=idx_dtype)
448 indices = np.empty(nnz, dtype=idx_dtype)
449 data = np.empty(nnz, dtype=upcast(self.dtype))
451 bsr_tocsr(M // R, # n_brow
452 N // C, # n_bcol
453 R, C,
454 self.indptr.astype(idx_dtype, copy=False),
455 self.indices.astype(idx_dtype, copy=False),
456 self.data,
457 indptr,
458 indices,
459 data)
460 return self._csr_container((data, indices, indptr), shape=self.shape)
462 tocsr.__doc__ = spmatrix.tocsr.__doc__
464 def tocsc(self, copy=False):
465 return self.tocsr(copy=False).tocsc(copy=copy)
467 tocsc.__doc__ = spmatrix.tocsc.__doc__
469 def tocoo(self, copy=True):
470 """Convert this matrix to COOrdinate format.
472 When copy=False the data array will be shared between
473 this matrix and the resultant coo_matrix.
474 """
476 M,N = self.shape
477 R,C = self.blocksize
479 indptr_diff = np.diff(self.indptr)
480 if indptr_diff.dtype.itemsize > np.dtype(np.intp).itemsize:
481 # Check for potential overflow
482 indptr_diff_limited = indptr_diff.astype(np.intp)
483 if np.any(indptr_diff_limited != indptr_diff):
484 raise ValueError("Matrix too big to convert")
485 indptr_diff = indptr_diff_limited
487 row = (R * np.arange(M//R)).repeat(indptr_diff)
488 row = row.repeat(R*C).reshape(-1,R,C)
489 row += np.tile(np.arange(R).reshape(-1,1), (1,C))
490 row = row.reshape(-1)
492 col = (C * self.indices).repeat(R*C).reshape(-1,R,C)
493 col += np.tile(np.arange(C), (R,1))
494 col = col.reshape(-1)
496 data = self.data.reshape(-1)
498 if copy:
499 data = data.copy()
501 return self._coo_container(
502 (data, (row, col)), shape=self.shape
503 )
505 def toarray(self, order=None, out=None):
506 return self.tocoo(copy=False).toarray(order=order, out=out)
508 toarray.__doc__ = spmatrix.toarray.__doc__
510 def transpose(self, axes=None, copy=False):
511 if axes is not None:
512 raise ValueError(("Sparse matrices do not support "
513 "an 'axes' parameter because swapping "
514 "dimensions is the only logical permutation."))
516 R, C = self.blocksize
517 M, N = self.shape
518 NBLK = self.nnz//(R*C)
520 if self.nnz == 0:
521 return self._bsr_container((N, M), blocksize=(C, R),
522 dtype=self.dtype, copy=copy)
524 indptr = np.empty(N//C + 1, dtype=self.indptr.dtype)
525 indices = np.empty(NBLK, dtype=self.indices.dtype)
526 data = np.empty((NBLK, C, R), dtype=self.data.dtype)
528 bsr_transpose(M//R, N//C, R, C,
529 self.indptr, self.indices, self.data.ravel(),
530 indptr, indices, data.ravel())
532 return self._bsr_container((data, indices, indptr),
533 shape=(N, M), copy=copy)
535 transpose.__doc__ = spmatrix.transpose.__doc__
537 ##############################################################
538 # methods that examine or modify the internal data structure #
539 ##############################################################
541 def eliminate_zeros(self):
542 """Remove zero elements in-place."""
544 if not self.nnz:
545 return # nothing to do
547 R,C = self.blocksize
548 M,N = self.shape
550 mask = (self.data != 0).reshape(-1,R*C).sum(axis=1) # nonzero blocks
552 nonzero_blocks = mask.nonzero()[0]
554 self.data[:len(nonzero_blocks)] = self.data[nonzero_blocks]
556 # modifies self.indptr and self.indices *in place*
557 _sparsetools.csr_eliminate_zeros(M//R, N//C, self.indptr,
558 self.indices, mask)
559 self.prune()
561 def sum_duplicates(self):
562 """Eliminate duplicate matrix entries by adding them together
564 The is an *in place* operation
565 """
566 if self.has_canonical_format:
567 return
568 self.sort_indices()
569 R, C = self.blocksize
570 M, N = self.shape
572 # port of _sparsetools.csr_sum_duplicates
573 n_row = M // R
574 nnz = 0
575 row_end = 0
576 for i in range(n_row):
577 jj = row_end
578 row_end = self.indptr[i+1]
579 while jj < row_end:
580 j = self.indices[jj]
581 x = self.data[jj]
582 jj += 1
583 while jj < row_end and self.indices[jj] == j:
584 x += self.data[jj]
585 jj += 1
586 self.indices[nnz] = j
587 self.data[nnz] = x
588 nnz += 1
589 self.indptr[i+1] = nnz
591 self.prune() # nnz may have changed
592 self.has_canonical_format = True
594 def sort_indices(self):
595 """Sort the indices of this matrix *in place*
596 """
597 if self.has_sorted_indices:
598 return
600 R,C = self.blocksize
601 M,N = self.shape
603 bsr_sort_indices(M//R, N//C, R, C, self.indptr, self.indices, self.data.ravel())
605 self.has_sorted_indices = True
607 def prune(self):
608 """ Remove empty space after all non-zero elements.
609 """
611 R,C = self.blocksize
612 M,N = self.shape
614 if len(self.indptr) != M//R + 1:
615 raise ValueError("index pointer has invalid length")
617 bnnz = self.indptr[-1]
619 if len(self.indices) < bnnz:
620 raise ValueError("indices array has too few elements")
621 if len(self.data) < bnnz:
622 raise ValueError("data array has too few elements")
624 self.data = self.data[:bnnz]
625 self.indices = self.indices[:bnnz]
627 # utility functions
628 def _binopt(self, other, op, in_shape=None, out_shape=None):
629 """Apply the binary operation fn to two sparse matrices."""
631 # Ideally we'd take the GCDs of the blocksize dimensions
632 # and explode self and other to match.
633 other = self.__class__(other, blocksize=self.blocksize)
635 # e.g. bsr_plus_bsr, etc.
636 fn = getattr(_sparsetools, self.format + op + self.format)
638 R,C = self.blocksize
640 max_bnnz = len(self.data) + len(other.data)
641 idx_dtype = get_index_dtype((self.indptr, self.indices,
642 other.indptr, other.indices),
643 maxval=max_bnnz)
644 indptr = np.empty(self.indptr.shape, dtype=idx_dtype)
645 indices = np.empty(max_bnnz, dtype=idx_dtype)
647 bool_ops = ['_ne_', '_lt_', '_gt_', '_le_', '_ge_']
648 if op in bool_ops:
649 data = np.empty(R*C*max_bnnz, dtype=np.bool_)
650 else:
651 data = np.empty(R*C*max_bnnz, dtype=upcast(self.dtype,other.dtype))
653 fn(self.shape[0]//R, self.shape[1]//C, R, C,
654 self.indptr.astype(idx_dtype),
655 self.indices.astype(idx_dtype),
656 self.data,
657 other.indptr.astype(idx_dtype),
658 other.indices.astype(idx_dtype),
659 np.ravel(other.data),
660 indptr,
661 indices,
662 data)
664 actual_bnnz = indptr[-1]
665 indices = indices[:actual_bnnz]
666 data = data[:R*C*actual_bnnz]
668 if actual_bnnz < max_bnnz/2:
669 indices = indices.copy()
670 data = data.copy()
672 data = data.reshape(-1,R,C)
674 return self.__class__((data, indices, indptr), shape=self.shape)
676 # needed by _data_matrix
677 def _with_data(self,data,copy=True):
678 """Returns a matrix with the same sparsity structure as self,
679 but with different data. By default the structure arrays
680 (i.e. .indptr and .indices) are copied.
681 """
682 if copy:
683 return self.__class__((data,self.indices.copy(),self.indptr.copy()),
684 shape=self.shape,dtype=data.dtype)
685 else:
686 return self.__class__((data,self.indices,self.indptr),
687 shape=self.shape,dtype=data.dtype)
689# # these functions are used by the parent class
690# # to remove redudancy between bsc_matrix and bsr_matrix
691# def _swap(self,x):
692# """swap the members of x if this is a column-oriented matrix
693# """
694# return (x[0],x[1])
697def isspmatrix_bsr(x):
698 """Is x of a bsr_matrix type?
700 Parameters
701 ----------
702 x
703 object to check for being a bsr matrix
705 Returns
706 -------
707 bool
708 True if x is a bsr matrix, False otherwise
710 Examples
711 --------
712 >>> from scipy.sparse import bsr_matrix, isspmatrix_bsr
713 >>> isspmatrix_bsr(bsr_matrix([[5]]))
714 True
716 >>> from scipy.sparse import bsr_matrix, csr_matrix, isspmatrix_bsr
717 >>> isspmatrix_bsr(csr_matrix([[5]]))
718 False
719 """
720 from ._arrays import bsr_array
721 return isinstance(x, bsr_matrix) or isinstance(x, bsr_array)