//! Frontage-first block subdivision (spec §3.1, milestone 0.2).
//!
//! For each block face we:
//!
//! 1. Identify "real" corners — block-boundary vertices whose
//!    underlying graph node has degree ≥3, or degree 2 with a bend
//!    angle below 150°. Smooth continuations (collinear-roads case)
//!    are not real corners.
//! 2. At each real corner, build a corner parcel: a parallelogram
//!    extending `frontage_width` along each adjacent block edge and
//!    closing diagonally inward. The longer adjacent edge becomes the
//!    parcel's frontage; the other road-adjacent edge is reclassified
//!    as `Side`. (See §11 D9–D11.)
//! 3. For each block boundary edge, walk the *remaining* frontage
//!    (between the corner footprints, if any) at arc-length intervals
//!    of approximately `params.frontage_width`, extruding perpendicular
//!    to the block interior. No more bisector clipping — the corners
//!    take care of corner geometry up front.
//!
//! Determinism (invariant **I6**) is preserved: random jitter is
//! drawn from a `ChaCha8Rng` whose seed is mixed from `params.seed`
//! and the road-segment id; corner parcels are pure functions of
//! geometry.

use std::time::{Duration, Instant};

use glam::DVec2;
use rand::{Rng, SeedableRng};
use rand_chacha::ChaCha8Rng;
use slotmap::Key;

use crate::config::SubdivisionParams;
use crate::error::SubdivisionError;
use crate::geometry::{Polygon, EPS_GEOM};
use crate::network::blocks::{extract_blocks, Block};
use crate::network::graph::FaceId;
use crate::network::{NodeId, RoadGraph, RoadId};

use super::classify::classify_edges;
use super::regularize::regularize_parcel;
use super::{Parcel, ParcelSet};

/// Per-phase timing collected by [`subdivide_all_with_stats`].
///
/// Phases are non-overlapping wall-clock spans. `total` is measured
/// independently and may be slightly larger than the sum of the
/// phases due to overhead between phases.
#[derive(Debug, Clone, Default)]
pub struct SubdivisionStats {
    /// Number of bounded faces (blocks) extracted from the graph.
    pub block_count: usize,
    /// Total parcels produced across all blocks.
    pub parcel_count: usize,
    /// Time to (re)build the half-edge / DCEL topology if needed.
    pub topology_rebuild: Duration,
    /// Time to extract bounded faces into [`Block`] objects.
    pub block_extraction: Duration,
    /// Cumulative time spent in [`subdivide_block`] across all blocks.
    pub block_subdivide_total: Duration,
    /// Wall-clock total of [`subdivide_all_with_stats`].
    pub total: Duration,
}

impl SubdivisionStats {
    /// Average parcels generated per second over the whole call.
    #[must_use]
    pub fn parcels_per_second(&self) -> f64 {
        let secs = self.total.as_secs_f64();
        if secs > 0.0 {
            self.parcel_count as f64 / secs
        } else {
            0.0
        }
    }

    /// Average wall-clock time per parcel, in microseconds.
    #[must_use]
    pub fn time_per_parcel_us(&self) -> f64 {
        if self.parcel_count > 0 {
            (self.total.as_nanos() as f64 / self.parcel_count as f64) / 1_000.0
        } else {
            0.0
        }
    }
}

/// Subdivide every bounded face of `graph` into parcels.
///
/// # Errors
///
/// Returns a [`SubdivisionError`] if `params` is invalid, the graph
/// is non-planar, or an internal geometric failure occurs.
pub fn subdivide_all(
    graph: &RoadGraph,
    params: &SubdivisionParams,
) -> Result<ParcelSet, SubdivisionError> {
    subdivide_all_with_stats(graph, params).map(|(p, _)| p)
}

