oxigrid 0.1.2

//! Grid Vulnerability Assessment.
//!
//! Provides comprehensive vulnerability assessment for power grids including:
//!
//! - Bridge detection (single points of failure)
//! - Branch betweenness centrality
//! - N-k contingency simulation with load-loss estimation
//! - Buldyrev percolation-based resilience index
//! - Worst-case attack scenario identification

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

// ─── Error ────────────────────────────────────────────────────────────────────

/// Errors from grid vulnerability assessment.
#[derive(Debug, thiserror::Error)]
pub enum VulnError {
    /// Network has no buses.
    #[error("network is empty — add buses and branches first")]
    EmptyNetwork,
    /// Invalid element referenced.
    #[error("invalid element id {0}")]
    InvalidElement(usize),
    /// Configuration error.
    #[error("configuration error: {0}")]
    Config(String),
}

// ─── Configuration ────────────────────────────────────────────────────────────

/// Threat model for vulnerability analysis.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum ThreatModel {
    /// Most impactful single or multi-element attack.
    WorstCase,
    /// Random element failures.
    RandomAttack,
    /// Targeted at critical infrastructure.
    IntelligentAttack,
    /// Correlated geographic failures.
    WeatherDriven,
}

/// Configuration for the vulnerability assessor.
#[derive(Debug, Clone)]
pub struct VulnerabilityConfig {
    /// Number of buses in the network.
    pub n_buses: usize,
    /// Number of branches in the network.
    pub n_branches: usize,
    /// Threat model used for attack scenario generation.
    pub threat_model: ThreatModel,
    /// Minimum load loss \[MW\] to include in report.
    pub impact_threshold: f64,
    /// Maximum number of simultaneous attack targets (N in N-k).
    pub max_attack_targets: usize,
}

impl VulnerabilityConfig {
    /// Create a default config for a network with the given size.
    pub fn new(n_buses: usize, n_branches: usize) -> Self {
        Self {
            n_buses,
            n_branches,
            threat_model: ThreatModel::WorstCase,
            impact_threshold: 1.0,
            max_attack_targets: 3,
        }
    }
}

// ─── Element Types ─────────────────────────────────────────────────────────────

/// Type of network element.
#[derive(Debug, Clone, PartialEq, Eq)]
pub enum ElementType {
    /// A network bus / node.
    Bus,
    /// A transmission branch (line or transformer).
    Branch,
    /// A generating unit.
    Generator,
    /// A substation (aggregated bus group).
    Substation,
    /// A transmission corridor (set of parallel branches).
    TransmissionCorridor,
}

impl std::fmt::Display for ElementType {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        match self {
            ElementType::Bus => write!(f, "Bus"),
            ElementType::Branch => write!(f, "Branch"),
            ElementType::Generator => write!(f, "Generator"),
            ElementType::Substation => write!(f, "Substation"),
            ElementType::TransmissionCorridor => write!(f, "TransmissionCorridor"),
        }
    }
}

// ─── Results ──────────────────────────────────────────────────────────────────

/// A critical network element with its vulnerability metrics.
#[derive(Debug, Clone)]
pub struct CriticalElement {
    /// Type of the element.
    pub element_type: ElementType,
    /// Unique identifier of the element.
    pub element_id: usize,
    /// Criticality score in \[0, 1\].
    pub criticality_score: f64,
    /// Estimated load loss \[MW\] if this element fails.
    pub load_loss_if_failed_mw: f64,
    /// Number of additional element failures triggered by cascading.
    pub cascading_failures: usize,
    /// Estimated restoration time \[h\].
    pub restoration_time_h: f64,
}

/// A concrete attack scenario targeting one or more elements.
#[derive(Debug, Clone)]
pub struct AttackScenario {
    /// Targeted elements (type, id).
    pub targets: Vec<(ElementType, usize)>,
    /// Total load loss \[MW\] from this scenario.
    pub total_impact_mw: f64,
    /// Buses affected (islanded or de-energised).
    pub affected_buses: Vec<usize>,
    /// Human-readable cascading chain.
    pub cascading_chain: Vec<String>,
    /// Estimated probability of occurrence.
    pub probability: f64,
    /// Short description of the attack type.
    pub attack_type: String,
}

