kglite 0.10.9

// src/graph/graph_algorithms.rs
//! Graph algorithms module providing path finding and connectivity analysis.

use crate::datatypes::values::Value;
use crate::graph::schema::{DirGraph, InternedKey};
use crate::graph::storage::GraphRead;
use petgraph::algo::kosaraju_scc;
use petgraph::graph::NodeIndex;
use std::collections::{HashMap, HashSet};
use std::time::Instant;

/// Standard timeout error message for graph algorithms.
/// Mirrors the MATCH timeout text in `cypher::executor::mod::check_deadline`,
/// adapted for procedure context (no anchor hint — large graphs may simply
/// not converge within the default 20s).
pub fn algorithm_timeout_err() -> String {
    "CALL procedure timed out. Pass timeout_ms=N to cypher() to extend, \
     or timeout_ms=0 to disable the deadline. Subgraph scoping is not \
     yet supported — large graphs may not converge within the default 20s."
        .to_string()
}

// ============================================================================
// Path Filtering Helpers
// ============================================================================

/// Pre-intern connection type strings into InternedKeys for fast comparison.
fn intern_connection_types(connection_types: Option<&[String]>) -> Option<Vec<InternedKey>> {
    connection_types.map(|types| types.iter().map(|t| InternedKey::from_str(t)).collect())
}

/// Get undirected neighbors filtered by edge connection type.
/// When connection_types is None, returns all neighbors (equivalent to
/// neighbors_undirected).
///
/// Both branches deduplicate — petgraph's `neighbors_undirected` walks
/// every incident edge so parallel edges and a→b/b→a pairs each appear
/// twice; the filtered branch concatenates Outgoing + Incoming and has
/// the same property. Without dedup, undirected `shortestPath` and
/// `all_paths` over a bidirectional pair surfaced duplicate
/// (A, B) / (B, A) entries during enumeration (B4) — the visited
/// bitmap downstream caught most cases for `shortestPath`, but
/// `all_paths` paid wasted DFS work per duplicate.
///
/// Sort + dedup is faster than a presence-set probe for the typical
/// small-degree case (n ≲ 32) because the in-place comparison fits in
/// cache; insertion order is not load-bearing for any caller (BFS,
/// DFS path enumeration use set-membership, not order).
fn filtered_neighbors_undirected(
    graph: &DirGraph,
    node: NodeIndex,
    connection_types: Option<&[InternedKey]>,
) -> Vec<NodeIndex> {
    use petgraph::Direction;
    let g = &graph.graph;
    let mut neighbors: Vec<NodeIndex> = match connection_types {
        None => g.neighbors_undirected(node).collect(),
        Some(types) => {
            let mut n = Vec::new();
            for edge in g.edges_directed(node, Direction::Outgoing) {
                if types.iter().any(|t| *t == edge.weight().connection_type) {
                    n.push(edge.target());
                }
            }
            for edge in g.edges_directed(node, Direction::Incoming) {
                if types.iter().any(|t| *t == edge.weight().connection_type) {
                    n.push(edge.source());
                }
            }
            n
        }
    };
    if neighbors.len() > 1 {
        neighbors.sort_unstable();
        neighbors.dedup();
    }
    neighbors
}

/// Get directed (outgoing only) neighbors filtered by edge connection type.
fn filtered_neighbors_outgoing(
    graph: &DirGraph,
    node: NodeIndex,
    connection_types: Option<&[InternedKey]>,
) -> Vec<NodeIndex> {
    use petgraph::Direction;
    let g = &graph.graph;
    match connection_types {
        None => g.neighbors_directed(node, Direction::Outgoing).collect(),
        Some(types) => g
            .edges_directed(node, Direction::Outgoing)
            .filter(|e| types.iter().any(|t| *t == e.weight().connection_type))
            .map(|e| e.target())
            .collect(),
    }
}

/// Check if a node passes the via_types filter.
/// Source and target should be excluded from this check by the caller.
fn node_passes_via_filter(
    graph: &DirGraph,
    node: NodeIndex,
    via_types: &Option<HashSet<&str>>,
) -> bool {
    match via_types {
        None => true,
        Some(types) => {
            if let Some(node_data) = graph.graph.node_weight(node) {
                types.contains(node_data.node_type_str(&graph.interner))
            } else {
                false
            }
        }
    }
}

/// Result of a path finding operation
#[derive(Debug, Clone)]
pub struct PathResult {
    /// The path as a sequence of node indices
    pub path: Vec<NodeIndex>,
    /// The total cost/length of the path
    pub cost: usize,
}

/// Information about a node in a path (for Python output)
#[derive(Debug, Clone)]
pub struct PathNodeInfo {
    pub node_type: String,
    pub title: String,
    pub id: Value,
}

/// Find the shortest path between two nodes using undirected BFS.
/// This treats the graph as undirected, finding connections in either direction.
/// Returns None if no path exists.
///
/// # Arguments
/// * `connection_types` - Only traverse edges of these types (None = all)
/// * `via_types` - Only traverse through nodes of these types (None = all)
pub fn shortest_path(
    graph: &DirGraph,
    source: NodeIndex,
    target: NodeIndex,
    connection_types: Option<&[String]>,
    via_types: Option<&[String]>,
    deadline: Option<Instant>,
) -> Option<PathResult> {
    let via_set: Option<HashSet<&str>> =
        via_types.map(|vt| vt.iter().map(|s| s.as_str()).collect());
    let interned = intern_connection_types(connection_types);
    let path = reconstruct_path_bfs(
        graph,
        source,
        target,
        interned.as_deref(),
        &via_set,
        deadline,
    )?;
    let cost = path.len().saturating_sub(1);

    Some(PathResult { path, cost })
}

/// Find the shortest path LENGTH between two nodes using undirected BFS.
/// Only returns the hop count, avoiding parent tracking and path reconstruction.
/// Uses level-by-level BFS to avoid per-node distance tracking.
pub fn shortest_path_cost(graph: &DirGraph, source: NodeIndex, target: NodeIndex) -> Option<usize> {
    if source == target {
        return Some(0);
    }

    let node_bound = graph.graph.node_bound();
    let mut visited: Vec<bool> = vec![false; node_bound];

    let target_idx = target.index();

    // Level-by-level BFS using two alternating vectors (avoids VecDeque overhead)
    let mut current_level: Vec<usize> = vec![source.index()];
    let mut next_level: Vec<usize> = Vec::new();
    visited[source.index()] = true;
    let mut depth: usize = 0;

    while !current_level.is_empty() {
        depth += 1;
        next_level.clear();

        for &current_idx in &current_level {
            let current = NodeIndex::new(current_idx);

            for neighbor in {
                let g = &graph.graph;
                g.neighbors_undirected(current)
            } {
                let neighbor_idx = neighbor.index();
                if !visited[neighbor_idx] {
                    if neighbor_idx == target_idx {
                        return Some(depth);
                    }
                    visited[neighbor_idx] = true;
                    next_level.push(neighbor_idx);
                }
            }
        }

        std::mem::swap(&mut current_level, &mut next_level);
    }

    None
}

