Merge pull request #513 from AdaWorldAPI/claude/perturbation-sim-inertia-clam

feat(perturbation-sim): inertia §0 promotion gate + CAKES/CHAODA + witness standing-wave + H ingest

Author: AdaWorldAPI

crates/perturbation-sim/Cargo.toml

@@ -100,3 +100,15 @@ path = "examples/calibrate.rs"
 [[example]]
 name = "hhtl_grid"
 path = "examples/hhtl_grid.rs"
+
+[[example]]
+name = "chaoda"
+path = "examples/chaoda.rs"
+
+[[example]]
+name = "witness"
+path = "examples/witness.rs"
+
+[[example]]
+name = "inertia_ingest"
+path = "examples/inertia_ingest.rs"

crates/perturbation-sim/examples/chaoda.rs

@@ -0,0 +1,133 @@
+//! CAKES + CHAODA over the HHTL basins of a real grid core — the similarity
+//! (attraction) / anomaly (repulsion) pair from the CLAM family, scoring grid
+//! compartments the way CHAODA scores papillary muscles or terrain tiles.
+//!
+//! HHTL is the family basin ("where"); CAKES finds each basin's relatives ("who
+//! looks similar"); CHAODA scores how far a basin sits from its family ("why am I
+//! different") — the top-anomaly basin is the fail-first compartment.
+//!
+//! Run: cargo run --release --manifest-path crates/perturbation-sim/Cargo.toml \
+//!        --example chaoda -- /tmp/pypsa/buses.csv /tmp/pypsa/lines.csv ES
+//! (no args → a synthetic 8×8 grid with one deliberately weakened block).
+
+use perturbation_sim::{
+    anomaly_ranking, cakes_neighbors, chaoda_scores, resilience_basin_features, Edge, Grid,
+    HhtlKey, CHAODA_FLAG,
+};
+
+/// Synthetic 8×8 lattice with one corner block weakened (low-susceptance internal
+/// edges) — a planted fail-first compartment for the no-data demo.
+fn synthetic() -> Grid {
+    let (rows, cols) = (8usize, 8usize);
+    let id = |r: usize, c: usize| r * cols + c;
+    let mut e = Vec::new();
+    for r in 0..rows {
+        for c in 0..cols {
+            // The 3×3 top-left block is brittle (weak internal coupling).
+            let brittle = r < 3 && c < 3;
+            let b = if brittle { 0.05 } else { 1.0 };
+            if c + 1 < cols {
+                e.push(Edge::new(id(r, c), id(r, c + 1), b, 1.0));
+            }
+            if r + 1 < rows {
+                e.push(Edge::new(id(r, c), id(r + 1, c), b, 1.0));
+            }
+        }
+    }
+    Grid::new(rows * cols, e)
+}
+
+/// Deterministic per-bus inertia proxy (real `H` is not in PyPSA topology): a small
+/// fixed cycle, decoupled from wiring, so the buffer axis is an independent input —
+/// the structure (basins as families, the outlier as fail-first) holds regardless.
+fn proxy_inertia(n: usize) -> Vec<f64> {
+    (0..n).map(|i| 2.0 + (i % 5) as f64).collect()
+}
+
+fn fmt_key(k: &HhtlKey) -> String {
+    format!("{}.{}.{}", k.heel, k.hip, k.twig)
+}
+
+fn main() {
+    let args: Vec<String> = std::env::args().collect();
+    let grid = if args.len() >= 3 {
+        let buses = std::fs::read_to_string(&args[1]).expect("buses.csv");
+        let lines = std::fs::read_to_string(&args[2]).expect("lines.csv");
+        let cc = args.get(3).map(|s| s.as_str()).unwrap_or("ES");
+        let imp = perturbation_sim::from_pypsa_csv(&buses, &lines, Some(cc))
+            .expect("import")
+            .largest_component();
+        println!("grid: {cc} PyPSA core — {} buses", imp.grid.n);
+        imp.grid
+    } else {
+        let g = synthetic();
+        println!(
+            "grid: synthetic 8×8 with a planted brittle 3×3 block — {} buses",
+            g.n
+        );
+        g
+    };
+
+    let h = proxy_inertia(grid.n);
+    let (basins, rows) = resilience_basin_features(&grid, &h, 0.2);
+    let scores = chaoda_scores(&rows, 2);
+
+    println!("\n== HHTL family basins — CAKES (similar) + CHAODA (anomalous) ==");
+    println!(
+        "  {:>10}  {:>6}  {:>8}  {:>8}  {:>8}  flag",
+        "basin", "size", "λ₂n", "inertia", "CHAODA"
+    );
+    for (i, k) in basins.iter().enumerate() {
+        let flag = if scores[i] >= CHAODA_FLAG {
+            "◀ ANOMALY"
+        } else {
+            ""
+        };
+        println!(
+            "  {:>10}  {:>6.2}  {:>8.3}  {:>8.3}  {:>8.3}  {}",
+            fmt_key(k),
+            rows[i][1],
+            rows[i][0],
+            rows[i][2],
+            scores[i],
+            flag
+        );
+    }
+
+    // CHAODA: the fail-first compartment is the top anomaly.
+    let rank = anomaly_ranking(&rows, 2);
+    if let Some(&(top, score)) = rank.first() {
+        println!(
+            "\n  CHAODA → fail-first compartment: basin {} (score {:.3}{})\n  \
+             \"why am I different\" — the basin whose resilience profile deviates most\n  \
+             from its family. {} basins; flag threshold {CHAODA_FLAG}.",
+            fmt_key(&basins[top]),
+            score,
+            if score >= CHAODA_FLAG {
+                ", FLAGGED"
+            } else {
+                ""
+            },
+            basins.len()
+        );
+
+        // CAKES: the top-anomaly basin's nearest relatives ("who looks similar").
+        let nbrs = cakes_neighbors(&rows, top, 3);
+        let rel: Vec<String> = nbrs
+            .iter()
+            .map(|(i, d)| format!("{} (d={:.3})", fmt_key(&basins[*i]), d))
+            .collect();
+        println!(
+            "  CAKES → {}'s nearest relatives: [{}]\n  \
+             attraction vs repulsion: the family it resembles, and how far it still sits from them.",
+            fmt_key(&basins[top]),
+            rel.join(", ")
+        );
+    }
+
+    println!(
+        "\n  CAKES pulls in the similar; CHAODA pushes out the unusual.\n  \
+         Family basin (HHTL) + deviation-from-family (CHAODA) = the fail-first locator.\n  \
+         (CHAODA-lite kNN scorer; ndarray::clam's ClamTree ensemble is the gated production path.)"
+    );
+}

