u_nesting_core/

sa.rs

//! Simulated Annealing framework for optimization.
//!
//! # Architecture
//!
//! This module maintains its own SA loop rather than delegating to u-metaheur's
//! SA runner. u-metaheur uses `cost(&self, &Solution) -> f64` (immutable),
//! while u-nesting uses `evaluate(&self, &mut Solution)` (mutable) for
//! consistency with the GA/BRKGA evaluation pattern.
//!
//! Additionally, this module provides features not available in u-metaheur SA:
//! - `NeighborhoodOperator` enum (5 operators) with `available_operators()`
//! - `PermutationSolution` with built-in swap/relocate/inversion/rotation/chain
//! - Reheating with configurable threshold and factor
//! - Adaptive cooling schedule
//!
//! The rand 0.9 API is shared with u-metaheur for ecosystem compatibility.

use rand::prelude::*;
use std::sync::atomic::{AtomicBool, Ordering};
use std::sync::Arc;
use std::time::Duration;
use u_metaheur::TemperatureScale;

use crate::error::{check_above, check_at_least, check_range, Result};
use crate::result::StopReason;
use crate::timing::{expired, Timer};

#[cfg(feature = "serde")]
use serde::{Deserialize, Serialize};

/// Cooling schedule types for Simulated Annealing.
#[derive(Debug, Clone, Copy, Default)]
#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
pub enum CoolingSchedule {
    /// Geometric cooling: T_new = T * alpha (alpha typically 0.95-0.99).
    #[default]
    Geometric,
    /// Linear cooling: T_new = T - delta.
    Linear,
    /// Adaptive cooling: adjusts based on acceptance rate.
    Adaptive,
    /// Lundy-Mees: T_new = T / (1 + beta * T).
    LundyMees,
}

/// Configuration for Simulated Annealing.
#[derive(Debug, Clone)]
#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
pub struct SaConfig {
    /// What the temperatures are measured in. The default measures the
    /// objective's own scale first (see [`TemperatureScale::Calibrated`]), so one
    /// configuration serves any unit of fitness;
    /// [`with_initial_temp`](Self::with_initial_temp) sets an absolute
    /// temperature, [`with_relative_temperature`](Self::with_relative_temperature)
    /// a calibrated one.
    #[cfg_attr(feature = "serde", serde(default))]
    pub temperature_scale: TemperatureScale,
    /// Initial temperature, in units of `temperature_scale`.
    pub initial_temp: f64,
    /// Final (minimum) temperature, in the same unit.
    pub final_temp: f64,
    /// Cooling rate (alpha for Geometric, delta for Linear, beta for LundyMees).
    pub cooling_rate: f64,
    /// Number of iterations at each temperature level.
    pub iterations_per_temp: usize,
    /// Maximum total iterations (None = temperature-based stopping only).
    pub max_iterations: Option<u64>,
    /// Cooling schedule type.
    pub cooling_schedule: CoolingSchedule,
    /// Maximum time limit (None = unlimited).
    pub time_limit: Option<Duration>,
    /// Target fitness to stop early (None = run until temperature limit).
    pub target_fitness: Option<f64>,
    /// Enable reheating when stagnation detected.
    pub enable_reheating: bool,
    /// Stagnation threshold for reheating.
    pub reheat_threshold: u64,
    /// Reheat factor (multiplier for current temperature).
    pub reheat_factor: f64,
}

impl Default for SaConfig {
    fn default() -> Self {
        Self {
            temperature_scale: TemperatureScale::default(),
            initial_temp: 1.0,
            final_temp: 0.001,
            cooling_rate: 0.95,
            iterations_per_temp: 100,
            max_iterations: Some(100_000),
            cooling_schedule: CoolingSchedule::Geometric,
            time_limit: None,
            target_fitness: None,
            enable_reheating: false,
            reheat_threshold: 1000,
            reheat_factor: 2.0,
        }
    }
}

impl SaConfig {
    /// Creates a new configuration with default values.
    pub fn new() -> Self {
        Self::default()
    }

    /// Sets the initial temperature **in the unit of the fitness**
    /// ([`TemperatureScale::Absolute`]).
    pub fn with_initial_temp(mut self, temp: f64) -> Self {
        self.temperature_scale = TemperatureScale::Absolute;
        self.initial_temp = temp;
        self
    }