/// Subdivide every bounded face and return per-phase performance
/// statistics alongside the parcel set.
///
/// # Errors
///
/// Same as [`subdivide_all`].
pub fn subdivide_all_with_stats(
    graph: &RoadGraph,
    params: &SubdivisionParams,
) -> Result<(ParcelSet, SubdivisionStats), SubdivisionError> {
    params.validate()?;
    let mut stats = SubdivisionStats::default();
    let total_start = Instant::now();

    let topo_start = Instant::now();
    let mut working: RoadGraph;
    let g_ref = if graph.topology_valid {
        graph
    } else {
        working = graph.clone();
        working.rebuild_topology()?;
        &working
    };
    stats.topology_rebuild = topo_start.elapsed();

    let blocks_start = Instant::now();
    let blocks = extract_blocks(g_ref);
    stats.block_extraction = blocks_start.elapsed();
    stats.block_count = blocks.len();

    let mut parcels = ParcelSet::default();
    for (block_idx, block) in blocks.iter().enumerate() {
        let block_start = Instant::now();
        let new = subdivide_block(g_ref, block, params, block_idx)?;
        stats.block_subdivide_total += block_start.elapsed();
        for parcel in new {
            parcels.insert(parcel);
        }
    }
    stats.parcel_count = parcels.len();
    stats.total = total_start.elapsed();
    Ok((parcels, stats))
}