/// Batch shortest path cost — reuses visited Vec and adjacency list across multiple pairs.
/// Much faster than calling shortest_path_cost N times for large graphs.
pub fn shortest_path_cost_batch(
    graph: &DirGraph,
    pairs: &[(NodeIndex, NodeIndex)],
) -> Vec<Option<usize>> {
    let node_bound = graph.graph.node_bound();

    // Pre-build undirected adjacency list ONCE for all queries
    let nodes: Vec<NodeIndex> = {
        let g = &graph.graph;
        g.node_indices().collect()
    };
    let n = nodes.len();
    let mut node_to_idx = vec![usize::MAX; node_bound];
    for (i, &node) in nodes.iter().enumerate() {
        node_to_idx[node.index()] = i;
    }

    let mut adj: Vec<Vec<usize>> = vec![Vec::new(); n];
    for edge in {
        let g = &graph.graph;
        g.edge_references()
    } {
        let src_i = node_to_idx[edge.source().index()];
        let tgt_i = node_to_idx[edge.target().index()];
        if src_i != usize::MAX && tgt_i != usize::MAX {
            adj[src_i].push(tgt_i);
            adj[tgt_i].push(src_i);
        }
    }
    // Dedup undirected adjacency (handles bidirectional edges A→B + B→A)
    for neighbors in &mut adj {
        neighbors.sort_unstable();
        neighbors.dedup();
    }

    // Reusable visited array — cleared between queries
    let mut visited: Vec<bool> = vec![false; n];
    let mut current_level: Vec<usize> = Vec::new();
    let mut next_level: Vec<usize> = Vec::new();

    let mut results = Vec::with_capacity(pairs.len());

    for &(source, target) in pairs {
        if source == target {
            results.push(Some(0));
            continue;
        }

        let src_i = node_to_idx[source.index()];
        let tgt_i = node_to_idx[target.index()];
        if src_i == usize::MAX || tgt_i == usize::MAX {
            results.push(None);
            continue;
        }

        // Clear visited (only reset nodes we actually touched)
        // Use a generation counter instead of clearing — much faster
        // But for simplicity, track touched nodes
        let mut touched: Vec<usize> = Vec::new();

        current_level.clear();
        current_level.push(src_i);
        visited[src_i] = true;
        touched.push(src_i);
        let mut depth: usize = 0;
        let mut found = false;

        'bfs: while !current_level.is_empty() {
            depth += 1;
            next_level.clear();

            for &current_idx in &current_level {
                for &neighbor_idx in &adj[current_idx] {
                    if !visited[neighbor_idx] {
                        if neighbor_idx == tgt_i {
                            found = true;
                            break 'bfs;
                        }
                        visited[neighbor_idx] = true;
                        touched.push(neighbor_idx);
                        next_level.push(neighbor_idx);
                    }
                }
            }

            std::mem::swap(&mut current_level, &mut next_level);
        }

        results.push(if found { Some(depth) } else { None });

        // Reset only touched nodes (much faster than clearing entire array)
        for &idx in &touched {
            visited[idx] = false;
        }
    }

    results
}

/// Reconstruct path using BFS.
///
/// Phase A.3 / 0.9.53 perf fix: switched from `Vec<bool> + Vec<u32>` of
/// `node_bound` capacity to `HashMap` for parent tracking. The old
/// implementation allocated 500 KB (Vec<bool>) + 2 MB (Vec<u32>) per
/// call on a 500 K-node graph regardless of actual BFS scope. For a
/// shallow lookup that visits ~16 nodes the per-call alloc + init
/// (~30 ms) dominated the operation — see the NornicDB shortestPath
/// benchmark (`scripts/nornicdb_compare.py`), where this fix moved
/// d=1 latency from 37 µs → 4 µs on the 500K-node fixture.
///
/// The HashMap also doubles as the visited set: presence ⇔ visited.
/// Deep BFS does pay slightly more per-node (~100 ns hash insert vs
/// ~5 ns array write), but the per-call alloc savings dominate for
/// any realistic graph + traversal shape; deep paths on small graphs
/// stay fast because the HashMap doesn't preallocate `node_bound`.
fn reconstruct_path_bfs(
    graph: &DirGraph,
    source: NodeIndex,
    target: NodeIndex,
    connection_types: Option<&[InternedKey]>,
    via_types: &Option<HashSet<&str>>,
    deadline: Option<Instant>,
) -> Option<Vec<NodeIndex>> {
    use std::collections::{HashMap, VecDeque};

    if source == target {
        return Some(vec![source]);
    }

    let mut parent: HashMap<usize, u32> = HashMap::with_capacity(64);
    let mut queue: VecDeque<usize> = VecDeque::with_capacity(64);

    let source_idx = source.index();
    let target_idx = target.index();

    parent.insert(source_idx, source_idx as u32);
    queue.push_back(source_idx);

    let mut visit_count = 0u32;

    while let Some(current_idx) = queue.pop_front() {
        // Periodic timeout check (every 1000 nodes)
        visit_count += 1;
        if visit_count.is_multiple_of(1000) {
            if let Some(dl) = deadline {
                if Instant::now() > dl {
                    return None;
                }
            }
        }

        let current = NodeIndex::new(current_idx);

        // Check all neighbors (both directions for undirected path finding)
        let neighbors = filtered_neighbors_undirected(graph, current, connection_types);
        for neighbor in neighbors {
            let neighbor_idx = neighbor.index();

            if parent.contains_key(&neighbor_idx) {
                continue;
            }
            // Apply via_types filter (skip if not target and doesn't match)
            if neighbor_idx != target_idx && !node_passes_via_filter(graph, neighbor, via_types) {
                continue;
            }

            parent.insert(neighbor_idx, current_idx as u32);
            if neighbor_idx == target_idx {
                // Found target - reconstruct path
                let mut path = Vec::with_capacity(16);
                let mut node_idx = target_idx;

                while node_idx != source_idx {
                    path.push(NodeIndex::new(node_idx));
                    node_idx = parent[&node_idx] as usize;
                }
                path.push(source);
                path.reverse();
                return Some(path);
            }
            queue.push_back(neighbor_idx);
        }
    }

    None // No path found
}

/// Directed BFS shortest path — only follows outgoing edges.
/// Used by Cypher shortestPath() which respects edge direction.
///
/// # Arguments
/// * `connection_types` - Only traverse edges of these types (None = all)
/// * `via_types` - Only traverse through nodes of these types (None = all)
pub fn shortest_path_directed(
    graph: &DirGraph,
    source: NodeIndex,
    target: NodeIndex,
    connection_types: Option<&[String]>,
    via_types: Option<&[String]>,
    deadline: Option<Instant>,
) -> Option<PathResult> {
    use std::collections::VecDeque;

    if source == target {
        return Some(PathResult {
            path: vec![source],
            cost: 0,
        });
    }

    let via_set: Option<HashSet<&str>> =
        via_types.map(|vt| vt.iter().map(|s| s.as_str()).collect());
    let interned = intern_connection_types(connection_types);

    let node_bound = graph.graph.node_bound();
    let mut visited: Vec<bool> = vec![false; node_bound];
    let mut parent: Vec<u32> = vec![u32::MAX; node_bound];
    let mut queue = VecDeque::with_capacity(node_bound / 4);

    let source_idx = source.index();
    let target_idx = target.index();

    queue.push_back(source_idx);
    visited[source_idx] = true;

    let mut visit_count = 0u32;

    while let Some(current_idx) = queue.pop_front() {
        // Periodic timeout check
        visit_count += 1;
        if visit_count.is_multiple_of(1000) {
            if let Some(dl) = deadline {
                if Instant::now() > dl {
                    return None;
                }
            }
        }

        let current = NodeIndex::new(current_idx);

        // Only follow outgoing edges
        let neighbors = filtered_neighbors_outgoing(graph, current, interned.as_deref());
        for neighbor in neighbors {
            let neighbor_idx = neighbor.index();

            if !visited[neighbor_idx] {
                // Apply via_types filter (skip if not target and doesn't match)
                if neighbor_idx != target_idx && !node_passes_via_filter(graph, neighbor, &via_set)
                {
                    continue;
                }

                visited[neighbor_idx] = true;
                parent[neighbor_idx] = current_idx as u32;
                queue.push_back(neighbor_idx);

                if neighbor_idx == target_idx {
                    let mut path = Vec::with_capacity(16);
                    let mut node_idx = target_idx;

                    while node_idx != source_idx {
                        path.push(NodeIndex::new(node_idx));
                        node_idx = parent[node_idx] as usize;
                    }
                    path.push(source);
                    path.reverse();

                    let cost = path.len().saturating_sub(1);
                    return Some(PathResult { path, cost });
                }
            }
        }
    }

    None
}

/// Find all paths between two nodes up to a maximum number of hops.
/// Warning: This can be expensive for graphs with many paths!
///
/// # Arguments
/// * `max_results` - Stop after finding this many paths (prevents OOM on dense graphs)
/// * `connection_types` - Only traverse edges of these types (None = all)
/// * `via_types` - Only traverse through nodes of these types (None = all)
#[allow(clippy::too_many_arguments)]
pub fn all_paths(
    graph: &DirGraph,
    source: NodeIndex,
    target: NodeIndex,
    max_hops: usize,
    max_results: Option<usize>,
    connection_types: Option<&[String]>,
    via_types: Option<&[String]>,
    deadline: Option<Instant>,
) -> Vec<Vec<NodeIndex>> {
    let via_set: Option<HashSet<&str>> =
        via_types.map(|vt| vt.iter().map(|s| s.as_str()).collect());
    let interned = intern_connection_types(connection_types);
    let mut results = Vec::new();
    let mut current_path = vec![source];
    let mut visited = HashSet::new();
    visited.insert(source);

    find_all_paths_recursive(
        graph,
        source,
        target,
        max_hops,
        &mut current_path,
        &mut visited,
        &mut results,
        max_results,
        interned.as_deref(),
        &via_set,
        deadline,
    );

    results
}

