Exploring SciPy sparse array migration from sparse matrices

libpysal/weights/raster.py

@@ -270,7 +270,13 @@ 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))
+        #### Can be this one-liner after scipy >=1.12 is assured
+        # sw = sum(sparse.linalg.matrix_power(sw, x) for x in range(1, k + 1))
+        tmp = sw.copy()
+        for _ in range(k - 1):
+            tmp += tmp @ sw + sw
+        sw = tmp
+        ####
         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,42 @@ def higher_order_sp(
         )
 
     if lower_order:
-        wk = sum(w**x for x in range(1, k + 1))
         shortest_path = False
+        #### Can be this one-liner after scipy >=1.12 is assured
+        # wk = sum(sparse.linalg.matrix_power(w, k) for k in range(1, k+1))
+        wk = w.copy()
+        for _ in range(k - 1):
+            wk = wk @ w + w
+        ####
     else:
-        wk = w**k
+        #### Can be this one-liner after scipy >=1.12 is assured
+        # wk = sparse.linalg.matrix_power(w, k)
+        wk = w.copy()
+        x = 1
+        while 2 * x < k:
+            wk = wk @ wk
+            x *= 2
+        while x < k:
+            wk = wk @ w
+            x += 1
+        ####
 
     rk, ck = wk.nonzero()
     sk = set(zip(rk, ck, strict=True))
 
     if shortest_path:
         for j in range(1, k):
-            wj = w**j
+            #### Can be this one-liner after scipy >=1.12 is assured
+            # wj = sparse.linalg.matrix_power(w, j)
+            wj = w.copy()
+            x = 1
+            while 2 * x < j:
+                wj = wj @ wj
+                x *= 2
+            while x < j:
+                wj = wj @ w
+                x += 1
+            ####
             rj, cj = wj.nonzero()
             sj = set(zip(rj, cj, strict=True))
             sk.difference_update(sj)
@@ -1225,7 +1250,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