/// Subdivide a single block. Internal entry used both by
/// `subdivide_all` and by `apply_road_edit`'s regenerate path.
pub(crate) fn subdivide_block(
    graph: &RoadGraph,
    block: &Block,
    params: &SubdivisionParams,
    _block_idx: usize,
) -> Result<Vec<Parcel>, SubdivisionError> {
    let n = block.polygon.len();
    if n < 3 {
        return Ok(Vec::new());
    }
    let verts: Vec<DVec2> = block.polygon.vertices().to_vec();
    let face = block.face_id();

    // Edge lengths (parallel to verts).
    let edge_lens: Vec<f64> = (0..n)
        .map(|i| (verts[(i + 1) % n] - verts[i]).length())
        .collect();

    // Real-corner classification per vertex (D11). Acute corners
    // (interior angle < 60°) are skipped here — sliver-merge is
    // milestone 0.3 — so parcels there fall back to whatever the
    // regular walk produces (usually a malformed shape that triggers
    // I1 rejection during validation, but never a panic).
    let real_corner: Vec<bool> = (0..n)
        .map(|i| {
            if !is_real_corner(graph, block, i, params) {
                return false;
            }
            let v_prev = block.polygon.vertices()[(i + n - 1) % n];
            let v_curr = block.polygon.vertices()[i];
            let v_next = block.polygon.vertices()[(i + 1) % n];
            let t_in = (v_prev - v_curr).normalize_or_zero();
            let t_out = (v_next - v_curr).normalize_or_zero();
            let t_arrive = -t_in;
            let cross = t_arrive.x * t_out.y - t_arrive.y * t_out.x;
            let dot = t_arrive.dot(t_out);
            let turn = cross.atan2(dot);
            let interior = std::f64::consts::PI - turn;
            interior > 60.0_f64.to_radians()
        })
        .collect();

    // Per-edge depth cap (block half-depth bound).
    let depth_caps = edge_depth_caps(&block.polygon);

    let r_target = params.frontage_width;

    // First pass: decide which adjacent road wins frontage at each
    // real corner. This drives both the corner geometry below and the
    // per-edge walk bounds.
    let mut frontage_on_next: Vec<bool> = vec![false; n];
    for i in 0..n {
        if !real_corner[i] {
            continue;
        }
        let prev_len = edge_lens[(i + n - 1) % n];
        let next_len = edge_lens[i];
        let prev_road = block.roads[(i + n - 1) % n];
        let next_road = block.roads[i];
        frontage_on_next[i] = if (next_len - prev_len).abs() < EPS_GEOM {
            // tie → deterministic by RoadId
            next_road.data().as_ffi() <= prev_road.data().as_ffi()
        } else {
            next_len > prev_len
        };
    }

    let mut out: Vec<Parcel> = Vec::new();

    // Build corner parcels.
    for i in 0..n {
        if !real_corner[i] {
            continue;
        }
        let v_curr = verts[i];
        let v_prev = verts[(i + n - 1) % n];
        let v_next = verts[(i + 1) % n];

        let prev_road = block.roads[(i + n - 1) % n];
        let next_road = block.roads[i];

        let t_in = (v_prev - v_curr).normalize_or_zero();
        let t_out = (v_next - v_curr).normalize_or_zero();
        if t_in.length_squared() < EPS_GEOM * EPS_GEOM
            || t_out.length_squared() < EPS_GEOM * EPS_GEOM
        {
            continue;
        }

        let prev_len = edge_lens[(i + n - 1) % n];
        let next_len = edge_lens[i];
        let r = r_target
            .min(prev_len * 0.5)
            .min(next_len * 0.5)
            .max(params.min_frontage * 0.5);
        if r < params.min_frontage * 0.5 {
            continue;
        }

        // Inward normals.
        let n_out = DVec2::new(-t_out.y, t_out.x);
        let n_in = DVec2::new(t_in.y, -t_in.x);

        // Two flavors of corner parcel, both with frontage = R on the
        // *winning* road and side2 = depth along the other:
        //
        //   A. frontage on NEXT road: v0 → v0+t_out·R is frontage,
        //      v0+t_in·depth → v0 is the side along the prev road.
        //   B. frontage on PREV road: v0+t_in·R → v0 is frontage,
        //      v0 → v0+t_out·depth is the side along the next road.
        let (frontage_road, polygon, frontage_idx, consume_next, consume_prev) =
            if frontage_on_next[i] {
                let perp_depth = params.depth.min(prev_len * 0.5);
                let v0 = v_curr;
                let v1 = v_curr + t_out * r;
                let v2 = v_curr + t_out * r + n_out * perp_depth;
                let v3 = v_curr + t_in * perp_depth;
                let Ok(p) = Polygon::new(vec![v0, v1, v2, v3]) else {
                    continue;
                };
                (next_road, p, 0_usize, r, perp_depth)
            } else {
                let perp_depth = params.depth.min(next_len * 0.5);
                let v0 = v_curr;
                let v1 = v_curr + t_out * perp_depth;
                let v2 = v_curr + t_out * perp_depth + n_out * r;
                let v3 = v_curr + t_in * r;
                let Ok(p) = Polygon::new(vec![v0, v1, v2, v3]) else {
                    continue;
                };
                // Frontage edge in CCW order is v3 → v0, i.e. edge 3.
                (prev_road, p, 3_usize, perp_depth, r)
            };
        // The variables `consume_next` and `consume_prev` propagate
        // out via the per-vertex helpers below.
        let _ = (consume_next, consume_prev, n_in);

        // Clip to the block boundary so the parcel can never extend
        // past it. Required for non-rectangular blocks where the
        // R×depth rectangle would otherwise cross another block edge.
        let Some(polygon) = clip_polygon_to_block(&polygon, &block.polygon) else {
            continue;
        };
        if polygon.area() < params.min_area {
            continue;
        }
        let edge_kinds = classify_edges(polygon.len(), frontage_idx);

        let mut parcel = Parcel {
            polygon,
            edge_kinds,
            frontage_road,
            frontage_edge_index: frontage_idx,
            block: face,
            building: None,
        };
        regularize_parcel(&mut parcel, params);
        out.push(parcel);
    }

    // For each edge, walk between the corner footprints (if any) and
    // emit regular frontage parcels.
    for i in 0..n {
        let p = verts[i];
        let q = verts[(i + 1) % n];
        let road = block.roads[i];
        let edge_vec = q - p;
        let edge_len = edge_vec.length();
        if edge_len < params.min_frontage {
            continue;
        }
        let edge_dir = edge_vec / edge_len;
        let inward = DVec2::new(-edge_dir.y, edge_dir.x);

        // Walk-bounds. The corner at vertex i (start of this edge)
        // consumes the corner's `r` along this edge (it's the corner's
        // *next* side). The corner at vertex i+1 (end of this edge)
        // consumes the corner's `perp_depth` along this edge (it's
        // *that* corner's *prev* side, i.e., its Side2 runs along
        // here for a depth-equivalent distance).
        let consume_at_start = if real_corner[i] {
            // Edge i is the *next* edge of corner at vertex i.
            // - frontage_on_next[i] true  → consume R (frontage).
            // - frontage_on_next[i] false → consume depth (side2).
            if frontage_on_next[i] {
                corner_consume_short(&edge_lens, i, r_target, params)
            } else {
                corner_consume_long(&edge_lens, i, params)
            }
        } else {
            0.0
        };
        let next_v = (i + 1) % n;
        let consume_at_end = if real_corner[next_v] {
            // Edge i is the *prev* edge of corner at vertex (i+1).
            // - frontage_on_next[next_v] true  → consume depth (side2 along prev).
            // - frontage_on_next[next_v] false → consume R (frontage on prev).
            if frontage_on_next[next_v] {
                corner_consume_long(&edge_lens, next_v, params)
            } else {
                corner_consume_short(&edge_lens, next_v, r_target, params)
            }
        } else {
            0.0
        };
        let t_start = consume_at_start;
        let t_end = edge_len - consume_at_end;
        if t_end - t_start < params.min_frontage {
            continue;
        }

        let mut rng = rng_for_road(params.seed, road);
        let max_depth = depth_caps[i].min(params.depth + params.depth_variance.abs() + EPS_GEOM);
        let setback = params.setback.max(0.0);

        // Place split positions between t_start and t_end.
        let splits = split_positions_between(t_start, t_end, params, &mut rng);

        for k in 0..splits.len() - 1 {
            let t_a = splits[k];
            let t_b = splits[k + 1];
            if t_b - t_a < params.min_frontage - EPS_GEOM {
                continue;
            }
            let p_a = p + edge_dir * t_a;
            let p_b = p + edge_dir * t_b;
            let depth_jitter = params.depth_variance * (rng.gen::<f64>() * 2.0 - 1.0);
            let raw_depth = (params.depth + depth_jitter).max(params.min_frontage);
            let total_depth = (setback + raw_depth)
                .min(max_depth)
                .max(params.min_frontage * 0.5);
            if total_depth < params.min_frontage * 0.25 {
                continue;
            }
            let q_b = p_b + inward * total_depth;
            let q_a = p_a + inward * total_depth;

            let Ok(raw_poly) = Polygon::new(vec![p_a, p_b, q_b, q_a]) else {
                continue;
            };
            // Clip to block boundary so parcels can never escape the
            // block (matters at acute corners and on non-rectangular
            // blocks).
            let Some(polygon) = clip_polygon_to_block(&raw_poly, &block.polygon) else {
                continue;
            };
            if polygon.area() < params.min_area {
                continue;
            }
            // After clipping, find the frontage edge — it's the one
            // that lies on the road segment within tolerance.
            let frontage_idx = match find_frontage_edge_after_clip(&polygon, p, q) {
                Some(idx) => idx,
                None => continue,
            };
            let frontage_len = {
                let v = polygon.vertices();
                (v[(frontage_idx + 1) % v.len()] - v[frontage_idx]).length()
            };
            if frontage_len < params.min_frontage {
                continue;
            }
            let edge_kinds = classify_edges(polygon.len(), frontage_idx);
            let mut parcel = Parcel {
                polygon,
                edge_kinds,
                frontage_road: road,
                frontage_edge_index: frontage_idx,
                block: face,
                building: None,
            };
            regularize_parcel(&mut parcel, params);
            out.push(parcel);
        }
    }

    Ok(out)
}