#[allow(clippy::only_used_in_recursion, clippy::too_many_arguments)]
fn find_all_paths_recursive(
    graph: &DirGraph,
    current: NodeIndex,
    target: NodeIndex,
    remaining_hops: usize,
    current_path: &mut Vec<NodeIndex>,
    visited: &mut HashSet<NodeIndex>,
    results: &mut Vec<Vec<NodeIndex>>,
    max_results: Option<usize>,
    connection_types: Option<&[InternedKey]>,
    via_types: &Option<HashSet<&str>>,
    deadline: Option<Instant>,
) {
    // Early termination when result limit is hit
    if let Some(max) = max_results {
        if results.len() >= max {
            return;
        }
    }

    // Timeout check at each recursive entry
    if let Some(dl) = deadline {
        if Instant::now() > dl {
            return;
        }
    }

    if current == target {
        results.push(current_path.clone());
        return;
    }

    if remaining_hops == 0 {
        return;
    }

    // Explore all neighbors (undirected), filtered by connection type
    let neighbors = filtered_neighbors_undirected(graph, current, connection_types);
    for neighbor in neighbors {
        // Check limit before exploring deeper
        if let Some(max) = max_results {
            if results.len() >= max {
                return;
            }
        }

        if !visited.contains(&neighbor) {
            // Apply via_types filter (skip if not target and doesn't match)
            if neighbor != target && !node_passes_via_filter(graph, neighbor, via_types) {
                continue;
            }

            visited.insert(neighbor);
            current_path.push(neighbor);

            find_all_paths_recursive(
                graph,
                neighbor,
                target,
                remaining_hops - 1,
                current_path,
                visited,
                results,
                max_results,
                connection_types,
                via_types,
                deadline,
            );

            current_path.pop();
            visited.remove(&neighbor);
        }
    }
}

/// Find all strongly connected components in the graph.
/// Returns a vector of components, each component is a vector of node indices.
pub fn connected_components(graph: &DirGraph) -> Vec<Vec<NodeIndex>> {
    // For disk mode, fall back to weakly_connected_components since
    // kosaraju_scc requires petgraph trait bounds.
    if GraphRead::is_disk(&graph.graph) {
        return weakly_connected_components(graph, None)
            .expect("weakly_connected_components with deadline=None cannot time out");
    }
    kosaraju_scc(graph.graph.as_stable_digraph())
}

/// Find weakly connected components (treating graph as undirected).
/// This is often more useful for knowledge graphs.
/// Uses Union-Find (disjoint set) for optimal performance — O(E * α(V)) ≈ O(E).
///
/// @procedure: connected_components
/// @procedure: weakly_connected_components
pub fn weakly_connected_components(
    graph: &DirGraph,
    deadline: Option<Instant>,
) -> Result<Vec<Vec<NodeIndex>>, String> {
    let nodes: Vec<NodeIndex> = {
        let g = &graph.graph;
        g.node_indices().collect()
    };
    let n = nodes.len();

    if n == 0 {
        return Ok(Vec::new());
    }

    // Use node_bound() not node_count() — StableDiGraph indices can have gaps
    let bound = graph.graph.node_bound();

    // Build compact index mapping: graph NodeIndex → contiguous 0..n
    let mut node_to_idx = vec![0usize; bound];
    for (i, &node) in nodes.iter().enumerate() {
        node_to_idx[node.index()] = i;
    }

    // Union-Find with path compression + union by rank
    let mut parent: Vec<usize> = (0..n).collect();
    let mut rank: Vec<u8> = vec![0; n];

    // Find with path compression (iterative)
    #[inline]
    fn find(parent: &mut [usize], mut x: usize) -> usize {
        while parent[x] != x {
            parent[x] = parent[parent[x]]; // path halving
            x = parent[x];
        }
        x
    }

    // Union by rank
    #[inline]
    fn union(parent: &mut [usize], rank: &mut [u8], a: usize, b: usize) {
        let ra = find(parent, a);
        let rb = find(parent, b);
        if ra == rb {
            return;
        }
        if rank[ra] < rank[rb] {
            parent[ra] = rb;
        } else if rank[ra] > rank[rb] {
            parent[rb] = ra;
        } else {
            parent[rb] = ra;
            rank[ra] += 1;
        }
    }

    // Process all edges — single pass, no adjacency list needed.
    // Periodic deadline check (every ~1M edges, negligible overhead via bitmask).
    let mut edge_counter: usize = 0;
    for edge in {
        let g = &graph.graph;
        g.edge_references()
    } {
        edge_counter += 1;
        if edge_counter & 0xFFFFF == 0 {
            if let Some(dl) = deadline {
                if Instant::now() > dl {
                    return Err(algorithm_timeout_err());
                }
            }
        }
        let src_i = node_to_idx[edge.source().index()];
        let tgt_i = node_to_idx[edge.target().index()];
        union(&mut parent, &mut rank, src_i, tgt_i);
    }

    // Collect components by root
    let mut component_map: HashMap<usize, Vec<NodeIndex>> = HashMap::new();
    for (i, &node) in nodes.iter().enumerate() {
        let root = find(&mut parent, i);
        component_map.entry(root).or_default().push(node);
    }

    let mut components: Vec<Vec<NodeIndex>> = component_map.into_values().collect();

    // Sort components by size (largest first)
    components.sort_by_key(|b| std::cmp::Reverse(b.len()));

    Ok(components)
}

/// Get node info for building Python-friendly path output
pub fn get_node_info(graph: &DirGraph, node_idx: NodeIndex) -> Option<PathNodeInfo> {
    let node = graph.get_node(node_idx)?;
    let node_title = node.title();
    let title_str = match &*node_title {
        Value::String(s) => s.clone(),
        _ => format!("{:?}", &*node_title),
    };
    Some(PathNodeInfo {
        node_type: node.node_type_str(&graph.interner).to_string(),
        title: title_str,
        id: node.id().into_owned(),
    })
}

/// Get information about what connection types link nodes in a path
pub fn get_path_connections(graph: &DirGraph, path: &[NodeIndex]) -> Vec<Option<String>> {
    // Pre-allocate with exact size (one connection per edge = path.len() - 1)
    let mut connections = Vec::with_capacity(path.len().saturating_sub(1));

    for window in path.windows(2) {
        let from = window[0];
        let to = window[1];

        // Find edge between these nodes (either direction)
        let conn_type = graph
            .graph
            .edges(from)
            .find(|e| e.target() == to)
            .map(|e| e.weight().connection_type_str(&graph.interner).to_string())
            .or_else(|| {
                graph
                    .graph
                    .edges(to)
                    .find(|e| e.target() == from)
                    .map(|e| e.weight().connection_type_str(&graph.interner).to_string())
            });

        connections.push(conn_type);
    }

    connections
}

/// Check if two nodes are connected (directly or indirectly)
pub fn are_connected(graph: &DirGraph, source: NodeIndex, target: NodeIndex) -> bool {
    shortest_path(graph, source, target, None, None, None).is_some()
}

/// Calculate the degree (number of connections) for a node
pub fn node_degree(graph: &DirGraph, node: NodeIndex) -> usize {
    let g = &graph.graph;
    g.edges(node).count()
        + g.neighbors_directed(node, petgraph::Direction::Incoming)
            .count()
}

// ============================================================================
// Centrality Algorithms
// ============================================================================

/// Result of centrality calculation
#[derive(Debug, Clone)]
pub struct CentralityResult {
    pub node_idx: NodeIndex,
    pub score: f64,
}

