perturbation-sim: resilience arc — meta-hop → buffer → exploration → fail-first scorecard

crates/perturbation-sim/Cargo.toml

@@ -60,3 +60,35 @@ path = "examples/reinforce.rs"
 [[example]]
 name = "hhtl_resident"
 path = "examples/hhtl_resident.rs"
+
+[[example]]
+name = "meta_hops"
+path = "examples/meta_hops.rs"
+
+[[example]]
+name = "validate_mediators"
+path = "examples/validate_mediators.rs"
+
+[[example]]
+name = "rolling_floor"
+path = "examples/rolling_floor.rs"
+
+[[example]]
+name = "modifier_validity"
+path = "examples/modifier_validity.rs"
+
+[[example]]
+name = "resilience"
+path = "examples/resilience.rs"
+
+[[example]]
+name = "buffer"
+path = "examples/buffer.rs"
+
+[[example]]
+name = "explore"
+path = "examples/explore.rs"
+
+[[example]]
+name = "scorecard"
+path = "examples/scorecard.rs"

crates/perturbation-sim/METHODS.md

@@ -365,6 +365,52 @@ Reserve **electricity-specific** encoding for the **raw physical layer** (AC
 256-palette lossy, ±½ bucket — fine for a live screen/stream; compute exact
 stats on raw. `turbovec` is coarse — retrieval, not values.)
 
