powerio-prob 0.11.1

//! Solver-neutral LinDist3Flow OPF instance front end.
//!
//! This module owns formulation semantics that precede numerical coefficient
//! assembly: applicability, conductor-resolved radial topology, and the fixed
//! voltage phasor reference. Sparse matrices and affine/SOC preparation remain
//! in `powerio-matrix`.

use std::collections::{BTreeMap, BTreeSet};

use powerio_core::{Diagnostic, DiagnosticInfo, DiagnosticSeverity, Error};
use powerio_dist::{Configuration, DistLoadVoltageModel, MulticonductorNetwork};
use serde::{Deserialize, Serialize};

use super::McAcOpfInstance;
use crate::diagnostics::codes;
use crate::{MulticonductorOperatingPointQuantity, ObjectiveTerm};

/// Selection of the fixed phasors used to form LinDist3Flow coefficients.
#[derive(Clone, Copy, Debug, Default, PartialEq, Eq, Serialize, Deserialize)]
#[cfg_attr(feature = "schema", derive(schemars::JsonSchema))]
#[serde(rename_all = "snake_case")]
#[non_exhaustive]
pub enum LinDist3FlowReferencePolicy {
    /// Use the instance initial point when present, otherwise propagate source
    /// phasors through the no-load line topology.
    #[default]
    Auto,
    /// Require the instance to carry a complete voltage initial point.
    Explicit,
    /// Ignore any initial point and propagate source phasors without line drop.
    SourcePropagated,
}

/// Policy for network features outside the canonical affine model.
#[derive(Clone, Copy, Debug, Default, PartialEq, Eq, Serialize, Deserialize)]
#[cfg_attr(feature = "schema", derive(schemars::JsonSchema))]
#[serde(rename_all = "snake_case")]
#[non_exhaustive]
pub enum LinDist3FlowUnsupported {
    #[default]
    Reject,
    Lower,
    Approximate,
    Permissive,
}

/// Semantic choices used while creating a LinDist3Flow instance.
#[derive(Clone, Copy, Debug, Default, PartialEq, Eq, Serialize, Deserialize)]
#[cfg_attr(feature = "schema", derive(schemars::JsonSchema))]
#[non_exhaustive]
pub struct LinDist3FlowBuildOptions {
    pub reference_policy: LinDist3FlowReferencePolicy,
    pub unsupported: LinDist3FlowUnsupported,
    /// Require typed neutral-Kron provenance even when no explicit neutral
    /// remains in the network.
    pub require_neutral_provenance: bool,
}

impl LinDist3FlowBuildOptions {
    #[must_use]
    pub const fn with_reference_policy(mut self, policy: LinDist3FlowReferencePolicy) -> Self {
        self.reference_policy = policy;
        self
    }

    #[must_use]
    pub const fn with_unsupported(mut self, policy: LinDist3FlowUnsupported) -> Self {
        self.unsupported = policy;
        self
    }

    #[must_use]
    pub const fn with_required_neutral_provenance(mut self, required: bool) -> Self {
        self.require_neutral_provenance = required;
        self
    }
}

/// One retained voltage node, identified independently of dense row order.
#[derive(Clone, Debug, PartialEq, Eq, PartialOrd, Ord, Serialize, Deserialize)]
#[cfg_attr(feature = "schema", derive(schemars::JsonSchema))]
#[non_exhaustive]
pub struct LinDist3FlowNode {
    pub bus: String,
    pub terminal: String,
}

/// One line conductor oriented away from its conductor island's source.
#[derive(Clone, Debug, PartialEq, Eq, Serialize, Deserialize)]
#[cfg_attr(feature = "schema", derive(schemars::JsonSchema))]
#[non_exhaustive]
pub struct LinDist3FlowOrientedConductor {
    pub line: String,
    pub source_line_row: usize,
    pub conductor_position: usize,
    pub parent: LinDist3FlowNode,
    pub child: LinDist3FlowNode,
    /// Whether this direction is opposite the input line's from/to order.
    pub reversed: bool,
}

/// The source-rooted conductor forest used by the formulation.
#[derive(Clone, Debug, PartialEq, Eq, Serialize, Deserialize)]
#[cfg_attr(feature = "schema", derive(schemars::JsonSchema))]
#[non_exhaustive]
pub struct LinDist3FlowTopology {
    pub nodes: Vec<LinDist3FlowNode>,
    pub conductors: Vec<LinDist3FlowOrientedConductor>,
    pub roots: Vec<LinDist3FlowNode>,
    pub islands: Vec<Vec<LinDist3FlowNode>>,
}

/// Whether applicability checks permit instance construction.
#[derive(Clone, Copy, Debug, PartialEq, Eq, Serialize, Deserialize)]
#[cfg_attr(feature = "schema", derive(schemars::JsonSchema))]
#[serde(rename_all = "snake_case")]
#[non_exhaustive]
pub enum LinDist3FlowApplicabilityStatus {
    Applicable,
    Inapplicable,
}

