eunoia 0.14.0

//! Ellipse shape implementation.
//!
//! This module provides an ellipse implementation for use in Euler and Venn diagrams.
//! Ellipses are more flexible than circles and can represent elongated or rotated regions,
//! making them useful for more accurate area-proportional diagrams.
//!
//! # Representation
//!
//! An ellipse is defined by:
//! - **Center point**: The center of the ellipse
//! - **Semi-major axis** (a): Half the length of the longest diameter
//! - **Semi-minor axis** (b): Half the length of the shortest diameter  
//! - **Rotation**: Rotation angle in radians (counterclockwise from x-axis)
//!
//! # Area Calculations
//!
//! The module provides several specialized area computation methods:
//!
//! - **Sector area**: The area swept by a radius from the center through an angle θ
//! - **Segment area**: The area between the ellipse boundary and a chord connecting two points
//!
//! These primitives are essential for computing intersection areas in Euler diagrams.
//!
//! # Examples
//!
//! ```
//! use eunoia::geometry::shapes::Ellipse;
//! use eunoia::geometry::traits::Area;
//! use eunoia::geometry::primitives::Point;
//!
//! // Create an ellipse centered at origin
//! let ellipse = Ellipse::new(Point::new(0.0, 0.0), 5.0, 3.0, 0.0);
//!
//! // Compute total area
//! let area = ellipse.area();
//!
//! // Compute sector area for π/4 radians
//! let sector = ellipse.sector_area(std::f64::consts::PI / 4.0);
//! ```

use std::f64::consts::PI;

use nalgebra::Matrix2;

use crate::geometry::diagram::RegionMask;
use crate::geometry::primitives::Point;
use crate::geometry::projective::Conic;
use crate::geometry::shapes::{Circle, Polygon, Rectangle};
use crate::geometry::traits::{
    Area, BoundingBox, Centroid, Closed, DiagramShape, Perimeter, Polygonize,
};

/// An ellipse defined by center, semi-major and semi-minor axes, and rotation.
///
/// Ellipses are oval-shaped closed curves that generalize circles. They are particularly
/// useful in Euler diagrams when sets have elongated or directional relationships.
///
/// # Representation
///
/// The ellipse is represented in standard form with:
/// - A center point `(h, k)`
/// - Semi-major axis length `a` (≥ semi-minor axis)
/// - Semi-minor axis length `b` (> 0)
/// - Rotation angle `φ` in radians (counterclockwise from the positive x-axis)
///
/// The canonical equation in the ellipse's local coordinate system is:
/// ```text
/// (x/a)² + (y/b)² = 1
/// ```
///
/// # Examples
///
/// ```
/// use eunoia::geometry::shapes::Ellipse;
/// use eunoia::geometry::primitives::Point;
/// use eunoia::geometry::traits::Area;
///
/// // Create an ellipse at the origin with no rotation
/// let ellipse = Ellipse::new(Point::new(0.0, 0.0), 5.0, 3.0, 0.0);
/// assert!((ellipse.area() - 47.12).abs() < 0.01);
///
/// // Create a rotated ellipse
/// let rotated = Ellipse::new(
///     Point::new(2.0, 3.0),
///     4.0,
///     2.0,
///     std::f64::consts::PI / 4.0  // 45 degrees
/// );
/// ```
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct Ellipse {
    center: Point,
    semi_major: f64,
    semi_minor: f64,
    rotation: f64, // in radians
}

impl Ellipse {
    /// Creates a new ellipse with the specified parameters.
    ///
    /// # Arguments
    ///
    /// * `center` - The center point of the ellipse
    /// * `semi_major` - The semi-major axis length (must be > 0)
    /// * `semi_minor` - The semi-minor axis length (must be > 0)
    /// * `rotation` - Rotation angle in radians (counterclockwise from x-axis)
    ///
    /// # Panics
    ///
    /// Panics if either axis length is `<= 0`. Use [`Ellipse::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`.
    ///
    /// # Examples
    ///
    /// ```
    /// use eunoia::geometry::shapes::Ellipse;
    /// use eunoia::geometry::primitives::Point;
    ///
    /// let ellipse = Ellipse::new(Point::new(1.0, 2.0), 4.0, 3.0, 0.0);
    /// assert_eq!(ellipse.semi_major(), 4.0);
    /// assert_eq!(ellipse.semi_minor(), 3.0);
    /// ```
    pub fn new(center: Point, semi_major: f64, semi_minor: f64, rotation: f64) -> Self {
        assert!(
            semi_major > 0.0,
            "Ellipse semi_major must be > 0, got {}",
            semi_major
        );
        assert!(
            semi_minor > 0.0,
            "Ellipse semi_minor must be > 0, got {}",
            semi_minor
        );
        Self {
            center,
            semi_major,
            semi_minor,
            rotation,
        }
    }

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

    /// Creates a new ellipse from radius and log-aspect ratio parameterization.
    ///
    /// This parameterization is optimized for numerical optimization:
    /// - `radius` controls the overall size (geometric mean of axes)
    /// - `log_aspect` controls elongation on a log scale (more uniform landscape)
    /// - `log_aspect = 0.0` gives a circle (aspect_ratio = 1.0)
    /// - `log_aspect < 0` gives elongated ellipses
    /// - `log_aspect = -0.693` gives aspect_ratio ≈ 0.5 (2:1 elongation)
    ///
    /// # Arguments
    ///
    /// * `center` - The center point of the ellipse
    /// * `radius` - Geometric mean radius: sqrt(semi_major * semi_minor)
    /// * `log_aspect` - Log of aspect ratio (semi_minor/semi_major), clamped to [-1.609, 0]
    ///   which corresponds to aspect_ratio in [0.2, 1.0]
    /// * `rotation` - Rotation angle in radians
    ///
    /// # Examples
    ///
    /// ```
    /// use eunoia::geometry::shapes::Ellipse;
    /// use eunoia::geometry::primitives::Point;
    ///
    /// // Circle: log_aspect = 0.0
    /// let circle = Ellipse::from_radius_ratio(Point::new(0.0, 0.0), 2.0, 0.0, 0.0);
    /// assert!((circle.semi_major() - 2.0).abs() < 1e-10);
    /// assert!((circle.semi_minor() - 2.0).abs() < 1e-10);
    ///
    /// // Elongated ellipse: log_aspect = -0.693 (aspect ≈ 0.5)
    /// let ellipse = Ellipse::from_radius_ratio(Point::new(0.0, 0.0), 2.0, -0.693147, 0.0);
    /// // radius = sqrt(a * b) = 2.0, aspect = 0.5
    /// // => a = 2.828, b = 1.414
    /// assert!((ellipse.semi_major() - 2.828).abs() < 0.01);
    /// assert!((ellipse.semi_minor() - 1.414).abs() < 0.01);
    /// ```
    pub fn from_radius_ratio(center: Point, radius: f64, log_aspect: f64, rotation: f64) -> Self {
        let radius = radius.abs();

        // Clamp log_aspect to [-ln(5), 0] which gives aspect_ratio in [0.2, 1.0]
        let log_aspect = log_aspect.clamp(-1.6094379124341003, 0.0); // -ln(5) to 0
        let aspect_ratio = log_aspect.exp();

        // Normalize rotation to [0, π) since ellipse has 180° rotational symmetry
        let rotation = rotation.rem_euclid(PI);

        // Given: r = sqrt(a * b) and aspect = b/a
        // Solving: a = r / sqrt(aspect), b = r * sqrt(aspect)
        let semi_major = radius / aspect_ratio.sqrt();
        let semi_minor = radius * aspect_ratio.sqrt();

        Self::new(center, semi_major, semi_minor, rotation)
    }

    pub fn from_conic(conic: Conic) -> Option<Self> {
        let c = conic.matrix();

        // Extract algebraic coefficients from the symmetric matrix
        // Q = [ A   B/2  D/2 ]
        //     [ B/2 C    E/2 ]
        //     [ D/2 E/2  F   ]
        let m1 = c[(0, 0)];
        let m2 = 2.0 * c[(0, 1)];
        let m3 = c[(1, 1)];
        let m4 = 2.0 * c[(0, 2)];
        let m5 = 2.0 * c[(1, 2)];
        let m6 = c[(2, 2)];

        // Check ellipse condition (negative discriminant required)
        if m2 * m2 - 4.0 * m1 * m3 >= 0.0 {
            return None;
        }

        // Quadratic form matrix
        let m = Matrix2::new(m1, m2 / 2.0, m2 / 2.0, m3);

        // Solve for center:
        // M * [xc, yc]^T = -1/2 * [D, E]^T
        let rhs = nalgebra::Vector2::new(-m4 / 2.0, -m5 / 2.0);
        let center_vec = m.lu().solve(&rhs)?;
        let xc = center_vec[0];
        let yc = center_vec[1];

        // Compute translated constant term F̄
        let f_bar = m6 + m1 * xc * xc + m2 * xc * yc + m3 * yc * yc + m4 * xc + m5 * yc;

        if f_bar >= 0.0 {
            return None; // Not an ellipse
        }

        // Eigen-decompose the quadratic part
        let eig = m.symmetric_eigen();
        let lambda1 = eig.eigenvalues[0];
        let lambda2 = eig.eigenvalues[1];
        let v = eig.eigenvectors;

        // Identify major/minor axes:
        // smaller eigenvalue ⇒ bigger radius
        let (lambda_major, lambda_minor, vec_major) = if lambda1 < lambda2 {
            (lambda1, lambda2, v.column(0))
        } else {
            (lambda2, lambda1, v.column(1))
        };

        // Radii from: a² = -F̄ / λ
        let a2 = -f_bar / lambda_major;
        let b2 = -f_bar / lambda_minor;

        if a2 <= 0.0 || b2 <= 0.0 {
            return None;
        }

        let a = a2.sqrt();
        let b = b2.sqrt();

        // Rotation angle from the major-axis eigenvector
        let vx = vec_major[0];
        let vy = vec_major[1];
        let mut phi = vy.atan2(vx);

        // Normalize angle into [0, π)
        if phi < 0.0 {
            phi += std::f64::consts::PI;
        }

        Some(Ellipse {
            center: Point::new(xc, yc),
            semi_major: a,
            semi_minor: b,
            rotation: phi,
        })
    }

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

    /// Returns the semi-major axis length.
    pub fn semi_major(&self) -> f64 {
        self.semi_major
    }

    /// Returns the semi-minor axis length.
    pub fn semi_minor(&self) -> f64 {
        self.semi_minor
    }

    /// Returns the rotation angle in radians.
    pub fn rotation(&self) -> f64 {
        self.rotation
    }

    /// Computes the area of an elliptical sector from the center through angle θ.
    ///
    /// An elliptical sector is the region bounded by:
    /// - Two radii from the center to the ellipse boundary
    /// - The ellipse arc between those radii
    ///
    /// This uses the exact formula for elliptical sectors, which accounts for
    /// the non-uniform curvature of the ellipse (unlike circular sectors where
    /// area is simply ½r²θ).
    ///
    /// # Arguments
    ///
    /// * `theta` - The angle in radians from the semi-major axis
    ///
    /// # Returns
    ///
    /// The area of the sector from 0 to θ
    ///
    /// # Examples
    ///
    /// ```
    /// use eunoia::geometry::shapes::Ellipse;
    /// use eunoia::geometry::primitives::Point;
    /// use eunoia::geometry::traits::Area;
    ///
    /// let ellipse = Ellipse::new(Point::new(0.0, 0.0), 4.0, 2.0, 0.0);
    ///
    /// // Quarter sector (π/2 radians)
    /// let quarter = ellipse.sector_area(std::f64::consts::PI / 2.0);
    /// assert!((quarter - ellipse.area() / 4.0).abs() < 1e-10);
    /// ```
    pub fn sector_area(&self, theta: f64) -> f64 {
        let a = self.semi_major;
        let b = self.semi_minor;

        let num = (b - a) * (2.0 * theta).sin();
        let den = a + b + (b - a) * (2.0 * theta).cos();

        0.5 * a * b * (theta - num.atan2(den))
    }

    /// Computes the area of an elliptical sector between two angles.
    ///
    /// This is equivalent to `sector_area(theta1) - sector_area(theta0)` but
    /// handles counter-clockwise angle wrapping automatically.
    ///
    /// # Arguments
    ///
    /// * `theta0` - Starting angle in radians
    /// * `theta1` - Ending angle in radians
    ///
    /// # Returns
    ///
    /// The area of the sector from θ₀ to θ₁ (counter-clockwise)
    ///
    /// # Note
    ///
    /// If θ₁ < θ₀, the function adds 2π to θ₁ to compute the counter-clockwise
    /// sector that wraps through angle 0.
    ///
    /// # Examples
    ///
    /// ```
    /// use eunoia::geometry::shapes::Ellipse;
    /// use eunoia::geometry::primitives::Point;
    ///
    /// let ellipse = Ellipse::new(Point::new(0.0, 0.0), 3.0, 2.0, 0.0);
    ///
    /// // Sector from 45° to 135°
    /// let sector = ellipse.sector_area_between(
    ///     std::f64::consts::PI / 4.0,
    ///     3.0 * std::f64::consts::PI / 4.0
    /// );
    /// assert!(sector > 0.0);
    /// ```
    pub fn sector_area_between(&self, theta0: f64, theta1: f64) -> f64 {
        let t0 = theta0;
        let mut t1 = theta1;

        // Ensure CCW ordering.
        if t1 < t0 {
            t1 += 2.0 * PI;
        }

        self.sector_area(t1) - self.sector_area(t0)
    }

    /// Computes the area of an ellipse segment between two boundary points.
    ///
    /// An ellipse segment is the region bounded by:
    /// - The ellipse arc between two points on the boundary
    /// - The chord (straight line) connecting those two points
    ///
    /// This method automatically handles:
    /// - Coordinate system transformation to the ellipse's local frame
    /// - Counter-clockwise angle ordering
    /// - Minor arc (≤ 180°) vs major arc (> 180°) selection
    ///
    /// # Algorithm
    ///
    /// 1. Transform points to ellipse coordinate system (centered, unrotated)
    /// 2. Compute angles θ₀ and θ₁ for each point
    /// 3. Ensure counter-clockwise ordering (add 2π if needed)
    /// 4. Calculate segment as: sector - triangle (for minor arc)
    ///    or: total_area - complementary_sector + triangle (for major arc)
    ///
    /// # Arguments
    ///
    /// * `p0` - First boundary point (in world coordinates)
    /// * `p1` - Second boundary point (in world coordinates)
    ///
    /// # Returns
    ///
    /// The area of the segment. When the angular span > π, returns the **major arc**
    /// segment (the larger of the two possible segments).
    ///
    /// # Examples
    ///
    /// ```
    /// use eunoia::geometry::shapes::Ellipse;
    /// use eunoia::geometry::primitives::Point;
    /// use eunoia::geometry::traits::Area;
    ///
    /// let ellipse = Ellipse::new(Point::new(0.0, 0.0), 5.0, 3.0, 0.0);
    ///
    /// // Segment between points on opposite sides (semicircle)
    /// let p0 = Point::new(5.0, 0.0);   // Right end of major axis
    /// let p1 = Point::new(-5.0, 0.0);  // Left end of major axis
    /// let segment = ellipse.ellipse_segment(p0, p1);
    ///
    /// // Should be approximately half the ellipse area
    /// assert!((segment - ellipse.area() / 2.0).abs() < 1e-10);
    /// ```
    pub fn ellipse_segment(&self, p0: Point, p1: Point) -> f64 {
        // 1. Move into ellipse coordinate system.
        let p0 = p0.to_ellipse_frame(self);
        let p1 = p1.to_ellipse_frame(self);

        // 2. Compute angles
        let theta0 = p0.angle_from_origin();
        let mut theta1 = p1.angle_from_origin();

        // 3. Ensure CCW ordering.
        if theta1 < theta0 {
            theta1 += 2.0 * PI;
        }

        // 4. Triangle correction (signed area of parallelogram / 2).
        let triangle = 0.5 * (p1.x() * p0.y() - p0.x() * p1.y()).abs();

        // 5. Minor or major arc?
        if (theta1 - theta0) <= PI {
            self.sector_area(theta1) - self.sector_area(theta0) - triangle
        } else {
            self.area() - (self.sector_area(theta0 + 2.0 * PI) - self.sector_area(theta1))
                + triangle
        }
    }