    /// Sets the start and the floor as multiples of the fitness scale measured
    /// over `samples` neighbours of the initial solution
    /// ([`TemperatureScale::Calibrated`]). Each sample is a full evaluation, so a
    /// problem whose evaluation is expensive wants few (the measurement stops at
    /// the time limit).
    pub fn with_relative_temperature(
        mut self,
        initial: f64,
        final_temp: f64,
        samples: usize,
    ) -> Self {
        self.temperature_scale = TemperatureScale::Calibrated { samples };
        self.initial_temp = initial;
        self.final_temp = final_temp;
        self
    }

    /// Sets the final temperature.
    pub fn with_final_temp(mut self, temp: f64) -> Self {
        self.final_temp = temp;
        self
    }

    /// Sets the cooling rate.
    pub fn with_cooling_rate(mut self, rate: f64) -> Self {
        self.cooling_rate = rate;
        self
    }

    /// Checks every field against its range; the builders store what they are
    /// given.
    ///
    /// # Errors
    /// `ConfigError` naming the first field out of range: `cooling_rate`
    /// outside `[0.001, 0.9999]`, temperatures that are not finite and
    /// positive, `final_temp` above `initial_temp`, `iterations_per_temp`
    /// below 1, `reheat_factor` below 1.
    pub fn validate(&self) -> Result<()> {
        check_range("cooling_rate", self.cooling_rate, 0.001, 0.9999)?;
        check_above("initial_temp", self.initial_temp, 0.0)?;
        check_above("final_temp", self.final_temp, 0.0)?;
        check_range("final_temp", self.final_temp, 0.0, self.initial_temp)?;
        check_at_least("iterations_per_temp", self.iterations_per_temp as f64, 1.0)?;
        check_at_least("reheat_factor", self.reheat_factor, 1.0)?;
        if let TemperatureScale::Calibrated { samples } = self.temperature_scale {
            check_at_least("temperature_scale.samples", samples as f64, 1.0)?;
        }
        Ok(())
    }

    /// Sets the iterations per temperature level.
    pub fn with_iterations_per_temp(mut self, iterations: usize) -> Self {
        self.iterations_per_temp = iterations;
        self
    }

    /// Sets the maximum iterations.
    pub fn with_max_iterations(mut self, iterations: u64) -> Self {
        self.max_iterations = Some(iterations);
        self
    }

    /// Sets the cooling schedule.
    pub fn with_cooling_schedule(mut self, schedule: CoolingSchedule) -> Self {
        self.cooling_schedule = schedule;
        self
    }

    /// Sets the time limit.
    pub fn with_time_limit(mut self, duration: Duration) -> Self {
        self.time_limit = Some(duration);
        self
    }

    /// Sets the target fitness.
    pub fn with_target_fitness(mut self, fitness: f64) -> Self {
        self.target_fitness = Some(fitness);
        self
    }

    /// Enables reheating.
    pub fn with_reheating(mut self, threshold: u64, factor: f64) -> Self {
        self.enable_reheating = true;
        self.reheat_threshold = threshold;
        self.reheat_factor = factor;
        self
    }
}

/// Trait for solutions in Simulated Annealing.
pub trait SaSolution: Clone + Send + Sync {
    /// Returns the objective value (fitness) of this solution.
    /// Higher values are better (maximization).
    fn objective(&self) -> f64;

    /// Sets the objective value.
    fn set_objective(&mut self, value: f64);
}

/// Neighborhood operator types.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum NeighborhoodOperator {
    /// Swap two elements.
    Swap,
    /// Relocate an element to a new position.
    Relocate,
    /// Reverse a segment (2-opt style).
    Inversion,
    /// Rotate/change orientation of an element.
    Rotation,
    /// Chain swap (3-opt style).
    Chain,
    /// Flip an element's mirror flag (`allow_flip` support, 2D nesting only
    /// — a solution's `mirrors` vector is harmless-if-unused by problems
    /// with no mirror concept, same as an all-`false` `rotations` vector is
    /// for a problem with a single rotation option).
    MirrorFlip,
}

/// Trait for problem-specific SA operations.
pub trait SaProblem: Send + Sync {
    /// The solution type for this problem.
    type Solution: SaSolution;

    /// Creates an initial solution.
    fn initial_solution<R: Rng>(&self, rng: &mut R) -> Self::Solution;