/// Structured findings and topology evidence from the formulation front end.
#[derive(Clone, Debug, PartialEq, Serialize, Deserialize)]
#[non_exhaustive]
pub struct LinDist3FlowApplicability {
    pub status: LinDist3FlowApplicabilityStatus,
    pub diagnostics: Vec<Diagnostic>,
    pub roots: Vec<LinDist3FlowNode>,
    pub islands: Vec<Vec<LinDist3FlowNode>>,
    pub reference_provenance: Option<LinDist3FlowReferenceProvenance>,
    pub kron_reduced: bool,
    pub lowered: bool,
}

impl LinDist3FlowApplicability {
    #[must_use]
    pub const fn is_applicable(&self) -> bool {
        matches!(self.status, LinDist3FlowApplicabilityStatus::Applicable)
    }
}

/// Origin of the fixed voltage reference.
#[derive(Clone, Copy, Debug, PartialEq, Eq, Serialize, Deserialize)]
#[cfg_attr(feature = "schema", derive(schemars::JsonSchema))]
#[serde(rename_all = "snake_case")]
#[non_exhaustive]
pub enum LinDist3FlowReferenceProvenance {
    InitialPoint,
    SourcePropagated,
}

/// One reference voltage phasor in polar SI coordinates.
#[derive(Clone, Debug, PartialEq, Serialize, Deserialize)]
#[cfg_attr(feature = "schema", derive(schemars::JsonSchema))]
#[non_exhaustive]
pub struct LinDist3FlowReferenceVoltage {
    pub node: LinDist3FlowNode,
    /// Volts.
    pub magnitude: f64,
    /// Radians.
    pub angle: f64,
}

/// Complete immutable phasor reference in topology node order.
#[derive(Clone, Debug, PartialEq, Serialize, Deserialize)]
#[cfg_attr(feature = "schema", derive(schemars::JsonSchema))]
#[non_exhaustive]
pub struct LinDist3FlowReferenceState {
    pub provenance: LinDist3FlowReferenceProvenance,
    pub voltages: Vec<LinDist3FlowReferenceVoltage>,
}

impl LinDist3FlowReferenceState {
    #[must_use]
    pub fn voltage(&self, bus: &str, terminal: &str) -> Option<&LinDist3FlowReferenceVoltage> {
        self.voltages.iter().find(|voltage| {
            voltage.node.bus.eq_ignore_ascii_case(bus) && voltage.node.terminal == terminal
        })
    }
}

/// Matrix-free LinDist3Flow OPF instance.
#[derive(Clone, Debug)]
pub struct LinDist3FlowOpfInstance {
    base: McAcOpfInstance,
    topology: LinDist3FlowTopology,
    reference: LinDist3FlowReferenceState,
    applicability: LinDist3FlowApplicability,
    options: LinDist3FlowBuildOptions,
}

impl LinDist3FlowOpfInstance {
    /// Build a LinDist3Flow instance directly from a multiconductor network.
    ///
    /// # Errors
    /// As [`McAcOpfInstance::from_network`] and [`Self::from_mc_ac`].
    pub fn from_network(
        network: MulticonductorNetwork,
        options: LinDist3FlowBuildOptions,
    ) -> Result<Self, Error> {
        Self::from_mc_ac(McAcOpfInstance::from_network(network)?, options)
    }

    /// Compile formulation semantics from a multiconductor AC OPF instance.
    ///
    /// # Errors
    /// The applicability report contains an error, or the selected voltage
    /// reference is absent, incomplete, non-finite, or zero magnitude.
    pub fn from_mc_ac(
        base: McAcOpfInstance,
        options: LinDist3FlowBuildOptions,
    ) -> Result<Self, Error> {
        let (applicability, topology, reference) = assess(&base, options);
        if let Some(diagnostic) = applicability
            .diagnostics
            .iter()
            .find(|finding| finding.severity() == DiagnosticSeverity::Error)
        {
            return Err(error_from_diagnostic(diagnostic));
        }
        let topology = topology.ok_or_else(|| {
            Error::new(
                &codes::BUILD_LINDIST3FLOW_TOPOLOGY_INVALID,
                "the applicable assessment did not produce a conductor topology",
            )
        })?;
        let reference = reference.ok_or_else(|| {
            Error::new(
                &codes::BUILD_LINDIST3FLOW_REFERENCE_INVALID,
                "the applicable assessment did not produce a coefficient reference",
            )
        })?;
        Ok(Self {
            base,
            topology,
            reference,
            applicability,
            options,
        })
    }

    #[must_use]
    pub fn base_instance(&self) -> &McAcOpfInstance {
        &self.base
    }

    #[must_use]
    pub fn network(&self) -> &MulticonductorNetwork {
        self.base.network()
    }

