graphina 0.4.0-alpha.2

/*!
# Graph-level Metrics

Graph-level metrics and statistics.
*/

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

use crate::core::types::{BaseGraph, GraphConstructor, NodeId};
use petgraph::EdgeType;

/// Orders a pair of node indices as `(low, high)`, the canonical key form for the
/// undirected edge set used by triangle-counting metrics.
fn order_pair(a: usize, b: usize) -> (usize, usize) {
    if a <= b { (a, b) } else { (b, a) }
}

/// Computes the diameter of the graph (longest shortest path).
///
/// For disconnected graphs, returns None.
///
/// # Time Complexity
/// O(V * (V + E)) - Runs BFS from each node
///
/// # Example
///
/// ```rust
/// use graphina::core::types::Graph;
/// use graphina::metrics::diameter;
///
/// let mut g = Graph::<i32, f64>::new();
/// let n1 = g.add_node(1);
/// let n2 = g.add_node(2);
/// let n3 = g.add_node(3);
/// g.add_edge(n1, n2, 1.0);
/// g.add_edge(n2, n3, 1.0);
///
/// assert_eq!(diameter(&g), Some(2));
/// ```
pub fn diameter<A, W, Ty: GraphConstructor<A, W> + EdgeType>(
    graph: &BaseGraph<A, W, Ty>,
) -> Option<usize> {
    if graph.is_empty() {
        return None;
    }

    let mut max_distance = 0;

    for start_node in graph.node_ids() {
        let distances = bfs_distances(graph, start_node);

        // If any node is unreachable, graph is disconnected
        if distances.len() != graph.node_count() {
            return None;
        }

        if let Some(&max_dist) = distances.values().max() {
            if max_dist > max_distance {
                max_distance = max_dist;
            }
        }
    }

    Some(max_distance)
}

/// Computes the radius of the graph (minimum eccentricity).
///
/// The radius is the minimum over all nodes of their eccentricity
/// (maximum distance to any other node).
///
/// # Time Complexity
/// O(V * (V + E))
pub fn radius<A, W, Ty: GraphConstructor<A, W> + EdgeType>(
    graph: &BaseGraph<A, W, Ty>,
) -> Option<usize> {
    if graph.is_empty() {
        return None;
    }

    let mut min_eccentricity = usize::MAX;

    for start_node in graph.node_ids() {
        let distances = bfs_distances(graph, start_node);

        // If any node is unreachable, graph is disconnected
        if distances.len() != graph.node_count() {
            return None;
        }

        if let Some(&max_dist) = distances.values().max() {
            if max_dist < min_eccentricity {
                min_eccentricity = max_dist;
            }
        }
    }

    Some(min_eccentricity)
}

/// Computes the average clustering coefficient of the graph.
///
/// The clustering coefficient measures the degree to which nodes tend to
/// cluster together. Returns the average of all local clustering coefficients.
///
/// # Time Complexity
/// O(V * d²) where d is average degree
///
/// # Example
///
/// ```rust
/// use graphina::core::types::Graph;
/// use graphina::metrics::average_clustering_coefficient;
///
/// let mut g = Graph::<i32, f64>::new();
/// let n1 = g.add_node(1);
/// let n2 = g.add_node(2);
/// let n3 = g.add_node(3);
/// g.add_edge(n1, n2, 1.0);
/// g.add_edge(n2, n3, 1.0);
/// g.add_edge(n3, n1, 1.0); // Triangle
///
/// assert!((average_clustering_coefficient(&g) - 1.0).abs() < 0.001);
/// ```
pub fn average_clustering_coefficient<A, W, Ty: GraphConstructor<A, W> + EdgeType>(
    graph: &BaseGraph<A, W, Ty>,
) -> f64 {
    if graph.is_empty() {
        return 0.0;
    }

    let coefficients: Vec<f64> = graph
        .node_ids()
        .map(|node| super::node_metrics::clustering_coefficient(graph, node))
        .collect();

    coefficients.iter().sum::<f64>() / coefficients.len() as f64
}

