diff --git a/libpysal/weights/raster.py b/libpysal/weights/raster.py index ee373609d..973a9ef17 100644 --- a/libpysal/weights/raster.py +++ b/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) diff --git a/libpysal/weights/set_operations.py b/libpysal/weights/set_operations.py index 2bf3d6037..2b72dcf7f 100644 --- a/libpysal/weights/set_operations.py +++ b/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) diff --git a/libpysal/weights/tests/test_util.py b/libpysal/weights/tests/test_util.py index 9d1b7208e..247b5aa10 100644 --- a/libpysal/weights/tests/test_util.py +++ b/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] diff --git a/libpysal/weights/util.py b/libpysal/weights/util.py index 6154387f1..e436035a9 100644 --- a/libpysal/weights/util.py +++ b/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 diff --git a/libpysal/weights/weights.py b/libpysal/weights/weights.py index 8e582953f..9c125f1d7 100644 --- a/libpysal/weights/weights.py +++ b/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