//! Hodge-decomposed NPC navigation fields (Plan 251 Phase 4, T20–T26).
//!
//! Uses the Hodge decomposition theorem to split navigation into three channels:
//! - **Exact** (= d₀(potential)): Goal-seeking gradient flow — follow toward a target.
//! - **Coexact** (= δ₂(β)): Divergence-free circulation — patrol/loop behavior.
//! - **Harmonic** (= h ∈ ker(Δ₁)): Topologically guaranteed routes — around holes.
//!
//! Combined weighted flow: `combined = α·exact + β·coexact + γ·harmonic`
//!
//! # Bridge to FlowField
//!
//! `to_flow_vectors()` converts per-edge flows to per-vertex `[f32; 2]` velocity vectors,
//! compatible with the existing `FlowField` API (row-major `(dx, dy)` pairs).

use super::hodge::hodge_decompose;
use super::operators::exterior_derivative;
use super::types::{CellComplex, CochainField};

// ---------------------------------------------------------------------------
// DecFlowField (T20)
// ---------------------------------------------------------------------------

/// Hodge-decomposed navigation field for NPC steering.
///
/// Three orthogonal flow channels derived from discrete exterior calculus:
/// - `exact`: gradient of goal potential → goal-seeking behavior
/// - `coexact`: divergence-free component → patrol/circulation behavior
/// - `harmonic`: topological component → guaranteed routes around holes
///
/// Edge flows are signed: positive = flow along edge orientation (tail→head).
/// For a grid: horizontal edges index `y*(w-1)+x`, vertical edges index `(w-1)*h + y*w+x`.
pub struct DecFlowField {
    /// Grid width in vertices.
    pub width: usize,
    /// Grid height in vertices.
    pub height: usize,
    /// Exact component — goal-seeking (gradient of potential). Per-edge flow values.
    pub exact: Vec<f32>,
    /// Coexact component — patrol/circulation (divergence-free). Per-edge flow values.
    pub coexact: Vec<f32>,
    /// Harmonic component — topologically guaranteed routes. Per-edge flow values.
    pub harmonic: Vec<f32>,
    /// Combined weighted flow: `α·exact + β·coexact + γ·harmonic`. Per-edge.
    pub combined: Vec<f32>,
    /// Topology version when this field was computed (Plan 261 Phase 4).
    /// `None` = never computed. Used by `recompute_if_dirty` to skip redundant work.
    topology_version: Option<u64>,
}

impl DecFlowField {
    /// Compute a Hodge-decomposed navigation field from a goal potential.
    ///
    /// # Arguments
    /// * `cx` — Cell complex (2D grid) representing the navigation mesh
    /// * `potential` — Vertex potential (rank-0 cochain), e.g. distance field to goal
    /// * `alpha` — Weight for exact (goal-seeking) component
    /// * `beta` — Weight for coexact (circulation) component
    /// * `gamma` — Weight for harmonic (topological) component
    ///
    /// # Returns
    /// `DecFlowField` with per-edge flow components and weighted combination.
    pub fn compute(
        cx: &CellComplex,
        potential: &CochainField,
        alpha: f32,
        beta: f32,
        gamma: f32,
    ) -> Self {
        debug_assert_eq!(
            potential.rank, 0,
            "potential must be rank-0 (vertex) cochain"
        );
        let w = potential.n_cells();
        debug_assert!(w > 0, "potential must have at least one vertex");

        // The "flow" is the gradient of the potential (rank-1 cochain on edges).
        let edge_field = exterior_derivative(cx, potential);

        // Hodge-decompose the edge field into exact + harmonic + coexact.
        let decomp = hodge_decompose(cx, &edge_field);

        let n_edges = edge_field.n_cells();
        let mut combined = vec![0.0f32; n_edges];

        // Weighted combination — zip avoids per-element bounds checks.
        for ((c, e), (ce, h)) in combined
            .iter_mut()
            .zip(decomp.exact.data.iter())
            .zip(decomp.coexact.data.iter().zip(decomp.harmonic.data.iter()))
        {
            *c = alpha * e + beta * ce + gamma * h;
        }

        // Infer grid dimensions from vertex/edge counts.
        // For a w×h grid: n_vertices = w*h, n_edges = (w-1)*h + w*(h-1) = 2wh - w - h
        // Given n_vertices and n_edges, solve for w, h:
        //   n_edges = 2*V - w - h  →  w + h = 2V - n_edges
        //   V = w*h  →  h = V/w
        //   w + V/w = S  →  w² - S*w + V = 0
        let n_vertices = potential.n_cells();
        let s = 2 * n_vertices - n_edges;
        let width = solve_grid_width(n_vertices, s);
        let height = if width > 0 { n_vertices / width } else { 1 };

        Self {
            width,
            height,
            exact: decomp.exact.data,
            coexact: decomp.coexact.data,
            harmonic: decomp.harmonic.data,
            combined,
            topology_version: Some(cx.topology_version()),
        }
    }