/// Computes the transitivity (global clustering coefficient) of the graph.
///
/// Measures the ratio of triangles to connected triples in the graph.
///
/// # Time Complexity
/// O(E^1.5) via degree-ordered forward triangle counting, down from O(V * d²).
pub fn transitivity<A, W, Ty: GraphConstructor<A, W> + EdgeType>(
    graph: &BaseGraph<A, W, Ty>,
) -> f64 {
    use petgraph::visit::NodeIndexable;

    // Adjacency test via a single Fx-hashed set of canonical (low, high) endpoint
    // pairs, so an edge lookup is one integer-tuple hash rather than an O(degree)
    // `contains_edge` call.
    let mut edge_set: HashSet<(usize, usize), rustc_hash::FxBuildHasher> =
        HashSet::with_capacity_and_hasher(graph.edge_count(), rustc_hash::FxBuildHasher);
    for (u, v, _w) in graph.edges() {
        let (lo, hi) = order_pair(u.index(), v.index());
        edge_set.insert((lo, hi));
    }

    // Rank each node by its neighbor count, keyed by the stable node index so ranks
    // are dense and comparable in O(1). A node's rank is the tuple (degree, index):
    // higher degree ranks higher, ties broken by index.
    let bound = graph.as_petgraph().node_bound();
    let mut degree = vec![0usize; bound];
    for node in graph.node_ids() {
        degree[node.index()] = graph.neighbors(node).count();
    }
    let higher_rank =
        |a: usize, v: usize| degree[a] > degree[v] || (degree[a] == degree[v] && a > v);

    // Connected triples are C(degree, 2) summed over all nodes, independent of the
    // triangle count.
    let mut triples = 0usize;
    for &k in &degree {
        if k >= 2 {
            triples += k * (k - 1) / 2;
        }
    }
    if triples == 0 {
        return 0.0;
    }

    // Count each triangle exactly once at its lowest-ranked vertex: from that vertex
    // both other members are higher ranked, so enumerating pairs among the
    // higher-ranked neighbors alone finds the triangle once and only once. On a
    // graph with n nodes each of degree d this examines about half the neighbors per
    // node, so roughly a quarter of the pairs the all-pairs enumeration would.
    let mut higher: Vec<usize> = Vec::new();
    let mut distinct_triangles = 0usize;
    for node in graph.node_ids() {
        let vi = node.index();
        higher.clear();
        higher.extend(
            graph
                .neighbors(node)
                .map(|nbr| nbr.index())
                .filter(|&nbr| higher_rank(nbr, vi)),
        );
        let k = higher.len();
        for i in 0..k {
            let ni = higher[i];
            for &other in higher.iter().skip(i + 1) {
                let (lo, hi) = order_pair(ni, other);
                if edge_set.contains(&(lo, hi)) {
                    distinct_triangles += 1;
                }
            }
        }
    }

    // Each distinct triangle closes three connected triples (one per apex), matching
    // the standard 3 * triangles / triples definition of transitivity.
    (3 * distinct_triangles) as f64 / triples as f64
}

/// Computes the average path length of the graph.
///
/// Returns the average shortest path length between all pairs of nodes.
/// For disconnected graphs, returns None.
///
/// # Time Complexity
/// O(V * (V + E))
pub fn average_path_length<A, W, Ty: GraphConstructor<A, W> + EdgeType>(
    graph: &BaseGraph<A, W, Ty>,
) -> Option<f64> {
    if graph.is_empty() {
        return None;
    }

    let mut total_distance = 0.0;
    let mut pair_count = 0;

    for start_node in graph.node_ids() {
        let distances = bfs_distances(graph, start_node);

        // If any node is unreachable, graph is disconnected
        if distances.len() != graph.node_count() {
            return None;
        }

        for &dist in distances.values() {
            if dist > 0 {
                total_distance += dist as f64;
                pair_count += 1;
            }
        }
    }

    if pair_count == 0 {
        return Some(0.0);
    }

    Some(total_distance / pair_count as f64)
}

