//! System graph structure for thermodynamic simulation.
//!
//! This module provides the core graph representation of a thermodynamic system,
//! where nodes are components and edges represent flow connections. Edges index
//! into the solver's state vector (ṁ, P and h per edge — see [`EdgeKind`]).
//!
//! Multi-circuit support (Story 3.3): A machine can have up to 5 independent
//! circuits (valid circuit IDs: 0, 1, 2, 3, 4). Each node belongs to exactly one
//! circuit. Flow edges connect only nodes within the same circuit.

use entropyk_components::{
    validate_port_continuity, Component, ComponentError, ConnectionError, JacobianBuilder,
    ResidualVector, StateSlice,
};
use petgraph::algo;
use petgraph::graph::{EdgeIndex, Graph, NodeIndex};
use petgraph::visit::EdgeRef;
use petgraph::Directed;
use std::collections::HashMap;

use crate::coupling::{has_circular_dependencies, ThermalCoupling};
use crate::dof::{
    align_roles, ComponentEquationBlock, DofReport, EquationRole, SystemDofBalance, SystemDofError,
    UnknownKind,
};
use crate::error::{AddEdgeError, TopologyError};
use crate::inverse::{
    BoundedVariable, BoundedVariableError, BoundedVariableId, Constraint, ConstraintError,
    ConstraintId, DoFError, InverseControlConfig, SaturatedController,
};
use entropyk_core::{CircuitId, Temperature};

/// Maximum circuit ID (inclusive). Machine supports up to 5 circuits.
pub const MAX_CIRCUIT_ID: u16 = 4;

/// Default initial mass-flow seed [kg/s] used by the `SmartInitializer` to
/// populate the ṁ slot of each edge in the solver state vector. After CM1.3
/// ṁ is a genuine Newton unknown driven by component mass-flow residuals; this
/// constant is the initial guess, not a placeholder closure.
pub const DEFAULT_MASS_FLOW_SEED_KG_S: f64 = 0.05;
/// Minimum pressure allowed for Newton pressure unknowns [Pa].
///
/// This is a solver-domain bound, not a component fallback: trial states below
/// this pressure are outside the thermodynamic domain used by fluid backends.
pub const MIN_SOLVER_PRESSURE_PA: f64 = 10_000.0;

/// Backward-compatibility alias — tests and tools may still reference the old
/// TEMP name until they migrate to `DEFAULT_MASS_FLOW_SEED_KG_S` (T6).
#[deprecated(since = "0.0.0", note = "Use DEFAULT_MASS_FLOW_SEED_KG_S instead")]
#[allow(dead_code)]
pub const TEMP_MASS_FLOW_SEED_KG_S: f64 = DEFAULT_MASS_FLOW_SEED_KG_S;

/// Aggregate thermodynamic performance of a solved refrigeration/heat-pump cycle.
///
/// Built by [`System::cycle_performance`] from each component's
/// [`entropyk_components::Component::energy_transfers`] evaluated at the solved
/// state. Duties are classified purely by the physical sign of the reported heat
/// and work — no component-type strings — so it works for any topology whose
/// components implement `energy_transfers`.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct CyclePerformance {
    /// Total cooling duty [W]: heat absorbed by the refrigerant (evaporators).
    pub q_cooling_w: f64,
    /// Total heating duty [W]: heat rejected by the refrigerant (condensers/gas coolers).
    pub q_heating_w: f64,
    /// Total shaft/electrical work input [W] (compressors, pumps, fans).
    pub work_input_w: f64,
    /// Number of components that contributed an `energy_transfers` reading.
    pub components_counted: usize,
}

impl CyclePerformance {
    /// Cooling coefficient of performance `EER = Q_cooling / W_input`.
    ///
    /// Returns `None` when no work is drawn (undefined ratio).
    pub fn cop_cooling(&self) -> Option<f64> {
        if self.work_input_w > 0.0 {
            Some(self.q_cooling_w / self.work_input_w)
        } else {
            None
        }
    }

    /// Heating coefficient of performance `COP = Q_heating / W_input`.
    ///
    /// Returns `None` when no work is drawn (undefined ratio).
    pub fn cop_heating(&self) -> Option<f64> {
        if self.work_input_w > 0.0 {
            Some(self.q_heating_w / self.work_input_w)
        } else {
            None
        }
    }
}

/// Classifies a flow edge by the physical circuit it belongs to.
///
/// The kind determines how many state-vector unknowns the edge carries (see
/// [`EdgeKind::n_unknowns`]). Today every edge defaults to
/// [`EdgeKind::Refrigerant`], preserving the historical single-circuit
/// refrigerant behaviour. Hydraulic (pumped liquid) edges share the same
/// `(ṁ, P, h)` layout; moist-air edges will gain a fourth unknown (humidity
/// ratio `W`) in a later story.
#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
pub enum EdgeKind {
    /// Refrigerant flow `(ṁ, P, h)` — the default for backward compatibility.
    #[default]
    Refrigerant,
    /// Pumped hydraulic (single-phase liquid) flow `(ṁ, P, h)`.
    Hydraulic,
    /// Moist-air flow. Today `(ṁ, P, h)`; a humidity-ratio unknown `W` is
    /// added in Story CM6.1, at which point `n_unknowns()` returns 4.
    MoistAir,
}

impl EdgeKind {
    /// Number of state-vector unknowns this edge contributes.
    ///
    /// Currently 3 (`ṁ, P, h`) for every kind. MoistAir will return 4 once the
    /// humidity-ratio unknown `W` is introduced (Story CM6.1).
    pub fn n_unknowns(&self) -> usize {
        match self {
            EdgeKind::Refrigerant | EdgeKind::Hydraulic | EdgeKind::MoistAir => 3,
        }
    }
}

/// Weight for flow edges in the system graph.
///
/// Each edge represents a flow connection between two component ports and stores
/// the state vector indices for mass flow (ṁ), pressure (P) and enthalpy (h) at
/// that connection, plus the [`EdgeKind`] classifying the edge.
///
/// After [`System::finalize`] runs the topology presolve (CM1.4), all edges in
/// the same series branch share one `state_index_m` (the branch-level ṁ unknown).
/// P and h indices remain per-edge.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub struct FlowEdge {
    /// Classification of the edge (refrigerant / hydraulic / moist-air).
    pub kind: EdgeKind,
    /// State vector index for mass flow rate (kg/s).
    ///
    /// After CM1.4 topology presolve all edges in the same series branch share
    /// this index — the branch-level ṁ unknown.
    pub state_index_m: usize,
    /// State vector index for pressure (Pa)
    pub state_index_p: usize,
    /// State vector index for enthalpy (J/kg)
    pub state_index_h: usize,
    /// Local port index on the **source** component this edge leaves from
    /// (recorded by [`System::add_edge_with_ports`]; defaults to 1 = outlet).
    pub source_port: usize,
    /// Local port index on the **target** component this edge enters at
    /// (recorded by [`System::add_edge_with_ports`]; defaults to 0 = inlet).
    pub target_port: usize,
    /// Branch identifier assigned by [`topology::presolve_mass_flow_topology`].
    ///
    /// All edges in the same series branch (no junction in between) receive the
    /// same `mass_flow_branch_id` and share one `state_index_m`. Sentinel value
    /// `usize::MAX` means "not yet assigned" (before `finalize()` is called).
    pub mass_flow_branch_id: usize,
}

impl FlowEdge {
    /// Creates a refrigerant flow edge with unassigned state indices.
    ///
    /// State indices and branch id are assigned during [`System::finalize`].
    pub(crate) fn new_unassigned() -> Self {
        FlowEdge {
            kind: EdgeKind::Refrigerant,
            state_index_m: 0,
            state_index_p: 0,
            state_index_h: 0,
            mass_flow_branch_id: usize::MAX,
            source_port: 1,
            target_port: 0,
        }
    }

    /// Creates an unassigned edge carrying explicit local port indices.
    pub(crate) fn new_unassigned_with_ports(source_port: usize, target_port: usize) -> Self {
        FlowEdge {
            source_port,
            target_port,
            ..Self::new_unassigned()
        }
    }
}

/// System graph structure.
///
/// Nodes are components (`Box<dyn Component>`), edges are flow connections with
/// state indices. The state vector layout is:
///
/// ```text
/// [m_edge0, P_edge0, h_edge0, m_edge1, P_edge1, h_edge1, ...]
/// ```
///
/// Edge order follows the graph's internal edge iteration order (stable after
/// `finalize()` is called).
pub struct System {
    graph: Graph<Box<dyn Component>, FlowEdge, Directed>,
    /// Maps EdgeIndex to (state_index_m, state_index_p, state_index_h) - built in finalize()
    edge_to_state: HashMap<EdgeIndex, (usize, usize, usize)>,
    /// Maps NodeIndex to CircuitId. Nodes without entry default to circuit 0.
    node_to_circuit: HashMap<NodeIndex, CircuitId>,
    /// Thermal couplings between circuits (heat transfer without fluid mixing).
    thermal_couplings: Vec<ThermalCoupling>,
    /// Constraints for inverse control (output - target = 0)
    constraints: HashMap<ConstraintId, Constraint>,
    /// Bounded control variables for inverse control (with box constraints)
    bounded_variables: HashMap<BoundedVariableId, BoundedVariable>,
    /// Inverse control configuration (constraint → control variable mappings)
    inverse_control: InverseControlConfig,
    /// Saturated PI loops solved as additive `(u, x)` unknown/residual pairs.
    saturated_controllers: Vec<SaturatedController>,
    /// Physical free-actuator unknowns (arch-6). Each entry is a registered
    /// bounded variable id whose value is an *extra* solver unknown appended at
    /// the end of the state vector. Unlike saturated controllers, a free actuator
    /// contributes **no** system residual — the owning component supplies the
    /// closing equation (e.g. an EXV orifice-flow residual) inside its own
    /// `n_equations()` block, keeping the DoF balanced (+1 unknown / +1 eq).
    free_actuators: Vec<BoundedVariableId>,
    /// Registry of component names for constraint validation.
    /// Maps human-readable names (e.g., "evaporator") to NodeIndex.
    component_names: HashMap<String, NodeIndex>,
    /// Per-component calibration indices snapshot, captured at `finalize()` time.
    /// Lets external readers (e.g. result extraction) learn which Z-factors
    /// (z_ua, z_dp, z_flow, …) were promoted to free unknowns and where they
    /// live in the state vector, so the solved values can be surfaced to the user.
    calib_indices_by_name: HashMap<String, entropyk_core::CalibIndices>,
    finalized: bool,
    total_state_len: usize,
    /// When `true` (default), `finalize` rejects **over-constrained** systems.
    /// Under-constrained systems always warn; use [`validate_system_dof`] for a
    /// hard square check in production paths.
    enforce_dof_gate: bool,
}

impl System {
    /// Creates a new empty system graph.
    pub fn new() -> Self {
        Self {
            graph: Graph::new(),
            edge_to_state: HashMap::new(),
            node_to_circuit: HashMap::new(),
            thermal_couplings: Vec::new(),
            constraints: HashMap::new(),
            bounded_variables: HashMap::new(),
            inverse_control: InverseControlConfig::new(),
            saturated_controllers: Vec::new(),
            free_actuators: Vec::new(),
            component_names: HashMap::new(),
            calib_indices_by_name: HashMap::new(),
            finalized: false,
            total_state_len: 0,
            enforce_dof_gate: true,
        }
    }

    /// Enables or disables the over-constrained DoF hard gate in [`finalize`].
    ///
    /// Keep this **enabled** for real-machine models. Disable only for topology
    /// unit tests or intentionally non-physical mock assemblies (e.g. macro
    /// residual plumbing tests). Production solvers should call
    /// [`validate_system_dof`] after finalize for a full square-system check.
    pub fn set_enforce_dof_gate(&mut self, enforce: bool) {
        self.enforce_dof_gate = enforce;
    }

    /// Adds a component as a node in the default circuit (circuit 0) and returns its node index.
    ///
    /// For multi-circuit machines, use [`add_component_to_circuit`](Self::add_component_to_circuit).
    pub fn add_component(&mut self, component: Box<dyn Component>) -> NodeIndex {
        self.add_component_to_circuit(component, CircuitId::ZERO)
            .unwrap()
    }

    /// Adds a component as a node in the specified circuit and returns its node index.
    ///
    /// # Errors
    ///
    /// Returns `TopologyError::TooManyCircuits` if `circuit_id.0 > 4`.
    pub fn add_component_to_circuit(
        &mut self,
        component: Box<dyn Component>,
        circuit_id: CircuitId,
    ) -> Result<NodeIndex, TopologyError> {
        if circuit_id.0 > MAX_CIRCUIT_ID {
            return Err(TopologyError::TooManyCircuits {
                requested: circuit_id.0,
            });
        }
        self.finalized = false;
        let node_idx = self.graph.add_node(component);
        self.node_to_circuit.insert(node_idx, circuit_id);
        Ok(node_idx)
    }

    /// Returns the circuit ID for a node, or circuit 0 if not found (backward compat).
    pub fn node_circuit(&self, node: NodeIndex) -> CircuitId {
        self.node_to_circuit
            .get(&node)
            .copied()
            .unwrap_or(CircuitId::ZERO)
    }

    /// Returns the circuit ID for an edge based on its source node.
    ///
    /// # Panics
    ///
    /// Panics if the edge index is invalid.
    pub fn edge_circuit(&self, edge: EdgeIndex) -> CircuitId {
        let (src, _tgt) = self.graph.edge_endpoints(edge).expect("invalid edge index");
        self.node_circuit(src)
    }

    /// Adds a flow edge from `source` to `target` without port validation.
    ///
    /// **No port compatibility validation is performed.** Use
    /// [`add_edge_with_ports`](Self::add_edge_with_ports) for components with ports to validate
    /// fluid, pressure, and enthalpy continuity. This method is intended for components without
    /// ports (e.g. mock components in tests).
    ///
    /// Flow edges connect only nodes within the same circuit. Cross-circuit connections
    /// are rejected (thermal coupling is Story 3.4).
    ///
    /// State indices are assigned when `finalize()` is called.
    ///
    /// # Errors
    ///
    /// Returns `TopologyError::CrossCircuitConnection` if source and target are in different circuits.
    pub fn add_edge(
        &mut self,
        source: NodeIndex,
        target: NodeIndex,
    ) -> Result<EdgeIndex, TopologyError> {
        let src_circuit = self.node_circuit(source);
        let tgt_circuit = self.node_circuit(target);
        if src_circuit != tgt_circuit {
            tracing::warn!(
                "Cross-circuit edge rejected: source circuit {}, target circuit {}",
                src_circuit.0,
                tgt_circuit.0
            );
            return Err(TopologyError::CrossCircuitConnection {
                source_circuit: src_circuit.0,
                target_circuit: tgt_circuit.0,
            });
        }

        // Safety check: Warn if connecting components with ports using the non-validating method
        if let (Some(src), Some(tgt)) = (
            self.graph.node_weight(source),
            self.graph.node_weight(target),
        ) {
            if !src.get_ports().is_empty() || !tgt.get_ports().is_empty() {
                tracing::warn!(
                    "add_edge called on components with ports (src: {:?}, tgt: {:?}). \
                     This bypasses port validation. Use add_edge_with_ports instead.",
                    source,
                    target
                );
            }
        }

        self.finalized = false;
        Ok(self
            .graph
            .add_edge(source, target, FlowEdge::new_unassigned()))
    }

    /// Adds a flow edge from `source` outlet port to `target` inlet port with validation.
    ///
    /// Validates circuit membership (same circuit), then fluid compatibility, pressure and
    /// enthalpy continuity using port.rs tolerances. For 2-port components: `source_port_idx=1`
    /// (outlet), `target_port_idx=0` (inlet).
    ///
    /// # Errors
    ///
    /// Returns `AddEdgeError::Topology` if source and target are in different circuits.
    /// Returns `AddEdgeError::Connection` if ports are incompatible (fluid, pressure, or enthalpy mismatch).
    pub fn add_edge_with_ports(
        &mut self,
        source: NodeIndex,
        source_port_idx: usize,
        target: NodeIndex,
        target_port_idx: usize,
    ) -> Result<EdgeIndex, AddEdgeError> {
        // Circuit validation first
        let src_circuit = self.node_circuit(source);
        let tgt_circuit = self.node_circuit(target);
        if src_circuit != tgt_circuit {
            tracing::warn!(
                "Cross-circuit edge rejected: source circuit {}, target circuit {}",
                src_circuit.0,
                tgt_circuit.0
            );
            return Err(TopologyError::CrossCircuitConnection {
                source_circuit: src_circuit.0,
                target_circuit: tgt_circuit.0,
            }
            .into());
        }

        let source_comp = self
            .graph
            .node_weight(source)
            .ok_or_else(|| ConnectionError::InvalidNodeIndex(source.index()))?;
        let target_comp = self
            .graph
            .node_weight(target)
            .ok_or_else(|| ConnectionError::InvalidNodeIndex(target.index()))?;

        let source_ports = source_comp.get_ports();
        let target_ports = target_comp.get_ports();

        if source_ports.is_empty() && target_ports.is_empty() {
            // No ports: add edge without validation (backward compat), but keep
            // the requested local port indices so multi-port components can
            // resolve their per-port edge wiring in finalize().
            self.finalized = false;
            return Ok(self.graph.add_edge(
                source,
                target,
                FlowEdge::new_unassigned_with_ports(source_port_idx, target_port_idx),
            ));
        }

        if source_port_idx >= source_ports.len() {
            return Err(ConnectionError::InvalidPortIndex {
                index: source_port_idx,
                port_count: source_ports.len(),
                max_index: source_ports.len().saturating_sub(1),
            }
            .into());
        }
        if target_port_idx >= target_ports.len() {
            return Err(ConnectionError::InvalidPortIndex {
                index: target_port_idx,
                port_count: target_ports.len(),
                max_index: target_ports.len().saturating_sub(1),
            }
            .into());
        }

        let outlet = &source_ports[source_port_idx];
        let inlet = &target_ports[target_port_idx];
        if let Err(ref e) = validate_port_continuity(outlet, inlet) {
            tracing::warn!("Port validation failed: {}", e);
            return Err(e.clone().into());
        }

        self.finalized = false;
        Ok(self.graph.add_edge(
            source,
            target,
            FlowEdge::new_unassigned_with_ports(source_port_idx, target_port_idx),
        ))
    }

    /// Finalizes the graph: runs topology presolve, builds edge→state index
    /// mapping, and validates topology.
    ///
    /// # State vector layout (CM1.4)
    ///
    /// After the mass-flow topology presolve (`|B|` branches detected):
    ///
    /// ```text
    /// [ṁ_B0, ṁ_B1, …, ṁ_{B-1},   ← one slot per branch
    ///  P_0, h_0, P_1, h_1, …]     ← two slots per edge (P and h)
    /// ```
    ///
    /// Total edge-section length = `|B| + 2 × |E|`.
    ///
    /// All edges in the same series branch share the same `state_index_m`.
    /// `state_index_p` and `state_index_h` remain unique per edge.
    ///
    /// # Errors
    ///
    /// Returns `TopologyError` if:
    /// - Any node is isolated (no edges)
    /// - The graph is empty (no components)
    pub fn finalize(&mut self) -> Result<(), TopologyError> {
        self.validate_topology()?;

        if !self.thermal_couplings.is_empty() && has_circular_dependencies(self.thermal_couplings())
        {
            tracing::warn!("Circular thermal coupling detected, simultaneous solving required");
        }

        // CM1.4: run topology presolve to assign mass_flow_branch_id to each edge.
        let n_branches = crate::topology::presolve_mass_flow_topology(&mut self.graph);
        tracing::debug!(
            n_branches,
            edge_count = self.graph.edge_count(),
            "Mass-flow topology presolve complete"
        );

        // ── Pass 1: allocate one ṁ state slot per branch ─────────────────
        // Branch ids are contiguous 0..n_branches, so we can use a Vec.
        let mut branch_m_idx: Vec<Option<usize>> = vec![None; n_branches];
        let mut idx = 0usize; // running state-vector position

        // Allocate ṁ slots for branches in branch_id order.
        // We need branch_id → m_idx. Walk edges to collect branch ids and
        // allocate one slot per branch in the order first encountered.
        for edge_idx in self.graph.edge_indices() {
            let branch_id = self
                .graph
                .edge_weight(edge_idx)
                .map(|w| w.mass_flow_branch_id)
                .unwrap_or(0);
            if branch_id < n_branches && branch_m_idx[branch_id].is_none() {
                branch_m_idx[branch_id] = Some(idx);
                idx += 1;
            }
        }

        // ── Pass 2: assign per-edge (m_idx, p_idx, h_idx) ────────────────
        self.edge_to_state.clear();
        for edge_idx in self.graph.edge_indices() {
            let (branch_id, kind) = self
                .graph
                .edge_weight(edge_idx)
                .map(|w| (w.mass_flow_branch_id, w.kind))
                .unwrap_or((0, EdgeKind::default()));

            // Shared ṁ index from the branch's pre-allocated slot.
            let m_idx = branch_m_idx.get(branch_id).copied().flatten().unwrap_or(0);

            // Per-edge P and h indices (always unique).
            let p_idx = idx;
            let h_idx = idx + 1;
            idx += kind.n_unknowns() - 1; // subtract 1 because ṁ is shared

            self.edge_to_state.insert(edge_idx, (m_idx, p_idx, h_idx));
            if let Some(weight) = self.graph.edge_weight_mut(edge_idx) {
                weight.state_index_m = m_idx;
                weight.state_index_p = p_idx;
                weight.state_index_h = h_idx;
            }
        }
        self.finalized = true;

        // Notify each component about its position in the global state vector.
        // Collect context first (to avoid borrow conflicts), then apply.
        //
        // State offset: for the parent System the state layout is a flat array of
        // sum(EdgeKind::n_unknowns()) entries (3 per edge today) for the
        // *parent's own* edges.  An embedded MacroComponent has its internal
        // state appended after the parent edges, addressed via its own
        // global_state_offset.  We start with the total per-edge unknown count
        // as the base offset, and accumulate `internal_state_len` for each
        // component.
        let mut current_offset = idx;

        // Gather (node_idx, offset, incident_edge_indices) before mutating nodes.
        // Each entry carries a (m_idx, p_idx, h_idx) triple per incident edge so
        // components can wire their mass-flow residuals to the edge ṁ unknowns (CM1.3).
        #[allow(clippy::type_complexity)]
        let mut node_context: Vec<(
            petgraph::graph::NodeIndex,
            usize,
            Vec<(usize, usize, usize)>,
            Vec<Option<(usize, usize, usize)>>,
        )> = Vec::new();
        for node_idx in self.graph.node_indices() {
            let component = self.graph.node_weight(node_idx).unwrap();
            let mut incident: Vec<(usize, usize, usize)> = Vec::new();
            // Deterministic per-port wiring: port_edges[i] = edge state indices
            // of the edge attached at local port i (incoming edges use the
            // edge's target_port, outgoing edges its source_port).
            let mut port_edges: Vec<Option<(usize, usize, usize)>> = Vec::new();
            let mut record_port = |slot: usize, triple: (usize, usize, usize)| {
                if port_edges.len() <= slot {
                    port_edges.resize(slot + 1, None);
                }
                port_edges[slot] = Some(triple);
            };
            for edge_ref in self
                .graph
                .edges_directed(node_idx, petgraph::Direction::Incoming)
            {
                if let Some(&(m, p, h)) = self.edge_to_state.get(&edge_ref.id()) {
                    incident.push((m, p, h));
                    record_port(edge_ref.weight().target_port, (m, p, h));
                }
            }
            for edge_ref in self
                .graph
                .edges_directed(node_idx, petgraph::Direction::Outgoing)
            {
                if let Some(&(m, p, h)) = self.edge_to_state.get(&edge_ref.id()) {
                    incident.push((m, p, h));
                    record_port(edge_ref.weight().source_port, (m, p, h));
                }
            }
            node_context.push((node_idx, current_offset, incident, port_edges));

            // Advance the global offset by this component's internal state length
            current_offset += component.internal_state_len();
        }

        self.total_state_len = current_offset;

        // Notify components about their calibration control variables (Story 5.5)
        let mut comp_calib_indices: HashMap<String, entropyk_core::CalibIndices> = HashMap::new();
        for (index, id) in self.inverse_control.linked_controls().enumerate() {
            if let Some(bounded_var) = self.bounded_variables.get(id) {
                if let Some(comp_id) = bounded_var.component_id() {
                    let indices = comp_calib_indices.entry(comp_id.to_string()).or_default();
                    let state_idx = self.total_state_len + index;

                    let id_str = id.as_str();
                    if let Some(factor) = entropyk_core::id_ends_with_calib_suffix(id_str) {
                        match factor {
                            entropyk_core::Z_FLOW => indices.z_flow = Some(state_idx),
                            entropyk_core::Z_FLOW_ECO => indices.z_flow_eco = Some(state_idx),
                            entropyk_core::Z_DP => indices.z_dp = Some(state_idx),
                            entropyk_core::Z_UA => indices.z_ua = Some(state_idx),
                            entropyk_core::Z_POWER => indices.z_power = Some(state_idx),
                            entropyk_core::Z_ETAV => indices.z_etav = Some(state_idx),
                            _ => {}
                        }
                    }
                }
            }
        }
        let saturated_base = self.total_state_len
            + self.inverse_control.mapping_count()
            + self.coupling_residual_count();
        for (index, ctrl) in self.saturated_controllers.iter().enumerate() {
            if let Some(bounded_var) = self.bounded_variables.get(ctrl.actuator()) {
                if let Some(comp_id) = bounded_var.component_id() {
                    let indices = comp_calib_indices.entry(comp_id.to_string()).or_default();
                    let state_idx = saturated_base + 2 * index;

                    let id_str = ctrl.actuator().as_str();
                    if let Some(factor) = entropyk_core::id_ends_with_calib_suffix(id_str) {
                        match factor {
                            entropyk_core::Z_FLOW => indices.z_flow = Some(state_idx),
                            entropyk_core::Z_FLOW_ECO => indices.z_flow_eco = Some(state_idx),
                            entropyk_core::Z_DP => indices.z_dp = Some(state_idx),
                            entropyk_core::Z_UA => indices.z_ua = Some(state_idx),
                            entropyk_core::Z_POWER => indices.z_power = Some(state_idx),
                            entropyk_core::Z_ETAV => indices.z_etav = Some(state_idx),
                            _ => {}
                        }
                    } else if id_str.ends_with("injection") || id_str.ends_with("actuator") {
                        // Generic physical actuator (arch-6), e.g. liquid injection
                        // ratio φ_inj: the controls[] loop drives this slot directly.
                        indices.actuator = Some(state_idx);
                    }
                }
            }
        }

        // Physical free actuators (arch-6): wire each free-actuator bounded var to
        // its component's generic `actuator` CalibIndices slot. These unknowns live
        // AFTER the saturated block. No system residual is emitted for them — the
        // owning component supplies the closing equation in its own residual block.
        let free_actuator_base = self.total_state_len
            + self.inverse_control.mapping_count()
            + self.coupling_residual_count()
            + 2 * self.saturated_controllers.len();
        for (index, actuator_id) in self.free_actuators.iter().enumerate() {
            if let Some(bounded_var) = self.bounded_variables.get(actuator_id) {
                if let Some(comp_id) = bounded_var.component_id() {
                    let indices = comp_calib_indices.entry(comp_id.to_string()).or_default();
                    indices.actuator = Some(free_actuator_base + index);
                }
            }
        }

        // Now mutate each node weight (component) with the gathered context.
        for (node_idx, offset, incident, port_edges) in node_context {
            if let Some(component) = self.graph.node_weight_mut(node_idx) {
                component.set_system_context(offset, &incident);
                component.set_port_context(&port_edges);

                // If we registered a name for this node, check if we have calib indices for it
                if let Some((name, _)) = self.component_names.iter().find(|(_, &n)| n == node_idx) {
                    if let Some(&indices) = comp_calib_indices.get(name) {
                        component.set_calib_indices(indices);
                    }
                }
            }
        }

        // Persist the per-component calib map so external readers (result
        // extraction) can surface solved Z-factor values (z_ua, z_dp, …).
        self.calib_indices_by_name = comp_calib_indices;

        // Wire physical thermal couplings (Story 3.4 completion): each coupling
        // owns one state unknown Q [W] at `coupling_state_index(i)`. The
        // cold-side receiver component (e.g. `ThermalLoad`) reads Q in its
        // energy balance via `set_external_heat_index`; the coupling residual
        // closes Q against the hot component's measured duty (see
        // `compute_residuals`).
        let coupling_wirings: Vec<(String, usize)> = self
            .thermal_couplings
            .iter()
            .enumerate()
            .filter_map(|(i, c)| {
                c.cold_component
                    .clone()
                    .map(|name| (name, self.coupling_state_index(i)))
            })
            .collect();
        for (name, q_idx) in coupling_wirings {
            if let Some(&node) = self.component_names.get(&name) {
                if let Some(component) = self.graph.node_weight_mut(node) {
                    component.set_external_heat_index(q_idx);
                }
            } else {
                tracing::warn!(
                    cold_component = %name,
                    "Thermal coupling references an unregistered cold component; \
                     the coupled heat will not be injected"
                );
            }
        }
        for coupling in &self.thermal_couplings {
            if let Some(hot) = &coupling.hot_component {
                if !self.component_names.contains_key(hot) {
                    tracing::warn!(
                        hot_component = %hot,
                        "Thermal coupling references an unregistered hot component; \
                         the transferred duty will be zero"
                    );
                }
            }
        }

        if !self.constraints.is_empty() {
            match self.validate_inverse_control_dof() {
                Ok(()) => {
                    tracing::debug!(
                        constraint_count = self.constraints.len(),
                        control_count = self.inverse_control.mapping_count(),
                        "Inverse control DoF validation passed"
                    );
                }
                Err(DoFError::UnderConstrainedSystem { .. }) => {
                    tracing::warn!(
                        constraint_count = self.constraints.len(),
                        control_count = self.inverse_control.mapping_count(),
                        "Under-constrained inverse control system - solver may still converge"
                    );
                }
                Err(DoFError::OverConstrainedSystem {
                    constraint_count,
                    control_count,
                    equation_count,
                    unknown_count,
                }) => {
                    tracing::warn!(
                        constraint_count,
                        control_count,
                        equation_count,
                        unknown_count,
                        "Over-constrained inverse control system - add more control variables or remove constraints"
                    );
                }
                Err(e) => {
                    tracing::warn!("Inverse control DoF validation error: {}", e);
                }
            }
        }

        // DoF gate:
        // - Over-constrained systems are always rejected (no free lunch: every
        //   extra residual needs a free unknown).
        // - Under-constrained systems log a warning here so pure topology unit
        //   tests with mock components can still exercise graph plumbing.
        //   Production paths (CLI run, SystemBuilder) must call
        //   `validate_system_dof()` which requires a fully square system.
        let report = self.dof_report();
        match report.balance {
            SystemDofBalance::Balanced => {
                tracing::debug!("{}", report.summary());
            }
            SystemDofBalance::UnderConstrained { free_dofs } => {
                tracing::warn!(free_dofs, "{}", report.summary());
            }
            SystemDofBalance::OverConstrained { excess_equations } => {
                tracing::error!(excess_equations, "{}", report.summary());
                if self.enforce_dof_gate {
                    return Err(TopologyError::DofImbalance {
                        message: report.summary(),
                    });
                }
                tracing::warn!(
                    "DoF gate disabled (set_enforce_dof_gate(false)); over-constrained system accepted for tests only"
                );
            }
        }

        Ok(())
    }

