tide-maxflow 0.1.0

//! # Tide Flow
//!
//! Implementation of the **Tide** max flow algorithm — a push-pull-relabel variant.
//!
//! The algorithm performs global sweeps (BFS relabeling, then global pull, then
//! global push) in each iteration ("tide"). It uses O(1) array-based data
//! structures with `i64` capacities for maximum performance.
//!
//! ## Algorithm Overview
//!
//! Each tide consists of three phases:
//! 1. **globalRelabel** — BFS from source and sink on the residual graph to
//!    compute exact vertex heights.
//! 2. **globalPull** — Right-fold over overflowing vertices (descending level),
//!    pulling flow on reverse edges.
//! 3. **globalPush** — Left-fold over overflowing vertices (ascending level),
//!    pushing flow on forward edges.
//!
//! Terminates when net flow and the overflowing set are unchanged between tides.
//!
//! ## Optimizations
//!
//! - Skips `globalRelabel` when no edge crossed a saturation boundary.
//! - Precomputed level buckets (levels are constant).
//! - Contiguous adjacency storage (CSR-like head/next arrays).
//! - Reused BFS buffers across tides.
//! - Allocation-free convergence check via excess hash.

/// Result of the Tide max flow computation.
#[derive(Debug, Clone)]
pub struct MaxFlowResult {
    /// Maximum flow value
    pub flow: i64,
    /// Number of tides (iterations)
    pub steps: usize,
    /// Flow on each edge (indexed by edge index; even = forward, odd = reverse)
    pub edge_flow: Vec<i64>,
}

/// Per-phase timing breakdown (in microseconds) for profiling.
#[derive(Debug, Clone, Default)]
pub struct PhaseTiming {
    /// Total time in globalRelabel (BFS), microseconds
    pub relabel_us: u64,
    /// Total time in globalPull, microseconds
    pub pull_us: u64,
    /// Total time in globalPush, microseconds
    pub push_us: u64,
    /// Total time in convergence checks, microseconds
    pub convergence_us: u64,
    /// Number of tides where globalRelabel was skipped
    pub relabel_skipped: usize,
}

/// Result with timing breakdown.
#[derive(Debug, Clone)]
pub struct TimedMaxFlowResult {
    /// Standard max flow result
    pub result: MaxFlowResult,
    /// Per-phase timing breakdown
    pub timing: PhaseTiming,
}

/// Internal state for the Tide algorithm.
///
/// Uses flat arrays with paired edges: forward edge `2*i`, reverse edge `2*i+1`.
/// Adjacency stored as CSR-like head[v]/adj_edges[] for cache locality.
/// Within each vertex's adjacency list, forward edges come first, then reverse.
struct TideState {
    n: usize,
    source: usize,
    sink: usize,

    // --- CSR-like adjacency ---
    /// head[v]..head[v+1] gives the range into adj_edges for vertex v
    head: Vec<usize>,
    /// Flat array of edge indices, grouped by source vertex.
    /// Forward edges (even) are stored before reverse edges (odd) within each group.
    adj_edges: Vec<usize>,
    /// fwd_end[v] is the end of forward edges for vertex v within adj_edges.
    /// Forward edges: head[v]..fwd_end[v], reverse: fwd_end[v]..head[v+1].
    fwd_end: Vec<usize>,

    // --- Edge data (flat, paired) ---
    /// Target vertex of each edge
    to: Vec<usize>,
    /// Capacity of each edge
    cap: Vec<i64>,
    /// Current flow on each edge (residual = cap[e] - flow[e])
    flow: Vec<i64>,

    // --- Vertex data ---
    /// Height (from globalRelabel BFS)
    height: Vec<i32>,
    /// Excess flow
    excess: Vec<i64>,
    /// Level (BFS distance from source in original graph, constant)
    #[allow(dead_code)]
    level: Vec<i32>,

    // --- Precomputed level buckets ---
    /// max_level value
    max_level: usize,
    /// level_buckets[l] = list of vertices at level l (excluding source/sink)
    level_buckets: Vec<Vec<usize>>,

    // --- Reusable BFS buffers ---
    /// Two flat vectors used as alternating current/next-level buffers
    bfs_cur: Vec<usize>,
    bfs_next: Vec<usize>,
    /// Distance buffers (reset to -1 each BFS)
    dist_buf1: Vec<i32>,
    dist_buf2: Vec<i32>,

    // --- State ---
    topology_changed: bool,
    steps: usize,
}