crates/perturbation-sim/examples/inertia_ingest.rs

@@ -0,0 +1,66 @@
+//! Probe 3 — the per-bus inertia (H) ingest path: parse a `bus,H` table, align it
+//! to the grid's bus order, disclose the proxy fallback, and feed the
+//! inertia-buffer column. No external data is bundled; this uses an inline fixture
+//! standing in for an ENTSO-E / ESIOS / TSO inertia export.
+//!
+//! Run: cargo run --release --manifest-path crates/perturbation-sim/Cargo.toml \
+//!        --example inertia_ingest
+
+use perturbation_sim::{
+    inertia_buffer_column, inertia_for_buses, parse_bus_inertia, proxy_inertia,
+};
+
+// Stand-in for an operator inertia export (real H is never committed).
+const INERTIA_CSV: &str = "\
+bus_id,inertia_h,source
+ES_NUC_1,6.5,nuclear
+ES_HYD_1,4.0,hydro
+ES_GAS_1,5.5,gas-synchronous
+ES_WND_1,0.8,wind-synthetic
+";
+
+fn main() {
+    // The grid's buses (would come from PypsaImport.bus_ids). ES_SOLAR_1 has no
+    // measured H → it falls back to the proxy value, and that is disclosed.
+    let bus_ids: Vec<String> = ["ES_NUC_1", "ES_HYD_1", "ES_GAS_1", "ES_WND_1", "ES_SOLAR_1"]
+        .iter()
+        .map(|s| s.to_string())
+        .collect();
+
+    let measured = parse_bus_inertia(INERTIA_CSV);
+    let fallback = 1.0; // low-inertia stand-in for an undocumented (likely inverter) bus
+    let (h, prov) = inertia_for_buses(&bus_ids, &measured, fallback);
+    let col = inertia_buffer_column(&h, 0.2);
+
+    println!("per-bus inertia ingest (measured H → grid order → buffer column)\n");
+    println!(
+        "  {:>12}  {:>8}  {:>10}  source",
+        "bus", "H (s)", "buffer_n"
+    );
+    for (i, id) in bus_ids.iter().enumerate() {
+        let src = if measured.contains_key(id) {
+            "measured"
+        } else {
+            "proxy"
+        };
+        println!("  {:>12}  {:>8.2}  {:>10.3}  {src}", id, h[i], col[i]);
+    }
+    println!(
+        "\n  provenance: {} measured, {} proxy (disclose the proxy fraction, like\n  \
+         PypsaImport's n_estimated_* counters).",
+        prov.measured, prov.proxy
+    );
+
+    // No-data path: a fully deterministic, topology-blind proxy field.
+    let demo = proxy_inertia(bus_ids.len(), 2.0, 6.0, 0xDEFA17);
+    println!(
+        "\n  no-data fallback — proxy_inertia(n=5, base=2, span=6, seed=0xDEFA17):\n  \
+         [{}]\n  decoupled from wiring on purpose (buffer ⊥ topology); deterministic, never\n  \
+         bundled. Wire measured H when an operator export is available; the column math\n  \
+         is identical either way.",
+        demo.iter()
+            .map(|h| format!("{h:.2}"))
+            .collect::<Vec<_>>()
+            .join(", ")
+    );
+}

