eunoia 0.16.0

//! Axis-aligned square shape implementation.
//!
//! Squares are parameterised by `[x, y, side]` (3 parameters per shape, the
//! same count as [`Circle`]). They are kept axis-aligned: rotation is not part
//! of the parameter vector. The 2D MDS phase of the fitter still works on
//! `(x, y)` pairs, with the canonical-direction overlap inversion provided by
//! [`Square::mds_target_distance`].
//!
//! The n-way intersection of axis-aligned squares is itself one axis-aligned
//! rectangle (the intersection of axis-aligned rectangles is axis-aligned),
//! so [`Square::compute_exclusive_regions`] is exact in closed form with no
//! polygon-clipping library required.
//!
//! [`Circle`]: crate::geometry::shapes::Circle

use std::collections::HashMap;
use std::f64::consts::PI;

use crate::geometry::diagram::{
    IntersectionPoint, RegionMask, discover_regions, mask_to_indices, to_exclusive_areas,
    to_exclusive_areas_and_gradients,
};
use crate::geometry::primitives::Point;
use crate::geometry::shapes::{Polygon, Rectangle};
use crate::geometry::traits::{
    Area, BoundingBox, Centroid, Closed, DiagramShape, Distance, ExclusiveRegionsAndGradient,
    Perimeter, Polygonize,
};

/// An axis-aligned square defined by a center point and side length.
///
/// # Examples
///
/// ```
/// use eunoia::geometry::shapes::Square;
/// use eunoia::geometry::traits::{Area, Closed};
/// use eunoia::geometry::primitives::Point;
///
/// let s1 = Square::new(Point::new(0.0, 0.0), 2.0);
/// let s2 = Square::new(Point::new(1.0, 0.0), 2.0);
///
/// let area = s1.area();
/// let overlap = s1.intersection_area(&s2);
/// ```
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct Square {
    center: Point,
    side: f64,
}

impl Square {
    /// Creates a new axis-aligned square with the given center and side length.
    ///
    /// # Panics
    ///
    /// Panics if `side <= 0`. Use [`Square::try_new`] to handle invalid
    /// input as a [`crate::error::DiagramError`] instead of a panic —
    /// bindings authors writing FFI wrappers should reach for `try_new`.
    pub fn new(center: Point, side: f64) -> Self {
        assert!(side > 0.0, "Square side must be > 0, got {}", side);
        Square { center, side }
    }

    /// Fallible constructor: returns
    /// [`crate::error::DiagramError::InvalidShapeParameter`] when `side <= 0`
    /// instead of panicking. Use this when constructing squares from
    /// untrusted input (e.g. across an FFI boundary).
    ///
    /// # Examples
    ///
    /// ```
    /// use eunoia::geometry::shapes::Square;
    /// use eunoia::geometry::primitives::Point;
    ///
    /// assert!(Square::try_new(Point::new(0.0, 0.0), 1.0).is_ok());
    /// assert!(Square::try_new(Point::new(0.0, 0.0), 0.0).is_err());
    /// assert!(Square::try_new(Point::new(0.0, 0.0), -2.0).is_err());
    /// ```
    pub fn try_new(center: Point, side: f64) -> Result<Self, crate::error::DiagramError> {
        if side > 0.0 {
            Ok(Square { center, side })
        } else {
            Err(crate::error::DiagramError::InvalidShapeParameter {
                shape: "Square",
                param: "side",
                value: side,
            })
        }
    }

    /// Returns the center point of the square.
    pub fn center(&self) -> Point {
        self.center
    }

    /// Returns the side length of the square.
    pub fn side(&self) -> f64 {
        self.side
    }

    /// Returns `(x_min, x_max, y_min, y_max)` for the square's axis-aligned bounds.
    pub fn bounds(&self) -> (f64, f64, f64, f64) {
        let h = self.side * 0.5;
        (
            self.center.x() - h,
            self.center.x() + h,
            self.center.y() - h,
            self.center.y() + h,
        )
    }

    /// View the square as the equivalent [`Rectangle`]. Used internally to
    /// delegate [`Closed`] operations.
    fn as_rectangle(&self) -> Rectangle {
        Rectangle::new(self.center, self.side, self.side)
    }
}

impl Area for Square {
    fn area(&self) -> f64 {
        self.side * self.side
    }
}

impl Centroid for Square {
    fn centroid(&self) -> Point {
        self.center
    }
}

impl Perimeter for Square {
    fn perimeter(&self) -> f64 {
        4.0 * self.side
    }
}

impl BoundingBox for Square {
    fn bounding_box(&self) -> Rectangle {
        self.as_rectangle()
    }
}

impl Distance for Square {
    fn distance(&self, other: &Self) -> f64 {
        self.as_rectangle().distance(&other.as_rectangle())
    }
}

impl Closed for Square {
    fn contains(&self, other: &Self) -> bool {
        self.as_rectangle().contains(&other.as_rectangle())
    }

    fn contains_point(&self, point: &Point) -> bool {
        self.as_rectangle().contains_point(point)
    }

    fn intersects(&self, other: &Self) -> bool {
        self.as_rectangle().intersects(&other.as_rectangle())
    }

    fn intersection_area(&self, other: &Self) -> f64 {
        let (ax0, ax1, ay0, ay1) = self.bounds();
        let (bx0, bx1, by0, by1) = other.bounds();
        let dx = (ax1.min(bx1) - ax0.max(bx0)).max(0.0);
        let dy = (ay1.min(by1) - ay0.max(by0)).max(0.0);
        dx * dy
    }

    /// Edge crossings between two axis-aligned squares.
    ///
    /// Each crossing is a horizontal edge of one square meeting a vertical
    /// edge of the other inside both squares. Returns 0–2 crossings on the
    /// shared boundary; coincident edges produce zero crossings (overlap is
    /// still witnessed by [`Closed::intersection_area`]).
    fn intersection_points(&self, other: &Self) -> Vec<Point> {
        let (ax0, ax1, ay0, ay1) = self.bounds();
        let (bx0, bx1, by0, by1) = other.bounds();
        let mut points = Vec::new();

        // Horizontal edges of `self` at y = ay0, ay1, x ∈ [ax0, ax1] cross
        // vertical edges of `other` at x = bx0, bx1, y ∈ [by0, by1].
        for &y in &[ay0, ay1] {
            for &x in &[bx0, bx1] {
                if x >= ax0 && x <= ax1 && y >= by0 && y <= by1 {
                    points.push(Point::new(x, y));
                }
            }
        }
        // And the symmetric case: horizontal edges of `other` × vertical
        // edges of `self`.
        for &y in &[by0, by1] {
            for &x in &[ax0, ax1] {
                if x >= bx0 && x <= bx1 && y >= ay0 && y <= ay1 {
                    points.push(Point::new(x, y));
                }
            }
        }

        // De-duplicate — corner-touching configurations register the same
        // point from both directions.
        points.sort_by(|p, q| {
            p.x()
                .partial_cmp(&q.x())
                .unwrap_or(std::cmp::Ordering::Equal)
                .then_with(|| {
                    p.y()
                        .partial_cmp(&q.y())
                        .unwrap_or(std::cmp::Ordering::Equal)
                })
        });
        points.dedup_by(|a, b| (a.x() - b.x()).abs() < 1e-12 && (a.y() - b.y()).abs() < 1e-12);
        points
    }
}