impl TideState {
    fn new(n: usize, source: usize, sink: usize, edges: &[(usize, usize, i64)]) -> Self {
        let m = edges.len();
        let total = 2 * m;

        let mut to = vec![0usize; total];
        let mut cap = vec![0i64; total];

        for (i, &(u, v, c)) in edges.iter().enumerate() {
            let fwd = 2 * i;
            let rev = 2 * i + 1;
            to[fwd] = v;
            to[rev] = u;
            cap[fwd] = c;
            cap[rev] = 0;
        }

        // Two-pass CSR construction with forward edges first per vertex.
        // Pass 1: count forward and reverse degree per vertex.
        let mut fwd_deg = vec![0u32; n];
        let mut rev_deg = vec![0u32; n];
        for (i, &(u, v, _)) in edges.iter().enumerate() {
            let _ = i;
            fwd_deg[u] += 1;
            rev_deg[v] += 1;
        }

        // Compute head offsets and fwd_end boundary
        let mut head = vec![0usize; n + 1];
        let mut fwd_end = vec![0usize; n];
        for v in 0..n {
            head[v + 1] = head[v] + (fwd_deg[v] + rev_deg[v]) as usize;
            fwd_end[v] = head[v] + fwd_deg[v] as usize;
        }
        let total_adj = head[n];
        let mut adj_edges = vec![0usize; total_adj];

        // Pass 2: fill CSR. Use fwd_deg/rev_deg as running insertion pointers.
        // Reset to use as write cursors: fwd starts at head[v], rev at fwd_end[v].
        let mut fwd_pos: Vec<usize> = (0..n).map(|v| head[v]).collect();
        let mut rev_pos: Vec<usize> = fwd_end.clone();

        for (i, &(u, v, _)) in edges.iter().enumerate() {
            let fwd = 2 * i;
            let rev = 2 * i + 1;
            adj_edges[fwd_pos[u]] = fwd;
            fwd_pos[u] += 1;
            adj_edges[rev_pos[v]] = rev;
            rev_pos[v] += 1;
        }
        drop(fwd_deg);
        drop(rev_deg);
        drop(fwd_pos);
        drop(rev_pos);

        let flow = vec![0i64; total];
        let height = vec![0i32; n];
        let excess = vec![0i64; n];

        // Compute levels: BFS from source on original graph using flat buffer
        let mut level = vec![-1i32; n];
        level[source] = 0;
        let mut bfs_cur = vec![source];
        let mut bfs_next = Vec::new();
        while !bfs_cur.is_empty() {
            for &v in &bfs_cur {
                // Only forward edges for level BFS
                for idx in head[v]..fwd_end[v] {
                    let e = adj_edges[idx];
                    let w = to[e];
                    if level[w] == -1 {
                        level[w] = level[v] + 1;
                        bfs_next.push(w);
                    }
                }
            }
            std::mem::swap(&mut bfs_cur, &mut bfs_next);
            bfs_next.clear();
        }

        // Precompute level buckets
        let max_level = level.iter().copied().filter(|&l| l >= 0).max().unwrap_or(0) as usize;
        let mut level_buckets: Vec<Vec<usize>> = vec![Vec::new(); max_level + 1];
        for v in 0..n {
            if v != source && v != sink && level[v] >= 0 {
                level_buckets[level[v] as usize].push(v);
            }
        }

        let dist_buf1 = vec![-1i32; n];
        let dist_buf2 = vec![-1i32; n];

        TideState {
            n,
            source,
            sink,
            head,
            adj_edges,
            fwd_end,
            to,
            cap,
            flow,
            height,
            excess,
            level,
            max_level,
            level_buckets,
            bfs_cur,
            bfs_next,
            dist_buf1,
            dist_buf2,
            topology_changed: true,
            steps: 0,
        }
    }

    /// Saturate all source edges to initialize the preflow.
    fn initialize(&mut self) {
        let source = self.source;
        for idx in self.head[source]..self.fwd_end[source] {
            let e = self.adj_edges[idx]; // guaranteed forward (even)
            let c = self.cap[e];
            if c > 0 {
                self.flow[e] = c;
                self.flow[e ^ 1] = -c;
                let v = self.to[e];
                self.excess[v] += c;
                self.excess[source] -= c;
            }
        }
        self.height[source] = self.n as i32;
    }

    /// Global relabel: BFS from sink (and source if needed) on residual graph.
    /// Uses flat level-synchronous BFS with separate forward/reverse CSR lists.
    fn global_relabel(&mut self) {
        let n = self.n;
        let source = self.source;
        let sink = self.sink;

        // Reset sink dist buffer
        let dist_sink = &mut self.dist_buf1;
        for x in dist_sink.iter_mut() {
            *x = -1;
        }
        dist_sink[sink] = 0;

        // Sink BFS using reverse-only CSR: at vertex v, iterate incoming forward
        // edges (stored as reverse edges in rev_edges). The reverse edge e (odd)
        // has e^1 = forward edge. To traverse backward on residual, we need
        // residual of the reverse edge from w→v, i.e. cap[e^1] - flow[e^1] > 0.
        // Wait — for BFS from sink, we go against the flow direction.
        // At vertex v in sink BFS, we want to reach neighbors w such that there
        // is residual capacity from w to v (so flow can reach sink via w→v).
        // That means: residual(w→v) > 0, i.e. cap[e_wv] - flow[e_wv] > 0.
        // Using the all-edges CSR at v: for each edge e from v, the reverse e^1
        // goes from to[e] back to v. Residual on e^1 = cap[e^1] - flow[e^1].
        // If that > 0, then to[e] can reach v on the residual graph.
        //
        // With separate CSR: at v, fwd_edges gives forward edges from v (we need
        // their reverses to check if someone can reach v), and rev_edges gives
        // reverse edges at v (original edges pointing TO v — their residual
        // gives backward reachability).
        //
        // Sink BFS wants: from v, find all w such that res(w→v) > 0.
        // Edge from w to v is a forward edge e (even), stored in fwd_edges of w
        // and as reverse edge e^1 (odd) in rev_edges of v.
        // res(w→v) via forward edge = cap[e] - flow[e] (the original forward edge).
        // The reverse edge at v is e^1 = e|1, and the forward is e^1^1 = e.
        // So: for rev edge r in rev_edges[v], fwd = r^1, w = to[r],
        //     res = cap[fwd] - flow[fwd]. But wait, that's residual on the
        //     original forward edge (w→v direction). For BFS from sink going
        //     backward, we want edges where flow can still be sent toward sink,
        //     meaning residual from w to v on the residual graph.
        //
        // Actually, let's think simply. In the residual graph, edge (w,v) exists
        // if cap[e_wv] - flow[e_wv] > 0 for any edge e_wv (forward or reverse).
        // BFS from sink traverses the REVERSE of residual edges.
        // So from v in BFS, we want w such that (w,v) is a residual edge.
        //
        // Using all-edges CSR at v: edge e from v→to[e]. The reverse e^1 is
        // from to[e]→v. If cap[e^1]-flow[e^1] > 0, then to[e] has residual to v.
        // This is the original code. Let's just use all-edges CSR for BFS
        // (it iterates all edges anyway) but with the flat buffer.
        self.bfs_cur.clear();
        self.bfs_cur.push(sink);
        let mut reached = 1usize;

        while !self.bfs_cur.is_empty() {
            self.bfs_next.clear();
            for i in 0..self.bfs_cur.len() {
                let v = self.bfs_cur[i];
                let d_next = dist_sink[v] + 1;
                for idx in self.head[v]..self.head[v + 1] {
                    let e = self.adj_edges[idx];
                    let w = self.to[e];
                    let rev = e ^ 1;
                    if dist_sink[w] == -1 && self.cap[rev] - self.flow[rev] > 0 {
                        dist_sink[w] = d_next;
                        self.bfs_next.push(w);
                        reached += 1;
                    }
                }
            }
            std::mem::swap(&mut self.bfs_cur, &mut self.bfs_next);
        }

        // Source BFS: only needed if sink BFS didn't reach all vertices
        let need_source_bfs = reached < n;
        if need_source_bfs {
            let dist_source = &mut self.dist_buf2;
            for x in dist_source.iter_mut() {
                *x = -1;
            }
            dist_source[source] = 0;

            self.bfs_cur.clear();
            self.bfs_cur.push(source);

            while !self.bfs_cur.is_empty() {
                self.bfs_next.clear();
                for i in 0..self.bfs_cur.len() {
                    let v = self.bfs_cur[i];
                    let d_next = dist_source[v] + 1;
                    for idx in self.head[v]..self.head[v + 1] {
                        let e = self.adj_edges[idx];
                        let w = self.to[e];
                        if dist_source[w] == -1 && self.cap[e] - self.flow[e] > 0 {
                            dist_source[w] = d_next;
                            self.bfs_next.push(w);
                        }
                    }
                }
                std::mem::swap(&mut self.bfs_cur, &mut self.bfs_next);
            }
        }

        // Set heights
        let n_i32 = n as i32;
        for v in 0..n {
            if v == source {
                self.height[v] = n_i32;
            } else if v == sink {
                self.height[v] = 0;
            } else if self.dist_buf1[v] >= 0 {
                self.height[v] = self.dist_buf1[v];
            } else if need_source_bfs && self.dist_buf2[v] >= 0 {
                self.height[v] = n_i32 + self.dist_buf2[v];
            }
        }
    }