    #[must_use]
    pub const fn topology(&self) -> &LinDist3FlowTopology {
        &self.topology
    }

    #[must_use]
    pub const fn reference(&self) -> &LinDist3FlowReferenceState {
        &self.reference
    }

    #[must_use]
    pub const fn applicability(&self) -> &LinDist3FlowApplicability {
        &self.applicability
    }

    #[must_use]
    pub const fn options(&self) -> LinDist3FlowBuildOptions {
        self.options
    }
}

fn finding(
    info: &'static DiagnosticInfo,
    message: impl Into<String>,
    target: Option<String>,
) -> Diagnostic {
    let mut diagnostic = Diagnostic::of(info, message);
    if let Some(target) = target {
        let _ = diagnostic.set_target(target);
    }
    diagnostic
}

fn error_from_diagnostic(diagnostic: &Diagnostic) -> Error {
    Error::new(
        diagnostic
            .registered_info()
            .expect("LinDist3Flow findings use registered codes"),
        diagnostic.message(),
    )
}

fn matrix_has_nonzero(matrix: &[Vec<f64>]) -> bool {
    matrix.iter().flatten().any(|value| *value != 0.0)
}

fn check_objective(instance: &McAcOpfInstance, diagnostics: &mut Vec<Diagnostic>) {
    let network = instance.network();
    let dispatch_cost = match instance.objective().terms() {
        [] => false,
        [ObjectiveTerm::ActivePowerDispatchCost] => true,
        terms => {
            diagnostics.push(finding(
                &codes::BUILD_LINDIST3FLOW_OBJECTIVE_UNSUPPORTED,
                format!(
                    "the strict LinDist3Flow slice accepts feasibility or exactly one active-power dispatch-cost term, but received {} term(s)",
                    terms.len()
                ),
                None,
            ));
            false
        }
    };
    if !dispatch_cost {
        return;
    }
    for (row, generator) in network.generators().iter().enumerate() {
        if generator.cost.is_none() {
            diagnostics.push(finding(
                &codes::BUILD_LINDIST3FLOW_COST_MISSING,
                format!(
                    "generator `{}` has no active-power dispatch cost; its coefficient is zero",
                    generator.name
                ),
                Some(format!("/generators/{row}/cost")),
            ));
        }
    }
    for (row, source) in network.sources().iter().enumerate() {
        if source.energy_cost_rate.is_none() {
            diagnostics.push(finding(
                &codes::BUILD_LINDIST3FLOW_COST_MISSING,
                format!(
                    "voltage source `{}` has no energy cost rate; its coefficient is zero",
                    source.name
                ),
                Some(format!("/sources/{row}/energy_cost_rate")),
            ));
        }
    }
}

fn supported_connection(configuration: Configuration, terminals: usize, channels: usize) -> bool {
    match configuration {
        Configuration::Wye => terminals == channels,
        Configuration::SinglePhase => terminals == channels || (terminals == 2 && channels == 1),
        Configuration::Delta => {
            (terminals == 2 && channels == 1) || (terminals == 3 && channels == 3)
        }
        _ => false,
    }
}

fn scalar_or_channels(values: &[f64], channels: usize) -> bool {
    matches!(values.len(), 1) || values.len() == channels
}

fn finite_scalar_or_channels(values: &[f64], channels: usize) -> bool {
    scalar_or_channels(values, channels) && values.iter().all(|value| value.is_finite())
}

fn positive_scalar_or_channels(values: &[f64], channels: usize) -> bool {
    scalar_or_channels(values, channels)
        && values.iter().all(|value| value.is_finite() && *value > 0.0)
}

fn valid_bounds(lower: Option<&[f64]>, upper: Option<&[f64]>, channels: usize) -> bool {
    match (lower, upper) {
        (None, None) => true,
        (Some(lower), Some(upper)) => {
            lower.len() == channels
                && upper.len() == channels
                && lower
                    .iter()
                    .zip(upper)
                    .all(|(lower, upper)| lower.is_finite() && upper.is_finite() && lower <= upper)
        }
        _ => false,
    }
}

fn square_finite(matrix: &[Vec<f64>], dimension: usize) -> bool {
    matrix.len() == dimension
        && matrix
            .iter()
            .all(|row| row.len() == dimension && row.iter().all(|value| value.is_finite()))
}

fn device_invalid(diagnostics: &mut Vec<Diagnostic>, message: impl Into<String>, target: String) {
    diagnostics.push(finding(
        &codes::BUILD_LINDIST3FLOW_DEVICE_INVALID,
        message,
        Some(target),
    ));
}