/// Computes the assortativity coefficient of the graph.
///
/// Measures the tendency of nodes to connect to others with similar degree.
/// Returns a value between -1 (disassortative) and 1 (assortative).
///
/// # Time Complexity
/// O(E)
pub fn assortativity<A, W, Ty: GraphConstructor<A, W> + EdgeType>(
    graph: &BaseGraph<A, W, Ty>,
) -> f64 {
    if graph.edge_count() == 0 {
        return 0.0;
    }

    let mut sum_jk = 0.0;
    let mut sum_j = 0.0;
    let mut sum_k = 0.0;
    let mut sum_j2 = 0.0;
    let mut sum_k2 = 0.0;
    // Degree assortativity (Newman) is the Pearson correlation over the joint
    // degree distribution of edge endpoints, which is symmetric: each
    // undirected edge contributes both orderings (j, k) and (k, j). Counting a
    // single ordering would give the two endpoints different means and yield a
    // different, direction-dependent coefficient.
    let m = (graph.edge_count() * 2) as f64;

    for (u, v, _) in graph.edges() {
        let j = graph.degree(u).unwrap_or(0) as f64;
        let k = graph.degree(v).unwrap_or(0) as f64;

        sum_jk += 2.0 * j * k;
        sum_j += j + k;
        sum_k += j + k;
        sum_j2 += j * j + k * k;
        sum_k2 += j * j + k * k;
    }

    let numerator = sum_jk / m - (sum_j / m) * (sum_k / m);
    let denominator =
        ((sum_j2 / m - (sum_j / m).powi(2)) * (sum_k2 / m - (sum_k / m).powi(2))).sqrt();

    if denominator == 0.0 {
        return 0.0;
    }

    numerator / denominator
}

/// Helper function: Computes BFS distances from a start node.
fn bfs_distances<A, W, Ty: GraphConstructor<A, W> + EdgeType>(
    graph: &BaseGraph<A, W, Ty>,
    start: NodeId,
) -> HashMap<NodeId, usize> {
    let mut distances = HashMap::new();
    let mut queue = VecDeque::new();

    distances.insert(start, 0);
    queue.push_back(start);

    while let Some(node) = queue.pop_front() {
        let dist = distances[&node];

        for neighbor in graph.neighbors(node) {
            if let std::collections::hash_map::Entry::Vacant(e) = distances.entry(neighbor) {
                e.insert(dist + 1);
                queue.push_back(neighbor);
            }
        }
    }

    distances
}

#[cfg(test)]
mod tests {

    #[test]
    fn test_assortativity_is_symmetric_newman_coefficient() {
        use crate::core::types::Graph;
        // Degree assortativity is the Pearson correlation over the symmetric joint
        // degree distribution of edge endpoints: each undirected edge contributes
        // both orderings. A star is perfectly disassortative (a degree-k hub joined
        // only to degree-1 leaves), so its coefficient is exactly -1.0. The earlier
        // implementation correlated a single edge ordering, giving the two endpoints
        // different means; on a star that collapsed the variance to zero and
        // returned 0.0 instead of -1.0.
        use crate::metrics::assortativity;

        let mut g: Graph<i32, f64> = Graph::new();
        let center = g.add_node(0);
        for i in 1..=4 {
            let leaf = g.add_node(i);
            g.add_edge(center, leaf, 1.0);
        }

        let r = assortativity(&g);
        assert!(
            (r - (-1.0)).abs() < 1e-9,
            "star degree assortativity should be -1.0, got {r}"
        );
    }
    use super::*;
    use crate::core::types::Graph;

    #[test]
    fn test_diameter() {
        let mut g = Graph::<i32, f64>::new();
        let n1 = g.add_node(1);
        let n2 = g.add_node(2);
        let n3 = g.add_node(3);

        g.add_edge(n1, n2, 1.0);
        g.add_edge(n2, n3, 1.0);

        assert_eq!(diameter(&g), Some(2));
    }

    #[test]
    fn test_diameter_disconnected() {
        let g = Graph::<i32, f64>::new();
        assert_eq!(diameter(&g), None);
    }

    #[test]
    fn test_radius() {
        let mut g = Graph::<i32, f64>::new();
        let n1 = g.add_node(1);
        let n2 = g.add_node(2);
        let n3 = g.add_node(3);

        g.add_edge(n1, n2, 1.0);
        g.add_edge(n2, n3, 1.0);

        assert_eq!(radius(&g), Some(1)); // Center node n2 has eccentricity 1
    }

