arcweight 0.3.0

//! Epsilon transition removal for weighted FSTs.
//!
//! This module implements epsilon removal, which eliminates $`\varepsilon`$-transitions
//! from a weighted FST while preserving the weighted language. The resulting FST has
//! no epsilon arcs but accepts the same input-output pairs with the same weights.
//!
//! # Mathematical Foundation
//!
//! For an FST $`T`$ with epsilon transitions, the epsilon-free equivalent $`T'`$
//! is constructed by computing epsilon closures:
//!
//! ```text
//! ε*(q) = {(p, w) : q →ε* p with total weight w}
//! ```
//!
//! For each non-epsilon arc $`q \xrightarrow{a:b/w} r`$, we add arcs
//! $`p \xrightarrow{a:b/w' \otimes w} r`$ for all $`(p, w') \in \varepsilon^*(q)`$.
//!
//! # Complexity
//!
//! | Operation | Time | Space |
//! |-----------|------|-------|
//! | Epsilon closure | $`O(V^2)`$ worst | $`O(V)`$ |
//! | Full algorithm | $`O(V^2 + VE)`$ | $`O(V^2)`$ |
//!
//! For sparse FSTs, the practical complexity is often much better.
//!
//! # Semiring Requirements
//!
//! Requires [`StarSemiring`] for handling epsilon cycles:
//!
//! | Semiring | Supported | Convergence |
//! |----------|-----------|-------------|
//! | [`TropicalWeight`] | Yes | Idempotent plus guarantees convergence |
//! | [`BooleanWeight`] | Yes | Trivial star: $`w^* = \bar{1}`$ |
//! | [`LogWeight`] | Yes | Converges for bounded cycles |
//! | [`ProbabilityWeight`] | Caution | May diverge for cycles with $`w \geq 1`$ |
//!
//! # Termination
//!
//! The algorithm terminates when:
//! - FST is epsilon-acyclic (no epsilon cycles)
//! - Semiring is k-closed ($`w^* = \bigoplus_{i=0}^{k} w^i`$ for some finite $`k`$)
//! - Epsilon cycles converge (idempotent semirings always do)
//!
//! # Algorithm
//!
//! 1. For each state $`q`$, compute epsilon closure using BFS
//! 2. For epsilon cycles, apply star operation to bound accumulated weight
//! 3. For each non-epsilon arc from closure states, create direct arc
//! 4. Propagate final weights through closures
//!
//! # Example
//!
//! ```rust
//! use arcweight::prelude::*;
//!
//! let mut fst = VectorFst::<BooleanWeight>::new();
//! let s0 = fst.add_state();
//! let s1 = fst.add_state();
//! let s2 = fst.add_state();
//! fst.set_start(s0);
//! fst.set_final(s2, BooleanWeight::one());
//!
//! // Path: s0 --ε--> s1 --a--> s2
//! fst.add_arc(s0, Arc::epsilon(BooleanWeight::one(), s1));
//! fst.add_arc(s1, Arc::new(1, 1, BooleanWeight::one(), s2));
//!
//! // Remove epsilon transitions
//! let no_eps: VectorFst<BooleanWeight> = remove_epsilons(&fst)?;
//!
//! // Now has direct arc: s0 --a--> s2
//! # Ok::<(), arcweight::Error>(())
//! ```
//!
//! # References
//!
//! \[1\] Mohri, M. 2002. Semiring frameworks and algorithms for shortest-distance
//!     problems. *Journal of Automata, Languages and Combinatorics* 7, 3, 321-350.
//!
//! \[2\] Mohri, M. 2009. Weighted automata algorithms. In *Handbook of Weighted
//!     Automata*, M. Droste, W. Kuich, and H. Vogler, Eds. Springer, 213-254.
//!     DOI: <https://doi.org/10.1007/978-3-642-01492-5_6>
//!
//! \[3\] Allauzen, C. and Mohri, M. 2006. Efficient algorithms for the composition
//!     of weighted finite-state transducers. In *Implementation and Application
//!     of Automata (CIAA 2006)*, Springer, 161-173.
//!
//! [`StarSemiring`]: crate::semiring::StarSemiring
//! [`TropicalWeight`]: crate::semiring::TropicalWeight
//! [`BooleanWeight`]: crate::semiring::BooleanWeight
//! [`LogWeight`]: crate::semiring::LogWeight
//! [`ProbabilityWeight`]: crate::semiring::ProbabilityWeight

use crate::arc::Arc;
use crate::fst::{Fst, MutableFst, StateId};
use crate::semiring::StarSemiring;
use crate::Result;
use rustc_hash::FxHashMap;
use std::collections::VecDeque;

