stat/spatial: use mat.RowNonZeroDoer

stat/spatial/spatial.go

@@ -12,7 +12,6 @@ import (
 )
 
 // TODO(kortschak): Implement weighted routines.
-// TODO(kortschak): Make use of banded matrices when they exist in mat.
 
 // GetisOrdGStar returns the Local Getis-Ord G*i statistic for element of of the
 // weighted data using the provided locality matrix. The returned value is a z-score.
@@ -43,12 +42,20 @@ func GetisOrdGStar(i int, data, weights []float64, locality mat.Matrix) float64
 	n := float64(len(data))
 	mean, std := stat.MeanStdDev(data, weights)
 	var dwd, dww, sw float64
+	if doer, ok := locality.(mat.RowNonZeroDoer); ok {
+		doer.DoRowNonZero(i, func(_, j int, w float64) {
+			sw += w
+			dwd += w * data[j]
+			dww += w * w
+		})
+	} else {
 		for j, v := range data {
 			w := locality.At(i, j)
 			sw += w
 			dwd += w * v
 			dww += w * w
 		}
+	}
 	s := std * math.Sqrt((n-1)/n)
 
 	return (dwd - mean*sw) / (s * math.Sqrt((n*dww-sw*sw)/(n-1)))
@@ -73,11 +80,20 @@ func GlobalMoransI(data, weights []float64, locality mat.Matrix) (i, v, z float6
 	}
 	mean := stat.Mean(data, nil)
 
