gridoxide_core/

jacobian.rs

//! Topology-derived Jacobian assembly: analyze the sparsity pattern once,
//! then refill values every iteration at fixed offsets.
//!
//! `solver::build_jacobian_triplets` rebuilds a fresh `Vec<(usize, usize,
//! f64)>` every Newton iteration, re-deriving row and column indices that are
//! by construction identical each time — the sparsity pattern is fixed for
//! the lifetime of a topology, which is the very property
//! `solver::LinearSolver` relies on to cache its symbolic factorization.
//! [`JacobianPattern`] hoists all of that out of the loop: the `(row, col)`
//! arrays and a per-entry "recipe" are computed once by
//! [`analyze`](JacobianPattern::analyze), and [`fill`](JacobianPattern::fill)
//! then writes only the `f64` values into a caller-owned buffer.
//!
//! Two payoffs:
//!
//! - **On the CPU**, per-iteration allocation and index arithmetic disappear.
//!   Jacobian assembly is 36–41% of iteration time on the cases measured.
//! - **On a GPU**, this *is* the kernel shape. [`Entry`] is a flat,
//!   fixed-size, branch-light record; `fill` is one independent write per
//!   entry into a preallocated array at a precomputed offset, with all reads
//!   gathers from small per-bus arrays. Assembly as one flat kernel is
//!   exactly this layout, extended to a batch by adding a scenario stride.
//!
//! [`JacobianPattern::fill`] is a transliteration of `build_jacobian_triplets`' H/N/M/L
//! formulas, not a rederivation. The one difference is how it reaches
//! `sin(θi − θk)` and `cos(θi − θk)`: from each bus's own sine and cosine,
//! once per bus, rather than per entry. That rounds differently, so the tests
//! below hold the two to the same positions exactly and to the same values
//! within a few units in the last place.

use num_complex::Complex;

use crate::network::{BusVoltages, YBusSparse};
use crate::types::{Bus, BusType};

/// Which of the Jacobian's eight distinct formulas an [`Entry`] evaluates.
///
/// ```text
/// J = [ H  N ]   H = dP/d_ang, N = dP/d_vmag
///     [ M  L ]   M = dQ/d_ang, L = dQ/d_vmag
/// ```
///
/// Diagonal and off-diagonal cases are separate variants because they read
/// different inputs, not merely different indices: the diagonal terms depend
/// on `p_calc`/`q_calc` at the bus, the off-diagonal ones on the angle
/// difference to a neighbor.
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
#[repr(u8)]
pub enum EntryKind {
    Hii = 0,
    Nii = 1,
    Hik = 2,
    Nik = 3,
    Mii = 4,
    Lii = 5,
    Mik = 6,
    Lik = 7,
}

/// One Jacobian nonzero's recipe. Deliberately flat and `Copy` — this is the
/// record a GPU kernel would read one-per-thread.
#[derive(Clone, Copy, Debug)]
pub struct Entry {
    pub kind: EntryKind,
    /// Row bus (physical index).
    pub i: u32,
    /// Column bus (physical index). Equals `i` for the diagonal kinds.
    pub k: u32,
    /// `Y[i][k]`, captured at analyze time. The Y-bus is fixed for the
    /// lifetime of a topology, so carrying the admittance here turns what
    /// would be a pointer-chase through `YBusSparse::row` into a value the
    /// entry already owns.
    pub y: Complex<f64>,
}

/// A topology's Jacobian sparsity pattern plus the per-entry recipes needed
/// to refill its values. Build once per topology via
/// [`analyze`](Self::analyze); call [`fill`](Self::fill) every iteration.
///
/// `rows`/`cols`/`entries` are parallel arrays of the same length, in the
/// exact emission order `solver::build_jacobian_triplets` used, so a cached
/// symbolic factorization keyed on that order stays valid — until
/// [`permute`](Self::permute) puts them in a sparse LU's own storage order.
pub struct JacobianPattern {
    pub n_unknowns: usize,
    rows: Vec<u32>,
    cols: Vec<u32>,
    entries: Vec<Entry>,
    /// Where each entry's value is in the array [`fill`](Self::fill)
    /// computes first, parallel to `entries`: the fill itself is then one
    /// gather per entry, with no branch on its kind.
    src: Vec<u32>,
    /// Row buses in non-slack order with their `Y[i][i]` (zero if the Y-bus
    /// has no diagonal there); bus `d` owns slots `4d..4d + 4` (H, N, M, L).
    diag: Vec<(u32, Complex<f64>)>,
    /// Off-diagonal `(i, k)` pairs with `i` a row bus, in Y-bus order; pair
    /// `p` owns slots `4 · diag.len() + 4p ..` (H, N, M, L).
    pairs: Vec<(u32, u32, Complex<f64>)>,
}