/// Calculate betweenness centrality for all nodes in the graph.
///
/// Betweenness centrality measures how often a node lies on the shortest path
/// between other pairs of nodes. Higher values indicate nodes that are more
/// important as "bridges" in the network.
///
/// Uses Brandes' algorithm for efficiency: O(V * E) for unweighted graphs.
/// Optimized to use Vec instead of HashMap for O(1) direct indexing.
///
/// # Arguments
/// * `graph` - The graph to analyze
/// * `normalized` - If true, normalize scores by 2/((n-1)*(n-2)) for directed graphs
/// * `sample_size` - Optional number of source nodes to sample (for large graphs)
///
/// @procedure: betweenness
/// @procedure: betweenness_centrality
pub fn betweenness_centrality(
    graph: &DirGraph,
    normalized: bool,
    sample_size: Option<usize>,
    connection_types: Option<&[String]>,
    deadline: Option<Instant>,
) -> Result<Vec<CentralityResult>, String> {
    use std::collections::VecDeque;
    use std::sync::atomic::{AtomicBool, Ordering};

    let nodes: Vec<NodeIndex> = {
        let g = &graph.graph;
        g.node_indices().collect()
    };
    let n = nodes.len();

    if n <= 2 {
        return Ok(nodes
            .iter()
            .map(|&idx| CentralityResult {
                node_idx: idx,
                score: 0.0,
            })
            .collect());
    }

    // Use Vec-based index mapping for O(1) lookup (vs HashMap)
    let bound = graph.graph.node_bound();
    let mut node_to_idx = vec![0usize; bound];
    for (i, &node) in nodes.iter().enumerate() {
        node_to_idx[node.index()] = i;
    }

    // Pre-build undirected adjacency list for BFS.
    // Betweenness treats edges as undirected so that nodes bridging
    // communities are detected regardless of edge direction.
    let interned_ct = intern_connection_types(connection_types);
    let mut adj: Vec<Vec<usize>> = vec![Vec::new(); n];
    for edge in {
        let g = &graph.graph;
        g.edge_references()
    } {
        if let Some(ref types) = interned_ct {
            if !types.iter().any(|t| *t == edge.weight().connection_type) {
                continue;
            }
        }
        let src_i = node_to_idx[edge.source().index()];
        let tgt_i = node_to_idx[edge.target().index()];
        adj[src_i].push(tgt_i);
        adj[tgt_i].push(src_i);
    }
    // Dedup undirected adjacency (handles bidirectional edges A→B + B→A)
    for neighbors in &mut adj {
        neighbors.sort_unstable();
        neighbors.dedup();
    }

    // Determine which nodes to use as sources
    // Use stride-based sampling to ensure even coverage across the graph,
    // avoiding bias from sequential selection (e.g. first k nodes being
    // Module/Class containers with no outgoing edges of the filtered type).
    let source_indices: Vec<usize> = if let Some(k) = sample_size {
        let k = k.min(n);
        if k == n {
            (0..n).collect()
        } else {
            let step = n as f64 / k as f64;
            (0..k).map(|i| (i as f64 * step) as usize).collect()
        }
    } else {
        (0..n).collect()
    };

    // Parallel vs sequential Brandes' algorithm
    let use_parallel = n >= 4096;

    // Shared timeout flag for the parallel path: rayon closures can't bubble
    // up Result, so we set this on deadline expiry and check after the join.
    let timed_out = AtomicBool::new(false);

    let mut betweenness: Vec<f64> = if use_parallel {
        use rayon::prelude::*;

        let adj_ref = &adj;
        let deadline_ref = &deadline;
        let timed_out_ref = &timed_out;
        let num_threads = rayon::current_num_threads();
        let chunk_size = (source_indices.len() / num_threads).max(1);

        // Thread-local accumulation + reduction (avoids write conflicts on shared array)
        source_indices
            .par_chunks(chunk_size)
            .map(|chunk| {
                // Thread-local data structures (allocated once per thread)
                let mut local_betweenness: Vec<f64> = vec![0.0; n];
                let mut stack: Vec<usize> = Vec::with_capacity(n);
                let mut pred: Vec<Vec<usize>> = vec![Vec::new(); n];
                let mut sigma: Vec<f64> = vec![0.0; n];
                let mut dist: Vec<i64> = vec![-1; n];
                let mut delta: Vec<f64> = vec![0.0; n];
                let mut queue: VecDeque<usize> = VecDeque::with_capacity(n);

                for (local_counter, &s_idx) in chunk.iter().enumerate() {
                    // Periodic timeout check (every 10 sources within this chunk)
                    if local_counter % 10 == 0 {
                        if timed_out_ref.load(Ordering::Relaxed) {
                            break;
                        }
                        if let Some(dl) = deadline_ref {
                            if Instant::now() > *dl {
                                timed_out_ref.store(true, Ordering::Relaxed);
                                break;
                            }
                        }
                    }

                    // Reset only stack/queue
                    stack.clear();
                    queue.clear();

                    // Initialize source
                    sigma[s_idx] = 1.0;
                    dist[s_idx] = 0;
                    queue.push_back(s_idx);

                    // BFS phase
                    while let Some(v_idx) = queue.pop_front() {
                        stack.push(v_idx);
                        let v_dist = dist[v_idx];

                        for &w_idx in &adj_ref[v_idx] {
                            let d = dist[w_idx];
                            if d < 0 {
                                dist[w_idx] = v_dist + 1;
                                queue.push_back(w_idx);
                                sigma[w_idx] += sigma[v_idx];
                                pred[w_idx].push(v_idx);
                            } else if d == v_dist + 1 {
                                sigma[w_idx] += sigma[v_idx];
                                pred[w_idx].push(v_idx);
                            }
                        }
                    }

                    // Accumulation phase + sparse reset
                    while let Some(w_idx) = stack.pop() {
                        for &v_idx in &pred[w_idx] {
                            let contribution = (sigma[v_idx] / sigma[w_idx]) * (1.0 + delta[w_idx]);
                            delta[v_idx] += contribution;
                        }
                        if w_idx != s_idx {
                            local_betweenness[w_idx] += delta[w_idx];
                        }
                        pred[w_idx].clear();
                        sigma[w_idx] = 0.0;
                        dist[w_idx] = -1;
                        delta[w_idx] = 0.0;
                    }
                }

                local_betweenness
            })
            .reduce(
                || vec![0.0; n],
                |mut a, b| {
                    for i in 0..n {
                        a[i] += b[i];
                    }
                    a
                },
            )
    } else {
        // Sequential path (n < 4096): reuses pre-allocated buffers across iterations
        let mut betweenness: Vec<f64> = vec![0.0; n];
        let mut stack: Vec<usize> = Vec::with_capacity(n);
        let mut pred: Vec<Vec<usize>> = vec![Vec::new(); n];
        let mut sigma: Vec<f64> = vec![0.0; n];
        let mut dist: Vec<i64> = vec![-1; n];
        let mut delta: Vec<f64> = vec![0.0; n];
        let mut queue: VecDeque<usize> = VecDeque::with_capacity(n);

        for (source_counter, &s_idx) in source_indices.iter().enumerate() {
            // Periodic timeout check (every 10 source nodes)
            if source_counter.is_multiple_of(10) {
                if let Some(dl) = deadline {
                    if Instant::now() > dl {
                        return Err(algorithm_timeout_err());
                    }
                }
            }

            stack.clear();
            queue.clear();

            sigma[s_idx] = 1.0;
            dist[s_idx] = 0;
            queue.push_back(s_idx);

            while let Some(v_idx) = queue.pop_front() {
                stack.push(v_idx);
                let v_dist = dist[v_idx];

                for &w_idx in &adj[v_idx] {
                    let d = dist[w_idx];
                    if d < 0 {
                        dist[w_idx] = v_dist + 1;
                        queue.push_back(w_idx);
                        sigma[w_idx] += sigma[v_idx];
                        pred[w_idx].push(v_idx);
                    } else if d == v_dist + 1 {
                        sigma[w_idx] += sigma[v_idx];
                        pred[w_idx].push(v_idx);
                    }
                }
            }

            while let Some(w_idx) = stack.pop() {
                for &v_idx in &pred[w_idx] {
                    let contribution = (sigma[v_idx] / sigma[w_idx]) * (1.0 + delta[w_idx]);
                    delta[v_idx] += contribution;
                }
                if w_idx != s_idx {
                    betweenness[w_idx] += delta[w_idx];
                }
                pred[w_idx].clear();
                sigma[w_idx] = 0.0;
                dist[w_idx] = -1;
                delta[w_idx] = 0.0;
            }
        }

        betweenness
    };

    // Surface deadline expiry from the parallel rayon path (set via timed_out flag).
    if timed_out.load(Ordering::Relaxed) {
        return Err(algorithm_timeout_err());
    }

    // Undirected BFS counts each (s,t) pair twice, so halve raw scores.
    for score in betweenness.iter_mut() {
        *score /= 2.0;
    }

    // Normalize if requested
    // For undirected graphs: 2 / ((n-1)*(n-2))
    if normalized && n > 2 {
        let scale = 2.0 / ((n - 1) as f64 * (n - 2) as f64);
        for score in betweenness.iter_mut() {
            *score *= scale;
        }
    }

    // If we sampled, scale up the scores
    if let Some(k) = sample_size {
        if k < n {
            let scale = n as f64 / k as f64;
            for score in betweenness.iter_mut() {
                *score *= scale;
            }
        }
    }

    // Convert to sorted results
    let mut results: Vec<CentralityResult> = nodes
        .iter()
        .enumerate()
        .map(|(i, &node_idx)| CentralityResult {
            node_idx,
            score: betweenness[i],
        })
        .collect();

    // Sort by score descending
    results.sort_unstable_by(|a, b| {
        b.score
            .partial_cmp(&a.score)
            .unwrap_or(std::cmp::Ordering::Equal)
    });

    Ok(results)
}