/// True iff vertex `i` of the block is a "real corner" — graph degree
/// ≥3, or degree 2 with a bend sharper than 150°. (D11.)
fn is_real_corner(graph: &RoadGraph, block: &Block, i: usize, _params: &SubdivisionParams) -> bool {
    let n = block.polygon.len();
    let half_edge_id = block.boundary_edges[i];
    let Some(he) = graph.half_edges.get(half_edge_id) else {
        return false;
    };
    let node_id = he.origin;
    let Some(node) = graph.nodes.get(node_id) else {
        return false;
    };
    if node.outgoing.len() >= 3 {
        return true;
    }
    // Degree-2: bend test.
    let v_prev = block.polygon.vertices()[(i + n - 1) % n];
    let v_curr = block.polygon.vertices()[i];
    let v_next = block.polygon.vertices()[(i + 1) % n];
    let in_dir = (v_curr - v_prev).normalize_or_zero();
    let out_dir = (v_next - v_curr).normalize_or_zero();
    if in_dir.length_squared() < EPS_GEOM * EPS_GEOM
        || out_dir.length_squared() < EPS_GEOM * EPS_GEOM
    {
        return false;
    }
    let cos_turn = in_dir.dot(out_dir).clamp(-1.0, 1.0);
    let turn_angle = cos_turn.acos();
    // turn > 30° means interior < 150° = "real corner".
    turn_angle > 30.0_f64.to_radians()
}

