gridoxide 0.1.0

use super::branch_flow::Terminal;
use super::network::{
    ShuntAdm, ShuntAdm3Ph, YBusSparse, branch_calc_param, half_open_branch_shunt,
    source_impedance_pu, source_impedance_pu_seq, tap_ratio_from_voltages,
    three_winding_star_params, transformer_admittances, transformer_admittances_ex,
    transformer_seq_params, transformer_tap,
};
use super::types::{Bus, BusType, Line, Line3Ph, Transformer, Transformer3PhSeq, ZipKind, ZipTerm};
use num_complex::Complex;
use serde::Deserialize;
use std::collections::HashMap;

pub mod measurement;
#[cfg(feature = "opf")]
pub mod opf;
pub mod se;
pub mod shortcircuit;

// ── Input structs ─────────────────────────────────────────────────────────────

#[derive(Clone, Deserialize)]
pub struct PgmInput {
    pub data: PgmData,
}

#[derive(Clone, Deserialize)]
pub struct PgmData {
    pub node: Vec<PgmNode>,
    #[serde(default)]
    pub line: Vec<PgmLine>,
    pub source: Vec<PgmSource>,
    #[serde(default)]
    pub sym_load: Vec<PgmSymLoad>,
    #[serde(default)]
    pub asym_load: Vec<PgmAsymLoad>,
    #[serde(default)]
    pub sym_gen: Vec<PgmSymGen>,
    #[serde(default)]
    pub asym_gen: Vec<PgmAsymGen>,
    #[serde(default)]
    pub shunt: Vec<PgmShunt>,
    #[serde(default)]
    pub transformer: Vec<PgmTransformer>,
    #[serde(default)]
    pub three_winding_transformer: Vec<PgmThreeWindingTransformer>,
    #[serde(default)]
    pub voltage_regulator: Vec<PgmVoltageRegulator>,
    #[serde(default)]
    pub link: Vec<PgmLink>,
    #[serde(default)]
    pub sym_voltage_sensor: Vec<PgmSymVoltageSensor>,
    #[serde(default)]
    pub sym_power_sensor: Vec<PgmSymPowerSensor>,
    #[serde(default)]
    pub asym_voltage_sensor: Vec<PgmAsymVoltageSensor>,
    #[serde(default)]
    pub asym_power_sensor: Vec<PgmAsymPowerSensor>,
    #[serde(default)]
    pub sym_current_sensor: Vec<PgmSymCurrentSensor>,
    #[serde(default)]
    pub asym_current_sensor: Vec<PgmAsymCurrentSensor>,
    /// Short-circuit fault locations. Read only by the short-circuit path
    /// (`crate::shortcircuit`); every other calculation ignores them, exactly
    /// as power-grid-model does.
    #[serde(default)]
    pub fault: Vec<PgmFault>,
}

#[derive(Clone, Deserialize)]
pub struct PgmNode {
    pub id: u64,
    pub u_rated: f64,
}

#[derive(Clone, Deserialize)]
pub struct PgmLine {
    pub id: u64,
    pub from_node: u64,
    pub to_node: u64,
    pub from_status: u8,
    pub to_status: u8,
    pub r1: f64,
    pub x1: f64,
    pub c1: f64,
    /// Positive-sequence shunt loss factor (tan δ). PGM forms the shunt
    /// admittance as `2πf·c1·(tan1 + j)` (`line.hpp`), so this is what gives a
    /// line's shunt a conductive part. Defaulted because not every fixture
    /// specifies it, and PGM treats an absent loss factor as zero.
    #[serde(default)]
    pub tan1: f64,
    /// Zero-sequence parameters, needed only for asymmetric calculations —
    /// power-grid-model marks them optional and many of its own symmetric
    /// fixtures omit them entirely.
    ///
    /// Defaulted to NaN rather than to zero or to the positive-sequence value:
    /// `r0 = x0 = 0` would make the zero-sequence admittance infinite, and
    /// `r0 = r1` is a modelling assumption that is wrong for real lines. NaN
    /// keeps the symmetric path (which never reads these) working while making
    /// any asymmetric use of an incomplete fixture loudly wrong instead of
    /// quietly plausible.
    #[serde(default = "nan")]
    pub r0: f64,
    #[serde(default = "nan")]
    pub x0: f64,
    #[serde(default = "nan")]
    pub c0: f64,
    /// Zero-sequence shunt loss factor. Parsed for completeness; the
    /// three-phase conversion (`pgm_to_3ph_network`) still models zero-sequence
    /// shunts as purely susceptive, since `Line3Ph` carries no conductance.
    #[serde(default)]
    pub tan0: f64,
}

/// power-grid-model's own `link` admittance, `1e8 + j1e8` per-unit — recorded
/// for reference, but **not** what gridoxide uses. See [`LINK_Y`].
///
/// Derived in its `common.hpp:83` as `1e6 / (base_power_3p / 10e3 / 10e3)` —
/// "1e6 siemens in a 10 kV network", evaluated once against a 1 MVA base and
/// frozen. It is a *per-unit* constant: the `10e3` is not a live voltage, so it
/// does not track the network's voltage level, and re-scaling it by voltage
/// would double-count what per-unit already removed. powsybl-open-loadflow
/// reaches the same order independently, clamping `|Z|` at `1e-8` p.u.
/// (`LfNetworkParameters.java:39`).
pub const PGM_LINK_Y: Complex<f64> = Complex::new(1e8, 1e8);

/// The admittance gridoxide stamps for a `link`.
///
/// An alias for [`topology::IDEAL_CONNECTION_Y`](crate::topology::IDEAL_CONNECTION_Y),
/// which documents the measurement behind the value and is shared with the
/// branches detected as ideal by impedance rather than declared as such — a
/// link and an undeclared jumper get the same treatment because they are the
/// same thing.
pub use crate::topology::IDEAL_CONNECTION_Y as LINK_Y;

fn one() -> f64 {
    1.0
}
/// A three-winding transformer's winding default — see `PgmThreeWindingTransformer`.
fn wye_n() -> u8 {
    crate::network::WYE_N
}
fn nan() -> f64 {
    f64::NAN
}
/// power-grid-model's documented `source.sk` default (`components.md`).
fn default_sk() -> f64 {
    1e10
}
/// power-grid-model's documented `source.rx_ratio` default (`components.md`).
fn default_rx_ratio() -> f64 {
    0.1
}
fn nan3() -> [f64; 3] {
    [f64::NAN; 3]
}
fn one_u8() -> u8 {
    1
}
/// power-grid-model's `FaultPhase::default_value`. Not a phase: a sentinel
/// meaning "whichever phases this fault type implies", resolved against the
/// fault type by `shortcircuit::FaultPhase::resolve`.
fn default_fault_phase() -> i8 {
    -1
}

/// Deserializes a number that power-grid-model may also write as the string
/// `"inf"`.
///
/// Its state-estimation fixtures use this for a measurement that is present but
/// carries no information (`inf-measurement-with-injection` and friends): an
/// infinite standard deviation is a zero weight, so the row exists structurally
/// and contributes nothing. JSON has no infinity literal, hence the string.
fn de_f64_or_inf<'de, D>(d: D) -> Result<f64, D::Error>
where
    D: serde::Deserializer<'de>,
{
    use serde::de::{Error, Unexpected};

    #[derive(Deserialize)]
    #[serde(untagged)]
    enum Raw {
        Num(f64),
        Str(String),
    }

    match Option::<Raw>::deserialize(d)? {
        None => Ok(f64::NAN),
        Some(Raw::Num(v)) => Ok(v),
        Some(Raw::Str(s)) => match s.trim() {
            "inf" | "+inf" | "Infinity" => Ok(f64::INFINITY),
            "-inf" | "-Infinity" => Ok(f64::NEG_INFINITY),
            "nan" | "NaN" => Ok(f64::NAN),
            other => Err(D::Error::invalid_value(
                Unexpected::Str(other),
                &"a number or the string \"inf\"",
            )),
        },
    }
}

#[derive(Clone, Deserialize)]
pub struct PgmSource {
    pub id: u64,
    pub node: u64,
    pub status: u8,
    /// Reference voltage, per-unit. Power-grid-model marks this required *only
    /// for power flow* — a state-estimation input has no reason to assert the
    /// source voltage, since that is precisely what is being estimated — so it
    /// falls back to 1.0 p.u. here.
    #[serde(default = "one")]
    pub u_ref: f64,
    /// Reference voltage angle, radians.
    ///
    /// **Read by the short-circuit path only.** `pgm_to_3ph_network` and
    /// `pgm_to_buses_and_branches` both build their virtual slack buses at
    /// angle zero (or, in the three-phase case, at the nominal phase angles)
    /// and ignore this field, which is what they have always done. No
    /// power-flow or state-estimation fixture in this tree sets it, so
    /// honouring it there would be an unvalidated behaviour change to those
    /// paths; the short-circuit fixtures *do* set it, which is why it is read
    /// here. Making the power-flow path honour it is a separate change that
    /// needs its own reference outputs.
    #[serde(default)]
    pub u_ref_angle: f64,
    /// Short-circuit power, VA. PGM's documented default is 1e10.
    #[serde(default = "default_sk")]
    pub sk: f64,
    /// R-to-X ratio. PGM's documented default is 0.1.
    #[serde(default = "default_rx_ratio")]
    pub rx_ratio: f64,
    #[serde(default = "one")]
    pub z01_ratio: f64,
}

/// PGM's ZIP-model type code: 0 = constant power, 1 = constant impedance
/// (S ∝ |V|²), 2 = constant current (S ∝ |V|).
fn default_load_type() -> u8 {
    0
}

#[derive(Clone, Deserialize)]
pub struct PgmSymLoad {
    pub id: u64,
    pub node: u64,
    pub status: u8,
    #[serde(rename = "type", default = "default_load_type")]
    pub load_type: u8,
    /// Specified power. power-grid-model marks this required *only for power
    /// flow*: a state-estimation input deliberately leaves it unset, because
    /// the appliance's power is what the estimator solves for. NaN when
    /// absent, and `pgm_to_network` contributes nothing to the bus injection
    /// for a non-finite value.
    #[serde(default = "nan")]
    pub p_specified: f64,
    #[serde(default = "nan")]
    pub q_specified: f64,
}

#[derive(Clone, Deserialize)]
pub struct PgmAsymLoad {
    pub id: u64,
    pub node: u64,
    pub status: u8,
    #[serde(rename = "type", default = "default_load_type")]
    pub load_type: u8,
    /// Per-phase specified power; see `PgmSymLoad::p_specified` for why this
    /// is optional.
    #[serde(default = "nan3")]
    pub p_specified: [f64; 3],
    #[serde(default = "nan3")]
    pub q_specified: [f64; 3],
}