+	doer, isDoer := locality.(mat.RowNonZeroDoer)
+
 	// Calculate Moran's I for the data.
 	var num, den, sum float64
 	for i, xi := range data {
 		zi := xi - mean
 		den += zi * zi
+		if isDoer {
+			doer.DoRowNonZero(i, func(_, j int, w float64) {
+				sum += w
+				zj := data[j] - mean
+				num += w * zi * zj
+			})
+		} else {
 			for j, xj := range data {
 				w := locality.At(i, j)
 				sum += w
@@ -85,6 +101,7 @@ func GlobalMoransI(data, weights []float64, locality mat.Matrix) (i, v, z float6
 				num += w * zi * zj
 			}
 		}
+	}
 	i = (float64(len(data)) / sum) * (num / den)
 
 	// Calculate Moran's E(I) for the data.
@@ -102,6 +119,18 @@ func GlobalMoransI(data, weights []float64, locality mat.Matrix) (i, v, z float6
 		var4 += v * v
 
 		var p2 float64
+		if isDoer {
+			doer.DoRowNonZero(i, func(i, j int, wij float64) {
+				wji := locality.At(j, i)
+
+				s0 += wij
+
+				v := wij + wji
+				s1 += v * v
+
+				p2 += v
+			})
+		} else {
 			for j := range data {
 				wij := locality.At(i, j)
 				wji := locality.At(j, i)
@@ -113,6 +142,7 @@ func GlobalMoransI(data, weights []float64, locality mat.Matrix) (i, v, z float6
 
 				p2 += v
 			}
+		}
 		s2 += p2 * p2
 	}
 	s1 *= 0.5

stat/spatial/spatial_example_test.go

@@ -38,6 +38,35 @@ func ExampleGlobalMoransI_linear() {
 	// Moran's I=0.1111 z-score=0.6335
 }
 
+func ExampleGlobalMoransI_banded() {
+	data := []float64{0, 0, 0, 1, 1, 1, 0, 1, 0, 0}
+
+	// The locality here describes spatial neighbor
+	// relationships.
+	// This example uses the band matrix representation
+	// to improve time and space efficiency.
+	locality := mat.NewBandDense(10, 10, 1, 1, []float64{
+		0, 0, 1,
+		1, 0, 1,
+		1, 0, 1,
+		1, 0, 1,
+		1, 0, 1,
+		1, 0, 1,
+		1, 0, 1,
+		1, 0, 1,
+		1, 0, 1,
+		1, 0, 0,
+	})
+
+	i, _, z := spatial.GlobalMoransI(data, nil, locality)
+
+	fmt.Printf("Moran's I=%.4v z-score=%.4v\n", i, z)
+
+	// Output:
+	//
+	// Moran's I=0.1111 z-score=0.6335
+}
+
 func ExampleGetisOrdGStar() {
 	data := []float64{0, 0, 0, 1, 1, 1, 0, 1, 0, 0}
 
@@ -73,3 +102,41 @@ func ExampleGetisOrdGStar() {
 	// v=0 G*i=-0.2673
 	// v=0 G*i=-1.225
 }
+
+func ExampleGetisOrd_band() {
+	data := []float64{0, 0, 0, 1, 1, 1, 0, 1, 0, 0}
+
+	// The locality here describes spatial neighbor
+	// relationships including self.
+	// This example uses the band matrix representation
+	// to improve time and space efficiency.
+	locality := mat.NewBandDense(10, 10, 1, 1, []float64{
+		0, 1, 1,
+		1, 1, 1,
+		1, 1, 1,
+		1, 1, 1,
+		1, 1, 1,
+		1, 1, 1,
+		1, 1, 1,
+		1, 1, 1,
+		1, 1, 1,
+		1, 1, 0,
+	})
+
+	for i, v := range data {
+		fmt.Printf("v=%v G*i=% .4v\n", v, spatial.GetisOrdGStar(i, data, nil, locality))
+	}
+
+	// Output:
+	//
+	// v=0 G*i=-1.225
+	// v=0 G*i=-1.604
+	// v=0 G*i=-0.2673
+	// v=1 G*i= 1.069
+	// v=1 G*i= 2.405
+	// v=1 G*i= 1.069
+	// v=0 G*i= 1.069
+	// v=1 G*i=-0.2673
+	// v=0 G*i=-0.2673
+	// v=0 G*i=-1.225
+}

stat/spatial/spatial_test.go

@@ -13,7 +13,7 @@ import (
 	"gonum.org/v1/gonum/mat"
 )
 
-func simpleAdjacency(n, wide int, diag bool) *mat.Dense {
+func simpleAdjacency(n, wide int, diag bool) mat.Matrix {
 	m := mat.NewDense(n, n, nil)
 	for i := 0; i < n; i++ {
 		for j := 1; j <= wide; j++ {
@@ -30,11 +30,28 @@ func simpleAdjacency(n, wide int, diag bool) *mat.Dense {
 	return m
 }
 
+func simpleAdjacencyBand(n, wide int, diag bool) mat.Matrix {
+	m := mat.NewBandDense(n, n, wide, wide, nil)
+	for i := 0; i < n; i++ {
+		for j := 1; j <= wide; j++ {
+			if j > i {
+				continue
+			}
+			m.SetBand(i-j, i, 1)
+			m.SetBand(i, i-j, 1)
+		}
+		if diag {
+			m.SetBand(i, i, 1)
+		}
+	}
+	return m
+}
+
 var spatialTests = []struct {
 	from, to float64
 	n, wide  int
 	fn       func(float64, int, *rand.Rand) float64
-	locality func(n, wide int, diag bool) *mat.Dense
+	locality func(n, wide int, diag bool) mat.Matrix
 
 	// Values for MoranI and z-score are obtained from
 	// an R reference implementation.
@@ -46,6 +63,7 @@ var spatialTests = []struct {
 	// of the plotted data.
 	wantSegs int
 }{
+	// Dense matrix locality.
 	{
 		from: 0, to: 1, n: 1000, wide: 1,
 		fn: func(_ float64, _ int, rnd *rand.Rand) float64 {
@@ -122,6 +140,84 @@ var spatialTests = []struct {
 		wantZ:      -31.559531064275987,
 		wantSegs:   0,
 	},
+
+	// Band matrix locality.
+	{
+		from: 0, to: 1, n: 1000, wide: 1,
+		fn: func(_ float64, _ int, rnd *rand.Rand) float64 {
+			return rnd.Float64()
+		},
+		locality: simpleAdjacencyBand,
+
+		wantMoranI: -0.0019631298955953233,
+		wantZ:      -0.03039477405151108,
+		wantSegs:   0,
+	},
+	{
+		from: -math.Pi / 2, to: 3 * math.Pi / 2, n: 1000, wide: 1,
+		fn: func(x float64, _ int, _ *rand.Rand) float64 {
+			y := math.Sin(x)
+			if math.Abs(y) > 0.5 {
+				y *= 1/math.Abs(y) - 1
+			}
+			return y * math.Sin(x*2)
+		},
+		locality: simpleAdjacencyBand,
+
+		wantMoranI: 1.0008149537991464,
+		wantZ:      31.648547078779092,
+		wantSegs:   4,
+	},
+	{
+		from: 0, to: 1, n: 1000, wide: 1,
+		fn: func(_ float64, _ int, rnd *rand.Rand) float64 {
+			return rnd.NormFloat64()
+		},
+		locality: simpleAdjacencyBand,
+
+		wantMoranI: 0.031195199553564902,
+		wantZ:      1.0171161514080056,
+		wantSegs:   0,
+	},
+	{
+		from: 0, to: 1, n: 1000, wide: 1,
+		fn: func(x float64, _ int, rnd *rand.Rand) float64 {
+			if rnd.Float64() < 0.5 {
+				return rnd.NormFloat64() + 5
+			}
+			return rnd.NormFloat64()
+		},
+		locality: simpleAdjacencyBand,
+
+		wantMoranI: -0.016245135637562223,
+		wantZ:      -0.48157993864993476,
+		wantSegs:   0,
+	},
+	{
+		from: 0, to: 1, n: 1000, wide: 1,
+		fn: func(x float64, i int, rnd *rand.Rand) float64 {
+			if i%2 == 0 {
+				return rnd.NormFloat64() + 5
+			}
+			return rnd.NormFloat64()
+		},
+		locality: simpleAdjacencyBand,
+
+		wantMoranI: -0.8565268969272998,
+		wantZ:      -27.027057520918113,
+		wantSegs:   0,
+	},
+	{
+		from: 0, to: 1, n: 1000, wide: 1,
+		fn: func(_ float64, i int, _ *rand.Rand) float64 {
+			return float64(i % 2)
+		},
+		locality: simpleAdjacencyBand,
+
+		wantMoranI: -1,
+		wantZ:      -31.559531064275987,
+		wantSegs:   0,
+	},
 }
 
 func TestGetisOrd(t *testing.T) {