/// Iteratively clip `parcel` by each inward half-plane of `block`.
/// Returns `None` if the result is empty or invalid. Correct for
/// convex blocks; may over-clip mild concavities (acceptable for
/// milestone 0.2).
fn clip_polygon_to_block(parcel: &Polygon, block: &Polygon) -> Option<Polygon> {
    let mut current = parcel.clone();
    let block_verts = block.vertices();
    let n = block_verts.len();
    for i in 0..n {
        let a = block_verts[i];
        let b = block_verts[(i + 1) % n];
        let edge = b - a;
        let len = edge.length();
        if len < EPS_GEOM {
            continue;
        }
        let dir = edge / len;
        let inward = DVec2::new(-dir.y, dir.x);
        current = current.clip_half_plane(a, inward)?;
    }
    Some(current)
}

/// Locate the polygon edge that lies on the road segment
/// `(road_a, road_b)` within tolerance.
fn find_frontage_edge_after_clip(polygon: &Polygon, road_a: DVec2, road_b: DVec2) -> Option<usize> {
    let road_dir = (road_b - road_a).normalize_or_zero();
    if road_dir.length_squared() < EPS_GEOM * EPS_GEOM {
        return None;
    }
    let road_normal = DVec2::new(-road_dir.y, road_dir.x);
    let v = polygon.vertices();
    let n = v.len();
    let tol = 100.0 * EPS_GEOM;
    let mut best: Option<(usize, f64)> = None;
    for i in 0..n {
        let p = v[i];
        let q = v[(i + 1) % n];
        let dp = (p - road_a).dot(road_normal).abs();
        let dq = (q - road_a).dot(road_normal).abs();
        if dp < tol && dq < tol {
            let edge_dir = (q - p).normalize_or_zero();
            if edge_dir.dot(road_dir).abs() > 0.5 {
                let len = (q - p).length();
                if best.is_none_or(|(_, blen)| len > blen) {
                    best = Some((i, len));
                }
            }
        }
    }
    best.map(|(i, _)| i)
}

/// "Short" corner consumption — equal to the corner radius `R`,
/// capped at half each adjacent edge.
fn corner_consume_short(
    edge_lens: &[f64],
    i: usize,
    r_target: f64,
    params: &SubdivisionParams,
) -> f64 {
    let n = edge_lens.len();
    let prev_len = edge_lens[(i + n - 1) % n];
    let next_len = edge_lens[i];
    r_target
        .min(prev_len * 0.5)
        .min(next_len * 0.5)
        .max(params.min_frontage * 0.5)
}

