oxigrid 0.1.2

//! Distribution network reconfiguration — minimum-loss switching.
//!
//! Distribution networks are typically operated radially (tree topology) for
//! protection simplicity, but are meshed in design with normally-open tie
//! switches.  Reconfiguration finds the optimal set of switch states that
//! minimises I²R losses while preserving radiality (connected spanning tree).
//!
//! # Algorithms
//!
//! * [`ReconfigAlgorithm::BranchExchange`] — Merlin-Back heuristic: iteratively
//!   close one tie switch, open the highest-loss branch in the created loop.
//! * [`ReconfigAlgorithm::SimulatedAnnealing`] — stochastic global search with
//!   geometric cooling; random LCG used (no `rand` crate).
//! * [`ReconfigAlgorithm::TabuSearch`] — deterministic neighbourhood search with
//!   a tabu list of the last 7 switch operations.
//! * [`ReconfigAlgorithm::Exhaustive`] — enumerate all radial spanning trees
//!   (only feasible when the number of tie switches is ≤ 10).
//!
//! # References
//! - Merlin & Back, "Search for a Minimum-Loss Operating Spanning Tree
//!   Configuration in an Urban Power Distribution System", PSCC 1975.
//! - Shirmohammadi & Hong, "Reconfiguration of Electric Distribution Networks
//!   for Resistive Line Losses Reduction", IEEE TPWRD, 1989.

use crate::error::{OxiGridError, Result};
use crate::network::topology::PowerNetwork;
use serde::{Deserialize, Serialize};
use std::collections::{HashSet, VecDeque};

// ─────────────────────────────────────────────────────────────────────────────
// Public data structures
// ─────────────────────────────────────────────────────────────────────────────

/// A switch in the distribution network.
///
/// Each switch is associated with a specific branch.  In the base (normal)
/// configuration, tie switches are open and section switches are closed.
#[derive(Debug, Clone, Serialize, Deserialize)]
pub struct Switch {
    /// Unique switch identifier.
    pub id: usize,
    /// Index into `PowerNetwork::branches` of the branch this switch controls.
    pub branch_idx: usize,
    /// `true` if this is a tie switch (normally open in base configuration).
    pub normally_open: bool,
    /// `false` means the switch is locked out (cannot be operated).
    pub can_operate: bool,
}

/// The state of all switches at a particular instant (a snapshot of the
/// network topology).
#[derive(Debug, Clone, Serialize, Deserialize)]
pub struct NetworkConfiguration {
    /// Switch IDs that are currently open (branch is de-energised).
    pub open_switches: Vec<usize>,
    /// Switch IDs that are currently closed (branch is energised).
    pub closed_switches: Vec<usize>,
    /// Whether the energised sub-graph forms a tree (no cycles).
    pub is_radial: bool,
    /// Number of isolated islands (connected components minus the main).
    pub n_islands: usize,
}

/// Algorithm selection for the reconfiguration solver.
#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)]
pub enum ReconfigAlgorithm {
    /// Merlin-Back branch-exchange heuristic (fast, locally optimal).
    BranchExchange,
    /// Simulated annealing — stochastic global search.
    SimulatedAnnealing,
    /// Tabu search — deterministic local search escaping local minima.
    TabuSearch,
    /// Exhaustive enumeration of radial spanning trees.
    ///
    /// Only practical when the number of tie switches is ≤ 10.
    Exhaustive,
}

/// Configuration parameters for the reconfiguration solver.
#[derive(Debug, Clone, Serialize, Deserialize)]
pub struct ReconfigConfig {
    /// All switches (tie and section) in the network.
    pub switches: Vec<Switch>,
    /// Maximum number of outer iterations / SA steps.
    pub max_iterations: usize,
    /// Algorithm to use.
    pub algorithm: ReconfigAlgorithm,
    /// Minimum acceptable bus voltage [p.u.].
    pub min_voltage_pu: f64,
    /// Maximum branch loading percentage.
    pub max_loading_pct: f64,
    /// Objective weight on total losses (MW).
    pub loss_weight: f64,
    /// Objective weight on maximum voltage deviation from 1.0 p.u.
    pub voltage_weight: f64,
}

impl Default for ReconfigConfig {
    fn default() -> Self {
        Self {
            switches: Vec::new(),
            max_iterations: 100,
            algorithm: ReconfigAlgorithm::BranchExchange,
            min_voltage_pu: 0.95,
            max_loading_pct: 100.0,
            loss_weight: 1.0,
            voltage_weight: 0.1,
        }
    }
}

/// A single switch operation (open or close) applied to a switch.
#[derive(Debug, Clone, Serialize, Deserialize)]
pub struct SwitchOperation {
    /// The switch that was operated.
    pub switch_id: usize,
    /// The action taken.
    pub action: SwitchAction,
}

/// The action taken on a switch.
#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)]
pub enum SwitchAction {
    Open,
    Close,
}

/// Result of the reconfiguration optimisation.
#[derive(Debug, Clone, Serialize, Deserialize)]
pub struct ReconfigResult {
    /// The optimal switch configuration found.
    pub best_configuration: NetworkConfiguration,
    /// Total I²R losses in the base (initial) configuration \[MW\].
    pub base_losses_mw: f64,
    /// Total I²R losses in the optimised configuration \[MW\].
    pub optimized_losses_mw: f64,
    /// Percentage reduction in losses: (base − opt) / base × 100.
    pub loss_reduction_pct: f64,
    /// List of switch operations that transform the base config to the best.
    pub switch_operations: Vec<SwitchOperation>,
    /// Number of iterations (or SA steps) actually performed.
    pub n_iterations: usize,
    /// Whether the algorithm converged within `max_iterations`.
    pub converged: bool,
    /// Minimum voltage seen at each bus across the optimised solution [p.u.].
    pub voltage_profile: Vec<f64>,
}

// ─────────────────────────────────────────────────────────────────────────────
// Main solver struct
// ─────────────────────────────────────────────────────────────────────────────

/// Distribution network reconfiguration solver.
pub struct NetworkReconfiguration {
    /// Solver configuration including switch list and algorithm choice.
    pub config: ReconfigConfig,
}