    #[test]
    fn test_average_clustering_coefficient() {
        let mut g = Graph::<i32, f64>::new();
        let n1 = g.add_node(1);
        let n2 = g.add_node(2);
        let n3 = g.add_node(3);

        g.add_edge(n1, n2, 1.0);
        g.add_edge(n2, n3, 1.0);
        g.add_edge(n3, n1, 1.0);

        assert!((average_clustering_coefficient(&g) - 1.0).abs() < 0.001);
    }

    #[test]
    fn test_transitivity() {
        let mut g = Graph::<i32, f64>::new();
        let n1 = g.add_node(1);
        let n2 = g.add_node(2);
        let n3 = g.add_node(3);

        g.add_edge(n1, n2, 1.0);
        g.add_edge(n2, n3, 1.0);
        g.add_edge(n3, n1, 1.0);

        assert!((transitivity(&g) - 1.0).abs() < 0.001);
    }

    #[test]
    fn test_transitivity_fractional() {
        // A triangle (0, 1, 2) with a pendant edge 2 - 3. There is one triangle,
        // counted once per apex (nodes 0, 1, and 2), giving 3 closed pairs. The
        // connected triples are: node 0 -> 1, node 1 -> 1, node 2 -> C(3, 2) = 3,
        // node 3 -> 0, for 5 total. Transitivity is therefore 3 / 5 = 0.6. This
        // guards the neighbor-set rewrite against miscounting on a hub node.
        let mut g = Graph::<i32, f64>::new();
        let n0 = g.add_node(0);
        let n1 = g.add_node(1);
        let n2 = g.add_node(2);
        let n3 = g.add_node(3);

        g.add_edge(n0, n1, 1.0);
        g.add_edge(n1, n2, 1.0);
        g.add_edge(n2, n0, 1.0);
        g.add_edge(n2, n3, 1.0);

        assert!((transitivity(&g) - 0.6).abs() < 1e-9);
    }

    #[test]
    fn test_transitivity_two_triangles_sharing_an_edge() {
        // Two triangles that share the edge (1, 2): {0, 1, 2} and {1, 2, 3}.
        // Degrees are 0 -> 2, 1 -> 3, 2 -> 3, and 3 -> 2, so the connected triples
        // are C(2, 2) + C(3, 2) + C(3, 2) + C(2, 2) = 1 + 3 + 3 + 1 = 8. There are
        // two distinct triangles, so the transitivity is 3 * 2 / 8 = 0.75. This pins
        // the count when a vertex belongs to more than one triangle, guarding the
        // degree-ordered forward counting against double counting or omission.
        let mut g = Graph::<i32, f64>::new();
        let n0 = g.add_node(0);
        let n1 = g.add_node(1);
        let n2 = g.add_node(2);
        let n3 = g.add_node(3);

        g.add_edge(n0, n1, 1.0);
        g.add_edge(n1, n2, 1.0);
        g.add_edge(n2, n0, 1.0);
        g.add_edge(n1, n3, 1.0);
        g.add_edge(n2, n3, 1.0);

        assert!((transitivity(&g) - 0.75).abs() < 1e-9);
    }

    #[test]
    fn test_average_path_length() {
        let mut g = Graph::<i32, f64>::new();
        let n1 = g.add_node(1);
        let n2 = g.add_node(2);
        let n3 = g.add_node(3);

        g.add_edge(n1, n2, 1.0);
        g.add_edge(n2, n3, 1.0);

        // Paths: 1->2 (1), 2->3 (1), 1->3 (2), avg = 4/3 ≈ 1.33
        let avg = average_path_length(&g).expect("Connected graph should have average path length");
        assert!((avg - 1.333).abs() < 0.01);
    }

    #[test]
    fn test_assortativity() {
        let mut g = Graph::<i32, f64>::new();
        let n1 = g.add_node(1);
        let n2 = g.add_node(2);
        let n3 = g.add_node(3);
        let n4 = g.add_node(4);

        // Create a simple graph
        g.add_edge(n1, n2, 1.0);
        g.add_edge(n2, n3, 1.0);
        g.add_edge(n3, n4, 1.0);

        // Just check it returns a value in valid range
        let assort = assortativity(&g);
        assert!((-1.0..=1.0).contains(&assort));
    }
}