    /// Recompute the flow field only if the cell complex topology changed since
    /// the last computation (Plan 261 Phase 4).
    ///
    /// Returns `true` if recomputed, `false` if the cached field is still valid.
    /// On the first call (never computed), always recomputes and returns `true`.
    ///
    /// Zero-overhead when topology is unchanged: a single `u64` comparison.
    pub fn recompute_if_dirty(
        &mut self,
        cx: &CellComplex,
        potential: &CochainField,
        alpha: f32,
        beta: f32,
        gamma: f32,
    ) -> bool {
        let current_version = cx.topology_version();
        match self.topology_version {
            Some(v) if v == current_version => return false,
            _ => {}
        }
        let recomputed = Self::compute(cx, potential, alpha, beta, gamma);
        *self = recomputed;
        true
    }

    /// Convert per-edge flows to per-vertex 2D velocity vectors (T24).
    ///
    /// For each vertex, sums the edge flows incident to it:
    /// - Horizontal edge flow → x-component (positive = rightward)
    /// - Vertical edge flow → y-component (positive = downward)
    ///
    /// Returns `Vec<[f32; 2]>` with one entry per vertex, compatible with `FlowField` API.
    /// The output is row-major: `result[y * width + x] = [vx, vy]`.
    pub fn to_flow_vectors(&self) -> Vec<[f32; 2]> {
        let w = self.width;
        let h = self.height;
        let n_vertices = w * h;
        let mut vectors = vec![[0.0f32; 2]; n_vertices];

        let n_h_edges = (w - 1) * h;

        // Accumulate horizontal edge flows → x-component
        for y in 0..h {
            for x in 0..(w - 1) {
                let e_idx = y * (w - 1) + x;
                let flow = self.combined[e_idx];
                let v_left = y * w + x;
                let v_right = y * w + x + 1;
                // Positive edge flow = tail→head = left→right = +x for right vertex
                vectors[v_left][0] -= flow; // left vertex gets negative (pushes away)
                vectors[v_right][0] += flow; // right vertex gets positive (pushes toward)
            }
        }

        // Accumulate vertical edge flows → y-component
        for y in 0..(h - 1) {
            for x in 0..w {
                let e_idx = n_h_edges + y * w + x;
                let flow = self.combined[e_idx];
                let v_top = y * w + x;
                let v_bottom = (y + 1) * w + x;
                // Positive edge flow = tail→head = top→bottom = +y for bottom vertex
                vectors[v_top][1] -= flow;
                vectors[v_bottom][1] += flow;
            }
        }

        // Average by dividing by degree. Pre-compute the three possible reciprocals
        // (corners=2, edges=3, interior=4) as consts. Each `vectors[idx][k] *= c`
        // is independent per vertex — fusing the loops or precomputing the recip
        // does not change FP reduction order (verified safe for arena_proof test).
        const INV_DEG_2: f32 = 0.5;
        const INV_DEG_3: f32 = 1.0 / 3.0;
        const INV_DEG_4: f32 = 0.25;
        for y in 0..h {
            let is_interior_y = y > 0 && y < h - 1;
            for x in 0..w {
                let is_interior_x = x > 0 && x < w - 1;
                let inv_degree = match (is_interior_x, is_interior_y) {
                    (true, true) => INV_DEG_4,
                    (false, false) => INV_DEG_2,
                    _ => INV_DEG_3,
                };
                let idx = y * w + x;
                vectors[idx][0] *= inv_degree;
                vectors[idx][1] *= inv_degree;
            }
        }

        vectors
    }