impl NetworkReconfiguration {
    /// Construct a solver with the supplied configuration.
    pub fn new(config: ReconfigConfig) -> Self {
        Self { config }
    }

    /// Solve the reconfiguration problem for `network`.
    ///
    /// Returns an error if the base configuration is infeasible (non-radial or
    /// disconnected).
    pub fn solve(&self, network: &PowerNetwork) -> Result<ReconfigResult> {
        // Build the base configuration from switch data
        let base_config = build_base_configuration(&self.config.switches, network)?;

        // Evaluate base losses
        let base_losses = self.evaluate_config(network, &base_config)?;

        let n_buses = network.bus_count();

        let (best_config, best_losses, n_iters, converged) = match self.config.algorithm {
            ReconfigAlgorithm::BranchExchange => {
                branch_exchange_solve(network, &base_config, &self.config.switches, &self.config)?
            }
            ReconfigAlgorithm::SimulatedAnnealing => {
                let (cfg, losses) = simulated_annealing(
                    network,
                    &base_config,
                    &self.config.switches,
                    &self.config,
                    42_u64,
                )?;
                (cfg, losses, self.config.max_iterations, true)
            }
            ReconfigAlgorithm::TabuSearch => {
                let (cfg, losses) =
                    tabu_search(network, &base_config, &self.config.switches, &self.config)?;
                (cfg, losses, self.config.max_iterations, true)
            }
            ReconfigAlgorithm::Exhaustive => {
                exhaustive_solve(network, &base_config, &self.config.switches, &self.config)?
            }
        };

        let loss_reduction_pct = if base_losses > 1e-12 {
            (base_losses - best_losses) / base_losses * 100.0
        } else {
            0.0
        };

        let switch_ops =
            compute_switch_operations(&base_config, &best_config, &self.config.switches);

        // Voltage profile: simple flat 1.0 estimate per bus (full power flow
        // integration would require running NR which is expensive here).
        let voltage_profile = vec![1.0_f64; n_buses];

        Ok(ReconfigResult {
            best_configuration: best_config,
            base_losses_mw: base_losses,
            optimized_losses_mw: best_losses,
            loss_reduction_pct,
            switch_operations: switch_ops,
            n_iterations: n_iters,
            converged,
            voltage_profile,
        })
    }

    /// Enumerate all valid radial configurations (for exhaustive search).
    ///
    /// A radial configuration is obtained by choosing exactly one switch to
    /// open from each fundamental loop created by the tie switches.
    pub fn enumerate_radial_configs(
        &self,
        network: &PowerNetwork,
    ) -> Result<Vec<NetworkConfiguration>> {
        let tie_switches: Vec<&Switch> = self
            .config
            .switches
            .iter()
            .filter(|s| s.normally_open)
            .collect();

        if tie_switches.len() > 10 {
            return Err(OxiGridError::InvalidParameter(
                "Exhaustive enumeration only supported for <= 10 tie switches".to_string(),
            ));
        }

        let n_tie = tie_switches.len();
        let n_configs = 1usize << n_tie; // 2^n_tie
        let mut configs = Vec::with_capacity(n_configs);

        for mask in 0..n_configs {
            // Each bit in `mask` indicates whether tie switch i is closed (1) or open (0)
            let open_sw: Vec<usize> = self
                .config
                .switches
                .iter()
                .filter_map(|s| {
                    if s.normally_open {
                        // index in tie_switches list
                        let tie_idx = tie_switches.iter().position(|t| t.id == s.id)?;
                        if (mask >> tie_idx) & 1 == 0 {
                            Some(s.id)
                        } else {
                            None
                        }
                    } else {
                        // section switch — normally closed; include in open set if bit not set
                        // For this enumeration we only toggle tie switches
                        None
                    }
                })
                .collect();

            let closed_sw: Vec<usize> = self
                .config
                .switches
                .iter()
                .filter_map(|s| {
                    if open_sw.contains(&s.id) {
                        None
                    } else {
                        Some(s.id)
                    }
                })
                .collect();

            let cfg = NetworkConfiguration {
                open_switches: open_sw,
                closed_switches: closed_sw,
                is_radial: false, // will be validated below
                n_islands: 0,
            };

            let valid = verify_radiality(network, &cfg, &self.config.switches);
            if valid {
                configs.push(NetworkConfiguration {
                    is_radial: true,
                    n_islands: 0,
                    ..cfg
                });
            }
        }

        Ok(configs)
    }

    /// Evaluate the I²R losses \[MW\] for a specific switch configuration.
    ///
    /// Uses a simple DC-style current estimate based on branch resistance
    /// and approximate branch currents derived from the network admittance
    /// structure and load data.
    pub fn evaluate_config(
        &self,
        network: &PowerNetwork,
        config: &NetworkConfiguration,
    ) -> Result<f64> {
        let open_branch_indices = config_to_open_branch_indices(config, &self.config.switches)?;
        compute_losses_estimate(network, &open_branch_indices)
    }
}

// ─────────────────────────────────────────────────────────────────────────────
// Loss sensitivity
// ─────────────────────────────────────────────────────────────────────────────

/// Compute loss sensitivity for each branch: ∂P_loss / ∂I_branch ≈ 2·R·I.
///
/// A higher value indicates that reducing flow on that branch has greater
/// impact on total losses.  Used to rank candidate branches for exchange.
///
/// # Arguments
/// * `network` — The power network.
/// * `v_mag`   — Bus voltage magnitudes [p.u.], length = n_buses.
/// * `v_ang`   — Bus voltage angles \[rad\], length = n_buses.
///
/// # Returns
/// Vector of length n_branches with ∂P_loss [MW / p.u. current].
pub fn compute_branch_loss_sensitivity(
    network: &PowerNetwork,
    v_mag: &[f64],
    v_ang: &[f64],
) -> Vec<f64> {
    let base_mva = network.base_mva;
    network
        .branches
        .iter()
        .map(|br| {
            if !br.status {
                return 0.0;
            }
            let fi = match network.bus_index(br.from_bus) {
                Ok(i) => i,
                Err(_) => return 0.0,
            };
            let ti = match network.bus_index(br.to_bus) {
                Ok(i) => i,
                Err(_) => return 0.0,
            };
            let vm_f = v_mag.get(fi).copied().unwrap_or(1.0);
            let vm_t = v_mag.get(ti).copied().unwrap_or(1.0);
            let va_f = v_ang.get(fi).copied().unwrap_or(0.0);
            let va_t = v_ang.get(ti).copied().unwrap_or(0.0);

            // Complex voltage difference
            let dv_re = vm_f * va_f.cos() - vm_t * va_t.cos();
            let dv_im = vm_f * va_f.sin() - vm_t * va_t.sin();

            // Series impedance magnitude squared
            let z_sq = br.r * br.r + br.x * br.x;
            if z_sq < 1e-15 {
                return 0.0;
            }

            // Approximate branch current magnitude [p.u.]
            let i_pu = (dv_re * dv_re + dv_im * dv_im).sqrt() / z_sq.sqrt();

            // ∂P_loss/∂I = 2·R·I, converted to MW
            2.0 * br.r * i_pu * base_mva
        })
        .collect()
}