crates/perturbation-sim/examples/witness.rs

@@ -0,0 +1,70 @@
+//! The witness arc as a standing wave (METHODS §11) — particle (pointer-chase) vs
+//! wave (Walsh-pyramid Parseval), on a real grid inertia field.
+//!
+//! Builds the per-bus inertia-buffer field, then reads three witness arcs both ways:
+//! the particle walk (`O(hops)` per arc) and the standing wave (one `field_spectrum`
+//! transform, then `O(N)` per arc). They agree to floating point — and the wave form
+//! amortizes one transform across all arcs.
+//!
+//! Run: cargo run --release --manifest-path crates/perturbation-sim/Cargo.toml \
+//!        --example witness
+
+use perturbation_sim::{
+    field_spectrum, inertia_buffer_column, witness_from_spectrum, witness_particle,
+};
+
+fn proxy_inertia(n: usize) -> Vec<f64> {
+    (0..n).map(|i| 2.0 + (i % 5) as f64).collect()
+}
+
+fn main() {
+    // The field: a per-bus inertia-buffer field over a 16-bus line (METHODS §11's
+    // "inertia field on power grids"; the promoted additive member as the carrier).
+    let n = 16;
+    let field = inertia_buffer_column(&proxy_inertia(n), 0.2);
+    let field: Vec<f64> = field.iter().map(|&x| x as f64).collect();
+
+    // Three witness arcs (Markov #1 reference chains over the buses):
+    //   - a single-hop witness (read one bus),
+    //   - a coarse dyadic arc (the first half — a Raumgewinn-scale reference),
+    //   - an alternating signed chain (a fine-scale infight reference).
+    let mut single = vec![0.0; n];
+    single[5] = 1.0;
+    let mut coarse = vec![0.0; n];
+    for c in coarse.iter_mut().take(n / 2) {
+        *c = 1.0;
+    }
+    let alt: Vec<f64> = (0..n)
+        .map(|i| if i % 2 == 0 { 1.0 } else { -1.0 })
+        .collect();
+    let arcs: [(&str, &[f64]); 3] = [
+        ("single-hop", &single),
+        ("coarse dyadic", &coarse),
+        ("alternating", &alt),
+    ];
+
+    // Wave view: transform the field ONCE; every arc is then an O(N) read off it.
+    let spectrum = field_spectrum(&field);
+
+    println!("witness arc as a standing wave — particle (walk) vs wave (Parseval)\n");
+    println!("  field: {n}-bus inertia-buffer field; one Walsh transform reused by all arcs\n");
+    println!(
+        "  {:>14}  {:>14}  {:>14}  {:>10}",
+        "arc", "particle", "wave", "Δ"
+    );
+    let mut max_err = 0.0_f64;
+    for (name, arc) in arcs {
+        let p = witness_particle(&field, arc);
+        let w = witness_from_spectrum(&spectrum, arc);
+        let d = (p - w).abs();
+        max_err = max_err.max(d);
+        println!("  {name:>14}  {p:>14.9}  {w:>14.9}  {d:>10.2e}");
+    }
+    println!(
+        "\n  max |particle − wave| = {max_err:.2e} (Parseval: Hᵀ H = N·I, exact up to fp).\n  \
+         particle = O(hops) pointer-chase per arc; wave = O(N log N) once + O(N) per arc.\n  \
+         The standing wave IS the witness arc — evaluated all at once, no chain walk.\n  \
+         (Demonstration in perturbation-sim; the contract witness_table evaluator is the\n  \
+         separate gated step — the SoA spine is additive-only behind the iron rules.)"
+    );
+}