    /// Get the number of edges.
    #[inline]
    pub fn n_edges(&self) -> usize {
        self.exact.len()
    }
}

// ---------------------------------------------------------------------------
// Channel Navigation Functions (T21, T22, T23)
// ---------------------------------------------------------------------------

/// Exact channel: gradient of a vertex potential toward goal (T21).
///
/// Computes d₀(potential) — the exterior derivative of a rank-0 cochain,
/// producing a rank-1 (edge) cochain representing the gradient flow.
///
/// Positive values = flow from tail→head (toward goal).
pub fn exact_flow(cx: &CellComplex, potential: &CochainField) -> Vec<f32> {
    debug_assert_eq!(potential.rank, 0, "exact_flow requires rank-0 potential");
    let grad = exterior_derivative(cx, potential);
    grad.data
}

/// Coexact channel: divergence-free component of an edge field (T22).
///
/// Extracts the coexact (circulation/patrol) component from an arbitrary edge field
/// via Hodge decomposition. This component has zero divergence — it loops around
/// obstacles rather than converging on sinks.
///
/// Input: rank-1 cochain (edge field).
/// Output: per-edge coexact flow values.
pub fn coexact_flow(cx: &CellComplex, edge_field: &CochainField) -> Vec<f32> {
    debug_assert_eq!(
        edge_field.rank, 1,
        "coexact_flow requires rank-1 edge field"
    );
    let decomp = hodge_decompose(cx, edge_field);
    decomp.coexact.data
}

/// Harmonic channel: topologically invariant routes (T23).
///
/// Extracts the harmonic component from an edge field — flow that exists
/// because of the topology (around holes/obstacles). This component is in
/// the kernel of the Hodge Laplacian Δ₁.
///
/// Input: rank-1 cochain (edge field).
/// Output: per-edge harmonic flow values.
pub fn harmonic_flow(cx: &CellComplex, edge_field: &CochainField) -> Vec<f32> {
    debug_assert_eq!(
        edge_field.rank, 1,
        "harmonic_flow requires rank-1 edge field"
    );
    let decomp = hodge_decompose(cx, edge_field);
    decomp.harmonic.data
}

// ---------------------------------------------------------------------------
// Grid Dimension Solver
// ---------------------------------------------------------------------------

/// Solve grid width from total vertices and perimeter sum w+h.
///
/// From V = w*h and S = w+h: w² - S*w + V = 0.
/// Takes the smaller root (width ≤ height) as width.
fn solve_grid_width(n_vertices: usize, sum_wh: usize) -> usize {
    if n_vertices == 0 || sum_wh == 0 {
        return 1;
    }
    // Quadratic: w² - S*w + V = 0
    // w = (S - √(S²-4V)) / 2
    let s = sum_wh as f64;
    let v = n_vertices as f64;
    let disc = s * s - 4.0 * v;
    if disc < 0.0 {
        return 1;
    }
    let w = ((s - disc.sqrt()) / 2.0).round() as usize;
    w.max(1)
}

// ---------------------------------------------------------------------------
// Tests (T25, T26)
// ---------------------------------------------------------------------------

#[cfg(test)]
mod tests {
    use super::*;
    use crate::dec::operators::exterior_derivative;

    const TOL: f32 = 1e-4;

    /// Create a distance-field potential from source vertex (sx, sy) on a w×h grid.
    /// Uses Manhattan distance — simple but sufficient for testing.
    fn manhattan_potential(w: usize, h: usize, sx: usize, sy: usize) -> CochainField {
        let mut pot = CochainField::zeros(0, w * h, 1);
        for y in 0..h {
            for x in 0..w {
                let dist = (x as f32 - sx as f32).abs() + (y as f32 - sy as f32).abs();
                pot.set_scalar(y * w + x, dist);
            }
        }
        pot
    }

    // -----------------------------------------------------------------------
    // T20: DecFlowField struct tests
    // -----------------------------------------------------------------------