/// "Long" corner consumption — equal to the parcel depth, capped at
/// half each adjacent edge.
fn corner_consume_long(edge_lens: &[f64], i: usize, params: &SubdivisionParams) -> f64 {
    let n = edge_lens.len();
    let prev_len = edge_lens[(i + n - 1) % n];
    let next_len = edge_lens[i];
    params
        .depth
        .min(prev_len * 0.5)
        .min(next_len * 0.5)
        .max(params.min_frontage * 0.5)
}

/// Per-edge depth cap (half the perpendicular distance from edge
/// midpoint to the next polygon edge along the inward normal).
fn edge_depth_caps(poly: &Polygon) -> Vec<f64> {
    let verts = poly.vertices();
    let n = verts.len();
    let mut caps = Vec::with_capacity(n);
    for i in 0..n {
        let p = verts[i];
        let q = verts[(i + 1) % n];
        let edge = q - p;
        let len = edge.length();
        if len < EPS_GEOM {
            caps.push(0.0);
            continue;
        }
        let dir = edge / len;
        let inward = DVec2::new(-dir.y, dir.x);
        let mid = (p + q) * 0.5;
        let mut min_t = f64::INFINITY;
        for j in 0..n {
            if j == i {
                continue;
            }
            let a = verts[j];
            let b = verts[(j + 1) % n];
            if let Some(t) = ray_segment_distance(mid, inward, a, b) {
                if t > EPS_GEOM && t < min_t {
                    min_t = t;
                }
            }
        }
        caps.push(min_t * 0.5);
    }
    caps
}

fn ray_segment_distance(origin: DVec2, dir: DVec2, a: DVec2, b: DVec2) -> Option<f64> {
    let v1 = origin - a;
    let v2 = b - a;
    let v3 = DVec2::new(-dir.y, dir.x);
    let denom = v2.dot(v3);
    if denom.abs() < EPS_GEOM {
        return None;
    }
    let t = v2.perp_dot(v1) / denom;
    let s = v1.dot(v3) / denom;
    if t >= 0.0 && (-EPS_GEOM..=1.0 + EPS_GEOM).contains(&s) {
        Some(t)
    } else {
        None
    }
}

/// Split positions between [`start`, `end`] along an edge. Begins
/// with `start` and ends with `end`.
fn split_positions_between(
    start: f64,
    end: f64,
    params: &SubdivisionParams,
    rng: &mut ChaCha8Rng,
) -> Vec<f64> {
    let mut splits: Vec<f64> = vec![start];
    let mut t = start;
    loop {
        let jitter = params.frontage_variance * (rng.gen::<f64>() * 2.0 - 1.0);
        let dt = (params.frontage_width + jitter).max(params.min_frontage);
        t += dt;
        if t >= end - params.min_frontage * 0.5 {
            break;
        }
        splits.push(t);
    }
    splits.push(end);
    splits
}

fn rng_for_road(seed: u64, road: RoadId) -> ChaCha8Rng {
    let road_bits = road.data().as_ffi();
    let mixed = seed
        .wrapping_mul(0x9E37_79B9_7F4A_7C15)
        .wrapping_add(road_bits)
        .wrapping_mul(0xBF58_476D_1CE4_E5B9);
    ChaCha8Rng::seed_from_u64(mixed)
}

#[allow(dead_code)]
fn _unused_node_id() -> NodeId {
    NodeId::null()
}