#[allow(clippy::too_many_lines)]
fn check_device_shapes(instance: &McAcOpfInstance, diagnostics: &mut Vec<Diagnostic>) {
    let network = instance.network();
    for (row, load) in network.loads().iter().enumerate() {
        let channels = load.p_nom.len();
        let mut valid = channels != 0
            && load.q_nom.len() == channels
            && load.p_nom.iter().all(|value| value.is_finite())
            && load.q_nom.iter().all(|value| value.is_finite())
            && supported_connection(load.configuration, load.terminal_map.len(), channels);
        valid &= match &load.voltage_model {
            DistLoadVoltageModel::ConstantPower { .. }
            | DistLoadVoltageModel::ConstantCurrent { .. }
            | DistLoadVoltageModel::Exponential { .. } => true,
            DistLoadVoltageModel::ConstantImpedance { v_nom } => {
                positive_scalar_or_channels(v_nom, channels)
            }
            DistLoadVoltageModel::Zip {
                v_nom,
                alpha_z,
                alpha_i,
                alpha_p,
                beta_z,
                beta_i,
                beta_p,
            } => {
                positive_scalar_or_channels(v_nom, channels)
                    && [alpha_z, alpha_i, alpha_p, beta_z, beta_i, beta_p]
                        .into_iter()
                        .all(|values| finite_scalar_or_channels(values, channels))
            }
            _ => false,
        };
        if !valid {
            device_invalid(
                diagnostics,
                format!(
                    "load `{}` has invalid channel, connection, nominal-power, or voltage-model dimensions",
                    load.name
                ),
                format!("/loads/{row}"),
            );
        }
    }

    for (row, generator) in network.generators().iter().enumerate() {
        let channels = generator.p_nom.len();
        let limits_valid = [generator.s_max.as_deref(), generator.i_max.as_deref()]
            .into_iter()
            .flatten()
            .all(|values| {
                values.len() == channels
                    && values.iter().all(|value| value.is_finite() && *value > 0.0)
            });
        let cost_valid = generator
            .cost
            .as_deref()
            .is_none_or(|values| finite_scalar_or_channels(values, channels));
        let valid = channels != 0
            && generator.q_nom.len() == channels
            && generator.p_nom.iter().all(|value| value.is_finite())
            && generator.q_nom.iter().all(|value| value.is_finite())
            && supported_connection(
                generator.configuration,
                generator.terminal_map.len(),
                channels,
            )
            && valid_bounds(
                generator.p_min.as_deref(),
                generator.p_max.as_deref(),
                channels,
            )
            && valid_bounds(
                generator.q_min.as_deref(),
                generator.q_max.as_deref(),
                channels,
            )
            && limits_valid
            && cost_valid;
        if !valid {
            device_invalid(
                diagnostics,
                format!(
                    "generator `{}` has invalid channel, connection, paired-bound, rating, or cost dimensions",
                    generator.name
                ),
                format!("/generators/{row}"),
            );
        }
    }

    for (row, shunt) in network.shunts().iter().enumerate() {
        let terminals = shunt.terminal_map.len();
        if terminals == 0
            || !square_finite(&shunt.g, terminals)
            || !square_finite(&shunt.b, terminals)
        {
            device_invalid(
                diagnostics,
                format!(
                    "shunt `{}` admittance matrices are not finite {terminals}x{terminals} arrays",
                    shunt.name
                ),
                format!("/shunts/{row}"),
            );
        }
    }

    for (row, source) in network.sources().iter().enumerate() {
        let channels = source.terminal_map.len();
        let valid = channels != 0
            && source.v_magnitude.len() == channels
            && source.v_angle.len() == channels
            && source
                .v_magnitude
                .iter()
                .all(|value| value.is_finite() && *value > 0.0)
            && source.v_angle.iter().all(|value| value.is_finite())
            && source
                .energy_cost_rate
                .as_deref()
                .is_none_or(|values| finite_scalar_or_channels(values, channels));
        if !valid {
            device_invalid(
                diagnostics,
                format!(
                    "voltage source `{}` has invalid terminal, phasor, or cost dimensions",
                    source.name
                ),
                format!("/sources/{row}"),
            );
        }
    }
}