/// Full vulnerability assessment result.
#[derive(Debug, Clone)]
pub struct VulnerabilityResult {
    /// Ranked list of critical elements (highest criticality first).
    pub critical_elements: Vec<CriticalElement>,
    /// Worst attack scenarios sorted by impact.
    pub worst_attack_scenarios: Vec<AttackScenario>,
    /// Network resilience index (0 = fragile, 1 = fully resilient).
    pub network_resilience_index: f64,
    /// Branch betweenness centrality: (branch\_id, score).
    pub betweenness_vulnerability: Vec<(usize, f64)>,
    /// Branch IDs that are bridges (single points of failure).
    pub single_points_of_failure: Vec<usize>,
    /// Network is N-k secure for this value of k.
    pub n_k_secure: usize,
    /// Human-readable hardening recommendations.
    pub recommendations: Vec<String>,
}

// ─── Assessor ─────────────────────────────────────────────────────────────────

/// Grid vulnerability assessor.
///
/// Accepts raw network topology (buses, branches, generators, loads, coordinates)
/// and computes a comprehensive [`VulnerabilityResult`].
pub struct GridVulnerabilityAssessor {
    config: VulnerabilityConfig,
    /// (from\_bus, to\_bus, rating\_mw)
    branches: Vec<(usize, usize, f64)>,
    /// (bus\_id, P\_mw)
    gen_buses: Vec<(usize, f64)>,
    /// (bus\_id, P\_mw)
    load_buses: Vec<(usize, f64)>,
    /// Geographic coordinates per bus.
    bus_coords: Vec<(f64, f64)>,
}

impl GridVulnerabilityAssessor {
    /// Create a new assessor with the given configuration.
    pub fn new(config: VulnerabilityConfig) -> Self {
        Self {
            config,
            branches: Vec::new(),
            gen_buses: Vec::new(),
            load_buses: Vec::new(),
            bus_coords: Vec::new(),
        }
    }

    /// Add a transmission branch from `from` to `to` with thermal rating \[MW\].
    pub fn add_branch(&mut self, from: usize, to: usize, rating_mw: f64) {
        self.branches.push((from, to, rating_mw));
    }

    /// Add a generator at `bus` with output \[MW\].
    pub fn add_generator(&mut self, bus: usize, p_mw: f64) {
        self.gen_buses.push((bus, p_mw));
    }

    /// Add a load at `bus` with demand \[MW\].
    pub fn add_load(&mut self, bus: usize, p_mw: f64) {
        self.load_buses.push((bus, p_mw));
    }

    /// Set geographic coordinates for each bus (indexed 0..n\_buses).
    pub fn set_coordinates(&mut self, coords: Vec<(f64, f64)>) {
        self.bus_coords = coords;
    }

