//! Region decomposition for diagram visualization.
//!
//! This module provides utilities for decomposing fitted shapes into
//! exclusive regions (one per set combination) for plotting.
//!
//! # Output structure
//!
//! [`RegionPolygons`] stores each region as a list of [`RegionPiece`]s, where
//! a *piece* is a single connected component with one outer boundary and
//! zero or more interior holes. Disconnected regions (e.g. an "A only"
//! region split by other sets cutting across A) produce multiple pieces.
//!
//! Outer rings are normalised to **CCW** (positive signed area) and holes
//! to **CW** (negative signed area) so renderers can draw each piece with
//! `fill-rule: nonzero` (the SVG default) without further bookkeeping.

use crate::geometry::diagram::{IntersectionPoint, discover_regions};
use crate::geometry::primitives::Point;
use crate::geometry::shapes::{Polygon, Rectangle};
use crate::geometry::traits::{Closed, DiagramShape, Polygonize};
use crate::plotting::clip::{ClipOperation, polygon_clip_many, polygon_difference};
use crate::spec::{Combination, DiagramSpec};
use std::collections::HashMap;

/// A connected component of a region: one outer boundary plus zero or more
/// interior holes (other regions or shapes cutting through this piece).
///
/// Orientation is normalised by the core library so consumers can render
/// each piece with the SVG default `fill-rule: nonzero` (or any equivalent
/// nonzero / "even-odd-of-rings" rule) **without classifying rings on their
/// own side**:
///
/// * `outer` — CCW (positive signed area).
/// * `holes` — CW (negative signed area), each strictly inside `outer`.
///
/// Use [`RegionPiece::area`] for the piece's net area (outer minus holes).
///
/// # Renderer recipe
///
/// SVG: emit one `<path>` per piece with `d` formed by concatenating the
/// outer ring's `M…Z` and each hole's `M…Z`. Default `fill-rule: nonzero`
/// fills only the donut/cookie shape because the rings have opposite
/// winding.
///
/// Canvas: `ctx.beginPath()`, walk the outer ring, walk each hole ring,
/// `ctx.fill()` (Canvas's default fill is also nonzero).
///
/// # Construction
///
/// `RegionPiece`s are produced by [`decompose_regions`] (and the higher-
/// level [`crate::Layout::region_polygons`]). Bindings authors should not
/// build them by hand: the topological classification that maps a flat
/// list of clipper-output rings to outer/holes pairs is not part of the
/// public API.
///
/// `#[non_exhaustive]`: an output type produced by [`decompose_regions`], not
/// constructed by downstream callers, so fields can be added non-breakingly.
#[derive(Debug, Clone)]
#[non_exhaustive]
pub struct RegionPiece {
    /// CCW outer boundary of the piece (always non-degenerate after
    /// [`decompose_regions`]).
    pub outer: Polygon,
    /// CW hole rings, each strictly inside `outer`. Empty when the piece
    /// is simply connected.
    pub holes: Vec<Polygon>,
}

impl RegionPiece {
    /// Net area of the piece — `outer` minus the sum of `holes`. Always
    /// non-negative because orientations are normalised.
    pub fn area(&self) -> f64 {
        let outer_area = self.outer.area();
        let hole_area: f64 = self.holes.iter().map(|h| h.area()).sum();
        (outer_area - hole_area).max(0.0)
    }
}

/// Collection of region pieces for each exclusive region in a diagram.
///
/// Each key is a [`Combination`] of set names, and the value is a list of
/// [`RegionPiece`]s that together represent that exclusive region. Multiple
/// pieces occur when a region is disconnected (e.g. an "A only" lobe split
/// by an intersection cutting across A); piece orientations and hole
/// containment are guaranteed by [`decompose_regions`] — see
/// [`RegionPiece`] for the rendering contract.
///
/// Pieces below the per-diagram sliver threshold (see [`decompose_regions`])
/// are filtered out before insertion; consumers can assume every retained
/// piece has positive net area.
#[derive(Debug, Clone)]
pub struct RegionPolygons {
    regions: HashMap<Combination, Vec<RegionPiece>>,
}

impl RegionPolygons {
    /// Creates a new empty `RegionPolygons` collection.
    ///
    /// `pub(crate)` on purpose: the only sound way to populate one is via
    /// [`decompose_regions`] (or [`crate::Layout::region_polygons`]), which
    /// guarantees the [`RegionPiece`] orientation and topology contract.
    /// Hand-built collections would let callers violate those invariants
    /// silently.
    pub(crate) fn new() -> Self {
        Self {
            regions: HashMap::new(),
        }
    }

    /// Adds pieces for a given region. `pub(crate)` to keep the orientation
    /// contract on [`RegionPiece`] enforceable — production code should use
    /// [`decompose_regions`].
    pub(crate) fn insert(&mut self, combination: Combination, pieces: Vec<RegionPiece>) {
        self.regions.insert(combination, pieces);
    }

    /// Wraps a pre-built `combination → pieces` map into a [`RegionPolygons`]
    /// without further validation.
    ///
    /// The caller is responsible for the [`RegionPiece`] orientation contract
    /// (outer rings CCW, hole rings CW, holes strictly inside their outer)
    /// and for region disjointness. **Most callers should use
    /// [`decompose_regions`]** (or [`crate::Layout::region_polygons`])
    /// instead — those guarantee the contract.
    ///
    /// This entry point exists for downstream bindings (WASM, R, Python,
    /// Julia) that already have a decomposed layout in hand and want to
    /// thread it back through APIs like [`crate::plotting::place_labels`]
    /// without re-running the fitter.
    pub fn from_map(regions: HashMap<Combination, Vec<RegionPiece>>) -> Self {
        Self { regions }
    }

    /// Returns the pieces for a given region, or `None` if the combination
    /// has no fitted (non-empty) area.
    pub fn get(&self, combination: &Combination) -> Option<&Vec<RegionPiece>> {
        self.regions.get(combination)
    }

    /// Iterates over every non-empty region and its pieces in canonical
    /// `(size, lexicographic)` order: singletons first, then pairs, then
    /// triples, etc.; ties broken by the lexicographic order of the
    /// (already-sorted) set names.
    ///
    /// The order is **deterministic** and independent of how the collection
    /// was built — bindings, snapshot tests, and CSV exports can rely on it.
    /// (The backing store is a `HashMap`, so the entries are collected and
    /// sorted on each call; the cost is negligible since the region count is
    /// at most the populated subset of `2^n`.)
    ///
    /// If you want the order to follow a spec's input order instead, use
    /// [`Self::iter_in_input_order`].
    pub fn iter(&self) -> impl Iterator<Item = (&Combination, &Vec<RegionPiece>)> {
        let mut entries: Vec<(&Combination, &Vec<RegionPiece>)> = self.regions.iter().collect();
        entries
            .sort_by(|(a, _), (b, _)| a.len().cmp(&b.len()).then_with(|| a.sets().cmp(b.sets())));
        entries.into_iter()
    }

    /// Iterates in canonical `(size, lexicographic)` order. This is an alias
    /// for [`Self::iter`], which now yields the same deterministic order;
    /// retained for discoverability and backwards compatibility.
    pub fn iter_sorted(&self) -> impl Iterator<Item = (&Combination, &Vec<RegionPiece>)> {
        self.iter()
    }

