add GeoGcd

* great circle distance utility

6 files changed, 155 insertions, 0 deletions

diff --git a/searchlib/CMakeLists.txt b/searchlib/CMakeLists.txt
index c9c44af3aa3..de46dfa6441 100644
--- a/searchlib/CMakeLists.txt
+++ b/searchlib/CMakeLists.txt
@@ -104,6 +104,7 @@ vespa_define_module(
     src/tests/bitvector
     src/tests/btree
     src/tests/common/bitvector
+    src/tests/common/geogcd
     src/tests/common/location
     src/tests/common/location_iterator
     src/tests/common/matching_elements
diff --git a/searchlib/src/tests/common/geogcd/CMakeLists.txt b/searchlib/src/tests/common/geogcd/CMakeLists.txt
new file mode 100644
index 00000000000..3721b609668
--- /dev/null
+++ b/searchlib/src/tests/common/geogcd/CMakeLists.txt
@@ -0,0 +1,9 @@
+# Copyright Yahoo. Licensed under the terms of the Apache 2.0 license. See LICENSE in the project root.
+vespa_add_executable(searchlib_geo_gcd_test_app TEST
+    SOURCES
+    geo_gcd_test.cpp
+    DEPENDS
+    searchlib
+    GTest::GTest
+)
+vespa_add_test(NAME searchlib_geo_gcd_test_app COMMAND searchlib_geo_gcd_test_app)
diff --git a/searchlib/src/tests/common/geogcd/geo_gcd_test.cpp b/searchlib/src/tests/common/geogcd/geo_gcd_test.cpp
new file mode 100644
index 00000000000..dd218c4ff04
--- /dev/null
+++ b/searchlib/src/tests/common/geogcd/geo_gcd_test.cpp
@@ -0,0 +1,63 @@
+// Copyright Yahoo. Licensed under the terms of the Apache 2.0 license. See LICENSE in the project root.
+
+#include <stdio.h>
+#include <vespa/searchlib/common/geo_gcd.h>
+#include <vespa/vespalib/gtest/gtest.h>
+
+using namespace search::common;
+
+struct Point {
+    const char *name;
+    double lat;
+    double lng;
+};
+
+constexpr size_t NUM = 9;
+
+const Point airports[NUM] = {
+    { "SFO", 37.61, -122.38},
+    { "LHR", 51.47, -0.46},
+    { "OSL", 60.20, 11.08},
+    { "GIG", -22.8, -43.25},
+    { "HKG", 22.31, 113.91},
+    { "TRD", 63.45, 10.92},
+    { "SYD", -33.95, 151.17},
+    { "LAX", 33.94, -118.41},
+    { "JFK", 40.64, -73.78},
+};
+
+const double exact_distances[NUM][NUM] = {
+    { 0, 5367, 5196, 6604, 6927, 5012, 7417, 337, 2586 },
+    { 5367, 0, 750, 5734, 5994, 928, 10573, 5456, 3451 },
+    { 5196, 750, 0, 6479, 5319, 226, 9888, 5345, 3687 },
+    { 6604, 5734, 6479, 0, 10989, 6623, 8414, 6294, 4786 },
+    { 6927, 5994, 5319, 10989, 0, 5240, 4581, 7260, 8072 },
+    { 5012, 928, 226, 6623, 5240, 0, 9782, 5171, 3611 },
+    { 7417, 10573, 9888, 8414, 4581, 9782, 0, 7488, 9950 },
+    { 337, 5456, 5345, 6294, 7260, 5171, 7488, 0, 2475 },
+    { 2586, 3451, 3687, 4786, 8072, 3611, 9950, 2475, 0 }
+};
+
+TEST(GeoGcdTest, computed_distances_seem_legit) {
+    for (size_t i = 0; i < NUM; ++i) {
+        const Point &from = airports[i];
+        printf("\n");
+        GeoGcd geo_from{from.lat, from.lng};
+        for (size_t j = 0; j < NUM; ++j) {
+            const Point &to = airports[j];
+            double km = geo_from.km_great_circle_distance(to.lat, to.lng);
+            double miles = km / 1.609344;
+            EXPECT_GE(miles, 0);
+            if (from.name == to.name) {
+                EXPECT_EQ(miles, 0.0);
+            } else {
+                double exact = exact_distances[i][j];
+                printf("Distance from %s to %s (in miles): %.1f [more exact would be %.1f]\n", from.name, to.name, miles, exact);
+                EXPECT_LT(miles*0.99, exact);
+                EXPECT_GT(miles*1.01, exact);
+            }
+        }
+    }
+}
+
+GTEST_MAIN_RUN_ALL_TESTS()
diff --git a/searchlib/src/vespa/searchlib/common/CMakeLists.txt b/searchlib/src/vespa/searchlib/common/CMakeLists.txt
index 73c8999520b..a7c8d56f11d 100644
--- a/searchlib/src/vespa/searchlib/common/CMakeLists.txt
+++ b/searchlib/src/vespa/searchlib/common/CMakeLists.txt
@@ -12,6 +12,7 @@ vespa_add_library(searchlib_common OBJECT
     featureset.cpp
     fileheadercontext.cpp
     flush_token.cpp
+    geo_gcd.cpp
     geo_location.cpp
     geo_location_parser.cpp
     geo_location_spec.cpp
diff --git a/searchlib/src/vespa/searchlib/common/geo_gcd.cpp b/searchlib/src/vespa/searchlib/common/geo_gcd.cpp
new file mode 100644
index 00000000000..a7fc870fc31
--- /dev/null
+++ b/searchlib/src/vespa/searchlib/common/geo_gcd.cpp
@@ -0,0 +1,50 @@
+// Copyright Yahoo. Licensed under the terms of the Apache 2.0 license. See LICENSE in the project root.
+
+#include "geo_gcd.h"
+
+namespace search::common {
+
+namespace {
+
+// in km, as defined by IUGG, see:
+// https://en.wikipedia.org/wiki/Earth_radius#Mean_radius
+static constexpr double earth_mean_radius = 6371.0088;
+
+static constexpr double degrees_to_radians = M_PI / 180.0;
+
+// with input in radians
+double greatCircleDistance(double theta_A, double phi_A,
+                           double theta_B, double phi_B)
+{
+    // convert to radians:
+    double theta_diff = theta_A - theta_B;
+    double phi_diff = phi_A - phi_B;
+    // haversines of differences:
+    double hav_theta = GeoGcd::haversine(theta_diff);
+    double hav_phi   = GeoGcd::haversine(phi_diff);
+    // haversine of central angle between the two points:
+    double hav_central_angle = hav_theta + cos(theta_A)*cos(theta_B)*hav_phi;
+    // sine of half the central angle:
+    double half_sine_diff = sqrt(hav_central_angle);
+    // distance in kilometers:
+    double d = 2 * asin(half_sine_diff) * earth_mean_radius;
+    return d;
+}
+
+}
+
+GeoGcd::GeoGcd(double lat, double lng)
+    : _latitude_radians(lat * degrees_to_radians),
+      _longitude_radians(lng * degrees_to_radians)
+{}
+
+
+double GeoGcd::km_great_circle_distance(double lat, double lng) const {
+    double theta_A = _latitude_radians;
+    double phi_A   = _longitude_radians;
+    double theta_B = lat * degrees_to_radians;
+    double phi_B   = lng * degrees_to_radians;
+    return  greatCircleDistance(theta_A, phi_A, theta_B, phi_B);
+}
+
+} // namespace search::common
diff --git a/searchlib/src/vespa/searchlib/common/geo_gcd.h b/searchlib/src/vespa/searchlib/common/geo_gcd.h
new file mode 100644
index 00000000000..acd10207057
--- /dev/null
+++ b/searchlib/src/vespa/searchlib/common/geo_gcd.h
@@ -0,0 +1,31 @@
+// Copyright Yahoo. Licensed under the terms of the Apache 2.0 license. See LICENSE in the project root.
+
+#pragma once
+
+#include <cstdint>
+#include <cmath>
+
+namespace search::common {
+
+/**
+ * An immutable struct for a (geo) location point,
+ * with methods for computing Great Circle Distance
+ * using the haversine formula
+ **/
+struct GeoGcd
+{
+    GeoGcd(double lat, double lng);
+
+    // haversine function:
+    static constexpr double haversine(double angle) {
+        double s = sin(0.5*angle);
+        return s*s;
+    }
+
+    double km_great_circle_distance(double lat, double lng) const;
+private:
+    double _latitude_radians;
+    double _longitude_radians;
+};
+
+} // namespace