u-nesting-d2 0.28.0

//! Shared utility functions for 2D nesting placement algorithms.
//!
//! This module consolidates common geometry operations used across
//! multiple nesting strategy implementations (GA, SA, BRKGA, ALNS, GDRR).

use crate::boundary::Boundary2D;
use crate::geometry::Geometry2D;
use crate::nfp::Nfp;
use std::sync::{Arc, Mutex};
use std::time::Duration;
use u_nesting_core::geometry::{Boundary, Geometry};
use u_nesting_core::timing::Timer;
use u_nesting_core::Placement;

/// Instance information for decoding placement orders.
///
/// Maps a flat instance index to the corresponding geometry and its
/// repetition number (when `quantity > 1`).
#[derive(Debug, Clone)]
pub struct InstanceInfo {
    /// Index into the geometries array.
    pub geometry_idx: usize,
    /// Instance number within this geometry's quantity.
    pub instance_num: usize,
}

/// Computes the centroid (arithmetic mean of vertices) of a polygon.
///
/// # Returns
///
/// `(0.0, 0.0)` for an empty polygon, otherwise the mean of all vertices.
pub fn polygon_centroid(polygon: &[(f64, f64)]) -> (f64, f64) {
    if polygon.is_empty() {
        return (0.0, 0.0);
    }

    let sum: (f64, f64) = polygon
        .iter()
        .fold((0.0, 0.0), |acc, &(x, y)| (acc.0 + x, acc.1 + y));
    let n = polygon.len() as f64;
    (sum.0 / n, sum.1 / n)
}

/// The no-fit polygons that keep `geom`, rotated by `rotation` and mirrored
/// when `mirror`, off every geometry already placed, each grown by `spacing`.
///
/// A no-fit polygon depends on the two shapes and the rotation between them,
/// not on where either lies, so one is computed between the placed geometry at
/// its origin and `geom` turned by the rotation relative to the placed one,
/// kept in `cache` under the two ids, and then turned by the placed rotation
/// and moved to the placed position. A search that places the same pieces in a
/// new order thousands of times finds nearly every pair in the cache; computing
/// each pair afresh is what a layout costs without it. Ids name shapes, so
/// they must be unique per shape within one cache (geometry ids are).
/// A pair whose no-fit polygon cannot be computed is left out.
pub(crate) fn nfps_against_placed(
    cache: &crate::nfp::NfpCache,
    placed: &[crate::nfp::PlacedGeometry],
    geom: &Geometry2D,
    rotation: f64,
    mirror: bool,
    spacing: f64,
) -> Vec<Nfp> {
    placed
        .iter()
        .filter_map(|p| {
            let relative = rotation - p.rotation;
            let at_origin = cache
                .get_or_compute_mirrored(
                    (
                        p.geometry.id().as_str(),
                        geom.id().as_str(),
                        relative,
                        p.mirrored,
                        mirror,
                    ),
                    || {
                        crate::nfp::compute_nfp_mirrored(
                            &p.geometry, geom, relative, p.mirrored, mirror,
                        )
                        // Grown before caching: an offset commutes with the
                        // rotation and translation below.
                        .map(|nfp| offset_nfp(&nfp, spacing))
                    },
                )
                .ok()?;
            Some(crate::nfp::translate_nfp(
                &crate::nfp::rotate_nfp(&at_origin, p.rotation),
                p.position,
            ))
        })
        .collect()
}

/// Grows an NFP by `spacing`, so that a reference point outside the result
/// keeps the two pieces at least `spacing` apart.
///
/// The NFP of pieces grown by `spacing / 2` each is the NFP grown by `spacing`
/// (the disc is symmetric), so this is the exact clearance region, computed as
/// a true polygon offset — see `polygon_ops::offset_polygon` for the accuracy
/// bound. Returns the NFP unchanged if `spacing <= 0.0`.
pub fn offset_nfp(nfp: &Nfp, spacing: f64) -> Nfp {
    if spacing <= 0.0 {
        return nfp.clone();
    }
    Nfp::from_polygons(
        nfp.polygons
            .iter()
            .flat_map(|polygon| crate::polygon_ops::offset_polygon(polygon, spacing))
            .collect(),
    )
}

