Coverage for /usr/lib/python3/dist-packages/scipy/sparse/_compressed.py: 11%
738 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1"""Base class for sparse matrix formats using compressed storage."""
2__all__ = []
4from warnings import warn
5import operator
7import numpy as np
8from scipy._lib._util import _prune_array
10from ._base import _spbase, issparse, SparseEfficiencyWarning
11from ._data import _data_matrix, _minmax_mixin
12from . import _sparsetools
13from ._sparsetools import (get_csr_submatrix, csr_sample_offsets, csr_todense,
14 csr_sample_values, csr_row_index, csr_row_slice,
15 csr_column_index1, csr_column_index2)
16from ._index import IndexMixin
17from ._sputils import (upcast, upcast_char, to_native, isdense, isshape,
18 getdtype, isscalarlike, isintlike, downcast_intp_index, get_sum_dtype, check_shape,
19 is_pydata_spmatrix)
22class _cs_matrix(_data_matrix, _minmax_mixin, IndexMixin):
23 """base matrix class for compressed row- and column-oriented matrices"""
25 def __init__(self, arg1, shape=None, dtype=None, copy=False):
26 _data_matrix.__init__(self)
28 if issparse(arg1):
29 if arg1.format == self.format and copy:
30 arg1 = arg1.copy()
31 else:
32 arg1 = arg1.asformat(self.format)
33 self._set_self(arg1)
35 elif isinstance(arg1, tuple):
36 if isshape(arg1):
37 # It's a tuple of matrix dimensions (M, N)
38 # create empty matrix
39 self._shape = check_shape(arg1)
40 M, N = self.shape
41 # Select index dtype large enough to pass array and
42 # scalar parameters to sparsetools
43 idx_dtype = self._get_index_dtype(maxval=max(M, N))
44 self.data = np.zeros(0, getdtype(dtype, default=float))
45 self.indices = np.zeros(0, idx_dtype)
46 self.indptr = np.zeros(self._swap((M, N))[0] + 1,
47 dtype=idx_dtype)
48 else:
49 if len(arg1) == 2:
50 # (data, ij) format
51 other = self.__class__(
52 self._coo_container(arg1, shape=shape, dtype=dtype)
53 )
54 self._set_self(other)
55 elif len(arg1) == 3:
56 # (data, indices, indptr) format
57 (data, indices, indptr) = arg1
59 # Select index dtype large enough to pass array and
60 # scalar parameters to sparsetools
61 maxval = None
62 if shape is not None:
63 maxval = max(shape)
64 idx_dtype = self._get_index_dtype((indices, indptr),
65 maxval=maxval,
66 check_contents=True)
68 self.indices = np.array(indices, copy=copy,
69 dtype=idx_dtype)
70 self.indptr = np.array(indptr, copy=copy, dtype=idx_dtype)
71 self.data = np.array(data, copy=copy, dtype=dtype)
72 else:
73 raise ValueError("unrecognized {}_matrix "
74 "constructor usage".format(self.format))
76 else:
77 # must be dense
78 try:
79 arg1 = np.asarray(arg1)
80 except Exception as e:
81 raise ValueError("unrecognized {}_matrix constructor usage"
82 "".format(self.format)) from e
83 self._set_self(self.__class__(
84 self._coo_container(arg1, dtype=dtype)
85 ))
87 # Read matrix dimensions given, if any
88 if shape is not None:
89 self._shape = check_shape(shape)
90 else:
91 if self.shape is None:
92 # shape not already set, try to infer dimensions
93 try:
94 major_dim = len(self.indptr) - 1
95 minor_dim = self.indices.max() + 1
96 except Exception as e:
97 raise ValueError('unable to infer matrix dimensions') from e
98 else:
99 self._shape = check_shape(self._swap((major_dim,
100 minor_dim)))
102 if dtype is not None:
103 self.data = self.data.astype(dtype, copy=False)
105 self.check_format(full_check=False)
107 def _getnnz(self, axis=None):
108 if axis is None:
109 return int(self.indptr[-1])
110 else:
111 if axis < 0:
112 axis += 2
113 axis, _ = self._swap((axis, 1 - axis))
114 _, N = self._swap(self.shape)
115 if axis == 0:
116 return np.bincount(downcast_intp_index(self.indices),
117 minlength=N)
118 elif axis == 1:
119 return np.diff(self.indptr)
120 raise ValueError('axis out of bounds')
122 _getnnz.__doc__ = _spbase._getnnz.__doc__
124 def _set_self(self, other, copy=False):
125 """take the member variables of other and assign them to self"""
127 if copy:
128 other = other.copy()
130 self.data = other.data
131 self.indices = other.indices
132 self.indptr = other.indptr
133 self._shape = check_shape(other.shape)
135 def check_format(self, full_check=True):
136 """check whether the matrix format is valid
138 Parameters
139 ----------
140 full_check : bool, optional
141 If `True`, rigorous check, O(N) operations. Otherwise
142 basic check, O(1) operations (default True).
143 """
144 # use _swap to determine proper bounds
145 major_name, minor_name = self._swap(('row', 'column'))
146 major_dim, minor_dim = self._swap(self.shape)
148 # index arrays should have integer data types
149 if self.indptr.dtype.kind != 'i':
150 warn("indptr array has non-integer dtype ({})"
151 "".format(self.indptr.dtype.name), stacklevel=3)
152 if self.indices.dtype.kind != 'i':
153 warn("indices array has non-integer dtype ({})"
154 "".format(self.indices.dtype.name), stacklevel=3)
156 idx_dtype = self._get_index_dtype((self.indptr, self.indices))
157 self.indptr = np.asarray(self.indptr, dtype=idx_dtype)
158 self.indices = np.asarray(self.indices, dtype=idx_dtype)
159 self.data = to_native(self.data)
161 # check array shapes
162 for x in [self.data.ndim, self.indices.ndim, self.indptr.ndim]:
163 if x != 1:
164 raise ValueError('data, indices, and indptr should be 1-D')
166 # check index pointer
167 if (len(self.indptr) != major_dim + 1):
168 raise ValueError("index pointer size ({}) should be ({})"
169 "".format(len(self.indptr), major_dim + 1))
170 if (self.indptr[0] != 0):
171 raise ValueError("index pointer should start with 0")
173 # check index and data arrays
174 if (len(self.indices) != len(self.data)):
175 raise ValueError("indices and data should have the same size")
176 if (self.indptr[-1] > len(self.indices)):
177 raise ValueError("Last value of index pointer should be less than "
178 "the size of index and data arrays")
180 self.prune()
182 if full_check:
183 # check format validity (more expensive)
184 if self.nnz > 0:
185 if self.indices.max() >= minor_dim:
186 raise ValueError("{} index values must be < {}"
187 "".format(minor_name, minor_dim))
188 if self.indices.min() < 0:
189 raise ValueError("{} index values must be >= 0"
190 "".format(minor_name))
191 if np.diff(self.indptr).min() < 0:
192 raise ValueError("index pointer values must form a "
193 "non-decreasing sequence")
195 # if not self.has_sorted_indices():
196 # warn('Indices were not in sorted order. Sorting indices.')
197 # self.sort_indices()
198 # assert(self.has_sorted_indices())
199 # TODO check for duplicates?
201 #######################
202 # Boolean comparisons #
203 #######################
205 def _scalar_binopt(self, other, op):
206 """Scalar version of self._binopt, for cases in which no new nonzeros
207 are added. Produces a new sparse array in canonical form.
208 """
209 self.sum_duplicates()
210 res = self._with_data(op(self.data, other), copy=True)
211 res.eliminate_zeros()
212 return res
214 def __eq__(self, other):
215 # Scalar other.
216 if isscalarlike(other):
217 if np.isnan(other):
218 return self.__class__(self.shape, dtype=np.bool_)
220 if other == 0:
221 warn("Comparing a sparse matrix with 0 using == is inefficient"
222 ", try using != instead.", SparseEfficiencyWarning,
223 stacklevel=3)
224 all_true = self.__class__(np.ones(self.shape, dtype=np.bool_))
225 inv = self._scalar_binopt(other, operator.ne)
226 return all_true - inv
227 else:
228 return self._scalar_binopt(other, operator.eq)
229 # Dense other.
230 elif isdense(other):
231 return self.todense() == other
232 # Pydata sparse other.
233 elif is_pydata_spmatrix(other):
234 return NotImplemented
235 # Sparse other.
236 elif issparse(other):
237 warn("Comparing sparse matrices using == is inefficient, try using"
238 " != instead.", SparseEfficiencyWarning, stacklevel=3)
239 # TODO sparse broadcasting
240 if self.shape != other.shape:
241 return False
242 elif self.format != other.format:
243 other = other.asformat(self.format)
244 res = self._binopt(other, '_ne_')
245 all_true = self.__class__(np.ones(self.shape, dtype=np.bool_))
246 return all_true - res
247 else:
248 return False
250 def __ne__(self, other):
251 # Scalar other.
252 if isscalarlike(other):
253 if np.isnan(other):
254 warn("Comparing a sparse matrix with nan using != is"
255 " inefficient", SparseEfficiencyWarning, stacklevel=3)
256 all_true = self.__class__(np.ones(self.shape, dtype=np.bool_))
257 return all_true
258 elif other != 0:
259 warn("Comparing a sparse matrix with a nonzero scalar using !="
260 " is inefficient, try using == instead.",
261 SparseEfficiencyWarning, stacklevel=3)
262 all_true = self.__class__(np.ones(self.shape), dtype=np.bool_)
263 inv = self._scalar_binopt(other, operator.eq)
264 return all_true - inv
265 else:
266 return self._scalar_binopt(other, operator.ne)
267 # Dense other.
268 elif isdense(other):
269 return self.todense() != other
270 # Pydata sparse other.
271 elif is_pydata_spmatrix(other):
272 return NotImplemented
273 # Sparse other.
274 elif issparse(other):
275 # TODO sparse broadcasting
276 if self.shape != other.shape:
277 return True
278 elif self.format != other.format:
279 other = other.asformat(self.format)
280 return self._binopt(other, '_ne_')
281 else:
282 return True
284 def _inequality(self, other, op, op_name, bad_scalar_msg):
285 # Scalar other.
286 if isscalarlike(other):
287 if 0 == other and op_name in ('_le_', '_ge_'):
288 raise NotImplementedError(" >= and <= don't work with 0.")
289 elif op(0, other):
290 warn(bad_scalar_msg, SparseEfficiencyWarning)
291 other_arr = np.empty(self.shape, dtype=np.result_type(other))
292 other_arr.fill(other)
293 other_arr = self.__class__(other_arr)
294 return self._binopt(other_arr, op_name)
295 else:
296 return self._scalar_binopt(other, op)
297 # Dense other.
298 elif isdense(other):
299 return op(self.todense(), other)
300 # Sparse other.
301 elif issparse(other):
302 # TODO sparse broadcasting
303 if self.shape != other.shape:
304 raise ValueError("inconsistent shapes")
305 elif self.format != other.format:
306 other = other.asformat(self.format)
307 if op_name not in ('_ge_', '_le_'):
308 return self._binopt(other, op_name)
310 warn("Comparing sparse matrices using >= and <= is inefficient, "
311 "using <, >, or !=, instead.", SparseEfficiencyWarning)
312 all_true = self.__class__(np.ones(self.shape, dtype=np.bool_))
313 res = self._binopt(other, '_gt_' if op_name == '_le_' else '_lt_')
314 return all_true - res
315 else:
316 raise ValueError("Operands could not be compared.")
318 def __lt__(self, other):
319 return self._inequality(other, operator.lt, '_lt_',
320 "Comparing a sparse matrix with a scalar "
321 "greater than zero using < is inefficient, "
322 "try using >= instead.")
324 def __gt__(self, other):
325 return self._inequality(other, operator.gt, '_gt_',
326 "Comparing a sparse matrix with a scalar "
327 "less than zero using > is inefficient, "
328 "try using <= instead.")
330 def __le__(self, other):
331 return self._inequality(other, operator.le, '_le_',
332 "Comparing a sparse matrix with a scalar "
333 "greater than zero using <= is inefficient, "
334 "try using > instead.")
336 def __ge__(self, other):
337 return self._inequality(other, operator.ge, '_ge_',
338 "Comparing a sparse matrix with a scalar "
339 "less than zero using >= is inefficient, "
340 "try using < instead.")
342 #################################
343 # Arithmetic operator overrides #
344 #################################
346 def _add_dense(self, other):
347 if other.shape != self.shape:
348 raise ValueError('Incompatible shapes ({} and {})'
349 .format(self.shape, other.shape))
350 dtype = upcast_char(self.dtype.char, other.dtype.char)
351 order = self._swap('CF')[0]
352 result = np.array(other, dtype=dtype, order=order, copy=True)
353 M, N = self._swap(self.shape)
354 y = result if result.flags.c_contiguous else result.T
355 csr_todense(M, N, self.indptr, self.indices, self.data, y)
356 return self._container(result, copy=False)
358 def _add_sparse(self, other):
359 return self._binopt(other, '_plus_')
361 def _sub_sparse(self, other):
362 return self._binopt(other, '_minus_')
364 def multiply(self, other):
365 """Point-wise multiplication by another matrix, vector, or
366 scalar.
367 """
368 # Scalar multiplication.
369 if isscalarlike(other):
370 return self._mul_scalar(other)
371 # Sparse matrix or vector.
372 if issparse(other):
373 if self.shape == other.shape:
374 other = self.__class__(other)
375 return self._binopt(other, '_elmul_')
376 # Single element.
377 elif other.shape == (1, 1):
378 return self._mul_scalar(other.toarray()[0, 0])
379 elif self.shape == (1, 1):
380 return other._mul_scalar(self.toarray()[0, 0])
381 # A row times a column.
382 elif self.shape[1] == 1 and other.shape[0] == 1:
383 return self._mul_sparse_matrix(other.tocsc())
384 elif self.shape[0] == 1 and other.shape[1] == 1:
385 return other._mul_sparse_matrix(self.tocsc())
386 # Row vector times matrix. other is a row.
387 elif other.shape[0] == 1 and self.shape[1] == other.shape[1]:
388 other = self._dia_container(
389 (other.toarray().ravel(), [0]),
390 shape=(other.shape[1], other.shape[1])
391 )
392 return self._mul_sparse_matrix(other)
393 # self is a row.
394 elif self.shape[0] == 1 and self.shape[1] == other.shape[1]:
395 copy = self._dia_container(
396 (self.toarray().ravel(), [0]),
397 shape=(self.shape[1], self.shape[1])
398 )
399 return other._mul_sparse_matrix(copy)
400 # Column vector times matrix. other is a column.
401 elif other.shape[1] == 1 and self.shape[0] == other.shape[0]:
402 other = self._dia_container(
403 (other.toarray().ravel(), [0]),
404 shape=(other.shape[0], other.shape[0])
405 )
406 return other._mul_sparse_matrix(self)
407 # self is a column.
408 elif self.shape[1] == 1 and self.shape[0] == other.shape[0]:
409 copy = self._dia_container(
410 (self.toarray().ravel(), [0]),
411 shape=(self.shape[0], self.shape[0])
412 )
413 return copy._mul_sparse_matrix(other)
414 else:
415 raise ValueError("inconsistent shapes")
417 # Assume other is a dense matrix/array, which produces a single-item
418 # object array if other isn't convertible to ndarray.
419 other = np.atleast_2d(other)
421 if other.ndim != 2:
422 return np.multiply(self.toarray(), other)
423 # Single element / wrapped object.
424 if other.size == 1:
425 return self._mul_scalar(other.flat[0])
426 # Fast case for trivial sparse matrix.
427 elif self.shape == (1, 1):
428 return np.multiply(self.toarray()[0, 0], other)
430 ret = self.tocoo()
431 # Matching shapes.
432 if self.shape == other.shape:
433 data = np.multiply(ret.data, other[ret.row, ret.col])
434 # Sparse row vector times...
435 elif self.shape[0] == 1:
436 if other.shape[1] == 1: # Dense column vector.
437 data = np.multiply(ret.data, other)
438 elif other.shape[1] == self.shape[1]: # Dense matrix.
439 data = np.multiply(ret.data, other[:, ret.col])
440 else:
441 raise ValueError("inconsistent shapes")
442 row = np.repeat(np.arange(other.shape[0]), len(ret.row))
443 col = np.tile(ret.col, other.shape[0])
444 return self._coo_container(
445 (data.view(np.ndarray).ravel(), (row, col)),
446 shape=(other.shape[0], self.shape[1]),
447 copy=False
448 )
449 # Sparse column vector times...
450 elif self.shape[1] == 1:
451 if other.shape[0] == 1: # Dense row vector.
452 data = np.multiply(ret.data[:, None], other)
453 elif other.shape[0] == self.shape[0]: # Dense matrix.
454 data = np.multiply(ret.data[:, None], other[ret.row])
455 else:
456 raise ValueError("inconsistent shapes")
457 row = np.repeat(ret.row, other.shape[1])
458 col = np.tile(np.arange(other.shape[1]), len(ret.col))
459 return self._coo_container(
460 (data.view(np.ndarray).ravel(), (row, col)),
461 shape=(self.shape[0], other.shape[1]),
462 copy=False
463 )
464 # Sparse matrix times dense row vector.
465 elif other.shape[0] == 1 and self.shape[1] == other.shape[1]:
466 data = np.multiply(ret.data, other[:, ret.col].ravel())
467 # Sparse matrix times dense column vector.
468 elif other.shape[1] == 1 and self.shape[0] == other.shape[0]:
469 data = np.multiply(ret.data, other[ret.row].ravel())
470 else:
471 raise ValueError("inconsistent shapes")
472 ret.data = data.view(np.ndarray).ravel()
473 return ret
475 ###########################
476 # Multiplication handlers #
477 ###########################
479 def _mul_vector(self, other):
480 M, N = self.shape
482 # output array
483 result = np.zeros(M, dtype=upcast_char(self.dtype.char,
484 other.dtype.char))
486 # csr_matvec or csc_matvec
487 fn = getattr(_sparsetools, self.format + '_matvec')
488 fn(M, N, self.indptr, self.indices, self.data, other, result)
490 return result
492 def _mul_multivector(self, other):
493 M, N = self.shape
494 n_vecs = other.shape[1] # number of column vectors
496 result = np.zeros((M, n_vecs),
497 dtype=upcast_char(self.dtype.char, other.dtype.char))
499 # csr_matvecs or csc_matvecs
500 fn = getattr(_sparsetools, self.format + '_matvecs')
501 fn(M, N, n_vecs, self.indptr, self.indices, self.data,
502 other.ravel(), result.ravel())
504 return result
506 def _mul_sparse_matrix(self, other):
507 M, K1 = self.shape
508 K2, N = other.shape
510 major_axis = self._swap((M, N))[0]
511 other = self.__class__(other) # convert to this format
513 idx_dtype = self._get_index_dtype((self.indptr, self.indices,
514 other.indptr, other.indices))
516 fn = getattr(_sparsetools, self.format + '_matmat_maxnnz')
517 nnz = fn(M, N,
518 np.asarray(self.indptr, dtype=idx_dtype),
519 np.asarray(self.indices, dtype=idx_dtype),
520 np.asarray(other.indptr, dtype=idx_dtype),
521 np.asarray(other.indices, dtype=idx_dtype))
523 idx_dtype = self._get_index_dtype((self.indptr, self.indices,
524 other.indptr, other.indices),
525 maxval=nnz)
527 indptr = np.empty(major_axis + 1, dtype=idx_dtype)
528 indices = np.empty(nnz, dtype=idx_dtype)
529 data = np.empty(nnz, dtype=upcast(self.dtype, other.dtype))
531 fn = getattr(_sparsetools, self.format + '_matmat')
532 fn(M, N, np.asarray(self.indptr, dtype=idx_dtype),
533 np.asarray(self.indices, dtype=idx_dtype),
534 self.data,
535 np.asarray(other.indptr, dtype=idx_dtype),
536 np.asarray(other.indices, dtype=idx_dtype),
537 other.data,
538 indptr, indices, data)
540 return self.__class__((data, indices, indptr), shape=(M, N))
542 def diagonal(self, k=0):
543 rows, cols = self.shape
544 if k <= -rows or k >= cols:
545 return np.empty(0, dtype=self.data.dtype)
546 fn = getattr(_sparsetools, self.format + "_diagonal")
547 y = np.empty(min(rows + min(k, 0), cols - max(k, 0)),
548 dtype=upcast(self.dtype))
549 fn(k, self.shape[0], self.shape[1], self.indptr, self.indices,
550 self.data, y)
551 return y
553 diagonal.__doc__ = _spbase.diagonal.__doc__
555 #####################
556 # Other binary ops #
557 #####################
559 def _maximum_minimum(self, other, npop, op_name, dense_check):
560 if isscalarlike(other):
561 if dense_check(other):
562 warn("Taking maximum (minimum) with > 0 (< 0) number results"
563 " to a dense matrix.", SparseEfficiencyWarning,
564 stacklevel=3)
565 other_arr = np.empty(self.shape, dtype=np.asarray(other).dtype)
566 other_arr.fill(other)
567 other_arr = self.__class__(other_arr)
568 return self._binopt(other_arr, op_name)
569 else:
570 self.sum_duplicates()
571 new_data = npop(self.data, np.asarray(other))
572 mat = self.__class__((new_data, self.indices, self.indptr),
573 dtype=new_data.dtype, shape=self.shape)
574 return mat
575 elif isdense(other):
576 return npop(self.todense(), other)
577 elif issparse(other):
578 return self._binopt(other, op_name)
579 else:
580 raise ValueError("Operands not compatible.")
582 def maximum(self, other):
583 return self._maximum_minimum(other, np.maximum,
584 '_maximum_', lambda x: np.asarray(x) > 0)
586 maximum.__doc__ = _spbase.maximum.__doc__
588 def minimum(self, other):
589 return self._maximum_minimum(other, np.minimum,
590 '_minimum_', lambda x: np.asarray(x) < 0)
592 minimum.__doc__ = _spbase.minimum.__doc__
594 #####################
595 # Reduce operations #
596 #####################
598 def sum(self, axis=None, dtype=None, out=None):
599 """Sum the matrix over the given axis. If the axis is None, sum
600 over both rows and columns, returning a scalar.
601 """
602 # The _spbase base class already does axis=0 and axis=1 efficiently
603 # so we only do the case axis=None here
604 if (not hasattr(self, 'blocksize') and
605 axis in self._swap(((1, -1), (0, 2)))[0]):
606 # faster than multiplication for large minor axis in CSC/CSR
607 res_dtype = get_sum_dtype(self.dtype)
608 ret = np.zeros(len(self.indptr) - 1, dtype=res_dtype)
610 major_index, value = self._minor_reduce(np.add)
611 ret[major_index] = value
612 ret = self._ascontainer(ret)
613 if axis % 2 == 1:
614 ret = ret.T
616 if out is not None and out.shape != ret.shape:
617 raise ValueError('dimensions do not match')
619 return ret.sum(axis=(), dtype=dtype, out=out)
620 # _spbase will handle the remaining situations when axis
621 # is in {None, -1, 0, 1}
622 else:
623 return _spbase.sum(self, axis=axis, dtype=dtype, out=out)
625 sum.__doc__ = _spbase.sum.__doc__
627 def _minor_reduce(self, ufunc, data=None):
628 """Reduce nonzeros with a ufunc over the minor axis when non-empty
630 Can be applied to a function of self.data by supplying data parameter.
632 Warning: this does not call sum_duplicates()
634 Returns
635 -------
636 major_index : array of ints
637 Major indices where nonzero
639 value : array of self.dtype
640 Reduce result for nonzeros in each major_index
641 """
642 if data is None:
643 data = self.data
644 major_index = np.flatnonzero(np.diff(self.indptr))
645 value = ufunc.reduceat(data,
646 downcast_intp_index(self.indptr[major_index]))
647 return major_index, value
649 #######################
650 # Getting and Setting #
651 #######################
653 def _get_intXint(self, row, col):
654 M, N = self._swap(self.shape)
655 major, minor = self._swap((row, col))
656 indptr, indices, data = get_csr_submatrix(
657 M, N, self.indptr, self.indices, self.data,
658 major, major + 1, minor, minor + 1)
659 return data.sum(dtype=self.dtype)
661 def _get_sliceXslice(self, row, col):
662 major, minor = self._swap((row, col))
663 if major.step in (1, None) and minor.step in (1, None):
664 return self._get_submatrix(major, minor, copy=True)
665 return self._major_slice(major)._minor_slice(minor)
667 def _get_arrayXarray(self, row, col):
668 # inner indexing
669 idx_dtype = self.indices.dtype
670 M, N = self._swap(self.shape)
671 major, minor = self._swap((row, col))
672 major = np.asarray(major, dtype=idx_dtype)
673 minor = np.asarray(minor, dtype=idx_dtype)
675 val = np.empty(major.size, dtype=self.dtype)
676 csr_sample_values(M, N, self.indptr, self.indices, self.data,
677 major.size, major.ravel(), minor.ravel(), val)
678 if major.ndim == 1:
679 return self._ascontainer(val)
680 return self.__class__(val.reshape(major.shape))
682 def _get_columnXarray(self, row, col):
683 # outer indexing
684 major, minor = self._swap((row, col))
685 return self._major_index_fancy(major)._minor_index_fancy(minor)
687 def _major_index_fancy(self, idx):
688 """Index along the major axis where idx is an array of ints.
689 """
690 idx_dtype = self.indices.dtype
691 indices = np.asarray(idx, dtype=idx_dtype).ravel()
693 _, N = self._swap(self.shape)
694 M = len(indices)
695 new_shape = self._swap((M, N))
696 if M == 0:
697 return self.__class__(new_shape, dtype=self.dtype)
699 row_nnz = self.indptr[indices + 1] - self.indptr[indices]
700 idx_dtype = self.indices.dtype
701 res_indptr = np.zeros(M+1, dtype=idx_dtype)
702 np.cumsum(row_nnz, out=res_indptr[1:])
704 nnz = res_indptr[-1]
705 res_indices = np.empty(nnz, dtype=idx_dtype)
706 res_data = np.empty(nnz, dtype=self.dtype)
707 csr_row_index(M, indices, self.indptr, self.indices, self.data,
708 res_indices, res_data)
710 return self.__class__((res_data, res_indices, res_indptr),
711 shape=new_shape, copy=False)
713 def _major_slice(self, idx, copy=False):
714 """Index along the major axis where idx is a slice object.
715 """
716 if idx == slice(None):
717 return self.copy() if copy else self
719 M, N = self._swap(self.shape)
720 start, stop, step = idx.indices(M)
721 M = len(range(start, stop, step))
722 new_shape = self._swap((M, N))
723 if M == 0:
724 return self.__class__(new_shape, dtype=self.dtype)
726 # Work out what slices are needed for `row_nnz`
727 # start,stop can be -1, only if step is negative
728 start0, stop0 = start, stop
729 if stop == -1 and start >= 0:
730 stop0 = None
731 start1, stop1 = start + 1, stop + 1
733 row_nnz = self.indptr[start1:stop1:step] - \
734 self.indptr[start0:stop0:step]
735 idx_dtype = self.indices.dtype
736 res_indptr = np.zeros(M+1, dtype=idx_dtype)
737 np.cumsum(row_nnz, out=res_indptr[1:])
739 if step == 1:
740 all_idx = slice(self.indptr[start], self.indptr[stop])
741 res_indices = np.array(self.indices[all_idx], copy=copy)
742 res_data = np.array(self.data[all_idx], copy=copy)
743 else:
744 nnz = res_indptr[-1]
745 res_indices = np.empty(nnz, dtype=idx_dtype)
746 res_data = np.empty(nnz, dtype=self.dtype)
747 csr_row_slice(start, stop, step, self.indptr, self.indices,
748 self.data, res_indices, res_data)
750 return self.__class__((res_data, res_indices, res_indptr),
751 shape=new_shape, copy=False)
753 def _minor_index_fancy(self, idx):
754 """Index along the minor axis where idx is an array of ints.
755 """
756 idx_dtype = self.indices.dtype
757 idx = np.asarray(idx, dtype=idx_dtype).ravel()
759 M, N = self._swap(self.shape)
760 k = len(idx)
761 new_shape = self._swap((M, k))
762 if k == 0:
763 return self.__class__(new_shape, dtype=self.dtype)
765 # pass 1: count idx entries and compute new indptr
766 col_offsets = np.zeros(N, dtype=idx_dtype)
767 res_indptr = np.empty_like(self.indptr)
768 csr_column_index1(k, idx, M, N, self.indptr, self.indices,
769 col_offsets, res_indptr)
771 # pass 2: copy indices/data for selected idxs
772 col_order = np.argsort(idx).astype(idx_dtype, copy=False)
773 nnz = res_indptr[-1]
774 res_indices = np.empty(nnz, dtype=idx_dtype)
775 res_data = np.empty(nnz, dtype=self.dtype)
776 csr_column_index2(col_order, col_offsets, len(self.indices),
777 self.indices, self.data, res_indices, res_data)
778 return self.__class__((res_data, res_indices, res_indptr),
779 shape=new_shape, copy=False)
781 def _minor_slice(self, idx, copy=False):
782 """Index along the minor axis where idx is a slice object.
783 """
784 if idx == slice(None):
785 return self.copy() if copy else self
787 M, N = self._swap(self.shape)
788 start, stop, step = idx.indices(N)
789 N = len(range(start, stop, step))
790 if N == 0:
791 return self.__class__(self._swap((M, N)), dtype=self.dtype)
792 if step == 1:
793 return self._get_submatrix(minor=idx, copy=copy)
794 # TODO: don't fall back to fancy indexing here
795 return self._minor_index_fancy(np.arange(start, stop, step))
797 def _get_submatrix(self, major=None, minor=None, copy=False):
798 """Return a submatrix of this matrix.
800 major, minor: None, int, or slice with step 1
801 """
802 M, N = self._swap(self.shape)
803 i0, i1 = _process_slice(major, M)
804 j0, j1 = _process_slice(minor, N)
806 if i0 == 0 and j0 == 0 and i1 == M and j1 == N:
807 return self.copy() if copy else self
809 indptr, indices, data = get_csr_submatrix(
810 M, N, self.indptr, self.indices, self.data, i0, i1, j0, j1)
812 shape = self._swap((i1 - i0, j1 - j0))
813 return self.__class__((data, indices, indptr), shape=shape,
814 dtype=self.dtype, copy=False)
816 def _set_intXint(self, row, col, x):
817 i, j = self._swap((row, col))
818 self._set_many(i, j, x)
820 def _set_arrayXarray(self, row, col, x):
821 i, j = self._swap((row, col))
822 self._set_many(i, j, x)
824 def _set_arrayXarray_sparse(self, row, col, x):
825 # clear entries that will be overwritten
826 self._zero_many(*self._swap((row, col)))
828 M, N = row.shape # matches col.shape
829 broadcast_row = M != 1 and x.shape[0] == 1
830 broadcast_col = N != 1 and x.shape[1] == 1
831 r, c = x.row, x.col
833 x = np.asarray(x.data, dtype=self.dtype)
834 if x.size == 0:
835 return
837 if broadcast_row:
838 r = np.repeat(np.arange(M), len(r))
839 c = np.tile(c, M)
840 x = np.tile(x, M)
841 if broadcast_col:
842 r = np.repeat(r, N)
843 c = np.tile(np.arange(N), len(c))
844 x = np.repeat(x, N)
845 # only assign entries in the new sparsity structure
846 i, j = self._swap((row[r, c], col[r, c]))
847 self._set_many(i, j, x)
849 def _setdiag(self, values, k):
850 if 0 in self.shape:
851 return
853 M, N = self.shape
854 broadcast = (values.ndim == 0)
856 if k < 0:
857 if broadcast:
858 max_index = min(M + k, N)
859 else:
860 max_index = min(M + k, N, len(values))
861 i = np.arange(max_index, dtype=self.indices.dtype)
862 j = np.arange(max_index, dtype=self.indices.dtype)
863 i -= k
865 else:
866 if broadcast:
867 max_index = min(M, N - k)
868 else:
869 max_index = min(M, N - k, len(values))
870 i = np.arange(max_index, dtype=self.indices.dtype)
871 j = np.arange(max_index, dtype=self.indices.dtype)
872 j += k
874 if not broadcast:
875 values = values[:len(i)]
877 self[i, j] = values
879 def _prepare_indices(self, i, j):
880 M, N = self._swap(self.shape)
882 def check_bounds(indices, bound):
883 idx = indices.max()
884 if idx >= bound:
885 raise IndexError('index (%d) out of range (>= %d)' %
886 (idx, bound))
887 idx = indices.min()
888 if idx < -bound:
889 raise IndexError('index (%d) out of range (< -%d)' %
890 (idx, bound))
892 i = np.array(i, dtype=self.indices.dtype, copy=False, ndmin=1).ravel()
893 j = np.array(j, dtype=self.indices.dtype, copy=False, ndmin=1).ravel()
894 check_bounds(i, M)
895 check_bounds(j, N)
896 return i, j, M, N
898 def _set_many(self, i, j, x):
899 """Sets value at each (i, j) to x
901 Here (i,j) index major and minor respectively, and must not contain
902 duplicate entries.
903 """
904 i, j, M, N = self._prepare_indices(i, j)
905 x = np.array(x, dtype=self.dtype, copy=False, ndmin=1).ravel()
907 n_samples = x.size
908 offsets = np.empty(n_samples, dtype=self.indices.dtype)
909 ret = csr_sample_offsets(M, N, self.indptr, self.indices, n_samples,
910 i, j, offsets)
911 if ret == 1:
912 # rinse and repeat
913 self.sum_duplicates()
914 csr_sample_offsets(M, N, self.indptr, self.indices, n_samples,
915 i, j, offsets)
917 if -1 not in offsets:
918 # only affects existing non-zero cells
919 self.data[offsets] = x
920 return
922 else:
923 warn("Changing the sparsity structure of a {}_matrix is expensive."
924 " lil_matrix is more efficient.".format(self.format),
925 SparseEfficiencyWarning, stacklevel=3)
926 # replace where possible
927 mask = offsets > -1
928 self.data[offsets[mask]] = x[mask]
929 # only insertions remain
930 mask = ~mask
931 i = i[mask]
932 i[i < 0] += M
933 j = j[mask]
934 j[j < 0] += N
935 self._insert_many(i, j, x[mask])
937 def _zero_many(self, i, j):
938 """Sets value at each (i, j) to zero, preserving sparsity structure.
940 Here (i,j) index major and minor respectively.
941 """
942 i, j, M, N = self._prepare_indices(i, j)
944 n_samples = len(i)
945 offsets = np.empty(n_samples, dtype=self.indices.dtype)
946 ret = csr_sample_offsets(M, N, self.indptr, self.indices, n_samples,
947 i, j, offsets)
948 if ret == 1:
949 # rinse and repeat
950 self.sum_duplicates()
951 csr_sample_offsets(M, N, self.indptr, self.indices, n_samples,
952 i, j, offsets)
954 # only assign zeros to the existing sparsity structure
955 self.data[offsets[offsets > -1]] = 0
957 def _insert_many(self, i, j, x):
958 """Inserts new nonzero at each (i, j) with value x
960 Here (i,j) index major and minor respectively.
961 i, j and x must be non-empty, 1d arrays.
962 Inserts each major group (e.g. all entries per row) at a time.
963 Maintains has_sorted_indices property.
964 Modifies i, j, x in place.
965 """
966 order = np.argsort(i, kind='mergesort') # stable for duplicates
967 i = i.take(order, mode='clip')
968 j = j.take(order, mode='clip')
969 x = x.take(order, mode='clip')
971 do_sort = self.has_sorted_indices
973 # Update index data type
974 idx_dtype = self._get_index_dtype((self.indices, self.indptr),
975 maxval=(self.indptr[-1] + x.size))
976 self.indptr = np.asarray(self.indptr, dtype=idx_dtype)
977 self.indices = np.asarray(self.indices, dtype=idx_dtype)
978 i = np.asarray(i, dtype=idx_dtype)
979 j = np.asarray(j, dtype=idx_dtype)
981 # Collate old and new in chunks by major index
982 indices_parts = []
983 data_parts = []
984 ui, ui_indptr = np.unique(i, return_index=True)
985 ui_indptr = np.append(ui_indptr, len(j))
986 new_nnzs = np.diff(ui_indptr)
987 prev = 0
988 for c, (ii, js, je) in enumerate(zip(ui, ui_indptr, ui_indptr[1:])):
989 # old entries
990 start = self.indptr[prev]
991 stop = self.indptr[ii]
992 indices_parts.append(self.indices[start:stop])
993 data_parts.append(self.data[start:stop])
995 # handle duplicate j: keep last setting
996 uj, uj_indptr = np.unique(j[js:je][::-1], return_index=True)
997 if len(uj) == je - js:
998 indices_parts.append(j[js:je])
999 data_parts.append(x[js:je])
1000 else:
1001 indices_parts.append(j[js:je][::-1][uj_indptr])
1002 data_parts.append(x[js:je][::-1][uj_indptr])
1003 new_nnzs[c] = len(uj)
1005 prev = ii
1007 # remaining old entries
1008 start = self.indptr[ii]
1009 indices_parts.append(self.indices[start:])
1010 data_parts.append(self.data[start:])
1012 # update attributes
1013 self.indices = np.concatenate(indices_parts)
1014 self.data = np.concatenate(data_parts)
1015 nnzs = np.empty(self.indptr.shape, dtype=idx_dtype)
1016 nnzs[0] = idx_dtype(0)
1017 indptr_diff = np.diff(self.indptr)
1018 indptr_diff[ui] += new_nnzs
1019 nnzs[1:] = indptr_diff
1020 self.indptr = np.cumsum(nnzs, out=nnzs)
1022 if do_sort:
1023 # TODO: only sort where necessary
1024 self.has_sorted_indices = False
1025 self.sort_indices()
1027 self.check_format(full_check=False)
1029 ######################
1030 # Conversion methods #
1031 ######################
1033 def tocoo(self, copy=True):
1034 major_dim, minor_dim = self._swap(self.shape)
1035 minor_indices = self.indices
1036 major_indices = np.empty(len(minor_indices), dtype=self.indices.dtype)
1037 _sparsetools.expandptr(major_dim, self.indptr, major_indices)
1038 row, col = self._swap((major_indices, minor_indices))
1040 return self._coo_container(
1041 (self.data, (row, col)), self.shape, copy=copy,
1042 dtype=self.dtype
1043 )
1045 tocoo.__doc__ = _spbase.tocoo.__doc__
1047 def toarray(self, order=None, out=None):
1048 if out is None and order is None:
1049 order = self._swap('cf')[0]
1050 out = self._process_toarray_args(order, out)
1051 if not (out.flags.c_contiguous or out.flags.f_contiguous):
1052 raise ValueError('Output array must be C or F contiguous')
1053 # align ideal order with output array order
1054 if out.flags.c_contiguous:
1055 x = self.tocsr()
1056 y = out
1057 else:
1058 x = self.tocsc()
1059 y = out.T
1060 M, N = x._swap(x.shape)
1061 csr_todense(M, N, x.indptr, x.indices, x.data, y)
1062 return out
1064 toarray.__doc__ = _spbase.toarray.__doc__
1066 ##############################################################
1067 # methods that examine or modify the internal data structure #
1068 ##############################################################
1070 def eliminate_zeros(self):
1071 """Remove zero entries from the matrix
1073 This is an *in place* operation.
1074 """
1075 M, N = self._swap(self.shape)
1076 _sparsetools.csr_eliminate_zeros(M, N, self.indptr, self.indices,
1077 self.data)
1078 self.prune() # nnz may have changed
1080 def __get_has_canonical_format(self):
1081 """Determine whether the matrix has sorted indices and no duplicates
1083 Returns
1084 - True: if the above applies
1085 - False: otherwise
1087 has_canonical_format implies has_sorted_indices, so if the latter flag
1088 is False, so will the former be; if the former is found True, the
1089 latter flag is also set.
1090 """
1092 # first check to see if result was cached
1093 if not getattr(self, '_has_sorted_indices', True):
1094 # not sorted => not canonical
1095 self._has_canonical_format = False
1096 elif not hasattr(self, '_has_canonical_format'):
1097 self.has_canonical_format = bool(
1098 _sparsetools.csr_has_canonical_format(
1099 len(self.indptr) - 1, self.indptr, self.indices))
1100 return self._has_canonical_format
1102 def __set_has_canonical_format(self, val):
1103 self._has_canonical_format = bool(val)
1104 if val:
1105 self.has_sorted_indices = True
1107 has_canonical_format = property(fget=__get_has_canonical_format,
1108 fset=__set_has_canonical_format)
1110 def sum_duplicates(self):
1111 """Eliminate duplicate matrix entries by adding them together
1113 This is an *in place* operation.
1114 """
1115 if self.has_canonical_format:
1116 return
1117 self.sort_indices()
1119 M, N = self._swap(self.shape)
1120 _sparsetools.csr_sum_duplicates(M, N, self.indptr, self.indices,
1121 self.data)
1123 self.prune() # nnz may have changed
1124 self.has_canonical_format = True
1126 def __get_sorted(self):
1127 """Determine whether the matrix has sorted indices
1129 Returns
1130 - True: if the indices of the matrix are in sorted order
1131 - False: otherwise
1133 """
1135 # first check to see if result was cached
1136 if not hasattr(self, '_has_sorted_indices'):
1137 self._has_sorted_indices = bool(
1138 _sparsetools.csr_has_sorted_indices(
1139 len(self.indptr) - 1, self.indptr, self.indices))
1140 return self._has_sorted_indices
1142 def __set_sorted(self, val):
1143 self._has_sorted_indices = bool(val)
1145 has_sorted_indices = property(fget=__get_sorted, fset=__set_sorted)
1147 def sorted_indices(self):
1148 """Return a copy of this matrix with sorted indices
1149 """
1150 A = self.copy()
1151 A.sort_indices()
1152 return A
1154 # an alternative that has linear complexity is the following
1155 # although the previous option is typically faster
1156 # return self.toother().toother()
1158 def sort_indices(self):
1159 """Sort the indices of this matrix *in place*
1160 """
1162 if not self.has_sorted_indices:
1163 _sparsetools.csr_sort_indices(len(self.indptr) - 1, self.indptr,
1164 self.indices, self.data)
1165 self.has_sorted_indices = True
1167 def prune(self):
1168 """Remove empty space after all non-zero elements.
1169 """
1170 major_dim = self._swap(self.shape)[0]
1172 if len(self.indptr) != major_dim + 1:
1173 raise ValueError('index pointer has invalid length')
1174 if len(self.indices) < self.nnz:
1175 raise ValueError('indices array has fewer than nnz elements')
1176 if len(self.data) < self.nnz:
1177 raise ValueError('data array has fewer than nnz elements')
1179 self.indices = _prune_array(self.indices[:self.nnz])
1180 self.data = _prune_array(self.data[:self.nnz])
1182 def resize(self, *shape):
1183 shape = check_shape(shape)
1184 if hasattr(self, 'blocksize'):
1185 bm, bn = self.blocksize
1186 new_M, rm = divmod(shape[0], bm)
1187 new_N, rn = divmod(shape[1], bn)
1188 if rm or rn:
1189 raise ValueError("shape must be divisible into {} blocks. "
1190 "Got {}".format(self.blocksize, shape))
1191 M, N = self.shape[0] // bm, self.shape[1] // bn
1192 else:
1193 new_M, new_N = self._swap(shape)
1194 M, N = self._swap(self.shape)
1196 if new_M < M:
1197 self.indices = self.indices[:self.indptr[new_M]]
1198 self.data = self.data[:self.indptr[new_M]]
1199 self.indptr = self.indptr[:new_M + 1]
1200 elif new_M > M:
1201 self.indptr = np.resize(self.indptr, new_M + 1)
1202 self.indptr[M + 1:].fill(self.indptr[M])
1204 if new_N < N:
1205 mask = self.indices < new_N
1206 if not np.all(mask):
1207 self.indices = self.indices[mask]
1208 self.data = self.data[mask]
1209 major_index, val = self._minor_reduce(np.add, mask)
1210 self.indptr.fill(0)
1211 self.indptr[1:][major_index] = val
1212 np.cumsum(self.indptr, out=self.indptr)
1214 self._shape = shape
1216 resize.__doc__ = _spbase.resize.__doc__
1218 ###################
1219 # utility methods #
1220 ###################
1222 # needed by _data_matrix
1223 def _with_data(self, data, copy=True):
1224 """Returns a matrix with the same sparsity structure as self,
1225 but with different data. By default the structure arrays
1226 (i.e. .indptr and .indices) are copied.
1227 """
1228 if copy:
1229 return self.__class__((data, self.indices.copy(),
1230 self.indptr.copy()),
1231 shape=self.shape,
1232 dtype=data.dtype)
1233 else:
1234 return self.__class__((data, self.indices, self.indptr),
1235 shape=self.shape, dtype=data.dtype)
1237 def _binopt(self, other, op):
1238 """apply the binary operation fn to two sparse matrices."""
1239 other = self.__class__(other)
1241 # e.g. csr_plus_csr, csr_minus_csr, etc.
1242 fn = getattr(_sparsetools, self.format + op + self.format)
1244 maxnnz = self.nnz + other.nnz
1245 idx_dtype = self._get_index_dtype((self.indptr, self.indices,
1246 other.indptr, other.indices),
1247 maxval=maxnnz)
1248 indptr = np.empty(self.indptr.shape, dtype=idx_dtype)
1249 indices = np.empty(maxnnz, dtype=idx_dtype)
1251 bool_ops = ['_ne_', '_lt_', '_gt_', '_le_', '_ge_']
1252 if op in bool_ops:
1253 data = np.empty(maxnnz, dtype=np.bool_)
1254 else:
1255 data = np.empty(maxnnz, dtype=upcast(self.dtype, other.dtype))
1257 fn(self.shape[0], self.shape[1],
1258 np.asarray(self.indptr, dtype=idx_dtype),
1259 np.asarray(self.indices, dtype=idx_dtype),
1260 self.data,
1261 np.asarray(other.indptr, dtype=idx_dtype),
1262 np.asarray(other.indices, dtype=idx_dtype),
1263 other.data,
1264 indptr, indices, data)
1266 A = self.__class__((data, indices, indptr), shape=self.shape)
1267 A.prune()
1269 return A
1271 def _divide_sparse(self, other):
1272 """
1273 Divide this matrix by a second sparse matrix.
1274 """
1275 if other.shape != self.shape:
1276 raise ValueError('inconsistent shapes')
1278 r = self._binopt(other, '_eldiv_')
1280 if np.issubdtype(r.dtype, np.inexact):
1281 # Eldiv leaves entries outside the combined sparsity
1282 # pattern empty, so they must be filled manually.
1283 # Everything outside of other's sparsity is NaN, and everything
1284 # inside it is either zero or defined by eldiv.
1285 out = np.empty(self.shape, dtype=self.dtype)
1286 out.fill(np.nan)
1287 row, col = other.nonzero()
1288 out[row, col] = 0
1289 r = r.tocoo()
1290 out[r.row, r.col] = r.data
1291 out = self._container(out)
1292 else:
1293 # integers types go with nan <-> 0
1294 out = r
1296 return out
1299def _process_slice(sl, num):
1300 if sl is None:
1301 i0, i1 = 0, num
1302 elif isinstance(sl, slice):
1303 i0, i1, stride = sl.indices(num)
1304 if stride != 1:
1305 raise ValueError('slicing with step != 1 not supported')
1306 i0 = min(i0, i1) # give an empty slice when i0 > i1
1307 elif isintlike(sl):
1308 if sl < 0:
1309 sl += num
1310 i0, i1 = sl, sl + 1
1311 if i0 < 0 or i1 > num:
1312 raise IndexError('index out of bounds: 0 <= %d < %d <= %d' %
1313 (i0, i1, num))
1314 else:
1315 raise TypeError('expected slice or scalar')
1317 return i0, i1