    /// Computes the lens-shaped intersection area between two ellipses with exactly two intersection points.
    ///
    /// Follows eulerr's polysegments: processes both arcs around the intersection perimeter.
    fn compute_lens_area(&self, other: &Self, p0: &Point, p1: &Point) -> f64 {
        // Two arcs: p1→p0 and p0→p1
        // For each arc, compute triangle + min(segment for common parents)

        // Arc 1: p1 → p0
        let seg1_10 = self.ellipse_segment(*p1, *p0);
        let seg2_10 = other.ellipse_segment(*p1, *p0);
        let triangle_10 = 0.5 * ((p0.x() + p1.x()) * (p0.y() - p1.y()));
        let arc1 = triangle_10 + seg1_10.min(seg2_10);

        // Arc 2: p0 → p1
        let seg1_01 = self.ellipse_segment(*p0, *p1);
        let seg2_01 = other.ellipse_segment(*p0, *p1);
        let triangle_01 = 0.5 * ((p1.x() + p0.x()) * (p1.y() - p0.y()));
        let arc2 = triangle_01 + seg1_01.min(seg2_01);

        arc1 + arc2
    }

    /// Computes intersection area for cases with 3 or 4 intersection points.
    ///
    /// This follows eulerr's polysegments algorithm.
    fn compute_multi_point_intersection_area(&self, other: &Self, points: &[Point]) -> f64 {
        if points.len() < 3 {
            return 0.0;
        }

        // Use the centroid of intersection points as a reference center
        let cx = points.iter().map(|p| p.x()).sum::<f64>() / points.len() as f64;
        let cy = points.iter().map(|p| p.y()).sum::<f64>() / points.len() as f64;

        // Sort points by angle around the centroid (match eulerr's atan2(x-cx, y-cy) order)
        let mut sorted_points: Vec<Point> = points.to_vec();
        sorted_points.sort_by(|a, b| {
            let angle_a = (a.x() - cx).atan2(a.y() - cy);
            let angle_b = (b.x() - cx).atan2(b.y() - cy);
            angle_a.partial_cmp(&angle_b).unwrap()
        });

        // Compute area using polysegments algorithm
        let mut total_area = 0.0;
        let n = sorted_points.len();

        for k in 0..n {
            let l = if k == 0 { n - 1 } else { k - 1 };

            let pi = sorted_points[k]; // current point (matches eulerr's i)
            let pj = sorted_points[l]; // previous point (matches eulerr's j)

            // Triangle contribution (shoelace formula) - matches eulerr line 90
            let triangle = 0.5 * ((pj.x() + pi.x()) * (pj.y() - pi.y()));

            // Segments from i to j (matches eulerr line 86: ellipse_segment(..., points[i], points[j]))
            let seg1 = self.ellipse_segment(pi, pj);
            let seg2 = other.ellipse_segment(pi, pj);

            // Add triangle + min(segment) as in eulerr line 89-91
            total_area += triangle + seg1.min(seg2);
        }

        total_area
    }
}

impl Area for Ellipse {
    /// Computes the total area of the ellipse using the formula A = πab.
    ///
    /// # Examples
    ///
    /// ```
    /// use eunoia::geometry::shapes::Ellipse;
    /// use eunoia::geometry::primitives::Point;
    /// use eunoia::geometry::traits::Area;
    ///
    /// let ellipse = Ellipse::new(Point::new(0.0, 0.0), 5.0, 3.0, 0.0);
    /// let area = ellipse.area();
    /// assert!((area - std::f64::consts::PI * 15.0).abs() < 1e-10);
    /// ```
    fn area(&self) -> f64 {
        PI * self.semi_major * self.semi_minor
    }
}

impl Perimeter for Ellipse {
    /// Computes an approximation of the ellipse perimeter.
    ///
    /// Uses Ramanujan's second approximation formula, which provides excellent
    /// accuracy for all ellipses:
    ///
    /// ```text
    /// P ≈ π(a + b)(1 + 3h / (10 + √(4 - 3h)))
    /// where h = ((a - b) / (a + b))²
    /// ```
    ///
    /// This approximation has a relative error of less than 0.01% for typical cases.
    ///
    /// # Examples
    ///
    /// ```
    /// use eunoia::geometry::shapes::Ellipse;
    /// use eunoia::geometry::primitives::Point;
    /// use eunoia::geometry::traits::Perimeter;
    ///
    /// let ellipse = Ellipse::new(Point::new(0.0, 0.0), 5.0, 3.0, 0.0);
    /// let perimeter = ellipse.perimeter();
    /// assert!(perimeter > 0.0);
    /// ```
    fn perimeter(&self) -> f64 {
        // Approximation using Ramanujan's second formula
        let a = self.semi_major;
        let b = self.semi_minor;
        let h = ((a - b).powi(2)) / ((a + b).powi(2));
        PI * (a + b) * (1.0 + (3.0 * h) / (10.0 + (4.0 - 3.0 * h).sqrt()))
    }
}

impl Centroid for Ellipse {
    /// Returns the centroid of the ellipse, which is its center point.
    fn centroid(&self) -> Point {
        self.center
    }
}

impl BoundingBox for Ellipse {
    /// Computes the axis-aligned bounding box that contains the ellipse.
    ///
    /// For a rotated ellipse, this computes the smallest axis-aligned rectangle
    /// that fully contains the ellipse. The calculation accounts for the rotation
    /// by projecting the ellipse axes onto the coordinate axes.
    ///
    /// # Examples
    ///
    /// ```
    /// use eunoia::geometry::shapes::Ellipse;
    /// use eunoia::geometry::primitives::Point;
    /// use eunoia::geometry::traits::BoundingBox;
    ///
    /// let ellipse = Ellipse::new(Point::new(0.0, 0.0), 5.0, 3.0, 0.0);
    /// let bbox = ellipse.bounding_box();
    /// // For unrotated ellipse, bbox dimensions equal 2*semi-major and 2*semi-minor
    /// ```
    fn bounding_box(&self) -> Rectangle {
        let cos_theta = self.rotation.cos();
        let sin_theta = self.rotation.sin();

        let width = 2.0
            * ((self.semi_major * cos_theta).powi(2) + (self.semi_minor * sin_theta).powi(2))
                .sqrt();
        let height = 2.0
            * ((self.semi_major * sin_theta).powi(2) + (self.semi_minor * cos_theta).powi(2))
                .sqrt();

        Rectangle::new(self.center, width, height)
    }
}

impl Closed for Ellipse {
    fn contains(&self, other: &Self) -> bool {
        // Quick check: if other's center is outside self, it can't be contained
        if !self.contains_point(&other.center) {
            return false;
        }

        // If other is larger than self (by area), it can't be contained
        if other.area() > self.area() {
            return false;
        }

        // Convert both ellipses to conics
        let c1 = Conic::from_ellipse(*self);
        let c2 = Conic::from_ellipse(*other);

        // Check for boundary intersections
        let intersection_points = c1.intersect_conic(&c2);

        // If no intersections found and center is inside, other is contained
        if intersection_points.is_empty() {
            return true;
        }

        // If intersections were found, verify they're not real by checking
        // the extreme points of other (ends of major/minor axes)
        let phi = other.rotation;
        let cos_phi = phi.cos();
        let sin_phi = phi.sin();

        // Check points at the ends of the semi-major axis
        let major_offset_x = other.semi_major * cos_phi;
        let major_offset_y = other.semi_major * sin_phi;
        let p1 = Point::new(
            other.center.x() + major_offset_x,
            other.center.y() + major_offset_y,
        );
        let p2 = Point::new(
            other.center.x() - major_offset_x,
            other.center.y() - major_offset_y,
        );

        // Check points at the ends of the semi-minor axis
        let minor_offset_x = other.semi_minor * -sin_phi;
        let minor_offset_y = other.semi_minor * cos_phi;
        let p3 = Point::new(
            other.center.x() + minor_offset_x,
            other.center.y() + minor_offset_y,
        );
        let p4 = Point::new(
            other.center.x() - minor_offset_x,
            other.center.y() - minor_offset_y,
        );

        // All extreme points must be inside or on self
        self.contains_point(&p1)
            && self.contains_point(&p2)
            && self.contains_point(&p3)
            && self.contains_point(&p4)
    }

    fn contains_point(&self, point: &Point) -> bool {
        // Transform point to ellipse's local coordinate system
        let dx = point.x() - self.center.x();
        let dy = point.y() - self.center.y();

        let cos_phi = self.rotation.cos();
        let sin_phi = self.rotation.sin();

        // Rotate point to align with ellipse axes
        let x_local = dx * cos_phi + dy * sin_phi;
        // Correct inverse rotation for y': y_local = -dx*sin(phi) + dy*cos(phi)
        let y_local = -dx * sin_phi + dy * cos_phi;

        // Check if point is inside using the ellipse equation
        (x_local * x_local) / (self.semi_major * self.semi_major)
            + (y_local * y_local) / (self.semi_minor * self.semi_minor)
            <= 1.0
    }

    fn intersects(&self, other: &Self) -> bool {
        // Quick rejection: if centers are too far apart, they can't intersect.
        // The farthest reach from an ellipse's center is `max(semi_major,
        // semi_minor)` — the field names don't enforce ordering (the optimizer
        // parameterises both axes independently in log-space, so either can
        // be larger), so we have to take the max of the two rather than
        // trusting `semi_major` alone.
        let center_distance = self.center.distance(&other.center);
        let max_reach_self = self.semi_major.max(self.semi_minor);
        let max_reach_other = other.semi_major.max(other.semi_minor);

        if center_distance > max_reach_self + max_reach_other {
            return false;
        }

        // Check if one contains the other (no intersection in that case)
        if self.contains(other) || other.contains(self) {
            return false;
        }

        // Convert to conics and check for intersection points
        let c1 = Conic::from_ellipse(*self);
        let c2 = Conic::from_ellipse(*other);

        let intersection_points = c1.intersect_conic(&c2);

        !intersection_points.is_empty()
    }

    fn intersection_area(&self, other: &Self) -> f64 {
        // Check if one ellipse contains the other
        if self.contains(other) {
            return other.area();
        }
        if other.contains(self) {
            return self.area();
        }

        // Get intersection points
        let points = self.intersection_points(other);
        let n = points.len();

        match n {
            0 => 0.0, // No intersection
            1 => {
                // Single point of tangency - negligible area
                0.0
            }
            2 => {
                // Two intersection points - compute the lens area
                self.compute_lens_area(other, &points[0], &points[1])
            }
            _ => {
                // Three or four intersection points
                self.compute_multi_point_intersection_area(other, &points)
            }
        }
    }

    fn intersection_points(&self, other: &Self) -> Vec<Point> {
        // Convert both ellipses to conics
        let c1 = Conic::from_ellipse(*self);
        let c2 = Conic::from_ellipse(*other);

        // Get homogeneous intersection points
        let homogeneous_points = c1.intersect_conic(&c2);

        // Convert to Cartesian points
        homogeneous_points
            .into_iter()
            .map(Point::from_homogeneous)
            .collect()
    }
}

/// One arc of an ellipse on the boundary of a region. The arc is traversed
/// CCW around the region; for intersections of convex sets this is also CCW
/// around the owning ellipse, so `delta_t > 0` always.
///
/// `t_start` is the canonical-frame parameter value at the arc's start
/// (`t = atan2(v·a, u·b)` where `(u, v)` is the point in the unrotated ellipse
/// frame). The endpoint can be recovered as `t_start + delta_t`; sin/cos of
/// that continuous value coincide with sin/cos of the wrapped atan2 endpoint.
#[derive(Debug, Clone, Copy)]
pub(crate) struct EllipseArc {
    pub ellipse: usize,
    pub t_start: f64,
    pub delta_t: f64,
}

/// World-frame point on ellipse `e` at canonical parameter `t`.
fn ellipse_point_at(e: &Ellipse, t: f64) -> Point {
    let cos_t = t.cos();
    let sin_t = t.sin();
    let u = e.semi_major * cos_t;
    let v = e.semi_minor * sin_t;
    let cos_phi = e.rotation.cos();
    let sin_phi = e.rotation.sin();
    Point::new(
        e.center.x() + u * cos_phi - v * sin_phi,
        e.center.y() + u * sin_phi + v * cos_phi,
    )
}

/// Canonical-frame parameter `t` for a world-frame point on (or near) `e`.
fn ellipse_param_for_point(e: &Ellipse, p: &Point) -> f64 {
    let local = p.to_ellipse_frame(e);
    // t such that (cos t, sin t) ∝ (local.x()/a, local.y()/b). Multiplying both
    // arguments by a·b > 0 is equivalent and avoids dividing by potentially
    // small axes:
    (local.y() * e.semi_major).atan2(local.x() * e.semi_minor)
}

/// Implicit-form value `(u/a)² + (v/b)²` for point `p` in ellipse `e`'s frame.
/// `< 1` strictly inside, `= 1` on the boundary, `> 1` outside.
fn ellipse_implicit_value(p: &Point, e: &Ellipse) -> f64 {
    let local = p.to_ellipse_frame(e);
    (local.x() / e.semi_major).powi(2) + (local.y() / e.semi_minor).powi(2)
}

/// Decide whether an arc whose midpoint is on `∂E_j` is owned by `j` for
/// the purpose of region-boundary contributions. The midpoint must be inside
/// every other mask ellipse; when its lies on another mask ellipse's
/// boundary (boundaries coincide), the smaller-index ellipse owns the arc to
/// avoid double-counting (the identical-ellipses case in particular would
/// otherwise emit one full-circle arc per ellipse).
///
/// Tolerance is scaled per comparator ellipse: a geometric near-tangency of
/// `~δ` translates to an implicit-value deviation of `~2δ/max(a,b)`, so we
/// pick `eps_l = 2·BOUNDARY_COINCIDENCE_GEOM_TOL / max(a_l, b_l)` (clamped
/// against floating-point noise) so the threshold corresponds to a roughly
/// constant geometric gap regardless of ellipse scale.
fn arc_midpoint_owned_by_j(j: usize, mid: &Point, indices: &[usize], ellipses: &[Ellipse]) -> bool {
    for &l in indices {
        if l == j {
            continue;
        }
        let el = &ellipses[l];
        let scale = el.semi_major.max(el.semi_minor);
        let eps = (2.0 * BOUNDARY_COINCIDENCE_GEOM_TOL / scale).max(IMPLICIT_VALUE_FP_TOL);
        let v = ellipse_implicit_value(mid, el);
        if v > 1.0 + eps {
            return false;
        }
        if (v - 1.0).abs() <= eps && l < j {
            return false;
        }
    }
    true
}

/// Geometric tolerance (in world-frame distance units) for treating two
/// shape boundaries as coincident at a probe point. Converted to a per-shape
/// implicit-value tolerance inside the tiebreaker.
const BOUNDARY_COINCIDENCE_GEOM_TOL: f64 = 1e-7;

/// Floor for the per-shape implicit-value tolerance, guarding against
/// floating-point noise on huge ellipses where the geometric tolerance would
/// otherwise map to a sub-ulp implicit-value threshold.
const IMPLICIT_VALUE_FP_TOL: f64 = 1e-12;

/// Build the CCW boundary arc list for an ellipse-mask region. Mirrors the
/// circle-side `region_boundary_arcs`, working per-ellipse: for each ellipse
/// in the mask, finds IPs on its boundary that sit inside every other mask
/// ellipse, sorts by canonical parameter `t`, and emits an arc for every
/// inter-IP segment whose midpoint sits inside every other mask ellipse.
pub(crate) fn region_boundary_arcs_ellipse(
    mask: RegionMask,
    ellipses: &[Ellipse],
    intersections: &[crate::geometry::diagram::IntersectionPoint],
    n_sets: usize,
) -> Vec<EllipseArc> {
    use crate::geometry::diagram::{adopters_to_mask, mask_to_indices};
    let count = mask.count_ones();
    if count == 0 {
        return Vec::new();
    }
    if count == 1 {
        let idx = mask.trailing_zeros() as usize;
        return vec![EllipseArc {
            ellipse: idx,
            t_start: 0.0,
            delta_t: 2.0 * PI,
        }];
    }

    let indices = mask_to_indices(mask, n_sets);
    let mut arcs: Vec<EllipseArc> = Vec::new();
    let two_pi = 2.0 * PI;
    let arc_eps = 1e-9;

    for &j in &indices {
        let ej = &ellipses[j];
        // IPs on ∂E_j (j is one of the parents) that sit inside every other
        // mask ellipse (mask ⊆ adopters).
        let mut j_ts: Vec<f64> = intersections
            .iter()
            .filter(|ip| {
                let (p1, p2) = ip.parents();
                if p1 != j && p2 != j {
                    return false;
                }
                let am = adopters_to_mask(ip.adopters());
                (mask & am) == mask
            })
            .map(|ip| ellipse_param_for_point(ej, ip.point()))
            .collect();

        if j_ts.is_empty() {
            // E_j has no IPs in the region: either it's wholly inside every
            // other mask ellipse (full-circle arc) or it's outside (no arc).
            // The probe lives on ∂E_j; if it also lies on ∂E_l for some
            // l ≠ j, the index tiebreaker hands the arc to the smallest-index
            // ellipse so identical / coincident-boundary ellipses don't all
            // emit duplicate full-circle arcs.
            let probe = ellipse_point_at(ej, 0.0);
            if arc_midpoint_owned_by_j(j, &probe, &indices, ellipses) {
                arcs.push(EllipseArc {
                    ellipse: j,
                    t_start: 0.0,
                    delta_t: two_pi,
                });
            }
            continue;
        }

        j_ts.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));
        let m = j_ts.len();
        for k in 0..m {
            let t_a = j_ts[k];
            let t_b = j_ts[(k + 1) % m];
            let mut delta = t_b - t_a;
            if delta <= 0.0 {
                delta += two_pi;
            }
            if delta < arc_eps || delta > two_pi - arc_eps {
                continue;
            }
            let t_mid = t_a + delta * 0.5;
            let mid = ellipse_point_at(ej, t_mid);
            if !arc_midpoint_owned_by_j(j, &mid, &indices, ellipses) {
                continue;
            }
            arcs.push(EllipseArc {
                ellipse: j,
                t_start: t_a,
                delta_t: delta,
            });
        }
    }

    arcs
}

