powerio_matrix/
pipeline.rs

1//! Orchestrates a single case → output directory.
2//!
3//! Given a parsed `Network`, builds the requested matrix family, writes
4//! `.mtx` files, and emits a `meta.json` sidecar describing what was
5//! produced. Used by both the `batch` CLI subcommand and the TUI's
6//! batch export screen.
7
8use std::path::{Path, PathBuf};
9
10use rand::SeedableRng;
11use rand::distr::{Distribution, StandardUniform};
12use rand_chacha::ChaCha8Rng;
13use serde::{Deserialize, Serialize};
14
15use crate::Result;
16use crate::indexed::IndexedNetwork;
17use crate::io::meta::{CaseMetadata, MatrixMetadata, write_meta_json};
18use crate::io::mtx::{write_mtx, write_vector_mtx};
19use crate::matrix::{
20    BuildOptions, MatrixStats, ZeroImpedanceRule, ZeroImpedanceSkips, build_adjacency,
21    build_bdoubleprime, build_bprime, build_lacpf, build_ybus, negate_into, sddm_check,
22    skipped_zero_impedance,
23};
24use crate::network::Network;
25
26#[derive(Debug, Clone, Copy, PartialEq, Eq, Serialize, Deserialize)]
27#[non_exhaustive]
28pub enum MatrixKind {
29    /// FDPF B' (shuntless, taps=1, shifts=0).
30    BPrime,
31    /// FDPF B'' (with shunts, taps; r=0 if BX scheme).
32    BDoublePrime,
33    /// `Re(Y_bus)` — full conductance matrix.
34    YbusG,
35    /// `-Im(Y_bus)` — full susceptance Laplacian (positive convention).
36    YbusB,
37    /// LACPF block: `[[G, -B], [-B, -G]]`, 2n × 2n indefinite.
38    Lacpf,
39    /// 0/1 bus adjacency matrix.
40    Adjacency,
41}
42
43impl MatrixKind {
44    pub const ALL: &'static [MatrixKind] = &[
45        Self::BPrime,
46        Self::BDoublePrime,
47        Self::YbusG,
48        Self::YbusB,
49        Self::Lacpf,
50        Self::Adjacency,
51    ];
52
53    pub fn slug(self) -> &'static str {
54        match self {
55            Self::BPrime => "bprime",
56            Self::BDoublePrime => "bdoubleprime",
57            Self::YbusG => "ybus_real",
58            Self::YbusB => "ybus_imag",
59            Self::Lacpf => "lacpf",
60            Self::Adjacency => "adjacency",
61        }
62    }
63
64    pub fn label(self) -> &'static str {
65        match self {
66            Self::BPrime => "B' (FDPF, shuntless)",
67            Self::BDoublePrime => "B'' (FDPF, with shunts)",
68            Self::YbusG => "Re(Y_bus)",
69            Self::YbusB => "-Im(Y_bus)",
70            Self::Lacpf => "LACPF block (2n×2n)",
71            Self::Adjacency => "adjacency (0/1)",
72        }
73    }
74}
75
76/// How to populate the RHS vector(s) emitted alongside each matrix.
77#[derive(Debug, Clone, Copy, PartialEq, Eq, Default, Serialize, Deserialize)]
78pub enum RhsKind {
79    #[default]
80    None,
81    /// Zero-mean Gaussian random (deterministic from `rng_seed`).
82    Random,
83    /// Power injections from the case: `b = (Pd, Qd) / baseMVA`.
84    Injection,
85}
86
87#[derive(Debug, Clone)]
88pub struct Pipeline {
89    pub matrices: Vec<MatrixKind>,
90    pub options: BuildOptions,
91    pub rhs: RhsKind,
92    pub rng_seed: u64,
93    pub source_file: Option<PathBuf>,
94}
95
96impl Default for Pipeline {
97    fn default() -> Self {
98        Self {
99            matrices: vec![MatrixKind::BPrime],
100            options: BuildOptions::default(),
101            rhs: RhsKind::None,
102            rng_seed: 0x00C0_FFEE,
103            source_file: None,
104        }
105    }
106}
107
108#[derive(Debug, Clone)]
109pub struct PipelineOutputs {
110    pub case_name: String,
111    pub files: Vec<PathBuf>,
112    pub metadata: CaseMetadata,
113}
114
115impl Pipeline {
116    pub fn run(&self, net: &Network, out_dir: impl AsRef<Path>) -> Result<PipelineOutputs> {
117        let out_dir = out_dir.as_ref();
118        std::fs::create_dir_all(out_dir)?;
119
120        let view = IndexedNetwork::new(net);
121
122        let mut files = Vec::new();
123        let mut matrices_meta = Vec::new();
124        let mut ybus_cache = None;
125
126        for &kind in &self.matrices {
127            let matrix_path = out_dir.join(format!("{}_{}.mtx", view.name(), kind.slug()));
128            let matrix = self.build_for_run(&view, kind, &mut ybus_cache)?;
129            write_mtx(&matrix, &matrix_path)?;
130            let stats = matrix_stats_for_kind(&matrix, &view, kind, &self.options);
131            let sddm = sddm_check(&matrix);
132            matrices_meta.push(MatrixMetadata {
133                kind: kind.slug().to_string(),
134                file: matrix_path
135                    .file_name()
136                    .and_then(|s| s.to_str())
137                    .unwrap_or("")
138                    .to_string(),
139                stats,
140                sddm,
141            });
142            files.push(matrix_path);
143
144            // RHS for matrices that take a RHS of length n (skip LACPF which is 2n).
145            if let Some(rhs) = self.build_rhs(&view, kind) {
146                let rhs_path = out_dir.join(format!("{}_{}_rhs.mtx", view.name(), kind.slug()));
147                write_vector_mtx(&rhs, &rhs_path)?;
148                files.push(rhs_path);
149            }
150        }
151
152        // Shunt vector as a sidecar (not always meaningful, but cheap).
153        let shunt_path = out_dir.join(format!("{}_shunt.mtx", view.name()));
154        let base = view.per_unit_base();
155        let shunt: Vec<f64> = view.bs().iter().map(|&b| b / base).collect();
156        write_vector_mtx(&shunt, &shunt_path)?;
157        files.push(shunt_path);
158
159        let metadata = CaseMetadata {
160            case_name: view.name().to_string(),
161            source_file: self
162                .source_file
163                .as_ref()
164                .and_then(|p| p.to_str())
165                .map(str::to_string),
166            source_sha256: self
167                .source_file
168                .as_ref()
169                .and_then(|p| std::fs::read(p).ok())
170                .map(|b| sha256_hex(&b)),
171            base_mva: view.base_mva(),
172            n_buses: view.n(),
173            n_branches: view.branches().len(),
174            build_options: self.options.clone(),
175            matrices: matrices_meta,
176            powerio_version: env!("CARGO_PKG_VERSION").to_string(),
177        };
178        let meta_path = out_dir.join(format!("{}_meta.json", view.name()));
179        write_meta_json(&metadata, &meta_path)?;
180        files.push(meta_path);
181
182        Ok(PipelineOutputs {
183            case_name: view.name().to_string(),
184            files,
185            metadata,
186        })
187    }
188
189    fn build_for_run(
190        &self,
191        case: &IndexedNetwork,
192        kind: MatrixKind,
193        ybus_cache: &mut Option<YbusCache>,
194    ) -> Result<sprs::CsMat<f64>> {
195        match kind {
196            MatrixKind::YbusG => take_ybus_g(case, &self.options, ybus_cache),
197            MatrixKind::YbusB => take_ybus_b(case, &self.options, ybus_cache),
198            _ => build_kind(case, kind, &self.options),
199        }
200    }
201
202    fn build_rhs(&self, case: &IndexedNetwork, kind: MatrixKind) -> Option<Vec<f64>> {
203        // No meaningful RHS for the 2n LACPF block or the structural adjacency.
204        if matches!(self.rhs, RhsKind::None)
205            || matches!(kind, MatrixKind::Lacpf | MatrixKind::Adjacency)
206        {
207            return None;
208        }
209        let n = case.n();
210        let v = match self.rhs {
211            RhsKind::Random => {
212                let mut rng = ChaCha8Rng::seed_from_u64(self.rng_seed.wrapping_add(kind as u64));
213                let dist = StandardUniform;
214                let mut v: Vec<f64> = (0..n)
215                    .map(|_| {
216                        let u: f64 = dist.sample(&mut rng);
217                        u - 0.5
218                    })
219                    .collect();
220                let mean = v.iter().sum::<f64>() / n as f64;
221                for x in &mut v {
222                    *x -= mean; // zero-mean for Laplacian compatibility
223                }
224                v
225            }
226            RhsKind::Injection => {
227                let base = case.per_unit_base();
228                match kind {
229                    MatrixKind::BPrime | MatrixKind::YbusG | MatrixKind::YbusB => {
230                        case.pd().iter().map(|&p| -p / base).collect()
231                    }
232                    MatrixKind::BDoublePrime => case.qd().iter().map(|&q| -q / base).collect(),
233                    MatrixKind::Lacpf | MatrixKind::Adjacency => unreachable!(),
234                }
235            }
236            RhsKind::None => unreachable!(),
237        };
238        Some(v)
239    }
240}
241
242struct YbusCache {
243    g: Option<sprs::CsMat<f64>>,
244    b: Option<sprs::CsMat<f64>>,
245}
246
247fn fill_ybus_cache(
248    view: &IndexedNetwork,
249    opts: &BuildOptions,
250    ybus_cache: &mut Option<YbusCache>,
251) -> Result<()> {
252    let parts = build_ybus(view, opts)?;
253    *ybus_cache = Some(YbusCache {
254        g: Some(parts.g),
255        b: Some(parts.b),
256    });
257    Ok(())
258}
259
260fn take_ybus_g(
261    view: &IndexedNetwork,
262    opts: &BuildOptions,
263    ybus_cache: &mut Option<YbusCache>,
264) -> Result<sprs::CsMat<f64>> {
265    if ybus_cache.as_ref().is_none_or(|c| c.g.is_none()) {
266        fill_ybus_cache(view, opts, ybus_cache)?;
267    }
268    Ok(ybus_cache
269        .as_mut()
270        .and_then(|c| c.g.take())
271        .expect("Ybus cache was just filled, so the real part must be present"))
272}
273
274fn take_ybus_b(
275    view: &IndexedNetwork,
276    opts: &BuildOptions,
277    ybus_cache: &mut Option<YbusCache>,
278) -> Result<sprs::CsMat<f64>> {
279    if ybus_cache.as_ref().is_none_or(|c| c.b.is_none()) {
280        fill_ybus_cache(view, opts, ybus_cache)?;
281    }
282    let b = ybus_cache
283        .as_mut()
284        .and_then(|c| c.b.take())
285        .expect("Ybus cache was just filled, so the imaginary part must be present");
286    Ok(negate_into(b))
287}
288
289/// Build the square matrix for one [`MatrixKind`] from an indexed network. The
290/// single dispatch shared by the [`Pipeline`], the `verify` CLI command, and the
291/// TUI inspect screen, so the `YbusB = -Im(Y_bus)` sign lives in one place.
292pub fn build_kind(
293    view: &IndexedNetwork,
294    kind: MatrixKind,
295    opts: &BuildOptions,
296) -> Result<sprs::CsMat<f64>> {
297    match kind {
298        MatrixKind::BPrime => build_bprime(view, opts),
299        MatrixKind::BDoublePrime => build_bdoubleprime(view, opts),
300        MatrixKind::YbusG => build_ybus(view, opts).map(|p| p.g),
301        MatrixKind::YbusB => build_ybus(view, opts).map(|p| negate_into(p.b)),
302        MatrixKind::Lacpf => build_lacpf(view, opts),
303        MatrixKind::Adjacency => build_adjacency(view),
304    }
305}
306
307pub fn zero_impedance_rule_for_kind(
308    kind: MatrixKind,
309    opts: &BuildOptions,
310) -> Option<ZeroImpedanceRule> {
311    match kind {
312        MatrixKind::BPrime => Some(match opts.scheme {
313            crate::matrix::Scheme::Bx => ZeroImpedanceRule::Series,
314            crate::matrix::Scheme::Xb => ZeroImpedanceRule::Reactance,
315        }),
316        MatrixKind::BDoublePrime => Some(match opts.scheme {
317            crate::matrix::Scheme::Bx => ZeroImpedanceRule::Reactance,
318            crate::matrix::Scheme::Xb => ZeroImpedanceRule::Series,
319        }),
320        MatrixKind::YbusG | MatrixKind::YbusB | MatrixKind::Lacpf => {
321            Some(ZeroImpedanceRule::Series)
322        }
323        MatrixKind::Adjacency => None,
324    }
325}
326
327pub fn zero_impedance_skips_for_kind(
328    view: &IndexedNetwork,
329    kind: MatrixKind,
330    opts: &BuildOptions,
331) -> ZeroImpedanceSkips {
332    if !opts.skip_zero_impedance {
333        return ZeroImpedanceSkips::default();
334    }
335    zero_impedance_rule_for_kind(kind, opts).map_or_else(ZeroImpedanceSkips::default, |rule| {
336        skipped_zero_impedance(view, rule)
337    })
338}
339
340pub fn matrix_stats_for_kind(
341    matrix: &sprs::CsMat<f64>,
342    view: &IndexedNetwork,
343    kind: MatrixKind,
344    opts: &BuildOptions,
345) -> MatrixStats {
346    MatrixStats::from_csr(matrix)
347        .with_zero_impedance_skips(zero_impedance_skips_for_kind(view, kind, opts))
348}
349
350fn sha256_hex(bytes: &[u8]) -> String {
351    use sha2::{Digest, Sha256};
352    use std::fmt::Write as _;
353    let mut hasher = Sha256::new();
354    hasher.update(bytes);
355    let digest = hasher.finalize();
356    let mut out = String::with_capacity(digest.len() * 2);
357    for byte in digest {
358        let _ = write!(out, "{byte:02x}");
359    }
360    out
361}
powerio_matrix/pipeline.rs

powerio_matrix/
pipeline.rs