// ─────────────────────────────────────────────────────────────────────────────
// Branch exchange algorithm (Merlin-Back)
// ─────────────────────────────────────────────────────────────────────────────

/// Run the Merlin-Back branch-exchange heuristic to completion.
fn branch_exchange_solve(
    network: &PowerNetwork,
    base_config: &NetworkConfiguration,
    switches: &[Switch],
    config: &ReconfigConfig,
) -> Result<(NetworkConfiguration, f64, usize, bool)> {
    let mut current_config = base_config.clone();
    let mut current_losses = compute_losses_estimate(
        network,
        &config_to_open_branch_indices(&current_config, switches)?,
    )?;

    let mut n_iters = 0;
    let mut improved = true;

    while improved && n_iters < config.max_iterations {
        improved = false;
        n_iters += 1;

        if let Some((new_config, new_losses)) =
            branch_exchange_step(network, &current_config, switches)?
        {
            if new_losses < current_losses - 1e-9 {
                current_config = new_config;
                current_losses = new_losses;
                improved = true;
            }
        }
    }

    let converged = n_iters < config.max_iterations;
    Ok((current_config, current_losses, n_iters, converged))
}

/// Perform one iteration of the branch-exchange heuristic.
///
/// For each tie switch (normally open):
///   1. Temporarily close it → creates one loop in the graph.
///   2. Identify all branches that form that loop.
///   3. Find the branch in the loop with the highest loss sensitivity.
///   4. If opening that branch (and closing the tie) reduces total losses
///      and maintains radiality, apply the exchange and return.
///
/// Returns `None` if no improving exchange is found.
fn branch_exchange_step(
    network: &PowerNetwork,
    config: &NetworkConfiguration,
    switches: &[Switch],
) -> Result<Option<(NetworkConfiguration, f64)>> {
    let n_buses = network.bus_count();
    let flat_v_mag = vec![1.0_f64; n_buses];
    let flat_v_ang = vec![0.0_f64; n_buses];
    let sensitivity = compute_branch_loss_sensitivity(network, &flat_v_mag, &flat_v_ang);

    let open_branch_indices = config_to_open_branch_indices(config, switches)?;
    let current_losses = compute_losses_estimate(network, &open_branch_indices)?;

    let mut best_improvement: Option<(NetworkConfiguration, f64)> = None;

    // Collect tie switches that are currently open and can be operated
    let tie_sw_ids: Vec<usize> = switches
        .iter()
        .filter(|s| s.normally_open && s.can_operate && config.open_switches.contains(&s.id))
        .map(|s| s.id)
        .collect();

    for &tie_id in &tie_sw_ids {
        let tie_sw = match switches.iter().find(|s| s.id == tie_id) {
            Some(s) => s,
            None => continue,
        };

        // Branches forming the loop when we close this tie switch
        let loop_branches =
            find_loop_branches_closed_tie(network, config, switches, tie_sw.branch_idx);

        if loop_branches.is_empty() {
            continue;
        }

        // Rank loop branches by loss sensitivity (descending) — best to open first
        let mut ranked = loop_branches.clone();
        ranked.sort_by(|&a, &b| {
            sensitivity
                .get(b)
                .copied()
                .unwrap_or(0.0)
                .partial_cmp(&sensitivity.get(a).copied().unwrap_or(0.0))
                .unwrap_or(std::cmp::Ordering::Equal)
        });

        for &branch_to_open in &ranked {
            // Find the section switch controlling this branch
            let sec_sw = switches
                .iter()
                .find(|s| s.branch_idx == branch_to_open && !s.normally_open && s.can_operate);

            let sec_sw_id = match sec_sw {
                Some(s) => s.id,
                None => {
                    // No operable section switch on this branch; try next
                    continue;
                }
            };

            // Skip if the section switch is already open (would disconnect)
            if config.open_switches.contains(&sec_sw_id) {
                continue;
            }

            // Build trial configuration: close tie, open section
            let mut new_open = config.open_switches.clone();
            new_open.retain(|&id| id != tie_id);
            new_open.push(sec_sw_id);
            new_open.sort_unstable();

            let mut new_closed = config.closed_switches.clone();
            new_closed.retain(|&id| id != sec_sw_id);
            new_closed.push(tie_id);
            new_closed.sort_unstable();

            let trial_config = NetworkConfiguration {
                open_switches: new_open,
                closed_switches: new_closed,
                is_radial: false,
                n_islands: 0,
            };

            // Verify radiality
            if !verify_radiality(network, &trial_config, switches) {
                continue;
            }

            let trial_open_branches = config_to_open_branch_indices(&trial_config, switches)?;
            let trial_losses = compute_losses_estimate(network, &trial_open_branches)?;

            if trial_losses < current_losses - 1e-9 {
                let improvement = current_losses - trial_losses;
                let is_better = best_improvement
                    .as_ref()
                    .map(|(_, prev_loss)| trial_losses < *prev_loss)
                    .unwrap_or(true);

                if is_better || improvement > 0.0 {
                    best_improvement = Some((
                        NetworkConfiguration {
                            is_radial: true,
                            n_islands: 0,
                            ..trial_config
                        },
                        trial_losses,
                    ));
                }
            }
            // Take the first improving exchange and break (Merlin-Back style)
            break;
        }
    }

    Ok(best_improvement)
}