#[derive(Clone, Deserialize)]
pub struct PgmSymGen {
    pub id: u64,
    pub node: u64,
    pub status: u8,
    #[serde(rename = "type", default = "default_load_type")]
    pub load_type: u8,
    /// Specified power. power-grid-model marks this required *only for power
    /// flow*: a state-estimation input deliberately leaves it unset, because
    /// the appliance's power is what the estimator solves for. NaN when
    /// absent, and `pgm_to_network` contributes nothing to the bus injection
    /// for a non-finite value.
    #[serde(default = "nan")]
    pub p_specified: f64,
    #[serde(default = "nan")]
    pub q_specified: f64,
}

/// PV-bus control: an active `voltage_regulator` pins its `regulated_object`
/// (a `sym_gen` id) bus's voltage magnitude to `u_ref`, letting Q float —
/// mirrors PGM's own `VoltageRegulator` component (`regulated_object` ==
/// `generator_id` in PGM's C++ `VoltageRegulator::calc_param()`).
///
/// `q_min`/`q_max` (VAr, PGM's own field names, optional — `NaN` if
/// omitted, matching PGM's own "unset" convention for this component) bound
/// the *bus's net* reactive injection, not any one generator's own gross
/// terminal output, since `Bus` (like PGM's own aggregated `q_specified`)
/// only tracks one netted P/Q per node — a real simplification versus a
/// bus with a PV generator *and* a significant co-located load, but exactly
/// consistent with how this parser already aggregates every other P/Q
/// quantity per bus. See [`outerloop::ReactiveLimits`](crate::outerloop::ReactiveLimits)
/// for where these get enforced (PV→PQ switching); plain `newton_raphson`
/// ignores them entirely, same as before.
#[derive(Clone, Deserialize)]
pub struct PgmVoltageRegulator {
    pub id: u64,
    pub regulated_object: u64,
    pub status: u8,
    pub u_ref: f64,
    #[serde(default = "nan")]
    pub q_min: f64,
    #[serde(default = "nan")]
    pub q_max: f64,
}

#[derive(Clone, Deserialize)]
pub struct PgmAsymGen {
    pub id: u64,
    pub node: u64,
    pub status: u8,
    #[serde(rename = "type", default = "default_load_type")]
    pub load_type: u8,
    /// Per-phase specified power; see `PgmSymLoad::p_specified` for why this
    /// is optional.
    #[serde(default = "nan3")]
    pub p_specified: [f64; 3],
    #[serde(default = "nan3")]
    pub q_specified: [f64; 3],
}

/// A short-circuit fault location, per power-grid-model's `fault` component
/// (`docs/user_manual/components.md`). A fault can only happen at a node.
///
/// **`r_f`/`x_f` default to zero, and zero means a *bolted* fault** — a dead
/// short, `y_fault = 1/0 = ∞` — not a fault that draws no current. That is
/// power-grid-model's documented default, it is what every one of its
/// `update_batch.json` scenarios drives, and it is the reading that makes the
/// component meaningful: a fault with no admittance would not be a fault.
/// `shortcircuit` handles the infinite-admittance case on its own code path
/// rather than letting an infinity into the matrix.
#[derive(Clone, Copy, Debug, Deserialize)]
pub struct PgmFault {
    pub id: u64,
    /// Whether the fault is active. Defaulted because power-grid-model's own
    /// batch scenarios omit it when only re-aiming an already-active fault.
    #[serde(default = "one_u8")]
    pub status: u8,
    /// `FaultType`: 0 three-phase, 1 single-phase-to-ground, 2 two-phase,
    /// 3 two-phase-to-ground. Interpreted by `shortcircuit::FaultType`.
    pub fault_type: i8,
    /// `FaultPhase`: 0 abc, 1 a, 2 b, 3 c, 4 ab, 5 ac, 6 bc, -1 default.
    #[serde(default = "default_fault_phase")]
    pub fault_phase: i8,
    /// The id of the faulted `node`.
    pub fault_object: u64,
    #[serde(default)]
    pub r_f: f64,
    #[serde(default)]
    pub x_f: f64,
}

#[derive(Clone, Deserialize)]
pub struct PgmShunt {
    pub id: u64,
    pub node: u64,
    pub status: u8,
    pub g1: f64,
    pub b1: f64,
    /// Zero-sequence admittance, needed only for asymmetric calculations and
    /// optional in power-grid-model. NaN when absent, for the same reason
    /// `PgmLine`'s zero-sequence fields are.
    #[serde(default = "nan")]
    pub g0: f64,
    #[serde(default = "nan")]
    pub b0: f64,
}

#[derive(Clone, Deserialize)]
pub struct PgmThreeWindingTransformer {
    pub id: u64,
    pub node_1: u64,
    pub node_2: u64,
    pub node_3: u64,
    pub status_1: u8,
    pub status_2: u8,
    pub status_3: u8,
    pub u1: f64,
    pub u2: f64,
    pub u3: f64,
    pub sn_1: f64,
    pub sn_2: f64,
    pub sn_3: f64,
    pub uk_12: f64,
    pub uk_13: f64,
    pub uk_23: f64,
    pub pk_12: f64,
    pub pk_13: f64,
    pub pk_23: f64,
    pub i0: f64,
    pub p0: f64,
    pub clock_12: i32,
    pub clock_13: i32,
    /// Winding configuration per side — `0` wye, `1` wye_n, `2` delta, `3`
    /// zigzag, `4` zigzag_n, as `network::WYE`..`ZIGZAG_N`.
    ///
    /// Read only by the phase domain, and unread for a long time because the
    /// symmetric star equivalent does not depend on them: the positive sequence
    /// is the same whatever the windings are, and the **zero** sequence is
    /// entirely determined by them. Defaulted to wye_n rather than required, so
    /// a document written for a symmetric calculation still parses; every
    /// committed three-winding fixture states all three.
    #[serde(default = "wye_n")]
    pub winding_1: u8,
    #[serde(default = "wye_n")]
    pub winding_2: u8,
    #[serde(default = "wye_n")]
    pub winding_3: u8,
    pub tap_side: u8,
    pub tap_pos: i32,
    pub tap_min: i32,
    pub tap_max: i32,
    pub tap_nom: i32,
    pub tap_size: f64,
}

#[derive(Deserialize, Clone)]
pub struct PgmTransformer {
    pub id: u64,
    pub from_node: u64,
    pub to_node: u64,
    pub from_status: u8,
    pub to_status: u8,
    pub u1: f64,
    pub u2: f64,
    pub sn: f64,
    pub uk: f64,
    pub pk: f64,
    pub i0: f64,
    pub p0: f64,
    #[serde(default)]
    pub winding_from: u8,
    #[serde(default)]
    pub winding_to: u8,
    pub clock: i32,
    pub tap_side: u8,
    pub tap_pos: i32,
    pub tap_min: i32,
    pub tap_max: i32,
    pub tap_nom: i32,
    pub tap_size: f64,
}

/// PGM's `sym_voltage_sensor`: a voltage measurement at a node.
///
/// Both `u_measured` and `u_angle_measured` are optional in practice — PGM's
/// own fixtures omit one or the other, and PGM reads an absent value as NaN
/// rather than as zero. They are defaulted to NaN here for the same reason, and
/// `measurement` skips any reading that isn't finite.
///
/// An angle-carrying voltage sensor is a phasor (PMU) measurement; one without
/// is an ordinary magnitude-only SCADA measurement.
#[derive(Clone, Deserialize)]
pub struct PgmSymVoltageSensor {
    pub id: u64,
    /// The measured `node`'s id.
    pub measured_object: u64,
    #[serde(default = "nan")]
    pub u_measured: f64,
    #[serde(default = "nan")]
    pub u_angle_measured: f64,
    /// Standard deviation of the magnitude error, in volts.
    #[serde(default = "nan")]
    pub u_sigma: f64,
    #[serde(default = "nan")]
    pub u_angle_sigma: f64,
}

/// PGM's `sym_power_sensor`: an active/reactive power measurement at one
/// terminal, of either a branch or an appliance.
///
/// `measured_terminal_type` is PGM's `MeasuredTerminalType`: 0 `branch_from`,
/// 1 `branch_to`, 2 `source`, 3 `shunt`, 4 `load`, 5 `generator`, 6/7/8
/// `branch3_1`/`_2`/`_3`, 9 `node`. It decides both what `measured_object`
/// refers to and the sign convention of the reading — see
/// [`measurement`](crate::measurement) for how each is mapped.
///
/// `power_sigma` applies to both components; `p_sigma`/`q_sigma` override it
/// per component when present (only two of power-grid-model's own fixtures use
/// them, but they are the more specific form).
#[derive(Clone, Deserialize)]
pub struct PgmSymPowerSensor {
    pub id: u64,
    pub measured_object: u64,
    pub measured_terminal_type: u8,
    #[serde(default = "nan", deserialize_with = "de_f64_or_inf")]
    pub p_measured: f64,
    #[serde(default = "nan", deserialize_with = "de_f64_or_inf")]
    pub q_measured: f64,
    #[serde(default = "nan", deserialize_with = "de_f64_or_inf")]
    pub power_sigma: f64,
    #[serde(default = "nan", deserialize_with = "de_f64_or_inf")]
    pub p_sigma: f64,
    #[serde(default = "nan", deserialize_with = "de_f64_or_inf")]
    pub q_sigma: f64,
}

/// PGM's `asym_voltage_sensor`: a per-phase voltage magnitude and angle.
///
/// The per-unit base differs from the symmetric sensor's. PGM divides by
/// `u_rated * u_scale`, with `u_scale` 1 for a symmetric sensor and `1/√3` for
/// an asymmetric one (`common.hpp`), because `u_measured` is line-to-line on the
/// first and line-to-neutral on the second. So an asymmetric reading is
/// `u_measured · √3 / u_rated` per unit.
///
/// `u_sigma` stays a single scalar covering all three phases; PGM has no
/// per-phase voltage sigma, and no `u_angle_sigma` at all (gridoxide's symmetric
/// struct carries one as an extension).
#[derive(Clone, Deserialize)]
pub struct PgmAsymVoltageSensor {
    pub id: u64,
    /// The measured `node`'s id.
    pub measured_object: u64,
    #[serde(default = "nan3")]
    pub u_measured: [f64; 3],
    #[serde(default = "nan3")]
    pub u_angle_measured: [f64; 3],
    /// Standard deviation of the magnitude error, in volts, line-to-neutral.
    #[serde(default = "nan", deserialize_with = "de_f64_or_inf")]
    pub u_sigma: f64,
}

