ruvector_mincut/algorithm/

approximate.rs

//! Approximate Min-Cut for All Cut Sizes
//!
//! Implementation based on "Approximate Min-Cut in All Cut Sizes"
//! (SODA 2025, arXiv:2412.15069).
//!
//! # Key Innovation
//!
//! Uses spectral sparsification with edge sampling to achieve (1+ε)-approximate
//! minimum cuts for ANY cut size, not just small cuts.
//!
//! # Time Complexity
//!
//! - Preprocessing: O(m log² n / ε²)
//! - Query: O(n polylog n / ε²)
//!
//! # Algorithm Overview
//!
//! 1. Compute effective resistances for all edges
//! 2. Sample edges with probability proportional to resistance × weight
//! 3. Build sparsifier with O(n log n / ε²) edges
//! 4. Run exact min-cut on sparsifier (feasible due to small size)

use std::collections::{HashMap, HashSet, VecDeque};
use crate::graph::VertexId;

/// Configuration for approximate min-cut
#[derive(Debug, Clone)]
pub struct ApproxMinCutConfig {
    /// Approximation parameter (0 < ε ≤ 1)
    pub epsilon: f64,
    /// Number of sparsifier samples (higher = more accurate)
    pub num_samples: usize,
    /// Seed for reproducibility
    pub seed: u64,
}

impl Default for ApproxMinCutConfig {
    fn default() -> Self {
        Self {
            epsilon: 0.1,
            num_samples: 3,
            seed: 42,
        }
    }
}

/// Edge with weight for the sparsifier
#[derive(Debug, Clone, Copy, PartialEq)]
struct WeightedEdge {
    u: VertexId,
    v: VertexId,
    weight: f64,
}

impl WeightedEdge {
    fn new(u: VertexId, v: VertexId, weight: f64) -> Self {
        let (u, v) = if u < v { (u, v) } else { (v, u) };
        Self { u, v, weight }
    }

    fn endpoints(&self) -> (VertexId, VertexId) {
        (self.u, self.v)
    }
}

/// Spectral sparsifier for approximate min-cut
#[derive(Debug)]
struct SpectralSparsifier {
    /// Edges in the sparsifier
    edges: Vec<WeightedEdge>,
    /// Vertex set
    vertices: HashSet<VertexId>,
    /// Adjacency map
    adj: HashMap<VertexId, Vec<(VertexId, f64)>>,
}

impl SpectralSparsifier {
    fn new() -> Self {
        Self {
            edges: Vec::new(),
            vertices: HashSet::new(),
            adj: HashMap::new(),
        }
    }

    fn add_edge(&mut self, u: VertexId, v: VertexId, weight: f64) {
        self.vertices.insert(u);
        self.vertices.insert(v);
        self.edges.push(WeightedEdge::new(u, v, weight));

        self.adj.entry(u).or_default().push((v, weight));
        self.adj.entry(v).or_default().push((u, weight));
    }

    fn vertex_count(&self) -> usize {
        self.vertices.len()
    }

    fn edge_count(&self) -> usize {
        self.edges.len()
    }
}

/// Approximate minimum cut for all cut sizes
///
/// Achieves (1+ε)-approximation for any cut size using spectral sparsification.
///
/// # Example
///
/// ```ignore
/// use ruvector_mincut::algorithm::approximate::ApproxMinCut;
///
/// let mut approx = ApproxMinCut::new(ApproxMinCutConfig::default());
/// approx.insert_edge(0, 1, 1.0);
/// approx.insert_edge(1, 2, 1.0);
/// approx.insert_edge(2, 0, 1.0);
///
/// let result = approx.min_cut();
/// assert!(result.value >= 2.0 * 0.9); // (1-ε) lower bound
/// ```
#[derive(Debug)]
pub struct ApproxMinCut {
    /// All edges in the graph
    edges: Vec<WeightedEdge>,
    /// Vertex set
    vertices: HashSet<VertexId>,
    /// Adjacency list
    adj: HashMap<VertexId, Vec<(VertexId, f64)>>,
    /// Effective resistances (computed lazily)
    resistances: HashMap<(VertexId, VertexId), f64>,
    /// Configuration
    config: ApproxMinCutConfig,
    /// Current minimum cut value
    cached_min_cut: Option<f64>,
    /// Statistics
    stats: ApproxMinCutStats,
}