    /// Iterates in spec input order: combinations are sorted by the
    /// position of their member sets in `set_names` (typically
    /// `spec.set_names()`), with smaller combinations before larger.
    /// Sets in a combination that don't appear in `set_names` are
    /// placed last by convention; if every combination's sets are present
    /// in `set_names`, this matches the order a user-facing summary
    /// would naturally take.
    ///
    /// # Examples
    ///
    /// ```
    /// use eunoia::{DiagramSpecBuilder, Fitter, InputType};
    /// use eunoia::geometry::shapes::Circle;
    ///
    /// let spec = DiagramSpecBuilder::new()
    ///     .set("A", 5.0)
    ///     .set("B", 3.0)
    ///     .intersection(&["A", "B"], 1.0)
    ///     .input_type(InputType::Exclusive)
    ///     .build()
    ///     .unwrap();
    ///
    /// let layout = Fitter::<Circle>::new(&spec).seed(0).fit().unwrap();
    /// let regions = layout.region_polygons(&spec, 64);
    ///
    /// // Pull out the keys in input order.
    /// let order: Vec<_> = regions
    ///     .iter_in_input_order(spec.set_names())
    ///     .map(|(c, _)| c.to_string())
    ///     .collect();
    /// // Singletons before the pair, and "A" before "B".
    /// assert_eq!(order, vec!["A", "B", "A&B"]);
    /// ```
    pub fn iter_in_input_order<'a, N: AsRef<str>>(
        &'a self,
        set_names: &'a [N],
    ) -> impl Iterator<Item = (&'a Combination, &'a Vec<RegionPiece>)> + 'a {
        // Build a name → index lookup once. Sets not found in `set_names`
        // sort after every known set (sentinel = `set_names.len()`).
        let lookup: HashMap<&str, usize> = set_names
            .iter()
            .enumerate()
            .map(|(i, name)| (name.as_ref(), i))
            .collect();
        let sentinel = set_names.len();
        let key = move |combo: &Combination| -> Vec<usize> {
            let mut idxs: Vec<usize> = combo
                .sets()
                .iter()
                .map(|s| lookup.get(s.as_str()).copied().unwrap_or(sentinel))
                .collect();
            idxs.sort_unstable();
            idxs
        };

        let mut entries: Vec<(&Combination, &Vec<RegionPiece>)> = self.regions.iter().collect();
        entries.sort_by(|(a, _), (b, _)| a.len().cmp(&b.len()).then_with(|| key(a).cmp(&key(b))));
        entries.into_iter()
    }

    /// Number of distinct non-empty regions.
    pub fn len(&self) -> usize {
        self.regions.len()
    }

    /// `true` when no regions are stored.
    pub fn is_empty(&self) -> bool {
        self.regions.is_empty()
    }

    /// Computes the total (net) area for each region — sum of piece areas,
    /// where each piece's area is its outer area minus its holes.
    pub fn areas(&self) -> HashMap<Combination, f64> {
        self.regions
            .iter()
            .map(|(combo, pieces)| {
                let area = pieces.iter().map(|p| p.area()).sum();
                (combo.clone(), area)
            })
            .collect()
    }

    /// Computes a hole-aware label anchor point for every non-empty region.
    ///
    /// For each region, this returns the *pole of inaccessibility* — the
    /// interior point that maximises the minimum distance to **every**
    /// boundary, both the outer ring and any holes (the Polylabel algorithm
    /// extended to hole-bearing shapes). When a region is disconnected into
    /// multiple [`RegionPiece`]s, the piece with the largest pole-to-boundary
    /// clearance wins — i.e. the piece with the most breathing room for a
    /// label, which is not always the piece with the largest area (a thin
    /// crescent can have large area but no spot wide enough to fit text).
    ///
    /// Regions with no pieces (or whose pieces are degenerate) are omitted.
    ///
    /// # When to use this vs [`crate::plotting::PlotData::region_anchors`]
    ///
    /// They compute the same thing. Use this method when you've built a
    /// `RegionPolygons` directly (e.g. via [`decompose_regions`]) and want
    /// region anchors without going through `PlotData`. Use
    /// [`PlotData::region_anchors`](crate::plotting::PlotData::region_anchors)
    /// if you're already constructing a [`crate::plotting::PlotData`] for
    /// rendering — the anchors are precomputed there.
    ///
    /// # Arguments
    ///
    /// * `precision` - Polylabel precision, in the same units as the polygon
    ///   coordinates. Smaller values yield more accurate anchors at higher
    ///   cost. A value of roughly 1% of the diagram's extent is typical
    ///   (e.g. `0.01` when coordinates are normalized to `[0, 1]`).
    ///
    /// # Examples
    ///
    /// ```
    /// use eunoia::{DiagramSpecBuilder, Fitter, InputType};
    /// use eunoia::geometry::shapes::Circle;
    ///
    /// let spec = DiagramSpecBuilder::new()
    ///     .set("A", 5.0)
    ///     .set("B", 3.0)
    ///     .intersection(&["A", "B"], 1.0)
    ///     .input_type(InputType::Exclusive)
    ///     .build()
    ///     .unwrap();
    ///
    /// let layout = Fitter::<Circle>::new(&spec).seed(42).fit().unwrap();
    /// let regions = layout.region_polygons(&spec, 64);
    /// let labels = regions.label_points(0.01);
    ///
    /// // One label point per non-empty region.
    /// for (combo, point) in &labels {
    ///     println!("{}: ({:.2}, {:.2})", combo, point.x(), point.y());
    /// }
    /// ```
    pub fn label_points(&self, precision: f64) -> HashMap<Combination, Point> {
        self.regions
            .iter()
            .filter_map(|(combo, pieces)| {
                let (point, _) = poi_with_holes(pieces, precision)?;
                Some((combo.clone(), point))
            })
            .collect()
    }

    /// **Raw** per-set label anchors — the hole-aware pole of inaccessibility
    /// of the union of every region containing the set. No fallback is
    /// applied: a set with no exclusive area of its own will land at the
    /// best POI of whatever pieces it appears in (which may be a thin
    /// intersection sliver), or be omitted entirely if it appears in no
    /// region.
    ///
    /// Most callers want
    /// [`PlotData::set_anchors`](crate::plotting::PlotData::set_anchors)
    /// instead, which adds the eulerr-style fallback chain (largest
    /// containing region → shape POI). Use this method only when you want
    /// the raw region-derived anchor without fallback — e.g. if you're
    /// implementing a different placement policy on top.
    ///
    /// Sets that are absent from every region are omitted from the result.
    ///
    /// # Arguments
    ///
    /// * `set_names` - The full list of set names in the diagram (typically
    ///   `spec.set_names()`).
    /// * `precision` - Polylabel precision, in the same units as the polygon
    ///   coordinates.
    ///
    /// # Examples
    ///
    /// ```
    /// use eunoia::{DiagramSpecBuilder, Fitter, InputType};
    /// use eunoia::geometry::shapes::Circle;
    ///
    /// let spec = DiagramSpecBuilder::new()
    ///     .set("A", 5.0)
    ///     .set("B", 3.0)
    ///     .intersection(&["A", "B"], 1.0)
    ///     .input_type(InputType::Exclusive)
    ///     .build()
    ///     .unwrap();
    ///
    /// let layout = Fitter::<Circle>::new(&spec).seed(42).fit().unwrap();
    /// let regions = layout.region_polygons(&spec, 64);
    /// let set_labels = regions.set_label_points(spec.set_names(), 0.01);
    ///
    /// assert!(set_labels.contains_key("A"));
    /// assert!(set_labels.contains_key("B"));
    /// ```
    pub fn set_label_points<N: AsRef<str>>(
        &self,
        set_names: &[N],
        precision: f64,
    ) -> HashMap<String, Point> {
        let mut result = HashMap::new();

        for name in set_names {
            let name = name.as_ref();
            // Flatten every piece (outer + holes) of every region containing
            // this set into a single ring list, union them, and reclassify
            // into pieces — gives the hole-aware POI of the unioned coverage.
            let mut owned: Vec<Polygon> = Vec::new();
            for (combo, pieces) in self.regions.iter() {
                if !combo.sets().iter().any(|s| s == name) {
                    continue;
                }
                for piece in pieces {
                    owned.push(piece.outer.clone());
                    owned.extend(piece.holes.iter().cloned());
                }
            }

            if owned.is_empty() {
                continue;
            }

            let mut merged = vec![owned.remove(0)];
            for p in owned {
                merged = polygon_clip_many(&merged, &p, ClipOperation::Union);
                if merged.is_empty() {
                    break;
                }
            }

            let pieces = classify_into_pieces(merged);
            if let Some((point, _)) = poi_with_holes(&pieces, precision) {
                result.insert(name.to_string(), point);
            }
        }

        result
    }
}