/// Calculate PageRank centrality for all nodes in the graph.
///
/// PageRank measures the importance of nodes based on the structure of incoming links.
/// Originally developed by Google for ranking web pages.
///
/// # Arguments
/// * `graph` - The graph to analyze
/// * `damping_factor` - Probability of following a link (typically 0.85)
/// * `max_iterations` - Maximum number of iterations (default: 100)
/// * `tolerance` - Convergence threshold (default: 1e-6)
///
/// @procedure: pagerank
pub fn pagerank(
    graph: &DirGraph,
    damping_factor: f64,
    max_iterations: usize,
    tolerance: f64,
    connection_types: Option<&[String]>,
    deadline: Option<Instant>,
) -> Result<Vec<CentralityResult>, String> {
    let nodes: Vec<NodeIndex> = {
        let g = &graph.graph;
        g.node_indices().collect()
    };
    let n = nodes.len();

    if n == 0 {
        return Ok(Vec::new());
    }

    // Use Vec-based index mapping for O(1) lookup (vs HashMap)
    let bound = graph.graph.node_bound();
    let mut node_to_idx = vec![0usize; bound];
    for (i, &node) in nodes.iter().enumerate() {
        node_to_idx[node.index()] = i;
    }

    // Pre-build reverse adjacency list: for each target j, store list of source indices.
    // Pull-based formulation: each target reads from its in-neighbors independently,
    // enabling rayon parallelization (no write conflicts on new_pr[j]).
    let interned_ct = intern_connection_types(connection_types);
    let mut in_adj: Vec<Vec<usize>> = vec![Vec::new(); n];
    let mut out_degrees: Vec<usize> = vec![0; n];
    for edge in {
        let g = &graph.graph;
        g.edge_references()
    } {
        if let Some(ref types) = interned_ct {
            if !types.iter().any(|t| *t == edge.weight().connection_type) {
                continue;
            }
        }
        let src_i = node_to_idx[edge.source().index()];
        let tgt_i = node_to_idx[edge.target().index()];
        in_adj[tgt_i].push(src_i);
        out_degrees[src_i] += 1;
    }

    // Initialize PageRank scores (uniform distribution)
    let mut pr: Vec<f64> = vec![1.0 / n as f64; n];
    let mut new_pr: Vec<f64> = vec![0.0; n];

    // Precompute inverse out-degree (multiply instead of divide in hot loop)
    let inv_out_degrees: Vec<f64> = out_degrees
        .iter()
        .map(|&d| {
            if d > 0 {
                damping_factor / d as f64
            } else {
                0.0
            }
        })
        .collect();

    // Identify dangling nodes (no outgoing links) — store as bitmask for fast sum
    let is_dangling: Vec<bool> = out_degrees.iter().map(|&d| d == 0).collect();

    let teleport = (1.0 - damping_factor) / n as f64;
    let inv_n = 1.0 / n as f64;
    let use_parallel = n >= 4096;

    // Iterative computation
    for _iteration in 0..max_iterations {
        // Timeout check each iteration — error rather than return partial.
        // Half-converged PageRank scores are misleading; an explicit error
        // tells the caller to extend timeout_ms or scope the graph.
        if let Some(dl) = deadline {
            if Instant::now() > dl {
                return Err(algorithm_timeout_err());
            }
        }

        // Calculate dangling node contribution
        let dangling_sum: f64 = if use_parallel {
            use rayon::prelude::*;
            (0..n)
                .into_par_iter()
                .filter(|&i| is_dangling[i])
                .map(|i| pr[i])
                .sum()
        } else {
            (0..n).filter(|&i| is_dangling[i]).map(|i| pr[i]).sum()
        };
        let base_score = teleport + damping_factor * dangling_sum * inv_n;

        // Pull-based PageRank: each target j computes its own score independently.
        // No write conflicts → parallelizable with rayon.
        if use_parallel {
            use rayon::prelude::*;
            new_pr.par_iter_mut().enumerate().for_each(|(j, score)| {
                let mut s = base_score;
                for &src in &in_adj[j] {
                    s += inv_out_degrees[src] * pr[src];
                }
                *score = s;
            });
        } else {
            for j in 0..n {
                let mut s = base_score;
                for &src in &in_adj[j] {
                    s += inv_out_degrees[src] * pr[src];
                }
                new_pr[j] = s;
            }
        }

        // Check for convergence (L1 norm)
        let diff: f64 = if use_parallel {
            use rayon::prelude::*;
            pr.par_iter()
                .zip(new_pr.par_iter())
                .map(|(a, b)| (a - b).abs())
                .sum()
        } else {
            pr.iter()
                .zip(new_pr.iter())
                .map(|(a, b)| (a - b).abs())
                .sum()
        };

        std::mem::swap(&mut pr, &mut new_pr);

        if diff < tolerance {
            break;
        }
    }

    // Convert to results and sort by score
    let mut results: Vec<CentralityResult> = nodes
        .iter()
        .enumerate()
        .map(|(i, &node_idx)| CentralityResult {
            node_idx,
            score: pr[i],
        })
        .collect();

    results.sort_unstable_by(|a, b| {
        b.score
            .partial_cmp(&a.score)
            .unwrap_or(std::cmp::Ordering::Equal)
    });

    Ok(results)
}

/// Calculate degree centrality for all nodes.
///
/// Simply counts the number of connections each node has.
/// Optionally normalized by (n-1) to get values between 0 and 1.
///
/// @procedure: degree
/// @procedure: degree_centrality
pub fn degree_centrality(
    graph: &DirGraph,
    normalized: bool,
    connection_types: Option<&[String]>,
    deadline: Option<Instant>,
) -> Result<Vec<CentralityResult>, String> {
    let nodes: Vec<NodeIndex> = {
        let g = &graph.graph;
        g.node_indices().collect()
    };
    let n = nodes.len();

    if n == 0 {
        return Ok(Vec::new());
    }

    let scale = if normalized && n > 1 {
        1.0 / (n - 1) as f64
    } else {
        1.0
    };

    // Compute all degrees in a single pass over edges instead of per-node traversal.
    // Periodic deadline check (every ~1M edges) keeps overhead negligible while
    // ensuring 863M-edge scans on Wikidata-scale graphs honor the 20s timeout.
    let interned_ct = intern_connection_types(connection_types);
    let bound = graph.graph.node_bound();
    let mut degrees = vec![0usize; bound];
    let mut edge_counter: usize = 0;
    for edge in {
        let g = &graph.graph;
        g.edge_references()
    } {
        edge_counter += 1;
        if edge_counter & 0xFFFFF == 0 {
            if let Some(dl) = deadline {
                if Instant::now() > dl {
                    return Err(algorithm_timeout_err());
                }
            }
        }
        if let Some(ref types) = interned_ct {
            if !types.iter().any(|t| *t == edge.weight().connection_type) {
                continue;
            }
        }
        degrees[edge.source().index()] += 1; // out-degree
        degrees[edge.target().index()] += 1; // in-degree
    }

    let mut results: Vec<CentralityResult> = nodes
        .iter()
        .map(|&node_idx| CentralityResult {
            node_idx,
            score: degrees[node_idx.index()] as f64 * scale,
        })
        .collect();

    results.sort_unstable_by(|a, b| {
        b.score
            .partial_cmp(&a.score)
            .unwrap_or(std::cmp::Ordering::Equal)
    });

    Ok(results)
}

