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

Author: dschult

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(sparse.linalg.matrix_power(sw, x) for x in range(1, k + 1))
+        #### Could be as follows after scipy >=1.12 is required
+        # 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/util.py

@@ -523,17 +523,42 @@ def higher_order_sp(
         )
 
     if lower_order:
-        wk = sum(sparse.linalg.matrix_power(w, k) for k in range(1, k+1))
         shortest_path = False
+        #### Could be as follows after scipy >=1.12 is required
+        # 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 = sparse.linalg.matrix_power(w, k)
+        #### Could be as follows after scipy >=1.12 is required
+        # 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 = sparse.linalg.matrix_power(w, j)
+            #### Could be as follows after scipy >=1.12 is required
+            # 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)