/// PGM's `asym_power_sensor`: per-phase active/reactive power at one terminal.
///
/// Per-phase powers are per-unit against `s_base / 3`, PGM's `base_power_1p`,
/// where the symmetric sensor's are against the three-phase `s_base`.
///
/// Unlike the voltage sensor, the per-component sigmas here *are* per-phase
/// (`RealValue<sym>` in PGM's `PowerSensorInput`), while the fallback
/// `power_sigma` remains one scalar for the whole sensor.
#[derive(Clone, Deserialize)]
pub struct PgmAsymPowerSensor {
    pub id: u64,
    pub measured_object: u64,
    pub measured_terminal_type: u8,
    #[serde(default = "nan3")]
    pub p_measured: [f64; 3],
    #[serde(default = "nan3")]
    pub q_measured: [f64; 3],
    #[serde(default = "nan", deserialize_with = "de_f64_or_inf")]
    pub power_sigma: f64,
    #[serde(default = "nan3")]
    pub p_sigma: [f64; 3],
    #[serde(default = "nan3")]
    pub q_sigma: [f64; 3],
}

/// PGM's `sym_current_sensor`: a current magnitude and angle at a branch
/// terminal.
///
/// `angle_measurement_type` is 0 for a *local* angle — the shift between the
/// terminal's voltage and its current, what a power-factor meter reads — and 1
/// for a *global* one, measured against the same reference voltage phasors use.
/// The two are different quantities, not two spellings of one; see
/// [`AngleFrame`](crate::measurement::AngleFrame).
///
/// Only branch terminals carry one: `measured_terminal_type` 0/1 and 6/7/8, and
/// never a `link`. Both sigmas are scalar, including on the asymmetric variant.
#[derive(Clone, Deserialize)]
pub struct PgmSymCurrentSensor {
    pub id: u64,
    pub measured_object: u64,
    pub measured_terminal_type: u8,
    #[serde(default)]
    pub angle_measurement_type: u8,
    #[serde(default = "nan", deserialize_with = "de_f64_or_inf")]
    pub i_measured: f64,
    #[serde(default = "nan", deserialize_with = "de_f64_or_inf")]
    pub i_angle_measured: f64,
    /// Standard deviation of the magnitude error, in amperes.
    #[serde(default = "nan", deserialize_with = "de_f64_or_inf")]
    pub i_sigma: f64,
    /// Standard deviation of the angle error, in radians.
    #[serde(default = "nan", deserialize_with = "de_f64_or_inf")]
    pub i_angle_sigma: f64,
}

/// PGM's `asym_current_sensor`: per-phase magnitudes and angles, scalar sigmas.
#[derive(Clone, Deserialize)]
pub struct PgmAsymCurrentSensor {
    pub id: u64,
    pub measured_object: u64,
    pub measured_terminal_type: u8,
    #[serde(default)]
    pub angle_measurement_type: u8,
    #[serde(default = "nan3")]
    pub i_measured: [f64; 3],
    #[serde(default = "nan3")]
    pub i_angle_measured: [f64; 3],
    #[serde(default = "nan", deserialize_with = "de_f64_or_inf")]
    pub i_sigma: f64,
    #[serde(default = "nan", deserialize_with = "de_f64_or_inf")]
    pub i_angle_sigma: f64,
}

// ── Output structs (used by integration tests) ────────────────────────────────

#[derive(Clone, Deserialize)]
pub struct PgmOutput<N> {
    pub data: PgmOutputData<N>,
}

#[derive(Clone, Deserialize)]
pub struct PgmOutputData<N> {
    pub node: Vec<N>,
}

/// A batch output document: one `PgmOutputData<N>` per scenario, in scenario order.
#[derive(Clone, Deserialize)]
pub struct PgmBatchOutput<N> {
    pub data: Vec<PgmOutputData<N>>,
}

#[derive(Clone, Deserialize)]
pub struct PgmNodeOutput {
    pub id: u64,
    pub u_pu: f64,
    pub u_angle: f64,
    #[serde(default)]
    pub u: f64,
}

/// PGM's per-branch output record, as found in `sym_output.json`'s `line`,
/// `transformer` and `link` arrays.
///
/// Powers are in W/var and currents in A — physical units, not per-unit, so
/// comparing against gridoxide's [`branch_flow`](crate::branch_flow) results
/// means scaling those by `s_base_va`. Only the flow fields are modelled;
/// `loading`/`i_from`/`s_from` are derived quantities gridoxide does not
/// compute.
#[derive(Debug, Deserialize)]
pub struct PgmLineOutput {
    pub id: u64,
    #[serde(default)]
    pub energized: u8,
    pub p_from: f64,
    pub q_from: f64,
    pub p_to: f64,
    pub q_to: f64,
}

#[derive(Debug, Deserialize)]
pub struct PgmNodeAsymOutput {
    pub id: u64,
    pub u_pu: [f64; 3],
    pub u_angle: [f64; 3],
}

// ── Short-circuit output structs ──────────────────────────────────────────────
//
// Short-circuit output is always asymmetric, whatever the fault type — so
// every quantity here is a three-element array, unlike the symmetric
// power-flow records above. Each component list is `#[serde(default)]`
// because power-grid-model's fixtures omit the ones a given network has none
// of.

#[derive(Debug, Deserialize)]
pub struct PgmScOutput {
    pub data: PgmScOutputData,
}

/// A batch short-circuit output: one entry per scenario, in scenario order.
#[derive(Debug, Deserialize)]
pub struct PgmScBatchOutput {
    pub data: Vec<PgmScOutputData>,
}

#[derive(Debug, Default, Deserialize)]
pub struct PgmScOutputData {
    #[serde(default)]
    pub node: Vec<PgmScNodeOutput>,
    #[serde(default)]
    pub fault: Vec<PgmScFaultOutput>,
    #[serde(default)]
    pub source: Vec<PgmScApplianceOutput>,
    #[serde(default)]
    pub line: Vec<PgmScBranchOutput>,
    #[serde(default)]
    pub transformer: Vec<PgmScBranchOutput>,
    #[serde(default)]
    pub link: Vec<PgmScBranchOutput>,
    #[serde(default)]
    pub shunt: Vec<PgmScApplianceOutput>,
}

// Every quantity is optional: power-grid-model's fixtures assert only the
// attributes each one is about, so a fixture may carry `u_pu` and no angle at
// all. `Option` rather than a zero default, so a test can tell "not asserted"
// from "asserted to be zero" — the two are very different when a bolted fault
// legitimately drives a voltage to zero.

#[derive(Debug, Deserialize)]
pub struct PgmScNodeOutput {
    pub id: u64,
    #[serde(default)]
    pub energized: Option<u8>,
    #[serde(default)]
    pub u_pu: Option<[f64; 3]>,
    #[serde(default)]
    pub u_angle: Option<[f64; 3]>,
    #[serde(default)]
    pub u: Option<[f64; 3]>,
}

#[derive(Debug, Deserialize)]
pub struct PgmScFaultOutput {
    pub id: u64,
    #[serde(default)]
    pub energized: Option<u8>,
    #[serde(default)]
    pub i_f: Option<[f64; 3]>,
    #[serde(default)]
    pub i_f_angle: Option<[f64; 3]>,
}

/// A `source` or `shunt` short-circuit record — a one-terminal appliance.
#[derive(Debug, Deserialize)]
pub struct PgmScApplianceOutput {
    pub id: u64,
    #[serde(default)]
    pub energized: Option<u8>,
    #[serde(default)]
    pub i: Option<[f64; 3]>,
    #[serde(default)]
    pub i_angle: Option<[f64; 3]>,
}

#[derive(Debug, Deserialize)]
pub struct PgmScBranchOutput {
    pub id: u64,
    #[serde(default)]
    pub energized: Option<u8>,
    #[serde(default)]
    pub i_from: Option<[f64; 3]>,
    #[serde(default)]
    pub i_from_angle: Option<[f64; 3]>,
    #[serde(default)]
    pub i_to: Option<[f64; 3]>,
    #[serde(default)]
    pub i_to_angle: Option<[f64; 3]>,
}

// ── Public helpers ────────────────────────────────────────────────────────────

/// power-grid-model's `link`: an ideal connection between two nodes, carrying
/// no attributes beyond its endpoints and their statuses.
///
/// Modelled as a branch rather than by merging its endpoints, because a `link`
/// has an identity in power-grid-model's output model — the schema carries a
/// `link` record with its own `p`/`q`/`i` flows, and its fixtures assert them.
/// Merging would delete the branch those numbers describe. See
/// [`crate::topology`] for the policy and
/// `docs/src/powerflow/zero_impedance_branches.md` for the alternatives.
#[derive(Clone, Deserialize)]
pub struct PgmLink {
    pub id: u64,
    pub from_node: u64,
    pub to_node: u64,
    pub from_status: u8,
    pub to_status: u8,
}

/// A line's total per-unit shunt admittance, `2πf·c1·(tan1 + j)·z_base`.
///
/// Matches power-grid-model's `line.hpp`:
/// `y1_shunt = 2π·f·c1/base_y·(tan1 + 1i)`. The conductive part comes entirely
/// from the loss factor, so it is zero for the many fixtures that leave `tan1`
/// unset — but not for the state-estimation fixtures, which specify it.
fn line_y_shunt(ln: &PgmLine, omega: f64, z_base: f64) -> Complex<f64> {
    Complex::new(ln.tan1, 1.0) * (omega * ln.c1 * z_base)
}

/// Builds the `Line` representing a branch with exactly one terminal connected.
///
/// PGM does not simply drop the series impedance here: an open-ended branch
/// still presents the near-end shunt half *in parallel with* the series
/// impedance feeding the far-end shunt half, i.e.
/// `y_sh/2 + 1/(1/y_s + 2/y_sh)` — [`half_open_branch_shunt`], the same
/// function [`branch_calc_param`] already applies to half-open transformers.
///
/// Modelling this as the bare shunt instead (which gridoxide did until the
/// branch-flow work) omits the small conductive path through the series
/// resistance: `tests/branch_flow_test.rs`'s `line` fixture has two half-open
/// lines, each carrying 0.68 W in PGM's own expected output, and their absence
/// showed up as a 1.36 W deficit on the healthy line feeding them.
///
/// The result is a self-loop `Line`, which `build_ybus` stamps as a pure
/// diagonal admittance. That representation loses which terminal was the
/// connected one — see [`branch_flow::line_params`](crate::branch_flow::line_params).
fn half_open_line(ln: &PgmLine, omega: f64, z_base: f64, idx: usize) -> Line {
    let y_shunt = line_y_shunt(ln, omega, z_base);
    let (r, x) = crate::topology::clamp_branch_impedance(ln.r1 / z_base, ln.x1 / z_base);
    let y_series = Complex::new(1.0, 0.0) / Complex::new(r, x);
    let y_eq = half_open_branch_shunt(y_series, y_shunt);
    Line {
        from: idx,
        to: idx,
        r: 0.0,
        x: 0.0,
        b_shunt: y_eq.im,
        g_shunt: y_eq.re,
    }
}