    /// Run the full vulnerability assessment.
    pub fn assess(&self) -> Result<VulnerabilityResult, VulnError> {
        if self.config.n_buses == 0 {
            return Err(VulnError::EmptyNetwork);
        }

        // 1. Bridge detection → single points of failure
        let spof = self.find_bridges();

        // 2. Branch betweenness centrality
        let betweenness = self.branch_betweenness();
        let mut betweenness_pairs: Vec<(usize, f64)> =
            betweenness.iter().copied().enumerate().collect();
        betweenness_pairs
            .sort_by(|a, b| b.1.partial_cmp(&a.1).unwrap_or(std::cmp::Ordering::Equal));

        // 3. Critical elements — combine SPOF status and betweenness
        let total_load: f64 = self.load_buses.iter().map(|(_, p)| p).sum();
        let n_branches = self.branches.len();

        let mut critical_elements: Vec<CriticalElement> = self
            .branches
            .iter()
            .enumerate()
            .map(|(i, (from, to, rating))| {
                let is_bridge = spof.contains(&i);
                let bw_score = betweenness.get(i).copied().unwrap_or(0.0);
                // Normalise betweenness to [0,1]
                let max_bw = betweenness
                    .iter()
                    .cloned()
                    .fold(f64::NEG_INFINITY, f64::max)
                    .max(1.0);
                let bw_norm = bw_score / max_bw;

                let criticality = if is_bridge {
                    0.5 + 0.5 * bw_norm
                } else {
                    0.5 * bw_norm
                };

                // Load loss: connectivity-based heuristic for branch removal
                let load_loss = self.simulate_contingency(&[(ElementType::Branch, i)]);

                // Restoration: bridges take longer (10-24 h), others 4-8 h
                let restoration_time_h = if is_bridge {
                    let _ = (from, to, rating); // referenced for completeness
                    14.0 + 10.0 * bw_norm
                } else {
                    4.0 + 4.0 * bw_norm
                };

                // Cascading: count branches that share a bus with this one
                let cascading = self
                    .branches
                    .iter()
                    .enumerate()
                    .filter(|(j, (f2, t2, _))| {
                        *j != i && (*f2 == *from || *f2 == *to || *t2 == *from || *t2 == *to)
                    })
                    .count();

                CriticalElement {
                    element_type: ElementType::Branch,
                    element_id: i,
                    criticality_score: criticality.clamp(0.0, 1.0),
                    load_loss_if_failed_mw: load_loss,
                    cascading_failures: cascading,
                    restoration_time_h,
                }
            })
            .filter(|ce| ce.load_loss_if_failed_mw >= self.config.impact_threshold)
            .collect();

        // Also add generators as potential critical elements
        for (idx, (bus, p_mw)) in self.gen_buses.iter().enumerate() {
            let load_loss = (*p_mw).min(total_load);
            if load_loss >= self.config.impact_threshold {
                critical_elements.push(CriticalElement {
                    element_type: ElementType::Generator,
                    element_id: idx,
                    criticality_score: (load_loss / total_load.max(1.0)).clamp(0.0, 1.0),
                    load_loss_if_failed_mw: load_loss,
                    cascading_failures: 0,
                    restoration_time_h: 8.0 + 4.0 * (load_loss / total_load.max(1.0)),
                    // suppress unused warning
                });
                let _ = bus;
            }
        }

        critical_elements.sort_by(|a, b| {
            b.criticality_score
                .partial_cmp(&a.criticality_score)
                .unwrap_or(std::cmp::Ordering::Equal)
        });

        // 4. Attack scenarios
        let worst_attack_scenarios =
            self.generate_attack_scenarios(&betweenness, &spof, n_branches);

        // 5. Resilience index (Buldyrev percolation)
        let network_resilience_index = self.resilience_index();

        // 6. N-k security level
        let n_k_secure = self.compute_nk_secure();

        // 7. Recommendations
        let recommendations = self.generate_recommendations(
            &critical_elements,
            &spof,
            network_resilience_index,
            n_k_secure,
        );

        Ok(VulnerabilityResult {
            critical_elements,
            worst_attack_scenarios,
            network_resilience_index,
            betweenness_vulnerability: betweenness_pairs,
            single_points_of_failure: spof,
            n_k_secure,
            recommendations,
        })
    }

    // ─── Internal helpers ───────────────────────────────────────────────────

    /// Find branches that are bridges using DFS low-link (Tarjan algorithm).
    ///
    /// Returns branch IDs (indices into `self.branches`) that are bridges.
    fn find_bridges(&self) -> Vec<usize> {
        let n = self.config.n_buses;
        if n == 0 || self.branches.is_empty() {
            return Vec::new();
        }

        // Build adjacency: node → Vec<(neighbor, branch_idx)>
        let mut adj: Vec<Vec<(usize, usize)>> = vec![Vec::new(); n];
        for (idx, (from, to, _)) in self.branches.iter().enumerate() {
            if *from < n && *to < n {
                adj[*from].push((*to, idx));
                adj[*to].push((*from, idx));
            }
        }

        let mut disc = vec![usize::MAX; n];
        let mut low = vec![usize::MAX; n];
        let mut visited = vec![false; n];
        let mut bridges = HashSet::new();
        let mut timer = 0usize;

        for start in 0..n {
            if !visited[start] {
                self.dfs_bridge(
                    start,
                    usize::MAX, // parent edge index
                    &adj,
                    &mut disc,
                    &mut low,
                    &mut visited,
                    &mut timer,
                    &mut bridges,
                );
            }
        }

        let mut result: Vec<usize> = bridges.into_iter().collect();
        result.sort_unstable();
        result
    }

    /// Recursive DFS for bridge detection.
    #[allow(clippy::too_many_arguments)]
    fn dfs_bridge(
        &self,
        u: usize,
        parent_edge: usize,
        adj: &[Vec<(usize, usize)>],
        disc: &mut Vec<usize>,
        low: &mut Vec<usize>,
        visited: &mut Vec<bool>,
        timer: &mut usize,
        bridges: &mut HashSet<usize>,
    ) {
        visited[u] = true;
        disc[u] = *timer;
        low[u] = *timer;
        *timer += 1;

        for &(v, edge_idx) in &adj[u] {
            if edge_idx == parent_edge {
                continue;
            }
            if visited[v] {
                low[u] = low[u].min(disc[v]);
            } else {
                self.dfs_bridge(v, edge_idx, adj, disc, low, visited, timer, bridges);
                low[u] = low[u].min(low[v]);
                if low[v] > disc[u] {
                    bridges.insert(edge_idx);
                }
            }
        }
    }