/// Find the set of branch indices that form the unique loop created by
/// closing tie switch at `tie_branch_idx` in the current configuration.
///
/// Algorithm: BFS from `from_bus` of the tie branch, traversing only
/// currently-energised branches, until we reach `to_bus`.  The path
/// returned is the loop body.
fn find_loop_branches_closed_tie(
    network: &PowerNetwork,
    config: &NetworkConfiguration,
    switches: &[Switch],
    tie_branch_idx: usize,
) -> Vec<usize> {
    let open_branch_set: HashSet<usize> = config_to_open_branch_indices_unchecked(config, switches);

    let tie_branch = match network.branches.get(tie_branch_idx) {
        Some(b) => b,
        None => return vec![],
    };

    let src_bus = tie_branch.from_bus;
    let dst_bus = tie_branch.to_bus;

    // Build adjacency: branch_idx → (bus_a, bus_b) for energised branches
    // (excluding the tie branch itself, which is currently open)
    let adj = build_adjacency_excluding(network, &open_branch_set, Some(tie_branch_idx));

    // BFS to find path from src_bus to dst_bus
    bfs_path(&adj, src_bus, dst_bus)
}

/// Build an adjacency map: bus_id → Vec<(neighbor_bus_id, branch_idx)>.
/// Excludes branches in `open_branches` and optionally `exclude_branch`.
fn build_adjacency_excluding(
    network: &PowerNetwork,
    open_branches: &HashSet<usize>,
    exclude_branch: Option<usize>,
) -> std::collections::HashMap<usize, Vec<(usize, usize)>> {
    let mut adj: std::collections::HashMap<usize, Vec<(usize, usize)>> =
        std::collections::HashMap::new();

    for (idx, branch) in network.branches.iter().enumerate() {
        if !branch.status {
            continue;
        }
        if open_branches.contains(&idx) {
            continue;
        }
        if exclude_branch == Some(idx) {
            continue;
        }
        adj.entry(branch.from_bus)
            .or_default()
            .push((branch.to_bus, idx));
        adj.entry(branch.to_bus)
            .or_default()
            .push((branch.from_bus, idx));
    }

    adj
}

/// BFS from `src` to `dst` through the adjacency map.
/// Returns the list of branch indices on the path, or empty if unreachable.
fn bfs_path(
    adj: &std::collections::HashMap<usize, Vec<(usize, usize)>>,
    src: usize,
    dst: usize,
) -> Vec<usize> {
    if src == dst {
        return vec![];
    }

    // BFS: queue holds (current_bus, path of branch indices to reach it)
    let mut visited: HashSet<usize> = HashSet::new();
    let mut queue: VecDeque<(usize, Vec<usize>)> = VecDeque::new();
    visited.insert(src);
    queue.push_back((src, vec![]));

    while let Some((bus, path)) = queue.pop_front() {
        if let Some(neighbors) = adj.get(&bus) {
            for &(nbr, branch_idx) in neighbors {
                if visited.contains(&nbr) {
                    continue;
                }
                let mut new_path = path.clone();
                new_path.push(branch_idx);
                if nbr == dst {
                    return new_path;
                }
                visited.insert(nbr);
                queue.push_back((nbr, new_path));
            }
        }
    }

    vec![]
}

// ─────────────────────────────────────────────────────────────────────────────
// Simulated annealing
// ─────────────────────────────────────────────────────────────────────────────

/// Simulated annealing reconfiguration.
///
/// Uses a linear congruential generator (LCG) for reproducible stochasticity
/// without the `rand` crate.
///
/// # SA parameters
/// - Initial temperature T₀ = 0.1 × initial_losses
/// - Final temperature T_f = 1 × 10⁻⁶
/// - Cooling factor α = 0.95 per outer iteration
/// - At each step: randomly toggle one operable switch, accept if Δ < 0 or
///   with probability exp(−Δ / T)
fn simulated_annealing(
    network: &PowerNetwork,
    initial_config: &NetworkConfiguration,
    switches: &[Switch],
    config: &ReconfigConfig,
    seed: u64,
) -> Result<(NetworkConfiguration, f64)> {
    let mut lcg = LcgRng::new(seed);

    let init_losses = compute_losses_estimate(
        network,
        &config_to_open_branch_indices(initial_config, switches)?,
    )?;

    let mut current_config = initial_config.clone();
    let mut current_losses = init_losses;
    let mut best_config = current_config.clone();
    let mut best_losses = current_losses;

    let t_init = 0.1 * init_losses.max(1e-3);
    let t_final = 1e-6_f64;
    let cooling = 0.95_f64;
    let mut temperature = t_init;

    let operable: Vec<&Switch> = switches.iter().filter(|s| s.can_operate).collect();
    if operable.is_empty() {
        return Ok((initial_config.clone(), init_losses));
    }

    let mut iter = 0;
    while temperature > t_final && iter < config.max_iterations {
        // Pick a random operable switch
        let sw_idx = (lcg.next_u64() as usize) % operable.len();
        let sw = operable[sw_idx];

        // Generate neighbour configuration by toggling this switch
        let (neighbour, valid) = toggle_switch(network, &current_config, switches, sw.id)?;
        if !valid {
            iter += 1;
            continue;
        }

        let neighbour_losses = compute_losses_estimate(
            network,
            &config_to_open_branch_indices(&neighbour, switches)?,
        )?;

        let delta = neighbour_losses - current_losses;

        // Metropolis acceptance criterion
        let accept = if delta < 0.0 {
            true
        } else {
            let prob = (-delta / temperature).exp();
            let r = lcg.next_f64();
            r < prob
        };

        if accept {
            current_config = neighbour;
            current_losses = neighbour_losses;

            if current_losses < best_losses {
                best_losses = current_losses;
                best_config = current_config.clone();
            }
        }

        temperature *= cooling;
        iter += 1;
    }

    Ok((best_config, best_losses))
}

// ─────────────────────────────────────────────────────────────────────────────
// Tabu search
// ─────────────────────────────────────────────────────────────────────────────