    /// Global push: ascending level order, forward edges only.
    /// Uses fwd_end boundary to iterate only forward edges.
    fn global_push(&mut self) {
        for lvl in 0..=self.max_level {
            for i in 0..self.level_buckets[lvl].len() {
                let v = self.level_buckets[lvl][i];
                if self.excess[v] <= 0 {
                    continue;
                }
                let h_target = self.height[v] - 1;
                for idx in self.head[v]..self.fwd_end[v] {
                    let e = self.adj_edges[idx]; // guaranteed even (forward)
                    let w = self.to[e];
                    let res_cap = self.cap[e] - self.flow[e];
                    if self.height[w] == h_target && res_cap > 0 {
                        let delta = self.excess[v].min(res_cap);

                        let old_bwd = self.cap[e ^ 1] - self.flow[e ^ 1] > 0;

                        self.flow[e] += delta;
                        self.flow[e ^ 1] -= delta;
                        self.excess[v] -= delta;
                        self.excess[w] += delta;

                        let new_fwd = self.cap[e] - self.flow[e] > 0;
                        let new_bwd = self.cap[e ^ 1] - self.flow[e ^ 1] > 0;
                        if !new_fwd || (old_bwd != new_bwd) {
                            self.topology_changed = true;
                        }

                        if self.excess[v] <= 0 {
                            break;
                        }
                    }
                }
            }
        }
    }

    /// Global pull: descending level order, reverse edges only.
    /// Uses fwd_end boundary to iterate only reverse edges.
    fn global_pull(&mut self) {
        for lvl in (0..=self.max_level).rev() {
            for i in 0..self.level_buckets[lvl].len() {
                let v = self.level_buckets[lvl][i];
                if self.excess[v] <= 0 {
                    continue;
                }
                let h_target = self.height[v] - 1;
                for idx in self.fwd_end[v]..self.head[v + 1] {
                    let e = self.adj_edges[idx]; // guaranteed odd (reverse)
                    let fwd = e ^ 1; // the forward edge
                    let u = self.to[e]; // = to[rev] = source of forward edge
                    let res_cap = self.cap[e] - self.flow[e]; // residual on reverse edge

                    if self.height[u] == h_target && res_cap > 0 {
                        let delta = self.excess[v].min(res_cap);

                        let old_fwd = self.cap[fwd] - self.flow[fwd] > 0;

                        self.flow[e] += delta;
                        self.flow[fwd] -= delta;
                        self.excess[v] -= delta;
                        self.excess[u] += delta;

                        let new_fwd = self.cap[fwd] - self.flow[fwd] > 0;
                        let new_bwd = self.cap[e] - self.flow[e] > 0;
                        if (old_fwd != new_fwd) || !new_bwd {
                            self.topology_changed = true;
                        }

                        if self.excess[v] <= 0 {
                            break;
                        }
                    }
                }
            }
        }
    }

    /// Local relabel: scan all residual neighbors of v, return
    /// min(neighbor height) + 1, or None if no residual neighbors.
    #[inline]
    fn local_relabel(&self, v: usize) -> Option<i32> {
        let mut min_h = i32::MAX;

        // Forward residual: v -> w
        for idx in self.head[v]..self.fwd_end[v] {
            let e = self.adj_edges[idx];
            if self.cap[e] - self.flow[e] > 0 {
                min_h = min_h.min(self.height[self.to[e]]);
            }
        }
        // Reverse residual: u -> v (stored as reverse edge at v)
        for idx in self.fwd_end[v]..self.head[v + 1] {
            let e = self.adj_edges[idx];
            if self.cap[e] - self.flow[e] > 0 {
                min_h = min_h.min(self.height[self.to[e]]);
            }
        }

        if min_h == i32::MAX {
            None
        } else {
            Some(min_h + 1)
        }
    }