/// Collect pairwise edge-crossings between axis-aligned squares for the
/// region-discovery pass. Mirrors [`crate::geometry::diagram::collect_intersections`]
/// (circles) and [`crate::geometry::shapes::ellipse::collect_intersections_ellipse`]
/// (ellipses) — same `IntersectionPoint` shape, with adopters populated by
/// scanning every other shape that contains the point.
fn collect_intersections_square(squares: &[Square], n_sets: usize) -> Vec<IntersectionPoint> {
    let mut intersections = Vec::new();
    for i in 0..n_sets {
        for j in (i + 1)..n_sets {
            let pts = squares[i].intersection_points(&squares[j]);
            for point in pts {
                let mut adopters = vec![i, j];
                for (k, sq) in squares.iter().enumerate().take(n_sets) {
                    if k != i && k != j && sq.contains_point(&point) {
                        adopters.push(k);
                    }
                }
                adopters.sort_unstable();
                intersections.push(IntersectionPoint::new(point, (i, j), adopters));
            }
        }
    }
    intersections
}

/// Companion to [`Square::compute_exclusive_regions`] that also returns the
/// analytical gradient of each exclusive area w.r.t. the flat parameter
/// vector `[x₀, y₀, s₀, x₁, …]`.
///
/// For each region the n-way intersection is one axis-aligned rectangle
/// with `dx = x_max − x_min` and `dy = y_max − y_min`; each of the four
/// extrema is a min/max over the mask. When a single shape strictly binds a
/// side, that side's contributions go to that shape:
///
/// ```text
/// ∂A/∂x_{i_L} = −dy   ∂A/∂x_{i_R} = +dy
/// ∂A/∂y_{i_B} = −dx   ∂A/∂y_{i_T} = +dx
/// ∂A/∂s_{i_L} = dy/2  ∂A/∂s_{i_R} = dy/2
/// ∂A/∂s_{i_B} = dx/2  ∂A/∂s_{i_T} = dx/2
/// ```
///
/// At a tie (multiple shapes share the extremum on a side, e.g. coincident
/// edges), the side's contribution is split equally among the tied shapes.
/// This matches the central-difference average of one-sided derivatives at
/// the non-smooth point and produces a valid subgradient. Contributions sum
/// when the same shape binds multiple sides (e.g. the 1-mask case yields
/// `∂A/∂s_i = 2s`, matching `A = s²`). When `dx ≤ 0` or `dy ≤ 0` the area
/// is clamped to 0 and the gradient for that region is the zero vector. The
/// shared inclusion-exclusion combiner [`to_exclusive_areas_and_gradients`]
/// already zeroes gradients for post-IE clamped-negative areas.
fn compute_exclusive_regions_with_gradient_squares(
    shapes: &[Square],
) -> ExclusiveRegionsAndGradient {
    let n_sets = shapes.len();
    let n_params = n_sets * 3;

    let intersections = collect_intersections_square(shapes, n_sets);
    let regions = discover_regions(shapes, &intersections, n_sets);

    let mut overlapping_areas: HashMap<RegionMask, f64> = HashMap::new();
    let mut overlapping_grads: HashMap<RegionMask, Vec<f64>> = HashMap::new();

    for &mask in &regions {
        let indices = mask_to_indices(mask, n_sets);

        let mut x_min = f64::NEG_INFINITY;
        let mut x_max = f64::INFINITY;
        let mut y_min = f64::NEG_INFINITY;
        let mut y_max = f64::INFINITY;
        for &i in &indices {
            let (a, b, c, d) = shapes[i].bounds();
            if a > x_min {
                x_min = a;
            }
            if b < x_max {
                x_max = b;
            }
            if c > y_min {
                y_min = c;
            }
            if d < y_max {
                y_max = d;
            }
        }

        let dx_raw = x_max - x_min;
        let dy_raw = y_max - y_min;
        let dx = dx_raw.max(0.0);
        let dy = dy_raw.max(0.0);
        overlapping_areas.insert(mask, dx * dy);

        let mut grad = vec![0.0; n_params];
        if dx_raw > 0.0 && dy_raw > 0.0 {
            // Collect tied binding indices per side. Bounds and the running
            // extremum are produced by arithmetic on input parameters, so
            // structural edge coincidences hit exact f64 equality.
            let mut tied_l: Vec<usize> = Vec::with_capacity(indices.len());
            let mut tied_r: Vec<usize> = Vec::with_capacity(indices.len());
            let mut tied_b: Vec<usize> = Vec::with_capacity(indices.len());
            let mut tied_t: Vec<usize> = Vec::with_capacity(indices.len());
            for &i in &indices {
                let (a, b, c, d) = shapes[i].bounds();
                #[allow(clippy::float_cmp)]
                {
                    if a == x_min {
                        tied_l.push(i);
                    }
                    if b == x_max {
                        tied_r.push(i);
                    }
                    if c == y_min {
                        tied_b.push(i);
                    }
                    if d == y_max {
                        tied_t.push(i);
                    }
                }
            }
            let w_l = 1.0 / tied_l.len() as f64;
            let w_r = 1.0 / tied_r.len() as f64;
            let w_b = 1.0 / tied_b.len() as f64;
            let w_t = 1.0 / tied_t.len() as f64;

            for &i in &tied_l {
                grad[3 * i] -= dy * w_l;
                grad[3 * i + 2] += dy * 0.5 * w_l;
            }
            for &i in &tied_r {
                grad[3 * i] += dy * w_r;
                grad[3 * i + 2] += dy * 0.5 * w_r;
            }
            for &i in &tied_b {
                grad[3 * i + 1] -= dx * w_b;
                grad[3 * i + 2] += dx * 0.5 * w_b;
            }
            for &i in &tied_t {
                grad[3 * i + 1] += dx * w_t;
                grad[3 * i + 2] += dx * 0.5 * w_t;
            }
        }
        overlapping_grads.insert(mask, grad);
    }

    to_exclusive_areas_and_gradients(&overlapping_areas, &overlapping_grads, n_params)
}

