Skip to content

perf: cover polygons by descending the geohash tree - #31

Merged
mailletf merged 8 commits into
mainfrom
perf/hierarchical-descent
Sep 16, 2026
Merged

mailletf merged 8 commits into
mainfrom
perf/hierarchical-descent

Conversation

@mailletf

@mailletf mailletf commented Aug 25, 2026 •

Copy link
Copy Markdown
Member

Problem

polygons_to_geohashes flood-filled at the target precision, testing every cell against the full polygon. Cost tracked the polygon's area, and each test walked every vertex.

Fix

Walk the geohash tree top-down. A cell disjoint from the polygon prunes its whole subtree; a cell the polygon contains yields all 32ⁿ descendants with no further geometry tests; only cells straddling the boundary subdivide. Cost now tracks the perimeter.

One PreparedGeometry::relate per visited cell answers intersects and contains together, and its R*-tree keeps each call proportional to the edges near that cell rather than the polygon's vertex count.

Both halves are needed. Measured separately, plain Contains inside the descent is ~6× slower than the old code, and relate without the descent is a net loss for the intersects-only case (943 ms → 1570 ms on verdun p9), because relate computes the full DE-9IM matrix where Intersects short-circuits.

Benchmarks

case before after
verdun p8 inner=false 48.4 ms 17.4 ms 2.8×
verdun p8 inner=true 86.4 ms 16.6 ms 5.2×
verdun p9 inner=false 943 ms 121 ms 7.8×
verdun p9 inner=true 2524 ms 121 ms 20.9×
verdun p10 inner=false 32.8 s 2.58 s 12.7×
verdun p10 inner=true 90.4 s 2.48 s 36.5×
whitehorse p7 inner=false 1472 ms 146 ms 10.1×
whitehorse p7 inner=true 3396 ms 141 ms 24.1×
whitehorse p8 inner=false 53.5 s 5.48 s 9.8×
whitehorse p8 inner=true 117 s 5.85 s 20.0×

Output sets are identical to the previous implementation in every case. Below p6 it is a wash — the gain scales with precision because interior cells stop costing anything.

Also

  • New ghbits module: a geohash of precision 1..=12 is 5·p interleaved bits and fits in a u64, so cells live in an FxHashSet<u64> and a cell's bbox is a couple of shifts rather than a base32 decode. Verified against the geohash crate across all 12 precisions, both bit parities, and the grid corners.
  • Precision is validated up front (1..=12) instead of failing later inside encode.
  • The descent seeds at the deepest level whose cover is still small, so it does not spend a relate call per cell on levels where the cell dwarfs the polygon and the R*-tree cannot prune.

⚠️ Breaking

Removes pub fn seed_interior_point_fast. The descent starts from the grid and never needs a point known to be inside the polygon, so the whole centroid → 24 probes → bbox centre → 4×4 grid → interior_point ladder has no callers. Deleting it also removes a failure mode rather than just dead code: the old walk silently skipped any polygon the ladder could not seed.

Tests

A brute-force oracle over both WKT fixtures checking every cell in the bounding box with no descent and no prepared geometry, plus a refinement-consistency check. The MultiPolygon regression tests from #30 still pass — the descent walks each part over the whole grid, so that class of bug is now structurally impossible.

Behaviour changes (deliberate)

  • inner=true keeps boundary-touching cells. Containment is now DE-9IM contains — the cell's interior inside the polygon, boundaries allowed to meet — for all polygons. The old fast path for hole-free polygons rejected any cell sharing an edge or vertex with the polygon boundary, so grid-aligned geometry gets more cells than before: a polygon exactly equal to one p5 cell now yields all 32 of its p6 children with inner=true, not just the 12 interior ones. This matches shapely's contains and how holed polygons were always treated; off the grid lines a boundary cell genuinely crosses the edge and is excluded either way, so non-aligned covers are unchanged. Pinned by test_inner_cover_keeps_boundary_touching_cells (Rust + Python).
  • Out-of-range coordinates now error. cover_range clamps grid indices, so a polygon past the antimeridian or a pole would have come back as a silently empty cover, and one straddling the domain edge as a half cover. The polygon's bounding box is validated up front and raises InvalidCoordinateRange, restoring the old flood-fill's contract.

@mailletf

mailletf commented Aug 25, 2026 •

Copy link
Copy Markdown
Member Author

Copilot AI left a comment

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Pull request overview

This pull request optimizes polygon-to-geohash coverage using hierarchical descent, packed geohash bits, and prepared geometry.

Changes:

  • Adds top-down pruning and subtree emission.
  • Adds packed geohash encoding and precision validation.
  • Removes obsolete interior-point seeding.
  • Updates dependencies and regression documentation.

Reviewed changes

Copilot reviewed 4 out of 5 changed files in this pull request and generated 1 comment.

