From bfc213459f7e5cfe7bddb819162365cb59d5d819 Mon Sep 17 00:00:00 2001 From: Alex Nodeland Date: Tue, 28 Jul 2026 11:10:06 -0400 Subject: [PATCH] Chebyshev scalarization + evolution-as-inference explorable (0.3.1) Chebyshev (inference/pareto.rs): ChebyshevScalarization - weighted-max norm over an ideal point, latent or fixed weight. Reaches non-convex Pareto front interiors where weighted-sum optima provably collapse to endpoints; pinned at the theorem level on the concave front f1=x, f2=1-x^2 (fixed-w contrast + weight sweep + latent-w conditional tracking). Docs now state the latent-w marginal tilt exp(-s*m(w)) honestly for both scalarizations. Explorable (crates/fugue-evo-wasm + docs): ExploreSmcInference steps the real inference layer one tempering rung per call (GaussianPrior program, twin-peaks factor likelihood, fugue smc_prior_particles/normalize/resample/ rejuvenate_particles/CrossoverKernel), streaming particles, ESS, resampling events, crossover swaps, and the running log-evidence; density_grid returns the exact tempered target for the heat overlay. Engine tests: determinism, ladder shape, analytic posterior mean + evidence vs grid quadrature, annealing concentration (distance-to-nearest-mode). Widget (viz/inference.js) follows the site conventions (wasm gate + notice, lazy init, seeded scrub, controls/canvas/instruction/readouts, reduced-motion batch); mounted on the "Evolution as Inference" architecture page; wasm crate gains ppl feature + direct fugue-ppl dep. Verified live in a browser against the built pkg (autoplay ladder, Replay, Reset, Step, readouts). Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01HmWLBR9Zq6hPG5o7URDuNJ --- CHANGELOG.md | 33 ++ Cargo.lock | 3 +- Cargo.toml | 2 +- crates/fugue-evo-wasm/Cargo.toml | 3 +- .../fugue-evo-wasm/src/explore_inference.rs | 442 ++++++++++++++++++ crates/fugue-evo-wasm/src/lib.rs | 2 + docs/book.toml | 1 + docs/src/architecture/fugue-integration.md | 14 + docs/viz/inference.js | 432 +++++++++++++++++ src/inference/pareto.rs | 292 +++++++++++- 10 files changed, 1211 insertions(+), 13 deletions(-) create mode 100644 crates/fugue-evo-wasm/src/explore_inference.rs create mode 100644 docs/viz/inference.js diff --git a/CHANGELOG.md b/CHANGELOG.md index 1f3d1a4..46f783d 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,39 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] +## [0.3.1] - 2026-07-28 + +### Added + +- **Chebyshev scalarization for non-convex Pareto fronts** + (`inference::pareto::ChebyshevScalarization`): the weighted-max norm + `max_i w_i·(f_i − z_i)` over an ideal point `z`, reaching every weakly + Pareto-optimal point — including the non-convex front regions where every + weighted-sum optimum provably collapses onto the endpoints. Supports a + latent weight (like `ParetoScalarization`) or a **fixed** weight + (`with_weight`) for uniform front sweeps. Pinned at the theorem level on + the concave front `f1 = x, f2 = 1 − x²`: fixed-w weighted-sum mass avoids + the interior (< 0.1) while fixed-w Chebyshev concentrates on the interior + front point `x* = (√5−1)/2`, and a weight sweep traces the whole front. +- **Honest marginal-tilt documentation** for the latent-weight Pareto models: + the `w`-marginal is tilted by `exp(−s·m(w))` (the scalarized optimum's + value), so the *conditional* `x | w` is what tracks the front; uniform + coverage comes from fixed-weight sweeps. New regression + `test_chebyshev_latent_weight_conditional_tracks_front` pins the + conditional property. +- **"Evolution as inference" explorable** (evo.fugue.run, Architecture → + Evolution as Inference): `ExploreSmcInference` in `fugue-evo-wasm` steps + the crate's real inference layer one tempering rung at a time — Gaussian + prior program, twin-peaks Boltzmann target, fugue's SMC primitives + (reweight / ESS-triggered systematic resampling / typed-MH rejuvenation / + crossover kernel), exact tempered-density heat recomputed each rung, + β-ladder up to annealed-optimizer territory, and live ESS / log-evidence / + swap readouts. The wasm crate now enables the `ppl` feature and depends on + `fugue-ppl` directly. Engine pinned by determinism, ladder-shape, + analytic-posterior-mean, analytic-evidence, and annealing-concentration + tests; the page verified live in a browser against the built wasm. + + ## [0.3.0] - 2026-07-28 **Inference-first.** fugue-evo's identity is now "an implementation of fugue diff --git a/Cargo.lock b/Cargo.lock index 24de9d9..70483c9 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -278,7 +278,7 @@ checksum = "3f9eec918d3f24069decb9af1554cad7c880e2da24a9afd88aca000531ab82c1" [[package]] name = "fugue-evo" -version = "0.3.0" +version = "0.3.1" dependencies = [ "approx", "bincode", @@ -304,6 +304,7 @@ version = "0.1.0" dependencies = [ "console_error_panic_hook", "fugue-evo", + "fugue-ppl", "getrandom 0.2.16", "js-sys", "rand", diff --git a/Cargo.toml b/Cargo.toml index b3eba10..b53d94c 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -12,7 +12,7 @@ lto = true [package] name = "fugue-evo" -version = "0.3.0" +version = "0.3.1" edition = "2021" authors = ["Alex Nodeland"] description = "An implementation of fugue for running evolutionary algorithms as Bayesian inference: priors and likelihoods as probabilistic programs, tempered SMC in trace space, annealed optimization, Pareto posteriors - plus a standalone classical EC toolkit" diff --git a/crates/fugue-evo-wasm/Cargo.toml b/crates/fugue-evo-wasm/Cargo.toml index 757cb48..4ffb0f7 100644 --- a/crates/fugue-evo-wasm/Cargo.toml +++ b/crates/fugue-evo-wasm/Cargo.toml @@ -18,7 +18,8 @@ crate-type = ["cdylib", "rlib"] [dependencies] # Classic EC layer without parallel/checkpoint/ppl features -fugue-evo = { path = "../..", default-features = false, features = ["std", "classic"] } +fugue-evo = { path = "../..", default-features = false, features = ["std", "classic", "ppl"] } +fugue-ppl = { path = "../../../fugue", version = "0.2.1" } # WASM bindings wasm-bindgen = "0.2" diff --git a/crates/fugue-evo-wasm/src/explore_inference.rs b/crates/fugue-evo-wasm/src/explore_inference.rs new file mode 100644 index 0000000..6681d8b --- /dev/null +++ b/crates/fugue-evo-wasm/src/explore_inference.rs @@ -0,0 +1,442 @@ +//! Incremental tempered-SMC explorable: evolution as inference, live. +//! +//! [`ExploreSmcInference`] drives the crate's *real* inference layer — a +//! [`GaussianPrior`] program, a bimodal fitness entering as a +//! [`FactorFitness`] likelihood, and fugue's SMC primitives +//! (`smc_prior_particles` / `normalize_particles` / `resample_particles` / +//! `rejuvenate_particles` / `CrossoverKernel`) — one tempering rung per +//! `step()`, streaming the particle population, ESS, resampling events, +//! accepted crossover swaps, and the running log-evidence as JSON. +//! +//! Module conventions match `explore.rs`: explicit `u64` seeds (a seed is a +//! replayable recording), every JS-supplied parameter clamped, and 2-D +//! genomes so positions map straight to the canvas. `density_grid` returns +//! the analytic **negative** log-target (lower = better, like the landscape +//! grids) so the same heat conventions apply and the particles can be seen +//! matching the exact tempered density at every β. + +use rand::rngs::StdRng; +use rand::SeedableRng; +use serde_json::json; +use wasm_bindgen::prelude::*; + +use fugue::{ + addr, effective_sample_size, normalize_particles, rejuvenate_particles, resample_particles, + smc_prior_particles, CrossoverKernel, Particle, PopulationKernel, ResamplingMethod, Trace, +}; +use fugue_evo::fitness::traits::Fitness; +use fugue_evo::genome::real_vector::RealVector; +use fugue_evo::genome::traits::RealValuedGenome; +use fugue_evo::inference::likelihood::FactorFitness; +use fugue_evo::inference::model::EvolutionModel; +use fugue_evo::inference::prior::GaussianPrior; + +/// Plot domain (matches the prior's ±2.25σ window). +const DOMAIN: (f64, f64) = (-4.5, 4.5); +/// Prior standard deviation per coordinate. +const PRIOR_STD: f64 = 2.0; + +/// The two modes of the twin-peaks fitness: (center, std, mixture weight). +const MODES: [([f64; 2], f64, f64); 2] = [([-1.6, -0.9], 0.55, 0.65), ([1.7, 1.2], 0.8, 0.35)]; + +/// Log-mixture-of-Gaussians fitness (higher is better): two unequal peaks, +/// so the Boltzmann posterior is visibly bimodal at β = 1 and annealing +/// (β > 1) shifts mass toward the sharper peak — "selection pressure is +/// conditioning; temperature is selection strength", drawable. +#[derive(Clone, Copy)] +struct TwinPeaks; + +fn twin_peaks_f(x: f64, y: f64) -> f64 { + let mut s = 0.0; + for (m, sd, w) in MODES { + let r2 = (x - m[0]).powi(2) + (y - m[1]).powi(2); + s += w * (-r2 / (2.0 * sd * sd)).exp(); + } + if s > 0.0 { + s.ln() + } else { + -1e12 // far outside the domain: crushed, never accepted + } +} + +impl Fitness for TwinPeaks { + type Genome = RealVector; + type Value = f64; + fn evaluate(&self, g: &RealVector) -> f64 { + let genes = g.genes(); + twin_peaks_f(genes[0], genes[1]) + } +} + +/// Prior log-density (up to a constant): i.i.d. `N(0, PRIOR_STD²)`. +fn log_prior(x: f64, y: f64) -> f64 { + -(x * x + y * y) / (2.0 * PRIOR_STD * PRIOR_STD) +} + +fn loglik(t: &Trace) -> f64 { + t.log_likelihood + t.log_factors +} + +fn log_sum_exp(v: &[f64]) -> f64 { + let m = v.iter().cloned().fold(f64::NEG_INFINITY, f64::max); + if !m.is_finite() { + return m; + } + m + v.iter().map(|x| (x - m).exp()).sum::().ln() +} + +/// Incremental likelihood-tempered SMC over the Boltzmann posterior +/// `π_β(x) ∝ p(x)·exp(β·f(x))`, one rung per `step()`. +#[wasm_bindgen] +pub struct ExploreSmcInference { + model: EvolutionModel>, + particles: Vec, + rng: StdRng, + rung: usize, + n_rungs: usize, + beta: f64, + beta_max: f64, + log_evidence: f64, + rejuv_steps: usize, + crossover: bool, + resampled_last: bool, + swaps_last: usize, +} + +#[wasm_bindgen] +impl ExploreSmcInference { + /// `pop_size` particles, a linear β ladder of `n_rungs` rungs from 0 to + /// `beta_max` (1.0 = the posterior; > 1 keeps annealing into optimizer + /// territory), optional population crossover, explicit seed. + #[wasm_bindgen(constructor)] + pub fn new( + pop_size: usize, + n_rungs: usize, + beta_max: f64, + crossover: bool, + seed: u64, + ) -> Result { + let pop_size = pop_size.clamp(32, 600); + let n_rungs = n_rungs.clamp(4, 60); + let beta_max = if beta_max.is_finite() { + beta_max.clamp(1.0, 8.0) + } else { + 1.0 + }; + let model = EvolutionModel::new(GaussianPrior::new(0.0, PRIOR_STD, 2), TwinPeaks); + let mut rng = StdRng::seed_from_u64(seed); + let model_fn = model.smc_model(); + let mut particles = smc_prior_particles(&mut rng, pop_size, &model_fn); + // β = 0 semantics: the ladder starts at the prior with uniform + // weights (smc_prior_particles sets β = 1 importance weights). + let inv_n = 1.0 / pop_size as f64; + for p in &mut particles { + p.weight = inv_n; + p.log_weight = inv_n.ln(); + } + drop(model_fn); + Ok(ExploreSmcInference { + model, + particles, + rng, + rung: 0, + n_rungs, + beta: 0.0, + beta_max, + log_evidence: 0.0, + rejuv_steps: 3, + crossover, + resampled_last: false, + swaps_last: 0, + }) + } + + /// Current state without stepping (for the initial paint). + pub fn snapshot(&self) -> String { + self.state_json() + } + + /// Advance one tempering rung: incremental reweight by `Δβ·loglik`, + /// evidence accumulation (β ≤ 1 portion only), normalize, ESS-triggered + /// systematic resampling, π_β-invariant MH rejuvenation, and (optionally) + /// a crossover sweep — all through fugue's real primitives. + pub fn step(&mut self) -> String { + if self.rung >= self.n_rungs { + return self.state_json(); + } + self.rung += 1; + let new_beta = self.beta_max * self.rung as f64 / self.n_rungs as f64; + let d_beta = new_beta - self.beta; + let n = self.particles.len() as f64; + + // Evidence increment covers only the β ≤ 1 portion of this rung + // (log Z is defined at the posterior; annealing past 1 adds nothing). + let d_cap = new_beta.min(1.0) - self.beta.min(1.0); + if d_cap > 0.0 { + let lw: Vec = self + .particles + .iter() + .map(|p| p.log_weight + d_cap * loglik(&p.trace)) + .collect(); + self.log_evidence += log_sum_exp(&lw); + } + + // (1) incremental reweight + normalize (keeping log_weight in sync — + // fugue's normalize_particles updates only the linear weights). + for p in &mut self.particles { + p.log_weight += d_beta * loglik(&p.trace); + } + normalize_particles(&mut self.particles); + for p in &mut self.particles { + p.log_weight = if p.weight > 0.0 { + p.weight.ln() + } else { + f64::NEG_INFINITY + }; + } + + // (2) ESS-triggered systematic resampling. + let ess = effective_sample_size(&self.particles); + self.resampled_last = ess < 0.5 * n; + if self.resampled_last { + self.particles = + resample_particles(&mut self.rng, &self.particles, ResamplingMethod::Systematic); + } + + // (3) π_β-invariant rejuvenation + optional crossover sweep. + let model_fn = self.model.smc_model(); + rejuvenate_particles( + &mut self.rng, + &mut self.particles, + &model_fn, + new_beta, + self.rejuv_steps, + ); + self.swaps_last = 0; + if self.crossover && self.particles.len() >= 2 { + let before: Vec> = self + .particles + .iter() + .map(|p| p.trace.get_f64(&addr!("gene", 0usize))) + .collect(); + let mut kernel = CrossoverKernel { + n_pairs: self.particles.len() / 2, + // Swap the x-coordinate block between the pair: a + // product-target Metropolis move that exchanges mode + // membership horizontally. + mask: Box::new(|_: &Trace, _: &Trace, _: &mut dyn rand::RngCore| { + vec![addr!("gene", 0usize)] + }), + }; + PopulationKernel::::sweep( + &mut kernel, + &mut self.rng as &mut dyn rand::RngCore, + &mut self.particles, + &model_fn, + new_beta, + ); + self.swaps_last = self + .particles + .iter() + .zip(&before) + .filter(|(p, b)| p.trace.get_f64(&addr!("gene", 0usize)) != **b) + .count(); + } + + self.beta = new_beta; + self.state_json() + } + + /// The analytic tempered target on an `nx × ny` grid over the plot + /// domain, as **negative** log-density `−(log p + β·f)` (lower = better, + /// matching the landscape-heat convention), row-major with `j` indexing + /// `y` upward and cell-centered sampling. + pub fn density_grid(&self, nx: usize, ny: usize) -> Vec { + let nx = nx.clamp(2, 400); + let ny = ny.clamp(2, 400); + let (lo, hi) = DOMAIN; + let mut out = Vec::with_capacity(nx * ny); + for j in 0..ny { + let y = lo + (hi - lo) * (j as f64 + 0.5) / ny as f64; + for i in 0..nx { + let x = lo + (hi - lo) * (i as f64 + 0.5) / nx as f64; + out.push(-(log_prior(x, y) + self.beta * twin_peaks_f(x, y))); + } + } + out + } + + /// Plot metadata: `{lo, hi, prior_std, modes: [[x, y]…]}`. + pub fn info(&self) -> String { + json!({ + "lo": DOMAIN.0, + "hi": DOMAIN.1, + "prior_std": PRIOR_STD, + "modes": MODES.iter().map(|(m, _, _)| vec![m[0], m[1]]).collect::>(), + }) + .to_string() + } + + #[wasm_bindgen(getter)] + pub fn generation(&self) -> usize { + self.rung + } + + fn state_json(&self) -> String { + let pts: Vec = self + .particles + .iter() + .map(|p| { + let x = p.trace.get_f64(&addr!("gene", 0usize)).unwrap_or(0.0); + let y = p.trace.get_f64(&addr!("gene", 1usize)).unwrap_or(0.0); + json!([x, y, p.weight]) + }) + .collect(); + json!({ + "rung": self.rung, + "n_rungs": self.n_rungs, + "beta": self.beta, + "done": self.rung >= self.n_rungs, + "ess": effective_sample_size(&self.particles), + "resampled": self.resampled_last, + "swaps": self.swaps_last, + "log_evidence": self.log_evidence, + "particles": pts, + }) + .to_string() + } +} + +#[cfg(all(test, not(target_arch = "wasm32")))] +mod tests { + use super::*; + + fn transcript(seed: u64) -> String { + let mut e = ExploreSmcInference::new(200, 10, 1.0, true, seed).unwrap(); + let mut out = e.snapshot(); + for _ in 0..10 { + out.push_str(&e.step()); + } + out + } + + #[test] + fn smc_inference_is_deterministic() { + assert_eq!(transcript(7), transcript(7)); + assert_ne!(transcript(7), transcript(8)); + } + + #[test] + fn smc_inference_ladder_completes_and_shapes_hold() { + let mut e = ExploreSmcInference::new(150, 8, 1.0, false, 11).unwrap(); + for r in 1..=8 { + let v: serde_json::Value = serde_json::from_str(&e.step()).unwrap(); + assert_eq!(v["rung"].as_u64().unwrap(), r); + assert_eq!(v["particles"].as_array().unwrap().len(), 150); + let ess = v["ess"].as_f64().unwrap(); + assert!((1.0..=150.0 + 1e-9).contains(&ess)); + let beta = v["beta"].as_f64().unwrap(); + assert!((beta - r as f64 / 8.0).abs() < 1e-12); + for p in v["particles"].as_array().unwrap() { + assert!(p[0].as_f64().unwrap().is_finite()); + assert!(p[1].as_f64().unwrap().is_finite()); + } + } + let v: serde_json::Value = serde_json::from_str(&e.step()).unwrap(); + assert!(v["done"].as_bool().unwrap()); + assert_eq!(v["rung"].as_u64().unwrap(), 8); + } + + /// The widget shows real inference: at β = 1 the particle population's + /// weighted mean must match the analytic posterior mean of the twin-peaks + /// Boltzmann target (computed by grid quadrature of the exact density). + #[test] + fn smc_inference_matches_analytic_posterior_mean() { + // Grid quadrature of π_1 ∝ exp(log_prior + f). + let n = 400; + let (lo, hi) = DOMAIN; + let (mut z, mut ex, mut ey) = (0.0f64, 0.0f64, 0.0f64); + for j in 0..n { + let y = lo + (hi - lo) * (j as f64 + 0.5) / n as f64; + for i in 0..n { + let x = lo + (hi - lo) * (i as f64 + 0.5) / n as f64; + let d = (log_prior(x, y) + twin_peaks_f(x, y)).exp(); + z += d; + ex += d * x; + ey += d * y; + } + } + let (ex, ey) = (ex / z, ey / z); + + let mut e = ExploreSmcInference::new(500, 16, 1.0, true, 42).unwrap(); + let mut last = String::new(); + for _ in 0..16 { + last = e.step(); + } + let v: serde_json::Value = serde_json::from_str(&last).unwrap(); + let (mut mx, mut my, mut tw) = (0.0f64, 0.0f64, 0.0f64); + for p in v["particles"].as_array().unwrap() { + let w = p[2].as_f64().unwrap(); + mx += w * p[0].as_f64().unwrap(); + my += w * p[1].as_f64().unwrap(); + tw += w; + } + assert!((tw - 1.0).abs() < 1e-6, "weights self-normalise"); + let (mx, my) = (mx / tw, my / tw); + assert!( + (mx - ex).abs() < 0.3 && (my - ey).abs() < 0.3, + "particle mean ({mx:.3}, {my:.3}) vs analytic posterior mean ({ex:.3}, {ey:.3})" + ); + + // Log-evidence agrees with the quadrature estimate of + // log ∫ p(x)·e^f dx / ∫ p(x) dx (both densities unnormalised the + // same way, so compare against the ratio). + let mut zp = 0.0f64; + for j in 0..n { + let y = lo + (hi - lo) * (j as f64 + 0.5) / n as f64; + for i in 0..n { + let x = lo + (hi - lo) * (i as f64 + 0.5) / n as f64; + zp += log_prior(x, y).exp(); + } + } + let analytic_log_z = (z / zp).ln(); + let got = v["log_evidence"].as_f64().unwrap(); + assert!( + (got - analytic_log_z).abs() < 0.25, + "log evidence {got:.3} vs analytic {analytic_log_z:.3}" + ); + } + + #[test] + fn smc_inference_annealing_concentrates() { + // β_max = 6 tightens the population onto the peaks: for a multimodal + // target the honest concentration measure is the mean distance to the + // NEAREST mode (total spread stays dominated by the inter-mode + // distance as long as both peaks keep any mass). + let run_mode_dist = |beta_max: f64| -> f64 { + let mut e = ExploreSmcInference::new(300, 24, beta_max, true, 5).unwrap(); + let mut last = String::new(); + for _ in 0..24 { + last = e.step(); + } + let v: serde_json::Value = serde_json::from_str(&last).unwrap(); + let pts = v["particles"].as_array().unwrap(); + pts.iter() + .map(|p| { + let (x, y) = (p[0].as_f64().unwrap(), p[1].as_f64().unwrap()); + MODES + .iter() + .map(|(m, _, _)| ((x - m[0]).powi(2) + (y - m[1]).powi(2)).sqrt()) + .fold(f64::INFINITY, f64::min) + }) + .sum::() + / pts.len() as f64 + }; + let posterior = run_mode_dist(1.0); + let annealed = run_mode_dist(6.0); + assert!( + annealed < 0.6 * posterior, + "annealed mean mode-distance {annealed:.3} should be well below posterior {posterior:.3}" + ); + } +} diff --git a/crates/fugue-evo-wasm/src/lib.rs b/crates/fugue-evo-wasm/src/lib.rs index 95e96b1..4353a61 100644 --- a/crates/fugue-evo-wasm/src/lib.rs +++ b/crates/fugue-evo-wasm/src/lib.rs @@ -25,6 +25,7 @@ use wasm_bindgen::prelude::*; mod config; mod error; mod explore; +mod explore_inference; mod fitness; mod interactive; mod optimizers; @@ -33,6 +34,7 @@ mod result; pub use config::*; pub use error::*; pub use explore::*; +pub use explore_inference::*; pub use fitness::*; pub use interactive::*; pub use optimizers::*; diff --git a/docs/book.toml b/docs/book.toml index 4e384df..e4eeb45 100644 --- a/docs/book.toml +++ b/docs/book.toml @@ -24,6 +24,7 @@ additional-js = [ "viz/nsga2.js", "viz/island.js", "viz/inline.js", + "viz/inference.js", "viz/playground.js", ] diff --git a/docs/src/architecture/fugue-integration.md b/docs/src/architecture/fugue-integration.md index 738b1c5..0be3edb 100644 --- a/docs/src/architecture/fugue-integration.md +++ b/docs/src/architecture/fugue-integration.md @@ -16,6 +16,20 @@ is not merely a mathematical analogy — it is a [fugue](https://fugue.run) program, and every sampler in the layer is fugue's own inference machinery run against that program. +
+ +*This is the whole idea, live: fugue-evo's real inference layer (compiled to +WASM) runs tempered SMC against a twin-peaks Boltzmann target. The yellow heat +is the exact tempered density `π_β ∝ p(x)·exp(β·f(x))`, recomputed each rung +as β climbs from the prior (β = 0, the blue rings) to the posterior (β = 1); +the green particles are the SMC population — reweighted, resampled when the +ESS drops (watch the ESS readout flash), rejuvenated by typed MH, and mixed +across the two modes by the crossover kernel (a product-target Metropolis +swap of x-coordinates). Push β MAX past 1 to watch inference become an +annealed optimizer, and drag the seed to replay a different history. The +log Z readout is the running evidence estimate — a model score no classic GA +can report.* + ## Priors are programs The `GenomePrior` trait replaces any notion of a built-in prior enum: diff --git a/docs/viz/inference.js b/docs/viz/inference.js new file mode 100644 index 0000000..07fcce0 --- /dev/null +++ b/docs/viz/inference.js @@ -0,0 +1,432 @@ +/* Evolution-as-inference explorable: tempered SMC over the Boltzmann + * posterior, run by the REAL fugue-evo inference layer compiled to wasm + * (ExploreSmcInference). One tempering rung per tick: the heat is the exact + * tempered density -(log p + beta*f) recomputed from the crate each rung, so + * the particle cloud can be watched matching the analytic target as beta + * climbs from the prior (beta = 0) to the posterior (beta = 1) and onward + * into annealed-optimizer territory. No JS fallback math. */ +(function () { + "use strict"; + if (!window.FugueViz) return; + var FV = window.FugueViz; + + var NOTICE = + "This figure runs the real fugue-evo crate compiled to WebAssembly — the wasm package isn't available in this build."; + var HEAT_N = 110; + + /* ---- local helpers (per-file convention) ---- */ + function el(tag, cls, parent) { + var e = document.createElement(tag); + if (cls) e.className = cls; + if (parent) parent.appendChild(e); + return e; + } + function showNotice(root, msg) { + var d = el("div", "fv-pg-notice", root); + d.textContent = msg; + } + function mkPlot(w, h, dom) { + var pad = { l: 40, r: 12, t: 10, b: 26 }; + var ix = pad.l, + iy = pad.t, + iw = Math.max(10, w - pad.l - pad.r), + ih = Math.max(10, h - pad.t - pad.b); + return { + ix: ix, + iy: iy, + iw: iw, + ih: ih, + sx: FV.scale(dom, [ix, ix + iw]), + sy: FV.scale(dom, [iy + ih, iy]), + }; + } + function dot(ctx, x, y, r, color, alpha) { + ctx.save(); + ctx.globalAlpha = alpha; + ctx.fillStyle = color; + ctx.beginPath(); + ctx.arc(x, y, r, 0, Math.PI * 2); + ctx.fill(); + ctx.restore(); + } + function diamond(ctx, x, y, s, color) { + ctx.save(); + ctx.globalAlpha = 0.22; + ctx.fillStyle = color; + ctx.beginPath(); + ctx.arc(x, y, 9, 0, Math.PI * 2); + ctx.fill(); + ctx.globalAlpha = 0.95; + ctx.beginPath(); + ctx.moveTo(x, y - s); + ctx.lineTo(x + s, y); + ctx.lineTo(x, y + s); + ctx.lineTo(x - s, y); + ctx.closePath(); + ctx.fill(); + ctx.restore(); + } + function fmtF(v) { + if (!isFinite(v)) return "—"; + var a = Math.abs(v); + if (a >= 1e4 || (a < 1e-3 && a > 0)) return v.toExponential(2); + return v.toFixed(3); + } + function buildHeat(off, grid, nx, ny, colHex) { + off.width = nx; + off.height = ny; + var octx = off.getContext("2d"); + var img = octx.createImageData(nx, ny); + var mn = Infinity, + mx = -Infinity; + for (var k = 0; k < grid.length; k++) { + if (grid[k] < mn) mn = grid[k]; + if (grid[k] > mx) mx = grid[k]; + } + var span = mx - mn || 1; + var col = hex(colHex); + for (var j = 0; j < ny; j++) { + for (var i = 0; i < nx; i++) { + var t = (grid[j * nx + i] - mn) / span; + var a = 0.6 * (1 - Math.sqrt(t)); + var p = ((ny - 1 - j) * nx + i) * 4; + img.data[p] = col[0]; + img.data[p + 1] = col[1]; + img.data[p + 2] = col[2]; + img.data[p + 3] = Math.round(255 * a); + } + } + octx.putImageData(img, 0, 0); + } + function hex(h) { + h = (h || "#f2cc60").replace("#", ""); + if (h.length === 3) + h = h[0] + h[0] + h[1] + h[1] + h[2] + h[2]; + return [ + parseInt(h.slice(0, 2), 16) || 242, + parseInt(h.slice(2, 4), 16) || 204, + parseInt(h.slice(4, 6), 16) || 96, + ]; + } + function seedControl(controls, value, onInput) { + var wrap = el("div", "fv-control", controls); + var lab = el("span", "fv-control-label", wrap); + lab.textContent = "SEED"; + var span = el("span", "fv-control-value", wrap); + FV.scrub(span, { min: 1, max: 999, step: 1, value: value, onInput: onInput }); + } + + /* ---- widget ---- */ + function smcInit(root, W) { + var params = { + seed: parseInt(root.getAttribute("data-seed") || "11", 10) || 11, + pop: parseInt(root.getAttribute("data-pop") || "240", 10) || 240, + rungs: parseInt(root.getAttribute("data-rungs") || "18", 10) || 18, + betamax: parseFloat(root.getAttribute("data-betamax") || "1") || 1, + crossover: (root.getAttribute("data-crossover") || "on") !== "off", + }; + + var engine = null, + cur = null, + info = null, + DOM = [-4.5, 4.5]; + var heatCv = document.createElement("canvas"); + var heatDirty = true; + + function advance() { + cur = JSON.parse(engine.step()); + heatDirty = true; // the tempered density changes with beta + } + function rebuild() { + engine = new W.ExploreSmcInference( + params.pop, + params.rungs, + params.betamax, + params.crossover, + BigInt(params.seed) + ); + info = JSON.parse(engine.info()); + DOM = [info.lo, info.hi]; + cur = JSON.parse(engine.snapshot()); + heatDirty = true; + if (loopApi && loopApi.reduced) { + // Reduced motion: run the whole ladder synchronously, show the end. + while (!cur.done) advance(); + } + } + + /* controls -> canvas -> instruction -> readouts */ + var controls = el("div", "fv-controls", root); + FV.slider(controls, { + label: "PARTICLES", + min: 48, + max: 480, + step: 16, + value: params.pop, + fmt: function (v) { + return String(Math.round(v)); + }, + onInput: function (v) { + params.pop = Math.round(v); + rebuild(); + renderReadouts(); + requestDraw(); + }, + }); + FV.slider(controls, { + label: "RUNGS", + min: 6, + max: 40, + step: 1, + value: params.rungs, + fmt: function (v) { + return String(Math.round(v)); + }, + onInput: function (v) { + params.rungs = Math.round(v); + rebuild(); + renderReadouts(); + requestDraw(); + }, + }); + FV.slider(controls, { + label: "β MAX", + min: 1, + max: 6, + step: 0.5, + value: params.betamax, + fmt: function (v) { + return "×" + v.toFixed(1); + }, + onInput: function (v) { + params.betamax = v; + rebuild(); + renderReadouts(); + requestDraw(); + }, + }); + FV.toggle(controls, { + label: "CROSSOVER", + value: params.crossover, + onChange: function (v) { + params.crossover = v; + rebuild(); + renderReadouts(); + requestDraw(); + }, + }); + seedControl(controls, params.seed, function (v) { + params.seed = v | 0; + rebuild(); + renderReadouts(); + requestDraw(); + }); + var btns = FV.buttons(controls, [ + { + label: "Play", + title: "Run the tempering ladder", + primary: true, + onClick: togglePlay, + }, + { + label: "Step", + title: "Advance one tempering rung", + onClick: function () { + loopApi.step(); + }, + }, + { + label: "Reset", + title: "Rebuild from the current seed", + onClick: function () { + rebuild(); + renderReadouts(); + requestDraw(); + if (!loopApi.playing) setPlayLabel("Play"); + }, + }, + ]); + + var cv = FV.canvas(root, { + height: 300, + onResize: function () { + draw(); + }, + }); + var ctx = cv.ctx; + + var instr = el("div", "fv-instruction", root); + instr.textContent = + "yellow heat = the exact tempered target πβ ∝ p(x)·exp(β·f(x)) · green dots = SMC particles (size = weight) · blue rings = the prior's 1σ/2σ · coral diamonds = the two fitness peaks"; + + var readouts = el("div", "fv-readouts", root); + var rBeta = FV.readout(readouts, { label: "β" }); + var rRung = FV.readout(readouts, { label: "RUNG" }); + var rEss = FV.readout(readouts, { label: "ESS" }); + var rZ = FV.readout(readouts, { label: "log Z" }); + var rSwap = FV.readout(readouts, { label: "SWAPS" }); + + function renderReadouts() { + if (!cur) return; + rBeta.set(cur.beta.toFixed(2), "post"); + rRung.set(cur.rung + "/" + cur.n_rungs); + rEss.set(String(Math.round(cur.ess)), cur.resampled ? "hot" : "post"); + rZ.set(fmtF(cur.log_evidence), "data"); + rSwap.set(params.crossover ? String(cur.swaps) : "off", "flow"); + } + + function draw() { + if (!cv) return; + cv.clear(); + var th = FV.theme(), + C = th.colors; + var plot = mkPlot(cv.w, cv.h, DOM); + if (heatDirty && engine) { + buildHeat(heatCv, engine.density_grid(HEAT_N, HEAT_N), HEAT_N, HEAT_N, C.data); + heatDirty = false; + } + ctx.drawImage(heatCv, plot.ix, plot.iy, plot.iw, plot.ih); + FV.axes(ctx, { + x: plot.ix, + y: plot.iy, + w: plot.iw, + h: plot.ih, + xscale: plot.sx, + yscale: plot.sy, + xlabel: "x", + ylabel: "y", + theme: th, + }); + if (!cur || !info) return; + ctx.save(); + ctx.beginPath(); + ctx.rect(plot.ix, plot.iy, plot.iw, plot.ih); + ctx.clip(); + + // Prior 1-sigma / 2-sigma rings (the beta = 0 starting law). + var cx = plot.sx(0), + cyp = plot.sy(0); + var r1 = Math.abs(plot.sx(info.prior_std) - plot.sx(0)); + ctx.save(); + ctx.strokeStyle = C.prior; + ctx.globalAlpha = 0.5; + ctx.setLineDash([]); + ctx.beginPath(); + ctx.arc(cx, cyp, r1, 0, Math.PI * 2); + ctx.stroke(); + ctx.setLineDash([4, 4]); + ctx.globalAlpha = 0.3; + ctx.beginPath(); + ctx.arc(cx, cyp, 2 * r1, 0, Math.PI * 2); + ctx.stroke(); + ctx.restore(); + + // Fitness peaks. + for (var m = 0; m < info.modes.length; m++) { + diamond(ctx, plot.sx(info.modes[m][0]), plot.sy(info.modes[m][1]), 5, C.hot); + } + + // Particles: radius and alpha carry the normalized weight. + var n = cur.particles.length; + for (var i = 0; i < n; i++) { + var p = cur.particles[i]; + var w = p[2]; + var r = Math.min(6, 1.6 + Math.sqrt(Math.max(0, w) * n) * 1.6); + var a = Math.min(0.9, 0.25 + w * n * 0.35); + dot(ctx, plot.sx(p[0]), plot.sy(p[1]), r, C.post, a); + } + ctx.restore(); + } + + var drawQueued = false; + function requestDraw() { + if (loopApi.playing || drawQueued) return; + drawQueued = true; + window.requestAnimationFrame(function () { + drawQueued = false; + draw(); + }); + } + + var pacer = FV.pace(2.5); + var loopApi = FV.loop( + root, + function (dt) { + if (!engine) return; + if (dt === 0) { + if (!cur.done) advance(); + renderReadouts(); + draw(); + return; + } + if (cur.done) { + loopApi.pause(); + setPlayLabel("Replay"); + draw(); + return; + } + var ticks = pacer(dt); + while (ticks-- > 0 && !cur.done) advance(); + renderReadouts(); + draw(); + }, + { autoplay: true } + ); + + function setPlayLabel(txt) { + btns.fvButtons["Play"].textContent = txt; + } + function togglePlay() { + if (loopApi.playing) { + loopApi.pause(); + setPlayLabel("Play"); + } else { + if (cur && cur.done) { + rebuild(); + renderReadouts(); + } + loopApi.play(); + if (loopApi.playing) setPlayLabel("Pause"); + } + } + if (loopApi.reduced) { + btns.fvButtons["Play"].disabled = true; + btns.fvButtons["Play"].title = + "Reduced motion is on — the completed ladder is shown"; + btns.fvButtons["Step"].disabled = true; + } + + FV.onThemeChange(function () { + heatDirty = true; + draw(); + }); + + rebuild(); + renderReadouts(); + draw(); + if (loopApi.playing) setPlayLabel("Pause"); + } + + FV.register("smc-inference", function (root) { + var p = FV.wasmReady || Promise.resolve(null); + p.then(function (W) { + if (!W) { + root.setAttribute("data-fugue-backend", "none"); + showNotice(root, NOTICE); + return; + } + root.setAttribute("data-fugue-backend", "wasm"); + try { + smcInit(root, W); + } catch (e) { + try { + console.error("[fugue-viz] smc-inference init failed", e); + } catch (e2) {} + showNotice( + root, + "smc-inference failed to initialize: " + (e && e.message ? e.message : e) + ); + } + }); + }); +})(); diff --git a/src/inference/pareto.rs b/src/inference/pareto.rs index 977f0a9..41229f8 100644 --- a/src/inference/pareto.rs +++ b/src/inference/pareto.rs @@ -11,17 +11,28 @@ //! π_β(x, w) ∝ p(x) · exp(−β · ⟨w, f(x)⟩) //! ``` //! -//! the joint posterior's marginal over genomes traces out the **Pareto -//! front**: each weight vector `w` selects a scalarized optimum on the front, -//! and integrating over `w` spreads the population across it. Each particle's -//! trace carries its own `w` (at `pareto#v{i}` stick-breaking sites), telling -//! you *where on the front* that particle lives — a posterior over the front, -//! with the usual inference dividends (uncertainty, evidence), which NSGA-II -//! cannot express. +//! the joint posterior spreads over front-adjacent configurations: each +//! weight vector `w` selects a scalarized optimum on the front, and each +//! particle's trace carries its own `w` (at `pareto#v{i}` stick-breaking +//! sites), telling you *where on the front* that particle lives — a posterior +//! over front positions, with the usual inference dividends (uncertainty, +//! evidence), which NSGA-II cannot express. //! -//! Weighted-sum scalarization recovers the convex part of the front; for -//! non-convex fronts a Chebyshev scalarization would be needed (future work, -//! same architecture). +//! **Marginal-tilt caveat (read this)**: in the latent-`w` model the +//! `w`-marginal is *not* uniform — it is tilted by `exp(−s·m(w))`, where +//! `m(w)` is the scalarized optimum's value at `w`, so weights whose optima +//! score better attract more mass, and high sharpness or heavy annealing +//! concentrates the population near the best-scoring front regions (often the +//! endpoints). The *conditional* `x | w` is what tracks the front. For +//! uniform front coverage, sweep **fixed** weights +//! ([`ChebyshevScalarization::with_weight`]) across a grid, or keep sharpness +//! moderate and read positions off `particle_weights`. +//! +//! [`ParetoScalarization`] uses weighted-sum scalarization, which recovers +//! the convex part of the front; [`ChebyshevScalarization`] uses the weighted +//! Chebyshev (weighted-max) norm, which reaches every (weakly) +//! Pareto-optimal point — including non-convex front regions where every +//! weighted-sum optimum collapses to the front's endpoints. use fugue::{addr, factor, Beta, Model, ModelExt, Trace}; @@ -99,6 +110,91 @@ where } } +/// The Chebyshev (weighted-max) scalarization likelihood with a latent +/// weight vector: +/// +/// ```text +/// w ~ Uniform(simplex) +/// π_β(x, w) ∝ p(x) · exp(−β · s · max_i w_i · (f_i(x) − z_i)) +/// ``` +/// +/// where `z` is the **ideal point** (a reference component-wise ≤ the +/// objective values of interest, e.g. per-objective minima or a slightly +/// optimistic estimate). Minimizing the weighted Chebyshev norm over `x` +/// reaches every weakly Pareto-optimal point as `w` varies over the simplex +/// (Miettinen 1999) — in particular the **non-convex** front regions where a +/// weighted sum's interior stationary point is a maximum and all its mass +/// collapses onto the front's endpoints. Use this when the front may be +/// non-convex; use [`ParetoScalarization`] when it is known convex (the +/// weighted sum is smoother). +#[derive(Clone)] +pub struct ChebyshevScalarization { + /// The multi-objective fitness (objectives **minimized**). + pub objectives: M, + /// Sharpness of the scalarized likelihood (see [`ParetoScalarization`]). + pub sharpness: f64, + /// The ideal/reference point `z` (one entry per objective). + pub ideal: Vec, + /// `None`: the weight is a latent site (subject to the marginal-tilt + /// caveat in the [module docs](self)). `Some(w)`: a fixed weight — the + /// posterior concentrates on that weight's own front point, which is the + /// mode to use for sweeping the front uniformly. + pub weight: Option>, +} + +impl ChebyshevScalarization { + /// Create a Chebyshev-scalarization likelihood with a **latent** weight. + pub fn new(objectives: M, sharpness: f64, ideal: Vec) -> Self { + Self { + objectives, + sharpness, + ideal, + weight: None, + } + } + + /// Fix the scalarization weight (front-sweeping mode): the posterior + /// targets this weight's own scalarized optimum — reaching interior + /// points of non-convex fronts that no weighted sum can select. + pub fn with_weight(mut self, weight: Vec) -> Self { + self.weight = Some(weight); + self + } +} + +impl GenomeLikelihood for ChebyshevScalarization +where + G: 'static, + M: MultiObjectiveFitness + Clone + Send + Sync + 'static, +{ + fn model(&self, genome: &G, beta: f64) -> Model<()> { + let objs = self.objectives.evaluate(genome); + let k = objs.len(); + let sharpness = self.sharpness; + let ideal = self.ideal.clone(); + if k == 0 { + return fugue::pure(()); + } + debug_assert_eq!(ideal.len(), k, "ideal point must match objective count"); + let cheby_factor = move |w: &[f64], objs: &[f64], ideal: &[f64]| -> Model<()> { + let cheby = w + .iter() + .zip(objs.iter().zip(ideal)) + .map(|(wi, (fi, zi))| wi * (fi - zi)) + .fold(f64::NEG_INFINITY, f64::max); + if cheby.is_finite() { + factor(-beta * sharpness * cheby) + } else { + factor(f64::NEG_INFINITY) + } + }; + match self.weight.clone() { + Some(w) => cheby_factor(&w, &objs, &ideal), + None => weight_model(k).bind(move |w| cheby_factor(&w, &objs, &ideal)), + } + } +} + /// Read a particle's weight vector back off its trace (the stick-breaking /// sites), i.e. *where on the front* the particle lives. Returns `None` when /// the sites are absent (e.g. a prior-only trace). @@ -209,6 +305,182 @@ mod tests { ); } + /// The non-convex-front contrast, at the theorem level. Objectives + /// (minimized) on x ∈ [0,1]: `f1 = x`, `f2 = 1 − x²`. Every x in [0,1] + /// is Pareto-optimal and the front `f2 = 1 − f1²` is CONCAVE, so for any + /// FIXED weight the weighted-sum scalarization `w·x + (1−w)(1−x²)` has + /// its interior stationary point as a MAXIMUM (second derivative + /// −2(1−w) < 0): its minimizers are always the endpoints, and interior + /// front points are unreachable. The Chebyshev scalarization's fixed-w + /// optimum is the interior crossing point `w·x = (1−w)(1−x²)` — for + /// w = 1/2, x* = (√5−1)/2 ≈ 0.618. We pin both facts. + #[test] + fn test_chebyshev_reaches_nonconvex_front_where_weighted_sum_cannot() { + #[derive(Clone)] + struct ConcaveFront; + impl MultiObjectiveFitness for ConcaveFront { + fn num_objectives(&self) -> usize { + 2 + } + fn evaluate(&self, g: &RealVector) -> Vec { + let x = g.genes()[0]; + vec![x, 1.0 - x * x] + } + } + + /// Test-local fixed-weight weighted-sum likelihood (the published + /// ParetoScalarization is latent-w only). + #[derive(Clone)] + struct FixedWeightSum { + w: f64, + sharpness: f64, + } + impl GenomeLikelihood for FixedWeightSum { + fn model(&self, g: &RealVector, beta: f64) -> Model<()> { + let objs = ConcaveFront.evaluate(g); + let s = self.w * objs[0] + (1.0 - self.w) * objs[1]; + factor(-beta * self.sharpness * s) + } + } + + let prior = || UniformBoxPrior::new(MultiBounds::new(vec![Bounds::new(0.0, 1.0)])); + let cfg = || EvoSmcConfig { + num_particles: 500, + ess_threshold: 0.5, + resampling: ResamplingMethod::Systematic, + rejuvenation_steps: 5, + crossover: Some(CrossoverConfig::default()), + }; + + let mut rng = StdRng::seed_from_u64(271828); + + // (a) Fixed w = 1/2, weighted sum: bimodal at the endpoints; the + // interior is a scalarization MAXIMUM and must be avoided. + let ws_model = EvolutionModel::from_likelihood( + prior(), + FixedWeightSum { + w: 0.5, + sharpness: 25.0, + }, + ); + let ws = EvolutionSMC::anneal(&mut rng, &ws_model, cfg(), 8.0, 8); + let ws_fn = ws_model.smc_model(); + let ws_interior: f64 = fugue::decode_particles(&ws.particles, &ws_fn) + .iter() + .filter(|(g, _)| (0.25..0.75).contains(&g.genes()[0])) + .map(|(_, w)| w) + .sum(); + assert!( + ws_interior < 0.1, + "fixed-w weighted sum put {ws_interior:.3} mass in the interior — impossible for a concave front" + ); + + // (b) Fixed w = 1/2, Chebyshev: concentrates on the interior front + // point x* = (√5 − 1)/2 ≈ 0.618 — the point weighted-sum cannot reach. + let x_star = (5.0f64.sqrt() - 1.0) / 2.0; + let ch_model = EvolutionModel::from_likelihood( + prior(), + ChebyshevScalarization::new(ConcaveFront, 25.0, vec![0.0, 0.0]) + .with_weight(vec![0.5, 0.5]), + ); + let ch = EvolutionSMC::anneal(&mut rng, &ch_model, cfg(), 8.0, 8); + let mean = ch.weighted_mean(0); + assert!( + (mean - x_star).abs() < 0.08, + "fixed-w Chebyshev posterior mean {mean:.3} should sit at the interior front point {x_star:.3}" + ); + let ch_fn = ch_model.smc_model(); + let ch_interior: f64 = fugue::decode_particles(&ch.particles, &ch_fn) + .iter() + .filter(|(g, _)| (0.25..0.75).contains(&g.genes()[0])) + .map(|(_, w)| w) + .sum(); + assert!( + ch_interior > 0.8, + "fixed-w Chebyshev interior mass {ch_interior:.3} — must reach the non-convex front interior" + ); + + // (c) Sweeping fixed weights traces the whole front, ends included. + for (w, lo, hi) in [(0.15, 0.75, 1.0), (0.5, 0.5, 0.75), (0.85, 0.1, 0.45)] { + let m = EvolutionModel::from_likelihood( + prior(), + ChebyshevScalarization::new(ConcaveFront, 25.0, vec![0.0, 0.0]) + .with_weight(vec![w, 1.0 - w]), + ); + let r = EvolutionSMC::anneal(&mut rng, &m, cfg(), 8.0, 8); + let mean = r.weighted_mean(0); + assert!( + (lo..=hi).contains(&mean), + "weight {w}: front point {mean:.3} outside expected band [{lo}, {hi}]" + ); + } + } + + /// Latent-weight Chebyshev: the CONDITIONAL x | w tracks the front even + /// though the w-marginal is tilted (module-docs caveat). Among particles + /// whose latent weight is interior (w₀ ∈ [0.35, 0.65]), most mass must + /// sit in the interior of the front — the region a weighted sum's + /// conditional never occupies. + #[test] + fn test_chebyshev_latent_weight_conditional_tracks_front() { + #[derive(Clone)] + struct ConcaveFront; + impl MultiObjectiveFitness for ConcaveFront { + fn num_objectives(&self) -> usize { + 2 + } + fn evaluate(&self, g: &RealVector) -> Vec { + let x = g.genes()[0]; + vec![x, 1.0 - x * x] + } + } + + let prior = UniformBoxPrior::new(MultiBounds::new(vec![Bounds::new(0.0, 1.0)])); + let model = EvolutionModel::from_likelihood( + prior, + ChebyshevScalarization::new(ConcaveFront, 8.0, vec![0.0, 0.0]), + ); + let mut rng = StdRng::seed_from_u64(314159); + // β = 1 posterior only — annealing would concentrate the tilted + // w-marginal onto the endpoints (see module docs). + let result = EvolutionSMC::run( + &mut rng, + &model, + EvoSmcConfig { + num_particles: 800, + ess_threshold: 0.5, + resampling: ResamplingMethod::Systematic, + rejuvenation_steps: 6, + crossover: Some(CrossoverConfig::default()), + }, + ); + let model_fn = model.smc_model(); + let decoded = fugue::decode_particles(&result.particles, &model_fn); + + let mut stratum_mass = 0.0; + let mut stratum_interior = 0.0; + for (p, (g, w)) in result.particles.iter().zip(&decoded) { + if let Some(wv) = particle_weights(&p.trace, 2) { + if (0.35..=0.65).contains(&wv[0]) { + stratum_mass += w; + let x = g.genes()[0]; + if (0.25..0.75).contains(&x) { + stratum_interior += w; + } + } + } + } + assert!( + stratum_mass > 0.02, + "interior-weight stratum carries only {stratum_mass:.4} mass — too depleted to test" + ); + let frac = stratum_interior / stratum_mass; + assert!( + frac > 0.5, + "interior-weight particles put only {frac:.2} of their mass on the front interior" + ); + } + /// Stick-breaking weights are a valid distribution over the simplex for /// k = 3: components positive, summing to 1, with symmetric means. #[test]