First part of SciPy sparse array migration from sparse matrices (#785)

* migration guide pass 1: make code work for both spmatrix and sparray

switch from isspmatrix_csr to issparse and format='csr'

* add conftest.txt to monkey-patch scipy sparse to fail on breaking * and **

* add test of row_st code to help coverage bot

* workaround scipy1.8 limitation. Could be reverted for scipy1.12

* Clean up comments and remove conftest.txt -- PR ready to remove Draft status

* bump scipy to 1.12

* cleanup

---------

Co-authored-by: Martin Fleischmann <martin@martinfleischmann.net>

Author: dschult

ci/310-oldest.yaml

@@ -10,7 +10,7 @@ dependencies:
   - packaging=22
   - pandas=1.4
   - requests=2.27
-  - scipy=1.8
+  - scipy=1.12
   - shapely=2.0.1
   # testing
   - codecov

libpysal/graph/base.py

@@ -3,8 +3,6 @@
 
 import numpy as np
 import pandas as pd
-from packaging.version import Version
-from scipy import __version__ as scipy_version
 from scipy import sparse
 
 from libpysal.weights import W
@@ -2000,9 +1998,6 @@ def higher_order(self, k=2, shortest_path=True, diagonal=False, lower_order=Fals
         <Graph of 85 nodes and 1176 nonzero edges indexed by
          [0, 1, 2, 3, 4, ...]>
         """
-        if not Version(scipy_version) >= Version("1.12.0"):
-            raise ImportError("Graph.higher_order() requires scipy>=1.12.0.")
-
         binary = self.transform("B")
         sp = binary.sparse
 

libpysal/graph/tests/test_base.py

@@ -7,8 +7,6 @@
 import numpy as np
 import pandas as pd
 import pytest
-from packaging.version import Version
-from scipy import __version__ as scipy_version
 from scipy import sparse
 
 from libpysal import graph, weights
@@ -1014,10 +1012,6 @@ def test_lag(self):
         with pytest.raises(ValueError, match="The length of `y`"):
             self.g_str.lag(list(range(1, 15)))
 
-    @pytest.mark.skipif(
-        Version(scipy_version) < Version("1.12.0"),
-        reason="sparse matrix power requires scipy>=1.12.0",
-    )
     def test_higher_order(self):
         cont = graph.Graph.build_contiguity(self.nybb)
         k2 = cont.higher_order(2)
@@ -1091,10 +1085,6 @@ def test_higher_order(self):
         lower = cont.higher_order(2, lower_order=True)
         assert lower == expected
 
-    @pytest.mark.skipif(
-        Version(scipy_version) < Version("1.12.0"),
-        reason="sparse matrix power requires scipy>=1.12.0",
-    )
     def test_higher_order_inclusive(self):  # GH738
         contig = graph.Graph.from_arrays(
             [0, 1, 2, 3, 3, 4, 4], [0, 3, 4, 1, 4, 2, 3], [0, 1, 1, 1, 1, 1, 1]

libpysal/graph/tests/test_raster.py

@@ -1,8 +1,6 @@
 import numpy as np
 import pandas as pd
 import pytest
-from packaging.version import Version
-from scipy import __version__ as scipy_version
 
 from libpysal import graph
 from libpysal.weights.raster import testDataArray as dummy_array  # noqa: N813
@@ -18,10 +16,6 @@ def setup_method(self):
         self.da4 = dummy_array((1, 1), missing_vals=False)
         self.da4.data = np.array([["test"]])
 
-    @pytest.mark.skipif(
-        Version(scipy_version) < Version("1.12.0"),
-        reason="sparse matrix power requires scipy>=1.12.0",
-    )
     def test_queen(self):
         g1 = graph.Graph.build_raster_contiguity(self.da1, rook=False, k=2, n_jobs=-1)
         assert g1[(1, -30.0, -180.0)].to_dict() == {

libpysal/weights/raster.py

@@ -270,7 +270,7 @@ def da2WSP(
     # then eliminate zeros from the data. This changes the
     # sparsity of the csr_matrix !!
     if k > 1 and not include_nodata:
-        sw = sum(sw**x for x in range(1, k + 1))
+        sw = sum(sparse.linalg.matrix_power(sw, x) for x in range(1, k + 1))
         sw.setdiag(0)
         sw.eliminate_zeros()
         sw.data[:] = np.ones_like(sw.data, dtype=np.int8)

libpysal/weights/set_operations.py

@@ -12,7 +12,7 @@
 import copy
 
 from numpy import ones
-from scipy.sparse import isspmatrix_csr
+from scipy.sparse import issparse
 
 from .weights import WSP, W
 
@@ -501,9 +501,9 @@ def w_clip(w1, w2, outSP=True, **kwargs):  # noqa: N803
     if not w1.id_order:
         w1.id_order = None
     id_order = w1.id_order
-    if not isspmatrix_csr(w1):
+    if not issparse(w1):
         w1 = w1.sparse
-    if not isspmatrix_csr(w2):
+    if not issparse(w2):
         w2 = w2.sparse
     w2.data = ones(w2.data.shape)
     wc = w1.multiply(w2)

libpysal/weights/tests/test_util.py

@@ -25,13 +25,21 @@ def test_lat2_w(self):
         assert w9[0] == {1: 1.0, 3: 1.0}
         assert w9[3] == {0: 1.0, 4: 1.0, 6: 1.0}
 
-    def test_lat2_sw(self):
-        w9 = util.lat2SW(3, 3)
+    @pytest.mark.parametrize("row_st", [True, False])
+    def test_lat2_sw(self, row_st):
+        w9 = util.lat2SW(3, 3, row_st=row_st)
         rows, cols = w9.shape
         n = rows * cols
         assert w9.nnz == 24
         pct_nonzero = w9.nnz / float(n)
         assert pct_nonzero == 0.29629629629629628
+        if row_st:
+            #### Can be this one-liner after scipy >=1.15 is assured
+            # w9 = (w9.T.multiply(w9.tocsr().count_nonzero(axis=1))).T
+            w9csr = w9.tocsr()
+            nnz_by_axis1 = np.diff(w9csr.indptr)
+            w9 = (w9csr.T.multiply(nnz_by_axis1)).T
+            ####
         data = w9.todense().tolist()
         assert data[0] == [0, 1, 0, 1, 0, 0, 0, 0, 0]
         assert data[1] == [1, 0, 1, 0, 1, 0, 0, 0, 0]

libpysal/weights/util.py

@@ -511,7 +511,7 @@ def higher_order_sp(
             w = w.sparse
         else:
             raise ValueError("Weights are not binary (0,1)")
-    elif scipy.sparse.isspmatrix_csr(w):
+    elif scipy.sparse.issparse(w) and w.format == "csr":
         if not np.unique(w.data) == np.array([1.0]):
             raise ValueError(
                 "Sparse weights matrix is not binary (0,1) weights matrix."
@@ -523,17 +523,17 @@ def higher_order_sp(
         )
 
     if lower_order:
-        wk = sum(w**x for x in range(1, k + 1))
         shortest_path = False
+        wk = sum(sparse.linalg.matrix_power(w, k) for k in range(1, k + 1))
     else:
-        wk = w**k
+        wk = sparse.linalg.matrix_power(w, k)
 
     rk, ck = wk.nonzero()
     sk = set(zip(rk, ck, strict=True))
 
     if shortest_path:
         for j in range(1, k):
-            wj = w**j
+            wj = sparse.linalg.matrix_power(w, j)
             rj, cj = wj.nonzero()
             sj = set(zip(rj, cj, strict=True))
             sk.difference_update(sj)
@@ -1225,7 +1225,7 @@ def lat2SW(nrows=3, ncols=5, criterion="rook", row_st=False):
     m = sparse.dia_matrix((data, offsets), shape=(n, n), dtype=np.int8)
     m = m + m.T
     if row_st:
-        m = sparse.spdiags(1.0 / m.sum(1).T, 0, *m.shape) * m
+        m = sparse.dia_matrix(((1.0 / m.sum(1).T), [0]), shape=m.shape) @ m
     m = m.tocsc()
     m.eliminate_zeros()
     return m

libpysal/weights/weights.py

@@ -721,7 +721,7 @@ def diagW2(self):
         trcW2
         """
         if "diagw2" not in self._cache:
-            self._diagW2 = (self.sparse * self.sparse).diagonal()
+            self._diagW2 = (self.sparse @ self.sparse).diagonal()
             self._cache["diagW2"] = self._diagW2
         return self._diagW2
 
@@ -734,7 +734,7 @@ def diagWtW(self):
         trcWtW
         """
         if "diagWtW" not in self._cache:
-            self._diagWtW = (self.sparse.transpose() * self.sparse).diagonal()
+            self._diagWtW = (self.sparse.transpose() @ self.sparse).diagonal()
             self._cache["diagWtW"] = self._diagWtW
         return self._diagWtW
 
@@ -757,7 +757,7 @@ def diagWtW_WW(self):
         if "diagWtW_WW" not in self._cache:
             wt = self.sparse.transpose()
             w = self.sparse
-            self._diagWtW_WW = (wt * w + w * w).diagonal()
+            self._diagWtW_WW = (wt @ w + w @ w).diagonal()
             self._cache["diagWtW_WW"] = self._diagWtW_WW
         return self._diagWtW_WW
 
@@ -1603,7 +1603,7 @@ def diagWtW_WW(self):
         if "diagWtW_WW" not in self._cache:
             wt = self.sparse.transpose()
             w = self.sparse
-            self._diagWtW_WW = (wt * w + w * w).diagonal()
+            self._diagWtW_WW = (wt @ w + w @ w).diagonal()
             self._cache["diagWtW_WW"] = self._diagWtW_WW
         return self._diagWtW_WW
 

pyproject.toml

@@ -32,7 +32,7 @@ dependencies = [
     "pandas>=1.4",
     "platformdirs>=2.0.2",
     "requests>=2.27",
-    "scipy>=1.8",
+    "scipy>=1.12",
     "shapely>=2.0.1",
     "scikit-learn>=1.1",
 ]