crates/perturbation-sim/src/buffer.rs

@@ -74,10 +74,57 @@ pub fn compartment_buffer(inertia_h: &[f64], df_band: f64) -> f64 {
     inertia_h.iter().map(|&h| impulse_buffer(h, df_band)).sum()
 }
 
+/// The `inertia_buffer` SoA column for a set of buses: each bus's impulse buffer
+/// ([`impulse_buffer`]) min-max **normalized to `[0,1]`** — the form the calibrated
+/// [`crate::INERTIA`] spec stores (normalized, not raw physical units; normalizing
+/// also lifts the axis clear of the ICC variance-underflow guard). This is the
+/// promoted additive value member, *computed*. It is orthogonal to topology by the
+/// HHTL-OGAR key/value split (topology is the GUID key; this is one more value) —
+/// the structural fact `buffer_is_independent_of_connectivity` witnesses. Degenerate
+/// (all-equal `H`, or empty) input yields all-zeros, never `NaN`.
+pub fn inertia_buffer_column(per_bus_h: &[f64], df_band: f64) -> Vec<f32> {
+    let raw: Vec<f64> = per_bus_h
+        .iter()
+        .map(|&h| impulse_buffer(h, df_band))
+        .collect();
+    let lo = raw.iter().copied().fold(f64::INFINITY, f64::min);
+    let hi = raw.iter().copied().fold(f64::NEG_INFINITY, f64::max);
+    let span = hi - lo;
+    raw.iter()
+        .map(|&r| {
+            if span > 0.0 {
+                ((r - lo) / span) as f32
+            } else {
+                0.0
+            }
+        })
+        .collect()
+}
+
 #[cfg(test)]
 mod tests {
     use super::*;
 
+    #[test]
+    fn inertia_column_is_normalized_and_monotone_in_h() {
+        // Per-bus inertia decoupled from any wiring: the column normalizes to [0,1]
+        // and preserves the H order (impulse_buffer is monotone in H at fixed df).
+        let h = [1.0_f64, 5.0, 3.0, 0.0, 8.0];
+        let col = inertia_buffer_column(&h, 0.2);
+        assert_eq!(col.len(), h.len());
+        for &c in &col {
+            assert!((0.0..=1.0).contains(&c), "normalized to [0,1]");
+        }
+        assert_eq!(col[3], 0.0, "H=0 is the min → 0.0");
+        assert_eq!(col[4], 1.0, "H=8 is the max → 1.0");
+        // Order preserved: H[0]=1 < H[2]=3 < H[1]=5 < H[4]=8.
+        assert!(col[0] < col[2] && col[2] < col[1] && col[1] < col[4]);
+        // Degenerate (all-equal H → no span) yields zeros, never NaN.
+        let flat = inertia_buffer_column(&[4.0, 4.0, 4.0], 0.2);
+        assert!(flat.iter().all(|&c| c == 0.0));
+        assert!(inertia_buffer_column(&[], 0.2).is_empty());
+    }
+
     #[test]
     fn buffer_grows_with_inertia() {
         // Twice the inertia ⇒ twice the impulse absorbed (the storage scaling).