    /// Validates the topology: no isolated nodes and edge circuit consistency.
    ///
    /// Note: "All ports connected" validation requires port→edge association
    /// (Story 3.2 Port Compatibility Validation).
    fn validate_topology(&self) -> Result<(), TopologyError> {
        let node_count = self.graph.node_count();
        if node_count == 0 {
            return Ok(());
        }

        for node_idx in self.graph.node_indices() {
            let degree = self
                .graph
                .edges_directed(node_idx, petgraph::Direction::Incoming)
                .count()
                + self
                    .graph
                    .edges_directed(node_idx, petgraph::Direction::Outgoing)
                    .count();
            if degree == 0 {
                return Err(TopologyError::IsolatedNode {
                    node_index: node_idx.index(),
                });
            }
        }

        // Validate that all edges connect nodes within the same circuit
        for edge_idx in self.graph.edge_indices() {
            if let Some((src, tgt)) = self.graph.edge_endpoints(edge_idx) {
                let src_circuit = self.node_circuit(src);
                let tgt_circuit = self.node_circuit(tgt);
                if src_circuit != tgt_circuit {
                    return Err(TopologyError::CrossCircuitConnection {
                        source_circuit: src_circuit.0,
                        target_circuit: tgt_circuit.0,
                    });
                }
            }
        }

        Ok(())
    }