/// Signed area of a polygon ring via the shoelace formula. Positive = CCW,
/// negative = CW.
fn signed_polygon_area(p: &Polygon) -> f64 {
    let v = p.vertices();
    let n = v.len();
    if n < 3 {
        return 0.0;
    }
    let mut s = 0.0;
    for i in 0..n {
        let j = (i + 1) % n;
        s += v[i].x() * v[j].y() - v[j].x() * v[i].y();
    }
    0.5 * s
}

pub(crate) fn point_in_polygon(p: &Point, poly: &Polygon) -> bool {
    let v = poly.vertices();
    let n = v.len();
    if n < 3 {
        return false;
    }
    let mut inside = false;
    let mut j = n - 1;
    for i in 0..n {
        let (xi, yi) = (v[i].x(), v[i].y());
        let (xj, yj) = (v[j].x(), v[j].y());
        let intersect = ((yi > p.y()) != (yj > p.y()))
            && (p.x() < (xj - xi) * (p.y() - yi) / (yj - yi + f64::EPSILON) + xi);
        if intersect {
            inside = !inside;
        }
        j = i;
    }
    inside
}

/// Reverse the vertex order of a polygon, flipping its winding.
fn reverse_polygon(p: &Polygon) -> Polygon {
    let mut v = p.vertices().to_vec();
    v.reverse();
    Polygon::new(v)
}

/// Classify a flat list of clipper-output rings into [`RegionPiece`]s by
/// topological containment. Each ring's depth is the number of other rings
/// strictly containing one of its vertices: even depth → outer, odd depth →
/// hole. Holes are assigned to the smallest enclosing outer.
///
/// Outputs are normalised so `outer` is CCW (positive signed area) and every
/// `hole` is CW (negative signed area). Degenerate rings (fewer than three
/// vertices, or near-zero area) are dropped.
///
/// We classify by containment rather than winding because i_overlay's output
/// orientation is not stable across operations — relying on the sign of the
/// signed area to identify outers vs holes silently breaks for chained
/// difference/union operations.
///
/// # When to use this
///
/// Reach for this when you've chained your own [`crate::plotting::polygon_clip`]
/// calls and have a flat `Vec<Polygon>` of rings that you'd like to render
/// with the same outer-and-holes contract that
/// [`decompose_regions`] guarantees. For the standard "fitted shapes →
/// per-region pieces" pipeline, prefer [`decompose_regions`] (or
/// [`crate::Layout::region_polygons`]) — it handles boundary discovery,
/// inclusion-exclusion, and sliver filtering for you.
///
/// # Examples
///
/// ```
/// use eunoia::geometry::primitives::Point;
/// use eunoia::geometry::shapes::Polygon;
/// use eunoia::plotting::classify_into_pieces;
///
/// // A square with a square hole — two rings, same direction here
/// // (classification is by containment, not winding).
/// let outer = Polygon::new(vec![
///     Point::new(0.0, 0.0),
///     Point::new(4.0, 0.0),
///     Point::new(4.0, 4.0),
///     Point::new(0.0, 4.0),
/// ]);
/// let hole = Polygon::new(vec![
///     Point::new(1.0, 1.0),
///     Point::new(3.0, 1.0),
///     Point::new(3.0, 3.0),
///     Point::new(1.0, 3.0),
/// ]);
///
/// let pieces = classify_into_pieces(vec![outer, hole]);
/// assert_eq!(pieces.len(), 1);
/// assert_eq!(pieces[0].holes.len(), 1);
/// // Net area: 16 (outer) - 4 (hole) = 12.
/// assert!((pieces[0].area() - 12.0).abs() < 1e-9);
/// ```
pub fn classify_into_pieces(rings: Vec<Polygon>) -> Vec<RegionPiece> {
    // Drop degenerate rings.
    let kept: Vec<(Polygon, f64)> = rings
        .into_iter()
        .filter_map(|p| {
            let area = signed_polygon_area(&p).abs();
            if p.vertices().len() < 3 || area < 1e-12 {
                None
            } else {
                Some((p, area))
            }
        })
        .collect();
    if kept.is_empty() {
        return Vec::new();
    }

    let n = kept.len();
    let mut containment_depth = vec![0usize; n];
    let mut smallest_container: Vec<Option<usize>> = vec![None; n];
    for i in 0..n {
        let probe = match kept[i].0.vertices().first() {
            Some(p) => *p,
            None => continue,
        };
        for (j, (other, area_j)) in kept.iter().enumerate() {
            if i == j {
                continue;
            }
            if point_in_polygon(&probe, other) {
                containment_depth[i] += 1;
                let cur_area = smallest_container[i]
                    .and_then(|k| kept.get(k).map(|(_, a)| *a))
                    .unwrap_or(f64::INFINITY);
                if *area_j < cur_area {
                    smallest_container[i] = Some(j);
                }
            }
        }
    }

    let mut outer_idx_to_piece_idx: Vec<Option<usize>> = vec![None; n];
    let mut pieces: Vec<RegionPiece> = Vec::new();
    for (i, (ring, _)) in kept.iter().enumerate() {
        if containment_depth[i].is_multiple_of(2) {
            outer_idx_to_piece_idx[i] = Some(pieces.len());
            // Normalise outer to CCW (positive signed area).
            let outer = if signed_polygon_area(ring) >= 0.0 {
                ring.clone()
            } else {
                reverse_polygon(ring)
            };
            pieces.push(RegionPiece {
                outer,
                holes: Vec::new(),
            });
        }
    }
    for (i, (ring, _)) in kept.iter().enumerate() {
        if containment_depth[i] % 2 == 1
            && let Some(parent_idx) = smallest_container[i].and_then(|k| outer_idx_to_piece_idx[k])
        {
            // Normalise holes to CW (negative signed area).
            let hole = if signed_polygon_area(ring) <= 0.0 {
                ring.clone()
            } else {
                reverse_polygon(ring)
            };
            pieces[parent_idx].holes.push(hole);
        }
    }
    pieces
}

/// Minimum Euclidean distance from `(px, py)` to any segment in any of the
/// given polygon rings. Returns `f64::INFINITY` when every ring is degenerate
/// (fewer than two vertices).
///
/// Shared between [`poi_with_holes`] and the inscribed-rectangle search in
/// [`crate::plotting::inscribed`]; both need the same per-point clearance
/// calculation against an outer ring plus its holes.
pub(crate) fn min_dist_to_rings(px: f64, py: f64, rings: &[&[Point]]) -> f64 {
    let mut best = f64::INFINITY;
    for ring in rings {
        let n = ring.len();
        if n < 2 {
            continue;
        }
        for i in 0..n {
            let a = ring[i];
            let b = ring[(i + 1) % n];
            let dx = b.x() - a.x();
            let dy = b.y() - a.y();
            let len2 = dx * dx + dy * dy;
            let t = if len2 > 0.0 {
                (((px - a.x()) * dx + (py - a.y()) * dy) / len2).clamp(0.0, 1.0)
            } else {
                0.0
            };
            let qx = a.x() + t * dx;
            let qy = a.y() + t * dy;
            let d = ((px - qx).powi(2) + (py - qy).powi(2)).sqrt();
            if d < best {
                best = d;
            }
        }
    }
    best
}