/// Returns a stable node-ID → 0-based-index map (sorted by node ID).
pub fn node_id_to_idx(input: &PgmInput) -> HashMap<u64, usize> {
    let mut ids: Vec<u64> = input.data.node.iter().map(|n| n.id).collect();
    ids.sort_unstable();
    ids.into_iter()
        .enumerate()
        .map(|(idx, id)| (id, idx))
        .collect()
}

/// Converts active `shunt` entries to per-unit self-admittances, for stamping
/// into a 1-phase Y-bus via `network::stamp_shunts`. Call before the owning
/// `PgmInput` is consumed by `pgm_to_buses_and_branches`.
pub fn pgm_shunts_1ph(
    input: &PgmInput,
    id_to_idx: &HashMap<u64, usize>,
    s_base_va: f64,
) -> Vec<ShuntAdm> {
    let id_to_u_rated: HashMap<u64, f64> =
        input.data.node.iter().map(|n| (n.id, n.u_rated)).collect();
    input
        .data
        .shunt
        .iter()
        .filter(|s| s.status != 0)
        .map(|s| {
            let u_rated = id_to_u_rated[&s.node];
            let base_y = s_base_va / (u_rated * u_rated);
            ShuntAdm {
                at: id_to_idx[&s.node],
                y: Complex::new(s.g1, s.b1) / base_y,
            }
        })
        .collect()
}

/// Converts active `shunt` entries to per-unit sequence self-admittances, for
/// stamping into a 3-phase Y-bus via `network::stamp_shunts_3ph`. `at` is the
/// physical node index (i.e. `id_to_idx` from `pgm_to_3ph_network`, not the
/// 3N-bus index — `stamp_shunts_3ph` multiplies by 3 internally).
pub fn pgm_shunts_3ph(
    input: &PgmInput,
    id_to_idx: &HashMap<u64, usize>,
    s_base_va: f64,
) -> Vec<ShuntAdm3Ph> {
    let id_to_u_rated: HashMap<u64, f64> =
        input.data.node.iter().map(|n| (n.id, n.u_rated)).collect();
    input
        .data
        .shunt
        .iter()
        .filter(|s| s.status != 0)
        .map(|s| {
            let u_rated = id_to_u_rated[&s.node];
            let base_y = s_base_va / (u_rated * u_rated);
            ShuntAdm3Ph {
                at: id_to_idx[&s.node],
                y1: Complex::new(s.g1, s.b1) / base_y,
                y0: Complex::new(s.g0, s.b0) / base_y,
            }
        })
        .collect()
}

/// Converts active `transformer` entries to sequence-domain branch admittance
/// parameters, for stamping into a 3-phase Y-bus via
/// `network::stamp_transformers_3ph`. Call before the owning `PgmInput` is
/// consumed by `pgm_to_3ph_network`. Only Dyn (`winding_from`=delta,
/// `winding_to`=wye_n) and YNyn (`winding_from`=`winding_to`=wye_n)
/// transformers are supported — see `network::transformer_seq_params`.
/// The three-phase sequence parameters of a `link`.
///
/// An ideal connection looks the same to all three sequences — a plain series
/// admittance, no shunt, no tap, no phase shift — so the three sequence
/// parameters are one `branch_calc_param` repeated. power-grid-model's own
/// `link.hpp` agrees literally: its asymmetric parameter is
/// `calc_param_y_asym(y_link, 0.0, y_link, 0.0, 1.0)`, the same admittance in
/// the positive and zero sequences with no shunt in either.
///
/// `None` for a link with both terminals open, which is not a branch at all.
pub fn link_seq_params(lk: &PgmLink, id_to_idx: &HashMap<u64, usize>) -> Option<Transformer3PhSeq> {
    if lk.from_status == 0 && lk.to_status == 0 {
        return None;
    }
    let y = branch_calc_param(
        LINK_Y,
        Complex::new(0.0, 0.0),
        Complex::new(1.0, 0.0),
        lk.from_status,
        lk.to_status,
    );
    Some(Transformer3PhSeq {
        from: id_to_idx[&lk.from_node],
        to: id_to_idx[&lk.to_node],
        y0: y,
        y1: y,
        y2: y,
    })
}

/// One leg of a three-winding transformer's star equivalent.
///
/// power-grid-model models a three-winding transformer as three two-winding
/// transformers to a shared internal star node
/// (`three_winding_transformer.hpp::convert_to_two_winding_transformers`), and
/// so does gridoxide. What the symmetric path never had to carry is each leg's
/// **winding configuration** — the zero sequence depends on it and the positive
/// sequence does not, so a balanced case cannot tell a right one from a wrong
/// one.
///
/// Assembled here once so the symmetric conversion and the phase-domain one
/// cannot disagree about what the legs are.
pub struct ThreeWindingLeg {
    /// The physical node this leg reaches, as a document id.
    pub node: u64,
    pub status: u8,
    pub y_series: Complex<f64>,
    pub y_shunt: Complex<f64>,
    pub tap: Complex<f64>,
    /// The winding on the *physical* side.
    pub winding_from: u8,
    /// The winding on the *star* side.
    pub winding_to: u8,
    pub sn: f64,
    pub uk: f64,
    pub clock: i32,
}

/// The three legs of a three-winding transformer's star equivalent.
///
/// Follows power-grid-model exactly, including the two asymmetries a reader
/// would not guess: T1 is `wye_n`/`wye_n` at clock 0 whatever `winding_1` is —
/// it is the leg the star node's base is *defined* by — while T2 and T3 carry
/// their own winding against `winding_1` and a **reversed** clock, `12 − clock`.
/// Only T1 carries the magnetizing branch; T2 and T3 get `i0 = p0 = 0`, so the
/// no-load loss is not counted three times.
pub fn three_winding_legs(
    t: &PgmThreeWindingTransformer,
    id_to_u_rated: &HashMap<u64, f64>,
    s_base_va: f64,
) -> [ThreeWindingLeg; 3] {
    let u1_rated = id_to_u_rated[&t.node_1];
    let u2_rated = id_to_u_rated[&t.node_2];
    let u3_rated = id_to_u_rated[&t.node_3];

    // Tap adjustment applies to exactly one of u1/u2/u3, matching `tap_side`.
    let tap_direction = if t.tap_max > t.tap_min { 1.0 } else { -1.0 };
    let delta = tap_direction * (t.tap_pos - t.tap_nom) as f64 * t.tap_size;
    let (u1_local, u2_local, u3_local) = match t.tap_side {
        0 => (t.u1 + delta, t.u2, t.u3),
        1 => (t.u1, t.u2 + delta, t.u3),
        _ => (t.u1, t.u2, t.u3 + delta),
    };

    let ((uk_t1, uk_t2, uk_t3), (pk_t1, pk_t2, pk_t3)) = three_winding_star_params(
        t.sn_1, t.sn_2, t.sn_3, t.uk_12, t.uk_13, t.uk_23, t.pk_12, t.pk_13, t.pk_23,
    );

    // All three legs' star-side nameplate is `u1_local` and their per-unit base
    // is pinned to `u1_rated` (physical, fixed), which is why
    // `transformer_admittances_ex` always takes `(u1_local, u1_rated)` whichever
    // leg is being built.
    let (t1_series, t1_shunt) = transformer_admittances_ex(
        u1_local, u1_rated, t.sn_1, uk_t1, pk_t1, t.i0, t.p0, s_base_va,
    );
    let (t2_series, t2_shunt) = transformer_admittances_ex(
        u1_local, u1_rated, t.sn_2, uk_t2, pk_t2, 0.0, 0.0, s_base_va,
    );
    let (t3_series, t3_shunt) = transformer_admittances_ex(
        u1_local, u1_rated, t.sn_3, uk_t3, pk_t3, 0.0, 0.0, s_base_va,
    );

    [
        ThreeWindingLeg {
            node: t.node_1,
            status: t.status_1,
            y_series: t1_series,
            y_shunt: t1_shunt,
            tap: tap_ratio_from_voltages(u1_local * u1_rated, u1_local * u1_rated, 0),
            winding_from: crate::network::WYE_N,
            winding_to: crate::network::WYE_N,
            sn: t.sn_1,
            uk: uk_t1,
            clock: 0,
        },
        ThreeWindingLeg {
            node: t.node_2,
            status: t.status_2,
            y_series: t2_series,
            y_shunt: t2_shunt,
            tap: tap_ratio_from_voltages(u2_local * u1_rated, u1_local * u2_rated, 12 - t.clock_12),
            winding_from: t.winding_2,
            winding_to: t.winding_1,
            sn: t.sn_2,
            uk: uk_t2,
            clock: 12 - t.clock_12,
        },
        ThreeWindingLeg {
            node: t.node_3,
            status: t.status_3,
            y_series: t3_series,
            y_shunt: t3_shunt,
            tap: tap_ratio_from_voltages(u3_local * u1_rated, u1_local * u3_rated, 12 - t.clock_13),
            winding_from: t.winding_3,
            winding_to: t.winding_1,
            sn: t.sn_3,
            uk: uk_t3,
            clock: 12 - t.clock_13,
        },
    ]
}