/// Tabu search reconfiguration.
///
/// Maintains a tabu list of the last 7 switch IDs that were operated.
/// At each iteration, evaluates all legal single-switch moves and selects
/// the best non-tabu neighbour.  If the best neighbour is globally better
/// than any seen so far (aspiration criterion), the tabu status is overridden.
fn tabu_search(
    network: &PowerNetwork,
    initial_config: &NetworkConfiguration,
    switches: &[Switch],
    config: &ReconfigConfig,
) -> Result<(NetworkConfiguration, f64)> {
    const TABU_SIZE: usize = 7;

    let init_losses = compute_losses_estimate(
        network,
        &config_to_open_branch_indices(initial_config, switches)?,
    )?;

    let mut current_config = initial_config.clone();
    let mut current_losses = init_losses;
    let mut best_config = current_config.clone();
    let mut best_losses = current_losses;

    let mut tabu_list: VecDeque<usize> = VecDeque::with_capacity(TABU_SIZE + 1);

    let operable: Vec<&Switch> = switches.iter().filter(|s| s.can_operate).collect();

    for _iter in 0..config.max_iterations {
        let mut best_candidate: Option<(NetworkConfiguration, f64, usize)> = None;

        for sw in &operable {
            let is_tabu = tabu_list.contains(&sw.id);

            let (neighbour, valid) = toggle_switch(network, &current_config, switches, sw.id)?;
            if !valid {
                continue;
            }

            let n_losses = compute_losses_estimate(
                network,
                &config_to_open_branch_indices(&neighbour, switches)?,
            )?;

            // Aspiration criterion: accept tabu move if it beats global best
            let aspiration = n_losses < best_losses;
            if is_tabu && !aspiration {
                continue;
            }

            let better = best_candidate
                .as_ref()
                .map(|(_, prev, _)| n_losses < *prev)
                .unwrap_or(true);

            if better {
                best_candidate = Some((neighbour, n_losses, sw.id));
            }
        }

        match best_candidate {
            None => break, // no improving non-tabu move
            Some((cfg, losses, sw_id)) => {
                current_config = cfg;
                current_losses = losses;

                if current_losses < best_losses {
                    best_losses = current_losses;
                    best_config = current_config.clone();
                }

                // Update tabu list
                tabu_list.push_back(sw_id);
                if tabu_list.len() > TABU_SIZE {
                    tabu_list.pop_front();
                }
            }
        }
    }

    Ok((best_config, best_losses))
}

// ─────────────────────────────────────────────────────────────────────────────
// Exhaustive search
// ─────────────────────────────────────────────────────────────────────────────

/// Exhaustive enumeration of all radial configurations.
fn exhaustive_solve(
    network: &PowerNetwork,
    base_config: &NetworkConfiguration,
    switches: &[Switch],
    config: &ReconfigConfig,
) -> Result<(NetworkConfiguration, f64, usize, bool)> {
    let solver = NetworkReconfiguration {
        config: ReconfigConfig {
            algorithm: ReconfigAlgorithm::Exhaustive,
            ..config.clone()
        },
    };

    let configs = solver.enumerate_radial_configs(network)?;

    let mut best_config = base_config.clone();
    let mut best_losses = f64::MAX;
    let n_iters = configs.len();

    for cfg in configs {
        let losses = solver.evaluate_config(network, &cfg)?;
        if losses < best_losses {
            best_losses = losses;
            best_config = cfg;
        }
    }

    if best_losses == f64::MAX {
        // No valid radial configs found; return base
        best_losses = compute_losses_estimate(
            network,
            &config_to_open_branch_indices(base_config, switches)?,
        )?;
        best_config = base_config.clone();
    }

    Ok((best_config, best_losses, n_iters, true))
}

// ─────────────────────────────────────────────────────────────────────────────
// Radiality verification
// ─────────────────────────────────────────────────────────────────────────────

/// Verify that the energised sub-graph (branches not in `open_switches`)
/// forms a spanning tree (connected, acyclic).
///
/// A graph on N buses is a tree iff it is connected and has exactly N−1 edges.
pub fn verify_radiality(
    network: &PowerNetwork,
    config: &NetworkConfiguration,
    switches: &[Switch],
) -> bool {
    let open_branches = config_to_open_branch_indices_unchecked(config, switches);
    let n_buses = network.bus_count();
    if n_buses == 0 {
        return true;
    }

    // Count energised branches
    let energised_count = network
        .branches
        .iter()
        .enumerate()
        .filter(|(idx, b)| b.status && !open_branches.contains(idx))
        .count();

    // For a spanning tree: exactly n_buses − 1 edges
    if energised_count != n_buses.saturating_sub(1) {
        return false;
    }

    // Check connectivity via BFS
    let open_set: HashSet<usize> = open_branches;
    let adj = build_adjacency_excluding(network, &open_set, None);

    let start_bus = match network.buses.first() {
        Some(b) => b.id,
        None => return true,
    };

    let mut visited: HashSet<usize> = HashSet::new();
    let mut queue = VecDeque::new();
    visited.insert(start_bus);
    queue.push_back(start_bus);

    while let Some(bus) = queue.pop_front() {
        if let Some(neighbors) = adj.get(&bus) {
            for &(nbr, _) in neighbors {
                if !visited.contains(&nbr) {
                    visited.insert(nbr);
                    queue.push_back(nbr);
                }
            }
        }
    }

    visited.len() == n_buses
}

// ─────────────────────────────────────────────────────────────────────────────
// Loss estimation
// ─────────────────────────────────────────────────────────────────────────────