impl JacobianPattern {
    /// Derives the pattern from a topology and a bus-type assignment.
    ///
    /// Must be recomputed if either changes — a different set of branches, or
    /// any bus switching `PV`↔`PQ` (which changes `n_unknowns` itself). This
    /// is the same validity condition `solver::PersistentSolver` documents
    /// for its cached factorization, and the two are reset together.
    pub fn analyze(buses: &[Bus], ybus: &YBusSparse) -> Self {
        let non_slack_idx: Vec<usize> = buses
            .iter()
            .filter(|b| !matches!(b.bus_type, BusType::Slack))
            .map(|b| b.idx)
            .collect();
        let pq_idx: Vec<usize> = buses
            .iter()
            .filter(|b| matches!(b.bus_type, BusType::PQ))
            .map(|b| b.idx)
            .collect();

        let n_angle = non_slack_idx.len();
        let n_unknowns = n_angle + pq_idx.len();

        let mut non_slack_pos: Vec<Option<usize>> = vec![None; buses.len()];
        for (pos, &i) in non_slack_idx.iter().enumerate() {
            non_slack_pos[i] = Some(pos);
        }
        let mut pq_pos: Vec<Option<usize>> = vec![None; buses.len()];
        for (pos, &i) in pq_idx.iter().enumerate() {
            pq_pos[i] = Some(pos);
        }

        let capacity: usize = non_slack_idx
            .iter()
            .map(|&i| ybus.row(i).len())
            .sum::<usize>()
            * 2
            + pq_idx.iter().map(|&i| ybus.row(i).len()).sum::<usize>() * 2;
        let mut rows = Vec::with_capacity(capacity);
        let mut cols = Vec::with_capacity(capacity);
        let mut entries = Vec::with_capacity(capacity);
        let mut src = Vec::with_capacity(capacity);
        // Slots: four per non-slack bus first (its diagonal H, N, M, L), then
        // four per off-diagonal pair. `slot` is the entry's offset within its
        // four, the kind's position in H, N, M, L.
        let pair_base = 4 * n_angle;
        let mut push = |row: usize,
                        col: usize,
                        kind: EntryKind,
                        i: usize,
                        k: usize,
                        y: Complex<f64>,
                        slot: usize| {
            rows.push(row as u32);
            cols.push(col as u32);
            entries.push(Entry {
                kind,
                i: i as u32,
                k: k as u32,
                y,
            });
            src.push(slot as u32);
        };
        let has_column = |k: usize| non_slack_pos[k].is_some() || pq_pos[k].is_some();

        // H and N blocks: rows = non_slack_idx (P-mismatch equations). Every
        // pair is met here first and numbered; a PQ row's M and L walk the
        // same neighbors in the same order below.
        let mut diag: Vec<(u32, Complex<f64>)> = non_slack_idx
            .iter()
            .map(|&i| (i as u32, Complex::new(0.0, 0.0)))
            .collect();
        let mut pairs: Vec<(u32, u32, Complex<f64>)> = Vec::new();
        let mut first_pair = vec![0usize; buses.len()];
        for (row_idx, &i) in non_slack_idx.iter().enumerate() {
            first_pair[i] = pairs.len();
            let d = 4 * row_idx;
            for &(k, y_ik) in ybus.row(i) {
                if k == i {
                    diag[row_idx].1 = y_ik;
                    push(row_idx, row_idx, EntryKind::Hii, i, i, y_ik, d);
                    if let Some(col_idx) = pq_pos[i] {
                        push(
                            row_idx,
                            n_angle + col_idx,
                            EntryKind::Nii,
                            i,
                            i,
                            y_ik,
                            d + 1,
                        );
                    }
                    continue;
                }
                if !has_column(k) {
                    continue;
                }
                let at = pair_base + 4 * pairs.len();
                pairs.push((i as u32, k as u32, y_ik));
                if let Some(col_idx) = non_slack_pos[k] {
                    push(row_idx, col_idx, EntryKind::Hik, i, k, y_ik, at);
                }
                if let Some(col_idx) = pq_pos[k] {
                    push(
                        row_idx,
                        n_angle + col_idx,
                        EntryKind::Nik,
                        i,
                        k,
                        y_ik,
                        at + 1,
                    );
                }
            }
        }

        // M and L blocks: rows = pq_idx (Q-mismatch equations).
        for (row_idx, &i) in pq_idx.iter().enumerate() {
            let d = 4 * non_slack_pos[i].expect("a PQ bus is always non-slack");
            let mut pair = first_pair[i];
            for &(k, y_ik) in ybus.row(i) {
                if k == i {
                    // Column is i's position among *all* non-slack buses, not
                    // just PQ ones — these differ once PV buses exist.
                    let col_idx = non_slack_pos[i].expect("a PQ bus is always non-slack");
                    push(
                        n_angle + row_idx,
                        col_idx,
                        EntryKind::Mii,
                        i,
                        i,
                        y_ik,
                        d + 2,
                    );
                    push(
                        n_angle + row_idx,
                        n_angle + row_idx,
                        EntryKind::Lii,
                        i,
                        i,
                        y_ik,
                        d + 3,
                    );
                    continue;
                }
                if !has_column(k) {
                    continue;
                }
                debug_assert_eq!((pairs[pair].0, pairs[pair].1), (i as u32, k as u32));
                let at = pair_base + 4 * pair;
                if let Some(col_idx) = non_slack_pos[k] {
                    push(
                        n_angle + row_idx,
                        col_idx,
                        EntryKind::Mik,
                        i,
                        k,
                        y_ik,
                        at + 2,
                    );
                }
                if let Some(col_idx) = pq_pos[k] {
                    push(
                        n_angle + row_idx,
                        n_angle + col_idx,
                        EntryKind::Lik,
                        i,
                        k,
                        y_ik,
                        at + 3,
                    );
                }
                pair += 1;
            }
        }

        Self {
            n_unknowns,
            rows,
            cols,
            entries,
            src,
            diag,
            pairs,
        }
    }