pub fn pgm_transformers_3ph(
    input: &PgmInput,
    id_to_idx: &HashMap<u64, usize>,
    s_base_va: f64,
) -> Vec<Transformer3PhSeq> {
    let id_to_u_rated: HashMap<u64, f64> =
        input.data.node.iter().map(|n| (n.id, n.u_rated)).collect();
    let n_nodes = input.data.node.len();
    input
        .data
        .transformer
        .iter()
        .map(|t| {
            let tap = transformer_tap(
                t.u1,
                t.u2,
                t.tap_side,
                t.tap_pos,
                t.tap_min,
                t.tap_max,
                t.tap_nom,
                t.tap_size,
                t.clock,
                id_to_u_rated[&t.from_node],
                id_to_u_rated[&t.to_node],
            );
            let (y_series, y_shunt) =
                transformer_admittances(t.u2, t.sn, t.uk, t.pk, t.i0, t.p0, s_base_va);
            let (y0, y1, y2) = transformer_seq_params(
                y_series,
                y_shunt,
                tap,
                t.from_status,
                t.to_status,
                t.winding_from,
                t.winding_to,
                t.sn,
                t.uk,
                s_base_va,
                t.clock,
            );
            Transformer3PhSeq {
                from: id_to_idx[&t.from_node],
                to: id_to_idx[&t.to_node],
                y0,
                y1,
                y2,
            }
        })
        // Links, after the transformers — the position `pgm_3ph_maps` assigns
        // them, and the one the short-circuit path already used. A link is a
        // branch like any other here; dropping it leaves whatever hangs off it
        // de-energized and absent from the answer.
        .chain(
            input
                .data
                .link
                .iter()
                .filter_map(|lk| link_seq_params(lk, id_to_idx)),
        )
        // Then the three-winding transformers, three legs each, to the star
        // node `pgm_to_3ph_network` appended for them. Same order as the
        // symmetric conversion, so the two agree about what a branch index is.
        .chain(
            input
                .data
                .three_winding_transformer
                .iter()
                .enumerate()
                .flat_map(|(i, t)| {
                    let star = n_nodes + i;
                    three_winding_legs(t, &id_to_u_rated, s_base_va)
                        .into_iter()
                        .map(move |leg| {
                            // The star side is always in service — it is an internal
                            // node, not a terminal anyone can open — so only the
                            // physical side's status varies, exactly as in
                            // power-grid-model's own two-winding equivalents.
                            let (y0, y1, y2) = crate::network::transformer_seq_params(
                                leg.y_series,
                                leg.y_shunt,
                                leg.tap,
                                leg.status,
                                1,
                                leg.winding_from,
                                leg.winding_to,
                                leg.sn,
                                leg.uk,
                                s_base_va,
                                leg.clock,
                            );
                            Transformer3PhSeq {
                                from: id_to_idx[&leg.node],
                                to: star,
                                y0,
                                y1,
                                y2,
                            }
                        })
                        .collect::<Vec<_>>()
                }),
        )
        .collect()
}

/// Converts a PGM input document to per-unit buses, lines, and transformers.
///
/// `s_base_va` — system base power in VA (e.g. 1e6 for 1 MVA).
/// `freq_hz`   — grid frequency in Hz (e.g. 50.0 or 60.0).
///
/// All PGM nodes become PQ buses. Each active source is modelled as a virtual
/// Slack bus appended after the physical nodes, connected via a source-impedance
/// `Line`. Transformers are returned as `Transformer` values in per-unit on the
/// system base; stamp them into a Y-bus with `network::stamp_transformers`.
/// A converted PGM network together with the object-ID maps needed to address
/// its pieces by their original PGM ids.
///
/// [`pgm_to_buses_and_branches`] returns only the three vectors, which is all
/// power flow ever needed: it addresses buses positionally and never has to name
/// a branch. Measurements do — a `sym_power_sensor` carries
/// `measured_object: 4`, a PGM object id — and the mapping from that id to a
/// branch index is *not* recoverable by arithmetic on input positions:
///
/// - a line with both terminals open is dropped entirely (no branch at all),
/// - each active `source` appends a virtual slack bus and an extra branch,
/// - each `three_winding_transformer` expands into three branches plus a star
///   bus,
///
/// so input position and branch index diverge as soon as any of those appear.
/// Hence these maps are recorded during conversion rather than reconstructed.
///
/// All branch indices are *flat*: lines first, then transformers, matching the
/// order [`branch_flow::branch_params`](crate::branch_flow::branch_params)
/// assembles them in.
#[derive(Clone, Debug)]
pub struct PgmNetwork {
    pub buses: Vec<Bus>,
    pub lines: Vec<Line>,
    pub transformers: Vec<Transformer>,
    /// Node id → bus index. Covers physical nodes only; star buses and virtual
    /// slack buses have no PGM id of their own.
    pub node_idx: HashMap<u64, usize>,
    /// `line`/`transformer` id → flat branch index. A line with both terminals
    /// open is absent, since it produces no branch.
    pub branch_idx: HashMap<u64, usize>,
    /// `three_winding_transformer` id → its three flat branch indices, in
    /// side-1/2/3 order.
    pub three_winding_branch_idx: HashMap<u64, [usize; 3]>,
    /// `source` id → the flat branch index of the virtual source-impedance
    /// branch gridoxide synthesizes for it. Inactive sources are absent.
    pub source_branch_idx: HashMap<u64, usize>,
    /// Appliance id (`sym_load`, `asym_load`, `sym_gen`, `asym_gen`, `shunt`,
    /// `source`) → the bus it is attached to. Includes inactive appliances: a
    /// sensor may reference one, and reporting a zero flow for a disconnected
    /// appliance is better than failing to resolve the id at all.
    pub appliance_bus: HashMap<u64, usize>,
    /// Buses whose net injection is identically zero, as a per-bus flag.
    ///
    /// A bus with no load and no generator injects nothing into the network,
    /// and that is *structural knowledge*, not a measurement — it holds exactly,
    /// with no uncertainty. State estimation can use it as a hard constraint
    /// rather than as a very-high-weight pseudo-measurement, which is what
    /// `se::constraints` does.
    ///
    /// Sources and shunts do not disqualify a bus here, because gridoxide
    /// models both structurally: a source's power arrives through its
    /// synthesized branch and a shunt sits on the Y-bus diagonal, so neither
    /// appears in `network::power_injections` at that bus. The virtual slack
    /// buses themselves *are* excluded — that is precisely where the source's
    /// unknown power enters the network.
    ///
    /// Note this is read off the *input document*, not off `Bus::p_spec`. A
    /// state-estimation document leaves `p_specified` unset, so an unmeasured
    /// load looks like zero injection in the converted network while being
    /// nothing of the sort.
    pub zero_injection: Vec<bool>,
    /// Branch ids that collapsed to a self-loop because exactly one terminal
    /// was connected, mapped to *which* PGM terminal is the live one.
    ///
    /// gridoxide represents such a branch as a single diagonal admittance
    /// (`from == to`), which by construction cannot distinguish its two ends,
    /// so the distinction is kept here instead. Use [`resolve_terminal`] rather
    /// than reading this directly.
    ///
    /// [`resolve_terminal`]: PgmNetwork::resolve_terminal
    pub half_open_terminal: HashMap<u64, Terminal>,
}

impl PgmNetwork {
    /// Maps a PGM branch id and the terminal a sensor measures to the branch
    /// index and terminal to evaluate with
    /// [`branch_flow::terminal_flow`](crate::branch_flow::terminal_flow).
    ///
    /// Returns `None` in the two cases where there is no flow to compute:
    ///
    /// - the branch has both terminals open, so no branch exists at all;
    /// - the requested terminal is the *open* end of a half-open branch, whose
    ///   flow is identically zero (which is what PGM reports for it too).
    ///
    /// For the live end of a half-open branch the returned terminal is always
    /// [`Terminal::From`], since the whole equivalent admittance sits there
    /// regardless of which PGM end was connected.
    pub fn resolve_terminal(&self, id: u64, terminal: Terminal) -> Option<(usize, Terminal)> {
        let &branch = self.branch_idx.get(&id)?;
        match self.half_open_terminal.get(&id) {
            Some(&live) if live == terminal => Some((branch, Terminal::From)),
            Some(_) => None,
            None => Some((branch, terminal)),
        }
    }

    /// The state-estimation model of this network, given its Y-bus and the
    /// shunts stamped into it — see [`SeNetwork::new`](crate::se::SeNetwork::new).
    pub fn se_network(&self, ybus: YBusSparse, shunts: &[ShuntAdm]) -> crate::se::SeNetwork {
        crate::se::SeNetwork::new(
            ybus,
            &self.lines,
            &self.transformers,
            shunts,
            &self.source_branch_idx,
            &self.zero_injection,
        )
    }
}