/// Estimate total I²R losses \[MW\] using a simplified load-flow approximation.
///
/// For each energised branch: P_loss ≈ R × (P²_load_downstream / V²_nominal).
/// This is a first-order radial distribution loss formula, suitable for
/// ranking configurations without a full iterative power-flow solution.
fn compute_losses_estimate(network: &PowerNetwork, open_branches: &HashSet<usize>) -> Result<f64> {
    let base_mva = network.base_mva;
    let n_buses = network.bus_count();

    if n_buses == 0 {
        return Ok(0.0);
    }

    // Compute approximate downstream real-power load for each branch via BFS
    // from each branch's to-bus, accumulating load in the subtree.
    // This is an O(branches × buses) pass suitable for distribution networks.

    // Build adjacency excluding open branches (directed: from→to of the branch
    // defines feeder flow direction).
    let adj = build_adjacency_excluding(network, open_branches, None);

    let mut total_losses = 0.0_f64;

    for (idx, branch) in network.branches.iter().enumerate() {
        if !branch.status || open_branches.contains(&idx) {
            continue;
        }

        // Compute subtree power downstream of this branch's to_bus
        let downstream_mw = downstream_load_mw(network, branch.to_bus, branch.from_bus, &adj);

        // I²R loss approximation: P_loss = R × (P_down/V_nom)²
        // With V_nom ≈ 1.0 p.u. and P in p.u.:
        let p_down_pu = downstream_mw / base_mva;
        let loss_pu = branch.r * p_down_pu * p_down_pu;
        total_losses += loss_pu * base_mva;
    }

    Ok(total_losses)
}

/// Compute total real-power load \[MW\] in the subtree rooted at `root_bus`,
/// not traversing back through `parent_bus` (to avoid double-counting in loops).
fn downstream_load_mw(
    network: &PowerNetwork,
    root_bus: usize,
    parent_bus: usize,
    adj: &std::collections::HashMap<usize, Vec<(usize, usize)>>,
) -> f64 {
    let mut total = 0.0_f64;
    let mut visited: HashSet<usize> = HashSet::new();
    visited.insert(parent_bus);
    visited.insert(root_bus);

    let mut queue = VecDeque::new();
    queue.push_back(root_bus);

    while let Some(bus_id) = queue.pop_front() {
        // Add load at this bus
        if let Ok(idx) = network.bus_index(bus_id) {
            total += network.buses[idx].pd.0;
        }

        if let Some(neighbors) = adj.get(&bus_id) {
            for &(nbr, _) in neighbors {
                if !visited.contains(&nbr) {
                    visited.insert(nbr);
                    queue.push_back(nbr);
                }
            }
        }
    }

    total
}

// ─────────────────────────────────────────────────────────────────────────────
// Helper utilities
// ─────────────────────────────────────────────────────────────────────────────

/// Build the base configuration from switch data:
/// - Tie switches (normally_open = true) → open
/// - Section switches (normally_open = false) → closed
fn build_base_configuration(
    switches: &[Switch],
    network: &PowerNetwork,
) -> Result<NetworkConfiguration> {
    let open_sw: Vec<usize> = switches
        .iter()
        .filter(|s| s.normally_open)
        .map(|s| s.id)
        .collect();

    let closed_sw: Vec<usize> = switches
        .iter()
        .filter(|s| !s.normally_open)
        .map(|s| s.id)
        .collect();

    let config = NetworkConfiguration {
        open_switches: open_sw,
        closed_switches: closed_sw,
        is_radial: false,
        n_islands: 0,
    };

    let radial = verify_radiality(network, &config, switches);

    Ok(NetworkConfiguration {
        is_radial: radial,
        n_islands: if radial { 0 } else { 1 },
        ..config
    })
}

/// Convert open switch IDs in a configuration to open branch indices.
fn config_to_open_branch_indices(
    config: &NetworkConfiguration,
    switches: &[Switch],
) -> Result<HashSet<usize>> {
    let mut set = HashSet::new();
    for &sw_id in &config.open_switches {
        match switches.iter().find(|s| s.id == sw_id) {
            Some(s) => {
                set.insert(s.branch_idx);
            }
            None => {
                return Err(OxiGridError::InvalidParameter(format!(
                    "Switch ID {sw_id} not found in switch list"
                )));
            }
        }
    }
    Ok(set)
}

/// Infallible version — silently skips unknown switch IDs.
fn config_to_open_branch_indices_unchecked(
    config: &NetworkConfiguration,
    switches: &[Switch],
) -> HashSet<usize> {
    config
        .open_switches
        .iter()
        .filter_map(|&sw_id| {
            switches
                .iter()
                .find(|s| s.id == sw_id)
                .map(|s| s.branch_idx)
        })
        .collect()
}

/// Toggle a switch: if it is open → close it; if closed → open it.
/// Returns `(new_config, is_radial)`.
fn toggle_switch(
    network: &PowerNetwork,
    config: &NetworkConfiguration,
    switches: &[Switch],
    sw_id: usize,
) -> Result<(NetworkConfiguration, bool)> {
    let mut new_open = config.open_switches.clone();
    let mut new_closed = config.closed_switches.clone();

    if new_open.contains(&sw_id) {
        // Open → Close
        new_open.retain(|&id| id != sw_id);
        new_closed.push(sw_id);
    } else {
        // Closed → Open
        new_closed.retain(|&id| id != sw_id);
        new_open.push(sw_id);
    }

    new_open.sort_unstable();
    new_closed.sort_unstable();

    let trial = NetworkConfiguration {
        open_switches: new_open,
        closed_switches: new_closed,
        is_radial: false,
        n_islands: 0,
    };

    let radial = verify_radiality(network, &trial, switches);

    Ok((
        NetworkConfiguration {
            is_radial: radial,
            n_islands: if radial { 0 } else { 1 },
            ..trial
        },
        radial,
    ))
}

/// Compute the list of switch operations that transform `base` into `target`.
fn compute_switch_operations(
    base: &NetworkConfiguration,
    target: &NetworkConfiguration,
    switches: &[Switch],
) -> Vec<SwitchOperation> {
    let mut ops = Vec::new();

    for sw in switches {
        let was_open = base.open_switches.contains(&sw.id);
        let is_open = target.open_switches.contains(&sw.id);

        if was_open && !is_open {
            ops.push(SwitchOperation {
                switch_id: sw.id,
                action: SwitchAction::Close,
            });
        } else if !was_open && is_open {
            ops.push(SwitchOperation {
                switch_id: sw.id,
                action: SwitchAction::Open,
            });
        }
    }

    ops
}

// ─────────────────────────────────────────────────────────────────────────────
// LCG random number generator (no `rand` crate)
// ─────────────────────────────────────────────────────────────────────────────

/// Linear congruential generator (Knuth / MMIX parameters).
struct LcgRng {
    state: u64,
}