fn check_supported_slice(
    instance: &McAcOpfInstance,
    options: LinDist3FlowBuildOptions,
    diagnostics: &mut Vec<Diagnostic>,
) {
    let network = instance.network();
    check_objective(instance, diagnostics);
    check_device_shapes(instance, diagnostics);
    if options.unsupported != LinDist3FlowUnsupported::Reject {
        diagnostics.push(finding(
            &codes::BUILD_LINDIST3FLOW_POLICY_UNAVAILABLE,
            format!(
                "the {:?} unsupported-data policy is reserved but not implemented yet",
                options.unsupported
            ),
            None,
        ));
    }

    let conventions = network.extras().get("bmopf_terminal_conventions");
    for (row, bus) in network.buses().iter().enumerate() {
        if bus.phase_indices(conventions).len() != bus.terminals.len() {
            diagnostics.push(finding(
                &codes::BUILD_LINDIST3FLOW_EXPLICIT_NEUTRAL,
                format!(
                    "bus `{}` still declares an explicit neutral conductor; apply neutral_kron_reduce first",
                    bus.id
                ),
                Some(format!("/buses/{row}/terminals")),
            ));
        }
    }
    if options.require_neutral_provenance && !network.extras().contains_key("powerio_neutral_kron")
    {
        diagnostics.push(finding(
            &codes::BUILD_LINDIST3FLOW_EXPLICIT_NEUTRAL,
            "the instance requires neutral-Kron provenance but the network carries none",
            None,
        ));
    }

    for (row, line) in network.lines().iter().enumerate() {
        let Some(code) = network.linecode(&line.linecode) else {
            continue;
        };
        if matrix_has_nonzero(&code.g_from)
            || matrix_has_nonzero(&code.b_from)
            || matrix_has_nonzero(&code.g_to)
            || matrix_has_nonzero(&code.b_to)
        {
            diagnostics.push(finding(
                &codes::BUILD_LINDIST3FLOW_UNSUPPORTED_COMPONENT,
                format!(
                    "line `{}` has pi shunt admittance; endpoint-shunt lowering is not implemented yet",
                    line.name
                ),
                Some(format!("/lines/{row}/linecode")),
            ));
        }
    }
    for (row, load) in network.loads().iter().enumerate() {
        let supported = match &load.voltage_model {
            DistLoadVoltageModel::ConstantPower { .. }
            | DistLoadVoltageModel::ConstantImpedance { .. } => true,
            DistLoadVoltageModel::Zip {
                alpha_i, beta_i, ..
            } => {
                alpha_i.iter().all(|value| value.abs() <= f64::EPSILON)
                    && beta_i.iter().all(|value| value.abs() <= f64::EPSILON)
            }
            DistLoadVoltageModel::ConstantCurrent { .. }
            | DistLoadVoltageModel::Exponential { .. }
            | _ => false,
        };
        if !supported {
            diagnostics.push(finding(
                &codes::BUILD_LINDIST3FLOW_UNSUPPORTED_COMPONENT,
                format!(
                    "load `{}` needs a current/exponential approximation outside the strict model",
                    load.name
                ),
                Some(format!("/loads/{row}/voltage_model")),
            ));
        }
    }

    for (family, count) in [
        ("switch", network.switches().len()),
        ("transformer", network.transformers().len()),
        ("capacitor", network.capacitors().len()),
        ("IBR", network.ibrs().len()),
        ("untyped object", network.untyped_objects().len()),
    ] {
        if count != 0 {
            diagnostics.push(finding(
                &codes::BUILD_LINDIST3FLOW_UNSUPPORTED_COMPONENT,
                format!(
                    "the initial strict LinDist3Flow slice does not yet compile {count} {family} record(s)"
                ),
                None,
            ));
        }
    }
}

/// Check the current strict LinDist3Flow boundary without constructing an
/// instance. Findings are stable diagnostics; callers can display all of them
/// before deciding whether to transform the network.
#[must_use]
pub fn check_lindist3flow_applicability(
    instance: &McAcOpfInstance,
    options: LinDist3FlowBuildOptions,
) -> LinDist3FlowApplicability {
    assess(instance, options).0
}

fn assess(
    instance: &McAcOpfInstance,
    options: LinDist3FlowBuildOptions,
) -> (
    LinDist3FlowApplicability,
    Option<LinDist3FlowTopology>,
    Option<LinDist3FlowReferenceState>,
) {
    let mut diagnostics = Vec::new();
    check_supported_slice(instance, options, &mut diagnostics);
    let topology = match build_topology(instance.network()) {
        Ok(topology) => Some(topology),
        Err(message) => {
            diagnostics.push(finding(
                &codes::BUILD_LINDIST3FLOW_TOPOLOGY_INVALID,
                message,
                None,
            ));
            None
        }
    };
    let reference = topology.as_ref().and_then(|topology| {
        match build_reference(instance, topology, options.reference_policy) {
            Ok(reference) => Some(reference),
            Err(error) => {
                diagnostics.extend(error.into_diagnostics());
                None
            }
        }
    });
    let status = if diagnostics
        .iter()
        .any(|finding| finding.severity() == DiagnosticSeverity::Error)
    {
        LinDist3FlowApplicabilityStatus::Inapplicable
    } else {
        LinDist3FlowApplicabilityStatus::Applicable
    };
    let roots = topology
        .as_ref()
        .map_or_else(Vec::new, |topology| topology.roots.clone());
    let islands = topology
        .as_ref()
        .map_or_else(Vec::new, |topology| topology.islands.clone());
    (
        LinDist3FlowApplicability {
            status,
            diagnostics,
            roots,
            islands,
            reference_provenance: reference.as_ref().map(|reference| reference.provenance),
            kron_reduced: instance
                .network()
                .extras()
                .contains_key("powerio_neutral_kron"),
            lowered: false,
        },
        topology,
        reference,
    )
}