    /// Reorders the entries so that entry `k` is the one that was at
    /// `order[k]`: from then on [`fill`](Self::fill) writes, and
    /// [`to_triplets`](Self::to_triplets) lists, in that order. What
    /// `lu::SparseLu::adopt_storage_order` returns, so the fill writes
    /// straight into the factorization's own layout.
    ///
    /// # Panics
    ///
    /// If `order` is not a permutation of `0..len()`.
    pub fn permute(&mut self, order: &[usize]) {
        assert_eq!(
            order.len(),
            self.entries.len(),
            "order must cover every entry"
        );
        let mut seen = vec![false; order.len()];
        for &i in order {
            assert!(
                !std::mem::replace(&mut seen[i], true),
                "order repeats entry {i}"
            );
        }
        self.rows = order.iter().map(|&i| self.rows[i]).collect();
        self.cols = order.iter().map(|&i| self.cols[i]).collect();
        self.entries = order.iter().map(|&i| self.entries[i]).collect();
        self.src = order.iter().map(|&i| self.src[i]).collect();
    }

    /// Number of Jacobian nonzeros — the length of `values` that
    /// [`fill`](Self::fill) expects.
    pub fn len(&self) -> usize {
        self.entries.len()
    }

    pub fn is_empty(&self) -> bool {
        self.entries.is_empty()
    }

    pub fn rows(&self) -> &[u32] {
        &self.rows
    }

    pub fn cols(&self) -> &[u32] {
        &self.cols
    }

    pub fn entries(&self) -> &[Entry] {
        &self.entries
    }