/// Remove epsilon (empty) transitions from an FST while preserving language semantics
///
/// Eliminates all epsilon transitions by computing epsilon closures and creating
/// equivalent direct transitions. The resulting FST accepts the same language
/// but with more efficient processing due to eliminated epsilon transitions.
///
/// # Algorithm Details
///
/// - **Epsilon Closure:** For each state, compute all states reachable via epsilon transitions
/// - **Direct Transitions:** Create direct arcs bypassing epsilon paths
/// - **Time Complexity:** O(|V| × (|V| + |E|)) = O(|V|² + |V| × |E|) in worst case
/// - **Space Complexity:** O(|V|²) for storing epsilon closures
/// - **Language Preservation:** L(remove_epsilons(T)) = L(T) exactly
///
/// # Mathematical Foundation
///
/// For an FST with epsilon transitions, the epsilon-free equivalent computes:
/// - **Epsilon Closure:** ε*(q) = {p : q →ε* p} (states reachable via epsilon paths)
/// - **Direct Arcs:** For each non-epsilon arc q →a p, add arcs r →a p for all r ∈ ε*(q)
/// - **Final Weights:** Combine final weights through epsilon closures
/// - **Weight Computation:** Uses semiring operations to combine path weights
///
/// # Algorithm Steps
///
/// 1. **Copy Structure:** Create new FST with same states as original
/// 2. **Epsilon Closure:** For each state, compute epsilon-reachable states with weights
/// 3. **Direct Arcs:** Add direct arcs bypassing epsilon transitions
/// 4. **Final Weight Update:** Propagate final weights through epsilon closures
/// 5. **Clean Result:** Result FST has no epsilon transitions
///
/// # Examples
///
/// ## Basic Epsilon Removal
///
/// ```rust
/// use arcweight::prelude::*;
///
/// // FST with epsilon transitions: start --a--> s1 --ε--> s2 --b--> final
/// let mut fst = VectorFst::<BooleanWeight>::new();
/// let s0 = fst.add_state();
/// let s1 = fst.add_state();
/// let s2 = fst.add_state();
/// let s3 = fst.add_state();
///
/// fst.set_start(s0);
/// fst.set_final(s3, BooleanWeight::one());
///
/// fst.add_arc(s0, Arc::new('a' as u32, 'a' as u32, BooleanWeight::one(), s1));
/// fst.add_arc(s1, Arc::epsilon(BooleanWeight::one(), s2)); // epsilon transition
/// fst.add_arc(s2, Arc::new('b' as u32, 'b' as u32, BooleanWeight::one(), s3));
///
/// // Remove epsilon transitions
/// let no_eps: VectorFst<BooleanWeight> = remove_epsilons(&fst)?;
///
/// // Result accepts "ab" directly without intermediate epsilon
/// assert_eq!(no_eps.num_states(), fst.num_states());
/// # Ok::<(), Box<dyn std::error::Error>>(())
/// ```
///
/// ## Multiple Epsilon Paths
///
/// ```rust
/// use arcweight::prelude::*;
///
/// // FST with multiple epsilon paths
/// let mut fst = VectorFst::<BooleanWeight>::new();
/// let s0 = fst.add_state();
/// let s1 = fst.add_state();
/// let s2 = fst.add_state();
/// let s3 = fst.add_state();
///
/// fst.set_start(s0);
/// fst.set_final(s3, BooleanWeight::one());
///
/// // Multiple epsilon paths: s0 --ε--> s1 --ε--> s2 --ε--> s3
/// fst.add_arc(s0, Arc::epsilon(BooleanWeight::one(), s1));
/// fst.add_arc(s1, Arc::epsilon(BooleanWeight::one(), s2));
/// fst.add_arc(s2, Arc::epsilon(BooleanWeight::one(), s3));
///
/// // Also direct non-epsilon arc
/// fst.add_arc(s1, Arc::new('x' as u32, 'x' as u32, BooleanWeight::one(), s3));
///
/// // Remove all epsilon chains
/// let cleaned: VectorFst<BooleanWeight> = remove_epsilons(&fst)?;
///
/// // Result has direct paths without epsilon steps
/// println!("Cleaned FST ready");
/// # Ok::<(), Box<dyn std::error::Error>>(())
/// ```
///
/// ## Weighted Epsilon Removal
///
/// ```rust
/// use arcweight::prelude::*;
///
/// // Weighted FST with epsilon transitions
/// let mut fst = VectorFst::<BooleanWeight>::new();
/// let s0 = fst.add_state();
/// let s1 = fst.add_state();
/// let s2 = fst.add_state();
///
/// fst.set_start(s0);
/// fst.set_final(s2, BooleanWeight::one());
///
/// // Path: s0 --a--> s1 --ε--> s2(final)
/// fst.add_arc(s0, Arc::new('a' as u32, 'a' as u32, BooleanWeight::one(), s1));
/// fst.add_arc(s1, Arc::epsilon(BooleanWeight::one(), s2));
///
/// // Epsilon removal creates direct final weight for s1
/// let result: VectorFst<BooleanWeight> = remove_epsilons(&fst)?;
///
/// // s1 now has final weight due to epsilon closure to s2
/// assert!(result.start().is_some());
/// # Ok::<(), Box<dyn std::error::Error>>(())
/// ```
///
/// ## Complex Epsilon Network
///
/// ```rust
/// use arcweight::prelude::*;
///
/// // Complex epsilon transition network
/// let mut fst = VectorFst::<BooleanWeight>::new();
/// let states: Vec<_> = (0..5).map(|_| fst.add_state()).collect();
///
/// fst.set_start(states[0]);
/// fst.set_final(states[4], BooleanWeight::one());
///
/// // Mix of epsilon and non-epsilon transitions
/// fst.add_arc(states[0], Arc::new('a' as u32, 'a' as u32, BooleanWeight::one(), states[1]));
/// fst.add_arc(states[1], Arc::epsilon(BooleanWeight::one(), states[2]));
/// fst.add_arc(states[1], Arc::epsilon(BooleanWeight::one(), states[3]));
/// fst.add_arc(states[2], Arc::new('b' as u32, 'b' as u32, BooleanWeight::one(), states[4]));
/// fst.add_arc(states[3], Arc::new('c' as u32, 'c' as u32, BooleanWeight::one(), states[4]));
///
/// // Remove epsilon network, creating direct alternatives
/// let simplified: VectorFst<BooleanWeight> = remove_epsilons(&fst)?;
///
/// // Result accepts "ab" and "ac" directly
/// assert_eq!(simplified.num_states(), fst.num_states());
/// # Ok::<(), Box<dyn std::error::Error>>(())
/// ```
///
/// ## Optimization Pipeline Integration
///
/// ```rust
/// use arcweight::prelude::*;
///
/// // Optimization pipeline with epsilon removal
/// fn optimize_fst(fst: &VectorFst<BooleanWeight>) -> Result<VectorFst<BooleanWeight>> {
///     // Step 1: Remove epsilon transitions for efficiency
///     let no_epsilon: VectorFst<BooleanWeight> = remove_epsilons(fst)?;
///     
///     // Step 2: Remove unreachable states
///     let connected: VectorFst<BooleanWeight> = connect(&no_epsilon)?;
///     
///     // Step 3: Determinize if needed (epsilon-free FSTs determinize better)
///     // Note: determinize requires DivisibleSemiring, not available for BooleanWeight
///     // let determinized: VectorFst<BooleanWeight> = determinize(&connected)?;
///     
///     Ok(connected)
/// }
///
/// // Create test FST
/// let mut test_fst = VectorFst::new();
/// let s0 = test_fst.add_state();
/// let s1 = test_fst.add_state();
/// test_fst.set_start(s0);
/// test_fst.set_final(s1, BooleanWeight::one());
/// test_fst.add_arc(s0, Arc::epsilon(BooleanWeight::one(), s1));
///
/// let optimized = optimize_fst(&test_fst)?;
/// println!("FST optimized");
/// # Ok::<(), Box<dyn std::error::Error>>(())
/// ```
///
/// # Use Cases
///
/// ## FST Preprocessing
/// - **Algorithm Preparation:** Remove epsilons before determinization or minimization
/// - **Performance Optimization:** Eliminate epsilon transitions for faster traversal
/// - **Canonical Form:** Create epsilon-free canonical representation
/// - **Memory Efficiency:** Reduce state space complexity
///
/// ## Regular Expression Processing
/// - **Regex Compilation:** Remove epsilons from regex-derived FSTs
/// - **Pattern Matching:** Optimize pattern matching automata
/// - **Text Processing:** Eliminate empty transitions in text processors
/// - **Lexical Analysis:** Clean up tokenizer automata
///
/// ## Natural Language Processing
/// - **Grammar Cleanup:** Remove epsilon productions from CFG-derived FSTs
/// - **Morphological Analysis:** Clean morphology automata
/// - **Phonological Rules:** Optimize phonology rule FSTs
/// - **Translation Models:** Simplify translation automata
///
/// ## Speech Processing
/// - **Pronunciation Models:** Remove epsilon paths in pronunciation FSTs
/// - **Acoustic Models:** Optimize acoustic model automata
/// - **Language Models:** Clean statistical language model FSTs
/// - **ASR Optimization:** Prepare FSTs for speech recognition
///
/// # Performance Characteristics
///
/// - **Time Complexity:** O(|V| × (|V| + |E|)) worst case for dense epsilon networks
/// - **Space Complexity:** O(|V|²) for storing epsilon closure information
/// - **Practical Performance:** Often much better than worst case for sparse graphs
/// - **Memory Access:** Sequential state processing improves cache locality
/// - **Convergence:** Guaranteed for k-closed semirings with proper star operation
/// - **Dense Graph Analysis:** For |E| = Θ(|V|²), complexity is O(|V|³)
///
/// # Mathematical Properties
///
/// Epsilon removal preserves fundamental FST properties:
/// - **Language Preservation:** L(remove_epsilons(T)) = L(T) exactly
/// - **Weight Preservation:** All path weights maintained through semiring operations
/// - **Structural Properties:** May change state connectivity but preserves semantics
/// - **Determinism:** Epsilon-free FSTs are easier to determinize efficiently
/// - **Compositionality:** remove_epsilons(T₁ ∘ T₂) relates to remove_epsilons(T₁) ∘ remove_epsilons(T₂)
///
/// # Implementation Details
///
/// The algorithm computes epsilon closures using breadth-first search with weight
/// accumulation. For each state, it finds all epsilon-reachable states and their
/// combined weights, then creates direct arcs that bypass epsilon paths.
///
/// Cycle handling relies on the semiring's plus operation to combine multiple
/// paths to the same state, ensuring convergence when the semiring is k-closed.
///
/// # Semiring Considerations
///
/// Different semirings have different epsilon removal characteristics:
/// - **Boolean Semiring:** Epsilon removal is always efficient and terminates
/// - **Tropical Semiring:** Idempotent plus ensures convergence
/// - **Probability Semiring:** May require careful handling of convergence
/// - **Log Semiring:** Convergence depends on cycle weight properties
///
/// # Optimization Opportunities
///
/// After epsilon removal, consider these optimizations:
/// - **Determinization:** Now more efficient without epsilon transitions
/// - **Minimization:** Can reduce states further after epsilon elimination
/// - **Connection:** Remove states that became unreachable
/// - **Topological Sort:** Order states for optimal processing
///
/// # Errors
///
/// Returns [`Error::Algorithm`](crate::Error::Algorithm) if:
/// - The input FST is invalid, corrupted, or malformed
/// - Memory allocation fails during epsilon closure computation
/// - Epsilon closure computation encounters infinite loops or non-convergent weights
/// - The semiring does not properly support required star operations
/// - Weight computation overflows or produces invalid results
///
/// # References
///
/// \[1\] Mohri, M. 2002. Semiring frameworks and algorithms for shortest-distance
///     problems. *Journal of Automata, Languages and Combinatorics* 7, 3, 321-350.
///
/// \[2\] Mohri, M. 2009. Weighted automata algorithms. In *Handbook of Weighted
///     Automata*, Springer, 213-254. DOI: <https://doi.org/10.1007/978-3-642-01492-5_6>
///
/// # See Also
///
/// - [`determinize`](crate::algorithms::determinize()) - Often applied after epsilon removal
/// - [`minimize`](crate::algorithms::minimize()) - State reduction after epsilon removal
/// - [`connect`](crate::algorithms::connect()) - Remove unreachable states
/// - [`StarSemiring`] - Required trait for epsilon cycles
pub fn remove_epsilons<W, F, M>(fst: &F) -> Result<M>
where
    W: StarSemiring,
    F: Fst<W>,
    M: MutableFst<W> + Default,
{
    let mut result = M::default();

    // copy states
    for _ in 0..fst.num_states() {
        result.add_state();
    }

    // set start
    if let Some(start) = fst.start() {
        result.set_start(start);
    }

    // compute epsilon closure for each state
    for state in fst.states() {
        let closure = compute_epsilon_closure(fst, state)?;

        // add non-epsilon arcs
        for arc in fst.arcs(state) {
            if !arc.is_epsilon() {
                result.add_arc(state, arc.clone());
            }
        }

        // Initialize final weight for this state if it has one
        let mut accumulated_final_weight = fst.final_weight(state).cloned();

        // add arcs from epsilon closure
        for &(closure_state, ref weight) in &closure {
            if closure_state != state {
                // add non-epsilon arcs from closure state
                for arc in fst.arcs(closure_state) {
                    if !arc.is_epsilon() {
                        result.add_arc(
                            state,
                            Arc::new(
                                arc.ilabel,
                                arc.olabel,
                                weight.times(&arc.weight),
                                arc.nextstate,
                            ),
                        );
                    }
                }

                // handle final weights from epsilon closure
                if let Some(final_weight) = fst.final_weight(closure_state) {
                    let propagated_weight = weight.times(final_weight);
                    accumulated_final_weight = match accumulated_final_weight {
                        Some(existing) => Some(existing.plus(&propagated_weight)),
                        None => Some(propagated_weight),
                    };
                }
            }
        }

        // Set the accumulated final weight if any
        if let Some(final_weight) = accumulated_final_weight {
            result.set_final(state, final_weight);
        }
    }

    Ok(result)
}