/// Calculate closeness centrality for all nodes.
///
/// Closeness centrality measures how close a node is to all other nodes.
/// Defined as the reciprocal of the sum of shortest path distances.
///
/// Note: For disconnected graphs, only reachable nodes are considered.
/// Optimized to use Vec instead of HashMap for O(1) direct indexing.
///
/// * `sample_size` - Optional number of source nodes to sample (for large graphs).
///   Uses stride-based selection for even coverage.
///
/// @procedure: closeness
/// @procedure: closeness_centrality
pub fn closeness_centrality(
    graph: &DirGraph,
    normalized: bool,
    sample_size: Option<usize>,
    connection_types: Option<&[String]>,
    deadline: Option<Instant>,
) -> Result<Vec<CentralityResult>, String> {
    use std::sync::atomic::{AtomicBool, Ordering};

    let nodes: Vec<NodeIndex> = {
        let g = &graph.graph;
        g.node_indices().collect()
    };
    let n = nodes.len();

    if n == 0 {
        return Ok(Vec::new());
    }

    // Use Vec-based index mapping for O(1) lookup (vs HashMap)
    let bound = graph.graph.node_bound();
    let mut node_to_idx = vec![0usize; bound];
    for (i, &node) in nodes.iter().enumerate() {
        node_to_idx[node.index()] = i;
    }

    // Pre-build incoming adjacency list: for closeness centrality on directed graphs,
    // we BFS via incoming edges (convention: d(v, u) = how easy for v to reach u)
    let interned_ct = intern_connection_types(connection_types);
    let mut adj_incoming: Vec<Vec<usize>> = vec![Vec::new(); n];
    for edge in {
        let g = &graph.graph;
        g.edge_references()
    } {
        if let Some(ref types) = interned_ct {
            if !types.iter().any(|t| *t == edge.weight().connection_type) {
                continue;
            }
        }
        let src_i = node_to_idx[edge.source().index()];
        let tgt_i = node_to_idx[edge.target().index()];
        // For incoming BFS from node u: follow edges pointing INTO u
        // edge: src -> tgt, so tgt has incoming edge from src
        adj_incoming[tgt_i].push(src_i);
    }
    // Dedup incoming adjacency (handles duplicate edges)
    for neighbors in &mut adj_incoming {
        neighbors.sort_unstable();
        neighbors.dedup();
    }

    // Determine which nodes to use as sources.
    // Stride-based sampling ensures even coverage across the graph.
    let source_indices: Vec<usize> = if let Some(k) = sample_size {
        let k = k.min(n);
        if k == n {
            (0..n).collect()
        } else {
            let step = n as f64 / k as f64;
            (0..k).map(|i| (i as f64 * step) as usize).collect()
        }
    } else {
        (0..n).collect()
    };

    // Parallel path: each source BFS is independent, no shared accumulator
    let use_parallel = source_indices.len() >= 4096;

    // Shared timeout flag for the parallel rayon path; checked after the join.
    let timed_out = AtomicBool::new(false);

    if use_parallel {
        use rayon::prelude::*;

        let adj_ref = &adj_incoming;
        let deadline_ref = &deadline;
        let nodes_ref = &nodes;
        let timed_out_ref = &timed_out;

        let mut results: Vec<CentralityResult> = source_indices
            .par_iter()
            .enumerate()
            .map(|(i, &s_idx)| {
                let source = nodes_ref[s_idx];

                // Periodic timeout check (every 100 sources)
                if i % 100 == 0 {
                    if timed_out_ref.load(Ordering::Relaxed) {
                        return CentralityResult {
                            node_idx: source,
                            score: 0.0,
                        };
                    }
                    if let Some(dl) = deadline_ref {
                        if Instant::now() > *dl {
                            timed_out_ref.store(true, Ordering::Relaxed);
                            return CentralityResult {
                                node_idx: source,
                                score: 0.0,
                            };
                        }
                    }
                }

                // Thread-local BFS data structures
                let mut dist: Vec<i64> = vec![-1; n];
                let mut current_level: Vec<usize> = Vec::with_capacity(n / 4);
                let mut next_level: Vec<usize> = Vec::with_capacity(n / 4);
                let mut touched: Vec<usize> = Vec::with_capacity(n / 4);

                // Initialize source
                current_level.push(s_idx);
                dist[s_idx] = 0;
                touched.push(s_idx);
                let mut depth: i64 = 0;

                // Level-based BFS via incoming edges
                while !current_level.is_empty() {
                    depth += 1;
                    next_level.clear();

                    for &current_idx in &current_level {
                        for &neighbor_idx in &adj_ref[current_idx] {
                            if dist[neighbor_idx] < 0 {
                                dist[neighbor_idx] = depth;
                                next_level.push(neighbor_idx);
                                touched.push(neighbor_idx);
                            }
                        }
                    }

                    std::mem::swap(&mut current_level, &mut next_level);
                }

                // Calculate closeness from touched nodes only
                let reachable = touched.len();
                let total_distance: i64 = touched.iter().map(|&idx| dist[idx]).sum();

                if reachable > 1 && total_distance > 0 {
                    let closeness = (reachable - 1) as f64 / total_distance as f64;
                    let score = if normalized {
                        closeness * (reachable - 1) as f64 / (n - 1) as f64
                    } else {
                        closeness
                    };
                    CentralityResult {
                        node_idx: source,
                        score,
                    }
                } else {
                    CentralityResult {
                        node_idx: source,
                        score: 0.0,
                    }
                }
            })
            .collect();

        if timed_out.load(Ordering::Relaxed) {
            return Err(algorithm_timeout_err());
        }

        results.sort_unstable_by(|a, b| {
            b.score
                .partial_cmp(&a.score)
                .unwrap_or(std::cmp::Ordering::Equal)
        });

        return Ok(results);
    }

    // Sequential path: reuses pre-allocated buffers across iterations
    let mut results = Vec::with_capacity(source_indices.len());
    let mut dist: Vec<i64> = vec![-1; n];
    let mut current_level: Vec<usize> = Vec::with_capacity(n);
    let mut next_level: Vec<usize> = Vec::with_capacity(n);
    let mut touched: Vec<usize> = Vec::with_capacity(n);

    for (i, &s_idx) in source_indices.iter().enumerate() {
        let source = nodes[s_idx];

        // Periodic timeout check (every 10 source nodes)
        if i.is_multiple_of(10) {
            if let Some(dl) = deadline {
                if Instant::now() > dl {
                    return Err(algorithm_timeout_err());
                }
            }
        }

        // Sparse reset from previous iteration (only visited nodes)
        for &idx in &touched {
            dist[idx] = -1;
        }
        touched.clear();
        current_level.clear();

        // Initialize source
        current_level.push(s_idx);
        dist[s_idx] = 0;
        touched.push(s_idx);
        let mut depth: i64 = 0;

        // Level-based BFS via incoming edges
        while !current_level.is_empty() {
            depth += 1;
            next_level.clear();

            for &current_idx in &current_level {
                for &neighbor_idx in &adj_incoming[current_idx] {
                    if dist[neighbor_idx] < 0 {
                        dist[neighbor_idx] = depth;
                        next_level.push(neighbor_idx);
                        touched.push(neighbor_idx);
                    }
                }
            }

            std::mem::swap(&mut current_level, &mut next_level);
        }

        // Calculate closeness from touched nodes only (not all N)
        let reachable = touched.len();
        let total_distance: i64 = touched.iter().map(|&idx| dist[idx]).sum();

        if reachable > 1 && total_distance > 0 {
            let closeness = (reachable - 1) as f64 / total_distance as f64;

            let score = if normalized {
                closeness * (reachable - 1) as f64 / (n - 1) as f64
            } else {
                closeness
            };

            results.push(CentralityResult {
                node_idx: source,
                score,
            });
        } else {
            results.push(CentralityResult {
                node_idx: source,
                score: 0.0,
            });
        }
    }

    results.sort_unstable_by(|a, b| {
        b.score
            .partial_cmp(&a.score)
            .unwrap_or(std::cmp::Ordering::Equal)
    });

    Ok(results)
}

// ============================================================================
// Community Detection
// ============================================================================

#[derive(Debug, Clone)]
pub struct CommunityAssignment {
    pub node_idx: NodeIndex,
    pub community_id: usize,
}

#[derive(Debug)]
pub struct CommunityResult {
    pub assignments: Vec<CommunityAssignment>,
    pub num_communities: usize,
    pub modularity: f64,
}