    /// Global pull with local relabeling.
    fn global_pull_relabel(&mut self) {
        for lvl in (0..=self.max_level).rev() {
            for i in 0..self.level_buckets[lvl].len() {
                let v = self.level_buckets[lvl][i];
                if self.excess[v] <= 0 {
                    continue;
                }

                let h_target = self.height[v] - 1;
                let mut pulled = false;
                for idx in self.fwd_end[v]..self.head[v + 1] {
                    let e = self.adj_edges[idx];
                    let fwd = e ^ 1;
                    let u = self.to[e];
                    let res_cap = self.cap[e] - self.flow[e];
                    if self.height[u] == h_target && res_cap > 0 {
                        let delta = self.excess[v].min(res_cap);
                        let old_fwd = self.cap[fwd] - self.flow[fwd] > 0;
                        self.flow[e] += delta;
                        self.flow[fwd] -= delta;
                        self.excess[v] -= delta;
                        self.excess[u] += delta;
                        let new_fwd = self.cap[fwd] - self.flow[fwd] > 0;
                        let new_bwd = self.cap[e] - self.flow[e] > 0;
                        if (old_fwd != new_fwd) || !new_bwd {
                            self.topology_changed = true;
                        }
                        pulled = true;
                        if self.excess[v] <= 0 {
                            break;
                        }
                    }
                }

                if !pulled && self.excess[v] > 0 {
                    if let Some(new_h) = self.local_relabel(v) {
                        if new_h > self.height[v] {
                            self.height[v] = new_h;
                            let h_target = new_h - 1;
                            for idx in self.fwd_end[v]..self.head[v + 1] {
                                let e = self.adj_edges[idx];
                                let fwd = e ^ 1;
                                let u = self.to[e];
                                let res_cap = self.cap[e] - self.flow[e];
                                if self.height[u] == h_target && res_cap > 0 {
                                    let delta = self.excess[v].min(res_cap);
                                    let old_fwd = self.cap[fwd] - self.flow[fwd] > 0;
                                    self.flow[e] += delta;
                                    self.flow[fwd] -= delta;
                                    self.excess[v] -= delta;
                                    self.excess[u] += delta;
                                    let new_fwd = self.cap[fwd] - self.flow[fwd] > 0;
                                    let new_bwd = self.cap[e] - self.flow[e] > 0;
                                    if (old_fwd != new_fwd) || !new_bwd {
                                        self.topology_changed = true;
                                    }
                                    if self.excess[v] <= 0 {
                                        break;
                                    }
                                }
                            }
                        }
                    }
                }
            }
        }
    }

    /// Global push with local relabeling.
    fn global_push_relabel(&mut self) {
        for lvl in 0..=self.max_level {
            for i in 0..self.level_buckets[lvl].len() {
                let v = self.level_buckets[lvl][i];
                if self.excess[v] <= 0 {
                    continue;
                }

                let h_target = self.height[v] - 1;
                let mut pushed = false;
                for idx in self.head[v]..self.fwd_end[v] {
                    let e = self.adj_edges[idx];
                    let w = self.to[e];
                    let res_cap = self.cap[e] - self.flow[e];
                    if self.height[w] == h_target && res_cap > 0 {
                        let delta = self.excess[v].min(res_cap);
                        let old_bwd = self.cap[e ^ 1] - self.flow[e ^ 1] > 0;
                        self.flow[e] += delta;
                        self.flow[e ^ 1] -= delta;
                        self.excess[v] -= delta;
                        self.excess[w] += delta;
                        let new_fwd = self.cap[e] - self.flow[e] > 0;
                        let new_bwd = self.cap[e ^ 1] - self.flow[e ^ 1] > 0;
                        if !new_fwd || (old_bwd != new_bwd) {
                            self.topology_changed = true;
                        }
                        pushed = true;
                        if self.excess[v] <= 0 {
                            break;
                        }
                    }
                }

                if !pushed && self.excess[v] > 0 {
                    if let Some(new_h) = self.local_relabel(v) {
                        if new_h > self.height[v] {
                            self.height[v] = new_h;
                            let h_target = new_h - 1;
                            for idx in self.head[v]..self.fwd_end[v] {
                                let e = self.adj_edges[idx];
                                let w = self.to[e];
                                let res_cap = self.cap[e] - self.flow[e];
                                if self.height[w] == h_target && res_cap > 0 {
                                    let delta = self.excess[v].min(res_cap);
                                    let old_bwd = self.cap[e ^ 1] - self.flow[e ^ 1] > 0;
                                    self.flow[e] += delta;
                                    self.flow[e ^ 1] -= delta;
                                    self.excess[v] -= delta;
                                    self.excess[w] += delta;
                                    let new_fwd = self.cap[e] - self.flow[e] > 0;
                                    let new_bwd = self.cap[e ^ 1] - self.flow[e ^ 1] > 0;
                                    if !new_fwd || (old_bwd != new_bwd) {
                                        self.topology_changed = true;
                                    }
                                    if self.excess[v] <= 0 {
                                        break;
                                    }
                                }
                            }
                        }
                    }
                }
            }
        }
    }