/// Signed clearance of `(px, py)` to a region piece. Positive when the point
/// lies inside `piece.outer` and outside every hole; negative otherwise.
/// Magnitude is the minimum Euclidean distance to any ring.
pub(crate) fn signed_clearance(px: f64, py: f64, piece: &RegionPiece) -> f64 {
    let mut all: Vec<&[Point]> = Vec::with_capacity(1 + piece.holes.len());
    all.push(piece.outer.vertices());
    for h in &piece.holes {
        all.push(h.vertices());
    }
    let dist = min_dist_to_rings(px, py, &all);
    let probe = Point::new(px, py);
    let in_outer = point_in_polygon(&probe, &piece.outer);
    let in_any_hole = piece.holes.iter().any(|h| point_in_polygon(&probe, h));
    if in_outer && !in_any_hole {
        dist
    } else {
        -dist
    }
}

/// Compute the pole of inaccessibility of a region described as a list of
/// already-classified [`RegionPiece`]s (one outer + holes per connected
/// component). Returns the (point, clearance) of the piece whose POI has
/// the largest minimum distance to its outer boundary and any of its holes.
///
/// Implemented in-house rather than via `polylabel-mini` because the latter
/// scores cells using distance to the exterior ring only — it ignores hole
/// boundaries when ranking candidates, so a point sitting right next to a
/// hole gets credited with the same clearance as a point far from any hole.
///
/// Returns `None` when there are no pieces with positive area or every
/// search-cell clearance comes out non-positive (degenerate input).
pub(crate) fn poi_with_holes(pieces: &[RegionPiece], precision: f64) -> Option<(Point, f64)> {
    let mut best_overall: Option<(Point, f64)> = None;
    for piece in pieces {
        let outer_verts = piece.outer.vertices();
        if outer_verts.len() < 3 {
            continue;
        }

        // Bounding box of the outer ring.
        let mut min_x = f64::INFINITY;
        let mut min_y = f64::INFINITY;
        let mut max_x = f64::NEG_INFINITY;
        let mut max_y = f64::NEG_INFINITY;
        for p in outer_verts {
            min_x = min_x.min(p.x());
            max_x = max_x.max(p.x());
            min_y = min_y.min(p.y());
            max_y = max_y.max(p.y());
        }
        let width = max_x - min_x;
        let height = max_y - min_y;
        if width <= 0.0 || height <= 0.0 {
            continue;
        }
        let cell_size = width.min(height);

        // Priority queue keyed by upper bound on a cell's clearance.
        #[derive(Copy, Clone)]
        struct Cell {
            x: f64,
            y: f64,
            h: f64,
            d: f64,  // signed clearance at the centre
            ub: f64, // upper bound on clearance anywhere in the cell
        }
        impl Eq for Cell {}
        impl PartialEq for Cell {
            fn eq(&self, other: &Self) -> bool {
                self.ub == other.ub
            }
        }
        impl PartialOrd for Cell {
            fn partial_cmp(&self, other: &Self) -> Option<std::cmp::Ordering> {
                Some(self.cmp(other))
            }
        }
        impl Ord for Cell {
            fn cmp(&self, other: &Self) -> std::cmp::Ordering {
                self.ub
                    .partial_cmp(&other.ub)
                    .unwrap_or(std::cmp::Ordering::Equal)
            }
        }

        let make_cell = |x: f64, y: f64, h: f64| -> Cell {
            let d = signed_clearance(x, y, piece);
            Cell {
                x,
                y,
                h,
                d,
                ub: d + h * std::f64::consts::SQRT_2,
            }
        };

        let mut queue: std::collections::BinaryHeap<Cell> = std::collections::BinaryHeap::new();
        let mut x = min_x;
        let h0 = cell_size / 2.0;
        while x < max_x {
            let mut y = min_y;
            while y < max_y {
                queue.push(make_cell(x + h0, y + h0, h0));
                y += cell_size;
            }
            x += cell_size;
        }
        let cx = (min_x + max_x) / 2.0;
        let cy = (min_y + max_y) / 2.0;
        let mut best = make_cell(cx, cy, 0.0);

        while let Some(cell) = queue.pop() {
            if cell.d > best.d {
                best = cell;
            }
            if cell.ub - best.d <= precision {
                continue;
            }
            let nh = cell.h / 2.0;
            queue.push(make_cell(cell.x - nh, cell.y - nh, nh));
            queue.push(make_cell(cell.x + nh, cell.y - nh, nh));
            queue.push(make_cell(cell.x - nh, cell.y + nh, nh));
            queue.push(make_cell(cell.x + nh, cell.y + nh, nh));
        }

        if best.d > 0.0 {
            let pt = Point::new(best.x, best.y);
            if best_overall.map(|(_, d)| best.d > d).unwrap_or(true) {
                best_overall = Some((pt, best.d));
            }
        }
    }
    best_overall
}

// Intentionally no `Default` impl: an empty `RegionPolygons` is not a
// meaningful default for a fitted diagram, and exposing one would let
// callers bypass the [`decompose_regions`] invariants. Construct via
// [`decompose_regions`] (or [`crate::Layout::region_polygons`]).