Show a summary per file
File Summary Review status
tests/test_geohasher.py Updates crescent regression documentation. No findings.
tests/conftest.py Updates crescent fixture documentation. No findings.
src/lib.rs Implements hierarchical coverage and packed geohashes. Moderate: Boundary-aligned polygon minima can miss west/south touching cells. Use boundary-inclusive lower-index calculation and add a regression test.
Cargo.toml Adds rustc-hash. No findings.
Cargo.lock Locks the new dependency. No findings.

💡 Add a code-review agent skill or configure MCP servers for context-aware, tailored reviews. Learn more in the docs.

Comment thread src/lib.rs Outdated
Comment on lines +207 to +220
let idx = |v: f64, span: f64, max: u64| -> u64 {
let raw = (v / span).floor();
if raw < 0.0 {
0
} else {
(raw as u64).min(max)
}
};
(
idx(bbox.min().x + 180.0, lon_span, lon_max),
idx(bbox.max().x + 180.0, lon_span, lon_max),
idx(bbox.min().y + 90.0, lat_span, lat_max),
idx(bbox.max().y + 90.0, lat_span, lat_max),
)

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Good catch — real bug, fixed in f924e36.

The asymmetry was the tell: floor on the upper bound already names the cell east/north of a grid-aligned maximum (the touching one), so only the lower bound was dropping them. A polygon equal to one p5 cell covered 4 of the 9 cells it meets. Switched the lower bound to ceil(v / span) - 1, clamped at zero, and added test_grid_aligned_polygon_reaches_touching_neighbours — it fails on the old code with exactly those 4 cells.

The seed range only has to be a superset; relate still decides what is kept, so the cost of an extra seed is one topological test.

@mailletf
mailletf changed the base branch from fix/multipolygon-dropped-parts to main August 26, 2026 17:59
@mailletf
mailletf force-pushed the perf/hierarchical-descent branch from f924e36 to 394db72 Compare September 1, 2026 18:19
polygons_to_geohashes flood-filled at the target precision, testing every
cell against the full polygon. Cost therefore tracked the polygon's *area*,
and each test walked every vertex.

Walk the geohash tree top-down instead. A cell disjoint from the polygon
prunes its whole subtree; a cell the polygon contains yields all 32^n
descendants with no further geometry tests; only cells straddling the
boundary are subdivided. Cost now tracks the polygon's perimeter.

One PreparedGeometry::relate per visited cell answers "intersects" and
"contains" together, and its R*-tree keeps each call proportional to the
edges near that cell rather than the polygon's vertex count.

Both halves are needed. Measured separately, plain Contains inside the
descent is ~6x slower than the old code, and relate without the descent is
a net loss for the intersects-only case (943ms -> 1570ms on verdun p9),
because relate computes the full DE-9IM matrix where Intersects
short-circuits. Together:

  verdun p8      inner=false    48.4ms ->  17.4ms    2.8x
  verdun p8      inner=true     86.4ms ->  16.6ms    5.2x
  verdun p9      inner=false     943ms ->   121ms    7.8x
  verdun p9      inner=true     2524ms ->   121ms   20.9x
  verdun p10     inner=false    32.8s  ->  2.58s    12.7x
  verdun p10     inner=true     90.4s  ->  2.48s    36.5x
  whitehorse p7  inner=false    1472ms ->   146ms   10.1x
  whitehorse p7  inner=true     3396ms ->   141ms   24.1x
  whitehorse p8  inner=false    53.5s  ->  5.48s     9.8x
  whitehorse p8  inner=true      117s  ->  5.85s    20.0x

Output sets are identical to the previous implementation in every case
above. Below p6 it is a wash; the gain scales with precision because
interior cells stop costing anything.

Supporting changes:

- New `ghbits` module: a geohash of precision 1..=12 is 5*p interleaved
  bits and fits in a u64, so cells live in an FxHashSet<u64> and a cell's
  bbox is a couple of shifts rather than a base32 decode. Verified against
  the geohash crate across all 12 precisions, both bit parities and the
  grid corners.
- Precision is now validated up front (1..=12, matching what the geohash
  crate accepts) instead of failing later inside encode.
- The descent seeds at the deepest level whose cover is still small, so it
  does not spend a relate call per cell on levels where the cell dwarfs the
  polygon and the R*-tree cannot prune.

Removes `seed_interior_point_fast`. The descent starts from the grid and
never needs a point known to be inside the polygon, so the whole
centroid -> 24 probes -> bbox centre -> 4x4 grid -> interior_point ladder
has no callers. Deleting it also removes a failure mode rather than just
dead code: the old walk silently skipped any polygon the ladder could not
find a seed for. This is a breaking change for Rust consumers of the rlib —
the function was `pub` — but it only ever existed to serve the flood fill.