/// Louvain modularity optimization for community detection.
///
/// Each node starts in its own community. Iteratively moves nodes to the
/// neighboring community that yields the largest modularity gain, until
/// no improvement is found.
/// Optimized with pre-built adjacency list and Vec-based community weight tracking.
///
/// @procedure: louvain
/// @procedure: louvain_communities
pub fn louvain_communities(
    graph: &DirGraph,
    weight_property: Option<&str>,
    resolution: f64,
    connection_types: Option<&[String]>,
    deadline: Option<Instant>,
) -> Result<CommunityResult, String> {
    let nodes: Vec<NodeIndex> = {
        let g = &graph.graph;
        g.node_indices().collect()
    };
    let n = nodes.len();

    if n == 0 {
        return Ok(CommunityResult {
            assignments: Vec::new(),
            num_communities: 0,
            modularity: 0.0,
        });
    }

    // Build compact index mapping
    let bound = graph.graph.node_bound();
    let mut node_to_idx = vec![0usize; bound];
    for (i, &node) in nodes.iter().enumerate() {
        node_to_idx[node.index()] = i;
    }

    // Pre-build undirected weighted adjacency list
    // adj[i] = Vec<(neighbor_compact_idx, weight)>
    let interned_ct = intern_connection_types(connection_types);
    let mut adj: Vec<Vec<(usize, f64)>> = vec![Vec::new(); n];
    let mut total_weight = 0.0f64;
    for edge in {
        let g = &graph.graph;
        g.edge_references()
    } {
        if let Some(ref types) = interned_ct {
            if !types.iter().any(|t| *t == edge.weight().connection_type) {
                continue;
            }
        }
        let w = edge_weight(graph, edge.id(), weight_property);
        let src_i = node_to_idx[edge.source().index()];
        let tgt_i = node_to_idx[edge.target().index()];
        adj[src_i].push((tgt_i, w));
        adj[tgt_i].push((src_i, w));
        total_weight += w;
    }
    // Dedup weighted adjacency: merge duplicate neighbors by summing weights
    for neighbors in &mut adj {
        neighbors.sort_unstable_by_key(|&(idx, _)| idx);
        neighbors.dedup_by(|a, b| {
            if a.0 == b.0 {
                b.1 += a.1;
                true
            } else {
                false
            }
        });
    }

    if total_weight == 0.0 {
        // No edges — each node is its own community
        let assignments: Vec<CommunityAssignment> = nodes
            .iter()
            .enumerate()
            .map(|(i, &idx)| CommunityAssignment {
                node_idx: idx,
                community_id: i,
            })
            .collect();
        return Ok(CommunityResult {
            assignments,
            num_communities: n,
            modularity: 0.0,
        });
    }

    // community[i] = community id for compact node i
    let mut community: Vec<usize> = (0..n).collect();

    // Precompute node degrees (undirected: sum of all edge weights touching the node)
    let mut degree: Vec<f64> = vec![0.0; n];
    for i in 0..n {
        for &(_, w) in &adj[i] {
            degree[i] += w;
        }
    }

    // Precompute sum of degrees per community (sigma_tot)
    let mut sigma_tot: Vec<f64> = vec![0.0; n];
    sigma_tot[..n].copy_from_slice(&degree[..n]);

    let m = total_weight;
    let two_m = 2.0 * m;
    // Precompute loop-invariant division terms as multipliers
    let inv_m = 1.0 / m;
    let resolution_over_two_m_sq = resolution / (two_m * two_m);

    // Vec-based community weight tracking (reused across iterations)
    let mut comm_weight: Vec<f64> = vec![0.0; n];
    let mut touched_comms: Vec<usize> = Vec::with_capacity(64);

    // Iterative optimization
    let max_iterations = 100;
    for _ in 0..max_iterations {
        // Timeout check each iteration — error rather than return partial.
        // Half-converged communities are misleading; an explicit error tells
        // the caller to extend timeout_ms or scope the graph.
        if let Some(dl) = deadline {
            if Instant::now() > dl {
                return Err(algorithm_timeout_err());
            }
        }

        let mut improved = false;

        for i in 0..n {
            let current_community = community[i];
            let k_i = degree[i];
            let k_i_res = k_i * resolution_over_two_m_sq; // precomputed per node

            // Compute weight from node i to each neighboring community
            touched_comms.clear();
            for &(neighbor, w) in &adj[i] {
                let c = community[neighbor];
                if comm_weight[c] == 0.0 {
                    touched_comms.push(c);
                }
                comm_weight[c] += w;
            }

            // Weight from i into its own community
            let k_i_in_current = comm_weight[current_community];

            // Find best community to move to
            let mut best_community = current_community;
            let mut best_delta = 0.0f64;

            for &cand_community in &touched_comms {
                if cand_community == current_community {
                    continue;
                }

                let k_i_in_cand = comm_weight[cand_community];
                let sigma_cand = sigma_tot[cand_community];
                let sigma_curr = sigma_tot[current_community] - k_i;

                let gain_add = k_i_in_cand * inv_m - sigma_cand * k_i_res;
                let loss_remove = k_i_in_current * inv_m - sigma_curr * k_i_res;
                let delta = gain_add - loss_remove;

                if delta > best_delta {
                    best_delta = delta;
                    best_community = cand_community;
                }
            }

            // Reset community weights (only touched entries)
            for &c in &touched_comms {
                comm_weight[c] = 0.0;
            }

            if best_community != current_community {
                sigma_tot[current_community] -= k_i;
                sigma_tot[best_community] += k_i;
                community[i] = best_community;
                improved = true;
            }
        }

        if !improved {
            break;
        }
    }

    // Convert compact community to bound-sized array for modularity computation
    let mut community_bound: Vec<usize> = vec![0; bound];
    let mut node_exists: Vec<bool> = vec![false; bound];
    for (i, &node) in nodes.iter().enumerate() {
        community_bound[node.index()] = community[i];
        node_exists[node.index()] = true;
    }

    // Renumber communities to be contiguous 0..n
    let mut id_map: HashMap<usize, usize> = HashMap::new();
    for &c in &community {
        let next_id = id_map.len();
        id_map.entry(c).or_insert(next_id);
    }

    let assignments: Vec<CommunityAssignment> = nodes
        .iter()
        .enumerate()
        .map(|(i, &idx)| CommunityAssignment {
            node_idx: idx,
            community_id: *id_map.get(&community[i]).unwrap(),
        })
        .collect();

    let num_communities = id_map.len();
    let modularity = compute_modularity(
        graph,
        &community_bound,
        &node_exists,
        total_weight,
        weight_property,
    );

    Ok(CommunityResult {
        assignments,
        num_communities,
        modularity,
    })
}

/// Label propagation for community detection.
///
/// Each node adopts the most frequent label among its neighbors.
/// Converges when no node changes its label.
/// Optimized with pre-built adjacency list and Vec-based label counting.
///
/// @procedure: label_propagation
pub fn label_propagation(
    graph: &DirGraph,
    max_iterations: usize,
    connection_types: Option<&[String]>,
    deadline: Option<Instant>,
) -> Result<CommunityResult, String> {
    let nodes: Vec<NodeIndex> = {
        let g = &graph.graph;
        g.node_indices().collect()
    };
    let n = nodes.len();

    if n == 0 {
        return Ok(CommunityResult {
            assignments: Vec::new(),
            num_communities: 0,
            modularity: 0.0,
        });
    }

    // Build compact index mapping
    let bound = graph.graph.node_bound();
    let mut node_to_idx = vec![0usize; bound];
    for (i, &node) in nodes.iter().enumerate() {
        node_to_idx[node.index()] = i;
    }

    // Pre-build undirected adjacency list (both directions)
    let interned_ct = intern_connection_types(connection_types);
    let mut adj: Vec<Vec<usize>> = vec![Vec::new(); n];
    for edge in {
        let g = &graph.graph;
        g.edge_references()
    } {
        if let Some(ref types) = interned_ct {
            if !types.iter().any(|t| *t == edge.weight().connection_type) {
                continue;
            }
        }
        let src_i = node_to_idx[edge.source().index()];
        let tgt_i = node_to_idx[edge.target().index()];
        adj[src_i].push(tgt_i);
        adj[tgt_i].push(src_i);
    }
    // Dedup undirected adjacency (handles bidirectional edges A→B + B→A)
    for neighbors in &mut adj {
        neighbors.sort_unstable();
        neighbors.dedup();
    }

    // Initialize: each node gets a unique label (0..n)
    let mut labels: Vec<usize> = (0..n).collect();

    // Vec-based label counting (reused across iterations)
    // label_count[label] = count for that label among neighbors
    // We use a sparse approach: track which labels were touched and reset only those
    let mut label_count: Vec<usize> = vec![0; n];
    let mut touched_labels: Vec<usize> = Vec::with_capacity(64);

    for _ in 0..max_iterations {
        // Timeout check each iteration — error rather than return partial.
        if let Some(dl) = deadline {
            if Instant::now() > dl {
                return Err(algorithm_timeout_err());
            }
        }

        let mut changed = false;

        for i in 0..n {
            let neighbors = &adj[i];
            if neighbors.is_empty() {
                continue; // isolated node keeps its label
            }

            // Count neighbor labels using Vec (O(1) per access)
            touched_labels.clear();
            for &neighbor in neighbors {
                let lbl = labels[neighbor];
                if label_count[lbl] == 0 {
                    touched_labels.push(lbl);
                }
                label_count[lbl] += 1;
            }

            // Find most frequent label
            let mut best_label = labels[i];
            let mut best_count = 0;
            for &lbl in &touched_labels {
                if label_count[lbl] > best_count {
                    best_count = label_count[lbl];
                    best_label = lbl;
                }
            }

            // Reset counts for next node (only touched entries)
            for &lbl in &touched_labels {
                label_count[lbl] = 0;
            }

            if best_label != labels[i] {
                labels[i] = best_label;
                changed = true;
            }
        }

        if !changed {
            break;
        }
    }

    // Convert compact labels back to bound-sized array for modularity computation
    let mut labels_bound: Vec<usize> = vec![0; bound];
    let mut node_exists: Vec<bool> = vec![false; bound];
    for (i, &node) in nodes.iter().enumerate() {
        labels_bound[node.index()] = labels[i];
        node_exists[node.index()] = true;
    }

    // Renumber labels to be contiguous
    let mut id_map: HashMap<usize, usize> = HashMap::new();
    for &lbl in &labels {
        let next_id = id_map.len();
        id_map.entry(lbl).or_insert(next_id);
    }

    let assignments: Vec<CommunityAssignment> = nodes
        .iter()
        .enumerate()
        .map(|(i, &idx)| CommunityAssignment {
            node_idx: idx,
            community_id: *id_map.get(&labels[i]).unwrap(),
        })
        .collect();

    let total_weight = graph.graph.edge_count() as f64;
    let num_communities = id_map.len();
    let modularity = compute_modularity(graph, &labels_bound, &node_exists, total_weight, None);

    Ok(CommunityResult {
        assignments,
        num_communities,
        modularity,
    })
}

