//! Comprehensive numerical validation tests for scirs2-graph
//!
//! This test suite validates the numerical accuracy of graph algorithms
//! against known reference values and mathematical properties.

use approx::assert_abs_diff_eq;
use rand::{rngs::StdRng, SeedableRng};
use scirs2_graph::*;
use std::collections::{HashMap, HashSet};

const TOLERANCE: f64 = 1e-6;
const LOOSE_TOLERANCE: f64 = 1e-4;

/// Load Zachary's Karate Club graph for testing
#[allow(dead_code)]
fn karate_club_graph() -> Graph<usize, f64> {
    let mut g = Graph::new();

    // Add nodes (34 members)
    for i in 0..34 {
        g.add_node(i);
    }

    // Add edges (original Zachary data)
    let edges = vec![
        (0, 1),
        (0, 2),
        (0, 3),
        (0, 4),
        (0, 5),
        (0, 6),
        (0, 7),
        (0, 8),
        (0, 10),
        (0, 11),
        (0, 12),
        (0, 13),
        (0, 17),
        (0, 19),
        (0, 21),
        (0, 31),
        (1, 2),
        (1, 3),
        (1, 7),
        (1, 13),
        (1, 17),
        (1, 19),
        (1, 21),
        (1, 30),
        (2, 3),
        (2, 7),
        (2, 8),
        (2, 9),
        (2, 13),
        (2, 27),
        (2, 28),
        (2, 32),
        (3, 7),
        (3, 12),
        (3, 13),
        (4, 6),
        (4, 10),
        (5, 6),
        (5, 10),
        (5, 16),
        (6, 16),
        (8, 30),
        (8, 32),
        (8, 33),
        (9, 33),
        (13, 33),
        (14, 32),
        (14, 33),
        (15, 32),
        (15, 33),
        (18, 32),
        (18, 33),
        (19, 33),
        (20, 32),
        (20, 33),
        (22, 32),
        (22, 33),
        (23, 25),
        (23, 27),
        (23, 29),
        (23, 32),
        (23, 33),
        (24, 25),
        (24, 27),
        (24, 31),
        (25, 31),
        (26, 29),
        (26, 33),
        (27, 33),
        (28, 31),
        (28, 33),
        (29, 32),
        (29, 33),
        (30, 32),
        (30, 33),
        (31, 32),
        (31, 33),
        (32, 33),
    ];

    for (u, v) in edges {
        g.add_edge(u, v, 1.0).unwrap();
    }

    g
}

/// Create a simple path graph for testing
#[allow(dead_code)]
fn simple_path_graph(n: usize) -> Graph<usize, f64> {
    let mut g = Graph::new();
    for i in 0..n {
        g.add_node(i);
    }
    for i in 0..n - 1 {
        g.add_edge(i, i + 1, 1.0).unwrap();
    }
    g
}

/// Create a complete graph for testing
#[allow(dead_code)]
fn simple_complete_graph(n: usize) -> Graph<usize, f64> {
    let mut g = Graph::new();
    for i in 0..n {
        g.add_node(i);
    }
    for i in 0..n {
        for j in i + 1..n {
            g.add_edge(i, j, 1.0).unwrap();
        }
    }
    g
}

#[test]
#[allow(dead_code)]
fn test_pagerank_numerical_properties() {
    let g = karate_club_graph();
    // Convert to DiGraph for pagerank
    let mut digraph = DiGraph::new();
    for i in 0..g.node_count() {
        digraph.add_node(i);
    }
    for edge in g.edges() {
        digraph
            .add_edge(edge.source, edge.target, edge.weight)
            .unwrap();
    }
    let pr = pagerank(&digraph, 0.85, 1e-8, 100);

    // Property 1: Sum should be 1.0
    let sum: f64 = pr.values().sum();
    assert_abs_diff_eq!(sum, 1.0, epsilon = TOLERANCE);

    // Property 2: All values should be positive
    for &value in pr.values() {
        assert!(value > 0.0);
    }

    // Property 3: Node 0 (instructor) should have high PageRank
    // Node 33 (Mr. Hi) should also have high PageRank
    assert!(pr[&0] > 0.02); // Instructor
    assert!(pr[&33] > 0.02); // Mr. Hi

    // Property 4: PageRank should correlate with degree
    let degrees: HashMap<usize, f64> = g
        .nodes()
        .map(|&n| (n, g.degree(&n).unwrap() as f64))
        .collect();

    let correlation = compute_correlation(&degrees, &pr);
    assert!(correlation > 0.5); // Should be positively correlated
}

#[test]
#[allow(dead_code)]
fn test_betweenness_centrality_properties() {
    let path = simple_path_graph(5);
    let bc = betweenness_centrality(&path, true).unwrap();

    // In a path graph, middle nodes should have highest betweenness
    assert!(bc[&2] > bc[&1]); // Node 2 > Node 1
    assert!(bc[&2] > bc[&3]); // Node 2 > Node 3
    assert!(bc[&1] > bc[&0]); // Node 1 > Node 0 (endpoint)
    assert!(bc[&3] > bc[&4]); // Node 3 > Node 4 (endpoint)

    // Endpoints should have zero betweenness
    assert_abs_diff_eq!(bc[&0], 0.0, epsilon = TOLERANCE);
    assert_abs_diff_eq!(bc[&4], 0.0, epsilon = TOLERANCE);
}

#[test]
#[allow(dead_code)]
fn test_closeness_centrality_properties() {
    let star = star_graph(5); // Center + 4 leaves
    let cc = closeness_centrality(&star).unwrap();

    // Center node (0) should have highest closeness centrality
    let center_closeness = cc[&0];
    for i in 1..5 {
        assert!(center_closeness > cc[&i]);
    }

    // All leaf nodes should have equal closeness
    for i in 2..5 {
        assert_abs_diff_eq!(cc[&1], cc[&i], epsilon = TOLERANCE);
    }
}

#[test]
#[allow(dead_code)]
fn test_eigenvector_centrality_properties() {
    let g = karate_club_graph();
    let ec = eigenvector_centrality(&g).unwrap();

    // Property 1: Should be normalized
    let norm: f64 = ec.values().map(|x| x * x).sum::<f64>().sqrt();
    assert_abs_diff_eq!(norm, 1.0, epsilon = TOLERANCE);

    // Property 2: All values should be non-negative for connected graph
    for &value in ec.values() {
        assert!(value >= 0.0);
    }

    // Property 3: Should correlate with degree for most graphs
    let degrees: HashMap<usize, f64> = g
        .nodes()
        .map(|&n| (n, g.degree(&n).unwrap() as f64))
        .collect();

    let correlation = compute_correlation(&degrees, &ec);
    assert!(correlation > 0.3);
}

#[test]
#[allow(dead_code)]
fn test_clustering_coefficient_known_values() {
    // Triangle: clustering coefficient should be 1.0 for all nodes
    let mut triangle = Graph::new();
    for i in 0..3 {
        triangle.add_node(i);
    }
    triangle.add_edge(0, 1, 1.0).unwrap();
    triangle.add_edge(1, 2, 1.0).unwrap();
    triangle.add_edge(2, 0, 1.0).unwrap();

    let cc = clustering_coefficient(&triangle).unwrap();
    for i in 0..3 {
        assert_abs_diff_eq!(cc[&i], 1.0, epsilon = TOLERANCE);
    }

    // Star graph: center has clustering 0, leaves undefined (no triangles possible)
    let star = star_graph(4);
    let cc_star = clustering_coefficient(&star).unwrap();
    assert_abs_diff_eq!(cc_star[&0], 0.0, epsilon = TOLERANCE); // Center
}

#[test]
#[allow(dead_code)]
fn test_shortest_path_properties() {
    let path = simple_path_graph(10);

    // Distance should equal number of edges in path
    match dijkstra_path(&path, &0, &9) {
        Ok((path_nodes, distance)) => {
            assert_eq!(path_nodes.len(), 10); // All nodes in path
            assert_abs_diff_eq!(distance, 9.0, epsilon = TOLERANCE); // 9 edges
        }
        Err(_) => panic!("Path should exist"),
    }

    // Distance should be symmetric for undirected graphs
    let dist1 = dijkstra_path(&path, &2, &7).unwrap().1;
    let dist2 = dijkstra_path(&path, &7, &2).unwrap().1;
    assert_abs_diff_eq!(dist1, dist2, epsilon = TOLERANCE);

    // Triangle inequality: d(a,c) <= d(a,b) + d(b,c)
    let d_ac = dijkstra_path(&path, &0, &5).unwrap().1;
    let d_ab = dijkstra_path(&path, &0, &3).unwrap().1;
    let d_bc = dijkstra_path(&path, &3, &5).unwrap().1;
    assert!(d_ac <= d_ab + d_bc + TOLERANCE);
}

#[test]
#[allow(dead_code)]
fn test_connected_components_correctness() {
    // Create disconnected graph: two triangles
    let mut g = Graph::new();
    for i in 0..6 {
        g.add_node(i);
    }

    // First triangle
    g.add_edge(0, 1, 1.0).unwrap();
    g.add_edge(1, 2, 1.0).unwrap();
    g.add_edge(2, 0, 1.0).unwrap();

    // Second triangle
    g.add_edge(3, 4, 1.0).unwrap();
    g.add_edge(4, 5, 1.0).unwrap();
    g.add_edge(5, 3, 1.0).unwrap();

    let components = connected_components(&g).unwrap();
    assert_eq!(components.len(), 2);

    // Each component should have 3 nodes
    for component in &components {
        assert_eq!(component.len(), 3);
    }

    // Verify partition property
    let all_nodes: HashSet<usize> = components.iter().flatten().copied().collect();
    assert_eq!(all_nodes.len(), 6);
    assert_eq!(all_nodes, (0..6).collect());
}

#[test]
#[allow(dead_code)]
fn test_minimum_spanning_tree_properties() {
    let complete = simple_complete_graph(5);
    let mst = minimum_spanning_tree(&complete).unwrap();

    // MST should have n-1 edges
    assert_eq!(mst.edge_count(), 4);

    // MST should be connected
    let mst_components = connected_components(&mst).unwrap();
    assert_eq!(mst_components.len(), 1);
    assert_eq!(mst_components[0].len(), 5);

    // MST should be acyclic (no simple test, but connected + n-1 edges guarantees it)
    assert_eq!(mst.node_count(), complete.node_count());
}

#[test]
#[allow(dead_code)]
fn test_community_detection_modularity() {
    let g = karate_club_graph();

    // Test multiple community detection algorithms
    let louvain_result = louvain_communities_result(&g).unwrap();
    let label_prop_result = label_propagation_result(&g).unwrap();

    // Both should find reasonable communities (2-6 is typical for karate club)
    assert!(louvain_result.num_communities >= 2);
    assert!(louvain_result.num_communities <= 6);
    assert!(label_prop_result.num_communities >= 2);
    assert!(label_prop_result.num_communities <= 6);

    // Modularity should be positive and reasonable
    assert!(louvain_result.modularity > 0.3);
    assert!(label_prop_result.modularity > 0.2);

    // Each node should be assigned to exactly one community
    let mut node_assignments = 0;
    for community in &louvain_result.communities {
        node_assignments += community.len();
    }
    assert_eq!(node_assignments, g.node_count());
}

#[test]
#[allow(dead_code)]
fn test_k_core_decomposition_properties() {
    let g = karate_club_graph();
    let k_cores = k_core_decomposition(&g).unwrap();

    // Every node should be assigned to some k-core
    let total_nodes: usize = k_cores.values().map(|nodes| nodes.len()).sum();
    assert_eq!(total_nodes, g.node_count());

    // k-core property: every node in k-core should have degree >= k within the k-core
    for (&k, nodes) in &k_cores {
        if k > 0 {
            let subgraph = subgraph(&g, nodes).unwrap();
            for &node in nodes {
                if let Some(remapped_node) = subgraph.nodes().find(|&&n| n == node) {
                    let degree = subgraph.degree(remapped_node).unwrap_or(0);
                    assert!(
                        degree >= k,
                        "Node {} in {}-core has degree {} < {}",
                        node,
                        k,
                        degree,
                        k
                    );
                }
            }
        }
    }
}

#[test]
#[allow(dead_code)]
fn test_graph_density_calculation() {
    // Empty graph
    let empty = Graph::<usize, f64>::new();
    assert_abs_diff_eq!(graph_density(&empty), 0.0, epsilon = TOLERANCE);

    // Single node
    let mut single = Graph::new();
    single.add_node(0);
    assert_abs_diff_eq!(graph_density(&single), 0.0, epsilon = TOLERANCE);

    // Complete graph: density should be 1.0
    let complete = simple_complete_graph(5);
    assert_abs_diff_eq!(graph_density(&complete), 1.0, epsilon = TOLERANCE);

    // Path graph: density = 2(n-1) / (n(n-1)) = 2/n for n > 1
    let path = simple_path_graph(10);
    let expected_density = 2.0 / 10.0;
    assert_abs_diff_eq!(graph_density(&path), expected_density, epsilon = TOLERANCE);
}

#[test]
#[allow(dead_code)]
fn test_random_walk_properties() {
    let g = karate_club_graph();
    let mut rng = StdRng::seed_from_u64(42);

    // Random walk should visit connected nodes
    let walk = random_walk(&g, &0, 100, &mut rng).unwrap();
    assert_eq!(walk.len(), 101); // Starting node + 100 steps

    // All nodes in walk should exist in graph
    for &node in &walk {
        assert!(g.has_node(&node));
    }

    // Consecutive nodes in walk should be connected
    for window in walk.windows(2) {
        let (current, next) = (window[0], window[1]);
        assert!(g.has_edge(&current, &next) || current == next);
    }

    // Stationary distribution should approximate degree distribution for long walks
    let long_walk = random_walk(&g, &0, 10000, &mut rng).unwrap();
    let mut visit_counts: HashMap<usize, usize> = HashMap::new();
    for &node in &long_walk {
        *visit_counts.entry(node).or_insert(0) += 1;
    }

    // Convert to probabilities
    let total_visits = long_walk.len() as f64;
    let visit_probs: HashMap<usize, f64> = visit_counts
        .iter()
        .map(|(&node, &count)| (node, count as f64 / total_visits))
        .collect();

    // Should correlate with degree (for simple random walk)
    let degrees: HashMap<usize, f64> = g
        .nodes()
        .map(|&n| (n, g.degree(&n).unwrap() as f64))
        .collect();

    // Normalize degrees to probabilities
    let total_degree: f64 = degrees.values().sum();
    let degree_probs: HashMap<usize, f64> = degrees
        .iter()
        .map(|(&node, &deg)| (node, deg / total_degree))
        .collect();

    let correlation = compute_correlation(&visit_probs, &degree_probs);
    assert!(correlation > 0.5);
}

#[test]
#[allow(dead_code)]
fn test_bipartite_detection() {
    // Create a bipartite graph
    let mut bipartite = Graph::new();
    for i in 0..6 {
        bipartite.add_node(i);
    }

    // Connect partition A (0,1,2) to partition B (3,4,5)
    for i in 0..3 {
        for j in 3..6 {
            bipartite.add_edge(i, j, 1.0).unwrap();
        }
    }

    let (is_bip, partition) = is_bipartite(&bipartite).unwrap();
    assert!(is_bip);

    if let Some((set_a, set_b)) = partition {
        // Each partition should have 3 nodes
        assert_eq!(set_a.len(), 3);
        assert_eq!(set_b.len(), 3);

        // Partitions should be disjoint and cover all nodes
        let union: HashSet<_> = set_a.union(&set_b).copied().collect();
        assert_eq!(union.len(), 6);
    }

    // Non-bipartite graph (triangle)
    let triangle = simple_complete_graph(3);
    let (is_bip_tri, _) = is_bipartite(&triangle).unwrap();
    assert!(!is_bip_tri);
}

#[test]
#[allow(dead_code)]
fn test_graph_isomorphism_simple_cases() {
    // Two identical paths should be isomorphic
    let path1 = simple_path_graph(5);
    let path2 = simple_path_graph(5);
    assert!(are_graphs_isomorphic(&path1, &path2).unwrap());

    // Path and cycle of same size should not be isomorphic
    let path = simple_path_graph(5);
    let cycle = cycle_graph(5);
    assert!(!are_graphs_isomorphic(&path, &cycle).unwrap());

    // Graphs with different number of nodes cannot be isomorphic
    let small_path = simple_path_graph(3);
    let large_path = simple_path_graph(5);
    assert!(!are_graphs_isomorphic(&small_path, &large_path).unwrap());
}

#[test]
#[allow(dead_code)]
fn test_numerical_stability_large_values() {
    // Test with large edge weights
    let mut g = Graph::new();
    for i in 0..5 {
        g.add_node(i);
    }
    g.add_edge(0, 1, 1e6).unwrap();
    g.add_edge(1, 2, 1e6).unwrap();
    g.add_edge(2, 3, 1e6).unwrap();
    g.add_edge(3, 4, 1e6).unwrap();

    // PageRank should still converge and sum to 1
    let pr = pagerank(&g, 0.85, Some(1e-8)).unwrap();
    let sum: f64 = pr.values().sum();
    assert_abs_diff_eq!(sum, 1.0, epsilon = TOLERANCE);

    // Shortest path should handle large weights correctly
    let (path, distance) = dijkstra_path(&g, &0, &4).unwrap();
    assert_eq!(path.len(), 5);
    assert_abs_diff_eq!(distance, 4e6, epsilon = 1.0); // Allow small numerical error
}

/// Helper function to compute Pearson correlation coefficient
#[allow(dead_code)]
fn compute_correlation(x: &HashMap<usize, f64>, y: &HashMap<usize, f64>) -> f64 {
    let common_keys: Vec<_> = x.keys().filter(|k| y.contains_key(k)).collect();
    if common_keys.is_empty() {
        return 0.0;
    }

    let n = common_keys.len() as f64;

    let sum_x: f64 = common_keys.iter().map(|&k| x[k]).sum();
    let sum_y: f64 = common_keys.iter().map(|&k| y[k]).sum();
    let sum_x2: f64 = common_keys.iter().map(|&k| x[k] * x[k]).sum();
    let sum_y2: f64 = common_keys.iter().map(|&k| y[k] * y[k]).sum();
    let sum_xy: f64 = common_keys.iter().map(|&k| x[k] * y[k]).sum();

    let numerator = n * sum_xy - sum_x * sum_y;
    let denominator = ((n * sum_x2 - sum_x * sum_x) * (n * sum_y2 - sum_y * sum_y)).sqrt();

    if denominator.abs() < 1e-10 {
        0.0
    } else {
        numerator / denominator
    }
}

#[test]
#[allow(dead_code)]
fn test_algorithm_determinism() {
    // Algorithms with fixed seeds should produce identical results
    let mut rng1 = StdRng::seed_from_u64(12345);
    let mut rng2 = StdRng::seed_from_u64(12345);

    let g1 = erdos_renyi_graph(100, 0.1, &mut rng1).unwrap();
    let g2 = erdos_renyi_graph(100, 0.1, &mut rng2).unwrap();

    // Graphs should be identical
    assert_eq!(g1.node_count(), g2.node_count());
    assert_eq!(g1.edge_count(), g2.edge_count());

    // PageRank should be identical
    let pr1 = pagerank(&g1, 0.85, Some(1e-8)).unwrap();
    let pr2 = pagerank(&g2, 0.85, Some(1e-8)).unwrap();

    for node in g1.nodes() {
        assert_abs_diff_eq!(pr1[node], pr2[node], epsilon = TOLERANCE);
    }
}

#[test]
#[allow(dead_code)]
fn test_edge_cases() {
    // Single node graph
    let mut single = Graph::new();
    single.add_node(42);

    // PageRank of single node should be 1.0
    let pr = pagerank(&single, 0.85, Some(1e-6)).unwrap();
    assert_abs_diff_eq!(pr[&42], 1.0, epsilon = TOLERANCE);

    // Connected components should find the single node
    let components = connected_components(&single).unwrap();
    assert_eq!(components.len(), 1);
    assert_eq!(components[0], vec![42]);

    // Clustering coefficient should be 0 (no triangles possible)
    let cc = clustering_coefficient(&single).unwrap();
    assert_abs_diff_eq!(cc[&42], 0.0, epsilon = TOLERANCE);

    // Two isolated nodes
    let mut isolated = Graph::new();
    isolated.add_node(1);
    isolated.add_node(2);

    let iso_components = connected_components(&isolated).unwrap();
    assert_eq!(iso_components.len(), 2);

    // No path should exist between isolated nodes
    assert!(dijkstra_path(&isolated, &1, &2).is_err());
}
