u_nesting_d2/

nester.rs

//! 2D nesting solver.

use crate::alns_nesting::run_alns_nesting;
use crate::boundary::Boundary2D;
use crate::brkga_nesting::run_brkga_nesting;
use crate::clamp_placement_to_boundary_with_margin;
use crate::ga_nesting::run_ga_nesting;
use crate::gdrr_nesting::run_gdrr_nesting;
use crate::geometry::Geometry2D;
use crate::nfp::{
    compute_ifp_with_margin_and_mirror, compute_nfp_mirrored, find_bottom_left_placement,
    rotate_nfp, translate_nfp, Nfp, NfpCache, PackingAxis, PlacedGeometry,
};
#[cfg(feature = "milp")]
use crate::nfp_cm_solver::run_nfp_cm_nesting;
use crate::sa_nesting::run_sa_nesting;
use crate::validate_and_filter_placements;
use u_nesting_core::alns::AlnsConfig;
use u_nesting_core::brkga::BrkgaConfig;
#[cfg(feature = "milp")]
use u_nesting_core::exact::ExactConfig;
use u_nesting_core::ga::GaConfig;
use u_nesting_core::gdrr::GdrrConfig;
use u_nesting_core::geometry::{Boundary, Geometry};
use u_nesting_core::sa::SaConfig;
use u_nesting_core::solver::{Config, ProgressCallback, ProgressInfo, Solver, Strategy};
use u_nesting_core::{Error, Placement, Result, SolveResult};

use crate::placement_utils::{hole_nfps, inset_boundary, nesting_fitness, offset_nfp};
use std::sync::atomic::{AtomicBool, Ordering};
use std::sync::Arc;
use u_nesting_core::timing::Timer;

/// The population the genetic searches (GA, BRKGA) run with when the caller
/// leaves [`Config::population_size`] unset. Each individual is decoded by
/// placing every part against the no-fit polygons of the parts placed before
/// it, so a generation costs one such layout per individual.
const SEARCH_POPULATION: usize = 30;

/// The generation budget of the genetic searches (GA, BRKGA) when the caller
/// leaves [`Config::max_generations`] unset.
const SEARCH_GENERATIONS: u32 = 50;

/// Mirror candidates to try for a geometry (`allow_flip` support): `[false,
/// true]` when mirroring is allowed, `[false]` otherwise — same shape as an
/// empty-vs-populated rotation list, so callers can nest a nested loop over
/// it exactly like `rotation_angles` without a branch at the call site.
///
/// `Geometry2D::validate()` currently rejects `allow_flip = true` outright
/// (mirroring the *strategy dispatch* is still incomplete — only BLF
/// enumerates mirror candidates so far), so in
/// every live code path today this returns `[false]`. Wired ahead of that
/// gate opening so BLF doesn't regress to the original silent-ignore bug the
/// moment it does.
fn mirror_candidates(geom: &Geometry2D) -> &'static [bool] {
    if geom.allow_flip() {
        &[false, true]
    } else {
        &[false]
    }
}

/// Returns `(placed_count, used_bounding_box_area)` for a solve result.
///
/// Used to compare two solutions on the same instance: a solution is better
/// when it places more pieces, or (tie) consumes less material *length*.
///
/// The second term is the extent along the boundary's **open (longer) axis** —
/// the material length consumed on an open-ended roll (see `used_bounding_box`
/// in `api_types`, which fixes the same convention). Comparing bounding-box
/// *area* instead is wrong for strip nesting: a tall, narrow column has a small
/// area yet a *longer* strip than a short, wide layout, so an area-based guard
/// would accept a metaheuristic result that packs the same pieces into a longer
/// roll than plain BLF — exactly the rotation-driven regression this floor
/// exists to prevent. Length is the objective consumers measure.
fn solution_quality(
    result: &SolveResult<f64>,
    geometries: &[Geometry2D],
    boundary: &Boundary2D,
) -> (usize, f64) {
    use std::collections::HashMap;
    let geom_map: HashMap<_, _> = geometries.iter().map(|g| (g.id().clone(), g)).collect();

    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 &result.placements {
        if let Some(geom) = geom_map.get(&p.geometry_id) {
            let x = p.position.first().copied().unwrap_or(0.0);
            let y = p.position.get(1).copied().unwrap_or(0.0);
            let rot = p.rotation.first().copied().unwrap_or(0.0);
            let (g_min, g_max) = geom.aabb_at_rotation(rot);
            min_x = min_x.min(x + g_min[0]);
            min_y = min_y.min(y + g_min[1]);
            max_x = max_x.max(x + g_max[0]);
            max_y = max_y.max(y + g_max[1]);
        }
    }

    if result.placements.is_empty() {
        return (0, f64::INFINITY);
    }

    // Material length = extent along the boundary's length axis.
    let strip_length = match PackingAxis::of(boundary) {
        PackingAxis::Y => max_y - min_y,
        PackingAxis::X => max_x - min_x,
    };
    (result.placements.len(), strip_length)
}

/// 2D nesting solver.
pub struct Nester2D {
    config: Config,
    cancelled: Arc<AtomicBool>,
    #[allow(dead_code)] // Will be used for caching in future optimization
    nfp_cache: NfpCache,
}

/// The spacing of the candidate grid the exact solver enumerates.
///
/// Every pair of candidates of two pieces becomes a conflict constraint, so the
/// model grows with the square of the candidate count. A fixed 1-unit grid put
/// roughly a thousand candidates on each piece -- a million constraints for two
/// pieces -- and the search returned nothing at all for anything past a single
/// piece, however long it was given. Deriving the step from the pieces keeps
/// the model at a size the search can actually finish: a quarter of the
/// smallest side any piece has, so a piece still has several distinct
/// positions along its own width.
#[cfg(feature = "milp")]
fn milp_grid_step(geometries: &[Geometry2D], boundary: &Boundary2D) -> f64 {
    let smallest_side = geometries
        .iter()
        .map(|g| {
            let (g_min, g_max) = g.aabb_at_rotation(0.0);
            (g_max[0] - g_min[0]).min(g_max[1] - g_min[1])
        })
        .fold(f64::INFINITY, f64::min);
    if !smallest_side.is_finite() || smallest_side <= 0.0 {
        return 1.0;
    }
    let (b_min, b_max) = boundary.aabb();
    let longest_side = (b_max[0] - b_min[0]).max(b_max[1] - b_min[1]);
    // Never finer than a unit, and never so fine that the sheet alone would
    // carry more than a few hundred steps in either direction.
    (smallest_side / 4.0).max(1.0).max(longest_side / 256.0)
}

/// The time budget the exact solver runs under.
///
/// `Config::time_limit_ms` of 0 means unlimited throughout this crate, but the
/// MIP search needs a number, so unlimited becomes a bound no real instance
/// reaches rather than a silent minute.
#[cfg(feature = "milp")]
fn milp_time_limit(requested_ms: u64) -> u64 {
    if requested_ms == 0 {
        u64::MAX / 2
    } else {
        requested_ms
    }
}

/// One rotation's landing spot in a Bottom-Left-Fill row, before the best of
/// them is chosen.
#[derive(Debug, Clone, Copy)]
struct BlfCandidate {
    rotation: f64,
    extent_across: f64,
    extent_along: f64,
    place_across: f64,
    place_along: f64,
    /// The candidate had to start a new row to fit.
    opened_row: bool,
}

impl Nester2D {
    /// Runs the configured strategy. A strategy 2D nesting does not provide
    /// (`ExtremePoint` is 3D only; the exact solvers need the `milp` feature)
    /// is refused rather than replaced by another.
    fn run_strategy(
        &self,
        geometries: &[Geometry2D],
        boundary: &Boundary2D,
    ) -> Result<SolveResult<f64>> {
        self.run_strategy_reporting(geometries, boundary, None)
    }

    /// [`Self::run_strategy`], reporting to `progress` where the strategy
    /// reports progress (bottom-left fill, NFP-guided, GA). The run is the
    /// same with or without a callback; the other strategies run without one.
    fn run_strategy_reporting(
        &self,
        geometries: &[Geometry2D],
        boundary: &Boundary2D,
        progress: Option<ProgressCallback>,
    ) -> Result<SolveResult<f64>> {
        match self.config.strategy {
            Strategy::BottomLeftFill => {
                self.bottom_left_fill_impl(geometries, boundary, progress.as_ref())
            }
            Strategy::NfpGuided => {
                self.nfp_guided_blf_impl(geometries, boundary, progress.as_ref())
            }
            Strategy::GeneticAlgorithm => self.genetic_algorithm(geometries, boundary, progress),
            Strategy::Brkga => self.brkga(geometries, boundary),
            Strategy::SimulatedAnnealing => self.simulated_annealing(geometries, boundary),
            Strategy::Gdrr => self.gdrr(geometries, boundary),
            Strategy::Alns => self.alns(geometries, boundary),
            #[cfg(feature = "milp")]
            Strategy::MilpExact => self.milp_exact(geometries, boundary),
            #[cfg(feature = "milp")]
            Strategy::HybridExact => self.hybrid_exact(geometries, boundary),
            Strategy::ExtremePoint => Err(Error::invalid_option(
                "strategy",
                "strategy ExtremePoint is 3D only; 2D nesting offers BottomLeftFill, \
                 NfpGuided, GeneticAlgorithm, Brkga, SimulatedAnnealing, Gdrr and Alns"
                    .into(),
            )),
            #[cfg(not(feature = "milp"))]
            Strategy::MilpExact | Strategy::HybridExact => Err(Error::invalid_option(
                "strategy",
                format!(
                    "strategy {:?} needs the `milp` feature, which this build does not have",
                    self.config.strategy
                ),
            )),
        }
    }

    /// Creates a new nester with the given configuration.
    pub fn new(config: Config) -> Self {
        Self {
            config,
            cancelled: Arc::new(AtomicBool::new(false)),
            nfp_cache: NfpCache::new(),
        }
    }

    /// Creates a nester with default configuration.
    pub fn default_config() -> Self {
        Self::new(Config::default())
    }

    /// The GA run the configuration asks for: the caller's population and
    /// generation budget, [`SEARCH_POPULATION`] / [`SEARCH_GENERATIONS`] for
    /// the ones left unset.
    fn ga_config(&self) -> GaConfig {
        let config = GaConfig::default()
            .with_population_size(self.config.population_size.unwrap_or(SEARCH_POPULATION))
            .with_max_generations(self.config.max_generations.unwrap_or(SEARCH_GENERATIONS))
            .with_crossover_rate(self.config.crossover_rate)
            .with_mutation_rate(self.config.mutation_rate)
            .with_elite_count(self.config.elite_count);
        match self.target_fitness() {
            Some(target) => config.with_target_fitness(target),
            None => config,
        }
    }

    /// The BRKGA run the configuration asks for, with the same population
    /// and generation budget as [`Self::ga_config`].
    fn brkga_config(&self) -> BrkgaConfig {
        let config = BrkgaConfig::default()
            .with_population_size(self.config.population_size.unwrap_or(SEARCH_POPULATION))
            .with_max_generations(self.config.max_generations.unwrap_or(SEARCH_GENERATIONS))
            .with_elite_fraction(0.2)
            .with_mutant_fraction(0.15)
            .with_elite_bias(0.7);
        match self.target_fitness() {
            Some(target) => config.with_target_fitness(target),
            None => config,
        }
    }

    /// The GDRR run for `time_limit_ms`, drawing from the caller's seed.
    fn gdrr_config(&self, time_limit_ms: u64) -> GdrrConfig {
        GdrrConfig {
            seed: self.config.seed,
            ..GdrrConfig::default()
                .with_max_iterations(1000) // Reduced from 5000 for faster execution
                .with_time_limit_ms(time_limit_ms)
                .with_ruin_ratio(0.1, 0.3) // Smaller ruin ratio for faster convergence
                .with_lahc_list_length(30) // Smaller list for faster convergence
        }
    }

    /// The ALNS run for `time_limit_ms`, drawing from the caller's seed.
    fn alns_config(&self, time_limit_ms: u64) -> AlnsConfig {
        AlnsConfig {
            seed: self.config.seed,
            ..AlnsConfig::default()
                .with_max_iterations(1000) // Reduced from 5000 for faster execution
                .with_time_limit_ms(time_limit_ms)
                .with_segment_size(50) // Smaller segments for faster adaptation
                .with_scores(33.0, 9.0, 13.0)
                .with_reaction_factor(0.15) // Slightly higher for faster adaptation
                .with_temperature(100.0, 0.999, 0.1) // Faster cooling
        }
    }

    /// Bottom-Left Fill algorithm implementation with rotation optimization.
    ///
    /// Intentionally does not enumerate mirror candidates (`allow_flip`):
    /// this is a pure AABB packer (row placement by bounding-box extent, no
    /// NFP/IFP shape awareness), and reflecting a polygon about an axis
    /// preserves its axis-aligned bounding-box width/height exactly — a
    /// mirrored candidate would always be bit-for-bit degenerate with its
    /// unmirrored counterpart here. See `nfp_guided_blf` for the strategy
    /// that actually benefits from mirroring.
    fn bottom_left_fill(
        &self,
        geometries: &[Geometry2D],
        boundary: &Boundary2D,
    ) -> Result<SolveResult<f64>> {
        self.bottom_left_fill_impl(geometries, boundary, None)
    }

    /// NFP-guided Bottom-Left Fill algorithm.
    ///
    /// Uses No-Fit Polygons to find optimal placement positions that minimize
    /// wasted space while ensuring no overlaps.
    fn nfp_guided_blf(
        &self,
        geometries: &[Geometry2D],
        boundary: &Boundary2D,
    ) -> Result<SolveResult<f64>> {
        self.nfp_guided_blf_impl(geometries, boundary, None)
    }