impl LcgRng {
    fn new(seed: u64) -> Self {
        Self {
            state: seed ^ 0x1234_5678_DEAD_BEEF,
        }
    }

    fn next_u64(&mut self) -> u64 {
        // LCG parameters from Knuth MMIX
        self.state = self
            .state
            .wrapping_mul(6_364_136_223_846_793_005)
            .wrapping_add(1_442_695_040_888_963_407);
        self.state
    }

    fn next_f64(&mut self) -> f64 {
        // Map to [0, 1)
        (self.next_u64() >> 11) as f64 / (1u64 << 53) as f64
    }
}

// ─────────────────────────────────────────────────────────────────────────────
// Tests
// ─────────────────────────────────────────────────────────────────────────────

#[cfg(test)]
mod tests {
    use super::*;
    use crate::network::branch::Branch;
    use crate::network::bus::{Bus, BusType};

    // ── helpers ──────────────────────────────────────────────────────────────

    /// Build a 5-bus radial feeder with one tie branch (bus 5 → bus 1).
    ///
    /// Topology:
    ///   Slack(1) — br0 — (2) — br1 — (3) — br2 — (4) — br3 — (5)
    ///    ↑                                                    |
    ///    └──────────────── br4 (tie) ─────────────────────────┘
    ///
    /// Branches:
    ///   br0: 1→2 (r=0.05)
    ///   br1: 2→3 (r=0.05)
    ///   br2: 3→4 (r=0.10)
    ///   br3: 4→5 (r=0.10)
    ///   br4: 5→1 (r=0.02) ← tie
    ///
    /// Switches:
    ///   sw0: br0, section
    ///   sw1: br1, section
    ///   sw2: br2, section
    ///   sw3: br3, section
    ///   sw4: br4, tie (normally open)
    fn five_bus_ring_network() -> (PowerNetwork, Vec<Switch>) {
        use crate::units::{Power, ReactivePower, Voltage};

        let mut net = PowerNetwork::new(100.0);

        let make_bus = |id: usize, bus_type: BusType, pd: f64| Bus {
            id,
            name: format!("Bus {id}"),
            bus_type,
            base_kv: Voltage(11.0),
            vm: 1.0,
            va: 0.0,
            pd: Power(pd),
            qd: ReactivePower(0.0),
            gs: 0.0,
            bs: 0.0,
            zone: None,
        };

        net.buses.push(make_bus(1, BusType::Slack, 0.0));
        net.buses.push(make_bus(2, BusType::PQ, 20.0));
        net.buses.push(make_bus(3, BusType::PQ, 20.0));
        net.buses.push(make_bus(4, BusType::PQ, 20.0));
        net.buses.push(make_bus(5, BusType::PQ, 20.0));

        let make_branch = |from: usize, to: usize, r: f64| Branch {
            from_bus: from,
            to_bus: to,
            r,
            x: 0.05,
            b: 0.0,
            rate_a: 0.0,
            rate_b: 0.0,
            rate_c: 0.0,
            tap: 0.0,
            shift: 0.0,
            status: true,
        };

        // br0..br3: section branches; br4: tie branch
        net.branches.push(make_branch(1, 2, 0.05)); // idx 0
        net.branches.push(make_branch(2, 3, 0.05)); // idx 1
        net.branches.push(make_branch(3, 4, 0.10)); // idx 2
        net.branches.push(make_branch(4, 5, 0.10)); // idx 3
        net.branches.push(make_branch(5, 1, 0.02)); // idx 4 — tie

        let switches = vec![
            Switch {
                id: 0,
                branch_idx: 0,
                normally_open: false,
                can_operate: true,
            },
            Switch {
                id: 1,
                branch_idx: 1,
                normally_open: false,
                can_operate: true,
            },
            Switch {
                id: 2,
                branch_idx: 2,
                normally_open: false,
                can_operate: true,
            },
            Switch {
                id: 3,
                branch_idx: 3,
                normally_open: false,
                can_operate: true,
            },
            Switch {
                id: 4,
                branch_idx: 4,
                normally_open: true,
                can_operate: true,
            },
        ];

        (net, switches)
    }

    // ── test_radiality_check ──────────────────────────────────────────────────

    #[test]
    fn test_radiality_check() {
        let (net, switches) = five_bus_ring_network();

        // Base configuration (tie switch 4 open) → radial tree
        let base_cfg = NetworkConfiguration {
            open_switches: vec![4],
            closed_switches: vec![0, 1, 2, 3],
            is_radial: false,
            n_islands: 0,
        };
        assert!(
            verify_radiality(&net, &base_cfg, &switches),
            "Base config should be radial"
        );

        // Close the tie switch (all closed) → creates a loop → NOT radial
        let loop_cfg = NetworkConfiguration {
            open_switches: vec![],
            closed_switches: vec![0, 1, 2, 3, 4],
            is_radial: false,
            n_islands: 0,
        };
        assert!(
            !verify_radiality(&net, &loop_cfg, &switches),
            "All-closed config should not be radial"
        );
    }

    // ── test_branch_exchange_basic ────────────────────────────────────────────

    #[test]
    fn test_branch_exchange_basic() {
        let (net, switches) = five_bus_ring_network();

        let config = ReconfigConfig {
            switches: switches.clone(),
            algorithm: ReconfigAlgorithm::BranchExchange,
            ..ReconfigConfig::default()
        };

        let solver = NetworkReconfiguration::new(config);
        let result = solver.solve(&net).expect("Branch exchange should succeed");

        assert!(
            result.loss_reduction_pct >= 0.0,
            "Loss reduction should be non-negative, got {:.4}%",
            result.loss_reduction_pct
        );
        assert!(result.best_configuration.is_radial, "Result must be radial");
    }

    // ── test_radiality_preserved ──────────────────────────────────────────────

    #[test]
    fn test_radiality_preserved() {
        let (net, switches) = five_bus_ring_network();

        for algo in [
            ReconfigAlgorithm::BranchExchange,
            ReconfigAlgorithm::SimulatedAnnealing,
            ReconfigAlgorithm::TabuSearch,
        ] {
            let config = ReconfigConfig {
                switches: switches.clone(),
                algorithm: algo,
                max_iterations: 50,
                ..ReconfigConfig::default()
            };

            let solver = NetworkReconfiguration::new(config);
            let result = solver.solve(&net).expect("Solve should succeed");

            assert!(
                result.best_configuration.is_radial,
                "{algo:?}: result configuration is not radial"
            );
        }
    }