/// The region pieces may occupy: the boundary moved inward by `margin`,
/// counter-clockwise.
///
/// A rectangle (or strip) is inset exactly on its bounding box. Any other
/// boundary is offset inward along its own edges — using its bounding box
/// instead would let a piece sit in the box's empty corners, closer than
/// `margin` to a slanted edge or outside the boundary altogether. If the
/// offset splits the boundary the largest piece is kept; if nothing is left,
/// the result is empty and nothing fits. Holes are not part of this region.
///
/// Placement code insets the boundary with this once and then computes the
/// inner-fit polygon with no further margin.
pub fn inset_boundary(boundary: &Boundary2D, margin: f64) -> Vec<(f64, f64)> {
    let rectangular =
        boundary.is_infinite() || (boundary.width().is_some() && boundary.height().is_some());
    if rectangular {
        let (b_min, b_max) = boundary.aabb();
        return vec![
            (b_min[0] + margin, b_min[1] + margin),
            (b_max[0] - margin, b_min[1] + margin),
            (b_max[0] - margin, b_max[1] - margin),
            (b_min[0] + margin, b_max[1] - margin),
        ];
    }
    let mut ring = crate::polygon_ops::offset_polygon(boundary.exterior(), -margin)
        .into_iter()
        .max_by(|a, b| {
            let area = |r: &[(f64, f64)]| crate::polygon_ops::signed_area2(r).abs();
            area(a).total_cmp(&area(b))
        })
        .unwrap_or_default();
    if crate::polygon_ops::signed_area2(&ring) < 0.0 {
        ring.reverse();
    }
    ring
}

/// Regions a piece's reference point must stay out of because of the
/// boundary's holes: for each hole, the no-fit polygon of the hole and the
/// piece (at `rotation`, mirrored if `mirror`), grown by `margin`.
///
/// Without these the placement search treats a hole as usable area and every
/// piece it puts there is later rejected by the boundary check.
pub fn hole_nfps(
    boundary: &Boundary2D,
    geometry: &Geometry2D,
    rotation: f64,
    mirror: bool,
    margin: f64,
) -> Vec<Nfp> {
    boundary
        .holes()
        .iter()
        .enumerate()
        .filter_map(|(i, hole)| {
            let stationary = Geometry2D::new(format!("_hole_{i}")).with_polygon(hole.clone());
            crate::nfp::compute_nfp_mirrored(&stationary, geometry, rotation, false, mirror).ok()
        })
        .map(|nfp| offset_nfp(&nfp, margin))
        .collect()
}

/// The best layout a search has produced, kept as the search goes.
#[derive(Debug, Clone)]
pub struct BestLayout {
    /// Fitness the layout scored.
    pub fitness: f64,
    /// The placements.
    pub placements: Vec<Placement<f64>>,
    /// Fraction of the boundary area the placements cover.
    pub utilization: f64,
}

/// What a search-based decoder shares with the run that owns it: when to stop
/// and the best layout seen so far.
///
/// A decoder that only watches the cancellation flag keeps placing past the
/// time limit until the whole layout is built, and a run that re-decodes its
/// best individual at the end pays for one more full layout after the limit.
/// Checking the deadline per placed piece and keeping the best decoded layout
/// removes both.
#[derive(Debug, Default)]
pub struct SearchState {
    budget: Option<(Timer, Duration)>,
    best: Arc<Mutex<Option<BestLayout>>>,
}

impl SearchState {
    /// Starts the clock for a search limited to `limit` (none = unlimited).
    pub fn with_time_limit(limit: Option<Duration>) -> Self {
        Self {
            budget: limit.map(|limit| (Timer::now(), limit)),
            best: Arc::default(),
        }
    }