That also clears the `items_after_test_module` clippy warning, since the
deleted function was the item sitting below `mod tests`.

Two test docstrings described the seed ladder and have been rewritten. The
crescent fixture itself stays: a thin C-shaped ring whose centroid and
bounding-box centre both fall in the hollow is still the sharpest test of a
cover, and it is exactly what the descent has to get right without a seed.

Tests add a brute-force oracle over both WKT fixtures that checks every
cell in the bounding box with no descent and no prepared geometry, plus a
refinement-consistency check. The multipolygon regression tests from the
previous commit still pass: the descent walks each part over the whole
grid, so that class of bug is now structurally impossible.

Squashed with:
- fix: seed the descent from cells the polygon only touches
- fix: reject out-of-range coordinates instead of covering nothing
- test: pin inner covers keeping boundary-touching cells

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
@mailletf
mailletf force-pushed the perf/hierarchical-descent branch from ed57150 to c9b999e Compare September 1, 2026 18:54
mailletf and others added 7 commits September 16, 2026 14:53
expand_geohash_set walked the grid through geohash::neighbors, which
allocates eight Strings per cell, and kept cells in a HashSet<String> hashed
with SipHash. The allocation churn dominated: the actual work is integer
arithmetic on a grid.

Pack each hash into its u64 form and walk with ghbits::neighbors, keeping
cells in an FxHashSet<u64> and rendering back to base32 once at the end.

Measured end to end from Python on whitehorse p6 (24,571 input cells),
output identical in every case:

  expansion_m       before     after
        0 m        15.1 ms    4.3 ms    3.5x
      700 m        16.8 ms    6.0 ms    2.8x
     3000 m        19.1 ms    6.4 ms    3.0x
    12000 m        31.4 ms    8.6 ms    3.7x

The zero-hop case also dropped a wasted pass. The old code scanned every
input cell for boundary membership before looking at n_hops, so an
expansion of 0 m paid the full O(N*8) neighbour scan to return its own
input. Input hashes are still validated in that case — packing does it — so
malformed input is rejected as before.

Ring expansion no longer needs that boundary pre-pass at all: cells already
in the set are never re-added, so interior cells fall out of the frontier
after the first pass on their own.

Two behaviour changes worth calling out:

- Latitude now clamps at the poles instead of wrapping. geohash::neighbors
  reports eight neighbours for a top-row cell such as "zzzzzzzzzzzz", three
  of which have wrapped over the pole to the opposite edge of the grid, so
  expanding a polar geography used to jump to the other hemisphere. Covered
  by test_ghbits_neighbors_clamp_at_poles and
  test_expand_at_the_pole_does_not_cross_over.
- A set mixing precisions is now an error rather than each cell expanding on
  its own grid. Every caller in this crate already validated uniform
  precision before calling in.

Squashed with:
- fix: reject out-of-range hash lengths in ghbits::pack

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
n_hops_for sizes an expansion against the smaller of a cell's height and
width. Cell width shrinks with cos(latitude), so the hop count climbs
steeply toward the poles. At precision 9, expanding 1 km needs 419 hops at
latitude 60, 1,206 at latitude 80, and around 1.5e8 in the top row. Each hop
adds a ring, so a polar input built a frontier of billions of cells and
never returned. If the narrowest extent ever reached zero the division gave
+inf, `as usize` saturated to usize::MAX, and the loop could not terminate
at all.

Refuse both cases with a message naming the offending hash, the hop count
and the cell's narrowest extent, so the caller can see whether to coarsen
the precision or shrink expansion_m. The cap is 10,000 hops: far past any
real expansion, far short of what a degenerate input produces. A 50 km
expansion at precision 6 still goes through.

Splits the sizing logic into n_hops_for_core, which returns Result<_,
String>, from the thin PyResult wrapper. Constructing a PyErr pulls in
Python symbols the Rust test binary does not link, so the policy was
otherwise unreachable from cargo test.

Also corrects two docstrings that said the hop count comes from a cell's
height; it has used min(height, width) since the function was written.

Squashed with:
- fix: size a group's expansion on its narrowest cell, not its first hash

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
polygon_to_geohashes took `_py: Python` and never used it, so the entire
cover walk ran with the GIL held. The walk touches no Python objects once
__geo_interface__ has been read, and at fine precisions it runs for a long
time — verdun at precision 10 takes minutes — blocking every other thread in
the interpreter for the duration. Every other heavy function in this module
already detaches.

Read the geometry under the GIL as before, then detach for the walk.

The test runs two concurrent calls and compares against one. Without the fix
it reports 2.06x, i.e. fully serialised.