/// Find strongly connected components in epsilon subgraph using Tarjan's algorithm
fn find_epsilon_sccs<W: StarSemiring, F: Fst<W>>(
    fst: &F,
) -> (
    FxHashMap<StateId, StateId>,
    FxHashMap<StateId, Vec<StateId>>,
) {
    let n = fst.num_states();
    let mut index_counter = 0;
    let mut stack = Vec::new();
    let mut indices = vec![None; n];
    let mut lowlinks = vec![None; n];
    let mut on_stack = vec![false; n];
    let mut scc_map = FxHashMap::default();
    let mut scc_id = 0;
    let mut scc_states: FxHashMap<StateId, Vec<StateId>> = FxHashMap::default();

    // Run DFS from each unvisited state
    for state_idx in 0..n {
        let state = state_idx as StateId;
        if indices[state_idx].is_none() {
            tarjan_epsilon_dfs(
                fst,
                state,
                &mut index_counter,
                &mut stack,
                &mut indices,
                &mut lowlinks,
                &mut on_stack,
                &mut scc_map,
                &mut scc_id,
                &mut scc_states,
            );
        }
    }

    (scc_map, scc_states)
}

/// Tarjan's DFS helper for epsilon subgraph
#[allow(clippy::too_many_arguments)]
fn tarjan_epsilon_dfs<W: StarSemiring, F: Fst<W>>(
    fst: &F,
    v: StateId,
    index_counter: &mut usize,
    stack: &mut Vec<StateId>,
    indices: &mut Vec<Option<usize>>,
    lowlinks: &mut Vec<Option<usize>>,
    on_stack: &mut Vec<bool>,
    scc_map: &mut FxHashMap<StateId, StateId>,
    scc_id: &mut StateId,
    scc_states: &mut FxHashMap<StateId, Vec<StateId>>,
) {
    let v_idx = v as usize;

    // Set the depth index for v
    indices[v_idx] = Some(*index_counter);
    lowlinks[v_idx] = Some(*index_counter);
    *index_counter += 1;
    stack.push(v);
    on_stack[v_idx] = true;

    // Consider successors of v (only epsilon arcs)
    for arc in fst.arcs(v) {
        if arc.is_epsilon() {
            let w = arc.nextstate;
            let w_idx = w as usize;
            if indices[w_idx].is_none() {
                // Successor w has not yet been visited; recurse on it
                tarjan_epsilon_dfs(
                    fst,
                    w,
                    index_counter,
                    stack,
                    indices,
                    lowlinks,
                    on_stack,
                    scc_map,
                    scc_id,
                    scc_states,
                );
                lowlinks[v_idx] = Some(lowlinks[v_idx].unwrap().min(lowlinks[w_idx].unwrap()));
            } else if on_stack[w_idx] {
                // Successor w is in stack and hence in the current SCC
                lowlinks[v_idx] = Some(lowlinks[v_idx].unwrap().min(indices[w_idx].unwrap()));
            }
        }
    }

    // If v is a root node, pop the stack and create an SCC
    if lowlinks[v_idx] == indices[v_idx] {
        let current_scc = *scc_id;
        *scc_id += 1;
        let mut scc = Vec::new();

        loop {
            let w = stack.pop().unwrap();
            let w_idx = w as usize;
            on_stack[w_idx] = false;
            scc_map.insert(w, current_scc);
            scc.push(w);
            if w == v {
                break;
            }
        }
        scc_states.insert(current_scc, scc);
    }
}