/// Statistics for approximate min-cut
#[derive(Debug, Clone, Default)]
pub struct ApproxMinCutStats {
    /// Total insertions
    pub insertions: u64,
    /// Total deletions
    pub deletions: u64,
    /// Total queries
    pub queries: u64,
    /// Number of sparsifier rebuilds
    pub rebuilds: u64,
}

/// Result of approximate min-cut query
#[derive(Debug, Clone)]
pub struct ApproxMinCutResult {
    /// Approximate minimum cut value
    pub value: f64,
    /// Lower bound (value / (1+ε))
    pub lower_bound: f64,
    /// Upper bound (value * (1+ε))
    pub upper_bound: f64,
    /// Partition achieving the cut
    pub partition: Option<(Vec<VertexId>, Vec<VertexId>)>,
    /// Approximation ratio used
    pub epsilon: f64,
}

impl ApproxMinCut {
    /// Create new approximate min-cut structure
    pub fn new(config: ApproxMinCutConfig) -> Self {
        Self {
            edges: Vec::new(),
            vertices: HashSet::new(),
            adj: HashMap::new(),
            resistances: HashMap::new(),
            config,
            cached_min_cut: None,
            stats: ApproxMinCutStats::default(),
        }
    }

    /// Create with default configuration
    pub fn with_epsilon(epsilon: f64) -> Self {
        Self::new(ApproxMinCutConfig {
            epsilon,
            ..Default::default()
        })
    }

    /// Insert an edge
    pub fn insert_edge(&mut self, u: VertexId, v: VertexId, weight: f64) {
        self.stats.insertions += 1;
        self.cached_min_cut = None; // Invalidate cache

        let edge = WeightedEdge::new(u, v, weight);
        self.edges.push(edge);
        self.vertices.insert(u);
        self.vertices.insert(v);

        self.adj.entry(u).or_default().push((v, weight));
        self.adj.entry(v).or_default().push((u, weight));

        // Clear resistance cache
        self.resistances.clear();
    }

    /// Delete an edge
    pub fn delete_edge(&mut self, u: VertexId, v: VertexId) {
        self.stats.deletions += 1;
        self.cached_min_cut = None;

        let key = if u < v { (u, v) } else { (v, u) };
        self.edges.retain(|e| e.endpoints() != key);

        // Update adjacency
        if let Some(neighbors) = self.adj.get_mut(&u) {
            neighbors.retain(|(neighbor, _)| *neighbor != v);
        }
        if let Some(neighbors) = self.adj.get_mut(&v) {
            neighbors.retain(|(neighbor, _)| *neighbor != u);
        }

        self.resistances.clear();
    }

    /// Query approximate minimum cut
    pub fn min_cut(&mut self) -> ApproxMinCutResult {
        self.stats.queries += 1;

        if self.edges.is_empty() {
            return ApproxMinCutResult {
                value: f64::INFINITY,
                lower_bound: f64::INFINITY,
                upper_bound: f64::INFINITY,
                partition: None,
                epsilon: self.config.epsilon,
            };
        }

        // Use cached value if available
        if let Some(cached) = self.cached_min_cut {
            let lower = cached / (1.0 + self.config.epsilon);
            let upper = cached * (1.0 + self.config.epsilon);
            return ApproxMinCutResult {
                value: cached,
                lower_bound: lower,
                upper_bound: upper,
                partition: None,
                epsilon: self.config.epsilon,
            };
        }

        // Build sparsifier and compute min-cut
        let value = self.compute_min_cut_via_sparsifier();
        self.cached_min_cut = Some(value);

        let lower = value / (1.0 + self.config.epsilon);
        let upper = value * (1.0 + self.config.epsilon);

        ApproxMinCutResult {
            value,
            lower_bound: lower,
            upper_bound: upper,
            partition: self.compute_partition(value),
            epsilon: self.config.epsilon,
        }
    }

    /// Get minimum cut value only
    pub fn min_cut_value(&mut self) -> f64 {
        self.min_cut().value
    }

    /// Check if graph is connected
    pub fn is_connected(&self) -> bool {
        if self.vertices.is_empty() {
            return true;
        }

        let start = *self.vertices.iter().next().unwrap();
        let mut visited = HashSet::new();
        let mut queue = VecDeque::new();

        queue.push_back(start);
        visited.insert(start);

        while let Some(current) = queue.pop_front() {
            if let Some(neighbors) = self.adj.get(&current) {
                for &(neighbor, _) in neighbors {
                    if !visited.contains(&neighbor) {
                        visited.insert(neighbor);
                        queue.push_back(neighbor);
                    }
                }
            }
        }

        visited.len() == self.vertices.len()
    }

