powerio_tx/

normalize.rs

//! The universal normalization shared by the PowerModels reader/writer and
//! [`BalancedNetwork::to_normalized`].
//!
//! Two things live here so there is one implementation of each:
//!
//! - **Per-unit scaling factors and the gen-cost rescale** ([`cost_to_pu`] /
//!   [`cost_from_pu`], [`DEG_TO_RAD`] / [`RAD_TO_DEG`], [`GEN_PU_KEYS`]). The
//!   PowerModels writer scales raw model values into its per-unit JSON; the
//!   reader inverts it; [`BalancedNetwork::to_normalized`] scales the same way into a new
//!   `BalancedNetwork`. The cost rescale is the one piece subtle enough that a second copy
//!   would drift, so it has a single home.
//! - **[`BalancedNetwork::to_normalized`]**: a derived, computation-ready form, per unit,
//!   radians, out of service filtered, source ID preserving, bus types canonicalized.

use std::collections::{HashMap, HashSet};

use crate::islands::{
    IslandPartition, IslandReferencePolicy, bus_list, kept_reference, largest_generator_bus,
};
use crate::merge::{BusMerge, BusMergeRule};
use crate::network::{
    BalancedNetwork, BalancedNetworkTables, Branch, Bus, BusId, BusType, GEN_EXTRA_KEYS, GenCost,
    Generator, Hvdc, Load, LoadVoltageModel, Shunt, SourceFormat, StaticVarCompensator, Storage,
    Switch, Transformer3W, TransformerControl, TransformerControlMode,
};
use crate::{Error, Result};

/// Degrees → radians. The per-unit convention stores angles in radians; the raw
/// model keeps MATPOWER degrees.
pub(crate) const DEG_TO_RAD: f64 = std::f64::consts::PI / 180.0;

/// Radians → degrees, the inverse of [`DEG_TO_RAD`], used when reading a per-unit
/// source back into the neutral degree model.
pub(crate) const RAD_TO_DEG: f64 = 180.0 / std::f64::consts::PI;

fn norm_transformer_control(
    control: &mut TransformerControl,
    base_mva: f64,
    map: &HashMap<BusId, BusId>,
) {
    control.controlled_bus = control
        .controlled_bus
        .and_then(|controlled_bus| remap(map, controlled_bus));
    // Siemens defines RMA/RMI as phase shift angles only for |COD| 3 and 5.
    // The values are unused for |COD| 0 and 4, so Fixed and DcLineQuantity
    // retain their source values instead of assigning units they do not have.
    if matches!(
        control.mode,
        TransformerControlMode::ActiveFlow | TransformerControlMode::AsymmetricActiveFlow
    ) {
        control.tap_min *= DEG_TO_RAD;
        control.tap_max *= DEG_TO_RAD;
    }
    // Siemens defines VMA/VMI as Mvar for |COD| 2 and MW for |COD| 3 and 5.
    // They are unused for |COD| 0 and 4.
    if matches!(
        control.mode,
        TransformerControlMode::ReactiveFlow
            | TransformerControlMode::ActiveFlow
            | TransformerControlMode::AsymmetricActiveFlow
    ) {
        control.band_min /= base_mva;
        control.band_max /= base_mva;
    }
    if let Some(angle) = &mut control.winding_connection_angle {
        *angle *= DEG_TO_RAD;
    }
}

/// The gen capability columns that are per-unitized (the ramp rates). The PQ-curve
/// points (`pc1`/`pc2`/`qc*`) and `apf` stay raw, exactly as PowerModels'
/// `make_per_unit!` leaves them, so a column is scaled in one place and can't drift
/// between the reader, the writer, and [`BalancedNetwork::to_normalized`].
pub(crate) const GEN_PU_KEYS: [&str; 4] = ["ramp_agc", "ramp_10", "ramp_30", "ramp_q"];

/// Default branch angle difference bound used by PowerModels parse time repair.
#[allow(clippy::approx_constant)]
pub const POWER_MODELS_ANGLE_BOUND_PAD: f64 = 1.0472;

/// What normalization does with a closed switch that joins two buses.
#[derive(Clone, Copy, Debug, Default, PartialEq, Eq)]
#[non_exhaustive]
pub enum ClosedSwitchPolicy {
    /// Merge the buses closed switches join before the pass, with
    /// [`BusMergeRule::closed_switches`], and report the merge under
    /// `CANONICALIZE.MERGE`. The matrix builders and solver preparations
    /// stamp branches only, so a normalized network carries no closed switch
    /// for them to leave apart. A closed switch on an isolated or undeclared
    /// bus is left to the filters that drop that bus. Returns
    /// [`Error::ClosedSwitch`] if a requested merge is refused, for example
    /// when it would join two windings of one three winding transformer.
    #[default]
    Merge,
    /// Refuse a network with such a switch with [`Error::ClosedSwitch`]
    /// (`BUILD.SWITCH.CLOSED`).
    Refuse,
    /// Keep closed switches as data, as normalization did before it merged
    /// them, for a consumer that models switches itself.
    Keep,
}

/// Options for [`BalancedNetwork::to_normalized_with_options`].
#[derive(Clone, Copy, Debug, PartialEq)]
pub struct NormalizeOptions {
    /// Clamp branch angle difference bounds to the interval PowerModels relaxations
    /// accept. Disabled by default so [`BalancedNetwork::to_normalized`] stays unchanged.
    pub clamp_angle_bounds: bool,
    /// Replacement magnitude, in radians, for clamped angle bounds.
    pub angle_bound_pad: f64,
}

impl Default for NormalizeOptions {
    fn default() -> Self {
        Self {
            clamp_angle_bounds: false,
            angle_bound_pad: POWER_MODELS_ANGLE_BOUND_PAD,
        }
    }
}

/// Options for [`BalancedNetwork::to_normalized_with`]: the
/// [`NormalizeOptions`] repairs and the normalization choices added after
/// them. Build one from its default or from [`NormalizeOptions`] and set the
/// rest with the `with_*` methods.
#[derive(Clone, Copy, Debug, Default, PartialEq)]
#[non_exhaustive]
pub struct NormalizeSettings {
    /// The repairs [`BalancedNetwork::to_normalized_with_options`] applies.
    pub options: NormalizeOptions,
    /// What happens to closed switches; by default they merge.
    pub closed_switches: ClosedSwitchPolicy,
    /// How the reference buses of a network with several islands are chosen.
    /// [`IslandReferencePolicy::Stated`], the default, keeps the behavior of
    /// [`BalancedNetwork::to_normalized`];
    /// [`IslandReferencePolicy::PerIsland`] gives every island one reference
    /// and leaves out the islands no in-service generator supplies.
    pub island_references: IslandReferencePolicy,
}

impl NormalizeSettings {
    /// `options` with the default choice for everything else.
    #[must_use]
    pub const fn new(options: NormalizeOptions) -> Self {
        Self {
            options,
            closed_switches: ClosedSwitchPolicy::Merge,
            island_references: IslandReferencePolicy::Stated,
        }
    }

    #[must_use]
    pub const fn with_options(mut self, options: NormalizeOptions) -> Self {
        self.options = options;
        self
    }

    #[must_use]
    pub const fn with_closed_switches(mut self, policy: ClosedSwitchPolicy) -> Self {
        self.closed_switches = policy;
        self
    }

    #[must_use]
    pub const fn with_island_references(mut self, policy: IslandReferencePolicy) -> Self {
        self.island_references = policy;
        self
    }
}

impl From<NormalizeOptions> for NormalizeSettings {
    fn from(options: NormalizeOptions) -> Self {
        Self::new(options)
    }
}

/// Output of [`BalancedNetwork::to_normalized_with_options`].
#[derive(Clone, Debug)]
pub struct NormalizedNetwork {
    pub network: BalancedNetwork,
    /// The pass's findings as structured records.
    pub diagnostics: Vec<crate::diagnostics::Diagnostic>,
    /// The same findings as `CODE: message` lines.
    pub warnings: Vec<String>,
}

/// What [`BalancedNetwork::to_normalized_with`] changed in the network before
/// the pass itself, beside the [`NormalizedNetwork`].
#[derive(Clone, Debug, Default)]
#[non_exhaustive]
pub struct NormalizeReport {
    /// The closed switch merge under [`ClosedSwitchPolicy::Merge`]: its bus
    /// map, removed switches, and row maps from the source. `None` when no
    /// closed switch joined two buses the pass keeps, or under another
    /// policy.
    pub bus_merge: Option<BusMerge>,
}

/// Row provenance for one normalize pass: for each dense position in the
/// [`IndexedNetwork`](crate::IndexedNetwork) view of the normalized network,
/// the row of the same element family in the source network. `None` marks an
/// element the pipeline synthesized, which has no source row.
///
/// Every field is positional over that view, so `buses[dense_index]` resolves
/// a matrix row back to its source element. Each length equals the matching
/// element table of the view: `buses` equals `view.n()`, `branches` equals
/// `view.branches().len()`.
///
/// The star lowering that the view applies to a 3-winding transformer appends
/// one bus, its star branches, and a magnetizing shunt. Those entries are
/// `None`. The lowering also consumes the transformer itself, so the view
/// holds none; `transformers_3w` therefore stays positional over the
/// normalized network's own list.
///
/// The map is valid only for the [`NormalizedNetwork`] returned beside it. A
/// later mutation of that network ([`BalancedNetwork::merge_bus`],
/// [`BalancedNetwork::reduce_zero_impedance`], [`BalancedNetwork::reduce_passthrough_buses`],
/// [`BalancedNetwork::subset`], or a hand edit) invalidates every entry; run the pass
/// again instead of patching the map.
#[derive(Clone, Debug)]
#[non_exhaustive]
pub struct NormalizeSourceRows {
    pub buses: Vec<Option<usize>>,
    pub loads: Vec<Option<usize>>,
    pub shunts: Vec<Option<usize>>,
    pub static_var_compensators: Vec<Option<usize>>,
    pub branches: Vec<Option<usize>>,
    pub switches: Vec<Option<usize>>,
    pub generators: Vec<Option<usize>>,
    pub storage: Vec<Option<usize>>,
    pub hvdc: Vec<Option<usize>>,
    pub transformers_3w: Vec<Option<usize>>,
}

