|
| 1 | +//! PROBE — CHAODA epicenter (`ndarray::hpc::clam`) on the real ES grid. |
| 2 | +//! **Gated behind `ndarray-simd`** (uses ndarray's production `ClamTree` + |
| 3 | +//! `anomaly_scores`, the real CLAM/CHAODA engine — not perturbation-sim's lite). |
| 4 | +//! |
| 5 | +//! The WHERE axis of the three-axis surge decomposition |
| 6 | +//! (`E-FAMILY-BASIN-WEYL-HOP-LOCAL-AT-CRISP-TIER`): a surge's fail-first |
| 7 | +//! compartment is the brittle seam — geometrically anomalous on the node |
| 8 | +//! manifold. CHAODA's LFD-based anomaly should flag it **statically, with no |
| 9 | +//! cascade simulation**. Hypothesis: the top-anomaly nodes concentrate on the |
| 10 | +//! HHTL seam — the inter-HEEL-basin cut, the crisp 2-split bottleneck (the |
| 11 | +//! "without family nodes" table's HIP 2-line seam, lines 46/150). |
| 12 | +//! |
| 13 | +//! Encoding: each node → a binary fingerprint = the **sign of its top-K Fiedler |
| 14 | +//! eigenvector coordinates**, packed to bytes (the spectral embedding, the same |
| 15 | +//! Laplacian spectrum the HHTL tiers come from). `ClamTree::build` (Hamming) + |
| 16 | +//! `anomaly_scores` → per-node LFD anomaly. We then compare the seam-cut |
| 17 | +//! endpoints' anomaly against the global mean and their share of the top anomalies. |
| 18 | +//! |
| 19 | +//! Honest gate: the epicenter↔anomaly claim is SUPPORTED only if the seam nodes' |
| 20 | +//! mean anomaly ≥ 1.30× the global mean AND seam nodes are over-represented in |
| 21 | +//! the top-quartile anomalies (lift ≥ 1.30). Reported as-is otherwise. |
| 22 | +//! |
| 23 | +//! Run: cargo run --manifest-path crates/perturbation-sim/Cargo.toml \ |
| 24 | +//! --features ndarray-simd --example chaoda_surge_epicenter \ |
| 25 | +//! -- /tmp/pypsa/buses.csv /tmp/pypsa/lines.csv ES |
| 26 | +
|
| 27 | +#[cfg(not(feature = "ndarray-simd"))] |
| 28 | +fn main() { |
| 29 | + eprintln!( |
| 30 | + "chaoda_surge_epicenter requires `--features ndarray-simd` \ |
| 31 | + (it uses ndarray::hpc::clam — the real CLAM/CHAODA engine)." |
| 32 | + ); |
| 33 | +} |
| 34 | + |
| 35 | +#[cfg(feature = "ndarray-simd")] |
| 36 | +fn main() { |
| 37 | + use ndarray::hpc::clam::ClamTree; |
| 38 | + use perturbation_sim::{from_pypsa_csv, hhtl_keys, symmetric_eigen, Grid}; |
| 39 | + use std::collections::HashSet; |
| 40 | + use std::fs; |
| 41 | + |
| 42 | + const K: usize = 64; // spectral-embedding bits per node (Fiedler coords 1..=K) |
| 43 | + |
| 44 | + let args: Vec<String> = std::env::args().collect(); |
| 45 | + let (bpath, lpath, country) = ( |
| 46 | + args.get(1) |
| 47 | + .map(String::as_str) |
| 48 | + .unwrap_or("/tmp/pypsa/buses.csv"), |
| 49 | + args.get(2) |
| 50 | + .map(String::as_str) |
| 51 | + .unwrap_or("/tmp/pypsa/lines.csv"), |
| 52 | + args.get(3).map(String::as_str).unwrap_or("ES"), |
| 53 | + ); |
| 54 | + let buses = fs::read_to_string(bpath).expect("read buses.csv"); |
| 55 | + let lines = fs::read_to_string(lpath).expect("read lines.csv"); |
| 56 | + let import = from_pypsa_csv(&buses, &lines, Some(country)) |
| 57 | + .expect("parse pypsa") |
| 58 | + .largest_component(); |
| 59 | + let grid: &Grid = &import.grid; |
| 60 | + let (n, m) = (grid.n, grid.edges.len()); |
| 61 | + |
| 62 | + // Spectral embedding: top-K Fiedler eigenvectors (skip j=0, the constant). |
| 63 | + let eig = symmetric_eigen(&grid.laplacian_of(&vec![true; m]), n); |
| 64 | + let k = K.min(n.saturating_sub(1)); |
| 65 | + let vecs: Vec<Vec<f64>> = (1..=k).map(|j| eig.eigenvector(j)).collect(); |
| 66 | + |
| 67 | + // Node fingerprint = sign bits of its K Fiedler coords, packed to bytes. |
| 68 | + let vec_len = k.div_ceil(8); |
| 69 | + let mut data = vec![0u8; n * vec_len]; |
| 70 | + for (node, fp) in data.chunks_mut(vec_len).enumerate() { |
| 71 | + for (j, v) in vecs.iter().enumerate() { |
| 72 | + if v[node] >= 0.0 { |
| 73 | + fp[j / 8] |= 1 << (j % 8); |
| 74 | + } |
| 75 | + } |
| 76 | + } |
| 77 | + |
| 78 | + // The real ndarray CLAM/CHAODA engine. |
| 79 | + let tree = ClamTree::build(&data, vec_len, 4); |
| 80 | + let scores = tree.anomaly_scores(&data, vec_len); |
| 81 | + |
| 82 | + // HHTL seam = the inter-HEEL-basin cut (the crisp 2-split bottleneck). |
| 83 | + let keys = hhtl_keys(grid); |
| 84 | + let mut seam_nodes: HashSet<usize> = HashSet::new(); |
| 85 | + for e in &grid.edges { |
| 86 | + if keys[e.from].heel != keys[e.to].heel { |
| 87 | + seam_nodes.insert(e.from); |
| 88 | + seam_nodes.insert(e.to); |
| 89 | + } |
| 90 | + } |
| 91 | + |
| 92 | + let mean = |idxs: &[usize]| -> f64 { |
| 93 | + if idxs.is_empty() { |
| 94 | + 0.0 |
| 95 | + } else { |
| 96 | + idxs.iter().map(|&i| scores[i].score).sum::<f64>() / idxs.len() as f64 |
| 97 | + } |
| 98 | + }; |
| 99 | + let all: Vec<usize> = (0..n).collect(); |
| 100 | + let seam: Vec<usize> = seam_nodes.iter().copied().collect(); |
| 101 | + let (anom_seam, anom_all) = (mean(&seam), mean(&all)); |
| 102 | + let anom_ratio = if anom_all > 0.0 { anom_seam / anom_all } else { 0.0 }; |
| 103 | + |
| 104 | + // Top-quartile anomalies: are seam nodes over-represented? |
| 105 | + let mut ranked: Vec<usize> = (0..n).collect(); |
| 106 | + ranked.sort_by(|&a, &b| { |
| 107 | + scores[b] |
| 108 | + .score |
| 109 | + .partial_cmp(&scores[a].score) |
| 110 | + .unwrap_or(std::cmp::Ordering::Equal) |
| 111 | + }); |
| 112 | + let top_n = (n / 4).max(1); |
| 113 | + let top: HashSet<usize> = ranked[..top_n].iter().copied().collect(); |
| 114 | + let seam_in_top = seam.iter().filter(|i| top.contains(i)).count(); |
| 115 | + let expected = top_n as f64 * seam.len() as f64 / n as f64; // null expectation |
| 116 | + let lift = if expected > 0.0 { |
| 117 | + seam_in_top as f64 / expected |
| 118 | + } else { |
| 119 | + 0.0 |
| 120 | + }; |
| 121 | + |
| 122 | + println!("PROBE — CHAODA surge epicenter (ndarray::hpc::clam, real {country} grid)"); |
| 123 | + println!(" grid: {n} buses, {m} lines; {k}-bit Fiedler fingerprints ({vec_len} B each)"); |
| 124 | + println!( |
| 125 | + " seam (inter-HEEL-basin cut): {} nodes; CLAM tree built, anomaly_scores computed", |
| 126 | + seam.len() |
| 127 | + ); |
| 128 | + println!("\n anomaly (LFD-derived, higher = more anomalous):"); |
| 129 | + println!(" mean anomaly, SEAM nodes : {anom_seam:.3}"); |
| 130 | + println!(" mean anomaly, ALL nodes : {anom_all:.3}"); |
| 131 | + println!(" ratio (seam / all) : {anom_ratio:.2}"); |
| 132 | + println!( |
| 133 | + " seam share of top-{top_n} anomalies: {seam_in_top}/{} (expected {expected:.1}, lift {lift:.2})", |
| 134 | + seam.len() |
| 135 | + ); |
| 136 | + |
| 137 | + let supported = !seam.is_empty() && anom_ratio >= 1.30 && lift >= 1.30; |
| 138 | + println!("\n VERDICT:"); |
| 139 | + if supported { |
| 140 | + println!( |
| 141 | + " [SUPPORTED] the seam (the brittle fail-first cut) IS a CHAODA anomaly — its mean \ |
| 142 | + LFD anomaly is {anom_ratio:.2}× the global mean and it is {lift:.2}× over-represented \ |
| 143 | + in the top-quartile anomalies. The surge EPICENTER is detectable statically from the \ |
| 144 | + manifold geometry, no cascade simulation — the WHERE axis closes on real data." |
| 145 | + ); |
| 146 | + } else { |
| 147 | + println!( |
| 148 | + " [NOT SUPPORTED] seam anomaly ratio {anom_ratio:.2} / top-quartile lift {lift:.2} \ |
| 149 | + do not clear 1.30 — on this grid the LFD anomaly does not single out the seam. Report \ |
| 150 | + honestly: CHAODA flags geometric outliers, which need not coincide with the \ |
| 151 | + spectral-cut bottleneck here. Do not promote." |
| 152 | + ); |
| 153 | + } |
| 154 | + println!( |
| 155 | + " NOTE: CHAODA anomaly = leaf-cluster LFD normalized over the tree; the seam is the \ |
| 156 | + inter-HEEL-basin cut (the crisp 2-split). This is the WHERE axis (epicenter), distinct \ |
| 157 | + from the Weyl HOW-MUCH and the family-basin HOW-it-spreads axes." |
| 158 | + ); |
| 159 | +} |
0 commit comments