    /// The NFP-guided placement both entry points run; `callback` hears of
    /// each placed piece and of the end of the run.
    fn nfp_guided_blf_impl(
        &self,
        geometries: &[Geometry2D],
        boundary: &Boundary2D,
        callback: Option<&ProgressCallback>,
    ) -> Result<SolveResult<f64>> {
        let start = Timer::now();
        let mut result = SolveResult::new();
        let mut placements = Vec::new();
        let mut placed_geometries: Vec<PlacedGeometry> = Vec::new();
        let report = |info: ProgressInfo| {
            if let Some(callback) = callback {
                callback(info);
            }
        };
        let total_pieces: usize = geometries.iter().map(|g| g.quantity()).sum();
        let mut placed_count = 0usize;
        report(
            ProgressInfo::new()
                .with_phase("NFP Placement")
                .with_items(0, total_pieces)
                .with_elapsed(0),
        );

        let margin = self.config.margin;
        let spacing = self.config.spacing;

        // Get boundary polygon with margin applied
        let boundary_polygon = inset_boundary(boundary, margin);

        let mut total_placed_area = 0.0;

        // Sampling step for grid search (adaptive based on geometry size)
        let sample_step = self.compute_sample_step(geometries);

        for geom in geometries {
            geom.validate()?;

            // Get allowed rotation angles
            // Never empty: `validate` refuses an empty angle list.
            let rotation_angles: Vec<f64> = geom.rotations();

            let mirror_candidates = mirror_candidates(geom);

            for instance in 0..geom.quantity() {
                if self.cancelled.load(Ordering::Relaxed) {
                    result.computation_time_ms = start.elapsed_ms();
                    report(
                        ProgressInfo::new()
                            .with_phase("Cancelled")
                            .with_items(placed_count, total_pieces)
                            .with_elapsed(result.computation_time_ms)
                            .finished(),
                    );
                    return Ok(result);
                }

                // Check time limit (0 = unlimited)
                if self.config.time_limit_ms > 0 && start.elapsed_ms() >= self.config.time_limit_ms
                {
                    result.boundaries_used = if placements.is_empty() { 0 } else { 1 };
                    result.utilization = total_placed_area / boundary.measure();
                    result.computation_time_ms = start.elapsed_ms();
                    result.placements = placements;
                    report(
                        ProgressInfo::new()
                            .with_phase("Time Limit Reached")
                            .with_items(placed_count, total_pieces)
                            .with_elapsed(result.computation_time_ms)
                            .finished(),
                    );
                    return Ok(result);
                }

                // Try each (rotation, mirror) candidate to find the best placement
                let mut best_placement: Option<(f64, f64, f64, bool)> = None; // (x, y, rotation, mirror)

                for &rotation in &rotation_angles {
                    for &mirror in mirror_candidates {
                        // Compute IFP for this candidate (with margin from boundary)
                        let ifp = match compute_ifp_with_margin_and_mirror(
                            &boundary_polygon,
                            geom,
                            rotation,
                            // The boundary polygon is already inset by `margin`.
                            0.0,
                            mirror,
                        ) {
                            Ok(ifp) => ifp,
                            Err(_) => continue,
                        };

                        if ifp.is_empty() {
                            continue;
                        }

                        // Compute NFPs with all placed geometries (using cache)
                        let mut nfps: Vec<Nfp> = Vec::new();
                        for placed in &placed_geometries {
                            // Use cache for NFP computation (between original geometries at origin)
                            // Key: (placed_geometry_id, current_geometry_id, rotation, mirror_a, mirror_b)
                            let cache_key = (
                                placed.geometry.id().as_str(),
                                geom.id().as_str(),
                                rotation - placed.rotation, // Relative rotation
                                placed.mirrored,
                                mirror,
                            );

                            // Compute NFP at origin and cache it (with relative rotation)
                            // NFP is computed between the placed geometry at origin (no rotation)
                            // and the new geometry with relative rotation applied.
                            // Formula: NFP_actual = translate(rotate(NFP_relative, placed.rotation), placed.position)
                            let nfp_at_origin =
                                match self.nfp_cache.get_or_compute_mirrored(cache_key, || {
                                    let placed_at_origin = placed.geometry.clone();
                                    compute_nfp_mirrored(
                                        &placed_at_origin,
                                        geom,
                                        rotation - placed.rotation,
                                        placed.mirrored,
                                        mirror,
                                    )
                                    // Grown by `spacing` before caching: an offset commutes with the
                                    // rotation and translation below, and `config` never changes for
                                    // this nester, so the cache never holds an NFP for another spacing.
                                    .map(|nfp| offset_nfp(&nfp, spacing))
                                }) {
                                    Ok(nfp) => nfp,
                                    Err(_) => continue,
                                };

                            // Transform NFP: first rotate by placed.rotation, then translate to placed.position
                            // This correctly accounts for the placed geometry's actual orientation
                            let rotated_nfp = rotate_nfp(&nfp_at_origin, placed.rotation);
                            let translated_nfp = translate_nfp(&rotated_nfp, placed.position);
                            nfps.push(translated_nfp);
                        }

                        // `spacing` separates pieces from each other, not from the boundary —
                        // clearance to the edge is `margin`, already applied to the boundary.

                        // Find the optimal valid placement (shortest strip first)
                        nfps.extend(hole_nfps(boundary, geom, rotation, mirror, margin));
                        let nfp_refs: Vec<&Nfp> = nfps.iter().collect();
                        if let Some((x, y)) = find_bottom_left_placement(
                            &ifp,
                            &nfp_refs,
                            sample_step,
                            PackingAxis::of(boundary),
                        ) {
                            // Compare with current best in packing order: shorter strip first
                            let is_better =
                                match best_placement {
                                    None => true,
                                    Some((best_x, best_y, _, _)) => PackingAxis::of(boundary)
                                        .precedes((x, y), (best_x, best_y), 1e-6),
                                };
                            if is_better {
                                best_placement = Some((x, y, rotation, mirror));
                            }
                        }
                    }
                }

                // Place the geometry at the best position found
                if let Some((x, y, rotation, mirror)) = best_placement {
                    // Clamp to ensure geometry stays within boundary. Must be
                    // mirror-aware: a mirrored candidate's true local extents
                    // are reflected, not just the unmirrored AABB — using the
                    // wrong one here shifts a already NFP-validated (x, y)
                    // into an overlap (found via `placement_has_no_overlap_with_mirroring`,
                    // `fuzz_robustness.rs`).
                    let geom_aabb = geom.aabb_at_rotation_mirrored(rotation, mirror);
                    let boundary_aabb = boundary.aabb();

                    if let Some((clamped_x, clamped_y)) = clamp_placement_to_boundary_with_margin(
                        x,
                        y,
                        geom_aabb,
                        boundary_aabb,
                        margin,
                    ) {
                        let placement = Placement::new_2d(
                            geom.id().clone(),
                            instance,
                            clamped_x,
                            clamped_y,
                            rotation,
                        )
                        .with_mirrored(mirror);

                        placements.push(placement);
                        placed_geometries.push(
                            PlacedGeometry::new(geom.clone(), (clamped_x, clamped_y), rotation)
                                .with_mirrored(mirror),
                        );
                        total_placed_area += geom.measure();
                        placed_count += 1;
                        report(
                            ProgressInfo::new()
                                .with_phase("NFP Placement")
                                .with_items(placed_count, total_pieces)
                                .with_utilization(total_placed_area / boundary.measure())
                                .with_elapsed(start.elapsed_ms()),
                        );
                    } else {
                        // Could not place - geometry doesn't fit
                        result.unplaced.push(geom.id().clone());
                    }
                } else {
                    // Could not place this instance
                    result.unplaced.push(geom.id().clone());
                }
            }
        }

        result.placements = placements;
        result.boundaries_used = 1;
        result.utilization = total_placed_area / boundary.measure();
        result.computation_time_ms = start.elapsed_ms();
        report(
            ProgressInfo::new()
                .with_phase("Complete")
                .with_items(placed_count, total_pieces)
                .with_utilization(result.utilization)
                .with_elapsed(result.computation_time_ms)
                .finished(),
        );

        Ok(result)
    }

    /// Computes an adaptive sample step based on geometry sizes.
    fn compute_sample_step(&self, geometries: &[Geometry2D]) -> f64 {
        if geometries.is_empty() {
            return 1.0;
        }

        // Use the smallest geometry dimension divided by 4 as sample step
        let mut min_dim = f64::INFINITY;
        for geom in geometries {
            let (g_min, g_max) = geom.aabb();
            let width = g_max[0] - g_min[0];
            let height = g_max[1] - g_min[1];
            min_dim = min_dim.min(width).min(height);
        }

        // Clamp sample step to reasonable range
        (min_dim / 4.0).clamp(0.5, 10.0)
    }

    /// Returns whichever of `meta` (a metaheuristic result) and a fresh
    /// Bottom-Left-Fill solve packs better.
    ///
    /// The search strategies (GA/BRKGA/SA/GDRR/ALNS) can converge to — or, on a
    /// short time limit, stop at — a solution worse
    /// than the deterministic greedy baseline. Guarding against this guarantees
    /// they never return a solution inferior to BLF: a metaheuristic that fails
    /// to beat the greedy floor simply returns the greedy solution. BLF is
    /// effectively free relative to a metaheuristic run.
    fn not_worse_than_baselines(
        &self,
        meta: SolveResult<f64>,
        geometries: &[Geometry2D],
        boundary: &Boundary2D,
        greedy: Option<SolveResult<f64>>,
    ) -> SolveResult<f64> {
        let baselines = self
            .bottom_left_fill(geometries, boundary)
            .ok()
            .into_iter()
            .chain(greedy);
        let quality = |r: &SolveResult<f64>| solution_quality(r, geometries, boundary);
        // A baseline wins if it places strictly more pieces, or ties on count
        // while consuming a strictly shorter strip length.
        meta.floored_at(baselines, |a, b| {
            let (a, b) = (quality(a), quality(b));
            a.0 > b.0 || (a.0 == b.0 && a.1 < b.1 - 1e-6)
        })
    }

    /// The time a search strategy may spend on one strip: a quarter of the
    /// limit (assuming up to about four strips), at least 5 s, never more than
    /// the limit itself; `default_ms` when there is no limit.
    /// The fitness a search stops at for `target_utilization`: every part
    /// placed, at that utilization ([`nesting_fitness`] of a complete layout).
    fn target_fitness(&self) -> Option<f64> {
        self.config
            .target_utilization
            .map(|target| nesting_fitness(1, 1, target))
    }

    fn search_budget_ms(&self, default_ms: u64) -> u64 {
        if self.config.time_limit_ms > 0 {
            (self.config.time_limit_ms / 4)
                .max(5000)
                .min(self.config.time_limit_ms)
        } else {
            default_ms
        }
    }

    /// Runs greedy no-fit-polygon placement within `budget_ms` as the starting
    /// point of a search, and returns it with the time left for the search.
    ///
    /// Searching from a random start can end, on a short limit or a large
    /// rotation set, with a layout the one greedy pass would have beaten. The
    /// search is seeded with this layout where it can be, and floored at it.
    fn greedy_start(
        &self,
        geometries: &[Geometry2D],
        boundary: &Boundary2D,
        budget_ms: u64,
    ) -> (Option<SolveResult<f64>>, u64) {
        let started = Timer::now();
        let greedy = Nester2D {
            config: Config {
                time_limit_ms: budget_ms,
                ..self.config.clone()
            },
            cancelled: self.cancelled.clone(),
            nfp_cache: NfpCache::new(),
        }
        .nfp_guided_blf(geometries, boundary)
        .ok();
        (
            greedy,
            budget_ms.saturating_sub(started.elapsed_ms()).max(1),
        )
    }

    /// Genetic Algorithm based nesting optimization.
    ///
    /// Uses GA to optimize placement order and rotations, with NFP-guided
    /// decoding for collision-free placements. `progress` hears of each
    /// generation; the run -- its time budget, greedy start and floor -- is
    /// the same with or without it.
    fn genetic_algorithm(
        &self,
        geometries: &[Geometry2D],
        boundary: &Boundary2D,
        progress: Option<ProgressCallback>,
    ) -> Result<SolveResult<f64>> {
        let (time_limit_ms, greedy) = {
            let budget = self.search_budget_ms(15000);
            let (greedy, remaining) = self.greedy_start(geometries, boundary, budget);
            (remaining, greedy)
        };

        let ga_config = self
            .ga_config()
            .with_time_limit(std::time::Duration::from_millis(time_limit_ms));

        let result = run_ga_nesting(
            geometries,
            boundary,
            &self.config,
            ga_config,
            self.cancelled.clone(),
            greedy.as_ref().map(|g| g.placements.as_slice()),
            progress,
        )?;

        Ok(self.not_worse_than_baselines(result, geometries, boundary, greedy))
    }

    /// BRKGA (Biased Random-Key Genetic Algorithm) based nesting optimization.
    ///
    /// Uses random-key encoding and biased crossover for robust optimization.
    fn brkga(&self, geometries: &[Geometry2D], boundary: &Boundary2D) -> Result<SolveResult<f64>> {
        let (time_limit_ms, greedy) = {
            let budget = self.search_budget_ms(15000);
            let (greedy, remaining) = self.greedy_start(geometries, boundary, budget);
            (remaining, greedy)
        };

        let brkga_config = self
            .brkga_config()
            .with_time_limit(std::time::Duration::from_millis(time_limit_ms));

        let result = run_brkga_nesting(
            geometries,
            boundary,
            &self.config,
            brkga_config,
            self.cancelled.clone(),
            greedy.as_ref().map(|g| g.placements.as_slice()),
        )?;

        Ok(self.not_worse_than_baselines(result, geometries, boundary, greedy))
    }

    /// Simulated Annealing based nesting optimization.
    ///
    /// Uses neighborhood operators to explore solution space with temperature-based
    /// acceptance probability.
    fn simulated_annealing(
        &self,
        geometries: &[Geometry2D],
        boundary: &Boundary2D,
    ) -> Result<SolveResult<f64>> {
        let (time_limit_ms, greedy) = {
            let budget = self.search_budget_ms(10000);
            let (greedy, remaining) = self.greedy_start(geometries, boundary, budget);
            (remaining, greedy)
        };

        let sa_config = SaConfig::default()
            .with_initial_temp(50.0) // Lower initial temp for faster convergence
            .with_final_temp(1.0) // Higher final temp to finish faster
            .with_cooling_rate(0.9) // Faster cooling (was 0.95)
            .with_iterations_per_temp(20) // Fewer iterations per temp (was 50)
            .with_max_iterations(500) // Much fewer max iterations (was 10000)
            .with_time_limit(std::time::Duration::from_millis(time_limit_ms));
        let sa_config = match self.target_fitness() {
            Some(target) => sa_config.with_target_fitness(target),
            None => sa_config,
        };

        let result = run_sa_nesting(
            geometries,
            boundary,
            &self.config,
            sa_config,
            self.cancelled.clone(),
            greedy.as_ref().map(|g| g.placements.as_slice()),
        )?;

        Ok(self.not_worse_than_baselines(result, geometries, boundary, greedy))
    }

    /// Goal-Driven Ruin and Recreate (GDRR) optimization.
    fn gdrr(&self, geometries: &[Geometry2D], boundary: &Boundary2D) -> Result<SolveResult<f64>> {
        let (time_limit, greedy) = {
            let budget = self.search_budget_ms(10000);
            let (greedy, remaining) = self.greedy_start(geometries, boundary, budget);
            (remaining, greedy)
        };
        let result = run_gdrr_nesting(
            geometries,
            boundary,
            &self.config,
            &self.gdrr_config(time_limit),
            self.cancelled.clone(),
        )?;

        Ok(self.not_worse_than_baselines(result, geometries, boundary, greedy))
    }

    /// Adaptive Large Neighborhood Search (ALNS) optimization.
    fn alns(&self, geometries: &[Geometry2D], boundary: &Boundary2D) -> Result<SolveResult<f64>> {
        self.alns_within(geometries, boundary, self.search_budget_ms(10000))
    }

    /// [`Self::alns`] in `budget_ms`, for a caller that has spent part of the
    /// limit already.
    fn alns_within(
        &self,
        geometries: &[Geometry2D],
        boundary: &Boundary2D,
        budget_ms: u64,
    ) -> Result<SolveResult<f64>> {
        let (greedy, time_limit) = self.greedy_start(geometries, boundary, budget_ms);
        let result = run_alns_nesting(
            geometries,
            boundary,
            &self.config,
            &self.alns_config(time_limit),
            self.cancelled.clone(),
        )?;

        Ok(self.not_worse_than_baselines(result, geometries, boundary, greedy))
    }

