RelateNG prepared mode rebuilds the B-side geometry every call: ~35 kB/call in a one-target x many-geometries sweep

[asinghvi17]

prepare(RelateNG(), a) caches the A side and relate_predicate(prep, pred, b) reuses it — that part works, and it is what we use downstream. But evaluate! builds a fresh RelateGeometry for b on every call:

GeometryOps.jl/src/methods/geom_relations/relateng/relate_ng.jl

Lines 272 to 280 in 4e0c0d5

	function evaluate!(alg::RelateNG, geom_a::RelateGeometry, b, predicate::TopologyPredicate,
	prep = nothing)
	#-- Java performs the envelope fast-exit before building the B
	#-- RelateGeometry, reading the envelope cached on the Geometry. Here the
	#-- RelateGeometry constructor is what caches the extents (one coordinate
	#-- pass either way), so B is built first and the check reads its extent.
	geom_b = RelateGeometry(geom_a.m, b; exact = geom_a.exact,
	boundary_rule = geom_a.boundary_rule)
	if !has_required_envelope_interaction(geom_a, get_extent(geom_b), predicate)

and there is no hook to reuse that work or its buffers across calls. In the one-target-versus-many-geometries sweep that prepared mode exists to serve, the B side is then the entire remaining cost.

Measured on main @ 4e0c0d5 (Apple M-series, Julia 1.11), one prepared Spherical() polygon target, pred_contains (the converse of a Within query). Full script in a comment below.

Per-call cost is flat in the target's size — prepared mode is doing its job:

target vertices	4	32	72	312	616	736
KiB/call	7.5	7.5	7.5	7.5	7.5	7.5
us/call	8.35	8.38	8.54	8.19	8.30	8.00

and linear in the query geometry's size, at ~1 KiB and ~0.8 us per vertex:

query vertices	3	6	16	32	64	128
KiB/call	4.5	7.5	17.7	33.8	66.1	132.0
us/call	6.23	7.86	15.29	27.79	55.50	102.75

For scale, prepare on those same targets costs 7.3 KiB / 13 us (4 verts) up to 899.5 KiB / 645 us (736 verts) — so against a large target a single 128-vertex query call is already ~15% of a full re-preparation, and against a small one it is several times the preparation.

A point query is the one cheap case, at 0.1 KiB / 0.01 us. That is a query-shape short-circuit, not a predicate one: pred_contains, pred_within and pred_intersects on the same 6-vertex ring all cost an identical 7.5 KiB.

Access pattern: one target × thousands of small cells. In DiscreteGlobalGrids' polygon-coverage traversal the same prepared outline is tested against every candidate cell of a DGGS level. Measured there, the predicate is 79% of wall time at H3 L7 and 81% at L9, and 90% of allocations (3.26 GiB total at L9). Removing the per-call rebuild is worth a measured ~3.7-4.2x end-to-end and ~6.6x memory, which would put a generic spherical DFS ahead of libh3's specialised polygonToCells + compactCells at every level tested. It is expensive enough today that we deliberately carry a cheap provenance-only "contained" label rather than ask the predicate per boundary leaf.

The B geometries genuinely differ call to call, so rebuilding something is unavoidable — but in this pattern they are all rings of the same small, fixed vertex count, so it is per-call scratch churn rather than retained state. A workspace/buffer argument (relate_predicate(prep, pred, b; workspace)) or an in-place RelateGeometry!(scratch, m, b; …) would cover it.

Separately: the clipping algorithms accept no prepared geometry at all — intersection(ConvexConvexSutherlandHodgman, a, b) re-walks both rings every call, with the same one-destination × many-sources access pattern in ConservativeRegridding.

Tested revs: measured on main @ 4e0c0d5; same behaviour on as/dualdfs-carry-extents @ a658e8b.