impl NormalizeSourceRows {
    /// Rows over a merged network, carried back to the network `merge`
    /// read. Buses resolve by id; a removed branch or switch has no merged
    /// row, and a shunt the merge added for a removed branch's charging has
    /// no source row. The other families keep their rows.
    fn through_merge(mut self, merge: &BusMerge) -> Self {
        let source_bus: HashMap<BusId, usize> = merge
            .source
            .buses()
            .iter()
            .enumerate()
            .map(|(row, bus)| (bus.id, row))
            .collect();
        let merged_buses = merge.network.buses();
        for row in self.buses.iter_mut().flatten() {
            *row = source_bus[&merged_buses[*row].id];
        }
        let inverse = |rows: &[Option<usize>]| {
            let mut inverse = vec![0; rows.iter().flatten().count()];
            for (source, merged) in rows.iter().enumerate() {
                if let Some(merged) = merged {
                    inverse[*merged] = source;
                }
            }
            inverse
        };
        let branches = inverse(&merge.branch_rows);
        for row in self.branches.iter_mut().flatten() {
            *row = branches[*row];
        }
        let switches = inverse(&merge.switch_rows);
        for row in self.switches.iter_mut().flatten() {
            *row = switches[*row];
        }
        let source_shunts = merge.source.shunts().len();
        for row in &mut self.shunts {
            if row.is_some_and(|row| row >= source_shunts) {
                *row = None;
            }
        }
        self
    }

    /// Grow the families the star lowering appends to so each length matches the
    /// lowered form of `net`. The appended entries have no source row. The
    /// lengths come from [`BalancedNetwork::lowered_lengths`], which counts them off the
    /// transformer records, so padding never builds the lowering itself.
    pub(crate) fn pad_to_lowered(&mut self, net: &BalancedNetwork) {
        let lengths = net.lowered_lengths();
        self.buses.resize(lengths.buses, None);
        self.branches.resize(lengths.branches, None);
        self.shunts.resize(lengths.shunts, None);
    }
}

#[derive(Clone, Copy, Debug, PartialEq, Eq)]
enum CostModel {
    Piecewise,
    Polynomial,
    Unknown,
}

impl From<u8> for CostModel {
    fn from(value: u8) -> Self {
        match value {
            1 => CostModel::Piecewise,
            2 => CostModel::Polynomial,
            _ => CostModel::Unknown,
        }
    }
}

/// Gen cost coefficients rescaled into the per-unit basis, trimmed to the length
/// the model implies (a polynomial keeps `ncost` coeffs; a piecewise curve keeps
/// `2·ncost` `(mw, cost)` values). MATPOWER pads every gencost row to the matrix
/// width with trailing zeros; the padding would make a polynomial read as a
/// higher-degree curve and mis-scale, so it is dropped here.
///
/// Polynomial (model 2): coeff `i` is the term `p^(k-1-i)`, so per unit scales it
/// by `base^(k-1-i)`. Piecewise (model 1): the MW breakpoints (even positions) are
/// divided by `base`; the cost ordinates (odd positions) stay. Any other model has
/// unknown coefficient semantics, so it passes through untouched — the exact
/// inverse of [`cost_from_pu`]'s own passthrough.
pub(crate) fn cost_to_pu(cost: &GenCost, base: f64) -> Vec<f64> {
    let mut coeffs = cost.coeffs.clone();
    scale_coeffs_to_pu(&mut coeffs, cost.ncost, cost.model, base);
    coeffs
}

/// [`cost_to_pu`] over a vector the caller already owns, so a rescale in place
/// keeps its allocation.
pub(crate) fn scale_coeffs_to_pu(coeffs: &mut Vec<f64>, ncost: usize, model: u8, base: f64) {
    match CostModel::from(model) {
        CostModel::Polynomial => {
            coeffs.truncate(ncost.min(coeffs.len()));
            let k = coeffs.len();
            // The exponent k-1-i is in [0, k-1]; a polynomial never has i32::MAX-many
            // terms, so the conversion can't fail (loud, not silent, if it ever did).
            for (i, c) in coeffs.iter_mut().enumerate() {
                *c *= base.powi(i32::try_from(k - 1 - i).expect("cost degree fits i32"));
            }
        }
        CostModel::Piecewise => {
            // saturating_mul: `ncost` comes from input (JSON deserializes it
            // unchecked), so an oversized count must clamp to the coefficient
            // length instead of overflowing.
            coeffs.truncate(ncost.saturating_mul(2).min(coeffs.len()));
            for c in coeffs.iter_mut().step_by(2) {
                *c /= base;
            }
        }
        CostModel::Unknown => {}
    }
}

/// Undo [`cost_to_pu`] for the neutral MW basis: a polynomial (model 2) divides
/// coeff `i` by `base^(k-1-i)`, a piecewise curve (model 1) multiplies its MW
/// breakpoints (even positions) by `base`. The exact inverse of [`cost_to_pu`] on
/// the trimmed coefficient vector — JSON-sourced coefficients arrive already
/// trimmed, so this does no trimming; other models pass through unchanged.
pub(crate) fn cost_from_pu(coeffs: &[f64], model: u8, base: f64) -> Vec<f64> {
    let k = coeffs.len();
    match CostModel::from(model) {
        CostModel::Polynomial => coeffs
            .iter()
            .enumerate()
            .map(|(i, &c)| c / base.powi(i32::try_from(k - 1 - i).expect("cost degree fits i32")))
            .collect(),
        CostModel::Piecewise => coeffs
            .iter()
            .enumerate()
            .map(|(i, &c)| if i % 2 == 0 { c * base } else { c })
            .collect(),
        CostModel::Unknown => coeffs.to_vec(),
    }
}

/// Map a source bus id to its surviving normalized id, or `None` if the bus was dropped.
fn remap(map: &HashMap<BusId, BusId>, id: BusId) -> Option<BusId> {
    map.get(&id).copied()
}

fn norm_loads(
    loads: &[Load],
    base: f64,
    map: &HashMap<BusId, BusId>,
) -> (Vec<Load>, Vec<Option<usize>>) {
    loads
        .iter()
        .enumerate()
        .filter(|(_, l)| l.in_service)
        .filter_map(|(row, l)| {
            Some((
                Load {
                    bus: remap(map, l.bus)?,
                    p: l.p / base,
                    q: l.q / base,
                    voltage_model: l
                        .voltage_model
                        .as_ref()
                        .map(|m| norm_load_voltage_model(m, base)),
                    ..l.clone()
                },
                Some(row),
            ))
        })
        .unzip()
}

fn norm_load_voltage_model(model: &LoadVoltageModel, base: f64) -> LoadVoltageModel {
    match model {
        LoadVoltageModel::ConstantPower => LoadVoltageModel::ConstantPower,
        LoadVoltageModel::Zip {
            p_constant_power,
            q_constant_power,
            p_constant_current,
            q_constant_current,
            p_constant_impedance,
            q_constant_impedance,
            v_nom,
            load_type,
            scaling,
        } => LoadVoltageModel::Zip {
            p_constant_power: p_constant_power / base,
            q_constant_power: q_constant_power / base,
            p_constant_current: p_constant_current / base,
            q_constant_current: q_constant_current / base,
            p_constant_impedance: p_constant_impedance / base,
            q_constant_impedance: q_constant_impedance / base,
            v_nom: *v_nom,
            load_type: *load_type,
            scaling: *scaling,
        },
        LoadVoltageModel::Exponential {
            p,
            q,
            v_nom,
            gamma_p,
            gamma_q,
        } => LoadVoltageModel::Exponential {
            p: p / base,
            q: q / base,
            v_nom: *v_nom,
            gamma_p: *gamma_p,
            gamma_q: *gamma_q,
        },
    }
}

fn norm_shunts(
    shunts: &[Shunt],
    base: f64,
    map: &HashMap<BusId, BusId>,
) -> (Vec<Shunt>, Vec<Option<usize>>) {
    shunts
        .iter()
        .enumerate()
        .filter(|(_, s)| s.in_service)
        .filter_map(|(row, s)| {
            let mut shunt = s.clone();
            shunt.bus = remap(map, s.bus)?;
            shunt.g = s.g / base;
            shunt.b = s.b / base;
            // Remap the switched-shunt control bus and drop it if its target was
            // filtered out, so the normalized network has no dangling reference.
            if let Some(c) = &mut shunt.control {
                c.control_bus = c.control_bus.and_then(|b| remap(map, b));
                for block in &mut c.blocks {
                    block.g /= base;
                    block.b /= base;
                }
            }
            Some((shunt, Some(row)))
        })
        .unzip()
}

fn norm_static_var_compensators(
    compensators: &[StaticVarCompensator],
    base: f64,
    map: &HashMap<BusId, BusId>,
) -> (Vec<StaticVarCompensator>, Vec<Option<usize>>) {
    compensators
        .iter()
        .enumerate()
        .filter(|(_, svc)| svc.in_service)
        .filter_map(|(row, svc)| {
            let mut normalized = svc.clone();
            normalized.bus = remap(map, svc.bus)?;
            normalized.p = svc.p / base;
            normalized.q = svc.q / base;
            Some((normalized, Some(row)))
        })
        .unzip()
}