    /// Compute branch betweenness centrality.
    ///
    /// For every source node, BFS shortest-path tree is computed; each branch
    /// on a path from source to all other nodes accumulates a fractional count.
    fn branch_betweenness(&self) -> Vec<f64> {
        let n = self.config.n_buses;
        let nb = self.branches.len();
        if n == 0 || nb == 0 {
            return vec![0.0; nb];
        }

        let mut adj: Vec<Vec<(usize, usize)>> = vec![Vec::new(); n];
        for (idx, (from, to, _)) in self.branches.iter().enumerate() {
            if *from < n && *to < n {
                adj[*from].push((*to, idx));
                adj[*to].push((*from, idx));
            }
        }

        let mut betweenness = vec![0.0f64; nb];

        for src in 0..n {
            // BFS to compute shortest path counts and predecessors
            let mut dist = vec![i64::MAX; n];
            let mut sigma = vec![0i64; n]; // number of shortest paths from src
            let mut pred: Vec<Vec<usize>> = vec![Vec::new(); n]; // predecessor nodes
            let mut queue = VecDeque::new();
            let mut order = Vec::new(); // BFS order

            dist[src] = 0;
            sigma[src] = 1;
            queue.push_back(src);

            while let Some(u) = queue.pop_front() {
                order.push(u);
                for &(v, _) in &adj[u] {
                    if dist[v] == i64::MAX {
                        dist[v] = dist[u] + 1;
                        queue.push_back(v);
                    }
                    if dist[v] == dist[u] + 1 {
                        sigma[v] += sigma[u];
                        pred[v].push(u);
                    }
                }
            }

            // Back-propagation of pair dependencies
            let mut delta = vec![0.0f64; n];
            while let Some(w) = order.pop() {
                for &u in &pred[w] {
                    let frac = if sigma[w] > 0 {
                        (sigma[u] as f64 / sigma[w] as f64) * (1.0 + delta[w])
                    } else {
                        0.0
                    };
                    delta[u] += frac;

                    // Find the branch between u and w, add contribution
                    for &(v, edge_idx) in &adj[u] {
                        if v == w && edge_idx < nb {
                            betweenness[edge_idx] += frac;
                        }
                    }
                }
            }
        }

        // Normalise (undirected: each pair counted twice)
        let norm = if n > 2 { 2.0 } else { 1.0 };
        for b in &mut betweenness {
            *b /= norm;
        }
        betweenness
    }

    /// Simulate the removal of `targets` and estimate load loss \[MW\].
    ///
    /// Uses a simple connectivity-based heuristic: load buses disconnected from
    /// all generator buses after branch/bus removals lose their full load.
    pub fn simulate_contingency(&self, targets: &[(ElementType, usize)]) -> f64 {
        let n = self.config.n_buses;
        if n == 0 {
            return 0.0;
        }

        // Build adjacency excluding removed branches / buses
        let removed_branches: HashSet<usize> = targets
            .iter()
            .filter_map(|(et, id)| {
                if *et == ElementType::Branch {
                    Some(*id)
                } else {
                    None
                }
            })
            .collect();
        let removed_buses: HashSet<usize> = targets
            .iter()
            .filter_map(|(et, id)| {
                if *et == ElementType::Bus {
                    Some(*id)
                } else {
                    None
                }
            })
            .collect();
        let removed_gens: HashSet<usize> = targets
            .iter()
            .filter_map(|(et, id)| {
                if *et == ElementType::Generator {
                    Some(*id)
                } else {
                    None
                }
            })
            .collect();

        let mut adj: Vec<Vec<usize>> = vec![Vec::new(); n];
        for (idx, (from, to, _)) in self.branches.iter().enumerate() {
            if removed_branches.contains(&idx) {
                continue;
            }
            if *from >= n || *to >= n {
                continue;
            }
            if removed_buses.contains(from) || removed_buses.contains(to) {
                continue;
            }
            adj[*from].push(*to);
            adj[*to].push(*from);
        }

        // Determine which generation buses still operate
        let active_gen_buses: HashSet<usize> = self
            .gen_buses
            .iter()
            .enumerate()
            .filter(|(idx, _)| !removed_gens.contains(idx))
            .map(|(_, (bus, _))| *bus)
            .filter(|b| !removed_buses.contains(b))
            .collect();

        // BFS connectivity: for each load bus, find if it can reach any gen bus
        let mut load_loss = 0.0f64;
        for (load_bus, p_mw) in &self.load_buses {
            if removed_buses.contains(load_bus) {
                load_loss += p_mw;
                continue;
            }
            // BFS from load_bus
            let mut visited = HashSet::new();
            let mut queue = VecDeque::new();
            queue.push_back(*load_bus);
            visited.insert(*load_bus);
            let mut connected_to_gen = false;

            'bfs: while let Some(node) = queue.pop_front() {
                if active_gen_buses.contains(&node) {
                    connected_to_gen = true;
                    break 'bfs;
                }
                for &nb in &adj[node] {
                    if !visited.contains(&nb) {
                        visited.insert(nb);
                        queue.push_back(nb);
                    }
                }
            }

            if !connected_to_gen {
                load_loss += p_mw;
            }
        }