    /// Generates a neighbor solution using the specified operator.
    fn neighbor<R: Rng>(
        &self,
        solution: &Self::Solution,
        operator: NeighborhoodOperator,
        rng: &mut R,
    ) -> Self::Solution;

    /// Evaluates the objective of a solution.
    fn evaluate(&self, solution: &mut Self::Solution);

    /// Returns available neighborhood operators for this problem.
    fn available_operators(&self) -> Vec<NeighborhoodOperator> {
        vec![
            NeighborhoodOperator::Swap,
            NeighborhoodOperator::Relocate,
            NeighborhoodOperator::Inversion,
        ]
    }

    /// Called after each temperature level (for progress reporting).
    fn on_temperature_change(
        &self,
        _temperature: f64,
        _iteration: u64,
        _best: &Self::Solution,
        _current: &Self::Solution,
    ) {
        // Default: do nothing
    }
}

/// Progress information during SA execution.
#[derive(Debug, Clone)]
pub struct SaProgress {
    /// Current temperature.
    pub temperature: f64,
    /// Current iteration number.
    pub iteration: u64,
    /// Best fitness so far.
    pub best_fitness: f64,
    /// Current fitness.
    pub current_fitness: f64,
    /// Acceptance rate at current temperature.
    pub acceptance_rate: f64,
    /// Elapsed time since start.
    pub elapsed: Duration,
    /// Whether the algorithm is still running.
    pub running: bool,
}

/// Result of a SA run.
#[derive(Debug, Clone)]
pub struct SaResult<S: SaSolution> {
    /// The best solution found.
    pub best: S,
    /// Final temperature reached, in the unit of the fitness.
    pub final_temperature: f64,
    /// The factor the configured temperatures were multiplied by: the measured
    /// fitness scale of a calibrated run, `1` for an absolute one.
    pub temperature_unit: f64,
    /// Total iterations performed.
    pub iterations: u64,
    /// Total elapsed time.
    pub elapsed: Duration,
    /// Why the run stopped: [`StopReason::Completed`] when the schedule
    /// cooled to its final temperature.
    pub stop_reason: StopReason,
    /// Number of reheats performed.
    pub reheat_count: u32,
    /// Fitness history (sampled at temperature changes).
    pub history: Vec<f64>,
}

/// Simulated Annealing runner.
pub struct SaRunner<P: SaProblem> {
    config: SaConfig,
    problem: P,
    cancelled: Arc<AtomicBool>,
}

impl<P: SaProblem> SaRunner<P> {
    /// Creates a new SA runner.
    ///
    /// # Errors
    /// As [`SaConfig::validate`].
    pub fn new(config: SaConfig, problem: P) -> Result<Self> {
        Self::with_cancellation(config, problem, Arc::new(AtomicBool::new(false)))
    }

    /// Creates a runner that stops when `cancelled` is set — the caller's own
    /// flag, observed directly, so no thread is needed to forward it.
    ///
    /// # Errors
    /// As [`SaConfig::validate`].
    pub fn with_cancellation(
        config: SaConfig,
        problem: P,
        cancelled: Arc<AtomicBool>,
    ) -> Result<Self> {
        config.validate()?;
        Ok(Self {
            config,
            problem,
            cancelled,
        })
    }

    /// Returns a handle to cancel the algorithm.
    pub fn cancel_handle(&self) -> Arc<AtomicBool> {
        self.cancelled.clone()
    }

    /// Runs the Simulated Annealing algorithm.
    pub fn run(&self) -> SaResult<P::Solution> {
        self.run_with_rng(&mut rand::rng())
    }

    /// Runs the Simulated Annealing algorithm with a specific RNG.
    pub fn run_with_rng<R: Rng>(&self, rng: &mut R) -> SaResult<P::Solution> {
        self.run_with_rng_and_progress::<R, fn(SaProgress)>(rng, None)
    }