#[derive(Clone, Copy)]
struct Edge {
    from: usize,
    to: usize,
    line: usize,
    conductor: usize,
}

struct UnionFind {
    parent: Vec<usize>,
}

impl UnionFind {
    fn new(n: usize) -> Self {
        Self {
            parent: (0..n).collect(),
        }
    }

    fn find(&mut self, node: usize) -> usize {
        let mut root = node;
        while self.parent[root] != root {
            root = self.parent[root];
        }
        let mut cursor = node;
        while self.parent[cursor] != root {
            let next = self.parent[cursor];
            self.parent[cursor] = root;
            cursor = next;
        }
        root
    }

    fn join(&mut self, left: usize, right: usize) -> bool {
        let left = self.find(left);
        let right = self.find(right);
        if left == right {
            return false;
        }
        self.parent[left.max(right)] = left.min(right);
        true
    }
}

fn node_key(bus: &str, terminal: &str) -> (String, String) {
    (bus.to_ascii_lowercase(), terminal.to_owned())
}

#[allow(clippy::too_many_lines)]
fn build_topology(network: &MulticonductorNetwork) -> Result<LinDist3FlowTopology, String> {
    let mut nodes = Vec::new();
    let mut positions = BTreeMap::new();
    let mut bus_ids = BTreeSet::new();
    let mut bus_positions = BTreeMap::new();
    for bus in network.buses() {
        if !bus_ids.insert(bus.id.to_ascii_lowercase()) {
            return Err(format!(
                "bus identity `{}` is duplicated case-insensitively",
                bus.id
            ));
        }
        bus_positions.insert(bus.id.to_ascii_lowercase(), bus_positions.len());
        let mut terminals = BTreeSet::new();
        for terminal in &bus.terminals {
            if !terminals.insert(terminal.clone()) {
                return Err(format!("bus `{}` repeats terminal `{terminal}`", bus.id));
            }
            let node = LinDist3FlowNode {
                bus: bus.id.clone(),
                terminal: terminal.clone(),
            };
            positions.insert(node_key(&node.bus, &node.terminal), nodes.len());
            nodes.push(node);
        }
    }
    if nodes.is_empty() {
        return Err("the network has no retained bus terminals".to_owned());
    }

    let mut forest = UnionFind::new(nodes.len());
    let mut bus_forest = UnionFind::new(bus_positions.len());
    let mut edges = Vec::new();
    for (line_row, line) in network.lines().iter().enumerate() {
        if line.terminal_map_from.len() != line.terminal_map_to.len()
            || line.terminal_map_from.is_empty()
        {
            return Err(format!(
                "line `{}` has empty or unequal terminal maps",
                line.name
            ));
        }
        let from_bus = *bus_positions
            .get(&line.bus_from.to_ascii_lowercase())
            .ok_or_else(|| {
                format!(
                    "line `{}` names unknown from bus `{}`",
                    line.name, line.bus_from
                )
            })?;
        let to_bus = *bus_positions
            .get(&line.bus_to.to_ascii_lowercase())
            .ok_or_else(|| {
                format!(
                    "line `{}` names unknown to bus `{}`",
                    line.name, line.bus_to
                )
            })?;
        bus_forest.join(from_bus, to_bus);
        for (conductor, (from_terminal, to_terminal)) in line
            .terminal_map_from
            .iter()
            .zip(&line.terminal_map_to)
            .enumerate()
        {
            let from = *positions
                .get(&node_key(&line.bus_from, from_terminal))
                .ok_or_else(|| {
                    format!(
                        "line `{}` names undeclared from terminal `{}/{from_terminal}`",
                        line.name, line.bus_from
                    )
                })?;
            let to = *positions
                .get(&node_key(&line.bus_to, to_terminal))
                .ok_or_else(|| {
                    format!(
                        "line `{}` names undeclared to terminal `{}/{to_terminal}`",
                        line.name, line.bus_to
                    )
                })?;
            if !forest.join(from, to) {
                return Err(format!(
                    "line `{}` conductor {} closes a cycle in the conductor-resolved graph",
                    line.name,
                    conductor + 1
                ));
            }
            edges.push(Edge {
                from,
                to,
                line: line_row,
                conductor,
            });
        }
    }

    let mut island_indices: BTreeMap<usize, Vec<usize>> = BTreeMap::new();
    for node in 0..nodes.len() {
        island_indices
            .entry(forest.find(node))
            .or_default()
            .push(node);
    }
    let mut source_nodes: BTreeMap<usize, Vec<usize>> = BTreeMap::new();
    let mut island_sources: BTreeMap<usize, Vec<&str>> = BTreeMap::new();
    for source in network.sources() {
        let bus = *bus_positions
            .get(&source.bus.to_ascii_lowercase())
            .ok_or_else(|| {
                format!(
                    "voltage source `{}` names unknown bus `{}`",
                    source.name, source.bus
                )
            })?;
        island_sources
            .entry(bus_forest.find(bus))
            .or_default()
            .push(&source.name);
        for terminal in &source.terminal_map {
            let node = *positions
                .get(&node_key(&source.bus, terminal))
                .ok_or_else(|| {
                    format!(
                        "voltage source `{}` names undeclared terminal `{}/{terminal}`",
                        source.name, source.bus
                    )
                })?;
            source_nodes
                .entry(forest.find(node))
                .or_default()
                .push(node);
        }
    }
    let mut physical_islands: BTreeMap<usize, Vec<&str>> = BTreeMap::new();
    for bus in network.buses() {
        let position = bus_positions[&bus.id.to_ascii_lowercase()];
        physical_islands
            .entry(bus_forest.find(position))
            .or_default()
            .push(&bus.id);
    }
    for (component, buses) in physical_islands {
        let sources: &[&str] = island_sources.get(&component).map_or(&[], Vec::as_slice);
        if sources.len() != 1 {
            return Err(format!(
                "physical island containing bus `{}` has {} voltage source records; exactly one is required",
                buses[0],
                sources.len()
            ));
        }
    }

    let mut ordered_islands = island_indices.into_values().collect::<Vec<_>>();
    ordered_islands.sort_by_key(|island| island[0]);
    let mut root_indices = Vec::with_capacity(ordered_islands.len());
    for island in &ordered_islands {
        let component = forest.find(island[0]);
        let roots: &[usize] = source_nodes.get(&component).map_or(&[], Vec::as_slice);
        if roots.len() != 1 {
            return Err(format!(
                "conductor island rooted at `{}/{}` contains {} fixed-voltage source terminals; exactly one is required",
                nodes[island[0]].bus,
                nodes[island[0]].terminal,
                roots.len()
            ));
        }
        root_indices.push(roots[0]);
    }

    let mut adjacency = vec![Vec::<(usize, usize)>::new(); nodes.len()];
    for (edge_index, edge) in edges.iter().enumerate() {
        adjacency[edge.from].push((edge_index, edge.to));
        adjacency[edge.to].push((edge_index, edge.from));
    }
    let mut directions = vec![None; edges.len()];
    for &root in &root_indices {
        let mut stack = vec![(root, usize::MAX)];
        while let Some((parent, incoming)) = stack.pop() {
            for &(edge_index, child) in &adjacency[parent] {
                if edge_index == incoming {
                    continue;
                }
                directions[edge_index] = Some((parent, child));
                stack.push((child, edge_index));
            }
        }
    }
    let mut line_directions = BTreeMap::new();
    for (edge, direction) in edges.iter().zip(&directions) {
        let (parent, _) = direction.expect("each source-rooted forest edge is visited");
        let reversed = parent != edge.from;
        if line_directions
            .insert(edge.line, reversed)
            .is_some_and(|previous| previous != reversed)
        {
            return Err(format!(
                "line `{}` is reached in conflicting directions across its coupled conductors",
                network.lines()[edge.line].name
            ));
        }
    }
    let conductors = edges
        .iter()
        .zip(directions)
        .map(|(edge, direction)| {
            let (parent, child) = direction.expect("each source-rooted forest edge is visited");
            LinDist3FlowOrientedConductor {
                line: network.lines()[edge.line].name.clone(),
                source_line_row: edge.line,
                conductor_position: edge.conductor,
                parent: nodes[parent].clone(),
                child: nodes[child].clone(),
                reversed: parent != edge.from,
            }
        })
        .collect();
    Ok(LinDist3FlowTopology {
        roots: root_indices
            .iter()
            .map(|&root| nodes[root].clone())
            .collect(),
        islands: ordered_islands
            .iter()
            .map(|island| island.iter().map(|&node| nodes[node].clone()).collect())
            .collect(),
        nodes,
        conductors,
    })
}