/// Clipped exclusive areas: each region's intersection with `container` is
/// itself an axis-aligned rectangle (axis-aligned squares + an axis-aligned
/// container rectangle still intersect to one axis-aligned rectangle). Mask
/// `0` is seeded with the container area; inclusion-exclusion produces
/// `complement = container.area − area(⋃ squares ∩ container)`.
pub(crate) fn compute_exclusive_regions_clipped_squares(
    shapes: &[Square],
    container: &Rectangle,
) -> HashMap<RegionMask, f64> {
    let n_sets = shapes.len();
    let intersections = collect_intersections_square(shapes, n_sets);
    let regions = discover_regions(shapes, &intersections, n_sets);

    let (cx_min, cx_max, cy_min, cy_max) = container.bounds();
    let mut overlapping_areas: HashMap<RegionMask, f64> = HashMap::new();
    overlapping_areas.insert(0, container.area());

    for &mask in &regions {
        let indices = mask_to_indices(mask, n_sets);
        let mut x_min = cx_min;
        let mut x_max = cx_max;
        let mut y_min = cy_min;
        let mut y_max = cy_max;
        for &i in &indices {
            let (a, b, c, d) = shapes[i].bounds();
            if a > x_min {
                x_min = a;
            }
            if b < x_max {
                x_max = b;
            }
            if c > y_min {
                y_min = c;
            }
            if d < y_max {
                y_max = d;
            }
        }
        let dx = (x_max - x_min).max(0.0);
        let dy = (y_max - y_min).max(0.0);
        overlapping_areas.insert(mask, dx * dy);
    }

    to_exclusive_areas(&overlapping_areas)
}

/// Gradient companion to [`compute_exclusive_regions_clipped_squares`].
///
/// Returns `(exclusive_areas, exclusive_grads)` where each gradient vector has
/// length `n_sets · 3 + 4` and is laid out as
/// `[x₀, y₀, s₀, …, x_c, y_c, u_c, v_c]` — the trailing four entries are the
/// container's optimizer encoding (`u = ln(w·h)`, `v = ln(w/h)`).
///
/// Per-region geometry mirrors [`compute_exclusive_regions_with_gradient_squares`]
/// extended with the container as one more potential binder per side. Tied
/// shapes plus the container split each side's contribution equally. Container
/// edges chain through `[x_c, y_c, w_c, h_c] → [x_c, y_c, u_c, v_c]` exactly
/// as in [`compute_exclusive_regions_clipped_with_gradient_rectangles`]. Mask
/// `0` is seeded with `container.area()` and a gradient
/// `[…, 0, 0, w·h, 0]`.
pub(crate) fn compute_exclusive_regions_clipped_with_gradient_squares(
    shapes: &[Square],
    container: &Rectangle,
) -> ExclusiveRegionsAndGradient {
    let n_sets = shapes.len();
    let n_params_per_shape = 3;
    let container_off = n_sets * n_params_per_shape;
    let n_params = container_off + 4;

    let intersections = collect_intersections_square(shapes, n_sets);
    let regions = discover_regions(shapes, &intersections, n_sets);

    let (cx_min, cx_max, cy_min, cy_max) = container.bounds();
    let container_area = container.area();

    let mut overlapping_areas: HashMap<RegionMask, f64> = HashMap::new();
    let mut overlapping_grads: HashMap<RegionMask, Vec<f64>> = HashMap::new();

    overlapping_areas.insert(0, container_area);
    let mut zero_grad = vec![0.0; n_params];
    zero_grad[container_off + 2] = container_area;
    overlapping_grads.insert(0, zero_grad);

    for &mask in &regions {
        let indices = mask_to_indices(mask, n_sets);

        let mut x_min = cx_min;
        let mut x_max = cx_max;
        let mut y_min = cy_min;
        let mut y_max = cy_max;
        for &i in &indices {
            let (a, b, c, d) = shapes[i].bounds();
            if a > x_min {
                x_min = a;
            }
            if b < x_max {
                x_max = b;
            }
            if c > y_min {
                y_min = c;
            }
            if d < y_max {
                y_max = d;
            }
        }

        let dx_raw = x_max - x_min;
        let dy_raw = y_max - y_min;
        let dx = dx_raw.max(0.0);
        let dy = dy_raw.max(0.0);
        overlapping_areas.insert(mask, dx * dy);

        let mut grad = vec![0.0; n_params];
        if dx_raw > 0.0 && dy_raw > 0.0 {
            let mut tied_l: Vec<usize> = Vec::with_capacity(indices.len());
            let mut tied_r: Vec<usize> = Vec::with_capacity(indices.len());
            let mut tied_b: Vec<usize> = Vec::with_capacity(indices.len());
            let mut tied_t: Vec<usize> = Vec::with_capacity(indices.len());
            for &i in &indices {
                let (a, b, c, d) = shapes[i].bounds();
                #[allow(clippy::float_cmp)]
                {
                    if a == x_min {
                        tied_l.push(i);
                    }
                    if b == x_max {
                        tied_r.push(i);
                    }
                    if c == y_min {
                        tied_b.push(i);
                    }
                    if d == y_max {
                        tied_t.push(i);
                    }
                }
            }
            #[allow(clippy::float_cmp)]
            let tied_l_c = cx_min == x_min;
            #[allow(clippy::float_cmp)]
            let tied_r_c = cx_max == x_max;
            #[allow(clippy::float_cmp)]
            let tied_b_c = cy_min == y_min;
            #[allow(clippy::float_cmp)]
            let tied_t_c = cy_max == y_max;

            let n_l = tied_l.len() + tied_l_c as usize;
            let n_r = tied_r.len() + tied_r_c as usize;
            let n_b = tied_b.len() + tied_b_c as usize;
            let n_t = tied_t.len() + tied_t_c as usize;

            let w_l = 1.0 / n_l as f64;
            let w_r = 1.0 / n_r as f64;
            let w_b = 1.0 / n_b as f64;
            let w_t = 1.0 / n_t as f64;

            // Square shape contributions: [x, y, side]; side derivative
            // collects from all four edges (each at half weight).
            for &i in &tied_l {
                grad[3 * i] -= dy * w_l;
                grad[3 * i + 2] += dy * 0.5 * w_l;
            }
            for &i in &tied_r {
                grad[3 * i] += dy * w_r;
                grad[3 * i + 2] += dy * 0.5 * w_r;
            }
            for &i in &tied_b {
                grad[3 * i + 1] -= dx * w_b;
                grad[3 * i + 2] += dx * 0.5 * w_b;
            }
            for &i in &tied_t {
                grad[3 * i + 1] += dx * w_t;
                grad[3 * i + 2] += dx * 0.5 * w_t;
            }

            // Container contributions: build geometric `[x, y, w, h]` then
            // chain into optimizer `[x, y, u, v]`.
            let mut da_dx_c = 0.0;
            let mut da_dy_c = 0.0;
            let mut da_dw_c = 0.0;
            let mut da_dh_c = 0.0;
            if tied_l_c {
                da_dx_c -= dy * w_l;
                da_dw_c += dy * 0.5 * w_l;
            }
            if tied_r_c {
                da_dx_c += dy * w_r;
                da_dw_c += dy * 0.5 * w_r;
            }
            if tied_b_c {
                da_dy_c -= dx * w_b;
                da_dh_c += dx * 0.5 * w_b;
            }
            if tied_t_c {
                da_dy_c += dx * w_t;
                da_dh_c += dx * 0.5 * w_t;
            }
            let w_c = container.width();
            let h_c = container.height();
            grad[container_off] = da_dx_c;
            grad[container_off + 1] = da_dy_c;
            grad[container_off + 2] = 0.5 * (da_dw_c * w_c + da_dh_c * h_c);
            grad[container_off + 3] = 0.5 * (da_dw_c * w_c - da_dh_c * h_c);
        }
        overlapping_grads.insert(mask, grad);
    }

    to_exclusive_areas_and_gradients(&overlapping_areas, &overlapping_grads, n_params)
}