    /// Hybrid Tide: alternates between full BFS tides and local-relabel tides.
    ///
    /// Structure: outer loop does full BFS tide, then inner loop does
    /// local-relabel tides until flow stops increasing. Terminates when
    /// a full BFS tide makes no progress.
    fn tide_hybrid(&mut self) -> MaxFlowResult {
        self.initialize();

        // Step budget: if hybrid exceeds this, fall back to standard tides
        // for the remainder. On small graphs (n<=30) standard takes <=50 steps,
        // so 10*n^2 is a generous hybrid budget.
        let step_budget = 10 * self.n * self.n;
        let mut bfs_stall_count = 0u32;

        loop {
            // Full BFS tide
            let (old_flow, old_hash) = self.convergence_state();
            self.global_relabel();
            self.topology_changed = false;
            self.global_pull();
            self.global_push();
            self.steps += 1;

            let (bfs_flow, bfs_hash) = self.convergence_state();

            // Terminate if this BFS tide made no progress at all
            if bfs_flow == old_flow && bfs_hash == old_hash {
                break;
            }

            // If flow didn't increase, count as BFS stall; bail after 3
            if bfs_flow <= old_flow {
                bfs_stall_count += 1;
                if bfs_stall_count >= 3 {
                    // Fall back to standard tides for the rest
                    break;
                }
            } else {
                bfs_stall_count = 0;
            }

            // If we've used too many steps in hybrid mode, stop local-relabeling
            if self.steps > step_budget {
                continue; // skip local-relabel, just do BFS tides
            }

            // Local-relabel tides: continue while flow keeps increasing.
            // Bail as soon as flow stalls for 2 consecutive tides.
            let mut prev_flow = bfs_flow;
            let mut stall_count = 0u32;
            loop {
                let (_, old_hash2) = self.convergence_state();
                self.topology_changed = false;
                self.global_pull_relabel();
                self.global_push_relabel();
                self.steps += 1;

                let (new_flow2, new_hash2) = self.convergence_state();

                // Full convergence: no change at all
                if new_flow2 == prev_flow && new_hash2 == old_hash2 {
                    break;
                }

                if new_flow2 > prev_flow {
                    // Flow increased — reset stall counter
                    prev_flow = new_flow2;
                    stall_count = 0;
                } else {
                    // Flow flat but excess moved — allow a few more tries
                    stall_count += 1;
                    if stall_count >= 2 {
                        break;
                    }
                }
            }
        }

        // If we broke out of the hybrid loop due to BFS stall,
        // finish with standard tides to guarantee correctness.
        loop {
            let (old_flow, old_hash) = self.convergence_state();
            self.global_relabel();
            self.topology_changed = false;
            self.global_pull();
            self.global_push();
            self.steps += 1;
            let (new_flow, new_hash) = self.convergence_state();
            if new_flow == old_flow && new_hash == old_hash {
                break;
            }
        }

        let (flow, _) = self.convergence_state();
        MaxFlowResult {
            flow,
            steps: self.steps,
            edge_flow: self.flow.clone(),
        }
    }

    /// Hybrid Tide with per-phase timing.
    fn tide_hybrid_timed(&mut self) -> TimedMaxFlowResult {
        use std::time::Instant;
        self.initialize();
        let mut timing = PhaseTiming::default();
        let step_budget = 10 * self.n * self.n;
        let mut bfs_stall_count = 0u32;

        loop {
            // Full BFS tide
            let tc0 = Instant::now();
            let (old_flow, old_hash) = self.convergence_state();
            let tc1 = Instant::now();
            timing.convergence_us += tc1.duration_since(tc0).as_micros() as u64;

            let tb0 = Instant::now();
            self.global_relabel();
            let tb1 = Instant::now();
            timing.relabel_us += tb1.duration_since(tb0).as_micros() as u64;
            self.topology_changed = false;

            let tp0 = Instant::now();
            self.global_pull();
            let tp1 = Instant::now();
            timing.pull_us += tp1.duration_since(tp0).as_micros() as u64;

            let ts0 = Instant::now();
            self.global_push();
            let ts1 = Instant::now();
            timing.push_us += ts1.duration_since(ts0).as_micros() as u64;
            self.steps += 1;

            let tc2 = Instant::now();
            let (bfs_flow, bfs_hash) = self.convergence_state();
            let tc3 = Instant::now();
            timing.convergence_us += tc3.duration_since(tc2).as_micros() as u64;

            if bfs_flow == old_flow && bfs_hash == old_hash {
                break;
            }

            if bfs_flow <= old_flow {
                bfs_stall_count += 1;
                if bfs_stall_count >= 3 {
                    break;
                }
            } else {
                bfs_stall_count = 0;
            }

            if self.steps > step_budget {
                continue;
            }

            // Local-relabel tides
            let mut prev_flow = bfs_flow;
            let mut stall_count = 0u32;
            loop {
                let tc4 = Instant::now();
                let (_, old_hash2) = self.convergence_state();
                let tc5 = Instant::now();
                timing.convergence_us += tc5.duration_since(tc4).as_micros() as u64;
                timing.relabel_skipped += 1;
                self.topology_changed = false;

                let t2 = Instant::now();
                self.global_pull_relabel();
                let t3 = Instant::now();
                timing.pull_us += t3.duration_since(t2).as_micros() as u64;

                let t4 = Instant::now();
                self.global_push_relabel();
                let t5 = Instant::now();
                timing.push_us += t5.duration_since(t4).as_micros() as u64;
                self.steps += 1;

                let tc6 = Instant::now();
                let (new_flow2, new_hash2) = self.convergence_state();
                let tc7 = Instant::now();
                timing.convergence_us += tc7.duration_since(tc6).as_micros() as u64;

                if new_flow2 == prev_flow && new_hash2 == old_hash2 {
                    break;
                }

                if new_flow2 > prev_flow {
                    prev_flow = new_flow2;
                    stall_count = 0;
                } else {
                    stall_count += 1;
                    if stall_count >= 2 {
                        break;
                    }
                }
            }
        }

        // Finish with standard tides to guarantee correctness
        loop {
            let tc0 = Instant::now();
            let (old_flow, old_hash) = self.convergence_state();
            let tc1 = Instant::now();
            timing.convergence_us += tc1.duration_since(tc0).as_micros() as u64;

            let tb0 = Instant::now();
            self.global_relabel();
            let tb1 = Instant::now();
            timing.relabel_us += tb1.duration_since(tb0).as_micros() as u64;
            self.topology_changed = false;

            let tp0 = Instant::now();
            self.global_pull();
            let tp1 = Instant::now();
            timing.pull_us += tp1.duration_since(tp0).as_micros() as u64;

            let ts0 = Instant::now();
            self.global_push();
            let ts1 = Instant::now();
            timing.push_us += ts1.duration_since(ts0).as_micros() as u64;
            self.steps += 1;

            let tc2 = Instant::now();
            let (new_flow, new_hash) = self.convergence_state();
            let tc3 = Instant::now();
            timing.convergence_us += tc3.duration_since(tc2).as_micros() as u64;

            if new_flow == old_flow && new_hash == old_hash {
                break;
            }
        }

        let (flow, _) = self.convergence_state();
        TimedMaxFlowResult {
            result: MaxFlowResult {
                flow,
                steps: self.steps,
                edge_flow: self.flow.clone(),
            },
            timing,
        }
    }