    /// MILP-based exact solver.
    ///
    /// Uses the NFP Covering Model formulation (`nfp_cm_solver`) — the other
    /// MILP module this crate carried (`milp_solver`, a continuous-position
    /// Big-M formulation) placed nothing for *any* input through this public
    /// entry point, including plain rectangles with no rotation or mirroring
    /// involved, and was covered by no test that actually went through
    /// `Strategy::MilpExact` end-to-end — every existing MILP test called
    /// straight into a solver module's own function. `nfp_cm_solver` is the
    /// one with working, tested coverage (including `allow_flip` support),
    /// so this now calls that instead; `milp_solver` was removed as
    /// orphaned (nothing else referenced it).
    #[cfg(feature = "milp")]
    fn milp_exact(
        &self,
        geometries: &[Geometry2D],
        boundary: &Boundary2D,
    ) -> Result<SolveResult<f64>> {
        // The caller's limit is the limit. It used to be raised to a minute,
        // so a request for a second ran for minutes -- the same shape as the
        // time contract the other strategies keep. 0 means unlimited here as
        // it does everywhere else in this crate.
        let exact_config = ExactConfig::default()
            .with_time_limit_ms(milp_time_limit(self.config.time_limit_ms))
            .with_max_items(15)
            .with_rotation_steps(4)
            .with_grid_step(milp_grid_step(geometries, boundary));

        run_nfp_cm_nesting(
            geometries,
            boundary,
            &self.config,
            &exact_config,
            self.cancelled.clone(),
        )
    }

    /// Hybrid exact solver: try MILP first, fallback to heuristic.
    #[cfg(feature = "milp")]
    fn hybrid_exact(
        &self,
        geometries: &[Geometry2D],
        boundary: &Boundary2D,
    ) -> Result<SolveResult<f64>> {
        let started = Timer::now();
        // Count total instances
        let total_instances: usize = geometries.iter().map(|g| g.quantity()).sum();

        // If small enough, try exact
        if total_instances <= 15 {
            // Half the caller's budget for the exact attempt, leaving the
            // rest for the heuristic this falls back to. No floor: see
            // `milp_time_limit`.
            let exact_config = ExactConfig::default()
                .with_time_limit_ms(milp_time_limit(self.config.time_limit_ms / 2))
                .with_max_items(15);

            let exact_result = run_nfp_cm_nesting(
                geometries,
                boundary,
                &self.config,
                &exact_config,
                self.cancelled.clone(),
            )?;

            // If got a good solution, return it
            if !exact_result.placements.is_empty() {
                return Ok(exact_result);
            }
        }

        // Fallback to ALNS (best heuristic), in what is left of the caller's
        // limit. It used to get its usual budget afresh after the exact attempt
        // had spent half the limit, so any limit under ten seconds ran over by
        // up to half again.
        let budget = match self.config.time_limit_ms {
            0 => self.search_budget_ms(10000),
            limit => self
                .search_budget_ms(10000)
                .min(limit.saturating_sub(started.elapsed_ms()).max(1)),
        };
        self.alns_within(geometries, boundary, budget)
    }

    /// The bottom-left fill both entry points run.
    ///
    /// Rows are filled across the boundary's width and advance along its length
    /// ([`PackingAxis::of`]), so a wide sheet is filled in columns and a tall one
    /// in rows: the layout uses as little of the length as it can. The row logic
    /// works in `(across, along)` coordinates and maps back to `(x, y)` only when
    /// a placement is recorded.
    fn bottom_left_fill_impl(
        &self,
        geometries: &[Geometry2D],
        boundary: &Boundary2D,
        callback: Option<&ProgressCallback>,
    ) -> Result<SolveResult<f64>> {
        let start = Timer::now();
        let mut result = SolveResult::new();
        let mut placements = Vec::new();
        let report = |info: ProgressInfo| {
            if let Some(callback) = callback {
                callback(info);
            }
        };

        let axis = PackingAxis::of(boundary);
        // (across, along) <-> (x, y)
        let to_ca = |x: f64, y: f64| match axis {
            PackingAxis::Y => (x, y),
            PackingAxis::X => (y, x),
        };
        let from_ca = |across: f64, along: f64| match axis {
            PackingAxis::Y => (across, along),
            PackingAxis::X => (along, across),
        };

        let (b_min, b_max) = boundary.aabb();
        let margin = self.config.margin;
        let spacing = self.config.spacing;
        let (min_across, min_along) = to_ca(b_min[0] + margin, b_min[1] + margin);
        let (max_across, max_along) = to_ca(b_max[0] - margin, b_max[1] - margin);

        // A hole is taken as its own AABB grown by `margin`, in (across, along).
        // This packer already reads the boundary as an AABB, so reading holes
        // the same way keeps it one kind of packer; it is conservative for a
        // non-rectangular hole, which blocks a little more than it occupies.
        // The alternative is what this replaces: placing over a hole and having
        // the placement filter drop the piece, which loses the position.
        let hole_boxes: Vec<(f64, f64, f64, f64)> = boundary
            .holes()
            .iter()
            .filter(|hole| !hole.is_empty())
            .map(|hole| {
                let (mut x0, mut y0) = (f64::INFINITY, f64::INFINITY);
                let (mut x1, mut y1) = (f64::NEG_INFINITY, f64::NEG_INFINITY);
                for &(x, y) in hole {
                    x0 = x0.min(x);
                    y0 = y0.min(y);
                    x1 = x1.max(x);
                    y1 = y1.max(y);
                }
                let (a0, l0) = to_ca(x0 - margin, y0 - margin);
                let (a1, l1) = to_ca(x1 + margin, y1 + margin);
                (a0.min(a1), l0.min(l1), a0.max(a1), l0.max(l1))
            })
            .collect();

        // The far edge of the hole a footprint runs into, if it runs into one.
        // Overlap is measured on open intervals: touching a hole's margin ring
        // exactly is clear, which is what `margin` means everywhere else.
        let blocking_hole_edge = |across: f64, along: f64, ea: f64, el: f64| -> Option<f64> {
            hole_boxes
                .iter()
                .filter(|(a0, l0, a1, l1)| {
                    across < *a1 && across + ea > *a0 && along < *l1 && along + el > *l0
                })
                .map(|(_, _, a1, _)| *a1)
                .fold(None, |far: Option<f64>, a1| {
                    Some(far.map_or(a1, |f| f.max(a1)))
                })
        };

        let mut cursor_across = min_across;
        let mut cursor_along = min_along;
        let mut row_depth = 0.0_f64;
        let mut total_placed_area = 0.0;
        let total_pieces: usize = geometries.iter().map(|g| g.quantity()).sum();
        let mut placed_count = 0usize;

        report(
            ProgressInfo::new()
                .with_phase("BLF Placement")
                .with_items(0, total_pieces)
                .with_elapsed(0),
        );

        for geom in geometries {
            geom.validate()?;

            // Never empty: `validate` refuses an empty angle list.
            let rotation_angles: Vec<f64> = geom.rotations();

            for instance in 0..geom.quantity() {
                if self.cancelled.load(Ordering::Relaxed) {
                    result.computation_time_ms = start.elapsed_ms();
                    report(
                        ProgressInfo::new()
                            .with_phase("Cancelled")
                            .with_items(placed_count, total_pieces)
                            .with_elapsed(result.computation_time_ms)
                            .finished(),
                    );
                    return Ok(result);
                }

                // Check time limit (0 = unlimited)
                if self.config.time_limit_ms > 0 && start.elapsed_ms() >= self.config.time_limit_ms
                {
                    result.boundaries_used = if placements.is_empty() { 0 } else { 1 };
                    result.utilization = total_placed_area / boundary.measure();
                    result.computation_time_ms = start.elapsed_ms();
                    result.placements = placements;
                    report(
                        ProgressInfo::new()
                            .with_phase("Time Limit Reached")
                            .with_items(placed_count, total_pieces)
                            .with_elapsed(result.computation_time_ms)
                            .finished(),
                    );
                    return Ok(result);
                }

                // The rotation that advances least: along the row if it still fits
                // there, otherwise along the length when it opens a new row.
                let mut best_fit: Option<BlfCandidate> = None;

                for &rotation in &rotation_angles {
                    let (g_min, g_max) = geom.aabb_at_rotation(rotation);
                    let (extent_across, extent_along) =
                        to_ca(g_max[0] - g_min[0], g_max[1] - g_min[1]);

                    if extent_across > max_across - min_across
                        || extent_along > max_along - min_along
                    {
                        continue;
                    }

                    let mut place_across = cursor_across;
                    let mut place_along = cursor_along;
                    // Depth of the row the candidate currently sits in. Once it
                    // opens a row, that row holds only this piece so far.
                    let mut depth = row_depth;
                    let mut opened_row = false;
                    // Each turn either opens a row (advancing `along` by at
                    // least the piece's own length once a row has been opened,
                    // so `max_along` ends it) or steps past one hole's far edge.
                    // The guard covers degenerate geometry rather than the
                    // normal case.
                    let mut turns = 0usize;
                    let max_turns = 4 * (hole_boxes.len() + 1) * (hole_boxes.len() + 2) + 16;
                    let fits = loop {
                        if place_across + extent_across > max_across {
                            place_across = min_across;
                            place_along += depth + spacing;
                            depth = extent_along;
                            opened_row = true;
                        }
                        if place_along + extent_along > max_along {
                            break false;
                        }
                        match blocking_hole_edge(
                            place_across,
                            place_along,
                            extent_across,
                            extent_along,
                        ) {
                            None => break true,
                            Some(past) => place_across = past + spacing,
                        }
                        turns += 1;
                        if turns > max_turns {
                            break false;
                        }
                    };
                    if !fits {
                        continue;
                    }

                    // Advancement a placement costs: along the length when it
                    // opens a row, across the row otherwise. A candidate is
                    // measured against the best so far with its own extents, so
                    // within one row the first rotation that fits is kept.
                    let cost = |c: &BlfCandidate| {
                        if c.opened_row {
                            c.place_along - min_along + c.extent_along
                        } else {
                            c.place_across - min_across + c.extent_across
                        }
                    };
                    let candidate = BlfCandidate {
                        rotation,
                        extent_across,
                        extent_along,
                        place_across,
                        place_along,
                        opened_row,
                    };
                    let better = best_fit
                        .as_ref()
                        .is_none_or(|best| cost(&candidate) < cost(best) - 1e-6);
                    if better {
                        best_fit = Some(candidate);
                    }
                }

                let Some(BlfCandidate {
                    rotation,
                    extent_across,
                    extent_along,
                    place_across,
                    place_along,
                    opened_row,
                }) = best_fit
                else {
                    result.unplaced.push(geom.id().clone());
                    continue;
                };
                if opened_row {
                    row_depth = 0.0;
                }

                let geom_aabb = geom.aabb_at_rotation(rotation);
                let (place_x, place_y) = from_ca(place_across, place_along);
                let Some((x, y)) = clamp_placement_to_boundary_with_margin(
                    place_x - geom_aabb.0[0],
                    place_y - geom_aabb.0[1],
                    geom_aabb,
                    (b_min, b_max),
                    margin,
                ) else {
                    result.unplaced.push(geom.id().clone());
                    continue;
                };

                placements.push(Placement::new_2d(
                    geom.id().clone(),
                    instance,
                    x,
                    y,
                    rotation,
                ));
                total_placed_area += geom.measure();
                placed_count += 1;

                // Continue from where the clamped piece actually sits.
                let (actual_across, actual_along) = to_ca(x + geom_aabb.0[0], y + geom_aabb.0[1]);
                cursor_across = actual_across + extent_across + spacing;
                cursor_along = actual_along;
                row_depth = row_depth.max(extent_along);

                report(
                    ProgressInfo::new()
                        .with_phase("BLF Placement")
                        .with_items(placed_count, total_pieces)
                        .with_utilization(total_placed_area / boundary.measure())
                        .with_elapsed(start.elapsed_ms()),
                );
            }
        }

        result.placements = placements;
        result.boundaries_used = 1;
        result.utilization = total_placed_area / boundary.measure();
        result.computation_time_ms = start.elapsed_ms();

        report(
            ProgressInfo::new()
                .with_phase("Complete")
                .with_items(placed_count, total_pieces)
                .with_utilization(result.utilization)
                .with_elapsed(result.computation_time_ms)
                .finished(),
        );

        Ok(result)
    }

    /// Solves nesting with automatic multi-strip support.
    ///
    /// When items don't fit in a single strip, automatically creates additional strips.
    /// Each placement's `boundary_index` indicates which strip it belongs to.
    /// Validates every input geometry once, at the solve entry point.
    ///
    /// Hoisting validation here (rather than relying on per-strategy calls)
    /// guarantees no dispatch path — including the progress/callback path and
    /// every metaheuristic — can bypass input rejection of degenerate,
    /// self-intersecting, or `allow_flip` geometries.
    ///
    /// An empty list is refused (`empty_input`): there is nothing to place,
    /// and a layout of no parts on one sheet answers a question nobody asked.
    fn validate_geometries(&self, geometries: &[Geometry2D]) -> Result<()> {
        use u_nesting_core::geometry::Geometry;
        if geometries.is_empty() {
            return Err(Error::empty_input("geometries"));
        }
        self.config.validate()?;
        for geom in geometries {
            geom.validate()?;
        }
        u_nesting_core::geometry::ensure_unique_ids(geometries)
    }