    // ── test_sa_finds_better_than_branch_exchange ─────────────────────────────

    #[test]
    fn test_sa_finds_better_than_branch_exchange() {
        let (net, switches) = five_bus_ring_network();

        let make_result = |algo: ReconfigAlgorithm| {
            let config = ReconfigConfig {
                switches: switches.clone(),
                algorithm: algo,
                max_iterations: 200,
                ..ReconfigConfig::default()
            };
            NetworkReconfiguration::new(config)
                .solve(&net)
                .expect("solve")
        };

        let be_result = make_result(ReconfigAlgorithm::BranchExchange);
        let sa_result = make_result(ReconfigAlgorithm::SimulatedAnnealing);

        // Both algorithms should find radial configurations with non-negative losses.
        // On a small 5-bus ring, SA may not always beat BE (there are only 2
        // distinct radial configs).  We verify that SA converges and that its
        // result is within 2× of the branch-exchange result.
        let be_losses = be_result.optimized_losses_mw;
        let sa_losses = sa_result.optimized_losses_mw;

        assert!(
            sa_result.best_configuration.is_radial,
            "SA must return a radial configuration"
        );
        assert!(
            sa_losses >= 0.0,
            "SA losses must be non-negative, got {sa_losses:.4}"
        );
        // The optimal losses across both configurations is min(be, sa);
        // neither should be more than 2× worse than the other.
        let min_losses = be_losses.min(sa_losses).max(1e-9);
        let max_losses = be_losses.max(sa_losses);
        assert!(
            max_losses <= min_losses * 2.0 + 1e-6,
            "SA losses {sa_losses:.4} and BE losses {be_losses:.4} should be within 2× of each other"
        );
    }

    // ── test_switch_operations_valid ──────────────────────────────────────────

    #[test]
    fn test_switch_operations_valid() {
        let (net, switches) = five_bus_ring_network();

        let config = ReconfigConfig {
            switches: switches.clone(),
            algorithm: ReconfigAlgorithm::BranchExchange,
            ..ReconfigConfig::default()
        };

        let solver = NetworkReconfiguration::new(config);
        let result = solver.solve(&net).expect("solve");

        // No duplicate switch IDs in operations
        let op_ids: Vec<usize> = result
            .switch_operations
            .iter()
            .map(|o| o.switch_id)
            .collect();
        let unique_ids: HashSet<usize> = op_ids.iter().copied().collect();
        assert_eq!(
            op_ids.len(),
            unique_ids.len(),
            "Duplicate switch in operations: {:?}",
            op_ids
        );

        // Every operated switch must exist in the switch list
        for op in &result.switch_operations {
            assert!(
                switches.iter().any(|s| s.id == op.switch_id),
                "Unknown switch ID {} in operations",
                op.switch_id
            );
        }
    }

    // ── test_loss_sensitivity_positive ───────────────────────────────────────

    #[test]
    fn test_loss_sensitivity_positive() {
        let (net, _) = five_bus_ring_network();
        let n = net.bus_count();
        let v_mag = vec![1.0_f64; n];
        let v_ang = vec![0.0_f64; n];

        let sens = compute_branch_loss_sensitivity(&net, &v_mag, &v_ang);

        assert_eq!(sens.len(), net.branch_count());
        for &s in &sens {
            assert!(
                s >= 0.0,
                "Loss sensitivity must be non-negative, got {s:.6}"
            );
        }
    }

    // ── test_tabu_search_convergence ─────────────────────────────────────────

    #[test]
    fn test_tabu_search_convergence() {
        let (net, switches) = five_bus_ring_network();

        let config = ReconfigConfig {
            switches: switches.clone(),
            algorithm: ReconfigAlgorithm::TabuSearch,
            max_iterations: 100,
            ..ReconfigConfig::default()
        };

        let solver = NetworkReconfiguration::new(config);
        let result = solver.solve(&net).expect("Tabu search should succeed");

        assert!(result.converged, "Tabu search should report converged");
        assert!(
            result.optimized_losses_mw >= 0.0,
            "Losses must be non-negative"
        );
        assert!(
            result.best_configuration.is_radial,
            "Best config must be radial"
        );
    }

    // ── test_exhaustive_enumerate ─────────────────────────────────────────────

    #[test]
    fn test_exhaustive_enumerate() {
        let (net, switches) = five_bus_ring_network();

        let config = ReconfigConfig {
            switches: switches.clone(),
            algorithm: ReconfigAlgorithm::Exhaustive,
            ..ReconfigConfig::default()
        };

        let solver = NetworkReconfiguration::new(config);
        let configs = solver
            .enumerate_radial_configs(&net)
            .expect("enumerate should succeed");

        // With 1 tie switch, there are 2 candidate combinations (open or closed),
        // and at least 1 should be radial (the base config)
        assert!(
            !configs.is_empty(),
            "Should find at least one radial config"
        );

        for cfg in &configs {
            assert!(
                verify_radiality(&net, cfg, &switches),
                "All enumerated configs must be radial"
            );
        }
    }

    // ── test_loss_sensitivity_with_voltage ────────────────────────────────────

    #[test]
    fn test_loss_sensitivity_with_voltage() {
        let (net, _) = five_bus_ring_network();
        let n = net.bus_count();

        // Introduce voltage gradient
        let v_mag: Vec<f64> = (0..n).map(|i| 1.0 - 0.01 * i as f64).collect();
        let v_ang: Vec<f64> = (0..n).map(|i| -0.02 * i as f64).collect();

        let sens = compute_branch_loss_sensitivity(&net, &v_mag, &v_ang);

        assert_eq!(sens.len(), net.branch_count());
        // All sensitivities must be finite and non-negative
        for &s in &sens {
            assert!(s.is_finite(), "Sensitivity must be finite");
            assert!(s >= 0.0, "Sensitivity must be non-negative");
        }
    }
}