    /// Writes this iteration's Jacobian values into `values`, which is
    /// resized to [`len`](Self::len) and otherwise fully overwritten.
    ///
    /// Every write is independent and lands at a fixed offset — no
    /// allocation, no index arithmetic, no branching beyond the eight-way
    /// `kind` dispatch. That independence is what makes this portable to one
    /// thread per entry.
    pub fn fill(&self, buses: &[Bus], p_calc: &[f64], q_calc: &[f64], values: &mut Vec<f64>) {
        values.clear();
        self.fill_into(buses, p_calc, q_calc, values);
    }

    /// [`fill`](Self::fill) without the clear — *appends* this scenario's
    /// values to `values`. What lets `bde::BlockDiagonal` concatenate B
    /// scenarios' blocks into one flat array for a block-diagonal solve.
    pub fn fill_into(&self, buses: &[Bus], p_calc: &[f64], q_calc: &[f64], values: &mut Vec<f64>) {
        self.fill_from(&BusVoltages::of(buses), p_calc, q_calc, values);
    }

    /// [`fill`](Self::fill) from voltages already taken, which the Newton
    /// loop shares with the power injections.
    pub fn fill_with(
        &self,
        v: &BusVoltages,
        p_calc: &[f64],
        q_calc: &[f64],
        values: &mut Vec<f64>,
    ) {
        values.clear();
        self.fill_from(v, p_calc, q_calc, values);
    }

    /// Appends, as [`fill_into`](Self::fill_into).
    ///
    /// In two passes. The first computes every value the Jacobian can hold,
    /// in the order the Y-bus stores them: for each non-slack bus its
    /// diagonal H, N, M and L, and for each off-diagonal pair its four,
    /// which share their sine and cosine — formed from the two buses' own
    /// rather than by calling the trigonometric functions per entry, which
    /// made `sincos` the largest cost of an iteration. The second copies
    /// them into the entries through the precomputed `src`, in whatever
    /// order the entries are in.
    ///
    /// The point of the split is that neither pass branches on an entry's
    /// kind. Once the entries are in a factorization's column order
    /// ([`permute`](Self::permute)) the kinds interleave unpredictably, and
    /// a per-entry `match` cost more than the arithmetic: this fill is
    /// several times faster than one, on the same values.
    fn fill_from(&self, v: &BusVoltages, p_calc: &[f64], q_calc: &[f64], values: &mut Vec<f64>) {
        let vm = &v.vm;
        let mut slots = Vec::with_capacity(4 * (self.diag.len() + self.pairs.len()));
        for &(i, y) in &self.diag {
            let i = i as usize;
            let (vm_i, g, b) = (vm[i], y.re, y.im);
            // `powi(2)` rather than `vm_i * vm_i` to stay textually
            // identical to `build_jacobian_triplets`, so the diagonal agrees
            // with it to the last bit.
            slots.extend([
                // H_ii = -Q_i - V_i^2 * B_ii
                -q_calc[i] - vm_i.powi(2) * b,
                // N_ii = P_i/V_i + V_i * G_ii
                p_calc[i] / vm_i + vm_i * g,
                // M_ii = P_i - V_i^2 * G_ii
                p_calc[i] - vm_i.powi(2) * g,
                // L_ii = Q_i/V_i - V_i * B_ii
                q_calc[i] / vm_i - vm_i * b,
            ]);
        }
        for &(i, k, y) in &self.pairs {
            let (i, k) = (i as usize, k as usize);
            let (vm_i, vm_k, g, b) = (vm[i], vm[k], y.re, y.im);
            let (s_i, c_i, s_k, c_k) = (v.sin[i], v.cos[i], v.sin[k], v.cos[k]);
            let sin = s_i * c_k - c_i * s_k;
            let cos = c_i * c_k + s_i * s_k;
            slots.extend([
                // H_ik = V_i * V_k * (G_ik * sin - B_ik * cos)
                vm_i * vm_k * (g * sin - b * cos),
                // N_ik = V_i * (G_ik * cos + B_ik * sin)
                vm_i * (g * cos + b * sin),
                // M_ik = -V_i * V_k * (G_ik * cos + B_ik * sin)
                -vm_i * vm_k * (g * cos + b * sin),
                // L_ik = V_i * (G_ik * sin - B_ik * cos)
                vm_i * (g * sin - b * cos),
            ]);
        }
        values.extend(self.src.iter().map(|&s| slots[s as usize]));
    }