    /// Positions are adjusted so that strip N items have x offset of N * strip_width.
    pub fn solve_multi_strip(
        &self,
        geometries: &[Geometry2D],
        boundary: &Boundary2D,
    ) -> Result<SolveResult<f64>> {
        boundary.validate()?;
        self.validate_geometries(geometries)?;
        self.cancelled.store(false, Ordering::Relaxed);

        let (b_min, b_max) = boundary.aabb();
        let strip_width = b_max[0] - b_min[0];

        let mut final_result = SolveResult::new();
        let mut remaining_geometries: Vec<Geometry2D> = geometries.to_vec();
        let mut strip_index = 0;
        let max_strips = 100; // Safety limit

        // Global per-geometry placed counter so instance indices stay unique across
        // strips (each strip re-numbers its own placements from 0).
        let mut placed_total: std::collections::HashMap<String, usize> =
            std::collections::HashMap::new();

        while !remaining_geometries.is_empty() && strip_index < max_strips {
            if self.cancelled.load(Ordering::Relaxed) {
                break;
            }

            // Solve on current strip
            let strip_result = self.run_strategy(&remaining_geometries, boundary)?;

            // Validate and filter out-of-bounds placements for this strip
            let strip_result = validate_and_filter_placements(
                strip_result,
                &remaining_geometries,
                boundary,
                self.config.margin,
            );

            if strip_result.placements.is_empty() {
                // No progress: every remaining geometry is individually too large for
                // an empty strip. Stop; the after-loop sweep records them as unplaced.
                break;
            }

            // Count instances of each geometry placed on this strip (instance-level,
            // not id-level) to drive remaining-quantity reduction and unique numbering.
            let mut strip_placed: std::collections::HashMap<String, usize> =
                std::collections::HashMap::new();

            // Adjust placements for this strip and add to final result.
            for mut placement in strip_result.placements {
                let gid = placement.geometry_id.clone();
                // Globally-unique instance index = already-placed of this id + running
                // count within this strip (each strip re-numbers its own from 0).
                let prior = placed_total.get(&gid).copied().unwrap_or(0);
                let in_strip = strip_placed.get(&gid).copied().unwrap_or(0);
                placement.instance = prior + in_strip;
                // Offset x position by strip_index * strip_width (global strip frame).
                if !placement.position.is_empty() {
                    placement.position[0] += strip_index as f64 * strip_width;
                }
                placement.boundary_index = strip_index;
                *strip_placed.entry(gid).or_insert(0) += 1;
                final_result.placements.push(placement);
            }

            // Reduce each geometry's remaining quantity by the instances placed this
            // strip. Fully-placed geometries drop out; partially-placed ones carry the
            // remainder to the next strip (fixes the prior id-level silent loss).
            for (gid, cnt) in &strip_placed {
                *placed_total.entry(gid.clone()).or_insert(0) += cnt;
            }
            remaining_geometries = remaining_geometries
                .into_iter()
                .filter_map(|g| {
                    let placed_here = strip_placed.get(g.id()).copied().unwrap_or(0);
                    let new_quantity = g.quantity().saturating_sub(placed_here);
                    if new_quantity == 0 {
                        None
                    } else {
                        Some(g.with_quantity(new_quantity))
                    }
                })
                .collect();

            strip_index += 1;
        }

        // Any geometry still remaining (too large to place, or hit the strip cap) is
        // unplaced at the instance level — record an id entry each (deduplicated below).
        for g in &remaining_geometries {
            final_result.unplaced.push(g.id().clone());
        }

        final_result.boundaries_used = strip_index;
        final_result.deduplicate_unplaced();
        // Authoritative instance-level request total (Σ quantity), mirroring solve().
        final_result.total_requested = geometries.iter().map(|g| g.quantity()).sum();

        // Calculate per-strip statistics for accurate utilization
        let (b_min, b_max) = boundary.aabb();
        let strip_height = b_max[1] - b_min[1]; // Height of each strip

        // Group placements by strip and calculate stats
        let mut strip_stats_map: std::collections::HashMap<usize, (f64, f64, usize)> =
            std::collections::HashMap::new(); // strip_index -> (max_x, piece_area, count)

        for placement in &final_result.placements {
            let strip_idx = placement.boundary_index;
            // Get the geometry to calculate its area and right edge
            if let Some(geom) = geometries.iter().find(|g| g.id() == &placement.geometry_id) {
                use u_nesting_core::geometry::Geometry;
                let piece_area = geom.measure();
                let rotation = placement.rotation.first().copied().unwrap_or(0.0);
                let (_g_min, g_max) = geom.aabb_at_rotation(rotation);
                // Position is where geometry's origin is placed
                // The actual right edge is position.x + g_max[0] (relative to origin)
                let local_x = placement.position[0] - (strip_idx as f64 * strip_width);
                let right_edge = local_x + g_max[0];

                let entry = strip_stats_map.entry(strip_idx).or_insert((0.0, 0.0, 0));
                entry.0 = entry.0.max(right_edge); // max_x (used_length)
                entry.1 += piece_area; // total piece area
                entry.2 += 1; // piece count
            }
        }

        // Convert to StripStats vec
        use u_nesting_core::result::StripStats;
        let mut strip_stats: Vec<StripStats> = strip_stats_map
            .into_iter()
            .map(|(idx, (used_length, piece_area, count))| StripStats {
                strip_index: idx,
                used_length,
                piece_area,
                piece_count: count,
                strip_width,  // Width of boundary (X dimension)
                strip_height, // Height of boundary (Y dimension, fixed)
            })
            .collect();
        strip_stats.sort_by_key(|s| s.strip_index);

        // Calculate accurate utilization
        // Material used = strip_height (fixed dimension) × used_length (consumed length)
        let total_piece_area: f64 = strip_stats.iter().map(|s| s.piece_area).sum();
        let total_material_used: f64 = strip_stats
            .iter()
            .map(|s| s.strip_height * s.used_length)
            .sum();

        final_result.strip_stats = strip_stats;
        final_result.total_piece_area = total_piece_area;
        final_result.total_material_used = total_material_used;

        if total_material_used > 0.0 {
            final_result.utilization = total_piece_area / total_material_used;
        }

        Ok(final_result)
    }
}

impl Nester2D {
    /// The solve both [`Solver`] entry points run: the same checks, strategy
    /// run and placement filtering, with `progress` told of the run where the
    /// strategy reports it.
    fn solve_reporting(
        &self,
        geometries: &[Geometry2D],
        boundary: &Boundary2D,
        progress: Option<ProgressCallback>,
    ) -> Result<SolveResult<f64>> {
        boundary.validate()?;
        self.validate_geometries(geometries)?;

        // Reset cancellation flag
        self.cancelled.store(false, Ordering::Relaxed);

        let initial_result = self.run_strategy_reporting(geometries, boundary, progress)?;

        // Validate all placements and remove any that are outside the boundary
        let mut result = validate_and_filter_placements(
            initial_result,
            geometries,
            boundary,
            self.config.margin,
        );

        // Remove duplicate entries from unplaced list
        result.deduplicate_unplaced();
        // Authoritative instance-level request total (Σ quantity), recorded once
        // at the top-level entry point where `geometries` is the full request.
        result.total_requested = geometries.iter().map(|g| g.quantity()).sum();
        Ok(result)
    }
}

impl Solver for Nester2D {
    type Geometry = Geometry2D;
    type Boundary = Boundary2D;
    type Scalar = f64;

    fn solve(
        &self,
        geometries: &[Self::Geometry],
        boundary: &Self::Boundary,
    ) -> Result<SolveResult<f64>> {
        self.solve_reporting(geometries, boundary, None)
    }

    fn solve_with_progress(
        &self,
        geometries: &[Self::Geometry],
        boundary: &Self::Boundary,
        callback: ProgressCallback,
    ) -> Result<SolveResult<f64>> {
        self.solve_reporting(geometries, boundary, Some(callback))
    }

    fn cancel(&self) {
        self.cancelled.store(true, Ordering::Relaxed);
    }
}

#[cfg(test)]
mod tests {
    #[cfg(feature = "milp")]
    mod milp_budget {
        use super::super::{milp_grid_step, milp_time_limit};
        use crate::{Boundary2D, Geometry2D};

        /// 0 means unlimited everywhere else in this crate, and the exact
        /// solver has to spell that as a number rather than as a minute.
        #[test]
        fn an_unlimited_request_does_not_become_a_minute() {
            assert_eq!(milp_time_limit(0), u64::MAX / 2);
        }

        /// The floor this replaces turned a one-second request into a
        /// one-minute run.
        #[test]
        fn a_short_request_is_kept_short() {
            assert_eq!(milp_time_limit(1_000), 1_000);
            assert_eq!(milp_time_limit(250), 250);
        }

        /// Candidates are enumerated on this grid and every pair of them
        /// becomes a constraint, so the step has to come from the pieces --
        /// a fixed unit grid is what made the model unsolvable.
        #[test]
        fn the_grid_step_follows_the_smallest_piece() {
            let boundary = Boundary2D::rectangle(400.0, 100.0);
            let pieces = [
                Geometry2D::rectangle("wide", 80.0, 40.0),
                Geometry2D::rectangle("small", 20.0, 60.0),
            ];
            // The smallest side across both pieces is 20, so a quarter of it.
            assert_eq!(milp_grid_step(&pieces, &boundary), 5.0);
        }

        /// A sheet far larger than its pieces would still carry an unbounded
        /// number of steps, so the sheet sets a floor of its own.
        #[test]
        fn a_huge_sheet_keeps_the_step_from_getting_too_fine() {
            let boundary = Boundary2D::rectangle(100_000.0, 100.0);
            let pieces = [Geometry2D::rectangle("small", 4.0, 4.0)];
            assert_eq!(milp_grid_step(&pieces, &boundary), 100_000.0 / 256.0);
        }

        /// Never finer than a unit, whatever the pieces are.
        #[test]
        fn the_step_never_goes_below_a_unit() {
            let boundary = Boundary2D::rectangle(10.0, 10.0);
            let pieces = [Geometry2D::rectangle("tiny", 0.5, 0.5)];
            assert_eq!(milp_grid_step(&pieces, &boundary), 1.0);
        }
    }

    use super::*;
    use crate::placement_utils::polygon_centroid;

    /// The result names placements by geometry id (and instance), and the
    /// placement check looks each geometry up by id -- two geometries with one
    /// id would be told apart by neither.
    #[test]
    fn a_geometry_id_given_twice_is_refused_by_every_entry_point() {
        let geometries = [
            Geometry2D::rectangle("P", 10.0, 10.0),
            Geometry2D::rectangle("Q", 10.0, 10.0),
            Geometry2D::rectangle("P", 30.0, 5.0),
        ];
        let boundary = Boundary2D::rectangle(100.0, 100.0);
        let nester = Nester2D::default_config();
        let errors = [
            nester.solve(&geometries, &boundary).expect_err("solve"),
            nester
                .solve_multi_strip(&geometries, &boundary)
                .expect_err("solve_multi_strip"),
            nester
                .solve_with_progress(&geometries, &boundary, Box::new(|_| {}))
                .expect_err("solve_with_progress"),
        ];
        for err in errors {
            let msg = err.to_string();
            assert!(msg.contains("'P'"), "{msg}");
            assert!(msg.contains("positions 0 and 2"), "{msg}");
        }
    }

    /// Nothing to place is a request to refuse, not an empty layout on one
    /// sheet.
    #[test]
    fn an_empty_list_of_geometries_is_refused_by_every_entry_point() {
        let boundary = Boundary2D::rectangle(100.0, 100.0);
        let nester = Nester2D::default_config();
        let errors = [
            nester.solve(&[], &boundary).expect_err("solve"),
            nester
                .solve_multi_strip(&[], &boundary)
                .expect_err("solve_multi_strip"),
            nester
                .solve_with_progress(&[], &boundary, Box::new(|_| {}))
                .expect_err("solve_with_progress"),
        ];
        for err in errors {
            assert_eq!(
                err,
                Error::EmptyInput {
                    parameter: "geometries".into()
                }
            );
        }
    }

    #[test]
    fn the_seed_reaches_the_searches_that_build_their_own_run() {
        // GA, BRKGA and SA read `config.seed` from the configuration they are
        // given; GDRR and ALNS run on a configuration of their own, which
        // used to leave the seed out.
        let seeded = Nester2D::new(Config::default().with_seed(7));
        assert_eq!(seeded.gdrr_config(1000).seed, Some(7));
        assert_eq!(seeded.alns_config(1000).seed, Some(7));
        let unseeded = Nester2D::default_config();
        assert_eq!(unseeded.gdrr_config(1000).seed, None);
        assert_eq!(unseeded.alns_config(1000).seed, None);
    }

    /// A population or generation budget the caller sets is the one the
    /// search runs with -- they used to be cut to 30 and 50, and BRKGA ran 50
    /// generations whatever was asked; left unset, the search runs with its
    /// own.
    #[test]
    fn a_set_budget_is_run_and_an_unset_one_is_the_strategy_s() {
        assert_eq!((SEARCH_POPULATION, SEARCH_GENERATIONS), (30, 50));
        let unset = Nester2D::default_config();
        assert_eq!(unset.ga_config().population_size, SEARCH_POPULATION);
        assert_eq!(unset.ga_config().max_generations, SEARCH_GENERATIONS);
        assert_eq!(unset.brkga_config().population_size, SEARCH_POPULATION);
        assert_eq!(unset.brkga_config().max_generations, SEARCH_GENERATIONS);

        let set = Nester2D::new(
            Config::default()
                .with_population_size(100)
                .with_max_generations(200),
        );
        assert_eq!(set.ga_config().population_size, 100);
        assert_eq!(set.ga_config().max_generations, 200);
        assert_eq!(set.brkga_config().population_size, 100);
        assert_eq!(set.brkga_config().max_generations, 200);

        let small = Nester2D::new(Config::default().with_population_size(3));
        assert_eq!(small.brkga_config().population_size, 3, "was raised to 4");
        assert!(small.brkga_config().validate().is_ok());
    }

    /// The generation budget reaches the running search: every report of the
    /// genetic search carries it as the total.
    #[test]
    fn the_generation_budget_reaches_the_running_search() {
        use std::sync::Mutex;
        let reports = |config: Config| {
            let seen = Arc::new(Mutex::new(Vec::new()));
            let sink = seen.clone();
            Nester2D::new(
                config
                    .with_strategy(Strategy::GeneticAlgorithm)
                    .with_seed(7),
            )
            .solve_with_progress(
                &[Geometry2D::rectangle("R", 10.0, 10.0).with_quantity(4)],
                &Boundary2D::rectangle(50.0, 50.0),
                Box::new(move |info| {
                    sink.lock()
                        .expect("no report panicked")
                        .push((info.iteration, info.total_iterations))
                }),
            )
            .expect("solve");
            let seen = seen.lock().expect("no report panicked").clone();
            seen
        };
        let set = reports(Config::default().with_max_generations(80));
        assert!(!set.is_empty(), "the search reports its progress");
        assert!(
            set.iter().all(|&(at, total)| total == 80 && at <= 80),
            "{set:?}"
        );
        let unset = reports(Config::default());
        assert!(!unset.is_empty());
        assert!(unset.iter().all(|&(_, total)| total == 50), "{unset:?}");
    }

    /// A progress callback only observes a solve: the same request with the
    /// same seed gives the same layout through `solve` and
    /// `solve_with_progress`. The searches end on their generation budget,
    /// not on the clock, so the comparison does not depend on the machine.
    /// The progress path used to give the genetic search the whole time
    /// limit (none at all when it was 0) where `solve` gives it a per-strip
    /// budget.
    #[test]
    fn a_progress_callback_does_not_change_the_solve() {
        let geometries = [
            Geometry2D::rectangle("A", 30.0, 12.0)
                .with_quantity(3)
                .with_rotations_deg(vec![0.0, 90.0]),
            Geometry2D::rectangle("B", 18.0, 18.0).with_quantity(2),
            Geometry2D::new("L")
                .with_polygon(vec![
                    (0.0, 0.0),
                    (20.0, 0.0),
                    (20.0, 8.0),
                    (8.0, 8.0),
                    (8.0, 20.0),
                    (0.0, 20.0),
                ])
                .with_quantity(2)
                .with_rotations_deg(vec![0.0, 90.0, 180.0, 270.0]),
        ];
        let boundary = Boundary2D::rectangle(60.0, 50.0);
        /// A placement as compared: id, instance, position, rotation, mirrored.
        type Laid = (String, usize, Vec<f64>, Vec<f64>, bool);
        let layout = |r: &SolveResult<f64>| -> Vec<Laid> {
            r.placements
                .iter()
                .map(|p| {
                    (
                        p.geometry_id.clone(),
                        p.instance,
                        p.position.clone(),
                        p.rotation.clone(),
                        p.mirrored,
                    )
                })
                .collect()
        };
        for strategy in [
            Strategy::BottomLeftFill,
            Strategy::NfpGuided,
            Strategy::GeneticAlgorithm,
            Strategy::Brkga,
        ] {
            for time_limit_ms in [0, 60_000] {
                let config = Config::default()
                    .with_strategy(strategy)
                    .with_seed(11)
                    .with_population_size(6)
                    .with_max_generations(4)
                    .with_time_limit(time_limit_ms);
                let plain = Nester2D::new(config.clone())
                    .solve(&geometries, &boundary)
                    .expect("solve");
                let reported = Nester2D::new(config)
                    .solve_with_progress(&geometries, &boundary, Box::new(|_| {}))
                    .expect("solve with progress");
                let case = format!("{strategy:?}, time_limit_ms {time_limit_ms}");
                assert!(!plain.placements.is_empty(), "{case}");
                assert_eq!(layout(&plain), layout(&reported), "{case}");
                assert_eq!(plain.utilization, reported.utilization, "{case}");
                assert_eq!(plain.generations, reported.generations, "{case}");
                assert_eq!(plain.unplaced, reported.unplaced, "{case}");
            }
        }
    }