    /// Compute a hash of the overflowing set for cheap convergence detection.
    /// Returns (net_flow, excess_hash).
    fn convergence_state(&self) -> (i64, u64) {
        // Sum flow on forward edges arriving at sink (iterate reverse edges at sink)
        let mut flow = 0i64;
        for idx in self.fwd_end[self.sink]..self.head[self.sink + 1] {
            let e = self.adj_edges[idx]; // reverse edge (odd)
            flow += self.flow[e ^ 1]; // flow on the forward edge
        }

        // Hash of overflowing set: XOR of (vertex * golden_ratio_prime + excess)
        let mut hash = 0u64;
        for v in 0..self.n {
            if v != self.source && v != self.sink && self.excess[v] > 0 {
                hash ^= (v as u64).wrapping_mul(0x9E3779B97F4A7C15)
                    ^ (self.excess[v] as u64).wrapping_mul(0x517CC1B727220A95);
            }
        }

        (flow, hash)
    }

    /// Adaptive Tide: starts with standard BFS tides, switches to hybrid
    /// (local-relabel) mode after `WARMUP` tides if the algorithm is still running.
    ///
    /// This gives the best of both worlds:
    /// - Graphs converging in few tides (layered, chains, grids) use only standard
    ///   tides with zero local-relabel overhead
    /// - Graphs needing many tides (bipartite, random, washington) get BFS-skip
    ///   savings from local relabeling
    fn tide_adaptive(&mut self) -> MaxFlowResult {
        const WARMUP: usize = 5;
        self.initialize();

        // Phase 1: standard BFS tides for warmup period
        for _ in 0..WARMUP {
            let (old_flow, old_hash) = self.convergence_state();
            if self.topology_changed {
                self.global_relabel();
            }
            self.topology_changed = false;
            self.global_pull();
            self.global_push();
            self.steps += 1;
            let (new_flow, new_hash) = self.convergence_state();
            if new_flow == old_flow && new_hash == old_hash {
                let (flow, _) = self.convergence_state();
                return MaxFlowResult {
                    flow,
                    steps: self.steps,
                    edge_flow: self.flow.clone(),
                };
            }
        }

        // Phase 2: hybrid mode — BFS tide then local-relabel tides
        let step_budget = 10 * self.n * self.n;
        let mut bfs_stall_count = 0u32;

        loop {
            // Full BFS tide
            let (old_flow, old_hash) = self.convergence_state();
            self.global_relabel();
            self.topology_changed = false;
            self.global_pull();
            self.global_push();
            self.steps += 1;

            let (bfs_flow, bfs_hash) = self.convergence_state();
            if bfs_flow == old_flow && bfs_hash == old_hash {
                break;
            }

            if bfs_flow <= old_flow {
                bfs_stall_count += 1;
                if bfs_stall_count >= 3 {
                    break;
                }
            } else {
                bfs_stall_count = 0;
            }

            if self.steps > step_budget {
                continue;
            }

            // Local-relabel tides
            let mut prev_flow = bfs_flow;
            let mut stall_count = 0u32;
            loop {
                let (_, old_hash2) = self.convergence_state();
                self.topology_changed = false;
                self.global_pull_relabel();
                self.global_push_relabel();
                self.steps += 1;

                let (new_flow2, new_hash2) = self.convergence_state();
                if new_flow2 == prev_flow && new_hash2 == old_hash2 {
                    break;
                }
                if new_flow2 > prev_flow {
                    prev_flow = new_flow2;
                    stall_count = 0;
                } else {
                    stall_count += 1;
                    if stall_count >= 2 {
                        break;
                    }
                }
            }
        }

        // Safety: finish with standard tides
        loop {
            let (old_flow, old_hash) = self.convergence_state();
            self.global_relabel();
            self.topology_changed = false;
            self.global_pull();
            self.global_push();
            self.steps += 1;
            let (new_flow, new_hash) = self.convergence_state();
            if new_flow == old_flow && new_hash == old_hash {
                break;
            }
        }

        let (flow, _) = self.convergence_state();
        MaxFlowResult {
            flow,
            steps: self.steps,
            edge_flow: self.flow.clone(),
        }
    }