fn valid_phasor(magnitude: f64, angle: f64) -> bool {
    magnitude.is_finite() && magnitude > 0.0 && angle.is_finite()
}

fn invalid_reference(message: impl Into<String>) -> Error {
    Error::new(&codes::BUILD_LINDIST3FLOW_REFERENCE_INVALID, message)
}

fn explicit_reference(
    instance: &McAcOpfInstance,
    topology: &LinDist3FlowTopology,
) -> Result<LinDist3FlowReferenceState, Error> {
    let point = instance.initial_point().ok_or_else(|| {
        invalid_reference("the explicit reference policy requires a voltage initial point")
    })?;
    if point
        .values(MulticonductorOperatingPointQuantity::TerminalVoltageMagnitude)
        .is_none()
        || point
            .values(MulticonductorOperatingPointQuantity::TerminalVoltageAngle)
            .is_none()
    {
        return Err(invalid_reference(
            "the initial point must contain both terminal voltage magnitude and angle columns",
        ));
    }
    let mut voltages = Vec::with_capacity(topology.nodes.len());
    for node in &topology.nodes {
        let magnitude = point
            .terminal_voltage_magnitude(&node.bus, &node.terminal)
            .ok_or_else(|| {
                invalid_reference(format!(
                    "the initial point has no voltage magnitude for `{}/{}`",
                    node.bus, node.terminal
                ))
            })?;
        let angle = point
            .terminal_voltage_angle(&node.bus, &node.terminal)
            .ok_or_else(|| {
                invalid_reference(format!(
                    "the initial point has no voltage angle for `{}/{}`",
                    node.bus, node.terminal
                ))
            })?;
        if !valid_phasor(magnitude, angle) {
            return Err(invalid_reference(format!(
                "the initial point voltage for `{}/{}` is zero or non-finite",
                node.bus, node.terminal
            )));
        }
        voltages.push(LinDist3FlowReferenceVoltage {
            node: node.clone(),
            magnitude,
            angle,
        });
    }
    Ok(LinDist3FlowReferenceState {
        provenance: LinDist3FlowReferenceProvenance::InitialPoint,
        voltages,
    })
}