(#87 asked for prepared geometry in general and is closed — prepare landed. This is the residual half: the unprepared side of a prepared query.) Related: #418, #408, #23.

[asinghvi17]

MWE

# fresh env
using Pkg
Pkg.activate(temp = true)
Pkg.add(url = "https://github.com/JuliaGeo/GeometryOps.jl.git", rev = "main")  # 4e0c0d5
Pkg.add("GeoInterface")

# MWE 4: per-call preparation cost of the B side in RelateNG prepared mode
import GeometryOps as GO, GeoInterface as GI
const US = GO.UnitSpherical
to_sphere = US.UnitSphereFromGeographic()

# one TARGET, prepared once
target = GI.Polygon([GI.LinearRing([to_sphere((60cosd(t), 60sind(t))) for t in 0:5:360])])
prep = GO.prepare(GO.RelateNG(; manifold = GO.Spherical()), target)

# many small 32-vertex "cells" tested against it
cell(lon, lat) = GI.Polygon([GI.LinearRing(
    [to_sphere((lon + 2cosd(t), lat + 2sind(t))) for t in range(0, 360; length = 33)])])
cells = [cell(lon, lat) for lon in -20:2:20, lat in -20:2:20] |> vec
centroids = [to_sphere((lon, lat)) for lon in -20:2:20, lat in -20:2:20] |> vec
println(length(cells), " cells of ",
    length(collect(GI.getpoint(GI.getexterior(cells[1])))), " vertices each")

# warm up
GO.relate_predicate(prep, GO.pred_contains(), cells[1])
GO.relate_predicate(prep, GO.pred_intersects(), cells[1])
GO.relate_predicate(prep, GO.pred_contains(), centroids[1])

for (name, pred, geoms) in [
    ("contains(prep, CELL)   [the Within query]", GO.pred_contains(), cells),
    ("intersects(prep, CELL)", GO.pred_intersects(), cells),
    ("contains(prep, POINT)  [Intersects fast path]", GO.pred_contains(), centroids),
]
    b = @allocated for g in geoms; GO.relate_predicate(prep, pred, g); end
    t = @elapsed for g in geoms; GO.relate_predicate(prep, pred, g); end
    println(rpad(name, 46), lpad(round(Int, b / length(geoms)), 8), " B/call  ",
        lpad(round(1e6 * t / length(geoms); digits = 2), 8), " us/call")
end

Actual output on main @ 4e0c0d5 (Apple M-series, Julia 1.11):

441 cells of 33 vertices each
contains(prep, CELL)   [the Within query]        35199 B/call     36.40 us/call
intersects(prep, CELL)                           34574 B/call     30.82 us/call
contains(prep, POINT)  [Intersects fast path]      249 B/call      0.07 us/call

The target is prepared once, outside the loop. The 35 kB is entirely the per-call B-side RelateGeometry build — ~1 kB and ~1 us per vertex of the query ring, versus 249 B for a point against the same prepared target.

The two access patterns

1. One target x many cells. DiscreteGlobalGrids' polygon coverage prepares the query outline once and then tests it against every candidate cell of a DGGS traversal — tens to hundreds of thousands of cells at a deep level, i.e. order 1 GB of transient allocation for a single query. This is the dominant cost of the traversal, and the reason we currently carry a cheap provenance-only "contained" label rather than asking the predicate per boundary leaf; the exact label would be nearly free if the B side were.

2. One cell x its own ring. The same rebuild happens for the small, fixed-size geometry a caller passes over and over with only its coordinates changed.

Both are the shape prepared mode is for, and in both the B geometries are rings of identical vertex count — so the scratch could be allocated once and refilled. relate_predicate(prep, pred, b; workspace), or an in-place RelateGeometry!(scratch, m, b; ...), would cover it.

[asinghvi17]

Full benchmark script

One harness for every number in the OP, so the rows are directly comparable.

# fresh env
using Pkg
Pkg.activate(temp = true)
Pkg.add(url = "https://github.com/JuliaGeo/GeometryOps.jl.git", rev = "main")  # 4e0c0d5
Pkg.add("GeoInterface")

# MWE 4c: one harness, isolate target-scaling vs cell-scaling vs cell POSITION
import GeometryOps as GO, GeoInterface as GI
const US = GO.UnitSpherical
to_sphere = US.UnitSphereFromGeographic()
poly(n, r, lon, lat) = GI.Polygon([GI.LinearRing(
    [to_sphere((lon + r*cosd(t), lat + r*sind(t))) for t in range(0, 360; length = n)])])
prep(t) = GO.prepare(GO.RelateNG(; manifold = GO.Spherical()), t)

# ONE harness for every measurement below.
function bench(p, pr, g; n = 300)
    GO.relate_predicate(p, pr, g)                       # warm up / compile
    b = @allocated for _ in 1:n; GO.relate_predicate(p, pr, g); end
    t = @elapsed for _ in 1:n; GO.relate_predicate(p, pr, g); end
    return (b / n / 1024, 1e6 * t / n)                  # KiB/call, us/call
end
row(lbl, x) = println("  ", rpad(lbl, 40), lpad(round(x[1]; digits = 1), 8), " KiB ",
                      lpad(round(x[2]; digits = 2), 9), " us")

const R = 40.0
# a 6-vertex cell DEEP INSIDE the target, and one STRADDLING its boundary
inside(r) = poly(7, 1.0, 3.0, 3.0)
straddle(r) = poly(7, 1.0, r, 0.0)

println("=== vary TARGET vertices; cell = 6 verts, DEEP INSIDE ===")
for nt in (5, 33, 73, 313, 617, 737)
    row("target $(nt-1) verts", bench(prep(poly(nt, R, 0.0, 0.0)), GO.pred_contains(), inside(R)))
end
println("=== vary TARGET vertices; cell = 6 verts, STRADDLING the boundary ===")
for nt in (5, 33, 73, 313, 617, 737)
    row("target $(nt-1) verts", bench(prep(poly(nt, R, 0.0, 0.0)), GO.pred_contains(), straddle(R)))
end
println("=== vary CELL vertices; target = 72 verts, cell DEEP INSIDE ===")
let p = prep(poly(73, R, 0.0, 0.0))
    for nc in (4, 7, 17, 33, 65, 129)
        row("cell $(nc-1) verts", bench(p, GO.pred_contains(), poly(nc, 1.0, 3.0, 3.0)))
    end
end
println("=== predicate comparison; target = 72 verts, cell = 6 verts INSIDE ===")
let p = prep(poly(73, R, 0.0, 0.0)), c = inside(R)
    row("pred_contains (the Within query)", bench(p, GO.pred_contains(), c))
    row("pred_within",                      bench(p, GO.pred_within(),   c))
    row("pred_intersects",                  bench(p, GO.pred_intersects(), c))
    row("pred_contains, POINT",             bench(p, GO.pred_contains(), to_sphere((3.0, 3.0))))
end
println("=== prepare(target) cost, for scale ===")
for nt in (5, 73, 617, 737)
    t = poly(nt, R, 0.0, 0.0); prep(t)
    b = @allocated prep(t); e = @elapsed prep(t)
    row("prepare, target $(nt-1) verts", (b/1024, 1e6*e))
end

Actual output on main @ 4e0c0d5:

=== vary TARGET vertices; cell = 6 verts, DEEP INSIDE ===
  target 4 verts                               7.5 KiB      8.35 us
  target 32 verts                              7.5 KiB      8.38 us
  target 72 verts                              7.5 KiB      8.54 us
  target 312 verts                             7.5 KiB      8.19 us
  target 616 verts                             7.5 KiB       8.3 us
  target 736 verts                             7.5 KiB       8.0 us
=== vary TARGET vertices; cell = 6 verts, STRADDLING the boundary ===
  target 4 verts                               7.5 KiB      9.09 us
  target 32 verts                              7.5 KiB      7.13 us
  target 72 verts                              7.5 KiB      7.51 us
  target 312 verts                             7.5 KiB      7.23 us
  target 616 verts                             7.5 KiB      7.85 us
  target 736 verts                             7.5 KiB      7.08 us
=== vary CELL vertices; target = 72 verts, cell DEEP INSIDE ===
  cell 3 verts                                 4.5 KiB      6.23 us
  cell 6 verts                                 7.5 KiB      7.86 us
  cell 16 verts                               17.7 KiB     15.29 us
  cell 32 verts                               33.8 KiB     27.79 us
  cell 64 verts                               66.1 KiB      55.5 us
  cell 128 verts                             132.0 KiB    102.75 us
=== predicate comparison; target = 72 verts, cell = 6 verts INSIDE ===
  pred_contains (the Within query)             7.5 KiB      8.58 us
  pred_within                                  7.5 KiB      7.33 us
  pred_intersects                              7.5 KiB      8.15 us
  pred_contains, POINT                         0.1 KiB      0.01 us
=== prepare(target) cost, for scale ===
  prepare, target 4 verts                      7.3 KiB     13.21 us
  prepare, target 72 verts                    87.9 KiB     67.46 us
  prepare, target 616 verts                  738.6 KiB     579.5 us
  prepare, target 736 verts                  899.5 KiB     645.38 us

Two hypotheses this rules out

Downstream measurements on real DGGS fixtures (a 5-vertex box, a 617-vertex California outline, a 734-vertex multipolygon) suggested the per-call cost was target-side and grew with the target's vertex count — 7.9 KiB, 263.1 KiB, 385.4 KiB respectively. That does not reproduce in isolation: the first table above varies the target across exactly that range, with the query held at 6 vertices, and the cost is flat at 7.5 KiB. The straddling variant rules out the obvious confounder (a boundary-adjacent cell touching more of the target's topology). The likely explanation for the downstream figures is that swapping the target also changes which cells the traversal visits and how many of them are boundary cells, so the per-call figure there is an average over a different cell population rather than a function of the target.

It also rules out an intersects-only short-circuit in RelateNG's prepared path: pred_intersects costs the same 7.5 KiB as pred_contains and pred_within. Downstream Intersects is cheap for a different reason — DiscreteGlobalGrids' own query code has a centroid fast-accept that passes a point to relate_predicate, and it is the point-shaped query that costs 0.1 KiB. So the cheap path already exists in RelateNG; it just is not reachable for a ring-shaped query.

That makes the fix query-side rather than target-side: reuse the B-side scratch, or give a small fixed-size ring the same kind of fast path a point already has.

[asinghvi17]

Related: #466, #467, #468 — same investigation (DiscreteGlobalGrids.jl conservative-regridding and polygon-coverage work). #466–#468 are all in the spherical _intersection_sutherland_hodgman, and a single tested patch on #466 covers #466, #467 and the first arm of #468; #469 is independent (RelateNG prepared mode).