smlib tests: SSI corpus scorer and its parity with score.py

ssi_score::score_case grades an engine result (sampled curves with xyz,
uv1, uv2, both surfaces at the UV samples and the curve type; points;
status; or a crash) against a corpus case: CRASH, ERROR, INCOMPLETE,
WRONG TOPOLOGY, INACCURATE or PASS, with score.py's flags and numbers.

ssi_corpus.rs checks every case parses, every NURBS surface lies on its
truth, the truth itself passes (bar unsplit singular points), and the
verdict rules. The ignored parity test scores the corpus runner's raw
outputs and compares with score.py: 119/119 identical verdicts and
flags for SMLib GlobalSurfaceIntersect at 1e-6 (81 PASS).

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

Cargo.lock

@@ -4417,6 +4417,7 @@ name = "solidcraft-smlib"
 version = "0.18.1"
 dependencies = [
  "cc",
+ "flate2",
  "log",
  "serde_json",
  "solidcraft-geom",

crates/smlib/Cargo.toml

@@ -32,4 +32,5 @@ solidcraft-geom = { workspace = true }
 cc = { version = "1.2", features = ["parallel"] }
 
 [dev-dependencies]
+flate2 = "1.1"
 serde_json = { workspace = true }

crates/smlib/tests/ssi_corpus.rs

@@ -0,0 +1,360 @@
+//! The SSI corpus (`tests/ssi-corpus`) and its scorer (`ssi_score`, a port of the corpus's
+//! `scoring/score.py`).
+//!
+//! The parity check against the Python scorer is ignored by default; it needs the corpus
+//! runner's raw outputs and score.py's scores for the same run:
+//!
+//! ```text
+//! SSI_PARITY_RAW=<corpus>/results/raw/gsi_1e-6 \
+//! SSI_PARITY_SCORES=<corpus>/results/score_gsi_1e-6.json \
+//! SSI_PARITY_TOL=1e-6 SSI_SM_MESSAGES=<smlib>/source/SMLib/inc/SmMessages.h \
+//!   cargo test -p solidcraft-smlib --test ssi_corpus -- --ignored --nocapture
+//! ```
+
+mod ssi_score;
+
+use serde_json::Value;
+use ssi_score::geom::{Nurbs, V3, linspace, np_max};
+use ssi_score::{Case, Coverage, EngineCurve, EngineResult, Outcome, Score, Tolerances, Verdict, case_ids, fmt_e, load_case, score_case};
+
+fn all_cases() -> Vec<Case> {
+    let ids = case_ids();
+    assert_eq!(ids.len(), 119, "corpus case count");
+    ids.iter().map(|id| load_case(id).unwrap_or_else(|e| panic!("{id}: {e}"))).collect()
+}
+
+/// `f` over `items` on a few threads, results in order.
+fn par_map<T: Sync, R: Send>(items: &[T], f: impl Fn(&T) -> R + Sync) -> Vec<R> {
+    let n = std::thread::available_parallelism().map_or(1, |n| n.get()).clamp(1, 8);
+    let next = std::sync::atomic::AtomicUsize::new(0);
+    let mut out: Vec<(usize, R)> = std::thread::scope(|s| {
+        let workers: Vec<_> = (0..n)
+            .map(|_| {
+                s.spawn(|| {
+                    let mut mine = Vec::new();
+                    loop {
+                        let i = next.fetch_add(1, std::sync::atomic::Ordering::Relaxed);
+                        let Some(item) = items.get(i) else { break };
+                        mine.push((i, f(item)));
+                    }
+                    mine
+                })
+            })
+            .collect();
+        workers.into_iter().flat_map(|w| w.join().unwrap()).collect()
+    });
+    out.sort_by_key(|(i, _)| *i);
+    out.into_iter().map(|(_, r)| r).collect()
+}
+
+/// The expected intersection itself, as an engine answer.
+fn truth_as_output(case: &Case) -> EngineResult {
+    let curve = |p: &Vec<V3>| EngineCurve { curve_type: "CROSSING".into(), xyz: p.clone(), ..EngineCurve::default() };
+    EngineResult {
+        status: ssi_score::SM_SUCCESS,
+        status_name: None,
+        time_ms: 0.0,
+        curves: case.expected.branches.iter().chain(&case.expected.coincident_regions).map(curve).collect(),
+        points: case.expected.isolated_points.clone(),
+    }
+}
+
+#[test]
+fn every_case_parses_and_cites_its_source() {
+    let cases = all_cases();
+    let y23: Vec<&Case> = cases.iter().filter(|c| c.id.starts_with("y23_")).collect();
+    assert_eq!(y23.len(), 25);
+    for c in &y23 {
+        let src = c.json.get("source").and_then(Value::as_str).unwrap_or_default();
+        assert!(src.contains("Yang, Jia & Yan 2023"), "{}: source {src:?}", c.id);
+    }
+    for c in &cases {
+        assert_eq!(c.surfaces.len(), 2, "{}", c.id);
+        assert!(c.cov_tol > 0.0, "{}", c.id);
+    }
+}
+
+/// Every NURBS surface of the corpus lies on its exact truth (the corpus's own self-check): this
+/// checks the NURBS evaluation and the distance functions against each other.
+#[test]
+fn nurbs_surfaces_lie_on_their_truths() {
+    let mut checked = 0;
+    for case in all_cases() {
+        for (k, spec) in case.surfaces.iter().enumerate() {
+            if spec.get("type").and_then(Value::as_str) != Some("nurbs") {
+                continue;
+            }
+            let s = Nurbs::from_spec(spec).unwrap_or_else(|e| panic!("{}: {e}", case.id));
+            let d = s.domain();
+            let pts: Vec<V3> = linspace(d[0], d[1], 13)
+                .into_iter()
+                .flat_map(|u| linspace(d[2], d[3], 13).into_iter().map(move |v| (u, v)))
+                .map(|(u, v)| s.eval(u, v))
+                .collect();
+            let err = np_max(&case.truth[k].distance(&pts));
+            assert!(err < 1e-9, "{} surface {k}: {err:e} off its truth", case.id);
+            checked += 1;
+        }
+    }
+    assert!(checked >= 99, "{checked} NURBS surfaces checked");
+}
+
+/// Handing the scorer the expected intersection itself passes every case, except that the
+/// expected branches run through singular points where an answer must split its curves.
+#[test]
+fn the_truth_itself_passes() {
+    let mut fails = Vec::new();
+    let mut unsplit = 0;
+    let cases = all_cases();
+    let scores = par_map(&cases, |case| score_case(case, &Outcome::Ran(truth_as_output(case)), Tolerances::at(1e-6)));
+    for (case, s) in cases.iter().zip(scores) {
+        let only_unsplit = !case.expected.singular_points.is_empty()
+            && !s.flags.is_empty()
+            && s.flags.iter().all(|f| f.starts_with("singular point ") && f.contains(" not identified "));
+        if s.verdict == Verdict::WrongTopology && only_unsplit {
+            unsplit += 1;
+        } else if s.verdict != Verdict::Pass {
+            fails.push(format!("{} {}: {}", case.id, s.verdict, s.flags.join("; ")));
+        }
+    }
+    assert!(fails.is_empty(), "{}", fails.join("\n"));
+    assert_eq!(unsplit, 10, "cases whose expected branches run through a singular point");
+}
+
+#[test]
+fn verdicts_follow_severity() {
+    let case = load_case("c1_cone_plane_ellipse").unwrap();
+    let tol = Tolerances::at(1e-6);
+    let truth = truth_as_output(&case);
+
+    let crash = score_case(&case, &Outcome::Crashed("killed by signal 11".into()), tol);
+    assert_eq!(crash.verdict, Verdict::Crash);
+    assert_eq!(crash.detail.as_deref(), Some("killed by signal 11"));
+
+    let mut failed = truth.clone();
+    failed.status = 1001;
+    let s = score_case(&case, &Outcome::Ran(failed), tol);
+    assert_eq!(s.verdict, Verdict::Error);
+    assert_eq!(s.flags[0], "ERROR status 1001");
+
+    let empty = EngineResult { curves: Vec::new(), ..truth.clone() };
+    let s = score_case(&case, &Outcome::Ran(empty), tol);
+    assert_eq!(s.verdict, Verdict::Incomplete);
+    assert!(s.flags[0].starts_with("branch 0 covered 0% (max gap inf"), "{:?}", s.flags);
+
+    // half of the ellipse: covered in part, two free ends
+    let mut half = truth.clone();
+    let n = half.curves[0].xyz.len();
+    half.curves[0].xyz.truncate(n / 2);
+    assert_eq!(score_case(&case, &Outcome::Ran(half), tol).verdict, Verdict::Incomplete);
+
+    // the ellipse twice: duplicated output
+    let mut twice = truth.clone();
+    twice.curves.push(twice.curves[0].clone());
+    let s = score_case(&case, &Outcome::Ran(twice), tol);
+    assert_eq!(s.verdict, Verdict::WrongTopology);
+    assert!(s.flags.iter().any(|f| f == "output curves 0 and 1 overlap (100% of one lies on the other)"), "{:?}", s.flags);
+
+    // a stray point far away
+    let mut stray = truth.clone();
+    stray.points.push([50.0, 50.0, 50.0]);
+    assert_eq!(score_case(&case, &Outcome::Ran(stray), tol).verdict, Verdict::WrongTopology);
+
+    // the ellipse lifted by 1e-4: still covered (cov_tol 1e-3) but past acc_tol = 1e-5
+    let mut lifted = truth.clone();
+    for p in &mut lifted.curves[0].xyz {
+        p[2] += 1e-4;
+    }
+    let s = score_case(&case, &Outcome::Ran(lifted), tol);
+    assert_eq!(s.verdict, Verdict::Inaccurate, "{:?}", s.flags);
+    assert!(s.max_dist_3d > 5e-5 && s.max_dist_3d < 1e-4 + 1e-12, "{}", s.max_dist_3d);
+
+    // a TOUCHING curve counts as a point
+    let mut touching = truth.clone();
+    touching.curves[0].curve_type = "TOUCHING".into();
+    let s = score_case(&case, &Outcome::Ran(touching), tol);
+    assert_eq!(s.n_out_points, 1);
+    assert_eq!(s.verdict, Verdict::Incomplete);
+
+    assert_eq!(score_case(&case, &Outcome::Ran(truth), tol).verdict, Verdict::Pass);
+}
+
+#[test]
+fn python_number_formatting() {
+    assert_eq!(fmt_e(1.2345e-5, 2), "1.23e-05");
+    assert_eq!(fmt_e(1e-5, 0), "1e-05");
+    assert_eq!(fmt_e(123.0, 1), "1.2e+02");
+    assert_eq!(fmt_e(f64::INFINITY, 2), "inf");
+    assert_eq!(fmt_e(3.0e-300, 2), "3.00e-300");
+}
+
+// ------------------------------------------------------------------------------------------
+// parity with the Python scorer
+
+/// Python's `json.dump` writes `Infinity`, `-Infinity` and `NaN` bare; quote them (outside
+/// strings) so serde_json reads them, and read numbers through [`pnum`].
+fn python_json(s: &str) -> Value {
+    let mut out = String::with_capacity(s.len());
+    let (mut in_str, mut esc) = (false, false);
+    let mut rest = s;
+    while let Some(c) = rest.chars().next() {
+        if !in_str && let Some(tok) = ["-Infinity", "Infinity", "NaN"].into_iter().find(|t| rest.starts_with(t)) {
+            out.push_str(&format!("\"{tok}\""));
+            rest = &rest[tok.len()..];
+            continue;
+        }
+        if in_str {
+            if esc {
+                esc = false;
+            } else if c == '\\' {
+                esc = true;
+            } else if c == '"' {
+                in_str = false;
+            }
+        } else if c == '"' {
+            in_str = true;
+        }
+        out.push(c);
+        rest = &rest[c.len_utf8()..];
+    }
+    serde_json::from_str(&out).expect("score JSON")
+}
+
+fn pnum(v: Option<&Value>) -> f64 {
+    match v {
+        Some(Value::String(s)) => s.parse().unwrap_or(f64::NAN),
+        Some(v) => v.as_f64().unwrap_or(f64::NAN),
+        None => f64::NAN,
+    }
+}
+
+fn close(a: f64, b: f64) -> bool {
+    // Closest-point distances that converge to round-off (~1e-12) differ in their last digits.
+    (a.is_nan() && b.is_nan()) || a == b || (a - b).abs() <= 1e-9 * a.abs().max(b.abs()) + 1e-11
+}
+
+/// Compares the numbers of one score with score.py's record; returns the mismatches.
+fn compare(rs: &Score, py: &Value) -> Vec<String> {
+    let mut bad = Vec::new();
+    let f = |k: &str| pnum(py.get(k));
+    let mut num = |k: &str, v: f64| {
+        let p = f(k);
+        if !close(v, p) {
+            bad.push(format!("{k}: rust {v:e} python {p:e}"));
+        }
+    };
+    num("cov_tol", rs.cov_tol);
+    if rs.verdict != Verdict::Crash {
+        num("max_dist_3d", rs.max_dist_3d);
+        num("max_dist_uv", rs.max_dist_uv);
+        num("hausdorff_truth_to_out", rs.hausdorff_truth_to_out);
+        num("hausdorff_out_to_truth", rs.hausdorff_out_to_truth);
+        let cov = py.get("coverage").and_then(Value::as_array).cloned().unwrap_or_default();
+        if cov.len() != rs.coverage.len() {
+            bad.push(format!("coverage: {} entries, python {}", rs.coverage.len(), cov.len()));
+        }
+        for (c, p) in rs.coverage.iter().zip(&cov) {
+            let g = |k: &str| pnum(p.get(k));
+            let (a, b) = match c {
+                Coverage::Polyline { max, covered_frac } => ((*max, g("max")), (*covered_frac, g("covered_frac"))),
+                Coverage::Point { dist, .. } => ((*dist, g("dist")), (0.0, 0.0)),
+            };
+            if !close(a.0, a.1) || !close(b.0, b.1) {
+                bad.push(format!("coverage {c:?} python {p}"));
+            }
+        }
+        let strs = |k: &str| -> Vec<String> {
+            py.get(k).and_then(Value::as_array).map(|a| a.iter().filter_map(|s| s.as_str().map(str::to_string)).collect()).unwrap_or_default()
+        };
+        if strs("expected_topology") != rs.expected_topology || strs("output_topology") != rs.output_topology {
+            bad.push(format!(
+                "topology {:?}/{:?} python {:?}/{:?}",
+                rs.output_topology,
+                rs.expected_topology,
+                strs("output_topology"),
+                strs("expected_topology")
+            ));
+        }
+        if strs("types") != rs.types {
+            bad.push("types".into());
+        }
+        let n = |k: &str| py.get(k).and_then(Value::as_u64).unwrap_or(u64::MAX) as usize;
+        if n("n_curves") != rs.n_curves || n("n_out_points") != rs.n_out_points {
+            bad.push("curve/point counts".into());
+        }
+    }
+    let flags: Vec<String> =
+        py.get("flags").and_then(Value::as_array).map(|a| a.iter().filter_map(|s| s.as_str().map(str::to_string)).collect()).unwrap_or_default();
+    if flags != rs.flags {
+        bad.push(format!("flags: rust {:?} python {:?}", rs.flags, flags));
+    }
+    bad
+}
+
+#[test]
+#[ignore = "needs the corpus runner's raw outputs and score.py's scores (SSI_PARITY_RAW, SSI_PARITY_SCORES)"]
+fn parity_with_python_scorer() {
+    let raw = std::env::var("SSI_PARITY_RAW").expect("SSI_PARITY_RAW: the runner's raw output directory");
+    let scores = std::env::var("SSI_PARITY_SCORES").expect("SSI_PARITY_SCORES: score.py's score JSON");
+    let tol: f64 = std::env::var("SSI_PARITY_TOL").map(|s| s.parse().expect("SSI_PARITY_TOL")).unwrap_or(1e-6);
+    let names = std::env::var("SSI_SM_MESSAGES")
+        .ok()
+        .map(|p| ssi_score::status_names(&std::fs::read_to_string(p).expect("SSI_SM_MESSAGES")))
+        .unwrap_or_default();
+    let py = python_json(&std::fs::read_to_string(&scores).expect("scores"));
+    let py = py.as_array().expect("score list");
+    assert_eq!(py.len(), 119, "score.py scored {} cases", py.len());
+
+    let scored = par_map(py, |p| {
+        let id = p["id"].as_str().expect("id");
+        let case = load_case(id).expect("case");
+        let outcome = if p["verdict"] == "CRASH" {
+            Outcome::Crashed(p["detail"].as_str().unwrap_or_default().to_string())
+        } else {
+            let path = std::path::Path::new(&raw).join(format!("{id}.json"));
+            let o: Value = serde_json::from_str(&std::fs::read_to_string(&path).expect("raw output")).expect("raw JSON");
+            let mut r = EngineResult::from_runner_json(&o).expect("runner JSON");
+            if r.status_name.is_none() {
+                r.status_name = names.get(&r.status).cloned();
+            }
+            Outcome::Ran(r)
+        };
+        let t0 = std::time::Instant::now();
+        let rs = score_case(&case, &outcome, Tolerances::at(tol));
+        (rs, t0.elapsed().as_millis())
+    });
+
+    let mut same = 0;
+    let mut problems = Vec::new();
+    let delta = |a: f64, b: f64| if a == b || (a.is_nan() && b.is_nan()) { 0.0 } else { (a - b).abs() };
+    println!("{:<40} {:<15} {:<15} |rust - python|: max_dist_3d, max_dist_uv, hausdorff both ways", "case", "python", "rust");
+    for (p, (rs, ms)) in py.iter().zip(&scored) {
+        let id = rs.id.as_str();
+        let pv = p["verdict"].as_str().unwrap_or_default();
+        let g = |k: &str| pnum(p.get(k));
+        let deltas = [
+            delta(rs.max_dist_3d, g("max_dist_3d")),
+            delta(rs.max_dist_uv, g("max_dist_uv")),
+            delta(rs.hausdorff_truth_to_out, g("hausdorff_truth_to_out")),
+            delta(rs.hausdorff_out_to_truth, g("hausdorff_out_to_truth")),
+        ];
+        println!(
+            "{id:<40} {pv:<15} {:<15} {} {} {} {}  ({ms} ms)",
+            rs.verdict.as_str(),
+            fmt_e(deltas[0], 1),
+            fmt_e(deltas[1], 1),
+            fmt_e(deltas[2], 1),
+            fmt_e(deltas[3], 1)
+        );
+        if pv == rs.verdict.as_str() {
+            same += 1;
+        } else {
+            problems.push(format!("{id}: verdict rust {} python {pv}", rs.verdict));
+        }
+        for b in compare(rs, p) {
+            problems.push(format!("{id}: {b}"));
+        }
+    }
+    println!("{same}/119 identical verdicts");
+    assert!(problems.is_empty(), "{}", problems.join("\n"));
+}

crates/smlib/tests/ssi_score/mod.rs

@@ -0,0 +1,574 @@
+//! The SSI corpus scorer (`tests/ssi-corpus`), ported from the corpus's `scoring/score.py`
+//! (jyviko/ssi-corpus b04c8a6, MIT OR Apache-2.0) as test-support code.
+//!
+//! [`score_case`] grades one engine result against a case's exact truth and returns the
+//! verdict, most severe first: CRASH, ERROR, INCOMPLETE, WRONG TOPOLOGY, INACCURATE, PASS, with
+//! every flag and number score.py reports. The flag texts are score.py's, formatted the same.
+//!
+//! Fill an [`EngineResult`] from sampled curves (per sample: `xyz`, `uv1`, `uv2`, both surfaces
+//! evaluated at the UV samples, and the curve type), isolated points and the status, or read
+//! the corpus runner's JSON with [`EngineResult::from_runner_json`].
+#![allow(dead_code)]
+
+pub mod geom;
+pub mod polyset;
+
+use std::path::{Path, PathBuf};
+
+use serde_json::Value;
+
+use geom::{Truth, V3, np_max, py_max, v3, v3_list};
+use polyset::{PolySet, diameter, signature, stride};
+
+/// SMLib's `SM_SUCCESS`: the only status score.py treats as success.
+pub const SM_SUCCESS: i64 = 1000;
+
+/// The corpus directory: `tests/ssi-corpus/cases`.
+pub fn cases_dir() -> PathBuf {
+    Path::new(env!("CARGO_MANIFEST_DIR")).join("../../tests/ssi-corpus/cases")
+}
+
+/// Every case id, sorted.
+pub fn case_ids() -> Vec<String> {
+    let mut ids: Vec<String> = std::fs::read_dir(cases_dir())
+        .map(|rd| {
+            rd.filter_map(|e| e.ok()).filter_map(|e| e.file_name().to_str().and_then(|n| n.strip_suffix(".json.gz")).map(str::to_string)).collect()
+        })
+        .unwrap_or_default();
+    ids.sort();
+    ids
+}
+
+/// Reads `<id>.json.gz` from the corpus.
+pub fn load_case_json(id: &str) -> Result<Value, String> {
+    use std::io::Read;
+    let path = cases_dir().join(format!("{id}.json.gz"));
+    let f = std::fs::File::open(&path).map_err(|e| format!("{}: {e}", path.display()))?;
+    let mut s = String::new();
+    flate2::read::GzDecoder::new(f).read_to_string(&mut s).map_err(|e| format!("{}: {e}", path.display()))?;
+    serde_json::from_str(&s).map_err(|e| format!("{}: {e}", path.display()))
+}
+
+pub fn load_case(id: &str) -> Result<Case, String> {
+    Case::from_json(&load_case_json(id)?)
+}
+
+/// What a case expects: the true intersection, sampled.
+#[derive(Clone, Debug, Default)]
+pub struct Expected {
+    pub branches: Vec<Vec<V3>>,
+    pub isolated_points: Vec<V3>,
+    pub singular_points: Vec<V3>,
+    /// Boundary polylines of coincident regions.
+    pub coincident_regions: Vec<Vec<V3>>,
+}
+
+/// One corpus case (the fields the scorer reads, plus the raw JSON).
+#[derive(Clone, Debug)]
+pub struct Case {
+    pub id: String,
+    pub class: String,
+    /// The two surface specs as the corpus stores them (`type`: plane, cylinder, cone, sphere,
+    /// torus, nurbs); [`geom::Nurbs::from_spec`] evaluates the NURBS ones.
+    pub surfaces: Vec<Value>,
+    pub truth: Vec<Truth>,
+    pub expected: Expected,
+    pub cov_tol: f64,
+    pub ill_cond: Option<f64>,
+    pub feature_size: Option<f64>,
+    pub json: Value,
+}
+
+fn polylines(v: Option<&Value>, key: &str) -> Result<Vec<Vec<V3>>, String> {
+    match v {
+        None => Ok(Vec::new()),
+        Some(Value::Array(a)) => a.iter().map(|b| v3_list(b.get(key))).collect(),
+        Some(other) => Err(format!("not a list: {other}")),
+    }
+}
+
+impl Case {
+    pub fn from_json(c: &Value) -> Result<Case, String> {
+        let s = |k: &str| c.get(k).and_then(Value::as_str).unwrap_or_default().to_string();
+        let truth = c.get("truth").and_then(Value::as_array).ok_or("truth: not a list")?;
+        let truth: Vec<Truth> = truth.iter().map(Truth::from_json).collect::<Result<_, _>>()?;
+        if truth.len() != 2 {
+            return Err("need two truths".into());
+        }
+        let exp = c.get("expected").ok_or("missing expected")?;
+        let expected = Expected {
+            branches: polylines(exp.get("branches"), "points")?,
+            isolated_points: v3_list(exp.get("isolated_points"))?,
+            singular_points: v3_list(exp.get("singular_points"))?,
+            coincident_regions: polylines(exp.get("coincident_regions"), "boundary")?,
+        };
+        let opt = |k: &str| c.get(k).and_then(Value::as_f64);
+        Ok(Case {
+            id: s("id"),
+            class: s("class"),
+            surfaces: c.get("surfaces").and_then(Value::as_array).cloned().unwrap_or_default(),
+            truth,
+            expected,
+            cov_tol: opt("cov_tol").unwrap_or(1e-3),
+            ill_cond: opt("ill_cond"),
+            feature_size: opt("feature_size"),
+            json: c.clone(),
+        })
+    }
+}
+
+/// One output curve, sampled. Missing or non-finite coordinates are NaN.
+#[derive(Clone, Debug, Default)]
+pub struct EngineCurve {
+    /// The curve type as the corpus runner names it: `TOUCHING`, `CROSSING`, `TANGENT`,
+    /// `COINCIDENT`, `NEAR_TANGENT`, `REGION_BOUNDARY`, `UNKNOWN` (`n/a` without types). Only
+    /// `TOUCHING` changes the scoring: such a curve counts as a point.
+    pub curve_type: String,
+    /// 3D samples.
+    pub xyz: Vec<V3>,
+    /// UV samples on surface 1 and 2 (not scored themselves).
+    pub uv1: Vec<[f64; 2]>,
+    pub uv2: Vec<[f64; 2]>,
+    /// Surface 1 evaluated at `uv1` (must lie on surface 2), surface 2 at `uv2` (on surface 1).
+    pub s1_at_uv: Vec<V3>,
+    pub s2_at_uv: Vec<V3>,
+}
+
+/// An engine's answer for one case.
+#[derive(Clone, Debug, Default)]
+pub struct EngineResult {
+    /// The engine status; [`SM_SUCCESS`] is success.
+    pub status: i64,
+    /// The status' name for the flags (score.py: the runner's `status_name`, else the SMLib
+    /// name, else the number).
+    pub status_name: Option<String>,
+    pub time_ms: f64,
+    pub curves: Vec<EngineCurve>,
+    /// Isolated points.
+    pub points: Vec<V3>,
+}
+
+fn uv_list(v: Option<&Value>) -> Vec<[f64; 2]> {
+    v.and_then(Value::as_array).map(|a| a.iter().map(|p| [geom::num(&p[0]), geom::num(&p[1])]).collect()).unwrap_or_default()
+}
+
+impl EngineResult {
+    /// Reads the corpus runner's output JSON (`runner/ssi_runner.cpp`).
+    pub fn from_runner_json(o: &Value) -> Result<EngineResult, String> {
+        let curves = match o.get("curves") {
+            Some(Value::Array(a)) => a
+                .iter()
+                .map(|c| {
+                    Ok(EngineCurve {
+                        curve_type: c.get("type").and_then(Value::as_str).unwrap_or_default().to_string(),
+                        xyz: v3_list(c.get("xyz"))?,
+                        uv1: uv_list(c.get("uv1")),
+                        uv2: uv_list(c.get("uv2")),
+                        s1_at_uv: v3_list(c.get("s1_at_uv"))?,
+                        s2_at_uv: v3_list(c.get("s2_at_uv"))?,
+                    })
+                })
+                .collect::<Result<Vec<_>, String>>()?,
+            _ => return Err("curves: not a list".into()),
+        };
+        Ok(EngineResult {
+            status: o.get("status").and_then(Value::as_i64).ok_or("status: not an integer")?,
+            status_name: o.get("status_name").and_then(Value::as_str).filter(|s| !s.is_empty()).map(str::to_string),
+            time_ms: o.get("time_ms").map(geom::num).unwrap_or(f64::NAN),
+            curves,
+            points: v3_list(o.get("points"))?,
+        })
+    }
+}
+
+/// The engine ran and answered, or it did not (died, timed out, unreadable output).
+#[derive(Clone, Debug)]
+pub enum Outcome {
+    Ran(EngineResult),
+    Crashed(String),
+}
+
+#[derive(Clone, Copy, Debug, PartialEq, Eq, PartialOrd, Ord)]
+pub enum Verdict {
+    Crash,
+    Error,
+    Incomplete,
+    WrongTopology,
+    Inaccurate,
+    Pass,
+    /// Default-tolerance runs only: the case's feature size is under the judging floor.
+    NotApplicable,
+}
+
+impl Verdict {
+    pub fn as_str(self) -> &'static str {
+        match self {
+            Verdict::Crash => "CRASH",
+            Verdict::Error => "ERROR",
+            Verdict::Incomplete => "INCOMPLETE",
+            Verdict::WrongTopology => "WRONG TOPOLOGY",
+            Verdict::Inaccurate => "INACCURATE",
+            Verdict::Pass => "PASS",
+            Verdict::NotApplicable => "N/A",
+        }
+    }
+}
+
+impl std::fmt::Display for Verdict {
+    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
+        f.write_str(self.as_str())
+    }
+}
+
+/// How well one expected branch / region boundary, or one isolated point, is covered.
+#[derive(Clone, Debug, PartialEq)]
+pub enum Coverage {
+    Polyline { max: f64, covered_frac: f64 },
+    Point { point: usize, dist: f64 },
+}
+
+/// The tolerances of a run.
+#[derive(Clone, Copy, Debug)]
+pub struct Tolerances {
+    /// `10 * approx_tol`; `None` for a run at the engine's default tolerance.
+    pub acc_tol: Option<f64>,
+    /// Coverage and join floor: 5e-3 for default-tolerance runs, else 0.
+    pub cov_floor: f64,
+}
+
+impl Tolerances {
+    /// A run at approximation tolerance `tol` (score.py: `acc_tol = 10 * tol`, no floor).
+    pub fn at(tol: f64) -> Tolerances {
+        Tolerances { acc_tol: Some(10.0 * tol), cov_floor: 0.0 }
+    }
+
+    /// A run at the engine's default tolerance.
+    pub fn default_tol() -> Tolerances {
+        Tolerances { acc_tol: None, cov_floor: 5e-3 }
+    }
+}
+
+/// The score of one case (score.py's result record).
+#[derive(Clone, Debug)]
+pub struct Score {
+    pub id: String,
+    pub class: String,
+    pub verdict: Verdict,
+    pub flags: Vec<String>,
+    /// The effective coverage tolerance.
+    pub cov_tol: f64,
+    /// The crash description (CRASH only).
+    pub detail: Option<String>,
+    pub status: Option<i64>,
+    pub status_name: Option<String>,
+    pub time_ms: Option<f64>,
+    pub n_curves: usize,
+    pub n_out_points: usize,
+    pub types: Vec<String>,
+    pub max_dist_3d: f64,
+    pub max_dist_uv: f64,
+    pub coverage: Vec<Coverage>,
+    pub hausdorff_truth_to_out: f64,
+    pub hausdorff_out_to_truth: f64,
+    pub expected_topology: Vec<String>,
+    pub output_topology: Vec<String>,
+}
+
+/// Python's `"%.<prec>e" % x`.
+pub fn fmt_e(x: f64, prec: usize) -> String {
+    if x.is_nan() {
+        return "nan".into();
+    }
+    if x.is_infinite() {
+        return if x > 0.0 { "inf".into() } else { "-inf".into() };
+    }
+    let s = format!("{x:.prec$e}");
+    let Some((m, e)) = s.split_once('e') else { return s };
+    let (sign, digits) = match e.strip_prefix('-') {
+        Some(d) => ('-', d),
+        None => ('+', e),
+    };
+    format!("{m}e{sign}{digits:0>2}")
+}
+
+/// Python's `"%.0f" % x`.
+fn fmt_f0(x: f64) -> String {
+    format!("{x:.0}")
+}
+
+/// Python's `str(list_of_str)`.
+fn py_list(v: &[String]) -> String {
+    format!("[{}]", v.iter().map(|s| format!("'{s}'")).collect::<Vec<_>>().join(", "))
+}
+
+fn finite(p: &[V3]) -> Vec<V3> {
+    p.iter().copied().filter(|q| q.iter().all(|c| c.is_finite())).collect()
+}
+
+/// Scores one engine outcome against a case (score.py's `score_case`).
+pub fn score_case(case: &Case, outcome: &Outcome, tol: Tolerances) -> Score {
+    let exp = &case.expected;
+    let cov_tol = case.cov_tol;
+    let mut cov_eff = cov_tol;
+    if let (Some(ill), Some(acc)) = (case.ill_cond, tol.acc_tol)
+        && ill != 0.0
+        && acc != 0.0
+    {
+        cov_eff = py_max(cov_tol, geom::py_min(0.5, 2.0 * acc * ill));
+    }
+    cov_eff = py_max(cov_eff, tol.cov_floor);
+    let join_tol = py_max(cov_tol, tol.cov_floor);
+    let mut r = Score {
+        id: case.id.clone(),
+        class: case.class.clone(),
+        verdict: Verdict::Crash,
+        flags: Vec::new(),
+        cov_tol: cov_eff,
+        detail: None,
+        status: None,
+        status_name: None,
+        time_ms: None,
+        n_curves: 0,
+        n_out_points: 0,
+        types: Vec::new(),
+        max_dist_3d: 0.0,
+        max_dist_uv: 0.0,
+        coverage: Vec::new(),
+        hausdorff_truth_to_out: 0.0,
+        hausdorff_out_to_truth: 0.0,
+        expected_topology: Vec::new(),
+        output_topology: Vec::new(),
+    };
+    let out = match outcome {
+        Outcome::Crashed(why) => {
+            r.detail = Some(why.clone());
+            return r;
+        }
+        Outcome::Ran(out) => out,
+    };
+    let status_name = out.status_name.clone().unwrap_or_else(|| out.status.to_string());
+    r.status = Some(out.status);
+    r.status_name = Some(status_name.clone());
+    r.time_ms = Some(out.time_ms);
+
+    let mut out_curves: Vec<Vec<V3>> = Vec::new();
+    let mut out_points: Vec<V3> = Vec::new();
+    let mut nonfinite = 0;
+    for c in &out.curves {
+        let p = finite(&c.xyz);
+        nonfinite += c.xyz.len() - p.len();
+        if p.is_empty() {
+            continue;
+        }
+        if c.curve_type == "TOUCHING" || diameter(&p) <= cov_tol {
+            let n = p.len() as f64;
+            let s = p.iter().fold([0.0; 3], |a, q| [a[0] + q[0], a[1] + q[1], a[2] + q[2]]);
+            out_points.push([s[0] / n, s[1] / n, s[2] / n]);
+        } else {
+            out_curves.push(p);
+        }
+        r.types.push(c.curve_type.clone());
+    }
+    out_points.extend(out.points.iter().copied());
+    r.n_curves = out.curves.len();
+    r.n_out_points = out_points.len();
+
+    // accuracy against the exact surfaces
+    let worst = |pts: &[V3]| -> f64 {
+        let mut m = np_max(&case.truth[0].distance(pts));
+        for t in &case.truth[1..] {
+            m = py_max(m, np_max(&t.distance(pts)));
+        }
+        m
+    };
+    let (mut acc3d, mut accuv) = (0.0_f64, 0.0_f64);
+    for c in &out.curves {
+        let p = finite(&c.xyz);
+        if !p.is_empty() {
+            acc3d = py_max(acc3d, worst(&p));
+        }
+        // S1(uv1) must lie on surface 2, S2(uv2) on surface 1.
+        for (k, q) in [(1, &c.s1_at_uv), (0, &c.s2_at_uv)] {
+            let q = finite(q);
+            if !q.is_empty() {
+                accuv = py_max(accuv, np_max(&case.truth[k].distance(&q)));
+            }
+        }
+    }
+    if !out_points.is_empty() {
+        acc3d = py_max(acc3d, worst(&out_points));
+    }
+    r.max_dist_3d = acc3d;
+    r.max_dist_uv = accuv;
+
+    let truth_polys: Vec<&[V3]> = exp.branches.iter().map(Vec::as_slice).chain(exp.coincident_regions.iter().map(Vec::as_slice)).collect();
+    let truth_pts = &exp.isolated_points;
+    let expected_any = !truth_polys.is_empty() || !truth_pts.is_empty();
+
+    let ok_status = out.status == SM_SUCCESS;
+    if !ok_status {
+        if !expected_any && out.curves.is_empty() && out_points.is_empty() {
+            r.flags.push(format!("status {status_name} with empty output (expected empty)"));
+        } else {
+            r.flags.push(format!("ERROR status {status_name}"));
+        }
+    }
+
+    // coverage truth -> output
+    let curve_refs: Vec<&[V3]> = out_curves.iter().map(Vec::as_slice).collect();
+    let outset = PolySet::new(&curve_refs, &out_points);
+    let mut missed = Vec::new();
+    for (k, p) in truth_polys.iter().enumerate() {
+        let d = outset.query(&stride(p, p.len() / 1000));
+        let frac = d.iter().filter(|&&x| x <= cov_eff).count() as f64 / d.len() as f64;
+        let dmax = np_max(&d);
+        r.coverage.push(Coverage::Polyline { max: dmax, covered_frac: frac });
+        if dmax > cov_eff {
+            let kind = if k >= exp.branches.len() { "region boundary".to_string() } else { format!("branch {k}") };
+            missed.push(format!("{kind} covered {}% (max gap {})", fmt_f0(100.0 * frac), fmt_e(dmax, 2)));
+        }
+    }
+    for (k, &p) in truth_pts.iter().enumerate() {
+        let d = outset.query1(p);
+        r.coverage.push(Coverage::Point { point: k, dist: d });
+        if d > cov_eff {
+            missed.push(format!("isolated point {k} not found (nearest output {})", fmt_e(d, 2)));
+        }
+    }
+    let mut h = 0.0;
+    for (i, c) in r.coverage.iter().enumerate() {
+        let x = match c {
+            Coverage::Polyline { max, .. } => *max,
+            Coverage::Point { dist, .. } => *dist,
+        };
+        h = if i == 0 { x } else { py_max(h, x) };
+    }
+    r.hausdorff_truth_to_out = h;
+
+    // output -> truth (spurious output)
+    let tset = PolySet::new(&truth_polys, truth_pts);
+    let mut spurious = Vec::new();
+    let mut h_out = 0.0;
+    for (k, p) in out_curves.iter().enumerate() {
+        let d = tset.query(&stride(p, p.len() / 1000));
+        let dmax = np_max(&d);
+        h_out = py_max(h_out, dmax);
+        if dmax > cov_eff {
+            let off = d.iter().filter(|&&x| x > cov_eff).count() as f64 / d.len() as f64;
+            spurious.push(format!("output curve {k} off the true intersection by up to {} ({}% of it)", fmt_e(dmax, 2), fmt_f0(100.0 * off)));
+        }
+    }
+    for (k, &p) in out_points.iter().enumerate() {
+        let d = tset.query1(p);
+        h_out = py_max(h_out, d);
+        if d > cov_eff {
+            spurious.push(format!("output point {k} off the true intersection by {}", fmt_e(d, 2)));
+        }
+    }
+    r.hausdorff_out_to_truth = h_out;
+
+    // topology signature
+    let truth_owned: Vec<Vec<V3>> = truth_polys.iter().map(|p| p.to_vec()).collect();
+    let sig_exp = signature(&truth_owned, truth_pts, join_tol);
+    let sig_out = signature(&out_curves, &out_points, join_tol);
+    let mut topo = Vec::new();
+    if nonfinite > 0 {
+        topo.push(format!("{nonfinite} non-finite output sample(s)"));
+    }
+    // duplicated / overlapping output curves
+    let osets: Vec<PolySet> = out_curves.iter().map(|c| PolySet::new(&[c], &[])).collect();
+    for i in 0..out_curves.len() {
+        for j in i + 1..out_curves.len() {
+            let (ci, cj) = (&out_curves[i], &out_curves[j]);
+            let di = osets[j].query(&stride(ci, ci.len() / 500));
+            let dj = osets[i].query(&stride(cj, cj.len() / 500));
+            let fi = di.iter().filter(|&&x| x <= join_tol).count() as f64 / di.len() as f64;
+            let fj = dj.iter().filter(|&&x| x <= join_tol).count() as f64 / dj.len() as f64;
+            let frac = py_max(fi, fj);
+            if frac > 0.2 {
+                topo.push(format!("output curves {i} and {j} overlap ({}% of one lies on the other)", fmt_f0(100.0 * frac)));
+            }
+        }
+    }
+    if sig_exp != sig_out {
+        topo.push(format!("topology {}, expected {}", py_list(&sig_out), py_list(&sig_exp)));
+    }
+    for (k, &s) in exp.singular_points.iter().enumerate() {
+        let ends = out_curves.iter().filter_map(|c| c.first()).chain(out_curves.iter().filter_map(|c| c.last())).chain(out_points.iter());
+        let mut d = f64::INFINITY;
+        for (n, &e) in ends.enumerate() {
+            let x = geom::norm(geom::sub(e, s));
+            d = if n == 0 { x } else { geom::py_min(d, x) };
+        }
+        if d > cov_eff {
+            topo.push(format!("singular point {k} not identified (no curve end/point within {}; nearest {})", fmt_e(cov_eff, 1), fmt_e(d, 2)));
+        }
+    }
+    let mut inacc = Vec::new();
+    if let Some(acc) = tol.acc_tol {
+        if acc3d > acc {
+            inacc.push(format!("3D output up to {} from a surface (acc_tol {})", fmt_e(acc3d, 2), fmt_e(acc, 0)));
+        }
+        if accuv > acc {
+            inacc.push(format!("surface at output UV up to {} from the other surface", fmt_e(accuv, 2)));
+        }
+    }
+
+    let has_output = !out.curves.is_empty() || !out_points.is_empty();
+    let verdict = if !ok_status && (expected_any || has_output) {
+        Verdict::Error
+    } else if !missed.is_empty() {
+        Verdict::Incomplete
+    } else if !spurious.is_empty() || !topo.is_empty() {
+        Verdict::WrongTopology
+    } else if !inacc.is_empty() {
+        Verdict::Inaccurate
+    } else {
+        Verdict::Pass
+    };
+    r.flags.extend(missed);
+    r.flags.extend(spurious);
+    r.flags.extend(topo);
+    r.flags.extend(inacc);
+    r.expected_topology = sig_exp;
+    r.output_topology = sig_out;
+    r.verdict = verdict;
+    if tol.cov_floor != 0.0
+        && let Some(fs) = case.feature_size
+        && fs < tol.cov_floor
+    {
+        r.flags.insert(0, format!("feature size {} below the {} judging floor of default-tolerance runs", fmt_e(fs, 1), fmt_e(tol.cov_floor, 0)));
+        r.verdict = Verdict::NotApplicable;
+    }
+    r
+}
+
+/// Parses `#define SM_NAME 1234L` lines of SMLib's `SmMessages.h` (first name per number), for
+/// status names in flags as score.py prints them.
+pub fn status_names(sm_messages_h: &str) -> std::collections::HashMap<i64, String> {
+    let mut names = std::collections::HashMap::new();
+    for line in sm_messages_h.lines() {
+        let Some(rest) = line.strip_prefix("#define") else { continue };
+        if !rest.starts_with(char::is_whitespace) {
+            continue;
+        }
+        let mut it = rest.split_whitespace();
+        let (Some(name), Some(val)) = (it.next(), it.next()) else { continue };
+        let Some(digits) = val.strip_suffix('L') else { continue };
+        if !name.starts_with("SM_") || !name[3..].chars().all(|c| c.is_ascii_alphanumeric() || c == '_') || name.len() == 3 {
+            continue;
+        }
+        if digits.is_empty() || !digits.chars().all(|c| c.is_ascii_digit()) {
+            continue;
+        }
+        if let Ok(n) = digits.parse::<i64>() {
+            names.entry(n).or_insert_with(|| name.to_string());
+        }
+    }
+    names
+}
+
+/// `v3` re-exported for callers building points from JSON.
+pub fn point(v: &Value) -> Result<V3, String> {
+    v3(v)
+}