/// Area of a region from its CCW boundary arc list, via
/// `A = (1/2) ∮ (x dy − y dx)`. For arc on ellipse k:
/// ```text
/// (x dy − y dx) dt =
///     [ x_c·(b cos φ cos t − a sin φ sin t)
///     + y_c·(a cos φ sin t + b sin φ cos t)
///     + a·b ] dt
/// ```
pub(crate) fn area_from_boundary_arcs_ellipse(arcs: &[EllipseArc], ellipses: &[Ellipse]) -> f64 {
    let mut total = 0.0;
    for arc in arcs {
        let e = &ellipses[arc.ellipse];
        let xc = e.center.x();
        let yc = e.center.y();
        let a = e.semi_major;
        let b = e.semi_minor;
        let cphi = e.rotation.cos();
        let sphi = e.rotation.sin();
        let t_a = arc.t_start;
        let t_b = t_a + arc.delta_t;
        let s_b_minus_s_a = t_b.sin() - t_a.sin();
        let c_b_minus_c_a = t_b.cos() - t_a.cos();
        // ∫(b cos φ cos t − a sin φ sin t) dt = b cos φ (sin t_b − sin t_a)
        //                                       + a sin φ (cos t_b − cos t_a)
        let int_x = b * cphi * s_b_minus_s_a + a * sphi * c_b_minus_c_a;
        // ∫(a cos φ sin t + b sin φ cos t) dt = −a cos φ (cos t_b − cos t_a)
        //                                       + b sin φ (sin t_b − sin t_a)
        let int_y = -a * cphi * c_b_minus_c_a + b * sphi * s_b_minus_s_a;
        total += xc * int_x + yc * int_y + a * b * arc.delta_t;
    }
    0.5 * total
}

/// Accumulate the gradient of a region's overlapping area into `grad`, where
/// `grad` is a length-`5 · n_sets` flat vector laid out as
/// `[x₀, y₀, u₀, v₀, φ₀, x₁, y₁, …]` with `u = ln(a)`, `v = ln(b)`. Each arc
/// on ellipse k contributes (via boundary velocity `dA/dθ = ∮ (v·n) ds`,
/// then chained through `a = exp(u)`, `b = exp(v)`):
///
/// ```text
/// ∂/∂x_c += b cos φ (sin t_b − sin t_a) + a sin φ (cos t_b − cos t_a)
/// ∂/∂y_c += −a cos φ (cos t_b − cos t_a) + b sin φ (sin t_b − sin t_a)
/// ∂/∂u   += a · [(b/2) Δt + (b/4) (sin 2t_b − sin 2t_a)]
/// ∂/∂v   += b · [(a/2) Δt − (a/4) (sin 2t_b − sin 2t_a)]
/// ∂/∂φ   += (a² − b²)/4 · (cos 2t_a − cos 2t_b)
/// ```
pub(crate) fn accumulate_region_overlap_gradient_ellipse(
    arcs: &[EllipseArc],
    ellipses: &[Ellipse],
    grad: &mut [f64],
) {
    for arc in arcs {
        let e = &ellipses[arc.ellipse];
        let a = e.semi_major;
        let b = e.semi_minor;
        let cphi = e.rotation.cos();
        let sphi = e.rotation.sin();
        let t_a = arc.t_start;
        let t_b = t_a + arc.delta_t;
        let s_a = t_a.sin();
        let s_b = t_b.sin();
        let c_a = t_a.cos();
        let c_b = t_b.cos();
        let s2_a = (2.0 * t_a).sin();
        let s2_b = (2.0 * t_b).sin();
        let c2_a = (2.0 * t_a).cos();
        let c2_b = (2.0 * t_b).cos();
        let off = arc.ellipse * 5;
        grad[off] += b * cphi * (s_b - s_a) + a * sphi * (c_b - c_a);
        grad[off + 1] += -a * cphi * (c_b - c_a) + b * sphi * (s_b - s_a);
        // ∂A/∂u = ∂A/∂a · da/du = (∂A/∂a) · a
        grad[off + 2] += a * (0.5 * b * arc.delta_t + 0.25 * b * (s2_b - s2_a));
        // ∂A/∂v = ∂A/∂b · db/dv = (∂A/∂b) · b
        grad[off + 3] += b * (0.5 * a * arc.delta_t - 0.25 * a * (s2_b - s2_a));
        grad[off + 4] += 0.25 * (a * a - b * b) * (c2_a - c2_b);
    }
}

/// Collect IPs between every pair of ellipses, with adopters computed via the
/// implicit-form check used by [`Ellipse::compute_exclusive_regions`] (small
/// tolerance to capture boundary points).
pub(crate) fn collect_intersections_ellipse(
    ellipses: &[Ellipse],
) -> Vec<crate::geometry::diagram::IntersectionPoint> {
    use crate::geometry::diagram::IntersectionPoint;
    let n_sets = ellipses.len();
    let mut intersections: Vec<IntersectionPoint> = Vec::new();
    for i in 0..n_sets {
        for j in (i + 1)..n_sets {
            let pts = ellipses[i].intersection_points(&ellipses[j]);
            for p in pts {
                let adopters: Vec<usize> = (0..n_sets)
                    .filter(|&k| {
                        let local = p.to_ellipse_frame(&ellipses[k]);
                        let val = (local.x() / ellipses[k].semi_major).powi(2)
                            + (local.y() / ellipses[k].semi_minor).powi(2);
                        val <= 1.0 + 1e-9
                    })
                    .collect();
                intersections.push(IntersectionPoint::new(p, (i, j), adopters));
            }
        }
    }
    intersections
}

/// Compute the per-mask exclusive areas of `ellipses` clipped to an
/// axis-aligned `container` rectangle, including the all-zeros region (mask
/// `0`) representing `container ∖ ⋃ ellipses` (the complement target).
///
/// Mirrors `circle::compute_exclusive_regions_clipped`. For each overlap
/// region `R_m = ⋂_{i ∈ m} E_i`, computes `area(R_m ∩ container)` via Green's
/// theorem with the boundary decomposed into:
///
/// - sub-arcs of `∂R_m` (arcs on each `∂E_i`) that lie inside the container, and
/// - sub-segments of `∂container` (axis-aligned edges) that lie inside `R_m`.
///
/// Mask `0` is seeded with `area(container)` so the inclusion-exclusion pass
/// produces `container.area − area(⋃ ellipses ∩ container)` for free.
pub fn compute_exclusive_regions_clipped_ellipse(
    ellipses: &[Ellipse],
    container: &Rectangle,
) -> std::collections::HashMap<RegionMask, f64> {
    use crate::geometry::diagram::{discover_regions, to_exclusive_areas};

    let n_sets = ellipses.len();
    let intersections = collect_intersections_ellipse(ellipses);
    let regions = discover_regions(ellipses, &intersections, n_sets);

    let mut overlapping_areas = std::collections::HashMap::new();
    overlapping_areas.insert(0, container.area());

    for &mask in &regions {
        let arcs = region_boundary_arcs_ellipse(mask, ellipses, &intersections, n_sets);
        let area = area_from_clipped_arcs_ellipse(&arcs, ellipses, container, mask, n_sets);
        overlapping_areas.insert(mask, area);
    }

    to_exclusive_areas(&overlapping_areas)
}

/// Gradient-aware companion to [`compute_exclusive_regions_clipped_ellipse`].
///
/// Returns `(exclusive_areas, exclusive_grads)` where each gradient vector has
/// length `n_sets · 5 + 4` and is laid out as
/// `[x₀, y₀, u₀, v₀, φ₀, x₁, …, x_c, y_c, u, v]` with `u_k = ln(a_k)`,
/// `v_k = ln(b_k)`. The trailing four entries are the container's optimizer
/// encoding (`u = ln(w·h)`, `v = ln(w/h)`).
///
/// The shape-param gradient comes from the same boundary-velocity identity
/// `dA/dθ = ∮_∂R (v_θ · n) ds` as the unclipped path
/// (`accumulate_region_overlap_gradient_ellipse`), restricted to inside-container
/// sub-arcs. Box-edge sub-segments contribute to the container-param block
/// only — the container moves rigidly, identical to the circle case.
/// Endpoints where an arc meets a box edge cancel pairwise (same point, same
/// `(v · n)`, opposite ds direction across the two boundary pieces), so this
/// builds a consistent gradient without tracking trim-point motion.
///
/// Mask `0` (the complement) is seeded with `container.area()` and a gradient
/// `∂(w·h)/∂(x_c, y_c, u, v) = [0, 0, w·h, 0]` so inclusion-exclusion produces
/// `complement = container.area − area(⋃ ellipses ∩ container)` along with its
/// matching gradient.
pub(crate) fn compute_exclusive_regions_clipped_ellipse_with_gradient(
    ellipses: &[Ellipse],
    container: &Rectangle,
) -> (
    std::collections::HashMap<RegionMask, f64>,
    std::collections::HashMap<RegionMask, Vec<f64>>,
) {
    use crate::geometry::diagram::{discover_regions, to_exclusive_areas_and_gradients};

    let n_sets = ellipses.len();
    let n_params = n_sets * 5 + 4;
    let intersections = collect_intersections_ellipse(ellipses);
    let regions = discover_regions(ellipses, &intersections, n_sets);

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

    // Seed mask 0 with `container.area()` and its gradient. ∂(w·h)/∂u = w·h
    // (chain through `u = ln(w·h)`); other partials are zero.
    let container_area = container.area();
    overlapping_areas.insert(0, container_area);
    let mut zero_grad = vec![0.0; n_params];
    zero_grad[n_sets * 5 + 2] = container_area;
    overlapping_grads.insert(0, zero_grad);

    for &mask in &regions {
        let arcs = region_boundary_arcs_ellipse(mask, ellipses, &intersections, n_sets);
        let mut grad = vec![0.0; n_params];
        let area = area_and_gradient_from_clipped_arcs_ellipse(
            &arcs, ellipses, container, mask, n_sets, &mut grad,
        );
        overlapping_areas.insert(mask, area);
        overlapping_grads.insert(mask, grad);
    }

    to_exclusive_areas_and_gradients(&overlapping_areas, &overlapping_grads, n_params)
}

/// Compute `area(R_m ∩ container)` from the CCW ellipse-boundary arcs of
/// `R_m` and the axis-aligned `container`. Mirrors `circle::area_from_clipped_arcs`:
/// arc contributions split at box-edge crossings, plus container-edge
/// contributions (intersection of per-ellipse interior-x / interior-y ranges).
fn area_from_clipped_arcs_ellipse(
    arcs: &[EllipseArc],
    ellipses: &[Ellipse],
    container: &Rectangle,
    mask: RegionMask,
    n_sets: usize,
) -> f64 {
    use crate::geometry::diagram::mask_to_indices;

    let (x_min, x_max, y_min, y_max) = container.bounds();
    if x_max <= x_min || y_max <= y_min {
        return 0.0;
    }

    let mut total = 0.0; // ∫ (x dy − y dx); area = 0.5 · total.

    // 1. Arc contributions.
    for arc in arcs {
        total += clipped_ellipse_arc_integral(arc, ellipses, x_min, x_max, y_min, y_max);
    }

    // 2. Container-edge contributions, each oriented CCW around the container
    // interior. For each edge we find the inside-`R_m` sub-interval (the
    // intersection of the per-ellipse projections onto that edge), then add
    // the closed-form line integral.
    let indices = mask_to_indices(mask, n_sets);

    // Bottom edge: y = y_min, walk x from x_min → x_max. ∫(−y dx) = −y_min·Δx.
    if let Some((a, b)) =
        ellipse_horizontal_edge_inside_interval(y_min, x_min, x_max, &indices, ellipses)
    {
        total += -y_min * (b - a);
    }
    // Right edge: x = x_max, walk y from y_min → y_max. ∫(x dy) = x_max·Δy.
    if let Some((a, b)) =
        ellipse_vertical_edge_inside_interval(x_max, y_min, y_max, &indices, ellipses)
    {
        total += x_max * (b - a);
    }
    // Top edge: y = y_max, walk x from x_max → x_min. ∫(−y dx) = +y_max·Δx
    // because dx is negative when walking right→left.
    if let Some((a, b)) =
        ellipse_horizontal_edge_inside_interval(y_max, x_min, x_max, &indices, ellipses)
    {
        total += y_max * (b - a);
    }
    // Left edge: x = x_min, walk y from y_max → y_min. ∫(x dy) = −x_min·Δy
    // because dy is negative when walking top→bottom.
    if let Some((a, b)) =
        ellipse_vertical_edge_inside_interval(x_min, y_min, y_max, &indices, ellipses)
    {
        total += -x_min * (b - a);
    }

    0.5 * total
}

/// Gradient-aware companion to [`area_from_clipped_arcs_ellipse`]. Computes
/// `area(R_m ∩ container)` and accumulates `∂area / ∂θ` into `grad`
/// (length `n_sets · 5 + 4`, layout `[x₀, y₀, u₀, v₀, φ₀, …, x_c, y_c, u, v]`
/// with `u_k = ln(a_k)`, `v_k = ln(b_k)` and the trailing four entries the
/// container's optimizer encoding).
///
/// The arc and box-edge decompositions are identical to the area-only path;
/// the gradient is produced inline using the per-piece boundary-velocity
/// formulas — see the function docs of
/// [`compute_exclusive_regions_clipped_ellipse_with_gradient`] for the
/// derivation summary.
fn area_and_gradient_from_clipped_arcs_ellipse(
    arcs: &[EllipseArc],
    ellipses: &[Ellipse],
    container: &Rectangle,
    mask: RegionMask,
    n_sets: usize,
    grad: &mut [f64],
) -> f64 {
    use crate::geometry::diagram::mask_to_indices;

    let (x_min, x_max, y_min, y_max) = container.bounds();
    if x_max <= x_min || y_max <= y_min {
        return 0.0;
    }

    let w = x_max - x_min;
    let h = y_max - y_min;
    let container_off = n_sets * 5;

    let mut total = 0.0; // ∫ (x dy − y dx); area = 0.5 · total.

    // 1. Arc contributions — area + shape-param gradient.
    for arc in arcs {
        total += clipped_ellipse_arc_integral_with_gradient(
            arc, ellipses, x_min, x_max, y_min, y_max, grad,
        );
    }

    // 2. Container-edge contributions. The container moves rigidly, so the
    // box-edge gradient identities are identical to the circle case — see
    // `circle::area_and_gradient_from_clipped_arcs` for the derivation. The
    // only ellipse-specific piece is the inside-region sub-interval.
    let indices = mask_to_indices(mask, n_sets);

    if let Some((a, b)) =
        ellipse_horizontal_edge_inside_interval(y_min, x_min, x_max, &indices, ellipses)
    {
        let l = b - a;
        // Area: bottom edge walked left→right; ∫(−y dx) = −y_min·L.
        total += -y_min * l;
        // Gradient: outward normal (0, −1); v·n = −∂y_min/∂θ.
        // ∂y_min/∂(x_c, y_c, u, v) = (0, 1, −h/4, h/4)
        grad[container_off + 1] -= l;
        grad[container_off + 2] += (h / 4.0) * l;
        grad[container_off + 3] -= (h / 4.0) * l;
    }
    if let Some((a, b)) =
        ellipse_vertical_edge_inside_interval(x_max, y_min, y_max, &indices, ellipses)
    {
        let l = b - a;
        total += x_max * l;
        // Right edge: outward normal (+1, 0); v·n = ∂x_max/∂θ.
        // ∂x_max/∂(x_c, y_c, u, v) = (1, 0, w/4, w/4)
        grad[container_off] += l;
        grad[container_off + 2] += (w / 4.0) * l;
        grad[container_off + 3] += (w / 4.0) * l;
    }
    if let Some((a, b)) =
        ellipse_horizontal_edge_inside_interval(y_max, x_min, x_max, &indices, ellipses)
    {
        let l = b - a;
        total += y_max * l;
        // Top edge: outward normal (0, +1); v·n = ∂y_max/∂θ.
        // ∂y_max/∂(x_c, y_c, u, v) = (0, 1, h/4, −h/4)
        grad[container_off + 1] += l;
        grad[container_off + 2] += (h / 4.0) * l;
        grad[container_off + 3] -= (h / 4.0) * l;
    }
    if let Some((a, b)) =
        ellipse_vertical_edge_inside_interval(x_min, y_min, y_max, &indices, ellipses)
    {
        let l = b - a;
        total += -x_min * l;
        // Left edge: outward normal (−1, 0); v·n = −∂x_min/∂θ.
        // ∂x_min/∂(x_c, y_c, u, v) = (1, 0, −w/4, −w/4)
        grad[container_off] -= l;
        grad[container_off + 2] += (w / 4.0) * l;
        grad[container_off + 3] += (w / 4.0) * l;
    }

    0.5 * total
}