/// Converts a PGM input document into gridoxide's own network types, keeping
/// the object-ID maps — see [`PgmNetwork`] for why those cannot be derived
/// afterwards.
///
/// [`pgm_to_buses_and_branches`] is the same conversion with the maps dropped.
pub fn pgm_to_network(input: PgmInput, s_base_va: f64, freq_hz: f64) -> PgmNetwork {
    let id_to_idx = node_id_to_idx(&input);
    let id_to_u_rated: HashMap<u64, f64> =
        input.data.node.iter().map(|n| (n.id, n.u_rated)).collect();

    // Accumulate per-node net injection from active loads (load = negative
    // injection) and gens (generation = positive injection). Constant-power
    // (type 0) entries fold into the flat p_inj/q_inj maps; constant-current/
    // -impedance (type 1/2) entries become voltage-dependent ZIP terms on Bus.
    let mut p_inj: HashMap<u64, f64> = HashMap::new();
    let mut q_inj: HashMap<u64, f64> = HashMap::new();
    let mut zip_map: HashMap<u64, Vec<ZipTerm>> = HashMap::new();
    let mut accumulate = |node: u64, load_type: u8, p: f64, q: f64, sign: f64| {
        // An appliance with no specified power contributes nothing. PGM marks
        // `p_specified`/`q_specified` required only for power flow, and a
        // state-estimation input leaves them unset precisely because the
        // appliance's power is an unknown — letting NaN through here would
        // poison the bus injection and, from there, the whole Y-bus solve.
        if !p.is_finite() && !q.is_finite() {
            return;
        }
        let p = if p.is_finite() { p } else { 0.0 };
        let q = if q.is_finite() { q } else { 0.0 };
        let s = Complex::new(p, q) * sign / s_base_va;
        match load_type {
            1 => zip_map.entry(node).or_default().push(ZipTerm {
                s_const: s,
                kind: ZipKind::ConstImpedance,
            }),
            2 => zip_map.entry(node).or_default().push(ZipTerm {
                s_const: s,
                kind: ZipKind::ConstCurrent,
            }),
            _ => {
                *p_inj.entry(node).or_insert(0.0) += s.re;
                *q_inj.entry(node).or_insert(0.0) += s.im;
            }
        }
    };
    for load in &input.data.sym_load {
        if load.status == 0 {
            continue;
        }
        accumulate(
            load.node,
            load.load_type,
            load.p_specified,
            load.q_specified,
            -1.0,
        );
    }
    for sgen in &input.data.sym_gen {
        if sgen.status == 0 {
            continue;
        }
        accumulate(
            sgen.node,
            sgen.load_type,
            sgen.p_specified,
            sgen.q_specified,
            1.0,
        );
    }
    // Asymmetric (3-phase) loads/gens are folded into their total three-phase
    // power for the symmetric (positive-sequence) equivalent.
    for load in &input.data.asym_load {
        if load.status == 0 {
            continue;
        }
        let p: f64 = load.p_specified.iter().sum();
        let q: f64 = load.q_specified.iter().sum();
        accumulate(load.node, load.load_type, p, q, -1.0);
    }
    for agen in &input.data.asym_gen {
        if agen.status == 0 {
            continue;
        }
        let p: f64 = agen.p_specified.iter().sum();
        let q: f64 = agen.q_specified.iter().sum();
        accumulate(agen.node, agen.load_type, p, q, 1.0);
    }

    // Physical PQ buses. Nodes with no path to any active source (PGM calls
    // these "de-energized") aren't special-cased here — every active source
    // already becomes its own real `Slack` bus wired in via a genuine Line
    // (below), so `network::connected_components`/`classify` downstream
    // (run inside the solver entry points themselves) correctly discovers
    // such a node's component has no reference bus and pins it to a V=0
    // placeholder without this function needing to duplicate that graph walk.
    let n_nodes = input.data.node.len();
    let mut sorted_ids: Vec<u64> = input.data.node.iter().map(|n| n.id).collect();
    sorted_ids.sort_unstable();
    let mut opt_buses = vec![None::<Bus>; n_nodes];
    for id in &sorted_ids {
        let idx = id_to_idx[id];
        opt_buses[idx] = Some(Bus {
            idx,
            bus_type: BusType::PQ,
            voltage_mag: 1.0,
            voltage_ang: 0.0,
            p_spec: *p_inj.get(id).unwrap_or(&0.0),
            q_spec: *q_inj.get(id).unwrap_or(&0.0),
            q_min: -f64::INFINITY,
            q_max: f64::INFINITY,
            u_rated: id_to_u_rated[id],
            zip_terms: zip_map.remove(id).unwrap_or_default(),
        });
    }
    let mut buses: Vec<Bus> = opt_buses.into_iter().map(|b| b.unwrap()).collect();

    // PV buses: an active `voltage_regulator` pins its regulated `sym_gen`'s
    // bus voltage magnitude to `u_ref`, letting Q float (see `PgmVoltageRegulator`).
    let sym_gen_node: HashMap<u64, u64> =
        input.data.sym_gen.iter().map(|g| (g.id, g.node)).collect();
    for vr in &input.data.voltage_regulator {
        if vr.status == 0 {
            continue;
        }
        let Some(&node) = sym_gen_node.get(&vr.regulated_object) else {
            continue;
        };
        let idx = id_to_idx[&node];
        if buses[idx].bus_type == BusType::PQ {
            buses[idx].bus_type = BusType::PV;
            buses[idx].voltage_mag = vr.u_ref;
            buses[idx].q_min = if vr.q_min.is_nan() {
                -f64::INFINITY
            } else {
                vr.q_min / s_base_va
            };
            buses[idx].q_max = if vr.q_max.is_nan() {
                f64::INFINITY
            } else {
                vr.q_max / s_base_va
            };
        }
    }

    // Lines. PGM's c1 is the *total* shunt capacitance; build_ybus splits b_shunt/2
    // per end, matching PGM's y_shunt/2. Half-open cases become self-loop shunts.
    let omega = 2.0 * std::f64::consts::PI * freq_hz;
    let mut lines: Vec<Line> = Vec::new();
    // Recorded as we go; see `PgmNetwork` for why these can't be rebuilt later.
    let mut branch_idx: HashMap<u64, usize> = HashMap::new();
    let mut three_winding_branch_idx: HashMap<u64, [usize; 3]> = HashMap::new();
    let mut source_branch_idx: HashMap<u64, usize> = HashMap::new();
    let mut half_open_terminal: HashMap<u64, Terminal> = HashMap::new();
    for ln in &input.data.line {
        match (ln.from_status, ln.to_status) {
            (1, 1) => {
                let z_base = id_to_u_rated[&ln.from_node].powi(2) / s_base_va;
                let y_shunt = line_y_shunt(ln, omega, z_base);
                // A line short enough to be a jumper would otherwise put an
                // unbounded admittance into the Y-bus; see
                // `topology::ZERO_IMPEDANCE_THRESHOLD`.
                let (r, x) =
                    crate::topology::clamp_branch_impedance(ln.r1 / z_base, ln.x1 / z_base);
                branch_idx.insert(ln.id, lines.len());
                lines.push(Line {
                    from: id_to_idx[&ln.from_node],
                    to: id_to_idx[&ln.to_node],
                    r,
                    x,
                    b_shunt: y_shunt.im,
                    g_shunt: y_shunt.re,
                });
            }
            (1, 0) => {
                let z_base = id_to_u_rated[&ln.from_node].powi(2) / s_base_va;
                let idx = id_to_idx[&ln.from_node];
                branch_idx.insert(ln.id, lines.len());
                half_open_terminal.insert(ln.id, Terminal::From);
                lines.push(half_open_line(ln, omega, z_base, idx));
            }
            (0, 1) => {
                let z_base = id_to_u_rated[&ln.to_node].powi(2) / s_base_va;
                let idx = id_to_idx[&ln.to_node];
                branch_idx.insert(ln.id, lines.len());
                half_open_terminal.insert(ln.id, Terminal::To);
                lines.push(half_open_line(ln, omega, z_base, idx));
            }
            // Both terminals open: no branch is created, so this line
            // deliberately gets no `branch_idx` entry.
            _ => {}
        }
    }

    // Transformers — convert physical-unit PGM parameters to system pu.
    let mut transformers: Vec<Transformer> = Vec::new();
    // Positions *within* `transformers`. They can't be turned into flat branch
    // indices yet: the source loop below still appends to `lines`, so the
    // offset (`lines.len()`) isn't final until the very end of this function.
    let mut transformer_pos: HashMap<u64, usize> = HashMap::new();
    let mut three_winding_pos: HashMap<u64, [usize; 3]> = HashMap::new();
    for t in &input.data.transformer {
        let tap = transformer_tap(
            t.u1,
            t.u2,
            t.tap_side,
            t.tap_pos,
            t.tap_min,
            t.tap_max,
            t.tap_nom,
            t.tap_size,
            t.clock,
            id_to_u_rated[&t.from_node],
            id_to_u_rated[&t.to_node],
        );
        let (y_series, y_shunt) =
            transformer_admittances(t.u2, t.sn, t.uk, t.pk, t.i0, t.p0, s_base_va);
        transformer_pos.insert(t.id, transformers.len());
        transformers.push(Transformer {
            from: id_to_idx[&t.from_node],
            to: id_to_idx[&t.to_node],
            from_status: t.from_status,
            to_status: t.to_status,
            y_series,
            y_shunt,
            tap,
        });
    }

    // Links: ideal connections, stamped as branches rather than merged. See
    // `PgmLink` for why, and `link_admittance` for the value.
    for ln in &input.data.link {
        if ln.from_status == 0 && ln.to_status == 0 {
            continue;
        }
        transformer_pos.insert(ln.id, transformers.len());
        transformers.push(Transformer {
            from: id_to_idx[&ln.from_node],
            to: id_to_idx[&ln.to_node],
            from_status: ln.from_status,
            to_status: ln.to_status,
            y_series: LINK_Y,
            y_shunt: Complex::new(0.0, 0.0),
            tap: Complex::new(1.0, 0.0),
        });
    }

    // Three-winding transformers — modelled as three virtual 2-winding
    // `Transformer` legs from each physical node to a synthesized internal
    // star bus (PGM's star-equivalent model: T1/T2/T3 all connect to a
    // shared dummy node whose per-unit base tracks side 1's rated voltage).
    // Star buses are appended right after the physical nodes and before
    // source virtual-slack buses, so the source loop's index offset below
    // must account for them via `n_3wdg`.
    for (i, t) in input.data.three_winding_transformer.iter().enumerate() {
        let star_idx = n_nodes + i;
        let legs = three_winding_legs(t, &id_to_u_rated, s_base_va);

        buses.push(Bus {
            idx: star_idx,
            bus_type: BusType::PQ,
            voltage_mag: 1.0,
            voltage_ang: 0.0,
            p_spec: 0.0,
            q_spec: 0.0,
            q_min: -f64::INFINITY,
            q_max: f64::INFINITY,
            u_rated: id_to_u_rated[&t.node_1],
            zip_terms: Vec::new(),
        });

        three_winding_pos.insert(
            t.id,
            [
                transformers.len(),
                transformers.len() + 1,
                transformers.len() + 2,
            ],
        );
        for leg in &legs {
            transformers.push(Transformer {
                from: id_to_idx[&leg.node],
                to: star_idx,
                from_status: leg.status,
                to_status: 1,
                y_series: leg.y_series,
                y_shunt: leg.y_shunt,
                tap: leg.tap,
            });
        }
    }
    let n_3wdg = input.data.three_winding_transformer.len();

    // Virtual Slack bus + source-impedance Line for each active source.
    for (i, src) in input
        .data
        .source
        .iter()
        .filter(|s| s.status != 0)
        .enumerate()
    {
        let virtual_idx = n_nodes + n_3wdg + i;
        let (r_s, x_s) = source_impedance_pu(src.sk, src.rx_ratio, s_base_va);
        buses.push(Bus {
            idx: virtual_idx,
            bus_type: BusType::Slack,
            voltage_mag: src.u_ref,
            voltage_ang: 0.0,
            p_spec: 0.0,
            q_spec: 0.0,
            q_min: -f64::INFINITY,
            q_max: f64::INFINITY,
            u_rated: id_to_u_rated[&src.node],
            zip_terms: Vec::new(),
        });
        source_branch_idx.insert(src.id, lines.len());
        lines.push(Line {
            from: virtual_idx,
            to: id_to_idx[&src.node],
            r: r_s,
            x: x_s,
            b_shunt: 0.0,
            g_shunt: 0.0,
        });
    }

    // `lines` is final now, so transformer positions can become flat indices.
    let n_lines = lines.len();
    branch_idx.extend(
        transformer_pos
            .into_iter()
            .map(|(id, pos)| (id, n_lines + pos)),
    );
    three_winding_branch_idx.extend(
        three_winding_pos
            .into_iter()
            .map(|(id, legs)| (id, legs.map(|pos| n_lines + pos))),
    );

    // Appliance → bus. Built in one pass at the end rather than inside the
    // accumulation closure above, which only sees active constant-power
    // appliances; a sensor may reference an inactive or ZIP-modelled one.
    let mut appliance_bus: HashMap<u64, usize> = HashMap::new();
    for (id, node) in input
        .data
        .sym_load
        .iter()
        .map(|a| (a.id, a.node))
        .chain(input.data.asym_load.iter().map(|a| (a.id, a.node)))
        .chain(input.data.sym_gen.iter().map(|a| (a.id, a.node)))
        .chain(input.data.asym_gen.iter().map(|a| (a.id, a.node)))
        .chain(input.data.shunt.iter().map(|a| (a.id, a.node)))
        .chain(input.data.source.iter().map(|a| (a.id, a.node)))
    {
        if let Some(&idx) = id_to_idx.get(&node) {
            appliance_bus.insert(id, idx);
        }
    }

    // Zero-injection buses: everything except the ones carrying an active load
    // or generator, and except the virtual slack buses where source power
    // enters. Star buses of three-winding transformers qualify, which is the
    // textbook example of the constraint being worth having.
    let mut zero_injection = vec![true; buses.len()];
    for (node, status) in input
        .data
        .sym_load
        .iter()
        .map(|a| (a.node, a.status))
        .chain(input.data.asym_load.iter().map(|a| (a.node, a.status)))
        .chain(input.data.sym_gen.iter().map(|a| (a.node, a.status)))
        .chain(input.data.asym_gen.iter().map(|a| (a.node, a.status)))
    {
        if status != 0 {
            if let Some(&idx) = id_to_idx.get(&node) {
                zero_injection[idx] = false;
            }
        }
    }
    for &branch in source_branch_idx.values() {
        zero_injection[lines[branch].from] = false;
    }

    PgmNetwork {
        buses,
        lines,
        transformers,
        zero_injection,
        node_idx: id_to_idx,
        branch_idx,
        three_winding_branch_idx,
        source_branch_idx,
        appliance_bus,
        half_open_terminal,
    }
}