+## 11. Why the pyramid is a *computation* substrate: witness arc as a standing wave
+
+The perturbation pyramid is not only an outage model — it is the **computation
+kernel** for evaluating a Mailbox-SoA **witness arc** as a *standing wave* instead
+of by *pointer chasing*. This is the bridge the operator named ("use Mailbox SoA
+and witness > pointer chasing as a standing wave, but it needs the perturbation
+pyramid for computation"), and it lands on architecture that already exists in
+`lance-graph-contract`.
+
+**The two views of the same arc** (already documented in
+`contract::witness_table` and `contract::soa_view`):
+
+- **Particle view (pointer chasing).** The `CausalEdge64` W-slot → witness chain
+  is a Markov #1 reference arc; walking it means dereferencing one witness per hop
+  across the SPO store — `O(hops)` dependent loads, each a cache miss. This is the
+  *discrete, addressable, exact* witness pointer (`witness_table.rs`: "the chain of
+  W-references across edges forms a Markov chain … walked backwards without
+  dereferencing the full SPO store on every hop"; `soa_view.rs`: "the *particle*").
+
+- **Wave view (standing wave).** The same arc, evaluated *all at once* as a bipolar
+  **interference field** over the HHTL tiers — the windowed/resonance reading
+  (`soa_view.rs`: "the windowed … *wave*"). No per-hop dereference: the field is a
+  function of the address tree, not of a pointer walk.
+
+**The pyramid is what makes the wave computable** — and it is exactly
+`timing::meta_cascade_phase` (§4.7 of `PAPER.md`):
+
+| Standing-wave ingredient | Pyramid mechanism (this crate) | OGAR canon |
+|---|---|---|
+| phase (±1) composes between levels | `phase_{i+1} = phaseᵢ · sign` (XOR/multiply) | "sign side = `vsa_bind` = XOR" |
+| magnitude bundles into the field | `field_k = Σ phaseᵢ·magnitudeᵢ` (running sum) | "magnitude side = `vsa_bundle` = add" |
+| the field is the Walsh–Hadamard of the tree | `sketch::fwht` / `walsh_pyramid_energy` (`ndarray::simd::wht_f32`) | "Bipolar-phase pyramid — Walsh-Hadamard on VSA" |
+| tier = one meta-hop, `O(tiers)` not `O(hops)` | 4 HHTL tiers, `tier = level >> 2` | "perturbation pyramid for computation" |
+
+So the witness arc that the particle view walks in `O(hops)` dependent loads, the
+wave view reads in `O(tiers)` (≤ 4) tier lookups — *because the deterministic
+bipolar phase is generated from the address, never stored* (OGAR "DETERMINISTIC
+PHASE"). The standing wave is the witness arc; the perturbation pyramid is its
+evaluator; the SoA columns (`FingerprintColumns`/`EdgeColumn`) are its backing
+store. **Scope (honest):** the bridge is *structural* — `perturbation-sim`
+demonstrates the pyramid/phase/inertia field on power grids; wiring it as the
+actual `witness`-arc evaluator in the contract is a separate, gated step (the
+witness/SoA types are the cognitive spine — additive only, behind the iron rules).
+It is recorded here as the computation-substrate connection, not yet as shipped
+contract code. CONJECTURE [H].
+
 ## Anti-dilution table — the distinctions to never collapse
 
 | Do NOT conflate | Because |

crates/perturbation-sim/examples/buffer.rs

@@ -0,0 +1,252 @@
+//! The buffer axis on the real ES core — impulse storage before yield, and the
+//! proof it is INDEPENDENT of the resistive (Kirchhoff) axis the modifier was
+//! confounded with.
+//!
+//! Per Cheeger compartment we read two resilience axes. STEADY / resistive: λ₂,
+//! Kf, mean pairwise resistance (topology) — §4.10. TRANSIENT / buffer: the
+//! sudden imbalance the compartment's inertia absorbs before a Ketchup yield
+//! (storage). Size-normalized, mean-buffer (per node = mean inertia,
+//! topology-FREE) and mean-resistance (per pair, topology) are independent by
+//! construction — the deconfound: `1/λ₂` was the dominant Kirchhoff term, the
+//! buffer is not.
+//!
+//! Then the Kugelstoßpendel demo: a sudden per-node impulse is metered against
+//! each compartment's buffer/node (mean inertia headroom); the thinnest-buffer
+//! compartment yields first (the Ketchup collapse) even when it is resistively
+//! healthy (high λ₂) — the low-inertia failure mode a resistance-only screen misses.
+//!
+//! Inertia `H` is NOT in the PyPSA topology, so it is an ILLUSTRATIVE per-node
+//! scenario field here (deterministic, topology-independent, a ~third of buses
+//! marked renewable-rich = low inertia). The *structure* (buffer ⊥ connectivity,
+//! low-buffer-yields-first) holds for any real `H`; feed measured inertia to calibrate.
+//!
+//! Run: cargo run --release --manifest-path crates/perturbation-sim/Cargo.toml \
+//!        --example buffer -- /tmp/pypsa/buses.csv /tmp/pypsa/lines.csv ES
+
+use perturbation_sim::{
+    cheeger_sweep, compartment_buffer, ketchup_yield, spearman, symmetric_eigen, Edge, Grid,
+    Resilience,
+};
+use std::collections::HashMap;
+
+fn synthetic(rows: usize, cols: usize) -> Grid {
+    let id = |r: usize, c: usize| r * cols + c;
+    let mut e = Vec::new();
+    for r in 0..rows {
+        for c in 0..cols {
+            if c + 1 < cols {
+                e.push(Edge::new(id(r, c), id(r, c + 1), 1.0, 1.0));
+            }
+            if r + 1 < rows {
+                e.push(Edge::new(id(r, c), id(r + 1, c), 1.0, 1.0));
+            }
+        }
+    }
+    Grid::new(rows * cols, e)
+}
+
+fn induced(grid: &Grid, members: &[usize]) -> Grid {
+    let mut remap = HashMap::new();
+    for (i, &m) in members.iter().enumerate() {
+        remap.insert(m, i);
+    }
+    let edges = grid
+        .edges
+        .iter()
+        .filter_map(|e| match (remap.get(&e.from), remap.get(&e.to)) {
+            (Some(&a), Some(&b)) => Some(Edge::new(a, b, e.susceptance, e.limit)),
+            _ => None,
+        })
+        .collect();
+    Grid::new(members.len(), edges)
+}
+
+fn bisect(grid: &Grid, members: &[usize]) -> Option<(Vec<usize>, Vec<usize>)> {
+    if members.len() < 6 {
+        return None;
+    }
+    let sub = induced(grid, members);
+    let c = cheeger_sweep(&sub, &vec![true; sub.edges.len()]);
+    let (mut a, mut b) = (Vec::new(), Vec::new());
+    for (i, &m) in members.iter().enumerate() {
+        if c.partition[i] {
+            a.push(m);
+        } else {
+            b.push(m);
+        }
+    }
+    if a.is_empty() || b.is_empty() {
+        None
+    } else {
+        Some((a, b))
+    }
+}
+
+/// Deterministic topology-independent inertia field (ILLUSTRATIVE, see header):
+/// most buses synchronous H∈[3,7] s; ~⅓ marked renewable-rich at H=1 s.
+fn inertia_scenario(n: usize) -> Vec<f64> {
+    (0..n)
+        .map(|i| {
+            let mut z = (i as u64).wrapping_mul(0x9E37_79B9_7F4A_7C15);
+            z = (z ^ (z >> 30)).wrapping_mul(0xBF58_476D_1CE4_E5B9);
+            let u = ((z ^ (z >> 31)) >> 11) as f64 / (1u64 << 53) as f64;
+            if u < 0.33 {
+                1.0 // renewable-rich: low inertia, thin buffer
+            } else {
+                3.0 + 4.0 * u // synchronous: H∈[3,7]
+            }
+        })
+        .collect()
+}
+
+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, {} lines",
+            imp.grid.n,
+            imp.grid.edges.len()
+        );
+        imp.grid
+    } else {
+        let g = synthetic(10, 10);
+        println!("grid: synthetic 10×10 — {} buses", g.n);
+        g
+    };
+    let n = grid.n;
+    let df_band = 0.2; // Hz protection band
+    let inertia = inertia_scenario(n);
+
+    // Depth-3 Cheeger compartments.
+    let mut basins: Vec<Vec<usize>> = vec![(0..n).collect()];
+    for _ in 0..3 {
+        let mut next = Vec::new();
+        for b in &basins {
+            match bisect(&grid, b) {
+                Some((x, y)) => {
+                    next.push(x);
+                    next.push(y);
+                }
+                None => next.push(b.clone()),
+            }
+        }
+        basins = next;
+    }
+    let comps: Vec<&Vec<usize>> = basins
+        .iter()
+        .filter(|b| b.len() >= 4 && induced(&grid, b).edges.len() >= 3)
+        .collect();
+
+    // Per compartment: resistive axis (λ₂, mean R) + buffer axis (total & per-node).
+    println!(
+        "\n== Two resilience axes per compartment ({} compartments) ==",
+        comps.len()
+    );
+    println!(
+        "  {:>4} {:>6} {:>12} {:>12} {:>12} {:>12}",
+        "comp", "nodes", "λ₂ (steady)", "mean R", "buffer Σ", "buffer/node"
+    );
+    let (mut mean_r, mut mean_buf, mut lam2s): (Vec<f64>, Vec<f64>, Vec<f64>) =
+        (Vec::new(), Vec::new(), Vec::new());
+    for (ci, b) in comps.iter().enumerate() {
+        let sub = induced(&grid, b);
+        let cert = Resilience::from_eigenvalues(
+            &symmetric_eigen(&sub.laplacian_of(&vec![true; sub.edges.len()]), sub.n).values,
+            1e-9,
+        );
+        let h_local: Vec<f64> = b.iter().map(|&node| inertia[node]).collect();
+        let buf = compartment_buffer(&h_local, df_band);
+        let buf_per_node = buf / b.len() as f64;
+        println!(
+            "  {:>4} {:>6} {:>12.3e} {:>12.3e} {:>12.3e} {:>12.4}",
+            ci,
+            b.len(),
+            cert.lambda2,
+            cert.mean_resistance(),
+            buf,
+            buf_per_node
+        );
+        mean_r.push(cert.mean_resistance());
+        mean_buf.push(buf_per_node);
+        lam2s.push(cert.lambda2);
+    }
+
+    // Deconfound: size-normalized buffer (= mean inertia, topology-free) vs mean
+    // resistance (topology). Report with Jirak |ρ|√n — at n compartments these are
+    // not significant, consistent with the structural independence.
+    if comps.len() >= 3 {
+        let nn = comps.len() as f64;
+        let rho_rb = spearman(&mean_r, &mean_buf);
+        let rho_lb = spearman(&lam2s, &mean_buf);
+        println!("\n== Deconfound: is the buffer independent of the resistive axis? ==");
+        println!(
+            "  Spearman(mean resistance, buffer/node) = {rho_rb:+.3}  (|ρ|√n = {:.2})",
+            rho_rb.abs() * nn.sqrt()
+        );
+        println!(
+            "  Spearman(λ₂,               buffer/node) = {rho_lb:+.3}  (|ρ|√n = {:.2})",
+            rho_lb.abs() * nn.sqrt()
+        );
+        println!(
+            "  → both below the ~2 Jirak floor ⇒ no significant coupling: buffer (storage,\n    \
+             set by inertia) is a SEPARATE axis from connectivity (λ₂/Kirchhoff). Contrast\n    \
+             the `1/λ₂`↔Kirchhoff confound, which was near +1.0 (definitional). buffer/node\n    \
+             is mean inertia — topology-free by construction; any residual ρ is n={} noise.",
+            comps.len()
+        );
+    }
+
+    // Kugelstoßpendel: a sudden impulse hits each node of a compartment; the
+    // compartment yields where its PER-NODE buffer (mean inertia headroom) is
+    // thinnest. Per-node buffer varies independently of λ₂, so the compartment
+    // that yields need not be the resistively-weakest — that is the point.
+    let impulse = 0.032; // per-node strike, fraction of nominal (illustrative)
+    println!(
+        "\n== Kugelstoßpendel: per-node impulse {impulse} vs each compartment's buffer/node =="
+    );
+    // (compartment, per-node buffer, that compartment's λ₂), sorted thinnest first.
+    let mut weakest: Vec<(usize, f64, f64)> = (0..comps.len())
+        .map(|ci| (ci, mean_buf[ci], lam2s[ci]))
+        .collect();
+    weakest.sort_by(|a, b| a.1.partial_cmp(&b.1).unwrap());
+    let mut first_yield = None;
+    for (ci, buf, l2) in &weakest {
+        let y = ketchup_yield(impulse, *buf);
+        if y.yielded && first_yield.is_none() {
+            first_yield = Some((*ci, *l2));
+        }
+        println!(
+            "  comp {ci:>2}: buffer/node={buf:.4}  λ₂={l2:.2e}  headroom={:+.0}%  {}",
+            100.0 * y.headroom,
+            if y.yielded {
+                "YIELDS → Ketchup seed"
+            } else {
+                "holds (elastic)"
+            }
+        );
+    }
+    match first_yield {
+        Some((ci, l2)) => {
+            let resistively_healthy = lam2s.iter().filter(|&&x| x < l2).count();
+            println!(
+                "\n  → compartment {ci} yields FIRST (thinnest local buffer) at λ₂={l2:.2e},\n    \
+                 which is MORE connected than {resistively_healthy}/{} compartments — a\n    \
+                 resistance-only screen ranks it safe, the buffer axis flags it. That is the\n    \
+                 low-inertia failure mode (renewable-rich node, 28 Apr 2025).",
+                comps.len()
+            );
+        }
+        None => println!("\n  → all weakest nodes hold this impulse elastically (raise the strike to find the yield)."),
+    }
+    println!(
+        "\n  buffer = transient storage (inertia·band/f₀, the swing buffer); the yield is\n  \
+         the sharp non-Newtonian threshold. Inertia field ILLUSTRATIVE — feed real per-bus\n  \
+         H to calibrate; impulse magnitude should come from the flow redistribution (LODF)."
+    );
+}