/// Sum the Green's-theorem arc integral over the inside-container sub-arcs of
/// `arc`, splitting at crossings with the four box edges. Mirrors
/// `circle::clipped_arc_integral`, with crossings found by solving
/// `A cos t + B sin t = C` for the rotated-ellipse parameterisation.
fn clipped_ellipse_arc_integral(
    arc: &EllipseArc,
    ellipses: &[Ellipse],
    x_min: f64,
    x_max: f64,
    y_min: f64,
    y_max: f64,
) -> f64 {
    let e = &ellipses[arc.ellipse];
    let xc = e.center.x();
    let yc = e.center.y();
    let a = e.semi_major;
    let b = e.semi_minor;
    let cphi = e.rotation.cos();
    let sphi = e.rotation.sin();

    let t_a = arc.t_start;
    let t_b = t_a + arc.delta_t; // continuous (not wrapped)

    let mut crossings: Vec<f64> = Vec::with_capacity(8);

    // Vertical edges x = X: a cos φ · cos t + (−b sin φ) · sin t = X − xc.
    for &x_edge in &[x_min, x_max] {
        push_axis_crossings(&mut crossings, a * cphi, -b * sphi, x_edge - xc, t_a, t_b);
    }
    // Horizontal edges y = Y: a sin φ · cos t + b cos φ · sin t = Y − yc.
    for &y_edge in &[y_min, y_max] {
        push_axis_crossings(&mut crossings, a * sphi, b * cphi, y_edge - yc, t_a, t_b);
    }

    let mut breaks: Vec<f64> = Vec::with_capacity(crossings.len() + 2);
    breaks.push(t_a);
    breaks.extend(crossings);
    breaks.push(t_b);
    breaks.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));

    let mut sum = 0.0;
    for w in breaks.windows(2) {
        let ta = w[0];
        let tb = w[1];
        let delta = tb - ta;
        if delta <= 1e-12 {
            continue;
        }
        let mid = 0.5 * (ta + tb);
        let mid_pt = ellipse_point_at(e, mid);
        let mx = mid_pt.x();
        let my = mid_pt.y();
        // Generous geometric tolerance to keep boundary-aligned midpoints
        // (e.g. arcs that exit and immediately re-enter via a corner) on the
        // "inside" side. Mismatched signs here would lose a slice of area.
        let tol = 1e-9;
        let inside =
            mx >= x_min - tol && mx <= x_max + tol && my >= y_min - tol && my <= y_max + tol;
        if !inside {
            continue;
        }

        // Integrand of `x dy − y dx` for the rotated ellipse, integrated
        // analytically from `ta` to `tb` (matches `area_from_boundary_arcs_ellipse`).
        let s_diff = tb.sin() - ta.sin();
        let c_diff = tb.cos() - ta.cos();
        let int_x = b * cphi * s_diff + a * sphi * c_diff;
        let int_y = -a * cphi * c_diff + b * sphi * s_diff;
        sum += xc * int_x + yc * int_y + a * b * delta;
    }
    sum
}

/// Gradient-aware companion to [`clipped_ellipse_arc_integral`]. Returns the
/// same area-integrand contribution to `2A` while also accumulating shape-param
/// gradient entries for the owning ellipse into `grad` (length
/// `≥ 5 · n_sets + 4`, layout `[x₀, y₀, u₀, v₀, φ₀, …, x_c, y_c, u, v]`).
///
/// Per inside-container sub-arc on ellipse `k` with continuous endpoints
/// `(t_a, t_b)`, `t_b > t_a`, the boundary-velocity identity (derived in
/// [`accumulate_region_overlap_gradient_ellipse`]) gives:
///
/// ```text
/// ∂A/∂x_k += b cos φ (sin t_b − sin t_a) + a sin φ (cos t_b − cos t_a)
/// ∂A/∂y_k += −a cos φ (cos t_b − cos t_a) + b sin φ (sin t_b − sin t_a)
/// ∂A/∂u_k += a · [(b/2) Δt + (b/4) (sin 2t_b − sin 2t_a)]
/// ∂A/∂v_k += b · [(a/2) Δt − (a/4) (sin 2t_b − sin 2t_a)]
/// ∂A/∂φ_k += (a² − b²)/4 · (cos 2t_a − cos 2t_b)
/// ```
///
/// Trim-point motion at box-edge endpoints is captured by the matching
/// box-edge integral (same point, opposite-sign `ds`), so the per-sub-arc
/// formula is self-contained.
fn clipped_ellipse_arc_integral_with_gradient(
    arc: &EllipseArc,
    ellipses: &[Ellipse],
    x_min: f64,
    x_max: f64,
    y_min: f64,
    y_max: f64,
    grad: &mut [f64],
) -> f64 {
    let e = &ellipses[arc.ellipse];
    let xc = e.center.x();
    let yc = e.center.y();
    let a = e.semi_major;
    let b = e.semi_minor;
    let cphi = e.rotation.cos();
    let sphi = e.rotation.sin();

    let t_a = arc.t_start;
    let t_b = t_a + arc.delta_t;

    let mut crossings: Vec<f64> = Vec::with_capacity(8);

    for &x_edge in &[x_min, x_max] {
        push_axis_crossings(&mut crossings, a * cphi, -b * sphi, x_edge - xc, t_a, t_b);
    }
    for &y_edge in &[y_min, y_max] {
        push_axis_crossings(&mut crossings, a * sphi, b * cphi, y_edge - yc, t_a, t_b);
    }

    let mut breaks: Vec<f64> = Vec::with_capacity(crossings.len() + 2);
    breaks.push(t_a);
    breaks.extend(crossings);
    breaks.push(t_b);
    breaks.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));

    let off = arc.ellipse * 5;
    let mut sum = 0.0;
    for w in breaks.windows(2) {
        let ta = w[0];
        let tb = w[1];
        let delta = tb - ta;
        if delta <= 1e-12 {
            continue;
        }
        let mid = 0.5 * (ta + tb);
        let mid_pt = ellipse_point_at(e, mid);
        let mx = mid_pt.x();
        let my = mid_pt.y();
        let tol = 1e-9;
        let inside =
            mx >= x_min - tol && mx <= x_max + tol && my >= y_min - tol && my <= y_max + tol;
        if !inside {
            continue;
        }

        let s_a = ta.sin();
        let s_b = tb.sin();
        let c_a = ta.cos();
        let c_b = tb.cos();
        let s_diff = s_b - s_a;
        let c_diff = c_b - c_a;
        let int_x = b * cphi * s_diff + a * sphi * c_diff;
        let int_y = -a * cphi * c_diff + b * sphi * s_diff;
        // 2A integrand (matches the area-only path).
        sum += xc * int_x + yc * int_y + a * b * delta;

        // Boundary-velocity gradient on ellipse k. Sub-arcs inherit the
        // parent arc's CCW orientation (delta > 0), so no sign flip needed.
        let s2_a = (2.0 * ta).sin();
        let s2_b = (2.0 * tb).sin();
        let c2_a = (2.0 * ta).cos();
        let c2_b = (2.0 * tb).cos();
        grad[off] += int_x;
        grad[off + 1] += int_y;
        // ∂A/∂u = ∂A/∂a · da/du = (∂A/∂a) · a
        grad[off + 2] += a * (0.5 * b * delta + 0.25 * b * (s2_b - s2_a));
        // ∂A/∂v = ∂A/∂b · db/dv = (∂A/∂b) · b
        grad[off + 3] += b * (0.5 * a * delta - 0.25 * a * (s2_b - s2_a));
        grad[off + 4] += 0.25 * (a * a - b * b) * (c2_a - c2_b);
    }
    sum
}

/// Solve `A cos t + B sin t = C` and push every periodic copy of each
/// solution that lands strictly inside `(t_a, t_b)` into `out`. Skips when
/// `R = √(A² + B²)` is zero or `|C| > R`.
fn push_axis_crossings(out: &mut Vec<f64>, big_a: f64, big_b: f64, big_c: f64, t_a: f64, t_b: f64) {
    let r_sq = big_a * big_a + big_b * big_b;
    if r_sq <= 0.0 {
        return;
    }
    let r = r_sq.sqrt();
    let ratio = big_c / r;
    if ratio.abs() > 1.0 {
        return;
    }
    // R cos(t − δ) = C, with δ = atan2(B, A).
    let delta = big_b.atan2(big_a);
    let alpha = ratio.clamp(-1.0, 1.0).acos();
    collect_periodic_in_range_ellipse(out, delta + alpha, t_a, t_b);
    collect_periodic_in_range_ellipse(out, delta - alpha, t_a, t_b);
}

/// For a candidate solution `cand`, push every periodic copy `cand + 2πk`
/// strictly inside `(t_a, t_b)` into `out`. Arc spans can reach 2π so at most
/// one copy per shift offset is in range; iterating `k ∈ {-1, 0, 1}` covers
/// every case for `delta_t ≤ 2π`.
fn collect_periodic_in_range_ellipse(out: &mut Vec<f64>, cand: f64, t_a: f64, t_b: f64) {
    let two_pi = 2.0 * PI;
    let tol = 1e-12;
    for k in -1..=1 {
        let p = cand + (k as f64) * two_pi;
        if p > t_a + tol && p < t_b - tol {
            out.push(p);
        }
    }
}

/// Inside-`R_m` sub-interval of a horizontal edge `y = const` over
/// `x ∈ [x_min, x_max]`, where `R_m = ⋂_{i ∈ indices} E_i`. Returns `None` if
/// the edge does not touch every ellipse in `indices`. With `indices` empty,
/// the whole `[x_min, x_max]` is "inside" (intersection over an empty
/// constraint set is the universe), which is exactly what mask `0` (the
/// complement) wants.
///
/// For each ellipse the inside-x range at fixed `y` is found by substituting
/// `dy = y − y_c` into the rotated-ellipse implicit form
/// `A·dx² + B·dx + C = 0` and taking the two real roots.
fn ellipse_horizontal_edge_inside_interval(
    y: f64,
    x_min: f64,
    x_max: f64,
    indices: &[usize],
    ellipses: &[Ellipse],
) -> Option<(f64, f64)> {
    let mut lo = x_min;
    let mut hi = x_max;
    for &i in indices {
        let e = &ellipses[i];
        let dy = y - e.center.y();
        let cphi = e.rotation.cos();
        let sphi = e.rotation.sin();
        let inv_a2 = 1.0 / (e.semi_major * e.semi_major);
        let inv_b2 = 1.0 / (e.semi_minor * e.semi_minor);
        let big_a = cphi * cphi * inv_a2 + sphi * sphi * inv_b2;
        let big_b = 2.0 * dy * sphi * cphi * (inv_a2 - inv_b2);
        let big_c = dy * dy * (sphi * sphi * inv_a2 + cphi * cphi * inv_b2) - 1.0;
        let disc = big_b * big_b - 4.0 * big_a * big_c;
        if disc < 0.0 {
            return None;
        }
        let sqrt_disc = disc.sqrt();
        let dx_lo = (-big_b - sqrt_disc) / (2.0 * big_a);
        let dx_hi = (-big_b + sqrt_disc) / (2.0 * big_a);
        let disk_lo = e.center.x() + dx_lo;
        let disk_hi = e.center.x() + dx_hi;
        if disk_lo > lo {
            lo = disk_lo;
        }
        if disk_hi < hi {
            hi = disk_hi;
        }
        if hi <= lo {
            return None;
        }
    }
    if hi > lo {
        Some((lo, hi))
    } else {
        None
    }
}

/// Inside-`R_m` sub-interval of a vertical edge `x = const` over
/// `y ∈ [y_min, y_max]`. Mirror of [`ellipse_horizontal_edge_inside_interval`].
fn ellipse_vertical_edge_inside_interval(
    x: f64,
    y_min: f64,
    y_max: f64,
    indices: &[usize],
    ellipses: &[Ellipse],
) -> Option<(f64, f64)> {
    let mut lo = y_min;
    let mut hi = y_max;
    for &i in indices {
        let e = &ellipses[i];
        let dx = x - e.center.x();
        let cphi = e.rotation.cos();
        let sphi = e.rotation.sin();
        let inv_a2 = 1.0 / (e.semi_major * e.semi_major);
        let inv_b2 = 1.0 / (e.semi_minor * e.semi_minor);
        let big_a = sphi * sphi * inv_a2 + cphi * cphi * inv_b2;
        let big_b = 2.0 * dx * sphi * cphi * (inv_a2 - inv_b2);
        let big_c = dx * dx * (cphi * cphi * inv_a2 + sphi * sphi * inv_b2) - 1.0;
        let disc = big_b * big_b - 4.0 * big_a * big_c;
        if disc < 0.0 {
            return None;
        }
        let sqrt_disc = disc.sqrt();
        let dy_lo = (-big_b - sqrt_disc) / (2.0 * big_a);
        let dy_hi = (-big_b + sqrt_disc) / (2.0 * big_a);
        let disk_lo = e.center.y() + dy_lo;
        let disk_hi = e.center.y() + dy_hi;
        if disk_lo > lo {
            lo = disk_lo;
        }
        if disk_hi < hi {
            hi = disk_hi;
        }
        if hi <= lo {
            return None;
        }
    }
    if hi > lo {
        Some((lo, hi))
    } else {
        None
    }
}

impl DiagramShape for Ellipse {
    fn compute_exclusive_regions(shapes: &[Self]) -> std::collections::HashMap<RegionMask, f64> {
        use crate::geometry::diagram::{discover_regions, to_exclusive_areas};
        use std::collections::HashMap;

        let n_sets = shapes.len();
        let intersections = collect_intersections_ellipse(shapes);
        // Sparse discovery: walk only regions that can geometrically be
        // non-empty. Full enumeration over `1..=(1<<n_sets)-1` is intractable
        // beyond a handful of sets — at n=17 it's 131,071 boundary-arc passes
        // per loss eval, which is the original eulerr issue #89 problem.
        let regions = discover_regions(shapes, &intersections, n_sets);

        let mut overlapping_areas: HashMap<RegionMask, f64> = HashMap::new();
        for &mask in &regions {
            let arcs = region_boundary_arcs_ellipse(mask, shapes, &intersections, n_sets);
            overlapping_areas.insert(mask, area_from_boundary_arcs_ellipse(&arcs, shapes));
        }

        to_exclusive_areas(&overlapping_areas)
    }

