osm-katana/benches/clip_containment.rs
//! **What a region clip costs per node, three ways — measured, not argued.**
//!
//! Cutting Sweden out of the planet is `nodes × vertices`: 10.58 billion nodes
//! against relation 52822's 14 602-vertex border, where the downloaded outline it
//! replaced had 146. Three approaches attack different factors, and this bench
//! reports all of them over the SAME real node coordinates:
//!
//! 1. **Simplify the polygon** — attacks `vertices`. Douglas–Peucker at several
//! `eps`, each reported with its vertex count, the area it changes, and the
//! 43-point named-place verdict. An `eps` that loses one probe is not a result.
//! 2. **Fan the test across cores** — attacks `nodes`. Reported, not measured
//! here, because it is already true: `clip_keeps` runs inside the gatling
//! pass-1 workers (`reader::read_pass1_mmap_planned` → `par::par_map`;
//! `xml_reader`'s bz2/XML pass-1 codecs call it in `transform`), measured at
//! 29.9 of 32 cores busy on a Sweden clip (`bench_history.json`, 2026-06-18).
//! 3. **A raster acceleration mask** — attacks the cost of one test, and changes
//! no geometry. [`osm_katana::polygon_mask::PolygonMask`].
//!
//! The corpus is REAL: node coordinates sampled out of a `.osm.bz2` with
//! `xml_scan::sample_node_coords`, because how many nodes land in the region's
//! bbox — and how many of those sit near the boundary — is a fact about the
//! planet, not something a generator can be trusted to imitate.
//!
//! ```text
//! OSM_KATANA_CLIP_POLY=~/git/korp/assets/sweden.poly \
//! OSM_KATANA_CLIP_NODES=~/scratch/.../slice.bz2 \
//! OSM_KATANA_CLIP_STRIDE=8 \
//! cargo bench -p osm-katana --bench clip_containment
//! ```
//! Without `OSM_KATANA_CLIP_NODES` it falls back to the in-repo liechtenstein
//! fixture, which is honest but tiny — say which corpus any quoted number used.
use std::path::PathBuf;
use std::time::Instant;
use osm_katana::Polygon;
const DEFAULT_NODES: &str = concat!(
env!("CARGO_MANIFEST_DIR"),
"/../lbzip2/test_data/liechtenstein.osm.bz2"
);
/// The 43 named places the clip polygon is judged at — 24 that ARE Sweden, 19
/// that are not. Same set korp's `bundled_sweden_poly_contains_the_islands…`
/// guard uses, because a simplification is only acceptable if it survives the
/// check that caught the 48-vertex ring.
const INSIDE: &[(&str, f64, f64)] = &[
("Visby, Gotland", 18.29, 57.64),
("Borgholm, Öland", 16.66, 56.88),
("Fårö", 19.10, 57.94),
("Gotska Sandön", 19.26, 58.38),
("Kalmar", 16.36, 56.66),
("Karlskrona", 15.59, 56.16),
("Sandhamn", 18.92, 59.29),
("Gräsö", 18.45, 60.40),
("Holmön", 20.86, 63.77),
("Smygehuk", 13.36, 55.34),
("Treriksröset", 20.5486, 69.0596),
("Strömstad", 11.17, 58.94),
("Stockholm", 18.07, 59.33),
("Malmö", 13.00, 55.60),
("Kiruna", 20.22, 67.85),
("Göteborg", 11.97, 57.71),
("Umeå", 20.26, 63.83),
("Luleå", 22.15, 65.58),
("Östersund", 14.64, 63.18),
("Karlstad", 13.50, 59.38),
("Uppsala", 17.64, 59.86),
("Visingsö", 14.36, 58.05),
("Ven", 12.69, 55.91),
("Ystad", 13.82, 55.43),
];
const OUTSIDE: &[(&str, f64, f64)] = &[
("Helsingør, DK", 12.61, 56.03),
("Copenhagen, DK", 12.57, 55.68),
("Bornholm, DK", 14.90, 55.13),
("Aalborg, DK", 9.92, 57.05),
("Mariehamn, FI", 19.94, 60.10),
("Turku, FI", 22.27, 60.45),
("Helsinki, FI", 24.94, 60.17),
("Vaasa, FI", 21.62, 63.10),
("Kolari, FI", 23.79, 67.33),
("Oslo, NO", 10.75, 59.91),
("Narvik, NO", 17.43, 68.44),
("Trondheim, NO", 10.40, 63.43),
("Kirkenes, NO", 30.05, 69.73),
("Fredrikstad, NO", 10.93, 59.22),
("Gdańsk, PL", 18.65, 54.35),
("Riga, LV", 24.11, 56.95),
("Tallinn, EE", 24.75, 59.44),
("Klaipėda, LT", 21.14, 55.71),
("Rostock, DE", 12.14, 54.09),
];
/// How many of the 43 named places a polygon gets wrong.
fn probes_wrong(poly: &Polygon) -> (usize, Vec<String>) {
let mut wrong = Vec::new();
for &(n, lon, lat) in INSIDE {
if !poly.contains_exact(lon, lat) {
wrong.push(format!("{n} dropped OUT"));
}
}
for &(n, lon, lat) in OUTSIDE {
if poly.contains_exact(lon, lat) {
wrong.push(format!("{n} pulled IN"));
}
}
(wrong.len(), wrong)
}
/// Absolute area of all rings (shoelace, in square degrees) — a scale-free way to
/// say how much shape a simplification moved, with no external truth set needed.
fn area(poly: &Polygon) -> f64 {
poly.rings()
.iter()
.map(|r| {
let mut a = 0.0;
for w in r.windows(2) {
a += w[0].0 * w[1].1 - w[1].0 * w[0].1;
}
(a / 2.0).abs()
})
.sum()
}
/// **A debug build measures nothing here.** `opt-level = 0` leaves the bounds checks,
/// the divisions and the edge walk unoptimised, so every arm slows by roughly an order
/// of magnitude — but *unevenly*, which silently reorders the very comparison this bench
/// exists to make. The numbers would still print, and they would still look like results.
///
/// Refuse instead. A bench that cannot tell you it was run wrong belongs to the same
/// family as a guard that can never report red.
#[cfg(debug_assertions)]
fn main() {
eprintln!(
"clip_containment: REFUSING to run under debug_assertions — every figure this \
bench prints would be an artefact of opt-level=0, and the ARMS would not even \
slow down evenly. Run it as `cargo bench -p osm-katana --bench clip_containment` \
(the bench profile inherits release), or `cargo test --release` for the \
equivalence test next to it."
);
std::process::exit(2);
}
#[cfg(not(debug_assertions))]
fn main() {
let poly_path = std::env::var("OSM_KATANA_CLIP_POLY")
.unwrap_or_else(|_| format!("{}/git/korp/assets/sweden.poly", std::env::var("HOME").unwrap_or_default()));
let nodes_path = PathBuf::from(
std::env::var("OSM_KATANA_CLIP_NODES").unwrap_or_else(|_| DEFAULT_NODES.to_string()),
);
let stride: usize = std::env::var("OSM_KATANA_CLIP_STRIDE")
.ok()
.and_then(|s| s.parse().ok())
.unwrap_or(1);
let workers = std::thread::available_parallelism()
.map(|n| n.get())
.unwrap_or(2)
.saturating_sub(1);
let text = match std::fs::read_to_string(&poly_path) {
Ok(t) => t,
Err(e) => {
eprintln!("clip_containment: cannot read {poly_path}: {e}");
std::process::exit(1);
}
};
let full = Polygon::from_poly_str(&text).expect("the clip polygon parses");
let bb = full.bbox();
eprintln!(
"clip_containment: polygon {poly_path} — {} ring(s), {} vertices, bbox \
{:.4}..{:.4} / {:.4}..{:.4}",
full.ring_count(),
full.vertex_count(),
bb.min_lon,
bb.max_lon,
bb.min_lat,
bb.max_lat,
);
// ── Corpus: real node coordinates ────────────────────────────────────────
// Sweeping the planet's node section for a sample costs a full pass, so the
// sample is cached to `OSM_KATANA_CLIP_CORPUS` (raw little-endian f64 pairs)
// and every later arm reads it back. One sweep, many measurements.
let corpus_cache = std::env::var("OSM_KATANA_CLIP_CORPUS")
.ok()
.map(PathBuf::from);
let t0 = Instant::now();
let nodes: Vec<(f64, f64)> = match corpus_cache.as_ref().filter(|p| p.exists()) {
Some(p) => {
let raw = std::fs::read(p).expect("read cached corpus");
raw.chunks_exact(16)
.map(|c| {
(
f64::from_le_bytes(c[..8].try_into().expect("8 bytes")),
f64::from_le_bytes(c[8..].try_into().expect("8 bytes")),
)
})
.collect()
}
None => {
let sampled = osm_katana::xml_scan::sample_node_coords(&nodes_path, workers, stride)
.expect("sample node coordinates");
if let Some(p) = corpus_cache.as_ref() {
let mut buf = Vec::with_capacity(sampled.len() * 16);
for &(lon, lat) in &sampled {
buf.extend_from_slice(&lon.to_le_bytes());
buf.extend_from_slice(&lat.to_le_bytes());
}
std::fs::write(p, &buf).expect("write corpus cache");
eprintln!("clip_containment: cached corpus → {}", p.display());
}
sampled
}
};
let sample_s = t0.elapsed().as_secs_f64();
assert!(
!nodes.is_empty(),
"the node corpus is empty — nothing would be measured"
);
let in_bbox = nodes
.iter()
.filter(|&&(lon, lat)| bb.contains(lon, lat))
.count();
eprintln!(
"clip_containment: corpus {} — {} nodes sampled (stride {stride}) in {sample_s:.1} s; \
{in_bbox} ({:.3} %) land inside the polygon's bbox and pay for a ray-cast",
nodes_path.display(),
nodes.len(),
100.0 * in_bbox as f64 / nodes.len() as f64,
);
println!(
"\n{:<34} {:>9} {:>8} {:>12} {:>12} {:>12} {:>9} {:>12}",
"variant",
"vertices",
"wrong43",
"wall_s",
"ns_per_node",
"kept",
"area_err",
"nodes_moved",
);
let full_area = area(&full);
// (label, ns/node, nodes whose verdict moved) for the extrapolation table.
let mut rows: Vec<(String, f64, usize)> = Vec::new();
// ── Baseline: today's clip, exact ray-cast, no mask ──────────────────────
// Every other row's verdicts are compared against THIS one, node by node.
let t0 = Instant::now();
let truth: Vec<bool> = nodes
.iter()
.map(|&(lon, lat)| full.contains_exact(lon, lat))
.collect();
let baseline_s = t0.elapsed().as_secs_f64();
let baseline_ns = baseline_s * 1e9 / nodes.len() as f64;
let baseline_kept = truth.iter().filter(|&&k| k).count();
{
let (w, _) = probes_wrong(&full);
println!(
"{:<34} {:>9} {:>8} {:>12.3} {:>12.1} {:>12} {:>9} {:>12}",
"exact (today's clip)",
full.vertex_count(),
w,
baseline_s,
baseline_ns,
baseline_kept,
"-",
0,
);
rows.push(("exact (today's clip)".into(), baseline_ns, 0));
}
/// How many real nodes a variant admits or drops that the full-fidelity clip
/// does not — the measure that matters, because it is counted in the data the
/// clip actually processes rather than at 43 hand-picked points.
let moved = |verdicts: &[bool]| -> usize {
verdicts.iter().zip(&truth).filter(|(a, b)| a != b).count()
};
// ── 3. Mask, at several grid resolutions. Geometry untouched. ────────────
for steps in [256usize, 512, 1024, 2048] {
let t0 = Instant::now();
let mask = osm_katana::polygon_mask::PolygonMask::build(full.rings(), bb, steps);
let build_s = t0.elapsed().as_secs_f64();
let (bcell, icell, ocell) = mask.cell_counts();
let t1 = Instant::now();
let verdicts: Vec<bool> = nodes
.iter()
.map(|&(lon, lat)| {
mask.lookup(lon, lat)
.unwrap_or_else(|| full.contains_exact(lon, lat))
})
.collect();
let s = t1.elapsed().as_secs_f64();
let kept = verdicts.iter().filter(|&&k| k).count();
let ns = s * 1e9 / nodes.len() as f64;
let m = moved(&verdicts);
// A mask that moved a single node is a defect, not a trade-off.
assert_eq!(
m, 0,
"mask {steps}² disagreed with the exact ray-cast on {m} real nodes"
);
println!(
"{:<34} {:>9} {:>8} {:>12.3} {:>12.1} {:>12} {:>9} {:>12}",
format!(
"mask {steps}² (+{:.0} ms, {} KiB)",
build_s * 1e3,
mask.bytes() / 1024
),
full.vertex_count(),
0,
s,
ns,
kept,
"0",
m,
);
rows.push((format!("mask {steps}²"), ns, m));
eprintln!(
" mask {steps}²: {:.2} % boundary cells, {:.2} % inside, {:.2} % outside; \
{:.2}× vs exact",
100.0 * bcell as f64 / (steps * steps) as f64,
100.0 * icell as f64 / (steps * steps) as f64,
100.0 * ocell as f64 / (steps * steps) as f64,
baseline_ns / ns,
);
}
// ── CEILING: where do the mask's nanoseconds actually go? ────────────────
// Before anyone builds a cleverer index, decompose the masked path into the
// three things it does, over the same corpus:
// A bbox reject only — the floor every node pays, mask or not
// B bbox + cell lookup — A plus the array index (WRONG verdicts
// at the boundary; a timing probe only)
// C the real masked path — B plus the ray-cast for cells an edge
// crosses
// C − B is the entire budget any sparse/stree index could compete for. If it
// is a sliver, the answer to "can Ragnar's stree be the mask" is "there is
// nothing left to win", and that is a complete answer.
{
let steps = 1024usize;
let mask = osm_katana::polygon_mask::PolygonMask::build(full.rings(), bb, steps);
let t = Instant::now();
let mut acc = 0usize;
for &(lon, lat) in &nodes {
acc += usize::from(bb.contains(lon, lat));
}
let a_s = t.elapsed().as_secs_f64();
let in_bbox_count = acc;
let t = Instant::now();
let mut fell_through = 0usize;
let mut acc_b = 0usize;
for &(lon, lat) in &nodes {
match mask.lookup(lon, lat) {
Some(v) => acc_b += usize::from(v),
None => fell_through += 1,
}
}
let b_s = t.elapsed().as_secs_f64();
let t = Instant::now();
let mut acc_c = 0usize;
for &(lon, lat) in &nodes {
acc_c += usize::from(
mask.lookup(lon, lat)
.unwrap_or_else(|| full.contains_exact(lon, lat)),
);
}
let c_s = t.elapsed().as_secs_f64();
let n = nodes.len() as f64;
let (a_ns, b_ns, c_ns) = (a_s * 1e9 / n, b_s * 1e9 / n, c_s * 1e9 / n);
println!(
"\nCEILING — where the mask 1024² nanoseconds go ({} corpus nodes, acc {acc_b}/{acc_c}):\n\
{:<40} {:>10} {:>10}\n\
{:<40} {:>10.2} {:>9.1}%\n\
{:<40} {:>10.2} {:>9.1}%\n\
{:<40} {:>10.2} {:>9.1}%\n\
{:<40} {:>10.2} {:>9.1}%",
nodes.len(),
"component",
"ns/node",
"share",
"A bbox reject (floor)",
a_ns,
100.0 * a_ns / c_ns,
// NOT "the array index": this is everything calling the mask adds
// over a bare inlined bbox compare, paid by EVERY node — the call
// itself, the coordinate divisions, and the loop no longer
// vectorising. The array index proper is measured per in-bbox query
// in the STREE AS MASK table below, at ~2.3 ns, i.e. ~0.05 ns/node
// amortised. Mislabelling this as the lookup would send the next
// reader optimising the wrong term.
"B− A mask call, every node",
b_ns - a_ns,
100.0 * (b_ns - a_ns) / c_ns,
"C− B boundary ray-cast <- the budget",
c_ns - b_ns,
100.0 * (c_ns - b_ns) / c_ns,
"C total (today's masked clip)",
c_ns,
100.0,
);
println!(
" {fell_through} of {} nodes ({:.4} % of all, {:.3} % of the {in_bbox_count} in-bbox) \
reach a ray-cast; each costs ~{:.1} µs over {} vertices.\n \
A perfect index that removed EVERY boundary ray-cast could reach {:.2} ns/node \
= {:.1} s of planet clip wall, against {:.2} ns / {:.1} s today.",
nodes.len(),
100.0 * fell_through as f64 / n,
100.0 * fell_through as f64 / in_bbox_count.max(1) as f64,
(c_s - b_s) * 1e6 / fell_through.max(1) as f64,
full.vertex_count(),
b_ns,
b_ns * 10_581_367_185.0 / 1e9 / workers as f64,
c_ns,
c_ns * 10_581_367_185.0 / 1e9 / workers as f64,
);
}
// ── Can Ragnar's static search tree be the mask? ─────────────────────────
// The dense mask is a direct array index; a tree is a search. The only place
// a tree could win is the SPARSE form: at a much finer grid, hold just the
// BOUNDARY cells as sorted keys and ask "is this cell one of them", so fewer
// nodes reach a ray-cast without paying dense memory for the fine grid.
//
// Measured against the dense array on the SAME queries — the in-bbox nodes,
// which are the only ones that ever reach the mask at all. Everything else is
// rejected by the bbox test both ways, so amortising over the whole corpus
// would hide the comparison rather than make it.
if std::env::var("OSM_KATANA_CLIP_STREE").is_ok() {
use znippy_zoomies::stree::STree64;
let in_bbox: Vec<(f64, f64)> = nodes
.iter()
.copied()
.filter(|&(lon, lat)| bb.contains(lon, lat))
.collect();
println!(
"\nSTREE AS MASK — {} in-bbox queries (the only nodes that reach a mask)\n\
{:<18} {:>10} {:>12} {:>11} {:>12} {:>12} {:>12}",
in_bbox.len(),
"grid",
"bnd cells",
"index bytes",
"build_ms",
"dense ns/q",
"stree ns/q",
"batch ns/q",
);
for steps in [1024usize, 4096, 16384] {
let t = Instant::now();
let mask = osm_katana::polygon_mask::PolygonMask::build(full.rings(), bb, steps);
let dense_build_ms = t.elapsed().as_secs_f64() * 1e3;
let keys = mask.boundary_cells();
if keys.is_empty() {
continue;
}
// The queries, precomputed once so neither arm pays for the divisions.
let q: Vec<i64> = in_bbox
.iter()
.filter_map(|&(lon, lat)| mask.cell_key(lon, lat))
.collect();
// Dense arm: one array index per query, on the SAME precomputed keys
// the stree arm gets, so neither pays for the other's arithmetic.
// `dense_in` is printed — a loop whose result is unused is deleted by
// the compiler and reports a speed that never happened.
let t = Instant::now();
let mut dense_in = 0usize;
let mut dense_bnd = 0usize;
for &k in &q {
match mask.lookup_key(k) {
Some(true) => dense_in += 1,
Some(false) => {}
None => dense_bnd += 1,
}
}
let dense_ns = t.elapsed().as_secs_f64() * 1e9 / q.len() as f64;
// stree arm: build over the sorted boundary keys, then query.
let t = Instant::now();
let tree = STree64::new(&keys);
let stree_build_ms = t.elapsed().as_secs_f64() * 1e3;
let t = Instant::now();
let mut hits = 0usize;
for &k in &q {
hits += usize::from(tree.find_exact(k).is_some());
}
let stree_ns = t.elapsed().as_secs_f64() * 1e9 / q.len() as f64;
// The batch path — where the software pipelining and prefetching pay,
// and the shape a clip actually has (millions of nodes in bulk).
let t = Instant::now();
let mut bhits = 0usize;
for chunk in q.chunks(16) {
if chunk.len() == 16 {
let mut arr = [0i64; 16];
arr.copy_from_slice(chunk);
for r in tree.batch_prefetch::<16>(&arr) {
bhits += usize::from(r.is_some());
}
} else {
for &k in chunk {
bhits += usize::from(tree.find_exact(k).is_some());
}
}
}
let batch_ns = t.elapsed().as_secs_f64() * 1e9 / q.len() as f64;
assert_eq!(hits, bhits, "batch and scalar stree disagree on membership");
println!(
"{:<18} {:>10} {:>12} {:>11} {:>12.1} {:>12.1} {:>12.1}",
format!("{steps}²"),
keys.len(),
keys.len() * 8,
format!("{dense_build_ms:.0}/{stree_build_ms:.0}"),
dense_ns,
stree_ns,
batch_ns,
);
assert_eq!(
dense_bnd, hits,
"the stree membership set and the dense mask disagree about which \
queries are on the boundary",
);
eprintln!(
" {steps}²: dense mask {} MiB vs stree keys {:.2} MiB; dense_in {dense_in}; \
{hits} of {} queries land in a boundary cell ({:.3} %)",
mask.bytes() / (1024 * 1024),
keys.len() as f64 * 8.0 / (1024.0 * 1024.0),
q.len(),
100.0 * hits as f64 / q.len() as f64,
);
}
}
// ── 1. Simplification sweep. Geometry CHANGES — probe every step. ────────
for eps in [0.0001_f64, 0.0005, 0.001, 0.005, 0.01, 0.05] {
let simp = match full.simplify(eps) {
Ok(p) => p,
Err(e) => {
println!("{:<34} simplify({eps}) failed: {e}", "");
continue;
}
};
let (w, detail) = probes_wrong(&simp);
let t1 = Instant::now();
let verdicts: Vec<bool> = nodes
.iter()
.map(|&(lon, lat)| simp.contains_exact(lon, lat))
.collect();
let s = t1.elapsed().as_secs_f64();
let kept = verdicts.iter().filter(|&&k| k).count();
let ns = s * 1e9 / nodes.len() as f64;
let area_err = 100.0 * (area(&simp) - full_area).abs() / full_area;
println!(
"{:<34} {:>9} {:>8} {:>12.3} {:>12.1} {:>12} {:>8.3}% {:>12}",
format!("simplified eps={eps}"),
simp.vertex_count(),
w,
s,
ns,
kept,
area_err,
moved(&verdicts),
);
rows.push((
format!("simplified eps={eps} ({} v)", simp.vertex_count()),
ns,
moved(&verdicts),
));
if !detail.is_empty() {
eprintln!(" eps={eps} BREAKS: {}", detail.join("; "));
}
}
// ── 1+3 composed: simplify AND mask, at the largest eps that stays 43/43 ──
for eps in [0.0001_f64, 0.001] {
if let Ok(simp) = full.simplify(eps) {
let (w, _) = probes_wrong(&simp);
let mask =
osm_katana::polygon_mask::PolygonMask::build(simp.rings(), simp.bbox(), 1024);
let t1 = Instant::now();
let verdicts: Vec<bool> = nodes
.iter()
.map(|&(lon, lat)| {
mask.lookup(lon, lat)
.unwrap_or_else(|| simp.contains_exact(lon, lat))
})
.collect();
let s = t1.elapsed().as_secs_f64();
let kept = verdicts.iter().filter(|&&k| k).count();
let ns = s * 1e9 / nodes.len() as f64;
println!(
"{:<34} {:>9} {:>8} {:>12.3} {:>12.1} {:>12} {:>9} {:>12}",
format!("simplified eps={eps} + mask 1024²"),
simp.vertex_count(),
w,
s,
ns,
kept,
"-",
moved(&verdicts),
);
rows.push((format!("simplified eps={eps} + mask"), ns, moved(&verdicts)));
}
}
// ── Extrapolation to a planet clip ───────────────────────────────────────
// 10 581 367 185 nodes is this planet's measured node count (bench_history,
// 2026-06-18). The clip test runs inside the gatling pass-1 workers, so the
// wall-clock share is the per-node cost times the node count over the worker
// count — stated as an extrapolation, not as a measured planet clip.
const PLANET_NODES: f64 = 10_581_367_185.0;
println!(
"\nExtrapolated to a {PLANET_NODES:.0}-node planet across {workers} gatling workers.\n\
This is ARITHMETIC on the per-node numbers above — not a measured re-cut, which is\n\
a deliberate non-goal here. `kept` above scales to {:.0} M nodes, against the\n\
103.5 M a real Sweden clip produced, so the sampled corpus is representative.\n\
{:<34} {:>12} {:>12} {:>10} {:>12}",
baseline_kept as f64 * PLANET_NODES / nodes.len() as f64 / 1e6,
"variant",
"clip_cpu_s",
"clip_wall_s",
"vs_today",
"nodes_moved",
);
for (label, ns, moved_n) in &rows {
let cpu = ns * PLANET_NODES / 1e9;
println!(
"{:<34} {:>12.0} {:>12.1} {:>9.1}× {:>12}",
label,
cpu,
cpu / workers as f64,
baseline_ns / ns,
moved_n,
);
}
}