fn norm_branches(
    branches: &[Branch],
    base: f64,
    map: &HashMap<BusId, BusId>,
) -> (Vec<Branch>, Vec<Option<usize>>) {
    branches
        .iter()
        .enumerate()
        .filter(|(_, br)| br.in_service)
        .filter_map(|(row, br)| {
            let mut branch = br.clone();
            branch.from = remap(map, br.from)?;
            branch.to = remap(map, br.to)?;
            branch.rate_a = br.rate_a / base;
            branch.rate_b = br.rate_b / base;
            branch.rate_c = br.rate_c / base;
            for set in &mut branch.rating_sets {
                set.rate_mva /= base;
            }
            branch.tap = br.calc_effective_tap();
            branch.shift = br.shift * DEG_TO_RAD;
            branch.angmin = br.angmin * DEG_TO_RAD;
            branch.angmax = br.angmax * DEG_TO_RAD;
            if let Some(s) = &mut branch.solution {
                s.pf /= base;
                s.qf /= base;
                s.pt /= base;
                s.qt /= base;
            }
            // Remap the regulated-bus reference through the id map and drop it
            // if its target was filtered out (out of service / isolated), so the
            // normalized network has no dangling control reference.
            if let Some(c) = &mut branch.control {
                norm_transformer_control(c, base, map);
            }
            Some((branch, Some(row)))
        })
        .unzip()
}

fn validate_normalize_options(options: &NormalizeOptions) -> Result<()> {
    if options.clamp_angle_bounds
        && (!options.angle_bound_pad.is_finite()
            || options.angle_bound_pad <= 0.0
            || options.angle_bound_pad >= std::f64::consts::FRAC_PI_2)
    {
        return Err(Error::InvalidNormalizeOption {
            field: "angle_bound_pad",
            value: options.angle_bound_pad,
        });
    }
    Ok(())
}

#[allow(clippy::float_cmp)] // exact replacement accounting, not numerical equivalence
fn clamp_angle_bounds(
    branches: &mut [Branch],
    pad: f64,
    warnings: &mut crate::diagnostics::Diagnostics,
) {
    for (idx, br) in branches.iter_mut().enumerate() {
        let old_min = br.angmin;
        let old_max = br.angmax;
        let mut changes = Vec::new();

        let (corrected_min, corrected_max) =
            correct_angle_difference_bounds_with_pad(old_min, old_max, pad, 0.0);

        if old_min != corrected_min {
            br.angmin = corrected_min;
            changes.push(format!("angmin {old_min} -> {}", br.angmin));
        }
        if old_max != corrected_max {
            br.angmax = corrected_max;
            changes.push(format!("angmax {old_max} -> {}", br.angmax));
        }

        if !changes.is_empty() {
            warnings.push(
                &crate::diagnostics::codes::CANONICALIZE_NORMALIZE_BOUNDS_CLAMPED,
                format!(
                    "branch {idx} angle difference bounds clamped: {}",
                    changes.join(", ")
                ),
            );
        }
    }
}

/// Correct one branch angle difference interval using the same ±60 degree
/// rule as PowerModels' `correct_voltage_angle_differences!`.
///
/// A lower bound at or below −90 degrees becomes −60 degrees, an upper bound
/// at or above 90 degrees becomes 60 degrees, and the MATPOWER 0/0 spelling
/// becomes ±60 degrees. If correcting one side would invert the interval, the
/// result is also ±60 degrees. Inputs and outputs are radians.
///
/// This is the pad centered on zero, which suits an unshifted branch only;
/// [`AngleDifferenceBounds::ShiftCenteredPad`] centers it on the branch's
/// phase shift.
#[must_use]
pub fn correct_angle_difference_bounds(angle_min: f64, angle_max: f64) -> (f64, f64) {
    correct_angle_difference_bounds_with_pad(
        angle_min,
        angle_max,
        POWER_MODELS_ANGLE_BOUND_PAD,
        0.0,
    )
}

/// One branch angle difference interval as the source states it, for a
/// consumer that applies no pad. Inputs and outputs are radians.
///
/// The interval constrains `theta_from - theta_to`, as MATPOWER and
/// PowerModels define it. The MATPOWER `0`/`0` spelling states no
/// constraint. A lower side at or below −360 degrees and an upper side at or
/// above 360 degrees also denote absent limits. Absent sides become negative
/// or positive infinity, respectively: unwrapped DC angle differences can
/// exceed a full turn. Every other stated bound passes through unchanged.
#[must_use]
pub fn stated_angle_difference_bounds(angle_min: f64, angle_max: f64) -> (f64, f64) {
    let full = std::f64::consts::TAU;
    if angle_min == 0.0 && angle_max == 0.0 {
        return (f64::NEG_INFINITY, f64::INFINITY);
    }
    let angle_min = if angle_min <= -full {
        f64::NEG_INFINITY
    } else {
        angle_min
    };
    let angle_max = if angle_max >= full {
        f64::INFINITY
    } else {
        angle_max
    };
    (angle_min, angle_max)
}

/// Which angle difference bounds an OPF or power flow preparation carries.
///
/// Every bound constrains `theta_from - theta_to`, the MATPOWER and
/// PowerModels convention, which is not net of a branch's phase shift.
///
/// The assembly options select one with `with_angle_difference_bounds`. Left
/// unselected, their `correct_angle_difference_bounds` flag decides as it
/// always has: `true`, the default, is [`Self::PowerModelsPad`] and `false`
/// is [`Self::Stated`].
#[derive(Clone, Copy, Debug, PartialEq, Eq, serde::Serialize, serde::Deserialize)]
#[serde(rename_all = "snake_case")]
#[non_exhaustive]
pub enum AngleDifferenceBounds {
    /// PowerModels' ±60 degree pad exactly as its
    /// `correct_voltage_angle_differences!` applies it
    /// ([`correct_angle_difference_bounds`]): a side the source leaves
    /// unconstrained or unusable (at or beyond ±90 degrees, the `0`/`0`
    /// spelling, or an inverted window) becomes ±60 degrees around zero, and
    /// a usable stated side is kept. A branch that shifts phase by more than
    /// 60 degrees cannot meet this window; select
    /// [`Self::ShiftCenteredPad`] or [`Self::Stated`] for such a case.
    PowerModelsPad,
    /// The same pad centered on the branch's phase shift, so a phase
    /// shifting branch keeps its operating point inside the window. An
    /// unshifted branch gets exactly [`Self::PowerModelsPad`].
    ShiftCenteredPad,
    /// The bounds the source states ([`stated_angle_difference_bounds`]):
    /// an absent side becomes negative or positive infinity.
    Stated,
    /// No angle difference bound: every interval is `(−∞, +∞)`.
    None,
}

impl AngleDifferenceBounds {
    /// The selection the legacy `correct_angle_difference_bounds` flag
    /// stands for: [`Self::PowerModelsPad`] when `true`, [`Self::Stated`]
    /// when `false`.
    #[must_use]
    pub const fn from_correct_flag(correct: bool) -> Self {
        if correct {
            Self::PowerModelsPad
        } else {
            Self::Stated
        }
    }

    /// Whether this selection pads unconstrained sides, which is what the
    /// preparations' `correct_angle_difference_bounds` flag reports.
    #[must_use]
    pub const fn is_padded(self) -> bool {
        matches!(self, Self::PowerModelsPad | Self::ShiftCenteredPad)
    }

    /// The interval for one branch, in radians, from its stated `angmin`
    /// and `angmax` and the phase shift the preparation models.
    #[must_use]
    pub fn calc_bounds(self, angle_min: f64, angle_max: f64, shift: f64) -> (f64, f64) {
        match self {
            Self::PowerModelsPad => correct_angle_difference_bounds(angle_min, angle_max),
            Self::ShiftCenteredPad => correct_angle_difference_bounds_with_pad(
                angle_min,
                angle_max,
                POWER_MODELS_ANGLE_BOUND_PAD,
                shift,
            ),
            Self::Stated => stated_angle_difference_bounds(angle_min, angle_max),
            Self::None => (f64::NEG_INFINITY, f64::INFINITY),
        }
    }
}

/// The PowerModels rule with the pad centered on `center`, the branch's
/// phase shift in radians: PowerModels itself centers it on zero.
#[allow(clippy::float_cmp)] // the `0`/`0` spelling is an exact sentinel
fn correct_angle_difference_bounds_with_pad(
    angle_min: f64,
    angle_max: f64,
    pad: f64,
    center: f64,
) -> (f64, f64) {
    let (mut low, mut high) = (angle_min, angle_max);
    if low <= -std::f64::consts::FRAC_PI_2 {
        low = center - pad;
    }
    if high >= std::f64::consts::FRAC_PI_2 {
        high = center + pad;
    }
    if angle_min == 0.0 && angle_max == 0.0 || low > high {
        return (center - pad, center + pad);
    }
    (low, high)
}

fn norm_gens(
    gens: &[Generator],
    base: f64,
    map: &HashMap<BusId, BusId>,
) -> (Vec<Generator>, Vec<Option<usize>>) {
    gens.iter()
        .enumerate()
        .filter(|(_, g)| g.in_service)
        .filter_map(|(row, g)| {
            let mut generator = g.clone();
            generator.bus = remap(map, g.bus)?;
            generator.pg = g.pg / base;
            generator.qg = g.qg / base;
            generator.pmax = g.pmax / base;
            generator.pmin = g.pmin / base;
            generator.qmax = g.qmax / base;
            generator.qmin = g.qmin / base;
            if let Some(c) = &mut generator.cost {
                scale_coeffs_to_pu(&mut c.coeffs, c.ncost, c.model, base);
            }
            // `GenCaps` is indexed by `GEN_EXTRA_KEYS`, so the two zip exactly.
            for (cap, key) in generator.caps.iter_mut().zip(GEN_EXTRA_KEYS) {
                if GEN_PU_KEYS.contains(&key)
                    && let Some(v) = cap
                {
                    *v /= base;
                }
            }
            // Remap the regulated bus through the same id map; drop it if its
            // target was filtered out so the normalized form stays consistent.
            generator.regulated_bus = g.regulated_bus.and_then(|b| remap(map, b));
            Some((generator, Some(row)))
        })
        .unzip()
}