    fn compute_exclusive_regions_with_gradient(
        shapes: &[Self],
    ) -> Option<crate::geometry::traits::ExclusiveRegionsAndGradient> {
        use crate::geometry::diagram::{discover_regions, to_exclusive_areas_and_gradients};
        use std::collections::HashMap;

        let n_sets = shapes.len();
        let n_params = n_sets * 5;
        let intersections = collect_intersections_ellipse(shapes);
        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 arcs = region_boundary_arcs_ellipse(mask, shapes, &intersections, n_sets);
            overlapping_areas.insert(mask, area_from_boundary_arcs_ellipse(&arcs, shapes));
            let mut grad = vec![0.0; n_params];
            accumulate_region_overlap_gradient_ellipse(&arcs, shapes, &mut grad);
            overlapping_grads.insert(mask, grad);
        }
        Some(to_exclusive_areas_and_gradients(
            &overlapping_areas,
            &overlapping_grads,
            n_params,
        ))
    }

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

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

    fn optimizer_params_from_circle(x: f64, y: f64, radius: f64) -> Vec<f64> {
        // Semi-axes are stored in log space (`u = ln(a)`, `v = ln(b)`) so the
        // unbounded LM solver stays on the positive-axis manifold without a
        // clamp; `from_optimizer_params` exponentiates back. For a circle,
        // `a = b = r`, so `u = v = ln(r)`. Rotation is unconstrained
        // (periodic).
        let log_radius = radius.ln();
        vec![x, y, log_radius, log_radius, 0.0]
    }

    fn mds_target_distance(
        area_i: f64,
        area_j: f64,
        target_overlap: f64,
    ) -> Result<f64, crate::error::DiagramError> {
        // The MDS warm-start treats every ellipse as a circle of equal area
        // (`r = sqrt(area/π)`), so the distance inversion matches Circle's.
        Circle::mds_target_distance(area_i, area_j, target_overlap)
    }

    fn n_params() -> usize {
        5 // x, y, ln(semi_major), ln(semi_minor), rotation
    }

    fn from_params(params: &[f64]) -> Self {
        assert_eq!(
            params.len(),
            5,
            "Ellipse requires 5 parameters: x, y, semi_major, semi_minor, rotation"
        );
        Ellipse::new(
            Point::new(params[0], params[1]),
            params[2],
            params[3],
            params[4],
        )
    }

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

    fn from_optimizer_params(params: &[f64]) -> Self {
        assert_eq!(
            params.len(),
            5,
            "Ellipse requires 5 parameters: x, y, ln(semi_major), ln(semi_minor), rotation"
        );

        // Floor at `f64::MIN_POSITIVE` so that very negative `u` / `v`
        // values (which `.exp()` underflows to 0) still produce a strictly
        // positive semi-axis — `Ellipse::new` asserts `> 0`. The optimizer
        // should not drive here in practice; this is purely defensive
        // against runaway log-space steps.
        Ellipse::new(
            Point::new(params[0], params[1]),
            params[2].exp().max(f64::MIN_POSITIVE),
            params[3].exp().max(f64::MIN_POSITIVE),
            params[4],
        )
    }

    fn to_optimizer_params(&self) -> Vec<f64> {
        vec![
            self.center.x(),
            self.center.y(),
            self.semi_major.ln(),
            self.semi_minor.ln(),
            self.rotation,
        ]
    }

    fn canonical_venn_layout(n: usize) -> Option<Vec<Self>> {
        let params: &[CanonicalVennParams] = match n {
            1 => &VENN_N1,
            2 => &VENN_N2,
            3 => &VENN_N3,
            4 => &VENN_N4,
            5 => &VENN_N5,
            _ => return None,
        };
        Some(
            params
                .iter()
                .map(|&(h, k, a, b, phi)| Ellipse::new(Point::new(h, k), a, b, phi))
                .collect(),
        )
    }
}

/// Ellipse parameters as `(h, k, a, b, phi)` for the canonical Venn arrangements.
type CanonicalVennParams = (f64, f64, f64, f64, f64);

// Canonical Venn-diagram arrangements taken verbatim from the eulerr R
// package's `venn_spec` data table, which in turn reproduces the published
// constructions:
//
// - n = 1, 2, 3: rotationally symmetric circles.
// - n = 4: Venn (1880) 4-ellipse arrangement; see Wilkinson (2012),
//   *JCGS* "Exact Rotational Symmetry in Venn Diagrams".
// - n = 5: Grünbaum (1975) symmetric 5-ellipse Venn diagram.
const VENN_N1: [CanonicalVennParams; 1] = [(0.0, 0.0, 1.0, 1.0, 0.0)];

const VENN_N2: [CanonicalVennParams; 2] = [(-0.5, 0.0, 1.0, 1.0, 1.0), (0.5, 0.0, 1.0, 1.0, 1.0)];

const VENN_N3: [CanonicalVennParams; 3] = [
    (-0.42, -0.36, 1.05, 1.05, 3.76),
    (0.42, -0.36, 1.05, 1.05, 3.76),
    (0.00, 0.36, 1.05, 1.05, 3.76),
];

const VENN_N4: [CanonicalVennParams; 4] = [
    (-0.8, 0.0, 1.2, 2.0, PI / 4.0),
    (0.8, 0.0, 1.2, 2.0, -PI / 4.0),
    (0.0, 1.0, 1.2, 2.0, PI / 4.0),
    (0.0, 1.0, 1.2, 2.0, -PI / 4.0),
];

const VENN_N5: [CanonicalVennParams; 5] = [
    (0.176, 0.096, 1.0, 0.6, 0.000),
    (-0.037, 0.197, 1.0, 0.6, 1.257),
    (-0.198, 0.026, 1.0, 0.6, 2.513),
    (-0.086, -0.181, 1.0, 0.6, 3.770),
    (0.145, -0.137, 1.0, 0.6, 5.027),
];