    #[test]
    fn test_simple_nesting() {
        let geometries = vec![
            Geometry2D::rectangle("R1", 20.0, 10.0).with_quantity(3),
            Geometry2D::rectangle("R2", 15.0, 15.0).with_quantity(2),
        ];

        let boundary = Boundary2D::rectangle(100.0, 50.0);
        let nester = Nester2D::default_config();

        let result = nester.solve(&geometries, &boundary).unwrap();

        assert!(result.utilization > 0.0);
        assert!(result.placements.len() <= 5); // 3 + 2 = 5 pieces
    }

    #[test]
    fn test_placement_within_bounds() {
        let geometries = vec![Geometry2D::rectangle("R1", 10.0, 10.0).with_quantity(4)];

        let boundary = Boundary2D::rectangle(50.0, 50.0);
        let config = Config::default().with_margin(5.0).with_spacing(2.0);
        let nester = Nester2D::new(config);

        let result = nester.solve(&geometries, &boundary).unwrap();

        // All pieces should be placed
        assert_eq!(result.placements.len(), 4);
        assert!(result.unplaced.is_empty());

        // Verify placements are within bounds (with margin)
        for p in &result.placements {
            assert!(p.position[0] >= 5.0);
            assert!(p.position[1] >= 5.0);
        }
    }

    #[test]
    fn test_nfp_guided_basic() {
        let geometries = vec![
            Geometry2D::rectangle("R1", 20.0, 10.0).with_quantity(2),
            Geometry2D::rectangle("R2", 15.0, 15.0).with_quantity(1),
        ];

        let boundary = Boundary2D::rectangle(100.0, 50.0);
        let config = Config::default().with_strategy(Strategy::NfpGuided);
        let nester = Nester2D::new(config);

        let result = nester.solve(&geometries, &boundary).unwrap();

        assert!(result.utilization > 0.0);
        assert_eq!(result.placements.len(), 3); // 2 + 1 = 3 pieces
        assert!(result.unplaced.is_empty());
    }

    #[test]
    fn test_nfp_guided_with_spacing() {
        let geometries = vec![Geometry2D::rectangle("R1", 10.0, 10.0).with_quantity(4)];

        let boundary = Boundary2D::rectangle(50.0, 50.0);
        let config = Config::default()
            .with_strategy(Strategy::NfpGuided)
            .with_margin(2.0)
            .with_spacing(3.0);
        let nester = Nester2D::new(config);

        let result = nester.solve(&geometries, &boundary).unwrap();

        // All pieces should be placed
        assert_eq!(result.placements.len(), 4);
        assert!(result.unplaced.is_empty());

        // Utilization should be positive
        assert!(result.utilization > 0.0);
    }

    #[test]
    fn test_nfp_guided_no_overlap() {
        let geometries = vec![Geometry2D::rectangle("R1", 20.0, 20.0).with_quantity(3)];

        let boundary = Boundary2D::rectangle(100.0, 100.0);
        let config = Config::default().with_strategy(Strategy::NfpGuided);
        let nester = Nester2D::new(config);

        let result = nester.solve(&geometries, &boundary).unwrap();

        assert_eq!(result.placements.len(), 3);

        // Verify no overlaps between placements
        for i in 0..result.placements.len() {
            for j in (i + 1)..result.placements.len() {
                let p1 = &result.placements[i];
                let p2 = &result.placements[j];

                // Simple AABB overlap check for rectangles
                let r1_min_x = p1.position[0];
                let r1_max_x = p1.position[0] + 20.0;
                let r1_min_y = p1.position[1];
                let r1_max_y = p1.position[1] + 20.0;

                let r2_min_x = p2.position[0];
                let r2_max_x = p2.position[0] + 20.0;
                let r2_min_y = p2.position[1];
                let r2_max_y = p2.position[1] + 20.0;

                // Check no overlap (with small tolerance for floating point)
                let overlaps_x = r1_min_x < r2_max_x - 0.01 && r1_max_x > r2_min_x + 0.01;
                let overlaps_y = r1_min_y < r2_max_y - 0.01 && r1_max_y > r2_min_y + 0.01;

                assert!(
                    !(overlaps_x && overlaps_y),
                    "Placements {} and {} overlap",
                    i,
                    j
                );
            }
        }
    }

    /// Returns true if any edge of `a` crosses any edge of `b` (world-space
    /// polygons). Sufficient for a placement-correctness regression check —
    /// the placement pipeline builds pieces edge-to-edge, so any real overlap
    /// between two placed pieces shows up as a boundary crossing.
    fn polygons_overlap(a: &[(f64, f64)], b: &[(f64, f64)]) -> bool {
        for i in 0..a.len() {
            let a1 = a[i];
            let a2 = a[(i + 1) % a.len()];
            for j in 0..b.len() {
                let b1 = b[j];
                let b2 = b[(j + 1) % b.len()];
                if crate::polygon_ops::segments_intersect(a1, a2, b1, b2) {
                    return true;
                }
            }
        }
        false
    }

    #[test]
    fn test_mirror_candidates_reflects_allow_flip() {
        let plain = Geometry2D::rectangle("R", 10.0, 10.0);
        assert_eq!(mirror_candidates(&plain), &[false]);

        let flippable = Geometry2D::rectangle("R", 10.0, 10.0).with_flip(true);
        assert_eq!(mirror_candidates(&flippable), &[false, true]);
    }

    /// Phase 1 (`allow_flip`/mirroring), placement-level correctness.
    ///
    /// Written while `Geometry2D::validate()` still rejected `allow_flip =
    /// true` unconditionally (that per-geometry gate was baked into the
    /// strategy function itself, so calling `nfp_guided_blf` directly
    /// couldn't bypass it either) — this replicated `nfp_guided_blf`'s own
    /// placement sequence one level down, using the exact primitives it
    /// calls (`compute_ifp_with_margin_and_mirror`, `compute_nfp_mirrored`,
    /// `find_bottom_left_placement`) — none of which call `.validate()` —
    /// to prove the thing that could actually go wrong: NFP-based collision
    /// avoidance between an unmirrored piece and a mirrored one.
    ///
    /// The gate is now open (Phase 4) and `nfp_guided_blf` *is* reachable
    /// end-to-end with `allow_flip = true`
    /// via `Nester2D::solve` — see
    /// `integration_tests::allow_flip_is_accepted_across_strategies` and
    /// `fuzz_robustness::placement_has_no_overlap_with_mirroring` for that
    /// path. This test stays as a focused, deterministic primitive-level
    /// regression; it is no longer the only way to exercise this.
    #[test]
    fn test_mirror_aware_placement_avoids_overlap() {
        // Nonzero spacing, matching how `nfp_guided_blf` actually calls these
        // primitives (`offset_nfp` by `self.config.spacing`) —
        // at zero spacing, BLF's own bottom-left search can legitimately
        // return pieces exactly edge-touching (a "kissing" placement is a
        // valid zero-gap packing, not an overlap), which would make this
        // test's boundary-crossing overlap check spuriously fail on a
        // shared-edge placement instead of the collision it's meant to catch.
        let spacing = 1.0;
        let geom = Geometry2D::l_shape("L", 30.0, 20.0, 20.0, 10.0);
        let boundary_polygon = vec![(0.0, 0.0), (65.0, 0.0), (65.0, 45.0), (0.0, 45.0)];
        let cache = NfpCache::new();

        // Place the first instance unmirrored, at the boundary's IFP origin.
        let ifp1 =
            compute_ifp_with_margin_and_mirror(&boundary_polygon, &geom, 0.0, 0.0, false).unwrap();
        let (x1, y1) = find_bottom_left_placement(&ifp1, &[], 1.0, PackingAxis::X)
            .expect("first piece must fit");
        let placed1 = PlacedGeometry::new(geom.clone(), (x1, y1), 0.0).with_mirrored(false);

        // Place the second instance MIRRORED, avoiding the first.
        let ifp2 =
            compute_ifp_with_margin_and_mirror(&boundary_polygon, &geom, 0.0, 0.0, true).unwrap();
        let nfp_at_origin = cache
            .get_or_compute_mirrored(("L", "L", 0.0, false, true), || {
                compute_nfp_mirrored(&placed1.geometry, &geom, 0.0, false, true)
            })
            .unwrap();
        let translated_nfp = translate_nfp(&nfp_at_origin, placed1.position);
        let expanded_nfp = offset_nfp(&translated_nfp, spacing);
        let (x2, y2) = find_bottom_left_placement(&ifp2, &[&expanded_nfp], 1.0, PackingAxis::X)
            .expect("mirrored second piece must fit avoiding the first");
        let placed2 = PlacedGeometry::new(geom.clone(), (x2, y2), 0.0).with_mirrored(true);

        assert_ne!(
            (x1, y1),
            (x2, y2),
            "mirrored placement must actually avoid the first piece's position"
        );
        assert!(
            !polygons_overlap(
                &placed1.translated_exterior(),
                &placed2.translated_exterior()
            ),
            "unmirrored piece at {:?} and mirrored piece at {:?} must not overlap",
            (x1, y1),
            (x2, y2)
        );
    }

    #[test]
    fn test_nfp_guided_utilization() {
        // Perfect fit: 4 rectangles of 25x25 in a 100x50 boundary
        let geometries = vec![Geometry2D::rectangle("R1", 25.0, 25.0).with_quantity(4)];

        let boundary = Boundary2D::rectangle(100.0, 50.0);
        let config = Config::default().with_strategy(Strategy::NfpGuided);
        let nester = Nester2D::new(config);

        let result = nester.solve(&geometries, &boundary).unwrap();

        // All pieces should be placed
        assert_eq!(result.placements.len(), 4);

        // Utilization should be 50% (4 * 625 = 2500 / 5000)
        assert!(result.utilization > 0.45);
    }

    #[test]
    fn test_polygon_centroid() {
        // Test the centroid calculation
        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() < 0.01);
        assert!((cy - 5.0).abs() < 0.01);

