Coverage for /home/martinb/.local/share/virtualenvs/camcops/lib/python3.6/site-packages/scipy/sparse/bsr.py : 15%

Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
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, to_native, upcast, get_index_dtype,
15 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 >>> bsr_matrix((3, 4), dtype=np.int8).toarray()
94 array([[0, 0, 0, 0],
95 [0, 0, 0, 0],
96 [0, 0, 0, 0]], dtype=int8)
98 >>> row = np.array([0, 0, 1, 2, 2, 2])
99 >>> col = np.array([0, 2, 2, 0, 1, 2])
100 >>> data = np.array([1, 2, 3 ,4, 5, 6])
101 >>> bsr_matrix((data, (row, col)), shape=(3, 3)).toarray()
102 array([[1, 0, 2],
103 [0, 0, 3],
104 [4, 5, 6]])
106 >>> indptr = np.array([0, 2, 3, 6])
107 >>> indices = np.array([0, 2, 2, 0, 1, 2])
108 >>> data = np.array([1, 2, 3, 4, 5, 6]).repeat(4).reshape(6, 2, 2)
109 >>> bsr_matrix((data,indices,indptr), shape=(6, 6)).toarray()
110 array([[1, 1, 0, 0, 2, 2],
111 [1, 1, 0, 0, 2, 2],
112 [0, 0, 0, 0, 3, 3],
113 [0, 0, 0, 0, 3, 3],
114 [4, 4, 5, 5, 6, 6],
115 [4, 4, 5, 5, 6, 6]])
117 """
118 format = 'bsr'
120 def __init__(self, arg1, shape=None, dtype=None, copy=False, blocksize=None):
121 _data_matrix.__init__(self)
123 if isspmatrix(arg1):
124 if isspmatrix_bsr(arg1) and copy:
125 arg1 = arg1.copy()
126 else:
127 arg1 = arg1.tobsr(blocksize=blocksize)
128 self._set_self(arg1)
130 elif isinstance(arg1,tuple):
131 if isshape(arg1):
132 # it's a tuple of matrix dimensions (M,N)
133 self._shape = check_shape(arg1)
134 M,N = self.shape
135 # process blocksize
136 if blocksize is None:
137 blocksize = (1,1)
138 else:
139 if not isshape(blocksize):
140 raise ValueError('invalid blocksize=%s' % blocksize)
141 blocksize = tuple(blocksize)
142 self.data = np.zeros((0,) + blocksize, getdtype(dtype, default=float))
144 R,C = blocksize
145 if (M % R) != 0 or (N % C) != 0:
146 raise ValueError('shape must be multiple of blocksize')
148 # Select index dtype large enough to pass array and
149 # scalar parameters to sparsetools
150 idx_dtype = get_index_dtype(maxval=max(M//R, N//C, R, C))
151 self.indices = np.zeros(0, dtype=idx_dtype)
152 self.indptr = np.zeros(M//R + 1, dtype=idx_dtype)
154 elif len(arg1) == 2:
155 # (data,(row,col)) format
156 from .coo import coo_matrix
157 self._set_self(coo_matrix(arg1, dtype=dtype).tobsr(blocksize=blocksize))
159 elif len(arg1) == 3:
160 # (data,indices,indptr) format
161 (data, indices, indptr) = arg1
163 # Select index dtype large enough to pass array and
164 # scalar parameters to sparsetools
165 maxval = 1
166 if shape is not None:
167 maxval = max(shape)
168 if blocksize is not None:
169 maxval = max(maxval, max(blocksize))
170 idx_dtype = get_index_dtype((indices, indptr), maxval=maxval, check_contents=True)
172 self.indices = np.array(indices, copy=copy, dtype=idx_dtype)
173 self.indptr = np.array(indptr, copy=copy, dtype=idx_dtype)
174 self.data = np.array(data, copy=copy, dtype=getdtype(dtype, data))
175 else:
176 raise ValueError('unrecognized bsr_matrix constructor usage')
177 else:
178 # must be dense
179 try:
180 arg1 = np.asarray(arg1)
181 except Exception:
182 raise ValueError("unrecognized form for"
183 " %s_matrix constructor" % self.format)
184 from .coo import coo_matrix
185 arg1 = coo_matrix(arg1, dtype=dtype).tobsr(blocksize=blocksize)
186 self._set_self(arg1)
188 if shape is not None:
189 self._shape = check_shape(shape)
190 else:
191 if self.shape is None:
192 # shape not already set, try to infer dimensions
193 try:
194 M = len(self.indptr) - 1
195 N = self.indices.max() + 1
196 except Exception:
197 raise ValueError('unable to infer matrix dimensions')
198 else:
199 R,C = self.blocksize
200 self._shape = check_shape((M*R,N*C))
202 if self.shape is None:
203 if shape is None:
204 # TODO infer shape here
205 raise ValueError('need to infer shape')
206 else:
207 self._shape = check_shape(shape)
209 if dtype is not None:
210 self.data = self.data.astype(dtype, copy=False)
212 self.check_format(full_check=False)
214 def check_format(self, full_check=True):
215 """check whether the matrix format is valid
217 *Parameters*:
218 full_check:
219 True - rigorous check, O(N) operations : default
220 False - basic check, O(1) operations
222 """
223 M,N = self.shape
224 R,C = self.blocksize
226 # index arrays should have integer data types
227 if self.indptr.dtype.kind != 'i':
228 warn("indptr array has non-integer dtype (%s)"
229 % self.indptr.dtype.name)
230 if self.indices.dtype.kind != 'i':
231 warn("indices array has non-integer dtype (%s)"
232 % self.indices.dtype.name)
234 idx_dtype = get_index_dtype((self.indices, self.indptr))
235 self.indptr = np.asarray(self.indptr, dtype=idx_dtype)
236 self.indices = np.asarray(self.indices, dtype=idx_dtype)
237 self.data = to_native(self.data)
239 # check array shapes
240 if self.indices.ndim != 1 or self.indptr.ndim != 1:
241 raise ValueError("indices, and indptr should be 1-D")
242 if self.data.ndim != 3:
243 raise ValueError("data should be 3-D")
245 # check index pointer
246 if (len(self.indptr) != M//R + 1):
247 raise ValueError("index pointer size (%d) should be (%d)" %
248 (len(self.indptr), M//R + 1))
249 if (self.indptr[0] != 0):
250 raise ValueError("index pointer should start with 0")
252 # check index and data arrays
253 if (len(self.indices) != len(self.data)):
254 raise ValueError("indices and data should have the same size")
255 if (self.indptr[-1] > len(self.indices)):
256 raise ValueError("Last value of index pointer should be less than "
257 "the size of index and data arrays")
259 self.prune()
261 if full_check:
262 # check format validity (more expensive)
263 if self.nnz > 0:
264 if self.indices.max() >= N//C:
265 raise ValueError("column index values must be < %d (now max %d)" % (N//C, self.indices.max()))
266 if self.indices.min() < 0:
267 raise ValueError("column index values must be >= 0")
268 if np.diff(self.indptr).min() < 0:
269 raise ValueError("index pointer values must form a "
270 "non-decreasing sequence")
272 # if not self.has_sorted_indices():
273 # warn('Indices were not in sorted order. Sorting indices.')
274 # self.sort_indices(check_first=False)
276 def _get_blocksize(self):
277 return self.data.shape[1:]
278 blocksize = property(fget=_get_blocksize)
280 def getnnz(self, axis=None):
281 if axis is not None:
282 raise NotImplementedError("getnnz over an axis is not implemented "
283 "for BSR format")
284 R,C = self.blocksize
285 return int(self.indptr[-1] * R * C)
287 getnnz.__doc__ = spmatrix.getnnz.__doc__
289 def __repr__(self):
290 format = _formats[self.getformat()][1]
291 return ("<%dx%d sparse matrix of type '%s'\n"
292 "\twith %d stored elements (blocksize = %dx%d) in %s format>" %
293 (self.shape + (self.dtype.type, self.nnz) + self.blocksize +
294 (format,)))
296 def diagonal(self, k=0):
297 rows, cols = self.shape
298 if k <= -rows or k >= cols:
299 return np.empty(0, dtype=self.data.dtype)
300 R, C = self.blocksize
301 y = np.zeros(min(rows + min(k, 0), cols - max(k, 0)),
302 dtype=upcast(self.dtype))
303 _sparsetools.bsr_diagonal(k, rows // R, cols // C, R, C,
304 self.indptr, self.indices,
305 np.ravel(self.data), y)
306 return y
308 diagonal.__doc__ = spmatrix.diagonal.__doc__
310 ##########################
311 # NotImplemented methods #
312 ##########################
314 def __getitem__(self,key):
315 raise NotImplementedError
317 def __setitem__(self,key,val):
318 raise NotImplementedError
320 ######################
321 # Arithmetic methods #
322 ######################
324 @np.deprecate(message="BSR matvec is deprecated in SciPy 0.19.0. "
325 "Use * operator instead.")
326 def matvec(self, other):
327 """Multiply matrix by vector."""
328 return self * other
330 @np.deprecate(message="BSR matmat is deprecated in SciPy 0.19.0. "
331 "Use * operator instead.")
332 def matmat(self, other):
333 """Multiply this sparse matrix by other matrix."""
334 return self * other
336 def _add_dense(self, other):
337 return self.tocoo(copy=False)._add_dense(other)
339 def _mul_vector(self, other):
340 M,N = self.shape
341 R,C = self.blocksize
343 result = np.zeros(self.shape[0], dtype=upcast(self.dtype, other.dtype))
345 bsr_matvec(M//R, N//C, R, C,
346 self.indptr, self.indices, self.data.ravel(),
347 other, result)
349 return result
351 def _mul_multivector(self,other):
352 R,C = self.blocksize
353 M,N = self.shape
354 n_vecs = other.shape[1] # number of column vectors
356 result = np.zeros((M,n_vecs), dtype=upcast(self.dtype,other.dtype))
358 bsr_matvecs(M//R, N//C, n_vecs, R, C,
359 self.indptr, self.indices, self.data.ravel(),
360 other.ravel(), result.ravel())
362 return result
364 def _mul_sparse_matrix(self, other):
365 M, K1 = self.shape
366 K2, N = other.shape
368 R,n = self.blocksize
370 # convert to this format
371 if isspmatrix_bsr(other):
372 C = other.blocksize[1]
373 else:
374 C = 1
376 from .csr import isspmatrix_csr
378 if isspmatrix_csr(other) and n == 1:
379 other = other.tobsr(blocksize=(n,C), copy=False) # lightweight conversion
380 else:
381 other = other.tobsr(blocksize=(n,C))
383 idx_dtype = get_index_dtype((self.indptr, self.indices,
384 other.indptr, other.indices))
386 bnnz = csr_matmat_maxnnz(M//R, N//C,
387 self.indptr.astype(idx_dtype),
388 self.indices.astype(idx_dtype),
389 other.indptr.astype(idx_dtype),
390 other.indices.astype(idx_dtype))
392 idx_dtype = get_index_dtype((self.indptr, self.indices,
393 other.indptr, other.indices),
394 maxval=bnnz)
395 indptr = np.empty(self.indptr.shape, dtype=idx_dtype)
396 indices = np.empty(bnnz, dtype=idx_dtype)
397 data = np.empty(R*C*bnnz, dtype=upcast(self.dtype,other.dtype))
399 bsr_matmat(bnnz, M//R, N//C, R, C, n,
400 self.indptr.astype(idx_dtype),
401 self.indices.astype(idx_dtype),
402 np.ravel(self.data),
403 other.indptr.astype(idx_dtype),
404 other.indices.astype(idx_dtype),
405 np.ravel(other.data),
406 indptr,
407 indices,
408 data)
410 data = data.reshape(-1,R,C)
412 # TODO eliminate zeros
414 return bsr_matrix((data,indices,indptr),shape=(M,N),blocksize=(R,C))
416 ######################
417 # Conversion methods #
418 ######################
420 def tobsr(self, blocksize=None, copy=False):
421 """Convert this matrix into Block Sparse Row Format.
423 With copy=False, the data/indices may be shared between this
424 matrix and the resultant bsr_matrix.
426 If blocksize=(R, C) is provided, it will be used for determining
427 block size of the bsr_matrix.
428 """
429 if blocksize not in [None, self.blocksize]:
430 return self.tocsr().tobsr(blocksize=blocksize)
431 if copy:
432 return self.copy()
433 else:
434 return self
436 def tocsr(self, copy=False):
437 M, N = self.shape
438 R, C = self.blocksize
439 nnz = self.nnz
440 idx_dtype = get_index_dtype((self.indptr, self.indices),
441 maxval=max(nnz, N))
442 indptr = np.empty(M + 1, dtype=idx_dtype)
443 indices = np.empty(nnz, dtype=idx_dtype)
444 data = np.empty(nnz, dtype=upcast(self.dtype))
446 bsr_tocsr(M // R, # n_brow
447 N // C, # n_bcol
448 R, C,
449 self.indptr.astype(idx_dtype, copy=False),
450 self.indices.astype(idx_dtype, copy=False),
451 self.data,
452 indptr,
453 indices,
454 data)
455 from .csr import csr_matrix
456 return csr_matrix((data, indices, indptr), shape=self.shape)
458 tocsr.__doc__ = spmatrix.tocsr.__doc__
460 def tocsc(self, copy=False):
461 return self.tocsr(copy=False).tocsc(copy=copy)
463 tocsc.__doc__ = spmatrix.tocsc.__doc__
465 def tocoo(self, copy=True):
466 """Convert this matrix to COOrdinate format.
468 When copy=False the data array will be shared between
469 this matrix and the resultant coo_matrix.
470 """
472 M,N = self.shape
473 R,C = self.blocksize
475 indptr_diff = np.diff(self.indptr)
476 if indptr_diff.dtype.itemsize > np.dtype(np.intp).itemsize:
477 # Check for potential overflow
478 indptr_diff_limited = indptr_diff.astype(np.intp)
479 if np.any(indptr_diff_limited != indptr_diff):
480 raise ValueError("Matrix too big to convert")
481 indptr_diff = indptr_diff_limited
483 row = (R * np.arange(M//R)).repeat(indptr_diff)
484 row = row.repeat(R*C).reshape(-1,R,C)
485 row += np.tile(np.arange(R).reshape(-1,1), (1,C))
486 row = row.reshape(-1)
488 col = (C * self.indices).repeat(R*C).reshape(-1,R,C)
489 col += np.tile(np.arange(C), (R,1))
490 col = col.reshape(-1)
492 data = self.data.reshape(-1)
494 if copy:
495 data = data.copy()
497 from .coo import coo_matrix
498 return coo_matrix((data,(row,col)), shape=self.shape)
500 def toarray(self, order=None, out=None):
501 return self.tocoo(copy=False).toarray(order=order, out=out)
503 toarray.__doc__ = spmatrix.toarray.__doc__
505 def transpose(self, axes=None, copy=False):
506 if axes is not None:
507 raise ValueError(("Sparse matrices do not support "
508 "an 'axes' parameter because swapping "
509 "dimensions is the only logical permutation."))
511 R, C = self.blocksize
512 M, N = self.shape
513 NBLK = self.nnz//(R*C)
515 if self.nnz == 0:
516 return bsr_matrix((N, M), blocksize=(C, R),
517 dtype=self.dtype, copy=copy)
519 indptr = np.empty(N//C + 1, dtype=self.indptr.dtype)
520 indices = np.empty(NBLK, dtype=self.indices.dtype)
521 data = np.empty((NBLK, C, R), dtype=self.data.dtype)
523 bsr_transpose(M//R, N//C, R, C,
524 self.indptr, self.indices, self.data.ravel(),
525 indptr, indices, data.ravel())
527 return bsr_matrix((data, indices, indptr),
528 shape=(N, M), copy=copy)
530 transpose.__doc__ = spmatrix.transpose.__doc__
532 ##############################################################
533 # methods that examine or modify the internal data structure #
534 ##############################################################
536 def eliminate_zeros(self):
537 """Remove zero elements in-place."""
539 if not self.nnz:
540 return # nothing to do
542 R,C = self.blocksize
543 M,N = self.shape
545 mask = (self.data != 0).reshape(-1,R*C).sum(axis=1) # nonzero blocks
547 nonzero_blocks = mask.nonzero()[0]
549 self.data[:len(nonzero_blocks)] = self.data[nonzero_blocks]
551 # modifies self.indptr and self.indices *in place*
552 _sparsetools.csr_eliminate_zeros(M//R, N//C, self.indptr,
553 self.indices, mask)
554 self.prune()
556 def sum_duplicates(self):
557 """Eliminate duplicate matrix entries by adding them together
559 The is an *in place* operation
560 """
561 if self.has_canonical_format:
562 return
563 self.sort_indices()
564 R, C = self.blocksize
565 M, N = self.shape
567 # port of _sparsetools.csr_sum_duplicates
568 n_row = M // R
569 nnz = 0
570 row_end = 0
571 for i in range(n_row):
572 jj = row_end
573 row_end = self.indptr[i+1]
574 while jj < row_end:
575 j = self.indices[jj]
576 x = self.data[jj]
577 jj += 1
578 while jj < row_end and self.indices[jj] == j:
579 x += self.data[jj]
580 jj += 1
581 self.indices[nnz] = j
582 self.data[nnz] = x
583 nnz += 1
584 self.indptr[i+1] = nnz
586 self.prune() # nnz may have changed
587 self.has_canonical_format = True
589 def sort_indices(self):
590 """Sort the indices of this matrix *in place*
591 """
592 if self.has_sorted_indices:
593 return
595 R,C = self.blocksize
596 M,N = self.shape
598 bsr_sort_indices(M//R, N//C, R, C, self.indptr, self.indices, self.data.ravel())
600 self.has_sorted_indices = True
602 def prune(self):
603 """ Remove empty space after all non-zero elements.
604 """
606 R,C = self.blocksize
607 M,N = self.shape
609 if len(self.indptr) != M//R + 1:
610 raise ValueError("index pointer has invalid length")
612 bnnz = self.indptr[-1]
614 if len(self.indices) < bnnz:
615 raise ValueError("indices array has too few elements")
616 if len(self.data) < bnnz:
617 raise ValueError("data array has too few elements")
619 self.data = self.data[:bnnz]
620 self.indices = self.indices[:bnnz]
622 # utility functions
623 def _binopt(self, other, op, in_shape=None, out_shape=None):
624 """Apply the binary operation fn to two sparse matrices."""
626 # Ideally we'd take the GCDs of the blocksize dimensions
627 # and explode self and other to match.
628 other = self.__class__(other, blocksize=self.blocksize)
630 # e.g. bsr_plus_bsr, etc.
631 fn = getattr(_sparsetools, self.format + op + self.format)
633 R,C = self.blocksize
635 max_bnnz = len(self.data) + len(other.data)
636 idx_dtype = get_index_dtype((self.indptr, self.indices,
637 other.indptr, other.indices),
638 maxval=max_bnnz)
639 indptr = np.empty(self.indptr.shape, dtype=idx_dtype)
640 indices = np.empty(max_bnnz, dtype=idx_dtype)
642 bool_ops = ['_ne_', '_lt_', '_gt_', '_le_', '_ge_']
643 if op in bool_ops:
644 data = np.empty(R*C*max_bnnz, dtype=np.bool_)
645 else:
646 data = np.empty(R*C*max_bnnz, dtype=upcast(self.dtype,other.dtype))
648 fn(self.shape[0]//R, self.shape[1]//C, R, C,
649 self.indptr.astype(idx_dtype),
650 self.indices.astype(idx_dtype),
651 self.data,
652 other.indptr.astype(idx_dtype),
653 other.indices.astype(idx_dtype),
654 np.ravel(other.data),
655 indptr,
656 indices,
657 data)
659 actual_bnnz = indptr[-1]
660 indices = indices[:actual_bnnz]
661 data = data[:R*C*actual_bnnz]
663 if actual_bnnz < max_bnnz/2:
664 indices = indices.copy()
665 data = data.copy()
667 data = data.reshape(-1,R,C)
669 return self.__class__((data, indices, indptr), shape=self.shape)
671 # needed by _data_matrix
672 def _with_data(self,data,copy=True):
673 """Returns a matrix with the same sparsity structure as self,
674 but with different data. By default the structure arrays
675 (i.e. .indptr and .indices) are copied.
676 """
677 if copy:
678 return self.__class__((data,self.indices.copy(),self.indptr.copy()),
679 shape=self.shape,dtype=data.dtype)
680 else:
681 return self.__class__((data,self.indices,self.indptr),
682 shape=self.shape,dtype=data.dtype)
684# # these functions are used by the parent class
685# # to remove redudancy between bsc_matrix and bsr_matrix
686# def _swap(self,x):
687# """swap the members of x if this is a column-oriented matrix
688# """
689# return (x[0],x[1])
692def isspmatrix_bsr(x):
693 """Is x of a bsr_matrix type?
695 Parameters
696 ----------
697 x
698 object to check for being a bsr matrix
700 Returns
701 -------
702 bool
703 True if x is a bsr matrix, False otherwise
705 Examples
706 --------
707 >>> from scipy.sparse import bsr_matrix, isspmatrix_bsr
708 >>> isspmatrix_bsr(bsr_matrix([[5]]))
709 True
711 >>> from scipy.sparse import bsr_matrix, csr_matrix, isspmatrix_bsr
712 >>> isspmatrix_bsr(csr_matrix([[5]]))
713 False
714 """
715 return isinstance(x, bsr_matrix)