/// Decomposes fitted shapes into exclusive [`RegionPiece`]s, one entry per
/// non-empty set combination.
///
/// # Algorithm
///
/// 1. Convert each shape to a polygon at `n_vertices` resolution.
/// 2. Discover candidate region masks sparsely from the actual fitted
///    geometry (via the shape-boundary intersection points), instead of
///    enumerating the full `2^n - 1` power set. Only regions that can
///    geometrically be non-empty are considered, so the cost scales with
///    the number of geometrically real regions and large `n` stays
///    tractable on sparse layouts.
/// 3. For each candidate, intersect the in-set polygons and subtract the
///    out-of-set polygons via repeated polygon clipping.
/// 4. Run the resulting flat ring list through the topological piece
///    classifier — each piece gets one CCW outer + zero or more CW holes,
///    with orientations normalised.
/// 5. Drop sliver pieces (see "Sliver filtering" below) and store the rest.
///
/// # Output guarantees
///
/// For every retained piece:
///
/// * `outer.area() > 0`, vertex count ≥ 3, winding is CCW.
/// * Each `hole` is strictly inside `outer`, with CW winding.
/// * Net area (`outer` minus `holes`) is positive and at least
///   `1e-3 ×` the largest piece in the diagram.
///
/// Renderers can therefore use SVG / Canvas defaults (`fill-rule: nonzero`)
/// and don't need to do their own ring classification — see [`RegionPiece`].
///
/// # Sliver filtering
///
/// 200-vertex polygonization (the typical default) produces tiny rounding
/// artifacts along seams between two shapes — e.g. a "B-only" region of
/// area ~0.02% of the largest piece when B is geometrically nested in A
/// but their polygonal boundaries don't agree exactly. These slivers
/// would otherwise mislead label anchors, produce duplicate strokes, and
/// expose binding authors to noise. The threshold is **net piece area
/// below `1e-3 ×` the largest single piece in the diagram**; smaller
/// pieces are dropped before insertion. This is scale-invariant (large
/// and small diagrams behave the same) and stable against `n_vertices`.
///
/// Driven by the actual fitted shapes, not the spec, so unexpected
/// overlaps the optimiser produces are still visualised.
///
/// # Complement region
///
/// When `spec.complement()` is `Some(_)` and a fitted `container` is passed,
/// the result also includes the **complement region** — the area inside
/// the container but outside every fitted shape — keyed on the empty
/// [`Combination`] (`Combination::new(&[])`). This mirrors how
/// `PreprocessedSpec::exclusive_areas` treats mask 0 as a first-class
/// region in the loss path. Renderers can identify it with
/// [`Combination::is_empty`] (or via the `""` string from
/// [`Combination::Display`](Combination)).
///
/// The complement piece is computed via
/// [`crate::plotting::polygon_difference`] (a single-pass multi-shape
/// difference, so outer/hole pairing is preserved when the union of
/// shapes is multi-ring inside the container) and runs through the same
/// piece classifier and sliver filter as every other region.
///
/// # Arguments
///
/// * `shapes` - The fitted diagram shapes (one per set).
/// * `set_names` - Set names in the same order as `shapes`.
/// * `spec` - The diagram specification. Read for `spec.complement()`; its
///   set-area entries are otherwise reserved for future use such as
///   spec-driven hole hints.
/// * `container` - The jointly-fitted container rectangle from
///   `Layout::container()`, or `None` when the spec has no complement.
///   Required for emitting the complement region; ignored otherwise.
/// * `n_vertices` - Resolution used when polygonizing each shape. Higher =
///   smoother boundaries at higher clipping cost. The eulerr default is 200.
///
/// # Examples
///
/// ```
/// use eunoia::{DiagramSpecBuilder, Fitter, InputType};
/// use eunoia::geometry::shapes::Circle;
/// use eunoia::plotting::decompose_regions;
///
/// let spec = DiagramSpecBuilder::new()
///     .set("A", 5.0)
///     .set("B", 3.0)
///     .intersection(&["A", "B"], 1.0)
///     .input_type(InputType::Exclusive)
///     .build()
///     .unwrap();
///
/// let layout = Fitter::<Circle>::new(&spec).fit().unwrap();
/// let shapes: Vec<_> = spec.set_names().iter()
///     .map(|name| *layout.shape_for_set(name).unwrap())
///     .collect();
///
/// let regions = decompose_regions(&shapes, spec.set_names(), &spec, layout.container(), 64);
/// ```
pub fn decompose_regions<S: DiagramShape + Polygonize>(
    shapes: &[S],
    set_names: &[String],
    spec: &DiagramSpec,
    container: Option<&Rectangle>,
    n_vertices: usize,
) -> RegionPolygons {
    decompose_regions_with(
        shapes,
        set_names,
        spec,
        container,
        n_vertices,
        DEFAULT_SLIVER_THRESHOLD,
    )
}

/// Default sliver threshold used by [`decompose_regions`] when called
/// without an explicit override. Mirrors `PlotOptions::default().sliver_threshold`.
pub const DEFAULT_SLIVER_THRESHOLD: f64 = 1e-3;

/// Like [`decompose_regions`] but with explicit control over the sliver
/// filter. `sliver_threshold` is the minimum net piece area relative to
/// the largest piece in the diagram; pieces below it are dropped. Set to
/// `0.0` to disable filtering.
pub fn decompose_regions_with<S: DiagramShape + Polygonize>(
    shapes: &[S],
    set_names: &[String],
    spec: &DiagramSpec,
    container: Option<&Rectangle>,
    n_vertices: usize,
    sliver_threshold: f64,
) -> RegionPolygons {
    if shapes.is_empty() {
        return RegionPolygons::new();
    }

    // Convert all shapes to polygons
    let shape_polygons: Vec<Polygon> = shapes.iter().map(|s| s.polygonize(n_vertices)).collect();

    let mut raw: Vec<(Combination, Vec<RegionPiece>)> = Vec::new();

    // Complement region (mask 0). Only emitted when the spec asked for it
    // *and* the layout supplied a fitted container — no inferring a
    // bounding box from the shapes, since that would silently change the
    // complement's size depending on the diagram extent. The piece is
    // computed in one pass via `polygon_difference` so outer/hole pairing
    // survives when the union of shapes is multi-ring inside the container.
    if spec.complement().is_some()
        && let Some(container_rect) = container
    {
        let container_poly = container_rect.polygonize(n_vertices);
        let rings = polygon_difference(&container_poly, &shape_polygons);
        let pieces = classify_into_pieces(rings);
        if !pieces.is_empty() {
            raw.push((Combination::new(&[]), pieces));
        }
    }

    let n = shapes.len();

    // Sparse region discovery: walk only candidate masks the actual fitted
    // geometry can populate. Avoids the 2^n power-set scan that the previous
    // implementation required and that put a hard practical ceiling on `n`.
    let intersections = collect_intersections_generic(shapes, n);
    let mut region_masks = discover_regions(shapes, &intersections, n);
    region_masks.sort_unstable();

    for mask in region_masks {
        let set_indices_in_combo: Vec<usize> = (0..n).filter(|&i| (mask >> i) & 1 == 1).collect();

        if set_indices_in_combo.is_empty() {
            continue;
        }

        // Start with the first shape in the combination
        let mut current_polygons = vec![shape_polygons[set_indices_in_combo[0]].clone()];

        // Intersect with remaining shapes that should be present
        for &idx in &set_indices_in_combo[1..] {
            current_polygons = polygon_clip_many(
                &current_polygons,
                &shape_polygons[idx],
                ClipOperation::Intersection,
            );

            if current_polygons.is_empty() {
                break;
            }
        }

        // Skip if intersection is empty
        if current_polygons.is_empty() {
            continue;
        }

        // Subtract all shapes that should NOT be present
        for (idx, _) in shapes.iter().enumerate() {
            if !set_indices_in_combo.contains(&idx) {
                current_polygons = polygon_clip_many(
                    &current_polygons,
                    &shape_polygons[idx],
                    ClipOperation::Difference,
                );

                if current_polygons.is_empty() {
                    break;
                }
            }
        }

        // Classify the clipper output (flat list of rings) into connected
        // pieces with explicit outer + holes structure. Skip empty.
        let pieces = classify_into_pieces(current_polygons);
        if pieces.is_empty() {
            continue;
        }

        let combo_sets: Vec<&str> = set_indices_in_combo
            .iter()
            .map(|&i| set_names[i].as_str())
            .collect();
        let combination = Combination::new(&combo_sets);

        raw.push((combination, pieces));
    }

    // Sliver filtering: drop pieces below `sliver_threshold` of the largest
    // piece anywhere in the diagram. Eliminates tiny shards produced by
    // 200-vertex polygonization along seams between geometrically tangent
    // shapes, which would otherwise mislead label anchors and stroke
    // rendering. `sliver_threshold == 0.0` disables filtering.
    let max_piece_area = raw
        .iter()
        .flat_map(|(_, pieces)| pieces.iter().map(|p| p.area()))
        .fold(0.0_f64, f64::max);
    let min_keep = if sliver_threshold > 0.0 {
        max_piece_area * sliver_threshold
    } else {
        0.0
    };

    let mut result = RegionPolygons::new();
    for (combo, pieces) in raw {
        let kept: Vec<RegionPiece> = pieces
            .into_iter()
            .filter(|p| p.area() >= min_keep)
            .collect();
        if !kept.is_empty() {
            result.insert(combo, kept);
        }
    }

    result
}