impl Polygonize for Ellipse {
    fn polygonize(&self, n_vertices: usize) -> Polygon {
        use std::f64::consts::PI;

        let n_vertices = n_vertices.max(3);
        let mut vertices = Vec::with_capacity(n_vertices);

        let cos_rot = self.rotation.cos();
        let sin_rot = self.rotation.sin();

        for i in 0..n_vertices {
            let angle = 2.0 * PI * (i as f64) / (n_vertices as f64);

            // Point in local ellipse frame
            let local_x = self.semi_major * angle.cos();
            let local_y = self.semi_minor * angle.sin();

            // Rotate and translate to world frame
            let x = self.center.x() + local_x * cos_rot - local_y * sin_rot;
            let y = self.center.y() + local_x * sin_rot + local_y * cos_rot;

            vertices.push(Point::new(x, y));
        }

        Polygon::new(vertices)
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::geometry::shapes::Circle;

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

    #[test]
    fn test_ellipse_try_new_accepts_positive() {
        let e = Ellipse::try_new(Point::new(0.0, 0.0), 4.0, 3.0, 0.5).unwrap();
        assert_eq!(e.semi_major(), 4.0);
        assert_eq!(e.semi_minor(), 3.0);
    }

    #[test]
    fn test_ellipse_try_new_rejects_non_positive_axes() {
        for (a, b, expected_param) in [
            (0.0, 3.0, "semi_major"),
            (-1.0, 3.0, "semi_major"),
            (4.0, 0.0, "semi_minor"),
            (4.0, -2.0, "semi_minor"),
        ] {
            let err = Ellipse::try_new(Point::new(0.0, 0.0), a, b, 0.0).unwrap_err();
            match err {
                crate::error::DiagramError::InvalidShapeParameter { shape, param, .. } => {
                    assert_eq!(shape, "Ellipse");
                    assert_eq!(param, expected_param);
                }
                _ => panic!("expected InvalidShapeParameter"),
            }
        }
    }

    #[test]
    #[should_panic(expected = "Ellipse semi_major must be > 0")]
    fn test_ellipse_new_panics_on_zero_semi_major() {
        let _ = Ellipse::new(Point::new(0.0, 0.0), 0.0, 3.0, 0.0);
    }

    #[test]
    #[should_panic(expected = "Ellipse semi_minor must be > 0")]
    fn test_ellipse_new_panics_on_zero_semi_minor() {
        let _ = Ellipse::new(Point::new(0.0, 0.0), 4.0, 0.0, 0.0);
    }

    #[test]
    fn test_sector_area_between_quarter() {
        let ellipse = Ellipse::new(Point::new(0.0, 0.0), 5.0, 3.0, 0.0);

        // Quarter sector from 0 to π/2
        let sector_between = ellipse.sector_area_between(0.0, PI / 2.0);
        let expected = ellipse.sector_area(PI / 2.0) - ellipse.sector_area(0.0);

        assert!(approx_eq(sector_between, expected));
        assert!(approx_eq(sector_between, ellipse.area() / 4.0));
    }

    #[test]
    fn test_sector_area_between_with_ccw_adjustment() {
        let ellipse = Ellipse::new(Point::new(0.0, 0.0), 4.0, 2.0, 0.0);

        // When theta1 < theta0, it should add 2π to theta1 for CCW ordering
        // Going from 3π/4 to π/4 CCW (wrapping around through 0)
        let theta0 = 3.0 * PI / 4.0;
        let theta1 = PI / 4.0;

        let sector_between = ellipse.sector_area_between(theta0, theta1);

        // This should be equivalent to going from 3π/4 to (π/4 + 2π)
        let expected = ellipse.sector_area(theta1 + 2.0 * PI) - ellipse.sector_area(theta0);

        assert!(
            approx_eq(sector_between, expected),
            "Expected {}, got {}",
            expected,
            sector_between
        );

        // The sector should be positive and cover most of the ellipse
        assert!(sector_between > 0.0);
        assert!(sector_between > ellipse.area() / 2.0);
    }

    #[test]
    fn test_sector_area_between_same_angle() {
        let ellipse = Ellipse::new(Point::new(1.0, 2.0), 3.0, 2.5, 0.3);

        // Same angle should give zero sector area
        let angle = PI / 3.0;
        let sector = ellipse.sector_area_between(angle, angle);

        assert!(approx_eq(sector, 0.0));
    }

    #[test]
    fn test_sector_area_between_full_rotation() {
        let ellipse = Ellipse::new(Point::new(0.0, 0.0), 5.0, 3.0, 0.0);

        // Full rotation: from 0 to 2π
        let sector = ellipse.sector_area_between(0.0, 2.0 * PI);

        assert!(approx_eq(sector, ellipse.area()));
    }

    #[test]
    fn test_sector_full_ellipse() {
        let ellipse = Ellipse::new(Point::new(0.0, 0.0), 3.0, 2.0, 0.0);
        let full_angle = 2.0 * PI;
        let sector = ellipse.sector_area(full_angle);
        assert!(approx_eq(sector, ellipse.area()));
    }

    #[test]
    fn test_sector_half_ellipse() {
        let ellipse = Ellipse::new(Point::new(0.0, 0.0), 4.0, 3.0, 0.0);
        let half_angle = PI;
        let sector = ellipse.sector_area(half_angle);
        assert!(approx_eq(sector, ellipse.area() / 2.0));
    }

    #[test]
    fn test_sector_quarter_ellipse() {
        let ellipse = Ellipse::new(Point::new(0.0, 0.0), 5.0, 3.0, 0.0);
        let quarter_angle = PI / 2.0;
        let sector = ellipse.sector_area(quarter_angle);
        assert!(approx_eq(sector, ellipse.area() / 4.0));
    }

    #[test]
    fn test_sector_zero_angle() {
        let ellipse = Ellipse::new(Point::new(1.0, 1.0), 3.0, 2.0, 0.5);
        let sector = ellipse.sector_area(0.0);
        assert!(approx_eq(sector, 0.0));
    }

    #[test]
    fn test_sector_circular_ellipse_matches_circle() {
        // Create a circle-like ellipse (semi_major == semi_minor)
        let radius = 3.0;
        let ellipse = Ellipse::new(Point::new(0.0, 0.0), radius, radius, 0.0);
        let circle = Circle::new(Point::new(0.0, 0.0), radius);

        let angle = PI / 3.0;
        let ellipse_sector = ellipse.sector_area(angle);
        let circle_sector = circle.sector_area(angle);

        assert!(approx_eq(ellipse_sector, circle_sector));
    }

    #[test]
    fn test_sector_circular_ellipse_multiple_angles() {
        let radius = 2.0;
        let ellipse = Ellipse::new(Point::new(0.0, 0.0), radius, radius, 0.0);
        let circle = Circle::new(Point::new(0.0, 0.0), radius);

        let test_angles = vec![
            0.0,
            PI / 6.0,
            PI / 4.0,
            PI / 2.0,
            PI,
            3.0 * PI / 2.0,
            2.0 * PI,
        ];

        for angle in test_angles {
            let ellipse_sector = ellipse.sector_area(angle);
            let circle_sector = circle.sector_area(angle);
            assert!(
                approx_eq(ellipse_sector, circle_sector),
                "Mismatch at angle {}: ellipse={}, circle={}",
                angle,
                ellipse_sector,
                circle_sector
            );
        }
    }

    #[test]
    fn test_ellipse_segment_semicircle() {
        let radius = 2.0;
        let ellipse = Ellipse::new(Point::new(0.0, 0.0), radius, radius, 0.0);

        // Points at opposite ends of diameter (along x-axis)
        let p0 = Point::new(-radius, 0.0);
        let p1 = Point::new(radius, 0.0);

        let segment = ellipse.ellipse_segment(p0, p1);
        // Semicircle segment equals half the circle area
        assert!(approx_eq(segment, ellipse.area() / 2.0));
    }

    #[test]
    fn test_ellipse_segment_circular_matches_circle_segment() {
        let radius = 3.0;
        let ellipse = Ellipse::new(Point::new(0.0, 0.0), radius, radius, 0.0);
        let circle = Circle::new(Point::new(0.0, 0.0), radius);

        // Test angle
        let angle = PI / 4.0;

        // Points on circle boundary
        let p0 = Point::new(radius, 0.0);
        let p1 = Point::new(radius * angle.cos(), radius * angle.sin());

        let ellipse_segment = ellipse.ellipse_segment(p0, p1);
        let circle_segment = circle.segment_area_from_angle(angle);

        assert!(
            approx_eq(ellipse_segment, circle_segment),
            "Circular ellipse segment should match circle segment: ellipse={}, circle={}",
            ellipse_segment,
            circle_segment
        );
    }

    #[test]
    fn test_ellipse_segment_small_angle() {
        let ellipse = Ellipse::new(Point::new(0.0, 0.0), 3.0, 2.0, 0.0);

        let angle = PI / 12.0; // 15 degrees
        let p0 = Point::new(ellipse.semi_major(), 0.0);
        let p1 = Point::new(
            ellipse.semi_major() * angle.cos(),
            ellipse.semi_minor() * angle.sin(),
        );

        let segment = ellipse.ellipse_segment(p0, p1);

        // Segment should be positive and less than the sector
        assert!(segment > 0.0);
        let sector = ellipse.sector_area(angle);
        assert!(segment < sector);
    }

    #[test]
    fn test_ellipse_segment_zero_angle() {
        let ellipse = Ellipse::new(Point::new(0.0, 0.0), 4.0, 2.5, 0.0);

        // Same point twice should give zero area
        let p = Point::new(ellipse.semi_major(), 0.0);
        let segment = ellipse.ellipse_segment(p, p);

        assert!(approx_eq(segment, 0.0));
    }

    #[test]
    fn test_ellipse_segment_with_known_values() {
        // Reference values computed for ellipse with semi_major=5.0, semi_minor=3.0
        // Total area: 47.12388980384689

        let ellipse = Ellipse::new(Point::new(0.0, 0.0), 5.0, 3.0, 0.0);

        // Quarter ellipse segment (from θ=0 to θ=π/2)
        let p0 = Point::new(ellipse.semi_major(), 0.0);
        let p1 = Point::new(0.0, ellipse.semi_minor());
        let segment_quarter = ellipse.ellipse_segment(p0, p1);
        let expected_quarter = 4.280972450961723;

        assert!(
            approx_eq(segment_quarter, expected_quarter),
            "Quarter segment: expected {}, got {}",
            expected_quarter,
            segment_quarter
        );

        // Half ellipse segment (from θ=0 to θ=π)
        let p2 = Point::new(-ellipse.semi_major(), 0.0);
        let segment_half = ellipse.ellipse_segment(p0, p2);
        let expected_half = 23.561944901923447;

        assert!(
            approx_eq(segment_half, expected_half),
            "Half segment: expected {}, got {}",
            expected_half,
            segment_half
        );

        // Half segment should also equal half the ellipse area
        assert!(approx_eq(segment_half, ellipse.area() / 2.0));
    }

    #[test]
    fn test_ellipse_segment_additional_angles() {
        let ellipse = Ellipse::new(Point::new(0.0, 0.0), 5.0, 3.0, 0.0);

        // 30 degree segment (π/6)
        let angle = PI / 6.0;
        let p0 = Point::new(ellipse.semi_major(), 0.0);
        let p1 = Point::new(
            ellipse.semi_major() * angle.cos(),
            ellipse.semi_minor() * angle.sin(),
        );
        let segment_30deg = ellipse.ellipse_segment(p0, p1);
        let expected_30deg = 0.17699081698724095;

        assert!(
            approx_eq(segment_30deg, expected_30deg),
            "30° segment: expected {}, got {}",
            expected_30deg,
            segment_30deg
        );

        // 60 degree segment (π/3)
        let angle = PI / 3.0;
        let p2 = Point::new(
            ellipse.semi_major() * angle.cos(),
            ellipse.semi_minor() * angle.sin(),
        );
        let segment_60deg = ellipse.ellipse_segment(p0, p2);
        let expected_60deg = 1.3587911055911919;

        assert!(
            approx_eq(segment_60deg, expected_60deg),
            "60° segment: expected {}, got {}",
            expected_60deg,
            segment_60deg
        );

        // Additional sanity checks
        assert!(segment_30deg < segment_60deg);
        assert!(segment_60deg < ellipse.area() / 4.0);
    }

    #[test]
    fn test_ellipse_segment_270_degrees() {
        let ellipse = Ellipse::new(Point::new(0.0, 0.0), 5.0, 3.0, 0.0);

        // 270 degree segment (3π/2)
        // When theta1 - theta0 > π, the algorithm computes the MAJOR arc segment.
        // This is the correct behavior: going CCW from 0° to 270° covers 270° of arc,
        // which is the larger of the two possible segments between these points.
        let angle = 3.0 * PI / 2.0;
        let p0 = Point::new(ellipse.semi_major(), 0.0);
        let p1 = Point::new(
            ellipse.semi_major() * angle.cos(),
            ellipse.semi_minor() * angle.sin(),
        );
        let segment_270deg = ellipse.ellipse_segment(p0, p1);
        let expected_270deg = 42.84291735288517;

        assert!(
            approx_eq(segment_270deg, expected_270deg),
            "270° segment: expected {}, got {}",
            expected_270deg,
            segment_270deg
        );

        // Verify this is indeed the major arc (> half the ellipse area)
        assert!(segment_270deg > ellipse.area() / 2.0);

        // The complementary segment (minor arc) would be the small 90° piece
        let complementary = ellipse.area() - segment_270deg;
        assert!(complementary > 0.0);
        assert!(complementary < ellipse.area() / 4.0);
    }

    #[test]
    fn test_contains_point_center_and_boundary() {
        let ellipse = Ellipse::new(Point::new(2.0, 3.0), 5.0, 3.0, 0.0);

        // Center should be inside
        assert!(ellipse.contains_point(&Point::new(2.0, 3.0)));

        // Points on the semi-major axis (horizontal)
        assert!(ellipse.contains_point(&Point::new(7.0, 3.0))); // right edge
        assert!(ellipse.contains_point(&Point::new(-3.0, 3.0))); // left edge

        // Points on the semi-minor axis (vertical)
        assert!(ellipse.contains_point(&Point::new(2.0, 6.0))); // top edge
        assert!(ellipse.contains_point(&Point::new(2.0, 0.0))); // bottom edge

        // Point clearly outside
        assert!(!ellipse.contains_point(&Point::new(10.0, 10.0)));
    }

    #[test]
    fn test_contains_point_rotated_ellipse() {
        // Create an ellipse rotated 45 degrees
        let ellipse = Ellipse::new(Point::new(0.0, 0.0), 4.0, 2.0, PI / 4.0);

        // Center should be inside
        assert!(ellipse.contains_point(&Point::new(0.0, 0.0)));

        // Point on the rotated semi-major axis
        let dist = 4.0 / (2.0_f64).sqrt(); // 4 * cos(45°) = 4 * sin(45°)
        assert!(ellipse.contains_point(&Point::new(dist, dist)));
        assert!(ellipse.contains_point(&Point::new(-dist, -dist)));

        // Point clearly outside
        assert!(!ellipse.contains_point(&Point::new(5.0, 0.0)));
        assert!(!ellipse.contains_point(&Point::new(0.0, 5.0)));

        // Point inside but not on axis
        assert!(ellipse.contains_point(&Point::new(0.5, 0.5)));
    }

    #[test]
    fn test_contains_ellipse_concentric() {
        // Larger ellipse should contain smaller concentric one
        let outer = Ellipse::new(Point::new(0.0, 0.0), 5.0, 3.0, 0.0);
        let inner = Ellipse::new(Point::new(0.0, 0.0), 3.0, 2.0, 0.0);

        assert!(outer.contains(&inner));
        assert!(!inner.contains(&outer));
    }

    #[test]
    fn test_contains_ellipse_offset_centers() {
        // Large ellipse at origin
        let outer = Ellipse::new(Point::new(0.0, 0.0), 10.0, 8.0, 0.0);

        // Small ellipse offset but still inside
        let inner = Ellipse::new(Point::new(2.0, 1.0), 2.0, 1.5, 0.0);

        assert!(outer.contains(&inner));

        // Small ellipse offset too far (outside)
        let outside = Ellipse::new(Point::new(9.0, 0.0), 2.0, 1.5, 0.0);

        assert!(!outer.contains(&outside));
    }

    #[test]
    fn test_contains_ellipse_rotated() {
        // Test with rotated ellipses
        let outer = Ellipse::new(Point::new(0.0, 0.0), 6.0, 4.0, PI / 6.0);
        let inner = Ellipse::new(Point::new(0.0, 0.0), 2.0, 1.5, PI / 4.0);

        assert!(outer.contains(&inner));
        assert!(!inner.contains(&outer));
    }

    #[test]
    fn test_contains_ellipse_intersecting() {
        // Overlapping but not containing
        let e1 = Ellipse::new(Point::new(0.0, 0.0), 5.0, 3.0, 0.0);
        let e2 = Ellipse::new(Point::new(3.0, 0.0), 4.0, 2.5, 0.0);

        assert!(!e1.contains(&e2));
        assert!(!e2.contains(&e1));
    }

    #[test]
    fn test_intersects_overlapping() {
        // Two ellipses that intersect
        let e1 = Ellipse::new(Point::new(0.0, 0.0), 5.0, 3.0, 0.0);
        let e2 = Ellipse::new(Point::new(3.0, 0.0), 4.0, 2.5, 0.0);

        assert!(e1.intersects(&e2));
        assert!(e2.intersects(&e1));
    }

    #[test]
    fn test_intersects_separate() {
        // Two ellipses that don't touch
        let e1 = Ellipse::new(Point::new(0.0, 0.0), 2.0, 1.5, 0.0);
        let e2 = Ellipse::new(Point::new(10.0, 0.0), 3.0, 2.0, 0.0);

        assert!(!e1.intersects(&e2));
        assert!(!e2.intersects(&e1));
    }

    #[test]
    fn test_intersects_contained() {
        // One ellipse contains the other - no intersection
        let outer = Ellipse::new(Point::new(0.0, 0.0), 5.0, 3.0, 0.0);
        let inner = Ellipse::new(Point::new(0.0, 0.0), 3.0, 2.0, 0.0);

        assert!(!outer.intersects(&inner));
        assert!(!inner.intersects(&outer));
    }

    #[test]
    fn test_intersects_touching() {
        // Two ellipses that barely touch (externally tangent)
        let e1 = Ellipse::new(Point::new(0.0, 0.0), 3.0, 2.0, 0.0);
        let e2 = Ellipse::new(Point::new(6.0, 0.0), 3.0, 2.0, 0.0);

        // They should intersect (tangent point counts as intersection)
        assert!(e1.intersects(&e2));
    }

    #[test]
    fn test_intersects_when_semi_minor_greater_than_semi_major() {
        // Regression for a fitter bug where the optimizer produces an ellipse
        // with `semi_minor > semi_major` (the field names are labels, not an
        // ordering invariant — `from_optimizer_params` just exponentiates the two log
        // params). The old quick-reject in `intersects` used `semi_major` as
        // the max reach, so two such ellipses whose actual long axes overlap
        // could be falsely reported as non-intersecting, breaking
        // `find_clusters` and tripping `normalize_layout`'s region-equality
        // assert. The reproducer below comes from `eulerape_3_set` seed=42.
        let e0 = Ellipse::new(
            Point::new(36.86134450095962, 78.5649458016757),
            (2.9862175109548414f64).exp(),
            (4.125140737010668f64).exp(),
            0.01270551956797799,
        );
        let e1 = Ellipse::new(
            Point::new(86.13182443212173, 20.46219955112993),
            (4.000229123217676f64).exp(),
            (3.0943704188347185f64).exp(),
            -0.023783971647556423,
        );

        assert!(
            e0.intersection_area(&e1) > 0.0,
            "sanity: ellipses do overlap by area",
        );
        assert!(
            e0.intersects(&e1),
            "intersects must agree with intersection_area on overlap",
        );
    }

    #[test]
    fn test_intersection_points_overlapping() {
        // Two ellipses that intersect at multiple points
        let e1 = Ellipse::new(Point::new(0.0, 0.0), 5.0, 3.0, 0.0);
        let e2 = Ellipse::new(Point::new(3.0, 0.0), 4.0, 2.5, 0.0);

        let points = e1.intersection_points(&e2);

        // Should have 2 or 4 intersection points (depending on overlap)
        assert!(!points.is_empty());
        assert!(points.len() <= 4);

        // Note: Due to numerical precision, points on the boundary may not satisfy
        // contains_point with the <= 1.0 check, so we just verify they exist
    }

    #[test]
    fn test_intersection_points_separate() {
        // Two ellipses that don't intersect
        let e1 = Ellipse::new(Point::new(0.0, 0.0), 2.0, 1.5, 0.0);
        let e2 = Ellipse::new(Point::new(10.0, 0.0), 3.0, 2.0, 0.0);

        let points = e1.intersection_points(&e2);

        // Note: The conic intersection algorithm may return spurious points
        // for widely separated ellipses due to numerical issues
        // In practice, these should be filtered or we should check intersects() first
        assert!(points.len() <= 4);
    }

    #[test]
    fn test_intersection_points_contained() {
        // One ellipse contains the other - no boundary intersection
        let outer = Ellipse::new(Point::new(0.0, 0.0), 5.0, 3.0, 0.0);
        let inner = Ellipse::new(Point::new(0.0, 0.0), 3.0, 2.0, 0.0);

        let points = outer.intersection_points(&inner);

        // Should have no intersection points (or handle numerical artifacts)
        // The contains check would return true, so we expect 0 or duplicate points
        assert!(points.len() <= 4);
    }

    #[test]
    fn test_intersection_points_tangent() {
        // Two circles (special case of ellipses) that are externally tangent
        let e1 = Ellipse::new(Point::new(0.0, 0.0), 2.0, 2.0, 0.0); // circle
        let e2 = Ellipse::new(Point::new(4.0, 0.0), 2.0, 2.0, 0.0); // circle

        let points = e1.intersection_points(&e2);

        // Should have 1 or 2 intersection points (tangent at one point, possibly duplicated)
        assert!(!points.is_empty());

        // The tangent point should be at (2.0, 0.0)
        if !points.is_empty() {
            let first_point = &points[0];
            assert!((first_point.x() - 2.0).abs() < 1e-6);
            assert!(first_point.y().abs() < 1e-6);
        }
    }

    #[test]
    fn test_intersection_area_no_overlap() {
        // Two ellipses that don't overlap
        let e1 = Ellipse::new(Point::new(0.0, 0.0), 2.0, 1.5, 0.0);
        let e2 = Ellipse::new(Point::new(10.0, 0.0), 2.0, 1.5, 0.0);

        let area = e1.intersection_area(&e2);
        assert!(approx_eq(area, 0.0), "Expected 0, got {}", area);
    }

    #[test]
    fn test_intersection_area_one_contains_other() {
        // Larger ellipse contains smaller one
        let e1 = Ellipse::new(Point::new(0.0, 0.0), 5.0, 4.0, 0.0);
        let e2 = Ellipse::new(Point::new(0.0, 0.0), 2.0, 1.5, 0.0);

        let area = e1.intersection_area(&e2);
        let expected = e2.area();

        assert!(
            (area - expected).abs() < 1e-8,
            "Expected {}, got {}",
            expected,
            area
        );

        // Symmetric test
        let area_reverse = e2.intersection_area(&e1);
        assert!(
            (area_reverse - expected).abs() < 1e-8,
            "Expected {}, got {}",
            expected,
            area_reverse
        );
    }

    #[test]
    fn test_intersection_area_identical_ellipses() {
        let e1 = Ellipse::new(Point::new(1.0, 2.0), 3.0, 2.0, PI / 4.0);
        let e2 = Ellipse::new(Point::new(1.0, 2.0), 3.0, 2.0, PI / 4.0);

        let area = e1.intersection_area(&e2);
        let expected = e1.area();

        assert!(
            (area - expected).abs() < 1e-8,
            "Expected {}, got {}",
            expected,
            area
        );
    }

    #[test]
    fn test_intersection_area_circles_partial_overlap() {
        // Two identical circles with partial overlap
        let radius = 3.0;
        let e1 = Ellipse::new(Point::new(0.0, 0.0), radius, radius, 0.0);
        let e2 = Ellipse::new(Point::new(4.0, 0.0), radius, radius, 0.0);

        let area = e1.intersection_area(&e2);

        // For circles, we can use the analytical formula to verify
        let c1 = Circle::new(Point::new(0.0, 0.0), radius);
        let c2 = Circle::new(Point::new(4.0, 0.0), radius);
        let expected = c1.intersection_area(&c2);

        assert!(
            (area - expected).abs() < 1e-6,
            "Expected {}, got {}",
            expected,
            area
        );
    }

    #[test]
    fn test_intersection_area_tangent_ellipses() {
        // Two circles that are externally tangent (touch at one point)
        let e1 = Ellipse::new(Point::new(0.0, 0.0), 2.0, 2.0, 0.0);
        let e2 = Ellipse::new(Point::new(4.0, 0.0), 2.0, 2.0, 0.0);

        let area = e1.intersection_area(&e2);

        // Tangent circles should have zero intersection area
        assert!(area < 1e-6, "Expected ~0, got {}", area);
    }

    #[test]
    fn test_intersection_area_symmetric() {
        // Test that intersection_area is symmetric
        let e1 = Ellipse::new(Point::new(0.0, 0.0), 4.0, 3.0, 0.0);
        let e2 = Ellipse::new(Point::new(3.0, 2.0), 3.5, 2.5, PI / 6.0);

        let area1 = e1.intersection_area(&e2);
        let area2 = e2.intersection_area(&e1);

        assert!(
            (area1 - area2).abs() < 1e-8,
            "Areas should be symmetric: {} vs {}",
            area1,
            area2
        );
    }

    #[test]
    fn test_intersection_area_bounds() {
        // Intersection area should be bounded by the smaller ellipse area
        let e1 = Ellipse::new(Point::new(0.0, 0.0), 5.0, 4.0, 0.0);
        let e2 = Ellipse::new(Point::new(2.0, 1.0), 3.0, 2.5, PI / 8.0);

        let area = e1.intersection_area(&e2);
        let min_area = e1.area().min(e2.area());

        assert!(
            area >= 0.0,
            "Intersection area should be non-negative: {}",
            area
        );
        assert!(
            area <= min_area + 1e-6,
            "Intersection area {} should not exceed smaller ellipse area {}",
            area,
            min_area
        );
    }

    #[test]
    fn test_intersection_area_rotated_ellipses() {
        // Two rotated ellipses with partial overlap
        let e1 = Ellipse::new(Point::new(0.0, 0.0), 4.0, 2.5, 0.0);
        let e2 = Ellipse::new(Point::new(3.0, 0.0), 4.0, 2.5, PI / 2.0);

        let area = e1.intersection_area(&e2);

        // Should have some intersection
        assert!(area > 0.0, "Expected positive area, got {}", area);
        assert!(
            area < e1.area() && area < e2.area(),
            "Area {} should be less than both ellipse areas",
            area
        );
    }

    #[test]
    fn test_intersection_area_nearly_identical() {
        // Two nearly identical ellipses (slightly offset)
        let e1 = Ellipse::new(Point::new(0.0, 0.0), 3.0, 2.0, 0.0);
        let e2 = Ellipse::new(Point::new(0.01, 0.01), 3.0, 2.0, 0.0);

        let area = e1.intersection_area(&e2);
        let expected = e1.area();

        // Should be very close to full area
        assert!(
            (area - expected).abs() < 0.1,
            "Expected ~{}, got {}",
            expected,
            area
        );
    }

    // Helper function for Monte Carlo validation
    fn monte_carlo_intersection_area(
        e1: &Ellipse,
        e2: &Ellipse,
        n_samples: usize,
        seed: u64,
    ) -> f64 {
        use rand::rngs::StdRng;
        use rand::{Rng, SeedableRng};

        let mut rng = StdRng::seed_from_u64(seed);
        let mut in_both = 0;

        // Bounding box for sampling - intersection of both ellipse bounding boxes
        let bbox1 = e1.bounding_box();
        let bbox2 = e2.bounding_box();
        let (min1, max1) = bbox1.to_points();
        let (min2, max2) = bbox2.to_points();

        let x_min = min1.x().max(min2.x());
        let x_max = max1.x().min(max2.x());
        let y_min = min1.y().max(min2.y());
        let y_max = max1.y().min(max2.y());

        if x_min >= x_max || y_min >= y_max {
            return 0.0; // No bounding box overlap
        }

        let bbox_area = (x_max - x_min) * (y_max - y_min);

        for _ in 0..n_samples {
            let x = x_min + (x_max - x_min) * rng.random::<f64>();
            let y = y_min + (y_max - y_min) * rng.random::<f64>();
            let p = Point::new(x, y);

            if e1.contains_point(&p) && e2.contains_point(&p) {
                in_both += 1;
            }
        }

        (in_both as f64 / n_samples as f64) * bbox_area
    }

    #[test]
    fn test_intersection_area_monte_carlo_circles() {
        // Validate circle intersection against Monte Carlo sampling
        let radius = 3.0;
        let e1 = Ellipse::new(Point::new(0.0, 0.0), radius, radius, 0.0);
        let e2 = Ellipse::new(Point::new(4.0, 0.0), radius, radius, 0.0);

        let exact = e1.intersection_area(&e2);
        let mc_estimate = monte_carlo_intersection_area(&e1, &e2, 100_000, 42);

        let error = (exact - mc_estimate).abs() / exact;
        assert!(
            error < 0.05,
            "Monte Carlo and exact should agree within 5%: exact={}, mc={}, error={:.1}%",
            exact,
            mc_estimate,
            error * 100.0
        );
    }

    #[test]
    fn test_intersection_area_monte_carlo_rotated_ellipses() {
        // Validate rotated ellipse intersection against Monte Carlo
        let e1 = Ellipse::new(Point::new(0.0, 0.0), 4.0, 2.5, 0.0);
        let e2 = Ellipse::new(Point::new(3.0, 0.0), 4.0, 2.5, PI / 2.0);

        let exact = e1.intersection_area(&e2);
        let mc_estimate = monte_carlo_intersection_area(&e1, &e2, 100_000, 42);

        let error = (exact - mc_estimate).abs() / exact.max(mc_estimate);
        assert!(
            error < 0.05,
            "Monte Carlo and exact should agree within 5%: exact={}, mc={}, error={:.1}%",
            exact,
            mc_estimate,
            error * 100.0
        );
    }

    #[test]
    fn test_intersection_area_monte_carlo_general_ellipses() {
        // Validate general ellipse intersection
        let e1 = Ellipse::new(Point::new(1.0, 0.5), 3.5, 2.0, PI / 6.0);
        let e2 = Ellipse::new(Point::new(3.0, 1.5), 3.0, 2.5, PI / 4.0);

        let exact = e1.intersection_area(&e2);
        let mc_estimate = monte_carlo_intersection_area(&e1, &e2, 150_000, 42);

        // For smaller intersections, allow slightly larger relative error
        if exact > 0.5 && mc_estimate > 0.5 {
            let error = (exact - mc_estimate).abs() / exact.max(mc_estimate);
            assert!(
                error < 0.06,
                "Monte Carlo and exact should agree within 6%: exact={}, mc={}, error={:.1}%",
                exact,
                mc_estimate,
                error * 100.0
            );
        } else {
            // Both should be small
            assert!(
                exact < 1.0 && mc_estimate < 1.0,
                "Both should be small: exact={}, mc={}",
                exact,
                mc_estimate
            );
        }
    }

    #[test]
    fn test_intersection_area_monte_carlo_contained() {
        // Validate containment case
        let e1 = Ellipse::new(Point::new(0.0, 0.0), 5.0, 4.0, 0.0);
        let e2 = Ellipse::new(Point::new(0.5, 0.5), 2.0, 1.5, PI / 8.0);

        let exact = e1.intersection_area(&e2);
        let expected = e2.area(); // Should be the smaller ellipse

        // Verify containment
        assert!(e1.contains(&e2), "e1 should contain e2");
        assert!(
            (exact - expected).abs() < 1e-6,
            "Intersection should equal smaller ellipse area: exact={}, expected={}",
            exact,
            expected
        );

        // Monte Carlo validation
        let mc_estimate = monte_carlo_intersection_area(&e1, &e2, 100_000, 42);

        let error = (exact - mc_estimate).abs() / exact;
        assert!(
            error < 0.03,
            "Monte Carlo should agree with exact (contained case): exact={}, mc={}, error={:.1}%",
            exact,
            mc_estimate,
            error * 100.0
        );
    }

    // ------------------------
    // Radius+ratio parameterization tests
    // ------------------------

    #[test]
    fn test_from_radius_ratio_circle() {
        // aspect_ratio = 1.0 should give a circle
        let e = Ellipse::from_radius_ratio(Point::new(1.0, 2.0), 3.0, 1.0, 0.5);
        assert!((e.semi_major() - 3.0).abs() < 1e-10);
        assert!((e.semi_minor() - 3.0).abs() < 1e-10);
        assert_eq!(e.center(), Point::new(1.0, 2.0));
        assert_eq!(e.rotation(), 0.5);
    }

    #[test]
    fn test_from_radius_ratio_elongated() {
        // log_aspect = ln(0.5) ≈ -0.693 means aspect_ratio = 0.5
        let radius = 2.0;
        let log_aspect = -std::f64::consts::LN_2; // ln(0.5)
        let e = Ellipse::from_radius_ratio(Point::new(0.0, 0.0), radius, log_aspect, 0.0);

        // radius = sqrt(a * b), aspect = exp(log_aspect) = 0.5
        // => a = radius / sqrt(aspect), b = radius * sqrt(aspect)
        let aspect = log_aspect.exp();
        let expected_a = radius / aspect.sqrt();
        let expected_b = radius * aspect.sqrt();

        assert!((e.semi_major() - expected_a).abs() < 1e-10);
        assert!((e.semi_minor() - expected_b).abs() < 1e-10);

        // Verify geometric mean
        let geom_mean = (e.semi_major() * e.semi_minor()).sqrt();
        assert!((geom_mean - radius).abs() < 1e-10);
    }

    #[test]
    fn test_from_radius_ratio_preserves_area() {
        // For fixed radius, changing log_aspect changes shape but preserves area
        let radius = 3.0;
        let expected_area = PI * radius * radius;

        for log_aspect in [-1.6, -0.693, -0.223, 0.0] {
            let e = Ellipse::from_radius_ratio(Point::new(0.0, 0.0), radius, log_aspect, 0.0);
            assert!((e.area() - expected_area).abs() < 1e-9);
        }
    }

    #[test]
    fn test_from_radius_ratio_clamping() {
        // log_aspect should be clamped to [-ln(5), 0] which is [-1.609, 0]
        let e_low = Ellipse::from_radius_ratio(Point::new(0.0, 0.0), 2.0, -5.0, 0.0);
        let e_high = Ellipse::from_radius_ratio(Point::new(0.0, 0.0), 2.0, 2.0, 0.0);

        // Both should be valid ellipses (clamped to valid range)
        assert!(e_low.semi_major() > 0.0);
        assert!(e_low.semi_minor() > 0.0);
        assert!(e_high.semi_major() > 0.0);
        assert!(e_high.semi_minor() > 0.0);

        // e_low should be at minimum aspect (most elongated)
        // log_aspect = -1.609 => aspect = 0.2
        assert!((e_low.semi_minor() / e_low.semi_major() - 0.2).abs() < 0.01);

        // e_high should be clamped to 0 (circle)
        assert!((e_high.semi_minor() / e_high.semi_major() - 1.0).abs() < 1e-10);
    }

    #[test]
    fn test_diagram_shape_to_params_is_geometric() {
        // Regression: `to_params` is the *geometric* encoding (linear
        // semi-axes), not the optimizer encoding. FFI consumers (e.g. the
        // eulerr R shim) rely on this — see eulerr#179.
        use crate::geometry::traits::DiagramShape;

        let e = Ellipse::new(Point::new(0.5, -1.5), 2.0, 0.5, 0.25);
        let params = e.to_params();
        assert_eq!(params.len(), 5);
        assert_eq!(params[0], 0.5); // x
        assert_eq!(params[1], -1.5); // y
        assert_eq!(params[2], 2.0); // semi_major (linear, NOT ln)
        assert_eq!(params[3], 0.5); // semi_minor (linear, NOT ln)
        assert_eq!(params[4], 0.25); // rotation

        // Round-trip through the geometric pair preserves the shape.
        let back = Ellipse::from_params(&params);
        assert_eq!(back.center(), e.center());
        assert!((back.semi_major() - e.semi_major()).abs() < 1e-12);
        assert!((back.semi_minor() - e.semi_minor()).abs() < 1e-12);
        assert!((back.rotation() - e.rotation()).abs() < 1e-12);
    }

    #[test]
    fn test_diagram_shape_optimizer_params_from_circle() {
        use crate::geometry::traits::DiagramShape;

        // optimizer_params_from_circle stores semi-axes in log space.
        let params = Ellipse::optimizer_params_from_circle(1.0, 2.0, 3.0);
        let log3 = 3.0_f64.ln();
        assert_eq!(params.len(), 5);
        assert_eq!(params[0], 1.0); // x
        assert_eq!(params[1], 2.0); // y
        assert!((params[2] - log3).abs() < 1e-12); // ln(semi_major)
        assert!((params[3] - log3).abs() < 1e-12); // ln(semi_minor)
        assert_eq!(params[4], 0.0); // rotation

        // from_optimizer_params should reconstruct a circle (exp roundtrip).
        let e = Ellipse::from_optimizer_params(&params);
        assert!((e.semi_major() - 3.0).abs() < 1e-10);
        assert!((e.semi_minor() - 3.0).abs() < 1e-10);
    }

    // ------------------------
    // Exclusive region tests
    // ------------------------

    #[test]
    fn test_exclusive_regions_two_ellipses_partial_overlap() {
        use crate::geometry::traits::DiagramShape;

        let e1 = Ellipse::new(Point::new(0.0, 0.0), 4.0, 2.5, 0.0);
        let e2 = Ellipse::new(Point::new(3.0, 1.0), 3.5, 2.0, PI / 6.0);

        let areas = Ellipse::compute_exclusive_regions(&[e1, e2]);
        let a_only = areas.get(&(1usize << 0)).copied().unwrap_or(0.0);
        let b_only = areas.get(&(1usize << 1)).copied().unwrap_or(0.0);
        let both = areas
            .get(&((1usize << 0) | (1usize << 1)))
            .copied()
            .unwrap_or(0.0);

        // Basic sanity: non-negative and sums to union
        assert!(a_only >= 0.0 && b_only >= 0.0 && both >= 0.0);
        let union = a_only + b_only + both;
        // Union must be <= sum of individual ellipse areas
        assert!(union <= e1.area() + e2.area() + 1e-6);

        // Intersection area should be close to direct computation
        let exact = e1.intersection_area(&e2);
        let error = (both - exact).abs() / exact.max(1e-6);
        assert!(
            error < 0.1,
            "Error too large: both={}, exact={}, error={:.2}%",
            both,
            exact,
            error * 100.0
        );
    }

    #[test]
    fn test_exclusive_regions_two_ellipses_contained() {
        use crate::geometry::traits::DiagramShape;

        let outer = Ellipse::new(Point::new(0.0, 0.0), 6.0, 4.0, 0.2);
        let inner = Ellipse::new(Point::new(0.5, 0.3), 2.0, 1.5, -0.3);

        let areas = Ellipse::compute_exclusive_regions(&[outer, inner]);
        let inner_only = areas.get(&(1usize << 1)).copied().unwrap_or(0.0);
        let both = areas
            .get(&((1usize << 0) | (1usize << 1)))
            .copied()
            .unwrap_or(0.0);

        // When contained, exclusive intersection should equal inner area
        assert!((both - inner.area()).abs() < 1e-3);
        // Inner exclusive should be ~0
        assert!(inner_only < 1e-6);
    }

    #[test]
    fn test_exclusive_regions_three_ellipses_basic() {
        use crate::geometry::traits::DiagramShape;

        let e1 = Ellipse::new(Point::new(-2.0, 0.0), 3.0, 2.0, 0.2);
        let e2 = Ellipse::new(Point::new(2.0, 0.5), 3.2, 2.1, -0.1);
        let e3 = Ellipse::new(Point::new(0.0, 1.5), 2.8, 1.8, 0.5);

        let areas = Ellipse::compute_exclusive_regions(&[e1, e2, e3]);
        let mask_all = (1usize << 0) | (1usize << 1) | (1usize << 2);
        let all = areas.get(&mask_all).copied().unwrap_or(0.0);

        // Non-negative and bounded by smallest ellipse area
        let min_area = e1.area().min(e2.area()).min(e3.area());
        assert!(all >= 0.0 && all <= min_area + 1e-6);

        // Union should not exceed sum of individuals
        let union: f64 = areas.values().sum();
        assert!(union <= e1.area() + e2.area() + e3.area() + 1e-3);
    }

    // Helper function for Monte Carlo validation of pairwise intersection
    fn monte_carlo_pairwise_intersection(
        e1: &Ellipse,
        e2: &Ellipse,
        n_samples: usize,
        seed: u64,
    ) -> f64 {
        use rand::rngs::StdRng;
        use rand::{Rng, SeedableRng};

        let mut rng = StdRng::seed_from_u64(seed);
        let mut in_both = 0;

        // Bounding box for sampling - intersection of both ellipse bounding boxes
        let bbox1 = e1.bounding_box();
        let bbox2 = e2.bounding_box();
        let (min1, max1) = bbox1.to_points();
        let (min2, max2) = bbox2.to_points();

        let x_min = min1.x().max(min2.x());
        let x_max = max1.x().min(max2.x());
        let y_min = min1.y().max(min2.y());
        let y_max = max1.y().min(max2.y());

        if x_min >= x_max || y_min >= y_max {
            return 0.0; // No bounding box overlap
        }

        let bbox_area = (x_max - x_min) * (y_max - y_min);

        for _ in 0..n_samples {
            let x = x_min + (x_max - x_min) * rng.random::<f64>();
            let y = y_min + (y_max - y_min) * rng.random::<f64>();
            let p = Point::new(x, y);

            if e1.contains_point(&p) && e2.contains_point(&p) {
                in_both += 1;
            }
        }

        (in_both as f64 / n_samples as f64) * bbox_area
    }

    #[test]
    #[ignore = "slow stochastic validation"]
    fn test_random_ellipse_intersections_monte_carlo() {
        use rand::rngs::StdRng;
        use rand::{Rng, SeedableRng};

        const N_TESTS: usize = 100;
        const MC_SAMPLES: usize = 50_000;
        const MAX_RELATIVE_ERROR: f64 = 0.10; // 10% tolerance

        let mut failures = Vec::new();

        for seed in 0..N_TESTS {
            let mut rng = StdRng::seed_from_u64(seed as u64);

            // Random ellipse 1
            let x1 = rng.random_range(-5.0..5.0);
            let y1 = rng.random_range(-5.0..5.0);
            let a1 = rng.random_range(0.5..3.0);
            let b1 = rng.random_range(0.5..3.0);
            let rot1 = rng.random_range(0.0..PI);

            // Random ellipse 2 (overlapping region)
            let x2 = x1 + rng.random_range(-2.0..2.0);
            let y2 = y1 + rng.random_range(-2.0..2.0);
            let a2 = rng.random_range(0.5..3.0);
            let b2 = rng.random_range(0.5..3.0);
            let rot2 = rng.random_range(0.0..PI);

            let e1 = Ellipse::new(Point::new(x1, y1), a1, b1, rot1);
            let e2 = Ellipse::new(Point::new(x2, y2), a2, b2, rot2);

            // Compute exact intersection area
            let exact = e1.intersection_area(&e2);

            // Compute Monte Carlo estimate
            let mc_estimate = monte_carlo_pairwise_intersection(&e1, &e2, MC_SAMPLES, seed as u64);

            // Skip cases with very small intersection (harder to validate)
            if exact < 0.01 && mc_estimate < 0.01 {
                continue;
            }

            // Compute relative error
            let max_val = exact.max(mc_estimate);
            let error = (exact - mc_estimate).abs() / max_val;

            if error > MAX_RELATIVE_ERROR {
                let n_points = e1.intersection_points(&e2).len();
                failures.push((seed, exact, mc_estimate, error, n_points));
            }
        }

        assert!(
            failures.is_empty(),
            "Monte Carlo validation failed for {} of {} cases: {:?}",
            failures.len(),
            N_TESTS,
            failures.iter().take(10).collect::<Vec<_>>()
        );
    }

    #[test]
    fn test_three_ellipse_complete_overlap() {
        use crate::geometry::traits::DiagramShape;

        // Test case: A=B=C=1, A&B&C=1 (all three completely overlap)
        // Three identical circles
        let e1 = Ellipse::new(Point::new(0.0, 0.0), 0.564, 0.564, 0.0); // area ≈ 1.0
        let e2 = Ellipse::new(Point::new(0.0, 0.0), 0.564, 0.564, 0.0);
        let e3 = Ellipse::new(Point::new(0.0, 0.0), 0.564, 0.564, 0.0);

        let areas = Ellipse::compute_exclusive_regions(&[e1, e2, e3]);

        let mask_all = (1usize << 0) | (1usize << 1) | (1usize << 2);
        let all_three = areas.get(&mask_all).copied().unwrap_or(0.0);

        // When all three are identical and overlapping, the 3-way intersection should be ~1.0
        assert!(
            (all_three - 1.0).abs() < 0.01,
            "Complete overlap should give area ~1.0, got {}",
            all_three
        );
    }

    #[test]
    fn test_three_ellipse_partial_overlap_monte_carlo() {
        use rand::rngs::StdRng;
        use rand::{Rng, SeedableRng};

        // Test a specific 3-ellipse configuration with Monte Carlo
        let seed = 42u64;
        let mut rng = StdRng::seed_from_u64(seed);

        let x1 = 0.0;
        let y1 = 0.0;
        let a1 = 1.5;
        let b1 = 1.0;
        let rot1 = 0.0;

        let x2 = 1.0;
        let y2 = 0.0;
        let a2 = 1.5;
        let b2 = 1.0;
        let rot2 = PI / 6.0;

        let x3 = 0.5;
        let y3 = 1.0;
        let a3 = 1.5;
        let b3 = 1.0;
        let rot3 = -PI / 6.0;

        let e1 = Ellipse::new(Point::new(x1, y1), a1, b1, rot1);
        let e2 = Ellipse::new(Point::new(x2, y2), a2, b2, rot2);
        let e3 = Ellipse::new(Point::new(x3, y3), a3, b3, rot3);

        // Monte Carlo estimate for 3-way intersection
        let n_samples = 100_000;
        let mut in_all_three = 0;

        // Find bounding box
        let bbox1 = e1.bounding_box();
        let bbox2 = e2.bounding_box();
        let bbox3 = e3.bounding_box();

        let (min1, max1) = bbox1.to_points();
        let (min2, max2) = bbox2.to_points();
        let (min3, max3) = bbox3.to_points();

        let x_min = min1.x().max(min2.x()).max(min3.x());
        let x_max = max1.x().min(max2.x()).min(max3.x());
        let y_min = min1.y().max(min2.y()).max(min3.y());
        let y_max = max1.y().min(max2.y()).min(max3.y());

        if x_min >= x_max || y_min >= y_max {
            return;
        }

        let bbox_area = (x_max - x_min) * (y_max - y_min);

        for _ in 0..n_samples {
            let x = x_min + (x_max - x_min) * rng.random::<f64>();
            let y = y_min + (y_max - y_min) * rng.random::<f64>();
            let p = Point::new(x, y);

            if e1.contains_point(&p) && e2.contains_point(&p) && e3.contains_point(&p) {
                in_all_three += 1;
            }
        }

        let mc_area = (in_all_three as f64 / n_samples as f64) * bbox_area;

        // Compute using compute_exclusive_regions
        use crate::geometry::traits::DiagramShape;
        let areas = Ellipse::compute_exclusive_regions(&[e1, e2, e3]);
        let mask_all = (1usize << 0) | (1usize << 1) | (1usize << 2);
        let exact_area = areas.get(&mask_all).copied().unwrap_or(0.0);

        if mc_area < 0.01 && exact_area < 0.01 {
            return;
        }

        let error = (exact_area - mc_area).abs() / mc_area.max(exact_area);

        assert!(
            error < 0.15,
            "Three-ellipse Monte Carlo error too large: {:.1}%",
            error * 100.0
        );
    }

    // ===== Container clipping (compute_exclusive_regions_clipped_ellipse) =====

    /// An axis-aligned ellipse fully inside the container should keep its
    /// full area; the complement region is the rest of the container.
    #[test]
    fn clipped_single_ellipse_inside_container() {
        // a = 2, b = 1 → area = 2π.
        let ellipses = vec![Ellipse::new(Point::new(0.0, 0.0), 2.0, 1.0, 0.0)];
        let container = Rectangle::new(Point::new(0.0, 0.0), 10.0, 6.0); // 60
        let areas = compute_exclusive_regions_clipped_ellipse(&ellipses, &container);

        let ellipse_area = areas.get(&0b1).copied().unwrap_or(0.0);
        let complement = areas.get(&0).copied().unwrap_or(0.0);
        let expected_ellipse = 2.0 * PI;

        assert!(
            (ellipse_area - expected_ellipse).abs() < 1e-9,
            "ellipse area {} ≠ 2π",
            ellipse_area
        );
        assert!(
            (complement - (60.0 - expected_ellipse)).abs() < 1e-9,
            "complement {} ≠ 60 − 2π",
            complement
        );
    }

    /// A rotated ellipse fully inside the container should still keep its
    /// full area πab.
    #[test]
    fn clipped_rotated_ellipse_inside_container() {
        let ellipses = vec![Ellipse::new(
            Point::new(0.0, 0.0),
            3.0,
            2.0,
            std::f64::consts::FRAC_PI_4,
        )];
        // Bounding box of rotated (a=3,b=2) is at most 3√2 ≈ 4.24 in each
        // dimension; 12×12 container gives plenty of room.
        let container = Rectangle::new(Point::new(0.0, 0.0), 12.0, 12.0);
        let areas = compute_exclusive_regions_clipped_ellipse(&ellipses, &container);

        let ellipse_area = areas.get(&0b1).copied().unwrap_or(0.0);
        let complement = areas.get(&0).copied().unwrap_or(0.0);
        let expected_ellipse = 6.0 * PI;

        assert!(
            (ellipse_area - expected_ellipse).abs() < 1e-9,
            "rotated ellipse area {} ≠ 6π",
            ellipse_area
        );
        assert!(
            (complement - (144.0 - expected_ellipse)).abs() < 1e-9,
            "complement {} ≠ 144 − 6π",
            complement
        );
    }

    /// An ellipse entirely outside the container should contribute zero; the
    /// complement equals the full container area.
    #[test]
    fn clipped_single_ellipse_outside_container() {
        let ellipses = vec![Ellipse::new(Point::new(100.0, 100.0), 2.0, 1.0, 0.3)];
        let container = Rectangle::new(Point::new(0.0, 0.0), 4.0, 4.0); // 16
        let areas = compute_exclusive_regions_clipped_ellipse(&ellipses, &container);

        let ellipse_area = areas.get(&0b1).copied().unwrap_or(0.0);
        let complement = areas.get(&0).copied().unwrap_or(0.0);

        assert!(
            ellipse_area.abs() < 1e-9,
            "ellipse area {} ≠ 0",
            ellipse_area
        );
        assert!(
            (complement - 16.0).abs() < 1e-9,
            "complement {} ≠ 16",
            complement
        );
    }

    /// Container fully inside an ellipse: clipped area equals container area;
    /// complement is zero.
    #[test]
    fn clipped_container_fully_inside_ellipse() {
        // Big ellipse covering the entire container (and then some).
        let ellipses = vec![Ellipse::new(Point::new(0.0, 0.0), 100.0, 80.0, 0.0)];
        let container = Rectangle::new(Point::new(0.0, 0.0), 4.0, 6.0);
        let areas = compute_exclusive_regions_clipped_ellipse(&ellipses, &container);

        let inside = areas.get(&0b1).copied().unwrap_or(0.0);
        let complement = areas.get(&0).copied().unwrap_or(0.0);

        assert!(
            (inside - 24.0).abs() < 1e-9,
            "ellipse-clipped {} ≠ container area",
            inside
        );
        assert!(complement.abs() < 1e-9, "complement {} ≠ 0", complement);
    }

    /// Sanity: with no shapes, the complement equals the container area.
    #[test]
    fn clipped_no_ellipses_complement_equals_container_area() {
        let ellipses: Vec<Ellipse> = Vec::new();
        let container = Rectangle::new(Point::new(1.0, 2.0), 3.0, 4.0);
        let areas = compute_exclusive_regions_clipped_ellipse(&ellipses, &container);

        let complement = areas.get(&0).copied().unwrap_or(0.0);
        assert!(
            (complement - 12.0).abs() < 1e-9,
            "complement {} ≠ 12",
            complement
        );
    }

    /// Two overlapping ellipses both fully inside the container: the per-mask
    /// exclusive areas sum to the union area; complement = box − union.
    #[test]
    fn clipped_two_ellipses_inside_sum_to_union_area() {
        // Two identical axis-aligned ellipses (a=1.5, b=1) horizontally
        // separated. Both inside an 8×6 box.
        let ellipses = vec![
            Ellipse::new(Point::new(-0.7, 0.0), 1.5, 1.0, 0.0),
            Ellipse::new(Point::new(0.7, 0.0), 1.5, 1.0, 0.0),
        ];
        let container = Rectangle::new(Point::new(0.0, 0.0), 8.0, 6.0); // 48
        let areas = compute_exclusive_regions_clipped_ellipse(&ellipses, &container);

        let a_only = areas.get(&0b01).copied().unwrap_or(0.0);
        let b_only = areas.get(&0b10).copied().unwrap_or(0.0);
        let a_and_b = areas.get(&0b11).copied().unwrap_or(0.0);
        let complement = areas.get(&0).copied().unwrap_or(0.0);

        // Union via per-mask sum.
        let union = a_only + b_only + a_and_b;
        // Reference: 2·area − lens area (intersection_area is exact for
        // ellipses).
        let lens = ellipses[0].intersection_area(&ellipses[1]);
        let expected_union = 2.0 * (1.5 * PI) - lens;

        assert!(
            (union - expected_union).abs() < 1e-9,
            "union via per-mask sum {} ≠ {}",
            union,
            expected_union
        );
        assert!(
            (complement - (48.0 - expected_union)).abs() < 1e-9,
            "complement {} ≠ box − union",
            complement
        );
    }

    /// Ellipse clipped by exactly one edge: validate against a Monte-Carlo
    /// reference. Closed-form for a rotated-ellipse single-edge cut is
    /// painful, so MC is more practical here.
    #[test]
    fn clipped_single_ellipse_one_edge_matches_monte_carlo() {
        // Rotated ellipse at origin; container's top edge cuts through.
        let ellipses = vec![Ellipse::new(Point::new(0.0, 0.0), 2.0, 1.0, 0.4)];
        let container = Rectangle::new(Point::new(0.0, -2.0), 100.0, 4.5);
        // Container bounds: x ∈ [−50, 50], y ∈ [−4.25, 0.25]. Top edge at
        // y = 0.25 cuts the rotated ellipse.
        let areas = compute_exclusive_regions_clipped_ellipse(&ellipses, &container);

        let kept = areas.get(&0b1).copied().unwrap_or(0.0);

        // Monte-Carlo reference: sample uniformly in the bounding box of the
        // rotated ellipse (which fits in [-2.2, 2.2]² for a=2, b=1).
        let n = 400_000;
        let mut hits = 0;
        let mut state: u64 = 0xDEAD_BEEF_CAFE_BABE;
        let next_u01 = |s: &mut u64| -> f64 {
            *s = s
                .wrapping_mul(6364136223846793005)
                .wrapping_add(1442695040888963407);
            ((*s >> 11) as f64) / ((1u64 << 53) as f64)
        };
        for _ in 0..n {
            let x = -2.2 + 4.4 * next_u01(&mut state);
            let y = -2.2 + 4.4 * next_u01(&mut state);
            let in_ellipse = ellipses[0].contains_point(&Point::new(x, y));
            let in_box = (-50.0..=50.0).contains(&x) && (-4.25..=0.25).contains(&y);
            if in_ellipse && in_box {
                hits += 1;
            }
        }
        let bbox_area = 4.4 * 4.4;
        let mc_area = bbox_area * (hits as f64) / (n as f64);
        // 400k samples → stderr ≈ √(p(1-p)/n) · bbox_area; about 0.02
        // empirically.
        assert!(
            (kept - mc_area).abs() < 0.05,
            "clipped ellipse {} vs MC {}",
            kept,
            mc_area
        );
    }

    // ===== Clipped exclusive areas: analytical gradient (S4) =====

    /// Pack `(ellipses, container)` into the same flat parameter layout the
    /// optimiser uses: per-ellipse `[x, y, ln a, ln b, φ]` blocks, then
    /// container `[x_c, y_c, ln(area), ln(ratio)]`.
    fn pack_clipped_params_ellipse(ellipses: &[Ellipse], container: &Rectangle) -> Vec<f64> {
        let mut p = Vec::with_capacity(ellipses.len() * 5 + 4);
        for e in ellipses {
            p.extend(e.to_optimizer_params());
        }
        p.extend(container.to_optimizer_params());
        p
    }

    /// Inverse of `pack_clipped_params_ellipse`: decode a flat parameter slice
    /// back into ellipses plus container.
    fn unpack_clipped_params_ellipse(p: &[f64], n_sets: usize) -> (Vec<Ellipse>, Rectangle) {
        let ellipses: Vec<Ellipse> = (0..n_sets)
            .map(|i| Ellipse::from_optimizer_params(&p[5 * i..5 * (i + 1)]))
            .collect();
        let container = Rectangle::from_optimizer_params(&p[5 * n_sets..5 * n_sets + 4]);
        (ellipses, container)
    }

    /// FD vs analytical gradient comparison for the clipped per-mask area
    /// helper. `params` is laid out as in `pack_clipped_params_ellipse`.
    fn assert_clipped_gradient_matches_fd_ellipse(
        ellipses: &[Ellipse],
        container: &Rectangle,
        h: f64,
        tol: f64,
        label: &str,
    ) {
        let (areas, grads) =
            compute_exclusive_regions_clipped_ellipse_with_gradient(ellipses, container);
        let n_sets = ellipses.len();
        let n_params = n_sets * 5 + 4;
        let p0 = pack_clipped_params_ellipse(ellipses, 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 (ep, kp) = unpack_clipped_params_ellipse(&plus, n_sets);
                let (em, km) = unpack_clipped_params_ellipse(&minus, n_sets);
                let ap = compute_exclusive_regions_clipped_ellipse(&ep, &kp)
                    .get(mask)
                    .copied()
                    .unwrap_or(0.0);
                let am = compute_exclusive_regions_clipped_ellipse(&em, &km)
                    .get(mask)
                    .copied()
                    .unwrap_or(0.0);
                fd[i] = (ap - am) / (2.0 * h);
            }
            // Check the analytical area matches the value-only path.
            let direct = compute_exclusive_regions_clipped_ellipse(ellipses, 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
            );
        }
    }