    /// Runs the Simulated Annealing algorithm with a specific RNG, reporting
    /// each step to `progress_callback` when one is given. The
    /// callback only observes: the run is the same with or without it.
    pub fn run_with_rng_and_progress<R: Rng, F>(
        &self,
        rng: &mut R,
        progress_callback: Option<F>,
    ) -> SaResult<P::Solution>
    where
        F: Fn(SaProgress),
    {
        let start = Timer::now();
        let report = |progress: SaProgress| {
            if let Some(ref callback) = progress_callback {
                callback(progress);
            }
        };
        let mut history = Vec::new();

        // Initialize
        let mut current = self.problem.initial_solution(rng);
        self.problem.evaluate(&mut current);
        let mut best = current.clone();
        let mut best_fitness = best.objective();

        // Temperatures are configured in units of `scale` (1 when absolute); a
        // calibrated run measures it on neighbours of the start, which are not
        // iterations and stop at the time limit.
        let operators = self.problem.available_operators();
        let scale = match self.config.temperature_scale {
            TemperatureScale::Absolute => 1.0,
            TemperatureScale::Calibrated { samples } => {
                let from = current.objective();
                let mut worsenings = Vec::with_capacity(samples);
                for _ in 0..samples {
                    if expired(&start, self.config.time_limit)
                        || self.cancelled.load(Ordering::Relaxed)
                    {
                        break;
                    }
                    let operator = operators[rng.random_range(0..operators.len())];
                    let mut neighbor = self.problem.neighbor(&current, operator, rng);
                    self.problem.evaluate(&mut neighbor);
                    worsenings.push(from - neighbor.objective());
                }
                TemperatureScale::cost_scale(worsenings)
            }
        };
        let initial_temp = self.config.initial_temp * scale;
        let final_temp = self.config.final_temp * scale;
        let mut temperature = initial_temp;
        let mut iteration = 0u64;
        let reached = |fitness: f64| self.config.target_fitness.is_some_and(|t| fitness >= t);
        let mut reheat_count = 0u32;
        let mut stagnation_count = 0u64;

        let temp_delta = if matches!(self.config.cooling_schedule, CoolingSchedule::Linear) {
            (initial_temp - final_temp)
                / (self.config.max_iterations.unwrap_or(10000) as f64
                    / self.config.iterations_per_temp as f64)
        } else {
            0.0
        };

        // For adaptive cooling
        let mut accepted_count = 0usize;
        let mut total_count = 0usize;

        // Each exit names its reason; the target is judged first, so a target
        // met by the last step is reported as met.
        let stop_reason = loop {
            if reached(best_fitness) {
                break StopReason::TargetReached;
            }
            if temperature <= final_temp {
                break StopReason::Completed;
            }
            if self
                .config
                .max_iterations
                .is_some_and(|max| iteration >= max)
            {
                break StopReason::IterationLimit;
            }
            if self.cancelled.load(Ordering::Relaxed) {
                break StopReason::Cancelled;
            }
            if expired(&start, self.config.time_limit) {
                break StopReason::TimeLimit;
            }

            // Iterations at this temperature — each evaluation can be a full
            // placement, so the time limit is checked per iteration, not per step
            for _ in 0..self.config.iterations_per_temp {
                if expired(&start, self.config.time_limit) {
                    break;
                }
                iteration += 1;
                total_count += 1;

                // Select random operator
                let operator = operators[rng.random_range(0..operators.len())];

                // Generate neighbor
                let mut neighbor = self.problem.neighbor(&current, operator, rng);
                self.problem.evaluate(&mut neighbor);

                let current_obj = current.objective();
                let neighbor_obj = neighbor.objective();
                let delta = neighbor_obj - current_obj;

                // Accept or reject
                let accept = if delta >= 0.0 {
                    // Better solution - always accept
                    true
                } else {
                    // Worse solution - accept with probability exp(delta/T)
                    let probability = (delta / temperature).exp();
                    rng.random::<f64>() < probability
                };

                if accept {
                    accepted_count += 1;
                    current = neighbor;

                    // Update best
                    if current.objective() > best_fitness {
                        best = current.clone();
                        best_fitness = best.objective();
                        stagnation_count = 0;
                    } else {
                        stagnation_count += 1;
                    }
                } else {
                    stagnation_count += 1;
                }

                // A report per step: one step of a large problem can take a
                // while, and a caller that cancels from its callback has to
                // be heard between them.
                report(SaProgress {
                    temperature,
                    iteration,
                    best_fitness,
                    current_fitness: current.objective(),
                    acceptance_rate: accepted_count as f64 / total_count.max(1) as f64,
                    elapsed: start.elapsed(),
                    running: true,
                });

                // Check max iterations inside inner loop
                if let Some(max) = self.config.max_iterations {
                    if iteration >= max {
                        break;
                    }
                }
            }

            // Record history
            history.push(best_fitness);

            // Callback
            self.problem
                .on_temperature_change(temperature, iteration, &best, &current);

            // Reheating check
            if self.config.enable_reheating && stagnation_count >= self.config.reheat_threshold {
                temperature *= self.config.reheat_factor;
                temperature = temperature.min(initial_temp);
                stagnation_count = 0;
                reheat_count += 1;
            }

            // Cool down
            temperature = self.cool_down(
                temperature,
                temp_delta,
                (final_temp, scale),
                accepted_count,
                total_count,
            );

            // Reset adaptive counters
            accepted_count = 0;
            total_count = 0;
        };

        // Final history entry
        history.push(best_fitness);
        report(SaProgress {
            temperature,
            iteration,
            best_fitness,
            current_fitness: current.objective(),
            acceptance_rate: 0.0,
            elapsed: start.elapsed(),
            running: false,
        });

        SaResult {
            best,
            final_temperature: temperature,
            temperature_unit: scale,
            iterations: iteration,
            elapsed: start.elapsed(),
            stop_reason,
            reheat_count,
            history,
        }
    }

