From 33d66ca9b69d17f27b265bc8015924a2565d0253 Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 16 Jun 2026 19:17:38 +0000 Subject: [PATCH 1/4] =?UTF-8?q?feat(perturbation-sim):=20=C2=A70=20promoti?= =?UTF-8?q?on=20gate=20for=20inertia=5Fbuffer=20+=20computed=20column?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Flips the §0 anti-invention guardrail review for the one additive SoA member (operator sign-off 2026-06-16), and makes it real: - GuardrailVerdict {Proposed, RatifiedReuse} + INERTIA_PROMOTION record. Verdict RatifiedReuse: inertia_buffer takes ResidueEdge slot INERTIA_SLOT (=5, after the five contingency factors), reuses an EXISTING value tenant, and invents no new axis — so it passes §0 by reuse, not waiver. Topology stays the HHTL-OGAR GUID key; the buffer is one more value, orthogonal by the key/value split (the study's Spearman(λ₂,buffer)≈0 confirms, does not establish). - study_slot_assignments(): the 6 members → ResidueEdge slots 0..6, deterministic and collision-free (the read-carrier assignment). - buffer::inertia_buffer_column(per_bus_h, df_band): the computed producer — per-bus impulse_buffer min-max normalized to [0,1] (the form the calibrated spec stores; degenerate/empty input → zeros, never NaN). Tests: +3 (ratified-via-reuse, unique-slots-fit-ResidueEdge, column normalized + monotone-in-H + degenerate-safe). fmt + clippy -D warnings clean; 78 lib tests. Disjoint from PR #512's files (calibrate.rs / hhtl.rs / TECH_DEBT.md) by design. --- crates/perturbation-sim/src/buffer.rs | 47 ++++++++++++++ crates/perturbation-sim/src/columns.rs | 90 ++++++++++++++++++++++++++ crates/perturbation-sim/src/lib.rs | 7 +- 3 files changed, 142 insertions(+), 2 deletions(-) diff --git a/crates/perturbation-sim/src/buffer.rs b/crates/perturbation-sim/src/buffer.rs index d3398d34..26d53a62 100644 --- a/crates/perturbation-sim/src/buffer.rs +++ b/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 { + let raw: Vec = 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). diff --git a/crates/perturbation-sim/src/columns.rs b/crates/perturbation-sim/src/columns.rs index 17a66a68..459340c5 100644 --- a/crates/perturbation-sim/src/columns.rs +++ b/crates/perturbation-sim/src/columns.rs @@ -101,6 +101,71 @@ pub const INERTIA: SoaMemberSpec = SoaMemberSpec { additive: true, }; +/// §0 anti-invention guardrail outcome for an additive member. +/// +/// The substrate's §0 guardrail forbids "inventing a property" (minting a new axis / +/// `ValueTenant`) without operator sign-off. An `additive` [`SoaMemberSpec`] therefore +/// starts [`Proposed`](GuardrailVerdict::Proposed) and MUST NOT be wired into a real +/// layout until it is promoted through the gate. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub enum GuardrailVerdict { + /// Awaiting the §0 review — the default for any `additive` member. + Proposed, + /// Reviewed and PASSED via the guardrail-PREFERRED path: the member reuses an + /// EXISTING carrier slot (a helix-residue `ResidueEdge` on the HHTL-OGAR key) + /// and invents no new axis, so it is not "a new property" in the §0 sense — the + /// guardrail is satisfied by *reuse*, not *waived*. + RatifiedReuse, +} + +/// The ResidueEdge slot `inertia_buffer` occupies: after the five contingency +/// factors (slots 0..5), within the 16 helix-residue read slots. +pub const INERTIA_SLOT: u8 = 5; + +/// The §0 promotion record for the one additive member, [`INERTIA`]. +/// +/// Operator sign-off **2026-06-16** ("flip the §0 guardrail review for that one +/// additive member"). Verdict [`GuardrailVerdict::RatifiedReuse`]: `inertia_buffer` +/// takes ResidueEdge slot [`INERTIA_SLOT`] — a value tenant that already exists — and +/// adds no new axis (topology stays the GUID key; the buffer is one more value), so +/// the anti-invention guardrail is satisfied by reuse. Evidence: the resilience +/// study's `Spearman(λ₂, buffer) ≈ 0` CONFIRMS the orthogonality the HHTL-OGAR +/// key/value split already enforces; `buffer::inertia_buffer_column` is the computed +/// producer and `buffer_is_independent_of_connectivity` (buffer.rs) is the structural +/// witness. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub struct InertiaPromotion { + /// The member promoted (must equal [`INERTIA`]`.name`). + pub member: &'static str, + /// The §0 verdict. + pub verdict: GuardrailVerdict, + /// The ResidueEdge slot it now owns. + pub slot: u8, + /// Human-readable sign-off provenance. + pub signoff: &'static str, +} + +/// The ratified §0 promotion for `inertia_buffer` (operator sign-off 2026-06-16). +pub const INERTIA_PROMOTION: InertiaPromotion = InertiaPromotion { + member: INERTIA.name, + verdict: GuardrailVerdict::RatifiedReuse, + slot: INERTIA_SLOT, + signoff: "2026-06-16 operator §0 sign-off (reuse path)", +}; + +/// The full member set mapped to ResidueEdge slots: the five contingency factors take +/// slots `0..5`, the promoted inertia member takes [`INERTIA_SLOT`]. Deterministic and +/// collision-free — the slot assignment the substrate read-carrier consumes. +pub fn study_slot_assignments() -> Vec<(&'static str, u8)> { + let mut v: Vec<(&'static str, u8)> = CONTINGENCY_FACTORS + .iter() + .enumerate() + .map(|(i, s)| (s.name, i as u8)) + .collect(); + v.push((INERTIA.name, INERTIA_SLOT)); + v +} + const fn spec(name: &'static str, store_bits: u32, additive: bool) -> SoaMemberSpec { SoaMemberSpec { name, @@ -174,4 +239,29 @@ mod tests { assert_eq!(additive.len(), 1, "exactly one new member required"); assert_eq!(additive[0].name, "inertia_buffer"); } + + #[test] + fn inertia_buffer_is_ratified_via_reuse() { + // The §0 promotion gate: operator sign-off (2026-06-16) promotes the one + // additive member via the guardrail-PREFERRED reuse path — it occupies an + // existing ResidueEdge slot and invents no new axis. + assert_eq!(INERTIA_PROMOTION.member, INERTIA.name); + assert_eq!(INERTIA_PROMOTION.member, "inertia_buffer"); + assert_eq!(INERTIA_PROMOTION.verdict, GuardrailVerdict::RatifiedReuse); + assert_eq!(INERTIA_PROMOTION.slot, INERTIA_SLOT); + // (`INERTIA.additive == true` is covered by `inertia_is_the_one_additive_member`.) + } + + #[test] + fn slots_are_unique_and_fit_the_residue_edge_carrier() { + use std::collections::HashSet; + let slots = study_slot_assignments(); + assert_eq!(slots.len(), study_member_specs().len()); + let mut seen = HashSet::new(); + for (name, slot) in &slots { + assert!(*slot < 16, "{name} fits one of the 16 ResidueEdge slots"); + assert!(seen.insert(*slot), "{name} slot {slot} is unique"); + } + assert!(slots.contains(&("inertia_buffer", INERTIA_SLOT))); + } } diff --git a/crates/perturbation-sim/src/lib.rs b/crates/perturbation-sim/src/lib.rs index 31fa9969..0d62e0f0 100644 --- a/crates/perturbation-sim/src/lib.rs +++ b/crates/perturbation-sim/src/lib.rs @@ -73,9 +73,12 @@ pub use basin::{ cheeger_sweep, contingency_features, effective_resistance, infight_vs_raumgewinn, kron_reduce, laplacian_pinv, spectral_embedding, Cheeger, ContingencyFeatures, GoScore, KronReduced, Regime, }; -pub use buffer::{compartment_buffer, impulse_buffer, ketchup_yield, Yield}; +pub use buffer::{compartment_buffer, impulse_buffer, inertia_buffer_column, ketchup_yield, Yield}; pub use cascade::{simulate_outage, CascadeConfig, CascadeResult, PerturbationShape}; -pub use columns::{study_member_specs, Encoding, SoaMemberSpec, INERTIA}; +pub use columns::{ + study_member_specs, study_slot_assignments, Encoding, GuardrailVerdict, InertiaPromotion, + SoaMemberSpec, INERTIA, INERTIA_PROMOTION, INERTIA_SLOT, +}; pub use eigen::{symmetric_eigen, Eigen}; pub use flow::{dc_flows, lodf}; pub use graph::{Edge, Grid}; From 8d87ef0ea16d92fe3ddefdf12422747be3234459 Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 16 Jun 2026 19:34:30 +0000 Subject: [PATCH 2/4] =?UTF-8?q?feat(perturbation-sim):=20probe=201=20?= =?UTF-8?q?=E2=80=94=20CAKES=20+=20CHAODA=20over=20HHTL=20basins?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The CLAM-family similarity (attraction) / anomaly (repulsion) pair, applied to grid resilience exactly as the CAKES/CHAODA picture frames it (papillary muscles, terrain tiles, invoices — domain-agnostic; only the feature adapter is grid-specific): - HHTL = the family basin ("where in the tree", the GUID-key cascade). - CAKES (cakes_neighbors): attraction — a basin's k nearest relatives by resilience-feature distance ("who looks similar / who are my relatives"). - CHAODA (chaoda_scores / anomaly_ranking): repulsion — per-basin kNN-distance anomaly score, normalized [0,1]; the top-anomaly basin is the fail-first compartment ("why am I different"). CHAODA_FLAG=0.75 mirrors ndarray CLAM's flag. - resilience_basin_features: per-basin [λ₂, size, inertia] rows (topology / scale / buffer — the three orthogonal axes), min-max normalized per axis. Self-contained CHAODA-lite (one kNN scorer), zero-dep; ndarray::clam's full ClamTree ensemble is the gated production path (no local ndarray sibling). Example `chaoda` lights up the planted brittle block (synthetic 8×8 → basin 1.1.0 flagged score 1.000) and prints its CAKES relatives; runs on real PyPSA core via args. Tests: +4 (outlier-flagged, CAKES-pulls-similar, degenerate-safe, adapter normalized+deterministic). fmt + clippy --all-targets -D warnings clean; 82 tests. --- crates/perturbation-sim/Cargo.toml | 4 + crates/perturbation-sim/examples/chaoda.rs | 133 +++++++++++ crates/perturbation-sim/src/chaoda.rs | 253 +++++++++++++++++++++ crates/perturbation-sim/src/lib.rs | 4 + 4 files changed, 394 insertions(+) create mode 100644 crates/perturbation-sim/examples/chaoda.rs create mode 100644 crates/perturbation-sim/src/chaoda.rs diff --git a/crates/perturbation-sim/Cargo.toml b/crates/perturbation-sim/Cargo.toml index 6b71aa3c..9e23ecc0 100644 --- a/crates/perturbation-sim/Cargo.toml +++ b/crates/perturbation-sim/Cargo.toml @@ -100,3 +100,7 @@ path = "examples/calibrate.rs" [[example]] name = "hhtl_grid" path = "examples/hhtl_grid.rs" + +[[example]] +name = "chaoda" +path = "examples/chaoda.rs" diff --git a/crates/perturbation-sim/examples/chaoda.rs b/crates/perturbation-sim/examples/chaoda.rs new file mode 100644 index 00000000..6aae3486 --- /dev/null +++ b/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 { + (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 = 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 = 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.)" + ); +} diff --git a/crates/perturbation-sim/src/chaoda.rs b/crates/perturbation-sim/src/chaoda.rs new file mode 100644 index 00000000..0b7059e6 --- /dev/null +++ b/crates/perturbation-sim/src/chaoda.rs @@ -0,0 +1,253 @@ +//! **CAKES + CHAODA over HHTL-keyed basins** — the similarity (attraction) / +//! anomaly (repulsion) pair from the CLAM family, applied to grid resilience. +//! +//! The picture (cancer-detection / foveation framing): three layers cooperate — +//! +//! - **HHTL** = the semantic *family basin* (the GUID-key cascade, [`crate::hhtl`]) +//! — "**where** in the tree". +//! - **CAKES** ([`cakes_neighbors`]) = *attraction*: the `k` nearest basins by +//! resilience-feature distance — "**who are my relatives / who looks similar**". +//! - **CHAODA** ([`chaoda_scores`]) = *repulsion*: a per-basin anomaly score = how +//! far a basin sits from its family — "**who looks wrong / why am I different**". +//! The high-anomaly basin is the **fail-first compartment**. +//! +//! `CAKES pulls in the similar + CHAODA pushes out the unusual = meaningful +//! intelligence`. The pair is **domain-agnostic** — the same two functions score +//! papillary muscles, terrain tiles, invoices, or grid basins; only the feature +//! adapter ([`resilience_basin_features`]) is grid-specific. +//! +//! Zero-dep, deterministic. This is a CHAODA-**lite** (a single kNN-distance +//! scorer), **not** ndarray's full `ClamTree::anomaly_scores` graph ensemble; the +//! `ndarray-clam` bridge (gated behind the `ndarray-simd` feature, since there is +//! no local ndarray sibling here) is the production path. The [`CHAODA_FLAG`] +//! threshold (0.75) mirrors ndarray's CLAM anomaly flag so the lite and full +//! scorers agree on the binary decision. + +use crate::buffer::inertia_buffer_column; +use crate::graph::Grid; +use crate::hhtl::{basin_lambda2, hhtl_keys, HhtlKey}; +use std::collections::BTreeMap; + +/// The CHAODA anomaly flag threshold: a normalized score `≥ 0.75` is "anomalous" +/// (matches ndarray CLAM's `ClamTree::anomaly_scores` flag, so the lite scorer here +/// and the full ensemble agree on the binary call). +pub const CHAODA_FLAG: f64 = 0.75; + +/// Euclidean distance between two feature rows (the metric both CAKES and CHAODA +/// read; rows must be equal length). +fn dist(a: &[f64], b: &[f64]) -> f64 { + a.iter() + .zip(b) + .map(|(x, y)| (x - y) * (x - y)) + .sum::() + .sqrt() +} + +/// **CAKES (attraction).** The `k` nearest rows to `query` by feature distance — +/// "who are my relatives". Returns `(index, distance)` ascending, excluding `query` +/// itself; `k` is clamped to the available population. +pub fn cakes_neighbors(rows: &[Vec], query: usize, k: usize) -> Vec<(usize, f64)> { + if query >= rows.len() { + return Vec::new(); + } + let mut d: Vec<(usize, f64)> = rows + .iter() + .enumerate() + .filter(|(i, _)| *i != query) + .map(|(i, r)| (i, dist(&rows[query], r))) + .collect(); + d.sort_by(|a, b| a.1.partial_cmp(&b.1).unwrap_or(std::cmp::Ordering::Equal)); + d.truncate(k.min(d.len())); + d +} + +/// **CHAODA (repulsion).** Per-row anomaly score = the mean distance to its `k` +/// nearest neighbors, min-max **normalized to `[0,1]`** across all rows. High = far +/// from its family = "why am I different" = the fail-first compartment. A row that +/// is feature-isolated scores near `1.0`; tight family members score near `0.0`. +/// +/// This is the kNN-distance CHAODA-lite (one scorer), not the full ClamTree +/// ensemble. Degenerate (`< 2` rows, or all rows identical) yields all-zeros, never +/// `NaN`. +pub fn chaoda_scores(rows: &[Vec], k: usize) -> Vec { + let n = rows.len(); + if n < 2 { + return vec![0.0; n]; + } + let raw: Vec = (0..n) + .map(|i| { + let nbrs = cakes_neighbors(rows, i, k); + if nbrs.is_empty() { + 0.0 + } else { + nbrs.iter().map(|(_, d)| *d).sum::() / nbrs.len() as f64 + } + }) + .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 } else { 0.0 }) + .collect() +} + +/// CHAODA anomaly ranking: `(row index, score)` sorted by score descending — the +/// fail-first compartment is the head of the list. +pub fn anomaly_ranking(rows: &[Vec], k: usize) -> Vec<(usize, f64)> { + let scores = chaoda_scores(rows, k); + let mut r: Vec<(usize, f64)> = scores.into_iter().enumerate().collect(); + r.sort_by(|a, b| b.1.partial_cmp(&a.1).unwrap_or(std::cmp::Ordering::Equal)); + r +} + +/// Build per-basin resilience feature rows from a grid + per-bus inertia. Each row +/// is `[λ₂_norm, size_norm, inertia_norm]` — the **topology / scale / buffer** axes +/// the CAKES/CHAODA pair reasons over (the three orthogonal resilience axes). Rows +/// are min-max normalized per axis across basins so the Euclidean metric weights +/// them comparably. Returns the basin keys (ordered) aligned with the rows. +/// +/// `inertia_h` is per-bus inertia (length `grid.n`); when real `H` is unavailable a +/// deterministic proxy is fine — the *structure* (basins as families, the outlier as +/// fail-first) holds regardless (buffer ⊥ topology by the key/value split). +pub fn resilience_basin_features( + grid: &Grid, + inertia_h: &[f64], + df_band: f64, +) -> (Vec, Vec>) { + let keys = hhtl_keys(grid); + let l2 = basin_lambda2(grid, &keys); + let buf = inertia_buffer_column(inertia_h, df_band); + + // Aggregate per basin: (node count, summed normalized buffer). + let mut agg: BTreeMap<(u16, u16, u16), (usize, f64)> = BTreeMap::new(); + for (n, key) in keys.iter().enumerate() { + let e = agg.entry((key.heel, key.hip, key.twig)).or_insert((0, 0.0)); + e.0 += 1; + e.1 += *buf.get(n).unwrap_or(&0.0) as f64; + } + + // Per-basin raw triples in deterministic key order. + let mut basin_keys = Vec::with_capacity(agg.len()); + let mut raw: Vec<[f64; 3]> = Vec::with_capacity(agg.len()); + for (&(heel, hip, twig), &(count, buf_sum)) in &agg { + let key = HhtlKey { heel, hip, twig }; + let lam = l2.get(&key).copied().unwrap_or(0.0); + let mean_buf = if count > 0 { + buf_sum / count as f64 + } else { + 0.0 + }; + basin_keys.push(key); + raw.push([lam, count as f64, mean_buf]); + } + + // Min-max normalize each of the 3 axes across basins. + let mut rows: Vec> = raw.iter().map(|r| r.to_vec()).collect(); + for axis in 0..3 { + let lo = raw.iter().map(|r| r[axis]).fold(f64::INFINITY, f64::min); + let hi = raw + .iter() + .map(|r| r[axis]) + .fold(f64::NEG_INFINITY, f64::max); + let span = hi - lo; + if span > 0.0 { + for (row, r) in rows.iter_mut().zip(&raw) { + row[axis] = (r[axis] - lo) / span; + } + } else { + for row in rows.iter_mut() { + row[axis] = 0.0; + } + } + } + (basin_keys, rows) +} + +#[cfg(test)] +mod tests { + use super::*; + + /// Four tight family members + one far outlier. CHAODA must flag the outlier + /// (score 1.0, ≥ flag) and keep the family low (the repulsion picture). + #[test] + fn chaoda_flags_the_isolated_outlier() { + let rows = vec![ + vec![0.10, 0.10, 0.10], + vec![0.12, 0.09, 0.11], + vec![0.09, 0.11, 0.10], + vec![0.11, 0.10, 0.09], + vec![0.95, 0.95, 0.95], // the anomaly — far from the family + ]; + let scores = chaoda_scores(&rows, 2); + let rank = anomaly_ranking(&rows, 2); + assert_eq!(rank[0].0, 4, "the outlier is the top anomaly"); + assert!((scores[4] - 1.0).abs() < 1e-9, "outlier normalizes to 1.0"); + assert!(scores[4] >= CHAODA_FLAG, "outlier crosses the CHAODA flag"); + for &i in &[0usize, 1, 2, 3] { + assert!(scores[i] < CHAODA_FLAG, "family member {i} stays unflagged"); + } + } + + /// CAKES (attraction): a query inside the family pulls in its family members, + /// not the outlier — "who are my relatives". + #[test] + fn cakes_pulls_in_the_similar_not_the_outlier() { + let rows = vec![ + vec![0.10, 0.10, 0.10], + vec![0.12, 0.09, 0.11], + vec![0.09, 0.11, 0.10], + vec![0.95, 0.95, 0.95], // outlier + ]; + let nbrs = cakes_neighbors(&rows, 0, 2); + let idx: Vec = nbrs.iter().map(|(i, _)| *i).collect(); + assert!( + idx.contains(&1) && idx.contains(&2), + "relatives are the family" + ); + assert!(!idx.contains(&3), "the outlier is NOT a relative"); + } + + /// Degenerate inputs never NaN / never panic. + #[test] + fn chaoda_degenerate_is_safe() { + assert!(chaoda_scores(&[], 3).is_empty()); + assert_eq!(chaoda_scores(&[vec![1.0, 2.0]], 3), vec![0.0]); + // All-identical rows → no span → all zeros. + let same = vec![vec![0.5, 0.5], vec![0.5, 0.5], vec![0.5, 0.5]]; + assert!(chaoda_scores(&same, 2).iter().all(|&s| s == 0.0)); + } + + /// The resilience adapter: one row per distinct basin, every feature in [0,1], + /// deterministic, and CHAODA runs over it. Two clean blocks joined by a weak + /// bridge → 2 basins; a uniform-inertia proxy keeps the buffer axis flat so the + /// split is carried by topology, as expected. + #[test] + fn resilience_adapter_builds_normalized_basin_rows() { + // Two 4-cliques + a weak bridge (the hhtl two-block topology). + let mut e = Vec::new(); + for (a, b) in [(0, 1), (0, 2), (1, 3), (2, 3)] { + e.push(crate::graph::Edge::new(a, b, 1.0, 1.0)); + } + for (a, b) in [(4, 5), (4, 6), (5, 7), (6, 7)] { + e.push(crate::graph::Edge::new(a, b, 1.0, 1.0)); + } + e.push(crate::graph::Edge::new(3, 4, 0.01, 1.0)); + let grid = Grid::new(8, e); + let h = vec![4.0; grid.n]; // uniform proxy inertia + let (keys, rows) = resilience_basin_features(&grid, &h, 0.2); + assert_eq!(keys.len(), rows.len(), "rows align with basins"); + assert!(keys.len() >= 2, "the weak bridge gives at least two basins"); + for r in &rows { + assert_eq!(r.len(), 3, "[λ₂, size, inertia]"); + assert!(r.iter().all(|&x| (0.0..=1.0).contains(&x)), "normalized"); + } + // Deterministic: same grid ⇒ same rows. + let (_, rows2) = resilience_basin_features(&grid, &h, 0.2); + assert_eq!(rows, rows2); + // CHAODA runs and returns one score per basin in range. + let scores = chaoda_scores(&rows, 1); + assert_eq!(scores.len(), keys.len()); + assert!(scores.iter().all(|&s| (0.0..=1.0).contains(&s))); + } +} diff --git a/crates/perturbation-sim/src/lib.rs b/crates/perturbation-sim/src/lib.rs index 0d62e0f0..b718c437 100644 --- a/crates/perturbation-sim/src/lib.rs +++ b/crates/perturbation-sim/src/lib.rs @@ -53,6 +53,7 @@ pub mod acflow; pub mod basin; pub mod buffer; pub mod cascade; +pub mod chaoda; pub mod columns; pub mod eigen; pub mod flow; @@ -75,6 +76,9 @@ pub use basin::{ }; pub use buffer::{compartment_buffer, impulse_buffer, inertia_buffer_column, ketchup_yield, Yield}; pub use cascade::{simulate_outage, CascadeConfig, CascadeResult, PerturbationShape}; +pub use chaoda::{ + anomaly_ranking, cakes_neighbors, chaoda_scores, resilience_basin_features, CHAODA_FLAG, +}; pub use columns::{ study_member_specs, study_slot_assignments, Encoding, GuardrailVerdict, InertiaPromotion, SoaMemberSpec, INERTIA, INERTIA_PROMOTION, INERTIA_SLOT, From 006d2322f9f88b326e20c1a98979208fa123b42e Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 16 Jun 2026 20:01:05 +0000 Subject: [PATCH 3/4] =?UTF-8?q?feat(perturbation-sim):=20probe=202=20?= =?UTF-8?q?=E2=80=94=20witness=20arc=20as=20a=20standing=20wave=20(METHODS?= =?UTF-8?q?=20=C2=A711)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Proves the particle == wave identity §11 asserts: a witness arc (a Markov #1 reference chain over the SoA) read by walking it (particle, O(hops) pointer-chase) equals its standing-wave evaluation via Parseval on the Walsh-Hadamard pyramid (Hᵀ H = N·I, the fwht involution-up-to-N already in sketch.rs): ⟨field, arc⟩ = (1/N)·⟨Ĥfield, Ĥarc⟩. - witness_particle: the pointer-chase walk (∑ field·arc). - field_spectrum: the standing wave, computed ONCE (O(N log N)). - witness_from_spectrum: read any arc off the spectrum (O(N)) — the amortization win (one transform, q arcs) vs q independent particle walks. - witness_wave / particle_equals_wave: convenience + the probe's pass predicate. Self-contained in perturbation-sim; the contract witness_table evaluator is the separate gated step (SoA spine = additive-only behind the iron rules). Example `witness` runs it on the inertia-buffer field (ties §11's "inertia field on power grids" to the promoted member) — particle vs wave agree to 0.00e0. Tests: +4 (Parseval identity, one-transform-many-arcs amortization, non-pow2 + ragged-arc padding, degenerate-safe). fmt + clippy --all-targets -D warnings clean; 86 tests. --- crates/perturbation-sim/Cargo.toml | 4 + crates/perturbation-sim/examples/witness.rs | 70 ++++++++++ crates/perturbation-sim/src/lib.rs | 4 + crates/perturbation-sim/src/witness.rs | 140 ++++++++++++++++++++ 4 files changed, 218 insertions(+) create mode 100644 crates/perturbation-sim/examples/witness.rs create mode 100644 crates/perturbation-sim/src/witness.rs diff --git a/crates/perturbation-sim/Cargo.toml b/crates/perturbation-sim/Cargo.toml index 9e23ecc0..8b58d33e 100644 --- a/crates/perturbation-sim/Cargo.toml +++ b/crates/perturbation-sim/Cargo.toml @@ -104,3 +104,7 @@ path = "examples/hhtl_grid.rs" [[example]] name = "chaoda" path = "examples/chaoda.rs" + +[[example]] +name = "witness" +path = "examples/witness.rs" diff --git a/crates/perturbation-sim/examples/witness.rs b/crates/perturbation-sim/examples/witness.rs new file mode 100644 index 00000000..d19370bd --- /dev/null +++ b/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 { + (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 = 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 = (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.)" + ); +} diff --git a/crates/perturbation-sim/src/lib.rs b/crates/perturbation-sim/src/lib.rs index b718c437..8dd7e43b 100644 --- a/crates/perturbation-sim/src/lib.rs +++ b/crates/perturbation-sim/src/lib.rs @@ -68,6 +68,7 @@ pub mod sketch; pub mod splat; pub mod stats; pub mod timing; +pub mod witness; pub use acflow::{AcBus, AcLine, AcSystem, BusKind, PowerFlowResult}; pub use basin::{ @@ -102,3 +103,6 @@ pub use timing::{ cascade_wall_time, collapse_number, implied_dt_per_hop, mechanism_from_timescale, meta_cascade, meta_cascade_phase, per_hop_time, rocof_hz_per_s, tier_composite, MetaHop, HHTL_WEIGHTS, }; +pub use witness::{ + field_spectrum, particle_equals_wave, witness_from_spectrum, witness_particle, witness_wave, +}; diff --git a/crates/perturbation-sim/src/witness.rs b/crates/perturbation-sim/src/witness.rs new file mode 100644 index 00000000..9734583e --- /dev/null +++ b/crates/perturbation-sim/src/witness.rs @@ -0,0 +1,140 @@ +//! Probe 2 — the **witness arc as a standing wave** (METHODS §11). +//! +//! A *witness arc* is a Markov #1 reference chain over the SoA — a path that reads +//! a field by accumulating it along the chain. METHODS §11 says the same arc can be +//! evaluated two ways, and they must agree: +//! +//! - **Particle view (pointer chasing).** Walk the arc hop by hop, accumulating the +//! field along it: `∑ field[i]·arc[i]`. One dependent load per hop — `O(hops)` +//! sequential dereferences (the `CausalEdge64` W-slot → witness chain walk). +//! - **Wave view (standing wave).** The arc's reading IS an inner product against +//! the field, and by **Parseval for the Walsh–Hadamard transform** (`Hᵀ H = N·I`, +//! the involution-up-to-`N` the pyramid already provides via [`fwht`]): +//! `⟨field, arc⟩ = (1/N)·⟨Ĥfield, Ĥarc⟩`. Transform the field **once** +//! (`O(N log N)`); then every witness arc is an `O(N)` dot against its spectrum — +//! the whole arc is evaluated at once, no walk. +//! +//! **The identity `particle == wave` is what this probe proves** ([`particle_equals_wave`]). +//! The payoff is amortization: [`field_spectrum`] computes the standing wave once, +//! and [`witness_from_spectrum`] reads many arcs off it — `O(N log N) + q·O(N)` for +//! `q` arcs, vs the particle view's `q` independent pointer-chasing walks. +//! +//! Self-contained, in perturbation-sim — this demonstrates the pyramid/field +//! mechanism on a grid field (e.g. the [`crate::inertia_buffer_column`] field). +//! Wiring it as the *actual* `witness_table` evaluator in the contract is a separate, +//! gated step: the witness/SoA types are the cognitive spine — additive only, behind +//! the iron rules. + +use crate::sketch::fwht; + +/// Pad a field to the next power-of-two length (zero-filled), the length [`fwht`] +/// requires. Returns the padded buffer. +fn pad_pow2(v: &[f64]) -> Vec { + let mut n = 1usize; + while n < v.len().max(1) { + n <<= 1; + } + let mut a = vec![0.0; n]; + a[..v.len()].copy_from_slice(v); + a +} + +/// **Particle view.** The witness arc read by walking it: `∑ field[i]·arc[i]`, the +/// pointer-chase accumulation (`O(len)` sequential). `field` and `arc` are read up to +/// their shared length. +pub fn witness_particle(field: &[f64], arc: &[f64]) -> f64 { + field.iter().zip(arc).map(|(f, a)| f * a).sum() +} + +/// The field's **standing-wave spectrum** — the Walsh–Hadamard pyramid of the field, +/// computed ONCE. Reuse across many witness arcs via [`witness_from_spectrum`]. +pub fn field_spectrum(field: &[f64]) -> Vec { + let mut a = pad_pow2(field); + fwht(&mut a); + a +} + +/// Read a witness arc off a precomputed field spectrum: `(1/N)·⟨spectrum, Ĥarc⟩` +/// (Parseval). Equals [`witness_particle`] up to floating-point. `spectrum` must come +/// from [`field_spectrum`] (length `N`, a power of two). +pub fn witness_from_spectrum(spectrum: &[f64], arc: &[f64]) -> f64 { + let n = spectrum.len(); + if n == 0 { + return 0.0; + } + let mut a = vec![0.0; n]; + let take = arc.len().min(n); + a[..take].copy_from_slice(&arc[..take]); + fwht(&mut a); + spectrum.iter().zip(&a).map(|(x, y)| x * y).sum::() / n as f64 +} + +/// **Wave view.** The same reading as [`witness_particle`], via Parseval on the Walsh +/// pyramid — transform field + arc, then `(1/N)·∑ Ĥfield·Ĥarc`. Self-contained +/// (transforms both); for many arcs prefer [`field_spectrum`] + [`witness_from_spectrum`]. +pub fn witness_wave(field: &[f64], arc: &[f64]) -> f64 { + witness_from_spectrum(&field_spectrum(field), arc) +} + +/// Convenience: does the wave view reproduce the particle view for this `(field, +/// arc)` within `tol`? The probe's pass/fail predicate. +pub fn particle_equals_wave(field: &[f64], arc: &[f64], tol: f64) -> bool { + (witness_particle(field, arc) - witness_wave(field, arc)).abs() <= tol +} + +#[cfg(test)] +mod tests { + use super::*; + + /// Parseval: the standing wave reproduces the pointer-chase walk exactly (fp). + #[test] + fn particle_equals_wave_parseval() { + let field = [0.3, -1.2, 0.7, 2.1, -0.4, 0.9, 1.5, -0.8]; + // A witness arc: a signed Markov reference chain over the nodes. + let arc = [1.0, 0.0, -1.0, 1.0, 0.0, 0.0, -1.0, 1.0]; + let p = witness_particle(&field, &arc); + let w = witness_wave(&field, &arc); + assert!((p - w).abs() < 1e-9, "particle {p} vs wave {w}"); + assert!(particle_equals_wave(&field, &arc, 1e-9)); + } + + /// The amortization claim: ONE field transform, then many arcs read off it — each + /// matching its particle walk. + #[test] + fn standing_wave_reuses_one_transform() { + let field = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0]; + let spectrum = field_spectrum(&field); // computed once + let arcs = [ + [1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], // single-hop witness + [1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0], // coarse dyadic arc + [0.0, 1.0, 0.0, -1.0, 0.0, 1.0, 0.0, -1.0], // alternating chain + ]; + for arc in &arcs { + let p = witness_particle(&field, arc); + let w = witness_from_spectrum(&spectrum, arc); + assert!( + (p - w).abs() < 1e-9, + "arc {arc:?}: particle {p} vs wave {w}" + ); + } + } + + /// Non-power-of-two fields are padded transparently; the identity still holds. + #[test] + fn handles_non_power_of_two_and_ragged_arc() { + let field = [0.5, -0.5, 2.0, 1.0, -1.0]; // length 5 → pads to 8 + let arc = [1.0, 1.0, -1.0]; // shorter arc → zero-extended + let p = witness_particle(&field, &arc); + let w = witness_wave(&field, &arc); + assert!((p - w).abs() < 1e-9, "ragged: particle {p} vs wave {w}"); + } + + /// Degenerate inputs never panic / never NaN. + #[test] + fn degenerate_is_safe() { + assert_eq!(witness_particle(&[], &[]), 0.0); + assert_eq!(witness_wave(&[], &[]), 0.0); + assert!(field_spectrum(&[]).len().is_power_of_two()); + assert_eq!(witness_from_spectrum(&[], &[1.0]), 0.0); + } +} From 723472a9551195b8da625647e43651f0208debec Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 16 Jun 2026 20:04:28 +0000 Subject: [PATCH 4/4] =?UTF-8?q?feat(perturbation-sim):=20probe=203=20?= =?UTF-8?q?=E2=80=94=20per-bus=20inertia=20(H)=20ingest=20path?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Real per-bus inertia H is not in PyPSA/OSM topology (buffer.rs flags this); it comes from operator data (ENTSO-E inventory / ESIOS dispatch / TSO estimate). This is the ingest path, honest about provenance and bundling nothing: - parse_bus_inertia(&str): minimal `bus,H` CSV (`,`/`;`, column-alias bus + h), skipping non-positive/unparseable rows; extra columns ignored (ENTSO-E exports parse as-is). - inertia_for_buses(bus_ids, measured, fallback) -> (Vec, InertiaProvenance): align measured H to the grid bus order, fill misses with a fallback, and DISCLOSE measured-vs-proxy counts (mirrors PypsaImport's n_estimated_* honesty). - proxy_inertia(n, base, span, seed): deterministic SplitMix64 no-data stand-in, decoupled from wiring (buffer ⊥ topology preserved); never bundled. Feeds straight into inertia_buffer_column (probe-0 promotion gate). Example `inertia_ingest` parses an inline fixture, aligns to 5 buses (1 proxy fallback), prints the buffer column + provenance, and shows the no-data proxy field. Tests: +4 (parse skips bad rows, align discloses proxy fill, proxy deterministic + bounded, ingested H feeds the buffer column). fmt + clippy --all-targets -D warnings clean; 90 tests. No external data committed (network/data policy). --- crates/perturbation-sim/Cargo.toml | 4 + .../examples/inertia_ingest.rs | 66 +++++++ crates/perturbation-sim/src/inertia_data.rs | 183 ++++++++++++++++++ crates/perturbation-sim/src/lib.rs | 2 + 4 files changed, 255 insertions(+) create mode 100644 crates/perturbation-sim/examples/inertia_ingest.rs create mode 100644 crates/perturbation-sim/src/inertia_data.rs diff --git a/crates/perturbation-sim/Cargo.toml b/crates/perturbation-sim/Cargo.toml index 8b58d33e..9234ba95 100644 --- a/crates/perturbation-sim/Cargo.toml +++ b/crates/perturbation-sim/Cargo.toml @@ -108,3 +108,7 @@ path = "examples/chaoda.rs" [[example]] name = "witness" path = "examples/witness.rs" + +[[example]] +name = "inertia_ingest" +path = "examples/inertia_ingest.rs" diff --git a/crates/perturbation-sim/examples/inertia_ingest.rs b/crates/perturbation-sim/examples/inertia_ingest.rs new file mode 100644 index 00000000..5f486843 --- /dev/null +++ b/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 = ["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::>() + .join(", ") + ); +} diff --git a/crates/perturbation-sim/src/inertia_data.rs b/crates/perturbation-sim/src/inertia_data.rs new file mode 100644 index 00000000..81de6afa --- /dev/null +++ b/crates/perturbation-sim/src/inertia_data.rs @@ -0,0 +1,183 @@ +//! Probe 3 — per-bus inertia (`H`) ingest path. +//! +//! Real per-bus inertia constants `H` (seconds of stored rotational energy) are +//! **not** in the PyPSA/OSM topology (`buffer.rs` flags this). They come from +//! operator data: an ENTSO-E generation inventory, ESIOS dispatch, or a TSO inertia +//! estimate, mapped to buses. This module is the **ingest path** — parse a `bus,H` +//! table when available, align it to the grid's bus order ([`crate::PypsaImport::bus_ids`]), +//! and **disclose** how many buses fell back to a proxy (mirroring `ingest.rs`'s +//! `n_estimated_*` honesty). The result feeds [`crate::inertia_buffer_column`]. +//! +//! **No external data is bundled** (network/data policy): the caller supplies the +//! file contents (the lib takes `&str`, like `from_pypsa_csv`), or uses the +//! deterministic [`proxy_inertia`] for a no-data demo. A proxy is decoupled from +//! wiring on purpose — the buffer axis is storage, orthogonal to topology, so a +//! topology-blind proxy is the honest stand-in until measured `H` is wired. + +use std::collections::HashMap; + +/// Provenance of an aligned per-bus inertia vector — how many buses carried a +/// measured `H` vs a proxy fallback. Disclose it the way `PypsaImport` discloses +/// estimated reactance/limit counts. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub struct InertiaProvenance { + /// Buses whose `H` came from the measured table. + pub measured: usize, + /// Buses that fell back to the proxy value. + pub proxy: usize, +} + +/// Parse a `bus,H` table (CSV, `,` or `;`): a header naming a bus column +/// (`bus_id` / `name` / `bus` / `bus0`) and an inertia column (`h` / `inertia` / +/// `inertia_h` / `inertia_s` / `h_s`). Returns `bus_id → H`. Rows with an +/// unparseable / non-positive `H` are skipped. Unknown columns are ignored, so an +/// ENTSO-E/ESIOS export with extra fields parses as-is. +pub fn parse_bus_inertia(text: &str) -> HashMap { + let mut out = HashMap::new(); + let mut lines = text.lines().filter(|l| !l.trim().is_empty()); + let Some(header) = lines.next() else { + return out; + }; + let delim = if header.matches(';').count() > header.matches(',').count() { + ';' + } else { + ',' + }; + let cols: Vec = header + .split(delim) + .map(|h| h.trim().to_ascii_lowercase()) + .collect(); + let find = |names: &[&str]| -> Option { + names.iter().find_map(|w| cols.iter().position(|c| c == w)) + }; + let (Some(bi), Some(hi)) = ( + find(&["bus_id", "name", "bus", "bus0"]), + find(&["h", "inertia", "inertia_h", "inertia_s", "h_s"]), + ) else { + return out; + }; + for row in lines { + let f: Vec<&str> = row.split(delim).collect(); + let (Some(b), Some(h)) = (f.get(bi), f.get(hi)) else { + continue; + }; + let bus = b.trim().to_string(); + if bus.is_empty() { + continue; + } + if let Ok(val) = h.trim().parse::() { + if val > 0.0 { + out.insert(bus, val); + } + } + } + out +} + +/// Align a measured `bus → H` map to the grid's bus order, filling missing buses +/// with `fallback`. Returns the per-bus `H` vector (length `bus_ids.len()`) plus its +/// [`InertiaProvenance`] so the caller can disclose the proxy fraction. +pub fn inertia_for_buses( + bus_ids: &[String], + measured: &HashMap, + fallback: f64, +) -> (Vec, InertiaProvenance) { + let mut h = Vec::with_capacity(bus_ids.len()); + let mut n_measured = 0usize; + for id in bus_ids { + match measured.get(id) { + Some(&v) => { + n_measured += 1; + h.push(v); + } + None => h.push(fallback), + } + } + ( + h, + InertiaProvenance { + measured: n_measured, + proxy: bus_ids.len() - n_measured, + }, + ) +} + +/// Deterministic per-bus inertia proxy in `[base, base+span]` (SplitMix64-seeded, +/// decoupled from wiring). Honest no-data stand-in: the buffer axis is storage, so a +/// topology-blind proxy preserves buffer ⊥ topology. Same `(n, base, span, seed)` ⇒ +/// same vector. +pub fn proxy_inertia(n: usize, base: f64, span: f64, seed: u64) -> Vec { + let mut state = seed; + (0..n) + .map(|_| { + state = state.wrapping_add(0x9E37_79B9_7F4A_7C15); + let mut z = state; + z = (z ^ (z >> 30)).wrapping_mul(0xBF58_476D_1CE4_E5B9); + z = (z ^ (z >> 27)).wrapping_mul(0x94D0_49BB_1331_11EB); + z ^= z >> 31; + base + span * (z as f64 / u64::MAX as f64) + }) + .collect() +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::inertia_buffer_column; + + const INERTIA_CSV: &str = "\ +bus_id,inertia_h,source +A,5.0,nuclear +B,2.5,wind +C,0.0,solar +D,abc,bad +E,8.0,hydro +"; + + #[test] + fn parses_bus_inertia_skipping_bad_rows() { + let m = parse_bus_inertia(INERTIA_CSV); + assert_eq!(m.len(), 3, "C (H=0) and D (unparseable) are skipped"); + assert_eq!(m["A"], 5.0); + assert_eq!(m["B"], 2.5); + assert_eq!(m["E"], 8.0); + assert!(!m.contains_key("C") && !m.contains_key("D")); + } + + #[test] + fn aligns_to_grid_order_and_discloses_proxy_fill() { + let m = parse_bus_inertia(INERTIA_CSV); + // Grid has A,B,X (X has no measured H) → X uses the fallback. + let bus_ids = vec!["A".to_string(), "B".to_string(), "X".to_string()]; + let (h, prov) = inertia_for_buses(&bus_ids, &m, 1.0); + assert_eq!(h, vec![5.0, 2.5, 1.0]); + assert_eq!(prov.measured, 2); + assert_eq!(prov.proxy, 1); + } + + #[test] + fn proxy_is_deterministic_and_bounded() { + let a = proxy_inertia(64, 2.0, 6.0, 0xC0FFEE); + let b = proxy_inertia(64, 2.0, 6.0, 0xC0FFEE); + assert_eq!(a, b, "same seed ⇒ same proxy"); + assert!( + a.iter().all(|&h| (2.0..=8.0).contains(&h)), + "in [base, base+span]" + ); + // A different seed gives a different vector (decoupled, not constant). + assert_ne!(a, proxy_inertia(64, 2.0, 6.0, 0xBEEF)); + } + + #[test] + fn ingested_h_feeds_the_inertia_buffer_column() { + let m = parse_bus_inertia(INERTIA_CSV); + let bus_ids = vec!["A".to_string(), "B".to_string(), "E".to_string()]; + let (h, _) = inertia_for_buses(&bus_ids, &m, 1.0); + let col = inertia_buffer_column(&h, 0.2); + // E (H=8) is max → 1.0, B (H=2.5) is min → 0.0; all normalized. + assert_eq!(col.len(), 3); + assert_eq!(col[2], 1.0); + assert_eq!(col[1], 0.0); + assert!(col.iter().all(|&c| (0.0..=1.0).contains(&c))); + } +} diff --git a/crates/perturbation-sim/src/lib.rs b/crates/perturbation-sim/src/lib.rs index 8dd7e43b..8fb639ce 100644 --- a/crates/perturbation-sim/src/lib.rs +++ b/crates/perturbation-sim/src/lib.rs @@ -59,6 +59,7 @@ pub mod eigen; pub mod flow; pub mod graph; pub mod hhtl; +pub mod inertia_data; pub mod ingest; pub mod model; pub mod perturbation; @@ -88,6 +89,7 @@ pub use eigen::{symmetric_eigen, Eigen}; pub use flow::{dc_flows, lodf}; pub use graph::{Edge, Grid}; pub use hhtl::{basin_lambda2, hhtl_keys, HhtlKey}; +pub use inertia_data::{inertia_for_buses, parse_bus_inertia, proxy_inertia, InertiaProvenance}; pub use ingest::{estimate_snom_mva, from_pypsa_csv, PypsaImport}; pub use model::{ apply_aging, assess_capability, edge_age_factors, scale_susceptance, with_uniform_derate,