    /// Whether the time limit has passed.
    pub fn out_of_time(&self) -> bool {
        self.budget
            .as_ref()
            .is_some_and(|(start, limit)| start.elapsed() > *limit)
    }

    /// Records a decoded layout if it beats the best one so far.
    ///
    /// Layouts are offered from parallel evaluations in whatever order the
    /// threads finish, so a tie on fitness is broken by the layout itself
    /// (`layout_order`), not by which arrived first -- otherwise a seeded
    /// run could return a different layout of the same fitness each time.
    pub fn offer(&self, fitness: f64, placements: &[Placement<f64>], utilization: f64) {
        let mut best = self
            .best
            .lock()
            .unwrap_or_else(|poisoned| poisoned.into_inner());
        if best.as_ref().is_none_or(|b| {
            fitness > b.fitness
                || (fitness == b.fitness && layout_order(placements, &b.placements).is_lt())
        }) {
            *best = Some(BestLayout {
                fitness,
                placements: placements.to_vec(),
                utilization,
            });
        }
    }

    /// A handle to the best layout, readable after the decoder is moved into a runner.
    pub fn best_handle(&self) -> Arc<Mutex<Option<BestLayout>>> {
        Arc::clone(&self.best)
    }
}

/// A total order on layouts that does not depend on how they were produced:
/// placement by placement, by geometry id, instance, position, rotation and
/// mirror flag.
fn layout_order(a: &[Placement<f64>], b: &[Placement<f64>]) -> std::cmp::Ordering {
    let floats = |x: &[f64], y: &[f64]| {
        x.iter()
            .zip(y)
            .map(|(p, q)| p.total_cmp(q))
            .find(|o| o.is_ne())
            .unwrap_or(x.len().cmp(&y.len()))
    };
    a.iter()
        .zip(b)
        .map(|(p, q)| {
            p.geometry_id
                .cmp(&q.geometry_id)
                .then(p.instance.cmp(&q.instance))
                .then_with(|| floats(&p.position, &q.position))
                .then_with(|| floats(&p.rotation, &q.rotation))
                .then(p.mirrored.cmp(&q.mirrored))
        })
        .find(|o| o.is_ne())
        .unwrap_or(a.len().cmp(&b.len()))
}

/// Takes the best layout out of a handle from [`SearchState::best_handle`].
pub fn take_best(handle: &Mutex<Option<BestLayout>>) -> Option<BestLayout> {
    handle
        .lock()
        .unwrap_or_else(|poisoned| poisoned.into_inner())
        .take()
}

/// Per-instance rotation index and mirror flag that reproduce `placements`
/// when pieces are placed in input order (each geometry's instances in turn).
///
/// This is how a search starts from the greedy layout: placing the instances
/// in that order with these genes makes the same choices the greedy pass made.
/// An instance the layout does not contain gets rotation 0, unmirrored.
pub fn seed_genes(
    geometries: &[Geometry2D],
    placements: &[Placement<f64>],
) -> (Vec<usize>, Vec<bool>) {
    let mut rotations = Vec::new();
    let mut mirrors = Vec::new();
    for geometry in geometries {
        let angles = geometry.rotations();
        for instance in 0..geometry.quantity() {
            let placed = placements
                .iter()
                .find(|p| &p.geometry_id == geometry.id() && p.instance == instance);
            let (rotation, mirrored) = placed.map_or((0, false), |p| {
                let angle = p.rotation.first().copied().unwrap_or(0.0);
                let turn = std::f64::consts::TAU;
                let index = angles
                    .iter()
                    .enumerate()
                    .min_by(|(_, a), (_, b)| {
                        let gap = |x: f64| {
                            (x - angle)
                                .rem_euclid(turn)
                                .min((angle - x).rem_euclid(turn))
                        };
                        gap(**a).total_cmp(&gap(**b))
                    })
                    .map_or(0, |(i, _)| i);
                (index, p.mirrored)
            });
            rotations.push(rotation);
            mirrors.push(mirrored);
        }
    }
    (rotations, mirrors)
}