    /// Apply cooling schedule.
    fn cool_down(
        &self,
        current_temp: f64,
        delta: f64,
        (final_temp, scale): (f64, f64),
        accepted: usize,
        total: usize,
    ) -> f64 {
        match self.config.cooling_schedule {
            CoolingSchedule::Geometric => current_temp * self.config.cooling_rate,
            CoolingSchedule::Linear => (current_temp - delta).max(final_temp),
            CoolingSchedule::Adaptive => {
                // Adjust cooling rate based on acceptance rate
                let acceptance_rate = if total > 0 {
                    accepted as f64 / total as f64
                } else {
                    0.5
                };

                // If acceptance rate is high, cool faster; if low, cool slower
                let adjusted_rate = if acceptance_rate > 0.5 {
                    self.config.cooling_rate * 0.95 // Cool faster
                } else if acceptance_rate < 0.1 {
                    self.config.cooling_rate.powf(0.5) // Cool slower (sqrt)
                } else {
                    self.config.cooling_rate
                };

                current_temp * adjusted_rate
            }
            CoolingSchedule::LundyMees => {
                // T_new = T / (1 + beta * T)
                // `cooling_rate` has the unit of 1 / temperature, so per `scale`.
                current_temp / (1.0 + self.config.cooling_rate / scale * current_temp)
            }
        }
    }
}

/// One of `count` options other than `current`, uniformly (`count` > 1).
pub fn another_option<R: Rng>(current: usize, count: usize, rng: &mut R) -> usize {
    let pick = rng.random_range(0..count - 1);
    if pick >= current.min(count - 1) {
        pick + 1
    } else {
        pick
    }
}

/// Permutation-based solution for SA.
#[derive(Debug, Clone)]
pub struct PermutationSolution {
    /// The permutation (indices).
    pub sequence: Vec<usize>,
    /// Rotation/orientation index per item (by item, not by sequence
    /// position — a reorder moves items, not their rotations).
    pub rotations: Vec<usize>,
    /// How many rotation/orientation options each item has: `rotations[i]`
    /// is always below `rotation_counts[i]`.
    pub rotation_counts: std::sync::Arc<[usize]>,
    /// Mirror flag per item (`allow_flip` support, 2D nesting only — see
    /// [`NeighborhoodOperator::MirrorFlip`]).
    pub mirrors: Vec<bool>,
    /// Cached objective value.
    objective: f64,
}

impl PermutationSolution {
    /// The identity order, every item at its first rotation; one item per
    /// entry of `rotation_counts` (each item's number of options).
    pub fn new(rotation_counts: impl Into<std::sync::Arc<[usize]>>) -> Self {
        let rotation_counts = rotation_counts.into();
        let size = rotation_counts.len();
        Self {
            sequence: (0..size).collect(),
            rotations: vec![0; size],
            rotation_counts,
            mirrors: vec![false; size],
            objective: f64::NEG_INFINITY,
        }
    }