impl DiagramShape for Square {
    fn compute_exclusive_regions(shapes: &[Self]) -> HashMap<RegionMask, f64> {
        let n_sets = shapes.len();
        let intersections = collect_intersections_square(shapes, n_sets);
        let regions = discover_regions(shapes, &intersections, n_sets);

        let mut overlapping_areas: HashMap<RegionMask, f64> = HashMap::new();
        for &mask in &regions {
            let indices = mask_to_indices(mask, n_sets);
            // The n-way intersection of axis-aligned squares is one
            // axis-aligned rectangle (or empty).
            let mut x_min = f64::NEG_INFINITY;
            let mut x_max = f64::INFINITY;
            let mut y_min = f64::NEG_INFINITY;
            let mut y_max = f64::INFINITY;
            for &i in &indices {
                let (a, b, c, d) = shapes[i].bounds();
                if a > x_min {
                    x_min = a;
                }
                if b < x_max {
                    x_max = b;
                }
                if c > y_min {
                    y_min = c;
                }
                if d < y_max {
                    y_max = d;
                }
            }
            let dx = (x_max - x_min).max(0.0);
            let dy = (y_max - y_min).max(0.0);
            overlapping_areas.insert(mask, dx * dy);
        }

        to_exclusive_areas(&overlapping_areas)
    }

    fn optimizer_params_from_circle(x: f64, y: f64, radius: f64) -> Vec<f64> {
        // Map the circle warm-start (area = π·r²) to a square of equal area
        // (`side = r·√π`). The optimiser refines from there; this just
        // minimises the global rescaling needed before overlap targets bite.
        vec![x, y, radius * PI.sqrt()]
    }

    fn mds_target_distance(
        area_i: f64,
        area_j: f64,
        target_overlap: f64,
    ) -> Result<f64, crate::error::DiagramError> {
        let s_i = area_i.sqrt();
        let s_j = area_j.sqrt();
        let half_sum = 0.5 * (s_i + s_j);

        // Inversion direction: |dx| = |dy|. With axis-aligned squares,
        //   overlap = ((s_i + s_j)/2 − |dx|) · ((s_i + s_j)/2 − |dy|),
        // so under |dx| = |dy| = d/√2 we get
        //   overlap = (half_sum − d/√2)²
        //         d = √2 · (half_sum − √overlap), clipped to non-negative.
        if target_overlap <= 0.0 {
            // Disjoint — use the diagonal-aligned tangency distance, which
            // is the smallest center distance for which the squares stop
            // touching when arranged at 45°.
            return Ok(half_sum * 2.0_f64.sqrt());
        }
        let root = target_overlap.sqrt();
        if root > half_sum {
            // Asked for more overlap than the smaller square can possibly
            // contribute (would require negative distance). Treat as full
            // containment along the canonical direction.
            return Ok(0.0);
        }
        let d = 2.0_f64.sqrt() * (half_sum - root);
        Ok(d.max(0.0))
    }

    fn n_params() -> usize {
        3 // x, y, side
    }

    fn from_params(params: &[f64]) -> Self {
        debug_assert_eq!(params.len(), 3, "Square requires 3 parameters: x, y, side");
        Square::new(
            Point::new(params[0], params[1]),
            params[2].max(f64::MIN_POSITIVE),
        )
    }

    fn to_params(&self) -> Vec<f64> {
        vec![self.center.x(), self.center.y(), self.side]
    }

    fn compute_exclusive_regions_with_gradient(
        shapes: &[Self],
    ) -> Option<ExclusiveRegionsAndGradient> {
        Some(compute_exclusive_regions_with_gradient_squares(shapes))
    }

    fn compute_exclusive_regions_clipped(
        shapes: &[Self],
        container: &Rectangle,
    ) -> Option<HashMap<RegionMask, f64>> {
        Some(compute_exclusive_regions_clipped_squares(shapes, container))
    }

    fn compute_exclusive_regions_clipped_with_gradient(
        shapes: &[Self],
        container: &Rectangle,
    ) -> Option<ExclusiveRegionsAndGradient> {
        Some(compute_exclusive_regions_clipped_with_gradient_squares(
            shapes, container,
        ))
    }

    /// Canonical axis-aligned Venn arrangements for `n ∈ {1, 2, 3}`.
    ///
    /// `n ≥ 4` returns `None`: there is no axis-aligned-square arrangement
    /// that opens all `2ⁿ − 1` regions for n ≥ 4. The footprint roughly
    /// matches `Ellipse::canonical_venn_layout` (radius ~1) so callers
    /// rescaling by a spec's mean circle radius land in the right magnitude.
    fn canonical_venn_layout(n: usize) -> Option<Vec<Self>> {
        // Layouts chosen so every one of the `2ⁿ − 1` non-empty regions has
        // strictly positive area; verified by the topology test in
        // `crate::venn` and by `test_canonical_venn_layout_topology` below.
        let centers_and_side: &[((f64, f64), f64)] = match n {
            1 => &[((0.0, 0.0), 2.0)],
            2 => &[((-0.4, 0.0), 1.0), ((0.4, 0.0), 1.0)],
            // n=3 footprint matches the existing N3 ellipse-Venn vertices
            // (equilateral-triangle-ish, circumradius ~0.55), with side=1.0
            // so all three pairwise overlaps and the central 3-way region
            // remain open.
            3 => &[
                ((0.0, 0.36), 1.0),
                ((0.42, -0.36), 1.0),
                ((-0.42, -0.36), 1.0),
            ],
            _ => return None,
        };
        Some(
            centers_and_side
                .iter()
                .map(|&((x, y), s)| Square::new(Point::new(x, y), s))
                .collect(),
        )
    }
}