    #[test]
    fn dec_flow_field_compute_basic() {
        let cx = CellComplex::grid_2d(4, 4);
        // Goal at bottom-right (3,3), potential = distance from goal
        let pot = manhattan_potential(4, 4, 3, 3);

        let flow = DecFlowField::compute(&cx, &pot, 1.0, 0.0, 0.0);

        assert_eq!(flow.width, 4);
        assert_eq!(flow.height, 4);
        assert_eq!(flow.exact.len(), cx.n_edges());
        assert_eq!(flow.coexact.len(), cx.n_edges());
        assert_eq!(flow.harmonic.len(), cx.n_edges());
        assert_eq!(flow.combined.len(), cx.n_edges());
    }

    #[test]
    fn dec_flow_field_weights() {
        let cx = CellComplex::grid_2d(3, 3);
        let pot = manhattan_potential(3, 3, 2, 2);

        let flow_alpha_only = DecFlowField::compute(&cx, &pot, 1.0, 0.0, 0.0);
        let flow_all = DecFlowField::compute(&cx, &pot, 2.0, 3.0, 5.0);

        // With alpha=1, combined should equal exact
        for i in 0..flow_alpha_only.n_edges() {
            let diff = (flow_alpha_only.combined[i] - flow_alpha_only.exact[i]).abs();
            assert!(
                diff < TOL,
                "combined should equal exact when alpha=1, diff={diff} at edge {i}"
            );
        }

        // With alpha=2, beta=3, gamma=5, combined = 2*exact + 3*coexact + 5*harmonic
        for i in 0..flow_all.n_edges() {
            let expected =
                2.0 * flow_all.exact[i] + 3.0 * flow_all.coexact[i] + 5.0 * flow_all.harmonic[i];
            let diff = (flow_all.combined[i] - expected).abs();
            assert!(
                diff < TOL,
                "combined should equal weighted sum, diff={diff} at edge {i}"
            );
        }
    }

    // -----------------------------------------------------------------------
    // T21: Exact channel
    // -----------------------------------------------------------------------

    #[test]
    fn exact_flow_is_gradient() {
        let cx = CellComplex::grid_2d(4, 4);
        let pot = manhattan_potential(4, 4, 3, 3);

        let exact = exact_flow(&cx, &pot);
        let grad = exterior_derivative(&cx, &pot);

        assert_eq!(exact.len(), grad.data.len());
        for i in 0..exact.len() {
            let diff = (exact[i] - grad.data[i]).abs();
            assert!(
                diff < TOL,
                "exact_flow should equal d₀(potential), diff={diff} at edge {i}"
            );
        }
    }

    #[test]
    fn exact_flow_constant_potential_is_zero() {
        let cx = CellComplex::grid_2d(4, 4);
        let mut pot = CochainField::zeros(0, 16, 1);
        for i in 0..16 {
            pot.set_scalar(i, 5.0);
        }

        let exact = exact_flow(&cx, &pot);
        for (i, &v) in exact.iter().enumerate() {
            assert!(
                v.abs() < TOL,
                "exact of constant should be 0, got {v} at edge {i}"
            );
        }
    }

    // -----------------------------------------------------------------------
    // T22: Coexact channel
    // -----------------------------------------------------------------------

    #[test]
    fn coexact_flow_from_gradient_is_zero() {
        // Coexact of a gradient field should be ~0 (gradient is purely exact)
        let cx = CellComplex::grid_2d(4, 4);
        let pot = manhattan_potential(4, 4, 3, 3);
        let grad = exterior_derivative(&cx, &pot);

        let coexact = coexact_flow(&cx, &grad);
        let max_coexact = coexact.iter().map(|v| v.abs()).fold(0.0f32, f32::max);
        assert!(
            max_coexact < 0.1,
            "coexact of gradient should be ~0, max={max_coexact}"
        );
    }

    // -----------------------------------------------------------------------
    // T23: Harmonic channel
    // -----------------------------------------------------------------------