#[allow(clippy::too_many_lines)]
fn propagated_reference(
    instance: &McAcOpfInstance,
    topology: &LinDist3FlowTopology,
) -> Result<LinDist3FlowReferenceState, Error> {
    let mut positions = BTreeMap::new();
    for (position, node) in topology.nodes.iter().enumerate() {
        positions.insert(node_key(&node.bus, &node.terminal), position);
    }
    let mut values = vec![None; topology.nodes.len()];
    for source in instance.network().sources() {
        if source.v_magnitude.len() != source.terminal_map.len()
            || source.v_angle.len() != source.terminal_map.len()
        {
            return Err(invalid_reference(format!(
                "voltage source `{}` values do not align with its terminal map",
                source.name
            )));
        }
        for (terminal_position, terminal) in source.terminal_map.iter().enumerate() {
            let node = *positions
                .get(&node_key(&source.bus, terminal))
                .ok_or_else(|| {
                    invalid_reference(format!(
                        "voltage source `{}` names unknown terminal `{}/{terminal}`",
                        source.name, source.bus
                    ))
                })?;
            let magnitude = source.v_magnitude[terminal_position];
            let angle = source.v_angle[terminal_position];
            if !valid_phasor(magnitude, angle) {
                return Err(invalid_reference(format!(
                    "voltage source `{}` terminal `{terminal}` is zero or non-finite",
                    source.name
                )));
            }
            values[node] = Some((magnitude, angle));
        }
    }

    let mut children = vec![Vec::new(); topology.nodes.len()];
    for conductor in &topology.conductors {
        let parent = positions[&node_key(&conductor.parent.bus, &conductor.parent.terminal)];
        let child = positions[&node_key(&conductor.child.bus, &conductor.child.terminal)];
        children[parent].push(child);
    }
    for root in &topology.roots {
        let root = positions[&node_key(&root.bus, &root.terminal)];
        let mut stack = vec![root];
        while let Some(parent) = stack.pop() {
            let value = values[parent].ok_or_else(|| {
                invalid_reference(format!(
                    "source root `{}/{}` has no phasor",
                    topology.nodes[parent].bus, topology.nodes[parent].terminal
                ))
            })?;
            for &child in &children[parent] {
                values[child] = Some(value);
                stack.push(child);
            }
        }
    }
    let voltages = topology
        .nodes
        .iter()
        .cloned()
        .zip(values)
        .map(|(node, value)| {
            let (magnitude, angle) = value.ok_or_else(|| {
                invalid_reference(format!(
                    "node `{}/{}` was not reached from a source",
                    node.bus, node.terminal
                ))
            })?;
            Ok(LinDist3FlowReferenceVoltage {
                node,
                magnitude,
                angle,
            })
        })
        .collect::<Result<Vec<_>, Error>>()?;
    Ok(LinDist3FlowReferenceState {
        provenance: LinDist3FlowReferenceProvenance::SourcePropagated,
        voltages,
    })
}

fn build_reference(
    instance: &McAcOpfInstance,
    topology: &LinDist3FlowTopology,
    policy: LinDist3FlowReferencePolicy,
) -> Result<LinDist3FlowReferenceState, Error> {
    match policy {
        LinDist3FlowReferencePolicy::Auto if instance.initial_point().is_some() => {
            explicit_reference(instance, topology)
        }
        LinDist3FlowReferencePolicy::Explicit => explicit_reference(instance, topology),
        LinDist3FlowReferencePolicy::Auto | LinDist3FlowReferencePolicy::SourcePropagated => {
            propagated_reference(instance, topology)
        }
    }
}