    /// Appends an identity block in this pattern's own layout: `1.0` at every
    /// stored diagonal position, `0.0` everywhere else.
    ///
    /// This is how `bde::BlockDiagonal` masks a converged or diverged
    /// scenario out of a batch. Dropping its block outright would change the
    /// stacked matrix's sparsity pattern and invalidate the cached symbolic
    /// factorization — the whole reason batching is cheap. Writing an
    /// identity into the *same* stored positions keeps the pattern
    /// bit-identical while making that scenario's update exactly zero (with a
    /// zero right-hand side), and identity is perfectly conditioned so it
    /// cannot degrade the factorization.
    ///
    /// Every diagonal position is structurally present: `Hii` supplies
    /// `(r, r)` for every angle row and `Lii` supplies it for every
    /// magnitude row, so no fill-in is required to write this.
    pub fn fill_identity_into(&self, values: &mut Vec<f64>) {
        values.reserve(self.entries.len());
        for k in 0..self.entries.len() {
            values.push(if self.rows[k] == self.cols[k] {
                1.0
            } else {
                0.0
            });
        }
    }

    /// Convenience for callers that still want `(row, col, value)` triplets,
    /// which is the shape every `solver::LinearSolver` backend consumes.
    pub fn to_triplets(&self, values: &[f64]) -> Vec<(usize, usize, f64)> {
        self.rows
            .iter()
            .zip(&self.cols)
            .zip(values)
            .map(|((&r, &c), &v)| (r as usize, c as usize, v))
            .collect()
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::network::{build_ybus, power_injections};
    use crate::types::Line;

    /// Builds a small mixed-type network: bus 0 Slack, bus 1 PV, buses 2-4 PQ.
    /// The PV bus matters — it is what makes a PQ row's `M_ii` column index
    /// differ from its row index, the one place the offset bookkeeping is
    /// genuinely easy to get wrong.
    fn mixed_network() -> (Vec<Bus>, YBusSparse) {
        let bus = |idx, bus_type, vm, p, q| Bus {
            idx,
            bus_type,
            voltage_mag: vm,
            voltage_ang: 0.0,
            p_spec: p,
            q_spec: q,
            q_min: -9.0,
            q_max: 9.0,
            u_rated: 0.0,
            zip_terms: Vec::new(),
        };
        let buses = vec![
            bus(0, BusType::Slack, 1.06, 0.0, 0.0),
            bus(1, BusType::PV, 1.04, 0.4, 0.0),
            bus(2, BusType::PQ, 1.0, -0.6, -0.2),
            bus(3, BusType::PQ, 1.0, -0.3, -0.1),
            bus(4, BusType::PQ, 1.0, -0.5, -0.15),
        ];
        let line = |from, to, r, x, b| Line {
            from,
            to,
            r,
            x,
            b_shunt: b,
            g_shunt: 0.0,
        };
        let lines = vec![
            line(0, 1, 0.02, 0.06, 0.03),
            line(0, 2, 0.08, 0.24, 0.025),
            line(1, 2, 0.06, 0.18, 0.02),
            line(1, 3, 0.06, 0.18, 0.02),
            line(2, 4, 0.04, 0.12, 0.015),
            line(3, 4, 0.01, 0.03, 0.01),
        ];
        let ybus = build_ybus(buses.len(), &lines, &[]).finish();
        (buses, ybus)
    }

    /// Values agree with the reference to rounding: the fill forms the
    /// angle difference's sine and cosine from the buses' own, which rounds
    /// differently from evaluating them at the difference.
    fn assert_same_value(got: f64, expected: f64, what: &str) {
        assert!(
            (got - expected).abs() <= 1e-13 * expected.abs().max(1.0),
            "{what}: {got} vs {expected}"
        );
    }

    /// The precomputed-offset assembler must reproduce the reference
    /// implementation's `(row, col, value)` sequence — same order and same
    /// indices exactly, values to rounding (see `assert_same_value`).
    #[test]
    fn pattern_matches_reference() {
        let (mut buses, ybus) = mixed_network();
        let pattern = JacobianPattern::analyze(&buses, &ybus);
        let mut values = Vec::new();

        // Several distinct operating points: a flat start, then perturbed
        // angles and magnitudes, so agreement can't be an artifact of every
        // sin/cos being evaluated at zero.
        for step in 0..4 {
            for (j, b) in buses.iter_mut().enumerate() {
                b.voltage_ang = -0.01 * (step * (j + 1)) as f64;
                if b.bus_type == BusType::PQ {
                    b.voltage_mag = 1.0 - 0.013 * step as f64 * (j + 1) as f64;
                }
            }
            let (p_calc, q_calc) = power_injections(&buses, &ybus);

            let expected = jacobian_triplets_reference(&buses, &ybus, &p_calc, &q_calc);
            pattern.fill(&buses, &p_calc, &q_calc, &mut values);
            let got = pattern.to_triplets(&values);

            assert_eq!(got.len(), expected.len(), "nonzero count at step {step}");
            for (idx, (g, e)) in got.iter().zip(&expected).enumerate() {
                assert_eq!(
                    (g.0, g.1),
                    (e.0, e.1),
                    "entry {idx} position at step {step}"
                );
                assert_same_value(g.2, e.2, &format!("entry {idx} value at step {step}"));
            }
        }
    }

    /// The same gate on the repo's committed 3-bus fixture, which has no PV
    /// bus — exercises the all-PQ path where row and column indices coincide.
    #[test]
    fn pattern_matches_reference_on_network_json() {
        let mut path = std::path::PathBuf::from(env!("CARGO_MANIFEST_DIR"));
        // The fixture is the workspace root's, two levels above this crate.
        path.push("../../tests/data/network.json");
        let raw = std::fs::read_to_string(path).expect("read network.json");
        let network: crate::json::NetworkData = serde_json::from_str(&raw).expect("parse");
        let buses = network.buses;
        let ybus = build_ybus(buses.len(), &network.lines, &[]).finish();

        let (p_calc, q_calc) = power_injections(&buses, &ybus);
        let expected = jacobian_triplets_reference(&buses, &ybus, &p_calc, &q_calc);

        let pattern = JacobianPattern::analyze(&buses, &ybus);
        let mut values = Vec::new();
        pattern.fill(&buses, &p_calc, &q_calc, &mut values);
        let got = pattern.to_triplets(&values);

        assert_eq!(got.len(), expected.len());
        for (idx, (g, e)) in got.iter().zip(&expected).enumerate() {
            assert_eq!((g.0, g.1), (e.0, e.1), "entry {idx} position");
            assert_same_value(g.2, e.2, &format!("entry {idx} value"));
        }
    }

    /// `permute` reorders what `fill` writes and `to_triplets` lists, and
    /// nothing else.
    #[test]
    fn permute_reorders_the_entries() {
        let (buses, ybus) = mixed_network();
        let (p_calc, q_calc) = power_injections(&buses, &ybus);
        let mut pattern = JacobianPattern::analyze(&buses, &ybus);
        let mut values = Vec::new();
        pattern.fill(&buses, &p_calc, &q_calc, &mut values);
        let before = pattern.to_triplets(&values);

        let order: Vec<usize> = (0..pattern.len()).rev().collect();
        pattern.permute(&order);
        pattern.fill(&buses, &p_calc, &q_calc, &mut values);
        let after = pattern.to_triplets(&values);
        let expected: Vec<_> = order.iter().map(|&i| before[i]).collect();
        assert_eq!(after, expected);
    }

    /// `fill` must fully overwrite its buffer, so reusing one across
    /// iterations (the entire point) can never leak stale values.
    #[test]
    fn fill_overwrites_a_reused_buffer() {
        let (buses, ybus) = mixed_network();
        let pattern = JacobianPattern::analyze(&buses, &ybus);
        let (p_calc, q_calc) = power_injections(&buses, &ybus);

        let mut fresh = Vec::new();
        pattern.fill(&buses, &p_calc, &q_calc, &mut fresh);

        let mut dirty = vec![f64::NAN; pattern.len() * 3];
        pattern.fill(&buses, &p_calc, &q_calc, &mut dirty);

        assert_eq!(dirty.len(), pattern.len());
        assert_eq!(dirty, fresh);
    }

    /// Derives the index bookkeeping `build_jacobian_triplets` needs and calls
    /// it, so a caller holding only `(buses, ybus, p_calc, q_calc)` can obtain
    /// the Jacobian triplets directly.
    ///
    /// Exists as the **reference oracle** for `jacobian::JacobianPattern`, whose
    /// precomputed-offset assembler must agree with this to rounding (see
    /// `jacobian`'s own tests). Kept deliberately naive — it re-derives
    /// everything from scratch on every call, which is exactly the cost
    /// `JacobianPattern` was written to remove, and exactly what makes it a
    /// trustworthy oracle.
    fn jacobian_triplets_reference(
        buses: &[Bus],
        ybus: &YBusSparse,
        p_calc: &[f64],
        q_calc: &[f64],
    ) -> Vec<(usize, usize, f64)> {
        let non_slack_idx: Vec<usize> = buses
            .iter()
            .filter(|b| !matches!(b.bus_type, BusType::Slack))
            .map(|b| b.idx)
            .collect();
        let pq_idx: Vec<usize> = buses
            .iter()
            .filter(|b| matches!(b.bus_type, BusType::PQ))
            .map(|b| b.idx)
            .collect();

        let n_angle = non_slack_idx.len();
        let mut non_slack_pos: Vec<Option<usize>> = vec![None; buses.len()];
        for (pos, &i) in non_slack_idx.iter().enumerate() {
            non_slack_pos[i] = Some(pos);
        }
        let mut pq_pos: Vec<Option<usize>> = vec![None; buses.len()];
        for (pos, &i) in pq_idx.iter().enumerate() {
            pq_pos[i] = Some(pos);
        }
        let cap = jacobian_triplet_capacity(ybus, &non_slack_idx, &pq_idx);

        build_jacobian_triplets(
            buses,
            ybus,
            &non_slack_idx,
            &pq_idx,
            &non_slack_pos,
            &pq_pos,
            n_angle,
            p_calc,
            q_calc,
            cap,
        )
    }

    /// Upper bound on `build_jacobian_triplets`' output length for a given
    /// topology — computed once per solve (not per iteration) and reused via
    /// `Vec::with_capacity`, since the sparsity pattern, and hence this bound,
    /// never changes across a solve's iterations while `build_jacobian_triplets`
    /// itself starts a fresh `Vec` from scratch every iteration. Each Y-bus
    /// neighbor contributes at most 2 triplets to the H/N block (one per
    /// non-slack row) and at most 2 to the M/L block (one per PQ row).
    fn jacobian_triplet_capacity(
        ybus: &YBusSparse,
        non_slack_idx: &[usize],
        pq_idx: &[usize],
    ) -> usize {
        let non_slack_degree: usize = non_slack_idx.iter().map(|&i| ybus.row(i).len()).sum();
        let pq_degree: usize = pq_idx.iter().map(|&i| ybus.row(i).len()).sum();
        2 * non_slack_degree + 2 * pq_degree
    }

    /// The **reference implementation** of Jacobian assembly, retained as the
    /// oracle `jacobian::JacobianPattern` is checked against. The
    /// Newton loop itself no longer calls this — it refills values at
    /// precomputed offsets instead — so it lives here, beside the tests it serves.
    /// Kept deliberately straightforward: its value is being obviously correct,
    /// not being fast.
    ///
    /// Assembles the Newton-Raphson Jacobian's nonzero entries as `(row, col,
    /// value)` triplets, walking only each unknown bus's actual Y-bus neighbors
    /// (`ybus.row(i)`) instead of the full cross product of unknown-bus indices.
    /// Every topological neighbor is always included regardless of its computed
    /// value (even if it happens to evaluate to exactly zero at some iteration),
    /// so the resulting sparsity *pattern* is identical across calls — required
    /// for `RealSparseSystem`'s cached symbolic factorization to stay valid.
    ///
    /// Block structure (H/N/M/L), matching the original dense Jacobian:
    /// ```text
    /// J = [ H  N ]   H = dP/d_ang, N = dP/d_vmag
    ///     [ M  L ]   M = dQ/d_ang, L = dQ/d_vmag
    /// ```
    #[allow(clippy::too_many_arguments)]
    fn build_jacobian_triplets(
        buses: &[Bus],
        ybus: &YBusSparse,
        non_slack_idx: &[usize],
        pq_idx: &[usize],
        non_slack_pos: &[Option<usize>],
        pq_pos: &[Option<usize>],
        n_angle: usize,
        p_calc: &[f64],
        q_calc: &[f64],
        triplet_capacity: usize,
    ) -> Vec<(usize, usize, f64)> {
        let vm: Vec<f64> = buses.iter().map(|b| b.voltage_mag).collect();
        let va: Vec<f64> = buses.iter().map(|b| b.voltage_ang).collect();
        let mut triplets = Vec::with_capacity(triplet_capacity);

        // H and N blocks: rows = non_slack_idx (P-mismatch equations).
        for (row_idx, &i) in non_slack_idx.iter().enumerate() {
            for &(k, y_ik) in ybus.row(i) {
                if k == i {
                    // H_ii = -Q_i - V_i^2 * B_ii
                    triplets.push((row_idx, row_idx, -q_calc[i] - vm[i].powi(2) * y_ik.im));
                    // N_ii = P_i/V_i + V_i * G_ii (only if i has a magnitude unknown)
                    if let Some(col_idx) = pq_pos[i] {
                        triplets.push((
                            row_idx,
                            n_angle + col_idx,
                            p_calc[i] / vm[i] + vm[i] * y_ik.re,
                        ));
                    }
                    continue;
                }
                let angle_ik = va[i] - va[k];
                if let Some(col_idx) = non_slack_pos[k] {
                    // H_ik = V_i * V_k * (G_ik * sin(d_ik) - B_ik * cos(d_ik))
                    let h_ik =
                        vm[i] * vm[k] * (y_ik.re * angle_ik.sin() - y_ik.im * angle_ik.cos());
                    triplets.push((row_idx, col_idx, h_ik));
                }
                if let Some(col_idx) = pq_pos[k] {
                    // N_ik = V_i * (G_ik * cos(d_ik) + B_ik * sin(d_ik))
                    let n_ik = vm[i] * (y_ik.re * angle_ik.cos() + y_ik.im * angle_ik.sin());
                    triplets.push((row_idx, n_angle + col_idx, n_ik));
                }
            }
        }

        // M and L blocks: rows = pq_idx (Q-mismatch equations).
        for (row_idx, &i) in pq_idx.iter().enumerate() {
            for &(k, y_ik) in ybus.row(i) {
                if k == i {
                    // M_ii = P_i - V_i^2 * G_ii; column is i's position among
                    // *all* non-slack buses, not just PQ ones (may differ from
                    // row_idx once PV buses exist).
                    let col_idx = non_slack_pos[i].expect("a PQ bus is always non-slack");
                    triplets.push((
                        n_angle + row_idx,
                        col_idx,
                        p_calc[i] - vm[i].powi(2) * y_ik.re,
                    ));
                    // L_ii = Q_i/V_i - V_i * B_ii
                    triplets.push((
                        n_angle + row_idx,
                        n_angle + row_idx,
                        q_calc[i] / vm[i] - vm[i] * y_ik.im,
                    ));
                    continue;
                }
                let angle_ik = va[i] - va[k];
                if let Some(col_idx) = non_slack_pos[k] {
                    // M_ik = -V_i * V_k * (G_ik * cos(d_ik) + B_ik * sin(d_ik))
                    let m_ik =
                        -vm[i] * vm[k] * (y_ik.re * angle_ik.cos() + y_ik.im * angle_ik.sin());
                    triplets.push((n_angle + row_idx, col_idx, m_ik));
                }
                if let Some(col_idx) = pq_pos[k] {
                    // L_ik = V_i * (G_ik * sin(d_ik) - B_ik * cos(d_ik))
                    let l_ik = vm[i] * (y_ik.re * angle_ik.sin() - y_ik.im * angle_ik.cos());
                    triplets.push((n_angle + row_idx, n_angle + col_idx, l_ik));
                }
            }
        }

        triplets
    }
}