    /// Adaptive Tide with per-phase timing.
    fn tide_adaptive_timed(&mut self) -> TimedMaxFlowResult {
        use std::time::Instant;
        const WARMUP: usize = 5;
        self.initialize();
        let mut timing = PhaseTiming::default();

        // Phase 1: standard BFS tides for warmup period
        for _ in 0..WARMUP {
            let t0 = Instant::now();
            let (old_flow, old_hash) = self.convergence_state();
            let t1 = Instant::now();
            timing.convergence_us += t1.duration_since(t0).as_micros() as u64;

            if self.topology_changed {
                let t2 = Instant::now();
                self.global_relabel();
                let t3 = Instant::now();
                timing.relabel_us += t3.duration_since(t2).as_micros() as u64;
            } else {
                timing.relabel_skipped += 1;
            }
            self.topology_changed = false;

            let t4 = Instant::now();
            self.global_pull();
            let t5 = Instant::now();
            timing.pull_us += t5.duration_since(t4).as_micros() as u64;

            let t6 = Instant::now();
            self.global_push();
            let t7 = Instant::now();
            timing.push_us += t7.duration_since(t6).as_micros() as u64;
            self.steps += 1;

            let t8 = Instant::now();
            let (new_flow, new_hash) = self.convergence_state();
            let t9 = Instant::now();
            timing.convergence_us += t9.duration_since(t8).as_micros() as u64;

            if new_flow == old_flow && new_hash == old_hash {
                let (flow, _) = self.convergence_state();
                return TimedMaxFlowResult {
                    result: MaxFlowResult {
                        flow,
                        steps: self.steps,
                        edge_flow: self.flow.clone(),
                    },
                    timing,
                };
            }
        }

        // Phase 2: hybrid mode
        let step_budget = 10 * self.n * self.n;
        let mut bfs_stall_count = 0u32;

        loop {
            let tc0 = Instant::now();
            let (old_flow, old_hash) = self.convergence_state();
            let tc1 = Instant::now();
            timing.convergence_us += tc1.duration_since(tc0).as_micros() as u64;

            let tb0 = Instant::now();
            self.global_relabel();
            let tb1 = Instant::now();
            timing.relabel_us += tb1.duration_since(tb0).as_micros() as u64;
            self.topology_changed = false;

            let tp0 = Instant::now();
            self.global_pull();
            let tp1 = Instant::now();
            timing.pull_us += tp1.duration_since(tp0).as_micros() as u64;

            let ts0 = Instant::now();
            self.global_push();
            let ts1 = Instant::now();
            timing.push_us += ts1.duration_since(ts0).as_micros() as u64;
            self.steps += 1;

            let tc2 = Instant::now();
            let (bfs_flow, bfs_hash) = self.convergence_state();
            let tc3 = Instant::now();
            timing.convergence_us += tc3.duration_since(tc2).as_micros() as u64;

            if bfs_flow == old_flow && bfs_hash == old_hash {
                break;
            }

            if bfs_flow <= old_flow {
                bfs_stall_count += 1;
                if bfs_stall_count >= 3 {
                    break;
                }
            } else {
                bfs_stall_count = 0;
            }

            if self.steps > step_budget {
                continue;
            }

            // Local-relabel tides
            let mut prev_flow = bfs_flow;
            let mut stall_count = 0u32;
            loop {
                let tc4 = Instant::now();
                let (_, old_hash2) = self.convergence_state();
                let tc5 = Instant::now();
                timing.convergence_us += tc5.duration_since(tc4).as_micros() as u64;
                timing.relabel_skipped += 1;
                self.topology_changed = false;

                let t2 = Instant::now();
                self.global_pull_relabel();
                let t3 = Instant::now();
                timing.pull_us += t3.duration_since(t2).as_micros() as u64;

                let t4 = Instant::now();
                self.global_push_relabel();
                let t5 = Instant::now();
                timing.push_us += t5.duration_since(t4).as_micros() as u64;
                self.steps += 1;

                let tc6 = Instant::now();
                let (new_flow2, new_hash2) = self.convergence_state();
                let tc7 = Instant::now();
                timing.convergence_us += tc7.duration_since(tc6).as_micros() as u64;

                if new_flow2 == prev_flow && new_hash2 == old_hash2 {
                    break;
                }
                if new_flow2 > prev_flow {
                    prev_flow = new_flow2;
                    stall_count = 0;
                } else {
                    stall_count += 1;
                    if stall_count >= 2 {
                        break;
                    }
                }
            }
        }

        // Safety: finish with standard tides
        loop {
            let tc0 = Instant::now();
            let (old_flow, old_hash) = self.convergence_state();
            let tc1 = Instant::now();
            timing.convergence_us += tc1.duration_since(tc0).as_micros() as u64;

            let tb0 = Instant::now();
            self.global_relabel();
            let tb1 = Instant::now();
            timing.relabel_us += tb1.duration_since(tb0).as_micros() as u64;
            self.topology_changed = false;

            let tp0 = Instant::now();
            self.global_pull();
            let tp1 = Instant::now();
            timing.pull_us += tp1.duration_since(tp0).as_micros() as u64;

            let ts0 = Instant::now();
            self.global_push();
            let ts1 = Instant::now();
            timing.push_us += ts1.duration_since(ts0).as_micros() as u64;
            self.steps += 1;

            let tc2 = Instant::now();
            let (new_flow, new_hash) = self.convergence_state();
            let tc3 = Instant::now();
            timing.convergence_us += tc3.duration_since(tc2).as_micros() as u64;

            if new_flow == old_flow && new_hash == old_hash {
                break;
            }
        }

        let (flow, _) = self.convergence_state();
        TimedMaxFlowResult {
            result: MaxFlowResult {
                flow,
                steps: self.steps,
                edge_flow: self.flow.clone(),
            },
            timing,
        }
    }

    /// Run the Tide algorithm to completion.
    fn tide(&mut self) -> MaxFlowResult {
        self.initialize();

        loop {
            let (old_flow, old_hash) = self.convergence_state();

            if self.topology_changed {
                self.global_relabel();
            }
            self.topology_changed = false;

            self.global_pull();
            self.global_push();
            self.steps += 1;

            let (new_flow, new_hash) = self.convergence_state();

            if new_flow == old_flow && new_hash == old_hash {
                break;
            }
        }

        let (flow, _) = self.convergence_state();
        MaxFlowResult {
            flow,
            steps: self.steps,
            edge_flow: self.flow.clone(),
        }
    }

    /// Run the Tide algorithm with per-phase timing instrumentation.
    fn tide_timed(&mut self) -> TimedMaxFlowResult {
        use std::time::Instant;

        self.initialize();
        let mut timing = PhaseTiming::default();

        loop {
            let t0 = Instant::now();
            let (old_flow, old_hash) = self.convergence_state();
            let t1 = Instant::now();
            timing.convergence_us += t1.duration_since(t0).as_micros() as u64;

            if self.topology_changed {
                let t2 = Instant::now();
                self.global_relabel();
                let t3 = Instant::now();
                timing.relabel_us += t3.duration_since(t2).as_micros() as u64;
            } else {
                timing.relabel_skipped += 1;
            }
            self.topology_changed = false;

            let t4 = Instant::now();
            self.global_pull();
            let t5 = Instant::now();
            timing.pull_us += t5.duration_since(t4).as_micros() as u64;

            let t6 = Instant::now();
            self.global_push();
            let t7 = Instant::now();
            timing.push_us += t7.duration_since(t6).as_micros() as u64;

            self.steps += 1;

            let t8 = Instant::now();
            let (new_flow, new_hash) = self.convergence_state();
            let t9 = Instant::now();
            timing.convergence_us += t9.duration_since(t8).as_micros() as u64;

            if new_flow == old_flow && new_hash == old_hash {
                break;
            }
        }

        let (flow, _) = self.convergence_state();
        TimedMaxFlowResult {
            result: MaxFlowResult {
                flow,
                steps: self.steps,
                edge_flow: self.flow.clone(),
            },
            timing,
        }
    }
}