    /// Two ellipses fully inside a wide container. No box edge clips any
    /// ellipse — only the complement region touches the box. Verifies that
    /// shape-param gradients reproduce the unclipped gradient and that the
    /// container gradient is purely (0, 0, container_area, 0).
    #[test]
    fn clipped_grad_two_ellipses_inside_box_matches_fd() {
        let ellipses = vec![
            Ellipse::new(Point::new(-0.7, 0.05), 1.1, 0.8, 0.0),
            Ellipse::new(Point::new(0.6, -0.04), 0.95, 0.6, 0.0),
        ];
        let container = Rectangle::new(Point::new(0.0, 0.0), 6.0, 5.0);
        assert_clipped_gradient_matches_fd_ellipse(
            &ellipses,
            &container,
            1e-6,
            1e-5,
            "two_ellipses_inside",
        );
    }

    /// Two ellipses each clipped by a different box edge. Exercises the
    /// box-edge boundary contribution to the container-param gradient and
    /// the trimmed-arc contribution to shape params.
    #[test]
    fn clipped_grad_two_ellipses_each_touching_an_edge_matches_fd() {
        let ellipses = vec![
            Ellipse::new(Point::new(1.4, 0.05), 0.9, 0.7, 0.0),
            Ellipse::new(Point::new(0.0, -0.07), 0.85, 0.6, 0.0),
        ];
        let container = Rectangle::new(Point::new(0.0, 0.0), 3.6, 3.0);
        assert_clipped_gradient_matches_fd_ellipse(
            &ellipses,
            &container,
            1e-6,
            1e-4,
            "two_ellipses_clipped_edge",
        );
    }