    /// A random order, each item at a rotation drawn uniformly from its own
    /// options.
    pub fn random<R: Rng>(
        rotation_counts: impl Into<std::sync::Arc<[usize]>>,
        rng: &mut R,
    ) -> Self {
        let rotation_counts = rotation_counts.into();
        let size = rotation_counts.len();
        let mut sequence: Vec<usize> = (0..size).collect();
        sequence.shuffle(rng);

        let rotations: Vec<usize> = rotation_counts
            .iter()
            .map(|&k| rng.random_range(0..k.max(1)))
            .collect();

        let mirrors: Vec<bool> = (0..size).map(|_| rng.random()).collect();

        Self {
            sequence,
            rotations,
            rotation_counts,
            mirrors,
            objective: f64::NEG_INFINITY,
        }
    }

    /// Returns the length of the sequence.
    pub fn len(&self) -> usize {
        self.sequence.len()
    }

    /// Returns true if empty.
    pub fn is_empty(&self) -> bool {
        self.sequence.is_empty()
    }

    /// Applies swap operator: swaps two elements in sequence.
    pub fn apply_swap<R: Rng>(&self, rng: &mut R) -> Self {
        let mut result = self.clone();
        if result.sequence.len() < 2 {
            return result;
        }

        let i = rng.random_range(0..result.sequence.len());
        let j = rng.random_range(0..result.sequence.len());
        result.sequence.swap(i, j);
        result.objective = f64::NEG_INFINITY;
        result
    }

    /// Applies relocate operator: moves an element to a new position.
    pub fn apply_relocate<R: Rng>(&self, rng: &mut R) -> Self {
        let mut result = self.clone();
        if result.sequence.len() < 2 {
            return result;
        }

        let from = rng.random_range(0..result.sequence.len());
        let to = rng.random_range(0..result.sequence.len());

        if from != to {
            let elem = result.sequence.remove(from);
            let insert_pos = if to > from { to - 1 } else { to };
            result
                .sequence
                .insert(insert_pos.min(result.sequence.len()), elem);
        }

        result.objective = f64::NEG_INFINITY;
        result
    }

    /// Applies inversion operator: reverses a segment.
    pub fn apply_inversion<R: Rng>(&self, rng: &mut R) -> Self {
        let mut result = self.clone();
        let n = result.sequence.len();
        if n < 2 {
            return result;
        }

        let (mut p1, mut p2) = (rng.random_range(0..n), rng.random_range(0..n));
        if p1 > p2 {
            std::mem::swap(&mut p1, &mut p2);
        }

        result.sequence[p1..=p2].reverse();
        result.objective = f64::NEG_INFINITY;
        result
    }

    /// Applies rotation operator: moves one item that has more than one
    /// option to another of its options, drawn uniformly. An item with one
    /// option is never chosen, so the move always changes the solution.
    pub fn apply_rotation<R: Rng>(&self, rng: &mut R) -> Self {
        let mut result = self.clone();
        let turnable: Vec<usize> = (0..result.rotations.len())
            .filter(|&i| result.rotation_counts.get(i).copied().unwrap_or(1) > 1)
            .collect();
        if turnable.is_empty() {
            return result;
        }

        let idx = turnable[rng.random_range(0..turnable.len())];
        result.rotations[idx] =
            another_option(result.rotations[idx], result.rotation_counts[idx], rng);
        result.objective = f64::NEG_INFINITY;
        result
    }

    /// Applies mirror-flip operator: flips one element's mirror bit
    /// (`allow_flip` support). Same shape as `apply_rotation` — unconditional
    /// on the gene, the consuming problem masks it against its own
    /// per-item mirror-eligibility (e.g. `Geometry2D::allow_flip()`).
    pub fn apply_mirror_flip<R: Rng>(&self, rng: &mut R) -> Self {
        let mut result = self.clone();
        if result.mirrors.is_empty() {
            return result;
        }

        let idx = rng.random_range(0..result.mirrors.len());
        result.mirrors[idx] = !result.mirrors[idx];
        result.objective = f64::NEG_INFINITY;
        result
    }

    /// Applies chain operator: 3-opt style move.
    pub fn apply_chain<R: Rng>(&self, rng: &mut R) -> Self {
        let mut result = self.clone();
        let n = result.sequence.len();
        if n < 4 {
            // Fall back to swap for small sequences
            return self.apply_swap(rng);
        }

        // Select three distinct positions
        let mut positions: Vec<usize> = (0..n).collect();
        positions.shuffle(rng);
        let mut selected: Vec<usize> = positions.into_iter().take(3).collect();
        selected.sort();

        let (p1, p2, p3) = (selected[0], selected[1], selected[2]);

        // Rotate segments: [0..p1] [p1..p2] [p2..p3] [p3..n]
        // New order: [0..p1] [p2..p3] [p1..p2] [p3..n]
        let seg1: Vec<usize> = result.sequence[..p1].to_vec();
        let seg2: Vec<usize> = result.sequence[p1..p2].to_vec();
        let seg3: Vec<usize> = result.sequence[p2..p3].to_vec();
        let seg4: Vec<usize> = result.sequence[p3..].to_vec();

        result.sequence = [seg1, seg3, seg2, seg4].concat();
        result.objective = f64::NEG_INFINITY;
        result
    }
}