/// The frame a strip is measured in: the boundary's length axis, where the strip
/// starts along it, and its extent across it.
pub fn strip_frame(boundary: &Boundary2D) -> (crate::nfp::PackingAxis, f64, f64) {
    let axis = crate::nfp::PackingAxis::of(boundary);
    let (b_min, b_max) = boundary.aabb();
    match axis {
        crate::nfp::PackingAxis::X => (axis, b_min[0], b_max[1] - b_min[1]),
        crate::nfp::PackingAxis::Y => (axis, b_min[1], b_max[0] - b_min[0]),
    }
}

/// How densely a layout fills the strip it uses: the area of the placed pieces
/// over the strip's cross extent times the length the pieces reach along the
/// boundary's length axis ([`crate::nfp::PackingAxis::of`]), counted from the boundary's
/// start edge. `0.0` when nothing is placed.
///
/// The area over the whole boundary is the same for every layout that places
/// the same pieces, so it cannot tell a short layout from a long one; this can,
/// and it is the density a strip is measured by.
pub fn layout_density(
    geometries: &[Geometry2D],
    placements: &[Placement<f64>],
    boundary: &Boundary2D,
) -> f64 {
    let (axis, start, cross) = strip_frame(boundary);
    let by_id: std::collections::HashMap<&str, &Geometry2D> =
        geometries.iter().map(|g| (g.id().as_str(), g)).collect();
    let (mut area, mut end) = (0.0, f64::NEG_INFINITY);
    for p in placements {
        let Some(geom) = by_id.get(p.geometry_id.as_str()) else {
            continue;
        };
        let rotation = p.rotation.first().copied().unwrap_or(0.0);
        let (_, hi) = geom.aabb_at_rotation_mirrored(rotation, p.mirrored);
        let k = match axis {
            crate::nfp::PackingAxis::X => 0,
            crate::nfp::PackingAxis::Y => 1,
        };
        area += geom.measure();
        end = end.max(p.position.get(k).copied().unwrap_or(0.0) + hi[k]);
    }
    let length = end - start;
    if area <= 0.0 || length <= 0.0 || cross <= 0.0 {
        return 0.0;
    }
    (area / (length * cross)).min(1.0)
}