fn norm_switches(
    switches: &[Switch],
    base: f64,
    map: &HashMap<BusId, BusId>,
) -> (Vec<Switch>, Vec<Option<usize>>) {
    switches
        .iter()
        .enumerate()
        .filter_map(|(row, s)| {
            let switch = Switch {
                from: remap(map, s.from)?,
                to: remap(map, s.to)?,
                thermal_rating: s.thermal_rating.map(|v| v / base),
                pf: s.pf.map(|v| v / base),
                qf: s.qf.map(|v| v / base),
                pt: s.pt.map(|v| v / base),
                qt: s.qt.map(|v| v / base),
                ..s.clone()
            };
            Some((switch, Some(row)))
        })
        .unzip()
}

fn norm_storage(
    storage: &[Storage],
    base: f64,
    map: &HashMap<BusId, BusId>,
) -> (Vec<Storage>, Vec<Option<usize>>) {
    storage
        .iter()
        .enumerate()
        .filter(|(_, s)| s.in_service)
        .filter_map(|(row, s)| {
            // ps/qs stay raw (PowerModels' make_per_unit! leaves the dispatch
            // setpoint alone); the energy, ratings, limits, and losses scale.
            let unit = Storage {
                bus: remap(map, s.bus)?,
                energy: s.energy / base,
                energy_rating: s.energy_rating / base,
                charge_rating: s.charge_rating / base,
                discharge_rating: s.discharge_rating / base,
                thermal_rating: s.thermal_rating / base,
                qmin: s.qmin / base,
                qmax: s.qmax / base,
                p_loss: s.p_loss / base,
                q_loss: s.q_loss / base,
                ..s.clone()
            };
            Some((unit, Some(row)))
        })
        .unzip()
}

fn norm_hvdc(
    hvdc: &[Hvdc],
    base: f64,
    map: &HashMap<BusId, BusId>,
) -> (Vec<Hvdc>, Vec<Option<usize>>) {
    hvdc.iter()
        .enumerate()
        .filter(|(_, d)| d.in_service)
        .filter_map(|(row, d)| {
            // No sign flip: the writer's Pt/Qf/Qt negation is a PowerModels output
            // convention, not part of per-unit normalization. The aggregate
            // pmin/pmax stay raw, matching make_per_unit!.
            let mut link = d.clone();
            link.from = remap(map, d.from)?;
            link.to = remap(map, d.to)?;
            link.pf = d.pf / base;
            link.pt = d.pt / base;
            link.qf = d.qf / base;
            link.qt = d.qt / base;
            link.qminf = d.qminf / base;
            link.qmaxf = d.qmaxf / base;
            link.qmint = d.qmint / base;
            link.qmaxt = d.qmaxt / base;
            link.loss0 = d.loss0 / base;
            if let Some(c) = &mut link.cost {
                scale_coeffs_to_pu(&mut c.coeffs, c.ncost, c.model, base);
            }
            Some((link, Some(row)))
        })
        .unzip()
}

fn norm_transformers_3w(
    xfmrs: &[Transformer3W],
    base: f64,
    map: &HashMap<BusId, BusId>,
    warnings: &mut crate::diagnostics::Diagnostics,
) -> (Vec<Transformer3W>, Vec<Option<usize>>) {
    let mut taken_out = Vec::new();
    let normalized = xfmrs
        .iter()
        .enumerate()
        .filter(|(_, t)| t.in_service)
        .filter_map(|(row, t)| {
            // A winding whose bus this pass drops goes out of service, as PSS/E
            // takes out a winding on a disconnected bus; the other two stay
            // coupled through the star point. Its bus field then names a
            // kept winding's bus, so the record stays reference consistent,
            // and the source row gives the original. A transformer left with
            // fewer than two windings couples nothing and is dropped. Phase
            // shifts and the star angle go to radians; winding ratings go per
            // unit; the pairwise impedances are already per unit on the system
            // base.
            let mut live = t.winding_in_service();
            for (state, winding) in live.iter_mut().zip(&t.windings) {
                *state &= remap(map, winding.bus).is_some();
            }
            let kept_bus = t
                .windings
                .iter()
                .zip(live)
                .find(|(_, on)| *on)
                .map(|(winding, _)| winding.bus)?;
            if live.iter().filter(|&&on| on).count() < 2 {
                return None;
            }
            let mut windings = t.windings.clone();
            for w in &mut windings {
                w.bus = remap(map, w.bus).unwrap_or(kept_bus);
                if let Some(control) = &mut w.control {
                    norm_transformer_control(control, base, map);
                }
                w.shift *= DEG_TO_RAD;
                w.rate_a /= base;
                w.rate_b /= base;
                w.rate_c /= base;
            }
            let mut transformer = Transformer3W {
                windings,
                star_va: t.star_va * DEG_TO_RAD,
                ..t.clone()
            };
            for winding in 0..3 {
                let mut sets = t.winding_rating_sets(winding);
                for set in &mut sets {
                    set.rate_mva /= base;
                }
                transformer.set_winding_rating_sets(winding, sets);
            }
            if live != t.winding_in_service() {
                transformer.set_winding_in_service(live);
                taken_out.push(row);
            }
            Some((transformer, Some(row)))
        })
        .unzip();
    if !taken_out.is_empty() {
        let names = taken_out
            .iter()
            .take(5)
            .map(|&row| {
                xfmrs[row]
                    .uid
                    .clone()
                    .unwrap_or_else(|| format!("transformers_3w:{row}"))
            })
            .collect::<Vec<_>>()
            .join(", ");
        warnings.push(
            &crate::diagnostics::codes::CANONICALIZE_NORMALIZE_WINDING_TAKEN_OUT,
            format!(
                "{} three winding transformer(s) have a winding on a bus normalization drops; \
                 that winding was taken out of service and the other two kept ({names})",
                taken_out.len()
            ),
        );
    }
    normalized
}

/// No reference survived the bus type pass: anchor the slack at the largest
/// pmax in-service generator's bus and record the designation on the coded
/// channel, or refuse when there is no generator to anchor it.
fn designate_reference(
    buses: &mut [Bus],
    generators: &[Generator],
    warnings: &mut crate::diagnostics::Diagnostics,
) -> Result<()> {
    let slack = generators
        .iter()
        .max_by(|a, b| {
            // A NaN pmax must never win the slack: map it below every real
            // bound so the choice stays deterministic (an unbounded +Inf
            // pmax still wins, as the largest capacity).
            let key = |p: f64| if p.is_nan() { f64::NEG_INFINITY } else { p };
            key(a.pmax).total_cmp(&key(b.pmax))
        })
        .map(|g| g.bus)
        .ok_or(Error::NoReferenceBus)?;
    if let Some(b) = buses.iter_mut().find(|b| b.id == slack) {
        b.kind = BusType::Ref;
        warnings.push(
            &crate::diagnostics::codes::CANONICALIZE_NORMALIZE_REFERENCE_DESIGNATED,
            format!(
                "the case states no reference bus that survives normalization; bus {slack} \
                 hosts the largest pmax in-service generator and was designated the slack"
            ),
        );
    }
    Ok(())
}

/// The buses of the islands no in-service generator supplies, each island
/// reported as left out.
fn unsupplied_island_buses(
    network: &BalancedNetwork,
    partition: &IslandPartition,
    warnings: &mut crate::diagnostics::Diagnostics,
) -> HashSet<BusId> {
    let mut unsupplied = HashSet::new();
    for island in partition
        .islands
        .iter()
        .filter(|island| !island.is_supplied())
    {
        warnings.record(
            network.unsupplied_island_finding(island, "was left out of the normalized network"),
        );
        unsupplied.extend(island.buses.iter().copied());
    }
    unsupplied
}

/// Canonical bus types: a bus hosting an in-service generator keeps `Ref` if
/// the file marked it `Ref`, else becomes `Pv`; a gen-less bus is `Pq`. A
/// `Pq` bus that becomes `Pv` is reported, since its generator now holds the
/// bus voltage.
fn canonicalize_bus_types(
    buses: &mut [Bus],
    generators: &[Generator],
    warnings: &mut crate::diagnostics::Diagnostics,
) {
    let gen_buses: HashSet<BusId> = generators.iter().map(|g| g.bus).collect();
    let mut retyped_pq = Vec::new();
    for b in buses.iter_mut() {
        if b.kind == BusType::Pq && gen_buses.contains(&b.id) {
            retyped_pq.push(b.id);
        }
        b.kind = match (gen_buses.contains(&b.id), b.kind) {
            (true, BusType::Ref) => BusType::Ref,
            (true, _) => BusType::Pv,
            (false, _) => BusType::Pq,
        };
    }
    if !retyped_pq.is_empty() {
        warnings.push(
            &crate::diagnostics::codes::CANONICALIZE_NORMALIZE_GENERATOR_BUS_RETYPED,
            format!(
                "{} PQ bus(es) host an in-service generator and were retyped PV, so the \
                 generator holds the bus voltage ({})",
                retyped_pq.len(),
                bus_list(&retyped_pq)
            ),
        );
    }
}

