From a16de16e1ee447af9939e1416c03f338eeceb7e3 Mon Sep 17 00:00:00 2001 From: Vince Buffalo Date: Sat, 11 Jul 2026 14:25:09 -0700 Subject: [PATCH] feat(survey): evaluate LHS points in greedy-nearest order MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Replace the flat index-order `par_iter` sweep with a batched greedy schedule: 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. The payoff is a live signal. The progress bar's best-loglik metric now climbs fast and plateaus early when the box holds a basin, which is an at-a-glance read on whether the survey bounds are placed well, and the good region is the part that fills in first — easier to watch, and quick to notice a mis-placed box. This reorders WHEN points run, never WHICH: every point is still evaluated (no early stop), each point's result is keyed on (seed, point_id) via `derive_point_seed` / `mix_cell_seed` rather than a shared RNG stream, and the final sort by (loglik desc, point_id asc) is order-independent. So `landscape.tsv` is byte-identical to the old sweep — no golden, expected-output, or CAS-identity impact. The distance metric reuses `EstimatedParam::to_transformed`, so points on log/logit params are compared on the same scale the sampler and landscape geometry live on; it only orders evaluation and never touches a loglik. `survey_greedy_order_invariant.rs` pins the invariant from outside: the same survey at `--parallel 1` vs `--parallel 4` walks two different greedy orders and must produce byte-identical `landscape.tsv`. If a future change makes a point's result order-dependent, it goes red. --- rust/crates/cli/src/survey.rs | 163 +++++++++++++----- .../tests/survey_greedy_order_invariant.rs | 155 +++++++++++++++++ 2 files changed, 279 insertions(+), 39 deletions(-) create mode 100644 rust/crates/cli/tests/survey_greedy_order_invariant.rs diff --git a/rust/crates/cli/src/survey.rs b/rust/crates/cli/src/survey.rs index 07586c37..85bfce99 100644 --- a/rust/crates/cli/src/survey.rs +++ b/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(), - ¶ms, &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, - ¶ms, &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(), + ¶ms, &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, + ¶ms, &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> = 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> = (0..n).map(|_| None).collect(); + let mut done = vec![false; n]; + let mut best_point: Option = 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 = if n_done == 0 { + (0..seed_size).collect() + } else { + let mut remaining: Vec = + (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::>() + }; let rows: Vec = 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, diff --git a/rust/crates/cli/tests/survey_greedy_order_invariant.rs b/rust/crates/cli/tests/survey_greedy_order_invariant.rs new file mode 100644 index 00000000..6d1d09df --- /dev/null +++ b/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 { + 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", ¶llel.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()); +}