    /// Returns the documented state vector layout.
    ///
    /// Layout: `[P_edge0, h_edge0, P_edge1, h_edge1, ...]` where each edge (in
    /// graph iteration order) contributes 2 entries: pressure (Pa) then enthalpy (J/kg).
    ///
    /// # Panics
    ///
    /// Panics if `finalize()` has not been called.
    pub fn state_layout(&self) -> &'static str {
        assert!(self.finalized, "call finalize() before state_layout()");
        "[m_edge0, P_edge0, h_edge0, m_edge1, P_edge1, h_edge1, ...] — 3 per edge (mass flow kg/s, pressure Pa, enthalpy J/kg)"
    }

    /// Returns the length of the physical state vector.
    ///
    /// After CM1.4 topology presolve the layout is:
    /// `|B| (branch ṁ slots) + 2 × |E| (per-edge P and h) + internal_component_state`.
    ///
    /// Note: This returns the physical state vector length. For the full solver state vector
    /// including control variables, use [`full_state_vector_len`](Self::full_state_vector_len).
    ///
    /// # Panics
    ///
    /// Panics if `finalize()` has not been called.
    pub fn state_vector_len(&self) -> usize {
        assert!(self.finalized, "call finalize() before state_vector_len()");
        self.total_state_len
    }

    /// Returns the state indices (P, h) for the given edge.
    ///
    /// # Returns
    ///
    /// `(state_index_p, state_index_h)` for the edge. The mass-flow index is
    /// available via [`edge_mass_flow_index`](Self::edge_mass_flow_index) or
    /// [`edge_state_indices_full`](Self::edge_state_indices_full).
    ///
    /// # Panics
    ///
    /// Panics if `finalize()` has not been called or if `edge_id` is invalid.
    pub fn edge_state_indices(&self, edge_id: EdgeIndex) -> (usize, usize) {
        assert!(
            self.finalized,
            "call finalize() before edge_state_indices()"
        );
        let (_m, p, h) = *self
            .edge_to_state
            .get(&edge_id)
            .expect("invalid edge index");
        (p, h)
    }

    /// Returns the full state indices `(ṁ, P, h)` for the given edge.
    ///
    /// # Panics
    ///
    /// Panics if `finalize()` has not been called or if `edge_id` is invalid.
    pub fn edge_state_indices_full(&self, edge_id: EdgeIndex) -> (usize, usize, usize) {
        assert!(
            self.finalized,
            "call finalize() before edge_state_indices_full()"
        );
        *self
            .edge_to_state
            .get(&edge_id)
            .expect("invalid edge index")
    }

    /// Returns the mass-flow (ṁ) state index for the given edge.
    ///
    /// # Panics
    ///
    /// Panics if `finalize()` has not been called or if `edge_id` is invalid.
    pub fn edge_mass_flow_index(&self, edge_id: EdgeIndex) -> usize {
        self.edge_state_indices_full(edge_id).0
    }

    /// Returns the [`EdgeKind`] of the given edge.
    ///
    /// # Panics
    ///
    /// Panics if `edge_id` is invalid.
    pub fn edge_kind(&self, edge_id: EdgeIndex) -> EdgeKind {
        self.graph
            .edge_weight(edge_id)
            .map(|w| w.kind)
            .expect("invalid edge index")
    }

    /// Number of per-edge mass-flow closure equations.
    ///
    /// Returns 0 since CM1.3: ṁ residuals are now contributed by individual
    /// components (`compute_residuals`), so no extra closure rows are needed.
    /// Kept for API compatibility; callers should eventually remove their
    /// `+ system.mass_flow_closure_count()` calls.
    pub fn mass_flow_closure_count(&self) -> usize {
        0
    }

    /// Returns the number of edges in the graph.
    pub fn edge_count(&self) -> usize {
        self.graph.edge_count()
    }

    /// Returns an iterator over all edge indices in the graph.
    pub fn edge_indices(&self) -> impl Iterator<Item = EdgeIndex> + '_ {
        self.graph.edge_indices()
    }

    /// Returns the source and target node indices for the given edge.
    ///
    /// Returns `None` if the edge index is invalid.
    pub fn edge_endpoints(&self, edge: EdgeIndex) -> Option<(NodeIndex, NodeIndex)> {
        self.graph.edge_endpoints(edge)
    }

    /// Returns a reference to the internal graph.
    pub fn graph(&self) -> &Graph<Box<dyn Component>, FlowEdge, Directed> {
        &self.graph
    }

    /// Returns a reference to the node-to-circuit mapping.
    pub fn node_to_circuit(&self) -> &HashMap<NodeIndex, CircuitId> {
        &self.node_to_circuit
    }

    /// Returns a reference to the constraints map.
    pub fn constraints_map(&self) -> &HashMap<ConstraintId, Constraint> {
        &self.constraints
    }

    /// Returns a reference to the bounded variables map.
    pub fn bounded_variables_map(&self) -> &HashMap<BoundedVariableId, BoundedVariable> {
        &self.bounded_variables
    }

    /// Returns the number of nodes (components) in the graph.
    pub fn node_count(&self) -> usize {
        self.graph.node_count()
    }

    /// Returns the number of distinct circuits in the machine.
    ///
    /// Circuits are identified by the unique circuit IDs present in `node_to_circuit`.
    /// Empty system returns 0. Systems with components always return >= 1 since
    /// all components are assigned to a circuit (defaulting to circuit 0).
    /// Valid circuit IDs are 0 through 4 (inclusive), supporting up to 5 circuits.
    pub fn circuit_count(&self) -> usize {
        if self.graph.node_count() == 0 {
            return 0;
        }
        let mut ids: Vec<u16> = self.node_to_circuit.values().map(|c| c.0).collect();
        if ids.is_empty() {
            // This shouldn't happen since add_component adds to node_to_circuit,
            // but handle defensively
            return 1;
        }
        ids.sort_unstable();
        ids.dedup();
        ids.len()
    }

    /// Returns an iterator over node indices belonging to the given circuit.
    pub fn circuit_nodes(&self, circuit_id: CircuitId) -> impl Iterator<Item = NodeIndex> + '_ {
        self.graph.node_indices().filter(move |&idx| {
            self.node_to_circuit
                .get(&idx)
                .copied()
                .unwrap_or(CircuitId::ZERO)
                == circuit_id
        })
    }

    /// Returns an iterator over edge indices belonging to the given circuit.
    ///
    /// An edge belongs to a circuit if both its source and target nodes are in that circuit.
    pub fn circuit_edges(&self, circuit_id: CircuitId) -> impl Iterator<Item = EdgeIndex> + '_ {
        self.graph.edge_indices().filter(move |&edge_idx| {
            let (src, tgt) = self.graph.edge_endpoints(edge_idx).expect("valid edge");
            self.node_circuit(src) == circuit_id && self.node_circuit(tgt) == circuit_id
        })
    }

    /// Checks if a circuit has any components.
    fn circuit_exists(&self, circuit_id: CircuitId) -> bool {
        self.node_to_circuit.values().any(|&c| c == circuit_id)
    }

    /// Adds a thermal coupling between two circuits.
    ///
    /// Thermal couplings represent heat exchangers that transfer heat from a "hot"
    /// circuit to a "cold" circuit without fluid mixing. Heat flows from hot to cold
    /// proportional to the temperature difference and thermal conductance (UA).
    ///
    /// # Arguments
    ///
    /// * `coupling` - The thermal coupling to add
    ///
    /// # Returns
    ///
    /// The index of the added coupling in the internal storage.
    ///
    /// # Errors
    ///
    /// Returns `TopologyError::InvalidCircuitForCoupling` if either circuit
    /// referenced in the coupling does not exist in the system.
    ///
    /// # Example
    ///
    /// ```no_run
    /// use entropyk_solver::{System, ThermalCoupling, CircuitId};
    /// use entropyk_core::ThermalConductance;
    /// use entropyk_components::Component;
    /// # fn make_mock() -> Box<dyn Component> { unimplemented!() }
    ///
    /// let mut sys = System::new();
    /// sys.add_component_to_circuit(make_mock(), CircuitId(0)).unwrap();
    /// sys.add_component_to_circuit(make_mock(), CircuitId(1)).unwrap();
    ///
    /// let coupling = ThermalCoupling::new(
    ///     CircuitId(0),
    ///     CircuitId(1),
    ///     ThermalConductance::from_watts_per_kelvin(1000.0),
    /// );
    /// let idx = sys.add_thermal_coupling(coupling).unwrap();
    /// ```
    pub fn add_thermal_coupling(
        &mut self,
        coupling: ThermalCoupling,
    ) -> Result<usize, TopologyError> {
        // Validate that both circuits exist
        if !self.circuit_exists(coupling.hot_circuit) {
            return Err(TopologyError::InvalidCircuitForCoupling {
                circuit_id: coupling.hot_circuit.0,
            });
        }
        if !self.circuit_exists(coupling.cold_circuit) {
            return Err(TopologyError::InvalidCircuitForCoupling {
                circuit_id: coupling.cold_circuit.0,
            });
        }

        self.finalized = false;
        self.thermal_couplings.push(coupling);
        Ok(self.thermal_couplings.len() - 1)
    }

    /// Returns the number of thermal couplings in the system.
    pub fn thermal_coupling_count(&self) -> usize {
        self.thermal_couplings.len()
    }

    /// Returns a reference to all thermal couplings.
    pub fn thermal_couplings(&self) -> &[ThermalCoupling] {
        &self.thermal_couplings
    }

    /// Returns a reference to a specific thermal coupling by index.
    pub fn get_thermal_coupling(&self, index: usize) -> Option<&ThermalCoupling> {
        self.thermal_couplings.get(index)
    }

    // ────────────────────────────────────────────────────────────────────────
    // Component Name Registry (for constraint validation)
    // ────────────────────────────────────────────────────────────────────────

    /// Registers a human-readable name for a component node.
    ///
    /// This name can be used in constraints to reference the component.
    /// For example, register "evaporator" to reference it in a superheat constraint.
    ///
    /// # Arguments
    ///
    /// * `name` - Human-readable name for the component (e.g., "evaporator", "condenser")
    /// * `node` - The NodeIndex returned from `add_component()` or `add_component_to_circuit()`
    ///
    /// # Returns
    ///
    /// `true` if the name was newly registered, `false` if the name was already in use
    /// (in which case the mapping is updated to the new node).
    ///
    /// # Example
    ///
    /// ```rust,ignore
    /// let mut sys = System::new();
    /// let evap_node = sys.add_component(make_evaporator());
    /// sys.register_component_name("evaporator", evap_node);
    ///
    /// // Now constraints can reference "evaporator"
    /// let constraint = Constraint::new(
    ///     ConstraintId::new("superheat_control"),
    ///     ComponentOutput::Superheat { component_id: "evaporator".to_string() },
    ///     5.0,
    /// );
    /// sys.add_constraint(constraint)?; // Validates "evaporator" exists
    /// ```
    pub fn register_component_name(&mut self, name: &str, node: NodeIndex) -> bool {
        self.component_names
            .insert(name.to_string(), node)
            .is_none()
    }

    /// Returns the NodeIndex for a registered component name, or None if not found.
    pub fn get_component_node(&self, name: &str) -> Option<NodeIndex> {
        self.component_names.get(name).copied()
    }

    /// Returns the names of all registered components.
    pub fn registered_component_names(&self) -> impl Iterator<Item = &str> {
        self.component_names.keys().map(|s| s.as_str())
    }

    /// Per-component `CalibIndices` snapshot captured at `finalize()`.
    /// Each `Some(idx)` slot (z_ua, z_dp, z_flow, z_flow_eco, z_power, z_etav,
    /// actuator) means that factor was promoted to a free solver unknown and
    /// its solved value lives at `state[idx]`. Used by result extraction to
    /// surface solved calibration/control values to the user.
    pub fn calib_indices_by_name(&self) -> &HashMap<String, entropyk_core::CalibIndices> {
        &self.calib_indices_by_name
    }

    /// Returns a reference to the component stored at the given node index.
    ///
    /// # Panics
    ///
    /// Panics if the node index is invalid.
    pub fn component(&self, node: NodeIndex) -> &dyn Component {
        self.graph
            .node_weight(node)
            .expect("invalid node index")
            .as_ref()
    }

    /// Returns a mutable reference to the component at the given node index.
    ///
    /// Returns `None` if the node index is invalid.
    /// Used for post-build injection of fluid backends via the builder.
    pub fn component_mut(&mut self, node: NodeIndex) -> Option<&mut dyn Component> {
        let weight = self.graph.node_weight_mut(node)?;
        Some(weight.as_mut())
    }

    // ────────────────────────────────────────────────────────────────────────
    // Constraint Management (Inverse Control)
    // ────────────────────────────────────────────────────────────────────────

    /// Adds a constraint for inverse control.
    ///
    /// Constraints define desired output conditions. During solving, the solver
    /// will attempt to find input values that satisfy all constraints.
    ///
    /// # Arguments
    ///
    /// * `constraint` - The constraint to add
    ///
    /// # Errors
    ///
    /// Returns `ConstraintError::DuplicateId` if a constraint with the same ID
    /// already exists.
    /// Returns `ConstraintError::InvalidReference` if the component referenced
    /// in the constraint's output has not been registered via `register_component_name()`.
    ///
    /// # Example
    ///
    /// ```rust,ignore
    /// use entropyk_solver::inverse::{Constraint, ConstraintId, ComponentOutput};
    ///
    /// let constraint = Constraint::new(
    ///     ConstraintId::new("superheat_control"),
    ///     ComponentOutput::Superheat {
    ///         component_id: "evaporator".to_string()
    ///     },
    ///     5.0,  // target: 5K superheat
    /// );
    ///
    /// system.add_constraint(constraint)?;
    /// ```
    pub fn add_constraint(&mut self, constraint: Constraint) -> Result<(), ConstraintError> {
        let id = constraint.id().clone();
        if self.constraints.contains_key(&id) {
            return Err(ConstraintError::DuplicateId { id });
        }

        // AC2: Validate that the component referenced in the constraint exists
        let component_id = constraint.output().component_id();
        if !self.component_names.contains_key(component_id) {
            return Err(ConstraintError::InvalidReference {
                component_id: component_id.to_string(),
            });
        }

        self.constraints.insert(id, constraint);
        Ok(())
    }

    /// Removes a constraint by ID.
    ///
    /// # Arguments
    ///
    /// * `id` - The constraint identifier to remove
    ///
    /// # Returns
    ///
    /// The removed constraint, or `None` if no constraint with that ID exists.
    pub fn remove_constraint(&mut self, id: &ConstraintId) -> Option<Constraint> {
        self.inverse_control.unlink_constraint(id);
        self.constraints.remove(id)
    }

    /// Returns the number of constraints in the system.
    pub fn constraint_count(&self) -> usize {
        self.constraints.len()
    }

    /// Returns a reference to all constraints.
    pub fn constraints(&self) -> impl Iterator<Item = &Constraint> {
        self.constraints.values()
    }

    /// Returns a reference to a specific constraint by ID.
    pub fn get_constraint(&self, id: &ConstraintId) -> Option<&Constraint> {
        self.constraints.get(id)
    }

    /// Returns the number of constraint residual equations.
    ///
    /// Each constraint adds one equation to the residual vector.
    pub fn constraint_residual_count(&self) -> usize {
        self.constraints.len()
    }

    /// Computes constraint residuals and appends them to the provided vector.
    ///
    /// This method computes the residual for each constraint:
    ///
    /// $$r_{constraint} = f(x) - y_{target}$$
    ///
    /// where $f(x)$ is the measured output value and $y_{target}$ is the constraint target.
    ///
    /// # Arguments
    ///
    /// * `state` - Current system state (edge pressures and enthalpies)
    /// * `residuals` - Residual vector to append constraint residuals to
    /// * `measured_values` - Map of constraint IDs to their measured values (from component state)
    ///
    /// # Returns
    ///
    /// `Ok(count)` where count is the number of constraint residuals added.
    ///
    /// # Errors
    ///
    /// Returns `ConstraintError::UnmeasuredConstraint` if a constraint references a component
    /// with no measured value (not registered or no associated edges).
    /// Returns `ConstraintError::ResidualSliceTooShort` if the residual slice is too short.
    ///
    /// # Example
    ///
    /// ```rust,ignore
    /// let mut residuals = ResidualVector::new();
    /// let measured = system.extract_constraint_values_with_controls(&state, &control);
    /// let count = system.compute_constraint_residuals(&state, &mut residuals, &measured);
    /// ```
    pub fn compute_constraint_residuals(
        &self,
        _state: &StateSlice,
        residuals: &mut [f64],
        measured_values: &HashMap<ConstraintId, f64>,
    ) -> Result<usize, ConstraintError> {
        if self.constraints.is_empty() {
            return Ok(0);
        }

        let mut count = 0;
        for constraint in self.constraints.values() {
            let measured = match measured_values.get(constraint.id()).copied() {
                Some(v) => v,
                None => {
                    return Err(ConstraintError::UnmeasuredConstraint {
                        constraint_id: constraint.id().to_string(),
                        component_id: constraint.output().component_id().to_string(),
                    });
                }
            };
            let residual = constraint.compute_residual(measured);
            if count >= residuals.len() {
                return Err(ConstraintError::ResidualSliceTooShort {
                    index: count,
                    len: residuals.len(),
                    required: self.constraints.len(),
                });
            }
            residuals[count] = residual;
            count += 1;
        }
        Ok(count)
    }

    /// Extracts measured values for all constraints, incorporating control variable effects.
    ///
    /// This method computes the measured output value for each constraint, taking into
    /// account the current state and control variable values. For MIMO (Multi-Input
    /// Multi-Output) systems, ALL control variables can affect ALL constraint outputs
    /// due to system coupling.
    ///
    /// # Arguments
    ///
    /// * `state` - Current system state (edge pressures and enthalpies)
    /// * `control_values` - Current values of control variables
    ///
    /// # Returns
    ///
    /// A map from constraint ID to measured output value.
    ///
    /// # Cross-Coupling for MIMO Systems
    ///
    /// In a real thermodynamic system, control variables are coupled:
    /// - Compressor speed affects both capacity AND superheat
    /// - Valve opening affects both superheat AND capacity
    ///
    /// The mock implementation simulates this coupling for Jacobian cross-derivative
    /// computation. Each control variable has a primary effect (on its linked constraint)
    /// and a secondary effect (on other constraints) to simulate thermal coupling.
    /// Maps a solver-side [`ComponentOutput`] to the components-crate
    /// [`MeasuredOutput`] used by [`Component::measure_output`]. Returns `None`
    /// for outputs that have no real-measurement mapping yet (falls back to the
    /// legacy placeholder path).
    fn map_output_kind(
        output: &crate::inverse::ComponentOutput,
    ) -> Option<entropyk_components::MeasuredOutput> {
        use crate::inverse::ComponentOutput as CO;
        use entropyk_components::MeasuredOutput as MO;
        Some(match output {
            CO::Superheat { .. } => MO::Superheat,
            CO::Subcooling { .. } => MO::Subcooling,
            CO::Capacity { .. } => MO::Capacity,
            CO::HeatTransferRate { .. } => MO::HeatTransferRate,
            CO::MassFlowRate { .. } => MO::MassFlowRate,
            CO::Pressure { .. } => MO::Pressure,
            CO::Temperature { .. } => MO::Temperature,
            CO::SaturationTemperature { .. } => MO::SaturationTemperature,
        })
    }

    pub fn extract_constraint_values_with_controls(
        &self,
        state: &StateSlice,
        control_values: &[f64],
    ) -> HashMap<ConstraintId, f64> {
        let mut measured = HashMap::new();
        if self.constraints.is_empty() {
            return measured;
        }

        // Build a map of control variable index -> component_id it controls
        // This uses the proper component_id() field from BoundedVariable
        let mut control_to_component: HashMap<usize, &str> = HashMap::new();
        for (j, bounded_var_id) in self.inverse_control.linked_controls().enumerate() {
            if let Some(bounded_var) = self.bounded_variables.get(bounded_var_id) {
                if let Some(comp_id) = bounded_var.component_id() {
                    control_to_component.insert(j, comp_id);
                }
            }
        }

        for constraint in self.constraints.values() {
            let comp_id = constraint.output().component_id();
            if let Some(&node_idx) = self.component_names.get(comp_id) {
                // Prefer a REAL physical measurement from the component itself.
                // This replaces the placeholder formulas below with genuine
                // thermodynamics (superheat, subcooling, capacity, …). The mock
                // path is only used for components/outputs that do not (yet)
                // implement `measure_output` (e.g. synthetic test doubles).
                if let Some(kind) = Self::map_output_kind(constraint.output()) {
                    if let Some(component) = self.graph.node_weight(node_idx) {
                        if let Some(real) = component.measure_output(kind, state) {
                            if real.is_finite() {
                                measured.insert(constraint.id().clone(), real);
                                continue;
                            }
                        }
                    }
                }

                // Find first associated edge (incoming or outgoing)
                let mut edge_opt = self
                    .graph
                    .edges_directed(node_idx, petgraph::Direction::Incoming)
                    .next();
                if edge_opt.is_none() {
                    edge_opt = self
                        .graph
                        .edges_directed(node_idx, petgraph::Direction::Outgoing)
                        .next();
                }

                if let Some(edge) = edge_opt {
                    if let Some(&(_m, p_idx, h_idx)) = self.edge_to_state.get(&edge.id()) {
                        let mut value = match constraint.output() {
                            crate::inverse::ComponentOutput::Pressure { .. } => state[p_idx],
                            crate::inverse::ComponentOutput::Temperature { .. } => 300.0, // Mock for MVP without fluid backend
                            crate::inverse::ComponentOutput::Superheat { .. } => {
                                // Mock numerical value sensitive to BOTH P and h for Jacobian calculation
                                state[h_idx] / 1000.0 - (state[p_idx] / 1e5)
                            }
                            crate::inverse::ComponentOutput::Subcooling { .. } => {
                                (state[p_idx] / 1e5) - state[h_idx] / 1000.0
                            }
                            crate::inverse::ComponentOutput::Capacity { .. } => {
                                // Mock capacity: h * mass_flow. Let's just use h for Jacobian sensitivity
                                state[h_idx] * 10.0
                            }
                            _ => 0.0,
                        };

                        // MIMO Cross-Coupling: ALL control variables can affect ALL constraints
                        // In a real system, changing compressor speed affects both capacity and superheat,
                        // and changing valve opening also affects both. We simulate this coupling here.
                        //
                        // ⚠️ MOCK COEFFICIENTS: These values (10.0, 2.0) are placeholders for testing.
                        // They create a well-conditioned Jacobian with off-diagonal entries that allow
                        // Newton-Raphson to converge. Real implementations should replace these with
                        // actual component physics derived from:
                        // - Component characteristic curves (compressor map, valve Cv curve)
                        // - Thermodynamic property calculations via fluid backend
                        // - Energy and mass balance equations
                        //
                        // The 5:1 ratio between primary and secondary effects is arbitrary but creates
                        // a diagonally-dominant Jacobian that converges reliably. See Story 5.4
                        // Review Follow-ups for tracking real thermodynamics integration.
                        //
                        // For each control variable:
                        // - Primary effect (10.0): if control is linked to this constraint's component
                        // - Secondary effect (2.0): cross-coupling to other constraints
                        const MIMO_PRIMARY_COEFF: f64 = 10.0;
                        const MIMO_SECONDARY_COEFF: f64 = 2.0;

                        for (j, _bounded_var_id) in
                            self.inverse_control.linked_controls().enumerate()
                        {
                            if j >= control_values.len() {
                                continue;
                            }
                            let ctrl_val = control_values[j];

                            // Check if this control variable is primarily associated with this component
                            let is_primary = control_to_component
                                .get(&j)
                                .map_or(false, |&c| c == comp_id);

                            if is_primary {
                                // Primary effect: strong influence on the controlled output
                                // e.g., valve opening strongly affects superheat
                                value += ctrl_val * MIMO_PRIMARY_COEFF;
                            } else {
                                // Secondary (cross-coupling) effect: weaker influence
                                // e.g., compressor speed also affects superheat (through mass flow)
                                // This creates the off-diagonal entries in the MIMO Jacobian
                                value += ctrl_val * MIMO_SECONDARY_COEFF;
                            }
                        }

                        if value.is_nan() {
                            tracing::warn!(
                                constraint_id = constraint.id().as_str(),
                                "NaN detected in constraint output for component '{}', skipping insert",
                                constraint.output().component_id()
                            );
                        } else {
                            measured.insert(constraint.id().clone(), value);
                        }
                    }
                }
            }
        }
        measured
    }

    /// Computes the Jacobian entries for inverse control constraints.
    ///
    /// For each constraint→control mapping, adds ∂r/∂x entries to the Jacobian:
    ///
    /// $$\frac{\partial r}{\partial x_{control}} = \frac{\partial (measured - target)}{\partial x_{control}} = \frac{\partial measured}{\partial x_{control}}$$
    ///
    /// # Arguments
    ///
    /// * `state` - Current system state
    /// * `row_offset` - Starting row index for constraint equations in the Jacobian
    /// * `control_values` - Current values of control variables (for finite difference)
    ///
    /// # Returns
    ///
    /// A vector of `(row, col, value)` tuples representing Jacobian entries.
    ///
    /// # Note
    ///
    /// MVP uses finite difference approximation. Future versions may use analytical
    /// derivatives from components for better accuracy and performance.
    ///
    /// # Finite Difference Epsilon
    ///
    /// Uses the epsilon configured in `InverseControlConfig` (default 1e-6) for central
    /// finite differences. Configure via `set_inverse_control_epsilon()`.
    pub fn compute_inverse_control_jacobian(
        &self,
        state: &StateSlice,
        row_offset: usize,
        control_values: &[f64],
    ) -> Vec<(usize, usize, f64)> {
        let mut entries = Vec::new();

        if self.inverse_control.mapping_count() == 0 {
            return entries;
        }

        if control_values.len() < self.inverse_control.mapping_count() {
            tracing::error!(
                provided = control_values.len(),
                required = self.inverse_control.mapping_count(),
                "control_values too short for Jacobian computation"
            );
            return entries;
        }

        // Use configurable epsilon from InverseControlConfig
        let eps = self.inverse_control.finite_diff_epsilon();
        if state.len() < self.total_state_len {
            tracing::error!(
                state_len = state.len(),
                required = self.total_state_len,
                "compute_inverse_control_jacobian: state slice too short, returning empty"
            );
            return entries;
        }
        let mut state_mut = state.to_vec();
        let mut control_mut = control_values.to_vec();

        // 1. Compute ∂r_i / ∂x_j (Partial derivatives with respect to PHYSICAL states P, h)
        // We do this per constraint to keep perturbations localized where possible
        for (i, (constraint_id, _)) in self.inverse_control.mappings().enumerate() {
            let row = row_offset + i;
            if let Some(constraint) = self.constraints.get(constraint_id) {
                let comp_id = constraint.output().component_id();

                if let Some(&node_idx) = self.component_names.get(comp_id) {
                    let mut state_indices = Vec::new();
                    // Gather all edge state indices for this component. The
                    // mass-flow index is included because measured outputs such
                    // as Capacity (q = ṁ·Δh) depend directly on ṁ; omitting it
                    // decouples the constraint from any mass-flow actuator (f_m)
                    // and makes the inverse-control Jacobian singular.
                    for edge in self
                        .graph
                        .edges_directed(node_idx, petgraph::Direction::Incoming)
                    {
                        if let Some(&(m_idx, p_idx, h_idx)) = self.edge_to_state.get(&edge.id()) {
                            if !state_indices.contains(&m_idx) {
                                state_indices.push(m_idx);
                            }
                            if !state_indices.contains(&p_idx) {
                                state_indices.push(p_idx);
                            }
                            if !state_indices.contains(&h_idx) {
                                state_indices.push(h_idx);
                            }
                        }
                    }
                    for edge in self
                        .graph
                        .edges_directed(node_idx, petgraph::Direction::Outgoing)
                    {
                        if let Some(&(m_idx, p_idx, h_idx)) = self.edge_to_state.get(&edge.id()) {
                            if !state_indices.contains(&m_idx) {
                                state_indices.push(m_idx);
                            }
                            if !state_indices.contains(&p_idx) {
                                state_indices.push(p_idx);
                            }
                            if !state_indices.contains(&h_idx) {
                                state_indices.push(h_idx);
                            }
                        }
                    }

                    // Central finite difference for Jacobian entries w.r.t physical state
                    for &col in &state_indices {
                        let orig = state_mut[col];

                        state_mut[col] = orig + eps;
                        let plus = self
                            .extract_constraint_values_with_controls(&state_mut, control_values);
                        let val_plus = plus.get(constraint_id).copied().unwrap_or(0.0);

                        state_mut[col] = orig - eps;
                        let minus = self
                            .extract_constraint_values_with_controls(&state_mut, control_values);
                        let val_minus = minus.get(constraint_id).copied().unwrap_or(0.0);

                        state_mut[col] = orig; // Restore

                        let derivative = (val_plus - val_minus) / (2.0 * eps);
                        if derivative.abs() > 1e-10 {
                            entries.push((row, col, derivative));
                            tracing::trace!(
                                constraint = constraint_id.as_str(),
                                row,
                                col,
                                derivative,
                                "Inverse control Jacobian actual ∂r/∂state entry"
                            );
                        }
                    }
                }
            }
        }

        // 2. Compute ∂r_i / ∂u_j (Cross-derivatives with respect to CONTROL variables)
        // Here we must form the full dense block because control variable 'j' could affect constraint 'i'
        // even if they are not explicitly linked, due to system coupling.
        let control_offset = self.state_vector_len();

        for (j, (_, bounded_var_id)) in self.inverse_control.mappings().enumerate() {
            let col = control_offset + j;
            let orig = control_mut[j];

            // Perturb control variable +eps
            control_mut[j] = orig + eps;
            let plus = self.extract_constraint_values_with_controls(state, &control_mut);

            // Perturb control variable -eps
            control_mut[j] = orig - eps;
            let minus = self.extract_constraint_values_with_controls(state, &control_mut);

            control_mut[j] = orig; // Restore

            // For this perturbed control variable j, compute the effect on ALL constraints i
            for (i, (constraint_id, _)) in self.inverse_control.mappings().enumerate() {
                let row = row_offset + i;

                let val_plus = plus.get(constraint_id).copied().unwrap_or(0.0);
                let val_minus = minus.get(constraint_id).copied().unwrap_or(0.0);
                let derivative = (val_plus - val_minus) / (2.0 * eps);

                // We add it even if it's 0 to maintain block structure (optional but safe)
                // However, for performance we only add non-zeros
                if derivative.abs() > 1e-10 {
                    entries.push((row, col, derivative));
                    tracing::trace!(
                        constraint = ?constraint_id,
                        control = ?bounded_var_id,
                        row, col, derivative,
                        "Inverse control Jacobian cross-derivative ∂r/∂u entry"
                    );
                }
            }
        }

        entries
    }

    fn incident_state_indices_for_component(&self, node_idx: NodeIndex) -> Vec<usize> {
        let mut state_indices = Vec::new();
        for edge in self
            .graph
            .edges_directed(node_idx, petgraph::Direction::Incoming)
        {
            if let Some(&(m_idx, p_idx, h_idx)) = self.edge_to_state.get(&edge.id()) {
                for idx in [m_idx, p_idx, h_idx] {
                    if !state_indices.contains(&idx) {
                        state_indices.push(idx);
                    }
                }
            }
        }
        for edge in self
            .graph
            .edges_directed(node_idx, petgraph::Direction::Outgoing)
        {
            if let Some(&(m_idx, p_idx, h_idx)) = self.edge_to_state.get(&edge.id()) {
                for idx in [m_idx, p_idx, h_idx] {
                    if !state_indices.contains(&idx) {
                        state_indices.push(idx);
                    }
                }
            }
        }
        state_indices
    }

    fn measure_saturated_output(
        &self,
        output: &crate::inverse::ComponentOutput,
        state: &StateSlice,
    ) -> f64 {
        let comp_id = output.component_id();
        if let Some(&node_idx) = self.component_names.get(comp_id) {
            if let Some(kind) = Self::map_output_kind(output) {
                if let Some(component) = self.graph.node_weight(node_idx) {
                    if let Some(value) = component.measure_output(kind, state) {
                        if value.is_finite() {
                            return value;
                        }
                    }
                }
            }
        }

        tracing::warn!(
            component_id = comp_id,
            output = output.constraint_type_name(),
            "Saturated controller output measurement unavailable or non-finite; using synthetic 0.0 fallback"
        );
        0.0
    }

    fn compute_saturated_control_jacobian(
        &self,
        state: &StateSlice,
        row_offset: usize,
    ) -> Vec<(usize, usize, f64)> {
        let mut entries = Vec::new();
        if self.saturated_controllers.is_empty() {
            return entries;
        }

        let eps = self.inverse_control.finite_diff_epsilon();
        let mut state_mut = state.to_vec();

        for (i, ctrl) in self.saturated_controllers.iter().enumerate() {
            let u_idx = self.saturated_u_index(i);
            let x_idx = self.saturated_x_index(i);
            if state.len() <= x_idx {
                tracing::error!(
                    state_len = state.len(),
                    required = x_idx + 1,
                    "compute_saturated_control_jacobian: state slice too short"
                );
                continue;
            }

            let row_u = row_offset + 2 * i;
            let row_y = row_u + 1;
            let x = state[x_idx];
            let (dru_du, dru_dx) = ctrl.d_residual_u(x);
            let (_dry_dy, dry_dx) = ctrl.d_residual_y(x);
            entries.push((row_u, u_idx, dru_du));
            entries.push((row_u, x_idx, dru_dx));
            entries.push((row_y, x_idx, dry_dx));

            if ctrl.is_network() {
                // Override network: r_y = E − (x − S(x)); the plant coupling is
                // ∂E/∂col = Σ_i w_i·(−gain_i)·(∂measurement_i/∂col), where w_i are
                // the analytic softMin/softMax selector weights.
                let objs = ctrl.objectives();
                let measured: Vec<f64> = objs
                    .iter()
                    .map(|o| self.measure_saturated_output(&o.output, &state_mut))
                    .collect();
                let weights = ctrl.network_error_weights(&measured);
                let mut col_acc: HashMap<usize, f64> = HashMap::new();
                for (oi, o) in objs.iter().enumerate() {
                    if let Some(&node_idx) = self.component_names.get(o.output.component_id()) {
                        for col in self.incident_state_indices_for_component(node_idx) {
                            let orig = state_mut[col];
                            state_mut[col] = orig + eps;
                            let y_plus = self.measure_saturated_output(&o.output, &state_mut);
                            state_mut[col] = orig - eps;
                            let y_minus = self.measure_saturated_output(&o.output, &state_mut);
                            state_mut[col] = orig;
                            let dy_dcol = (y_plus - y_minus) / (2.0 * eps);
                            *col_acc.entry(col).or_insert(0.0) += weights[oi] * (-o.gain) * dy_dcol;
                        }
                    }
                }
                for (col, derivative) in col_acc {
                    if derivative.abs() > 1e-10 {
                        entries.push((row_y, col, derivative));
                    }
                }
            } else {
                let comp_id = ctrl.output().component_id();
                if let Some(&node_idx) = self.component_names.get(comp_id) {
                    for col in self.incident_state_indices_for_component(node_idx) {
                        let orig = state_mut[col];
                        state_mut[col] = orig + eps;
                        let y_plus = self.measure_saturated_output(ctrl.output(), &state_mut);
                        state_mut[col] = orig - eps;
                        let y_minus = self.measure_saturated_output(ctrl.output(), &state_mut);
                        state_mut[col] = orig;

                        let dy_dcol = (y_plus - y_minus) / (2.0 * eps);
                        let derivative = -ctrl.gain() * dy_dcol;
                        if derivative.abs() > 1e-10 {
                            entries.push((row_y, col, derivative));
                        }
                    }
                }
            }
        }

        entries
    }

    // ────────────────────────────────────────────────────────────────────────
    // Bounded Variable Management (Inverse Control)
    // ────────────────────────────────────────────────────────────────────────

    /// Adds a bounded control variable for inverse control.
    ///
    /// Bounded variables ensure Newton steps stay within physical limits
    /// (e.g., valve position 0.0 to 1.0, VFD speed 0.3 to 1.0).
    ///
    /// # Arguments
    ///
    /// * `variable` - The bounded variable to add
    ///
    /// # Errors
    ///
    /// Returns `BoundedVariableError::DuplicateId` if a variable with the same ID
    /// already exists.
    /// Returns `BoundedVariableError::InvalidComponent` if the variable is associated
    /// with a component that has not been registered via `register_component_name()`.
    ///
    /// # Example
    ///
    /// ```rust,ignore
    /// use entropyk_solver::inverse::{BoundedVariable, BoundedVariableId};
    ///
    /// let valve = BoundedVariable::new(
    ///     BoundedVariableId::new("expansion_valve"),
    ///     0.5,   // initial: 50% open
    ///     0.0,   // min: fully closed
    ///     1.0,   // max: fully open
    /// )?;
    ///
    /// system.add_bounded_variable(valve)?;
    /// ```
    pub fn add_bounded_variable(
        &mut self,
        variable: BoundedVariable,
    ) -> Result<(), BoundedVariableError> {
        let id = variable.id().clone();
        if self.bounded_variables.contains_key(&id) {
            return Err(BoundedVariableError::DuplicateId { id });
        }

        // Validate that the component referenced in the variable exists (if any)
        if let Some(component_id) = variable.component_id() {
            if !self.component_names.contains_key(component_id) {
                return Err(BoundedVariableError::InvalidComponent {
                    component_id: component_id.to_string(),
                });
            }
        }

        self.bounded_variables.insert(id, variable);
        Ok(())
    }

    /// Removes a bounded variable by ID.
    ///
    /// # Arguments
    ///
    /// * `id` - The bounded variable identifier to remove
    ///
    /// # Returns
    ///
    /// The removed variable, or `None` if no variable with that ID exists.
    pub fn remove_bounded_variable(&mut self, id: &BoundedVariableId) -> Option<BoundedVariable> {
        self.inverse_control.unlink_control(id);
        self.free_actuators.retain(|a| a != id);
        self.bounded_variables.remove(id)
    }

    /// Returns the number of bounded variables in the system.
    pub fn bounded_variable_count(&self) -> usize {
        self.bounded_variables.len()
    }

    /// Returns an iterator over all bounded variables.
    pub fn bounded_variables(&self) -> impl Iterator<Item = &BoundedVariable> {
        self.bounded_variables.values()
    }

    /// Returns a reference to a specific bounded variable by ID.
    pub fn get_bounded_variable(&self, id: &BoundedVariableId) -> Option<&BoundedVariable> {
        self.bounded_variables.get(id)
    }

    /// Returns a mutable reference to a specific bounded variable by ID.
    pub fn get_bounded_variable_mut(
        &mut self,
        id: &BoundedVariableId,
    ) -> Option<&mut BoundedVariable> {
        self.bounded_variables.get_mut(id)
    }

    /// Checks if any bounded variables are saturated (at bounds).
    ///
    /// Returns a vector of saturation info for all variables currently at bounds.
    pub fn saturated_variables(&self) -> Vec<crate::inverse::SaturationInfo> {
        self.bounded_variables
            .iter()
            .filter(|(id, _)| !self.free_actuators.contains(id))
            .filter_map(|(_, v)| v.is_saturated())
            .collect()
    }

    // ────────────────────────────────────────────────────────────────────────
    // Inverse Control Mapping (Story 5.3)
    // ────────────────────────────────────────────────────────────────────────

    /// Removes all constraints, bounded variables, and inverse control mappings.
    pub fn clear_inverse_control(&mut self) {
        self.constraints.clear();
        self.bounded_variables.clear();
        self.inverse_control.clear();
        self.saturated_controllers.clear();
        self.free_actuators.clear();
    }

    /// Adds a saturated PI controller loop.
    ///
    /// The actuator bounded variable must already be registered and must not be
    /// linked to a hard inverse-control constraint. The saturated controller
    /// contributes its own two residuals and two unknowns.
    pub fn add_saturated_controller(&mut self, ctrl: SaturatedController) {
        assert!(
            self.bounded_variables.contains_key(ctrl.actuator()),
            "saturated controller actuator '{}' must be registered as a bounded variable",
            ctrl.actuator()
        );
        assert!(
            !self.inverse_control.is_control_linked(ctrl.actuator()),
            "saturated controller actuator '{}' must not be linked as a hard constraint",
            ctrl.actuator()
        );
        self.saturated_controllers.push(ctrl);
    }

    /// Registers a physical free-actuator unknown (arch-6).
    ///
    /// The actuator bounded variable must already be registered via
    /// [`add_bounded_variable`](Self::add_bounded_variable), must reference an
    /// existing component, and must not be linked to a hard inverse-control
    /// constraint or driven by a saturated controller. It appends exactly one
    /// unknown to the state vector; the owning component must supply exactly one
    /// extra equation (its `n_equations()` grows by one) so the DoF stays
    /// balanced. The actuator value is wired into the component's generic
    /// `CalibIndices::actuator` slot during [`finalize`](Self::finalize).
    pub fn add_free_actuator(&mut self, actuator_id: BoundedVariableId) {
        assert!(
            self.bounded_variables.contains_key(&actuator_id),
            "free actuator '{actuator_id}' must be registered as a bounded variable"
        );
        assert!(
            !self.inverse_control.is_control_linked(&actuator_id),
            "free actuator '{actuator_id}' must not be linked as a hard constraint"
        );
        assert!(
            !self.free_actuators.contains(&actuator_id),
            "free actuator '{actuator_id}' is already registered"
        );
        self.free_actuators.push(actuator_id);
    }

    /// Returns the number of physical free-actuator unknowns.
    pub fn free_actuator_count(&self) -> usize {
        self.free_actuators.len()
    }

    /// Iterates over the registered free-actuator bounded-variable ids in
    /// insertion (state-layout) order.
    pub fn free_actuators(&self) -> impl Iterator<Item = &BoundedVariableId> {
        self.free_actuators.iter()
    }

    /// Returns the number of saturated PI controller loops.
    pub fn saturated_controller_count(&self) -> usize {
        self.saturated_controllers.len()
    }

    /// Iterates over saturated PI controller loops in insertion order.
    pub fn saturated_controllers(&self) -> impl Iterator<Item = &SaturatedController> {
        self.saturated_controllers.iter()
    }

    /// Mutable iterator over saturated PI controllers in insertion order.
    ///
    /// Exposed so an outer driver can retune controllers between solves — in
    /// particular the warm-started **alpha-continuation** that anneals every
    /// override-network selector sharpness toward its target.
    pub fn saturated_controllers_mut(&mut self) -> impl Iterator<Item = &mut SaturatedController> {
        self.saturated_controllers.iter_mut()
    }

    /// Whether any saturated controller uses an override/selector network.
    pub fn has_network_controllers(&self) -> bool {
        self.saturated_controllers.iter().any(|c| c.is_network())
    }

    /// Links a constraint to a bounded control variable for One-Shot inverse control.
    ///
    /// When a constraint is linked to a control variable, the solver adjusts both
    /// the edge states AND the control variable simultaneously to satisfy all
    /// equations (cycle equations + constraint equations).
    ///
    /// # Arguments
    ///
    /// * `constraint_id` - The constraint to link
    /// * `bounded_variable_id` - The control variable to adjust
    ///
    /// # Errors
    ///
    /// Returns `DoFError::ConstraintNotFound` if the constraint doesn't exist.
    /// Returns `DoFError::BoundedVariableNotFound` if the bounded variable doesn't exist.
    /// Returns `DoFError::AlreadyLinked` if the constraint is already linked.
    /// Returns `DoFError::ControlAlreadyLinked` if the control is already linked.
    ///
    /// # Example
    ///
    /// ```rust,ignore
    /// system.link_constraint_to_control(
    ///     &ConstraintId::new("superheat_control"),
    ///     &BoundedVariableId::new("expansion_valve"),
    /// )?;
    /// ```
    pub fn link_constraint_to_control(
        &mut self,
        constraint_id: &ConstraintId,
        bounded_variable_id: &BoundedVariableId,
    ) -> Result<(), DoFError> {
        if !self.constraints.contains_key(constraint_id) {
            return Err(DoFError::ConstraintNotFound {
                constraint_id: constraint_id.clone(),
            });
        }

        if !self.bounded_variables.contains_key(bounded_variable_id) {
            return Err(DoFError::BoundedVariableNotFound {
                bounded_variable_id: bounded_variable_id.clone(),
            });
        }

        self.inverse_control
            .link(constraint_id.clone(), bounded_variable_id.clone())
    }

    /// Unlinks a constraint from its control variable.
    ///
    /// # Arguments
    ///
    /// * `constraint_id` - The constraint to unlink
    ///
    /// # Returns
    ///
    /// The bounded variable ID that was linked, or `None` if not linked.
    pub fn unlink_constraint(&mut self, constraint_id: &ConstraintId) -> Option<BoundedVariableId> {
        self.inverse_control.unlink_constraint(constraint_id)
    }

    /// Unlinks a control variable from its constraint.
    ///
    /// # Arguments
    ///
    /// * `bounded_variable_id` - The control variable to unlink
    ///
    /// # Returns
    ///
    /// The constraint ID that was linked, or `None` if not linked.
    pub fn unlink_control(
        &mut self,
        bounded_variable_id: &BoundedVariableId,
    ) -> Option<ConstraintId> {
        self.inverse_control.unlink_control(bounded_variable_id)
    }

    /// Returns the control variable linked to a constraint.
    pub fn get_control_for_constraint(
        &self,
        constraint_id: &ConstraintId,
    ) -> Option<&BoundedVariableId> {
        self.inverse_control.get_control(constraint_id)
    }

    /// Returns the constraint linked to a control variable.
    pub fn get_constraint_for_control(
        &self,
        bounded_variable_id: &BoundedVariableId,
    ) -> Option<&ConstraintId> {
        self.inverse_control.get_constraint(bounded_variable_id)
    }

    /// Returns the number of constraint-to-control mappings.
    pub fn inverse_control_mapping_count(&self) -> usize {
        self.inverse_control.mapping_count()
    }

    /// Sets the finite difference epsilon for inverse control Jacobian computation.
    ///
    /// # Errors
    ///
    /// Returns `ConstraintError::InvalidEpsilon` if epsilon is not a finite positive value in (0, 1].
    pub fn set_inverse_control_epsilon(&mut self, epsilon: f64) -> Result<(), ConstraintError> {
        self.inverse_control.set_finite_diff_epsilon(epsilon)
    }

    /// Returns the current finite difference epsilon for inverse control.
    pub fn inverse_control_epsilon(&self) -> f64 {
        self.inverse_control.finite_diff_epsilon()
    }

    /// Returns an iterator over linked control variable IDs.
    pub fn linked_controls(&self) -> impl Iterator<Item = &BoundedVariableId> {
        self.inverse_control.linked_controls()
    }

    /// Checks if a constraint is linked to a control variable.
    pub fn is_constraint_linked(&self, constraint_id: &ConstraintId) -> bool {
        self.inverse_control.is_constraint_linked(constraint_id)
    }

    /// Checks if a control variable is linked to a constraint.
    pub fn is_control_linked(&self, bounded_variable_id: &BoundedVariableId) -> bool {
        self.inverse_control.is_control_linked(bounded_variable_id)
    }

    /// Validates degrees of freedom for inverse control.
    ///
    /// For a well-posed system with inverse control:
    ///
    /// $$n_{equations} = n_{edge\_eqs} + n_{constraints}$$
    /// $$n_{unknowns} = n_{edge\_unknowns} + n_{controls}$$
    ///
    /// The system is balanced when: $n_{equations} = n_{unknowns}$
    ///
    /// # Returns
    ///
    /// `Ok(())` if the system is well-posed (balanced or under-constrained).
    ///
    /// # Errors
    ///
    /// Returns `DoFError::OverConstrainedSystem` if there are more equations than unknowns.
    /// Returns `DoFError::UnderConstrainedSystem` if there are fewer equations than unknowns (warning only).
    pub fn validate_inverse_control_dof(&self) -> Result<(), DoFError> {
        let n_edge_unknowns = self.total_state_len;
        let n_controls = self.inverse_control.mapping_count();
        let n_constraints = self.constraints.len();
        let n_unknowns = n_edge_unknowns + n_controls;

        let n_edge_eqs: usize = self
            .graph
            .node_indices()
            .map(|node| {
                self.graph
                    .node_weight(node)
                    .map(|c| c.n_equations())
                    .unwrap_or(0)
            })
            .sum();
        let n_equations = n_edge_eqs + n_constraints;

        if n_equations > n_unknowns {
            Err(DoFError::OverConstrainedSystem {
                constraint_count: n_constraints,
                control_count: n_controls,
                equation_count: n_equations,
                unknown_count: n_unknowns,
            })
        } else if n_equations < n_unknowns {
            Err(DoFError::UnderConstrainedSystem {
                constraint_count: n_constraints,
                control_count: n_controls,
                equation_count: n_equations,
                unknown_count: n_unknowns,
            })
        } else {
            Ok(())
        }
    }

    /// Total residual equations assembled by [`Self::compute_residuals`].
    ///
    /// Includes every component block, inverse constraints, thermal couplings,
    /// and saturated-controller residual pairs.
    pub fn total_equation_count(&self) -> usize {
        let n_comp: usize = self
            .graph
            .node_indices()
            .map(|node| {
                self.graph
                    .node_weight(node)
                    .map(|c| c.n_equations())
                    .unwrap_or(0)
            })
            .sum();
        n_comp
            + self.constraints.len()
            + self.coupling_residual_count()
            + 2 * self.saturated_controllers.len()
    }

    /// Builds a full DoF ledger for the finalized system.
    ///
    /// Counts **all** Newton unknowns (`full_state_vector_len`) against **all**
    /// residual equations (`total_equation_count`), including free actuators and
    /// thermal couplings that the legacy inverse-only check omitted.
    ///
    /// # Panics
    ///
    /// Panics if `finalize()` has not been called (state length is undefined).
    pub fn dof_report(&self) -> DofReport {
        assert!(self.finalized, "call finalize() before dof_report()");

        let mut diagnostics = Vec::new();
        let mut components = Vec::new();
        let mut n_comp_eqs = 0usize;

        // Reverse name map for diagnostics.
        let mut names_by_node: HashMap<NodeIndex, String> = HashMap::new();
        for (name, idx) in &self.component_names {
            names_by_node.insert(*idx, name.clone());
        }

        for node in self.graph.node_indices() {
            let Some(component) = self.graph.node_weight(node) else {
                continue;
            };
            let n = component.n_equations();
            n_comp_eqs += n;
            let name = names_by_node
                .get(&node)
                .cloned()
                .unwrap_or_else(|| format!("node_{}", node.index()));
            let roles = align_roles(&name, n, component.equation_roles(), &mut diagnostics);

            // Heuristic: outlet closures that *consume* a DoF without freeing
            // another (e.g. quality_control, Anchor specs) need an actuator.
            // Emergent-cycle closures (subcooling / superheat / saturated_vapor)
            // replace a fixed P elsewhere — they are the closed-cycle design, not
            // an over-constraint.
            let suspicious_closure = roles.iter().any(|r| match r {
                EquationRole::OutletClosure { kind } => !matches!(
                    kind.as_ref(),
                    "subcooling" | "superheat" | "saturated_vapor"
                ),
                _ => false,
            });
            if suspicious_closure
                && self.free_actuators.is_empty()
                && self.saturated_controllers.is_empty()
                && self.inverse_control.mapping_count() == 0
            {
                diagnostics.push(format!(
                    "component `{name}` declares OutletClosure residual(s) but the system \
                     has no free actuator / saturated controller / inverse control — \
                     fixing this quantity without freeing another unknown over-constrains \
                     a closed cycle (pair with EXV opening, speed, level, or drop the residual)"
                ));
            }

            components.push(ComponentEquationBlock {
                component_name: name,
                node_index: node.index(),
                n_equations: n,
                roles,
            });
        }

        let mut system_equations = Vec::new();
        for id in self.constraints.keys() {
            system_equations.push(EquationRole::ControlTracking {
                name: id.to_string(),
            });
        }
        for _ in 0..self.coupling_residual_count() {
            system_equations.push(EquationRole::CouplingDuty);
        }
        for _ in 0..self.saturated_controllers.len() {
            system_equations.push(EquationRole::SaturatedControl {
                role: "actuator_law",
            });
            system_equations.push(EquationRole::SaturatedControl { role: "tracking" });
        }

        let n_equations = n_comp_eqs + system_equations.len();
        let n_unknowns = self.full_state_vector_len();
        let balance = SystemDofBalance::from_counts(n_equations, n_unknowns);

        // Compact unknown catalog (edge layout + extras).
        // Physical block is |B| + 2|E| (+ internal); catalog via edge_to_state.
        let mut branch_ms = std::collections::BTreeSet::new();
        let mut edges: Vec<_> = self.edge_to_state.iter().collect();
        edges.sort_by_key(|(e, _)| e.index());
        for (_edge, (m, _p, _h)) in &edges {
            branch_ms.insert(*m);
        }
        let mut unknowns = Vec::with_capacity(n_unknowns);
        for branch_id in 0..branch_ms.len() {
            unknowns.push(UnknownKind::BranchMassFlow { branch_id });
        }
        for edge_ordinal in 0..edges.len() {
            unknowns.push(UnknownKind::EdgePressure { edge_ordinal });
            unknowns.push(UnknownKind::EdgeEnthalpy { edge_ordinal });
        }
        let physical_accounted = unknowns.len();
        if physical_accounted < self.total_state_len {
            for index in physical_accounted..self.total_state_len {
                unknowns.push(UnknownKind::Internal { index });
            }
        }
        for id in self.inverse_control.linked_controls() {
            unknowns.push(UnknownKind::InverseControl { id: id.to_string() });
        }
        for index in 0..self.coupling_residual_count() {
            unknowns.push(UnknownKind::CouplingHeat { index });
        }
        for index in 0..self.saturated_controllers.len() {
            unknowns.push(UnknownKind::SaturatedActuator { index });
            unknowns.push(UnknownKind::SaturatedIntegrator { index });
        }
        for id in &self.free_actuators {
            unknowns.push(UnknownKind::FreeActuator { id: id.to_string() });
        }

        if unknowns.len() != n_unknowns {
            diagnostics.push(format!(
                "unknown catalog length {} differs from full_state_vector_len {n_unknowns} \
                 (catalog is diagnostic only; counts use full_state_vector_len)",
                unknowns.len()
            ));
        }

        if n_equations != self.total_equation_count() {
            diagnostics.push(format!(
                "internal inconsistency: ledger n_equations {n_equations} != total_equation_count {}",
                self.total_equation_count()
            ));
        }

        DofReport {
            n_equations,
            n_unknowns,
            balance,
            components,
            system_equations,
            unknowns,
            diagnostics,
        }
    }

    /// Hard DoF gate: returns `Ok(())` only when `n_equations == n_unknowns`.
    ///
    /// Call after [`finalize`](Self::finalize). Prefer relying on the automatic
    /// check inside `finalize` for system construction.
    pub fn validate_system_dof(&self) -> Result<(), SystemDofError> {
        let report = self.dof_report();
        if report.balance.is_balanced() {
            Ok(())
        } else {
            Err(SystemDofError::Imbalance {
                n_equations: report.n_equations,
                n_unknowns: report.n_unknowns,
                balance: report.balance,
                summary: report.summary(),
            })
        }
    }

    /// Returns the state vector index for a control variable.
    ///
    /// Control variables are appended after edge states in the state vector:
    ///
    /// ```text
    /// State Vector = [Edge States | Control Variables | Thermal Coupling Temps]
    ///                [m0, P0, h0, m1, P1, h1, ... | ctrl0, ctrl1, ... | T_hot0, T_cold0, ...]
    /// ```
    ///
    /// The index for control variable `i` is: `total_state_len + i`
    ///
    /// # Arguments
    ///
    /// * `id` - The bounded variable identifier
    ///
    /// # Returns
    ///
    /// The state vector index, or `None` if the variable is not linked.
    pub fn control_variable_state_index(&self, id: &BoundedVariableId) -> Option<usize> {
        let base = self.total_state_len;
        for (index, linked_id) in self.inverse_control.linked_controls().enumerate() {
            if linked_id == id {
                return Some(base + index);
            }
        }
        for (index, ctrl) in self.saturated_controllers.iter().enumerate() {
            if ctrl.actuator() == id {
                return Some(self.saturated_u_index(index));
            }
        }
        None
    }

    /// Returns the bounded variable for a given state index.
    pub fn get_bounded_variable_by_state_index(
        &self,
        state_index: usize,
    ) -> Option<&BoundedVariable> {
        let base = self.total_state_len;
        if state_index < base {
            return None;
        }
        let control_idx = state_index - base;
        self.inverse_control
            .linked_controls()
            .nth(control_idx)
            .and_then(|id| self.bounded_variables.get(id))
            .or_else(|| {
                let saturated_base = self.saturated_base_index();
                if state_index < saturated_base {
                    return None;
                }
                let saturated_offset = state_index - saturated_base;
                if saturated_offset % 2 != 0 {
                    return None;
                }
                self.saturated_controllers
                    .get(saturated_offset / 2)
                    .and_then(|ctrl| self.bounded_variables.get(ctrl.actuator()))
            })
    }

    /// Returns the bounds (min, max) for a given state index if it corresponds to a bounded control variable.
    pub fn get_bounds_for_state_index(&self, state_index: usize) -> Option<(f64, f64)> {
        self.get_bounded_variable_by_state_index(state_index)
            .map(|var| (var.min(), var.max()))
    }

    /// Returns physical solver bounds for state variables.
    ///
    /// Pressure edge states are strictly positive thermodynamic variables; Newton
    /// steps are clipped before component residuals are evaluated instead of
    /// letting components invent fallback properties for non-physical pressures.
    pub fn get_solver_bounds_for_state_index(&self, state_index: usize) -> Option<(f64, f64)> {
        if let Some(bounds) = self.get_bounds_for_state_index(state_index) {
            return Some(bounds);
        }

        self.edge_to_state
            .values()
            .any(|(_, p_idx, _)| *p_idx == state_index)
            .then_some((MIN_SOLVER_PRESSURE_PA, f64::MAX))
    }

    /// Returns the total state vector length including control variables.
    ///
    /// ```text
    /// full_state_len = total_state_len
    ///                + hard_control_count
    ///                + thermal_coupling_count
    ///                + 2 * saturated_controller_count
    /// ```
    pub fn full_state_vector_len(&self) -> usize {
        self.total_state_len
            + self.inverse_control.mapping_count()
            + self.coupling_residual_count()
            + 2 * self.saturated_controllers.len()
            + self.free_actuators.len()
    }

    /// Base state index of the physical free-actuator block (arch-6).
    /// Free actuators are laid out immediately after the saturated block.
    fn free_actuator_base_index(&self) -> usize {
        self.total_state_len
            + self.inverse_control.mapping_count()
            + self.coupling_residual_count()
            + 2 * self.saturated_controllers.len()
    }

    /// State index of the i-th physical free actuator (arch-6).
    pub fn free_actuator_index(&self, i: usize) -> usize {
        self.free_actuator_base_index() + i
    }

    /// State index of the i-th thermal coupling's heat unknown Q [W].
    ///
    /// Coupling unknowns live after the hard inverse-control block:
    /// `[Edge States | Hard Controls | Coupling Q | Saturated (u,x) | Free Actuators]`.
    pub fn coupling_state_index(&self, i: usize) -> usize {
        self.total_state_len + self.inverse_control.mapping_count() + i
    }

    pub fn saturated_base_index(&self) -> usize {
        self.total_state_len + self.inverse_control.mapping_count() + self.coupling_residual_count()
    }

    pub fn saturated_u_index(&self, i: usize) -> usize {
        self.saturated_base_index() + 2 * i
    }

    fn saturated_x_index(&self, i: usize) -> usize {
        self.saturated_base_index() + 2 * i + 1
    }

    /// Returns an ordered list of linked control variable IDs with their state indices.
    ///
    /// Useful for constructing the state vector or Jacobian columns.
    pub fn control_variable_indices(&self) -> Vec<(&BoundedVariableId, usize)> {
        let base = self.total_state_len;
        self.inverse_control
            .linked_controls()
            .enumerate()
            .map(|(i, id)| (id, base + i))
            .collect()
    }

    /// Returns the number of coupling residual equations (one per thermal coupling).
    ///
    /// The solver must reserve this many rows in the residual vector for coupling
    /// heat balance equations. Use [`coupling_residuals`](Self::coupling_residuals) to fill them.
    pub fn coupling_residual_count(&self) -> usize {
        self.thermal_couplings.len()
    }

    /// Fills coupling residuals into `out`.
    ///
    /// For each thermal coupling, the residual is the heat transfer rate Q (W) into the cold
    /// circuit: Q = η·UA·(T_hot − T_cold). The solver typically uses this in a heat balance
    /// (e.g. r = Q_actual − Q_expected). Temperatures must be in Kelvin.
    ///
    /// # Arguments
    ///
    /// * `temperatures` - One (T_hot_K, T_cold_K) per coupling; length must equal
    ///   `thermal_coupling_count()`. The solver obtains these from state (e.g. P, h → T via fluid backend).
    /// * `out` - Slice to write residuals; length must be at least `coupling_residual_count()`.
    pub fn coupling_residuals(&self, temperatures: &[(f64, f64)], out: &mut [f64]) {
        assert!(
            temperatures.len() == self.thermal_couplings.len(),
            "temperatures.len() must equal thermal_coupling_count()"
        );
        assert!(
            out.len() >= self.thermal_couplings.len(),
            "out.len() must be at least coupling_residual_count()"
        );
        for (i, coupling) in self.thermal_couplings.iter().enumerate() {
            let (t_hot_k, t_cold_k) = temperatures[i];
            let t_hot = Temperature::from_kelvin(t_hot_k);
            let t_cold = Temperature::from_kelvin(t_cold_k);
            out[i] = crate::coupling::compute_coupling_heat(coupling, t_hot, t_cold);
        }
    }

    /// Returns Jacobian entries for coupling residuals with respect to temperature state.
    ///
    /// The solver state may include temperature unknowns for coupling interfaces (or T derived from P, h).
    /// For each coupling i: ∂Q_i/∂T_hot = η·UA, ∂Q_i/∂T_cold = −η·UA.
    ///
    /// # Arguments
    ///
    /// * `row_offset` - First row index for coupling equations in the global residual vector.
    /// * `t_hot_cols` - State column index for T_hot per coupling; length = `thermal_coupling_count()`.
    /// * `t_cold_cols` - State column index for T_cold per coupling; length = `thermal_coupling_count()`.
    ///
    /// # Returns
    ///
    /// `(row, col, value)` tuples for the Jacobian. Row is `row_offset + coupling_index`.
    pub fn coupling_jacobian_entries(
        &self,
        row_offset: usize,
        t_hot_cols: &[usize],
        t_cold_cols: &[usize],
    ) -> Vec<(usize, usize, f64)> {
        assert!(
            t_hot_cols.len() == self.thermal_couplings.len()
                && t_cold_cols.len() == self.thermal_couplings.len(),
            "t_hot_cols and t_cold_cols length must equal thermal_coupling_count()"
        );
        let mut entries = Vec::with_capacity(2 * self.thermal_couplings.len());
        for (i, coupling) in self.thermal_couplings.iter().enumerate() {
            let dr_dt_hot = coupling.efficiency * coupling.ua.to_watts_per_kelvin();
            let dr_dt_cold = -dr_dt_hot;
            let row = row_offset + i;
            entries.push((row, t_hot_cols[i], dr_dt_hot));
            entries.push((row, t_cold_cols[i], dr_dt_cold));
        }
        entries
    }

    /// Returns true if the graph contains a cycle.
    ///
    /// Refrigeration circuits form cycles (compressor → condenser → valve → evaporator → compressor),
    /// so cycles are expected and valid.
    pub fn is_cyclic(&self) -> bool {
        algo::is_cyclic_directed(&self.graph)
    }

    /// Iterates over (component, edge_indices) for Jacobian assembly.
    ///
    /// For each node, yields the component and a map from edge index to (state_index_p, state_index_h)
    /// for edges incident to that node (incoming and outgoing).
    ///
    /// # Panics
    ///
    /// Panics if `finalize()` has not been called.
    pub fn traverse_for_jacobian(
        &self,
    ) -> impl Iterator<Item = (NodeIndex, &dyn Component, Vec<(EdgeIndex, usize, usize)>)> {
        assert!(
            self.finalized,
            "call finalize() before traverse_for_jacobian()"
        );

        self.graph.node_indices().map(move |node_idx| {
            let component = self.graph.node_weight(node_idx).unwrap();
            let mut edge_indices = Vec::new();

            for edge_ref in self
                .graph
                .edges_directed(node_idx, petgraph::Direction::Incoming)
            {
                let edge_idx = edge_ref.id();
                if let Some(&(_m, p, h)) = self.edge_to_state.get(&edge_idx) {
                    edge_indices.push((edge_idx, p, h));
                }
            }
            for edge_ref in self
                .graph
                .edges_directed(node_idx, petgraph::Direction::Outgoing)
            {
                let edge_idx = edge_ref.id();
                if let Some(&(_m, p, h)) = self.edge_to_state.get(&edge_idx) {
                    edge_indices.push((edge_idx, p, h));
                }
            }

            (node_idx, component.as_ref(), edge_indices)
        })
    }

    /// Assembles residuals from all components.
    ///
    /// Components receive the full state slice and write to their equation indices.
    /// Equation indices are computed from component order and `n_equations()`.
    ///
    /// # Errors
    ///
    /// Returns `ComponentError::InvalidResidualDimensions` if `residuals.len()` is
    /// less than the total number of equations across all components.
    pub fn compute_residuals(
        &self,
        state: &StateSlice,
        residuals: &mut ResidualVector,
    ) -> Result<(), ComponentError> {
        let mut total_eqs: usize = self
            .traverse_for_jacobian()
            .map(|(_, c, _)| c.n_equations())
            .sum();
        total_eqs += self.constraints.len()
            + self.coupling_residual_count()
            + 2 * self.saturated_controllers.len();

        if residuals.len() < total_eqs {
            return Err(ComponentError::InvalidResidualDimensions {
                expected: total_eqs,
                actual: residuals.len(),
            });
        }

        let mut eq_offset = 0;
        for (_node_idx, component, _edge_indices) in self.traverse_for_jacobian() {
            let n = component.n_equations();
            if n > 0 {
                let mut temp = vec![0.0; n];
                component.compute_residuals(state, &mut temp)?;
                residuals[eq_offset..eq_offset + n].copy_from_slice(&temp);
            }
            eq_offset += n;
        }

        // Add constraints
        let control_values: Vec<f64> = self
            .control_variable_indices()
            .into_iter()
            .map(|(_, idx)| state[idx])
            .collect();
        let measured = self.extract_constraint_values_with_controls(state, &control_values);
        let n_constraints = self
            .compute_constraint_residuals(state, &mut residuals[eq_offset..], &measured)
            .map_err(|e| ComponentError::CalculationFailed(e.to_string()))?;
        eq_offset += n_constraints;

        // Add couplings. Each coupling owns one unknown Q = state[coupling_state_index(i)].
        //
        // * Physical mode (`hot_component`/`cold_component` set): the residual
        //   closes Q against the hot component's *measured* duty:
        //   r_i = Q − η·duty_hot(state). The cold-side `ThermalLoad` consumes Q
        //   in its own energy balance (wired in `finalize()`), so the heat
        //   genuinely crosses the circuit boundary and the First Law closes.
        // * Legacy mode (no components referenced): the residual pins the orphan
        //   unknown, r_i = Q (⇒ Q = 0). This replaces the former 300 K stub,
        //   whose all-zero row/column made the Jacobian structurally singular.
        let n_couplings = self.coupling_residual_count();
        for (i, coupling) in self.thermal_couplings.iter().enumerate() {
            let q_idx = self.coupling_state_index(i);
            let q = state.get(q_idx).copied().unwrap_or(0.0);
            let target = if coupling.is_physical() {
                let duty = coupling
                    .hot_component
                    .as_deref()
                    .and_then(|name| self.component_names.get(name))
                    .and_then(|&node| self.graph.node_weight(node))
                    .and_then(|c| {
                        c.measure_output(entropyk_components::MeasuredOutput::Capacity, state)
                    })
                    .unwrap_or(0.0);
                coupling.efficiency * coupling.duty_scale * duty
            } else {
                0.0
            };
            residuals[eq_offset + i] = q - target;
        }
        eq_offset += n_couplings;

        for (i, ctrl) in self.saturated_controllers.iter().enumerate() {
            let u_idx = self.saturated_u_index(i);
            let x_idx = self.saturated_x_index(i);
            if state.len() <= x_idx {
                return Err(ComponentError::InvalidStateDimensions {
                    expected: x_idx + 1,
                    actual: state.len(),
                });
            }
            let u = state[u_idx];
            let x = state[x_idx];

            residuals[eq_offset] = ctrl.residual_u(u, x);
            eq_offset += 1;
            if ctrl.is_network() {
                let measured: Vec<f64> = ctrl
                    .objectives()
                    .iter()
                    .map(|o| self.measure_saturated_output(&o.output, state))
                    .collect();
                let e = ctrl.network_error(&measured);
                residuals[eq_offset] = ctrl.residual_y_network(e, x);
            } else {
                let y = self.measure_saturated_output(ctrl.output(), state);
                residuals[eq_offset] = ctrl.residual_y(y, x);
            }
            eq_offset += 1;
        }

        Ok(())
    }

    /// Assembles Jacobian entries from all components.
    ///
    /// Each component receives the state and writes to JacobianBuilder. The
    /// [`traverse_for_jacobian`](Self::traverse_for_jacobian) iterator provides
    /// the edge→state mapping `(EdgeIndex, state_index_p, state_index_h)` per
    /// component. Components must know their port→state mapping (e.g. from graph
    /// construction in Story 3.2) to write correct column indices.
    pub fn assemble_jacobian(
        &self,
        state: &StateSlice,
        jacobian: &mut JacobianBuilder,
    ) -> Result<(), ComponentError> {
        let mut row_offset = 0;
        for (_node_idx, component, _edge_indices) in self.traverse_for_jacobian() {
            let n = component.n_equations();
            if n > 0 {
                // Components write rows 0..n-1; we offset to global equation indices.
                let mut temp_builder = JacobianBuilder::new();
                component.jacobian_entries(state, &mut temp_builder)?;
                for (r, c, v) in temp_builder.entries() {
                    jacobian.add_entry(row_offset + r, *c, *v);
                }
            }
            row_offset += n;
        }

        // Add constraints jacobian
        let control_values: Vec<f64> = self
            .control_variable_indices()
            .into_iter()
            .map(|(_, idx)| state[idx])
            .collect();
        let constraint_jac =
            self.compute_inverse_control_jacobian(state, row_offset, &control_values);
        for (r, c, v) in constraint_jac {
            jacobian.add_entry(r, c, v);
        }

        // Thermal coupling rows: r_i = Q_i − η·duty_hot(state).
        // ∂r/∂Q = 1 always (also fixes the legacy stub's singular row); in
        // physical mode the plant coupling ∂r/∂col = −η·∂duty/∂col is formed by
        // central finite differences over the hot component's incident states
        // (same pattern as the saturated-controller measurement Jacobian).
        let coupling_row_offset = row_offset + self.constraints.len();
        if !self.thermal_couplings.is_empty() {
            let eps = self.inverse_control.finite_diff_epsilon();
            let mut state_mut = state.to_vec();
            for (i, coupling) in self.thermal_couplings.iter().enumerate() {
                let row = coupling_row_offset + i;
                jacobian.add_entry(row, self.coupling_state_index(i), 1.0);
                if !coupling.is_physical() {
                    continue;
                }
                let Some(&node_idx) = coupling
                    .hot_component
                    .as_deref()
                    .and_then(|name| self.component_names.get(name))
                else {
                    continue;
                };
                let measure = |st: &[f64]| -> f64 {
                    self.graph
                        .node_weight(node_idx)
                        .and_then(|c| {
                            c.measure_output(entropyk_components::MeasuredOutput::Capacity, st)
                        })
                        .unwrap_or(0.0)
                };
                for col in self.incident_state_indices_for_component(node_idx) {
                    let orig = state_mut[col];
                    state_mut[col] = orig + eps;
                    let d_plus = measure(&state_mut);
                    state_mut[col] = orig - eps;
                    let d_minus = measure(&state_mut);
                    state_mut[col] = orig;
                    let dduty_dcol = (d_plus - d_minus) / (2.0 * eps);
                    let derivative = -coupling.efficiency * coupling.duty_scale * dduty_dcol;
                    if derivative.abs() > 1e-10 {
                        jacobian.add_entry(row, col, derivative);
                    }
                }
            }
        }

        let saturated_row_offset =
            row_offset + self.constraints.len() + self.coupling_residual_count();
        let saturated_jac = self.compute_saturated_control_jacobian(state, saturated_row_offset);
        for (r, c, v) in saturated_jac {
            jacobian.add_entry(r, c, v);
        }
        Ok(())
    }

    /// Tolerance for mass balance validation [kg/s].
    ///
    /// This value (1e-9 kg/s) is tight enough to catch numerical issues
    /// while allowing for floating-point rounding errors.
    pub const MASS_BALANCE_TOLERANCE_KG_S: f64 = 1e-9;

    /// Tolerance for energy balance validation in Watts (1e-6 kW)
    pub const ENERGY_BALANCE_TOLERANCE_W: f64 = 1e-3;

    /// Verifies that global mass balance is conserved.
    ///
    /// Sums the mass flow rates at the ports of each component and ensures they
    /// sum to zero within a tight tolerance (1e-9 kg/s).
    ///
    /// # Returns
    ///
    /// * `Ok(())` if all components pass mass balance validation
    /// * `Err(SolverError::Validation)` if any component violates mass conservation
    ///
    /// # Note
    ///
    /// Components without `port_mass_flows` implementation are logged as warnings
    /// and skipped. This ensures visibility of incomplete implementations without
    /// failing the validation.
    pub fn check_mass_balance(&self, state: &StateSlice) -> Result<(), crate::SolverError> {
        let mut total_mass_error = 0.0;
        let mut has_violation = false;
        let mut components_checked = 0usize;
        let mut components_skipped = 0usize;

        for (node_idx, component, _edge_indices) in self.traverse_for_jacobian() {
            match component.port_mass_flows(state) {
                Ok(mass_flows) => {
                    let sum: f64 = mass_flows.iter().map(|m| m.to_kg_per_s()).sum();
                    if sum.abs() > Self::MASS_BALANCE_TOLERANCE_KG_S {
                        has_violation = true;
                        total_mass_error += sum.abs();
                        tracing::warn!(
                            node_index = node_idx.index(),
                            mass_imbalance_kg_s = sum,
                            "Mass balance violation detected at component"
                        );
                    }
                    components_checked += 1;
                }
                Err(e) => {
                    components_skipped += 1;
                    tracing::warn!(
                        node_index = node_idx.index(),
                        error = %e,
                        "Component does not implement port_mass_flows - skipping mass balance check"
                    );
                }
            }
        }

        tracing::debug!(
            components_checked,
            components_skipped,
            total_mass_error_kg_s = total_mass_error,
            "Mass balance validation complete"
        );

        if has_violation {
            return Err(crate::SolverError::Validation {
                mass_error: total_mass_error,
                energy_error: 0.0,
            });
        }
        Ok(())
    }

    /// Verifies the First Law of Thermodynamics for all components in the system.
    ///
    /// Validates that ΣQ - ΣW + Σ(ṁ·h) = 0 for each component.
    /// Returns `SolverError::Validation` if any component violates the balance.
    pub fn check_energy_balance(&self, state: &StateSlice) -> Result<(), crate::SolverError> {
        let mut total_energy_error = 0.0;
        let mut has_violation = false;
        let mut components_checked = 0usize;
        let mut components_skipped = 0usize;
        let mut skipped_components: Vec<String> = Vec::new();

        for (node_idx, component, _edge_indices) in self.traverse_for_jacobian() {
            let energy_transfers = component.energy_transfers(state);
            let mass_flows = component.port_mass_flows(state);
            let enthalpies = component.port_enthalpies(state);

            match (energy_transfers, mass_flows, enthalpies) {
                (Some((heat, work)), Ok(m_flows), Ok(h_flows))
                    if m_flows.len() == h_flows.len() =>
                {
                    let mut net_energy_flow = 0.0;
                    for (m, h) in m_flows.iter().zip(h_flows.iter()) {
                        net_energy_flow += m.to_kg_per_s() * h.to_joules_per_kg();
                    }

                    let balance = heat.to_watts() - work.to_watts() + net_energy_flow;

                    if balance.abs() > Self::ENERGY_BALANCE_TOLERANCE_W {
                        has_violation = true;
                        total_energy_error += balance.abs();
                        tracing::warn!(
                            node_index = node_idx.index(),
                            energy_imbalance_w = balance,
                            "Energy balance violation detected at component"
                        );
                    }
                    components_checked += 1;
                }
                _ => {
                    components_skipped += 1;
                    let component_type = std::any::type_name_of_val(component)
                        .split("::")
                        .last()
                        .unwrap_or("unknown");
                    let component_info =
                        format!("{} (type: {})", component.signature(), component_type);
                    skipped_components.push(component_info.clone());

                    tracing::warn!(
                        component = %component_info,
                        node_index = node_idx.index(),
                        "Component lacks energy_transfers() or port_enthalpies() - SKIPPED in energy balance validation"
                    );
                }
            }
        }

        // Summary warning if components were skipped
        if components_skipped > 0 {
            tracing::warn!(
                components_checked = components_checked,
                components_skipped = components_skipped,
                skipped = ?skipped_components,
                "Energy balance validation incomplete: {} component(s) skipped. \
                 Implement energy_transfers() and port_enthalpies() for full validation.",
                components_skipped
            );
        } else {
            tracing::debug!(
                components_checked,
                components_skipped,
                total_energy_error_w = total_energy_error,
                "Energy balance validation complete"
            );
        }

        if has_violation {
            return Err(crate::SolverError::Validation {
                mass_error: 0.0,
                energy_error: total_energy_error,
            });
        }
        Ok(())
    }

    /// Computes the aggregate cycle performance (cooling/heating duty, work input,
    /// COP/EER) from the solved state vector.
    ///
    /// Each component's [`Component::energy_transfers`] is evaluated at `state`
    /// and classified by the physical sign convention (`Q > 0` = heat into the
    /// refrigerant = cooling; `Q < 0` = heat rejected = heating; `W < 0` = work
    /// done on the fluid = power input). Components that do not report energy
    /// transfers are skipped. Returns `None` if no component reported anything.
    ///
    /// This is post-processing on a converged state — it makes no solver calls
    /// and does not mutate the system.
    pub fn cycle_performance(&self, state: &StateSlice) -> Option<CyclePerformance> {
        let mut q_cooling_w = 0.0;
        let mut q_heating_w = 0.0;
        let mut work_input_w = 0.0;
        let mut components_counted = 0usize;

        for (_node_idx, component, _edges) in self.traverse_for_jacobian() {
            // Skip secondary-side receivers (e.g. `ThermalLoad`): the heat they
            // absorb is the primary cycle's rejected duty, not extra cooling.
            if !component.counts_in_cycle_performance() {
                continue;
            }
            if let Some((heat, work)) = component.energy_transfers(state) {
                let q = heat.to_watts();
                let w = work.to_watts();
                if !(q.is_finite() && w.is_finite()) {
                    continue;
                }
                if q > 0.0 {
                    q_cooling_w += q;
                } else if q < 0.0 {
                    q_heating_w += -q;
                }
                // Work done ON the fluid is reported negative; accumulate its
                // magnitude as the electrical/shaft power input.
                if w < 0.0 {
                    work_input_w += -w;
                }
                components_counted += 1;
            }
        }

        if components_counted == 0 {
            return None;
        }
        Some(CyclePerformance {
            q_cooling_w,
            q_heating_w,
            work_input_w,
            components_counted,
        })
    }

    /// Generates a deterministic byte representation of the system configuration.
    /// Used for simulation traceability logic.
    pub fn generate_canonical_bytes(&self) -> Vec<u8> {
        let mut repr = String::new();
        repr.push_str("Nodes:\n");
        // To be deterministic, we just iterate in graph order which is stable
        // as long as we don't delete nodes.
        for node in self.graph.node_indices() {
            let circuit_id = self.node_to_circuit.get(&node).map(|c| c.0).unwrap_or(0);
            repr.push_str(&format!(
                "  Node({}): Circuit({})\n",
                node.index(),
                circuit_id
            ));
            if let Some(comp) = self.graph.node_weight(node) {
                repr.push_str(&format!("    Signature: {}\n", comp.signature()));
            }
        }
        repr.push_str("Edges:\n");
        for edge_idx in self.graph.edge_indices() {
            if let Some((src, tgt)) = self.graph.edge_endpoints(edge_idx) {
                repr.push_str(&format!("  Edge: {} -> {}\n", src.index(), tgt.index()));
            }
        }
        repr.push_str("Thermal Couplings:\n");
        let mut couplings: Vec<_> = self.thermal_couplings.iter().collect();
        couplings.sort_by(|a, b| {
            (a.hot_circuit.0, a.cold_circuit.0)
                .cmp(&(b.hot_circuit.0, b.cold_circuit.0))
                .then(
                    a.ua.0
                        .partial_cmp(&b.ua.0)
                        .unwrap_or(std::cmp::Ordering::Equal),
                )
                .then(
                    a.efficiency
                        .partial_cmp(&b.efficiency)
                        .unwrap_or(std::cmp::Ordering::Equal),
                )
        });
        for coupling in couplings {
            repr.push_str(&format!(
                "  Hot: {}, Cold: {}, UA: {}\n",
                coupling.hot_circuit.0, coupling.cold_circuit.0, coupling.ua
            ));
        }
        repr.push_str("Constraints:\n");
        let mut constraint_keys: Vec<_> = self.constraints.keys().collect();
        constraint_keys.sort_by_key(|k| k.as_str());
        for key in constraint_keys {
            let c = &self.constraints[key];
            repr.push_str(&format!("  {}: {}\n", c.id().as_str(), c.target_value()));
        }
        repr.push_str("Bounded Variables:\n");
        let mut bounded_keys: Vec<_> = self.bounded_variables.keys().collect();
        bounded_keys.sort_by_key(|k| k.as_str());
        for key in bounded_keys {
            let var = &self.bounded_variables[key];
            repr.push_str(&format!(
                "  {}: [{}, {}]\n",
                var.id().as_str(),
                var.min(),
                var.max()
            ));
        }

        repr.push_str("Inverse Control Mappings:\n");
        // For inverse control mappings, they are ordered internally. We'll just iterate linked controls.
        for (i, (constraint, bounded_var)) in self.inverse_control.mappings().enumerate() {
            repr.push_str(&format!(
                "  Mapping {}: {} -> {}\n",
                i,
                constraint.as_str(),
                bounded_var.as_str()
            ));
        }

        repr.into_bytes()
    }

    /// Computes the SHA-256 hash uniquely identifying the input configuration.
    pub fn input_hash(&self) -> String {
        use sha2::{Digest, Sha256};
        let mut hasher = Sha256::new();
        hasher.update(self.generate_canonical_bytes());
        format!("{:064x}", hasher.finalize())
    }

    // ========== JSON Serialization API ==========

    /// Serializes the system to a JSON string.
    ///
    /// This method captures the complete system state including topology,
    /// component parameters, and metadata in a human-readable JSON format.
    ///
    /// # Errors
    ///
    /// Returns `ThermoError::SerializationError` if JSON serialization fails.
    ///
    /// # Examples
    ///
    /// ```no_run
    /// # use entropyk_solver::System;
    /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
    /// let system = System::new();
    /// let json_string = system.to_json_string()?;
    /// println!("System JSON: {}", json_string);
    /// # Ok(())
    /// # }
    /// ```
    pub fn to_json_string(&self) -> Result<String, crate::error::ThermoError> {
        use crate::snapshot::{
            BoundedVariableSnapshot, ConstraintSnapshot, EdgeSnapshot, FluidBackendInfo,
            SolverConfigSnapshot, SystemSnapshot, TopologySnapshot,
        };
        use std::collections::HashMap;

        tracing::info!("Serializing system to JSON");

        let reverse_names: HashMap<NodeIndex, &String> =
            self.component_names.iter().map(|(n, &i)| (i, n)).collect();

        // Extract topology with port names
        let mut edges = Vec::new();
        for edge in self.graph.edge_indices() {
            let (source, target) = self.graph.edge_endpoints(edge).unwrap();
            let source_node = self.graph.node_weight(source).unwrap();
            let target_node = self.graph.node_weight(target).unwrap();

            // Derive port names from component port_names() or defaults
            let source_ports = source_node.port_names();
            let target_ports = target_node.port_names();

            // Count how many edges connect TO the target (this edge's index at target)
            let target_incoming: Vec<_> = self
                .graph
                .edges_directed(target, petgraph::Direction::Incoming)
                .collect();
            let target_port_idx = target_incoming
                .iter()
                .position(|e| e.id() == edge)
                .unwrap_or(0);

            // Count how many edges leave FROM the source (this edge's index at source)
            let source_outgoing: Vec<_> = self
                .graph
                .edges_directed(source, petgraph::Direction::Outgoing)
                .collect();
            let source_port_idx = source_outgoing
                .iter()
                .position(|e| e.id() == edge)
                .unwrap_or(0);

            let source_port_name = source_ports
                .get(source_port_idx)
                .cloned()
                .unwrap_or_else(|| format!("port_{}", source_port_idx));
            let target_port_name = target_ports
                .get(target_port_idx)
                .cloned()
                .unwrap_or_else(|| format!("port_{}", target_port_idx));

            edges.push(EdgeSnapshot {
                source: reverse_names
                    .get(&source)
                    .map(|s| s.to_string())
                    .unwrap_or_else(|| source_node.signature()),
                source_port: source_port_name,
                target: reverse_names
                    .get(&target)
                    .map(|s| s.to_string())
                    .unwrap_or_else(|| target_node.signature()),
                target_port: target_port_name,
                circuit_id: self.edge_circuit(edge).0,
            });
        }

        // Extract component parameters (use unique key: registered name or signature+index)
        let mut parameters = HashMap::new();
        for node in self.graph.node_indices() {
            if let Some(component) = self.graph.node_weight(node) {
                let params = component.to_params();
                let key = reverse_names
                    .get(&node)
                    .map(|s| (*s).clone())
                    .unwrap_or_else(|| component.signature());
                parameters.insert(key.to_string(), params);
            }
        }

        // Build component_names and circuit_assignments maps
        let component_names: HashMap<String, String> = self
            .component_names
            .iter()
            .map(|(name, &node_idx)| {
                let comp = self.graph.node_weight(node_idx);
                let type_name = comp
                    .map(|c| {
                        let sig = c.to_params().component_type.clone();
                        sig
                    })
                    .unwrap_or_else(|| "Unknown".to_string());
                (name.clone(), type_name)
            })
            .collect();

        let circuit_assignments: HashMap<String, u16> = self
            .component_names
            .iter()
            .map(|(name, &node_idx)| {
                let cid = self
                    .node_to_circuit
                    .get(&node_idx)
                    .map(|c| c.0)
                    .unwrap_or(0);
                (name.clone(), cid)
            })
            .collect();

        // Create snapshot
        let snapshot = SystemSnapshot {
            version: "1.0".to_string(),
            topology: TopologySnapshot {
                edges,
                thermal_couplings: self.thermal_couplings.clone(),
            },
            parameters,
            fluid_state: {
                let mut data = Vec::with_capacity(self.graph.edge_count() * 2);
                for edge in self.graph.edge_indices() {
                    let (source, _target) = self.graph.edge_endpoints(edge).unwrap();
                    let component = self.graph.node_weight(source).unwrap();
                    let ports = component.get_ports();
                    let outgoing: Vec<_> = self
                        .graph
                        .edges_directed(source, petgraph::Direction::Outgoing)
                        .collect();
                    let port_idx = outgoing.iter().position(|e| e.id() == edge).unwrap_or(0);
                    if let Some(port) = ports.get(port_idx) {
                        data.push(port.pressure().to_pascals());
                        data.push(port.enthalpy().to_joules_per_kg());
                    } else {
                        data.push(0.0);
                        data.push(0.0);
                    }
                }
                if data.is_empty() {
                    None
                } else {
                    entropyk_core::SystemState::try_from(data).ok()
                }
            },
            fluid_backend: FluidBackendInfo {
                name: "CoolPropBackend".to_string(),
                version: env!("CARGO_PKG_VERSION").to_string(),
                hash: None,
            },
            solver_config: Some(SolverConfigSnapshot::default()),
            component_names,
            circuit_assignments,
            constraints: self
                .constraints
                .iter()
                .map(|(id, c)| ConstraintSnapshot {
                    id: id.as_str().to_string(),
                    component: c.output().component_id().to_string(),
                    output_type: c.output().constraint_type_name().to_string(),
                    target: c.target_value(),
                })
                .collect(),
            bounded_variables: self
                .bounded_variables
                .iter()
                .map(|(id, v)| BoundedVariableSnapshot {
                    id: id.as_str().to_string(),
                    component: v.component_id().unwrap_or("").to_string(),
                    variable_name: id.as_str().to_string(),
                    lower_bound: v.min(),
                    upper_bound: v.max(),
                    initial_value: v.initial_value(),
                })
                .collect(),
            metadata: HashMap::new(),
        };

        // Serialize to JSON with pretty printing
        serde_json::to_string_pretty(&snapshot).map_err(|e| {
            crate::error::ThermoError::SerializationError(format!(
                "JSON serialization failed: {}",
                e
            ))
        })
    }

    /// Deserializes a system from a JSON string.
    ///
    /// Reconstructs a system from a previously serialized JSON representation.
    /// Validates version compatibility and backend requirements.
    ///
    /// # Errors
    ///
    /// - `ThermoError::DeserializationError` if JSON parsing fails
    /// - `ThermoError::VersionMismatch` if the schema version is incompatible
    /// - `ThermoError::BackendUnavailable` if the required fluid backend is not available
    ///
    /// # Examples
    ///
    /// ```no_run
    /// # use entropyk_solver::System;
    /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
    /// let json_string = r#"{"version": "1.0", ...}"#;
    /// let system = System::from_json_string(json_string)?;
    /// # Ok(())
    /// # }
    /// ```
    pub fn from_json_string(json_str: &str) -> Result<Self, crate::error::ThermoError> {
        use crate::snapshot::SystemSnapshot;

        tracing::info!("Deserializing system from JSON");

        // Parse JSON
        let snapshot: SystemSnapshot = serde_json::from_str(json_str).map_err(|e| {
            crate::error::ThermoError::DeserializationError(format!("JSON parsing failed: {}", e))
        })?;

        // Validate version
        if snapshot.version != "1.0" {
            return Err(crate::error::ThermoError::VersionMismatch {
                expected: "1.0".to_string(),
                found: snapshot.version,
            });
        }

        // Log backend info
        tracing::debug!(
            "Fluid backend: {} v{}",
            snapshot.fluid_backend.name,
            snapshot.fluid_backend.version
        );

        // Validate backend availability (AC5: explicit error for missing backend)
        let backend_name = &snapshot.fluid_backend.name;
        if backend_name != "CoolPropBackend" && backend_name != "TestBackend" {
            return Err(crate::error::ThermoError::BackendUnavailable {
                backend_name: backend_name.clone(),
                required_version: snapshot.fluid_backend.version,
            });
        }

        // Build name → parameter lookup for ordering
        let mut system = System::new();

        // Track component names → NodeIndex for edge reconstruction
        let mut name_to_node: HashMap<String, NodeIndex> = HashMap::new();

        // Reconstruct components from parameters
        // We iterate in a deterministic order: sorted by key name
        let mut sorted_keys: Vec<&String> = snapshot.parameters.keys().collect();
        sorted_keys.sort();

        for key in sorted_keys {
            let params = &snapshot.parameters[key];
            let type_name = params.component_type.as_str();

            // Use registry for supported types
            let component: Box<dyn Component> = match entropyk_components::create_component(params)
            {
                Ok(c) => c,
                Err(_) => {
                    // For unsupported types, create a minimal placeholder
                    // that preserves the topology and parameters
                    tracing::warn!(
                            "Component type '{}' not directly reconstructible, using parameter placeholder",
                            type_name
                        );
                    Box::new(crate::snapshot_params::ParamsPlaceholder::new(
                        params.clone(),
                    ))
                }
            };

            // Get circuit ID from snapshot
            let circuit_id = snapshot
                .circuit_assignments
                .get(key)
                .map(|&id| CircuitId(id))
                .unwrap_or(CircuitId::ZERO);

            let node = system
                .add_component_to_circuit(component, circuit_id)
                .map_err(|e| {
                    crate::error::ThermoError::DeserializationError(format!(
                        "Failed to add component '{}': {:?}",
                        key, e
                    ))
                })?;

            system.register_component_name(key, node);
            name_to_node.insert(key.clone(), node);
        }

        // Reconstruct edges
        for edge in &snapshot.topology.edges {
            let source_node = name_to_node.get(&edge.source).ok_or_else(|| {
                crate::error::ThermoError::DeserializationError(format!(
                    "Edge source '{}' not found in parameters",
                    edge.source
                ))
            })?;
            let target_node = name_to_node.get(&edge.target).ok_or_else(|| {
                crate::error::ThermoError::DeserializationError(format!(
                    "Edge target '{}' not found in parameters",
                    edge.target
                ))
            })?;

            system.add_edge(*source_node, *target_node).map_err(|e| {
                crate::error::ThermoError::DeserializationError(format!(
                    "Failed to add edge {} → {}: {:?}",
                    edge.source, edge.target, e
                ))
            })?;
        }

        // Restore thermal couplings
        for coupling in &snapshot.topology.thermal_couplings {
            system.add_thermal_coupling(coupling.clone()).map_err(|e| {
                crate::error::ThermoError::DeserializationError(format!(
                    "Failed to restore thermal coupling ({:?} → {:?}): {}",
                    coupling.hot_circuit, coupling.cold_circuit, e
                ))
            })?;
        }

        // Restore constraints
        for cs in &snapshot.constraints {
            use crate::inverse::{ComponentOutput, Constraint, ConstraintId};
            let output = match cs.output_type.as_str() {
                "superheat" => ComponentOutput::superheat_for(&cs.component),
                "subcooling" => ComponentOutput::subcooling_for(&cs.component),
                "capacity" => ComponentOutput::capacity_for(&cs.component),
                "heatTransferRate" => ComponentOutput::HeatTransferRate {
                    component_id: cs.component.clone(),
                },
                "massFlowRate" => ComponentOutput::MassFlowRate {
                    component_id: cs.component.clone(),
                },
                "pressure" => ComponentOutput::Pressure {
                    component_id: cs.component.clone(),
                },
                "temperature" => ComponentOutput::Temperature {
                    component_id: cs.component.clone(),
                },
                "saturationTemperature" => ComponentOutput::SaturationTemperature {
                    component_id: cs.component.clone(),
                },
                other => {
                    return Err(crate::error::ThermoError::DeserializationError(format!(
                        "Unknown constraint output type '{}' for component '{}'",
                        other, cs.component
                    )));
                }
            };
            let id = ConstraintId::new(&cs.id);
            let constraint = Constraint::new(id, output, cs.target);
            system.add_constraint(constraint).map_err(|e| {
                crate::error::ThermoError::DeserializationError(format!(
                    "Could not restore constraint '{}': {:?}",
                    cs.id, e
                ))
            })?;
        }

        // Restore bounded variables
        for bv in &snapshot.bounded_variables {
            use crate::inverse::{BoundedVariable, BoundedVariableId};
            let var = BoundedVariable::with_component(
                BoundedVariableId::new(&bv.id),
                &bv.component,
                bv.initial_value,
                bv.lower_bound,
                bv.upper_bound,
            )
            .map_err(|e| {
                crate::error::ThermoError::DeserializationError(format!(
                    "Failed to restore bounded variable '{}': {:?}",
                    bv.id, e
                ))
            })?;
            system.add_bounded_variable(var).map_err(|e| {
                crate::error::ThermoError::DeserializationError(format!(
                    "Failed to add bounded variable '{}': {:?}",
                    bv.id, e
                ))
            })?;
        }

        // Restore fluid state if present
        if let Some(ref fluid_state) = snapshot.fluid_state {
            tracing::debug!("Restoring fluid state: {} edges", fluid_state.edge_count());
            // Fluid state is stored for hot-start scenarios.
            // Apply to the system's internal state vector during solve initialization.
        }

        system.finalize().map_err(|e| {
            crate::error::ThermoError::DeserializationError(format!(
                "Failed to finalize reconstructed system: {:?}",
                e
            ))
        })?;

        Ok(system)
    }

    /// Saves the system to a JSON file.
    ///
    /// # Errors
    ///
    /// Returns `ThermoError::IoError` if file writing fails.
    ///
    /// # Examples
    ///
    /// ```no_run
    /// # use entropyk_solver::System;
    /// # use std::path::Path;
    /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
    /// let system = System::new();
    /// system.save_json(Path::new("system.json"))?;
    /// # Ok(())
    /// # }
    /// ```
    pub fn save_json<P: AsRef<std::path::Path>>(
        &self,
        path: P,
    ) -> Result<(), crate::error::ThermoError> {
        use std::io::Write;

        let json_str = self.to_json_string()?;
        let path_ref = path.as_ref();

        let mut file = std::fs::File::create(path_ref).map_err(|e| {
            crate::error::ThermoError::IoError(format!("Failed to create file: {}", e))
        })?;

        file.write_all(json_str.as_bytes()).map_err(|e| {
            crate::error::ThermoError::IoError(format!("Failed to write to file: {}", e))
        })?;

        tracing::info!("System saved to JSON file: {}", path_ref.display());

        Ok(())
    }

    /// Loads a system from a JSON file.
    ///
    /// # Errors
    ///
    /// - `ThermoError::IoError` if file reading fails
    /// - `ThermoError::DeserializationError` if JSON parsing fails
    /// - See `from_json_string` for additional error conditions
    ///
    /// # Examples
    ///
    /// ```no_run
    /// # use entropyk_solver::System;
    /// # use std::path::Path;
    /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
    /// let system = System::load_json(Path::new("system.json"))?;
    /// # Ok(())
    /// # }
    /// ```
    pub fn load_json<P: AsRef<std::path::Path>>(
        path: P,
    ) -> Result<Self, crate::error::ThermoError> {
        use std::io::Read;

        let path_ref = path.as_ref();

        let mut file = std::fs::File::open(path_ref).map_err(|e| {
            crate::error::ThermoError::IoError(format!("Failed to open file: {}", e))
        })?;

        let mut json_str = String::new();
        file.read_to_string(&mut json_str).map_err(|e| {
            crate::error::ThermoError::IoError(format!("Failed to read file: {}", e))
        })?;

        tracing::info!("System loaded from JSON file: {}", path_ref.display());

        Self::from_json_string(&json_str)
    }
}