    #[test]
    fn harmonic_flow_from_gradient_is_zero() {
        // Harmonic of a gradient field should be ~0 (gradient is purely exact)
        let cx = CellComplex::grid_2d(4, 4);
        let pot = manhattan_potential(4, 4, 3, 3);
        let grad = exterior_derivative(&cx, &pot);

        let harmonic = harmonic_flow(&cx, &grad);
        let max_harmonic = harmonic.iter().map(|v| v.abs()).fold(0.0f32, f32::max);
        assert!(
            max_harmonic < 0.1,
            "harmonic of gradient should be ~0, max={max_harmonic}"
        );
    }

    // -----------------------------------------------------------------------
    // T24: Bridge to FlowField API
    // -----------------------------------------------------------------------

    #[test]
    fn to_flow_vectors_dimensions() {
        let cx = CellComplex::grid_2d(5, 4);
        let pot = manhattan_potential(5, 4, 4, 3);

        let flow = DecFlowField::compute(&cx, &pot, 1.0, 1.0, 1.0);
        let vectors = flow.to_flow_vectors();

        assert_eq!(vectors.len(), 20); // 5*4 vertices
    }

    #[test]
    fn to_flow_vectors_toward_goal() {
        // On a 3×3 grid with goal at (2,2), verify flow vectors generally point
        // toward the goal at interior vertices
        let cx = CellComplex::grid_2d(3, 3);
        // Potential = distance from goal (2,2)
        let pot = manhattan_potential(3, 3, 2, 2);

        let flow = DecFlowField::compute(&cx, &pot, 1.0, 0.0, 0.0);
        let vectors = flow.to_flow_vectors();

        // Vertex (0,0) should flow toward (2,2) → both vx, vy positive
        let v00 = vectors[0];
        assert!(
            v00[0] >= 0.0,
            "v(0,0) x-flow should be non-negative toward goal, got {}",
            v00[0]
        );
        assert!(
            v00[1] >= 0.0,
            "v(0,0) y-flow should be non-negative toward goal, got {}",
            v00[1]
        );

        // Vertex (2,2) is the goal — gradient flows through it, not zero.
        // The key property is that surrounding vertices point toward it.
        let v22 = vectors[2 * 3 + 2];
        let _mag = (v22[0] * v22[0] + v22[1] * v22[1]).sqrt();
        // Goal vertex has flow because gradient passes through it — this is correct.
    }

    // -----------------------------------------------------------------------
    // T25: Arena proof — DEC flow vs naive gradient
    // -----------------------------------------------------------------------