impl SaSolution for PermutationSolution {
    fn objective(&self) -> f64 {
        self.objective
    }

    fn set_objective(&mut self, value: f64) {
        self.objective = value;
    }
}

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

    #[test]
    fn a_rotation_move_turns_only_items_with_options_within_their_own() {
        use rand::rngs::StdRng;
        use rand::SeedableRng;
        let mut rng = StdRng::seed_from_u64(7);
        let counts: Vec<usize> = vec![1, 3, 1, 8];
        let start = PermutationSolution::random(counts.clone(), &mut rng);
        for (r, k) in start.rotations.iter().zip(&counts) {
            assert!(r < k);
        }
        let mut reached = vec![std::collections::HashSet::new(); counts.len()];
        let mut current = start;
        for _ in 0..2000 {
            let next = current.apply_rotation(&mut rng);
            let changed: Vec<usize> = (0..counts.len())
                .filter(|&i| next.rotations[i] != current.rotations[i])
                .collect();
            assert_eq!(changed.len(), 1, "a move changes exactly one item");
            assert!(counts[changed[0]] > 1, "a fixed item was turned");
            for (i, &r) in next.rotations.iter().enumerate() {
                assert!(r < counts[i]);
                reached[i].insert(r);
            }
            current = next;
        }
        // Every option of every item is reachable: 3 and 8 of them.
        assert_eq!(reached[1].len(), 3);
        assert_eq!(reached[3].len(), 8);
        assert_eq!(reached[0].len(), 1);
        // Nothing to turn: the move leaves the solution as it is.
        let fixed = PermutationSolution::new(vec![1; 4]);
        assert_eq!(fixed.apply_rotation(&mut rng).rotations, vec![0; 4]);
    }

    struct SimpleMaxProblem {
        size: usize,
    }

    impl SaProblem for SimpleMaxProblem {
        type Solution = PermutationSolution;

        fn initial_solution<R: Rng>(&self, rng: &mut R) -> Self::Solution {
            PermutationSolution::random(vec![1; self.size], rng)
        }

        fn neighbor<R: Rng>(
            &self,
            solution: &Self::Solution,
            operator: NeighborhoodOperator,
            rng: &mut R,
        ) -> Self::Solution {
            match operator {
                NeighborhoodOperator::Swap => solution.apply_swap(rng),
                NeighborhoodOperator::Relocate => solution.apply_relocate(rng),
                NeighborhoodOperator::Inversion => solution.apply_inversion(rng),
                NeighborhoodOperator::Rotation => solution.apply_rotation(rng),
                NeighborhoodOperator::Chain => solution.apply_chain(rng),
                NeighborhoodOperator::MirrorFlip => solution.apply_mirror_flip(rng),
            }
        }

        fn evaluate(&self, solution: &mut Self::Solution) {
            // Maximize: sequence should be in ascending order
            // Fitness = negative sum of inversions
            let mut inversions = 0i64;
            for i in 0..solution.sequence.len() {
                for j in (i + 1)..solution.sequence.len() {
                    if solution.sequence[i] > solution.sequence[j] {
                        inversions += 1;
                    }
                }
            }
            solution.set_objective(-inversions as f64);
        }
    }

    #[test]
    fn test_sa_basic() {
        let config = SaConfig::default()
            .with_initial_temp(100.0)
            .with_final_temp(0.1)
            .with_cooling_rate(0.9)
            .with_iterations_per_temp(50)
            .with_max_iterations(5000);

        let problem = SimpleMaxProblem { size: 10 };
        let runner = SaRunner::new(config, problem).expect("valid config");
        let result = runner.run();

        // Should find something reasonably good (fewer inversions)
        assert!(result.best.objective() > -20.0);
        assert!(result.iterations > 0);
    }

    /// `SimpleMaxProblem` with every objective multiplied by `unit`.
    struct ScaledMax {
        inner: SimpleMaxProblem,
        unit: f64,
    }

    impl SaProblem for ScaledMax {
        type Solution = PermutationSolution;

        fn initial_solution<R: Rng>(&self, rng: &mut R) -> Self::Solution {
            self.inner.initial_solution(rng)
        }

        fn neighbor<R: Rng>(
            &self,
            solution: &Self::Solution,
            operator: NeighborhoodOperator,
            rng: &mut R,
        ) -> Self::Solution {
            self.inner.neighbor(solution, operator, rng)
        }

        fn evaluate(&self, solution: &mut Self::Solution) {
            self.inner.evaluate(solution);
            solution.set_objective(solution.objective() * self.unit);
        }
    }

    /// A calibrated temperature is relative to the objective: the same
    /// configuration on the same problem with every objective multiplied by a
    /// power of two (exact in floating point) follows the same trajectory,
    /// which an absolute temperature cannot do.
    #[test]
    fn a_calibrated_run_does_not_depend_on_the_unit_of_the_objective() {
        use rand::rngs::StdRng;
        let config = SaConfig::default()
            .with_relative_temperature(1.0, 0.01, 20)
            .with_iterations_per_temp(10)
            .with_max_iterations(1500);
        let run = |unit: f64| {
            let problem = ScaledMax {
                inner: SimpleMaxProblem { size: 12 },
                unit,
            };
            let runner = SaRunner::new(config.clone(), problem).expect("valid config");
            runner.run_with_rng(&mut StdRng::seed_from_u64(9))
        };
        let (small, large) = (run(1.0), run(1024.0));
        assert_eq!(small.best.sequence, large.best.sequence);
        assert_eq!(small.iterations, large.iterations);
        assert_eq!(small.best.objective() * 1024.0, large.best.objective());
        assert_eq!(small.temperature_unit * 1024.0, large.temperature_unit);
    }

    /// An absolute temperature is used as configured.
    #[test]
    fn an_absolute_temperature_is_in_the_unit_of_the_objective() {
        let config = SaConfig::default()
            .with_initial_temp(10.0)
            .with_final_temp(0.5)
            .with_max_iterations(200);
        let runner = SaRunner::new(config, SimpleMaxProblem { size: 6 }).expect("valid config");
        let result = runner.run();
        assert_eq!(result.temperature_unit, 1.0);
        assert!(result.final_temperature <= 10.0);
    }

    #[test]
    fn test_cooling_schedules() {
        let problem = SimpleMaxProblem { size: 5 };

        for schedule in [
            CoolingSchedule::Geometric,
            CoolingSchedule::Linear,
            CoolingSchedule::Adaptive,
            CoolingSchedule::LundyMees,
        ] {
            let config = SaConfig::default()
                .with_cooling_schedule(schedule)
                .with_max_iterations(1000);

            let runner = SaRunner::new(config, problem.clone()).expect("valid config");
            let result = runner.run();

            // Should complete without panic
            assert!(result.iterations > 0);
        }
    }

    #[test]
    fn test_neighborhood_operators() {
        let mut rng = rand::rng();
        let solution = PermutationSolution::random(vec![4; 10], &mut rng);

        // Test all operators produce valid permutations
        let swap = solution.apply_swap(&mut rng);
        let relocate = solution.apply_relocate(&mut rng);
        let inversion = solution.apply_inversion(&mut rng);
        let rotation = solution.apply_rotation(&mut rng);
        let chain = solution.apply_chain(&mut rng);

        for sol in [&swap, &relocate, &inversion, &rotation, &chain] {
            let mut sorted = sol.sequence.clone();
            sorted.sort();
            assert_eq!(sorted, (0..10).collect::<Vec<_>>());
        }
    }

    #[test]
    fn test_reheating() {
        let config = SaConfig::default()
            .with_initial_temp(10.0)
            .with_final_temp(0.1)
            .with_max_iterations(500)
            .with_reheating(50, 1.5);

        let problem = SimpleMaxProblem { size: 8 };
        let runner = SaRunner::new(config, problem).expect("valid config");
        let result = runner.run();

        // Should complete (reheating may or may not trigger)
        assert!(result.iterations > 0);
    }

    impl Clone for SimpleMaxProblem {
        fn clone(&self) -> Self {
            Self { size: self.size }
        }
    }
}