/// The reference buses of the normalized tables: one per island of
/// `partition` under [`IslandReferencePolicy::PerIsland`], and the largest
/// generator's bus when no reference survives at all.
fn assign_normalized_references(
    buses: &mut [Bus],
    generators: &[Generator],
    partition: Option<&IslandPartition>,
    warnings: &mut crate::diagnostics::Diagnostics,
) -> Result<()> {
    if let Some(partition) = partition {
        assign_normalized_island_references(buses, generators, partition, warnings);
    }
    if !buses.iter().any(|b| b.kind == BusType::Ref) {
        designate_reference(buses, generators, warnings)?;
    }
    Ok(())
}

/// [`IslandReferencePolicy::PerIsland`] over the normalized tables: each
/// supplied island of `partition` keeps exactly one reference bus. The bus
/// types are already canonical, so a reference hosts a generator and a demoted
/// one becomes PV. This picks the same buses as
/// [`BalancedNetwork::assign_island_references`] on the raw network.
fn assign_normalized_island_references(
    buses: &mut [Bus],
    generators: &[Generator],
    partition: &IslandPartition,
    warnings: &mut crate::diagnostics::Diagnostics,
) {
    let mut island_generators: HashMap<usize, Vec<usize>> = HashMap::new();
    for (row, generator) in generators.iter().enumerate() {
        if let Some(island) = partition.island_of(generator.bus) {
            island_generators.entry(island).or_default().push(row);
        }
    }
    let mut island_references: HashMap<usize, Vec<BusId>> = HashMap::new();
    for bus in buses.iter().filter(|bus| bus.kind == BusType::Ref) {
        if let Some(island) = partition.island_of(bus.id) {
            island_references.entry(island).or_default().push(bus.id);
        }
    }
    let mut designated = HashSet::new();
    let mut demoted = HashSet::new();
    for (position, island) in partition.islands.iter().enumerate() {
        let Some(rows) = island_generators.get(&position) else {
            continue;
        };
        let mut references = island_references.remove(&position).unwrap_or_default();
        references.sort_unstable();
        match references.as_slice() {
            [_] => {}
            [] => {
                if let Some(bus) = largest_generator_bus(generators, rows) {
                    warnings.push(
                        &crate::diagnostics::codes::CANONICALIZE_ISLAND_REFERENCE_DESIGNATED,
                        format!(
                            "the island of {} bus(es) ({}) has no reference bus that survives \
                             normalization; bus {bus} hosts its largest pmax in-service \
                             generator and was designated the reference",
                            island.buses.len(),
                            bus_list(&island.buses)
                        ),
                    );
                    designated.insert(bus);
                }
            }
            many => {
                let kept = kept_reference(many, generators, rows).unwrap_or(many[0]);
                let others: Vec<BusId> = many.iter().copied().filter(|&b| b != kept).collect();
                warnings.push(
                    &crate::diagnostics::codes::CANONICALIZE_ISLAND_REFERENCE_DEMOTED,
                    format!(
                        "the island of {} bus(es) keeps {} reference buses; bus {kept} hosts the \
                         most generation and stays the reference, and bus(es) {} became PV",
                        island.buses.len(),
                        many.len(),
                        bus_list(&others)
                    ),
                );
                demoted.extend(others);
            }
        }
    }
    for bus in buses.iter_mut() {
        if designated.contains(&bus.id) {
            bus.kind = BusType::Ref;
        } else if demoted.contains(&bus.id) {
            bus.kind = BusType::Pv;
        }
    }
}

impl BalancedNetwork {
    /// A normalized, computation-ready copy of this network. The raw `BalancedNetwork` is
    /// kept lossless (MATPOWER units, 1-based sparse ids, out-of-service elements
    /// retained); `to_normalized` derives the form a solver or ML pipeline wants:
    ///
    /// - **Per unit** (÷`base_mva`): gen `pg/qg/pmax/pmin/qmax/qmin` and the ramp
    ///   caps (`GEN_PU_KEYS`); load `p/q`; shunt `g/b`; branch `rate_a/b/c`;
    ///   storage energy/ratings/limits/losses; HVDC `pf/pt/qf/qt`, reactive limits,
    ///   `loss0`; gen-cost coefficients (`cost_to_pu`). Storage `ps/qs` and HVDC
    ///   aggregate `pmin/pmax` stay raw, matching the PowerModels per-unit
    ///   convention. Voltages, impedances, tap, and `loss1` are already
    ///   dimensionless.
    /// - **Radians**: bus `va`; branch `shift/angmin/angmax`.
    /// - **Tap**: `0 → 1.0` (an explicit `1` is kept).
    /// - **Filtered**: drop buses typed isolated (`BusType::Isolated`) and every
    ///   out-of-service element, then drop any element left referencing a dropped
    ///   bus. A bus orphaned by the out-of-service filter (no in-service branch,
    ///   but not typed isolated) is kept — its load is real — and surfaces as its
    ///   own island, which the grounding check reports if it has no reference.
    /// - **Closed switches**: the buses a closed switch joins merge into one
    ///   ([`ClosedSwitchPolicy::Merge`], reported under `CANONICALIZE.MERGE`),
    ///   whose survivor is the reference bus, else a generator bus, else a
    ///   regulated bus, smallest id first. The matrix builders model no
    ///   switches, so this keeps their buses together;
    ///   [`BalancedNetwork::to_normalized_with`] keeps or refuses them instead.
    /// - **IDs**: kept buses retain their source bus ids, and every surviving
    ///   endpoint stays in the same id space. Consumers that need dense rows should
    ///   use [`IndexedNetwork`](crate::IndexedNetwork), which derives `[0, n)`
    ///   indices without destroying source ids.
    /// - **Bus types**: a bus hosting a surviving generator keeps `REF` if the file
    ///   marked it `REF`, otherwise becomes `PV`; a generator-less bus is `PQ` (so a
    ///   generator-less `REF` is demoted). A `PQ` bus retyped `PV` is reported as
    ///   `CANONICALIZE.NORMALIZE.GENERATOR_BUS_RETYPED`. The file's `REF` buses are
    ///   kept, several included, and the consumer picks the slack. Only when no
    ///   reference bus survives is the largest-`pmax` in-service generator's bus
    ///   promoted to `REF`;
    ///   [`NormalizeSettings::island_references`] applies that rule per island
    ///   instead.
    ///
    /// This is a derived product, not a source for write-back: `source` is dropped
    /// and `source_format` is [`SourceFormat::Normalized`], so writing it serializes
    /// the per-unit/radian model instead of echoing the raw bytes, and a consumer
    /// can tell it apart from a raw in-memory network.
    ///
    /// Scope is the universal canonicalization only. It does not synthesize a
    /// missing `rate_a` or restrict the gen-cost model — those are solver
    /// preparation choices a consumer applies on top. Use
    /// [`BalancedNetwork::to_normalized_with_options`] for the opt in PowerModels angle
    /// bound repair. The cost *rescale* is
    /// universal and lives here; the model *restriction* does not.
    ///
    /// # Errors
    /// [`Error::InvalidBaseMva`] if `base_mva` is not a positive, finite number
    /// (every per-unit divisor), so a malformed base can't silently poison the
    /// whole network with `NaN`/`Inf` or sign-flipped values.
    /// [`Error::NoReferenceBus`] if no reference bus can be established — no `REF`
    /// survives and there is no in-service generator to anchor one.
    /// [`Error::ClosedSwitch`] if a requested switch merge cannot complete,
    /// such as a connection between windings of one three winding transformer.
    pub fn to_normalized(&self) -> Result<BalancedNetwork> {
        Ok(self
            .to_normalized_with_options(&NormalizeOptions::default())?
            .network)
    }

    /// Like [`BalancedNetwork::to_normalized`], with opt in solver preparation repairs
    /// that report fidelity warnings.
    pub fn to_normalized_with_options(
        &self,
        options: &NormalizeOptions,
    ) -> Result<NormalizedNetwork> {
        Ok(self.normalize_inner(&NormalizeSettings::new(*options))?.0)
    }

    /// Like [`BalancedNetwork::to_normalized_with_options`], under every
    /// choice [`NormalizeSettings`] states, also returning what the pass
    /// changed before normalizing ([`NormalizeReport`]), such as the closed
    /// switch merge.
    ///
    /// # Errors
    /// As [`BalancedNetwork::to_normalized`], and [`Error::ClosedSwitch`]
    /// under [`ClosedSwitchPolicy::Refuse`] or an incomplete
    /// [`ClosedSwitchPolicy::Merge`].
    pub fn to_normalized_with(
        &self,
        settings: &NormalizeSettings,
    ) -> Result<(NormalizedNetwork, NormalizeReport)> {
        let (normalized, _, report) = self.normalize_inner(settings)?;
        Ok((normalized, report))
    }

    /// Like [`BalancedNetwork::to_normalized_with_options`], also returning the
    /// [`NormalizeSourceRows`] row provenance.
    ///
    /// The rows are positional over the
    /// [`IndexedNetwork`](crate::IndexedNetwork) view of the returned network,
    /// which is the index space a matrix row or a solver table row lives in.
    /// That view star-lowers a 3-winding transformer, so it holds more buses,
    /// branches, and shunts than the returned [`NormalizedNetwork`] does;
    /// indexing the returned network by a row position is out of bounds on any
    /// case that carries one. Resolve a row through the view:
    ///
    /// ```
    /// # use powerio_tx::{IndexedNetwork, BalancedNetwork, NormalizeOptions};
    /// # fn f(raw: &BalancedNetwork) -> powerio_tx::Result<()> {
    /// let (normalized, rows) = raw.to_normalized_with_source_rows(&NormalizeOptions::default())?;
    /// let view = IndexedNetwork::new(&normalized.network);
    /// for (dense, source) in rows.buses.iter().enumerate() {
    ///     let bus = &view.network().buses()[dense];
    ///     // `source` is `None` for the synthetic star bus the view appended.
    ///     let _ = (bus, source);
    /// }
    /// # Ok(())
    /// # }
    /// ```
    #[doc(hidden)]
    pub fn to_normalized_with_source_rows(
        &self,
        options: &NormalizeOptions,
    ) -> Result<(NormalizedNetwork, NormalizeSourceRows)> {
        let (normalized, mut rows, _) = self.normalize_inner(&NormalizeSettings::new(*options))?;
        rows.pad_to_lowered(&normalized.network);
        Ok((normalized, rows))
    }