/// Computes the nesting fitness score from placement results.
///
/// The fitness function prioritizes placement count (weight 100) over
/// utilization (weight 10), ensuring all-placed solutions always rank
/// higher than partial solutions.
///
/// # Arguments
///
/// * `placed_count` - Number of successfully placed instances
/// * `total_count` - Total number of instances to place
/// * `utilization` - Fraction of boundary area used (0.0..1.0)
///
/// # Returns
///
/// Fitness score in range [0, 110] where higher is better.
pub fn nesting_fitness(placed_count: usize, total_count: usize, utilization: f64) -> f64 {
    let placement_ratio = placed_count as f64 / total_count.max(1) as f64;
    placement_ratio * 100.0 + utilization * 10.0
}

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

    #[test]
    fn density_tells_a_short_layout_from_a_long_one() {
        // Two 10 x 10 squares on a 100 x 20 strip: side by side they use 20 of
        // its length, one behind the other 50 of it. The area over the whole
        // boundary is the same for both.
        let g = [Geometry2D::rectangle("s", 10.0, 10.0).with_quantity(2)];
        let strip = Boundary2D::rectangle(100.0, 20.0);
        let at = |i: usize, x: f64, y: f64| Placement::new_2d("s".to_string(), i, x, y, 0.0);
        let short = [at(0, 0.0, 0.0), at(1, 0.0, 10.0)];
        let long = [at(0, 0.0, 0.0), at(1, 40.0, 0.0)];
        assert!((layout_density(&g, &short, &strip) - 200.0 / (10.0 * 20.0)).abs() < 1e-12);
        assert!((layout_density(&g, &long, &strip) - 200.0 / (50.0 * 20.0)).abs() < 1e-12);
        assert_eq!(layout_density(&g, &[], &strip), 0.0);
        // On a tall strip the length is y.
        let tall = Boundary2D::rectangle(20.0, 100.0);
        let stacked = [at(0, 0.0, 0.0), at(1, 10.0, 0.0)];
        assert!((layout_density(&g, &stacked, &tall) - 1.0).abs() < 1e-12);
    }

    #[test]
    fn test_polygon_centroid_empty() {
        assert_eq!(polygon_centroid(&[]), (0.0, 0.0));
    }

    #[test]
    fn test_polygon_centroid_square() {
        let square = vec![(0.0, 0.0), (10.0, 0.0), (10.0, 10.0), (0.0, 10.0)];
        let (cx, cy) = polygon_centroid(&square);
        assert!((cx - 5.0).abs() < 1e-10);
        assert!((cy - 5.0).abs() < 1e-10);
    }

    #[test]
    fn offset_nfp_with_no_spacing_is_the_nfp() {
        let nfp = Nfp::from_polygons(vec![vec![
            (0.0, 0.0),
            (10.0, 0.0),
            (10.0, 10.0),
            (0.0, 10.0),
        ]]);
        assert_eq!(offset_nfp(&nfp, 0.0).polygons, nfp.polygons);
    }

    #[test]
    fn offset_nfp_moves_every_edge_by_the_spacing() {
        // An NFP whose corners sit on the diagonals: moving vertices away from
        // the centre shifted its edges by only 0.71 × spacing.
        let nfp = Nfp::from_polygons(vec![vec![
            (-300.0, -300.0),
            (300.0, -300.0),
            (300.0, 300.0),
            (-300.0, 300.0),
        ]]);
        let grown = offset_nfp(&nfp, 50.0);
        assert_eq!(grown.polygons.len(), 1);
        let ring = &grown.polygons[0];
        let max_x = ring.iter().map(|p| p.0).fold(f64::MIN, f64::max);
        let max_y = ring.iter().map(|p| p.1).fold(f64::MIN, f64::max);
        assert!((350.0..350.1).contains(&max_x), "right edge at {max_x}");
        assert!((350.0..350.1).contains(&max_y), "top edge at {max_y}");
    }

    #[test]
    fn inset_boundary_applies_the_margin_once_to_a_rectangle() {
        let rect = inset_boundary(&Boundary2D::rectangle(1000.0, 500.0), 50.0);
        assert_eq!(
            rect,
            vec![(50.0, 50.0), (950.0, 50.0), (950.0, 450.0), (50.0, 450.0)]
        );
    }

    #[test]
    fn inset_boundary_follows_a_slanted_edge() {
        // The bounding box of this triangle is 100 × 100; its hypotenuse is the
        // edge a box inset would ignore.
        let triangle = Boundary2D::new(vec![(0.0, 0.0), (100.0, 0.0), (0.0, 100.0)]);
        let ring = inset_boundary(&triangle, 10.0);
        assert!(ring.len() >= 3);
        for &(x, y) in &ring {
            let to_hypotenuse = (100.0 - x - y) / 2f64.sqrt();
            assert!(x >= 10.0 - 1e-6 && y >= 10.0 - 1e-6 && to_hypotenuse >= 10.0 - 1e-6);
        }
    }

    #[test]
    fn test_nesting_fitness_all_placed() {
        let f = nesting_fitness(10, 10, 0.85);
        assert!((f - 108.5).abs() < 1e-10);
    }

    #[test]
    fn test_nesting_fitness_partial() {
        let f = nesting_fitness(5, 10, 0.40);
        // 0.5 * 100 + 0.4 * 10 = 54.0
        assert!((f - 54.0).abs() < 1e-10);
    }

    #[test]
    fn test_nesting_fitness_none_placed() {
        let f = nesting_fitness(0, 10, 0.0);
        assert!((f - 0.0).abs() < 1e-10);
    }

    #[test]
    fn test_nesting_fitness_empty_total() {
        let f = nesting_fitness(0, 0, 0.0);
        assert!((f - 0.0).abs() < 1e-10);
    }
}