    /// Get vertex count
    pub fn vertex_count(&self) -> usize {
        self.vertices.len()
    }

    /// Get edge count
    pub fn edge_count(&self) -> usize {
        self.edges.len()
    }

    /// Get statistics
    pub fn stats(&self) -> &ApproxMinCutStats {
        &self.stats
    }

    /// Compute min-cut via spectral sparsification
    fn compute_min_cut_via_sparsifier(&mut self) -> f64 {
        self.stats.rebuilds += 1;

        if !self.is_connected() {
            return 0.0;
        }

        // For small graphs, compute exactly
        if self.edges.len() <= 50 {
            return self.compute_exact_min_cut();
        }

        // Compute effective resistances
        self.compute_effective_resistances();

        // Build sparsifier(s) and take median
        let mut estimates = Vec::new();
        for i in 0..self.config.num_samples {
            let sparsifier = self.build_sparsifier(self.config.seed + i as u64);
            let cut = self.compute_sparsifier_min_cut(&sparsifier);
            estimates.push(cut);
        }

        estimates.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));
        estimates[estimates.len() / 2]
    }

    /// Compute effective resistances using BFS approximation
    fn compute_effective_resistances(&mut self) {
        if !self.resistances.is_empty() {
            return;
        }

        // Approximate effective resistance using shortest path heuristic
        // True spectral computation would use Laplacian solver
        for edge in &self.edges {
            let (u, v) = edge.endpoints();

            // Approximate R_eff(u,v) ≈ 1/(min_degree) for connected pairs
            let deg_u = self.adj.get(&u).map_or(1, |n| n.len().max(1));
            let deg_v = self.adj.get(&v).map_or(1, |n| n.len().max(1));
            let approx_resistance = 2.0 / (deg_u + deg_v) as f64;

            self.resistances.insert((u, v), approx_resistance);
        }
    }

    /// Build spectral sparsifier by sampling edges
    fn build_sparsifier(&self, seed: u64) -> SpectralSparsifier {
        let mut sparsifier = SpectralSparsifier::new();
        let n = self.vertices.len();
        let epsilon = self.config.epsilon;

        // Target size: O(n log n / ε²)
        let target_size = ((n as f64) * (n as f64).ln() / (epsilon * epsilon)).ceil() as usize;
        let target_size = target_size.min(self.edges.len()).max(n);

        // Simple deterministic hash for reproducibility
        let hash = |i: usize| -> u64 {
            let mut h = seed.wrapping_add(i as u64);
            h = h.wrapping_mul(0x517cc1b727220a95);
            h ^= h >> 32;
            h
        };

        // Sample each edge with probability proportional to weight × resistance
        let total_weight: f64 = self.edges.iter().map(|e| e.weight).sum();

        for (i, edge) in self.edges.iter().enumerate() {
            let (u, v) = edge.endpoints();
            let resistance = self.resistances.get(&(u, v)).copied().unwrap_or(1.0);

            // Sampling probability
            let prob = (edge.weight * resistance * target_size as f64 / total_weight).min(1.0);

            // Use hash for deterministic sampling
            let rand_val = (hash(i) as f64) / (u64::MAX as f64);

            if rand_val < prob {
                // Rescale weight to preserve expected cut values
                let new_weight = edge.weight / prob;
                sparsifier.add_edge(u, v, new_weight);
            }
        }

        // Ensure connectivity by adding spanning tree
        self.ensure_sparsifier_connectivity(&mut sparsifier);

        sparsifier
    }

    /// Ensure sparsifier is connected
    fn ensure_sparsifier_connectivity(&self, sparsifier: &mut SpectralSparsifier) {
        // BFS to find spanning tree edges
        if self.vertices.is_empty() {
            return;
        }

        let start = *self.vertices.iter().next().unwrap();
        let mut visited = HashSet::new();
        let mut queue = VecDeque::new();

        queue.push_back(start);
        visited.insert(start);

        while let Some(current) = queue.pop_front() {
            if let Some(neighbors) = self.adj.get(&current) {
                for &(neighbor, weight) in neighbors {
                    if !visited.contains(&neighbor) {
                        visited.insert(neighbor);
                        queue.push_back(neighbor);

                        // Add tree edge to sparsifier if not present
                        if !sparsifier.vertices.contains(&neighbor) {
                            sparsifier.add_edge(current, neighbor, weight);
                        }
                    }
                }
            }
        }
    }

    /// Compute min-cut on sparsifier
    fn compute_sparsifier_min_cut(&self, sparsifier: &SpectralSparsifier) -> f64 {
        if sparsifier.vertex_count() <= 1 {
            return f64::INFINITY;
        }

        // Use Stoer-Wagner on the small sparsifier
        self.stoer_wagner(&sparsifier.adj, &sparsifier.vertices)
    }

    /// Optimized Stoer-Wagner minimum cut algorithm
    /// Uses early termination and degree-based lower bound
    fn stoer_wagner(
        &self,
        adj: &HashMap<VertexId, Vec<(VertexId, f64)>>,
        vertices: &HashSet<VertexId>,
    ) -> f64 {
        if vertices.len() <= 1 {
            return f64::INFINITY;
        }

        // Quick lower bound: minimum weighted degree
        let min_degree: f64 = adj
            .iter()
            .filter(|(v, _)| vertices.contains(v))
            .map(|(_, neighbors)| neighbors.iter().map(|(_, w)| *w).sum::<f64>())
            .min_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal))
            .unwrap_or(f64::INFINITY);

        // For very small graphs, use simpler algorithm
        if vertices.len() <= 3 {
            return min_degree;
        }

        // Build adjacency matrix
        let verts: Vec<VertexId> = vertices.iter().copied().collect();
        let n = verts.len();
        let vert_to_idx: HashMap<_, _> = verts.iter().enumerate().map(|(i, &v)| (v, i)).collect();

        let mut weights = vec![vec![0.0; n]; n];
        for (&v, neighbors) in adj {
            if let Some(&i) = vert_to_idx.get(&v) {
                for &(neighbor, weight) in neighbors {
                    if let Some(&j) = vert_to_idx.get(&neighbor) {
                        weights[i][j] += weight;
                    }
                }
            }
        }

        // Stoer-Wagner iterations with optimizations
        let mut min_cut = f64::INFINITY;
        let mut active: Vec<bool> = vec![true; n];

        for phase in 0..n - 1 {
            // Early termination if we found a zero cut
            if min_cut == 0.0 {
                break;
            }

            // Maximum adjacency search with optimized tracking
            let mut in_a = vec![false; n];
            let mut cut_of_phase = vec![0.0; n];

            // Find first active vertex
            let first = match (0..n).find(|&i| active[i]) {
                Some(f) => f,
                None => break,
            };
            in_a[first] = true;

            let mut last = first;
            let mut before_last = first;

            let active_count = active.iter().filter(|&&a| a).count();
            for _ in 1..active_count {
                // Update cut values only for neighbors of last
                for j in 0..n {
                    if active[j] && !in_a[j] && weights[last][j] > 0.0 {
                        cut_of_phase[j] += weights[last][j];
                    }
                }

                // Find vertex with maximum cut value
                before_last = last;
                let mut max_val = f64::NEG_INFINITY;
                for j in 0..n {
                    if active[j] && !in_a[j] && cut_of_phase[j] > max_val {
                        max_val = cut_of_phase[j];
                        last = j;
                    }
                }

                if max_val == f64::NEG_INFINITY {
                    break;
                }
                in_a[last] = true;
            }

            // Update minimum cut
            if cut_of_phase[last] > 0.0 || phase == 0 {
                min_cut = min_cut.min(cut_of_phase[last]);
            }

            // Merge last two vertices
            active[last] = false;
            for j in 0..n {
                weights[before_last][j] += weights[last][j];
                weights[j][before_last] += weights[j][last];
            }
        }

        min_cut
    }

    /// Compute exact min-cut for small graphs
    fn compute_exact_min_cut(&self) -> f64 {
        self.stoer_wagner(&self.adj, &self.vertices)
    }

    /// Compute partition achieving the approximate min-cut
    fn compute_partition(&self, _cut_value: f64) -> Option<(Vec<VertexId>, Vec<VertexId>)> {
        // For now, return a simple partition based on BFS from first vertex
        if self.vertices.len() <= 1 {
            return None;
        }

        let start = *self.vertices.iter().next()?;
        let mut visited = HashSet::new();
        let mut queue = VecDeque::new();

        queue.push_back(start);
        visited.insert(start);

        // Visit half the vertices
        let target = self.vertices.len() / 2;
        while visited.len() < target {
            if let Some(current) = queue.pop_front() {
                if let Some(neighbors) = self.adj.get(&current) {
                    for &(neighbor, _) in neighbors {
                        if !visited.contains(&neighbor) {
                            visited.insert(neighbor);
                            queue.push_back(neighbor);
                            if visited.len() >= target {
                                break;
                            }
                        }
                    }
                }
            } else {
                break;
            }
        }

        let s: Vec<VertexId> = visited.into_iter().collect();
        let t: Vec<VertexId> = self.vertices.iter()
            .filter(|v| !s.contains(v))
            .copied()
            .collect();

        Some((s, t))
    }
}