/// Collect pairwise boundary intersection points across arbitrary closed
/// shapes. Uses only the `Closed` trait, so it works for any `DiagramShape`
/// implementor — including future shapes that don't have a hand-rolled
/// `collect_intersections_*` helper.
///
/// The resulting `IntersectionPoint` list is what `discover_regions` consumes
/// to figure out which region masks the geometry can actually populate.
fn collect_intersections_generic<S: Closed>(shapes: &[S], n_sets: usize) -> Vec<IntersectionPoint> {
    let mut intersections = Vec::new();

    for i in 0..n_sets {
        for j in (i + 1)..n_sets {
            let pts = shapes[i].intersection_points(&shapes[j]);
            for point in pts {
                let mut adopters = vec![i, j];
                for (k, shape) in shapes.iter().enumerate() {
                    if k != i && k != j && shape.contains_point(&point) {
                        adopters.push(k);
                    }
                }
                adopters.sort_unstable();

                intersections.push(IntersectionPoint::new(point, (i, j), adopters));
            }
        }
    }

    intersections
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::fitter::Fitter;
    use crate::geometry::primitives::Bounds;
    use crate::geometry::shapes::Circle;
    use crate::geometry::traits::BoundingBox;
    use crate::spec::{DiagramSpecBuilder, InputType};

    #[test]
    fn test_decompose_disjoint_circles_skips_pairwise() {
        // Two disjoint circles. The previous power-set scanner would still
        // walk the {A,B} mask and only discover emptiness via polygon
        // clipping. The sparse path should skip {A,B} entirely (no boundary
        // intersection points, no containment), keeping cost proportional to
        // the number of geometrically real regions.
        let circles = [
            Circle::new(Point::new(0.0, 0.0), 1.0),
            Circle::new(Point::new(10.0, 0.0), 1.0),
        ];
        let set_names = vec!["A".to_string(), "B".to_string()];
        let spec = DiagramSpecBuilder::new()
            .set("A", std::f64::consts::PI)
            .set("B", std::f64::consts::PI)
            .input_type(InputType::Exclusive)
            .build()
            .unwrap();

        let regions = decompose_regions(&circles, &set_names, &spec, None, 64);

        // Exactly the two singleton regions, never the pair.
        assert_eq!(regions.len(), 2, "expected only A-only and B-only regions");
        assert!(regions.get(&Combination::new(&["A"])).is_some());
        assert!(regions.get(&Combination::new(&["B"])).is_some());
        assert!(
            regions.get(&Combination::new(&["A", "B"])).is_none(),
            "disjoint pair should not appear in region polygons"
        );
    }

    #[test]
    fn test_decompose_many_disjoint_circles_scales_sparsely() {
        // 20 disjoint circles in a row. Power-set decomposition would walk
        // 2^20 - 1 ≈ 1M masks; the sparse path walks only the 20 singletons.
        // This test asserts both correctness (only singletons appear) and
        // implicitly that the call returns in well under a second.
        let n = 20;
        let circles: Vec<Circle> = (0..n)
            .map(|i| Circle::new(Point::new(10.0 * i as f64, 0.0), 1.0))
            .collect();
        let set_names: Vec<String> = (0..n).map(|i| format!("S{i}")).collect();

        let mut builder = DiagramSpecBuilder::new();
        for name in &set_names {
            builder = builder.set(name.clone(), std::f64::consts::PI);
        }
        let spec = builder.input_type(InputType::Exclusive).build().unwrap();

        let regions = decompose_regions(&circles, &set_names, &spec, None, 32);
        assert_eq!(regions.len(), n, "expected one region per disjoint circle");
    }

    #[test]
    fn test_decompose_two_squares() {
        // The plotting path is generic over any `DiagramShape + Polygonize`.
        // Square implements both, so polygon decomposition into per-region
        // polygons should work without any plotting-side changes.
        use crate::geometry::shapes::Square;

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

        let layout = Fitter::<Square>::new(&spec).seed(42).fit().unwrap();
        let shapes: Vec<Square> = spec
            .set_names()
            .iter()
            .map(|name| *layout.shape_for_set(name).unwrap())
            .collect();

        // n_vertices is ignored by Square::polygonize (always 4 corners).
        let regions = decompose_regions(&shapes, spec.set_names(), &spec, None, 64);
        assert!(!regions.is_empty(), "no regions decomposed");
        for (combo, polys) in regions.iter() {
            assert!(!polys.is_empty(), "Region {combo:?} should have polygons");
        }
    }

    #[test]
    fn test_decompose_two_circles() {
        let spec = DiagramSpecBuilder::new()
            .set("A", 5.0)
            .set("B", 3.0)
            .intersection(&["A", "B"], 1.0)
            .input_type(InputType::Exclusive)
            .build()
            .unwrap();

        let layout = Fitter::<Circle>::new(&spec).seed(42).fit().unwrap();

        let shapes: Vec<Circle> = spec
            .set_names()
            .iter()
            .map(|name| *layout.shape_for_set(name).unwrap())
            .collect();

        let regions = decompose_regions(&shapes, spec.set_names(), &spec, None, 64);

        // Should have regions for A-only, B-only, and A&B
        assert!(regions.len() >= 2); // At least A-only and B-only

        // Check that we have some polygons
        for (combo, polys) in regions.iter() {
            assert!(!polys.is_empty(), "Region {combo:?} should have polygons");
        }
    }

    #[test]
    fn test_decompose_three_circles() {
        let spec = DiagramSpecBuilder::new()
            .set("A", 4.0)
            .set("B", 4.0)
            .set("C", 4.0)
            .intersection(&["A", "B"], 1.0)
            .intersection(&["B", "C"], 1.0)
            .intersection(&["A", "C"], 1.0)
            .intersection(&["A", "B", "C"], 0.5)
            .input_type(InputType::Exclusive)
            .build()
            .unwrap();

        let layout = Fitter::<Circle>::new(&spec).seed(123).fit().unwrap();

        let shapes: Vec<Circle> = spec
            .set_names()
            .iter()
            .map(|name| *layout.shape_for_set(name).unwrap())
            .collect();

        let regions = decompose_regions(&shapes, spec.set_names(), &spec, None, 64);

        // Should have multiple regions
        assert!(regions.len() >= 3);
    }

    #[test]
    fn test_label_points_two_circles() {
        // Two overlapping circles should yield 3 non-empty regions: A-only,
        // B-only, A&B. Each should get a label point inside the region itself.
        let spec = DiagramSpecBuilder::new()
            .set("A", 5.0)
            .set("B", 3.0)
            .intersection(&["A", "B"], 1.0)
            .input_type(InputType::Exclusive)
            .build()
            .unwrap();

        let layout = Fitter::<Circle>::new(&spec).seed(42).fit().unwrap();
        let shapes: Vec<Circle> = spec
            .set_names()
            .iter()
            .map(|name| *layout.shape_for_set(name).unwrap())
            .collect();

        let regions = decompose_regions(&shapes, spec.set_names(), &spec, None, 64);
        let labels = regions.label_points(0.01);

        // Every non-empty region in `regions` should appear in `labels`.
        for combo in regions.iter().map(|(c, _)| c) {
            assert!(
                labels.contains_key(combo),
                "Missing label point for region {combo:?}"
            );
        }

        // Each label point must sit inside the axis-aligned bounding box of
        // its region's largest piece's outer ring.
        for (combo, pieces) in regions.iter() {
            let label = labels.get(combo).unwrap();
            let largest = pieces
                .iter()
                .max_by(|a, b| a.area().partial_cmp(&b.area()).unwrap())
                .unwrap();
            let (mut min_x, mut min_y) = (f64::INFINITY, f64::INFINITY);
            let (mut max_x, mut max_y) = (f64::NEG_INFINITY, f64::NEG_INFINITY);
            for v in largest.outer.vertices() {
                min_x = min_x.min(v.x());
                min_y = min_y.min(v.y());
                max_x = max_x.max(v.x());
                max_y = max_y.max(v.y());
            }
            assert!(
                label.x() >= min_x - 1e-9
                    && label.x() <= max_x + 1e-9
                    && label.y() >= min_y - 1e-9
                    && label.y() <= max_y + 1e-9,
                "Label for {:?} at ({:.3}, {:.3}) is outside its region's bounding box [{:.3}, {:.3}] x [{:.3}, {:.3}]",
                combo,
                label.x(),
                label.y(),
                min_x,
                max_x,
                min_y,
                max_y
            );
        }
    }

    #[test]
    fn test_set_label_points_two_circles() {
        let spec = DiagramSpecBuilder::new()
            .set("A", 5.0)
            .set("B", 3.0)
            .intersection(&["A", "B"], 1.0)
            .input_type(InputType::Exclusive)
            .build()
            .unwrap();

        let layout = Fitter::<Circle>::new(&spec).seed(42).fit().unwrap();
        let shapes: Vec<Circle> = spec
            .set_names()
            .iter()
            .map(|name| *layout.shape_for_set(name).unwrap())
            .collect();

        let regions = decompose_regions(&shapes, spec.set_names(), &spec, None, 128);
        let set_labels = regions.set_label_points(spec.set_names(), 0.01);

        assert_eq!(set_labels.len(), 2);
        for name in ["A", "B"] {
            let label = set_labels.get(name).expect("missing set label");
            // The set label must lie inside the corresponding circle.
            let circle = layout.shape_for_set(name).unwrap();
            let dx = label.x() - circle.center().x();
            let dy = label.y() - circle.center().y();
            let r = circle.radius();
            assert!(
                dx * dx + dy * dy <= r * r + 1e-6,
                "set label for {} at ({:.3}, {:.3}) is outside circle (center=({:.3}, {:.3}), r={:.3})",
                name,
                label.x(),
                label.y(),
                circle.center().x(),
                circle.center().y(),
                r,
            );
        }
    }

    #[test]
    fn test_set_label_points_skips_absent_sets() {
        let regions = RegionPolygons::new();
        let names = vec!["A".to_string()];
        assert!(regions.set_label_points(&names, 0.01).is_empty());
    }

    #[test]
    fn test_label_points_empty() {
        let empty = RegionPolygons::new();
        assert!(empty.label_points(0.01).is_empty());
    }

    #[test]
    fn test_label_points_skips_zero_area_regions() {
        // A region composed only of a degenerate (zero-area) outer should be
        // omitted from the label map.
        let mut regions = RegionPolygons::new();
        let degenerate = Polygon::new(vec![
            Point::new(0.0, 0.0),
            Point::new(1.0, 0.0),
            Point::new(0.5, 0.0), // collinear → zero area
        ]);
        let piece = RegionPiece {
            outer: degenerate,
            holes: Vec::new(),
        };
        regions.insert(Combination::new(&["X"]), vec![piece]);

        assert!(regions.label_points(0.01).is_empty());
    }

    #[test]
    fn test_region_areas() {
        let spec = DiagramSpecBuilder::new()
            .set("A", 5.0)
            .set("B", 3.0)
            .intersection(&["A", "B"], 1.0)
            .input_type(InputType::Exclusive)
            .build()
            .unwrap();

        let layout = Fitter::<Circle>::new(&spec).seed(42).fit().unwrap();

        let shapes: Vec<Circle> = spec
            .set_names()
            .iter()
            .map(|name| *layout.shape_for_set(name).unwrap())
            .collect();

        let regions = decompose_regions(&shapes, spec.set_names(), &spec, None, 128);
        let areas = regions.areas();

        // Total area should be close to sum of fitted values
        let total_area: f64 = areas.values().sum();
        let expected_total: f64 = spec.exclusive_areas().values().sum();

        // Allow some tolerance due to polygonization
        assert!(
            (total_area - expected_total).abs() < 0.5,
            "Total area {total_area:.3} should be close to expected {expected_total:.3}"
        );
    }

    #[test]
    fn test_iter_sorted_is_deterministic_by_size_then_lex() {
        let mut regions = RegionPolygons::new();
        let dummy_piece = || RegionPiece {
            outer: Polygon::new(vec![
                Point::new(0.0, 0.0),
                Point::new(1.0, 0.0),
                Point::new(1.0, 1.0),
            ]),
            holes: Vec::new(),
        };
        // Insert in scrambled order to make sure iter_sorted re-sorts.
        regions.insert(Combination::new(&["A", "B"]), vec![dummy_piece()]);
        regions.insert(Combination::new(&["B"]), vec![dummy_piece()]);
        regions.insert(Combination::new(&["A"]), vec![dummy_piece()]);
        regions.insert(Combination::new(&["A", "B", "C"]), vec![dummy_piece()]);
        regions.insert(Combination::new(&["A", "C"]), vec![dummy_piece()]);

        let order: Vec<String> = regions.iter_sorted().map(|(c, _)| c.to_string()).collect();
        assert_eq!(order, vec!["A", "B", "A&B", "A&C", "A&B&C"]);
    }

    #[test]
    fn test_iter_is_deterministic_and_matches_iter_sorted() {
        let mut regions = RegionPolygons::new();
        let dummy_piece = || RegionPiece {
            outer: Polygon::new(vec![
                Point::new(0.0, 0.0),
                Point::new(1.0, 0.0),
                Point::new(1.0, 1.0),
            ]),
            holes: Vec::new(),
        };
        // Insert in scrambled order; iter() must impose the canonical order
        // regardless of how (or in what hash-seeded order) entries are stored.
        regions.insert(Combination::new(&["A", "B"]), vec![dummy_piece()]);
        regions.insert(Combination::new(&["B"]), vec![dummy_piece()]);
        regions.insert(Combination::new(&["A"]), vec![dummy_piece()]);
        regions.insert(Combination::new(&["A", "B", "C"]), vec![dummy_piece()]);
        regions.insert(Combination::new(&["A", "C"]), vec![dummy_piece()]);

        // Canonical (size, lexicographic) order: singletons first, then pairs,
        // then triples; ties broken alphabetically.
        let order: Vec<String> = regions.iter().map(|(c, _)| c.to_string()).collect();
        assert_eq!(order, vec!["A", "B", "A&B", "A&C", "A&B&C"]);

        // iter() and iter_sorted() are the same order.
        let sorted: Vec<String> = regions.iter_sorted().map(|(c, _)| c.to_string()).collect();
        assert_eq!(order, sorted);
    }

    #[test]
    fn test_iter_in_input_order_respects_set_names_order() {
        let mut regions = RegionPolygons::new();
        let dummy_piece = || RegionPiece {
            outer: Polygon::new(vec![
                Point::new(0.0, 0.0),
                Point::new(1.0, 0.0),
                Point::new(1.0, 1.0),
            ]),
            holes: Vec::new(),
        };
        regions.insert(Combination::new(&["B"]), vec![dummy_piece()]);
        regions.insert(Combination::new(&["A"]), vec![dummy_piece()]);
        regions.insert(Combination::new(&["C"]), vec![dummy_piece()]);
        regions.insert(Combination::new(&["A", "B"]), vec![dummy_piece()]);
        regions.insert(Combination::new(&["A", "C"]), vec![dummy_piece()]);
        regions.insert(Combination::new(&["B", "C"]), vec![dummy_piece()]);

        // Spec defined in order C, A, B — iteration should follow that.
        let names = vec!["C".to_string(), "A".to_string(), "B".to_string()];
        let order: Vec<String> = regions
            .iter_in_input_order(&names)
            .map(|(c, _)| c.to_string())
            .collect();
        // Singletons first, in input order: C, A, B.
        // Pairs sorted by their canonical-form member indices in `names`:
        //   "A&C"  → indices (1, 0)  → sorted (0, 1)
        //   "B&C"  → indices (2, 0)  → sorted (0, 2)
        //   "A&B"  → indices (1, 2)  → sorted (1, 2)
        assert_eq!(order, vec!["C", "A", "B", "A&C", "B&C", "A&B"]);
    }

    #[test]
    fn test_decompose_with_zero_spec_area() {
        // Regression test: When spec has zero area for a set (e.g., C),
        // but the fitted layout has non-zero area, we should still generate
        // polygons for all possible regions, even if after subtraction some
        // regions are empty.
        let spec = DiagramSpecBuilder::new()
            .set("A", 3.0)
            .set("B", 5.0)
            .intersection(&["A", "B", "C"], 1.0)
            .input_type(InputType::Exclusive)
            .build()
            .unwrap();

        // Spec should have C with zero area
        let c_combo = crate::spec::Combination::new(&["C"]);
        assert!(
            spec.exclusive_areas().get(&c_combo).copied().unwrap_or(0.0) < 1e-10,
            "Spec should have zero area for C-only"
        );

        // Fit the layout - this will create a non-zero ellipse for C
        use crate::geometry::shapes::Ellipse;
        let layout = Fitter::<Ellipse>::new(&spec).seed(1).fit().unwrap();

        let shapes: Vec<Ellipse> = spec
            .set_names()
            .iter()
            .map(|name| *layout.shape_for_set(name).unwrap())
            .collect();

        // Decompose regions
        let regions = decompose_regions(&shapes, spec.set_names(), &spec, None, 64);

        // For this particular configuration, C is fully contained within
        // A&B&C intersection, so there won't be a C-only region. But we
        // should verify that regions involving C are present.
        let abc_combo = crate::spec::Combination::new(&["A", "B", "C"]);
        let abc_polygons = regions.get(&abc_combo);
        assert!(abc_polygons.is_some(), "Should have polygons for A&B&C");

        // Verify total area is reasonable
        let total_area: f64 = regions.areas().values().sum();

        // The total should be close to sum of individual shapes minus overlaps
        assert!(
            total_area > 5.0,
            "Total area should be substantial, got {total_area:.3}"
        );
    }

    #[test]
    fn test_complement_emitted_under_empty_combination_when_set() {
        // With a complement target, the fitter produces a container, and
        // `decompose_regions` should emit the complement piece keyed on
        // the empty Combination (mirroring how mask 0 is treated in the
        // loss path).
        let spec = DiagramSpecBuilder::new()
            .set("A", 25.0)
            .set("B", 25.0)
            .intersection(&["A", "B"], 5.0)
            .complement(20.0) // universe = 75
            .input_type(InputType::Exclusive)
            .build()
            .unwrap();

        let layout = Fitter::<Circle>::new(&spec).seed(42).fit().unwrap();
        let container = layout.container().expect("container fitted");

        let shapes: Vec<Circle> = spec
            .set_names()
            .iter()
            .map(|name| *layout.shape_for_set(name).unwrap())
            .collect();

        let regions = decompose_regions(&shapes, spec.set_names(), &spec, Some(container), 128);

        // The complement region appears under the empty combination.
        let empty = Combination::new(&[]);
        let pieces = regions
            .get(&empty)
            .expect("complement region missing under empty Combination");
        assert!(!pieces.is_empty());

        // Net complement area should be close to the spec target. We polygonise
        // shapes (so accuracy is bounded by `n_vertices`); allow 5% slack.
        let complement_area: f64 = pieces.iter().map(|p| p.area()).sum();
        assert!(
            (complement_area - 20.0).abs() / 20.0 < 0.05,
            "complement area {complement_area:.3} should be ≈ 20.0 (within 5%)"
        );

        // We don't pin down the topology (single piece with a hole vs.
        // several disjoint pieces) because the fitter is free to pack
        // shapes against the container boundary, which fragments the
        // complement. The area + label-anchor checks already cover what
        // bindings actually consume.

        // Label points must include an anchor inside the complement region.
        let labels = regions.label_points(0.05);
        let anchor = labels.get(&empty).expect("complement label anchor missing");
        // Anchor must lie inside the container.
        let Bounds {
            x_min: xmin,
            x_max: xmax,
            y_min: ymin,
            y_max: ymax,
        } = container.bounds();
        assert!(anchor.x() >= xmin && anchor.x() <= xmax);
        assert!(anchor.y() >= ymin && anchor.y() <= ymax);
    }

    #[test]
    fn test_complement_with_strictly_nested_shapes_yields_outer_plus_holes() {
        // Synthetic layout: two shapes strictly inside a generously sized
        // container, far from each other and from the container edges.
        // The complement must therefore be a single piece (the container
        // outline) with two holes (one per shape).
        use crate::geometry::primitives::Point;
        use crate::geometry::shapes::Rectangle;

        let spec = DiagramSpecBuilder::new()
            .set("A", std::f64::consts::PI)
            .set("B", std::f64::consts::PI)
            .complement(96.0)
            .input_type(InputType::Exclusive)
            .build()
            .unwrap();

        let shapes = [
            Circle::new(Point::new(-3.0, 0.0), 1.0),
            Circle::new(Point::new(3.0, 0.0), 1.0),
        ];
        // 10×10 box, two unit circles ⇒ complement ≈ 100 − 2π ≈ 93.7.
        let container = Rectangle::new(Point::new(0.0, 0.0), 10.0, 10.0);

        let regions = decompose_regions(&shapes, spec.set_names(), &spec, Some(&container), 128);

        let pieces = regions
            .get(&Combination::new(&[]))
            .expect("complement region missing");
        assert_eq!(pieces.len(), 1, "complement should be a single piece");
        assert_eq!(
            pieces[0].holes.len(),
            2,
            "expected exactly two holes, one per nested shape"
        );

        let net = pieces[0].area();
        let expected = 100.0 - 2.0 * std::f64::consts::PI;
        assert!(
            (net - expected).abs() < 0.1,
            "net complement area {net:.4} should be ≈ {expected:.4}"
        );
    }

    #[test]
    fn test_no_complement_region_when_spec_has_none() {
        // Without `.complement(...)`, the empty Combination must never be
        // emitted, even if a container ref is somehow passed.
        let spec = DiagramSpecBuilder::new()
            .set("A", 5.0)
            .set("B", 3.0)
            .intersection(&["A", "B"], 1.0)
            .input_type(InputType::Exclusive)
            .build()
            .unwrap();

        let layout = Fitter::<Circle>::new(&spec).seed(42).fit().unwrap();
        let shapes: Vec<Circle> = spec
            .set_names()
            .iter()
            .map(|name| *layout.shape_for_set(name).unwrap())
            .collect();

        let regions = decompose_regions(&shapes, spec.set_names(), &spec, None, 64);
        assert!(regions.get(&Combination::new(&[])).is_none());
    }

    #[test]
    fn test_complement_region_in_plot_data_shows_up_in_anchors_and_areas() {
        // End-to-end via `Layout::plot_data`: the complement combo string ""
        // must appear in `region_anchors`, `region_areas`, and
        // `pieces_for("")`.
        use crate::plotting::PlotOptions;

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

        let layout = Fitter::<Circle>::new(&spec).seed(42).fit().unwrap();
        let plot = layout.plot_data(&spec, PlotOptions::default());

        let key = Combination::new(&[]).to_string();
        assert_eq!(key, "", "empty combination must serialise to empty string");
        assert!(
            plot.region_anchors.contains_key(&key),
            "region_anchors should include the empty combo for complement"
        );
        assert!(
            plot.region_areas.contains_key(&key),
            "region_areas should include the empty combo for complement"
        );
        assert!(
            plot.pieces_for(&key).is_some(),
            "pieces_for(\"\") should return the complement region"
        );
    }
}