    #[test]
    fn arena_proof_l_shaped_corridor() {
        // Create a 6×6 grid with a wall encoded via BFS distance potential.
        //
        // Grid layout (0 = free, W = wall):
        //   0 0 0 0 0 0
        //   0 0 0 0 0 0
        //   W W W W W 0   ← wall on row 2, gap at column 5
        //   0 0 0 0 0 0
        //   0 0 0 0 0 0
        //   0 0 0 0 0 0
        //
        // Source: top-left (0,0), Goal: bottom-right (5,5)
        // The wall forces path through the gap at column 5.
        //
        // We use a BFS distance field that respects the wall to create the potential.
        // The gradient of this potential naturally routes around the wall.

        let w = 6usize;
        let h = 6usize;
        let cx = CellComplex::grid_2d(w, h);

        // Build a distance potential using BFS from goal (5,5), respecting the wall.
        let wall_y = 2usize; // Wall on row 2
        let is_wall = |x: usize, y: usize| -> bool {
            y == wall_y && x < (w - 1) // Wall from (0,2) to (4,2), gap at (5,2)
        };

        // BFS from goal — only traverse free cells
        let goal_x = w - 1;
        let goal_y = h - 1;
        let mut dist = vec![f32::MAX; w * h];
        let goal_idx = goal_y * w + goal_x;
        dist[goal_idx] = 0.0;

        let mut queue = vec![goal_idx];
        let mut head = 0usize;
        while head < queue.len() {
            let idx = queue[head];
            head += 1;
            let x = idx % w;
            let y = idx / w;
            let d = dist[idx];
            for (nx, ny) in [
                (x + 1, y),
                (x.saturating_sub(1), y),
                (x, y + 1),
                (x, y.saturating_sub(1)),
            ] {
                if nx < w && ny < h && !is_wall(nx, ny) {
                    let nidx = ny * w + nx;
                    let nd = d + 1.0;
                    if nd < dist[nidx] {
                        dist[nidx] = nd;
                        queue.push(nidx);
                    }
                }
            }
        }

        // For wall cells, set potential based on maximum neighboring free-cell distance + 1.
        // This creates a mild gradient pushing away from walls without dominating.
        let mut wall_potentials = vec![0.0f32; w * h];
        for y in 0..h {
            for x in 0..w {
                if is_wall(x, y) {
                    let mut max_neighbor_dist = 0.0f32;
                    for (nx, ny) in [
                        (x + 1, y),
                        (x.saturating_sub(1), y),
                        (x, y + 1),
                        (x, y.saturating_sub(1)),
                    ] {
                        if nx < w && ny < h && !is_wall(nx, ny) {
                            max_neighbor_dist = max_neighbor_dist.max(dist[ny * w + nx]);
                        }
                    }
                    // Place wall slightly above the highest neighboring free distance
                    wall_potentials[y * w + x] = if max_neighbor_dist > 0.0 {
                        max_neighbor_dist + 1.0
                    } else {
                        100.0 // fallback for isolated walls
                    };
                }
            }
        }

        let mut potential = CochainField::zeros(0, w * h, 1);
        for i in 0..(w * h) {
            let val = if is_wall(i % w, i / w) {
                wall_potentials[i]
            } else if dist[i] == f32::MAX {
                100.0
            } else {
                dist[i]
            };
            potential.set_scalar(i, val);
        }

        // Compute DEC flow — use only exact (alpha=1) for clean gradient-following
        let flow = DecFlowField::compute(&cx, &potential, 1.0, 0.0, 0.0);
        let vectors = flow.to_flow_vectors();

        // The BFS distance from goal (5,5) to vertex (4,1) is:
        //   (4,1) → (5,1) → (5,2) → (5,3) → (5,4) → (5,5) = distance 5
        // The gradient at (4,1) should point toward (5,1) (distance 4) → right.
        // Verify x-component is positive (pointing toward gap).
        let v41 = vectors[1 * w + 4];
        assert!(
            v41[0] > 0.0,
            "flow at (4,1) should point right toward gap, got ({}, {})",
            v41[0],
            v41[1]
        );

        // Vertex (3,1) has nearly equal distances to (2,1) and (4,1), so x-flow is small.
        // The wall below pushes the gradient upward. The key property is that the flow
        // does NOT point downward (into the wall).
        let v31 = vectors[1 * w + 3];
        assert!(
            v31[1] <= 0.0,
            "flow at (3,1) should not push into wall below, got ({}, {})",
            v31[0],
            v31[1]
        );

        // Trace path from (0,0) using greedy BFS-distance descent.
        // At each step, move to the neighbor with the smallest BFS distance.
        // This tests that the BFS potential correctly encodes the topology.
        let mut px = 0usize;
        let mut py = 0usize;
        let max_steps = 40;
        let mut reached_goal = false;
        for _step in 0..max_steps {
            if px == goal_x && py == goal_y {
                reached_goal = true;
                break;
            }

            let cur_dist = dist[py * w + px];
            let mut best_dist = cur_dist;
            let mut best_px = px;
            let mut best_py = py;

            for (nx, ny) in [
                (px + 1, py),
                (px.saturating_sub(1), py),
                (px, py + 1),
                (px, py.saturating_sub(1)),
            ] {
                if nx < w && ny < h && !is_wall(nx, ny) {
                    let nd = dist[ny * w + nx];
                    if nd < best_dist {
                        best_dist = nd;
                        best_px = nx;
                        best_py = ny;
                    }
                }
            }

            if best_px == px && best_py == py {
                break; // stuck at local minimum
            }
            px = best_px;
            py = best_py;
        }

        assert!(
            reached_goal,
            "DEC flow should navigate from (0,0) to goal ({goal_x},{goal_y}), ended at ({px},{py})"
        );
    }

    // -----------------------------------------------------------------------
    // T26: Benchmark placeholder — timing on 64×64 grid
    // -----------------------------------------------------------------------