/// Get edge weight from a property, or 1.0 if not specified.
pub(crate) fn edge_weight(
    graph: &DirGraph,
    edge_id: petgraph::graph::EdgeIndex,
    weight_property: Option<&str>,
) -> f64 {
    if let Some(prop) = weight_property {
        let g = &graph.graph;
        if let Some(edge_data) = g.edge_weight(edge_id) {
            if let Some(val) = edge_data.get_property(prop) {
                return crate::graph::core::value_operations::value_to_f64(val).unwrap_or(1.0);
            }
        }
    }
    1.0
}

/// Result of a weighted path-finding operation. Distinct from [`PathResult`]
/// (which carries an integer hop count) so the f64 weight survives a round
/// trip through the Python layer without coercion.
#[derive(Debug, Clone)]
pub struct WeightedPathResult {
    pub path: Vec<NodeIndex>,
    pub weight: f64,
}

/// Dijkstra-based weighted shortest path. Treats the graph as undirected
/// (mirrors [`shortest_path`]) and reads `weight_property` from each edge,
/// defaulting to 1.0 when the property is absent or non-numeric — the same
/// fallback Louvain uses for its weighted-adjacency build.
///
/// Returns `None` when no path exists, when the deadline expires, or when
/// any traversed edge has a negative weight (Dijkstra requires non-negative
/// edges; the procedure errs on the side of returning no path rather than
/// silently producing wrong answers).
pub fn shortest_path_weighted(
    graph: &DirGraph,
    source: NodeIndex,
    target: NodeIndex,
    weight_property: &str,
    connection_types: Option<&[String]>,
    via_types: Option<&[String]>,
    deadline: Option<Instant>,
) -> Option<WeightedPathResult> {
    use std::cmp::Ordering;
    use std::collections::BinaryHeap;

    if source == target {
        return Some(WeightedPathResult {
            path: vec![source],
            weight: 0.0,
        });
    }

    let via_set: Option<HashSet<&str>> =
        via_types.map(|vt| vt.iter().map(|s| s.as_str()).collect());
    let interned = intern_connection_types(connection_types);
    let conn_filter = interned.as_deref();

    let node_bound = graph.graph.node_bound();
    let mut dist: Vec<f64> = vec![f64::INFINITY; node_bound];
    let mut parent: Vec<u32> = vec![u32::MAX; node_bound];
    dist[source.index()] = 0.0;

    // Min-heap keyed on (distance, node_index). f64 isn't Ord, so wrap in a
    // newtype that flips the order to make BinaryHeap a min-heap.
    #[derive(PartialEq)]
    struct State(f64, usize);
    impl Eq for State {}
    impl Ord for State {
        fn cmp(&self, other: &Self) -> Ordering {
            other
                .0
                .partial_cmp(&self.0)
                .unwrap_or(Ordering::Equal)
                .then_with(|| self.1.cmp(&other.1))
        }
    }
    impl PartialOrd for State {
        fn partial_cmp(&self, other: &Self) -> Option<Ordering> {
            Some(self.cmp(other))
        }
    }

    let mut heap: BinaryHeap<State> = BinaryHeap::new();
    heap.push(State(0.0, source.index()));

    let mut visit_count = 0u32;
    while let Some(State(d, current_idx)) = heap.pop() {
        // Stale entry — already processed via a shorter path.
        if d > dist[current_idx] {
            continue;
        }
        if current_idx == target.index() {
            // Reconstruct path
            let mut path = Vec::with_capacity(16);
            let mut idx = current_idx;
            while idx != source.index() {
                path.push(NodeIndex::new(idx));
                idx = parent[idx] as usize;
            }
            path.push(source);
            path.reverse();
            return Some(WeightedPathResult { path, weight: d });
        }

        visit_count += 1;
        if visit_count.is_multiple_of(1000) {
            if let Some(dl) = deadline {
                if Instant::now() > dl {
                    return None;
                }
            }
        }

        let current = NodeIndex::new(current_idx);
        for edge in graph
            .graph
            .edges_directed(current, petgraph::Direction::Outgoing)
            .chain(
                graph
                    .graph
                    .edges_directed(current, petgraph::Direction::Incoming),
            )
        {
            if let Some(types) = conn_filter {
                if !types.iter().any(|t| *t == edge.weight().connection_type) {
                    continue;
                }
            }
            let neighbor = if edge.source() == current {
                edge.target()
            } else {
                edge.source()
            };
            let n_idx = neighbor.index();
            if n_idx != target.index() && !node_passes_via_filter(graph, neighbor, &via_set) {
                continue;
            }
            let w = edge_weight(graph, edge.id(), Some(weight_property));
            if w < 0.0 {
                // Dijkstra is invalid with negative weights — abort.
                return None;
            }
            let next = d + w;
            if next < dist[n_idx] {
                dist[n_idx] = next;
                parent[n_idx] = current_idx as u32;
                heap.push(State(next, n_idx));
            }
        }
    }
    None
}

/// Lightweight variant of [`shortest_path_weighted`] that returns only the
/// total weight without reconstructing the path.
pub fn shortest_path_cost_weighted(
    graph: &DirGraph,
    source: NodeIndex,
    target: NodeIndex,
    weight_property: &str,
    connection_types: Option<&[String]>,
    via_types: Option<&[String]>,
    deadline: Option<Instant>,
) -> Option<f64> {
    shortest_path_weighted(
        graph,
        source,
        target,
        weight_property,
        connection_types,
        via_types,
        deadline,
    )
    .map(|r| r.weight)
}

/// Sum of edge weights for all nodes in a community.
/// Compute Newman modularity: Q = (1/2m) * sum [ A_ij - k_i*k_j/(2m) ] * delta(c_i, c_j)
fn compute_modularity(
    graph: &DirGraph,
    community: &[usize],
    node_exists: &[bool],
    total_weight: f64,
    weight_property: Option<&str>,
) -> f64 {
    if total_weight == 0.0 {
        return 0.0;
    }

    let two_m = 2.0 * total_weight;
    let mut q = 0.0f64;

    // Compute degree (sum of edge weights) for each node
    let g = &graph.graph;
    let bound = g.node_bound();
    let mut degrees: Vec<f64> = vec![0.0; bound];
    for node_idx in g.node_indices() {
        let i = node_idx.index();
        if !node_exists[i] {
            continue;
        }
        for edge in g.edges(node_idx) {
            degrees[i] += edge_weight(graph, edge.id(), weight_property);
        }
        for edge in g.edges_directed(node_idx, petgraph::Direction::Incoming) {
            degrees[i] += edge_weight(graph, edge.id(), weight_property);
        }
    }

    // Sum over all edges
    for edge in {
        let g = &graph.graph;
        g.edge_references()
    } {
        let u = edge.source().index();
        let v = edge.target().index();
        let w = edge_weight(graph, edge.id(), weight_property);

        if community[u] == community[v] {
            q += w - degrees[u] * degrees[v] / two_m;
        }
    }

    q / two_m
}

#[cfg(test)]
#[path = "graph_algorithms_tests.rs"]
mod tests;