    /// The pass itself, after what `settings` asks of closed switches. The
    /// rows it gives cover the normalized network before the star lowering,
    /// so only [`Self::to_normalized_with_source_rows`] pays for the lowered
    /// lengths.
    fn normalize_inner(
        &self,
        settings: &NormalizeSettings,
    ) -> Result<(NormalizedNetwork, NormalizeSourceRows, NormalizeReport)> {
        let options = &settings.options;
        validate_normalize_options(options)?;
        let (kept, left) = self.closed_switch_rows();
        if kept.is_empty() || settings.closed_switches == ClosedSwitchPolicy::Keep {
            let (normalized, rows) = self.normalize_pass(settings)?;
            return Ok((normalized, rows, NormalizeReport::default()));
        }
        if settings.closed_switches == ClosedSwitchPolicy::Refuse {
            let switch = &self.switches()[kept[0]];
            return Err(Error::ClosedSwitch {
                row: kept[0],
                from: switch.from,
                to: switch.to,
                closed: kept.len(),
            });
        }
        let merge = self.merge_normalization_switches(&left)?;
        // A winding-pair refusal leaves the selected switch in the result.
        // Report the source row, even when earlier switches were removed or
        // their endpoints were remapped by a successful partial merge.
        let mut remaining = kept
            .into_iter()
            .filter(|&row| merge.switch_rows[row].is_some());
        if let Some(row) = remaining.next() {
            let switch = &self.switches()[row];
            return Err(Error::ClosedSwitch {
                row,
                from: switch.from,
                to: switch.to,
                closed: 1 + remaining.count(),
            });
        }
        let (mut normalized, rows) = merge.network.normalize_pass(settings)?;
        let rows = rows.through_merge(&merge);
        let mut diagnostics = merge.diagnostics.clone();
        diagnostics.append(&mut normalized.diagnostics);
        normalized.warnings = crate::diagnostics::render_diagnostics(&diagnostics);
        normalized.diagnostics = diagnostics;
        Ok((
            normalized,
            rows,
            NormalizeReport {
                bus_merge: Some(merge),
            },
        ))
    }

    /// Merge only switches whose endpoints survive bus filtering, then
    /// compose the merge's switch provenance back to the original table.
    fn merge_normalization_switches(&self, left: &[usize]) -> Result<BusMerge> {
        if left.is_empty() {
            return self.merge_buses(&BusMergeRule::closed_switches());
        }
        let mut eligible = self.clone();
        let excluded: HashSet<_> = left.iter().copied().collect();
        let mut source_rows = Vec::new();
        let mut row = 0;
        eligible.switches_mut().retain_mut(|switch| {
            let keep = !excluded.contains(&row);
            if keep {
                source_rows.push(row);
                // Keep anonymous source identities in merge diagnostics even
                // though filtering changes this temporary table's row numbers.
                if switch.uid.is_none() {
                    switch.uid = Some(format!("switches:{row}"));
                }
            }
            row += 1;
            keep
        });
        // Opening an ineligible switch is insufficient: an undeclared bus
        // still leaves a dangling reference for merge_buses to validate.
        let mut merge = eligible.merge_buses(&BusMergeRule::closed_switches())?;
        let mut switch_rows = vec![None; self.switches().len()];
        for (&source, &merged) in source_rows.iter().zip(&merge.switch_rows) {
            switch_rows[source] = merged;
            if let Some(row) = merged {
                merge.network.switches_mut()[row]
                    .uid
                    .clone_from(&self.switches()[source].uid);
            }
        }
        for removed in &mut merge.removed_switches {
            removed.row = source_rows[removed.row];
        }
        merge.switch_rows = switch_rows;
        merge.source = self.clone();
        Ok(merge)
    }

    /// Closed switches joining distinct buses the pass keeps, and all
    /// switches on isolated or undeclared buses, which the pass drops.
    fn closed_switch_rows(&self) -> (Vec<usize>, Vec<usize>) {
        let kept_bus: HashMap<BusId, bool> = self
            .buses()
            .iter()
            .map(|bus| (bus.id, bus.kind != BusType::Isolated))
            .collect();
        let kept = |bus: BusId| kept_bus.get(&bus).copied().unwrap_or(false);
        self.switches().iter().enumerate().fold(
            (Vec::new(), Vec::new()),
            |(mut eligible, mut left), (row, switch)| {
                if !kept(switch.from) || !kept(switch.to) {
                    left.push(row);
                } else if switch.closed && switch.from != switch.to {
                    eligible.push(row);
                }
                (eligible, left)
            },
        )
    }

