util/neighbour: add basic nearest-neighbour package

Change-Id: I1aeb02f5de6d9ac47e5127cb1c0d5c2f4ea6457a

util/neighbour/nearest.go

@@ -0,0 +1,71 @@
+// Package neighbour contains a basic non-optimized nearest-neighbour
+// algorithm implementation for terrestrial GPS coordinates.
+package neighbour
+
+import (
+	"math"
+
+	"golang.org/x/exp/slices"
+)
+
+// Location is a latitude/longitude pair representing a point on Earth.
+type Location struct {
+	Latitude  float64
+	Longitude float64
+}
+
+// Distance calculates the great-circle distance between two points using the
+// Haversine formula, in kilometers.
+//
+// This is also known as the "as-the-crow-flies" distance.
+func (l Location) Distance(other Location) float64 {
+	// For the following variable definitions:
+	//    φ is latitude ("phi")
+	//    λ is longitude ("lambda")
+	//    R is earth’s radius (mean radius = 6,371km)
+	//
+	// We can calculate the distance using the haversine formula as such:
+	//    a = sin²((φB - φA)/2) + cos φA * cos φB * sin²((λB - λA)/2)
+	//    c = 2 * atan2( √a, √(1−a) )
+	//    d = R * c
+
+	// Convert our latitude/longitude to radians, since the various math
+	// functions take radians but latitude/longitude are in degrees.
+	lat1, lon1 := degreesToRadians(l.Latitude), degreesToRadians(l.Longitude)
+	lat2, lon2 := degreesToRadians(other.Latitude), degreesToRadians(other.Longitude)
+
+	deltaPhi := lat2 - lat1
+	deltaLambda := lon2 - lon1
+
+	// Haversine
+	a := math.Pow(math.Sin(deltaPhi/2), 2) +
+		math.Cos(lat1)*math.Cos(lat2)*math.Pow(math.Sin(deltaLambda/2), 2)
+	c := 2 * math.Atan2(math.Sqrt(a), math.Sqrt(1-a))
+
+	// Return the distance in km.
+	const earthRadiusKm = 6371
+	return c * earthRadiusKm
+}
+
+func degreesToRadians(d float64) float64 {
+	return d * math.Pi / 180
+}
+
+// Neighbours returns the nearest n neighbours to the provided point from the
+// set of candidates.
+func Neighbours(point Location, n int, candidates []Location) []Location {
+	// Calculate all distances up-front to avoid recalculating during the
+	// sort below.
+	distances := map[Location]float64{}
+	for _, candidate := range candidates {
+		distances[candidate] = candidate.Distance(point)
+	}
+
+	// Sort the candidates slice by their distance to the provided point.
+	candidates = slices.Clone(candidates)
+	slices.SortFunc(candidates, func(a, b Location) bool {
+		return distances[a] < distances[b]
+	})
+
+	return candidates[:n]
+}

util/neighbour/nearest_test.go

@@ -0,0 +1,51 @@
+package neighbour
+
+import (
+	"math"
+	"reflect"
+	"testing"
+)
+
+func TestHaversine(t *testing.T) {
+	one := Location{51.510357, -0.116773}  // King's College, London
+	two := Location{38.889931, -77.009003} // The White House
+
+	dist := one.Distance(two)
+	want := 5897.658
+
+	if math.Abs(want-dist) > 0.001 {
+		t.Fatalf("distance mismatch; got %v, want %v", dist, want)
+	}
+}
+
+func TestNeighbours(t *testing.T) {
+	// Provincial capitals
+	capitals := []Location{
+		{48.4283182, -123.3649533}, // Victoria, BC, Canada
+		{60.721571, -135.054932},   // Whitehorse, YT, Canada
+		{53.5462055, -113.491241},  // Edmonton, AB, Canada
+		{62.4540807, -114.377385},  // Yellowknife, NT, Canada
+		{50.44876, -104.61731},     // Regina, SK, Canada
+		{49.8955367, -97.1384584},  // Winnipeg, MB, Canada
+		{63.74944, -68.521857},     // Iqaluit, NU, Canada
+		{43.6534817, -79.3839347},  // Toronto, ON, Canada
+		{45.5031824, -73.5698065},  // Montreal, QC, Canada
+		{45.94780155, -66.6534707}, // Fredericton, NB, Canada
+		{44.648618, -63.5859487},   // Halifax, NS, Canada
+		{46.234953, -63.132935},    // Charlottetown, PE, Canada
+		{47.5614705, -52.7126162},  // St. John’s, NL, Canada
+	}
+
+	// Thunder Bay, Ontario, Canada
+	point := Location{48.382221, -89.246109}
+	nearest := Neighbours(point, 4, capitals)
+	want := []Location{
+		{49.8955367, -97.1384584}, // Winnipeg, MB, Canada
+		{43.6534817, -79.3839347}, // Toronto, ON, Canada
+		{50.44876, -104.61731},    // Regina, SK, Canada
+		{45.5031824, -73.5698065}, // Montreal, QC, Canada
+	}
+	if !reflect.DeepEqual(nearest, want) {
+		t.Errorf("nearest points mismatch\ngot: %v\nwant: %v", nearest, want)
+	}
+}