/// The [`pgm_to_network`] conversion with the object-ID maps dropped — the
/// shape power flow has always used, kept so its many call sites stay
/// unchanged.
pub fn pgm_to_buses_and_branches(
    input: PgmInput,
    s_base_va: f64,
    freq_hz: f64,
) -> (Vec<Bus>, Vec<Line>, Vec<Transformer>) {
    let net = pgm_to_network(input, s_base_va, freq_hz);
    (net.buses, net.lines, net.transformers)
}

/// Converts a PGM input document (with `asym_load`) into a 3N-bus expanded
/// representation suitable for three-phase power flow.
///
/// Physical node `k` (sorted by PGM node ID) maps to buses at indices
/// `3k`, `3k+1`, `3k+2` for phases a, b, c respectively.  Each active
/// source adds one virtual Slack node appended after all physical nodes,
/// whose three phase buses carry fixed voltages `u_ref ∠ 0°`, `u_ref ∠ -120°`,
/// `u_ref ∠ +120°`.
///
/// Returns `(buses, lines_3ph, id_to_physical_idx)`.  Pass `buses.len() / 3`
/// as `n` to `build_ybus_3ph`.
/// The object-ID maps a three-phase network needs beyond its buses and branches.
///
/// [`pgm_to_3ph_network`] returns only what a power flow needs — buses, lines
/// and a node map — because a power flow never has to answer "which branch
/// terminal does this sensor name". State estimation does, so this carries the
/// same maps [`PgmNetwork`] does for the scalar case.
///
/// **Every index here is a *physical* one.** Node `k` is buses `3k..3k+3` and
/// branch `b` is terminals `3b..3b+3`, matching
/// [`SeNetwork::from_3ph`](crate::se::SeNetwork::from_3ph). Storing the physical
/// index and expanding at the point of use keeps one map serving all three
/// phases.
pub struct PgmNetwork3Ph {
    /// Node id → physical node index. Physical nodes only; the virtual slack
    /// buses behind sources have no PGM id.
    pub node_idx: HashMap<u64, usize>,
    /// `three_winding_transformer` id → its three legs' flat branch indices,
    /// in side order. A sensor names the transformer and a side; this is what
    /// turns that into a branch.
    pub three_winding_idx: HashMap<u64, [usize; 3]>,
    /// `line` id → physical flat branch index.
    ///
    /// The flat order is the one [`pgm_to_3ph_network`] produces and
    /// `SeNetwork::from_3ph` consumes: the document's own lines that yield a
    /// branch, then one synthesized branch per active source, then the
    /// transformers. A line with both terminals open yields no branch and is
    /// absent, which is why this cannot be derived from document position.
    pub branch_idx: HashMap<u64, usize>,
    /// `source` id → the physical index of its synthesized impedance branch.
    pub source_branch_idx: HashMap<u64, usize>,
    /// Appliance id → the physical node index it sits on.
    pub appliance_bus: HashMap<u64, usize>,
    /// Per *physical node*, whether its net injection is identically zero.
    pub zero_injection: Vec<bool>,
    /// Branch ids that collapsed to a self-loop because exactly one terminal was
    /// connected, mapped to which PGM terminal is the live one — see
    /// [`PgmNetwork::half_open_terminal`].
    pub half_open_terminal: HashMap<u64, Terminal>,
    /// Physical nodes including the virtual slack bus per active source.
    pub n_nodes: usize,
}

/// Builds [`PgmNetwork3Ph`] for a document, without consuming it.
///
/// Deliberately separate from [`pgm_to_3ph_network`] rather than folded into it:
/// that function consumes its input and is on the power-flow path, which has no
/// use for any of this.
///
/// Returns an error for the components the three-phase conversion does not
/// model, rather than letting them vanish. `pgm_to_3ph_network` drops all three
/// silently, which is survivable for a power flow whose caller chose the
/// document, and not for an estimate whose answer would quietly describe a
/// different network.
pub fn pgm_3ph_maps(input: &PgmInput) -> PgmNetwork3Ph {
    // `voltage_regulator` is deliberately *ignored* rather than refused, and it
    // used to be refused. These maps serve state estimation and nothing else,
    // and a voltage regulator has no meaning in an estimate: it pins a
    // generator's bus to `u_ref` and lets Q float, which is a *boundary
    // condition* for a power flow, while an estimate treats every bus voltage
    // as an unknown to be recovered from measurements. There are no PV buses
    // here for it to create.
    //
    // Refusing it therefore made an estimable document inestimable for a reason
    // that could not have changed the estimate — and inconsistently, since the
    // symmetric path reads the same document without comment. `pgm_to_network`
    // does set the bus type, and the estimator ignores it for exactly this
    // reason: every bus carries a magnitude unknown whatever its type says.
    // There used to be a winding-pair guard here, because
    // `transformer_seq_params` panicked outside Dyn and YNyn and a panic
    // reached from a document is not a diagnosis. That function now implements
    // PGM's general zero-sequence algorithm and has no unsupported pairs left
    // to guard against, so the check went with it.

    let node_idx = node_id_to_idx(input);
    let n_physical = input.data.node.len();

    // Lines, in the order `pgm_to_3ph_network` emits them: a branch per line
    // with at least one live terminal, skipping the fully-open ones.
    let mut branch_idx = HashMap::new();
    let mut half_open_terminal = HashMap::new();
    let mut next = 0usize;
    for ln in &input.data.line {
        match (ln.from_status, ln.to_status) {
            (1, 1) => {}
            (1, 0) => {
                half_open_terminal.insert(ln.id, Terminal::From);
            }
            (0, 1) => {
                half_open_terminal.insert(ln.id, Terminal::To);
            }
            _ => continue,
        }
        branch_idx.insert(ln.id, next);
        next += 1;
    }

    // Then one synthesized branch per active source, then the transformers.
    let mut source_branch_idx = HashMap::new();
    for src in input.data.source.iter().filter(|s| s.status != 0) {
        source_branch_idx.insert(src.id, next);
        next += 1;
    }
    for t in &input.data.transformer {
        branch_idx.insert(t.id, next);
        next += 1;
    }
    // Then the links, which `pgm_transformers_3ph` appends after the
    // transformers. A fully-open link is not a branch and takes no index, the
    // same rule a fully-open line follows above.
    for lk in &input.data.link {
        if lk.from_status == 0 && lk.to_status == 0 {
            continue;
        }
        branch_idx.insert(lk.id, next);
        next += 1;
    }
    // Then the three-winding transformers, three legs each. A sensor names the
    // transformer and a *side*, so all three indices are recorded against the
    // one id — `three_winding_idx` is what resolves a side to its leg, exactly
    // as the symmetric `three_winding_pos` does.
    let mut three_winding_idx = HashMap::new();
    for t in &input.data.three_winding_transformer {
        three_winding_idx.insert(t.id, [next, next + 1, next + 2]);
        next += 3;
    }

    let mut appliance_bus = HashMap::new();
    for (id, node) in input
        .data
        .sym_load
        .iter()
        .map(|a| (a.id, a.node))
        .chain(input.data.asym_load.iter().map(|a| (a.id, a.node)))
        .chain(input.data.sym_gen.iter().map(|a| (a.id, a.node)))
        .chain(input.data.asym_gen.iter().map(|a| (a.id, a.node)))
        .chain(input.data.shunt.iter().map(|a| (a.id, a.node)))
        .chain(input.data.source.iter().map(|a| (a.id, a.node)))
    {
        if let Some(&k) = node_idx.get(&node) {
            appliance_bus.insert(id, k);
        }
    }

    // A node with any load or generator injects something; sources and shunts
    // are structural here exactly as in the scalar case. Virtual slack buses are
    // excluded, being where a source's own unknown power enters.
    let n_sources = input.data.source.iter().filter(|s| s.status != 0).count();
    let n_3wdg = input.data.three_winding_transformer.len();
    let mut zero_injection = vec![true; n_physical + n_3wdg + n_sources];
    // A three-winding transformer's star node is internal: nothing is attached
    // to it, so it injects nothing, and saying so is real information rather
    // than an assumption — it is the one place a hard zero-injection constraint
    // is unarguable. The virtual buses behind sources are the opposite case,
    // and are where a source's own unknown power enters.
    for flag in zero_injection.iter_mut().skip(n_physical + n_3wdg) {
        *flag = false;
    }
    for node in input
        .data
        .sym_load
        .iter()
        .map(|a| a.node)
        .chain(input.data.asym_load.iter().map(|a| a.node))
        .chain(input.data.sym_gen.iter().map(|a| a.node))
        .chain(input.data.asym_gen.iter().map(|a| a.node))
    {
        if let Some(&k) = node_idx.get(&node) {
            zero_injection[k] = false;
        }
    }

    PgmNetwork3Ph {
        node_idx,
        branch_idx,
        three_winding_idx,
        source_branch_idx,
        appliance_bus,
        zero_injection,
        half_open_terminal,
        n_nodes: n_physical + n_3wdg + n_sources,
    }
}

