feat(survey): evaluate LHS points in greedy-nearest order

rust/crates/cli/src/survey.rs

@@ -457,50 +457,113 @@ pub fn cmd_survey(a: &crate::args::SurveyArgs) {
             }))
     };
 
-    // Progress: one overall bar over the LHS points, ticked from the parallel
-    // sweep (`Task` is `Send + Sync`), with the best loglik found so far as the
-    // researcher metric. Honors `--progress none/plain`.
+    // Progress: one overall bar over the LHS points, with the best loglik found
+    // so far as the researcher metric. Honors `--progress none/plain`.
     let bar = crate::progress::Reporter::new().task(n_points as u64, "survey", "pts");
-    let best = std::sync::Mutex::new(f64::NEG_INFINITY);
-    let sweep = || lhs_starts.par_iter().enumerate()
-        .map(|(point_id, draw)| {
-            // Build the full parameter vector: base_params overwritten
-            // at each estimated index. Fixed params are already baked
-            // into base_params (resolve_survey_inputs).
-            let mut params = resolved.base_params.clone();
-            for spec in draw {
-                params[spec.index] = spec.initial;
-            }
-            let row = match eval_method {
-                SurveyEvalMethod::Pfilter => eval_point_pfilter(
-                    &process, obs_model.as_ref(),
-                    &params, &resolved.estimated, draw,
-                    a.eval_particles, a.eval_replicates,
-                    smc_dt, t_start, a.seed, point_id,
-                ),
-                SurveyEvalMethod::Auto => unreachable!(
-                    "Auto resolved before parallel eval loop"),
-                SurveyEvalMethod::Simulate => eval_point_simulate(
-                    &resolved.compiled, &obs_model, &obs_times,
-                    &params, &resolved.estimated, draw,
-                    smc_dt, point_id,
-                ),
+
+    // Evaluate a single LHS point. A point's result is a pure function of
+    // (seed, point_id) — `derive_point_seed` / `mix_cell_seed` key the RNG on
+    // the point index, not on evaluation order — so the greedy schedule below
+    // reorders WHEN points run without changing any row. `bar.inc` is the only
+    // shared-state touch; the running best is folded in single-threaded after
+    // each batch.
+    let eval_one = |point_id: usize| -> LandscapeRow {
+        let draw = &lhs_starts[point_id];
+        // Build the full parameter vector: base_params overwritten at each
+        // estimated index. Fixed params are already baked into base_params.
+        let mut params = resolved.base_params.clone();
+        for spec in draw {
+            params[spec.index] = spec.initial;
+        }
+        let row = match eval_method {
+            SurveyEvalMethod::Pfilter => eval_point_pfilter(
+                &process, obs_model.as_ref(),
+                &params, &resolved.estimated, draw,
+                a.eval_particles, a.eval_replicates,
+                smc_dt, t_start, a.seed, point_id,
+            ),
+            SurveyEvalMethod::Auto => unreachable!(
+                "Auto resolved before parallel eval loop"),
+            SurveyEvalMethod::Simulate => eval_point_simulate(
+                &resolved.compiled, &obs_model, &obs_times,
+                &params, &resolved.estimated, draw,
+                smc_dt, point_id,
+            ),
+        };
+        bar.inc(1);
+        row
+    };
+
+    // Greedy-nearest evaluation order (progressive, never early-stopping).
+    // Seed a coarse spread of points, then each round evaluate the unevaluated
+    // points nearest — in transform-normalized parameter space — to the current
+    // best-loglik point, a batch at a time so every core stays busy. Every
+    // point is still evaluated: this reorders WHEN, never WHICH, so
+    // `landscape.tsv` is byte-identical to an index-order sweep (the final sort
+    // by (loglik desc, point_id asc) is order-independent, and each point's RNG
+    // is keyed on point_id). The payoff is a live signal: the best-loglik metric
+    // climbs fast and plateaus early when the box holds a basin — an
+    // at-a-glance read on whether the survey bounds are placed well — and the
+    // good region is the part that fills in first.
+    let normed: Vec<Vec<f64>> = lhs_starts.iter()
+        .map(|draw| draw.iter().map(|ep| survey_norm_coord(ep, ep.initial)).collect())
+        .collect();
+    let greedy = || {
+        let n = lhs_starts.len();
+        let n_threads = rayon::current_num_threads().max(1);
+        let batch_size = n_threads;
+        // Coarse global spread before drilling: √n, at least one full core-width.
+        let seed_size = n_threads.max((n as f64).sqrt().ceil() as usize).min(n);
+
+        let mut result: Vec<Option<LandscapeRow>> = (0..n).map(|_| None).collect();
+        let mut done = vec![false; n];
+        let mut best_point: Option<usize> = None;
+        let mut running_best = f64::NEG_INFINITY;
+        let mut n_done = 0usize;
+
+        while n_done < n {
+            // This round's batch: the seed spread first, then the unevaluated
+            // points closest to the running-best point. Ties (and the
+            // no-finite-best-yet case) fall back to ascending point_id, so the
+            // schedule is deterministic given the point set.
+            let batch: Vec<usize> = if n_done == 0 {
+                (0..seed_size).collect()
+            } else {
+                let mut remaining: Vec<usize> =
+                    (0..n).filter(|&i| !done[i]).collect();
+                if let Some(bp) = best_point {
+                    remaining.sort_by(|&x, &y| {
+                        let dx = survey_dist2(&normed[bp], &normed[x]);
+                        let dy = survey_dist2(&normed[bp], &normed[y]);
+                        dx.partial_cmp(&dy)
+                            .unwrap_or(std::cmp::Ordering::Equal)
+                            .then(x.cmp(&y))
+                    });
+                }
+                remaining.truncate(batch_size);
+                remaining
             };
-            bar.inc(1);
-            if row.loglik.is_finite() {
-                if let Ok(mut b) = best.lock() {
-                    if row.loglik > *b {
-                        *b = row.loglik;
-                        bar.set(crate::progress::best_ll(*b));
-                    }
+            let batch_rows: Vec<(usize, LandscapeRow)> = batch.par_iter()
+                .map(|&pid| (pid, eval_one(pid)))
+                .collect();
+            for (pid, row) in batch_rows {
+                if row.loglik.is_finite() && row.loglik > running_best {
+                    running_best = row.loglik;
+                    best_point = Some(pid);
+                    bar.set(crate::progress::best_ll(running_best));
                 }
+                result[pid] = Some(row);
+                done[pid] = true;
+                n_done += 1;
             }
-            row
-        })
-        .collect();
+        }
+        result.into_iter()
+            .map(|o| o.expect("every survey point evaluated"))
+            .collect::<Vec<LandscapeRow>>()
+    };
     let rows: Vec<LandscapeRow> = match &survey_pool {
-        Some(pool) => pool.install(sweep),
-        None => sweep(),
+        Some(pool) => pool.install(greedy),
+        None => greedy(),
     };
     bar.finish();
 
@@ -998,6 +1061,28 @@ struct LandscapeRow {
     n_replicates: usize,
 }
 
+/// Normalize a natural-scale parameter value to `[0, 1]` on the parameter's
+/// *transform* scale — the same scale the sampler and the landscape geometry
+/// live on (log for rates, logit for probabilities, linear otherwise). Used
+/// only to order survey evaluation (nearest-to-best first); it never touches a
+/// `loglik`, so the exact metric is a heuristic, not a result.
+fn survey_norm_coord(ep: &EstimatedParam, x: f64) -> f64 {
+    let zlo = ep.to_transformed(ep.lower);
+    let zhi = ep.to_transformed(ep.upper);
+    let denom = zhi - zlo;
+    if !denom.is_finite() || denom.abs() < 1e-12 {
+        return 0.0; // degenerate / near-fixed dimension contributes no distance
+    }
+    ((ep.to_transformed(x) - zlo) / denom).clamp(0.0, 1.0)
+}
+
+/// Squared Euclidean distance between two transform-normalized coordinate
+/// vectors (see [`survey_norm_coord`]). Squared because ordering only needs the
+/// monotone comparison — the sqrt would be wasted.
+fn survey_dist2(a: &[f64], b: &[f64]) -> f64 {
+    a.iter().zip(b).map(|(x, y)| (x - y) * (x - y)).sum()
+}
+
 #[allow(clippy::too_many_arguments)]
 fn eval_point_pfilter(
     process: &ChainBinomialProcess,

rust/crates/cli/tests/survey_greedy_order_invariant.rs

@@ -0,0 +1,155 @@
+//! `camdl survey` evaluates LHS points in a greedy-nearest order (drill toward
+//! the running-best point) so the best-loglik metric converges fast and the
+//! interesting region fills in first. That reordering must be *purely* a
+//! scheduling change: every point is still evaluated, and each point's result
+//! is keyed on `(seed, point_id)`, so `landscape.tsv` must be **byte-identical**
+//! regardless of the order points happen to run in.
+//!
+//! This test pins that invariant the only way it can be observed from outside:
+//! run the same survey twice with different `--parallel` throttles. A different
+//! thread budget means a different batch width, hence a different greedy
+//! evaluation order — and yet the two landscape files must match exactly. If a
+//! future change makes a point's result depend on when it ran (a shared RNG
+//! stream, an order-dependent accumulator), this test goes red.
+//!
+//! `--eval simulate` keeps it deterministic and sub-second (no particle filter).
+
+use std::path::{Path, PathBuf};
+use std::process::Command;
+
+fn camdl_bin() -> PathBuf {
+    let manifest = std::env::var("CARGO_MANIFEST_DIR")
+        .expect("CARGO_MANIFEST_DIR set under cargo test");
+    let p = Path::new(&manifest).join("../../target/release/camdl");
+    assert!(
+        p.exists(),
+        "release camdl binary missing: {} - run `make build-rust` or `make test` (gh#105)",
+        p.display()
+    );
+    p
+}
+
+fn camdlc_bin() -> PathBuf {
+    let manifest = std::env::var("CARGO_MANIFEST_DIR")
+        .expect("CARGO_MANIFEST_DIR set under cargo test");
+    let p = Path::new(&manifest).join("../../../ocaml/_build/default/bin/camdlc.exe");
+    assert!(p.exists(), "camdlc.exe missing: {} - run `make build-ocaml`", p.display());
+    p
+}
+
+struct Tmp(PathBuf);
+impl Tmp { fn path(&self) -> &Path { &self.0 } }
+impl Drop for Tmp { fn drop(&mut self) { let _ = std::fs::remove_dir_all(&self.0); } }
+fn tempdir(tag: &str) -> Tmp {
+    let ns = std::time::SystemTime::now()
+        .duration_since(std::time::UNIX_EPOCH).unwrap().as_nanos();
+    let base = std::env::temp_dir().join(format!(
+        "camdl_survey_order_{}_{}_{}", tag, std::process::id(), ns));
+    std::fs::create_dir_all(&base).unwrap();
+    Tmp(base)
+}
+
+/// A small deterministic SIR + dataset. Three estimated params (one fixed) so
+/// the transform-normalized distance metric that drives the greedy order has a
+/// nontrivial space to walk.
+fn write_fixture(dir: &Path) -> (PathBuf, PathBuf) {
+    let camdlc = camdlc_bin();
+    let src = r#"
+time_unit = 'days
+compartments { S, I, R }
+parameters {
+  beta  : rate  in [0.001, 5.0]
+  gamma : rate  in [0.01, 1.0]
+  N0    : count in [100, 10000]
+}
+transitions {
+  infection : S --> I @ beta * S * I / N0
+  recovery  : I --> R @ gamma * I
+}
+observations {
+  cases {
+    columns       { time : time, cases : count }
+    projected  = prevalence(I)
+    emit_schedule = every 1 'days
+    cases ~ poisson(rate = projected)
+  }
+}
+init { S = 999  I = 1 }
+simulate { from = 0 'days  to = 6 'days }
+"#;
+    let model_path = dir.join("sir.camdl");
+    std::fs::write(&model_path, src).unwrap();
+    let ir_path = dir.join("sir.ir.json");
+    let out = Command::new(&camdlc).arg(&model_path).output().unwrap();
+    assert!(out.status.success(),
+        "camdlc failed: {}", String::from_utf8_lossy(&out.stderr));
+    std::fs::write(&ir_path, &out.stdout).unwrap();
+
+    let data_path = dir.join("cases.tsv");
+    std::fs::write(&data_path,
+        "time\tcases\n1\t2\n2\t4\n3\t8\n4\t6\n5\t4\n6\t2\n").unwrap();
+
+    (ir_path, data_path)
+}
+
+/// The one `surveys/…` leaf holding `landscape.tsv`.
+fn find_landscape(root: &Path) -> PathBuf {
+    let mut stack = vec![root.join("surveys")];
+    while let Some(dir) = stack.pop() {
+        if dir.join("landscape.tsv").is_file() {
+            return dir.join("landscape.tsv");
+        }
+        if let Ok(entries) = std::fs::read_dir(&dir) {
+            for e in entries.flatten() {
+                let p = e.path();
+                if p.is_dir() { stack.push(p); }
+            }
+        }
+    }
+    panic!("no landscape.tsv under {}", root.display());
+}
+
+/// Run a survey with a given `--parallel` throttle; return the landscape bytes.
+fn run_survey(bin: &Path, ir: &Path, data: &Path, out_root: &Path, parallel: u32) -> Vec<u8> {
+    let out = Command::new(bin)
+        .env("CAMDL_SKIP_VERSION_CHECK", "1")
+        .args([
+            "survey", &ir.to_string_lossy(),
+            "--data", &data.to_string_lossy(),
+            "--estimate", "beta=0.001:5.0",
+            "--estimate", "gamma=0.01:1.0",
+            "--fixed", "N0=1000",
+            "--eval", "simulate",
+            "--n-points", "250",
+            "--seed", "1",
+            "--parallel", &parallel.to_string(),
+            "--output", &out_root.to_string_lossy(),
+        ])
+        .output().expect("spawn camdl survey");
+    assert!(out.status.success(),
+        "camdl survey (--parallel {}) must exit 0.\nstderr:\n{}",
+        parallel, String::from_utf8_lossy(&out.stderr));
+    std::fs::read(find_landscape(out_root)).expect("read landscape.tsv")
+}
+
+#[test]
+fn survey_landscape_is_byte_identical_across_parallelism() {
+    let bin = camdl_bin();
+    let tmp = tempdir("inv");
+    let (ir, data) = write_fixture(tmp.path());
+
+    // Single-threaded (batch width 1) vs multi-threaded (wider batches) exercise
+    // two very different greedy evaluation orders over the same 250 points.
+    let serial   = run_survey(&bin, &ir, &data, &tmp.path().join("p1"), 1);
+    let parallel = run_survey(&bin, &ir, &data, &tmp.path().join("p4"), 4);
+
+    // Non-vacuous: header + 250 data rows.
+    let n_lines = serial.iter().filter(|&&b| b == b'\n').count();
+    assert!(n_lines >= 250,
+        "expected ~250 landscape rows, got {} lines — fixture/eval may be broken", n_lines);
+
+    assert!(serial == parallel,
+        "survey landscape.tsv must be byte-identical regardless of evaluation \
+         order (greedy scheduling is order-only); serial vs parallel differ.\n\
+         serial {} bytes, parallel {} bytes", serial.len(), parallel.len());
+}