    /// Three overlapping ellipses inside the box, all axis-aligned.
    /// Exercises 3-way intersection regions plus the complement.
    #[test]
    fn clipped_grad_three_ellipses_inside_box_matches_fd() {
        let ellipses = vec![
            Ellipse::new(Point::new(-0.5, -0.3), 1.0, 0.7, 0.0),
            Ellipse::new(Point::new(0.5, -0.3), 1.0, 0.7, 0.0),
            Ellipse::new(Point::new(0.0, 0.55), 1.0, 0.7, 0.0),
        ];
        let container = Rectangle::new(Point::new(0.0, 0.05), 3.5, 3.5);
        assert_clipped_gradient_matches_fd_ellipse(
            &ellipses,
            &container,
            1e-6,
            1e-4,
            "three_ellipses",
        );
    }

    /// Two rotated ellipses with the box clipping the right-most ellipse.
    /// Exercises non-zero rotation in both the arc-gradient terms and the
    /// box-edge intersection helpers.
    #[test]
    fn clipped_grad_rotated_ellipses_one_clipped_matches_fd() {
        let ellipses = vec![
            Ellipse::new(Point::new(-0.4, -0.1), 1.0, 0.6, 0.5),
            Ellipse::new(Point::new(0.7, 0.05), 0.95, 0.55, -0.3),
        ];
        // Right edge cuts through the right ellipse; left ellipse fully inside.
        let container = Rectangle::new(Point::new(0.0, 0.0), 2.6, 2.4);
        assert_clipped_gradient_matches_fd_ellipse(
            &ellipses,
            &container,
            1e-6,
            1e-4,
            "rotated_one_clipped",
        );
    }
}