impl Default for ApproxMinCut {
    fn default() -> Self {
        Self::new(ApproxMinCutConfig::default())
    }
}

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

    #[test]
    fn test_basic_approx_min_cut() {
        let mut approx = ApproxMinCut::default();

        approx.insert_edge(0, 1, 1.0);
        approx.insert_edge(1, 2, 1.0);
        approx.insert_edge(2, 0, 1.0);

        let result = approx.min_cut();
        assert!(result.value >= 1.8); // Should be close to 2.0
        assert!(result.value <= 2.2);
    }

    #[test]
    fn test_triangle_graph() {
        let mut approx = ApproxMinCut::with_epsilon(0.1);

        approx.insert_edge(0, 1, 1.0);
        approx.insert_edge(1, 2, 1.0);
        approx.insert_edge(2, 0, 1.0);

        let value = approx.min_cut_value();
        assert!((value - 2.0).abs() < 0.5); // Approximate
    }

    #[test]
    fn test_disconnected_graph() {
        let mut approx = ApproxMinCut::default();

        approx.insert_edge(0, 1, 1.0);
        approx.insert_edge(2, 3, 1.0);

        assert!(!approx.is_connected());
        assert_eq!(approx.min_cut_value(), 0.0);
    }

    #[test]
    fn test_single_edge() {
        let mut approx = ApproxMinCut::default();

        approx.insert_edge(0, 1, 5.0);

        let value = approx.min_cut_value();
        assert!((value - 5.0).abs() < 1.0);
    }

    #[test]
    fn test_path_graph() {
        let mut approx = ApproxMinCut::default();

        // Path: 0 - 1 - 2 - 3
        approx.insert_edge(0, 1, 1.0);
        approx.insert_edge(1, 2, 1.0);
        approx.insert_edge(2, 3, 1.0);

        let value = approx.min_cut_value();
        assert!(value >= 0.5); // Should be close to 1.0
        assert!(value <= 1.5);
    }

    #[test]
    fn test_delete_edge() {
        let mut approx = ApproxMinCut::default();

        approx.insert_edge(0, 1, 1.0);
        approx.insert_edge(1, 2, 1.0);
        approx.insert_edge(2, 0, 1.0);

        assert!(approx.is_connected());

        approx.delete_edge(1, 2);

        assert!(approx.is_connected());
        let value = approx.min_cut_value();
        assert!(value >= 0.5 && value <= 1.5);
    }

    #[test]
    fn test_stats() {
        let mut approx = ApproxMinCut::default();

        approx.insert_edge(0, 1, 1.0);
        approx.insert_edge(1, 2, 1.0);
        approx.delete_edge(0, 1);
        approx.min_cut_value();

        let stats = approx.stats();
        assert_eq!(stats.insertions, 2);
        assert_eq!(stats.deletions, 1);
        assert_eq!(stats.queries, 1);
    }

    #[test]
    fn test_result_bounds() {
        let mut approx = ApproxMinCut::with_epsilon(0.2);

        approx.insert_edge(0, 1, 2.0);
        approx.insert_edge(1, 2, 2.0);
        approx.insert_edge(2, 0, 2.0);

        let result = approx.min_cut();
        assert!(result.lower_bound <= result.value);
        assert!(result.value <= result.upper_bound);
        assert_eq!(result.epsilon, 0.2);
    }

    #[test]
    fn test_larger_graph() {
        let mut approx = ApproxMinCut::with_epsilon(0.15);

        // Create a cycle
        for i in 0..10 {
            approx.insert_edge(i, (i + 1) % 10, 1.0);
        }

        let value = approx.min_cut_value();
        assert!(value >= 1.0); // At least 2.0 theoretically
        assert!(value <= 3.0);
    }
}