        load_loss
    }

    /// Compute network resilience using Buldyrev percolation approach.
    ///
    /// Estimates the fraction of load that survives random branch failures by
    /// computing the giant connected component fraction under incremental removal.
    fn resilience_index(&self) -> f64 {
        let n = self.config.n_buses;
        let nb = self.branches.len();
        if n == 0 || nb == 0 {
            return 0.0;
        }

        // Build initial adjacency list
        let mut adj: Vec<HashSet<usize>> = vec![HashSet::new(); n];
        for (from, to, _) in &self.branches {
            if *from < n && *to < n {
                adj[*from].insert(*to);
                adj[*to].insert(*from);
            }
        }

        // Giant component size without any removals
        let gcc_full = self.giant_component_size(&adj, n);
        if gcc_full == 0 {
            return 0.0;
        }

        // Area under the percolation curve (integrate GCC size / GCC_full over removals)
        // Use LCG for deterministic pseudorandom order
        let mut rng = 6364136223846793005u64;
        let mut order: Vec<usize> = (0..nb).collect();

        // Fisher-Yates shuffle with LCG
        for i in (1..nb).rev() {
            rng = rng
                .wrapping_mul(6364136223846793005u64)
                .wrapping_add(1442695040888963407u64);
            let j = (rng >> 33) as usize % (i + 1);
            order.swap(i, j);
        }

        // Sample 10 evenly spaced removal fractions
        let sample_points = 10usize;
        let mut area = 0.0f64;
        let mut current_adj: Vec<HashSet<usize>> = adj.clone();

        // For resilience, we integrate the GCC curve: higher area = more resilient
        area += 1.0f64; // at 0 removals = full network (fraction = 1.0)

        for k in 1..=sample_points {
            // Remove k/sample_points fraction of branches
            let remove_up_to = k * nb / sample_points;
            let remove_from = (k - 1) * nb / sample_points;
            for &branch_idx in &order[remove_from..remove_from.max(remove_up_to).min(order.len())] {
                let (from, to, _) = self.branches[branch_idx];
                if from < n && to < n {
                    current_adj[from].remove(&to);
                    current_adj[to].remove(&from);
                }
            }
            let gcc = self.giant_component_size(&current_adj, n);
            let fraction = gcc as f64 / gcc_full as f64;
            area += fraction;
        }

        // Normalise to [0,1]: area / (sample_points + 1 steps)
        (area / (sample_points + 1) as f64).clamp(0.0, 1.0)
    }

    /// Compute size of the giant connected component.
    fn giant_component_size(&self, adj: &[HashSet<usize>], n: usize) -> usize {
        if n == 0 {
            return 0;
        }
        let mut visited = vec![false; n];
        let mut max_size = 0usize;

        for start in 0..n {
            if visited[start] {
                continue;
            }
            let mut size = 0usize;
            let mut queue = VecDeque::new();
            queue.push_back(start);
            visited[start] = true;

            while let Some(u) = queue.pop_front() {
                size += 1;
                for &v in &adj[u] {
                    if v < n && !visited[v] {
                        visited[v] = true;
                        queue.push_back(v);
                    }
                }
            }
            if size > max_size {
                max_size = size;
            }
        }
        max_size
    }

    /// Generate attack scenarios based on the configured threat model.
    fn generate_attack_scenarios(
        &self,
        betweenness: &[f64],
        bridges: &[usize],
        n_branches: usize,
    ) -> Vec<AttackScenario> {
        let mut scenarios = Vec::new();
        let max_k = self.config.max_attack_targets.min(n_branches);

        // N-1 scenarios for all branches (worst case)
        let mut n1_scores: Vec<(usize, f64)> = (0..n_branches)
            .map(|i| {
                let loss = self.simulate_contingency(&[(ElementType::Branch, i)]);
                (i, loss)
            })
            .filter(|(_, loss)| *loss >= self.config.impact_threshold)
            .collect();
        n1_scores.sort_by(|a, b| b.1.partial_cmp(&a.1).unwrap_or(std::cmp::Ordering::Equal));

        for (branch_idx, loss) in n1_scores.iter().take(5) {
            let (from, to, _) = self.branches[*branch_idx];
            let is_bridge = bridges.contains(branch_idx);
            let bw = betweenness.get(*branch_idx).copied().unwrap_or(0.0);
            let prob = match self.config.threat_model {
                ThreatModel::RandomAttack => 1.0 / n_branches as f64,
                ThreatModel::WorstCase | ThreatModel::IntelligentAttack => {
                    0.05 * (bw / betweenness.iter().cloned().fold(0.0f64, f64::max).max(1.0))
                }
                ThreatModel::WeatherDriven => 0.02,
            };
            scenarios.push(AttackScenario {
                targets: vec![(ElementType::Branch, *branch_idx)],
                total_impact_mw: *loss,
                affected_buses: vec![from, to],
                cascading_chain: vec![format!(
                    "Branch {branch_idx} ({from}→{to}) removed{}",
                    if is_bridge { " [BRIDGE]" } else { "" }
                )],
                probability: prob.clamp(0.0, 1.0),
                attack_type: "N-1 Branch Outage".to_string(),
            });
        }

        // N-k scenarios: combine top-betweenness branches
        if max_k >= 2 && !n1_scores.is_empty() {
            let top_k: Vec<usize> = n1_scores.iter().take(max_k).map(|(idx, _)| *idx).collect();
            if top_k.len() >= 2 {
                let targets: Vec<(ElementType, usize)> =
                    top_k.iter().map(|&i| (ElementType::Branch, i)).collect();
                let loss = self.simulate_contingency(&targets);
                if loss >= self.config.impact_threshold {
                    let chain: Vec<String> = top_k
                        .iter()
                        .map(|&i| {
                            let (f, t, _) = self.branches[i];
                            format!("Branch {i} ({f}→{t})")
                        })
                        .collect();
                    scenarios.push(AttackScenario {
                        targets,
                        total_impact_mw: loss,
                        affected_buses: top_k
                            .iter()
                            .flat_map(|&i| {
                                let (f, t, _) = self.branches[i];
                                vec![f, t]
                            })
                            .collect::<HashSet<_>>()
                            .into_iter()
                            .collect(),
                        cascading_chain: chain,
                        probability: 0.001,
                        attack_type: format!("N-{} Coordinated Attack", top_k.len()),
                    });
                }
            }
        }

        scenarios.sort_by(|a, b| {
            b.total_impact_mw
                .partial_cmp(&a.total_impact_mw)
                .unwrap_or(std::cmp::Ordering::Equal)
        });
        scenarios
    }

    /// Determine the N-k security level by incremental contingency testing.
    fn compute_nk_secure(&self) -> usize {
        let nb = self.branches.len();
        if nb == 0 {
            return 0;
        }

        // N-1: check all single branch outages
        let n1_ok = (0..nb).all(|i| {
            self.simulate_contingency(&[(ElementType::Branch, i)]) < self.config.impact_threshold
        });
        if !n1_ok {
            return 0;
        }

        // N-2: try a representative sample of pairs (limit to 100 pairs)
        let pairs: Vec<(usize, usize)> = {
            let mut ps = Vec::new();
            'outer: for i in 0..nb {
                for j in (i + 1)..nb {
                    ps.push((i, j));
                    if ps.len() >= 100 {
                        break 'outer;
                    }
                }
            }
            ps
        };

        let n2_ok = pairs.iter().all(|&(i, j)| {
            self.simulate_contingency(&[(ElementType::Branch, i), (ElementType::Branch, j)])
                < self.config.impact_threshold
        });

        if n2_ok {
            2
        } else {
            1
        }
    }

    /// Generate hardening recommendations.
    fn generate_recommendations(
        &self,
        critical_elements: &[CriticalElement],
        bridges: &[usize],
        resilience: f64,
        n_k_secure: usize,
    ) -> Vec<String> {
        let mut recs = Vec::new();

        if !bridges.is_empty() {
            recs.push(format!(
                "Install redundant paths for {} bridge branch(es) to eliminate single points of failure.",
                bridges.len()
            ));
        }

        if resilience < 0.5 {
            recs.push(
                "Network resilience is low (<0.5). Consider meshing the topology with additional tie lines."
                    .to_string(),
            );
        } else if resilience < 0.75 {
            recs.push(
                "Network resilience is moderate. Adding 1-2 strategic tie lines would improve robustness."
                    .to_string(),
            );
        }

        if n_k_secure == 0 {
            recs.push(
                "Network is NOT N-1 secure. Immediate corrective actions required.".to_string(),
            );
        } else if n_k_secure == 1 {
            recs.push(
                "Network is N-1 secure but not N-2. Consider N-2 hardening for critical corridors."
                    .to_string(),
            );
        }

        let top_critical = critical_elements.iter().take(3);
        for ce in top_critical {
            if ce.criticality_score > 0.7 {
                recs.push(format!(
                    "Priority hardening: {} #{} (criticality={:.2}, load loss if failed={:.1} MW).",
                    ce.element_type, ce.element_id, ce.criticality_score, ce.load_loss_if_failed_mw
                ));
            }
        }

        if recs.is_empty() {
            recs.push(
                "Network appears resilient. Maintain N-2 monitoring and annual review.".to_string(),
            );
        }

        recs
    }
}

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

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

    fn simple_chain() -> GridVulnerabilityAssessor {
        // 3-bus chain: 0-1-2, bus 1 is a load bus fed by gen at 0 and 2
        // Branch 0: 0-1 (bridge), Branch 1: 1-2 (bridge)
        let config = VulnerabilityConfig::new(3, 2);
        let mut a = GridVulnerabilityAssessor::new(config);
        a.add_branch(0, 1, 100.0);
        a.add_branch(1, 2, 100.0);
        a.add_generator(0, 50.0);
        a.add_load(1, 30.0);
        a
    }

    fn simple_ring() -> GridVulnerabilityAssessor {
        // 4-bus ring: 0-1-2-3-0, no bridges
        let config = VulnerabilityConfig {
            n_buses: 4,
            n_branches: 4,
            threat_model: ThreatModel::WorstCase,
            impact_threshold: 0.1,
            max_attack_targets: 2,
        };
        let mut a = GridVulnerabilityAssessor::new(config);
        a.add_branch(0, 1, 100.0);
        a.add_branch(1, 2, 100.0);
        a.add_branch(2, 3, 100.0);
        a.add_branch(3, 0, 100.0);
        a.add_generator(0, 100.0);
        a.add_load(2, 50.0);
        a
    }

    // Test 1: Bridge detection on a chain graph (both branches are bridges)
    #[test]
    fn test_single_point_of_failure_chain() {
        let a = simple_chain();
        let result = a.assess().expect("assessment ok");
        // Both branches in a chain are bridges
        assert!(
            !result.single_points_of_failure.is_empty(),
            "Chain graph must have bridges (SPOFs): got {:?}",
            result.single_points_of_failure
        );
        assert!(
            result.single_points_of_failure.contains(&0),
            "Branch 0 must be a bridge"
        );
        assert!(
            result.single_points_of_failure.contains(&1),
            "Branch 1 must be a bridge"
        );
    }

    // Test 2: Ring graph has no bridges
    #[test]
    fn test_ring_no_bridges() {
        let a = simple_ring();
        let result = a.assess().expect("assessment ok");
        assert!(
            result.single_points_of_failure.is_empty(),
            "Ring graph has no bridges: got {:?}",
            result.single_points_of_failure
        );
    }

    // Test 3: N-1 worst contingency — load loss detected
    #[test]
    fn test_n1_worst_contingency_identified() {
        let a = simple_chain();
        let result = a.assess().expect("assessment ok");
        // Should have at least one attack scenario since branch removals disconnect load
        assert!(
            !result.worst_attack_scenarios.is_empty(),
            "Should have N-1 attack scenarios"
        );
        let max_impact = result
            .worst_attack_scenarios
            .iter()
            .map(|s| s.total_impact_mw)
            .fold(f64::NEG_INFINITY, f64::max);
        assert!(
            max_impact > 0.0,
            "Worst contingency must have positive load loss: {max_impact}"
        );
    }

    // Test 4: Betweenness — highly connected middle bus branches should score highest
    #[test]
    fn test_betweenness_centrality() {
        // 5-bus star: all branches connect to central node 0
        let config = VulnerabilityConfig {
            n_buses: 5,
            n_branches: 4,
            threat_model: ThreatModel::WorstCase,
            impact_threshold: 0.1,
            max_attack_targets: 2,
        };
        let mut a = GridVulnerabilityAssessor::new(config);
        a.add_branch(0, 1, 100.0);
        a.add_branch(0, 2, 100.0);
        a.add_branch(0, 3, 100.0);
        a.add_branch(0, 4, 100.0);
        a.add_generator(0, 200.0);
        a.add_load(1, 10.0);
        a.add_load(2, 10.0);
        a.add_load(3, 10.0);
        a.add_load(4, 10.0);

        let result = a.assess().expect("assessment ok");
        // All branches in a star connect through hub, so they should all have equal high betweenness
        // Verify betweenness list is non-empty
        assert!(
            !result.betweenness_vulnerability.is_empty(),
            "Betweenness must be computed"
        );
        let max_bw = result
            .betweenness_vulnerability
            .iter()
            .map(|(_, s)| *s)
            .fold(f64::NEG_INFINITY, f64::max);
        assert!(
            max_bw > 0.0,
            "At least one branch must have positive betweenness: {max_bw}"
        );
    }

    // Test 5: Resilience index — fully connected ring approaches 1.0
    #[test]
    fn test_resilience_fully_connected() {
        // Dense 5-bus complete graph
        let config = VulnerabilityConfig {
            n_buses: 4,
            n_branches: 6,
            threat_model: ThreatModel::WorstCase,
            impact_threshold: 0.1,
            max_attack_targets: 2,
        };
        let mut a = GridVulnerabilityAssessor::new(config);
        // K4 complete graph: 6 edges
        a.add_branch(0, 1, 100.0);
        a.add_branch(0, 2, 100.0);
        a.add_branch(0, 3, 100.0);
        a.add_branch(1, 2, 100.0);
        a.add_branch(1, 3, 100.0);
        a.add_branch(2, 3, 100.0);
        a.add_generator(0, 100.0);
        a.add_load(3, 50.0);

        let result = a.assess().expect("assessment ok");
        assert!(
            result.network_resilience_index > 0.5,
            "Fully connected graph should have resilience > 0.5: got {:.4}",
            result.network_resilience_index
        );
    }

    // Test 6: Recommendations generated for critical elements
    #[test]
    fn test_recommendations_generated() {
        let a = simple_chain();
        let result = a.assess().expect("assessment ok");
        assert!(
            !result.recommendations.is_empty(),
            "Should generate recommendations for a chain with bridges"
        );
        // The recommendation should mention bridge or redundant paths
        let has_bridge_rec = result.recommendations.iter().any(|r| {
            r.to_lowercase().contains("bridge")
                || r.to_lowercase().contains("redundant")
                || r.to_lowercase().contains("n-1")
                || r.to_lowercase().contains("single point")
        });
        assert!(
            has_bridge_rec,
            "Recommendation must mention bridge hardening. Got: {:?}",
            result.recommendations
        );
    }

    // Test 7: Empty network returns error
    #[test]
    fn test_empty_network_error() {
        let config = VulnerabilityConfig::new(0, 0);
        let a = GridVulnerabilityAssessor::new(config);
        assert!(a.assess().is_err(), "Empty network must return error");
    }

    // Test 8: Resilience index for single-branch network
    #[test]
    fn test_single_branch_resilience() {
        let config = VulnerabilityConfig::new(2, 1);
        let mut a = GridVulnerabilityAssessor::new(config);
        a.add_branch(0, 1, 100.0);
        a.add_generator(0, 50.0);
        a.add_load(1, 30.0);
        let result = a.assess().expect("assessment ok");
        // Single branch = single bridge, low resilience
        assert!(
            result.single_points_of_failure.contains(&0),
            "Single branch must be a bridge"
        );
        assert!(
            result.network_resilience_index <= 1.0,
            "Resilience must be <= 1.0"
        );
    }
}