impl Polygonize for Square {
    /// Returns the four corners as a CCW polygon. `n_vertices` is ignored —
    /// a square has exactly four vertices.
    fn polygonize(&self, _n_vertices: usize) -> Polygon {
        let (x_min, x_max, y_min, y_max) = self.bounds();
        Polygon::new(vec![
            Point::new(x_min, y_min),
            Point::new(x_max, y_min),
            Point::new(x_max, y_max),
            Point::new(x_min, y_max),
        ])
    }
}

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

    const EPSILON: f64 = 1e-10;

    fn approx_eq(a: f64, b: f64) -> bool {
        (a - b).abs() < EPSILON
    }

    #[test]
    fn test_new_and_accessors() {
        let s = Square::new(Point::new(1.0, 2.0), 3.0);
        assert_eq!(s.center().x(), 1.0);
        assert_eq!(s.center().y(), 2.0);
        assert_eq!(s.side(), 3.0);
    }

    #[test]
    fn test_try_new_accepts_positive() {
        let s = Square::try_new(Point::new(0.0, 0.0), 1.0).unwrap();
        assert_eq!(s.side(), 1.0);
    }

    #[test]
    fn test_try_new_rejects_zero_and_negative() {
        for bad in [0.0, -1.0] {
            let err = Square::try_new(Point::new(0.0, 0.0), bad).unwrap_err();
            assert!(matches!(
                err,
                crate::error::DiagramError::InvalidShapeParameter {
                    shape: "Square",
                    param: "side",
                    ..
                }
            ));
        }
    }

    #[test]
    #[should_panic(expected = "Square side must be > 0")]
    fn test_new_panics_on_zero_side() {
        let _ = Square::new(Point::new(0.0, 0.0), 0.0);
    }

    #[test]
    fn test_area_perimeter_centroid() {
        let s = Square::new(Point::new(0.0, 0.0), 4.0);
        assert!(approx_eq(s.area(), 16.0));
        assert!(approx_eq(s.perimeter(), 16.0));
        assert_eq!(s.centroid(), Point::new(0.0, 0.0));
    }

    #[test]
    fn test_bounds_and_bounding_box() {
        let s = Square::new(Point::new(2.0, 3.0), 4.0);
        let (x0, x1, y0, y1) = s.bounds();
        assert!(approx_eq(x0, 0.0));
        assert!(approx_eq(x1, 4.0));
        assert!(approx_eq(y0, 1.0));
        assert!(approx_eq(y1, 5.0));
        let bb = s.bounding_box();
        assert!(approx_eq(bb.area(), 16.0));
    }

    #[test]
    fn test_contains_point_inside_outside_boundary() {
        let s = Square::new(Point::new(0.0, 0.0), 2.0);
        assert!(s.contains_point(&Point::new(0.0, 0.0)));
        assert!(s.contains_point(&Point::new(1.0, 1.0))); // corner
        assert!(s.contains_point(&Point::new(1.0, 0.0))); // edge midpoint
        assert!(!s.contains_point(&Point::new(1.5, 0.0)));
    }

    #[test]
    fn test_contains_square_strict_equal_partial() {
        let outer = Square::new(Point::new(0.0, 0.0), 4.0);
        let inner = Square::new(Point::new(0.0, 0.0), 2.0);
        let same = Square::new(Point::new(0.0, 0.0), 4.0);
        let partial = Square::new(Point::new(2.0, 0.0), 4.0);

        assert!(outer.contains(&inner));
        assert!(!inner.contains(&outer));
        assert!(outer.contains(&same)); // edges coincident → contained
        assert!(!outer.contains(&partial));
    }

    #[test]
    fn test_intersects_disjoint_touching_partial_nested() {
        let a = Square::new(Point::new(0.0, 0.0), 2.0);
        let disjoint = Square::new(Point::new(5.0, 0.0), 2.0);
        let touching = Square::new(Point::new(2.0, 0.0), 2.0); // share an edge
        let partial = Square::new(Point::new(1.0, 0.0), 2.0);
        let nested = Square::new(Point::new(0.0, 0.0), 1.0);

        assert!(!a.intersects(&disjoint));
        assert!(a.intersects(&touching));
        assert!(a.intersects(&partial));
        assert!(a.intersects(&nested));
    }

    #[test]
    fn test_intersection_area_disjoint_partial_nested() {
        let a = Square::new(Point::new(0.0, 0.0), 2.0);
        let disjoint = Square::new(Point::new(5.0, 0.0), 2.0);
        let partial = Square::new(Point::new(1.0, 0.0), 2.0); // 1×2 overlap
        let nested = Square::new(Point::new(0.0, 0.0), 1.0); // 1×1 inside

        assert!(approx_eq(a.intersection_area(&disjoint), 0.0));
        assert!(approx_eq(a.intersection_area(&partial), 2.0));
        assert!(approx_eq(a.intersection_area(&nested), 1.0));
    }

    #[test]
    fn test_intersection_points_partial_overlap() {
        // Two unit squares overlapping in a 1x1 patch in the corner: the
        // shared boundary contributes two crossings.
        let a = Square::new(Point::new(0.0, 0.0), 2.0);
        let b = Square::new(Point::new(2.0, 2.0), 2.0); // overlaps at upper-right corner
        let pts = a.intersection_points(&b);
        // The overlap region's corners on the shared boundary lie at
        // (1, 1) (single shared corner), and the boundary crossings at
        // (1, 1) only. With the overlap at exactly the corner, expect a
        // single shared point.
        assert!(!pts.is_empty());
    }

    #[test]
    fn test_compute_exclusive_regions_two_disjoint() {
        let a = Square::new(Point::new(0.0, 0.0), 2.0);
        let b = Square::new(Point::new(10.0, 0.0), 2.0);
        let regions = Square::compute_exclusive_regions(&[a, b]);
        assert!(approx_eq(regions[&0b01], 4.0));
        assert!(approx_eq(regions[&0b10], 4.0));
        assert_eq!(regions.get(&0b11).copied().unwrap_or(0.0), 0.0);
    }

    #[test]
    fn test_compute_exclusive_regions_two_partial() {
        let a = Square::new(Point::new(0.0, 0.0), 2.0);
        let b = Square::new(Point::new(1.0, 0.0), 2.0);
        // a∩b is a 1×2 rectangle of area 2; exclusive areas are each 4 − 2 = 2.
        let regions = Square::compute_exclusive_regions(&[a, b]);
        assert!(approx_eq(regions[&0b01], 2.0));
        assert!(approx_eq(regions[&0b10], 2.0));
        assert!(approx_eq(regions[&0b11], 2.0));
    }

    #[test]
    fn test_compute_exclusive_regions_three_way_grid() {
        // Three side-2 squares with centers (0,0), (1,0), (0.5, 1). Each pair
        // overlaps in a 1×2 strip; A∩B∩C is the rectangle x ∈ [0.5−ε, 0.5+ε]
        // (because A and B meet at x ∈ [0,1], C at y ∈ [0,2]; A∩B has
        // y ∈ [-1,1], so A∩B∩C = (x ∈ [0,1]) ∩ (y ∈ [0,1]) = 1.
        let a = Square::new(Point::new(0.0, 0.0), 2.0);
        let b = Square::new(Point::new(1.0, 0.0), 2.0);
        let c = Square::new(Point::new(0.5, 1.0), 2.0);
        let regions = Square::compute_exclusive_regions(&[a, b, c]);

        // Triple-intersection rectangle: x ∈ [0,1], y ∈ [0,1] → area 1.
        assert!(
            approx_eq(regions[&0b111], 1.0),
            "triple ∩ area = {}, expected 1.0",
            regions[&0b111]
        );
    }

    #[test]
    fn test_compute_exclusive_regions_nested() {
        let outer = Square::new(Point::new(0.0, 0.0), 4.0);
        let inner = Square::new(Point::new(0.0, 0.0), 2.0);
        let regions = Square::compute_exclusive_regions(&[outer, inner]);
        assert!(approx_eq(regions[&0b11], 4.0)); // inner is fully inside outer
        assert!(approx_eq(regions[&0b01], 16.0 - 4.0)); // outer-only
        assert_eq!(regions.get(&0b10).copied().unwrap_or(0.0), 0.0);
    }

    /// Central-difference reference for `compute_exclusive_regions`. Returns a
    /// HashMap matching the analytical layout: per region, a length-`3·n_sets`
    /// gradient vector ordered `[x₀, y₀, s₀, x₁, …]`.
    fn fd_exclusive_region_gradients(shapes: &[Square], h: f64) -> HashMap<RegionMask, Vec<f64>> {
        let n_sets = shapes.len();
        let n_params = n_sets * 3;
        let base = Square::compute_exclusive_regions(shapes);

        let mut grads: HashMap<RegionMask, Vec<f64>> =
            base.keys().map(|&m| (m, vec![0.0; n_params])).collect();

        for i in 0..n_sets {
            for k in 0..3 {
                let perturb = |delta: f64| -> HashMap<RegionMask, f64> {
                    let mut copy: Vec<Square> = shapes.to_vec();
                    let (cx, cy, side) =
                        (copy[i].center().x(), copy[i].center().y(), copy[i].side());
                    let new = match k {
                        0 => Square::new(Point::new(cx + delta, cy), side),
                        1 => Square::new(Point::new(cx, cy + delta), side),
                        2 => Square::new(Point::new(cx, cy), side + delta),
                        _ => unreachable!(),
                    };
                    copy[i] = new;
                    Square::compute_exclusive_regions(&copy)
                };
                let plus = perturb(h);
                let minus = perturb(-h);
                for (&mask, g) in grads.iter_mut() {
                    let p = plus.get(&mask).copied().unwrap_or(0.0);
                    let m = minus.get(&mask).copied().unwrap_or(0.0);
                    g[3 * i + k] = (p - m) / (2.0 * h);
                }
            }
        }
        grads
    }

    /// Compare analytical and FD gradient maps mask-by-mask, indexing both by
    /// the analytical map's masks. Tolerance `tol` is per-component absolute.
    fn assert_grad_matches_fd(
        analytical: &HashMap<RegionMask, Vec<f64>>,
        fd: &HashMap<RegionMask, Vec<f64>>,
        tol: f64,
    ) {
        for (&mask, ag) in analytical.iter() {
            let fg = fd
                .get(&mask)
                .expect("FD missing mask present in analytical");
            assert_eq!(ag.len(), fg.len(), "param count mismatch for mask {mask:b}");
            for (k, (&a, &f)) in ag.iter().zip(fg.iter()).enumerate() {
                assert!(
                    (a - f).abs() < tol,
                    "mask {mask:b} param {k}: analytical={a} fd={f} (tol={tol})"
                );
            }
        }
    }

    #[test]
    fn test_gradient_single_square_matches_2s_on_side() {
        // A = s²; ∂A/∂s = 2s, ∂A/∂x = ∂A/∂y = 0.
        let s = 3.0;
        let sq = Square::new(Point::new(1.5, -2.0), s);
        let (areas, grads) = Square::compute_exclusive_regions_with_gradient(&[sq]).unwrap();
        assert!(approx_eq(areas[&0b1], s * s));
        let g = &grads[&0b1];
        assert_eq!(g.len(), 3);
        assert!(approx_eq(g[0], 0.0));
        assert!(approx_eq(g[1], 0.0));
        assert!(approx_eq(g[2], 2.0 * s));
    }

    #[test]
    fn test_gradient_two_squares_partial_overlap_matches_fd() {
        // 1×2 overlap rectangle. Edge-coincident on y=±1, so the gradient
        // exercises the tie-split path; central FD truncation degrades from
        // O(h²) to O(h) at those kinks, hence a 1e-5 tolerance.
        let a = Square::new(Point::new(0.0, 0.0), 2.0);
        let b = Square::new(Point::new(1.0, 0.0), 2.0);
        let (_, grads) = Square::compute_exclusive_regions_with_gradient(&[a, b]).unwrap();
        let fd = fd_exclusive_region_gradients(&[a, b], 1e-6);
        assert_grad_matches_fd(&grads, &fd, 1e-5);
    }

    #[test]
    fn test_gradient_three_squares_overlap_matches_fd() {
        // Triple-overlap configuration with multiple edge ties.
        let a = Square::new(Point::new(0.0, 0.0), 2.0);
        let b = Square::new(Point::new(1.0, 0.0), 2.0);
        let c = Square::new(Point::new(0.5, 1.0), 2.0);
        let (_, grads) = Square::compute_exclusive_regions_with_gradient(&[a, b, c]).unwrap();
        let fd = fd_exclusive_region_gradients(&[a, b, c], 1e-6);
        assert_grad_matches_fd(&grads, &fd, 1e-5);
    }

    #[test]
    fn test_gradient_generic_no_ties_matches_fd_tightly() {
        // Generic configuration with all distinct edges so no tie-split
        // applies; central FD is O(h²) accurate and the gradient should
        // match to ~1e-9 with h=1e-5.
        let a = Square::new(Point::new(0.0, 0.0), 2.3);
        let b = Square::new(Point::new(1.1, 0.4), 1.7);
        let c = Square::new(Point::new(0.6, 1.2), 2.1);
        let (_, grads) = Square::compute_exclusive_regions_with_gradient(&[a, b, c]).unwrap();
        let fd = fd_exclusive_region_gradients(&[a, b, c], 1e-5);
        assert_grad_matches_fd(&grads, &fd, 1e-7);
    }

    #[test]
    fn test_gradient_disjoint_pair_is_zero_on_intersection() {
        let a = Square::new(Point::new(0.0, 0.0), 2.0);
        let b = Square::new(Point::new(10.0, 0.0), 2.0);
        let (_, grads) = Square::compute_exclusive_regions_with_gradient(&[a, b]).unwrap();
        // The {a,b} mask might be absent (sparse discovery) or present with
        // zero area + zero gradient.
        if let Some(g) = grads.get(&0b11) {
            for &v in g {
                assert!(approx_eq(v, 0.0), "expected zero on disjoint pair, got {v}");
            }
        }
        // Singletons must agree with FD.
        let fd = fd_exclusive_region_gradients(&[a, b], 1e-5);
        for &mask in &[0b01_usize, 0b10_usize] {
            let ag = grads.get(&mask).expect("singleton missing");
            let fg = fd.get(&mask).expect("FD singleton missing");
            for (k, (&a, &f)) in ag.iter().zip(fg.iter()).enumerate() {
                assert!(
                    (a - f).abs() < 1e-6,
                    "mask {mask:b} param {k}: analytical={a} fd={f}"
                );
            }
        }
    }

    #[test]
    fn test_gradient_nested_matches_fd() {
        // Inner square fully inside outer; the inner-only exclusive area is 0.
        // Edges are not coincident, so this is a smooth point — tight FD parity.
        let outer = Square::new(Point::new(0.0, 0.0), 4.0);
        let inner = Square::new(Point::new(0.0, 0.0), 2.0);
        let (areas, grads) =
            Square::compute_exclusive_regions_with_gradient(&[outer, inner]).unwrap();
        let fd = fd_exclusive_region_gradients(&[outer, inner], 1e-5);
        assert!(approx_eq(areas[&0b11], 4.0));
        assert!(approx_eq(areas[&0b01], 12.0));
        assert_grad_matches_fd(&grads, &fd, 1e-7);
    }

    #[test]
    fn test_optimizer_params_round_trip() {
        let s = Square::new(Point::new(1.5, -2.0), 3.5);
        let p = s.to_optimizer_params();
        let back = Square::from_optimizer_params(&p);
        assert_eq!(s, back);
    }

    #[test]
    fn test_optimizer_params_from_circle_equal_area() {
        // optimizer_params_from_circle should produce a square of the same
        // area as the seed circle (πr²).
        let r = 2.0;
        let p = Square::optimizer_params_from_circle(0.0, 0.0, r);
        let s = Square::from_optimizer_params(&p);
        assert!(approx_eq(s.area(), PI * r * r));
    }

    #[test]
    fn test_mds_target_distance_zero_overlap_is_diagonal_tangency() {
        // Two unit-area squares (s = 1) with target overlap = 0 should be
        // placed √2 · half_sum = √2 · 1 apart along the diagonal.
        let area = 1.0;
        let d = Square::mds_target_distance(area, area, 0.0).unwrap();
        let s = area.sqrt();
        let expected = 2.0_f64.sqrt() * s;
        assert!(approx_eq(d, expected), "d = {d}, expected {expected}");
    }

    #[test]
    fn test_mds_target_distance_full_overlap_is_zero() {
        // Asking for full overlap of equal-area squares saturates: distance 0.
        let area = 4.0;
        let full_overlap = area; // both squares fully overlapping
        let d = Square::mds_target_distance(area, area, full_overlap).unwrap();
        assert!(approx_eq(d, 0.0));
    }

    #[test]
    fn test_polygonize_returns_4_ccw_vertices_with_correct_area() {
        let s = Square::new(Point::new(0.0, 0.0), 2.0);
        let p = s.polygonize(0);
        assert_eq!(p.vertices().len(), 4);
        // CCW shoelace gives positive area = 4.
        let v = p.vertices();
        let mut shoelace = 0.0;
        for i in 0..4 {
            let j = (i + 1) % 4;
            shoelace += v[i].x() * v[j].y() - v[j].x() * v[i].y();
        }
        assert!(approx_eq(0.5 * shoelace, 4.0));
    }

    #[test]
    fn test_fitter_end_to_end_two_partial_overlap() {
        use crate::{DiagramSpecBuilder, Fitter, InputType};

        // Two sets with sizes 4 (= side 2) and 4, overlap 2. The
        // diagonal-direction MDS inversion places centers √2·(2 − √2) apart
        // (≈ 0.83), and the final-stage optimizer should refine to a layout
        // whose exclusive areas match the spec to high precision.
        let spec = DiagramSpecBuilder::new()
            .set("A", 4.0)
            .set("B", 4.0)
            .intersection(&["A", "B"], 2.0)
            .input_type(InputType::Exclusive)
            .build()
            .unwrap();

        let layout = Fitter::<Square>::new(&spec).seed(42).fit().unwrap();
        // Existence check only: we want to confirm the full
        // Fitter<Square>::fit pipeline runs end-to-end (MDS init via
        // mds_target_distance, final-stage optimisation, layout
        // construction). The numeric loss is sensitive to the non-smooth
        // `(s - |dx|)·(s - |dy|)` overlap; we do not pin a specific basin
        // at this stage.
        let fitted = layout.fitted();
        assert!(
            fitted.values().all(|&v| v.is_finite()),
            "non-finite fitted areas in {fitted:?}"
        );
        assert!(
            layout.loss().is_finite(),
            "non-finite loss {}",
            layout.loss()
        );
        assert_eq!(layout.shapes().len(), 2);
    }

    #[test]
    fn test_distance_between_squares() {
        let a = Square::new(Point::new(0.0, 0.0), 2.0);
        let b_overlap = Square::new(Point::new(1.0, 0.0), 2.0);
        let b_far = Square::new(Point::new(5.0, 0.0), 2.0);
        assert!(approx_eq(a.distance(&b_overlap), 0.0));
        assert!(approx_eq(a.distance(&b_far), 3.0)); // center gap 5 minus two half-widths 1+1
    }

    fn assert_square(actual: &Square, x: f64, y: f64, side: f64) {
        assert!(
            approx_eq(actual.center().x(), x),
            "center.x: {} vs {}",
            actual.center().x(),
            x
        );
        assert!(
            approx_eq(actual.center().y(), y),
            "center.y: {} vs {}",
            actual.center().y(),
            y
        );
        assert!(
            approx_eq(actual.side(), side),
            "side: {} vs {}",
            actual.side(),
            side
        );
    }

    #[test]
    fn test_canonical_venn_layout_n1() {
        let shapes = Square::canonical_venn_layout(1).unwrap();
        assert_eq!(shapes.len(), 1);
        assert_square(&shapes[0], 0.0, 0.0, 2.0);
    }

    #[test]
    fn test_canonical_venn_layout_n2() {
        let shapes = Square::canonical_venn_layout(2).unwrap();
        assert_eq!(shapes.len(), 2);
        assert_square(&shapes[0], -0.4, 0.0, 1.0);
        assert_square(&shapes[1], 0.4, 0.0, 1.0);
    }

    #[test]
    fn test_canonical_venn_layout_n3() {
        let shapes = Square::canonical_venn_layout(3).unwrap();
        assert_eq!(shapes.len(), 3);
        assert_square(&shapes[0], 0.0, 0.36, 1.0);
        assert_square(&shapes[1], 0.42, -0.36, 1.0);
        assert_square(&shapes[2], -0.42, -0.36, 1.0);
    }

    #[test]
    fn test_canonical_venn_layout_unsupported() {
        assert!(Square::canonical_venn_layout(0).is_none());
        assert!(Square::canonical_venn_layout(4).is_none());
        assert!(Square::canonical_venn_layout(5).is_none());
    }

    // ----------------------------------------------------------------
    // Container clipping (S5)
    // ----------------------------------------------------------------

    #[test]
    fn test_clipped_no_squares_complement_equals_container_area() {
        let container = Rectangle::new(Point::new(0.0, 0.0), 5.0, 4.0);
        let areas = compute_exclusive_regions_clipped_squares(&[], &container);
        assert!(approx_eq(areas[&0], 20.0));
    }

    #[test]
    fn test_clipped_single_square_inside_container() {
        let s = Square::new(Point::new(0.0, 0.0), 2.0);
        let container = Rectangle::new(Point::new(0.0, 0.0), 5.0, 4.0);
        let areas = compute_exclusive_regions_clipped_squares(&[s], &container);
        assert!(approx_eq(areas[&0b1], 4.0));
        assert!(approx_eq(areas[&0], 16.0));
    }

    #[test]
    fn test_clipped_single_square_partial_clip() {
        // Square side 4 at (1, 0): bounds x ∈ [-1, 3], y ∈ [-1, 3]. Container
        // 4×4 at origin: bounds [-2, 2] × [-2, 2]. Overlap:
        // x ∈ [-1, 2] (= 3) × y ∈ [-1, 2] (= 3) = 9.
        let s = Square::new(Point::new(1.0, 1.0), 4.0);
        let container = Rectangle::new(Point::new(0.0, 0.0), 4.0, 4.0);
        let areas = compute_exclusive_regions_clipped_squares(&[s], &container);
        assert!(approx_eq(areas[&0b1], 9.0));
        assert!(approx_eq(areas[&0], 16.0 - 9.0));
    }

    #[test]
    fn test_clipped_container_fully_inside_square() {
        let s = Square::new(Point::new(0.0, 0.0), 100.0);
        let container = Rectangle::new(Point::new(0.0, 0.0), 4.0, 3.0);
        let areas = compute_exclusive_regions_clipped_squares(&[s], &container);
        assert!(approx_eq(areas[&0b1], 12.0));
        assert!(approx_eq(areas[&0], 0.0));
    }

    /// Pack `(squares, container)` into the same flat parameter layout the
    /// optimiser uses: per-square `[x, y, side]` blocks, then container
    /// `[x_c, y_c, u_c, v_c]`.
    fn pack_clipped_square_params(squares: &[Square], container: &Rectangle) -> Vec<f64> {
        let mut p = Vec::with_capacity(squares.len() * 3 + 4);
        for s in squares {
            p.push(s.center().x());
            p.push(s.center().y());
            p.push(s.side());
        }
        p.extend(container.to_optimizer_params());
        p
    }

    fn unpack_clipped_square_params(p: &[f64], n_sets: usize) -> (Vec<Square>, Rectangle) {
        let squares: Vec<Square> = (0..n_sets)
            .map(|i| {
                Square::new(
                    Point::new(p[3 * i], p[3 * i + 1]),
                    p[3 * i + 2].max(f64::MIN_POSITIVE),
                )
            })
            .collect();
        let container = Rectangle::from_optimizer_params(&p[3 * n_sets..3 * n_sets + 4]);
        (squares, container)
    }

    fn assert_clipped_square_gradient_matches_fd(
        squares: &[Square],
        container: &Rectangle,
        h: f64,
        tol: f64,
        label: &str,
    ) {
        let (areas, grads) =
            compute_exclusive_regions_clipped_with_gradient_squares(squares, container);
        let n_sets = squares.len();
        let n_params = n_sets * 3 + 4;
        let p0 = pack_clipped_square_params(squares, container);

        for (mask, analytic) in &grads {
            assert_eq!(
                analytic.len(),
                n_params,
                "{}: gradient length {} ≠ expected {}",
                label,
                analytic.len(),
                n_params
            );
            let mut fd = vec![0.0; n_params];
            for i in 0..n_params {
                let mut plus = p0.clone();
                let mut minus = p0.clone();
                plus[i] += h;
                minus[i] -= h;
                let (sp, kp) = unpack_clipped_square_params(&plus, n_sets);
                let (sm, km) = unpack_clipped_square_params(&minus, n_sets);
                let ap = compute_exclusive_regions_clipped_squares(&sp, &kp)
                    .get(mask)
                    .copied()
                    .unwrap_or(0.0);
                let am = compute_exclusive_regions_clipped_squares(&sm, &km)
                    .get(mask)
                    .copied()
                    .unwrap_or(0.0);
                fd[i] = (ap - am) / (2.0 * h);
            }
            let direct = compute_exclusive_regions_clipped_squares(squares, container)
                .get(mask)
                .copied()
                .unwrap_or(0.0);
            let analytic_area = areas.get(mask).copied().unwrap_or(0.0);
            assert!(
                (analytic_area - direct).abs() < 1e-12,
                "{}: mask {:b} area {} vs direct {} mismatch",
                label,
                mask,
                analytic_area,
                direct
            );
            let diff_norm: f64 = analytic
                .iter()
                .zip(fd.iter())
                .map(|(a, b)| (a - b).powi(2))
                .sum::<f64>()
                .sqrt();
            let fd_norm: f64 = fd.iter().map(|b| b * b).sum::<f64>().sqrt();
            let rel = if fd_norm > 1e-9 {
                diff_norm / fd_norm
            } else {
                diff_norm
            };
            assert!(
                rel < tol,
                "{}: mask {:b} analytic vs FD mismatch (rel={:.3e}, |fd|={:.3e})\n  analytic={:?}\n  fd      ={:?}",
                label,
                mask,
                rel,
                fd_norm,
                analytic,
                fd
            );
        }
    }

    #[test]
    fn test_clipped_grad_two_squares_inside_box_matches_fd() {
        let squares = vec![
            Square::new(Point::new(-0.6, 0.07), 1.4),
            Square::new(Point::new(0.62, -0.04), 1.5),
        ];
        let container = Rectangle::new(Point::new(0.0, 0.0), 6.0, 5.0);
        assert_clipped_square_gradient_matches_fd(
            &squares,
            &container,
            1e-5,
            1e-7,
            "two_squares_inside",
        );
    }

    #[test]
    fn test_clipped_grad_two_squares_one_clipped_matches_fd() {
        let squares = vec![
            Square::new(Point::new(0.0, 0.07), 1.5),
            Square::new(Point::new(1.4, -0.03), 1.5),
        ];
        let container = Rectangle::new(Point::new(0.0, 0.0), 3.6, 3.0);
        assert_clipped_square_gradient_matches_fd(
            &squares,
            &container,
            1e-5,
            1e-6,
            "two_squares_one_clipped",
        );
    }

    #[test]
    fn test_clipped_grad_three_squares_inside_box_matches_fd() {
        let squares = vec![
            Square::new(Point::new(-0.5, -0.3), 1.4),
            Square::new(Point::new(0.5, -0.3), 1.5),
            Square::new(Point::new(0.0, 0.55), 1.6),
        ];
        let container = Rectangle::new(Point::new(0.0, 0.05), 3.5, 3.5);
        assert_clipped_square_gradient_matches_fd(
            &squares,
            &container,
            1e-5,
            1e-7,
            "three_squares_inside",
        );
    }

    #[test]
    fn fit_two_squares_with_complement_runs_to_finite_loss() {
        use crate::{DiagramSpecBuilder, Fitter, InputType};

        let spec = DiagramSpecBuilder::new()
            .set("A", 4.0)
            .set("B", 4.0)
            .intersection(&["A", "B"], 1.0)
            .complement(20.0)
            .input_type(InputType::Exclusive)
            .build()
            .unwrap();

        let layout = Fitter::<Square>::new(&spec).seed(42).fit().unwrap();
        let fitted = layout.fitted();
        assert!(
            fitted.values().all(|&v| v.is_finite()),
            "non-finite fitted areas in {fitted:?}"
        );
        assert!(
            layout.loss().is_finite(),
            "non-finite loss {}",
            layout.loss()
        );
        let container = layout
            .container()
            .expect("complement spec carries container");
        // Universe = 4 + 4 + 1 + 20 = 29.
        assert!(
            (container.area() - 29.0).abs() / 29.0 < 0.1,
            "container area {} should be near universe 29",
            container.area()
        );
    }
}