//! 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).
//!
//! 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::error::{AddEdgeError, TopologyError};
use crate::inverse::{
    BoundedVariable, BoundedVariableError, BoundedVariableId, Constraint, ConstraintError,
    ConstraintId, DoFError, InverseControlConfig,
};
use entropyk_core::{CircuitId, Temperature};

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

/// 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 pressure (P) and enthalpy (h) at that connection.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub struct FlowEdge {
    /// State vector index for pressure (Pa)
    pub state_index_p: usize,
    /// State vector index for enthalpy (J/kg)
    pub state_index_h: usize,
}

/// System graph structure.
///
/// Nodes are components (`Box<dyn Component>`), edges are flow connections with
/// state indices. The state vector layout is:
///
/// ```text
/// [P_edge0, h_edge0, 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_p, state_index_h) - built in finalize()
    edge_to_state: HashMap<EdgeIndex, (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,
    /// Registry of component names for constraint validation.
    /// Maps human-readable names (e.g., "evaporator") to NodeIndex.
    component_names: HashMap<String, NodeIndex>,
    finalized: bool,
    total_state_len: usize,
}

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(),
            component_names: HashMap::new(),
            finalized: false,
            total_state_len: 0,
        }
    }

    /// 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 {
                state_index_p: 0,
                state_index_h: 0,
            },
        ))
    }

    /// 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)
            self.finalized = false;
            return Ok(self.graph.add_edge(
                source,
                target,
                FlowEdge {
                    state_index_p: 0,
                    state_index_h: 0,
                },
            ));
        }

        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 {
                state_index_p: 0,
                state_index_h: 0,
            },
        ))
    }

    /// Finalizes the graph: builds edge→state index mapping and validates topology.
    ///
    /// # State vector layout
    ///
    /// The state vector has length `2 * edge_count`. For each edge (in graph iteration order):
    /// - `state[2*i]` = pressure at edge i (Pa)
    /// - `state[2*i + 1]` = enthalpy at edge i (J/kg)
    ///
    /// # 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");
        }

        self.edge_to_state.clear();
        let mut idx = 0usize;
        for edge_idx in self.graph.edge_indices() {
            let (p_idx, h_idx) = (idx, idx + 1);
            self.edge_to_state.insert(edge_idx, (p_idx, h_idx));
            if let Some(weight) = self.graph.edge_weight_mut(edge_idx) {
                weight.state_index_p = p_idx;
                weight.state_index_h = h_idx;
            }
            idx += 2;
        }
        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
        // 2 * edge_count entries 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 `2 * edge_count` as
        // the base offset, and accumulate `internal_state_len` for each component.
        let mut current_offset = 2 * self.graph.edge_count();

        // Gather (node_idx, offset, incident_edge_indices) before mutating nodes.
        #[allow(clippy::type_complexity)]
        let mut node_context: Vec<(
            petgraph::graph::NodeIndex,
            usize,
            Vec<(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)> = Vec::new();
            for edge_ref in self
                .graph
                .edges_directed(node_idx, petgraph::Direction::Incoming)
            {
                if let Some(&ph) = self.edge_to_state.get(&edge_ref.id()) {
                    incident.push(ph);
                }
            }
            for edge_ref in self
                .graph
                .edges_directed(node_idx, petgraph::Direction::Outgoing)
            {
                if let Some(&ph) = self.edge_to_state.get(&edge_ref.id()) {
                    incident.push(ph);
                }
            }
            node_context.push((node_idx, current_offset, incident));

            // 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 id_str.ends_with("f_m") || id_str == "f_m" {
                        indices.f_m = Some(state_idx);
                    } else if id_str.ends_with("f_dp") || id_str == "f_dp" {
                        indices.f_dp = Some(state_idx);
                    } else if id_str.ends_with("f_ua") || id_str == "f_ua" {
                        indices.f_ua = Some(state_idx);
                    } else if id_str.ends_with("f_power") || id_str == "f_power" {
                        indices.f_power = Some(state_idx);
                    } else if id_str.ends_with("f_etav") || id_str == "f_etav" {
                        indices.f_etav = Some(state_idx);
                    }
                }
            }
        }

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

                // 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);
                    }
                }
            }
        }

        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);
                }
            }
        }

        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()");
        "[P_edge0, h_edge0, P_edge1, h_edge1, ...] — 2 per edge (pressure Pa, enthalpy J/kg)"
    }

    /// Returns the length of the state vector: `2 * edge_count + internal_components_length`.
    ///
    /// 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.
    ///
    /// # 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()"
        );
        *self
            .edge_to_state
            .get(&edge_id)
            .expect("invalid edge index")
    }

    /// 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 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())
    }

    // ────────────────────────────────────────────────────────────────────────
    // 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.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
    ///
    /// The number of constraint residuals added.
    ///
    /// # 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>,
    ) -> usize {
        if self.constraints.is_empty() {
            return 0;
        }

        let mut count = 0;
        for constraint in self.constraints.values() {
            let measured = measured_values
                .get(constraint.id())
                .copied()
                .unwrap_or_else(|| {
                    tracing::warn!(
                        constraint_id = constraint.id().as_str(),
                        "No measured value for constraint, using zero residual"
                    );
                    constraint.target_value()
                });
            let residual = constraint.compute_residual(measured);
            if count < residuals.len() {
                residuals[count] = residual;
            }
            count += 1;
        }
        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.
    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) {
                // 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(&(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;
                            }
                        }

                        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;
        }

        // Use configurable epsilon from InverseControlConfig
        let eps = self.inverse_control.finite_diff_epsilon();
        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
                    for edge in self
                        .graph
                        .edges_directed(node_idx, petgraph::Direction::Incoming)
                    {
                        if let Some(&(p_idx, h_idx)) = self.edge_to_state.get(&edge.id()) {
                            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(&(p_idx, h_idx)) = self.edge_to_state.get(&edge.id()) {
                            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
    }

    // ────────────────────────────────────────────────────────────────────────
    // 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.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
            .values()
            .filter_map(|v| v.is_saturated())
            .collect()
    }

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

    /// 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.
    ///
    /// # Panics
    ///
    /// Panics if epsilon is non-positive.
    pub fn set_inverse_control_epsilon(&mut self, epsilon: f64) {
        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(())
        }
    }

    /// 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]
    ///                [P0, h0, P1, h1, ... | ctrl0, ctrl1, ... | T_hot0, T_cold0, ...]
    /// ```
    ///
    /// The index for control variable `i` is: `2 * edge_count + 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> {
        if !self.inverse_control.is_control_linked(id) {
            return None;
        }

        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);
            }
        }
        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))
    }

    /// 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 the total state vector length including control variables.
    ///
    /// ```text
    /// full_state_len = 2 * edge_count + control_variable_count + 2 * thermal_coupling_count
    /// ```
    pub fn full_state_vector_len(&self) -> usize {
        self.total_state_len
            + self.inverse_control.mapping_count()
            + 2 * self.thermal_couplings.len()
    }

    /// 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(&(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(&(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();

        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);
        eq_offset += n_constraints;

        // Add couplings
        let n_couplings = self.coupling_residual_count();
        if n_couplings > 0 {
            // MVP: We need real temperatures in K from components, using dummy 300K for now
            let temps = vec![(300.0, 300.0); n_couplings];
            self.coupling_residuals(&temps, &mut residuals[eq_offset..eq_offset + n_couplings]);
        }
        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);
        }
        row_offset += self.constraints.len();

        // Add couplings jacobian (MVP: missing T indices)
        let _ = row_offset; // avoid unused warning
        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(())
    }

    /// 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");
        for coupling in &self.thermal_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())
    }
}

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 })
    }

    /// 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();
        assert_eq!(sys.state_vector_len(), 4); // 2 edges * 2 = 4
    }

    #[test]
    fn test_edge_indices_contiguous() {
        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 mut indices: Vec<usize> = Vec::new();
        for edge_idx in sys.edge_indices() {
            let (p, h) = sys.edge_state_indices(edge_idx);
            indices.push(p);
            indices.push(h);
        }

        let n = sys.edge_count();
        assert_eq!(indices.len(), 2 * n);
        let expected: Vec<usize> = (0..2 * n).collect();
        assert_eq!(indices, expected, "indices should be 0..2n");
    }

    #[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"));

        let state_len = sys.state_vector_len();
        assert_eq!(state_len, 4);

        let mut state = vec![0.0; state_len];
        state[0] = 1e5; // P_edge0
        state[1] = 250000.0; // h_edge0
        state[2] = 5e5; // P_edge1
        state[3] = 300000.0; // h_edge1

        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
        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;
        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 = 2
        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 2 * edge_count = 2
        let idx = system.control_variable_state_index(&BoundedVariableId::new("valve"));
        assert!(idx.is_some());
        assert_eq!(idx.unwrap(), 2);

        // 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 with 1 equation each and one edge
        let n0 = system.add_component(make_mock(1));
        let n1 = system.add_component(make_mock(1));
        system.register_component_name("evaporator", n0);
        system.add_edge(n0, n1).unwrap();
        system.finalize().unwrap();

        // n_edge_eqs = 2 (each mock component has 1 equation)
        // n_edge_unknowns = 2 * 1 = 2
        // Without constraints: 2 eqs = 2 unknowns (balanced)

        // 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 = 2 + 1 = 3
        // n_unknowns = 2 + 1 = 3
        // 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 with 1 equation each and one edge
        let n0 = system.add_component(make_mock(1));
        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 = 2 + 2 = 4
        // n_unknowns = 2 + 1 = 3
        // 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, 4);
            assert_eq!(unknown_count, 3);
        }
    }

    #[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 with 1 equation each and one edge
        let n0 = system.add_component(make_mock(1));
        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 = 2 + 1 = 3
        // n_unknowns = 2 + 1 = 3 (only linked controls count)
        // 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: 2 * 1 = 2
        assert_eq!(system.full_state_vector_len(), 2);

        // 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: 2, control vars: 1, thermal couplings: 0
        // Total: 2 + 1 + 0 = 3
        assert_eq!(system.full_state_vector_len(), 3);
    }

    #[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, 2); // 2 * edge_count = 2
    }

    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 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
        );
    }
}