#[allow(dead_code)]
fn _unused_face_id() -> FaceId {
    FaceId::null()
}

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

    fn rectangle_graph(w: f64, h: f64) -> RoadGraph {
        let mut g = RoadGraph::new();
        let a = g.add_node(DVec2::new(0.0, 0.0));
        let b = g.add_node(DVec2::new(w, 0.0));
        let c = g.add_node(DVec2::new(w, h));
        let d = g.add_node(DVec2::new(0.0, h));
        g.add_road(a, b).unwrap();
        g.add_road(b, c).unwrap();
        g.add_road(c, d).unwrap();
        g.add_road(d, a).unwrap();
        g.rebuild_topology().unwrap();
        g
    }

    #[test]
    fn rectangle_subdivides_into_many_parcels() {
        let g = rectangle_graph(200.0, 100.0);
        let params = SubdivisionParams::default();
        let set = subdivide_all(&g, &params).unwrap();
        assert!(set.len() > 8, "got {} parcels", set.len());
    }

    #[test]
    fn rectangle_parcels_have_one_frontage_each() {
        let g = rectangle_graph(200.0, 100.0);
        let params = SubdivisionParams::default();
        let set = subdivide_all(&g, &params).unwrap();
        for (_, p) in set.iter() {
            let frontage_count = p
                .edge_kinds()
                .iter()
                .filter(|k| matches!(k, super::super::EdgeKind::Frontage))
                .count();
            assert_eq!(frontage_count, 1, "I2 violation");
        }
    }

    #[test]
    fn rectangle_parcels_do_not_overlap() {
        let g = rectangle_graph(200.0, 100.0);
        let params = SubdivisionParams::default();
        let set = subdivide_all(&g, &params).unwrap();
        let total_area: f64 = set.iter().map(|(_, p)| p.area()).sum();
        let block_area = 200.0 * 100.0;
        assert!(
            total_area <= block_area + 1e-3,
            "parcel areas {total_area} > block area {block_area}",
        );
    }

    #[test]
    fn deterministic_with_fixed_seed() {
        let g = rectangle_graph(200.0, 100.0);
        let params = SubdivisionParams::default();
        let a = subdivide_all(&g, &params).unwrap();
        let b = subdivide_all(&g, &params).unwrap();
        assert_eq!(a.len(), b.len());
        for ((_, p1), (_, p2)) in a.iter().zip(b.iter()) {
            assert_eq!(p1.vertices().len(), p2.vertices().len());
            for (v1, v2) in p1.vertices().iter().zip(p2.vertices().iter()) {
                assert!((v1.x - v2.x).abs() < 1e-12);
                assert!((v1.y - v2.y).abs() < 1e-12);
            }
        }
    }

    #[test]
    fn invalid_params_rejected() {
        let g = rectangle_graph(50.0, 50.0);
        let params = SubdivisionParams {
            frontage_width: -1.0,
            ..SubdivisionParams::default()
        };
        let err = subdivide_all(&g, &params).unwrap_err();
        assert!(matches!(err, SubdivisionError::InvalidParams(_)));
    }

    #[test]
    fn rectangle_has_corner_parcels() {
        let g = rectangle_graph(200.0, 100.0);
        let params = SubdivisionParams::default();
        let set = subdivide_all(&g, &params).unwrap();
        // Parcel at the (0, 0) corner: its polygon should contain the
        // origin (or a vertex right at it) and have 4 vertices.
        let near_origin = set
            .iter()
            .find(|(_, p)| p.vertices().iter().any(|v| v.length_squared() < 1.0))
            .map(|(_, p)| p.vertices().len());
        assert_eq!(
            near_origin,
            Some(4),
            "expected a 4-vertex corner parcel near the origin"
        );
    }

    #[test]
    fn stats_report_nonzero_phases() {
        let g = rectangle_graph(200.0, 100.0);
        let params = SubdivisionParams::default();
        let (set, stats) = subdivide_all_with_stats(&g, &params).unwrap();
        assert_eq!(stats.parcel_count, set.len());
        assert!(stats.block_count >= 1);
        assert!(stats.total > Duration::ZERO);
        assert!(stats.parcels_per_second() > 0.0);
    }
}