    /// One normalization pass over this network as it stands.
    fn normalize_pass(
        &self,
        settings: &NormalizeSettings,
    ) -> Result<(NormalizedNetwork, NormalizeSourceRows)> {
        let options = &settings.options;
        validate_normalize_options(options)?;
        self.check_base_mva()?;
        let base = self.base_mva();
        let mut warnings = crate::diagnostics::Diagnostics::new();

        // Per island references: an island no in-service generator supplies
        // cannot be solved, so its buses stay out like isolated ones and the
        // element filters below drop what sits on them. The partition follows
        // what this pass keeps: a three winding transformer keeps its live
        // windings when another winding sits on an isolated bus.
        let partition = (settings.island_references == IslandReferencePolicy::PerIsland)
            .then(|| self.calc_islands());
        let unsupplied = partition
            .as_ref()
            .map(|partition| unsupplied_island_buses(self, partition, &mut warnings))
            .unwrap_or_default();

        // Kept buses keep their original `kind` for now (the reference scan below
        // reads it) and their source ids. Isolated buses are dropped.
        let mut id_map: HashMap<BusId, BusId> = HashMap::with_capacity(self.buses().len());
        let mut buses: Vec<Bus> = Vec::with_capacity(self.buses().len());
        // The pass keeps only elements that came from a source row, so each row
        // here is `Some`; the `None` entries appear later, when
        // `pad_to_lowered` extends the map over what the star lowering appends.
        let mut bus_rows: Vec<Option<usize>> = Vec::with_capacity(self.buses().len());
        for (row, b) in self.buses().iter().enumerate() {
            if b.kind == BusType::Isolated || unsupplied.contains(&b.id) {
                continue;
            }
            id_map.insert(b.id, b.id);
            buses.push(Bus {
                va: b.va * DEG_TO_RAD,
                ..b.clone()
            });
            bus_rows.push(Some(row));
        }
        let (loads, load_rows) = norm_loads(self.loads(), base, &id_map);
        let (shunts, shunt_rows) = norm_shunts(self.shunts(), base, &id_map);
        let (static_var_compensators, static_var_compensator_rows) =
            norm_static_var_compensators(self.static_var_compensators(), base, &id_map);
        let (mut branches, branch_rows) = norm_branches(self.branches(), base, &id_map);
        if options.clamp_angle_bounds {
            clamp_angle_bounds(&mut branches, options.angle_bound_pad, &mut warnings);
        }
        let (switches, switch_rows) = norm_switches(self.switches(), base, &id_map);
        let (generators, generator_rows) = norm_gens(self.generators(), base, &id_map);
        let (storage, storage_rows) = norm_storage(self.storage(), base, &id_map);
        let (hvdc, hvdc_rows) = norm_hvdc(self.hvdc(), base, &id_map);
        let (transformers_3w, transformer_3w_rows) =
            norm_transformers_3w(self.transformers_3w(), base, &id_map, &mut warnings);
        let source_rows = NormalizeSourceRows {
            buses: bus_rows,
            loads: load_rows,
            shunts: shunt_rows,
            static_var_compensators: static_var_compensator_rows,
            branches: branch_rows,
            switches: switch_rows,
            generators: generator_rows,
            storage: storage_rows,
            hvdc: hvdc_rows,
            transformers_3w: transformer_3w_rows,
        };

        // Bus types. Multiple file `Ref` buses are kept as-is unless the per
        // island policy keeps one per island, and only when no `Ref` survives
        // is the largest-pmax generator's bus promoted.
        canonicalize_bus_types(&mut buses, &generators, &mut warnings);
        assign_normalized_references(&mut buses, &generators, partition.as_ref(), &mut warnings)?;
        // The other silent semantic decision this gateway announces: a
        // solver-ready copy whose cost objective is identically zero.
        if !generators.is_empty() && generators.iter().all(|g| g.cost.is_none()) {
            warnings.push(
                &crate::diagnostics::codes::CANONICALIZE_NORMALIZE_GEN_COST_ABSENT,
                format!(
                    "the case has {} in-service generator(s) and no cost data; any cost \
                     objective built from it is identically zero",
                    generators.len()
                ),
            );
        }

        let mut net = BalancedNetwork::from_tables(BalancedNetworkTables {
            name: self.name().clone(),
            base_mva: base,
            base_frequency: self.base_frequency(),
            geo: self.geo().clone(),
            case_metadata: self.case_metadata().clone(),
            detailed_connectivity: self.detailed_connectivity().clone(),
            generated_uids: self.generated_uids().clone(),
            buses: buses.into(),
            loads: loads.into(),
            shunts: shunts.into(),
            static_var_compensators: static_var_compensators.into(),
            branches: branches.into(),
            switches: switches.into(),
            generators: generators.into(),
            storage: storage.into(),
            hvdc: hvdc.into(),
            transformers_3w: transformers_3w.into(),
            // Areas (interchange schedule, per-area swing) are interchange metadata,
            // not part of the per unit electrical view, so they are not carried.
            areas: Vec::new().into(),
            solver: None,
            source_format: SourceFormat::Normalized,
        });
        // Project source topology after all electrical filters have run, including
        // unsupplied islands, isolated buses, and out-of-service equipment.
        crate::operations::filter_connectivity(self, &mut net);
        // The filter drops every reference to a dropped bus by
        // construction, so the result is reference-consistent. Assert it in
        // debug builds to catch a future regression in the filtering logic.
        debug_assert!(
            net.validate().is_ok(),
            "to_normalized produced a dangling reference"
        );
        Ok((
            NormalizedNetwork {
                network: net,
                warnings: warnings.lines(),
                diagnostics: warnings.into_records(),
            },
            source_rows,
        ))
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn normalization_merges_closed_switches_unless_told_otherwise() {
        use crate::network::Switch;
        let mut buses: Vec<Bus> = (1..=4)
            .map(|id| Bus::new(BusId(id), BusType::Pq, 230.0))
            .collect();
        buses[0].kind = BusType::Ref;
        buses[3].kind = BusType::Isolated;
        let mut net = BalancedNetwork::in_memory(
            "switched",
            100.0,
            buses,
            vec![Branch::new(BusId(1), BusId(2), 0.0, 0.1)],
        );
        // A closed switch 2-3, an open one 1-3, and a closed one onto the
        // isolated bus 4, which the pass drops with that bus.
        net.switches_mut().extend([
            Switch::new(BusId(2), BusId(3), true),
            Switch::new(BusId(1), BusId(3), false),
            Switch::new(BusId(2), BusId(4), true),
        ]);
        net.generators_mut().push(Generator::new(BusId(1)));
        net.loads_mut().push(Load::new(BusId(3), 10.0, 0.0));

        // The default merges, reports the merge, and maps rows to the source.
        let (merged, rows) = net
            .to_normalized_with_source_rows(&NormalizeOptions::default())
            .unwrap();
        assert_eq!(merged.network.buses().len(), 2);
        assert_eq!(merged.network.loads()[0].bus, BusId(2));
        crate::IndexedNetwork::new(&merged.network)
            .check_closed_switches()
            .unwrap();
        assert!(
            merged
                .diagnostics
                .iter()
                .any(|d| d.code() == "CANONICALIZE.MERGE.CLOSED_SWITCH")
        );
        // The open switch survives at merged row 0 and maps back to source row 1.
        assert_eq!(rows.switches, [Some(1)]);
        assert_eq!(rows.buses, [Some(0), Some(1)]);
        assert_eq!(rows.branches, [Some(0)]);
        let plain = net.to_normalized().unwrap();
        assert_eq!(plain.buses(), merged.network.buses());
        assert_eq!(plain.switches(), merged.network.switches());

        let (_, report) = net
            .to_normalized_with(&NormalizeSettings::default())
            .unwrap();
        let merge = report.bus_merge.unwrap();
        assert_eq!(merge.survivor(BusId(3)), BusId(2));
        assert_eq!(merge.survivor(BusId(4)), BusId(4));

        // Keep leaves the switches as data, as normalization used to.
        let (kept, report) = net
            .to_normalized_with(
                &NormalizeSettings::default().with_closed_switches(ClosedSwitchPolicy::Keep),
            )
            .unwrap();
        assert!(report.bus_merge.is_none());
        assert_eq!(kept.network.buses().len(), 3);
        let error = crate::IndexedNetwork::new(&kept.network)
            .check_closed_switches()
            .unwrap_err();
        assert!(error.to_string().starts_with("1 closed switch(es)"));

        // Refuse states the switch the pass would keep.
        let error = net
            .to_normalized_with(
                &NormalizeSettings::default().with_closed_switches(ClosedSwitchPolicy::Refuse),
            )
            .unwrap_err();
        assert!(matches!(
            error,
            Error::ClosedSwitch {
                row: 0,
                closed: 1,
                ..
            }
        ));
        let message = error.to_string();
        for call in [
            "BusMergeRule::closed_switches()",
            "closed_switches=\"merge\"",
            "--merge-buses switches",
        ] {
            assert!(message.contains(call), "{message}");
        }

        // A network with no closed switch normalizes as it always has.
        net.switches_mut()[0].closed = false;
        net.switches_mut()[2].closed = false;
        let (plain, report) = net
            .to_normalized_with(&NormalizeSettings::default())
            .unwrap();
        assert!(report.bus_merge.is_none());
        assert_eq!(plain.network.buses().len(), 3);
    }
    use crate::network::GeneratorEnergySource;

    fn approx(a: f64, b: f64) -> bool {
        (a - b).abs() < 1e-9
    }

    fn angle_bound_fixture() -> BalancedNetwork {
        let path = std::path::PathBuf::from(env!("CARGO_MANIFEST_DIR"))
            .join("../tests/data/angle_bounds_clamp.m");
        crate::parse_file(path, None).unwrap().network
    }

    #[test]
    fn transformer_control_normalization_follows_psse_field_units_for_every_mode() {
        // Siemens PSS/E 35.4.1 defines RMA/RMI as ratios for |COD| 1/2,
        // degrees for |COD| 3/5, and unused for |COD| 0/4. VMA/VMI are p.u.
        // voltage for |COD| 1, Mvar for |COD| 2, MW for |COD| 3/5, and unused
        // for |COD| 0/4. Preserve unused values rather than inventing units.
        let map = HashMap::from([(BusId(7), BusId(7))]);
        let cases = [
            (TransformerControlMode::Fixed, false, false),
            (TransformerControlMode::Voltage, false, false),
            (TransformerControlMode::ReactiveFlow, false, true),
            (TransformerControlMode::ActiveFlow, true, true),
            (TransformerControlMode::DcLineQuantity, false, false),
            (TransformerControlMode::AsymmetricActiveFlow, true, true),
        ];

        for (mode, angle_limits, power_band) in cases {
            let mut control = TransformerControl::new(mode);
            control.controlled_bus = Some(BusId(7));
            control.tap_min = -10.0;
            control.tap_max = 20.0;
            control.band_min = -50.0;
            control.band_max = 75.0;
            control.winding_connection_angle =
                (mode == TransformerControlMode::AsymmetricActiveFlow).then_some(30.0);

            norm_transformer_control(&mut control, 100.0, &map);

            let tap_scale = if angle_limits { DEG_TO_RAD } else { 1.0 };
            let band_scale = if power_band { 0.01 } else { 1.0 };
            assert!(approx(control.tap_min, -10.0 * tap_scale), "{mode:?}");
            assert!(approx(control.tap_max, 20.0 * tap_scale), "{mode:?}");
            assert!(approx(control.band_min, -50.0 * band_scale), "{mode:?}");
            assert!(approx(control.band_max, 75.0 * band_scale), "{mode:?}");
            assert_eq!(control.controlled_bus, Some(BusId(7)), "{mode:?}");
            if mode == TransformerControlMode::AsymmetricActiveFlow {
                assert!(approx(
                    control.winding_connection_angle.unwrap(),
                    30.0 * DEG_TO_RAD
                ));
            } else {
                assert_eq!(control.winding_connection_angle, None, "{mode:?}");
            }
        }
    }

    #[test]
    fn angle_bound_clamp_is_opt_in_and_matches_powermodels_rules() {
        let net = angle_bound_fixture();

        let plain = net.to_normalized().unwrap();
        assert!(approx(plain.branches()[0].angmin, -std::f64::consts::TAU));
        assert!(approx(plain.branches()[0].angmax, std::f64::consts::TAU));
        assert!(approx(plain.branches()[1].angmin, 0.0));
        assert!(approx(plain.branches()[1].angmax, 0.0));
        assert!(approx(plain.branches()[3].angmin, -120.0 * DEG_TO_RAD));
        assert!(approx(plain.branches()[3].angmax, -100.0 * DEG_TO_RAD));
        assert!(approx(plain.branches()[4].angmin, 100.0 * DEG_TO_RAD));
        assert!(approx(plain.branches()[4].angmax, 120.0 * DEG_TO_RAD));

        let out = net
            .to_normalized_with_options(&NormalizeOptions {
                clamp_angle_bounds: true,
                ..NormalizeOptions::default()
            })
            .unwrap();
        // The fixture also carries no gencost, so the costless-case warning
        // rides beside the clamp lines; hold the clamp set on its own code.
        let clamps: Vec<&String> = out
            .warnings
            .iter()
            .filter(|w| w.contains("BOUNDS_CLAMPED"))
            .collect();
        assert_eq!(clamps.len(), 4, "{:?}", out.warnings);
        assert!(clamps[0].contains("branch 0"));
        assert!(clamps[1].contains("branch 1"));
        assert!(clamps[2].contains("branch 3"));
        assert!(clamps[3].contains("branch 4"));

        let branches = &out.network.branches();
        assert!(approx(branches[0].angmin, -POWER_MODELS_ANGLE_BOUND_PAD));
        assert!(approx(branches[0].angmax, POWER_MODELS_ANGLE_BOUND_PAD));
        assert!(approx(branches[1].angmin, -POWER_MODELS_ANGLE_BOUND_PAD));
        assert!(approx(branches[1].angmax, POWER_MODELS_ANGLE_BOUND_PAD));
        assert!(approx(branches[2].angmin, -30.0 * DEG_TO_RAD));
        assert!(approx(branches[2].angmax, 30.0 * DEG_TO_RAD));
        assert!(approx(branches[3].angmin, -POWER_MODELS_ANGLE_BOUND_PAD));
        assert!(approx(branches[3].angmax, POWER_MODELS_ANGLE_BOUND_PAD));
        assert!(approx(branches[4].angmin, -POWER_MODELS_ANGLE_BOUND_PAD));
        assert!(approx(branches[4].angmax, POWER_MODELS_ANGLE_BOUND_PAD));
        assert!(branches.iter().all(|br| br.angmin <= br.angmax));
    }

    #[test]
    fn angle_bound_clamp_rejects_invalid_pad() {
        let net = angle_bound_fixture();
        let err = net
            .to_normalized_with_options(&NormalizeOptions {
                clamp_angle_bounds: true,
                angle_bound_pad: std::f64::consts::FRAC_PI_2,
            })
            .unwrap_err();
        assert!(matches!(
            err,
            Error::InvalidNormalizeOption {
                field: "angle_bound_pad",
                ..
            }
        ));
    }

    #[test]
    fn to_normalized_drops_a_control_bus_whose_target_was_filtered_out() {
        use crate::network::{Extras, ShuntBlock, SwitchedShuntControl, SwitchedShuntMode};

        let mkbus = |id: usize, kind: BusType| Bus {
            id: BusId(id),
            kind,
            vm: 1.0,
            va: 0.0,
            base_kv: 230.0,
            vmax: 1.1,
            vmin: 0.9,
            evhi: None,
            evlo: None,
            area: 1,
            zone: 1,
            name: None,
            uid: None,
            location: None,
            extras: Extras::new(),
        };
        let branch = Branch {
            name: None,
            from: BusId(1),
            to: BusId(2),
            r: 0.0,
            x: 0.1,
            b: 0.0,
            charging: None,
            rate_a: 0.0,
            rate_b: 0.0,
            rate_c: 0.0,
            rating_sets: Vec::new(),
            current_ratings: None,
            tap: 0.0,
            shift: 0.0,
            in_service: true,
            angmin: -360.0,
            angmax: 360.0,
            control: None,
            solution: None,
            uid: None,
            route: None,
            extras: Extras::new(),
        };
        // Bus 3 is isolated, so to_normalized drops it.
        let mut net = BalancedNetwork::in_memory(
            "n",
            100.0,
            vec![
                mkbus(1, BusType::Ref),
                mkbus(2, BusType::Pq),
                mkbus(3, BusType::Isolated),
            ],
            vec![branch],
        );
        net.generators_mut().push(Generator {
            bus: BusId(1),
            energy_source: GeneratorEnergySource::default(),
            pg: 10.0,
            qg: 0.0,
            pmax: 100.0,
            pmin: 0.0,
            qmax: 50.0,
            qmin: -50.0,
            vg: 1.0,
            mbase: 100.0,
            in_service: true,
            cost: None,
            caps: Default::default(),
            voltage_regulation_on: true,
            regulating_terminal: None,
            regulated_bus: None,
            active_power_control: None,
            uid: None,
        });
        // A switched shunt on bus 2 whose control bus is the (dropped) isolated bus 3.
        net.shunts_mut().push(Shunt {
            bus: BusId(2),
            g: 0.0,
            b: 10.0,
            in_service: true,
            section_count: None,
            control: Some(SwitchedShuntControl {
                mode: SwitchedShuntMode::Discrete,
                vhigh: 1.05,
                vlow: 0.95,
                control_bus: Some(BusId(3)),
                regulating_terminal: None,
                rmpct: 100.0,
                blocks: vec![ShuntBlock::with_admittance(2, 4.0, 20.0)],
            }),
            uid: None,
            extras: Extras::new(),
        });

        let norm = net.to_normalized().unwrap();
        norm.validate().unwrap();
        let c = norm.shunts()[0].control.as_ref().expect("control retained");
        assert_eq!(
            c.control_bus, None,
            "a control bus pointing at a filtered-out isolated bus is dropped, not left dangling"
        );
        assert!(approx(c.blocks[0].g, 0.04));
        assert!(approx(c.blocks[0].b, 0.2));
    }

    #[test]
    fn normalized_slack_tiebreak_ignores_nan_pmax() {
        use crate::network::Extras;

        let mkbus = |id: usize| Bus {
            id: BusId(id),
            kind: BusType::Pq,
            vm: 1.0,
            va: 0.0,
            base_kv: 230.0,
            vmax: 1.1,
            vmin: 0.9,
            evhi: None,
            evlo: None,
            area: 1,
            zone: 1,
            name: None,
            uid: None,
            location: None,
            extras: Extras::new(),
        };
        let mkgen = |bus: usize, pmax: f64| Generator {
            bus: BusId(bus),
            energy_source: GeneratorEnergySource::default(),
            pg: 0.0,
            qg: 0.0,
            pmax,
            pmin: 0.0,
            qmax: 0.0,
            qmin: 0.0,
            vg: 1.0,
            mbase: 100.0,
            in_service: true,
            cost: None,
            caps: Default::default(),
            voltage_regulation_on: true,
            regulating_terminal: None,
            regulated_bus: None,
            active_power_control: None,
            uid: None,
        };
        let mut net = BalancedNetwork::in_memory("n", 100.0, vec![mkbus(1), mkbus(2)], Vec::new());
        *net.generators_mut() = vec![mkgen(1, f64::NAN), mkgen(2, 10.0)];
        let norm = net.to_normalized().unwrap();

        assert_eq!(
            norm.buses().iter().find(|b| b.id == BusId(1)).unwrap().kind,
            BusType::Pv
        );
        assert_eq!(
            norm.buses().iter().find(|b| b.id == BusId(2)).unwrap().kind,
            BusType::Ref
        );
    }

    #[test]
    fn cost_to_pu_polynomial_scales_and_trims() {
        // Model 2: the coeff of p^j scales by base^j; MATPOWER's trailing-zero
        // padding (beyond ncost) is dropped.
        let cost = GenCost {
            model: 2,
            startup: 0.0,
            shutdown: 0.0,
            ncost: 2,
            coeffs: vec![24.035, -403.5, 0.0, 0.0, 0.0, 0.0],
        };
        let out = cost_to_pu(&cost, 100.0);
        assert_eq!(out.len(), 2, "padding dropped");
        assert!(approx(out[0], 2403.5)); // 24.035 · 100^1
        assert!(approx(out[1], -403.5)); // -403.5 · 100^0
    }

    #[test]
    fn cost_to_pu_piecewise_scales_mw_only_and_trims() {
        // Model 1: MW breakpoints (even positions) ÷ base; cost ordinates (odd) raw.
        let cost = GenCost {
            model: 1,
            startup: 0.0,
            shutdown: 0.0,
            ncost: 4,
            coeffs: vec![
                0.0, 0.0, 100.0, 2500.0, 200.0, 5500.0, 250.0, 7250.0, 0.0, 0.0,
            ],
        };
        let out = cost_to_pu(&cost, 100.0);
        assert_eq!(out.len(), 8, "trimmed to 2·ncost, padding dropped");
        assert!(
            approx(out[0], 0.0)
                && approx(out[2], 1.0)
                && approx(out[4], 2.0)
                && approx(out[6], 2.5)
        );
        assert!(
            approx(out[1], 0.0)
                && approx(out[3], 2500.0)
                && approx(out[5], 5500.0)
                && approx(out[7], 7250.0)
        );
    }

    #[test]
    fn cost_rescale_round_trips() {
        // c2 p² + c1 p + c0 with base 100: per unit then back is the identity.
        let cost = GenCost {
            model: 2,
            startup: 0.0,
            shutdown: 0.0,
            ncost: 3,
            coeffs: vec![0.11, 5.0, 150.0],
        };
        let pu = cost_to_pu(&cost, 100.0);
        // p^2 coeff scales by 100^2, p^1 by 100, constant unchanged.
        assert!((pu[0] - 0.11 * 100.0 * 100.0).abs() < 1e-9);
        assert!((pu[1] - 5.0 * 100.0).abs() < 1e-9);
        assert!((pu[2] - 150.0).abs() < 1e-9);
        let back = cost_from_pu(&pu, 2, 100.0);
        for (a, b) in back.iter().zip(&cost.coeffs) {
            assert!((a - b).abs() < 1e-9);
        }
    }

    #[test]
    fn cost_rescale_passes_through_unknown_model() {
        // A model outside {1,2} has unknown coefficient semantics, so neither
        // direction may touch it; to_pu and from_pu must both be the identity,
        // or the round trip silently corrupts a curve we don't understand.
        let cost = GenCost {
            model: 0,
            startup: 0.0,
            shutdown: 0.0,
            ncost: 2,
            coeffs: vec![3.0, 7.0, 9.0],
        };
        let pu = cost_to_pu(&cost, 100.0);
        assert_eq!(pu, cost.coeffs, "to_pu must not scale an unknown model");
        let back = cost_from_pu(&pu, cost.model, 100.0);
        assert_eq!(back, cost.coeffs, "from_pu must not scale an unknown model");
    }

    #[test]
    fn cost_rescale_round_trips_piecewise() {
        // Model 1: cost_from_pu multiplies the MW breakpoints back by base and
        // leaves the cost ordinates, the exact inverse of cost_to_pu's even/odd
        // split. (cost_to_pu trims, cost_from_pu doesn't, so feed a trimmed row.)
        let cost = GenCost {
            model: 1,
            startup: 0.0,
            shutdown: 0.0,
            ncost: 4,
            coeffs: vec![0.0, 0.0, 100.0, 2500.0, 200.0, 5500.0, 250.0, 7250.0],
        };
        let pu = cost_to_pu(&cost, 100.0);
        let back = cost_from_pu(&pu, 1, 100.0);
        for (a, b) in back.iter().zip(&cost.coeffs) {
            assert!((a - b).abs() < 1e-9, "{a} != {b}");
        }
    }
}