/// Compute cycle weight for an SCC in epsilon subgraph
/// Returns the star of the cycle weight if the SCC is cyclic
fn compute_scc_cycle_weight<W: StarSemiring, F: Fst<W>>(
    fst: &F,
    scc_states: &[StateId],
) -> Option<W> {
    if scc_states.is_empty() {
        return None;
    }

    // For single state, check for self-loop
    if scc_states.len() == 1 {
        let state = scc_states[0];
        for arc in fst.arcs(state) {
            if arc.is_epsilon() && arc.nextstate == state {
                // Self-loop: apply star operation
                return Some(arc.weight.star());
            }
        }
        return None; // No self-loop, not a cycle
    }

    // For multi-state SCC, find a simple cycle
    // Use the first state as starting point
    let start = scc_states[0];
    find_cycle_weight_from_state(fst, start, scc_states).map(|w| w.star())
}

/// Find the weight of a cycle starting from a given state within an SCC
fn find_cycle_weight_from_state<W: StarSemiring, F: Fst<W>>(
    fst: &F,
    start: StateId,
    scc_states: &[StateId],
) -> Option<W> {
    // Use BFS to find shortest cycle back to start
    let mut queue = VecDeque::new();
    // Track (state, accumulated_weight, path_length)
    let mut visited: FxHashMap<(StateId, usize), W> = FxHashMap::default();

    // Start BFS from start state's neighbors
    for arc in fst.arcs(start) {
        if arc.is_epsilon() && scc_states.contains(&arc.nextstate) {
            let next = arc.nextstate;
            if next == start {
                // Direct self-loop (should have been caught in single-state case, but handle it)
                return Some(arc.weight.clone());
            }
            queue.push_back((next, arc.weight.clone(), 1));
            visited.insert((next, 1), arc.weight.clone());
        }
    }

    // Limit search depth to prevent infinite loops.
    // A simple cycle can have at most |SCC| edges, but complex graphs may need
    // more exploration to find all cycle paths. Use 2*|SCC| as a conservative limit.
    let max_depth = scc_states.len() * 2;

    while let Some((state, weight, depth)) = queue.pop_front() {
        if depth > max_depth {
            continue;
        }

        for arc in fst.arcs(state) {
            if arc.is_epsilon() && scc_states.contains(&arc.nextstate) {
                let next = arc.nextstate;
                let next_weight = weight.times(&arc.weight);
                let next_depth = depth + 1;

                // If we've returned to start, we found a cycle
                if next == start {
                    return Some(next_weight);
                }

                // Continue BFS if we haven't seen this state at this depth with better weight
                let key = (next, next_depth);
                let should_explore = match visited.get(&key) {
                    Some(existing) => {
                        let combined = existing.plus(&next_weight);
                        if combined != *existing {
                            visited.insert(key, combined.clone());
                            true
                        } else {
                            false
                        }
                    }
                    None => {
                        visited.insert(key, next_weight.clone());
                        true
                    }
                };

                if should_explore && next_depth <= max_depth {
                    queue.push_back((next, visited[&key].clone(), next_depth));
                }
            }
        }
    }

    None
}