        let triangle = vec![(0.0, 0.0), (6.0, 0.0), (3.0, 6.0)];
        let (cx, cy) = polygon_centroid(&triangle);
        assert!((cx - 3.0).abs() < 0.01);
        assert!((cy - 2.0).abs() < 0.01);
    }

    #[test]
    fn test_ga_strategy_basic() {
        let geometries = vec![
            Geometry2D::rectangle("R1", 20.0, 10.0).with_quantity(2),
            Geometry2D::rectangle("R2", 15.0, 15.0).with_quantity(2),
        ];

        let boundary = Boundary2D::rectangle(100.0, 50.0);
        let config = Config::default().with_strategy(Strategy::GeneticAlgorithm);
        let nester = Nester2D::new(config);

        let result = nester.solve(&geometries, &boundary).unwrap();

        assert!(result.utilization > 0.0);
        assert!(!result.placements.is_empty());
        // GA should report generations and fitness
        assert!(result.generations.is_some());
        assert!(result.best_fitness.is_some());
        assert!(result.strategy == Some("GeneticAlgorithm".to_string()));
    }

    #[test]
    fn test_ga_strategy_all_placed() {
        // Easy case: 4 small rectangles in large boundary
        let geometries = vec![Geometry2D::rectangle("R1", 20.0, 20.0).with_quantity(4)];

        let boundary = Boundary2D::rectangle(100.0, 100.0);
        let config = Config::default().with_strategy(Strategy::GeneticAlgorithm);
        let nester = Nester2D::new(config);

        let result = nester.solve(&geometries, &boundary).unwrap();

        // All 4 pieces should fit
        assert_eq!(result.placements.len(), 4);
        assert!(result.unplaced.is_empty());
    }

    #[test]
    fn test_brkga_strategy_basic() {
        let geometries = vec![
            Geometry2D::rectangle("R1", 20.0, 10.0).with_quantity(2),
            Geometry2D::rectangle("R2", 15.0, 15.0).with_quantity(2),
        ];

        let boundary = Boundary2D::rectangle(100.0, 50.0);
        let config = Config::default().with_strategy(Strategy::Brkga);
        let nester = Nester2D::new(config);

        let result = nester.solve(&geometries, &boundary).unwrap();

        assert!(result.utilization > 0.0);
        assert!(!result.placements.is_empty());
        // BRKGA should report generations and fitness
        assert!(result.generations.is_some());
        assert!(result.best_fitness.is_some());
        assert!(result.strategy == Some("BRKGA".to_string()));
    }

    #[test]
    fn test_brkga_strategy_all_placed() {
        // Easy case: 4 small rectangles in large boundary
        let geometries = vec![Geometry2D::rectangle("R1", 20.0, 20.0).with_quantity(4)];

        let boundary = Boundary2D::rectangle(100.0, 100.0);
        // Use longer time limit to ensure BRKGA converges on all platforms
        let config = Config::default()
            .with_strategy(Strategy::Brkga)
            .with_time_limit(30000); // 30 seconds
        let nester = Nester2D::new(config);

        let result = nester.solve(&geometries, &boundary).unwrap();

        // BRKGA is stochastic; expect at least 3 of 4 pieces placed
        // (4 x 20x20 = 1600 area in 10000 boundary = 16% utilization, easy case)
        assert!(
            result.placements.len() >= 3,
            "Expected at least 3 placements, got {}",
            result.placements.len()
        );
    }

    #[test]
    fn test_gdrr_strategy_basic() {
        let geometries = vec![
            Geometry2D::rectangle("R1", 20.0, 10.0).with_quantity(2),
            Geometry2D::rectangle("R2", 15.0, 15.0).with_quantity(2),
        ];

        let boundary = Boundary2D::rectangle(100.0, 50.0);
        let config = Config::default().with_strategy(Strategy::Gdrr);
        let nester = Nester2D::new(config);

        let result = nester.solve(&geometries, &boundary).unwrap();

        assert!(result.utilization > 0.0);
        assert!(!result.placements.is_empty());
        // GDRR should report iterations and fitness
        assert!(result.iterations.is_some());
        assert!(result.best_fitness.is_some());
        assert!(result.strategy == Some("GDRR".to_string()));
    }

    #[test]
    fn test_gdrr_strategy_all_placed() {
        // Easy case: 4 small rectangles in large boundary
        let geometries = vec![Geometry2D::rectangle("R1", 20.0, 20.0).with_quantity(4)];

        let boundary = Boundary2D::rectangle(100.0, 100.0);
        let config = Config::default().with_strategy(Strategy::Gdrr);
        let nester = Nester2D::new(config);

        let result = nester.solve(&geometries, &boundary).unwrap();

        // All 4 pieces should fit
        assert_eq!(result.placements.len(), 4);
        assert!(result.unplaced.is_empty());
    }

    #[test]
    fn test_alns_strategy_basic() {
        let geometries = vec![
            Geometry2D::rectangle("R1", 20.0, 10.0).with_quantity(2),
            Geometry2D::rectangle("R2", 15.0, 15.0).with_quantity(2),
        ];

        let boundary = Boundary2D::rectangle(100.0, 50.0);
        let config = Config::default().with_strategy(Strategy::Alns);
        let nester = Nester2D::new(config);

        let result = nester.solve(&geometries, &boundary).unwrap();

        assert!(result.utilization > 0.0);
        assert!(!result.placements.is_empty());
        // ALNS should report iterations and fitness
        assert!(result.iterations.is_some());
        assert!(result.best_fitness.is_some());
        assert!(result.strategy == Some("ALNS".to_string()));
    }

    #[test]
    fn test_alns_strategy_all_placed() {
        // Easy case: 4 small rectangles in large boundary
        let geometries = vec![Geometry2D::rectangle("R1", 20.0, 20.0).with_quantity(4)];

        let boundary = Boundary2D::rectangle(100.0, 100.0);
        let config = Config::default().with_strategy(Strategy::Alns);
        let nester = Nester2D::new(config);

        let result = nester.solve(&geometries, &boundary).unwrap();

        // All 4 pieces should fit
        assert_eq!(result.placements.len(), 4);
        assert!(result.unplaced.is_empty());
    }

    #[test]
    fn test_blf_rotation_optimization() {
        // Test that BLF uses rotation to optimize placement
        // A 30x10 rectangle can fit better in a narrow strip when rotated 90 degrees
        let geometries = vec![Geometry2D::rectangle("R1", 30.0, 10.0)
                .with_rotations(vec![0.0, std::f64::consts::FRAC_PI_2]) // 0 and 90 degrees
                .with_quantity(3)];

        // Strip that's 35 wide: 30x10 won't fit two side-by-side at 0 deg
        // But two 10x30 (rotated 90 deg) can fit vertically in 95 height
        let boundary = Boundary2D::rectangle(35.0, 95.0);
        let nester = Nester2D::default_config();

        let result = nester.solve(&geometries, &boundary).unwrap();

        // All 3 pieces should be placed (by rotating)
        assert_eq!(
            result.placements.len(),
            3,
            "All pieces should be placed with rotation optimization"
        );
        assert!(result.unplaced.is_empty());
    }

    #[test]
    fn test_blf_selects_best_rotation() {
        // Verify BLF selects optimal rotation, not just the first one
        let geometries = vec![Geometry2D::rectangle("R1", 40.0, 10.0)
                .with_rotations(vec![0.0, std::f64::consts::FRAC_PI_2]) // 0 and 90 degrees
                .with_quantity(2)];

        // In a 45x50 boundary:
        // - At 0 deg: 40x10, only one fits horizontally (40 < 45), next row needed
        // - At 90 deg: 10x40, two can fit side-by-side (10+10 < 45) in one row
        let boundary = Boundary2D::rectangle(45.0, 50.0);
        let nester = Nester2D::default_config();

        let result = nester.solve(&geometries, &boundary).unwrap();

        assert_eq!(result.placements.len(), 2);
        assert!(result.unplaced.is_empty());
    }

    #[test]
    fn test_progress_callback_blf() {
        use std::sync::atomic::{AtomicUsize, Ordering};
        use std::sync::Arc;

        let geometries = vec![Geometry2D::rectangle("R1", 10.0, 10.0).with_quantity(4)];
        let boundary = Boundary2D::rectangle(50.0, 50.0);
        let config = Config::default().with_strategy(Strategy::BottomLeftFill);
        let nester = Nester2D::new(config);

        let callback_count = Arc::new(AtomicUsize::new(0));
        let callback_count_clone = callback_count.clone();
        let last_items_placed = Arc::new(AtomicUsize::new(0));
        let last_items_placed_clone = last_items_placed.clone();

        let callback: ProgressCallback = Box::new(move |info| {
            callback_count_clone.fetch_add(1, Ordering::Relaxed);
            last_items_placed_clone.store(info.items_placed, Ordering::Relaxed);
        });

        let result = nester
            .solve_with_progress(&geometries, &boundary, callback)
            .unwrap();

        // Verify callback was called (at least once per piece + initial + final)
        let count = callback_count.load(Ordering::Relaxed);
        assert!(
            count >= 5,
            "Expected at least 5 callbacks (1 initial + 4 pieces + 1 final), got {}",
            count
        );

        // Verify final items_placed
        let final_placed = last_items_placed.load(Ordering::Relaxed);
        assert_eq!(final_placed, 4, "Should report 4 items placed");

        // Verify result
        assert_eq!(result.placements.len(), 4);
    }

    #[test]
    fn test_progress_callback_nfp() {
        use std::sync::atomic::{AtomicUsize, Ordering};
        use std::sync::Arc;

        let geometries = vec![Geometry2D::rectangle("R1", 10.0, 10.0).with_quantity(2)];
        let boundary = Boundary2D::rectangle(50.0, 50.0);
        let config = Config::default().with_strategy(Strategy::NfpGuided);
        let nester = Nester2D::new(config);

        let callback_count = Arc::new(AtomicUsize::new(0));
        let callback_count_clone = callback_count.clone();

        let callback: ProgressCallback = Box::new(move |info| {
            callback_count_clone.fetch_add(1, Ordering::Relaxed);
            assert!(info.items_placed <= info.total_items);
        });

        let result = nester
            .solve_with_progress(&geometries, &boundary, callback)
            .unwrap();

        // Verify callback was called
        let count = callback_count.load(Ordering::Relaxed);
        assert!(count >= 3, "Expected at least 3 callbacks, got {}", count);

        // Verify result
        assert_eq!(result.placements.len(), 2);
    }

    #[test]
    fn test_time_limit_honored() {
        // Create many geometries to ensure BLF takes measurable time
        let geometries: Vec<Geometry2D> = (0..100)
            .map(|i| Geometry2D::rectangle(format!("R{}", i), 5.0, 5.0))
            .collect();
        let boundary = Boundary2D::rectangle(1000.0, 1000.0);

        // Set a very short time limit (1ms) to ensure timeout
        let config = Config::default()
            .with_strategy(Strategy::BottomLeftFill)
            .with_time_limit(1);
        let nester = Nester2D::new(config);

        let result = nester.solve(&geometries, &boundary).unwrap();

        // With such a short time limit, we may not place all items
        // The test verifies that the solver respects the time limit
        assert!(
            result.computation_time_ms <= 100, // Allow some margin for overhead
            "Computation took too long: {}ms (expected <= 100ms with 1ms limit)",
            result.computation_time_ms
        );
    }

    #[test]
    fn test_time_limit_zero_unlimited() {
        // time_limit_ms = 0 means unlimited
        let geometries = vec![Geometry2D::rectangle("R1", 10.0, 10.0).with_quantity(4)];
        let boundary = Boundary2D::rectangle(50.0, 50.0);

        let config = Config::default()
            .with_strategy(Strategy::BottomLeftFill)
            .with_time_limit(0); // Unlimited
        let nester = Nester2D::new(config);

        let result = nester.solve(&geometries, &boundary).unwrap();

        // Should place all items (no early exit)
        assert_eq!(result.placements.len(), 4);
    }

    #[test]
    fn test_blf_bounds_clamping() {
        // Test that BLF correctly clamps placements within boundary
        // Create a shape with non-zero g_min (similar to Gear shape)
        // Gear-like: x ranges from 5 to 95 (width=90), y from 5 to 95 (height=90)
        let gear_like = Geometry2D::new("gear")
            .with_polygon(vec![
                (50.0, 5.0), // Bottom
                (65.0, 15.0),
                (77.0, 18.0),
                (80.0, 32.0),
                (95.0, 50.0), // Right
                (80.0, 68.0),
                (77.0, 82.0),
                (65.0, 85.0),
                (50.0, 95.0), // Top
                (35.0, 85.0),
                (23.0, 82.0),
                (20.0, 68.0),
                (5.0, 50.0), // Left (min_x = 5)
                (20.0, 32.0),
                (23.0, 18.0),
                (35.0, 15.0),
            ])
            .with_quantity(1);

        // Boundary is 100x100
        let boundary = Boundary2D::rectangle(100.0, 100.0);

        let config = Config::default().with_strategy(Strategy::BottomLeftFill);
        let nester = Nester2D::new(config);

        let result = nester
            .solve(std::slice::from_ref(&gear_like), &boundary)
            .unwrap();

        assert_eq!(result.placements.len(), 1);
        let placement = &result.placements[0];

        // Origin position
        let origin_x = placement.position[0];
        let origin_y = placement.position[1];

        // Get rotation from placement (2D rotation is a single value in Vec)
        let rotation = placement.rotation.first().copied().unwrap_or(0.0);

        // Get AABB at rotation
        let (g_min, g_max) = gear_like.aabb_at_rotation(rotation);

        // Actual geometry bounds after placement
        let actual_min_x = origin_x + g_min[0];
        let actual_max_x = origin_x + g_max[0];
        let actual_min_y = origin_y + g_min[1];
        let actual_max_y = origin_y + g_max[1];

        // All edges should be within boundary [0, 100]
        assert!(
            actual_min_x >= 0.0,
            "Left edge {} should be >= 0",
            actual_min_x
        );
        assert!(
            actual_max_x <= 100.0,
            "Right edge {} should be <= 100",
            actual_max_x
        );
        assert!(
            actual_min_y >= 0.0,
            "Bottom edge {} should be >= 0",
            actual_min_y
        );
        assert!(
            actual_max_y <= 100.0,
            "Top edge {} should be <= 100",
            actual_max_y
        );
    }

    #[test]
    fn test_blf_bounds_clamping_many_pieces() {
        // Test BLF bounds clamping with many pieces to trigger row overflow
        // This mimics the actual failing case from test_blf.py
        let gear_like = Geometry2D::new("gear")
            .with_polygon(vec![
                (50.0, 5.0),
                (65.0, 15.0),
                (77.0, 18.0),
                (80.0, 32.0),
                (95.0, 50.0),
                (80.0, 68.0),
                (77.0, 82.0),
                (65.0, 85.0),
                (50.0, 95.0),
                (35.0, 85.0),
                (23.0, 82.0),
                (20.0, 68.0),
                (5.0, 50.0),
                (20.0, 32.0),
                (23.0, 18.0),
                (35.0, 15.0),
            ])
            .with_quantity(13); // Same as Gear (shape 8) in test_blf.py

        // Boundary is 500x500 like the test
        let boundary = Boundary2D::rectangle(500.0, 500.0);

        let config = Config::default().with_strategy(Strategy::BottomLeftFill);
        let nester = Nester2D::new(config);

        let result = nester
            .solve(std::slice::from_ref(&gear_like), &boundary)
            .unwrap();

        // Check that ALL placements are within bounds
        for (i, placement) in result.placements.iter().enumerate() {
            let origin_x = placement.position[0];
            let origin_y = placement.position[1];
            let rotation = placement.rotation.first().copied().unwrap_or(0.0);

            let (g_min, g_max) = gear_like.aabb_at_rotation(rotation);

            let actual_min_x = origin_x + g_min[0];
            let actual_max_x = origin_x + g_max[0];
            let actual_min_y = origin_y + g_min[1];
            let actual_max_y = origin_y + g_max[1];

            assert!(
                actual_min_x >= -0.01,
                "Piece {}: Left edge {} should be >= 0",
                i,
                actual_min_x
            );
            assert!(
                actual_max_x <= 500.01,
                "Piece {}: Right edge {} should be <= 500",
                i,
                actual_max_x
            );
            assert!(
                actual_min_y >= -0.01,
                "Piece {}: Bottom edge {} should be >= 0",
                i,
                actual_min_y
            );
            assert!(
                actual_max_y <= 500.01,
                "Piece {}: Top edge {} should be <= 500",
                i,
                actual_max_y
            );
        }
    }

    #[test]
    fn test_blf_bounds_trace() {
        // Debug test: trace through BLF to understand why clamping doesn't work
        let gear = Geometry2D::new("gear").with_polygon(vec![
            (50.0, 5.0),
            (65.0, 15.0),
            (77.0, 18.0),
            (80.0, 32.0),
            (95.0, 50.0),
            (80.0, 68.0),
            (77.0, 82.0),
            (65.0, 85.0),
            (50.0, 95.0),
            (35.0, 85.0),
            (23.0, 82.0),
            (20.0, 68.0),
            (5.0, 50.0),
            (20.0, 32.0),
            (23.0, 18.0),
            (35.0, 15.0),
        ]);

        // Verify AABB
        let (g_min, g_max) = gear.aabb();
        println!("Gear AABB: min={:?}, max={:?}", g_min, g_max);
        assert!(
            (g_min[0] - 5.0).abs() < 0.01,
            "g_min[0] should be 5, got {}",
            g_min[0]
        );
        assert!(
            (g_max[0] - 95.0).abs() < 0.01,
            "g_max[0] should be 95, got {}",
            g_max[0]
        );

        // Verify valid origin range for 500x500 boundary
        let b_max_x = 500.0;
        let margin = 0.0;
        let max_valid_x = b_max_x - margin - g_max[0];
        println!(
            "max_valid_x = {} - {} - {} = {}",
            b_max_x, margin, g_max[0], max_valid_x
        );
        assert!(
            (max_valid_x - 405.0).abs() < 0.01,
            "max_valid_x should be 405, got {}",
            max_valid_x
        );

        // Run BLF and check the result
        let boundary = Boundary2D::rectangle(500.0, 500.0);
        let config = Config::default().with_strategy(Strategy::BottomLeftFill);
        let nester = Nester2D::new(config);

        let result = nester
            .solve(&[gear.clone().with_quantity(1)], &boundary)
            .unwrap();

        assert_eq!(result.placements.len(), 1);
        let p = &result.placements[0];
        let origin_x = p.position[0];
        let rotation = p.rotation.first().copied().unwrap_or(0.0);

        let (g_min_r, g_max_r) = gear.aabb_at_rotation(rotation);
        let actual_max_x = origin_x + g_max_r[0];

        println!("Placement: origin_x={}, rotation={}", origin_x, rotation);
        println!(
            "At rotation {}: g_min={:?}, g_max={:?}",
            rotation, g_min_r, g_max_r
        );
        println!(
            "Actual max x: {} + {} = {}",
            origin_x, g_max_r[0], actual_max_x
        );

        assert!(
            actual_max_x <= 500.01,
            "Geometry exceeds boundary: max_x={} > 500",
            actual_max_x
        );
    }

    #[test]
    fn test_blf_bounds_many_pieces_direct() {
        // Test with many pieces to trigger the boundary violation
        let gear = Geometry2D::new("gear")
            .with_polygon(vec![
                (50.0, 5.0),
                (65.0, 15.0),
                (77.0, 18.0),
                (80.0, 32.0),
                (95.0, 50.0),
                (80.0, 68.0),
                (77.0, 82.0),
                (65.0, 85.0),
                (50.0, 95.0),
                (35.0, 85.0),
                (23.0, 82.0),
                (20.0, 68.0),
                (5.0, 50.0),
                (20.0, 32.0),
                (23.0, 18.0),
                (35.0, 15.0),
            ])
            .with_quantity(25); // Many pieces

        let boundary = Boundary2D::rectangle(500.0, 500.0);
        let config = Config::default().with_strategy(Strategy::BottomLeftFill);
        let nester = Nester2D::new(config);

        let result = nester
            .solve(std::slice::from_ref(&gear), &boundary)
            .unwrap();

        println!("Placed {} pieces", result.placements.len());

        // Check all placements
        for (i, p) in result.placements.iter().enumerate() {
            let origin_x = p.position[0];
            let origin_y = p.position[1];
            let rotation = p.rotation.first().copied().unwrap_or(0.0);

            let (g_min_r, g_max_r) = gear.aabb_at_rotation(rotation);

            let actual_min_x = origin_x + g_min_r[0];
            let actual_max_x = origin_x + g_max_r[0];
            let actual_min_y = origin_y + g_min_r[1];
            let actual_max_y = origin_y + g_max_r[1];

            println!(
                "Piece {}: origin=({:.1}, {:.1}), rot={:.2}, bounds=[{:.1},{:.1}]x[{:.1},{:.1}]",
                i,
                origin_x,
                origin_y,
                rotation,
                actual_min_x,
                actual_max_x,
                actual_min_y,
                actual_max_y
            );

            assert!(
                actual_max_x <= 500.01,
                "Piece {}: Right edge {} > 500",
                i,
                actual_max_x
            );
            assert!(
                actual_max_y <= 500.01,
                "Piece {}: Top edge {} > 500",
                i,
                actual_max_y
            );
        }
    }

    #[test]
    fn test_blf_bounds_multi_strip() {
        // Test with solve_multi_strip which is what benchmark runner uses
        let gear = Geometry2D::new("gear")
            .with_polygon(vec![
                (50.0, 5.0),
                (65.0, 15.0),
                (77.0, 18.0),
                (80.0, 32.0),
                (95.0, 50.0),
                (80.0, 68.0),
                (77.0, 82.0),
                (65.0, 85.0),
                (50.0, 95.0),
                (35.0, 85.0),
                (23.0, 82.0),
                (20.0, 68.0),
                (5.0, 50.0),
                (20.0, 32.0),
                (23.0, 18.0),
                (35.0, 15.0),
            ])
            .with_quantity(50); // Many pieces to force multiple strips

        let boundary = Boundary2D::rectangle(500.0, 500.0);
        let config = Config::default().with_strategy(Strategy::BottomLeftFill);
        let nester = Nester2D::new(config);

        // Use solve_multi_strip like benchmark runner does
        let result = nester
            .solve_multi_strip(std::slice::from_ref(&gear), &boundary)
            .unwrap();

        println!(
            "Placed {} pieces across {} strips",
            result.placements.len(),
            result.boundaries_used
        );

        // Check all placements - within their respective strips
        let strip_width = 500.0;
        for (i, p) in result.placements.iter().enumerate() {
            let origin_x = p.position[0];
            let origin_y = p.position[1];
            let rotation = p.rotation.first().copied().unwrap_or(0.0);
            let strip_idx = p.boundary_index;

            // Calculate local position within strip
            let local_x = origin_x - (strip_idx as f64 * strip_width);

            let (_g_min_r, g_max_r) = gear.aabb_at_rotation(rotation);

            let local_max_x = local_x + g_max_r[0];
            let local_max_y = origin_y + g_max_r[1];

            println!(
                "Piece {}: strip={}, origin=({:.1}, {:.1}), local_x={:.1}, rot={:.2}, local_max_x={:.1}",
                i, strip_idx, origin_x, origin_y, local_x, rotation, local_max_x
            );

            assert!(
                local_max_x <= 500.01,
                "Piece {}: In strip {}, local right edge {:.1} > 500",
                i,
                strip_idx,
                local_max_x
            );
            assert!(
                local_max_y <= 500.01,
                "Piece {}: Top edge {:.1} > 500",
                i,
                local_max_y
            );
        }
    }

    #[test]
    fn test_blf_bounds_mixed_shapes() {
        // Replicate test_blf.py with all 9 shapes
        let shapes = vec![
            // Shape 0: Rounded rectangle (demand 2)
            Geometry2D::new("shape0")
                .with_polygon(vec![
                    (0.0, 0.0),
                    (180.0, 0.0),
                    (195.0, 15.0),
                    (200.0, 50.0),
                    (200.0, 150.0),
                    (195.0, 185.0),
                    (180.0, 200.0),
                    (20.0, 200.0),
                    (5.0, 185.0),
                    (0.0, 150.0),
                    (0.0, 50.0),
                    (5.0, 15.0),
                ])
                .with_quantity(2),
            // Shape 1: Circular-ish (demand 4)
            Geometry2D::new("shape1")
                .with_polygon(vec![
                    (60.0, 0.0),
                    (85.0, 7.0),
                    (104.0, 25.0),
                    (118.0, 50.0),
                    (120.0, 60.0),
                    (118.0, 70.0),
                    (104.0, 95.0),
                    (85.0, 113.0),
                    (60.0, 120.0),
                    (35.0, 113.0),
                    (16.0, 95.0),
                    (2.0, 70.0),
                    (0.0, 60.0),
                    (2.0, 50.0),
                    (16.0, 25.0),
                    (35.0, 7.0),
                ])
                .with_quantity(4),
            // Shape 2: L-shape (demand 6)
            Geometry2D::new("shape2")
                .with_polygon(vec![
                    (0.0, 0.0),
                    (80.0, 0.0),
                    (80.0, 20.0),
                    (20.0, 20.0),
                    (20.0, 80.0),
                    (0.0, 80.0),
                ])
                .with_quantity(6),
            // Shape 3: Triangle (demand 6)
            Geometry2D::new("shape3")
                .with_polygon(vec![(0.0, 0.0), (70.0, 0.0), (0.0, 70.0)])
                .with_quantity(6),
            // Shape 4: Rectangle (demand 4)
            Geometry2D::new("shape4")
                .with_polygon(vec![(0.0, 0.0), (120.0, 0.0), (120.0, 60.0), (0.0, 60.0)])
                .with_quantity(4),
            // Shape 5: Hexagon (demand 8)
            Geometry2D::new("shape5")
                .with_polygon(vec![
                    (15.0, 0.0),
                    (45.0, 0.0),
                    (60.0, 26.0),
                    (45.0, 52.0),
                    (15.0, 52.0),
                    (0.0, 26.0),
                ])
                .with_quantity(8),
            // Shape 6: T-shape (demand 4)
            Geometry2D::new("shape6")
                .with_polygon(vec![
                    (0.0, 0.0),
                    (90.0, 0.0),
                    (90.0, 12.0),
                    (55.0, 12.0),
                    (55.0, 60.0),
                    (35.0, 60.0),
                    (35.0, 12.0),
                    (0.0, 12.0),
                ])
                .with_quantity(4),
            // Shape 7: Rounded square (demand 3)
            Geometry2D::new("shape7")
                .with_polygon(vec![
                    (0.0, 10.0),
                    (10.0, 0.0),
                    (70.0, 0.0),
                    (80.0, 10.0),
                    (80.0, 70.0),
                    (70.0, 80.0),
                    (10.0, 80.0),
                    (0.0, 70.0),
                ])
                .with_quantity(3),
            // Shape 8: Gear (demand 13) - the problematic shape
            Geometry2D::new("shape8_gear")
                .with_polygon(vec![
                    (50.0, 5.0),
                    (65.0, 15.0),
                    (77.0, 18.0),
                    (80.0, 32.0),
                    (95.0, 50.0),
                    (80.0, 68.0),
                    (77.0, 82.0),
                    (65.0, 85.0),
                    (50.0, 95.0),
                    (35.0, 85.0),
                    (23.0, 82.0),
                    (20.0, 68.0),
                    (5.0, 50.0),
                    (20.0, 32.0),
                    (23.0, 18.0),
                    (35.0, 15.0),
                ])
                .with_quantity(13),
        ];

        // Total: 2+4+6+6+4+8+4+3+13 = 50 pieces
        let boundary = Boundary2D::rectangle(500.0, 500.0);
        let config = Config::default().with_strategy(Strategy::BottomLeftFill);
        let nester = Nester2D::new(config);

        let result = nester.solve_multi_strip(&shapes, &boundary).unwrap();

        println!(
            "Placed {} pieces across {} strips",
            result.placements.len(),
            result.boundaries_used
        );

        // Check placements for Gear (shape8) specifically
        let strip_width = 500.0;
        let gear_aabb = shapes[8].aabb();
        println!("Gear AABB: min={:?}, max={:?}", gear_aabb.0, gear_aabb.1);

        let mut violations = Vec::new();
        for p in &result.placements {
            if p.geometry_id.as_str().starts_with("shape8") {
                let origin_x = p.position[0];
                let _origin_y = p.position[1];
                let rotation = p.rotation.first().copied().unwrap_or(0.0);
                let strip_idx = p.boundary_index;
                let local_x = origin_x - (strip_idx as f64 * strip_width);

                let (_g_min_r, g_max_r) = shapes[8].aabb_at_rotation(rotation);
                let local_max_x = local_x + g_max_r[0];

                println!(
                    "{}: strip={}, local_x={:.1}, rot={:.2}, local_max_x={:.1}",
                    p.geometry_id, strip_idx, local_x, rotation, local_max_x
                );

                if local_max_x > 500.01 {
                    violations.push((p.geometry_id.clone(), strip_idx, local_x, local_max_x));
                }
            }
        }

        assert!(
            violations.is_empty(),
            "Found {} Gear pieces exceeding boundary: {:?}",
            violations.len(),
            violations
        );
    }

    #[test]
    fn test_blf_bounds_expanded_like_benchmark() {
        // Replicate EXACTLY how benchmark runner creates geometries:
        // Each piece is a separate Geometry2D with quantity=1
        // (vertices, demand, allowed_rotations_deg)
        type ShapeDef = (Vec<(f64, f64)>, usize, Vec<f64>);
        let shape_defs: Vec<ShapeDef> = vec![
            (
                vec![
                    (0.0, 0.0),
                    (180.0, 0.0),
                    (195.0, 15.0),
                    (200.0, 50.0),
                    (200.0, 150.0),
                    (195.0, 185.0),
                    (180.0, 200.0),
                    (20.0, 200.0),
                    (5.0, 185.0),
                    (0.0, 150.0),
                    (0.0, 50.0),
                    (5.0, 15.0),
                ],
                2,
                vec![0.0, 90.0, 180.0, 270.0],
            ),
            (
                vec![
                    (60.0, 0.0),
                    (85.0, 7.0),
                    (104.0, 25.0),
                    (118.0, 50.0),
                    (120.0, 60.0),
                    (118.0, 70.0),
                    (104.0, 95.0),
                    (85.0, 113.0),
                    (60.0, 120.0),
                    (35.0, 113.0),
                    (16.0, 95.0),
                    (2.0, 70.0),
                    (0.0, 60.0),
                    (2.0, 50.0),
                    (16.0, 25.0),
                    (35.0, 7.0),
                ],
                4,
                vec![0.0, 45.0, 90.0, 135.0],
            ),
            (
                vec![
                    (0.0, 0.0),
                    (80.0, 0.0),
                    (80.0, 20.0),
                    (20.0, 20.0),
                    (20.0, 80.0),
                    (0.0, 80.0),
                ],
                6,
                vec![0.0, 90.0, 180.0, 270.0],
            ),
            (
                vec![(0.0, 0.0), (70.0, 0.0), (0.0, 70.0)],
                6,
                vec![0.0, 90.0, 180.0, 270.0],
            ),
            (
                vec![(0.0, 0.0), (120.0, 0.0), (120.0, 60.0), (0.0, 60.0)],
                4,
                vec![0.0, 90.0],
            ),
            (
                vec![
                    (15.0, 0.0),
                    (45.0, 0.0),
                    (60.0, 26.0),
                    (45.0, 52.0),
                    (15.0, 52.0),
                    (0.0, 26.0),
                ],
                8,
                vec![0.0, 60.0, 120.0],
            ),
            (
                vec![
                    (0.0, 0.0),
                    (90.0, 0.0),
                    (90.0, 12.0),
                    (55.0, 12.0),
                    (55.0, 60.0),
                    (35.0, 60.0),
                    (35.0, 12.0),
                    (0.0, 12.0),
                ],
                4,
                vec![0.0, 90.0, 180.0, 270.0],
            ),
            (
                vec![
                    (0.0, 10.0),
                    (10.0, 0.0),
                    (70.0, 0.0),
                    (80.0, 10.0),
                    (80.0, 70.0),
                    (70.0, 80.0),
                    (10.0, 80.0),
                    (0.0, 70.0),
                ],
                3,
                vec![0.0, 90.0],
            ),
            // Shape 8: Gear - with all 8 rotations
            (
                vec![
                    (50.0, 5.0),
                    (65.0, 15.0),
                    (77.0, 18.0),
                    (80.0, 32.0),
                    (95.0, 50.0),
                    (80.0, 68.0),
                    (77.0, 82.0),
                    (65.0, 85.0),
                    (50.0, 95.0),
                    (35.0, 85.0),
                    (23.0, 82.0),
                    (20.0, 68.0),
                    (5.0, 50.0),
                    (20.0, 32.0),
                    (23.0, 18.0),
                    (35.0, 15.0),
                ],
                13,
                vec![0.0, 45.0, 90.0, 135.0, 180.0, 225.0, 270.0, 315.0],
            ),
        ];

        // Expand like benchmark runner: each piece is separate geometry
        let mut geometries = Vec::new();
        let mut piece_id = 0;
        for (vertices, demand, rotations) in shape_defs.iter() {
            for _ in 0..*demand {
                let geom = Geometry2D::new(format!("piece_{}", piece_id))
                    .with_polygon(vertices.clone())
                    .with_rotations_deg(rotations.clone());
                geometries.push(geom);
                piece_id += 1;
            }
        }

        // Store gear AABB for checking
        let gear_geom = Geometry2D::new("gear_check").with_polygon(shape_defs[8].0.clone());
        let (gear_min, gear_max) = gear_geom.aabb();
        println!("Gear AABB: min={:?}, max={:?}", gear_min, gear_max);

        let boundary = Boundary2D::rectangle(500.0, 500.0);
        let config = Config::default().with_strategy(Strategy::BottomLeftFill);
        let nester = Nester2D::new(config);

        let result = nester.solve_multi_strip(&geometries, &boundary).unwrap();

        println!(
            "Placed {} pieces across {} strips",
            result.placements.len(),
            result.boundaries_used
        );

        // Check Gear placements (piece_37 to piece_49)
        let strip_width = 500.0;
        let mut violations = Vec::new();

        for p in &result.placements {
            let id_num: usize = p
                .geometry_id
                .as_str()
                .strip_prefix("piece_")
                .and_then(|s| s.parse().ok())
                .unwrap_or(0);

            // piece_37 to piece_49 are Gear shapes
            if (37..=49).contains(&id_num) {
                let origin_x = p.position[0];
                let rotation = p.rotation.first().copied().unwrap_or(0.0);
                let strip_idx = p.boundary_index;
                let local_x = origin_x - (strip_idx as f64 * strip_width);

                let (_, g_max_r) = gear_geom.aabb_at_rotation(rotation);
                let local_max_x = local_x + g_max_r[0];

                println!(
                    "{}: strip={}, local_x={:.1}, rot={:.2}, local_max_x={:.1}",
                    p.geometry_id, strip_idx, local_x, rotation, local_max_x
                );

                if local_max_x > 500.01 {
                    violations.push((p.geometry_id.clone(), strip_idx, local_x, local_max_x));
                }
            }
        }

        assert!(
            violations.is_empty(),
            "Found {} Gear pieces exceeding boundary: {:?}",
            violations.len(),
            violations
        );
    }

    /// Helper function to check if two AABBs overlap
    fn aabbs_overlap(
        a_min: [f64; 2],
        a_max: [f64; 2],
        b_min: [f64; 2],
        b_max: [f64; 2],
        tolerance: f64,
    ) -> bool {
        // Two AABBs overlap if they overlap on both axes
        let x_overlap = a_min[0] < b_max[0] - tolerance && a_max[0] > b_min[0] + tolerance;
        let y_overlap = a_min[1] < b_max[1] - tolerance && a_max[1] > b_min[1] + tolerance;
        x_overlap && y_overlap
    }

    /// Comprehensive test for all strategies - checks boundary and overlap violations
    #[test]
    fn test_all_strategies_boundary_and_overlap() {
        use std::collections::HashMap;

        // Create test shapes similar to demo
        let shapes = vec![
            Geometry2D::new("shape0")
                .with_polygon(vec![
                    (0.0, 0.0),
                    (180.0, 0.0),
                    (195.0, 15.0),
                    (200.0, 50.0),
                    (200.0, 150.0),
                    (195.0, 185.0),
                    (180.0, 200.0),
                    (20.0, 200.0),
                    (5.0, 185.0),
                    (0.0, 150.0),
                    (0.0, 50.0),
                    (5.0, 15.0),
                ])
                .with_rotations_deg(vec![0.0, 90.0, 180.0, 270.0])
                .with_quantity(2),
            Geometry2D::new("shape1_flange")
                .with_polygon(vec![
                    (60.0, 0.0),
                    (85.0, 7.0),
                    (104.0, 25.0),
                    (118.0, 50.0),
                    (120.0, 60.0),
                    (118.0, 70.0),
                    (104.0, 95.0),
                    (85.0, 113.0),
                    (60.0, 120.0),
                    (35.0, 113.0),
                    (16.0, 95.0),
                    (2.0, 70.0),
                    (0.0, 60.0),
                    (2.0, 50.0),
                    (16.0, 25.0),
                    (35.0, 7.0),
                ])
                .with_rotations_deg(vec![0.0, 45.0, 90.0, 135.0])
                .with_quantity(4),
            Geometry2D::new("shape2_lbracket")
                .with_polygon(vec![
                    (0.0, 0.0),
                    (80.0, 0.0),
                    (80.0, 20.0),
                    (20.0, 20.0),
                    (20.0, 80.0),
                    (0.0, 80.0),
                ])
                .with_rotations_deg(vec![0.0, 90.0, 180.0, 270.0])
                .with_quantity(6),
            Geometry2D::new("shape3_triangle")
                .with_polygon(vec![(0.0, 0.0), (70.0, 0.0), (0.0, 70.0)])
                .with_rotations_deg(vec![0.0, 90.0, 180.0, 270.0])
                .with_quantity(6),
            Geometry2D::new("shape4_rect")
                .with_polygon(vec![(0.0, 0.0), (120.0, 0.0), (120.0, 60.0), (0.0, 60.0)])
                .with_rotations_deg(vec![0.0, 90.0])
                .with_quantity(4),
            Geometry2D::new("shape5_hexagon")
                .with_polygon(vec![
                    (15.0, 0.0),
                    (45.0, 0.0),
                    (60.0, 26.0),
                    (45.0, 52.0),
                    (15.0, 52.0),
                    (0.0, 26.0),
                ])
                .with_rotations_deg(vec![0.0, 60.0, 120.0])
                .with_quantity(8),
            Geometry2D::new("shape6_tstiff")
                .with_polygon(vec![
                    (0.0, 0.0),
                    (90.0, 0.0),
                    (90.0, 12.0),
                    (55.0, 12.0),
                    (55.0, 60.0),
                    (35.0, 60.0),
                    (35.0, 12.0),
                    (0.0, 12.0),
                ])
                .with_rotations_deg(vec![0.0, 90.0, 180.0, 270.0])
                .with_quantity(4),
            Geometry2D::new("shape7_mount")
                .with_polygon(vec![
                    (0.0, 10.0),
                    (10.0, 0.0),
                    (70.0, 0.0),
                    (80.0, 10.0),
                    (80.0, 70.0),
                    (70.0, 80.0),
                    (10.0, 80.0),
                    (0.0, 70.0),
                ])
                .with_rotations_deg(vec![0.0, 90.0])
                .with_quantity(3),
            Geometry2D::new("shape8_gear")
                .with_polygon(vec![
                    (50.0, 5.0),
                    (65.0, 15.0),
                    (77.0, 18.0),
                    (80.0, 32.0),
                    (95.0, 50.0),
                    (80.0, 68.0),
                    (77.0, 82.0),
                    (65.0, 85.0),
                    (50.0, 95.0),
                    (35.0, 85.0),
                    (23.0, 82.0),
                    (20.0, 68.0),
                    (5.0, 50.0),
                    (20.0, 32.0),
                    (23.0, 18.0),
                    (35.0, 15.0),
                ])
                .with_rotations_deg(vec![0.0, 45.0, 90.0, 135.0, 180.0, 225.0, 270.0, 315.0])
                .with_quantity(13),
        ];

        // Build geometry lookup map
        let geom_map: HashMap<String, &Geometry2D> =
            shapes.iter().map(|g| (g.id().clone(), g)).collect();

        let boundary = Boundary2D::rectangle(500.0, 500.0);
        let strip_width = 500.0;

        // Test each strategy
        let strategies = vec![
            Strategy::BottomLeftFill,
            Strategy::NfpGuided,
            Strategy::GeneticAlgorithm,
            Strategy::Brkga,
            Strategy::SimulatedAnnealing,
            Strategy::Gdrr,
            Strategy::Alns,
        ];

        for strategy in strategies {
            println!("\n========== Testing {:?} ==========", strategy);

            let config = Config::default()
                .with_strategy(strategy)
                .with_time_limit(30000); // 30s max per strategy
            let nester = Nester2D::new(config);

            let result = match nester.solve_multi_strip(&shapes, &boundary) {
                Ok(r) => r,
                Err(e) => {
                    println!("  Strategy {:?} failed: {}", strategy, e);
                    continue;
                }
            };

            println!(
                "  Placed {} pieces across {} strips",
                result.placements.len(),
                result.boundaries_used
            );

            // Check 1: Boundary violations
            let mut boundary_violations = Vec::new();
            for p in &result.placements {
                // Find the base geometry ID (without instance suffix)
                let base_id = p.geometry_id.split('_').next().unwrap_or(&p.geometry_id);
                let full_id = if base_id.starts_with("shape") {
                    // Find matching geometry by checking all shape IDs
                    shapes
                        .iter()
                        .find(|g| p.geometry_id.starts_with(g.id()))
                        .map(|g| g.id().as_str())
                } else {
                    geom_map.get(&p.geometry_id).map(|g| g.id().as_str())
                };

                let geom = match full_id.and_then(|id| geom_map.get(id)) {
                    Some(g) => *g,
                    None => {
                        // Try to find by prefix match
                        match shapes.iter().find(|g| p.geometry_id.starts_with(g.id())) {
                            Some(g) => g,
                            None => {
                                println!(
                                    "  WARNING: Could not find geometry for {}",
                                    p.geometry_id
                                );
                                continue;
                            }
                        }
                    }
                };

                let origin_x = p.position[0];
                let origin_y = p.position[1];
                let rotation = p.rotation.first().copied().unwrap_or(0.0);
                let strip_idx = p.boundary_index;

                // Calculate local position within strip
                let local_x = origin_x - (strip_idx as f64 * strip_width);

                let (g_min, g_max) = geom.aabb_at_rotation(rotation);

                // Calculate actual bounds in local strip coordinates
                let local_min_x = local_x + g_min[0];
                let local_max_x = local_x + g_max[0];
                let local_min_y = origin_y + g_min[1];
                let local_max_y = origin_y + g_max[1];

                // Check boundary (with small tolerance)
                let tolerance = 0.1;
                if local_min_x < -tolerance
                    || local_max_x > 500.0 + tolerance
                    || local_min_y < -tolerance
                    || local_max_y > 500.0 + tolerance
                {
                    boundary_violations.push(format!(
                        "{} in strip {}: bounds ({:.1}, {:.1}) to ({:.1}, {:.1})",
                        p.geometry_id,
                        strip_idx,
                        local_min_x,
                        local_min_y,
                        local_max_x,
                        local_max_y
                    ));
                }
            }

            if !boundary_violations.is_empty() {
                println!("  BOUNDARY VIOLATIONS ({}):", boundary_violations.len());
                for v in &boundary_violations {
                    println!("    - {}", v);
                }
            }

            // Check 2: Overlaps (within same strip)
            let mut overlaps = Vec::new();
            let placements: Vec<_> = result.placements.iter().collect();

            for i in 0..placements.len() {
                for j in (i + 1)..placements.len() {
                    let p1 = placements[i];
                    let p2 = placements[j];

                    // Only check overlaps within the same strip
                    if p1.boundary_index != p2.boundary_index {
                        continue;
                    }

                    // Find geometries
                    let g1 = shapes.iter().find(|g| p1.geometry_id.starts_with(g.id()));
                    let g2 = shapes.iter().find(|g| p2.geometry_id.starts_with(g.id()));

                    let (g1, g2) = match (g1, g2) {
                        (Some(a), Some(b)) => (a, b),
                        _ => continue,
                    };

                    let strip_idx = p1.boundary_index;
                    let local_x1 = p1.position[0] - (strip_idx as f64 * strip_width);
                    let local_x2 = p2.position[0] - (strip_idx as f64 * strip_width);

                    let rot1 = p1.rotation.first().copied().unwrap_or(0.0);
                    let rot2 = p2.rotation.first().copied().unwrap_or(0.0);

                    let (g1_min, g1_max) = g1.aabb_at_rotation(rot1);
                    let (g2_min, g2_max) = g2.aabb_at_rotation(rot2);

                    let a_min = [local_x1 + g1_min[0], p1.position[1] + g1_min[1]];
                    let a_max = [local_x1 + g1_max[0], p1.position[1] + g1_max[1]];
                    let b_min = [local_x2 + g2_min[0], p2.position[1] + g2_min[1]];
                    let b_max = [local_x2 + g2_max[0], p2.position[1] + g2_max[1]];

                    if aabbs_overlap(a_min, a_max, b_min, b_max, 1.0) {
                        overlaps.push(format!(
                            "{} and {} in strip {}",
                            p1.geometry_id, p2.geometry_id, strip_idx
                        ));
                    }
                }
            }

            if !overlaps.is_empty() {
                println!("  OVERLAPS ({}):", overlaps.len());
                for o in overlaps.iter().take(10) {
                    println!("    - {}", o);
                }
                if overlaps.len() > 10 {
                    println!("    ... and {} more", overlaps.len() - 10);
                }
            }

            // Assert no boundary violations
            assert!(
                boundary_violations.is_empty(),
                "{:?}: Found {} boundary violations",
                strategy,
                boundary_violations.len()
            );

            println!("  ✓ All placements within boundary");
            println!("  ✓ No AABB overlaps detected");
        }
    }

    /// Regression guard for the multi-strip overflow distribution (ISSUE-20260621b).
    ///
    /// A single geometry with quantity 20 (100×100) cannot fit in one 300×300 sheet
    /// (9 per sheet). The prior id-level `retain` dropped a geometry from `remaining`
    /// as soon as *any* instance was placed, silently losing the rest. The fixed
    /// instance-level reduction must spread all 20 across sheets: 9 + 9 + 2 = 3 sheets,
    /// 0 unplaced, and every `(geometry_id, instance)` pair unique across all sheets.
    #[test]
    fn test_multi_strip_distributes_all_instances() {
        let geometries = vec![Geometry2D::rectangle("part", 100.0, 100.0).with_quantity(20)];
        let boundary = Boundary2D::rectangle(300.0, 300.0);
        let config = Config::default().with_strategy(Strategy::BottomLeftFill);
        let nester = Nester2D::new(config);

        let result = nester.solve_multi_strip(&geometries, &boundary).unwrap();

        // All 20 instances placed across exactly 3 sheets (9 + 9 + 2), none unplaced.
        assert_eq!(
            result.placements.len(),
            20,
            "all 20 instances must be placed"
        );
        assert_eq!(
            result.boundaries_used, 3,
            "20 of 100x100 in 300x300 => 3 sheets"
        );
        assert!(
            result.unplaced.is_empty(),
            "nothing should be unplaced, got {:?}",
            result.unplaced
        );
        // Instance-level request total is recorded (mirrors single-strip solve()).
        assert_eq!(result.total_requested, 20);

        // Every (geometry_id, instance) pair is globally unique — strips re-number
        // from 0 internally, so the multi-strip path must reassign global indices.
        let mut seen = std::collections::HashSet::new();
        for p in &result.placements {
            assert!(
                seen.insert((p.geometry_id.clone(), p.instance)),
                "duplicate (id, instance) = ({}, {}) across sheets",
                p.geometry_id,
                p.instance
            );
        }
        // Sheet indices are contiguous 0..3.
        let mut sheets: Vec<usize> = result.placements.iter().map(|p| p.boundary_index).collect();
        sheets.sort_unstable();
        sheets.dedup();
        assert_eq!(sheets, vec![0, 1, 2]);
    }

    /// When an item is genuinely too large for the sheet, the after-loop sweep must
    /// report it as unplaced (instance-level) rather than silently dropping it.
    #[test]
    fn test_multi_strip_oversized_reported_unplaced() {
        let geometries = vec![
            Geometry2D::rectangle("ok", 50.0, 50.0).with_quantity(2),
            Geometry2D::rectangle("toobig", 400.0, 400.0).with_quantity(3),
        ];
        let boundary = Boundary2D::rectangle(300.0, 300.0);
        let config = Config::default().with_strategy(Strategy::BottomLeftFill);
        let nester = Nester2D::new(config);

        let result = nester.solve_multi_strip(&geometries, &boundary).unwrap();

        assert_eq!(result.total_requested, 5, "2 + 3 instances requested");
        // The two 50x50 fit; the oversized geometry is reported unplaced (deduped id).
        assert_eq!(result.placements.len(), 2);
        assert!(
            result.unplaced.contains(&"toobig".to_string()),
            "oversized geometry must surface in unplaced, got {:?}",
            result.unplaced
        );
    }
}