Squashed with:
- test: make the GIL test robust off a multi-core runner

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
polygons_to_geohashes_handbrake had no callers outside benches/bench.rs,
where it served as the "oldfunc" comparison arm. It was also broken: a cell
failing the envelope test was recorded in neither the inner nor the outer
set and was not expanded from, so every neighbour that reached it queued it
again — up to eight times the work. Drop the function and its four
benchmark arms rather than fix a variant nothing calls.

That leaves several imports unused outside the test module. VecDeque and
geohash::neighbors are no longer needed at all in the library; Centroid was
the last thing the handbrake used; and Contains is now only reached by the
brute-force oracle, so it moves into `mod tests` alongside the comment
explaining why the oracle avoids the prepared-geometry path.

Also: expand_geohash_mapping_arrow reported a null geohash as "all geohashes
in a group must have the same precision", because a null reads back as an
empty string and so fails the length check. Detect the null first and name
its position.

Squashed with:
- docs: say which crate the neighbours oracle comes from

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
The tree has been drifting from rustfmt for a while, so `make check-rust`
fails on formatting alone and every code change has to dodge unrelated
reformatting in its diff. Run cargo fmt once, on its own, so later diffs
stay about their own change.

Pure formatting: no logic, no renames, no reordering beyond import sorting.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
For 500k precision-7 geohashes, the pure Rust core of decode_many_to_wkb
runs in under 4 ms, while the function called from Python takes 45 ms. Over
90% of the wall time is the boundary: converting list[str] into Vec<String>
on the way in, and building half a million Python bytes objects on the way
out. That is also why the num_threads knob does nothing here — it
parallelises the 4 ms.

Add decode_many_to_wkb_arrow and decode_many_to_ewkb_arrow, following
expand_geohash_mapping_arrow. Strings are read straight from the input's
Arrow buffers and the output is built into Arrow buffers, so no Python
object exists per row.

  N = 500,000, medians of 15 runs

  decode_many_to_wkb    (list -> list[bytes])    45.0 ms
  decode_many_to_wkb_arrow   (arrow -> arrow)     3.1 ms   14.5x
  decode_many_to_ewkb   (list -> list[bytes])    51.8 ms
  decode_many_to_ewkb_arrow  (arrow -> arrow)     3.2 ms   16.2x

Even when the data starts as a Python list and has to be converted, list ->
pa.array -> wkb_arrow is 10.5 ms, still 4.3x. Output is byte-identical to
the list API in every case.

Every row serializes to the same width — 93 bytes for WKB, 97 for EWKB — so
the values buffer is allocated once for the whole column and filled across
threads in place via par_chunks_exact_mut, with offsets computed from the
row width. That removes the per-row Vec the list API still pays.

Supporting changes:

- serialize_bbox now wraps write_bbox, which writes into a caller-provided
  slice. The allocating form stays for the existing list API.
- StringColumn accepts either Utf8 or LargeUtf8 input, so callers do not
  have to match the offset width the library happens to prefer.
- Null inputs produce null outputs rather than being rejected or silently
  decoded as an empty string.

The list-returning functions are untouched.

Squashed with:
- test: assert write_bbox fills its slot exactly
- test: put the n_hops_for tests back under their section

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Completes the Arrow bulk codec started for WKB output. encode_many,
decode_many and decode_many_exactly all spend most of their time building
Python objects — 500k strings, or 500k tuples — rather than doing geohash
work.

  N = 500,000 at precision 7, medians of 15 runs

  encode_many          (list -> list[str])       24.8 ms
  encode_many_arrow    (arrow -> arrow)           3.4 ms    7.3x
  decode_many          (list -> list[tuple])     52.8 ms
  decode_many_arrow    (arrow -> arrow)           1.9 ms   27.8x
  decode_many_exactly  (list -> list[tuple])     68.8 ms
  decode_many_exactly_arrow                       2.8 ms   24.6x

decode_many gains the most because a tuple per row is the most expensive
thing any of these functions does.

encode_many_arrow reuses the fixed-width column trick: every hash is exactly
`precision` characters, so the buffer is allocated once and filled across
threads. The decoders write one row-major cell buffer, so each row is
touched by exactly one thread, then split it into columns.

The decoders return a RecordBatch rather than an array of tuples — (lng,
lat) and (lng, lat, lng_err, lat_err) — which is the shape a dataframe or a
DuckDB table actually wants, and keeps each column contiguous.

Nulls propagate: a null geohash gives null coordinates, and a row is null if
either input coordinate is null. Precision is validated up front, so
encode_many_arrow reports a bad precision as such rather than failing per
row.

The list-returning functions are untouched.

Squashed with:
- perf: decode coordinates straight into the output columns

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
@mailletf
mailletf merged commit 7c3b53e into main Sep 16, 2026
4 checks passed
@mailletf
mailletf deleted the perf/hierarchical-descent branch September 16, 2026 20:02
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants