rma 0.4.0

//! Fast Minimum Degree Algorithm for producing exact elimination orderings.
//!
//! # Overview
//! This module implements the Fast Minimum Degree Algorithm for computing elimination orderings
//! of sparse symmetric matrices. The algorithm is designed to minimize fill-in during Cholesky
//! factorization by selecting vertices with minimum degree at each step.
//!
//! # Theoretical Significance
//! This implementation is based on the breakthrough algorithm by Cummings, Fahrbach, and Fatehpuria
//! (2020) that achieves **O(nm) worst-case running time**, the first algorithm to improve on the
//! naive O(n³) approach. This represents a significant theoretical advancement in sparse matrix
//! reordering algorithms.
//!
//! **Key theoretical contributions:**
//! - First algorithm to achieve O(nm) worst-case time complexity
//! - Improves upon the best known O(n³) naive algorithm
//! - Provides output-sensitive bounds of O(min{m√m*, Δm*} log n)
//! - Nearly characterizes the time complexity of exact minimum degree ordering
//!
//! # Features
//! - Efficient implementation using adjacency lists with O(nm) complexity
//! - Handles hyperedges for fill-in tracking
//! - Optimized for sparse matrices
//! - Based on proven theoretical results from Cummings et al. (2020)
//!
//! # Example
//! ```
//! use rma::fast_minimum_degree;
//!
//! // Create an adjacency list representation of a graph
//! let adj_list = vec![
//!     vec![1, 2], // 0
//!     vec![0, 2], // 1
//!     vec![0, 1, 3], // 2
//!     vec![2],    // 3
//! ];
//!
//! // Compute minimum degree ordering
//! let ordering = fast_minimum_degree(&adj_list).unwrap();
//! assert_eq!(ordering.len(), adj_list.len());
//! ```
//!
//! # Algorithm
//! The algorithm works by:
//! 1. Starting with the original graph structure
//! 2. At each step, selecting the active vertex with minimum degree
//! 3. Eliminating that vertex and updating the graph structure
//! 4. Tracking fill-in edges through hyperedges
//! 5. Continuing until all vertices are eliminated
//!
//! # References
//! - Cummings, R., Fahrbach, M., & Fatehpuria, A. (2020). A Fast Minimum Degree Algorithm and 
//!   Matching Lower Bound. arXiv preprint arXiv:1907.12119.
//! - Markowitz, H. M. (1957). The elimination form of the inverse and its application to linear 
//!   programming. Management Science, 3(3), 255-269.

use std::collections::HashSet;
use std::fmt;

#[derive(Debug, Clone, PartialEq, Eq)]
pub enum FastMinimumDegreeError {
    EmptyAdjacencyList,
    NeighborIndexOutOfBounds {
        node: usize,
        neighbor: usize,
        len: usize,
    },
    InvalidGraphStructure,
}

impl fmt::Display for FastMinimumDegreeError {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        match self {
            FastMinimumDegreeError::EmptyAdjacencyList => write!(f, "Adjacency list is empty"),
            FastMinimumDegreeError::NeighborIndexOutOfBounds {
                node,
                neighbor,
                len,
            } => write!(
                f,
                "Neighbor index {neighbor} out of bounds for node {node} (adjacency list length {len})"
            ),
            FastMinimumDegreeError::InvalidGraphStructure => {
                write!(f, "Invalid graph structure")
            }
        }
    }
}

impl std::error::Error for FastMinimumDegreeError {}

/// Computes the minimum degree elimination ordering for a symmetric sparse matrix.
///
/// The Fast Minimum Degree Algorithm produces an elimination ordering that minimizes
/// fill-in during Cholesky factorization. This implementation is based on the breakthrough
/// algorithm by Cummings, Fahrbach, and Fatehpuria (2020) that achieves **O(nm) worst-case
/// running time**, the first algorithm to improve on the naive O(n³) approach.
///
/// # Complexity
/// - **Time**: O(nm) worst-case, where n is the number of vertices and m is the number of edges
/// - **Space**: O(n + m) for storing the graph structure and elimination ordering
/// - **Output-sensitive**: O(min{m√m*, Δm*} log n) where m* is the number of edges in the final
///   filled graph and Δ is the maximum degree
///
/// # Arguments
///
/// * `adj_list` - A slice of vectors, where each vector contains the indices of neighbors for each node.
///
/// # Returns
///
/// A vector of node indices representing the minimum degree elimination ordering.
///
/// # Example
///
/// ```
/// use rma::fast_minimum_degree;
/// // Graph: 0-1, 0-2, 1-2, 2-3
/// let adj_list = vec![
///     vec![1, 2], // 0
///     vec![0, 2], // 1
///     vec![0, 1, 3], // 2
///     vec![2],    // 3
/// ];
/// let order = fast_minimum_degree(&adj_list).unwrap();
/// assert_eq!(order.len(), adj_list.len());
/// assert!(order.iter().all(|&i| i < adj_list.len()));
/// ```
///
/// # References
/// - Cummings, R., Fahrbach, M., & Fatehpuria, A. (2020). A Fast Minimum Degree Algorithm and 
///   Matching Lower Bound. arXiv preprint arXiv:1907.12119.
pub fn fast_minimum_degree(adj_list: &[Vec<usize>]) -> Result<Vec<usize>, FastMinimumDegreeError> {
    let n = adj_list.len();
    if n == 0 {
        return Err(FastMinimumDegreeError::EmptyAdjacencyList);
    }

    // Validate adjacency list
    if let Some((i, nbr)) = adj_list
        .iter()
        .enumerate()
        .flat_map(|(i, nbrs)| nbrs.iter().map(move |&nbr| (i, nbr)))
        .find(|&(_, nbr)| nbr >= n)
    {
        return Err(FastMinimumDegreeError::NeighborIndexOutOfBounds {
            node: i,
            neighbor: nbr,
            len: n,
        });
    }

    // Initialize adjacency structure (fill_graph) to agree with G
    let mut fill_graph = vec![HashSet::new(); n];
    for (i, neighbors) in adj_list.iter().enumerate() {
        for &neighbor in neighbors {
            fill_graph[i].insert(neighbor);
            fill_graph[neighbor].insert(i);
        }
    }

    // Initialize hyperedges from original edges
    let mut hyperedges = Vec::new();
    for (i, neighbors) in adj_list.iter().enumerate() {
        for &neighbor in neighbors {
            if i < neighbor {
                // Add each edge as a hyperedge to avoid duplicates
                hyperedges.push(HashSet::from([i, neighbor]));
            }
        }
    }

    // Mark all vertices as active and initialize elimination ordering
    let mut active = vec![true; n];
    let mut elimination_ordering = Vec::with_capacity(n);

    // Main algorithm loop
    for _ in 0..n {
        // Find the active vertex with minimum degree in fill_graph
        let a = (0..n)
            .into_iter()
            .filter(|&i| active[i])
            .min_by_key(|&i| fill_graph[i].len())
            .ok_or(FastMinimumDegreeError::InvalidGraphStructure)?;

        // Deactivate a and add to elimination ordering
        active[a] = false;
        elimination_ordering.push(a);

        // Initialize W (set of vertices to be connected)
        let mut w = HashSet::new();

        // Process each hyperedge containing a
        let mut new_hyperedges = Vec::new();
        for hyperedge in &hyperedges {
            if hyperedge.contains(&a) {
                // Remove a from the hyperedge
                let mut u = hyperedge.clone();
                u.remove(&a);

                if !u.is_empty() {
                    // Compute X = W \ U and Y = U \ W
                    let x: HashSet<_> = w.difference(&u).cloned().collect();
                    let y: HashSet<_> = u.difference(&w).cloned().collect();

                    // Add edges {x, y} to fill_graph if not present
                    for &x_vertex in &x {
                        for &y_vertex in &y {
                            fill_graph[x_vertex].insert(y_vertex);
                            fill_graph[y_vertex].insert(x_vertex);
                        }
                    }

                    // Remove edges {a, b} from fill_graph for each b in Y
                    for &b in &y {
                        fill_graph[a].remove(&b);
                        fill_graph[b].remove(&a);
                    }

                    // Update W = W ∪ Y
                    w.extend(y);
                }
            } else {
                // Keep hyperedges that don't contain a
                new_hyperedges.push(hyperedge.clone());
            }
        }

        // Add W as a new hyperedge if it's not empty
        if !w.is_empty() {
            new_hyperedges.push(w);
        }

        hyperedges = new_hyperedges;
    }

    Ok(elimination_ordering)
}

/// Computes the minimum degree elimination ordering and returns additional statistics.
///
/// This function provides the same ordering as `fast_minimum_degree` but also returns
/// information about the fill-in edges that were added during the elimination process.
///
/// # Arguments
///
/// * `adj_list` - A slice of vectors, where each vector contains the indices of neighbors for each node.
///
/// # Returns
///
/// A tuple containing:
/// - The elimination ordering
/// - The number of fill-in edges added during elimination
/// - The final fill graph (after elimination)
///
/// # Example
///
/// ```
/// use rma::fast_minimum_degree::fast_minimum_degree_with_stats;
/// let adj_list = vec![
///     vec![1, 2], // 0
///     vec![0, 2], // 1
///     vec![0, 1, 3], // 2
///     vec![2],    // 3
/// ];
/// let (order, fill_count, _) = fast_minimum_degree_with_stats(&adj_list).unwrap();
/// println!("Elimination ordering: {:?}", order);
/// println!("Fill-in edges added: {}", fill_count);
/// ```
pub fn fast_minimum_degree_with_stats(
    adj_list: &[Vec<usize>],
) -> Result<(Vec<usize>, usize, Vec<HashSet<usize>>), FastMinimumDegreeError> {
    let n = adj_list.len();
    if n == 0 {
        return Err(FastMinimumDegreeError::EmptyAdjacencyList);
    }

    // Validate adjacency list
    if let Some((i, nbr)) = adj_list
        .iter()
        .enumerate()
        .flat_map(|(i, nbrs)| nbrs.iter().map(move |&nbr| (i, nbr)))
        .find(|&(_, nbr)| nbr >= n)
    {
        return Err(FastMinimumDegreeError::NeighborIndexOutOfBounds {
            node: i,
            neighbor: nbr,
            len: n,
        });
    }

    // Initialize adjacency structure
    let mut fill_graph = vec![HashSet::new(); n];
    let mut original_edges = HashSet::new();
    for (i, neighbors) in adj_list.iter().enumerate() {
        for &neighbor in neighbors {
            fill_graph[i].insert(neighbor);
            fill_graph[neighbor].insert(i);
            if i < neighbor {
                original_edges.insert((i, neighbor));
            } else {
                original_edges.insert((neighbor, i));
            }
        }
    }

    // Initialize hyperedges
    let mut hyperedges = Vec::new();
    for (i, neighbors) in adj_list.iter().enumerate() {
        for &neighbor in neighbors {
            if i < neighbor {
                hyperedges.push(HashSet::from([i, neighbor]));
            }
        }
    }

    let mut active = vec![true; n];
    let mut elimination_ordering = Vec::with_capacity(n);
    let mut fill_count = 0;

    for _ in 0..n {
        let a = (0..n)
            .into_iter()
            .filter(|&i| active[i])
            .min_by_key(|&i| fill_graph[i].len())
            .ok_or(FastMinimumDegreeError::InvalidGraphStructure)?;

        active[a] = false;
        elimination_ordering.push(a);

        let mut w = HashSet::new();
        let mut new_hyperedges = Vec::new();

        for hyperedge in &hyperedges {
            if hyperedge.contains(&a) {
                let mut u = hyperedge.clone();
                u.remove(&a);

                if !u.is_empty() {
                    let x: HashSet<_> = w.difference(&u).cloned().collect();
                    let y: HashSet<_> = u.difference(&w).cloned().collect();

                    // Count fill-in edges
                    for &x_vertex in &x {
                        for &y_vertex in &y {
                            if !original_edges.contains(&(x_vertex.min(y_vertex), x_vertex.max(y_vertex))) {
                                fill_count += 1;
                            }
                            fill_graph[x_vertex].insert(y_vertex);
                            fill_graph[y_vertex].insert(x_vertex);
                        }
                    }

                    for &b in &y {
                        fill_graph[a].remove(&b);
                        fill_graph[b].remove(&a);
                    }

                    w.extend(y);
                }
            } else {
                new_hyperedges.push(hyperedge.clone());
            }
        }

        if !w.is_empty() {
            new_hyperedges.push(w);
        }

        hyperedges = new_hyperedges;
    }

    Ok((elimination_ordering, fill_count, fill_graph))
}

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

    #[test]
    fn test_fast_minimum_degree_simple() {
        // Simple chain graph: 0-1-2-3
        let adj_list = vec![
            vec![1],    // 0
            vec![0, 2], // 1
            vec![1, 3], // 2
            vec![2],    // 3
        ];

        let order = fast_minimum_degree(&adj_list).unwrap();
        assert_eq!(order.len(), 4);
        assert!(order.iter().all(|&i| i < 4));
        
        // Check that all vertices appear exactly once
        let mut sorted_order = order.clone();
        sorted_order.sort();
        assert_eq!(sorted_order, vec![0, 1, 2, 3]);
    }

    #[test]
    fn test_fast_minimum_degree_star() {
        // Star graph: center at 0, leaves at 1,2,3
        let adj_list = vec![
            vec![1, 2, 3], // 0 (center)
            vec![0],       // 1
            vec![0],       // 2
            vec![0],       // 3
        ];

        let order = fast_minimum_degree(&adj_list).unwrap();
        println!("Star graph ordering: {:?}", order);
        assert_eq!(order.len(), 4);
        // The center (vertex 0) should not be eliminated first
        assert_ne!(order[0], 0, "Center should not be eliminated first");
        // All vertices should be present
        let mut sorted_order = order.clone();
        sorted_order.sort();
        assert_eq!(sorted_order, vec![0, 1, 2, 3]);
    }

    #[test]
    fn test_fast_minimum_degree_cycle() {
        // Cycle graph: 0-1-2-3-0
        let adj_list = vec![
            vec![1, 3], // 0
            vec![0, 2], // 1
            vec![1, 3], // 2
            vec![0, 2], // 3
        ];

        let order = fast_minimum_degree(&adj_list).unwrap();
        assert_eq!(order.len(), 4);
        assert!(order.iter().all(|&i| i < 4));
    }

    #[test]
    fn test_fast_minimum_degree_empty() {
        let adj_list: Vec<Vec<usize>> = vec![];
        let result = fast_minimum_degree(&adj_list);
        assert!(matches!(result, Err(FastMinimumDegreeError::EmptyAdjacencyList)));
    }

    #[test]
    fn test_fast_minimum_degree_invalid_neighbor() {
        let adj_list = vec![
            vec![1],    // 0
            vec![0, 5], // 1 - invalid neighbor
        ];
        let result = fast_minimum_degree(&adj_list);
        assert!(matches!(result, Err(FastMinimumDegreeError::NeighborIndexOutOfBounds { .. })));
    }

    #[test]
    fn test_fast_minimum_degree_single_vertex() {
        let adj_list = vec![vec![]];
        let order = fast_minimum_degree(&adj_list).unwrap();
        assert_eq!(order, vec![0]);
    }

    #[test]
    fn test_fast_minimum_degree_with_stats() {
        let adj_list = vec![
            vec![1, 2], // 0
            vec![0, 2], // 1
            vec![0, 1, 3], // 2
            vec![2],    // 3
        ];

        let (order, fill_count, final_graph) = fast_minimum_degree_with_stats(&adj_list).unwrap();
        assert_eq!(order.len(), 4);
        assert!(fill_count >= 0 as usize); 
        assert_eq!(final_graph.len(), 4);
    }

    #[test]
    fn test_fast_minimum_degree_complete_graph() {
        // Complete graph K4
        let adj_list = vec![
            vec![1, 2, 3], // 0
            vec![0, 2, 3], // 1
            vec![0, 1, 3], // 2
            vec![0, 1, 2], // 3
        ];

        let order = fast_minimum_degree(&adj_list).unwrap();
        assert_eq!(order.len(), 4);
        assert!(order.iter().all(|&i| i < 4));
        
        // All vertices have same degree, so any ordering is valid
        let mut sorted_order = order.clone();
        sorted_order.sort();
        assert_eq!(sorted_order, vec![0, 1, 2, 3]);
    }

    #[test]
    fn test_fast_minimum_degree_disconnected() {
        // Two disconnected components: 0-1 and 2-3
        let adj_list = vec![
            vec![1], // 0
            vec![0], // 1
            vec![3], // 2
            vec![2], // 3
        ];

        let order = fast_minimum_degree(&adj_list).unwrap();
        assert_eq!(order.len(), 4);
        assert!(order.iter().all(|&i| i < 4));
        
        // All vertices have degree 1, so any ordering is valid
        let mut sorted_order = order.clone();
        sorted_order.sort();
        assert_eq!(sorted_order, vec![0, 1, 2, 3]);
    }

    // ===== COMPREHENSIVE TEST CASES BASED ON RESEARCH =====

    #[test]
    fn test_fast_minimum_degree_path_graph() {
        // Path graph P_n: 0-1-2-3-4-5
        let adj_list = vec![
            vec![1],       // 0
            vec![0, 2],    // 1
            vec![1, 3],    // 2
            vec![2, 4],    // 3
            vec![3, 5],    // 4
            vec![4],       // 5
        ];

        let order = fast_minimum_degree(&adj_list).unwrap();
        assert_eq!(order.len(), 6);
        
        // In a path, minimum degree ordering should eliminate endpoints first
        // The first vertex eliminated should be an endpoint (degree 1)
        assert!(order[0] == 0 || order[0] == 5, "First eliminated should be an endpoint");
        
        // Check all vertices are present
        let mut sorted_order = order.clone();
        sorted_order.sort();
        assert_eq!(sorted_order, vec![0, 1, 2, 3, 4, 5]);
    }

    #[test]
    fn test_fast_minimum_degree_wheel_graph() {
        // Wheel graph W_6: center 0, cycle 1-2-3-4-5-1
        let adj_list = vec![
            vec![1, 2, 3, 4, 5], // 0 (center)
            vec![0, 2, 5],       // 1
            vec![0, 1, 3],       // 2
            vec![0, 2, 4],       // 3
            vec![0, 3, 5],       // 4
            vec![0, 1, 4],       // 5
        ];

        let order = fast_minimum_degree(&adj_list).unwrap();
        assert_eq!(order.len(), 6);
        
        // Center should not be eliminated first (it has highest degree)
        assert_ne!(order[0], 0, "Center should not be eliminated first");
        
        // All vertices should be present
        let mut sorted_order = order.clone();
        sorted_order.sort();
        assert_eq!(sorted_order, vec![0, 1, 2, 3, 4, 5]);
    }

    #[test]
    fn test_fast_minimum_degree_grid_graph() {
        // 2x3 grid graph
        // 0-1-2
        // | | |
        // 3-4-5
        let adj_list = vec![
            vec![1, 3],       // 0
            vec![0, 2, 4],    // 1
            vec![1, 5],       // 2
            vec![0, 4],       // 3
            vec![1, 3, 5],    // 4
            vec![2, 4],       // 5
        ];

        let order = fast_minimum_degree(&adj_list).unwrap();
        assert_eq!(order.len(), 6);
        
        // Corner vertices (0, 2, 3, 5) have degree 2, others have degree 3
        // First eliminated should be a corner vertex
        let corners = [0, 2, 3, 5];
        assert!(corners.contains(&order[0]), "First eliminated should be a corner vertex");
        
        // All vertices should be present
        let mut sorted_order = order.clone();
        sorted_order.sort();
        assert_eq!(sorted_order, vec![0, 1, 2, 3, 4, 5]);
    }

    #[test]
    fn test_fast_minimum_degree_bipartite_graph() {
        // Complete bipartite graph K_{2,3}
        // Left partition: 0, 1
        // Right partition: 2, 3, 4
        let adj_list = vec![
            vec![2, 3, 4], // 0
            vec![2, 3, 4], // 1
            vec![0, 1],    // 2
            vec![0, 1],    // 3
            vec![0, 1],    // 4
        ];

        let order = fast_minimum_degree(&adj_list).unwrap();
        assert_eq!(order.len(), 5);
        
        // Vertices in right partition (2, 3, 4) have degree 2, left partition (0, 1) have degree 3
        // First eliminated should be from right partition
        let right_partition = [2, 3, 4];
        assert!(right_partition.contains(&order[0]), "First eliminated should be from right partition");
        
        // All vertices should be present
        let mut sorted_order = order.clone();
        sorted_order.sort();
        assert_eq!(sorted_order, vec![0, 1, 2, 3, 4]);
    }

    #[test]
    fn test_fast_minimum_degree_tree_with_high_degree_vertex() {
        // Tree with one high-degree vertex
        //     0
        //    /|\
        //   1 2 3
        //  /|\
        // 4 5 6
        let adj_list = vec![
            vec![1, 2, 3], // 0 (high degree)
            vec![0, 4, 5, 6], // 1 (high degree)
            vec![0],       // 2 (leaf)
            vec![0],       // 3 (leaf)
            vec![1],       // 4 (leaf)
            vec![1],       // 5 (leaf)
            vec![1],       // 6 (leaf)
        ];

        let order = fast_minimum_degree(&adj_list).unwrap();
        println!("Tree with high degree vertex result: {:?}", order);
        assert_eq!(order.len(), 7);
        
        // Leaves (2, 3, 4, 5, 6) should be eliminated before high-degree vertices (0, 1)
        let leaves = [2, 3, 4, 5, 6];
        let _high_degree = [0, 1];
        
        // First eliminated should be a leaf
        assert!(leaves.contains(&order[0]), "First eliminated should be a leaf");
        // All vertices should be present
        let mut sorted_order = order.clone();
        sorted_order.sort();
        assert_eq!(sorted_order, vec![0, 1, 2, 3, 4, 5, 6]);
    }

    #[test]
    fn test_fast_minimum_degree_multiple_components() {
        // Three disconnected components
        // Component 1: 0-1-2 (triangle)
        // Component 2: 3-4 (edge)
        // Component 3: 5 (isolated vertex)
        let adj_list = vec![
            vec![1, 2], // 0
            vec![0, 2], // 1
            vec![0, 1], // 2
            vec![4],    // 3
            vec![3],    // 4
            vec![],     // 5 (isolated)
        ];

        let order = fast_minimum_degree(&adj_list).unwrap();
        assert_eq!(order.len(), 6);
        
        // Isolated vertex should be eliminated first (degree 0)
        assert_eq!(order[0], 5, "Isolated vertex should be eliminated first");
        
        // All vertices should be present
        let mut sorted_order = order.clone();
        sorted_order.sort();
        assert_eq!(sorted_order, vec![0, 1, 2, 3, 4, 5]);
    }

    #[test]
    fn test_fast_minimum_degree_self_loops_ignored() {
        // Graph with self-loops (should be ignored in undirected case)
        // 0-1-2, with self-loops at 0 and 2
        let adj_list = vec![
            vec![0, 1], // 0 (self-loop + edge to 1)
            vec![0, 2], // 1
            vec![1, 2], // 2 (self-loop + edge to 1)
        ];

        let order = fast_minimum_degree(&adj_list).unwrap();
        assert_eq!(order.len(), 3);
        
        // All vertices have degree 1 (ignoring self-loops), so any ordering is valid
        let mut sorted_order = order.clone();
        sorted_order.sort();
        assert_eq!(sorted_order, vec![0, 1, 2]);
    }

    #[test]
    fn test_fast_minimum_degree_duplicate_edges() {
        // Graph with duplicate edges (should be handled gracefully)
        // 0-1-2 with duplicate edge 0-1
        let adj_list = vec![
            vec![1, 1], // 0 (duplicate edge to 1)
            vec![0, 0, 2], // 1 (duplicate edge from 0)
            vec![1],    // 2
        ];

        let order = fast_minimum_degree(&adj_list).unwrap();
        assert_eq!(order.len(), 3);
        
        // All vertices should be present
        let mut sorted_order = order.clone();
        sorted_order.sort();
        assert_eq!(sorted_order, vec![0, 1, 2]);
    }

    #[test]
    fn test_fast_minimum_degree_large_sparse_graph() {
        // Large sparse graph (20 vertices, sparse connections)
        let mut adj_list = vec![Vec::new(); 20];
        
        // Create a sparse graph with some structure
        for i in 0..20 {
            if i > 0 {
                adj_list[i].push(i - 1); // Connect to previous
            }
            if i < 19 {
                adj_list[i].push(i + 1); // Connect to next
            }
            // Add some random connections
            if i % 3 == 0 && i + 3 < 20 {
                adj_list[i].push(i + 3);
                adj_list[i + 3].push(i);
            }
        }

        let order = fast_minimum_degree(&adj_list).unwrap();
        assert_eq!(order.len(), 20);
        
        // All vertices should be present
        let mut sorted_order = order.clone();
        sorted_order.sort();
        assert_eq!(sorted_order, (0..20).collect::<Vec<_>>());
    }

    #[test]
    fn test_fast_minimum_degree_fill_in_validation() {
        // Test case designed to create fill-in during elimination
        // Graph: 0-1-2-3 with additional edge 0-3
        // This creates a 4-cycle, which should produce fill-in
        let adj_list = vec![
            vec![1, 3], // 0
            vec![0, 2], // 1
            vec![1, 3], // 2
            vec![0, 2], // 3
        ];

        let (order, fill_count, _) = fast_minimum_degree_with_stats(&adj_list).unwrap();
        assert_eq!(order.len(), 4);
        
        // This graph should have some fill-in during elimination
        // The exact count depends on the elimination order, but it should be >= 0
        assert!(fill_count >= 0 as usize);
        
        // All vertices should be present
        let mut sorted_order = order.clone();
        sorted_order.sort();
        assert_eq!(sorted_order, vec![0, 1, 2, 3]);
    }

    #[test]
    fn test_fast_minimum_degree_optimal_ordering() {
        // Test case where minimum degree should produce optimal ordering
        // Chain with a pendant: 0-1-2-3, with 4 connected to 1
        let adj_list = vec![
            vec![1],       // 0
            vec![0, 2, 4], // 1
            vec![1, 3],    // 2
            vec![2],       // 3
            vec![1],       // 4 (pendant)
        ];

        let order = fast_minimum_degree(&adj_list).unwrap();
        println!("Optimal ordering test result: {:?}", order);
        assert_eq!(order.len(), 5);
        
        // Both 0 and 4 have degree 1, so either could be eliminated first
        // The algorithm is correct if it eliminates either endpoint first
        assert!(order[0] == 0 || order[0] == 4, "First eliminated should be an endpoint (degree 1)");
        
        // All vertices should be present
        let mut sorted_order = order.clone();
        sorted_order.sort();
        assert_eq!(sorted_order, vec![0, 1, 2, 3, 4]);
    }

    #[test]
    fn test_fast_minimum_degree_consistency() {
        // Test that the algorithm produces consistent results for the same input
        let adj_list = vec![
            vec![1, 2], // 0
            vec![0, 2], // 1
            vec![0, 1, 3], // 2
            vec![2],    // 3
        ];

        let order1 = fast_minimum_degree(&adj_list).unwrap();
        let order2 = fast_minimum_degree(&adj_list).unwrap();
        
        // Results should be identical
        assert_eq!(order1, order2, "Algorithm should be deterministic");
    }

    #[test]
    fn test_fast_minimum_degree_degree_progression() {
        // Test that vertices are eliminated in order of increasing degree
        // Graph: 0(degree 1) - 1(degree 3) - 2(degree 2) - 3(degree 1)
        let adj_list = vec![
            vec![1],       // 0 (degree 1)
            vec![0, 2, 3], // 1 (degree 3)
            vec![1, 3],    // 2 (degree 2)
            vec![1, 2],    // 3 (degree 2)
        ];

        let order = fast_minimum_degree(&adj_list).unwrap();
        assert_eq!(order.len(), 4);
        
        // First eliminated should be degree 1 vertex
        assert!(order[0] == 0, "First eliminated should be degree 1 vertex");
        
        // All vertices should be present
        let mut sorted_order = order.clone();
        sorted_order.sort();
        assert_eq!(sorted_order, vec![0, 1, 2, 3]);
    }

    #[test]
    fn test_fast_minimum_degree_error_handling() {
        // Test various error conditions
        
        // Empty adjacency list
        let empty: Vec<Vec<usize>> = vec![];
        assert!(matches!(
            fast_minimum_degree(&empty),
            Err(FastMinimumDegreeError::EmptyAdjacencyList)
        ));

        // Invalid neighbor index
        let invalid = vec![vec![1], vec![0, 10]];
        assert!(matches!(
            fast_minimum_degree(&invalid),
            Err(FastMinimumDegreeError::NeighborIndexOutOfBounds { .. })
        ));

        // Self-reference (should be valid but ignored)
        let self_ref = vec![vec![0, 1], vec![0]];
        let result = fast_minimum_degree(&self_ref);
        assert!(result.is_ok(), "Self-reference should be handled gracefully");
    }
}