Skip to content

Commit 18f369c

Browse files
committed
update to graph-compatible code
1 parent 9ca259d commit 18f369c

File tree

2 files changed

+33
-19
lines changed

2 files changed

+33
-19
lines changed

esda/moran_local_mv.py

Lines changed: 10 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -23,7 +23,7 @@ def __init__(
2323
y : (N,1) array
2424
array of data that is the targeted "outcome" covariate
2525
to compute the multivariable Moran's I
26-
X : (N,3) array
26+
X : (N,p) array
2727
array of data that is used as "confounding factors"
2828
to account for their covariance with Y.
2929
W : (N,N) weights object
@@ -83,7 +83,7 @@ def __init__(
8383
if isinstance(W, Graph):
8484
W = W.transform("R")
8585
else:
86-
W.transform = "r" # TODO: as a function for graph
86+
W.transform = "r"
8787
y = y - y.mean()
8888
if unit_scale:
8989
y /= y.std()
@@ -154,7 +154,10 @@ def __init__(
154154
self.quads_ = self._uvquads_[:, 1]
155155

156156
def _make_data(self, z, X, W):
157-
Wz = lag_spatial(W, z)
157+
if isinstance(W, Graph):
158+
Wz = W.lag(z)
159+
else:
160+
Wz = lag_spatial(W, z)
158161
if X is not None:
159162
D = np.hstack((np.ones(z.shape), z, X))
160163
P = X.shape[1] + 1
@@ -368,7 +371,10 @@ def _crand(self):
368371
lisas = np.zeros((self.connectivity.n, self.permutations))
369372
n_1 = self.connectivity.n - 1
370373
prange = list(range(self.permutations))
371-
k = self.connectivity.cardinalities.max() + 1
374+
if isinstance(self.connectivity, Graph):
375+
k = self.connectivity.cardinalities.max() + 1
376+
else:
377+
k = self.connectivity.max_neighbors + 1
372378
nn = self.connectivity.n - 1
373379
rids = np.array([np.random.permutation(nn)[0:k] for i in prange])
374380
ids = np.arange(self.connectivity.n)

esda/tests/test_moran_local_mv.py

Lines changed: 23 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -3,6 +3,7 @@
33
import geopandas
44
import pytest
55
from libpysal.weights import Queen
6+
from libpysal.graph import Graph
67
from sklearn.linear_model import TheilSenRegressor
78
from esda.moran_local_mv import Partial_Moran_Local, Auxiliary_Moran_Local
89
from esda.moran import Moran_Local_BV
@@ -12,21 +13,24 @@
1213

1314
y = df.HH_INC.values.reshape(-1,1)
1415
X = df.HSG_VAL.values.reshape(-1,1)
15-
w = Queen.from_dataframe(df)
16-
w.transform = 'r'
16+
w_classic = Queen.from_dataframe(df)
17+
w_classic.transform = 'r'
18+
g = Graph.build_contiguity(df).transform("r")
1719

18-
def test_partial_runs():
20+
pytestmark = pytest.mark.parametrize("w", [g, w_classic])
21+
22+
def test_partial_runs(w):
1923
"""Check if the class computes successfully in a default configuration"""
2024
m = Partial_Moran_Local(y,X,w, permutations=1)
2125
# done, just check if it runs
2226

23-
def test_partial_accuracy():
27+
def test_partial_accuracy(w):
2428
"""Check if the class outputs expected results at a given seed"""
2529
numpy.random.seed(111221)
2630
m = Partial_Moran_Local(y,X,w, permutations=10)
2731
# compute result by hand
2832
zy = (y - y.mean())/y.std()
29-
wz = (w.sparse * zy)
33+
wz = (w.sparse @ zy)
3034
zx = (X - X.mean(axis=0))/X.std(axis=0)
3135
rho = numpy.corrcoef(zy.squeeze(), zx.squeeze())[0,1]
3236
left = zy - rho * zx
@@ -45,7 +49,7 @@ def test_partial_accuracy():
4549
is_odd_label = m.labels_ % 2
4650
numpy.testing.assert_equal(is_cluster, is_odd_label)
4751

48-
def test_partial_unscaled():
52+
def test_partial_unscaled(w):
4953
"""Check if the variance scaling behaves as expected"""
5054
m = Partial_Moran_Local(y,X,w, permutations=0, unit_scale=True)
5155
m2 = Partial_Moran_Local(y,X,w, permutations=0, unit_scale=False)
@@ -55,33 +59,35 @@ def test_partial_unscaled():
5559
assert s1y > s2y, "variance is incorrectly scaled for y"
5660
assert s1x < s2x, "variance is incorrectly scaled for x"
5761

58-
def test_partial_uvquads():
62+
def test_partial_uvquads(w):
5963
"""Check that the quadrant decisions vary correctly with the inputs"""
6064
m = Partial_Moran_Local(y,X,w, permutations=0, mvquads=False)
6165
bv = Moran_Local_BV(y,X,w,permutations=0)
6266
# TODO: this currently fails, and it should pass. I am probably mis-calculating the bivariate quadrants for this option, and need to correct the code.
6367
numpy.testing.assert_array_equal(m.quads_, bv.q)
6468

65-
def test_aux_runs():
69+
def test_aux_runs(w):
70+
print(type(w), w.transform)
6671
"""Check that the class completes successfully in a default configuration"""
67-
a = Auxiliary_Moran_Local(y,X,w, permutations=1)
72+
a = Auxiliary_Moran_Local(y,X, w, permutations=1)
6873
#done, just check if it runs
6974

70-
def test_aux_accuracy():
75+
def test_aux_accuracy(w):
76+
print(type(w), w.transform)
7177
"""Check that the class outputs expected values for a given seed"""
7278
numpy.random.seed(111221)
7379
a = Auxiliary_Moran_Local(y,X,w, permutations=10)
7480

7581
# compute result by hand
7682
zy = (y - y.mean())/y.std()
77-
wz = (w.sparse * zy)
83+
wz = (w.sparse @ zy)
7884
zx = (X - X.mean(axis=0))/X.std(axis=0)
79-
wzx = w.sparse * zx
85+
wzx = w.sparse @ zx
8086
rho = numpy.corrcoef(zy.squeeze(), zx.squeeze())[0,1]
8187
mean = zy * wz - rho * zx * wz - rho * zy * wzx + rho**2 * zx * wzx
8288
scale = (w.n-1) / (w.n * (1 - rho**2))
8389

84-
manual = (mean * scale).squeeze()
90+
manual = numpy.asarray(mean * scale).squeeze()
8591
# check values, may not be identical because of the
8692
# matrix inversion least squares estimator used in scikit
8793
numpy.testing.assert_allclose(manual, a.associations_)
@@ -94,15 +100,17 @@ def test_aux_accuracy():
94100
is_odd_label = (a.labels_ % 2).astype(bool)
95101
numpy.testing.assert_equal(is_cluster, is_odd_label)
96102

97-
def test_aux_unscaled():
103+
def test_aux_unscaled(w):
104+
print(type(w), w.transform)
98105
"""Check that the variance scaling behaves as expected"""
99106
a = Auxiliary_Moran_Local(y,X/10000,w, permutations=0, unit_scale=True)
100107
a2 = Auxiliary_Moran_Local(y,X,w, permutations=0, unit_scale=False)
101108
assert (a.partials_.std(axis=0) < a2.partials_.std(axis=0)).all(), (
102109
"variance is not scaled correctly in partial regression."
103110
)
104111

105-
def test_aux_transformer():
112+
def test_aux_transformer(w):
113+
print(type(w), w.transform)
106114
"""Check that an alternative regressor can be used to calculate y|X"""
107115
a = Auxiliary_Moran_Local(y,X,w, permutations=0, transformer=TheilSenRegressor)
108116
# done, should just complete

0 commit comments

Comments
 (0)