/// Compute the maximum flow in a network using the Tide algorithm.
///
/// # Arguments
/// * `n` - Number of vertices (0-indexed)
/// * `source` - Source vertex
/// * `sink` - Sink vertex
/// * `edges` - List of (from, to, capacity) tuples
///
/// # Returns
/// A `MaxFlowResult` containing the max flow value, step count, and per-edge flows.
///
/// # Example
/// ```
/// use tide_maxflow::max_flow;
/// let result = max_flow(4, 0, 3, &[(0, 1, 10), (0, 2, 10), (1, 3, 10), (2, 3, 10)]);
/// assert_eq!(result.flow, 20);
/// ```
pub fn max_flow(
    n: usize,
    source: usize,
    sink: usize,
    edges: &[(usize, usize, i64)],
) -> MaxFlowResult {
    let mut state = TideState::new(n, source, sink, edges);
    state.tide()
}

/// Compute max flow with per-phase timing instrumentation.
///
/// Same as `max_flow` but returns a `TimedMaxFlowResult` with timing breakdown
/// for each phase (relabel, pull, push, convergence check).
pub fn max_flow_timed(
    n: usize,
    source: usize,
    sink: usize,
    edges: &[(usize, usize, i64)],
) -> TimedMaxFlowResult {
    let mut state = TideState::new(n, source, sink, edges);
    state.tide_timed()
}

/// Compute max flow using the hybrid Tide algorithm.
///
/// Full BFS on the first tide, then local relabeling during pull/push
/// until flow stalls, then full BFS again. Terminates when a full BFS
/// tide makes no progress.
pub fn max_flow_hybrid(
    n: usize,
    source: usize,
    sink: usize,
    edges: &[(usize, usize, i64)],
) -> MaxFlowResult {
    let mut state = TideState::new(n, source, sink, edges);
    state.tide_hybrid()
}

/// Compute max flow using adaptive Tide.
///
/// Starts with standard BFS tides. If the algorithm hasn't converged after
/// 5 tides, switches to hybrid mode (local relabeling between BFS tides).
/// Best of both worlds: no overhead on easy graphs, BFS-skip savings on hard ones.
pub fn max_flow_adaptive(
    n: usize,
    source: usize,
    sink: usize,
    edges: &[(usize, usize, i64)],
) -> MaxFlowResult {
    let mut state = TideState::new(n, source, sink, edges);
    state.tide_adaptive()
}

/// Compute max flow using adaptive Tide with per-phase timing.
pub fn max_flow_adaptive_timed(
    n: usize,
    source: usize,
    sink: usize,
    edges: &[(usize, usize, i64)],
) -> TimedMaxFlowResult {
    let mut state = TideState::new(n, source, sink, edges);
    state.tide_adaptive_timed()
}

/// Compute max flow using the hybrid Tide algorithm with timing.
pub fn max_flow_hybrid_timed(
    n: usize,
    source: usize,
    sink: usize,
    edges: &[(usize, usize, i64)],
) -> TimedMaxFlowResult {
    let mut state = TideState::new(n, source, sink, edges);
    state.tide_hybrid_timed()
}

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

    #[test]
    fn test_simple_path() {
        let result = max_flow(3, 0, 2, &[(0, 1, 10), (1, 2, 5)]);
        assert_eq!(result.flow, 5);
    }

    #[test]
    fn test_parallel_paths() {
        let result = max_flow(4, 0, 3, &[(0, 1, 10), (1, 3, 10), (0, 2, 15), (2, 3, 15)]);
        assert_eq!(result.flow, 25);
    }

    #[test]
    fn test_diamond() {
        let result = max_flow(
            4,
            0,
            3,
            &[(0, 1, 10), (0, 2, 10), (1, 2, 1), (1, 3, 10), (2, 3, 10)],
        );
        assert_eq!(result.flow, 20);
    }

    #[test]
    fn test_bottleneck() {
        let result = max_flow(4, 0, 3, &[(0, 1, 100), (1, 2, 1), (2, 3, 100)]);
        assert_eq!(result.flow, 1);
    }

    #[test]
    fn test_no_path() {
        let result = max_flow(4, 0, 3, &[(0, 1, 10), (2, 3, 10)]);
        assert_eq!(result.flow, 0);
    }

    #[test]
    fn test_single_edge() {
        let result = max_flow(2, 0, 1, &[(0, 1, 42)]);
        assert_eq!(result.flow, 42);
    }

    #[test]
    fn test_multiple_paths_shared_edge() {
        let result = max_flow(
            5,
            0,
            4,
            &[(0, 1, 10), (0, 2, 10), (1, 3, 10), (2, 3, 10), (3, 4, 15)],
        );
        assert_eq!(result.flow, 15);
    }

    #[test]
    fn test_layered_graph() {
        let mut edges = Vec::new();
        for i in 0..3 {
            edges.push((0, i + 1, 10));
        }
        for i in 0..3 {
            for j in 0..3 {
                edges.push((i + 1, j + 4, 5));
            }
        }
        for j in 0..3 {
            edges.push((j + 4, 7, 10));
        }
        let result = max_flow(8, 0, 7, &edges);
        assert_eq!(result.flow, 30);
    }

    #[test]
    fn test_clrs() {
        let result = max_flow(
            6,
            0,
            5,
            &[
                (0, 1, 16),
                (0, 2, 13),
                (1, 2, 10),
                (1, 3, 12),
                (2, 1, 4),
                (2, 4, 14),
                (3, 2, 9),
                (3, 5, 20),
                (4, 3, 7),
                (4, 5, 4),
            ],
        );
        assert_eq!(result.flow, 23);
    }
}