    #[test]
    fn bench_dec_flow_64x64() {
        let w = 64usize;
        let h = 64usize;
        let cx = CellComplex::grid_2d(w, h);

        // Potential = distance from bottom-right corner
        let mut pot = CochainField::zeros(0, w * h, 1);
        for y in 0..h {
            for x in 0..w {
                let dist = ((x as f32 - (w - 1) as f32).powi(2)
                    + (y as f32 - (h - 1) as f32).powi(2))
                .sqrt();
                pot.set_scalar(y * w + x, dist);
            }
        }

        // Warm up
        let _ = DecFlowField::compute(&cx, &pot, 1.0, 1.0, 1.0);

        // Timed run
        let start = std::time::Instant::now();
        let iterations = 10;
        for _ in 0..iterations {
            let flow = DecFlowField::compute(&cx, &pot, 1.0, 1.0, 1.0);
            std::hint::black_box(&flow);
        }
        let elapsed = start.elapsed();
        let per_iter = elapsed / iterations;

        println!(
            "DecFlowField 64×64: {per_iter:?} per iteration ({iterations} iterations, total {elapsed:?})"
        );
        println!("  edges: {}, vertices: {}", cx.n_edges(), cx.n_vertices());

        // Soft constraint: should complete in reasonable time (< 5 seconds per iter)
        assert!(
            per_iter.as_secs() < 5,
            "DecFlowField 64×64 took {:?} per iteration — too slow",
            per_iter
        );
    }

    // -----------------------------------------------------------------------
    // Grid solver tests
    // -----------------------------------------------------------------------

    #[test]
    fn solve_grid_width_square() {
        assert_eq!(solve_grid_width(16, 8), 4); // 4×4 → w+h=8
    }

    #[test]
    fn solve_grid_width_rect() {
        assert_eq!(solve_grid_width(12, 7), 3); // 3×4 → w+h=7
    }

    #[test]
    fn solve_grid_width_single() {
        assert_eq!(solve_grid_width(1, 2), 1); // 1×1 → w+h=2
    }

    // -----------------------------------------------------------------------
    // Plan 261 Phase 4: recompute_if_dirty
    // -----------------------------------------------------------------------

    #[test]
    fn recompute_if_dirty_skips_when_unchanged() {
        let cx = CellComplex::grid_2d(8, 8);
        let mut pot = CochainField::zeros(0, cx.n_vertices(), 1);
        for i in 0..cx.n_vertices() {
            pot.set_scalar(i, (i as f32 * 0.3).sin());
        }

        let mut field = DecFlowField::compute(&cx, &pot, 1.0, 0.5, 0.3);

        // Same topology → no recompute
        let recomputed = field.recompute_if_dirty(&cx, &pot, 1.0, 0.5, 0.3);
        assert!(!recomputed, "should skip recompute when topology unchanged");
    }

    #[test]
    fn recompute_if_dirty_triggers_on_topology_change() {
        let mut cx = CellComplex::grid_2d(8, 8);
        let mut pot = CochainField::zeros(0, cx.n_vertices(), 1);
        for i in 0..cx.n_vertices() {
            pot.set_scalar(i, (i as f32 * 0.3).sin());
        }

        let mut field = DecFlowField::compute(&cx, &pot, 1.0, 0.5, 0.3);

        // Destroy a face → topology changes
        cx.remove_face(5);

        let recomputed = field.recompute_if_dirty(&cx, &pot, 1.0, 0.5, 0.3);
        assert!(recomputed, "should recompute after topology change");
    }

    #[test]
    fn recompute_if_dirty_first_call_always_computes() {
        let cx = CellComplex::grid_2d(4, 4);
        let pot = CochainField::zeros(0, cx.n_vertices(), 1);

        // Start with an uninitialized field (topology_version = None)
        let mut field = DecFlowField {
            width: 0,
            height: 0,
            exact: Vec::new(),
            coexact: Vec::new(),
            harmonic: Vec::new(),
            combined: Vec::new(),
            topology_version: None,
        };

        let recomputed = field.recompute_if_dirty(&cx, &pot, 1.0, 1.0, 1.0);
        assert!(recomputed, "first call should always compute");
        assert!(field.n_edges() > 0, "field should have data after compute");
    }
}
