From df1b842fcd484315d4ca441a36634d82a42b14b3 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Wed, 26 Nov 2025 23:46:34 -0500 Subject: [PATCH 01/62] area.c --- justfile | 45 +++++++++++++++++++++++++++ src/h3lib/lib/area.c | 74 ++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 119 insertions(+) create mode 100644 justfile create mode 100644 src/h3lib/lib/area.c diff --git a/justfile b/justfile new file mode 100644 index 0000000000..8ab59e9ced --- /dev/null +++ b/justfile @@ -0,0 +1,45 @@ +# TODO: Remove before landing + +init: purge + mkdir build + +build: + cd build; cmake -DCMAKE_BUILD_TYPE=Release ..; make + +profile: init + cd build; cmake -DCMAKE_BUILD_TYPE=RelWithDebInfo -DCMAKE_C_FLAGS='-fno-omit-frame-pointer' ..; make + # cd build; cmake -DCMAKE_BUILD_TYPE=Release ..; make + xcrun xctrace record \ + --template 'Time Profiler' \ + --output h3-prof.trace \ + --launch -- ./build/bin/benchmarkCellsToLinkedMultiPolygon + open h3-prof.trace + +purge: + rm -rf build + rm -rf *.trace + rm -rf .ipynb_checkpoints + +test: build + # ./build/bin/testH3CellAreaExhaustive + # ./build/bin/testEdgeCellsToPoly + ./build/bin/testDirectedEdge + # ./build/bin/testArea + # just test-slow + +time: + time ./build/bin/testArea + +test-fast: build + cd build; make test-fast + +test-slow: build + cd build; make test + +bench: build + ./build/bin/benchmarkCellsToLinkedMultiPolygon + # ./build/bin/benchmarkCellsToMultiPolygon + # ./build/bin/benchmarkDirectedEdge + +fuzz: + cd build; CC=clang cmake -DENABLE_LIBFUZZER=ON ..; make fuzzers diff --git a/src/h3lib/lib/area.c b/src/h3lib/lib/area.c new file mode 100644 index 0000000000..e14677f307 --- /dev/null +++ b/src/h3lib/lib/area.c @@ -0,0 +1,74 @@ +#include +#include +#include +#include +#include +#include +#include + +#include "adder.h" +#include "alloc.h" +#include "h3api.h" + + +void H3_EXPORT(destroyGeoLoop)(GeoLoop* loop) { + H3_MEMORY(free)(loop->verts); + loop->verts = NULL; +} + +static double cagnoli(LatLng x, LatLng y) { + x.lat = x.lat / 2.0 + M_PI / 4.0; + y.lat = y.lat / 2.0 + M_PI / 4.0; + + double sa = sin(x.lat) * sin(y.lat); + double ca = cos(x.lat) * cos(y.lat); + + double d = y.lng - x.lng; + double sd = sin(d); + double cd = cos(d); + + return -2.0 * atan2(sa * sd, sa * cd + ca); +} + +static double vertArea(LatLng *verts, int numVerts) { + Kahan k = {0.0, 0.0}; + + for (int i = 0; i < numVerts; i++) { + int j = (i + 1) % numVerts; + kadd(&k, cagnoli(verts[i], verts[j])); + } + + if (k.sum < 0.0) { + kadd(&k, 4.0 * M_PI); + } + return k.sum; +} + +H3Error H3_EXPORT(geoLoopArea)(GeoLoop loop, double *out) { + *out = vertArea(loop.verts, loop.numVerts); + return E_SUCCESS; +} + +/** + * Area of H3 cell in radians^2. + * + * The area is calculated by breaking the cell into spherical triangles and + * summing up their areas. Note that some H3 cells (hexagons and pentagons) + * are irregular, and have more than 6 or 5 sides. + * + * todo: optimize the computation by re-using the edges shared between triangles + * + * @param cell H3 cell + * @param out cell area in radians^2 + * @return E_SUCCESS on success, or an error code otherwise + */ +H3Error H3_EXPORT(cellAreaRads2)(H3Index cell, double *out) { + CellBoundary cb; + H3Error err = H3_EXPORT(cellToBoundary)(cell, &cb); + if (err) { + return err; + } + *out = vertArea(cb.verts, cb.numVerts); + + return E_SUCCESS; +} From b6d9ddb9c16ba6881f8d76f9c79329fb881407da Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Wed, 26 Nov 2025 23:58:13 -0500 Subject: [PATCH 02/62] adding in more stuff --- CMakeLists.txt | 5 ++- src/apps/testapps/testArea.c | 56 +++++++++++++++++++++++++++ src/h3lib/include/adder.h | 16 ++++++++ src/h3lib/include/h3api.h.in | 3 ++ src/h3lib/lib/latLng.c | 74 ------------------------------------ 5 files changed, 79 insertions(+), 75 deletions(-) create mode 100644 src/apps/testapps/testArea.c create mode 100644 src/h3lib/include/adder.h diff --git a/CMakeLists.txt b/CMakeLists.txt index 7eb6dbddb5..828c02bcef 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -170,6 +170,7 @@ set(LIB_SOURCE_FILES src/h3lib/include/constants.h src/h3lib/include/coordijk.h src/h3lib/include/algos.h + src/h3lib/include/adder.h src/h3lib/lib/h3Assert.c src/h3lib/lib/algos.c src/h3lib/lib/coordijk.c @@ -188,7 +189,8 @@ set(LIB_SOURCE_FILES src/h3lib/lib/iterators.c src/h3lib/lib/vertexGraph.c src/h3lib/lib/faceijk.c - src/h3lib/lib/baseCells.c) + src/h3lib/lib/baseCells.c + src/h3lib/lib/area.c) set(APP_SOURCE_FILES src/apps/applib/include/kml.h src/apps/applib/include/benchmark.h @@ -278,6 +280,7 @@ set(OTHER_SOURCE_FILES src/apps/testapps/testH3IteratorsInternal.c src/apps/testapps/testMathExtensionsInternal.c src/apps/testapps/testDescribeH3Error.c + src/apps/testapps/testArea.c src/apps/miscapps/cellToBoundaryHier.c src/apps/miscapps/cellToLatLngHier.c src/apps/miscapps/generateBaseCellNeighbors.c diff --git a/src/apps/testapps/testArea.c b/src/apps/testapps/testArea.c new file mode 100644 index 0000000000..1482e74176 --- /dev/null +++ b/src/apps/testapps/testArea.c @@ -0,0 +1,56 @@ +#include +#include +#include +#include + +#include "h3api.h" +#include "iterators.h" +#include "test.h" +#include "utility.h" +#include "adder.h" + + +static void do_res_sum(int res) { + Kahan k = {0,0}; + + for (IterCellsResolution iter = iterInitRes(res); iter.h; + iterStepRes(&iter)) { + double cellArea; + H3_EXPORT(cellAreaRads2)(iter.h, &cellArea); + + kadd(&k, cellArea); + } + + double area = k.sum; + double diff = fabs(area - 4 * M_PI); + printf("res: %d, diff: %e\n", res, diff); +} + +SUITE(test_for_area) { + TEST(some_area_test) { + // the numerical test is how close to 4*pi do we get when adding up all + // cells at smaller and smaller resolution? + + printf("\n"); + + do_res_sum(0); + do_res_sum(1); + do_res_sum(2); + do_res_sum(3); + do_res_sum(4); + // do_res_sum(5); + // do_res_sum(6); + // do_res_sum(7); + // do_res_sum(8); + + // res: 0, diff: 7.105427e-15 + // res: 1, diff: 3.907985e-14 + // res: 2, diff: 1.421085e-13 + // res: 3, diff: 4.689582e-13 + // res: 4, diff: 2.161826e-12 + // res: 5, diff: 1.543654e-12 + // res: 6, diff: 1.768363e-11 + // res: 7, diff: 1.490719e-11 + // res: 8, diff: 2.937917e-11 + } +} diff --git a/src/h3lib/include/adder.h b/src/h3lib/include/adder.h new file mode 100644 index 0000000000..60e19891d6 --- /dev/null +++ b/src/h3lib/include/adder.h @@ -0,0 +1,16 @@ +#ifndef ADDER_H +#define ADDER_H + +typedef struct { + double sum; // running total + double c; // compensation term +} Kahan; + +static inline void kadd(Kahan *k, double x) { + double y = x - k->c; + double t = k->sum + y; + k->c = (t - k->sum) - y; + k->sum = t; +} + +#endif diff --git a/src/h3lib/include/h3api.h.in b/src/h3lib/include/h3api.h.in index dfca733620..a3159cffee 100644 --- a/src/h3lib/include/h3api.h.in +++ b/src/h3lib/include/h3api.h.in @@ -348,6 +348,9 @@ DECLSPEC H3Error H3_EXPORT(cellsToLinkedMultiPolygon)(const H3Index *h3Set, DECLSPEC void H3_EXPORT(destroyLinkedMultiPolygon)(LinkedGeoPolygon *polygon); /** @} */ +DECLSPEC void H3_EXPORT(destroyGeoLoop)(GeoLoop* loop); +DECLSPEC H3Error H3_EXPORT(geoLoopArea)(GeoLoop loop, double *out) + /** @defgroup degsToRads degsToRads * Functions for degsToRads * @{ diff --git a/src/h3lib/lib/latLng.c b/src/h3lib/lib/latLng.c index f90d9ee2b4..13bf196e8e 100644 --- a/src/h3lib/lib/latLng.c +++ b/src/h3lib/lib/latLng.c @@ -353,80 +353,6 @@ H3Error H3_EXPORT(getNumCells)(int res, int64_t *out) { return E_SUCCESS; } -/** - * Surface area in radians^2 of spherical triangle on unit sphere. - * - * For the math, see: - * https://en.wikipedia.org/wiki/Spherical_trigonometry#Area_and_spherical_excess - * - * @param a length of triangle side A in radians - * @param b length of triangle side B in radians - * @param c length of triangle side C in radians - * - * @return area in radians^2 of triangle on unit sphere - */ -double triangleEdgeLengthsToArea(double a, double b, double c) { - double s = (a + b + c) * 0.5; - - a = (s - a) * 0.5; - b = (s - b) * 0.5; - c = (s - c) * 0.5; - s = s * 0.5; - - return 4 * atan(sqrt(tan(s) * tan(a) * tan(b) * tan(c))); -} - -/** - * Compute area in radians^2 of a spherical triangle, given its vertices. - * - * @param a vertex lat/lng in radians - * @param b vertex lat/lng in radians - * @param c vertex lat/lng in radians - * - * @return area of triangle on unit sphere, in radians^2 - */ -double triangleArea(const LatLng *a, const LatLng *b, const LatLng *c) { - return triangleEdgeLengthsToArea(H3_EXPORT(greatCircleDistanceRads)(a, b), - H3_EXPORT(greatCircleDistanceRads)(b, c), - H3_EXPORT(greatCircleDistanceRads)(c, a)); -} - -/** - * Area of H3 cell in radians^2. - * - * The area is calculated by breaking the cell into spherical triangles and - * summing up their areas. Note that some H3 cells (hexagons and pentagons) - * are irregular, and have more than 6 or 5 sides. - * - * todo: optimize the computation by re-using the edges shared between triangles - * - * @param cell H3 cell - * @param out cell area in radians^2 - * @return E_SUCCESS on success, or an error code otherwise - */ -H3Error H3_EXPORT(cellAreaRads2)(H3Index cell, double *out) { - LatLng c; - CellBoundary cb; - H3Error err = H3_EXPORT(cellToLatLng)(cell, &c); - if (err) { - return err; - } - err = H3_EXPORT(cellToBoundary)(cell, &cb); - if (NEVER(err)) { - // Uncoverable because cellToLatLng will have returned an error already - return err; - } - - double area = 0.0; - for (int i = 0; i < cb.numVerts; i++) { - int j = (i + 1) % cb.numVerts; - area += triangleArea(&cb.verts[i], &cb.verts[j], &c); - } - - *out = area; - return E_SUCCESS; -} - /** * Area of H3 cell in kilometers^2. * From 8080e568e290766f61c0f049fa762c332f08a705 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Thu, 27 Nov 2025 00:03:06 -0500 Subject: [PATCH 03/62] merp --- CMakeTests.cmake | 1 + justfile | 4 ++-- src/apps/testapps/testArea.c | 5 ++--- src/h3lib/include/h3api.h.in | 4 ++-- src/h3lib/lib/area.c | 3 +-- 5 files changed, 8 insertions(+), 9 deletions(-) diff --git a/CMakeTests.cmake b/CMakeTests.cmake index 1901a61fac..cb705fc94a 100644 --- a/CMakeTests.cmake +++ b/CMakeTests.cmake @@ -259,6 +259,7 @@ add_h3_test(testH3IteratorsInternal src/apps/testapps/testH3IteratorsInternal.c) add_h3_test(testMathExtensionsInternal src/apps/testapps/testMathExtensionsInternal.c) add_h3_test(testDescribeH3Error src/apps/testapps/testDescribeH3Error.c) +add_h3_test(testArea src/apps/testapps/testArea.c) add_h3_test_with_arg(testH3NeighborRotations src/apps/testapps/testH3NeighborRotations.c 0) diff --git a/justfile b/justfile index 8ab59e9ced..c9992e3ce8 100644 --- a/justfile +++ b/justfile @@ -23,8 +23,8 @@ purge: test: build # ./build/bin/testH3CellAreaExhaustive # ./build/bin/testEdgeCellsToPoly - ./build/bin/testDirectedEdge - # ./build/bin/testArea + # ./build/bin/testDirectedEdge + ./build/bin/testArea # just test-slow time: diff --git a/src/apps/testapps/testArea.c b/src/apps/testapps/testArea.c index 1482e74176..b103c950ec 100644 --- a/src/apps/testapps/testArea.c +++ b/src/apps/testapps/testArea.c @@ -3,15 +3,14 @@ #include #include +#include "adder.h" #include "h3api.h" #include "iterators.h" #include "test.h" #include "utility.h" -#include "adder.h" - static void do_res_sum(int res) { - Kahan k = {0,0}; + Kahan k = {0, 0}; for (IterCellsResolution iter = iterInitRes(res); iter.h; iterStepRes(&iter)) { diff --git a/src/h3lib/include/h3api.h.in b/src/h3lib/include/h3api.h.in index a3159cffee..042b778471 100644 --- a/src/h3lib/include/h3api.h.in +++ b/src/h3lib/include/h3api.h.in @@ -348,8 +348,8 @@ DECLSPEC H3Error H3_EXPORT(cellsToLinkedMultiPolygon)(const H3Index *h3Set, DECLSPEC void H3_EXPORT(destroyLinkedMultiPolygon)(LinkedGeoPolygon *polygon); /** @} */ -DECLSPEC void H3_EXPORT(destroyGeoLoop)(GeoLoop* loop); -DECLSPEC H3Error H3_EXPORT(geoLoopArea)(GeoLoop loop, double *out) +DECLSPEC void H3_EXPORT(destroyGeoLoop)(GeoLoop *loop); +DECLSPEC H3Error H3_EXPORT(geoLoopArea)(GeoLoop loop, double *out); /** @defgroup degsToRads degsToRads * Functions for degsToRads diff --git a/src/h3lib/lib/area.c b/src/h3lib/lib/area.c index e14677f307..08a227629c 100644 --- a/src/h3lib/lib/area.c +++ b/src/h3lib/lib/area.c @@ -10,8 +10,7 @@ #include "alloc.h" #include "h3api.h" - -void H3_EXPORT(destroyGeoLoop)(GeoLoop* loop) { +void H3_EXPORT(destroyGeoLoop)(GeoLoop *loop) { H3_MEMORY(free)(loop->verts); loop->verts = NULL; } From 5c8950ade3992c4c55a56e3f00261645917c1c33 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Thu, 27 Nov 2025 00:06:05 -0500 Subject: [PATCH 04/62] kahan accumulator --- src/apps/testapps/testH3CellAreaExhaustive.c | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/src/apps/testapps/testH3CellAreaExhaustive.c b/src/apps/testapps/testH3CellAreaExhaustive.c index 13261dae15..c0bb793de7 100644 --- a/src/apps/testapps/testH3CellAreaExhaustive.c +++ b/src/apps/testapps/testH3CellAreaExhaustive.c @@ -136,15 +136,16 @@ static void cell_area_assert(H3Index cell) { */ static void earth_area_test(int res, H3Error (*cell_area)(H3Index, double *), double target, double tol) { - double area = 0.0; - for (IterCellsResolution iter = iterInitRes(res); iter.h; - iterStepRes(&iter)) { + Kahan k = {0.0, 0.0}; + IterCellsResolution iter = iterInitRes(res); + + for (; iter.h; iterStepRes(&iter)) { double cellArea; t_assertSuccess((*cell_area)(iter.h, &cellArea)); - area += cellArea; + kadd(&k, cellArea); } - t_assert(fabs(area - target) < tol, + t_assert(fabs(k.sum - target) < tol, "sum of all cells should give earth area"); } From b82cab22a46877986dd70c6d9ab58e967ac54ed8 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Thu, 27 Nov 2025 11:18:35 -0500 Subject: [PATCH 05/62] neumaier_add --- src/apps/testapps/testArea.c | 15 ++++++--- src/apps/testapps/testH3CellAreaExhaustive.c | 1 + src/h3lib/include/adder.h | 32 +++++++++++++++++++- src/h3lib/lib/area.c | 2 +- 4 files changed, 44 insertions(+), 6 deletions(-) diff --git a/src/apps/testapps/testArea.c b/src/apps/testapps/testArea.c index b103c950ec..813b542d14 100644 --- a/src/apps/testapps/testArea.c +++ b/src/apps/testapps/testArea.c @@ -9,6 +9,13 @@ #include "test.h" #include "utility.h" +// TODO: turn this into a benchmark, report accuracy and time +// see if we get any benefit by avoiding the repeated sin/cos computation + +// TODO: document algos: https://en.wikipedia.org/wiki/Kahan_summation_algorithm#Further_enhancements +// TODO: better more general name for accumulator struct and functions +// TOOD: demonstrate the area alg works for global polygons. + static void do_res_sum(int res) { Kahan k = {0, 0}; @@ -20,7 +27,7 @@ static void do_res_sum(int res) { kadd(&k, cellArea); } - double area = k.sum; + double area = kresult(k); double diff = fabs(area - 4 * M_PI); printf("res: %d, diff: %e\n", res, diff); } @@ -37,9 +44,9 @@ SUITE(test_for_area) { do_res_sum(2); do_res_sum(3); do_res_sum(4); - // do_res_sum(5); - // do_res_sum(6); - // do_res_sum(7); + do_res_sum(5); + do_res_sum(6); + do_res_sum(7); // do_res_sum(8); // res: 0, diff: 7.105427e-15 diff --git a/src/apps/testapps/testH3CellAreaExhaustive.c b/src/apps/testapps/testH3CellAreaExhaustive.c index c0bb793de7..3f39fac7bb 100644 --- a/src/apps/testapps/testH3CellAreaExhaustive.c +++ b/src/apps/testapps/testH3CellAreaExhaustive.c @@ -23,6 +23,7 @@ #include +#include "adder.h" #include "constants.h" #include "h3api.h" #include "iterators.h" diff --git a/src/h3lib/include/adder.h b/src/h3lib/include/adder.h index 60e19891d6..a9763be3ed 100644 --- a/src/h3lib/include/adder.h +++ b/src/h3lib/include/adder.h @@ -6,11 +6,41 @@ typedef struct { double c; // compensation term } Kahan; -static inline void kadd(Kahan *k, double x) { +static inline void simple_add(Kahan *k, double x) { k->sum += x; } + +static inline void kahan_add(Kahan *k, double x) { double y = x - k->c; double t = k->sum + y; k->c = (t - k->sum) - y; k->sum = t; } +static inline void neumaier_add(Kahan *k, double x) { + double t = k->sum + x; + + if (fabs(k->sum) >= fabs(x)) { + k->c += (k->sum - t) + x; + } else { + k->c += (x - t) + k->sum; + } + + k->sum = t; +} + +static inline double simple_result(Kahan k) { return k.sum; } +static inline double kahan_result(Kahan k) { return k.sum; } +static inline double neumaier_result(Kahan k) { return k.sum + k.c; } + +static inline void kadd(Kahan *k, double x) { + // simple_add(k, x); + // kahan_add(k, x); + neumaier_add(k, x); +} + +static inline double kresult(Kahan k) { + // return simple_result(k); + // return kahan_result(k); + return neumaier_result(k); +} + #endif diff --git a/src/h3lib/lib/area.c b/src/h3lib/lib/area.c index 08a227629c..1eef4940e9 100644 --- a/src/h3lib/lib/area.c +++ b/src/h3lib/lib/area.c @@ -40,7 +40,7 @@ static double vertArea(LatLng *verts, int numVerts) { if (k.sum < 0.0) { kadd(&k, 4.0 * M_PI); } - return k.sum; + return kresult(k); } H3Error H3_EXPORT(geoLoopArea)(GeoLoop loop, double *out) { From ec26cb85f39618d649203399721e72fc5a8d0882 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Thu, 27 Nov 2025 11:41:43 -0500 Subject: [PATCH 06/62] notes --- src/apps/testapps/testArea.c | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/src/apps/testapps/testArea.c b/src/apps/testapps/testArea.c index 813b542d14..d927210ff8 100644 --- a/src/apps/testapps/testArea.c +++ b/src/apps/testapps/testArea.c @@ -12,15 +12,18 @@ // TODO: turn this into a benchmark, report accuracy and time // see if we get any benefit by avoiding the repeated sin/cos computation -// TODO: document algos: https://en.wikipedia.org/wiki/Kahan_summation_algorithm#Further_enhancements +// TODO: document algos: +// https://en.wikipedia.org/wiki/Kahan_summation_algorithm#Further_enhancements // TODO: better more general name for accumulator struct and functions // TOOD: demonstrate the area alg works for global polygons. +// TODO: is kadd faster without pointers? k = kadd(k, cellArea); + static void do_res_sum(int res) { Kahan k = {0, 0}; + IterCellsResolution iter = iterInitRes(res); - for (IterCellsResolution iter = iterInitRes(res); iter.h; - iterStepRes(&iter)) { + for (; iter.h; iterStepRes(&iter)) { double cellArea; H3_EXPORT(cellAreaRads2)(iter.h, &cellArea); From 74e277ff43312452ab8b6d0eaa922a2d38f4b579 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Thu, 27 Nov 2025 11:52:27 -0500 Subject: [PATCH 07/62] benchmark area --- CMakeLists.txt | 6 ++-- CMakeTests.cmake | 1 - justfile | 6 ++-- .../testArea.c => benchmarks/benchmarkArea.c} | 35 +++++-------------- 4 files changed, 16 insertions(+), 32 deletions(-) rename src/apps/{testapps/testArea.c => benchmarks/benchmarkArea.c} (56%) diff --git a/CMakeLists.txt b/CMakeLists.txt index 828c02bcef..7e4f88f7d9 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -280,7 +280,6 @@ set(OTHER_SOURCE_FILES src/apps/testapps/testH3IteratorsInternal.c src/apps/testapps/testMathExtensionsInternal.c src/apps/testapps/testDescribeH3Error.c - src/apps/testapps/testArea.c src/apps/miscapps/cellToBoundaryHier.c src/apps/miscapps/cellToLatLngHier.c src/apps/miscapps/generateBaseCellNeighbors.c @@ -320,7 +319,8 @@ set(OTHER_SOURCE_FILES src/apps/benchmarks/benchmarkDirectedEdge.c src/apps/benchmarks/benchmarkVertex.c src/apps/benchmarks/benchmarkIsValidCell.c - src/apps/benchmarks/benchmarkH3Api.c) + src/apps/benchmarks/benchmarkH3Api.c + src/apps/benchmarks/benchmarkArea.c) set(ALL_SOURCE_FILES ${LIB_SOURCE_FILES} ${APP_SOURCE_FILES} ${TEST_APP_SOURCE_FILES} @@ -675,6 +675,8 @@ if(BUILD_BENCHMARKS) src/apps/benchmarks/benchmarkPolygonToCells.c) add_h3_benchmark(benchmarkPolygonToCellsExperimental src/apps/benchmarks/benchmarkPolygonToCellsExperimental.c) + add_h3_benchmark(benchmarkArea + src/apps/benchmarks/benchmarkArea.c) if(ENABLE_REQUIRES_ALL_SYMBOLS) add_h3_benchmark(benchmarkPolygon src/apps/benchmarks/benchmarkPolygon.c) diff --git a/CMakeTests.cmake b/CMakeTests.cmake index cb705fc94a..1901a61fac 100644 --- a/CMakeTests.cmake +++ b/CMakeTests.cmake @@ -259,7 +259,6 @@ add_h3_test(testH3IteratorsInternal src/apps/testapps/testH3IteratorsInternal.c) add_h3_test(testMathExtensionsInternal src/apps/testapps/testMathExtensionsInternal.c) add_h3_test(testDescribeH3Error src/apps/testapps/testDescribeH3Error.c) -add_h3_test(testArea src/apps/testapps/testArea.c) add_h3_test_with_arg(testH3NeighborRotations src/apps/testapps/testH3NeighborRotations.c 0) diff --git a/justfile b/justfile index c9992e3ce8..5d7ff258c8 100644 --- a/justfile +++ b/justfile @@ -24,8 +24,9 @@ test: build # ./build/bin/testH3CellAreaExhaustive # ./build/bin/testEdgeCellsToPoly # ./build/bin/testDirectedEdge - ./build/bin/testArea + # ./build/bin/testArea # just test-slow + just test-fast time: time ./build/bin/testArea @@ -37,7 +38,8 @@ test-slow: build cd build; make test bench: build - ./build/bin/benchmarkCellsToLinkedMultiPolygon + ./build/bin/benchmarkArea + # ./build/bin/benchmarkCellsToLinkedMultiPolygon # ./build/bin/benchmarkCellsToMultiPolygon # ./build/bin/benchmarkDirectedEdge diff --git a/src/apps/testapps/testArea.c b/src/apps/benchmarks/benchmarkArea.c similarity index 56% rename from src/apps/testapps/testArea.c rename to src/apps/benchmarks/benchmarkArea.c index d927210ff8..1ffe6fa858 100644 --- a/src/apps/testapps/testArea.c +++ b/src/apps/benchmarks/benchmarkArea.c @@ -4,9 +4,9 @@ #include #include "adder.h" +#include "benchmark.h" #include "h3api.h" #include "iterators.h" -#include "test.h" #include "utility.h" // TODO: turn this into a benchmark, report accuracy and time @@ -35,31 +35,12 @@ static void do_res_sum(int res) { printf("res: %d, diff: %e\n", res, diff); } -SUITE(test_for_area) { - TEST(some_area_test) { - // the numerical test is how close to 4*pi do we get when adding up all - // cells at smaller and smaller resolution? +BEGIN_BENCHMARKS(); - printf("\n"); - - do_res_sum(0); - do_res_sum(1); - do_res_sum(2); - do_res_sum(3); - do_res_sum(4); - do_res_sum(5); - do_res_sum(6); - do_res_sum(7); - // do_res_sum(8); - - // res: 0, diff: 7.105427e-15 - // res: 1, diff: 3.907985e-14 - // res: 2, diff: 1.421085e-13 - // res: 3, diff: 4.689582e-13 - // res: 4, diff: 2.161826e-12 - // res: 5, diff: 1.543654e-12 - // res: 6, diff: 1.768363e-11 - // res: 7, diff: 1.490719e-11 - // res: 8, diff: 2.937917e-11 +BENCHMARK(directedEdgeToBoundary, 1, { + for (int i = 0; i < 7; i++) { + do_res_sum(i); } -} +}); + +END_BENCHMARKS(); From c915b5183011839bb211a305241454134ee8c759 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Thu, 27 Nov 2025 22:11:55 -0500 Subject: [PATCH 08/62] turn it to 11. --- src/apps/benchmarks/benchmarkArea.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/apps/benchmarks/benchmarkArea.c b/src/apps/benchmarks/benchmarkArea.c index 1ffe6fa858..ae006a09ef 100644 --- a/src/apps/benchmarks/benchmarkArea.c +++ b/src/apps/benchmarks/benchmarkArea.c @@ -38,7 +38,7 @@ static void do_res_sum(int res) { BEGIN_BENCHMARKS(); BENCHMARK(directedEdgeToBoundary, 1, { - for (int i = 0; i < 7; i++) { + for (int i = 0; i <= 8; i++) { do_res_sum(i); } }); From a74b9d949042e56873150cf71d3a9af828544386 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Thu, 27 Nov 2025 23:31:41 -0500 Subject: [PATCH 09/62] sadly, no improvement from reusing some trig --- src/apps/benchmarks/benchmarkArea.c | 2 +- src/h3lib/lib/area.c | 54 ++++++++++++++++++++++++++--- 2 files changed, 51 insertions(+), 5 deletions(-) diff --git a/src/apps/benchmarks/benchmarkArea.c b/src/apps/benchmarks/benchmarkArea.c index ae006a09ef..d246493941 100644 --- a/src/apps/benchmarks/benchmarkArea.c +++ b/src/apps/benchmarks/benchmarkArea.c @@ -38,7 +38,7 @@ static void do_res_sum(int res) { BEGIN_BENCHMARKS(); BENCHMARK(directedEdgeToBoundary, 1, { - for (int i = 0; i <= 8; i++) { + for (int i = 0; i <= 6; i++) { do_res_sum(i); } }); diff --git a/src/h3lib/lib/area.c b/src/h3lib/lib/area.c index 1eef4940e9..bd8b5551b5 100644 --- a/src/h3lib/lib/area.c +++ b/src/h3lib/lib/area.c @@ -15,7 +15,53 @@ void H3_EXPORT(destroyGeoLoop)(GeoLoop *loop) { loop->verts = NULL; } -static double cagnoli(LatLng x, LatLng y) { +typedef struct { + double s; + double c; + double lng; +} LatLngPre; + +static inline LatLngPre precompute_latlng(LatLng x) { + double lat = x.lat / 2.0 + M_PI / 4.0; + + LatLngPre llp = {.s = sin(lat), .c = cos(lat), .lng = x.lng}; + + return llp; +} + +static inline double cagnoli_pre(LatLngPre x, LatLngPre y) { + double sa = x.s * y.s; + double ca = x.c * y.c; + + double d = y.lng - x.lng; + double sd = sin(d); + double cd = cos(d); + + return -2.0 * atan2(sa * sd, sa * cd + ca); +} + +static inline double vertArea2(LatLng *verts, int numVerts) { + Kahan k = {0.0, 0.0}; + + LatLngPre x = precompute_latlng(verts[0]); + LatLngPre first = x; + LatLngPre y; + + for (int i = 1; i < numVerts; i++) { + y = precompute_latlng(verts[i]); + + kadd(&k, cagnoli_pre(x, y)); + x = y; + } + kadd(&k, cagnoli_pre(x, first)); + + if (k.sum < 0.0) { + kadd(&k, 4.0 * M_PI); + } + return kresult(k); +} + +static inline double cagnoli(LatLng x, LatLng y) { x.lat = x.lat / 2.0 + M_PI / 4.0; y.lat = y.lat / 2.0 + M_PI / 4.0; @@ -29,7 +75,7 @@ static double cagnoli(LatLng x, LatLng y) { return -2.0 * atan2(sa * sd, sa * cd + ca); } -static double vertArea(LatLng *verts, int numVerts) { +static inline double vertArea(LatLng *verts, int numVerts) { Kahan k = {0.0, 0.0}; for (int i = 0; i < numVerts; i++) { @@ -44,7 +90,7 @@ static double vertArea(LatLng *verts, int numVerts) { } H3Error H3_EXPORT(geoLoopArea)(GeoLoop loop, double *out) { - *out = vertArea(loop.verts, loop.numVerts); + *out = vertArea2(loop.verts, loop.numVerts); return E_SUCCESS; } @@ -67,7 +113,7 @@ H3Error H3_EXPORT(cellAreaRads2)(H3Index cell, double *out) { if (err) { return err; } - *out = vertArea(cb.verts, cb.numVerts); + *out = vertArea2(cb.verts, cb.numVerts); return E_SUCCESS; } From 98ecbee87cb17e308beb549854db55bdf1c5744e Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Thu, 27 Nov 2025 23:32:43 -0500 Subject: [PATCH 10/62] back to original --- src/h3lib/lib/area.c | 50 ++------------------------------------------ 1 file changed, 2 insertions(+), 48 deletions(-) diff --git a/src/h3lib/lib/area.c b/src/h3lib/lib/area.c index bd8b5551b5..ec7c39ecc4 100644 --- a/src/h3lib/lib/area.c +++ b/src/h3lib/lib/area.c @@ -15,52 +15,6 @@ void H3_EXPORT(destroyGeoLoop)(GeoLoop *loop) { loop->verts = NULL; } -typedef struct { - double s; - double c; - double lng; -} LatLngPre; - -static inline LatLngPre precompute_latlng(LatLng x) { - double lat = x.lat / 2.0 + M_PI / 4.0; - - LatLngPre llp = {.s = sin(lat), .c = cos(lat), .lng = x.lng}; - - return llp; -} - -static inline double cagnoli_pre(LatLngPre x, LatLngPre y) { - double sa = x.s * y.s; - double ca = x.c * y.c; - - double d = y.lng - x.lng; - double sd = sin(d); - double cd = cos(d); - - return -2.0 * atan2(sa * sd, sa * cd + ca); -} - -static inline double vertArea2(LatLng *verts, int numVerts) { - Kahan k = {0.0, 0.0}; - - LatLngPre x = precompute_latlng(verts[0]); - LatLngPre first = x; - LatLngPre y; - - for (int i = 1; i < numVerts; i++) { - y = precompute_latlng(verts[i]); - - kadd(&k, cagnoli_pre(x, y)); - x = y; - } - kadd(&k, cagnoli_pre(x, first)); - - if (k.sum < 0.0) { - kadd(&k, 4.0 * M_PI); - } - return kresult(k); -} - static inline double cagnoli(LatLng x, LatLng y) { x.lat = x.lat / 2.0 + M_PI / 4.0; y.lat = y.lat / 2.0 + M_PI / 4.0; @@ -90,7 +44,7 @@ static inline double vertArea(LatLng *verts, int numVerts) { } H3Error H3_EXPORT(geoLoopArea)(GeoLoop loop, double *out) { - *out = vertArea2(loop.verts, loop.numVerts); + *out = vertArea(loop.verts, loop.numVerts); return E_SUCCESS; } @@ -113,7 +67,7 @@ H3Error H3_EXPORT(cellAreaRads2)(H3Index cell, double *out) { if (err) { return err; } - *out = vertArea2(cb.verts, cb.numVerts); + *out = vertArea(cb.verts, cb.numVerts); return E_SUCCESS; } From db21eb0bc0ba056a474f873c038dc8145546a5ff Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Fri, 28 Nov 2025 00:12:41 -0500 Subject: [PATCH 11/62] simplify --- src/h3lib/lib/area.c | 30 ++++++++++++++++-------------- 1 file changed, 16 insertions(+), 14 deletions(-) diff --git a/src/h3lib/lib/area.c b/src/h3lib/lib/area.c index ec7c39ecc4..decc096afd 100644 --- a/src/h3lib/lib/area.c +++ b/src/h3lib/lib/area.c @@ -10,11 +10,6 @@ #include "alloc.h" #include "h3api.h" -void H3_EXPORT(destroyGeoLoop)(GeoLoop *loop) { - H3_MEMORY(free)(loop->verts); - loop->verts = NULL; -} - static inline double cagnoli(LatLng x, LatLng y) { x.lat = x.lat / 2.0 + M_PI / 4.0; y.lat = y.lat / 2.0 + M_PI / 4.0; @@ -29,22 +24,19 @@ static inline double cagnoli(LatLng x, LatLng y) { return -2.0 * atan2(sa * sd, sa * cd + ca); } -static inline double vertArea(LatLng *verts, int numVerts) { +H3Error H3_EXPORT(geoLoopArea)(GeoLoop loop, double *out) { Kahan k = {0.0, 0.0}; - for (int i = 0; i < numVerts; i++) { - int j = (i + 1) % numVerts; - kadd(&k, cagnoli(verts[i], verts[j])); + for (int i = 0; i < loop.numVerts; i++) { + int j = (i + 1) % loop.numVerts; + kadd(&k, cagnoli(loop.verts[i], loop.verts[j])); } if (k.sum < 0.0) { kadd(&k, 4.0 * M_PI); } - return kresult(k); -} -H3Error H3_EXPORT(geoLoopArea)(GeoLoop loop, double *out) { - *out = vertArea(loop.verts, loop.numVerts); + *out = kresult(k); return E_SUCCESS; } @@ -67,7 +59,17 @@ H3Error H3_EXPORT(cellAreaRads2)(H3Index cell, double *out) { if (err) { return err; } - *out = vertArea(cb.verts, cb.numVerts); + + GeoLoop loop = {.verts = cb.verts, .numVerts = cb.numVerts}; + err = H3_EXPORT(geoLoopArea)(loop, out); + if (err) { + return err; + } return E_SUCCESS; } + +void H3_EXPORT(destroyGeoLoop)(GeoLoop *loop) { + H3_MEMORY(free)(loop->verts); + loop->verts = NULL; +} From 04640ef8b648f90e8cdfcafce3ee7e80e770fbb5 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Fri, 28 Nov 2025 20:49:25 -0500 Subject: [PATCH 12/62] try adding constants --- src/h3lib/lib/area.c | 1 + 1 file changed, 1 insertion(+) diff --git a/src/h3lib/lib/area.c b/src/h3lib/lib/area.c index decc096afd..cec3c90ddd 100644 --- a/src/h3lib/lib/area.c +++ b/src/h3lib/lib/area.c @@ -9,6 +9,7 @@ #include "adder.h" #include "alloc.h" #include "h3api.h" +#include "constants.h" static inline double cagnoli(LatLng x, LatLng y) { x.lat = x.lat / 2.0 + M_PI / 4.0; From 2349ab4c333a35b9bd8aff560be60b0a3933e949 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Fri, 28 Nov 2025 20:53:27 -0500 Subject: [PATCH 13/62] format, we must --- src/h3lib/lib/area.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/h3lib/lib/area.c b/src/h3lib/lib/area.c index cec3c90ddd..203f039c9a 100644 --- a/src/h3lib/lib/area.c +++ b/src/h3lib/lib/area.c @@ -8,8 +8,8 @@ #include "adder.h" #include "alloc.h" -#include "h3api.h" #include "constants.h" +#include "h3api.h" static inline double cagnoli(LatLng x, LatLng y) { x.lat = x.lat / 2.0 + M_PI / 4.0; From 0dd6a5eec25bfb3a3ff93f934a9d52c36c4b11a7 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Fri, 28 Nov 2025 22:10:07 -0500 Subject: [PATCH 14/62] Adder docs --- src/apps/benchmarks/benchmarkArea.c | 17 ++--- src/apps/testapps/testH3CellAreaExhaustive.c | 6 +- src/h3lib/include/adder.h | 78 ++++++++++++++------ src/h3lib/lib/area.c | 10 +-- 4 files changed, 67 insertions(+), 44 deletions(-) diff --git a/src/apps/benchmarks/benchmarkArea.c b/src/apps/benchmarks/benchmarkArea.c index d246493941..69ce545cfb 100644 --- a/src/apps/benchmarks/benchmarkArea.c +++ b/src/apps/benchmarks/benchmarkArea.c @@ -9,29 +9,22 @@ #include "iterators.h" #include "utility.h" -// TODO: turn this into a benchmark, report accuracy and time -// see if we get any benefit by avoiding the repeated sin/cos computation - -// TODO: document algos: -// https://en.wikipedia.org/wiki/Kahan_summation_algorithm#Further_enhancements -// TODO: better more general name for accumulator struct and functions // TOOD: demonstrate the area alg works for global polygons. // TODO: is kadd faster without pointers? k = kadd(k, cellArea); static void do_res_sum(int res) { - Kahan k = {0, 0}; + Adder adder = {0, 0}; + double cellArea; IterCellsResolution iter = iterInitRes(res); for (; iter.h; iterStepRes(&iter)) { - double cellArea; H3_EXPORT(cellAreaRads2)(iter.h, &cellArea); - - kadd(&k, cellArea); + kadd(&adder, cellArea); } - double area = kresult(k); - double diff = fabs(area - 4 * M_PI); + double total_area = kresult(adder); + double diff = fabs(total_area - 4 * M_PI); printf("res: %d, diff: %e\n", res, diff); } diff --git a/src/apps/testapps/testH3CellAreaExhaustive.c b/src/apps/testapps/testH3CellAreaExhaustive.c index 3f39fac7bb..b464014b46 100644 --- a/src/apps/testapps/testH3CellAreaExhaustive.c +++ b/src/apps/testapps/testH3CellAreaExhaustive.c @@ -137,16 +137,16 @@ static void cell_area_assert(H3Index cell) { */ static void earth_area_test(int res, H3Error (*cell_area)(H3Index, double *), double target, double tol) { - Kahan k = {0.0, 0.0}; + Adder adder = {0.0, 0.0}; IterCellsResolution iter = iterInitRes(res); for (; iter.h; iterStepRes(&iter)) { double cellArea; t_assertSuccess((*cell_area)(iter.h, &cellArea)); - kadd(&k, cellArea); + kadd(&adder, cellArea); } - t_assert(fabs(k.sum - target) < tol, + t_assert(fabs(kresult(adder) - target) < tol, "sum of all cells should give earth area"); } diff --git a/src/h3lib/include/adder.h b/src/h3lib/include/adder.h index a9763be3ed..e391987964 100644 --- a/src/h3lib/include/adder.h +++ b/src/h3lib/include/adder.h @@ -1,46 +1,76 @@ #ifndef ADDER_H #define ADDER_H +/* +This is a header-only implementation of a "compensated summation" algorithm, +which allows us to add up sequences of floating point numbers with better +precision than naive summation, especially when the terms in the sum +can vary in magnitude. +See: https://en.wikipedia.org/wiki/Kahan_summation_algorithm + +This is useful when computing the area of (multi-)polygons, which can involve +adding many small terms to large aggregates. +For example, D3 uses this when computing polygonal areas via the `Kahan` class: +https://github.com/d3/d3-geo/blob/main/src/area.js + +There are a few potential algorithms we might consider for summation: + +1. Naive sum +2. Kahan summation +3. Neumaier summation +4. Other approaches like pairwise summation, or Pythons `fsum` + +We considered the first three for simplicity, and settled on Kahan summation, +as it acheives noticably better accuracy than naive summation and almost as good +accuracy as Neumaier, while only being slightly slower than naive, but slightly +faster and simpler than Neumaier. + +See also: https://github.com/python/cpython/issues/100425 for discussions of +tradeoffs between Kahan, Neumaier, and `fsum`. +*/ + typedef struct { double sum; // running total double c; // compensation term -} Kahan; +} Adder; -static inline void simple_add(Kahan *k, double x) { k->sum += x; } +static inline void simple_add(Adder *adder, double x) { adder->sum += x; } -static inline void kahan_add(Kahan *k, double x) { - double y = x - k->c; - double t = k->sum + y; - k->c = (t - k->sum) - y; - k->sum = t; +static inline void kahan_add(Adder *adder, double x) { + double y = x - adder->c; + double t = adder->sum + y; + adder->c = (t - adder->sum) - y; + adder->sum = t; } -static inline void neumaier_add(Kahan *k, double x) { - double t = k->sum + x; +static inline void neumaier_add(Adder *adder, double x) { + double t = adder->sum + x; - if (fabs(k->sum) >= fabs(x)) { - k->c += (k->sum - t) + x; + if (fabs(adder->sum) >= fabs(x)) { + adder->c += (adder->sum - t) + x; } else { - k->c += (x - t) + k->sum; + adder->c += (x - t) + adder->sum; } - k->sum = t; + adder->sum = t; } -static inline double simple_result(Kahan k) { return k.sum; } -static inline double kahan_result(Kahan k) { return k.sum; } -static inline double neumaier_result(Kahan k) { return k.sum + k.c; } +static inline double simple_result(Adder adder) { return adder.sum; } +static inline double kahan_result(Adder adder) { return adder.sum; } +static inline double neumaier_result(Adder adder) { + return adder.sum + adder.c; +} -static inline void kadd(Kahan *k, double x) { - // simple_add(k, x); - // kahan_add(k, x); - neumaier_add(k, x); +static inline void kadd(Adder *adder, double x) { + // simple_add(adder, x); + // kahan_add(adder, x); + neumaier_add(adder, x); } -static inline double kresult(Kahan k) { - // return simple_result(k); - // return kahan_result(k); - return neumaier_result(k); +static inline double kresult(Adder adder) { + // return simple_result(adder); + // return kahan_result(adder); + return neumaier_result(adder); } #endif diff --git a/src/h3lib/lib/area.c b/src/h3lib/lib/area.c index 203f039c9a..a10e337dcf 100644 --- a/src/h3lib/lib/area.c +++ b/src/h3lib/lib/area.c @@ -26,18 +26,18 @@ static inline double cagnoli(LatLng x, LatLng y) { } H3Error H3_EXPORT(geoLoopArea)(GeoLoop loop, double *out) { - Kahan k = {0.0, 0.0}; + Adder adder = {0.0, 0.0}; for (int i = 0; i < loop.numVerts; i++) { int j = (i + 1) % loop.numVerts; - kadd(&k, cagnoli(loop.verts[i], loop.verts[j])); + kadd(&adder, cagnoli(loop.verts[i], loop.verts[j])); } - if (k.sum < 0.0) { - kadd(&k, 4.0 * M_PI); + if (adder.sum < 0.0) { + kadd(&adder, 4.0 * M_PI); } - *out = kresult(k); + *out = kresult(adder); return E_SUCCESS; } From 7528234979b9e5492da5bc1effe0199b9b4df595 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Fri, 28 Nov 2025 22:25:10 -0500 Subject: [PATCH 15/62] bench --- src/apps/benchmarks/benchmarkArea.c | 7 ++----- src/h3lib/lib/area.c | 4 ++++ 2 files changed, 6 insertions(+), 5 deletions(-) diff --git a/src/apps/benchmarks/benchmarkArea.c b/src/apps/benchmarks/benchmarkArea.c index 69ce545cfb..87a05e1e0a 100644 --- a/src/apps/benchmarks/benchmarkArea.c +++ b/src/apps/benchmarks/benchmarkArea.c @@ -9,10 +9,6 @@ #include "iterators.h" #include "utility.h" -// TOOD: demonstrate the area alg works for global polygons. - -// TODO: is kadd faster without pointers? k = kadd(k, cellArea); - static void do_res_sum(int res) { Adder adder = {0, 0}; double cellArea; @@ -31,7 +27,8 @@ static void do_res_sum(int res) { BEGIN_BENCHMARKS(); BENCHMARK(directedEdgeToBoundary, 1, { - for (int i = 0; i <= 6; i++) { + int MAX_RES = 8; + for (int i = 0; i <= MAX_RES; i++) { do_res_sum(i); } }); diff --git a/src/h3lib/lib/area.c b/src/h3lib/lib/area.c index a10e337dcf..6af6632122 100644 --- a/src/h3lib/lib/area.c +++ b/src/h3lib/lib/area.c @@ -11,6 +11,10 @@ #include "constants.h" #include "h3api.h" +// TOOD: demonstrate the area alg works for global polygons. + +// TODO: is kadd faster without pointers? k = kadd(k, cellArea); + static inline double cagnoli(LatLng x, LatLng y) { x.lat = x.lat / 2.0 + M_PI / 4.0; y.lat = y.lat / 2.0 + M_PI / 4.0; From 4a4fb073b4eebcd9237a0489e8fea1b717ccee3b Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Fri, 28 Nov 2025 22:48:05 -0500 Subject: [PATCH 16/62] adder docs --- src/h3lib/include/adder.h | 35 ++++++++++++++++++----------------- 1 file changed, 18 insertions(+), 17 deletions(-) diff --git a/src/h3lib/include/adder.h b/src/h3lib/include/adder.h index e391987964..9b8ad25e07 100644 --- a/src/h3lib/include/adder.h +++ b/src/h3lib/include/adder.h @@ -2,15 +2,16 @@ #define ADDER_H /* -This is a header-only implementation of a "compensated summation" algorithm, -which allows us to add up sequences of floating point numbers with better -precision than naive summation, especially when the terms in the sum -can vary in magnitude. +Header-only implementation of a "compensated summation" algorithm (Kahan +summation), which allows us to add up sequences of floating-point numbers +with better precision than naive summation, especially when the terms in +the sum vary significantly in magnitude. See: https://en.wikipedia.org/wiki/Kahan_summation_algorithm -This is useful when computing the area of (multi-)polygons, which can involve -adding many small terms to large aggregates. -For example, D3 uses this when computing polygonal areas via the `Kahan` class: +This is useful when computing the area of (multi)polygons, which +often involves adding many small terms to a large aggregate. For example, +D3 uses an improved accuracy summation when computing polygonal areas via its +`Adder` class: https://github.com/d3/d3-geo/blob/main/src/area.js There are a few potential algorithms we might consider for summation: @@ -18,14 +19,14 @@ There are a few potential algorithms we might consider for summation: 1. Naive sum 2. Kahan summation 3. Neumaier summation -4. Other approaches like pairwise summation, or Pythons `fsum` +4. Other approaches like pairwise summation, or Python's `fsum` -We considered the first three for simplicity, and settled on Kahan summation, -as it acheives noticably better accuracy than naive summation and almost as good -accuracy as Neumaier, while only being slightly slower than naive, but slightly -faster and simpler than Neumaier. +We considered the first three for simplicity, and settled on Kahan summation: +it achieves noticeably better accuracy than naive summation and almost as +good accuracy as Neumaier, while being only slightly slower than naive and +slightly faster and simpler than Neumaier. -See also: https://github.com/python/cpython/issues/100425 for discussions of +See also: https://github.com/python/cpython/issues/100425 for discussion of tradeoffs between Kahan, Neumaier, and `fsum`. */ @@ -62,15 +63,15 @@ static inline double neumaier_result(Adder adder) { } static inline void kadd(Adder *adder, double x) { - // simple_add(adder, x); + simple_add(adder, x); // kahan_add(adder, x); - neumaier_add(adder, x); + // neumaier_add(adder, x); } static inline double kresult(Adder adder) { - // return simple_result(adder); + return simple_result(adder); // return kahan_result(adder); - return neumaier_result(adder); + // return neumaier_result(adder); } #endif From 6bd329904066292c54ee395de901b651c1d3fde3 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Fri, 28 Nov 2025 23:03:18 -0500 Subject: [PATCH 17/62] Settled on Kahan implementation --- src/apps/benchmarks/benchmarkArea.c | 2 +- src/apps/testapps/testH3CellAreaExhaustive.c | 2 +- src/h3lib/include/adder.h | 34 +------------------- src/h3lib/lib/area.c | 2 +- 4 files changed, 4 insertions(+), 36 deletions(-) diff --git a/src/apps/benchmarks/benchmarkArea.c b/src/apps/benchmarks/benchmarkArea.c index 87a05e1e0a..e9fa852d2d 100644 --- a/src/apps/benchmarks/benchmarkArea.c +++ b/src/apps/benchmarks/benchmarkArea.c @@ -19,7 +19,7 @@ static void do_res_sum(int res) { kadd(&adder, cellArea); } - double total_area = kresult(adder); + double total_area = adder.sum; double diff = fabs(total_area - 4 * M_PI); printf("res: %d, diff: %e\n", res, diff); } diff --git a/src/apps/testapps/testH3CellAreaExhaustive.c b/src/apps/testapps/testH3CellAreaExhaustive.c index b464014b46..4b0dc2b6b3 100644 --- a/src/apps/testapps/testH3CellAreaExhaustive.c +++ b/src/apps/testapps/testH3CellAreaExhaustive.c @@ -146,7 +146,7 @@ static void earth_area_test(int res, H3Error (*cell_area)(H3Index, double *), kadd(&adder, cellArea); } - t_assert(fabs(kresult(adder) - target) < tol, + t_assert(fabs(adder.sum - target) < tol, "sum of all cells should give earth area"); } diff --git a/src/h3lib/include/adder.h b/src/h3lib/include/adder.h index 9b8ad25e07..1a6fc79695 100644 --- a/src/h3lib/include/adder.h +++ b/src/h3lib/include/adder.h @@ -35,43 +35,11 @@ typedef struct { double c; // compensation term } Adder; -static inline void simple_add(Adder *adder, double x) { adder->sum += x; } - -static inline void kahan_add(Adder *adder, double x) { +static inline void kadd(Adder *adder, double x) { double y = x - adder->c; double t = adder->sum + y; adder->c = (t - adder->sum) - y; adder->sum = t; } -static inline void neumaier_add(Adder *adder, double x) { - double t = adder->sum + x; - - if (fabs(adder->sum) >= fabs(x)) { - adder->c += (adder->sum - t) + x; - } else { - adder->c += (x - t) + adder->sum; - } - - adder->sum = t; -} - -static inline double simple_result(Adder adder) { return adder.sum; } -static inline double kahan_result(Adder adder) { return adder.sum; } -static inline double neumaier_result(Adder adder) { - return adder.sum + adder.c; -} - -static inline void kadd(Adder *adder, double x) { - simple_add(adder, x); - // kahan_add(adder, x); - // neumaier_add(adder, x); -} - -static inline double kresult(Adder adder) { - return simple_result(adder); - // return kahan_result(adder); - // return neumaier_result(adder); -} - #endif diff --git a/src/h3lib/lib/area.c b/src/h3lib/lib/area.c index 6af6632122..e026ddd525 100644 --- a/src/h3lib/lib/area.c +++ b/src/h3lib/lib/area.c @@ -41,7 +41,7 @@ H3Error H3_EXPORT(geoLoopArea)(GeoLoop loop, double *out) { kadd(&adder, 4.0 * M_PI); } - *out = kresult(adder); + *out = adder.sum; return E_SUCCESS; } From 1b21bed993ef58a0ae24301fa91a8c0b6433c599 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Fri, 28 Nov 2025 23:07:42 -0500 Subject: [PATCH 18/62] Tighten tolerances in testH3CellAreaExhaustive.c due to compensated sum --- src/apps/testapps/testH3CellAreaExhaustive.c | 26 +++++++++----------- 1 file changed, 11 insertions(+), 15 deletions(-) diff --git a/src/apps/testapps/testH3CellAreaExhaustive.c b/src/apps/testapps/testH3CellAreaExhaustive.c index 4b0dc2b6b3..315eb1e91e 100644 --- a/src/apps/testapps/testH3CellAreaExhaustive.c +++ b/src/apps/testapps/testH3CellAreaExhaustive.c @@ -178,28 +178,24 @@ SUITE(h3CellAreaExhaustive) { double km2 = rads2 * EARTH_RADIUS_KM * EARTH_RADIUS_KM; double m2 = km2 * 1000 * 1000; - // Notice the drop in accuracy at resolution 1. - // I think this has something to do with Class II vs Class III - // resolutions. - earth_area_test(0, H3_EXPORT(cellAreaRads2), rads2, 1e-14); earth_area_test(0, H3_EXPORT(cellAreaKm2), km2, 1e-6); earth_area_test(0, H3_EXPORT(cellAreaM2), m2, 1e0); - earth_area_test(1, H3_EXPORT(cellAreaRads2), rads2, 1e-9); - earth_area_test(1, H3_EXPORT(cellAreaKm2), km2, 1e-1); - earth_area_test(1, H3_EXPORT(cellAreaM2), m2, 1e5); + earth_area_test(1, H3_EXPORT(cellAreaRads2), rads2, 1e-14); + earth_area_test(1, H3_EXPORT(cellAreaKm2), km2, 1e-6); + earth_area_test(1, H3_EXPORT(cellAreaM2), m2, 1e0); - earth_area_test(2, H3_EXPORT(cellAreaRads2), rads2, 1e-12); - earth_area_test(2, H3_EXPORT(cellAreaKm2), km2, 1e-5); + earth_area_test(2, H3_EXPORT(cellAreaRads2), rads2, 1e-14); + earth_area_test(2, H3_EXPORT(cellAreaKm2), km2, 1e-6); earth_area_test(2, H3_EXPORT(cellAreaM2), m2, 1e0); - earth_area_test(3, H3_EXPORT(cellAreaRads2), rads2, 1e-11); - earth_area_test(3, H3_EXPORT(cellAreaKm2), km2, 1e-3); - earth_area_test(3, H3_EXPORT(cellAreaM2), m2, 1e3); + earth_area_test(3, H3_EXPORT(cellAreaRads2), rads2, 1e-14); + earth_area_test(3, H3_EXPORT(cellAreaKm2), km2, 1e-6); + earth_area_test(3, H3_EXPORT(cellAreaM2), m2, 1e0); - earth_area_test(4, H3_EXPORT(cellAreaRads2), rads2, 1e-11); - earth_area_test(4, H3_EXPORT(cellAreaKm2), km2, 1e-3); - earth_area_test(4, H3_EXPORT(cellAreaM2), m2, 1e2); + earth_area_test(4, H3_EXPORT(cellAreaRads2), rads2, 1e-14); + earth_area_test(4, H3_EXPORT(cellAreaKm2), km2, 1e-6); + earth_area_test(4, H3_EXPORT(cellAreaM2), m2, 1e0); } } From 64853b3c011c6af823d5507cf6afc44406659513 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Fri, 28 Nov 2025 23:31:33 -0500 Subject: [PATCH 19/62] 6 --- src/apps/benchmarks/benchmarkArea.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/apps/benchmarks/benchmarkArea.c b/src/apps/benchmarks/benchmarkArea.c index e9fa852d2d..a29a771680 100644 --- a/src/apps/benchmarks/benchmarkArea.c +++ b/src/apps/benchmarks/benchmarkArea.c @@ -27,7 +27,7 @@ static void do_res_sum(int res) { BEGIN_BENCHMARKS(); BENCHMARK(directedEdgeToBoundary, 1, { - int MAX_RES = 8; + int MAX_RES = 6; for (int i = 0; i <= MAX_RES; i++) { do_res_sum(i); } From 3547c47d829e2bf611d8fb4f7313f7f87b373a62 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Fri, 28 Nov 2025 23:33:51 -0500 Subject: [PATCH 20/62] notes --- src/h3lib/include/adder.h | 6 ++++++ src/h3lib/lib/area.c | 2 -- 2 files changed, 6 insertions(+), 2 deletions(-) diff --git a/src/h3lib/include/adder.h b/src/h3lib/include/adder.h index 1a6fc79695..53e31ac5f8 100644 --- a/src/h3lib/include/adder.h +++ b/src/h3lib/include/adder.h @@ -28,6 +28,12 @@ slightly faster and simpler than Neumaier. See also: https://github.com/python/cpython/issues/100425 for discussion of tradeoffs between Kahan, Neumaier, and `fsum`. + +## Usage + +Adder adder = {0.0, 0.0}; // initialize +kadd(&adder, x); // add x to summation +out = adder.sum; // extract final sum */ typedef struct { diff --git a/src/h3lib/lib/area.c b/src/h3lib/lib/area.c index e026ddd525..0a43bcf32a 100644 --- a/src/h3lib/lib/area.c +++ b/src/h3lib/lib/area.c @@ -13,8 +13,6 @@ // TOOD: demonstrate the area alg works for global polygons. -// TODO: is kadd faster without pointers? k = kadd(k, cellArea); - static inline double cagnoli(LatLng x, LatLng y) { x.lat = x.lat / 2.0 + M_PI / 4.0; y.lat = y.lat / 2.0 + M_PI / 4.0; From 0c73921821d631af1d159ab7fc21000db9348612 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Sat, 29 Nov 2025 09:36:28 -0500 Subject: [PATCH 21/62] starting new tests --- CMakeLists.txt | 1 + CMakeTests.cmake | 1 + justfile | 3 ++- src/apps/testapps/testGeoLoopArea.c | 32 +++++++++++++++++++++++++++++ src/h3lib/lib/area.c | 2 ++ 5 files changed, 38 insertions(+), 1 deletion(-) create mode 100644 src/apps/testapps/testGeoLoopArea.c diff --git a/CMakeLists.txt b/CMakeLists.txt index 7e4f88f7d9..dcc7069a46 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -280,6 +280,7 @@ set(OTHER_SOURCE_FILES src/apps/testapps/testH3IteratorsInternal.c src/apps/testapps/testMathExtensionsInternal.c src/apps/testapps/testDescribeH3Error.c + src/apps/testapps/testGeoLoopArea.c src/apps/miscapps/cellToBoundaryHier.c src/apps/miscapps/cellToLatLngHier.c src/apps/miscapps/generateBaseCellNeighbors.c diff --git a/CMakeTests.cmake b/CMakeTests.cmake index 1901a61fac..3508ff23a3 100644 --- a/CMakeTests.cmake +++ b/CMakeTests.cmake @@ -259,6 +259,7 @@ add_h3_test(testH3IteratorsInternal src/apps/testapps/testH3IteratorsInternal.c) add_h3_test(testMathExtensionsInternal src/apps/testapps/testMathExtensionsInternal.c) add_h3_test(testDescribeH3Error src/apps/testapps/testDescribeH3Error.c) +add_h3_test(testGeoLoopArea src/apps/testapps/testGeoLoopArea.c) add_h3_test_with_arg(testH3NeighborRotations src/apps/testapps/testH3NeighborRotations.c 0) diff --git a/justfile b/justfile index 5d7ff258c8..8df37001ba 100644 --- a/justfile +++ b/justfile @@ -26,7 +26,8 @@ test: build # ./build/bin/testDirectedEdge # ./build/bin/testArea # just test-slow - just test-fast + # just test-fast + ./build/bin/testGeoLoopArea time: time ./build/bin/testArea diff --git a/src/apps/testapps/testGeoLoopArea.c b/src/apps/testapps/testGeoLoopArea.c new file mode 100644 index 0000000000..f423c1db75 --- /dev/null +++ b/src/apps/testapps/testGeoLoopArea.c @@ -0,0 +1,32 @@ +/* + * Copyright 2020-2021 Uber Technologies, Inc. + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +/** @file + * @brief tests H3 cell area functions on a few specific cases + * + * usage: `testH3CellArea` + */ + +#include + +#include "h3Index.h" +#include "test.h" + +SUITE(h3CellArea) { + TEST(specific_cell_area) { + t_assert(false, "should fail"); + } +} diff --git a/src/h3lib/lib/area.c b/src/h3lib/lib/area.c index 0a43bcf32a..bdbbc0f351 100644 --- a/src/h3lib/lib/area.c +++ b/src/h3lib/lib/area.c @@ -14,6 +14,7 @@ // TOOD: demonstrate the area alg works for global polygons. static inline double cagnoli(LatLng x, LatLng y) { + // https://github.com/d3/d3-geo/blob/8c53a90ae70c94bace73ecb02f2c792c649c86ba/src/area.js#L51-L70 x.lat = x.lat / 2.0 + M_PI / 4.0; y.lat = y.lat / 2.0 + M_PI / 4.0; @@ -74,5 +75,6 @@ H3Error H3_EXPORT(cellAreaRads2)(H3Index cell, double *out) { void H3_EXPORT(destroyGeoLoop)(GeoLoop *loop) { H3_MEMORY(free)(loop->verts); + loop->numVerts = 0; loop->verts = NULL; } From 31fd41644f3576efc1d22f97fdfc6ffc67a14627 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Sat, 29 Nov 2025 18:49:08 -0500 Subject: [PATCH 22/62] tests work --- src/apps/testapps/testGeoLoopArea.c | 141 +++++++++++++++++++++++++++- 1 file changed, 138 insertions(+), 3 deletions(-) diff --git a/src/apps/testapps/testGeoLoopArea.c b/src/apps/testapps/testGeoLoopArea.c index f423c1db75..985a48b324 100644 --- a/src/apps/testapps/testGeoLoopArea.c +++ b/src/apps/testapps/testGeoLoopArea.c @@ -24,9 +24,144 @@ #include "h3Index.h" #include "test.h" +#include "utility.h" -SUITE(h3CellArea) { - TEST(specific_cell_area) { - t_assert(false, "should fail"); +// GeoMultiPolygon alloc_global_mpoly() { +// const int numPolygons = 8; +// const int numVerts = 3; +// const LatLng verts[numPolygons][numVerts] = { +// {{M_PI_2, 0.0}, {0.0, 0.0}, {0.0, M_PI_2}}, +// {{M_PI_2, 0.0}, {0.0, M_PI_2}, {0.0, M_PI}}, +// {{M_PI_2, 0.0}, {0.0, M_PI}, {0.0, -M_PI_2}}, +// {{M_PI_2, 0.0}, {0.0, -M_PI_2}, {0.0, 0.0}}, +// {{-M_PI_2, 0.0}, {0.0, 0.0}, {0.0, -M_PI_2}}, +// {{-M_PI_2, 0.0}, {0.0, -M_PI_2}, {0.0, -M_PI}}, +// {{-M_PI_2, 0.0}, {0.0, -M_PI}, {0.0, M_PI_2}}, +// {{-M_PI_2, 0.0}, {0.0, M_PI_2}, {0.0, 0.0}}, +// }; + +// GeoMultiPolygon mpoly = { +// .numPolygons = numPolygons, +// .polygons = H3_MEMORY(malloc)(sizeof(GeoPolygon) * numPolygons), +// }; + +// for (int i = 0; i < numPolygons; i++) { +// GeoPolygon *poly = &mpoly.polygons[i]; +// poly->numHoles = 0; +// poly->holes = NULL; +// poly->geoloop.numVerts = numVerts; +// poly->geoloop.verts = H3_MEMORY(malloc)(sizeof(LatLng) * numVerts); + +// for (int j = 0; j < numVerts; j++) { +// poly->geoloop.verts[j] = verts[i][j]; +// } +// } + +// qsort(mpoly.polygons, numPolygons, sizeof(GeoPolygon), cmp_poly_area); + +// return mpoly; +// } + +#define TOL 1e-14 + +static void helper(LatLng *verts, int numVerts, double target_area) { + GeoLoop loop = {.verts = verts, .numVerts = numVerts}; + + double out; + t_assertSuccess(H3_EXPORT(geoLoopArea)(loop, &out)); + + t_assert(fabs(out - target_area) < TOL, "area should match"); +} + +SUITE(geoLoopArea) { + TEST(triangle) { + // GeoLoop representing a triangle covering 1/8 of the globe, with + // points ordered according to right-hand rule (counter-clockwise). + // Expected area: pi/2 + LatLng verts[] = {{M_PI_2, 0.0}, {0.0, 0.0}, {0.0, M_PI_2}}; + GeoLoop loop = {.verts = verts, .numVerts = ARRAY_SIZE(verts)}; + + double out; + t_assertSuccess(H3_EXPORT(geoLoopArea)(loop, &out)); + + t_assert(fabs(out - M_PI / 2.0) < 1e-14, "area should match"); + } + TEST(reverse_triangle) { + // Reversed the order of the points, so they are in clockwise order. + // This GeoLoop represents the whole globe minus the triangle above. + // Expected area: 7*pi/2 + LatLng verts[] = {{0.0, M_PI_2}, {0.0, 0.0}, {M_PI_2, 0.0}}; + GeoLoop loop = {.verts = verts, .numVerts = ARRAY_SIZE(verts)}; + + double out; + t_assertSuccess(H3_EXPORT(geoLoopArea)(loop, &out)); + + t_assert(fabs(out - 7 * M_PI / 2.0) < 1e-14, "area should match"); + } + + TEST(slice) { + // 1/4 slice of the globe, from north to south pole + // Expected area: pi + LatLng verts[] = {{M_PI_2, 0.0}, + {0.0, 0.0}, + {-M_PI_2, 0.0}, + {0.0, M_PI_2}, + {M_PI_2, 0.0}}; + GeoLoop loop = {.verts = verts, .numVerts = ARRAY_SIZE(verts)}; + + double out; + t_assertSuccess(H3_EXPORT(geoLoopArea)(loop, &out)); + + t_assert(fabs(out - M_PI) < 1e-14, "area should match"); + } + + TEST(reverse_slice) { + // 3/4 slice of the globe, from north to south pole, formed by reversing + // order of points from example above. + // Expected area: 3*pi + LatLng verts[] = {{M_PI_2, 0.0}, + {0.0, M_PI_2}, + {-M_PI_2, 0.0}, + {0.0, 0.0}, + {M_PI_2, 0.0}}; + GeoLoop loop = {.verts = verts, .numVerts = ARRAY_SIZE(verts)}; + + double out; + t_assertSuccess(H3_EXPORT(geoLoopArea)(loop, &out)); + + t_assert(fabs(out - 3 * M_PI) < 1e-14, "area should match"); + } + + TEST(hemisphere_east) { + // one hemisphere of globe. western or eastern? + // Expected area: 2*pi + LatLng verts[] = { + {M_PI_2, 0.0}, // north pole + {0.0, -M_PI_2}, // equator + {-M_PI_2, 0.0}, // south pole + {0.0, M_PI_2}, // equator + }; + GeoLoop loop = {.verts = verts, .numVerts = ARRAY_SIZE(verts)}; + + double out; + t_assertSuccess(H3_EXPORT(geoLoopArea)(loop, &out)); + + t_assert(fabs(out - 2 * M_PI) < 1e-14, "area should match"); + } + + TEST(hemisphere_north) { + // one hemisphere of globe. northern + // Expected area: 2*pi + LatLng verts[] = { + {0.0, -M_PI}, // equator + {0.0, -M_PI_2}, // equator + {0.0, 0.0}, // equator + {0.0, M_PI_2}, // equator + // {0.0, M_PI}, // equator don't need last one here. last + // point + // unnecessary, but doesn't harm anything. + }; + + helper(verts, ARRAY_SIZE(verts), 2 * M_PI); } } From 61702f13dab6573b3a0dc2a18bc29041c2993db1 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Sat, 29 Nov 2025 21:31:32 -0500 Subject: [PATCH 23/62] more tests --- src/apps/testapps/testGeoLoopArea.c | 194 +++++++++++++--------------- 1 file changed, 93 insertions(+), 101 deletions(-) diff --git a/src/apps/testapps/testGeoLoopArea.c b/src/apps/testapps/testGeoLoopArea.c index 985a48b324..a6663bd0a5 100644 --- a/src/apps/testapps/testGeoLoopArea.c +++ b/src/apps/testapps/testGeoLoopArea.c @@ -1,5 +1,5 @@ /* - * Copyright 2020-2021 Uber Technologies, Inc. + * Copyright 2025 Uber Technologies, Inc. * * Licensed under the Apache License, Version 2.0 (the "License"); * you may not use this file except in compliance with the License. @@ -15,9 +15,9 @@ */ /** @file - * @brief tests H3 cell area functions on a few specific cases + * @brief tests H3 GeoLoop area calculation * - * usage: `testH3CellArea` + * usage: `testGeoLoopArea` */ #include @@ -26,142 +26,134 @@ #include "test.h" #include "utility.h" -// GeoMultiPolygon alloc_global_mpoly() { -// const int numPolygons = 8; -// const int numVerts = 3; -// const LatLng verts[numPolygons][numVerts] = { -// {{M_PI_2, 0.0}, {0.0, 0.0}, {0.0, M_PI_2}}, -// {{M_PI_2, 0.0}, {0.0, M_PI_2}, {0.0, M_PI}}, -// {{M_PI_2, 0.0}, {0.0, M_PI}, {0.0, -M_PI_2}}, -// {{M_PI_2, 0.0}, {0.0, -M_PI_2}, {0.0, 0.0}}, -// {{-M_PI_2, 0.0}, {0.0, 0.0}, {0.0, -M_PI_2}}, -// {{-M_PI_2, 0.0}, {0.0, -M_PI_2}, {0.0, -M_PI}}, -// {{-M_PI_2, 0.0}, {0.0, -M_PI}, {0.0, M_PI_2}}, -// {{-M_PI_2, 0.0}, {0.0, M_PI_2}, {0.0, 0.0}}, -// }; - -// GeoMultiPolygon mpoly = { -// .numPolygons = numPolygons, -// .polygons = H3_MEMORY(malloc)(sizeof(GeoPolygon) * numPolygons), -// }; - -// for (int i = 0; i < numPolygons; i++) { -// GeoPolygon *poly = &mpoly.polygons[i]; -// poly->numHoles = 0; -// poly->holes = NULL; -// poly->geoloop.numVerts = numVerts; -// poly->geoloop.verts = H3_MEMORY(malloc)(sizeof(LatLng) * numVerts); - -// for (int j = 0; j < numVerts; j++) { -// poly->geoloop.verts[j] = verts[i][j]; -// } -// } - -// qsort(mpoly.polygons, numPolygons, sizeof(GeoPolygon), cmp_poly_area); - -// return mpoly; -// } - -#define TOL 1e-14 - -static void helper(LatLng *verts, int numVerts, double target_area) { +static void _compareArea(LatLng *verts, int numVerts, double target_area) { + double tol = 1e-14; GeoLoop loop = {.verts = verts, .numVerts = numVerts}; double out; t_assertSuccess(H3_EXPORT(geoLoopArea)(loop, &out)); - - t_assert(fabs(out - target_area) < TOL, "area should match"); + t_assert(fabs(out - target_area) < tol, "area should match"); } SUITE(geoLoopArea) { TEST(triangle) { // GeoLoop representing a triangle covering 1/8 of the globe, with // points ordered according to right-hand rule (counter-clockwise). - // Expected area: pi/2 - LatLng verts[] = {{M_PI_2, 0.0}, {0.0, 0.0}, {0.0, M_PI_2}}; - GeoLoop loop = {.verts = verts, .numVerts = ARRAY_SIZE(verts)}; - - double out; - t_assertSuccess(H3_EXPORT(geoLoopArea)(loop, &out)); + LatLng verts[] = { + {M_PI_2, 0.0}, + {0.0, 0.0}, + {0.0, M_PI_2}, + }; - t_assert(fabs(out - M_PI / 2.0) < 1e-14, "area should match"); + _compareArea(verts, ARRAY_SIZE(verts), M_PI / 2); } + TEST(reverse_triangle) { // Reversed the order of the points, so they are in clockwise order. // This GeoLoop represents the whole globe minus the triangle above. - // Expected area: 7*pi/2 - LatLng verts[] = {{0.0, M_PI_2}, {0.0, 0.0}, {M_PI_2, 0.0}}; - GeoLoop loop = {.verts = verts, .numVerts = ARRAY_SIZE(verts)}; - - double out; - t_assertSuccess(H3_EXPORT(geoLoopArea)(loop, &out)); + LatLng verts[] = { + {0.0, M_PI_2}, + {0.0, 0.0}, + {M_PI_2, 0.0}, + }; - t_assert(fabs(out - 7 * M_PI / 2.0) < 1e-14, "area should match"); + _compareArea(verts, ARRAY_SIZE(verts), 7 * M_PI / 2); } TEST(slice) { // 1/4 slice of the globe, from north to south pole - // Expected area: pi - LatLng verts[] = {{M_PI_2, 0.0}, - {0.0, 0.0}, - {-M_PI_2, 0.0}, - {0.0, M_PI_2}, - {M_PI_2, 0.0}}; - GeoLoop loop = {.verts = verts, .numVerts = ARRAY_SIZE(verts)}; - - double out; - t_assertSuccess(H3_EXPORT(geoLoopArea)(loop, &out)); - - t_assert(fabs(out - M_PI) < 1e-14, "area should match"); + LatLng verts[] = { + {M_PI_2, 0.0}, + {0.0, 0.0}, + {-M_PI_2, 0.0}, + {0.0, M_PI_2}, + }; + + _compareArea(verts, ARRAY_SIZE(verts), M_PI); } TEST(reverse_slice) { // 3/4 slice of the globe, from north to south pole, formed by reversing // order of points from example above. - // Expected area: 3*pi - LatLng verts[] = {{M_PI_2, 0.0}, - {0.0, M_PI_2}, - {-M_PI_2, 0.0}, - {0.0, 0.0}, - {M_PI_2, 0.0}}; - GeoLoop loop = {.verts = verts, .numVerts = ARRAY_SIZE(verts)}; - - double out; - t_assertSuccess(H3_EXPORT(geoLoopArea)(loop, &out)); - - t_assert(fabs(out - 3 * M_PI) < 1e-14, "area should match"); + LatLng verts[] = { + {M_PI_2, 0.0}, + {0.0, M_PI_2}, + {-M_PI_2, 0.0}, + {0.0, 0.0}, + }; + + _compareArea(verts, ARRAY_SIZE(verts), 3 * M_PI); } TEST(hemisphere_east) { // one hemisphere of globe. western or eastern? - // Expected area: 2*pi LatLng verts[] = { - {M_PI_2, 0.0}, // north pole - {0.0, -M_PI_2}, // equator - {-M_PI_2, 0.0}, // south pole - {0.0, M_PI_2}, // equator + {M_PI_2, 0.0}, + {0.0, -M_PI_2}, + {-M_PI_2, 0.0}, + {0.0, M_PI_2}, }; - GeoLoop loop = {.verts = verts, .numVerts = ARRAY_SIZE(verts)}; - - double out; - t_assertSuccess(H3_EXPORT(geoLoopArea)(loop, &out)); - t_assert(fabs(out - 2 * M_PI) < 1e-14, "area should match"); + _compareArea(verts, ARRAY_SIZE(verts), 2 * M_PI); } TEST(hemisphere_north) { // one hemisphere of globe. northern - // Expected area: 2*pi LatLng verts[] = { - {0.0, -M_PI}, // equator - {0.0, -M_PI_2}, // equator - {0.0, 0.0}, // equator - {0.0, M_PI_2}, // equator - // {0.0, M_PI}, // equator don't need last one here. last - // point - // unnecessary, but doesn't harm anything. + {0.0, -M_PI}, + {0.0, -M_PI_2}, + {0.0, 0.0}, + {0.0, M_PI_2}, }; - helper(verts, ARRAY_SIZE(verts), 2 * M_PI); + _compareArea(verts, ARRAY_SIZE(verts), 2 * M_PI); + } + + TEST(percentage_slice) { + // Demonstrate that edge arcs between points in a polygon or geoloop + // should be less than 180 degrees (M_PI radians). + // + // Create a slice from north pole to equator and back to the north pole + // that sweeps out an edge arc of t * M_PI radians along the equator, + // so it should have an area of t*M_PI for t \in [0,1]. + // + // However, notice a discontinuity at t = 1 (i.e., M_PI radians or 180 + // degrees), where expected area goes to (2 + t) * M_PI for 1 < t < 2. + // + // Recall that the area in stradians of the entire globe is 4*M_PI. + for (double t = 0; t <= 1.2; t += 0.01) { + LatLng verts[] = { + {M_PI_2, 0.0}, + {0.0, -M_PI_2}, + {0.0, t * M_PI - M_PI_2}, + }; + GeoLoop loop = {.verts = verts, .numVerts = ARRAY_SIZE(verts)}; + + double out; + t_assertSuccess(H3_EXPORT(geoLoopArea)(loop, &out)); + + double tol = 1e-13; + if (t < 0.99) { + // When t < 1, the largest angle in the triangle is less than + // 180 degrees + t_assert(fabs(out - t * M_PI) <= tol, "expected area"); + } else if (t > 1.01) { + // Discontinuity at t == 1. For t > 1, the triangle "flips", + // because the shortest geodesic path is on the other side of + // the globe. The triangle is now oriented in clockwise order, + // and the area computed is the area *outside* of the triangle, + // which starts at 3*pi. + t_assert(fabs(out - (2 + t) * M_PI) <= tol, "expected area"); + } + // Note that we avoid testing t == 1, since the triangle + // isn't well defined because there are many possible geodesic + // shortest paths when consecutive points are antipodal + // (180 degrees apart). + + // Also note that a large polygon with t > 1 is *still* + // represntable and we can compute its area accurately. + // We just need to add intermediate vertices so that + // no edge arc is greater than 180 degrees. + } } } From 4b33b0278a9cc0496275842f4e55176a7e0be5ac Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Sat, 29 Nov 2025 21:57:57 -0500 Subject: [PATCH 24/62] comments --- src/apps/testapps/testGeoLoopArea.c | 110 ++++++++++++++++++---------- 1 file changed, 70 insertions(+), 40 deletions(-) diff --git a/src/apps/testapps/testGeoLoopArea.c b/src/apps/testapps/testGeoLoopArea.c index a6663bd0a5..61b52d09c5 100644 --- a/src/apps/testapps/testGeoLoopArea.c +++ b/src/apps/testapps/testGeoLoopArea.c @@ -36,9 +36,18 @@ static void _compareArea(LatLng *verts, int numVerts, double target_area) { } SUITE(geoLoopArea) { - TEST(triangle) { - // GeoLoop representing a triangle covering 1/8 of the globe, with - // points ordered according to right-hand rule (counter-clockwise). + TEST(triangle_basic) { + /* + GeoLoop representing a triangle covering 1/8 of the globe, with + points ordered according to right-hand rule (counter-clockwise). + + The triangle starts at the north pole, moves down 90 degrees to the + equator, and then sweeps out 90 degrees along the equator + before returning to the north pole. + + The globe has an area of 4*pi radians, so this 1/8 triangle piece of + the globe should have area pi/2. + */ LatLng verts[] = { {M_PI_2, 0.0}, {0.0, 0.0}, @@ -48,9 +57,14 @@ SUITE(geoLoopArea) { _compareArea(verts, ARRAY_SIZE(verts), M_PI / 2); } - TEST(reverse_triangle) { - // Reversed the order of the points, so they are in clockwise order. - // This GeoLoop represents the whole globe minus the triangle above. + TEST(triangle_reversed) { + /* + Reverse the order of the points in the triangle from the previous test, + so that they are in clockwise order. + + Since the points are in clockwise order, GeoLoop represents the whole + globe minus the triangle above. + */ LatLng verts[] = { {0.0, M_PI_2}, {0.0, 0.0}, @@ -61,7 +75,11 @@ SUITE(geoLoopArea) { } TEST(slice) { - // 1/4 slice of the globe, from north to south pole + /* + Stitch two 1/8 triangles together, sharing an edge along the equator + to create a 1/4 slice of the globe, with vertices at the north and south + pole. + */ LatLng verts[] = { {M_PI_2, 0.0}, {0.0, 0.0}, @@ -72,9 +90,11 @@ SUITE(geoLoopArea) { _compareArea(verts, ARRAY_SIZE(verts), M_PI); } - TEST(reverse_slice) { - // 3/4 slice of the globe, from north to south pole, formed by reversing - // order of points from example above. + TEST(slice_reversed) { + /* + 3/4 slice of the globe, from north to south pole, formed by reversing + order of points from example above. + */ LatLng verts[] = { {M_PI_2, 0.0}, {0.0, M_PI_2}, @@ -86,19 +106,23 @@ SUITE(geoLoopArea) { } TEST(hemisphere_east) { - // one hemisphere of globe. western or eastern? + /* + Stitch 4 1/8 triangles together to cover the eastern hemisphere. + */ LatLng verts[] = { {M_PI_2, 0.0}, - {0.0, -M_PI_2}, + {0.0, 0}, {-M_PI_2, 0.0}, - {0.0, M_PI_2}, + {0.0, M_PI}, }; _compareArea(verts, ARRAY_SIZE(verts), 2 * M_PI); } TEST(hemisphere_north) { - // one hemisphere of globe. northern + /* + Stitch 4 1/8 triangles together to cover the northern hemisphere. + */ LatLng verts[] = { {0.0, -M_PI}, {0.0, -M_PI_2}, @@ -109,18 +133,20 @@ SUITE(geoLoopArea) { _compareArea(verts, ARRAY_SIZE(verts), 2 * M_PI); } - TEST(percentage_slice) { - // Demonstrate that edge arcs between points in a polygon or geoloop - // should be less than 180 degrees (M_PI radians). - // - // Create a slice from north pole to equator and back to the north pole - // that sweeps out an edge arc of t * M_PI radians along the equator, - // so it should have an area of t*M_PI for t \in [0,1]. - // - // However, notice a discontinuity at t = 1 (i.e., M_PI radians or 180 - // degrees), where expected area goes to (2 + t) * M_PI for 1 < t < 2. - // - // Recall that the area in stradians of the entire globe is 4*M_PI. + TEST(percentageSlice) { + /* + Demonstrate that edge arcs between points in a polygon or geoloop + should be less than 180 degrees (M_PI radians). + + Create a triangle from north pole to equator and back to the north pole + that sweeps out an edge arc of t * M_PI radians along the equator, + so it should have an area of t*M_PI for t \in [0,1]. + + However, there is a discontinuity at t = 1 (i.e., M_PI radians or 180 + degrees), where expected area goes to (2 + t) * M_PI for 1 < t < 2. + + Recall that the area in stradians of the entire globe is 4*M_PI. + */ for (double t = 0; t <= 1.2; t += 0.01) { LatLng verts[] = { {M_PI_2, 0.0}, @@ -138,22 +164,26 @@ SUITE(geoLoopArea) { // 180 degrees t_assert(fabs(out - t * M_PI) <= tol, "expected area"); } else if (t > 1.01) { - // Discontinuity at t == 1. For t > 1, the triangle "flips", - // because the shortest geodesic path is on the other side of - // the globe. The triangle is now oriented in clockwise order, - // and the area computed is the area *outside* of the triangle, - // which starts at 3*pi. + /* + Discontinuity at t == 1. For t > 1, the triangle "flips", + because the shortest geodesic path is on the other side of + the globe. The triangle is now oriented in clockwise order, + and the area computed is the area *outside* of the triangle, + which starts at 3*pi. + */ t_assert(fabs(out - (2 + t) * M_PI) <= tol, "expected area"); } - // Note that we avoid testing t == 1, since the triangle - // isn't well defined because there are many possible geodesic - // shortest paths when consecutive points are antipodal - // (180 degrees apart). - - // Also note that a large polygon with t > 1 is *still* - // represntable and we can compute its area accurately. - // We just need to add intermediate vertices so that - // no edge arc is greater than 180 degrees. + /* + Note that we avoid testing t == 1, since the triangle + isn't well defined because there are many possible geodesic + shortest paths when consecutive points are antipodal + (180 degrees apart). + + Also note that a large polygon with t > 1 is *still* + representable and we can compute its area accurately; + we just need to add intermediate vertices so that + no edge arc is greater than 180 degrees. + */ } } } From 26eb14376f1916f65d9996774168b66d1f732798 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Sat, 29 Nov 2025 22:45:06 -0500 Subject: [PATCH 25/62] docstring --- src/h3lib/lib/area.c | 38 +++++++++++++++++++++++++------------- 1 file changed, 25 insertions(+), 13 deletions(-) diff --git a/src/h3lib/lib/area.c b/src/h3lib/lib/area.c index bdbbc0f351..a1c00a8526 100644 --- a/src/h3lib/lib/area.c +++ b/src/h3lib/lib/area.c @@ -11,8 +11,6 @@ #include "constants.h" #include "h3api.h" -// TOOD: demonstrate the area alg works for global polygons. - static inline double cagnoli(LatLng x, LatLng y) { // https://github.com/d3/d3-geo/blob/8c53a90ae70c94bace73ecb02f2c792c649c86ba/src/area.js#L51-L70 x.lat = x.lat / 2.0 + M_PI / 4.0; @@ -28,6 +26,30 @@ static inline double cagnoli(LatLng x, LatLng y) { return -2.0 * atan2(sa * sd, sa * cd + ca); } +/** + * Area in radians^2 enclosed by vertices in GeoLoop. + * + * The GeoLoop should represent a simple curve with no self-intersections. + * Vertices should be ordered according to the "right hand rule". + * That is, if you are looking from outer space at a spherical polygon + * on the surface of the earth whose interior is contained within a hemisphere, + * then the vertices should be ordered counter-clockwise. The interior of the + * polygon is to the left of a person walking along the boundary of the polygon + * in the counter-clockwise direction. + * + * Note that GeoLoops do not need to repeat the first vertex at the end of the + * array to close the loop; this is done automatically. + * + * The area of the entire globe is 4*pi radians^2. If, for example, you have a + * small GeoLoop with area `a << 4*pi` and then reverse the order of the + * vertices, you produce GeoLoop with area `4*pi - a`, since by the right hand + * rule, the interior of the new loop's interior is the majority of the globe, + * or "everything except the original polygon". + * + * @param loop GeoLoop of boundary vertices in counter-clockwise order + * @param out loop area in radians^2, in interval [0, 4*pi] + * @return E_SUCCESS on success, or an error code otherwise + */ H3Error H3_EXPORT(geoLoopArea)(GeoLoop loop, double *out) { Adder adder = {0.0, 0.0}; @@ -47,11 +69,7 @@ H3Error H3_EXPORT(geoLoopArea)(GeoLoop loop, double *out) { /** * Area of H3 cell in radians^2. * - * The area is calculated by breaking the cell into spherical triangles and - * summing up their areas. Note that some H3 cells (hexagons and pentagons) - * are irregular, and have more than 6 or 5 sides. - * - * todo: optimize the computation by re-using the edges shared between triangles + * Uses `geoLoopArea` to compute cell area. * * @param cell H3 cell * @param out cell area in radians^2 @@ -72,9 +90,3 @@ H3Error H3_EXPORT(cellAreaRads2)(H3Index cell, double *out) { return E_SUCCESS; } - -void H3_EXPORT(destroyGeoLoop)(GeoLoop *loop) { - H3_MEMORY(free)(loop->verts); - loop->numVerts = 0; - loop->verts = NULL; -} From 06b2a6b8b6d9348311483d296b5af740d584623a Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Sat, 29 Nov 2025 23:08:53 -0500 Subject: [PATCH 26/62] better docstrings --- src/h3lib/lib/area.c | 27 ++++++++++++++++++++++----- 1 file changed, 22 insertions(+), 5 deletions(-) diff --git a/src/h3lib/lib/area.c b/src/h3lib/lib/area.c index a1c00a8526..db640bc8a8 100644 --- a/src/h3lib/lib/area.c +++ b/src/h3lib/lib/area.c @@ -11,8 +11,11 @@ #include "constants.h" #include "h3api.h" + +// Cagnoli contribution for edge arc x to y, following d3-geo’s +// area implementation: +// https://github.com/d3/d3-geo/blob/8c53a90ae70c94bace73ecb02f2c792c649c86ba/src/area.js#L51-L70 static inline double cagnoli(LatLng x, LatLng y) { - // https://github.com/d3/d3-geo/blob/8c53a90ae70c94bace73ecb02f2c792c649c86ba/src/area.js#L51-L70 x.lat = x.lat / 2.0 + M_PI / 4.0; y.lat = y.lat / 2.0 + M_PI / 4.0; @@ -31,20 +34,30 @@ static inline double cagnoli(LatLng x, LatLng y) { * * The GeoLoop should represent a simple curve with no self-intersections. * Vertices should be ordered according to the "right hand rule". - * That is, if you are looking from outer space at a spherical polygon + * That is, if you are looking from outer space at a spherical polygon loop * on the surface of the earth whose interior is contained within a hemisphere, * then the vertices should be ordered counter-clockwise. The interior of the - * polygon is to the left of a person walking along the boundary of the polygon + * loop is to the left of a person walking along the boundary of the polygon * in the counter-clockwise direction. * * Note that GeoLoops do not need to repeat the first vertex at the end of the * array to close the loop; this is done automatically. * + * The edge arcs between adjacent vertices are assumed to be the shortest + * geodesic path between them; that is, all arcs are interpreted to be less + * than 180 degrees or pi radians. + * Avoid arcs that are exactly pi (i.e, two antipodal vertices). + * "Large" polygon loops (e.g., that cannot be contained in a hemisphere) can + * still be constructed by using intermediate vertices with arcs less than + * 180 degrees, and the loop area will still be computed correctly. + * * The area of the entire globe is 4*pi radians^2. If, for example, you have a * small GeoLoop with area `a << 4*pi` and then reverse the order of the - * vertices, you produce GeoLoop with area `4*pi - a`, since by the right hand - * rule, the interior of the new loop's interior is the majority of the globe, + * vertices, you produce GeoLoop with area `4*pi - a`, since, by the right hand + * rule, the new loop's interior is the majority of the globe, * or "everything except the original polygon". + * Note that the area enclosed by the loop is determined by the vertex order; + * this function does **not** return `min(a, 4*pi - a)`. * * @param loop GeoLoop of boundary vertices in counter-clockwise order * @param out loop area in radians^2, in interval [0, 4*pi] @@ -58,6 +71,10 @@ H3Error H3_EXPORT(geoLoopArea)(GeoLoop loop, double *out) { kadd(&adder, cagnoli(loop.verts[i], loop.verts[j])); } + // The Cagnoli sum above yields a signed area, with the sign switching + // with the orientation of the vertices. Since we want our area to always be + // positive, we normalize into [0, 4*pi] by adding 4*pi when the signed + // area is negative. if (adder.sum < 0.0) { kadd(&adder, 4.0 * M_PI); } From f55e1b14366a6fbc085e3cffb6f20f2b68fb34e1 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Sat, 29 Nov 2025 23:10:54 -0500 Subject: [PATCH 27/62] Adder note --- src/h3lib/lib/area.c | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/h3lib/lib/area.c b/src/h3lib/lib/area.c index db640bc8a8..ac8f82a3d4 100644 --- a/src/h3lib/lib/area.c +++ b/src/h3lib/lib/area.c @@ -11,7 +11,6 @@ #include "constants.h" #include "h3api.h" - // Cagnoli contribution for edge arc x to y, following d3-geo’s // area implementation: // https://github.com/d3/d3-geo/blob/8c53a90ae70c94bace73ecb02f2c792c649c86ba/src/area.js#L51-L70 @@ -64,6 +63,8 @@ static inline double cagnoli(LatLng x, LatLng y) { * @return E_SUCCESS on success, or an error code otherwise */ H3Error H3_EXPORT(geoLoopArea)(GeoLoop loop, double *out) { + // Use `Adder` to improve numerical accuracy of the sum of many Cagnoli + // terms Adder adder = {0.0, 0.0}; for (int i = 0; i < loop.numVerts; i++) { From 0c939fce58de248d6424c1197effcb3a8d636c03 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Sat, 29 Nov 2025 23:20:11 -0500 Subject: [PATCH 28/62] don't need these imports --- justfile | 4 ++-- src/apps/testapps/testGeoLoopArea.c | 2 +- src/apps/testapps/testH3CellArea.c | 2 +- src/h3lib/lib/area.c | 7 ------- 4 files changed, 4 insertions(+), 11 deletions(-) diff --git a/justfile b/justfile index 8df37001ba..a9919cc8ce 100644 --- a/justfile +++ b/justfile @@ -26,8 +26,8 @@ test: build # ./build/bin/testDirectedEdge # ./build/bin/testArea # just test-slow - # just test-fast - ./build/bin/testGeoLoopArea + just test-fast + # ./build/bin/testGeoLoopArea time: time ./build/bin/testArea diff --git a/src/apps/testapps/testGeoLoopArea.c b/src/apps/testapps/testGeoLoopArea.c index 61b52d09c5..ffac0a505f 100644 --- a/src/apps/testapps/testGeoLoopArea.c +++ b/src/apps/testapps/testGeoLoopArea.c @@ -22,7 +22,7 @@ #include -#include "h3Index.h" +#include "h3api.h" #include "test.h" #include "utility.h" diff --git a/src/apps/testapps/testH3CellArea.c b/src/apps/testapps/testH3CellArea.c index 060e14d169..4189dfc7de 100644 --- a/src/apps/testapps/testH3CellArea.c +++ b/src/apps/testapps/testH3CellArea.c @@ -22,7 +22,7 @@ #include -#include "h3Index.h" +#include "h3api.h" #include "test.h" static const double areasKm2[] = { diff --git a/src/h3lib/lib/area.c b/src/h3lib/lib/area.c index ac8f82a3d4..ef1bb845c1 100644 --- a/src/h3lib/lib/area.c +++ b/src/h3lib/lib/area.c @@ -1,13 +1,6 @@ -#include #include -#include -#include -#include -#include -#include #include "adder.h" -#include "alloc.h" #include "constants.h" #include "h3api.h" From a292924a53d230ac6351a08c5997421292dfcbfd Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Sat, 29 Nov 2025 23:24:34 -0500 Subject: [PATCH 29/62] clean a few more imports --- src/apps/benchmarks/benchmarkArea.c | 3 --- src/h3lib/lib/area.c | 19 +++++++++++++++++++ 2 files changed, 19 insertions(+), 3 deletions(-) diff --git a/src/apps/benchmarks/benchmarkArea.c b/src/apps/benchmarks/benchmarkArea.c index a29a771680..c353ed4628 100644 --- a/src/apps/benchmarks/benchmarkArea.c +++ b/src/apps/benchmarks/benchmarkArea.c @@ -1,5 +1,3 @@ -#include -#include #include #include @@ -7,7 +5,6 @@ #include "benchmark.h" #include "h3api.h" #include "iterators.h" -#include "utility.h" static void do_res_sum(int res) { Adder adder = {0, 0}; diff --git a/src/h3lib/lib/area.c b/src/h3lib/lib/area.c index ef1bb845c1..73c5e819ac 100644 --- a/src/h3lib/lib/area.c +++ b/src/h3lib/lib/area.c @@ -1,3 +1,22 @@ +/* + * Copyright 2025 Uber Technologies, Inc. + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +/** @file area.c + * @brief Algorithms for computing areas of regions on a sphere (GeoLoop, cells, polygons, multipolygons, etc.) + */ + #include #include "adder.h" From 00f45a2f0f8ba7628d0d09853389caf204b6f188 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Sat, 29 Nov 2025 23:24:55 -0500 Subject: [PATCH 30/62] format --- src/h3lib/lib/area.c | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/h3lib/lib/area.c b/src/h3lib/lib/area.c index 73c5e819ac..88cec9565d 100644 --- a/src/h3lib/lib/area.c +++ b/src/h3lib/lib/area.c @@ -14,7 +14,8 @@ * limitations under the License. */ /** @file area.c - * @brief Algorithms for computing areas of regions on a sphere (GeoLoop, cells, polygons, multipolygons, etc.) + * @brief Algorithms for computing areas of regions on a sphere (GeoLoop, + * cells, polygons, multipolygons, etc.) */ #include From 39be84c2c6c32e4f307ca62ed95cf8f12aadc732 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Sat, 29 Nov 2025 23:30:05 -0500 Subject: [PATCH 31/62] needs constants --- src/apps/benchmarks/benchmarkArea.c | 1 + 1 file changed, 1 insertion(+) diff --git a/src/apps/benchmarks/benchmarkArea.c b/src/apps/benchmarks/benchmarkArea.c index c353ed4628..997e95e69c 100644 --- a/src/apps/benchmarks/benchmarkArea.c +++ b/src/apps/benchmarks/benchmarkArea.c @@ -5,6 +5,7 @@ #include "benchmark.h" #include "h3api.h" #include "iterators.h" +#include "constants.h" static void do_res_sum(int res) { Adder adder = {0, 0}; From 53e5be9ccbd3b7102699618f088be20cd3d3b646 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Sat, 29 Nov 2025 23:35:04 -0500 Subject: [PATCH 32/62] format --- src/apps/benchmarks/benchmarkArea.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/apps/benchmarks/benchmarkArea.c b/src/apps/benchmarks/benchmarkArea.c index 997e95e69c..c11f4842b1 100644 --- a/src/apps/benchmarks/benchmarkArea.c +++ b/src/apps/benchmarks/benchmarkArea.c @@ -3,9 +3,9 @@ #include "adder.h" #include "benchmark.h" +#include "constants.h" #include "h3api.h" #include "iterators.h" -#include "constants.h" static void do_res_sum(int res) { Adder adder = {0, 0}; From 139f483147413e0b1affcc60e642ec94d73b4893 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Sun, 30 Nov 2025 09:53:58 -0500 Subject: [PATCH 33/62] try fuzzer --- CMakeLists.txt | 2 ++ src/apps/fuzzers/README.md | 1 + src/apps/fuzzers/fuzzerGeoLoopArea.c | 39 ++++++++++++++++++++++++++++ src/h3lib/include/h3api.h.in | 1 - 4 files changed, 42 insertions(+), 1 deletion(-) create mode 100644 src/apps/fuzzers/fuzzerGeoLoopArea.c diff --git a/CMakeLists.txt b/CMakeLists.txt index dcc7069a46..d2bea788f9 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -310,6 +310,7 @@ set(OTHER_SOURCE_FILES src/apps/fuzzers/fuzzerCellToChildPos.c src/apps/fuzzers/fuzzerInternalAlgos.c src/apps/fuzzers/fuzzerInternalCoordIjk.c + src/apps/fuzzers/fuzzerGeoLoopArea.c src/apps/benchmarks/benchmarkPolygonToCells.c src/apps/benchmarks/benchmarkPolygonToCellsExperimental.c src/apps/benchmarks/benchmarkPolygon.c @@ -639,6 +640,7 @@ if(BUILD_FUZZERS) add_h3_fuzzer(fuzzerPolygonToCellsExperimentalNoHoles src/apps/fuzzers/fuzzerPolygonToCellsExperimentalNoHoles.c) add_h3_fuzzer(fuzzerCellToChildPos src/apps/fuzzers/fuzzerCellToChildPos.c) + add_h3_fuzzer(fuzzerGeoLoopArea src/apps/fuzzers/fuzzerGeoLoopArea.c) if(ENABLE_REQUIRES_ALL_SYMBOLS) add_h3_fuzzer(fuzzerInternalAlgos src/apps/fuzzers/fuzzerInternalAlgos.c) diff --git a/src/apps/fuzzers/README.md b/src/apps/fuzzers/README.md index 31135f0519..5b2a931d9f 100644 --- a/src/apps/fuzzers/README.md +++ b/src/apps/fuzzers/README.md @@ -64,6 +64,7 @@ The public API of H3 is covered in the following fuzzers: | originToDirectedEdges | [fuzzerDirectedEdge](./fuzzerDirectedEdge.c) | polygonToCells | [fuzzerPoylgonToCells](./fuzzerPolygonToCells.c) | polygonToCellsExperimental | [fuzzerPoylgonToCellsExperimental](./fuzzerPolygonToCellsExperimental.c) [fuzzerPoylgonToCellsExperimentalNoHoles](./fuzzerPolygonToCellsExperimentalNoHoles.c) +| geoLoopArea | [fuzzerGeoLoopArea](./fuzzerGeoLoopArea.c) | radsToDegs | Trivial | stringToH3 | [fuzzerIndexIO](./fuzzerIndexIO.c) | uncompactCells | [fuzzerCompact](./fuzzerCompact.c) diff --git a/src/apps/fuzzers/fuzzerGeoLoopArea.c b/src/apps/fuzzers/fuzzerGeoLoopArea.c new file mode 100644 index 0000000000..a5510d1519 --- /dev/null +++ b/src/apps/fuzzers/fuzzerGeoLoopArea.c @@ -0,0 +1,39 @@ +/* + * Copyright 2025 Uber Technologies, Inc. + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +/** @file + * @brief Fuzzer program for geoLoopArea + */ + +#include "aflHarness.h" +#include "h3api.h" +#include "utility.h" + +int LLVMFuzzerTestOneInput(const uint8_t *data, size_t size) { + // Interpret the input buffer as an array of LatLng vertices. + int numVerts = (int)(size / sizeof(LatLng)); + // Need at least 3 vertices to form a polygon. + if (numVerts < 3) { + return 0; + } + + GeoLoop loop = {.numVerts = numVerts, .verts = (LatLng *)data}; + double area; + H3_EXPORT(geoLoopArea)(loop, &area); + + return 0; +} + +AFL_HARNESS_MAIN(sizeof(LatLng) * 1024); diff --git a/src/h3lib/include/h3api.h.in b/src/h3lib/include/h3api.h.in index 042b778471..8fe85cbcb0 100644 --- a/src/h3lib/include/h3api.h.in +++ b/src/h3lib/include/h3api.h.in @@ -348,7 +348,6 @@ DECLSPEC H3Error H3_EXPORT(cellsToLinkedMultiPolygon)(const H3Index *h3Set, DECLSPEC void H3_EXPORT(destroyLinkedMultiPolygon)(LinkedGeoPolygon *polygon); /** @} */ -DECLSPEC void H3_EXPORT(destroyGeoLoop)(GeoLoop *loop); DECLSPEC H3Error H3_EXPORT(geoLoopArea)(GeoLoop loop, double *out); /** @defgroup degsToRads degsToRads From afa64f166931325f9112c614fb24f62e05144df0 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Sun, 30 Nov 2025 10:24:08 -0500 Subject: [PATCH 34/62] .mdx --- website/docs/api/regions.mdx | 90 ++++++++++++++++++++++++++++++++++++ 1 file changed, 90 insertions(+) diff --git a/website/docs/api/regions.mdx b/website/docs/api/regions.mdx index 048d62935e..0d8ede3871 100644 --- a/website/docs/api/regions.mdx +++ b/website/docs/api/regions.mdx @@ -647,4 +647,94 @@ This function exists for memory management and is not exposed. +## geoLoopArea + +Compute the spherical area, in square radians, that is enclosed by a closed +`GeoLoop`. The loop must be a simple polygon traced counter-clockwise as viewed +from above the Earth, and each edge is interpreted as the shortest geodesic +between consecutive vertices (i.e., no side spans more than π radians). Areas +are normalized to the [0, 4π] interval, so reversing the vertex order yields +the complementary area of the loop. + + + + +```c +H3Error geoLoopArea(GeoLoop loop, double *out); +``` + +The caller supplies the `GeoLoop` memory (for example, by filling `loop.verts` +with previously allocated `LatLng` vertices) and receives the enclosed area in +square radians in `out`. Returns 0 (`E_SUCCESS`) on success. + + + + +:::note + +This function is not exposed in the Java bindings. + +::: + + + + +:::note + +This function is not exposed in the JavaScript bindings. + +::: + + + + +:::note + +This function is not exposed in the Python bindings. + +TODO: direct users to create a polygon with no holes. + +::: + + + + +:::note + +This function is not exposed in the Go bindings. + +::: + + + + +:::note + +This function is not exposed in the DuckDB bindings. + +::: + + + + +:::note + +This function is not exposed in the CLI bindings. + +::: + + + From f855776e60da80714c673e9a316d29e4697d966f Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Sun, 30 Nov 2025 10:27:30 -0500 Subject: [PATCH 35/62] docs in h3api.h.in --- src/h3lib/include/h3api.h.in | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/src/h3lib/include/h3api.h.in b/src/h3lib/include/h3api.h.in index 8fe85cbcb0..19837bcb42 100644 --- a/src/h3lib/include/h3api.h.in +++ b/src/h3lib/include/h3api.h.in @@ -348,7 +348,21 @@ DECLSPEC H3Error H3_EXPORT(cellsToLinkedMultiPolygon)(const H3Index *h3Set, DECLSPEC void H3_EXPORT(destroyLinkedMultiPolygon)(LinkedGeoPolygon *polygon); /** @} */ +/** @defgroup geoLoopArea geoLoopArea + * Functions for geoLoopArea + * @{ + */ +/** + * @brief Area in radians^2 enclosed by vertices in a GeoLoop. + * + * Vertices should define a simple polygon ordered counter-clockwise according + * to the right-hand rule. Edge arcs follow the shortest geodesic path (less + * than pi radians), and the first vertex does not need to be repeated to close + * the loop. The returned area depends on vertex orientation and is in the + * interval [0, 4*pi]. + */ DECLSPEC H3Error H3_EXPORT(geoLoopArea)(GeoLoop loop, double *out); +/** @} */ /** @defgroup degsToRads degsToRads * Functions for degsToRads From e48b735105053541eff39352cdc3c08a9e9a6d35 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Sun, 30 Nov 2025 12:05:37 -0500 Subject: [PATCH 36/62] rads2 --- src/apps/fuzzers/fuzzerGeoLoopArea.c | 4 ++-- src/apps/testapps/testGeoLoopArea.c | 4 ++-- src/h3lib/include/h3api.h.in | 6 +++--- src/h3lib/lib/area.c | 6 +++--- website/docs/api/regions.mdx | 4 ++-- 5 files changed, 12 insertions(+), 12 deletions(-) diff --git a/src/apps/fuzzers/fuzzerGeoLoopArea.c b/src/apps/fuzzers/fuzzerGeoLoopArea.c index a5510d1519..88b96e1b7b 100644 --- a/src/apps/fuzzers/fuzzerGeoLoopArea.c +++ b/src/apps/fuzzers/fuzzerGeoLoopArea.c @@ -14,7 +14,7 @@ * limitations under the License. */ /** @file - * @brief Fuzzer program for geoLoopArea + * @brief Fuzzer program for geoLoopAreaRads2 */ #include "aflHarness.h" @@ -31,7 +31,7 @@ int LLVMFuzzerTestOneInput(const uint8_t *data, size_t size) { GeoLoop loop = {.numVerts = numVerts, .verts = (LatLng *)data}; double area; - H3_EXPORT(geoLoopArea)(loop, &area); + H3_EXPORT(geoLoopAreaRads2)(loop, &area); return 0; } diff --git a/src/apps/testapps/testGeoLoopArea.c b/src/apps/testapps/testGeoLoopArea.c index ffac0a505f..f043cc3e7d 100644 --- a/src/apps/testapps/testGeoLoopArea.c +++ b/src/apps/testapps/testGeoLoopArea.c @@ -31,7 +31,7 @@ static void _compareArea(LatLng *verts, int numVerts, double target_area) { GeoLoop loop = {.verts = verts, .numVerts = numVerts}; double out; - t_assertSuccess(H3_EXPORT(geoLoopArea)(loop, &out)); + t_assertSuccess(H3_EXPORT(geoLoopAreaRads2)(loop, &out)); t_assert(fabs(out - target_area) < tol, "area should match"); } @@ -156,7 +156,7 @@ SUITE(geoLoopArea) { GeoLoop loop = {.verts = verts, .numVerts = ARRAY_SIZE(verts)}; double out; - t_assertSuccess(H3_EXPORT(geoLoopArea)(loop, &out)); + t_assertSuccess(H3_EXPORT(geoLoopAreaRads2)(loop, &out)); double tol = 1e-13; if (t < 0.99) { diff --git a/src/h3lib/include/h3api.h.in b/src/h3lib/include/h3api.h.in index 19837bcb42..1220414f08 100644 --- a/src/h3lib/include/h3api.h.in +++ b/src/h3lib/include/h3api.h.in @@ -348,8 +348,8 @@ DECLSPEC H3Error H3_EXPORT(cellsToLinkedMultiPolygon)(const H3Index *h3Set, DECLSPEC void H3_EXPORT(destroyLinkedMultiPolygon)(LinkedGeoPolygon *polygon); /** @} */ -/** @defgroup geoLoopArea geoLoopArea - * Functions for geoLoopArea +/** @defgroup geoLoopAreaRads2 geoLoopAreaRads2 + * Functions for geoLoopAreaRads2 * @{ */ /** @@ -361,7 +361,7 @@ DECLSPEC void H3_EXPORT(destroyLinkedMultiPolygon)(LinkedGeoPolygon *polygon); * the loop. The returned area depends on vertex orientation and is in the * interval [0, 4*pi]. */ -DECLSPEC H3Error H3_EXPORT(geoLoopArea)(GeoLoop loop, double *out); +DECLSPEC H3Error H3_EXPORT(geoLoopAreaRads2)(GeoLoop loop, double *out); /** @} */ /** @defgroup degsToRads degsToRads diff --git a/src/h3lib/lib/area.c b/src/h3lib/lib/area.c index 88cec9565d..f94b68f306 100644 --- a/src/h3lib/lib/area.c +++ b/src/h3lib/lib/area.c @@ -75,7 +75,7 @@ static inline double cagnoli(LatLng x, LatLng y) { * @param out loop area in radians^2, in interval [0, 4*pi] * @return E_SUCCESS on success, or an error code otherwise */ -H3Error H3_EXPORT(geoLoopArea)(GeoLoop loop, double *out) { +H3Error H3_EXPORT(geoLoopAreaRads2)(GeoLoop loop, double *out) { // Use `Adder` to improve numerical accuracy of the sum of many Cagnoli // terms Adder adder = {0.0, 0.0}; @@ -100,7 +100,7 @@ H3Error H3_EXPORT(geoLoopArea)(GeoLoop loop, double *out) { /** * Area of H3 cell in radians^2. * - * Uses `geoLoopArea` to compute cell area. + * Uses `geoLoopAreaRads2` to compute cell area. * * @param cell H3 cell * @param out cell area in radians^2 @@ -114,7 +114,7 @@ H3Error H3_EXPORT(cellAreaRads2)(H3Index cell, double *out) { } GeoLoop loop = {.verts = cb.verts, .numVerts = cb.numVerts}; - err = H3_EXPORT(geoLoopArea)(loop, out); + err = H3_EXPORT(geoLoopAreaRads2)(loop, out); if (err) { return err; } diff --git a/website/docs/api/regions.mdx b/website/docs/api/regions.mdx index 0d8ede3871..7ccf0c1248 100644 --- a/website/docs/api/regions.mdx +++ b/website/docs/api/regions.mdx @@ -647,7 +647,7 @@ This function exists for memory management and is not exposed. -## geoLoopArea +## geoLoopAreaRads2 Compute the spherical area, in square radians, that is enclosed by a closed `GeoLoop`. The loop must be a simple polygon traced counter-clockwise as viewed @@ -672,7 +672,7 @@ the complementary area of the loop. ```c -H3Error geoLoopArea(GeoLoop loop, double *out); +H3Error geoLoopAreaRads2(GeoLoop loop, double *out); ``` The caller supplies the `GeoLoop` memory (for example, by filling `loop.verts` From 11b277a2b97282d52393d7c40bbdd40b468d94fb Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Sun, 30 Nov 2025 12:25:47 -0500 Subject: [PATCH 37/62] python plan --- website/docs/api/regions.mdx | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/website/docs/api/regions.mdx b/website/docs/api/regions.mdx index 7ccf0c1248..57711dfcde 100644 --- a/website/docs/api/regions.mdx +++ b/website/docs/api/regions.mdx @@ -702,9 +702,13 @@ This function is not exposed in the JavaScript bindings. :::note -This function is not exposed in the Python bindings. +This function isn't exposed directly, but can be used +via `h3.h3shape_area()`: -TODO: direct users to create a polygon with no holes. +```python +poly = h3.LatLngPoly([(90,0),(0,0),(0,90)]) +h3.h3shape_area(poly, unit='km^2') +``` ::: From 226da8cadffb27f0110e2d8dcf43d537d03e04f1 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Sun, 30 Nov 2025 12:32:43 -0500 Subject: [PATCH 38/62] show, don't tell --- src/apps/testapps/testGeoLoopArea.c | 23 ++++++++++++++++++----- 1 file changed, 18 insertions(+), 5 deletions(-) diff --git a/src/apps/testapps/testGeoLoopArea.c b/src/apps/testapps/testGeoLoopArea.c index f043cc3e7d..95c3c5f4aa 100644 --- a/src/apps/testapps/testGeoLoopArea.c +++ b/src/apps/testapps/testGeoLoopArea.c @@ -178,12 +178,25 @@ SUITE(geoLoopArea) { isn't well defined because there are many possible geodesic shortest paths when consecutive points are antipodal (180 degrees apart). - - Also note that a large polygon with t > 1 is *still* - representable and we can compute its area accurately; - we just need to add intermediate vertices so that - no edge arc is greater than 180 degrees. */ } } + + TEST(percentageSlice_large) { + /* + Continuing from the test above, note that a large polygon with + t > 1 is *still* representable and we can compute its area accurately; + we just need to add intermediate vertices so that + no edge arc is greater than 180 degrees. + */ + double t = 1.2; + LatLng verts[] = { + {M_PI_2, 0.0}, + {0.0, -M_PI_2}, + {0.0, 0.0}, // Extra vertex so every angle is < 180 degrees + {0.0, t * M_PI - M_PI_2}, + }; + + _compareArea(verts, ARRAY_SIZE(verts), t * M_PI); + } } From c274bff131a6bb129194e2d75610e5d8d460e5e9 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Sun, 30 Nov 2025 12:55:40 -0500 Subject: [PATCH 39/62] benchmark clean up --- src/apps/benchmarks/benchmarkArea.c | 20 +++++++++++++++----- 1 file changed, 15 insertions(+), 5 deletions(-) diff --git a/src/apps/benchmarks/benchmarkArea.c b/src/apps/benchmarks/benchmarkArea.c index c11f4842b1..cf412ce3b5 100644 --- a/src/apps/benchmarks/benchmarkArea.c +++ b/src/apps/benchmarks/benchmarkArea.c @@ -1,4 +1,5 @@ #include +#include #include #include "adder.h" @@ -7,7 +8,7 @@ #include "h3api.h" #include "iterators.h" -static void do_res_sum(int res) { +static void doResSum(int res, bool print) { Adder adder = {0, 0}; double cellArea; IterCellsResolution iter = iterInitRes(res); @@ -19,15 +20,24 @@ static void do_res_sum(int res) { double total_area = adder.sum; double diff = fabs(total_area - 4 * M_PI); - printf("res: %d, diff: %e\n", res, diff); + if (print) { + printf("res: %d, diff: %e\n", res, diff); + } } BEGIN_BENCHMARKS(); -BENCHMARK(directedEdgeToBoundary, 1, { - int MAX_RES = 6; +static int MAX_RES = 4; + +BENCHMARK(allCellsAtRes_print, 1, { + for (int i = 0; i <= MAX_RES; i++) { + doResSum(i, true); + } +}); + +BENCHMARK(allCellsAtRes_noprint, 100, { for (int i = 0; i <= MAX_RES; i++) { - do_res_sum(i); + doResSum(i, false); } }); From 43f19a949bc1728dec3e545d5240d33a075d0746 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Sun, 30 Nov 2025 14:04:19 -0500 Subject: [PATCH 40/62] never say never again --- src/h3lib/lib/area.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/h3lib/lib/area.c b/src/h3lib/lib/area.c index f94b68f306..a62f5b68bf 100644 --- a/src/h3lib/lib/area.c +++ b/src/h3lib/lib/area.c @@ -115,7 +115,7 @@ H3Error H3_EXPORT(cellAreaRads2)(H3Index cell, double *out) { GeoLoop loop = {.verts = cb.verts, .numVerts = cb.numVerts}; err = H3_EXPORT(geoLoopAreaRads2)(loop, out); - if (err) { + if (NEVER(err)) { return err; } From fffd72399e9956da2cd117fa0fe8eb997bff8c73 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Sun, 30 Nov 2025 14:07:25 -0500 Subject: [PATCH 41/62] ugh, would have been cooler if it compiled the first time --- src/h3lib/lib/area.c | 1 + 1 file changed, 1 insertion(+) diff --git a/src/h3lib/lib/area.c b/src/h3lib/lib/area.c index a62f5b68bf..9ef3b3c629 100644 --- a/src/h3lib/lib/area.c +++ b/src/h3lib/lib/area.c @@ -22,6 +22,7 @@ #include "adder.h" #include "constants.h" +#include "h3Assert.h" #include "h3api.h" // Cagnoli contribution for edge arc x to y, following d3-geo’s From efcbea0dae9b3b11f4abdf32a5c5e293ac479e5d Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Sun, 30 Nov 2025 14:25:51 -0500 Subject: [PATCH 42/62] remove justfile --- justfile | 48 ------------------------------------------------ 1 file changed, 48 deletions(-) delete mode 100644 justfile diff --git a/justfile b/justfile deleted file mode 100644 index a9919cc8ce..0000000000 --- a/justfile +++ /dev/null @@ -1,48 +0,0 @@ -# TODO: Remove before landing - -init: purge - mkdir build - -build: - cd build; cmake -DCMAKE_BUILD_TYPE=Release ..; make - -profile: init - cd build; cmake -DCMAKE_BUILD_TYPE=RelWithDebInfo -DCMAKE_C_FLAGS='-fno-omit-frame-pointer' ..; make - # cd build; cmake -DCMAKE_BUILD_TYPE=Release ..; make - xcrun xctrace record \ - --template 'Time Profiler' \ - --output h3-prof.trace \ - --launch -- ./build/bin/benchmarkCellsToLinkedMultiPolygon - open h3-prof.trace - -purge: - rm -rf build - rm -rf *.trace - rm -rf .ipynb_checkpoints - -test: build - # ./build/bin/testH3CellAreaExhaustive - # ./build/bin/testEdgeCellsToPoly - # ./build/bin/testDirectedEdge - # ./build/bin/testArea - # just test-slow - just test-fast - # ./build/bin/testGeoLoopArea - -time: - time ./build/bin/testArea - -test-fast: build - cd build; make test-fast - -test-slow: build - cd build; make test - -bench: build - ./build/bin/benchmarkArea - # ./build/bin/benchmarkCellsToLinkedMultiPolygon - # ./build/bin/benchmarkCellsToMultiPolygon - # ./build/bin/benchmarkDirectedEdge - -fuzz: - cd build; CC=clang cmake -DENABLE_LIBFUZZER=ON ..; make fuzzers From 6e955c9541a0fbc45743fa6a1c50acb3a90c12a8 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Sun, 30 Nov 2025 14:44:00 -0500 Subject: [PATCH 43/62] lighten the h3api.h.in description, as per usual in that file --- CHANGELOG.md | 5 +++++ src/h3lib/include/h3api.h.in | 8 +------- 2 files changed, 6 insertions(+), 7 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 415bcdf886..a08417da08 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,6 +6,11 @@ The public API of this library consists of the functions declared in file [h3api.h.in](./src/h3lib/include/h3api.h.in). ## [Unreleased] +### Added +- `geoLoopArea` function (#1101) + +### Changed +- `cellAreaRads2` uses `geoLoopArea` (#1101) ## [4.4.1] - 2025-11-11 ### Fixed diff --git a/src/h3lib/include/h3api.h.in b/src/h3lib/include/h3api.h.in index 1220414f08..de589d50fa 100644 --- a/src/h3lib/include/h3api.h.in +++ b/src/h3lib/include/h3api.h.in @@ -353,13 +353,7 @@ DECLSPEC void H3_EXPORT(destroyLinkedMultiPolygon)(LinkedGeoPolygon *polygon); * @{ */ /** - * @brief Area in radians^2 enclosed by vertices in a GeoLoop. - * - * Vertices should define a simple polygon ordered counter-clockwise according - * to the right-hand rule. Edge arcs follow the shortest geodesic path (less - * than pi radians), and the first vertex does not need to be repeated to close - * the loop. The returned area depends on vertex orientation and is in the - * interval [0, 4*pi]. + * @brief Area in radians^2 enclosed by vertices of a GeoLoop. */ DECLSPEC H3Error H3_EXPORT(geoLoopAreaRads2)(GeoLoop loop, double *out); /** @} */ From dfeb6d828ca7963f7fd3cd08082c5e6c138f065e Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Sun, 30 Nov 2025 18:23:08 -0500 Subject: [PATCH 44/62] simplify adder initialization --- src/apps/benchmarks/benchmarkArea.c | 2 +- src/apps/testapps/testH3CellAreaExhaustive.c | 2 +- src/h3lib/include/adder.h | 2 +- src/h3lib/lib/area.c | 2 +- 4 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/apps/benchmarks/benchmarkArea.c b/src/apps/benchmarks/benchmarkArea.c index cf412ce3b5..fe1ffa4e29 100644 --- a/src/apps/benchmarks/benchmarkArea.c +++ b/src/apps/benchmarks/benchmarkArea.c @@ -9,7 +9,7 @@ #include "iterators.h" static void doResSum(int res, bool print) { - Adder adder = {0, 0}; + Adder adder = {0}; double cellArea; IterCellsResolution iter = iterInitRes(res); diff --git a/src/apps/testapps/testH3CellAreaExhaustive.c b/src/apps/testapps/testH3CellAreaExhaustive.c index 315eb1e91e..5a8bd30aa6 100644 --- a/src/apps/testapps/testH3CellAreaExhaustive.c +++ b/src/apps/testapps/testH3CellAreaExhaustive.c @@ -137,7 +137,7 @@ static void cell_area_assert(H3Index cell) { */ static void earth_area_test(int res, H3Error (*cell_area)(H3Index, double *), double target, double tol) { - Adder adder = {0.0, 0.0}; + Adder adder = {0}; IterCellsResolution iter = iterInitRes(res); for (; iter.h; iterStepRes(&iter)) { diff --git a/src/h3lib/include/adder.h b/src/h3lib/include/adder.h index 53e31ac5f8..57682d4a78 100644 --- a/src/h3lib/include/adder.h +++ b/src/h3lib/include/adder.h @@ -31,7 +31,7 @@ tradeoffs between Kahan, Neumaier, and `fsum`. ## Usage -Adder adder = {0.0, 0.0}; // initialize +Adder adder = {0}; // initialize kadd(&adder, x); // add x to summation out = adder.sum; // extract final sum */ diff --git a/src/h3lib/lib/area.c b/src/h3lib/lib/area.c index 9ef3b3c629..5b0b2ac3cb 100644 --- a/src/h3lib/lib/area.c +++ b/src/h3lib/lib/area.c @@ -79,7 +79,7 @@ static inline double cagnoli(LatLng x, LatLng y) { H3Error H3_EXPORT(geoLoopAreaRads2)(GeoLoop loop, double *out) { // Use `Adder` to improve numerical accuracy of the sum of many Cagnoli // terms - Adder adder = {0.0, 0.0}; + Adder adder = {0}; for (int i = 0; i < loop.numVerts; i++) { int j = (i + 1) % loop.numVerts; From 6db33a30962f327c1b94307f3869ab703623790f Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Sun, 30 Nov 2025 20:33:07 -0500 Subject: [PATCH 45/62] drop numVerts < 3 --- justfile | 53 ++++++++++++++++++++++++++++ src/apps/fuzzers/fuzzerGeoLoopArea.c | 11 +++--- 2 files changed, 57 insertions(+), 7 deletions(-) create mode 100644 justfile diff --git a/justfile b/justfile new file mode 100644 index 0000000000..cffb5a6fde --- /dev/null +++ b/justfile @@ -0,0 +1,53 @@ +# TODO: Remove before landing + +init: purge + mkdir build + +build: + cd build; cmake -DCMAKE_BUILD_TYPE=Release ..; make + +profile: init + cd build; cmake -DCMAKE_BUILD_TYPE=RelWithDebInfo -DCMAKE_C_FLAGS='-fno-omit-frame-pointer' ..; make + # cd build; cmake -DCMAKE_BUILD_TYPE=Release ..; make + xcrun xctrace record \ + --template 'Time Profiler' \ + --output h3-prof.trace \ + --launch -- ./build/bin/benchmarkCellsToLinkedMultiPolygon + open h3-prof.trace + +purge: + rm -rf build + rm -rf *.trace + rm -rf .ipynb_checkpoints + +test: build + # ./build/bin/testH3CellAreaExhaustive + # ./build/bin/testEdgeCellsToPoly + # ./build/bin/testDirectedEdge + # ./build/bin/testArea + # just test-slow + just test-fast + # ./build/bin/testGeoLoopArea + +time: + time ./build/bin/testArea + +test-fast: build + cd build; make test-fast + +test-slow: build + cd build; make test + +bench: build + ./build/bin/benchmarkArea + # ./build/bin/benchmarkCellsToLinkedMultiPolygon + # ./build/bin/benchmarkCellsToMultiPolygon + # ./build/bin/benchmarkDirectedEdge + +fuzz: + cd build; CC=clang cmake -DENABLE_LIBFUZZER=ON ..; make fuzzers + +coverage: purge + mkdir build + cd build; cmake -DCMAKE_BUILD_TYPE=Debug -DENABLE_COVERAGE=ON ..; make; make coverage + open build/coverage/index.html diff --git a/src/apps/fuzzers/fuzzerGeoLoopArea.c b/src/apps/fuzzers/fuzzerGeoLoopArea.c index 88b96e1b7b..b13b69f488 100644 --- a/src/apps/fuzzers/fuzzerGeoLoopArea.c +++ b/src/apps/fuzzers/fuzzerGeoLoopArea.c @@ -23,13 +23,10 @@ int LLVMFuzzerTestOneInput(const uint8_t *data, size_t size) { // Interpret the input buffer as an array of LatLng vertices. - int numVerts = (int)(size / sizeof(LatLng)); - // Need at least 3 vertices to form a polygon. - if (numVerts < 3) { - return 0; - } - - GeoLoop loop = {.numVerts = numVerts, .verts = (LatLng *)data}; + GeoLoop loop = { + .numVerts = (int)(size / sizeof(LatLng)), + .verts = (LatLng *)data, + }; double area; H3_EXPORT(geoLoopAreaRads2)(loop, &area); From e2c20b0d96e2c96997be7be822cfb713cff6cf19 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Sun, 30 Nov 2025 20:42:05 -0500 Subject: [PATCH 46/62] degenerate loop tests --- src/apps/testapps/testGeoLoopArea.c | 23 +++++++++++++++++++++++ 1 file changed, 23 insertions(+) diff --git a/src/apps/testapps/testGeoLoopArea.c b/src/apps/testapps/testGeoLoopArea.c index 95c3c5f4aa..5de08002a5 100644 --- a/src/apps/testapps/testGeoLoopArea.c +++ b/src/apps/testapps/testGeoLoopArea.c @@ -199,4 +199,27 @@ SUITE(geoLoopArea) { _compareArea(verts, ARRAY_SIZE(verts), t * M_PI); } + + TEST(degenerateLoop2) { + // TODO: should we raise an error when given degenerate loops? + LatLng verts[] = { + {M_PI_2, 0.0}, + {0.0, -M_PI_2}, + }; + _compareArea(verts, ARRAY_SIZE(verts), 0.0); + } + + TEST(degenerateLoop1) { + // TODO: should we raise an error when given degenerate loops? + LatLng verts[] = { + {0.0, 0.0}, + }; + _compareArea(verts, ARRAY_SIZE(verts), 0.0); + } + + TEST(degenerateLoop0) { + // TODO: should we raise an error when given degenerate loops? + LatLng verts[] = {}; + _compareArea(verts, ARRAY_SIZE(verts), 0.0); + } } From eeb3046fc694ec2bd310edf225fac47343b630df Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Sun, 30 Nov 2025 20:51:44 -0500 Subject: [PATCH 47/62] try _compareArea(NULL, 0, 0.0); --- src/apps/testapps/testGeoLoopArea.c | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/apps/testapps/testGeoLoopArea.c b/src/apps/testapps/testGeoLoopArea.c index 5de08002a5..db331c71f8 100644 --- a/src/apps/testapps/testGeoLoopArea.c +++ b/src/apps/testapps/testGeoLoopArea.c @@ -219,7 +219,6 @@ SUITE(geoLoopArea) { TEST(degenerateLoop0) { // TODO: should we raise an error when given degenerate loops? - LatLng verts[] = {}; - _compareArea(verts, ARRAY_SIZE(verts), 0.0); + _compareArea(NULL, 0, 0.0); } } From 28105c3a51c833073bf22a16266d8abf14e088b8 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Mon, 1 Dec 2025 18:20:59 -0500 Subject: [PATCH 48/62] remove docs --- website/docs/api/regions.mdx | 96 ------------------------------------ 1 file changed, 96 deletions(-) diff --git a/website/docs/api/regions.mdx b/website/docs/api/regions.mdx index 1431cc55bd..ceebf05ade 100644 --- a/website/docs/api/regions.mdx +++ b/website/docs/api/regions.mdx @@ -661,99 +661,3 @@ This function exists for memory management and is not exposed. - -## geoLoopAreaRads2 - -Compute the spherical area, in square radians, that is enclosed by a closed -`GeoLoop`. The loop must be a simple polygon traced counter-clockwise as viewed -from above the Earth, and each edge is interpreted as the shortest geodesic -between consecutive vertices (i.e., no side spans more than π radians). Areas -are normalized to the [0, 4π] interval, so reversing the vertex order yields -the complementary area of the loop. - - - - -```c -H3Error geoLoopAreaRads2(GeoLoop loop, double *out); -``` - -The caller supplies the `GeoLoop` memory (for example, by filling `loop.verts` -with previously allocated `LatLng` vertices) and receives the enclosed area in -square radians in `out`. Returns 0 (`E_SUCCESS`) on success. - - - - -:::note - -This function is not exposed in the Java bindings. - -::: - - - - -:::note - -This function is not exposed in the JavaScript bindings. - -::: - - - - -:::note - -This function isn't exposed directly, but can be used -via `h3.h3shape_area()`: - -```python -poly = h3.LatLngPoly([(90,0),(0,0),(0,90)]) -h3.h3shape_area(poly, unit='km^2') -``` - -::: - - - - -:::note - -This function is not exposed in the Go bindings. - -::: - - - - -:::note - -This function is not exposed in the DuckDB bindings. - -::: - - - - -:::note - -This function is not exposed in the CLI bindings. - -::: - - - - From a3ee94dac936611b64b2d7bbdb0a4579e2233234 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Mon, 1 Dec 2025 18:22:57 -0500 Subject: [PATCH 49/62] drop H3_EXPORT --- src/apps/fuzzers/fuzzerGeoLoopArea.c | 2 +- src/apps/testapps/testGeoLoopArea.c | 4 ++-- src/h3lib/include/h3api.h.in | 2 +- src/h3lib/lib/area.c | 4 ++-- 4 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/apps/fuzzers/fuzzerGeoLoopArea.c b/src/apps/fuzzers/fuzzerGeoLoopArea.c index b13b69f488..652c7865cf 100644 --- a/src/apps/fuzzers/fuzzerGeoLoopArea.c +++ b/src/apps/fuzzers/fuzzerGeoLoopArea.c @@ -28,7 +28,7 @@ int LLVMFuzzerTestOneInput(const uint8_t *data, size_t size) { .verts = (LatLng *)data, }; double area; - H3_EXPORT(geoLoopAreaRads2)(loop, &area); + geoLoopAreaRads2(loop, &area); return 0; } diff --git a/src/apps/testapps/testGeoLoopArea.c b/src/apps/testapps/testGeoLoopArea.c index db331c71f8..dacdfa3265 100644 --- a/src/apps/testapps/testGeoLoopArea.c +++ b/src/apps/testapps/testGeoLoopArea.c @@ -31,7 +31,7 @@ static void _compareArea(LatLng *verts, int numVerts, double target_area) { GeoLoop loop = {.verts = verts, .numVerts = numVerts}; double out; - t_assertSuccess(H3_EXPORT(geoLoopAreaRads2)(loop, &out)); + t_assertSuccess(geoLoopAreaRads2(loop, &out)); t_assert(fabs(out - target_area) < tol, "area should match"); } @@ -156,7 +156,7 @@ SUITE(geoLoopArea) { GeoLoop loop = {.verts = verts, .numVerts = ARRAY_SIZE(verts)}; double out; - t_assertSuccess(H3_EXPORT(geoLoopAreaRads2)(loop, &out)); + t_assertSuccess(geoLoopAreaRads2(loop, &out)); double tol = 1e-13; if (t < 0.99) { diff --git a/src/h3lib/include/h3api.h.in b/src/h3lib/include/h3api.h.in index c8e7df2fee..bf13ff9d1c 100644 --- a/src/h3lib/include/h3api.h.in +++ b/src/h3lib/include/h3api.h.in @@ -355,7 +355,7 @@ DECLSPEC void H3_EXPORT(destroyLinkedMultiPolygon)(LinkedGeoPolygon *polygon); /** * @brief Area in radians^2 enclosed by vertices of a GeoLoop. */ -DECLSPEC H3Error H3_EXPORT(geoLoopAreaRads2)(GeoLoop loop, double *out); +DECLSPEC H3Error geoLoopAreaRads2(GeoLoop loop, double *out); /** @} */ /** @defgroup degsToRads degsToRads diff --git a/src/h3lib/lib/area.c b/src/h3lib/lib/area.c index 5b0b2ac3cb..f1499357a5 100644 --- a/src/h3lib/lib/area.c +++ b/src/h3lib/lib/area.c @@ -76,7 +76,7 @@ static inline double cagnoli(LatLng x, LatLng y) { * @param out loop area in radians^2, in interval [0, 4*pi] * @return E_SUCCESS on success, or an error code otherwise */ -H3Error H3_EXPORT(geoLoopAreaRads2)(GeoLoop loop, double *out) { +H3Error geoLoopAreaRads2(GeoLoop loop, double *out) { // Use `Adder` to improve numerical accuracy of the sum of many Cagnoli // terms Adder adder = {0}; @@ -115,7 +115,7 @@ H3Error H3_EXPORT(cellAreaRads2)(H3Index cell, double *out) { } GeoLoop loop = {.verts = cb.verts, .numVerts = cb.numVerts}; - err = H3_EXPORT(geoLoopAreaRads2)(loop, out); + err = geoLoopAreaRads2(loop, out); if (NEVER(err)) { return err; } From 571d72ea4f9ac90e452ecd127c3d922b9d62c38d Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Mon, 1 Dec 2025 18:28:34 -0500 Subject: [PATCH 50/62] move geoLoopAreaRads2 out of public API for now --- CMakeLists.txt | 1 + src/apps/fuzzers/fuzzerGeoLoopArea.c | 1 + src/apps/testapps/testGeoLoopArea.c | 1 + src/h3lib/include/adder.h | 16 ++++++++++++++++ src/h3lib/include/area.h | 16 ++++++++++++++++ src/h3lib/include/h3api.h.in | 10 ---------- src/h3lib/lib/area.c | 2 ++ 7 files changed, 37 insertions(+), 10 deletions(-) create mode 100644 src/h3lib/include/area.h diff --git a/CMakeLists.txt b/CMakeLists.txt index d2bea788f9..2ee8421547 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -171,6 +171,7 @@ set(LIB_SOURCE_FILES src/h3lib/include/coordijk.h src/h3lib/include/algos.h src/h3lib/include/adder.h + src/h3lib/include/area.h src/h3lib/lib/h3Assert.c src/h3lib/lib/algos.c src/h3lib/lib/coordijk.c diff --git a/src/apps/fuzzers/fuzzerGeoLoopArea.c b/src/apps/fuzzers/fuzzerGeoLoopArea.c index 652c7865cf..eb39842b7a 100644 --- a/src/apps/fuzzers/fuzzerGeoLoopArea.c +++ b/src/apps/fuzzers/fuzzerGeoLoopArea.c @@ -18,6 +18,7 @@ */ #include "aflHarness.h" +#include "area.h" #include "h3api.h" #include "utility.h" diff --git a/src/apps/testapps/testGeoLoopArea.c b/src/apps/testapps/testGeoLoopArea.c index dacdfa3265..fd9ba1b230 100644 --- a/src/apps/testapps/testGeoLoopArea.c +++ b/src/apps/testapps/testGeoLoopArea.c @@ -22,6 +22,7 @@ #include +#include "area.h" #include "h3api.h" #include "test.h" #include "utility.h" diff --git a/src/h3lib/include/adder.h b/src/h3lib/include/adder.h index 57682d4a78..dfebb0bcf1 100644 --- a/src/h3lib/include/adder.h +++ b/src/h3lib/include/adder.h @@ -1,3 +1,19 @@ +/* + * Copyright 2025 Uber Technologies, Inc. + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + #ifndef ADDER_H #define ADDER_H diff --git a/src/h3lib/include/area.h b/src/h3lib/include/area.h new file mode 100644 index 0000000000..22f43aeb2d --- /dev/null +++ b/src/h3lib/include/area.h @@ -0,0 +1,16 @@ +#ifndef AREA_H +#define AREA_H + +#include "h3api.h" + +/** @defgroup geoLoopAreaRads2 geoLoopAreaRads2 + * Functions for geoLoopAreaRads2 + * @{ + */ +/** + * @brief Area in radians^2 enclosed by vertices of a GeoLoop. + */ +H3Error geoLoopAreaRads2(GeoLoop loop, double *out); +/** @} */ + +#endif diff --git a/src/h3lib/include/h3api.h.in b/src/h3lib/include/h3api.h.in index bf13ff9d1c..f237dc3582 100644 --- a/src/h3lib/include/h3api.h.in +++ b/src/h3lib/include/h3api.h.in @@ -348,16 +348,6 @@ DECLSPEC H3Error H3_EXPORT(cellsToLinkedMultiPolygon)(const H3Index *h3Set, DECLSPEC void H3_EXPORT(destroyLinkedMultiPolygon)(LinkedGeoPolygon *polygon); /** @} */ -/** @defgroup geoLoopAreaRads2 geoLoopAreaRads2 - * Functions for geoLoopAreaRads2 - * @{ - */ -/** - * @brief Area in radians^2 enclosed by vertices of a GeoLoop. - */ -DECLSPEC H3Error geoLoopAreaRads2(GeoLoop loop, double *out); -/** @} */ - /** @defgroup degsToRads degsToRads * Functions for degsToRads * @{ diff --git a/src/h3lib/lib/area.c b/src/h3lib/lib/area.c index f1499357a5..e80d25da5f 100644 --- a/src/h3lib/lib/area.c +++ b/src/h3lib/lib/area.c @@ -18,6 +18,8 @@ * cells, polygons, multipolygons, etc.) */ +#include "area.h" + #include #include "adder.h" From 7d0f56a8488ccb447d6243bac80a925432f3173d Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Mon, 1 Dec 2025 18:31:23 -0500 Subject: [PATCH 51/62] slim down benchmark --- src/apps/benchmarks/benchmarkArea.c | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/apps/benchmarks/benchmarkArea.c b/src/apps/benchmarks/benchmarkArea.c index fe1ffa4e29..193f6e66f2 100644 --- a/src/apps/benchmarks/benchmarkArea.c +++ b/src/apps/benchmarks/benchmarkArea.c @@ -27,7 +27,7 @@ static void doResSum(int res, bool print) { BEGIN_BENCHMARKS(); -static int MAX_RES = 4; +static int MAX_RES = 3; BENCHMARK(allCellsAtRes_print, 1, { for (int i = 0; i <= MAX_RES; i++) { @@ -35,7 +35,7 @@ BENCHMARK(allCellsAtRes_print, 1, { } }); -BENCHMARK(allCellsAtRes_noprint, 100, { +BENCHMARK(allCellsAtRes_noprint, 10, { for (int i = 0; i <= MAX_RES; i++) { doResSum(i, false); } From 77278a6089e1a83222c317fc78a87b0735ace5d6 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Mon, 1 Dec 2025 18:50:26 -0500 Subject: [PATCH 52/62] clear these comments, maybe? --- src/h3lib/include/area.h | 8 -------- 1 file changed, 8 deletions(-) diff --git a/src/h3lib/include/area.h b/src/h3lib/include/area.h index 22f43aeb2d..bf9235eb48 100644 --- a/src/h3lib/include/area.h +++ b/src/h3lib/include/area.h @@ -3,14 +3,6 @@ #include "h3api.h" -/** @defgroup geoLoopAreaRads2 geoLoopAreaRads2 - * Functions for geoLoopAreaRads2 - * @{ - */ -/** - * @brief Area in radians^2 enclosed by vertices of a GeoLoop. - */ H3Error geoLoopAreaRads2(GeoLoop loop, double *out); -/** @} */ #endif From fa555e33a2b063f84bcdf8347e266cd1db6ffabf Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Mon, 1 Dec 2025 19:04:54 -0500 Subject: [PATCH 53/62] try adding area.h to APP_SOURCE_FILES --- CMakeLists.txt | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 2ee8421547..f1a127a5ab 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -200,7 +200,8 @@ set(APP_SOURCE_FILES src/apps/applib/include/aflHarness.h src/apps/applib/lib/kml.c src/apps/applib/lib/utility.c - src/apps/applib/lib/args.c) + src/apps/applib/lib/args.c + src/h3lib/include/area.h) set(TEST_APP_SOURCE_FILES src/apps/applib/include/test.h src/apps/applib/lib/test.c) set(EXAMPLE_SOURCE_FILES From 64e9ad88115645e1301b7c74a0ade94cd9150742 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Mon, 1 Dec 2025 19:14:12 -0500 Subject: [PATCH 54/62] can't fail if you don't try! Also, probably do not need to fuzz if no longer in the public API --- CMakeLists.txt | 5 +--- src/apps/fuzzers/README.md | 1 - src/apps/fuzzers/fuzzerGeoLoopArea.c | 37 ---------------------------- 3 files changed, 1 insertion(+), 42 deletions(-) delete mode 100644 src/apps/fuzzers/fuzzerGeoLoopArea.c diff --git a/CMakeLists.txt b/CMakeLists.txt index f1a127a5ab..6b4c8983ec 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -200,8 +200,7 @@ set(APP_SOURCE_FILES src/apps/applib/include/aflHarness.h src/apps/applib/lib/kml.c src/apps/applib/lib/utility.c - src/apps/applib/lib/args.c - src/h3lib/include/area.h) + src/apps/applib/lib/args.c) set(TEST_APP_SOURCE_FILES src/apps/applib/include/test.h src/apps/applib/lib/test.c) set(EXAMPLE_SOURCE_FILES @@ -312,7 +311,6 @@ set(OTHER_SOURCE_FILES src/apps/fuzzers/fuzzerCellToChildPos.c src/apps/fuzzers/fuzzerInternalAlgos.c src/apps/fuzzers/fuzzerInternalCoordIjk.c - src/apps/fuzzers/fuzzerGeoLoopArea.c src/apps/benchmarks/benchmarkPolygonToCells.c src/apps/benchmarks/benchmarkPolygonToCellsExperimental.c src/apps/benchmarks/benchmarkPolygon.c @@ -642,7 +640,6 @@ if(BUILD_FUZZERS) add_h3_fuzzer(fuzzerPolygonToCellsExperimentalNoHoles src/apps/fuzzers/fuzzerPolygonToCellsExperimentalNoHoles.c) add_h3_fuzzer(fuzzerCellToChildPos src/apps/fuzzers/fuzzerCellToChildPos.c) - add_h3_fuzzer(fuzzerGeoLoopArea src/apps/fuzzers/fuzzerGeoLoopArea.c) if(ENABLE_REQUIRES_ALL_SYMBOLS) add_h3_fuzzer(fuzzerInternalAlgos src/apps/fuzzers/fuzzerInternalAlgos.c) diff --git a/src/apps/fuzzers/README.md b/src/apps/fuzzers/README.md index a559075950..192d4a37d5 100644 --- a/src/apps/fuzzers/README.md +++ b/src/apps/fuzzers/README.md @@ -65,7 +65,6 @@ The public API of H3 is covered in the following fuzzers: | originToDirectedEdges | [fuzzerDirectedEdge](./fuzzerDirectedEdge.c) | polygonToCells | [fuzzerPoylgonToCells](./fuzzerPolygonToCells.c) | polygonToCellsExperimental | [fuzzerPoylgonToCellsExperimental](./fuzzerPolygonToCellsExperimental.c) [fuzzerPoylgonToCellsExperimentalNoHoles](./fuzzerPolygonToCellsExperimentalNoHoles.c) -| geoLoopArea | [fuzzerGeoLoopArea](./fuzzerGeoLoopArea.c) | radsToDegs | Trivial | stringToH3 | [fuzzerIndexIO](./fuzzerIndexIO.c) | uncompactCells | [fuzzerCompact](./fuzzerCompact.c) diff --git a/src/apps/fuzzers/fuzzerGeoLoopArea.c b/src/apps/fuzzers/fuzzerGeoLoopArea.c deleted file mode 100644 index eb39842b7a..0000000000 --- a/src/apps/fuzzers/fuzzerGeoLoopArea.c +++ /dev/null @@ -1,37 +0,0 @@ -/* - * Copyright 2025 Uber Technologies, Inc. - * - * Licensed under the Apache License, Version 2.0 (the "License"); - * you may not use this file except in compliance with the License. - * You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ -/** @file - * @brief Fuzzer program for geoLoopAreaRads2 - */ - -#include "aflHarness.h" -#include "area.h" -#include "h3api.h" -#include "utility.h" - -int LLVMFuzzerTestOneInput(const uint8_t *data, size_t size) { - // Interpret the input buffer as an array of LatLng vertices. - GeoLoop loop = { - .numVerts = (int)(size / sizeof(LatLng)), - .verts = (LatLng *)data, - }; - double area; - geoLoopAreaRads2(loop, &area); - - return 0; -} - -AFL_HARNESS_MAIN(sizeof(LatLng) * 1024); From a05939767c393eec9c04c0bfeb52b668a51819b0 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Mon, 1 Dec 2025 19:28:15 -0500 Subject: [PATCH 55/62] minor --- CHANGELOG.md | 2 +- website/docs/api/regions.mdx | 2 ++ 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 0e68810554..73f71e5dee 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -8,7 +8,7 @@ The public API of this library consists of the functions declared in file ## [Unreleased] ### Added - `reverseDirectedEdge` function (#1098) -- `geoLoopArea` function (#1101) +- (internal) `geoLoopArea` function (#1101) ### Changed - `cellAreaRads2` uses `geoLoopArea` (#1101) diff --git a/website/docs/api/regions.mdx b/website/docs/api/regions.mdx index ceebf05ade..a5bf4c8510 100644 --- a/website/docs/api/regions.mdx +++ b/website/docs/api/regions.mdx @@ -661,3 +661,5 @@ This function exists for memory management and is not exposed. + + From b24442f334cd7552348025ed9a627304b67b3a78 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Mon, 1 Dec 2025 19:54:42 -0500 Subject: [PATCH 56/62] remove justfile --- justfile | 53 ----------------------------------------------------- 1 file changed, 53 deletions(-) delete mode 100644 justfile diff --git a/justfile b/justfile deleted file mode 100644 index cffb5a6fde..0000000000 --- a/justfile +++ /dev/null @@ -1,53 +0,0 @@ -# TODO: Remove before landing - -init: purge - mkdir build - -build: - cd build; cmake -DCMAKE_BUILD_TYPE=Release ..; make - -profile: init - cd build; cmake -DCMAKE_BUILD_TYPE=RelWithDebInfo -DCMAKE_C_FLAGS='-fno-omit-frame-pointer' ..; make - # cd build; cmake -DCMAKE_BUILD_TYPE=Release ..; make - xcrun xctrace record \ - --template 'Time Profiler' \ - --output h3-prof.trace \ - --launch -- ./build/bin/benchmarkCellsToLinkedMultiPolygon - open h3-prof.trace - -purge: - rm -rf build - rm -rf *.trace - rm -rf .ipynb_checkpoints - -test: build - # ./build/bin/testH3CellAreaExhaustive - # ./build/bin/testEdgeCellsToPoly - # ./build/bin/testDirectedEdge - # ./build/bin/testArea - # just test-slow - just test-fast - # ./build/bin/testGeoLoopArea - -time: - time ./build/bin/testArea - -test-fast: build - cd build; make test-fast - -test-slow: build - cd build; make test - -bench: build - ./build/bin/benchmarkArea - # ./build/bin/benchmarkCellsToLinkedMultiPolygon - # ./build/bin/benchmarkCellsToMultiPolygon - # ./build/bin/benchmarkDirectedEdge - -fuzz: - cd build; CC=clang cmake -DENABLE_LIBFUZZER=ON ..; make fuzzers - -coverage: purge - mkdir build - cd build; cmake -DCMAKE_BUILD_TYPE=Debug -DENABLE_COVERAGE=ON ..; make; make coverage - open build/coverage/index.html From dd316e0a1b30c647014118304293f42131641508 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Tue, 2 Dec 2025 04:43:40 -0500 Subject: [PATCH 57/62] Adder adder = {}; --- src/apps/benchmarks/benchmarkArea.c | 2 +- src/apps/testapps/testH3CellAreaExhaustive.c | 2 +- src/h3lib/include/adder.h | 2 +- src/h3lib/lib/area.c | 2 +- 4 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/apps/benchmarks/benchmarkArea.c b/src/apps/benchmarks/benchmarkArea.c index 193f6e66f2..cbeba2d8b8 100644 --- a/src/apps/benchmarks/benchmarkArea.c +++ b/src/apps/benchmarks/benchmarkArea.c @@ -9,7 +9,7 @@ #include "iterators.h" static void doResSum(int res, bool print) { - Adder adder = {0}; + Adder adder = {}; double cellArea; IterCellsResolution iter = iterInitRes(res); diff --git a/src/apps/testapps/testH3CellAreaExhaustive.c b/src/apps/testapps/testH3CellAreaExhaustive.c index 5a8bd30aa6..784aa51bde 100644 --- a/src/apps/testapps/testH3CellAreaExhaustive.c +++ b/src/apps/testapps/testH3CellAreaExhaustive.c @@ -137,7 +137,7 @@ static void cell_area_assert(H3Index cell) { */ static void earth_area_test(int res, H3Error (*cell_area)(H3Index, double *), double target, double tol) { - Adder adder = {0}; + Adder adder = {}; IterCellsResolution iter = iterInitRes(res); for (; iter.h; iterStepRes(&iter)) { diff --git a/src/h3lib/include/adder.h b/src/h3lib/include/adder.h index dfebb0bcf1..9796ae520a 100644 --- a/src/h3lib/include/adder.h +++ b/src/h3lib/include/adder.h @@ -47,7 +47,7 @@ tradeoffs between Kahan, Neumaier, and `fsum`. ## Usage -Adder adder = {0}; // initialize +Adder adder = {}; // initialize kadd(&adder, x); // add x to summation out = adder.sum; // extract final sum */ diff --git a/src/h3lib/lib/area.c b/src/h3lib/lib/area.c index e80d25da5f..2394373e72 100644 --- a/src/h3lib/lib/area.c +++ b/src/h3lib/lib/area.c @@ -81,7 +81,7 @@ static inline double cagnoli(LatLng x, LatLng y) { H3Error geoLoopAreaRads2(GeoLoop loop, double *out) { // Use `Adder` to improve numerical accuracy of the sum of many Cagnoli // terms - Adder adder = {0}; + Adder adder = {}; for (int i = 0; i < loop.numVerts; i++) { int j = (i + 1) % loop.numVerts; From 0c7efce1e97f7c9af73fef8bf1d9c92b34d34a21 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Tue, 2 Dec 2025 07:41:12 -0500 Subject: [PATCH 58/62] degenerate loop note --- src/apps/testapps/testGeoLoopArea.c | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/src/apps/testapps/testGeoLoopArea.c b/src/apps/testapps/testGeoLoopArea.c index fd9ba1b230..a8331451d6 100644 --- a/src/apps/testapps/testGeoLoopArea.c +++ b/src/apps/testapps/testGeoLoopArea.c @@ -202,7 +202,8 @@ SUITE(geoLoopArea) { } TEST(degenerateLoop2) { - // TODO: should we raise an error when given degenerate loops? + // Note that `geoLoopAreaRads2()` works without error on degenerate + // loops, returning 0 area. LatLng verts[] = { {M_PI_2, 0.0}, {0.0, -M_PI_2}, @@ -211,7 +212,8 @@ SUITE(geoLoopArea) { } TEST(degenerateLoop1) { - // TODO: should we raise an error when given degenerate loops? + // Note that `geoLoopAreaRads2()` works without error on degenerate + // loops, returning 0 area. LatLng verts[] = { {0.0, 0.0}, }; @@ -219,7 +221,8 @@ SUITE(geoLoopArea) { } TEST(degenerateLoop0) { - // TODO: should we raise an error when given degenerate loops? + // Note that `geoLoopAreaRads2()` works without error on degenerate + // loops, returning 0 area. _compareArea(NULL, 0, 0.0); } } From ce65b2d04d3db0f80b47fea4801ce76808bf7f68 Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Wed, 10 Dec 2025 12:18:58 -0800 Subject: [PATCH 59/62] use same public/private underscore convention as iterators.h --- src/h3lib/include/adder.h | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/h3lib/include/adder.h b/src/h3lib/include/adder.h index 9796ae520a..715470dd67 100644 --- a/src/h3lib/include/adder.h +++ b/src/h3lib/include/adder.h @@ -53,14 +53,14 @@ out = adder.sum; // extract final sum */ typedef struct { - double sum; // running total - double c; // compensation term + double sum; // running total. public + double _c; // compensation term. private } Adder; static inline void kadd(Adder *adder, double x) { - double y = x - adder->c; + double y = x - adder->_c; double t = adder->sum + y; - adder->c = (t - adder->sum) - y; + adder->_c = (t - adder->sum) - y; adder->sum = t; } From 2dac0471b96c52467389a219afe037d8430031ee Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Tue, 16 Dec 2025 19:41:34 -0500 Subject: [PATCH 60/62] move other cellArea* functions to area.c --- src/h3lib/lib/area.c | 30 ++++++++++++++++++++++++++++++ src/h3lib/lib/latLng.c | 30 ------------------------------ 2 files changed, 30 insertions(+), 30 deletions(-) diff --git a/src/h3lib/lib/area.c b/src/h3lib/lib/area.c index 2394373e72..166059fcd9 100644 --- a/src/h3lib/lib/area.c +++ b/src/h3lib/lib/area.c @@ -124,3 +124,33 @@ H3Error H3_EXPORT(cellAreaRads2)(H3Index cell, double *out) { return E_SUCCESS; } + +/** + * Area of H3 cell in kilometers^2. + * + * @param cell H3 cell + * @param out cell area in kilometers^2 + * @return E_SUCCESS on success, or an error code otherwise + */ +H3Error H3_EXPORT(cellAreaKm2)(H3Index cell, double *out) { + H3Error err = H3_EXPORT(cellAreaRads2)(cell, out); + if (!err) { + *out *= EARTH_RADIUS_KM * EARTH_RADIUS_KM; + } + return err; +} + +/** + * Area of H3 cell in meters^2. + * + * @param cell H3 cell + * @param out cell area in meters^2 + * @return E_SUCCESS on success, or an error code otherwise + */ +H3Error H3_EXPORT(cellAreaM2)(H3Index cell, double *out) { + H3Error err = H3_EXPORT(cellAreaKm2)(cell, out); + if (!err) { + *out *= 1000 * 1000; + } + return err; +} diff --git a/src/h3lib/lib/latLng.c b/src/h3lib/lib/latLng.c index 13bf196e8e..082db77460 100644 --- a/src/h3lib/lib/latLng.c +++ b/src/h3lib/lib/latLng.c @@ -353,36 +353,6 @@ H3Error H3_EXPORT(getNumCells)(int res, int64_t *out) { return E_SUCCESS; } -/** - * Area of H3 cell in kilometers^2. - * - * @param cell H3 cell - * @param out cell area in kilometers^2 - * @return E_SUCCESS on success, or an error code otherwise - */ -H3Error H3_EXPORT(cellAreaKm2)(H3Index cell, double *out) { - H3Error err = H3_EXPORT(cellAreaRads2)(cell, out); - if (!err) { - *out = *out * EARTH_RADIUS_KM * EARTH_RADIUS_KM; - } - return err; -} - -/** - * Area of H3 cell in meters^2. - * - * @param cell H3 cell - * @param out cell area in meters^2 - * @return E_SUCCESS on success, or an error code otherwise - */ -H3Error H3_EXPORT(cellAreaM2)(H3Index cell, double *out) { - H3Error err = H3_EXPORT(cellAreaKm2)(cell, out); - if (!err) { - *out = *out * 1000 * 1000; - } - return err; -} - /** * Length of a directed edge in radians. * From 27911ae3d1c38bf76ac52c665533510ca86a7f8e Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Tue, 16 Dec 2025 19:45:29 -0500 Subject: [PATCH 61/62] Stitch two 1/4 triangles --- src/apps/testapps/testGeoLoopArea.c | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/apps/testapps/testGeoLoopArea.c b/src/apps/testapps/testGeoLoopArea.c index a8331451d6..49643ca245 100644 --- a/src/apps/testapps/testGeoLoopArea.c +++ b/src/apps/testapps/testGeoLoopArea.c @@ -108,11 +108,11 @@ SUITE(geoLoopArea) { TEST(hemisphere_east) { /* - Stitch 4 1/8 triangles together to cover the eastern hemisphere. + Stitch two 1/4 triangles together to cover the eastern hemisphere. */ LatLng verts[] = { {M_PI_2, 0.0}, - {0.0, 0}, + {0.0, 0.0}, {-M_PI_2, 0.0}, {0.0, M_PI}, }; From 49318002a17cef244380401d7ba4b242ec31c51b Mon Sep 17 00:00:00 2001 From: AJ Friend Date: Tue, 16 Dec 2025 19:53:49 -0500 Subject: [PATCH 62/62] add copyright headers --- src/apps/benchmarks/benchmarkArea.c | 15 +++++++++++++++ src/h3lib/include/area.h | 19 +++++++++++++++++++ 2 files changed, 34 insertions(+) diff --git a/src/apps/benchmarks/benchmarkArea.c b/src/apps/benchmarks/benchmarkArea.c index cbeba2d8b8..fa1b29e75e 100644 --- a/src/apps/benchmarks/benchmarkArea.c +++ b/src/apps/benchmarks/benchmarkArea.c @@ -1,3 +1,18 @@ +/* + * Copyright 2025 Uber Technologies, Inc. + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ #include #include #include diff --git a/src/h3lib/include/area.h b/src/h3lib/include/area.h index bf9235eb48..c561f53773 100644 --- a/src/h3lib/include/area.h +++ b/src/h3lib/include/area.h @@ -1,3 +1,22 @@ +/* + * Copyright 2025 Uber Technologies, Inc. + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +/** @file area.h + * @brief Area computation functions + */ + #ifndef AREA_H #define AREA_H