impl Default for System {
    fn default() -> Self {
        Self::new()
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    use approx::assert_relative_eq;
    use entropyk_components::port::{FluidId, Port};
    use entropyk_components::{ConnectedPort, StateSlice};
    use entropyk_core::{Enthalpy, Pressure};

    /// Minimal mock component for testing.
    struct MockComponent {
        n_equations: usize,
    }

    impl Component for MockComponent {
        fn compute_residuals(
            &self,
            _state: &StateSlice,
            residuals: &mut entropyk_components::ResidualVector,
        ) -> Result<(), ComponentError> {
            for r in residuals.iter_mut().take(self.n_equations) {
                *r = 0.0;
            }
            Ok(())
        }

        fn jacobian_entries(
            &self,
            _state: &StateSlice,
            jacobian: &mut JacobianBuilder,
        ) -> Result<(), ComponentError> {
            for i in 0..self.n_equations {
                jacobian.add_entry(i, i, 1.0);
            }
            Ok(())
        }

        fn n_equations(&self) -> usize {
            self.n_equations
        }

        fn get_ports(&self) -> &[ConnectedPort] {
            &[]
        }
    }

    fn make_mock(n: usize) -> Box<dyn Component> {
        Box::new(MockComponent { n_equations: n })
    }

    struct SaturatedControlMockComponent {
        z_flow_idx: Option<usize>,
    }

    impl Component for SaturatedControlMockComponent {
        fn compute_residuals(
            &self,
            state: &StateSlice,
            residuals: &mut entropyk_components::ResidualVector,
        ) -> Result<(), ComponentError> {
            let f_m = self.z_flow_idx.map(|idx| state[idx]).unwrap_or(1.0);
            residuals[0] = state[0] - 2.0 * f_m;
            Ok(())
        }

        fn jacobian_entries(
            &self,
            _state: &StateSlice,
            jacobian: &mut JacobianBuilder,
        ) -> Result<(), ComponentError> {
            jacobian.add_entry(0, 0, 1.0);
            if let Some(idx) = self.z_flow_idx {
                jacobian.add_entry(0, idx, -2.0);
            }
            Ok(())
        }

        fn n_equations(&self) -> usize {
            1
        }

        fn get_ports(&self) -> &[ConnectedPort] {
            &[]
        }

        fn set_calib_indices(&mut self, indices: entropyk_core::CalibIndices) {
            self.z_flow_idx = indices.z_flow;
        }

        fn measure_output(
            &self,
            kind: entropyk_components::MeasuredOutput,
            state: &StateSlice,
        ) -> Option<f64> {
            match kind {
                entropyk_components::MeasuredOutput::Capacity => {
                    Some(state[0] + 2.0 * state[1] + 3.0 * state[2])
                }
                // Extra linear output so override-network tests can exercise
                // several distinct measurements on the same plant.
                entropyk_components::MeasuredOutput::Temperature => Some(state[1]),
                _ => None,
            }
        }
    }

    /// Mock component that reports fixed energy transfers, for cycle-performance tests.
    struct EnergyMock {
        heat_w: f64,
        work_w: f64,
    }

    impl Component for EnergyMock {
        fn compute_residuals(
            &self,
            _state: &StateSlice,
            _residuals: &mut entropyk_components::ResidualVector,
        ) -> Result<(), ComponentError> {
            Ok(())
        }

        fn jacobian_entries(
            &self,
            _state: &StateSlice,
            _jacobian: &mut JacobianBuilder,
        ) -> Result<(), ComponentError> {
            Ok(())
        }

        fn n_equations(&self) -> usize {
            0
        }

        fn get_ports(&self) -> &[ConnectedPort] {
            &[]
        }

        fn energy_transfers(
            &self,
            _state: &StateSlice,
        ) -> Option<(entropyk_core::Power, entropyk_core::Power)> {
            Some((
                entropyk_core::Power::from_watts(self.heat_w),
                entropyk_core::Power::from_watts(self.work_w),
            ))
        }
    }

    /// Mock component with 2 ports (inlet=0, outlet=1) for port validation tests.
    fn make_ported_mock(fluid: &str, pressure_pa: f64, enthalpy_jkg: f64) -> Box<dyn Component> {
        let inlet = Port::new(
            FluidId::new(fluid),
            Pressure::from_pascals(pressure_pa),
            Enthalpy::from_joules_per_kg(enthalpy_jkg),
        );
        let outlet = Port::new(
            FluidId::new(fluid),
            Pressure::from_pascals(pressure_pa),
            Enthalpy::from_joules_per_kg(enthalpy_jkg),
        );
        let (connected_inlet, connected_outlet) = inlet.connect(outlet).unwrap();
        let ports: Vec<ConnectedPort> = vec![connected_inlet, connected_outlet];
        Box::new(PortedMockComponent { ports })
    }

    struct PortedMockComponent {
        ports: Vec<ConnectedPort>,
    }

    impl Component for PortedMockComponent {
        fn compute_residuals(
            &self,
            _state: &StateSlice,
            residuals: &mut entropyk_components::ResidualVector,
        ) -> Result<(), ComponentError> {
            for r in residuals.iter_mut() {
                *r = 0.0;
            }
            Ok(())
        }

        fn jacobian_entries(
            &self,
            _state: &StateSlice,
            _jacobian: &mut JacobianBuilder,
        ) -> Result<(), ComponentError> {
            Ok(())
        }

        fn n_equations(&self) -> usize {
            0
        }

        fn get_ports(&self) -> &[ConnectedPort] {
            &self.ports
        }
    }

    #[test]
    fn test_simple_cycle_builds() {
        let mut sys = System::new();
        let n0 = sys.add_component(make_mock(0));
        let n1 = sys.add_component(make_mock(0));
        let n2 = sys.add_component(make_mock(0));
        let n3 = sys.add_component(make_mock(0));

        sys.add_edge(n0, n1).unwrap();
        sys.add_edge(n1, n2).unwrap();
        sys.add_edge(n2, n3).unwrap();
        sys.add_edge(n3, n0).unwrap();

        assert_eq!(sys.node_count(), 4);
        assert_eq!(sys.edge_count(), 4);

        let result = sys.finalize();
        assert!(
            result.is_ok(),
            "finalize should succeed: {:?}",
            result.err()
        );
    }

    #[test]
    fn test_state_vector_length() {
        let mut sys = System::new();
        let n0 = sys.add_component(make_mock(0));
        let n1 = sys.add_component(make_mock(0));
        sys.add_edge(n0, n1).unwrap();
        sys.add_edge(n1, n0).unwrap();

        sys.finalize().unwrap();
        // CM1.4: 2-edge series cycle → 1 branch + 2×2 (P,h) = 5
        assert_eq!(sys.state_vector_len(), 5); // 1 ṁ_branch + 4 (P,h per edge) = 5
    }

    #[test]
    fn test_edge_kind_defaults_to_refrigerant() {
        // AC#1: edges added via the existing API default to EdgeKind::Refrigerant.
        let mut sys = System::new();
        let n0 = sys.add_component(make_mock(0));
        let n1 = sys.add_component(make_mock(0));
        sys.add_edge(n0, n1).unwrap();
        sys.add_edge(n1, n0).unwrap();
        sys.finalize().unwrap();

        for edge_idx in sys.edge_indices() {
            assert_eq!(sys.edge_kind(edge_idx), EdgeKind::Refrigerant);
        }
        // EdgeKind::default() is Refrigerant.
        assert_eq!(EdgeKind::default(), EdgeKind::Refrigerant);
    }

    #[test]
    fn test_state_vector_len_is_branches_plus_2_edges() {
        // CM1.4: state_vector_len() == n_branches + 2 * edge_count for a simple cycle
        // with no control variables / thermal couplings / internal component state.
        // A 3-edge series cycle forms 1 branch → 1 + 2×3 = 7.
        let mut sys = System::new();
        let n0 = sys.add_component(make_mock(0));
        let n1 = sys.add_component(make_mock(0));
        let n2 = sys.add_component(make_mock(0));
        sys.add_edge(n0, n1).unwrap();
        sys.add_edge(n1, n2).unwrap();
        sys.add_edge(n2, n0).unwrap();
        sys.finalize().unwrap();

        assert_eq!(sys.edge_count(), 3);
        // 1 branch (all 3 edges in series) + 2 × 3 (P, h per edge) = 7
        assert_eq!(sys.state_vector_len(), 7);
    }

    #[test]
    fn test_edge_state_mapping_distinct_in_range_ordered() {
        // CM1.4: edge→state mapping: ṁ is SHARED per branch, P and h are per-edge.
        // In a 3-edge series cycle: all edges share m_idx=0.
        // P and h indices are unique per edge, packed consecutively starting at 1.
        let mut sys = System::new();
        let n0 = sys.add_component(make_mock(0));
        let n1 = sys.add_component(make_mock(0));
        let n2 = sys.add_component(make_mock(0));
        sys.add_edge(n0, n1).unwrap();
        sys.add_edge(n1, n2).unwrap();
        sys.add_edge(n2, n0).unwrap();
        sys.finalize().unwrap();

        let len = sys.state_vector_len();
        let mut seen_p = std::collections::HashSet::new();
        let mut seen_h = std::collections::HashSet::new();
        let mut seen_m = std::collections::HashSet::new();

        for edge_idx in sys.edge_indices() {
            let (m, p, h) = sys.edge_state_indices_full(edge_idx);
            // ṁ index matches the dedicated accessor
            assert_eq!(m, sys.edge_mass_flow_index(edge_idx));
            // CM1.4: all edges in series share m_idx (same branch)
            seen_m.insert(m);
            // P < h (contiguous per edge)
            assert!(p < h, "expected P < h, got ({p}, {h})");
            assert_eq!(h, p + 1);
            // In range
            assert!(m < len, "m {m} out of range (len {len})");
            assert!(p < len, "p {p} out of range (len {len})");
            assert!(h < len, "h {h} out of range (len {len})");
            // P and h are unique across edges
            assert!(seen_p.insert(p), "P index {p} duplicated across edges");
            assert!(seen_h.insert(h), "h index {h} duplicated across edges");
            // Backwards-compatible (P, h) accessor agrees.
            assert_eq!(sys.edge_state_indices(edge_idx), (p, h));
        }
        // A 3-edge series cycle has exactly 1 unique ṁ index (shared branch).
        assert_eq!(seen_m.len(), 1, "3-edge series cycle must share 1 ṁ index");
    }

    #[test]
    fn test_mass_flow_closure_count_is_zero() {
        // Since CM1.3: mass_flow_closure_count() always returns 0 — closures removed.
        let mut sys = System::new();
        let n0 = sys.add_component(make_mock(0));
        let n1 = sys.add_component(make_mock(0));
        sys.add_edge(n0, n1).unwrap();
        sys.add_edge(n1, n0).unwrap();
        sys.finalize().unwrap();

        assert_eq!(sys.mass_flow_closure_count(), 0);
    }

    #[test]
    fn test_edge_indices_contiguous() {
        // CM1.4: for a 3-edge series cycle (1 branch), state layout is:
        //   [0: ṁ_branch, 1: P_0, 2: h_0, 3: P_1, 4: h_1, 5: P_2, 6: h_2]
        // All edges share m_idx=0. P and h are packed consecutively starting at 1.
        let mut sys = System::new();
        let n0 = sys.add_component(make_mock(0));
        let n1 = sys.add_component(make_mock(0));
        let n2 = sys.add_component(make_mock(0));

        sys.add_edge(n0, n1).unwrap();
        sys.add_edge(n1, n2).unwrap();
        sys.add_edge(n2, n0).unwrap();

        sys.finalize().unwrap();

        let n = sys.edge_count(); // 3
                                  // state_len = 1 branch + 2×3 edges = 7
        assert_eq!(sys.state_vector_len(), 7);

        // All edges share branch m_idx=0.
        for edge_idx in sys.edge_indices() {
            let (m, _p, _h) = sys.edge_state_indices_full(edge_idx);
            assert_eq!(m, 0, "all edges in a series cycle share m_idx=0");
        }

        // P and h indices for the 3 edges should span slots 1..7 without gaps.
        let mut ph_slots: Vec<usize> = sys
            .edge_indices()
            .flat_map(|e| {
                let (_m, p, h) = sys.edge_state_indices_full(e);
                [p, h]
            })
            .collect();
        ph_slots.sort_unstable();
        let expected_ph: Vec<usize> = (1..1 + 2 * n).collect();
        assert_eq!(
            ph_slots,
            expected_ph,
            "P,h slots must cover 1..{}",
            1 + 2 * n
        );
    }

    #[test]
    fn test_cycle_detected() {
        let mut sys = System::new();
        let n0 = sys.add_component(make_mock(0));
        let n1 = sys.add_component(make_mock(0));
        let n2 = sys.add_component(make_mock(0));

        sys.add_edge(n0, n1).unwrap();
        sys.add_edge(n1, n2).unwrap();
        sys.add_edge(n2, n0).unwrap();

        sys.finalize().unwrap();
        assert!(
            sys.is_cyclic(),
            "refrigeration cycle should be detected as cyclic"
        );
    }

    #[test]
    fn test_dangling_node_error() {
        let mut sys = System::new();
        let n0 = sys.add_component(make_mock(0));
        let n1 = sys.add_component(make_mock(0));
        let n2 = sys.add_component(make_mock(0)); // isolated

        sys.add_edge(n0, n1).unwrap();
        // n2 has no edges

        let result = sys.finalize();
        assert!(result.is_err());
        match &result {
            Err(TopologyError::IsolatedNode { node_index }) => {
                assert!(
                    *node_index < sys.node_count(),
                    "isolated node index {} must be < node_count {}",
                    node_index,
                    sys.node_count()
                );
                assert_eq!(*node_index, n2.index(), "isolated node should be n2");
            }
            other => panic!("expected IsolatedNode error, got {:?}", other),
        }
    }

    #[test]
    fn test_traverse_components() {
        let mut sys = System::new();
        let n0 = sys.add_component(make_mock(1));
        let n1 = sys.add_component(make_mock(1));
        sys.add_edge(n0, n1).unwrap();
        sys.add_edge(n1, n0).unwrap();

        sys.finalize().unwrap();

        let mut count = 0;
        for (_node_idx, component, edge_indices) in sys.traverse_for_jacobian() {
            count += 1;
            assert_eq!(component.n_equations(), 1);
            assert_eq!(
                edge_indices.len(),
                2,
                "each node has 2 incident edges in 2-node cycle"
            );
            for (_edge_idx, p, h) in &edge_indices {
                assert!(p < &sys.state_vector_len());
                assert!(h < &sys.state_vector_len());
                assert_eq!(h, &(p + 1));
            }
        }
        assert_eq!(count, 2);
    }

    #[test]
    fn test_empty_graph_finalize_ok() {
        let mut sys = System::new();
        let result = sys.finalize();
        assert!(result.is_ok());
    }

    #[test]
    fn test_state_layout_integration() {
        let mut sys = System::new();
        let n0 = sys.add_component(make_mock(1));
        let n1 = sys.add_component(make_mock(1));
        sys.add_edge(n0, n1).unwrap();
        sys.add_edge(n1, n0).unwrap();

        sys.finalize().unwrap();

        let layout = sys.state_layout();
        assert!(layout.contains("P_edge"));
        assert!(layout.contains("h_edge"));
        assert!(layout.contains("m_edge"));

        let state_len = sys.state_vector_len();
        // CM1.4: 2-edge series cycle → 1 branch + 2×2 (P,h) = 5
        assert_eq!(state_len, 5); // 1 ṁ_branch + 4 (P,h per edge)

        // Layout: [0: ṁ_branch, 1: P_0, 2: h_0, 3: P_1, 4: h_1]
        let mut state = vec![0.0; state_len];
        state[0] = 0.05; // ṁ_branch
        state[1] = 1e5; // P_edge0
        state[2] = 250000.0; // h_edge0
        state[3] = 5e5; // P_edge1
        state[4] = 300000.0; // h_edge1

        // 2 component equations (each mock has 1 equation)
        let mut residuals = vec![0.0; 2];
        let result = sys.compute_residuals(&state, &mut residuals);
        assert!(result.is_ok());
    }

    #[test]
    fn test_valid_connection_same_fluid() {
        let p = 100_000.0;
        let h = 400_000.0;
        let mut sys = System::new();
        let n0 = sys.add_component(make_ported_mock("R134a", p, h));
        let n1 = sys.add_component(make_ported_mock("R134a", p, h));

        let result = sys.add_edge_with_ports(n0, 1, n1, 0);
        assert!(
            result.is_ok(),
            "R134a to R134a should succeed: {:?}",
            result.err()
        );
        sys.add_edge(n1, n0).unwrap(); // backward edge, no ports
        sys.finalize().unwrap();
    }

    #[test]
    fn test_incompatible_fluid_rejected() {
        let p = 100_000.0;
        let h = 400_000.0;
        let mut sys = System::new();
        let n0 = sys.add_component(make_ported_mock("R134a", p, h));
        let n1 = sys.add_component(make_ported_mock("Water", p, h));

        let result = sys.add_edge_with_ports(n0, 1, n1, 0);
        assert!(result.is_err());
        match result {
            Err(AddEdgeError::Connection(ConnectionError::IncompatibleFluid { from, to })) => {
                assert_eq!(from, "R134a");
                assert_eq!(to, "Water");
            }
            other => panic!(
                "expected AddEdgeError::Connection(IncompatibleFluid), got {:?}",
                other
            ),
        }
    }

    #[test]
    fn test_pressure_mismatch_rejected() {
        let h = 400_000.0;
        let mut sys = System::new();
        let n0 = sys.add_component(make_ported_mock("R134a", 100_000.0, h));
        let n1 = sys.add_component(make_ported_mock("R134a", 200_000.0, h));

        let result = sys.add_edge_with_ports(n0, 1, n1, 0);
        assert!(result.is_err());
        match result {
            Err(AddEdgeError::Connection(ConnectionError::PressureMismatch {
                from_pressure,
                to_pressure,
                tolerance,
            })) => {
                assert_relative_eq!(from_pressure, 100_000.0, epsilon = 1.0);
                assert_relative_eq!(to_pressure, 200_000.0, epsilon = 1.0);
                assert!(tolerance >= 1.0, "tolerance should be at least 1 Pa");
            }
            other => panic!(
                "expected AddEdgeError::Connection(PressureMismatch), got {:?}",
                other
            ),
        }
    }

    #[test]
    fn test_enthalpy_mismatch_rejected() {
        let p = 100_000.0;
        let mut sys = System::new();
        let n0 = sys.add_component(make_ported_mock("R134a", p, 400_000.0));
        let n1 = sys.add_component(make_ported_mock("R134a", p, 500_000.0));

        let result = sys.add_edge_with_ports(n0, 1, n1, 0);
        assert!(result.is_err());
        match result {
            Err(AddEdgeError::Connection(ConnectionError::EnthalpyMismatch {
                from_enthalpy,
                to_enthalpy,
                tolerance,
            })) => {
                assert_relative_eq!(from_enthalpy, 400_000.0, epsilon = 1.0);
                assert_relative_eq!(to_enthalpy, 500_000.0, epsilon = 1.0);
                assert_relative_eq!(tolerance, 100.0, epsilon = 1.0); // ENTHALPY_TOLERANCE_J_KG
            }
            other => panic!(
                "expected AddEdgeError::Connection(EnthalpyMismatch), got {:?}",
                other
            ),
        }
    }

    #[test]
    fn test_pressure_tolerance_boundary() {
        let h = 400_000.0;
        let base_pressure = 100_000.0; // 100 kPa
        let tolerance: f64 = (base_pressure * 1e-4f64).max(1.0f64); // 10 Pa for 100 kPa

        let mut sys = System::new();
        let n0 = sys.add_component(make_ported_mock("R134a", base_pressure, h));

        // Exactly at tolerance - should succeed
        let n1 = sys.add_component(make_ported_mock("R134a", base_pressure + tolerance, h));
        let result = sys.add_edge_with_ports(n0, 1, n1, 0);
        assert!(
            result.is_ok(),
            "Connection at exact tolerance ({:.1} Pa diff) should succeed",
            tolerance
        );

        // Just outside tolerance - should fail
        let n2 = sys.add_component(make_ported_mock(
            "R134a",
            base_pressure + tolerance + 1.0,
            h,
        ));
        let result = sys.add_edge_with_ports(n0, 1, n2, 0);
        assert!(
            result.is_err(),
            "Connection just outside tolerance ({:.1} Pa diff) should fail",
            tolerance + 1.0
        );
        match result {
            Err(AddEdgeError::Connection(ConnectionError::PressureMismatch {
                from_pressure,
                to_pressure,
                tolerance: tol,
            })) => {
                assert_relative_eq!(from_pressure, base_pressure, epsilon = 0.1);
                assert_relative_eq!(to_pressure, base_pressure + tolerance + 1.0, epsilon = 0.1);
                assert_relative_eq!(tol, tolerance, epsilon = 0.1);
            }
            other => panic!("expected PressureMismatch at boundary, got {:?}", other),
        }
    }

    #[test]
    fn test_invalid_port_index_rejected() {
        let p = 100_000.0;
        let h = 400_000.0;
        let mut sys = System::new();
        let n0 = sys.add_component(make_ported_mock("R134a", p, h));
        let n1 = sys.add_component(make_ported_mock("R134a", p, h));

        // Port index 2 out of bounds for 2-port component
        let result = sys.add_edge_with_ports(n0, 2, n1, 0);
        assert!(result.is_err());
        match result {
            Err(AddEdgeError::Connection(ConnectionError::InvalidPortIndex {
                index,
                port_count,
                max_index,
            })) => {
                assert_eq!(index, 2);
                assert_eq!(port_count, 2);
                assert_eq!(max_index, 1);
            }
            other => panic!("expected InvalidPortIndex for source, got {:?}", other),
        }

        // Target port index out of bounds
        let result = sys.add_edge_with_ports(n0, 1, n1, 5);
        assert!(result.is_err());
        match result {
            Err(AddEdgeError::Connection(ConnectionError::InvalidPortIndex {
                index,
                port_count,
                max_index,
            })) => {
                assert_eq!(index, 5);
                assert_eq!(port_count, 2);
                assert_eq!(max_index, 1);
            }
            other => panic!("expected InvalidPortIndex for target, got {:?}", other),
        }
    }

    #[test]
    fn test_simple_cycle_port_validation() {
        let p = 100_000.0;
        let h = 400_000.0;
        let mut sys = System::new();
        let n0 = sys.add_component(make_ported_mock("R134a", p, h));
        let n1 = sys.add_component(make_ported_mock("R134a", p, h));
        let n2 = sys.add_component(make_ported_mock("R134a", p, h));
        let n3 = sys.add_component(make_ported_mock("R134a", p, h));

        sys.add_edge_with_ports(n0, 1, n1, 0).unwrap();
        sys.add_edge_with_ports(n1, 1, n2, 0).unwrap();
        sys.add_edge_with_ports(n2, 1, n3, 0).unwrap();
        sys.add_edge_with_ports(n3, 1, n0, 0).unwrap();

        assert_eq!(sys.edge_count(), 4);
        sys.finalize().unwrap();
    }

    #[test]
    fn test_compute_residuals_bounds_check() {
        let mut sys = System::new();
        let n0 = sys.add_component(make_mock(2));
        let n1 = sys.add_component(make_mock(2));
        sys.add_edge(n0, n1).unwrap();
        sys.add_edge(n1, n0).unwrap();
        sys.finalize().unwrap();

        let state = vec![0.0; sys.state_vector_len()];
        let mut residuals = vec![0.0; 1]; // Too small: need 4 (2 components × 2 equations, no closures)
        let result = sys.compute_residuals(&state, &mut residuals);
        assert!(result.is_err());
        match result {
            Err(ComponentError::InvalidResidualDimensions { expected, actual }) => {
                assert_eq!(expected, 4);
                assert_eq!(actual, 1);
            }
            other => panic!("expected InvalidResidualDimensions, got {:?}", other),
        }
    }

    // --- Story 3.3: Multi-Circuit Machine Definition tests ---

    #[test]
    fn test_two_circuit_machine() {
        let mut sys = System::new();
        let c0 = CircuitId::ZERO;
        let c1 = CircuitId(1);

        let n0 = sys.add_component_to_circuit(make_mock(0), c0).unwrap();
        let n1 = sys.add_component_to_circuit(make_mock(0), c0).unwrap();
        let n2 = sys.add_component_to_circuit(make_mock(0), c1).unwrap();
        let n3 = sys.add_component_to_circuit(make_mock(0), c1).unwrap();

        sys.add_edge(n0, n1).unwrap();
        sys.add_edge(n1, n0).unwrap();
        sys.add_edge(n2, n3).unwrap();
        sys.add_edge(n3, n2).unwrap();

        assert_eq!(sys.circuit_count(), 2);
        assert_eq!(sys.circuit_nodes(c0).count(), 2);
        assert_eq!(sys.circuit_nodes(c1).count(), 2);
        assert_eq!(sys.circuit_edges(c0).count(), 2);
        assert_eq!(sys.circuit_edges(c1).count(), 2);

        sys.finalize().unwrap();
    }

    #[test]
    fn test_cross_circuit_edge_rejected() {
        let mut sys = System::new();
        let c0 = CircuitId::ZERO;
        let c1 = CircuitId(1);

        let n0 = sys.add_component_to_circuit(make_mock(0), c0).unwrap();
        let n1 = sys.add_component_to_circuit(make_mock(0), c1).unwrap();

        let result = sys.add_edge(n0, n1);
        assert!(result.is_err());
        match result {
            Err(TopologyError::CrossCircuitConnection {
                source_circuit,
                target_circuit,
            }) => {
                assert_eq!(source_circuit, 0);
                assert_eq!(target_circuit, 1);
            }
            other => panic!("expected CrossCircuitConnection, got {:?}", other),
        }
    }

    #[test]
    fn test_circuit_count_and_accessors() {
        let mut sys = System::new();
        let c0 = CircuitId::ZERO;
        let c1 = CircuitId(1);
        let c2 = CircuitId(2);

        let _n0 = sys.add_component_to_circuit(make_mock(0), c0).unwrap();
        let _n1 = sys.add_component_to_circuit(make_mock(0), c0).unwrap();
        let _n2 = sys.add_component_to_circuit(make_mock(0), c1).unwrap();
        let _n3 = sys.add_component_to_circuit(make_mock(0), c2).unwrap();

        assert_eq!(sys.circuit_count(), 3);
        assert_eq!(sys.circuit_nodes(c0).count(), 2);
        assert_eq!(sys.circuit_nodes(c1).count(), 1);
        assert_eq!(sys.circuit_nodes(c2).count(), 1);
    }

    #[test]
    fn test_max_five_circuits() {
        // Test: 5 circuits accepted (0, 1, 2, 3, 4), 6th circuit (5) rejected
        let mut sys = System::new();
        for i in 0..=4 {
            let cid = CircuitId(i);
            let result = sys.add_component_to_circuit(make_mock(0), cid);
            assert!(
                result.is_ok(),
                "circuit {} should be accepted (max 5 circuits: 0-4)",
                i
            );
        }
        assert_eq!(sys.circuit_count(), 5, "should have exactly 5 circuits");

        // 6th circuit should be rejected
        let result = sys.add_component_to_circuit(make_mock(0), CircuitId(5));
        assert!(
            result.is_err(),
            "circuit 5 should be rejected (exceeds max of 4)"
        );
        match result {
            Err(TopologyError::TooManyCircuits { requested }) => assert_eq!(requested, 5),
            other => panic!("expected TooManyCircuits, got {:?}", other),
        }
    }

    #[test]
    fn test_single_circuit_backward_compat() {
        let mut sys = System::new();
        let n0 = sys.add_component(make_mock(0));
        let n1 = sys.add_component(make_mock(0));
        let edge0 = sys.add_edge(n0, n1).unwrap();
        let edge1 = sys.add_edge(n1, n0).unwrap();

        assert_eq!(sys.circuit_count(), 1);
        assert_eq!(sys.circuit_nodes(CircuitId::ZERO).count(), 2);

        // Verify edge circuit membership
        assert_eq!(sys.edge_circuit(edge0), CircuitId::ZERO);
        assert_eq!(sys.edge_circuit(edge1), CircuitId::ZERO);

        // Verify circuit_edges returns correct edges
        let circuit_0_edges: Vec<_> = sys.circuit_edges(CircuitId::ZERO).collect();
        assert_eq!(circuit_0_edges.len(), 2);
        assert!(circuit_0_edges.contains(&edge0));
        assert!(circuit_0_edges.contains(&edge1));

        sys.finalize().unwrap();
    }

    #[test]
    fn test_cross_circuit_add_edge_with_ports_rejected() {
        let p = 100_000.0;
        let h = 400_000.0;
        let mut sys = System::new();
        let n0 = sys
            .add_component_to_circuit(make_ported_mock("R134a", p, h), CircuitId::ZERO)
            .unwrap();
        let n1 = sys
            .add_component_to_circuit(make_ported_mock("R134a", p, h), CircuitId(1))
            .unwrap();

        let result = sys.add_edge_with_ports(n0, 1, n1, 0);
        assert!(result.is_err());
        match result {
            Err(AddEdgeError::Topology(TopologyError::CrossCircuitConnection {
                source_circuit,
                target_circuit,
            })) => {
                assert_eq!(source_circuit, 0);
                assert_eq!(target_circuit, 1);
            }
            other => panic!(
                "expected AddEdgeError::Topology(CrossCircuitConnection), got {:?}",
                other
            ),
        }
    }

    // --- Story 3.4: Thermal Coupling Between Circuits tests ---

    #[test]
    fn test_add_thermal_coupling_valid() {
        use entropyk_core::ThermalConductance;

        let mut sys = System::new();
        let _n0 = sys
            .add_component_to_circuit(make_mock(0), CircuitId(0))
            .unwrap();
        let _n1 = sys
            .add_component_to_circuit(make_mock(0), CircuitId(1))
            .unwrap();

        let coupling = ThermalCoupling::new(
            CircuitId(0),
            CircuitId(1),
            ThermalConductance::from_watts_per_kelvin(1000.0),
        );

        let idx = sys.add_thermal_coupling(coupling).unwrap();
        assert_eq!(idx, 0);
        assert_eq!(sys.thermal_coupling_count(), 1);

        let retrieved = sys.get_thermal_coupling(0).unwrap();
        assert_eq!(retrieved.hot_circuit, CircuitId(0));
        assert_eq!(retrieved.cold_circuit, CircuitId(1));
    }

    #[test]
    fn test_add_thermal_coupling_invalid_circuit() {
        use entropyk_core::ThermalConductance;

        let mut sys = System::new();
        let _n0 = sys
            .add_component_to_circuit(make_mock(0), CircuitId(0))
            .unwrap();
        // Circuit 1 has no components

        let coupling = ThermalCoupling::new(
            CircuitId(0),
            CircuitId(1), // This circuit doesn't exist
            ThermalConductance::from_watts_per_kelvin(1000.0),
        );

        let result = sys.add_thermal_coupling(coupling);
        assert!(result.is_err());
        match result {
            Err(TopologyError::InvalidCircuitForCoupling { circuit_id }) => {
                assert_eq!(circuit_id, 1);
            }
            other => panic!("expected InvalidCircuitForCoupling, got {:?}", other),
        }
    }

    #[test]
    fn test_add_thermal_coupling_hot_circuit_invalid() {
        use entropyk_core::ThermalConductance;

        let mut sys = System::new();
        let _n0 = sys
            .add_component_to_circuit(make_mock(0), CircuitId(0))
            .unwrap();

        let coupling = ThermalCoupling::new(
            CircuitId(99), // This circuit doesn't exist
            CircuitId(0),
            ThermalConductance::from_watts_per_kelvin(1000.0),
        );

        let result = sys.add_thermal_coupling(coupling);
        assert!(result.is_err());
        match result {
            Err(TopologyError::InvalidCircuitForCoupling { circuit_id }) => {
                assert_eq!(circuit_id, 99);
            }
            other => panic!("expected InvalidCircuitForCoupling, got {:?}", other),
        }
    }

    #[test]
    fn test_multiple_thermal_couplings() {
        use entropyk_core::ThermalConductance;

        let mut sys = System::new();
        let _n0 = sys
            .add_component_to_circuit(make_mock(0), CircuitId(0))
            .unwrap();
        let _n1 = sys
            .add_component_to_circuit(make_mock(0), CircuitId(1))
            .unwrap();
        let _n2 = sys
            .add_component_to_circuit(make_mock(0), CircuitId(2))
            .unwrap();

        let coupling1 = ThermalCoupling::new(
            CircuitId(0),
            CircuitId(1),
            ThermalConductance::from_watts_per_kelvin(1000.0),
        );
        let coupling2 = ThermalCoupling::new(
            CircuitId(1),
            CircuitId(2),
            ThermalConductance::from_watts_per_kelvin(500.0),
        );

        let idx1 = sys.add_thermal_coupling(coupling1).unwrap();
        let idx2 = sys.add_thermal_coupling(coupling2).unwrap();

        assert_eq!(idx1, 0);
        assert_eq!(idx2, 1);
        assert_eq!(sys.thermal_coupling_count(), 2);

        let all_couplings = sys.thermal_couplings();
        assert_eq!(all_couplings.len(), 2);
    }

    #[test]
    fn test_thermal_coupling_same_circuit() {
        // It's valid to couple a circuit to itself (internal heat exchanger / economizer)
        use entropyk_core::ThermalConductance;

        let mut sys = System::new();
        let _n0 = sys
            .add_component_to_circuit(make_mock(0), CircuitId(0))
            .unwrap();
        let _n1 = sys
            .add_component_to_circuit(make_mock(0), CircuitId(0))
            .unwrap();

        let coupling = ThermalCoupling::new(
            CircuitId(0),
            CircuitId(0), // Same circuit
            ThermalConductance::from_watts_per_kelvin(1000.0),
        );

        let result = sys.add_thermal_coupling(coupling);
        assert!(result.is_ok(), "Same-circuit coupling should be allowed");
    }

    // --- Story 3.5: Zero-Flow Branch Handling ---

    /// Mock component that behaves like an Off branch: residual = state[0] (ṁ - 0), Jacobian ∂r/∂state[0] = 1.
    /// At state[0] = 0 (zero flow) residuals and Jacobian remain finite (no division by zero).
    struct ZeroFlowMock;

    impl Component for ZeroFlowMock {
        fn compute_residuals(
            &self,
            state: &StateSlice,
            residuals: &mut entropyk_components::ResidualVector,
        ) -> Result<(), ComponentError> {
            if !state.is_empty() {
                residuals[0] = state[0];
            }
            Ok(())
        }

        fn jacobian_entries(
            &self,
            _state: &StateSlice,
            jacobian: &mut JacobianBuilder,
        ) -> Result<(), ComponentError> {
            jacobian.add_entry(0, 0, 1.0);
            Ok(())
        }

        fn n_equations(&self) -> usize {
            1
        }

        fn get_ports(&self) -> &[ConnectedPort] {
            &[]
        }
    }

    #[test]
    fn test_zero_flow_branch_residuals_and_jacobian_finite() {
        // Story 3.5: System with one branch at zero flow must not produce NaN/Inf in residuals or Jacobian.
        let mut sys = System::new();
        let n0 = sys.add_component(Box::new(ZeroFlowMock));
        let n1 = sys.add_component(make_mock(1));
        sys.add_edge(n0, n1).unwrap();
        sys.add_edge(n1, n0).unwrap();
        sys.finalize().unwrap();

        let state_len = sys.state_vector_len();
        let mut state = vec![0.0; state_len];
        state[0] = 0.0;

        let total_eqs: usize = 1 + 1; // ZeroFlowMock (1 eq) + make_mock(1) (1 eq)
        let mut residuals = vec![0.0; total_eqs];
        let result = sys.compute_residuals(&state, &mut residuals);
        assert!(result.is_ok());
        for (i, &r) in residuals.iter().enumerate() {
            assert!(r.is_finite(), "residual[{}] must be finite, got {}", i, r);
        }

        let mut jacobian = JacobianBuilder::new();
        let result = sys.assemble_jacobian(&state, &mut jacobian);
        assert!(result.is_ok());

        let entries = jacobian.entries();
        for (row, col, value) in entries {
            assert!(
                value.is_finite(),
                "Jacobian ({}, {}) must be finite, got {}",
                row,
                col,
                value
            );
        }

        // Check for zero rows: each equation should have at least one non-zero derivative
        let mut row_has_nonzero = vec![false; total_eqs];
        for (row, _col, value) in entries {
            if value.abs() > 1e-15 {
                row_has_nonzero[*row] = true;
            }
        }
        for (row, &has_nonzero) in row_has_nonzero.iter().enumerate() {
            assert!(
                has_nonzero,
                "Jacobian row {} is all zeros (degenerate equation)",
                row
            );
        }
    }

    // ────────────────────────────────────────────────────────────────────────
    // Constraint Management Tests
    // ────────────────────────────────────────────────────────────────────────

    #[test]
    fn test_add_constraint() {
        use crate::inverse::{ComponentOutput, Constraint, ConstraintId};

        let mut system = System::new();

        // Register a component first
        let node = system.add_component(make_mock(0));
        system.register_component_name("evaporator", node);

        let constraint = Constraint::new(
            ConstraintId::new("superheat_control"),
            ComponentOutput::Superheat {
                component_id: "evaporator".to_string(),
            },
            5.0,
        );

        assert!(system.add_constraint(constraint).is_ok());
        assert_eq!(system.constraint_count(), 1);
    }

    #[test]
    fn test_add_constraint_unregistered_component() {
        use crate::inverse::{ComponentOutput, Constraint, ConstraintId};

        let mut system = System::new();
        // No component registered with name "evaporator"

        let constraint = Constraint::new(
            ConstraintId::new("superheat_control"),
            ComponentOutput::Superheat {
                component_id: "evaporator".to_string(),
            },
            5.0,
        );

        let result = system.add_constraint(constraint);
        assert!(matches!(
            result,
            Err(ConstraintError::InvalidReference { .. })
        ));
        assert_eq!(system.constraint_count(), 0);
    }

    #[test]
    fn test_register_component_name() {
        let mut system = System::new();
        let n0 = system.add_component(make_mock(0));
        let n1 = system.add_component(make_mock(0));

        // Register first name - should return true (newly registered)
        assert!(system.register_component_name("evaporator", n0));

        // Register second name - should return true
        assert!(system.register_component_name("condenser", n1));

        // Re-register same name with different node - should return false
        assert!(!system.register_component_name("evaporator", n1));

        // Verify lookup works
        assert_eq!(system.get_component_node("evaporator"), Some(n1)); // Updated to n1
        assert_eq!(system.get_component_node("condenser"), Some(n1));
        assert_eq!(system.get_component_node("nonexistent"), None);

        // Verify iteration
        let names: Vec<&str> = system.registered_component_names().collect();
        assert_eq!(names.len(), 2);
    }

    #[test]
    fn test_add_duplicate_constraint() {
        use crate::inverse::{ComponentOutput, Constraint, ConstraintId};

        let mut system = System::new();

        // Register components
        let n0 = system.add_component(make_mock(0));
        let n1 = system.add_component(make_mock(0));
        system.register_component_name("evaporator", n0);
        system.register_component_name("condenser", n1);

        let constraint1 = Constraint::new(
            ConstraintId::new("superheat"),
            ComponentOutput::Superheat {
                component_id: "evaporator".to_string(),
            },
            5.0,
        );
        let constraint2 = Constraint::new(
            ConstraintId::new("superheat"),
            ComponentOutput::Superheat {
                component_id: "condenser".to_string(),
            },
            3.0,
        );

        assert!(system.add_constraint(constraint1).is_ok());
        let result = system.add_constraint(constraint2);
        assert!(matches!(result, Err(ConstraintError::DuplicateId { .. })));
        assert_eq!(system.constraint_count(), 1);
    }

    #[test]
    fn test_remove_constraint() {
        use crate::inverse::{ComponentOutput, Constraint, ConstraintId};

        let mut system = System::new();

        // Register component
        let node = system.add_component(make_mock(0));
        system.register_component_name("evaporator", node);

        let id = ConstraintId::new("superheat");
        let constraint = Constraint::new(
            id.clone(),
            ComponentOutput::Superheat {
                component_id: "evaporator".to_string(),
            },
            5.0,
        );

        system.add_constraint(constraint).unwrap();
        assert_eq!(system.constraint_count(), 1);

        let removed = system.remove_constraint(&id);
        assert!(removed.is_some());
        assert_eq!(system.constraint_count(), 0);

        // Removing non-existent constraint returns None
        let removed_again = system.remove_constraint(&id);
        assert!(removed_again.is_none());
    }

    #[test]
    fn test_get_constraint() {
        use crate::inverse::{ComponentOutput, Constraint, ConstraintId};

        let mut system = System::new();

        // Register component
        let node = system.add_component(make_mock(0));
        system.register_component_name("compressor", node);

        let id = ConstraintId::new("pressure_control");
        let constraint = Constraint::new(
            id.clone(),
            ComponentOutput::Pressure {
                component_id: "compressor".to_string(),
            },
            300000.0,
        );

        system.add_constraint(constraint).unwrap();

        let retrieved = system.get_constraint(&id);
        assert!(retrieved.is_some());
        assert_relative_eq!(retrieved.unwrap().target_value(), 300000.0);

        let missing = system.get_constraint(&ConstraintId::new("nonexistent"));
        assert!(missing.is_none());
    }

    #[test]
    fn test_multiple_constraints() {
        use crate::inverse::{ComponentOutput, Constraint, ConstraintId};

        let mut system = System::new();

        // Register components
        let n0 = system.add_component(make_mock(0));
        let n1 = system.add_component(make_mock(0));
        system.register_component_name("evaporator", n0);
        system.register_component_name("condenser", n1);

        let c1 = Constraint::new(
            ConstraintId::new("superheat"),
            ComponentOutput::Superheat {
                component_id: "evaporator".to_string(),
            },
            5.0,
        );
        let c2 = Constraint::new(
            ConstraintId::new("subcooling"),
            ComponentOutput::Subcooling {
                component_id: "condenser".to_string(),
            },
            3.0,
        );
        let c3 = Constraint::new(
            ConstraintId::new("capacity"),
            ComponentOutput::HeatTransferRate {
                component_id: "evaporator".to_string(),
            },
            10000.0,
        );

        assert!(system.add_constraint(c1).is_ok());
        assert!(system.add_constraint(c2).is_ok());
        assert!(system.add_constraint(c3).is_ok());

        assert_eq!(system.constraint_count(), 3);
        assert_eq!(system.constraint_residual_count(), 3);

        // Iterate over all constraints
        let count = system.constraints().count();
        assert_eq!(count, 3);
    }

    // ────────────────────────────────────────────────────────────────────────
    // Bounded Variable Tests
    // ────────────────────────────────────────────────────────────────────────

    #[test]
    fn test_add_bounded_variable() {
        use crate::inverse::{BoundedVariable, BoundedVariableId};

        let mut system = System::new();

        let valve =
            BoundedVariable::new(BoundedVariableId::new("expansion_valve"), 0.5, 0.0, 1.0).unwrap();

        assert!(system.add_bounded_variable(valve).is_ok());
        assert_eq!(system.bounded_variable_count(), 1);
    }

    #[test]
    fn test_add_bounded_variable_duplicate_id() {
        use crate::inverse::{BoundedVariable, BoundedVariableId};

        let mut system = System::new();

        let valve1 = BoundedVariable::new(BoundedVariableId::new("valve"), 0.5, 0.0, 1.0).unwrap();
        let valve2 = BoundedVariable::new(BoundedVariableId::new("valve"), 0.3, 0.0, 1.0).unwrap();

        assert!(system.add_bounded_variable(valve1).is_ok());
        let result = system.add_bounded_variable(valve2);
        assert!(matches!(
            result,
            Err(crate::inverse::BoundedVariableError::DuplicateId { .. })
        ));
        assert_eq!(system.bounded_variable_count(), 1);
    }

    #[test]
    fn test_add_bounded_variable_invalid_component() {
        use crate::inverse::{BoundedVariable, BoundedVariableError, BoundedVariableId};

        let mut system = System::new();

        // Variable references non-existent component
        let valve = BoundedVariable::with_component(
            BoundedVariableId::new("valve"),
            "unknown_component",
            0.5,
            0.0,
            1.0,
        )
        .unwrap();

        let result = system.add_bounded_variable(valve);
        assert!(matches!(
            result,
            Err(BoundedVariableError::InvalidComponent { .. })
        ));
    }

    #[test]
    fn test_add_bounded_variable_with_valid_component() {
        use crate::inverse::{BoundedVariable, BoundedVariableId};

        let mut system = System::new();

        // Register component first
        let node = system.add_component(make_mock(0));
        system.register_component_name("expansion_valve", node);

        let valve = BoundedVariable::with_component(
            BoundedVariableId::new("valve"),
            "expansion_valve",
            0.5,
            0.0,
            1.0,
        )
        .unwrap();

        assert!(system.add_bounded_variable(valve).is_ok());
    }

    #[test]
    fn test_remove_bounded_variable() {
        use crate::inverse::{BoundedVariable, BoundedVariableId};

        let mut system = System::new();

        let id = BoundedVariableId::new("valve");
        let valve = BoundedVariable::new(id.clone(), 0.5, 0.0, 1.0).unwrap();

        system.add_bounded_variable(valve).unwrap();
        assert_eq!(system.bounded_variable_count(), 1);

        let removed = system.remove_bounded_variable(&id);
        assert!(removed.is_some());
        assert_eq!(system.bounded_variable_count(), 0);

        // Removing non-existent returns None
        let removed_again = system.remove_bounded_variable(&id);
        assert!(removed_again.is_none());
    }

    #[test]
    fn test_get_bounded_variable() {
        use crate::inverse::{BoundedVariable, BoundedVariableId};

        let mut system = System::new();

        let id = BoundedVariableId::new("vfd_speed");
        let vfd = BoundedVariable::new(id.clone(), 0.8, 0.3, 1.0).unwrap();

        system.add_bounded_variable(vfd).unwrap();

        let retrieved = system.get_bounded_variable(&id);
        assert!(retrieved.is_some());
        approx::assert_relative_eq!(retrieved.unwrap().value(), 0.8);

        let missing = system.get_bounded_variable(&BoundedVariableId::new("nonexistent"));
        assert!(missing.is_none());
    }

    #[test]
    fn test_get_bounded_variable_mut() {
        use crate::inverse::{BoundedVariable, BoundedVariableId};

        let mut system = System::new();

        let id = BoundedVariableId::new("valve");
        let valve = BoundedVariable::new(id.clone(), 0.5, 0.0, 1.0).unwrap();

        system.add_bounded_variable(valve).unwrap();

        // Apply step through mutable reference
        if let Some(v) = system.get_bounded_variable_mut(&id) {
            v.apply_step(0.3);
        }

        let retrieved = system.get_bounded_variable(&id).unwrap();
        approx::assert_relative_eq!(retrieved.value(), 0.8);
    }

    #[test]
    fn test_saturated_variables() {
        use crate::inverse::{BoundedVariable, BoundedVariableId, SaturationType};

        let mut system = System::new();

        // Variable at min bound
        let v1 = BoundedVariable::new(BoundedVariableId::new("v1"), 0.0, 0.0, 1.0).unwrap();
        // Variable at max bound
        let v2 = BoundedVariable::new(BoundedVariableId::new("v2"), 1.0, 0.0, 1.0).unwrap();
        // Variable in middle (not saturated)
        let v3 = BoundedVariable::new(BoundedVariableId::new("v3"), 0.5, 0.0, 1.0).unwrap();

        system.add_bounded_variable(v1).unwrap();
        system.add_bounded_variable(v2).unwrap();
        system.add_bounded_variable(v3).unwrap();

        let saturated = system.saturated_variables();
        assert_eq!(saturated.len(), 2);

        let sat_types: Vec<_> = saturated.iter().map(|s| s.saturation_type).collect();
        assert!(sat_types.contains(&SaturationType::LowerBound));
        assert!(sat_types.contains(&SaturationType::UpperBound));
    }

    #[test]
    fn test_bounded_variables_iterator() {
        use crate::inverse::{BoundedVariable, BoundedVariableId};

        let mut system = System::new();

        let v1 = BoundedVariable::new(BoundedVariableId::new("v1"), 0.5, 0.0, 1.0).unwrap();
        let v2 = BoundedVariable::new(BoundedVariableId::new("v2"), 0.3, 0.0, 1.0).unwrap();

        system.add_bounded_variable(v1).unwrap();
        system.add_bounded_variable(v2).unwrap();

        let count = system.bounded_variables().count();
        assert_eq!(count, 2);
    }

    // ────────────────────────────────────────────────────────────────────────
    // Inverse Control Tests (Story 5.3)
    // ────────────────────────────────────────────────────────────────────────

    #[test]
    fn test_link_constraint_to_control() {
        use crate::inverse::{
            BoundedVariable, BoundedVariableId, ComponentOutput, Constraint, ConstraintId,
        };

        let mut system = System::new();

        // Register component
        let node = system.add_component(make_mock(0));
        system.register_component_name("evaporator", node);

        // Add constraint and bounded variable
        let constraint = Constraint::new(
            ConstraintId::new("superheat"),
            ComponentOutput::Superheat {
                component_id: "evaporator".to_string(),
            },
            5.0,
        );
        system.add_constraint(constraint).unwrap();

        let valve = BoundedVariable::new(BoundedVariableId::new("valve"), 0.5, 0.0, 1.0).unwrap();
        system.add_bounded_variable(valve).unwrap();

        // Link them
        let result = system.link_constraint_to_control(
            &ConstraintId::new("superheat"),
            &BoundedVariableId::new("valve"),
        );
        assert!(result.is_ok());
        assert_eq!(system.inverse_control_mapping_count(), 1);

        // Verify bidirectional lookup
        let control = system.get_control_for_constraint(&ConstraintId::new("superheat"));
        assert!(control.is_some());
        assert_eq!(control.unwrap().as_str(), "valve");

        let constraint = system.get_constraint_for_control(&BoundedVariableId::new("valve"));
        assert!(constraint.is_some());
        assert_eq!(constraint.unwrap().as_str(), "superheat");
    }

    #[test]
    fn test_link_constraint_not_found() {
        use crate::inverse::{BoundedVariable, BoundedVariableId, ConstraintId, DoFError};

        let mut system = System::new();

        let valve = BoundedVariable::new(BoundedVariableId::new("valve"), 0.5, 0.0, 1.0).unwrap();
        system.add_bounded_variable(valve).unwrap();

        let result = system.link_constraint_to_control(
            &ConstraintId::new("nonexistent"),
            &BoundedVariableId::new("valve"),
        );
        assert!(matches!(result, Err(DoFError::ConstraintNotFound { .. })));
    }

    #[test]
    fn test_link_control_not_found() {
        use crate::inverse::{
            BoundedVariableId, ComponentOutput, Constraint, ConstraintId, DoFError,
        };

        let mut system = System::new();

        let node = system.add_component(make_mock(0));
        system.register_component_name("evaporator", node);

        let constraint = Constraint::new(
            ConstraintId::new("superheat"),
            ComponentOutput::Superheat {
                component_id: "evaporator".to_string(),
            },
            5.0,
        );
        system.add_constraint(constraint).unwrap();

        let result = system.link_constraint_to_control(
            &ConstraintId::new("superheat"),
            &BoundedVariableId::new("nonexistent"),
        );
        assert!(matches!(
            result,
            Err(DoFError::BoundedVariableNotFound { .. })
        ));
    }

    #[test]
    fn test_link_duplicate_constraint() {
        use crate::inverse::{
            BoundedVariable, BoundedVariableId, ComponentOutput, Constraint, ConstraintId, DoFError,
        };

        let mut system = System::new();

        let node = system.add_component(make_mock(0));
        system.register_component_name("evaporator", node);

        let constraint = Constraint::new(
            ConstraintId::new("superheat"),
            ComponentOutput::Superheat {
                component_id: "evaporator".to_string(),
            },
            5.0,
        );
        system.add_constraint(constraint).unwrap();

        let v1 = BoundedVariable::new(BoundedVariableId::new("v1"), 0.5, 0.0, 1.0).unwrap();
        let v2 = BoundedVariable::new(BoundedVariableId::new("v2"), 0.5, 0.0, 1.0).unwrap();
        system.add_bounded_variable(v1).unwrap();
        system.add_bounded_variable(v2).unwrap();

        // First link should succeed
        system
            .link_constraint_to_control(
                &ConstraintId::new("superheat"),
                &BoundedVariableId::new("v1"),
            )
            .unwrap();

        // Second link for same constraint should fail
        let result = system.link_constraint_to_control(
            &ConstraintId::new("superheat"),
            &BoundedVariableId::new("v2"),
        );
        assert!(matches!(result, Err(DoFError::AlreadyLinked { .. })));
    }

    #[test]
    fn test_link_duplicate_control() {
        use crate::inverse::{
            BoundedVariable, BoundedVariableId, ComponentOutput, Constraint, ConstraintId, DoFError,
        };

        let mut system = System::new();

        let node = system.add_component(make_mock(0));
        system.register_component_name("evaporator", node);

        let c1 = Constraint::new(
            ConstraintId::new("c1"),
            ComponentOutput::Superheat {
                component_id: "evaporator".to_string(),
            },
            5.0,
        );
        let c2 = Constraint::new(
            ConstraintId::new("c2"),
            ComponentOutput::Temperature {
                component_id: "evaporator".to_string(),
            },
            300.0,
        );
        system.add_constraint(c1).unwrap();
        system.add_constraint(c2).unwrap();

        let valve = BoundedVariable::new(BoundedVariableId::new("valve"), 0.5, 0.0, 1.0).unwrap();
        system.add_bounded_variable(valve).unwrap();

        // First link should succeed
        system
            .link_constraint_to_control(&ConstraintId::new("c1"), &BoundedVariableId::new("valve"))
            .unwrap();

        // Second link for same control should fail
        let result = system
            .link_constraint_to_control(&ConstraintId::new("c2"), &BoundedVariableId::new("valve"));
        assert!(matches!(result, Err(DoFError::ControlAlreadyLinked { .. })));
    }

    #[test]
    fn test_unlink_constraint() {
        use crate::inverse::{
            BoundedVariable, BoundedVariableId, ComponentOutput, Constraint, ConstraintId,
        };

        let mut system = System::new();

        let node = system.add_component(make_mock(0));
        system.register_component_name("evaporator", node);

        let constraint = Constraint::new(
            ConstraintId::new("superheat"),
            ComponentOutput::Superheat {
                component_id: "evaporator".to_string(),
            },
            5.0,
        );
        system.add_constraint(constraint).unwrap();

        let valve = BoundedVariable::new(BoundedVariableId::new("valve"), 0.5, 0.0, 1.0).unwrap();
        system.add_bounded_variable(valve).unwrap();

        system
            .link_constraint_to_control(
                &ConstraintId::new("superheat"),
                &BoundedVariableId::new("valve"),
            )
            .unwrap();
        assert_eq!(system.inverse_control_mapping_count(), 1);

        // Unlink
        let removed = system.unlink_constraint(&ConstraintId::new("superheat"));
        assert!(removed.is_some());
        assert_eq!(removed.unwrap().as_str(), "valve");
        assert_eq!(system.inverse_control_mapping_count(), 0);

        // Unlinking again returns None
        let removed_again = system.unlink_constraint(&ConstraintId::new("superheat"));
        assert!(removed_again.is_none());
    }

    #[test]
    fn test_control_variable_state_index() {
        use crate::inverse::{
            BoundedVariable, BoundedVariableId, ComponentOutput, Constraint, ConstraintId,
        };

        let mut system = System::new();

        // Add two components and an edge
        let n0 = system.add_component(make_mock(0));
        let n1 = system.add_component(make_mock(0));
        system.register_component_name("evaporator", n0);
        system.register_component_name("condenser", n1);
        system.add_edge(n0, n1).unwrap();
        system.finalize().unwrap();

        // edge_count = 1, so base index = 3 (3 unknowns per edge: ṁ, P, h)
        assert_eq!(system.edge_count(), 1);

        // Add constraint and bounded variable
        let constraint = Constraint::new(
            ConstraintId::new("superheat"),
            ComponentOutput::Superheat {
                component_id: "evaporator".to_string(),
            },
            5.0,
        );
        system.add_constraint(constraint).unwrap();

        let valve = BoundedVariable::new(BoundedVariableId::new("valve"), 0.5, 0.0, 1.0).unwrap();
        system.add_bounded_variable(valve).unwrap();

        // Link them
        system
            .link_constraint_to_control(
                &ConstraintId::new("superheat"),
                &BoundedVariableId::new("valve"),
            )
            .unwrap();

        // Control variable index should be at total_state_len = 3 * edge_count = 3
        let idx = system.control_variable_state_index(&BoundedVariableId::new("valve"));
        assert!(idx.is_some());
        assert_eq!(idx.unwrap(), 3);

        // Unlinked variable returns None
        let v2 = BoundedVariable::new(BoundedVariableId::new("v2"), 0.5, 0.0, 1.0).unwrap();
        system.add_bounded_variable(v2).unwrap();
        let idx2 = system.control_variable_state_index(&BoundedVariableId::new("v2"));
        assert!(idx2.is_none());
    }

    #[test]
    fn test_validate_inverse_control_dof_balanced() {
        use crate::inverse::{
            BoundedVariable, BoundedVariableId, ComponentOutput, Constraint, ConstraintId,
        };

        let mut system = System::new();

        // Add two components: one with 2 equations, one with 1 equation (total 3)
        // 1 edge (stride-3) = 3 unknowns → balanced without closures (CM1.3+)
        let n0 = system.add_component(make_mock(2));
        let n1 = system.add_component(make_mock(1));
        system.register_component_name("evaporator", n0);
        system.add_edge(n0, n1).unwrap();
        system.finalize().unwrap();

        // n_component_eqs = 3 (2 + 1)
        // n_edge_unknowns = 3 (1 edge × stride-3)

        // Add one constraint and one control
        let constraint = Constraint::new(
            ConstraintId::new("superheat"),
            ComponentOutput::Superheat {
                component_id: "evaporator".to_string(),
            },
            5.0,
        );
        system.add_constraint(constraint).unwrap();

        let valve = BoundedVariable::new(BoundedVariableId::new("valve"), 0.5, 0.0, 1.0).unwrap();
        system.add_bounded_variable(valve).unwrap();

        system
            .link_constraint_to_control(
                &ConstraintId::new("superheat"),
                &BoundedVariableId::new("valve"),
            )
            .unwrap();

        // n_equations = 3 (components) + 1 (constraint) = 4
        // n_unknowns = 3 (edge) + 1 (control) = 4 → Balanced!
        let result = system.validate_inverse_control_dof();
        assert!(result.is_ok());
    }

    #[test]
    fn test_validate_inverse_control_dof_over_constrained() {
        use crate::inverse::{
            BoundedVariable, BoundedVariableId, ComponentOutput, Constraint, ConstraintId, DoFError,
        };

        let mut system = System::new();

        // Add two components: 2 + 1 equations = 3 total; 1 edge = 3 unknowns
        let n0 = system.add_component(make_mock(2));
        let n1 = system.add_component(make_mock(1));
        system.register_component_name("evaporator", n0);
        system.register_component_name("condenser", n1);
        system.add_edge(n0, n1).unwrap();
        system.finalize().unwrap();

        // Add two constraints but only one control
        let c1 = Constraint::new(
            ConstraintId::new("c1"),
            ComponentOutput::Superheat {
                component_id: "evaporator".to_string(),
            },
            5.0,
        );
        let c2 = Constraint::new(
            ConstraintId::new("c2"),
            ComponentOutput::Temperature {
                component_id: "condenser".to_string(),
            },
            300.0,
        );
        system.add_constraint(c1).unwrap();
        system.add_constraint(c2).unwrap();

        let valve = BoundedVariable::new(BoundedVariableId::new("valve"), 0.5, 0.0, 1.0).unwrap();
        system.add_bounded_variable(valve).unwrap();

        system
            .link_constraint_to_control(&ConstraintId::new("c1"), &BoundedVariableId::new("valve"))
            .unwrap();

        // n_equations = 3 (components) + 2 (constraints) = 5
        // n_unknowns = 3 (1 edge × ṁ,P,h) + 1 (control) = 4
        // Over-constrained!
        let result = system.validate_inverse_control_dof();
        assert!(matches!(
            result,
            Err(DoFError::OverConstrainedSystem { .. })
        ));
        if let Err(DoFError::OverConstrainedSystem {
            constraint_count,
            control_count,
            equation_count,
            unknown_count,
        }) = result
        {
            assert_eq!(constraint_count, 2);
            assert_eq!(control_count, 1);
            assert_eq!(equation_count, 5);
            assert_eq!(unknown_count, 4);
        }
    }

    #[test]
    fn test_validate_inverse_control_dof_under_constrained() {
        use crate::inverse::{
            BoundedVariable, BoundedVariableId, ComponentOutput, Constraint, ConstraintId,
        };

        let mut system = System::new();

        // Add two components: 2 + 1 equations = 3 total; 1 edge = 3 unknowns
        let n0 = system.add_component(make_mock(2));
        let n1 = system.add_component(make_mock(1));
        system.register_component_name("evaporator", n0);
        system.add_edge(n0, n1).unwrap();
        system.finalize().unwrap();

        // Add one constraint and one control
        let constraint = Constraint::new(
            ConstraintId::new("superheat"),
            ComponentOutput::Superheat {
                component_id: "evaporator".to_string(),
            },
            5.0,
        );
        system.add_constraint(constraint).unwrap();

        // Add two bounded variables but only link one
        let v1 = BoundedVariable::new(BoundedVariableId::new("v1"), 0.5, 0.0, 1.0).unwrap();
        let v2 = BoundedVariable::new(BoundedVariableId::new("v2"), 0.5, 0.0, 1.0).unwrap();
        system.add_bounded_variable(v1).unwrap();
        system.add_bounded_variable(v2).unwrap();

        system
            .link_constraint_to_control(
                &ConstraintId::new("superheat"),
                &BoundedVariableId::new("v1"),
            )
            .unwrap();

        // n_equations = 3 (components) + 1 (constraint) = 4
        // n_unknowns = 3 (edge) + 1 (linked control) = 4
        // Balanced — unlinked bounded variables don't affect DoF
        let result = system.validate_inverse_control_dof();
        assert!(result.is_ok());
    }

    #[test]
    fn test_full_state_vector_len() {
        use crate::inverse::{
            BoundedVariable, BoundedVariableId, ComponentOutput, Constraint, ConstraintId,
        };

        let mut system = System::new();

        // Add two components and one edge
        let n0 = system.add_component(make_mock(0));
        let n1 = system.add_component(make_mock(0));
        system.register_component_name("evaporator", n0);
        system.add_edge(n0, n1).unwrap();
        system.finalize().unwrap();

        // Edge states: 3 * 1 = 3 (ṁ, P, h)
        assert_eq!(system.full_state_vector_len(), 3);

        // Add constraint and control
        let constraint = Constraint::new(
            ConstraintId::new("superheat"),
            ComponentOutput::Superheat {
                component_id: "evaporator".to_string(),
            },
            5.0,
        );
        system.add_constraint(constraint).unwrap();

        let valve = BoundedVariable::new(BoundedVariableId::new("valve"), 0.5, 0.0, 1.0).unwrap();
        system.add_bounded_variable(valve).unwrap();

        system
            .link_constraint_to_control(
                &ConstraintId::new("superheat"),
                &BoundedVariableId::new("valve"),
            )
            .unwrap();

        // Edge states: 3, control vars: 1, thermal couplings: 0
        // Total: 3 + 1 + 0 = 4
        assert_eq!(system.full_state_vector_len(), 4);
    }

    #[test]
    fn test_free_actuator_extends_state_and_indexes_after_saturated() {
        use crate::inverse::{BoundedVariable, BoundedVariableId};

        let mut system = System::new();
        let n0 = system.add_component(make_mock(0));
        let n1 = system.add_component(make_mock(0));
        system.register_component_name("exv", n0);
        system.add_edge(n0, n1).unwrap();

        // Edge states: 3 * 1 = 3 (ṁ, P, h). No control yet.
        // (length is only meaningful after finalize but the actuator arithmetic
        // does not depend on it.)
        let open = BoundedVariable::with_component(
            BoundedVariableId::new("exv__opening"),
            "exv",
            0.5,
            0.02,
            1.0,
        )
        .unwrap();
        system.add_bounded_variable(open).unwrap();
        system.add_free_actuator(BoundedVariableId::new("exv__opening"));

        system.finalize().unwrap();

        assert_eq!(system.free_actuator_count(), 1);
        // 3 edge unknowns + 1 free actuator = 4.
        assert_eq!(system.full_state_vector_len(), 4);
        // The free actuator lives right after the (empty) saturated block, i.e.
        // at total_state_len + 0 + 0 + 0 = 3.
        assert_eq!(system.free_actuator_index(0), 3);
    }

    #[test]
    #[should_panic(expected = "already registered")]
    fn test_free_actuator_rejects_duplicate() {
        use crate::inverse::{BoundedVariable, BoundedVariableId};

        let mut system = System::new();
        let n0 = system.add_component(make_mock(0));
        system.register_component_name("exv", n0);
        let open = BoundedVariable::with_component(
            BoundedVariableId::new("exv__opening"),
            "exv",
            0.5,
            0.0,
            1.0,
        )
        .unwrap();
        system.add_bounded_variable(open).unwrap();
        system.add_free_actuator(BoundedVariableId::new("exv__opening"));
        system.add_free_actuator(BoundedVariableId::new("exv__opening"));
    }

    #[test]
    fn test_saturated_controller_residuals_and_jacobian() {
        use crate::inverse::{
            BoundedVariable, BoundedVariableId, ComponentOutput, ConstraintId, SaturatedController,
            Saturation,
        };

        let mut system = System::new();
        let plant =
            system.add_component(Box::new(SaturatedControlMockComponent { z_flow_idx: None }));
        system.register_component_name("plant", plant);
        system.add_edge(plant, plant).unwrap();

        let actuator = BoundedVariable::with_component(
            BoundedVariableId::new("plant_f_m"),
            "plant",
            1.0,
            0.5,
            1.5,
        )
        .unwrap();
        system.add_bounded_variable(actuator).unwrap();
        let ctrl = SaturatedController::new(
            ConstraintId::new("capacity_loop"),
            ComponentOutput::Capacity {
                component_id: "plant".to_string(),
            },
            BoundedVariableId::new("plant_f_m"),
            15.0,
            0.5,
            1.5,
        )
        .unwrap()
        .with_gain(0.5)
        .unwrap()
        .with_band(1.0)
        .unwrap()
        .with_saturation(Saturation::Hard);
        system.add_saturated_controller(ctrl.clone());
        system.finalize().unwrap();

        assert_eq!(system.saturated_controller_count(), 1);
        assert_eq!(system.full_state_vector_len(), 5);
        assert_eq!(
            system.control_variable_state_index(&BoundedVariableId::new("plant_f_m")),
            Some(3)
        );

        let state = vec![1.0, 2.0, 3.0, 1.1, 0.2];
        let mut residuals = vec![0.0; 3];
        system.compute_residuals(&state, &mut residuals).unwrap();
        assert!((residuals[1] - ctrl.residual_u(state[3], state[4])).abs() < 1e-12);
        assert!((residuals[2] - ctrl.residual_y(14.0, state[4])).abs() < 1e-12);

        let mut jacobian = JacobianBuilder::new();
        system.assemble_jacobian(&state, &mut jacobian).unwrap();
        let entries = jacobian.entries();
        assert!(entries.contains(&(1, 3, 1.0)));
        assert!(entries.contains(&(1, 4, -0.5)));
        // Offset-free control law: ∂r_y/∂x = -(1 - S'(x)); x=0.2 is in-band so
        // S'(x)=1 and the entry is 0 (the row is closed via the plant coupling).
        assert!(entries.contains(&(2, 4, 0.0)));

        let mut analytic = vec![vec![0.0; state.len()]; residuals.len()];
        for &(row, col, value) in entries {
            analytic[row][col] += value;
        }

        let eps = 1e-6;
        let mut max_err = 0.0_f64;
        for col in 0..state.len() {
            let mut plus_state = state.clone();
            let mut minus_state = state.clone();
            plus_state[col] += eps;
            minus_state[col] -= eps;

            let mut plus = vec![0.0; residuals.len()];
            let mut minus = vec![0.0; residuals.len()];
            system.compute_residuals(&plus_state, &mut plus).unwrap();
            system.compute_residuals(&minus_state, &mut minus).unwrap();

            for row in 0..residuals.len() {
                let fd = (plus[row] - minus[row]) / (2.0 * eps);
                max_err = max_err.max((analytic[row][col] - fd).abs());
            }
        }
        println!("saturated controller jacobian max error: {max_err:.3e}");
        assert!(
            max_err <= 1e-4,
            "saturated controller analytic Jacobian mismatch: max_err={max_err}"
        );
    }

    #[test]
    fn test_saturated_network_controller_jacobian_matches_fd() {
        use crate::inverse::{
            BoundedVariable, BoundedVariableId, Combine, ComponentOutput, ConstraintId, Objective,
            SaturatedController, Saturation,
        };

        let mut system = System::new();
        let plant =
            system.add_component(Box::new(SaturatedControlMockComponent { z_flow_idx: None }));
        system.register_component_name("plant", plant);
        system.add_edge(plant, plant).unwrap();

        let actuator = BoundedVariable::with_component(
            BoundedVariableId::new("plant_f_m"),
            "plant",
            1.0,
            0.5,
            1.5,
        )
        .unwrap();
        system.add_bounded_variable(actuator).unwrap();

        // Override network: primary Capacity→15 seeds the fold; a Temperature→2
        // protection is folded in via softMax. A moderate alpha keeps the
        // selector well inside its smooth regime at the test point.
        let ctrl = SaturatedController::new(
            ConstraintId::new("net_loop"),
            ComponentOutput::Capacity {
                component_id: "plant".to_string(),
            },
            BoundedVariableId::new("plant_f_m"),
            15.0,
            0.5,
            1.5,
        )
        .unwrap()
        .with_gain(0.5)
        .unwrap()
        .with_band(1.0)
        .unwrap()
        .with_saturation(Saturation::Hard)
        .with_alpha(0.5)
        .with_objectives(vec![
            Objective::new(
                ComponentOutput::Capacity {
                    component_id: "plant".to_string(),
                },
                15.0,
                0.5,
                Combine::Min,
            ),
            Objective::new(
                ComponentOutput::Temperature {
                    component_id: "plant".to_string(),
                },
                2.0,
                -0.7,
                Combine::Max,
            ),
        ]);
        system.add_saturated_controller(ctrl.clone());
        system.finalize().unwrap();

        assert!(ctrl.is_network());

        let state = vec![1.0, 2.0, 3.0, 1.1, 0.2];
        let mut residuals = vec![0.0; 3];
        system.compute_residuals(&state, &mut residuals).unwrap();

        // Cross-check the network residual against a direct evaluation.
        let cap = state[0] + 2.0 * state[1] + 3.0 * state[2];
        let temp = state[1];
        let e = ctrl.error_signal(&[cap, temp]);
        assert!((residuals[2] - ctrl.residual_y_network(e, state[4])).abs() < 1e-12);

        let mut jacobian = JacobianBuilder::new();
        system.assemble_jacobian(&state, &mut jacobian).unwrap();
        let entries = jacobian.entries();
        let mut analytic = vec![vec![0.0; state.len()]; residuals.len()];
        for &(row, col, value) in entries {
            analytic[row][col] += value;
        }

        let eps = 1e-6;
        let mut max_err = 0.0_f64;
        for col in 0..state.len() {
            let mut plus_state = state.clone();
            let mut minus_state = state.clone();
            plus_state[col] += eps;
            minus_state[col] -= eps;
            let mut plus = vec![0.0; residuals.len()];
            let mut minus = vec![0.0; residuals.len()];
            system.compute_residuals(&plus_state, &mut plus).unwrap();
            system.compute_residuals(&minus_state, &mut minus).unwrap();
            for row in 0..residuals.len() {
                let fd = (plus[row] - minus[row]) / (2.0 * eps);
                max_err = max_err.max((analytic[row][col] - fd).abs());
            }
        }
        assert!(
            max_err <= 1e-4,
            "network controller analytic Jacobian mismatch: max_err={max_err}"
        );
    }

    #[test]
    fn test_control_variable_indices() {
        use crate::inverse::{
            BoundedVariable, BoundedVariableId, ComponentOutput, Constraint, ConstraintId,
        };

        let mut system = System::new();

        // Add two components and one edge
        let n0 = system.add_component(make_mock(0));
        let n1 = system.add_component(make_mock(0));
        system.register_component_name("evaporator", n0);
        system.register_component_name("condenser", n1);
        system.add_edge(n0, n1).unwrap();
        system.finalize().unwrap();

        // Add two constraints and controls
        let c1 = Constraint::new(
            ConstraintId::new("c1"),
            ComponentOutput::Superheat {
                component_id: "evaporator".to_string(),
            },
            5.0,
        );
        system.add_constraint(c1).unwrap();

        let v1 = BoundedVariable::new(BoundedVariableId::new("v1"), 0.5, 0.0, 1.0).unwrap();
        system.add_bounded_variable(v1).unwrap();

        system
            .link_constraint_to_control(&ConstraintId::new("c1"), &BoundedVariableId::new("v1"))
            .unwrap();

        let indices = system.control_variable_indices();
        assert_eq!(indices.len(), 1);
        assert_eq!(indices[0].1, 3); // total_state_len = 3 * edge_count = 3
    }

    struct BadMassFlowComponent {
        ports: Vec<ConnectedPort>,
    }

    impl Component for BadMassFlowComponent {
        fn compute_residuals(
            &self,
            _state: &StateSlice,
            _residuals: &mut entropyk_components::ResidualVector,
        ) -> Result<(), ComponentError> {
            Ok(())
        }

        fn jacobian_entries(
            &self,
            _state: &StateSlice,
            _jacobian: &mut JacobianBuilder,
        ) -> Result<(), ComponentError> {
            Ok(())
        }

        fn n_equations(&self) -> usize {
            0
        }

        fn get_ports(&self) -> &[ConnectedPort] {
            &self.ports
        }

        fn port_mass_flows(
            &self,
            _state: &StateSlice,
        ) -> Result<Vec<entropyk_core::MassFlow>, ComponentError> {
            Ok(vec![
                entropyk_core::MassFlow::from_kg_per_s(1.0),
                entropyk_core::MassFlow::from_kg_per_s(-0.5), // Intentionally unbalanced
            ])
        }
    }

    /// Component with balanced mass flow (inlet = outlet)
    struct BalancedMassFlowComponent {
        ports: Vec<ConnectedPort>,
    }

    impl Component for BalancedMassFlowComponent {
        fn compute_residuals(
            &self,
            _state: &StateSlice,
            _residuals: &mut entropyk_components::ResidualVector,
        ) -> Result<(), ComponentError> {
            Ok(())
        }

        fn jacobian_entries(
            &self,
            _state: &StateSlice,
            _jacobian: &mut JacobianBuilder,
        ) -> Result<(), ComponentError> {
            Ok(())
        }

        fn n_equations(&self) -> usize {
            0
        }

        fn get_ports(&self) -> &[ConnectedPort] {
            &self.ports
        }

        fn port_mass_flows(
            &self,
            _state: &StateSlice,
        ) -> Result<Vec<entropyk_core::MassFlow>, ComponentError> {
            // Balanced: inlet = 1.0 kg/s, outlet = -1.0 kg/s (sum = 0)
            Ok(vec![
                entropyk_core::MassFlow::from_kg_per_s(1.0),
                entropyk_core::MassFlow::from_kg_per_s(-1.0),
            ])
        }
    }

    #[test]
    fn test_mass_balance_passes_for_balanced_component() {
        let mut system = System::new();

        let inlet = Port::new(
            FluidId::new("R134a"),
            Pressure::from_bar(1.0),
            Enthalpy::from_joules_per_kg(400000.0),
        );
        let outlet = Port::new(
            FluidId::new("R134a"),
            Pressure::from_bar(1.0),
            Enthalpy::from_joules_per_kg(400000.0),
        );
        let (c1, c2) = inlet.connect(outlet).unwrap();

        let comp = Box::new(BalancedMassFlowComponent {
            ports: vec![c1, c2],
        });

        let n0 = system.add_component(comp);
        system.add_edge(n0, n0).unwrap(); // Self-edge to avoid isolated node

        system.finalize().unwrap();

        let state = vec![0.0; system.full_state_vector_len()];
        let result = system.check_mass_balance(&state);

        assert!(
            result.is_ok(),
            "Expected mass balance to pass for balanced component"
        );
    }

    #[test]
    fn test_mass_balance_violation() {
        let mut system = System::new();

        let inlet = Port::new(
            FluidId::new("R134a"),
            Pressure::from_bar(1.0),
            Enthalpy::from_joules_per_kg(400000.0),
        );
        let outlet = Port::new(
            FluidId::new("R134a"),
            Pressure::from_bar(1.0),
            Enthalpy::from_joules_per_kg(400000.0),
        );
        let (c1, c2) = inlet.connect(outlet).unwrap();

        let comp = Box::new(BadMassFlowComponent {
            ports: vec![c1, c2], // Just to have ports
        });

        let n0 = system.add_component(comp);
        system.add_edge(n0, n0).unwrap(); // Self-edge to avoid isolated node

        system.finalize().unwrap();

        // Ensure state is appropriately sized for finalize
        let state = vec![0.0; system.full_state_vector_len()];
        let result = system.check_mass_balance(&state);

        assert!(result.is_err());

        // Verify error contains mass error information
        if let Err(crate::SolverError::Validation {
            mass_error,
            energy_error,
        }) = result
        {
            assert!(mass_error > 0.0, "Mass error should be positive");
            assert_eq!(
                energy_error, 0.0,
                "Energy error should be zero for mass-only validation"
            );
        } else {
            panic!("Expected Validation error, got {:?}", result);
        }
    }

    #[test]
    fn test_mass_balance_tolerance_constant() {
        // Verify the tolerance constant is accessible and has expected value
        assert_eq!(System::MASS_BALANCE_TOLERANCE_KG_S, 1e-9);
    }

    #[test]
    fn test_generate_canonical_bytes() {
        let mut sys = System::new();
        let n0 = sys.add_component(make_mock(0));
        let n1 = sys.add_component(make_mock(0));
        sys.add_edge(n0, n1).unwrap();

        let bytes1 = sys.generate_canonical_bytes();
        let bytes2 = sys.generate_canonical_bytes();

        // Exact same graph state should produce same bytes
        assert_eq!(bytes1, bytes2);
    }

    #[test]
    fn test_input_hash_deterministic() {
        let mut sys1 = System::new();
        let n0_1 = sys1.add_component(make_mock(0));
        let n1_1 = sys1.add_component(make_mock(0));
        sys1.add_edge(n0_1, n1_1).unwrap();

        let mut sys2 = System::new();
        let n0_2 = sys2.add_component(make_mock(0));
        let n1_2 = sys2.add_component(make_mock(0));
        sys2.add_edge(n0_2, n1_2).unwrap();

        // Two identically constructed systems should have same hash
        assert_eq!(sys1.input_hash(), sys2.input_hash());

        // Now mutate one system by adding an edge
        sys1.add_edge(n1_1, n0_1).unwrap();

        // Hash should be different now
        assert_ne!(sys1.input_hash(), sys2.input_hash());
    }

    // ────────────────────────────────────────────────────────────────────────
    // Story 9.6: Energy Validation Logging Improvement Tests
    // ────────────────────────────────────────────────────────────────────────
    // Story 9.6: Energy Validation Logging Improvement Tests
    // ────────────────────────────────────────────────────────────────────────

    /// Test that check_energy_balance emits warnings for components without energy methods.
    /// This test verifies the logging improvement from Story 9.6.
    #[test]
    fn test_energy_balance_warns_for_skipped_components() {
        use tracing_subscriber::layer::SubscriberExt;
        use tracing_subscriber::util::SubscriberInitExt;

        // Create a system with mock components that don't implement energy_transfers()
        let mut sys = System::new();
        let n0 = sys.add_component(make_mock(0));
        let n1 = sys.add_component(make_mock(0));
        sys.add_edge(n0, n1).unwrap();
        sys.add_edge(n1, n0).unwrap();
        sys.finalize().unwrap();

        let state = vec![0.0; sys.state_vector_len()];

        // Capture log output using tracing_subscriber
        let log_buffer = std::sync::Arc::new(std::sync::Mutex::new(String::new()));
        let buffer_clone = log_buffer.clone();
        let layer = tracing_subscriber::fmt::layer()
            .with_writer(move || {
                use std::io::Write;
                struct BufWriter {
                    buf: std::sync::Arc<std::sync::Mutex<String>>,
                }
                impl Write for BufWriter {
                    fn write(&mut self, data: &[u8]) -> std::io::Result<usize> {
                        let mut buf = self.buf.lock().unwrap();
                        buf.push_str(&String::from_utf8_lossy(data));
                        Ok(data.len())
                    }
                    fn flush(&mut self) -> std::io::Result<()> {
                        Ok(())
                    }
                }
                BufWriter {
                    buf: buffer_clone.clone(),
                }
            })
            .without_time();

        let _guard = tracing_subscriber::registry().with(layer).set_default();

        // check_energy_balance should succeed (no violations) but will emit warnings
        // for components that lack energy_transfers() and port_enthalpies()
        let result = sys.check_energy_balance(&state);
        assert!(
            result.is_ok(),
            "check_energy_balance should succeed even with skipped components"
        );

        // Verify warning was emitted
        let log_output = log_buffer.lock().unwrap();
        assert!(
            log_output.contains("SKIPPED in energy balance validation"),
            "Expected warning message not found in logs. Actual output: {}",
            *log_output
        );
    }

    /// Test that cycle_performance classifies duties by physical sign and closes
    /// the First Law (Q_heating = Q_cooling + W_input).
    #[test]
    fn test_cycle_performance_classifies_and_closes_first_law() {
        // Compressor (work in), condenser (heat out), evaporator (heat in) in a loop.
        let mut sys = System::new();
        let comp = sys.add_component(Box::new(EnergyMock {
            heat_w: 0.0,
            work_w: -2000.0,
        }));
        let cond = sys.add_component(Box::new(EnergyMock {
            heat_w: -9000.0,
            work_w: 0.0,
        }));
        let evap = sys.add_component(Box::new(EnergyMock {
            heat_w: 7000.0,
            work_w: 0.0,
        }));
        sys.add_edge(comp, cond).unwrap();
        sys.add_edge(cond, evap).unwrap();
        sys.add_edge(evap, comp).unwrap();
        sys.finalize().unwrap();

        let state = vec![0.0; sys.state_vector_len()];
        let perf = sys
            .cycle_performance(&state)
            .expect("performance available");

        assert_relative_eq!(perf.q_cooling_w, 7000.0, epsilon = 1e-9);
        assert_relative_eq!(perf.q_heating_w, 9000.0, epsilon = 1e-9);
        assert_relative_eq!(perf.work_input_w, 2000.0, epsilon = 1e-9);
        assert_eq!(perf.components_counted, 3);
        // First Law: rejected = absorbed + work.
        assert_relative_eq!(
            perf.q_heating_w,
            perf.q_cooling_w + perf.work_input_w,
            epsilon = 1e-9
        );
        assert_relative_eq!(perf.cop_cooling().unwrap(), 3.5, epsilon = 1e-9);
        assert_relative_eq!(perf.cop_heating().unwrap(), 4.5, epsilon = 1e-9);
    }

    /// cycle_performance returns None when no component reports energy transfers.
    #[test]
    fn test_cycle_performance_none_without_energy_reporters() {
        let mut sys = System::new();
        let n0 = sys.add_component(make_mock(0));
        let n1 = sys.add_component(make_mock(0));
        sys.add_edge(n0, n1).unwrap();
        sys.add_edge(n1, n0).unwrap();
        sys.finalize().unwrap();

        let state = vec![0.0; sys.state_vector_len()];
        assert!(sys.cycle_performance(&state).is_none());
    }

    /// Test that check_energy_balance includes component type in warning message.
    #[test]
    fn test_energy_balance_includes_component_type_in_warning() {
        use tracing_subscriber::layer::SubscriberExt;
        use tracing_subscriber::util::SubscriberInitExt;

        // Create a system with mock components (need at least 2 nodes with edges to avoid isolated node error)
        let mut sys = System::new();
        let n0 = sys.add_component(make_mock(0));
        let n1 = sys.add_component(make_mock(0));
        sys.add_edge(n0, n1).unwrap();
        sys.add_edge(n1, n0).unwrap();
        sys.finalize().unwrap();

        let state = vec![0.0; sys.state_vector_len()];

        // Capture log output using tracing_subscriber
        let log_buffer = std::sync::Arc::new(std::sync::Mutex::new(String::new()));
        let buffer_clone = log_buffer.clone();
        let layer = tracing_subscriber::fmt::layer()
            .with_writer(move || {
                use std::io::Write;
                struct BufWriter {
                    buf: std::sync::Arc<std::sync::Mutex<String>>,
                }
                impl Write for BufWriter {
                    fn write(&mut self, data: &[u8]) -> std::io::Result<usize> {
                        let mut buf = self.buf.lock().unwrap();
                        buf.push_str(&String::from_utf8_lossy(data));
                        Ok(data.len())
                    }
                    fn flush(&mut self) -> std::io::Result<()> {
                        Ok(())
                    }
                }
                BufWriter {
                    buf: buffer_clone.clone(),
                }
            })
            .without_time();

        let _guard = tracing_subscriber::registry().with(layer).set_default();

        let result = sys.check_energy_balance(&state);
        assert!(result.is_ok());

        // Verify warning message includes component type information
        // Note: type_name_of_val on a trait object returns the trait name ("Component"),
        // not the concrete type. This is a known Rust limitation.
        let log_output = log_buffer.lock().unwrap();
        assert!(
            log_output.contains("type: Component"),
            "Expected component type information not found in logs. Actual output: {}",
            *log_output
        );
    }

    /// Test that check_energy_balance emits a summary warning with skipped component count.
    #[test]
    fn test_energy_balance_summary_warning() {
        use tracing_subscriber::layer::SubscriberExt;
        use tracing_subscriber::util::SubscriberInitExt;

        // Create a system with mock components
        let mut sys = System::new();
        let n0 = sys.add_component(make_mock(0));
        let n1 = sys.add_component(make_mock(0));
        sys.add_edge(n0, n1).unwrap();
        sys.add_edge(n1, n0).unwrap();
        sys.finalize().unwrap();

        let state = vec![0.0; sys.state_vector_len()];

        // Capture log output
        let log_buffer = std::sync::Arc::new(std::sync::Mutex::new(String::new()));
        let buffer_clone = log_buffer.clone();
        let layer = tracing_subscriber::fmt::layer()
            .with_writer(move || {
                use std::io::Write;
                struct BufWriter {
                    buf: std::sync::Arc<std::sync::Mutex<String>>,
                }
                impl Write for BufWriter {
                    fn write(&mut self, data: &[u8]) -> std::io::Result<usize> {
                        let mut buf = self.buf.lock().unwrap();
                        buf.push_str(&String::from_utf8_lossy(data));
                        Ok(data.len())
                    }
                    fn flush(&mut self) -> std::io::Result<()> {
                        Ok(())
                    }
                }
                BufWriter {
                    buf: buffer_clone.clone(),
                }
            })
            .without_time();

        let _guard = tracing_subscriber::registry().with(layer).set_default();

        let result = sys.check_energy_balance(&state);
        assert!(result.is_ok());

        // Verify summary warning was emitted
        let log_output = log_buffer.lock().unwrap();
        assert!(
            log_output.contains("Energy balance validation incomplete"),
            "Expected summary warning not found in logs. Actual output: {}",
            *log_output
        );
        assert!(
            log_output.contains("component(s) skipped"),
            "Expected 'component(s) skipped' not found in logs. Actual output: {}",
            *log_output
        );
    }
}