/// Converts a document's `line` entries to per-unit three-phase branches.
///
/// Extracted from [`pgm_to_3ph_network`] so the short-circuit path
/// ([`pgm_to_3ph_sc_network`](shortcircuit::pgm_to_3ph_sc_network)) can reuse
/// the conversion without also picking up that function's injection
/// accumulation and virtual slack buses, neither of which a short-circuit
/// calculation wants. `pgm_transformers_3ph` and
/// `pgm_shunts_3ph` were already standalone; the line conversion was the last
/// piece still inlined.
///
/// A line with exactly one live terminal collapses to a self-loop carrying
/// only its own shunt, and a fully open line yields no branch at all — so the
/// returned order is *not* document order, which is why `PgmNetwork3Ph`
/// derives its `branch_idx` from the same rules rather than from position.
pub fn pgm_lines_3ph(
    input: &PgmInput,
    id_to_idx: &HashMap<u64, usize>,
    s_base_va: f64,
    freq_hz: f64,
) -> Vec<Line3Ph> {
    let id_to_u_rated: HashMap<u64, f64> =
        input.data.node.iter().map(|n| (n.id, n.u_rated)).collect();
    let two_pi_f = 2.0 * std::f64::consts::PI * freq_hz;

    let mut lines: Vec<Line3Ph> = Vec::new();
    for ln in &input.data.line {
        let u_rated_from = id_to_u_rated[&ln.from_node];
        let z_base = u_rated_from * u_rated_from / s_base_va;
        let b1 = two_pi_f * ln.c1 * z_base;
        let b0 = two_pi_f * ln.c0 * z_base;
        // power-grid-model's line shunt is `2πf·c·(tan δ + j)`, so the loss
        // factor turns into a conductance alongside the susceptance. `tan1`
        // and `tan0` are zero in almost every fixture, which is why this term
        // went unmodelled until a short-circuit fixture that sets them
        // (`dummy-test-line-into-itself`) made it visible.
        let g1 = b1 * ln.tan1;
        let g0 = b0 * ln.tan0;
        let r1_pu = ln.r1 / z_base;
        let x1_pu = ln.x1 / z_base;
        let r0_pu = ln.r0 / z_base;
        let x0_pu = ln.x0 / z_base;
        let from_phys = id_to_idx[&ln.from_node];
        let to_phys = id_to_idx[&ln.to_node];

        match (ln.from_status, ln.to_status) {
            (1, 1) => {
                lines.push(Line3Ph {
                    from: from_phys,
                    to: to_phys,
                    r1: r1_pu,
                    x1: x1_pu,
                    b1,
                    g1,
                    r0: r0_pu,
                    x0: x0_pu,
                    b0,
                    g0,
                });
            }
            (1, 0) => {
                lines.push(Line3Ph {
                    from: from_phys,
                    to: from_phys,
                    r1: 0.0,
                    x1: 0.0,
                    b1,
                    g1,
                    r0: 0.0,
                    x0: 0.0,
                    b0,
                    g0,
                });
            }
            (0, 1) => {
                lines.push(Line3Ph {
                    from: to_phys,
                    to: to_phys,
                    r1: 0.0,
                    x1: 0.0,
                    b1,
                    g1,
                    r0: 0.0,
                    x0: 0.0,
                    b0,
                    g0,
                });
            }
            _ => {}
        }
    }
    lines
}

pub fn pgm_to_3ph_network(
    input: PgmInput,
    s_base_va: f64,
    freq_hz: f64,
) -> (Vec<Bus>, Vec<Line3Ph>, HashMap<u64, usize>) {
    let id_to_idx = node_id_to_idx(&input);
    let id_to_u_rated: HashMap<u64, f64> =
        input.data.node.iter().map(|n| (n.id, n.u_rated)).collect();

    let n_nodes = input.data.node.len();
    let phase_ang = [
        0.0_f64,
        -2.0 * std::f64::consts::PI / 3.0,
        2.0 * std::f64::consts::PI / 3.0,
    ];

    // Per-node, per-phase net injection [phase_a, phase_b, phase_c] in p.u.,
    // plus per-phase ZIP terms for constant-current/-impedance loads and gens.
    // The phase-domain power flow equations naturally work in units of S_base/3
    // (per-phase base), so per-phase Watts must be divided by s_base_va/3.
    let mut p_inj: HashMap<u64, [f64; 3]> = HashMap::new();
    let mut q_inj: HashMap<u64, [f64; 3]> = HashMap::new();
    let mut zip_map: HashMap<u64, [Vec<ZipTerm>; 3]> = HashMap::new();
    let s_base_1ph = s_base_va / 3.0;
    let mut accumulate3 = |node: u64,
                           load_type: u8,
                           p: [f64; 3],
                           q: [f64; 3],
                           sign: f64,
                           base: f64| {
        match load_type {
            1 | 2 => {
                let kind = if load_type == 1 {
                    ZipKind::ConstImpedance
                } else {
                    ZipKind::ConstCurrent
                };
                let entry = zip_map
                    .entry(node)
                    .or_insert_with(|| [Vec::new(), Vec::new(), Vec::new()]);
                for ph in 0..3 {
                    entry[ph].push(ZipTerm {
                        s_const: Complex::new(p[ph], q[ph]) * sign / base,
                        kind,
                    });
                }
            }
            _ => {
                let pe = p_inj.entry(node).or_insert([0.0; 3]);
                let qe = q_inj.entry(node).or_insert([0.0; 3]);
                for ph in 0..3 {
                    pe[ph] += sign * p[ph] / base;
                    qe[ph] += sign * q[ph] / base;
                }
            }
        }
    };
    for load in &input.data.asym_load {
        if load.status == 0 {
            continue;
        }
        accumulate3(
            load.node,
            load.load_type,
            load.p_specified,
            load.q_specified,
            -1.0,
            s_base_1ph,
        );
    }
    for agen in &input.data.asym_gen {
        if agen.status == 0 {
            continue;
        }
        accumulate3(
            agen.node,
            agen.load_type,
            agen.p_specified,
            agen.q_specified,
            1.0,
            s_base_1ph,
        );
    }
    // sym_load/sym_gen p_specified is 3-phase total; each phase gets 1/3 of
    // total, and P_1ph_pu = (P_total/3) / (s_base/3) = P_total / s_base.
    for load in &input.data.sym_load {
        if load.status == 0 {
            continue;
        }
        let p = [load.p_specified; 3];
        let q = [load.q_specified; 3];
        accumulate3(load.node, load.load_type, p, q, -1.0, s_base_va);
    }
    for sgen in &input.data.sym_gen {
        if sgen.status == 0 {
            continue;
        }
        let p = [sgen.p_specified; 3];
        let q = [sgen.q_specified; 3];
        accumulate3(sgen.node, sgen.load_type, p, q, 1.0, s_base_va);
    }

    // Build 3N buses: physical node k → buses 3k, 3k+1, 3k+2.
    let mut sorted_ids: Vec<u64> = input.data.node.iter().map(|n| n.id).collect();
    sorted_ids.sort_unstable();

    let mut buses: Vec<Bus> = Vec::with_capacity(3 * n_nodes);
    for _ in 0..3 * n_nodes {
        buses.push(Bus {
            idx: 0,
            bus_type: BusType::PQ,
            voltage_mag: 1.0,
            voltage_ang: 0.0,
            p_spec: 0.0,
            q_spec: 0.0,
            q_min: -f64::INFINITY,
            q_max: f64::INFINITY,
            u_rated: 0.0,
            zip_terms: Vec::new(),
        });
    }
    for id in &sorted_ids {
        let phys = id_to_idx[id];
        let p_arr = p_inj.get(id).copied().unwrap_or([0.0; 3]);
        let q_arr = q_inj.get(id).copied().unwrap_or([0.0; 3]);
        let mut zip_arr = zip_map.remove(id).unwrap_or_default();
        for ph in 0..3 {
            let bus_idx = 3 * phys + ph;
            buses[bus_idx] = Bus {
                idx: bus_idx,
                bus_type: BusType::PQ,
                voltage_mag: 1.0,
                voltage_ang: phase_ang[ph],
                p_spec: p_arr[ph],
                q_spec: q_arr[ph],
                q_min: -f64::INFINITY,
                q_max: f64::INFINITY,
                u_rated: id_to_u_rated[id],
                zip_terms: std::mem::take(&mut zip_arr[ph]),
            };
        }
    }

    let mut lines = pgm_lines_3ph(&input, &id_to_idx, s_base_va, freq_hz);

    // Star nodes for the three-winding transformers, three buses each, between
    // the physical nodes and the sources' virtual ones — the same layout the
    // symmetric conversion uses, so a bus index means the same thing in both.
    // Their branches are built by `pgm_transformers_3ph`, which is where every
    // three-phase caller gets its transformers.
    for (i, t) in input.data.three_winding_transformer.iter().enumerate() {
        let star_phys = n_nodes + i;
        for ph in 0..3 {
            buses.push(Bus {
                idx: 3 * star_phys + ph,
                bus_type: BusType::PQ,
                voltage_mag: 1.0,
                voltage_ang: phase_ang[ph],
                p_spec: 0.0,
                q_spec: 0.0,
                q_min: -f64::INFINITY,
                q_max: f64::INFINITY,
                u_rated: id_to_u_rated[&t.node_1],
                zip_terms: Vec::new(),
            });
        }
    }
    let n_3wdg = input.data.three_winding_transformer.len();

    // Virtual Slack buses + source-impedance lines for each active source.
    for (i, src) in input
        .data
        .source
        .iter()
        .filter(|s| s.status != 0)
        .enumerate()
    {
        let virtual_phys = n_nodes + n_3wdg + i;
        let (r1_s, x1_s, r0_s, x0_s) =
            source_impedance_pu_seq(src.sk, src.rx_ratio, src.z01_ratio, s_base_va);

        for ph in 0..3 {
            let bus_idx = 3 * virtual_phys + ph;
            buses.push(Bus {
                idx: bus_idx,
                bus_type: BusType::Slack,
                voltage_mag: src.u_ref,
                voltage_ang: phase_ang[ph],
                p_spec: 0.0,
                q_spec: 0.0,
                q_min: -f64::INFINITY,
                q_max: f64::INFINITY,
                u_rated: id_to_u_rated[&src.node],
                zip_terms: Vec::new(),
            });
        }

        // A source's internal impedance is a pure series branch — no shunt at
        // all, hence no charging susceptance and no dielectric loss.
        lines.push(Line3Ph {
            from: virtual_phys,
            to: id_to_idx[&src.node],
            r1: r1_s,
            x1: x1_s,
            b1: 0.0,
            g1: 0.0,
            r0: r0_s,
            x0: x0_s,
            b0: 0.0,
            g0: 0.0,
        });
    }

    (buses, lines, id_to_idx)
}