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, run_ga_nesting_with_progress};
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, 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::{Placement, Result, SolveResult};

use crate::placement_utils::{expand_nfp, shrink_ifp};
use std::sync::atomic::{AtomicBool, Ordering};
use std::sync::Arc;
use u_nesting_core::timing::Timer;

/// 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, see `logs/ROADMAP.md` §12), 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 longer (open-roll) axis.
    let (b_min, b_max) = boundary.aabb();
    let bound_w = b_max[0] - b_min[0];
    let bound_h = b_max[1] - b_min[1];
    let strip_length = if bound_h >= bound_w {
        max_y - min_y
    } else {
        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,
}

impl Nester2D {
    /// 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())
    }

    /// 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>> {
        let start = Timer::now();
        let mut result = SolveResult::new();
        let mut placements = Vec::new();

        // Get boundary dimensions
        let (b_min, b_max) = boundary.aabb();
        let margin = self.config.margin;
        let spacing = self.config.spacing;

        let bound_min_x = b_min[0] + margin;
        let bound_min_y = b_min[1] + margin;
        let bound_max_x = b_max[0] - margin;
        let bound_max_y = b_max[1] - margin;

        let strip_width = bound_max_x - bound_min_x;
        let strip_height = bound_max_y - bound_min_y;

        // Simple row-based placement with rotation optimization
        let mut current_x = bound_min_x;
        let mut current_y = bound_min_y;
        let mut row_height = 0.0_f64;

        let mut total_placed_area = 0.0;

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

            // Get allowed rotation angles (default to 0 if none specified)
            let rotations = geom.rotations();
            let rotation_angles: Vec<f64> = if rotations.is_empty() {
                vec![0.0]
            } else {
                rotations
            };

            for instance in 0..geom.quantity() {
                if self.cancelled.load(Ordering::Relaxed) {
                    result.computation_time_ms = start.elapsed_ms();
                    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;
                    return Ok(result);
                }

                // Find the best rotation for current position
                let mut best_fit: Option<(f64, f64, f64, f64, f64, [f64; 2])> = None; // (rotation, width, height, x, y, g_min)

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

                    // Skip if geometry doesn't fit in boundary at all
                    if g_width > strip_width || g_height > strip_height {
                        continue;
                    }

                    // Calculate placement position for this rotation
                    let mut place_x = current_x;
                    let mut place_y = current_y;

                    // Check if piece fits in remaining row space
                    if place_x + g_width > bound_max_x {
                        // Would need to move to next row
                        place_x = bound_min_x;
                        place_y += row_height + spacing;
                    }

                    // Check if piece fits in boundary height
                    if place_y + g_height > bound_max_y {
                        continue; // This rotation doesn't fit
                    }

                    // Calculate score: prefer rotations that minimize wasted space
                    // Score = row advancement (lower is better)
                    let score = if place_x == bound_min_x && place_y > current_y {
                        // New row: score is based on new Y position
                        place_y - bound_min_y + g_height
                    } else {
                        // Same row: score is based on strip length advancement
                        place_x - bound_min_x + g_width
                    };

                    let is_better = match &best_fit {
                        None => true,
                        Some((_, _, _, _, _, _)) => {
                            // Prefer placements that don't start new rows
                            let best_score = if let Some((_, _, _, bx, by, _)) = best_fit {
                                if bx == bound_min_x && by > current_y {
                                    by - bound_min_y + g_height
                                } else {
                                    bx - bound_min_x + g_width
                                }
                            } else {
                                f64::INFINITY
                            };
                            score < best_score - 1e-6
                        }
                    };

                    if is_better {
                        best_fit = Some((rotation, g_width, g_height, place_x, place_y, g_min));
                    }
                }

                // Place the geometry with the best rotation
                if let Some((rotation, g_width, g_height, place_x, place_y, g_min)) = best_fit {
                    // Update row tracking if we moved to a new row
                    if place_x == bound_min_x && place_y > current_y {
                        row_height = 0.0;
                    }

                    // Compute origin position from AABB position
                    let origin_x = place_x - g_min[0];
                    let origin_y = place_y - g_min[1];

                    // Clamp to ensure geometry stays within boundary
                    let geom_aabb = geom.aabb_at_rotation(rotation);
                    let boundary_aabb = (b_min, b_max);

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

                        placements.push(placement);
                        total_placed_area += geom.measure();

                        // Update position for next piece
                        // Use actual clamped AABB position, not original place_x/place_y
                        let actual_place_x = clamped_x + g_min[0];
                        let actual_place_y = clamped_y + g_min[1];
                        current_x = actual_place_x + g_width + spacing;
                        current_y = actual_place_y;
                        row_height = row_height.max(g_height);
                    }
                } else {
                    // Can't place this piece with any rotation
                    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();

        Ok(result)
    }

    /// 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>> {
        let start = Timer::now();
        let mut result = SolveResult::new();
        let mut placements = Vec::new();
        let mut placed_geometries: Vec<PlacedGeometry> = Vec::new();

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

        // Get boundary polygon with margin applied
        let boundary_polygon = self.get_boundary_polygon_with_margin(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
            let rotations = geom.rotations();
            let rotation_angles: Vec<f64> = if rotations.is_empty() {
                vec![0.0]
            } else {
                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();
                    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;
                    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,
                            margin,
                            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,
                                    )
                                }) {
                                    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);
                            let expanded = self.expand_nfp(&translated_nfp, spacing);
                            nfps.push(expanded);
                        }

                        // Shrink IFP by spacing from boundary
                        let ifp_shrunk = self.shrink_ifp(&ifp, spacing);

                        // Find the optimal valid placement (minimize X for shorter strip)
                        let nfp_refs: Vec<&Nfp> = nfps.iter().collect();
                        if let Some((x, y)) =
                            find_bottom_left_placement(&ifp_shrunk, &nfp_refs, sample_step)
                        {
                            // Compare with current best: prefer smaller X (shorter strip), then smaller Y
                            let is_better = match best_placement {
                                None => true,
                                Some((best_x, best_y, _, _)) => {
                                    x < best_x - 1e-6 || (x < best_x + 1e-6 && y < 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();
                    } 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();

        Ok(result)
    }

    /// Gets the boundary polygon with margin applied.
    fn get_boundary_polygon_with_margin(
        &self,
        boundary: &Boundary2D,
        margin: f64,
    ) -> Vec<(f64, f64)> {
        let (b_min, b_max) = boundary.aabb();

        // Create a rectangular boundary polygon with margin
        vec![
            (b_min[0] + margin, b_min[1] + margin),
            (b_max[0] - margin, b_min[1] + margin),
            (b_max[0] - margin, b_max[1] - margin),
            (b_min[0] + margin, b_max[1] - margin),
        ]
    }

    /// 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)
    }

    /// Expands an NFP by the given spacing amount.
    fn expand_nfp(&self, nfp: &Nfp, spacing: f64) -> Nfp {
        expand_nfp(nfp, spacing)
    }

    /// Shrinks an IFP by the given spacing amount.
    fn shrink_ifp(&self, ifp: &Nfp, spacing: f64) -> Nfp {
        shrink_ifp(ifp, spacing)
    }

    /// Returns whichever of `meta` (a metaheuristic result) and a fresh
    /// Bottom-Left-Fill solve packs better.
    ///
    /// The stochastic strategies (GA/BRKGA/SA) can converge to 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_blf(
        &self,
        meta: SolveResult<f64>,
        geometries: &[Geometry2D],
        boundary: &Boundary2D,
    ) -> SolveResult<f64> {
        let blf = match self.bottom_left_fill(geometries, boundary) {
            Ok(b) => b,
            Err(_) => return meta,
        };
        let (meta_placed, meta_len) = solution_quality(&meta, geometries, boundary);
        let (blf_placed, blf_len) = solution_quality(&blf, geometries, boundary);
        // BLF wins if it places strictly more pieces, or ties on count while
        // consuming a strictly shorter strip length.
        if blf_placed > meta_placed || (blf_placed == meta_placed && blf_len < meta_len - 1e-6) {
            // The greedy layout is better, so its placements are returned — but the
            // metaheuristic *did* run. Preserve its search diagnostics (strategy
            // label, generation/fitness history) so a caller inspecting the result
            // still sees which strategy executed and how it converged, rather than a
            // bare BLF. Only the placements are floored, not the provenance.
            let mut floored = blf;
            floored.strategy = meta.strategy;
            floored.generations = meta.generations;
            floored.best_fitness = meta.best_fitness;
            floored.fitness_history = meta.fitness_history;
            floored.target_reached = meta.target_reached;
            floored
        } else {
            meta
        }
    }

    /// Genetic Algorithm based nesting optimization.
    ///
    /// Uses GA to optimize placement order and rotations, with NFP-guided
    /// decoding for collision-free placements.
    fn genetic_algorithm(
        &self,
        geometries: &[Geometry2D],
        boundary: &Boundary2D,
    ) -> Result<SolveResult<f64>> {
        // Configure GA with time limit for multi-strip scenarios
        let time_limit_ms = if self.config.time_limit_ms > 0 {
            // Use 1/4 of total time limit per strip to allow for multiple strips
            // Budget a quarter of the total per strip (assuming up to ~4 strips), but
            // never exceed the user's total limit: a single-strip solve must honor an
            // explicit short budget instead of being floored up to 5s.
            (self.config.time_limit_ms / 4)
                .max(5000)
                .min(self.config.time_limit_ms)
        } else {
            15000 // 15 seconds default per strip
        };

        let ga_config = GaConfig::default()
            .with_population_size(self.config.population_size.min(30)) // Limit population
            .with_max_generations(self.config.max_generations.min(50)) // Limit generations
            .with_crossover_rate(self.config.crossover_rate)
            .with_mutation_rate(self.config.mutation_rate)
            .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(),
        );

        Ok(self.not_worse_than_blf(result, geometries, boundary))
    }

    /// 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>> {
        // Configure BRKGA with time limit for multi-strip scenarios
        let time_limit_ms = if self.config.time_limit_ms > 0 {
            // Use 1/4 of total time limit per strip to allow for multiple strips
            // Budget a quarter of the total per strip (assuming up to ~4 strips), but
            // never exceed the user's total limit: a single-strip solve must honor an
            // explicit short budget instead of being floored up to 5s.
            (self.config.time_limit_ms / 4)
                .max(5000)
                .min(self.config.time_limit_ms)
        } else {
            15000 // 15 seconds default per strip
        };

        let brkga_config = BrkgaConfig::default()
            // Honor the caller's population_size (was hardcoded to 30, silently
            // ignoring config). Capped at 30 as a performance ceiling, matching
            // the GA path — each decode is an O(N²) NFP evaluation.
            .with_population_size(self.config.population_size.min(30))
            .with_max_generations(50) // Fewer generations
            .with_elite_fraction(0.2)
            .with_mutant_fraction(0.15)
            .with_elite_bias(0.7)
            .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(),
        );

        Ok(self.not_worse_than_blf(result, geometries, boundary))
    }

    /// 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>> {
        // Configure SA with faster defaults for multi-strip scenarios
        // Note: Each decode() call is O(N²) NFP computations, so we need fewer iterations
        let time_limit_ms = if self.config.time_limit_ms > 0 {
            // Use 1/4 of total time limit per strip to allow for multiple strips
            // Budget a quarter of the total per strip (assuming up to ~4 strips), but
            // never exceed the user's total limit: a single-strip solve must honor an
            // explicit short budget instead of being floored up to 5s.
            (self.config.time_limit_ms / 4)
                .max(5000)
                .min(self.config.time_limit_ms)
        } else {
            10000 // 10 seconds default per strip
        };

        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 result = run_sa_nesting(
            geometries,
            boundary,
            &self.config,
            sa_config,
            self.cancelled.clone(),
        );

        Ok(self.not_worse_than_blf(result, geometries, boundary))
    }

    /// Goal-Driven Ruin and Recreate (GDRR) optimization.
    fn gdrr(&self, geometries: &[Geometry2D], boundary: &Boundary2D) -> Result<SolveResult<f64>> {
        // Configure GDRR with faster defaults for multi-strip scenarios
        // Use user's time limit, default to 10s per strip if not specified
        let time_limit = if self.config.time_limit_ms > 0 {
            // Use 1/4 of total time limit per strip to allow for multiple strips
            // Budget a quarter of the total per strip (assuming up to ~4 strips), but
            // never exceed the user's total limit: a single-strip solve must honor an
            // explicit short budget instead of being floored up to 5s.
            (self.config.time_limit_ms / 4)
                .max(5000)
                .min(self.config.time_limit_ms)
        } else {
            10000 // 10 seconds default per strip
        };
        let gdrr_config = GdrrConfig::default()
            .with_max_iterations(1000) // Reduced from 5000 for faster execution
            .with_time_limit_ms(time_limit)
            .with_ruin_ratio(0.1, 0.3) // Smaller ruin ratio for faster convergence
            .with_lahc_list_length(30); // Smaller list for faster convergence

        let result = run_gdrr_nesting(
            geometries,
            boundary,
            &self.config,
            &gdrr_config,
            self.cancelled.clone(),
        );

        Ok(result)
    }

    /// Adaptive Large Neighborhood Search (ALNS) optimization.
    fn alns(&self, geometries: &[Geometry2D], boundary: &Boundary2D) -> Result<SolveResult<f64>> {
        // Configure ALNS with faster defaults for multi-strip scenarios
        // Use user's time limit, default to 10s per strip if not specified
        let time_limit = if self.config.time_limit_ms > 0 {
            // Use 1/4 of total time limit per strip to allow for multiple strips
            // Budget a quarter of the total per strip (assuming up to ~4 strips), but
            // never exceed the user's total limit: a single-strip solve must honor an
            // explicit short budget instead of being floored up to 5s.
            (self.config.time_limit_ms / 4)
                .max(5000)
                .min(self.config.time_limit_ms)
        } else {
            10000 // 10 seconds default per strip
        };
        let alns_config = AlnsConfig::default()
            .with_max_iterations(1000) // Reduced from 5000 for faster execution
            .with_time_limit_ms(time_limit)
            .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

        let result = run_alns_nesting(
            geometries,
            boundary,
            &self.config,
            &alns_config,
            self.cancelled.clone(),
        );

        Ok(result)
    }

    /// 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>> {
        let exact_config = ExactConfig::default()
            .with_time_limit_ms(self.config.time_limit_ms.max(60000))
            .with_max_items(15)
            .with_rotation_steps(4)
            .with_grid_step(1.0);

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

        Ok(result)
    }

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

        // If small enough, try exact
        if total_instances <= 15 {
            let exact_config = ExactConfig::default()
                .with_time_limit_ms((self.config.time_limit_ms / 2).max(30000))
                .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)
        self.alns(geometries, boundary)
    }

    /// Bottom-Left Fill with progress callback.
    ///
    /// Same AABB-only reasoning as `bottom_left_fill` — mirror candidates are
    /// intentionally not enumerated here, see that function's doc comment.
    fn bottom_left_fill_with_progress(
        &self,
        geometries: &[Geometry2D],
        boundary: &Boundary2D,
        callback: &ProgressCallback,
    ) -> Result<SolveResult<f64>> {
        let start = Timer::now();
        let mut result = SolveResult::new();
        let mut placements = Vec::new();

        // Get boundary dimensions
        let (b_min, b_max) = boundary.aabb();
        let margin = self.config.margin;
        let spacing = self.config.spacing;

        let bound_min_x = b_min[0] + margin;
        let bound_min_y = b_min[1] + margin;
        let bound_max_x = b_max[0] - margin;
        let bound_max_y = b_max[1] - margin;

        let strip_width = bound_max_x - bound_min_x;
        let strip_height = bound_max_y - bound_min_y;

        let mut current_x = bound_min_x;
        let mut current_y = bound_min_y;
        let mut row_height = 0.0_f64;
        let mut total_placed_area = 0.0;

        // Count total pieces for progress
        let total_pieces: usize = geometries.iter().map(|g| g.quantity()).sum();
        let mut placed_count = 0usize;

        // Initial progress callback
        callback(
            ProgressInfo::new()
                .with_phase("BLF Placement")
                .with_items(0, total_pieces)
                .with_elapsed(0),
        );

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

            let rotations = geom.rotations();
            let rotation_angles: Vec<f64> = if rotations.is_empty() {
                vec![0.0]
            } else {
                rotations
            };

            for instance in 0..geom.quantity() {
                if self.cancelled.load(Ordering::Relaxed) {
                    result.computation_time_ms = start.elapsed_ms();
                    callback(
                        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;
                    callback(
                        ProgressInfo::new()
                            .with_phase("Time Limit Reached")
                            .with_items(placed_count, total_pieces)
                            .with_elapsed(result.computation_time_ms)
                            .finished(),
                    );
                    return Ok(result);
                }

                let mut best_fit: Option<(f64, f64, f64, f64, f64, [f64; 2])> = None;

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

                    if g_width > strip_width || g_height > strip_height {
                        continue;
                    }

                    let mut place_x = current_x;
                    let mut place_y = current_y;

                    if place_x + g_width > bound_max_x {
                        place_x = bound_min_x;
                        place_y += row_height + spacing;
                    }

                    if place_y + g_height > bound_max_y {
                        continue;
                    }

                    let score = if place_x == bound_min_x && place_y > current_y {
                        place_y - bound_min_y + g_height
                    } else {
                        place_x - bound_min_x + g_width
                    };

                    let is_better = match &best_fit {
                        None => true,
                        Some((_, _, _, bx, by, _)) => {
                            let best_score = if *bx == bound_min_x && *by > current_y {
                                by - bound_min_y
                            } else {
                                bx - bound_min_x
                            };
                            score < best_score - 1e-6
                        }
                    };

                    if is_better {
                        best_fit = Some((rotation, g_width, g_height, place_x, place_y, g_min));
                    }
                }

                if let Some((rotation, g_width, g_height, place_x, place_y, g_min)) = best_fit {
                    if place_x == bound_min_x && place_y > current_y {
                        row_height = 0.0;
                    }

                    // Compute origin position from AABB position
                    let origin_x = place_x - g_min[0];
                    let origin_y = place_y - g_min[1];

                    // Clamp to ensure geometry stays within boundary
                    let geom_aabb = geom.aabb_at_rotation(rotation);
                    let boundary_aabb = (b_min, b_max);

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

                        placements.push(placement);
                        total_placed_area += geom.measure();
                        placed_count += 1;

                        current_x = place_x + g_width + spacing;
                        current_y = place_y;
                        row_height = row_height.max(g_height);

                        // Progress callback every piece
                        callback(
                            ProgressInfo::new()
                                .with_phase("BLF Placement")
                                .with_items(placed_count, total_pieces)
                                .with_utilization(total_placed_area / boundary.measure())
                                .with_elapsed(start.elapsed_ms()),
                        );
                    } else {
                        result.unplaced.push(geom.id().clone());
                    }
                } else {
                    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();

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

        Ok(result)
    }

    /// NFP-guided BLF with progress callback.
    fn nfp_guided_blf_with_progress(
        &self,
        geometries: &[Geometry2D],
        boundary: &Boundary2D,
        callback: &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 margin = self.config.margin;
        let spacing = self.config.spacing;
        let boundary_polygon = self.get_boundary_polygon_with_margin(boundary, margin);

        let mut total_placed_area = 0.0;
        let sample_step = self.compute_sample_step(geometries);

        // Count total pieces for progress
        let total_pieces: usize = geometries.iter().map(|g| g.quantity()).sum();
        let mut placed_count = 0usize;

        // Initial progress callback
        callback(
            ProgressInfo::new()
                .with_phase("NFP Placement")
                .with_items(0, total_pieces)
                .with_elapsed(0),
        );

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

            let rotations = geom.rotations();
            let rotation_angles: Vec<f64> = if rotations.is_empty() {
                vec![0.0]
            } else {
                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();
                    callback(
                        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;
                    callback(
                        ProgressInfo::new()
                            .with_phase("Time Limit Reached")
                            .with_items(placed_count, total_pieces)
                            .with_elapsed(result.computation_time_ms)
                            .finished(),
                    );
                    return Ok(result);
                }

                let mut best_placement: Option<(f64, f64, f64, bool)> = None;

                for &rotation in &rotation_angles {
                    for &mirror in mirror_candidates {
                        let ifp = match compute_ifp_with_margin_and_mirror(
                            &boundary_polygon,
                            geom,
                            rotation,
                            margin,
                            mirror,
                        ) {
                            Ok(ifp) => ifp,
                            Err(_) => continue,
                        };

                        if ifp.is_empty() {
                            continue;
                        }

                        let mut nfps: Vec<Nfp> = Vec::new();
                        for placed in &placed_geometries {
                            // Use cache for NFP computation
                            let cache_key = (
                                placed.geometry.id().as_str(),
                                geom.id().as_str(),
                                rotation - placed.rotation,
                                placed.mirrored,
                                mirror,
                            );

                            // Compute NFP at origin and cache it (with relative rotation)
                            // 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,
                                    )
                                }) {
                                    Ok(nfp) => nfp,
                                    Err(_) => continue,
                                };

                            // Transform NFP: first rotate by placed.rotation, then translate
                            let rotated_nfp = rotate_nfp(&nfp_at_origin, placed.rotation);
                            let translated_nfp = translate_nfp(&rotated_nfp, placed.position);
                            let expanded = self.expand_nfp(&translated_nfp, spacing);
                            nfps.push(expanded);
                        }

                        let ifp_shrunk = self.shrink_ifp(&ifp, spacing);
                        let nfp_refs: Vec<&Nfp> = nfps.iter().collect();

                        if let Some((x, y)) =
                            find_bottom_left_placement(&ifp_shrunk, &nfp_refs, sample_step)
                        {
                            let is_better = match best_placement {
                                None => true,
                                Some((best_x, best_y, _, _)) => {
                                    x < best_x - 1e-6 || (x < best_x + 1e-6 && y < best_y - 1e-6)
                                }
                            };
                            if is_better {
                                best_placement = Some((x, y, rotation, mirror));
                            }
                        }
                    }
                }

                if let Some((x, y, rotation, mirror)) = best_placement {
                    // Clamp to ensure geometry stays within boundary
                    // (mirror-aware — see the same fix in `nfp_guided_blf`).
                    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;

                        // Progress callback every piece
                        callback(
                            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 {
                        result.unplaced.push(geom.id().clone());
                    }
                } else {
                    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();

        // Final progress callback
        callback(
            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.
    fn validate_geometries(&self, geometries: &[Geometry2D]) -> Result<()> {
        use u_nesting_core::geometry::Geometry;
        for geom in geometries {
            geom.validate()?;
        }
        Ok(())
    }

    /// 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 = match self.config.strategy {
                Strategy::BottomLeftFill => self.bottom_left_fill(&remaining_geometries, boundary),
                Strategy::NfpGuided => self.nfp_guided_blf(&remaining_geometries, boundary),
                Strategy::GeneticAlgorithm => {
                    self.genetic_algorithm(&remaining_geometries, boundary)
                }
                Strategy::Brkga => self.brkga(&remaining_geometries, boundary),
                Strategy::SimulatedAnnealing => {
                    self.simulated_annealing(&remaining_geometries, boundary)
                }
                Strategy::Gdrr => self.gdrr(&remaining_geometries, boundary),
                Strategy::Alns => self.alns(&remaining_geometries, boundary),
                #[cfg(feature = "milp")]
                Strategy::MilpExact => self.milp_exact(&remaining_geometries, boundary),
                #[cfg(feature = "milp")]
                Strategy::HybridExact => self.hybrid_exact(&remaining_geometries, boundary),
                _ => self.nfp_guided_blf(&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);

            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 Solver for Nester2D {
    type Geometry = Geometry2D;
    type Boundary = Boundary2D;
    type Scalar = f64;

    fn solve(
        &self,
        geometries: &[Self::Geometry],
        boundary: &Self::Boundary,
    ) -> Result<SolveResult<f64>> {
        boundary.validate()?;
        self.validate_geometries(geometries)?;

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

        let initial_result = match self.config.strategy {
            Strategy::BottomLeftFill => self.bottom_left_fill(geometries, boundary),
            Strategy::NfpGuided => self.nfp_guided_blf(geometries, boundary),
            Strategy::GeneticAlgorithm => self.genetic_algorithm(geometries, boundary),
            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),
            _ => {
                // Fall back to NFP-guided BLF for other strategies
                log::warn!(
                    "Strategy {:?} not yet implemented, using NfpGuided",
                    self.config.strategy
                );
                self.nfp_guided_blf(geometries, boundary)
            }
        }?;

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

        // 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)
    }

    fn solve_with_progress(
        &self,
        geometries: &[Self::Geometry],
        boundary: &Self::Boundary,
        callback: ProgressCallback,
    ) -> Result<SolveResult<f64>> {
        boundary.validate()?;
        self.validate_geometries(geometries)?;

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

        let initial_result = match self.config.strategy {
            Strategy::BottomLeftFill => {
                self.bottom_left_fill_with_progress(geometries, boundary, &callback)?
            }
            Strategy::NfpGuided => {
                self.nfp_guided_blf_with_progress(geometries, boundary, &callback)?
            }
            Strategy::GeneticAlgorithm => {
                // Cap population/generations to match the non-progress
                // `genetic_algorithm` path. Without this the callback path ran the
                // full default 500 generations × 100 population (vs 50 × 30),
                // making a progress-driven solve dramatically slower for no quality
                // gain and letting it overrun a modest time budget.
                let mut ga_config = GaConfig::default()
                    .with_population_size(self.config.population_size.min(30))
                    .with_max_generations(self.config.max_generations.min(50))
                    .with_crossover_rate(self.config.crossover_rate)
                    .with_mutation_rate(self.config.mutation_rate);

                // Apply time limit if specified
                if self.config.time_limit_ms > 0 {
                    ga_config = ga_config.with_time_limit(std::time::Duration::from_millis(
                        self.config.time_limit_ms,
                    ));
                }

                let ga_result = run_ga_nesting_with_progress(
                    geometries,
                    boundary,
                    &self.config,
                    ga_config,
                    self.cancelled.clone(),
                    callback,
                );
                // Same BLF floor as the non-progress path: never return a
                // layout worse than deterministic bottom-left-fill. The
                // callback-driven entry points (FFI `solve_2d_with_callback`,
                // WASM/demo) route through here, so the guard must apply here too.
                self.not_worse_than_blf(ga_result, geometries, boundary)
            }
            // For other strategies, use basic progress reporting
            _ => {
                log::warn!(
                    "Strategy {:?} not yet implemented, using NfpGuided",
                    self.config.strategy
                );
                self.nfp_guided_blf_with_progress(geometries, boundary, &callback)?
            }
        };

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

        // 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)
    }

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

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

    #[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 (`expand_nfp`/`shrink_ifp` 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 ifp1_shrunk = shrink_ifp(&ifp1, spacing);
        let (x1, y1) =
            find_bottom_left_placement(&ifp1_shrunk, &[], 1.0).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 ifp2_shrunk = shrink_ifp(&ifp2, spacing);
        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 = expand_nfp(&translated_nfp, spacing);
        let (x2, y2) = find_bottom_left_placement(&ifp2_shrunk, &[&expanded_nfp], 1.0)
            .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
        );
    }
}