fn compute_epsilon_closure<W: StarSemiring, F: Fst<W>>(
    fst: &F,
    start: StateId,
) -> Result<Vec<(StateId, W)>> {
    // Find SCCs in epsilon subgraph for cycle detection
    let (scc_map, scc_states_map) = find_epsilon_sccs(fst);

    // Compute star weights for cyclic SCCs
    let mut scc_star_weights: FxHashMap<StateId, W> = FxHashMap::default();
    for (scc_id, states) in &scc_states_map {
        if let Some(star_weight) = compute_scc_cycle_weight(fst, states) {
            scc_star_weights.insert(*scc_id, star_weight);
        }
    }

    let mut closure = FxHashMap::default();
    let mut queue = VecDeque::new();

    queue.push_back((start, W::one()));
    closure.insert(start, W::one());

    while let Some((state, weight)) = queue.pop_front() {
        // Update closure with current weight
        let mut current_best = match closure.get(&state) {
            Some(existing) => {
                let combined = existing.plus(&weight);
                if combined != *existing {
                    // Weight improved, update closure
                    closure.insert(state, combined.clone());
                    combined
                } else {
                    // Weight didn't improve - for idempotent semirings, this means
                    // we've already seen this path. But we should still process arcs
                    // in case there are new arcs or the state was reached via a different path.
                    // However, to avoid infinite loops, we only process if this is a new
                    // state-weight combination we haven't fully explored.
                    // For now, always process arcs - the should_update check will prevent
                    // adding redundant states to queue
                    existing.clone()
                }
            }
            None => {
                // New state, add to closure
                closure.insert(state, weight.clone());
                weight.clone()
            }
        };

        // Handle cycles: if state is in a cyclic SCC, account for cycle traversal
        // Only apply if this SCC actually has a cycle (not just an SCC)
        let state_scc = scc_map.get(&state).copied();
        if let Some(scc_id) = state_scc {
            // Check if this is a cyclic SCC (has cycle weight computed)
            if let Some(star_weight) = scc_star_weights.get(&scc_id) {
                // We can traverse the cycle zero or more times
                // Weight with cycle: current_best * star(cycle_weight)
                let cycle_enhanced = current_best.times(star_weight);
                // Take the best of direct path and cycle path
                let best_with_cycle = current_best.plus(&cycle_enhanced);
                if best_with_cycle != current_best {
                    closure.insert(state, best_with_cycle.clone());
                    current_best = best_with_cycle;
                }
            }
        }

        // Follow epsilon transitions (normal BFS)
        // Use current_best which accounts for cycles
        for arc in fst.arcs(state) {
            if arc.is_epsilon() {
                let next_state = arc.nextstate;
                let next_weight = current_best.times(&arc.weight);

                let should_update = match closure.get(&next_state) {
                    Some(existing) => {
                        let combined = existing.plus(&next_weight);
                        if combined != *existing {
                            closure.insert(next_state, combined.clone());
                            true
                        } else {
                            false
                        }
                    }
                    None => {
                        closure.insert(next_state, next_weight.clone());
                        true
                    }
                };

                if should_update {
                    queue.push_back((next_state, closure[&next_state].clone()));
                }
            }
        }
    }

    // Convert FxHashMap to Vec for compatibility
    let mut result: Vec<(StateId, W)> = closure.into_iter().collect();
    result.sort_by_key(|(state, _)| *state);
    Ok(result)
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::prelude::*;
    use num_traits::One;

    #[test]
    fn test_remove_epsilons() {
        // Test epsilon removal with BooleanWeight
        let mut bool_fst = VectorFst::<BooleanWeight>::new();
        let s0 = bool_fst.add_state();
        let s1 = bool_fst.add_state();
        let s2 = bool_fst.add_state();

        bool_fst.set_start(s0);
        bool_fst.set_final(s2, BooleanWeight::one());

        bool_fst.add_arc(s0, Arc::epsilon(BooleanWeight::new(true), s1));
        bool_fst.add_arc(s1, Arc::new(1, 1, BooleanWeight::new(true), s2));

        let no_eps: VectorFst<BooleanWeight> =
            remove_epsilons::<BooleanWeight, VectorFst<BooleanWeight>, VectorFst<BooleanWeight>>(
                &bool_fst,
            )
            .unwrap();

        // Check no epsilon transitions remain
        for state in no_eps.states() {
            for arc in no_eps.arcs(state) {
                assert!(!arc.is_epsilon(), "Found epsilon arc: {arc:?}");
            }
        }

        // Should preserve the language
        assert!(no_eps.start().is_some());
    }

    #[test]
    fn test_remove_epsilons_none() {
        // Test epsilon removal with BooleanWeight
        let mut bool_fst = VectorFst::<BooleanWeight>::new();
        let s0 = bool_fst.add_state();
        let s1 = bool_fst.add_state();
        let s2 = bool_fst.add_state();

        bool_fst.set_start(s0);
        bool_fst.set_final(s2, BooleanWeight::one());

        bool_fst.add_arc(s0, Arc::epsilon(BooleanWeight::new(true), s1));
        bool_fst.add_arc(s1, Arc::new(1, 1, BooleanWeight::new(true), s2));

        let no_eps: VectorFst<BooleanWeight> =
            remove_epsilons::<BooleanWeight, VectorFst<BooleanWeight>, VectorFst<BooleanWeight>>(
                &bool_fst,
            )
            .unwrap();

        // Should preserve structure when no epsilons present
        // Note: Implementation might add states for proper structure
        assert!(no_eps.num_states() >= bool_fst.num_states() - 1); // Allow for structure changes
    }

    #[test]
    fn test_remove_epsilons_multiple_paths() {
        // Test epsilon removal with multiple epsilon paths to same state
        let mut fst = VectorFst::<BooleanWeight>::new();
        let s0 = fst.add_state();
        let s1 = fst.add_state();
        let s2 = fst.add_state();
        let s3 = fst.add_state();

        fst.set_start(s0);
        fst.set_final(s3, BooleanWeight::one());

        // Multiple epsilon paths from s0 to s2
        fst.add_arc(s0, Arc::epsilon(BooleanWeight::one(), s1));
        fst.add_arc(s1, Arc::epsilon(BooleanWeight::one(), s2));
        fst.add_arc(s0, Arc::epsilon(BooleanWeight::one(), s2)); // direct path

        // Non-epsilon arc from s2
        fst.add_arc(s2, Arc::new(1, 1, BooleanWeight::one(), s3));

        let result: VectorFst<BooleanWeight> = remove_epsilons(&fst).unwrap();

        // Should have direct arc from s0 to s3
        let arcs_from_start: Vec<_> = result.arcs(s0).collect();
        assert!(arcs_from_start
            .iter()
            .any(|a| a.ilabel == 1 && a.nextstate == s3));

        // No epsilon transitions should remain
        for state in result.states() {
            for arc in result.arcs(state) {
                assert!(!arc.is_epsilon());
            }
        }
    }

    #[test]
    fn test_remove_epsilons_cycles() {
        // Test epsilon removal with epsilon cycles
        let mut fst = VectorFst::<BooleanWeight>::new();
        let s0 = fst.add_state();
        let s1 = fst.add_state();
        let s2 = fst.add_state();
        let s3 = fst.add_state();

        fst.set_start(s0);
        fst.set_final(s3, BooleanWeight::one());

        // Create epsilon cycle: s0 -> s1 -> s2 -> s0
        fst.add_arc(s0, Arc::epsilon(BooleanWeight::one(), s1));
        fst.add_arc(s1, Arc::epsilon(BooleanWeight::one(), s2));
        fst.add_arc(s2, Arc::epsilon(BooleanWeight::one(), s0)); // cycle

        // Exit from cycle
        fst.add_arc(s1, Arc::new(1, 1, BooleanWeight::one(), s3));

        let result: VectorFst<BooleanWeight> = remove_epsilons(&fst).unwrap();

        // Should handle epsilon cycles correctly
        assert!(result.start().is_some());
        assert!(result.is_final(s3));

        // Should have path from start to final
        let arcs_from_start: Vec<_> = result.arcs(s0).collect();
        assert!(!arcs_from_start.is_empty());
    }

    #[test]
    fn test_remove_epsilons_to_final() {
        // Test epsilon transitions to final states
        let mut fst = VectorFst::<BooleanWeight>::new();
        let s0 = fst.add_state();
        let s1 = fst.add_state();
        let s2 = fst.add_state();

        fst.set_start(s0);
        fst.set_final(s2, BooleanWeight::one());

        // Path with epsilon to final
        fst.add_arc(s0, Arc::new(1, 1, BooleanWeight::one(), s1));
        fst.add_arc(s1, Arc::epsilon(BooleanWeight::one(), s2)); // epsilon to final

        let result: VectorFst<BooleanWeight> = remove_epsilons(&fst).unwrap();

        // s1 should become final since it has epsilon path to s2
        assert!(result.is_final(s1));
        assert!(result.is_final(s2));
    }

    #[test]
    fn test_remove_epsilons_all_epsilon() {
        // Test FST with only epsilon transitions
        let mut fst = VectorFst::<BooleanWeight>::new();
        let s0 = fst.add_state();
        let s1 = fst.add_state();
        let s2 = fst.add_state();

        fst.set_start(s0);
        fst.set_final(s2, BooleanWeight::one());

        // All transitions are epsilon
        fst.add_arc(s0, Arc::epsilon(BooleanWeight::one(), s1));
        fst.add_arc(s1, Arc::epsilon(BooleanWeight::one(), s2));

        let result: VectorFst<BooleanWeight> = remove_epsilons(&fst).unwrap();

        // Start state should become final (epsilon path to final)
        assert!(result.is_final(s0));

        // Should have no arcs (all were epsilon)
        let total_arcs: usize = result.states().map(|s| result.num_arcs(s)).sum();
        assert_eq!(total_arcs, 0);
    }

    #[test]
    fn test_remove_epsilons_mixed_paths() {
        // Test with mixed epsilon and non-epsilon paths
        let mut fst = VectorFst::<BooleanWeight>::new();
        let states: Vec<_> = (0..5).map(|_| fst.add_state()).collect();

        fst.set_start(states[0]);
        fst.set_final(states[4], BooleanWeight::one());

        // Mixed paths
        fst.add_arc(states[0], Arc::new(1, 1, BooleanWeight::one(), states[1]));
        fst.add_arc(states[1], Arc::epsilon(BooleanWeight::one(), states[2]));
        fst.add_arc(states[2], Arc::new(2, 2, BooleanWeight::one(), states[3]));
        fst.add_arc(states[3], Arc::epsilon(BooleanWeight::one(), states[4]));

        // Alternative path
        fst.add_arc(states[0], Arc::epsilon(BooleanWeight::one(), states[2]));

        let result: VectorFst<BooleanWeight> = remove_epsilons(&fst).unwrap();

        // Should have no epsilon transitions
        for state in result.states() {
            for arc in result.arcs(state) {
                assert!(!arc.is_epsilon());
            }
        }

        // Should preserve language
        assert!(result.is_final(states[4]));
    }

    #[test]
    fn test_epsilon_closure_computation() {
        // Test the epsilon closure computation directly
        let mut fst = VectorFst::<BooleanWeight>::new();
        let s0 = fst.add_state();
        let s1 = fst.add_state();
        let s2 = fst.add_state();
        let s3 = fst.add_state();

        fst.set_start(s0);
        fst.set_final(s3, BooleanWeight::one());

        // Epsilon transitions forming a DAG
        fst.add_arc(s0, Arc::epsilon(BooleanWeight::one(), s1));
        fst.add_arc(s0, Arc::epsilon(BooleanWeight::one(), s2));
        fst.add_arc(s1, Arc::epsilon(BooleanWeight::one(), s3));
        fst.add_arc(s2, Arc::epsilon(BooleanWeight::one(), s3));

        let closure = compute_epsilon_closure(&fst, s0).unwrap();

        // Should reach all states from s0
        let reached_states: Vec<_> = closure.iter().map(|(state, _)| *state).collect();
        assert!(reached_states.contains(&s0));
        assert!(reached_states.contains(&s1));
        assert!(reached_states.contains(&s2));
        assert!(reached_states.contains(&s3));
    }

    #[test]
    fn test_remove_epsilons_preserves_weights() {
        // Test that epsilon removal preserves path weights correctly
        let mut fst = VectorFst::<BooleanWeight>::new();
        let s0 = fst.add_state();
        let s1 = fst.add_state();
        let s2 = fst.add_state();
        let s3 = fst.add_state();

        fst.set_start(s0);
        fst.set_final(s3, BooleanWeight::one());

        // Path with weights
        fst.add_arc(s0, Arc::new(1, 1, BooleanWeight::new(true), s1));
        fst.add_arc(s1, Arc::epsilon(BooleanWeight::new(true), s2));
        fst.add_arc(s2, Arc::new(2, 2, BooleanWeight::new(true), s3));

        let result: VectorFst<BooleanWeight> = remove_epsilons(&fst).unwrap();

        // Should have combined the path correctly
        let has_direct_path = result.arcs(s1).any(|a| a.ilabel == 2 && a.nextstate == s3);
        assert!(has_direct_path);
    }

    #[test]
    fn test_remove_epsilons_empty_fst() {
        // Test epsilon removal on empty FST
        let fst = VectorFst::<BooleanWeight>::new();
        let result: VectorFst<BooleanWeight> = remove_epsilons(&fst).unwrap();

        assert_eq!(result.num_states(), 0);
        assert!(result.start().is_none());
    }

    #[test]
    fn test_remove_epsilons_single_state() {
        // Test epsilon removal on single-state FST
        let mut fst = VectorFst::<BooleanWeight>::new();
        let s0 = fst.add_state();

        fst.set_start(s0);
        fst.set_final(s0, BooleanWeight::one());

        // Self-loop epsilon
        fst.add_arc(s0, Arc::epsilon(BooleanWeight::one(), s0));

        let result: VectorFst<BooleanWeight> = remove_epsilons(&fst).unwrap();

        assert_eq!(result.num_states(), 1);
        assert!(result.is_final(s0));

        // Self-loop epsilon should be removed
        let self_loops: Vec<_> = result.arcs(s0).filter(|a| a.is_epsilon()).collect();
        assert!(self_loops.is_empty());
    }

    #[test]
    fn test_remove_epsilons_tropical_weight() {
        // Test that remove_epsilons works with TropicalWeight now that it implements StarSemiring
        let mut fst = VectorFst::<TropicalWeight>::new();
        let s0 = fst.add_state();
        let s1 = fst.add_state();
        let s2 = fst.add_state();
        let s3 = fst.add_state();

        fst.set_start(s0);
        fst.set_final(s3, TropicalWeight::new(1.0));

        // Add epsilon transitions with weights
        fst.add_arc(s0, Arc::epsilon(TropicalWeight::new(0.5), s1));
        fst.add_arc(s0, Arc::epsilon(TropicalWeight::new(0.3), s2));

        // Add regular transitions
        fst.add_arc(s1, Arc::new(1, 1, TropicalWeight::new(2.0), s3));
        fst.add_arc(s2, Arc::new(1, 1, TropicalWeight::new(1.5), s3));

        // Remove epsilons - should preserve weight semantics
        let result: VectorFst<TropicalWeight> = remove_epsilons(&fst).unwrap();

        // Verify no epsilons remain
        for state in result.states() {
            for arc in result.arcs(state) {
                assert!(!arc.is_epsilon(), "Epsilon arc found after removal");
            }
        }

        // Verify structure is preserved (should have direct paths from s0 to s3)
        assert!(result.num_states() >= fst.num_states());
        assert!(result.start().is_some());
    }

    #[test]
    fn test_epsilon_closure_with_self_loop() {
        // Test epsilon closure with self-loop (single state cycle)
        let mut fst = VectorFst::<TropicalWeight>::new();
        let s0 = fst.add_state();
        let s1 = fst.add_state();
        let s2 = fst.add_state();

        fst.set_start(s0);
        fst.set_final(s2, TropicalWeight::one());

        // Self-loop epsilon on s0
        fst.add_arc(s0, Arc::epsilon(TropicalWeight::new(0.5), s0));
        // Path from s0 to s1
        fst.add_arc(s0, Arc::epsilon(TropicalWeight::new(1.0), s1));
        // Path from s1 to s2
        fst.add_arc(s1, Arc::new(1, 1, TropicalWeight::new(2.0), s2));

        let closure = compute_epsilon_closure(&fst, s0).unwrap();

        // Should include s0 (self), s1 (via epsilon), and potentially s0 again via cycle
        let reached_states: Vec<_> = closure.iter().map(|(state, _)| *state).collect();
        assert!(reached_states.contains(&s0));
        assert!(reached_states.contains(&s1));

        // Remove epsilons and verify correctness
        let result: VectorFst<TropicalWeight> = remove_epsilons(&fst).unwrap();
        assert!(result.start().is_some());
        // Should have path from s0 to s2
        let arcs_from_s0: Vec<_> = result.arcs(s0).collect();
        assert!(arcs_from_s0
            .iter()
            .any(|a| a.ilabel == 1 && a.nextstate == s2));
    }

    #[test]
    fn test_epsilon_closure_with_multi_state_cycle() {
        // Test epsilon closure with multi-state cycle
        let mut fst = VectorFst::<TropicalWeight>::new();
        let s0 = fst.add_state();
        let s1 = fst.add_state();
        let s2 = fst.add_state();
        let s3 = fst.add_state();

        fst.set_start(s0);
        fst.set_final(s3, TropicalWeight::one());

        // Cycle: s0 -> s1 -> s2 -> s0 (all epsilon)
        fst.add_arc(s0, Arc::epsilon(TropicalWeight::new(0.5), s1));
        fst.add_arc(s1, Arc::epsilon(TropicalWeight::new(0.3), s2));
        fst.add_arc(s2, Arc::epsilon(TropicalWeight::new(0.2), s0));
        // Exit from cycle
        fst.add_arc(s1, Arc::new(1, 1, TropicalWeight::new(1.0), s3));

        let closure = compute_epsilon_closure(&fst, s0).unwrap();

        // Should reach all states in cycle
        let reached_states: Vec<_> = closure.iter().map(|(state, _)| *state).collect();
        assert!(reached_states.contains(&s0));
        assert!(reached_states.contains(&s1));
        assert!(reached_states.contains(&s2));

        // Remove epsilons
        let result: VectorFst<TropicalWeight> = remove_epsilons(&fst).unwrap();
        assert!(result.start().is_some());
        // Should have path from s0 to s3 (through cycle)
        let arcs_from_s0: Vec<_> = result.arcs(s0).collect();
        assert!(arcs_from_s0
            .iter()
            .any(|a| a.ilabel == 1 && a.nextstate == s3));
    }

    #[test]
    fn test_epsilon_closure_boolean_weight_cycle() {
        // Test epsilon closure with BooleanWeight (trivial star operation)
        let mut fst = VectorFst::<BooleanWeight>::new();
        let s0 = fst.add_state();
        let s1 = fst.add_state();
        let s2 = fst.add_state();

        fst.set_start(s0);
        fst.set_final(s2, BooleanWeight::one());

        // Cycle: s0 -> s1 -> s0
        fst.add_arc(s0, Arc::epsilon(BooleanWeight::one(), s1));
        fst.add_arc(s1, Arc::epsilon(BooleanWeight::one(), s0));
        // Exit
        fst.add_arc(s0, Arc::new(1, 1, BooleanWeight::one(), s2));

        let closure = compute_epsilon_closure(&fst, s0).unwrap();

        // Should reach s0 and s1
        let reached_states: Vec<_> = closure.iter().map(|(state, _)| *state).collect();
        assert!(reached_states.contains(&s0));
        assert!(reached_states.contains(&s1));

        let result: VectorFst<BooleanWeight> = remove_epsilons(&fst).unwrap();
        assert!(result.start().is_some());
        assert!(result.is_final(s2));
    }

    #[test]
    fn test_epsilon_closure_multiple_cycles() {
        // Test FST with multiple epsilon cycles
        let mut fst = VectorFst::<TropicalWeight>::new();
        let s0 = fst.add_state();
        let s1 = fst.add_state();
        let s2 = fst.add_state();
        let s3 = fst.add_state();
        let s4 = fst.add_state();

        fst.set_start(s0);
        fst.set_final(s4, TropicalWeight::one());

        // First cycle: s0 -> s1 -> s0
        fst.add_arc(s0, Arc::epsilon(TropicalWeight::new(0.5), s1));
        fst.add_arc(s1, Arc::epsilon(TropicalWeight::new(0.3), s0));
        // Second cycle: s2 -> s3 -> s2
        fst.add_arc(s2, Arc::epsilon(TropicalWeight::new(0.2), s3));
        fst.add_arc(s3, Arc::epsilon(TropicalWeight::new(0.4), s2));
        // Connect cycles
        fst.add_arc(s0, Arc::epsilon(TropicalWeight::new(1.0), s2));
        // Exit
        fst.add_arc(s2, Arc::new(1, 1, TropicalWeight::new(2.0), s4));

        let result: VectorFst<TropicalWeight> = remove_epsilons(&fst).unwrap();
        assert!(result.start().is_some());
        // Should have path from s0 to s4
        let arcs_from_s0: Vec<_> = result.arcs(s0).collect();
        assert!(arcs_from_s0
            .iter()
            .any(|a| a.ilabel == 1 && a.nextstate == s4));
    }

    #[test]
    fn test_epsilon_closure_acyclic_path() {
        // Test that acyclic epsilon paths still work correctly
        let mut fst = VectorFst::<TropicalWeight>::new();
        let s0 = fst.add_state();
        let s1 = fst.add_state();
        let s2 = fst.add_state();
        let s3 = fst.add_state();

        fst.set_start(s0);
        fst.set_final(s3, TropicalWeight::one());

        // Acyclic epsilon path: s0 -> s1 -> s2 (no cycles)
        fst.add_arc(s0, Arc::epsilon(TropicalWeight::new(0.5), s1));
        fst.add_arc(s1, Arc::epsilon(TropicalWeight::new(0.3), s2));
        fst.add_arc(s2, Arc::new(1, 1, TropicalWeight::new(1.0), s3));

        let closure = compute_epsilon_closure(&fst, s0).unwrap();

        // Should reach s0, s1, s2
        let reached_states: Vec<_> = closure.iter().map(|(state, _)| *state).collect();
        assert!(reached_states.contains(&s0));
        assert!(reached_states.contains(&s1));
        assert!(reached_states.contains(&s2));

        let result: VectorFst<TropicalWeight> = remove_epsilons(&fst).unwrap();
        // Should have direct path from s0 to s3
        let arcs_from_s0: Vec<_> = result.arcs(s0).collect();
        assert!(arcs_from_s0
            .iter()
            .any(|a| a.ilabel == 1 && a.nextstate == s3));
    }

    #[test]
    fn test_epsilon_closure_termination() {
        // Test that epsilon closure terminates even with cycles
        // This is especially important for non-idempotent semirings
        let mut fst = VectorFst::<TropicalWeight>::new();
        let s0 = fst.add_state();
        let s1 = fst.add_state();

        fst.set_start(s0);
        fst.set_final(s1, TropicalWeight::one());

        // Large cycle that would cause infinite loop without star operation
        fst.add_arc(s0, Arc::epsilon(TropicalWeight::new(0.1), s0)); // self-loop
        fst.add_arc(s0, Arc::new(1, 1, TropicalWeight::new(1.0), s1));

        // This should terminate quickly with star operation
        let closure = compute_epsilon_closure(&fst, s0).unwrap();
        assert!(!closure.is_empty());

        let result: VectorFst<TropicalWeight> = remove_epsilons(&fst).unwrap();
        assert!(result.start().is_some());
    }
}