1use rand::prelude::*;
20#[cfg(feature = "parallel")]
21use rayon::prelude::*;
22use std::sync::atomic::{AtomicBool, Ordering};
23use std::sync::Arc;
24use std::time::Duration;
25
26use crate::error::{check_above, check_at_least, check_range, Result};
27use crate::timing::{expired, Timer};
28
29#[cfg(feature = "serde")]
30use serde::{Deserialize, Serialize};
31
32#[derive(Debug, Clone, Copy, Default)]
34#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
35pub enum CoolingSchedule {
36 #[default]
38 Geometric,
39 Linear,
41 Adaptive,
43 LundyMees,
45}
46
47#[derive(Debug, Clone)]
49#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
50pub struct SaConfig {
51 pub initial_temp: f64,
53 pub final_temp: f64,
55 pub cooling_rate: f64,
57 pub iterations_per_temp: usize,
59 pub max_iterations: Option<u64>,
61 pub cooling_schedule: CoolingSchedule,
63 pub time_limit: Option<Duration>,
65 pub target_fitness: Option<f64>,
67 pub enable_reheating: bool,
69 pub reheat_threshold: u64,
71 pub reheat_factor: f64,
73}
74
75impl Default for SaConfig {
76 fn default() -> Self {
77 Self {
78 initial_temp: 1000.0,
79 final_temp: 0.001,
80 cooling_rate: 0.95,
81 iterations_per_temp: 100,
82 max_iterations: Some(100_000),
83 cooling_schedule: CoolingSchedule::Geometric,
84 time_limit: None,
85 target_fitness: None,
86 enable_reheating: false,
87 reheat_threshold: 1000,
88 reheat_factor: 2.0,
89 }
90 }
91}
92
93impl SaConfig {
94 pub fn new() -> Self {
96 Self::default()
97 }
98
99 pub fn with_initial_temp(mut self, temp: f64) -> Self {
101 self.initial_temp = temp;
102 self
103 }
104
105 pub fn with_final_temp(mut self, temp: f64) -> Self {
107 self.final_temp = temp;
108 self
109 }
110
111 pub fn with_cooling_rate(mut self, rate: f64) -> Self {
113 self.cooling_rate = rate;
114 self
115 }
116
117 pub fn validate(&self) -> Result<()> {
126 check_range("cooling_rate", self.cooling_rate, 0.001, 0.9999)?;
127 check_above("initial_temp", self.initial_temp, 0.0)?;
128 check_above("final_temp", self.final_temp, 0.0)?;
129 check_range("final_temp", self.final_temp, 0.0, self.initial_temp)?;
130 check_at_least("iterations_per_temp", self.iterations_per_temp as f64, 1.0)?;
131 check_at_least("reheat_factor", self.reheat_factor, 1.0)?;
132 Ok(())
133 }
134
135 pub fn with_iterations_per_temp(mut self, iterations: usize) -> Self {
137 self.iterations_per_temp = iterations;
138 self
139 }
140
141 pub fn with_max_iterations(mut self, iterations: u64) -> Self {
143 self.max_iterations = Some(iterations);
144 self
145 }
146
147 pub fn with_cooling_schedule(mut self, schedule: CoolingSchedule) -> Self {
149 self.cooling_schedule = schedule;
150 self
151 }
152
153 pub fn with_time_limit(mut self, duration: Duration) -> Self {
155 self.time_limit = Some(duration);
156 self
157 }
158
159 pub fn with_target_fitness(mut self, fitness: f64) -> Self {
161 self.target_fitness = Some(fitness);
162 self
163 }
164
165 pub fn with_reheating(mut self, threshold: u64, factor: f64) -> Self {
167 self.enable_reheating = true;
168 self.reheat_threshold = threshold;
169 self.reheat_factor = factor;
170 self
171 }
172}
173
174pub trait SaSolution: Clone + Send + Sync {
176 fn objective(&self) -> f64;
179
180 fn set_objective(&mut self, value: f64);
182}
183
184#[derive(Debug, Clone, Copy, PartialEq, Eq)]
186pub enum NeighborhoodOperator {
187 Swap,
189 Relocate,
191 Inversion,
193 Rotation,
195 Chain,
197 MirrorFlip,
202}
203
204pub trait SaProblem: Send + Sync {
206 type Solution: SaSolution;
208
209 fn initial_solution<R: Rng>(&self, rng: &mut R) -> Self::Solution;
211
212 fn neighbor<R: Rng>(
214 &self,
215 solution: &Self::Solution,
216 operator: NeighborhoodOperator,
217 rng: &mut R,
218 ) -> Self::Solution;
219
220 fn evaluate(&self, solution: &mut Self::Solution);
222
223 fn available_operators(&self) -> Vec<NeighborhoodOperator> {
225 vec![
226 NeighborhoodOperator::Swap,
227 NeighborhoodOperator::Relocate,
228 NeighborhoodOperator::Inversion,
229 ]
230 }
231
232 fn on_temperature_change(
234 &self,
235 _temperature: f64,
236 _iteration: u64,
237 _best: &Self::Solution,
238 _current: &Self::Solution,
239 ) {
240 }
242}
243
244#[derive(Debug, Clone)]
246pub struct SaProgress {
247 pub temperature: f64,
249 pub iteration: u64,
251 pub best_fitness: f64,
253 pub current_fitness: f64,
255 pub acceptance_rate: f64,
257 pub elapsed: Duration,
259 pub running: bool,
261}
262
263#[derive(Debug, Clone)]
265pub struct SaResult<S: SaSolution> {
266 pub best: S,
268 pub final_temperature: f64,
270 pub iterations: u64,
272 pub elapsed: Duration,
274 pub target_reached: bool,
276 pub reheat_count: u32,
278 pub history: Vec<f64>,
280}
281
282pub struct SaRunner<P: SaProblem> {
284 config: SaConfig,
285 problem: P,
286 cancelled: Arc<AtomicBool>,
287}
288
289impl<P: SaProblem> SaRunner<P> {
290 pub fn new(config: SaConfig, problem: P) -> Result<Self> {
295 Self::with_cancellation(config, problem, Arc::new(AtomicBool::new(false)))
296 }
297
298 pub fn with_cancellation(
304 config: SaConfig,
305 problem: P,
306 cancelled: Arc<AtomicBool>,
307 ) -> Result<Self> {
308 config.validate()?;
309 Ok(Self {
310 config,
311 problem,
312 cancelled,
313 })
314 }
315
316 pub fn cancel_handle(&self) -> Arc<AtomicBool> {
318 self.cancelled.clone()
319 }
320
321 pub fn run(&self) -> SaResult<P::Solution> {
323 self.run_with_rng(&mut rand::rng())
324 }
325
326 pub fn run_with_rng<R: Rng>(&self, rng: &mut R) -> SaResult<P::Solution> {
328 let start = Timer::now();
329 let mut history = Vec::new();
330
331 let mut current = self.problem.initial_solution(rng);
333 self.problem.evaluate(&mut current);
334 let mut best = current.clone();
335 let mut best_fitness = best.objective();
336
337 let mut temperature = self.config.initial_temp;
338 let mut iteration = 0u64;
339 let mut target_reached = false;
340 let mut reheat_count = 0u32;
341 let mut stagnation_count = 0u64;
342
343 let operators = self.problem.available_operators();
344 let temp_delta = if matches!(self.config.cooling_schedule, CoolingSchedule::Linear) {
345 (self.config.initial_temp - self.config.final_temp)
346 / (self.config.max_iterations.unwrap_or(10000) as f64
347 / self.config.iterations_per_temp as f64)
348 } else {
349 0.0
350 };
351
352 let mut accepted_count = 0usize;
354 let mut total_count = 0usize;
355
356 while temperature > self.config.final_temp {
357 if self.cancelled.load(Ordering::Relaxed) {
359 break;
360 }
361
362 if expired(&start, self.config.time_limit) {
364 break;
365 }
366
367 if let Some(max) = self.config.max_iterations {
369 if iteration >= max {
370 break;
371 }
372 }
373
374 if let Some(target) = self.config.target_fitness {
376 if best_fitness >= target {
377 target_reached = true;
378 break;
379 }
380 }
381
382 for _ in 0..self.config.iterations_per_temp {
385 if expired(&start, self.config.time_limit) {
386 break;
387 }
388 iteration += 1;
389 total_count += 1;
390
391 let operator = operators[rng.random_range(0..operators.len())];
393
394 let mut neighbor = self.problem.neighbor(¤t, operator, rng);
396 self.problem.evaluate(&mut neighbor);
397
398 let current_obj = current.objective();
399 let neighbor_obj = neighbor.objective();
400 let delta = neighbor_obj - current_obj;
401
402 let accept = if delta >= 0.0 {
404 true
406 } else {
407 let probability = (delta / temperature).exp();
409 rng.random::<f64>() < probability
410 };
411
412 if accept {
413 accepted_count += 1;
414 current = neighbor;
415
416 if current.objective() > best_fitness {
418 best = current.clone();
419 best_fitness = best.objective();
420 stagnation_count = 0;
421 } else {
422 stagnation_count += 1;
423 }
424 } else {
425 stagnation_count += 1;
426 }
427
428 if let Some(max) = self.config.max_iterations {
430 if iteration >= max {
431 break;
432 }
433 }
434 }
435
436 history.push(best_fitness);
438
439 self.problem
441 .on_temperature_change(temperature, iteration, &best, ¤t);
442
443 if self.config.enable_reheating && stagnation_count >= self.config.reheat_threshold {
445 temperature *= self.config.reheat_factor;
446 temperature = temperature.min(self.config.initial_temp);
447 stagnation_count = 0;
448 reheat_count += 1;
449 }
450
451 temperature = self.cool_down(temperature, temp_delta, accepted_count, total_count);
453
454 accepted_count = 0;
456 total_count = 0;
457 }
458
459 history.push(best_fitness);
461
462 SaResult {
463 best,
464 final_temperature: temperature,
465 iterations: iteration,
466 elapsed: start.elapsed(),
467 target_reached,
468 reheat_count,
469 history,
470 }
471 }
472
473 fn cool_down(&self, current_temp: f64, delta: f64, accepted: usize, total: usize) -> f64 {
475 match self.config.cooling_schedule {
476 CoolingSchedule::Geometric => current_temp * self.config.cooling_rate,
477 CoolingSchedule::Linear => (current_temp - delta).max(self.config.final_temp),
478 CoolingSchedule::Adaptive => {
479 let acceptance_rate = if total > 0 {
481 accepted as f64 / total as f64
482 } else {
483 0.5
484 };
485
486 let adjusted_rate = if acceptance_rate > 0.5 {
488 self.config.cooling_rate * 0.95 } else if acceptance_rate < 0.1 {
490 self.config.cooling_rate.powf(0.5) } else {
492 self.config.cooling_rate
493 };
494
495 current_temp * adjusted_rate
496 }
497 CoolingSchedule::LundyMees => {
498 current_temp / (1.0 + self.config.cooling_rate * current_temp)
500 }
501 }
502 }
503
504 #[cfg(feature = "parallel")]
517 pub fn run_parallel(&self, num_restarts: usize) -> SaResult<P::Solution>
518 where
519 P: Clone,
520 {
521 let num_restarts = num_restarts.max(1);
522
523 let results: Vec<SaResult<P::Solution>> = (0..num_restarts)
525 .into_par_iter()
526 .map(|_| {
527 let mut rng = rand::rng();
528 self.run_with_rng(&mut rng)
529 })
530 .collect();
531
532 results
534 .into_iter()
535 .max_by(|a, b| {
536 a.best
537 .objective()
538 .partial_cmp(&b.best.objective())
539 .unwrap_or(std::cmp::Ordering::Equal)
540 })
541 .expect("At least one result should exist")
542 }
543}
544
545#[derive(Debug, Clone)]
547pub struct PermutationSolution {
548 pub sequence: Vec<usize>,
550 pub rotations: Vec<usize>,
552 pub rotation_options: usize,
554 pub mirrors: Vec<bool>,
557 objective: f64,
559}
560
561impl PermutationSolution {
562 pub fn new(size: usize, rotation_options: usize) -> Self {
564 Self {
565 sequence: (0..size).collect(),
566 rotations: vec![0; size],
567 rotation_options,
568 mirrors: vec![false; size],
569 objective: f64::NEG_INFINITY,
570 }
571 }
572
573 pub fn random<R: Rng>(size: usize, rotation_options: usize, rng: &mut R) -> Self {
575 let mut sequence: Vec<usize> = (0..size).collect();
576 sequence.shuffle(rng);
577
578 let rotations: Vec<usize> = (0..size)
579 .map(|_| rng.random_range(0..rotation_options.max(1)))
580 .collect();
581
582 let mirrors: Vec<bool> = (0..size).map(|_| rng.random()).collect();
583
584 Self {
585 sequence,
586 rotations,
587 rotation_options,
588 mirrors,
589 objective: f64::NEG_INFINITY,
590 }
591 }
592
593 pub fn len(&self) -> usize {
595 self.sequence.len()
596 }
597
598 pub fn is_empty(&self) -> bool {
600 self.sequence.is_empty()
601 }
602
603 pub fn apply_swap<R: Rng>(&self, rng: &mut R) -> Self {
605 let mut result = self.clone();
606 if result.sequence.len() < 2 {
607 return result;
608 }
609
610 let i = rng.random_range(0..result.sequence.len());
611 let j = rng.random_range(0..result.sequence.len());
612 result.sequence.swap(i, j);
613 result.objective = f64::NEG_INFINITY;
614 result
615 }
616
617 pub fn apply_relocate<R: Rng>(&self, rng: &mut R) -> Self {
619 let mut result = self.clone();
620 if result.sequence.len() < 2 {
621 return result;
622 }
623
624 let from = rng.random_range(0..result.sequence.len());
625 let to = rng.random_range(0..result.sequence.len());
626
627 if from != to {
628 let elem = result.sequence.remove(from);
629 let insert_pos = if to > from { to - 1 } else { to };
630 result
631 .sequence
632 .insert(insert_pos.min(result.sequence.len()), elem);
633 }
634
635 result.objective = f64::NEG_INFINITY;
636 result
637 }
638
639 pub fn apply_inversion<R: Rng>(&self, rng: &mut R) -> Self {
641 let mut result = self.clone();
642 let n = result.sequence.len();
643 if n < 2 {
644 return result;
645 }
646
647 let (mut p1, mut p2) = (rng.random_range(0..n), rng.random_range(0..n));
648 if p1 > p2 {
649 std::mem::swap(&mut p1, &mut p2);
650 }
651
652 result.sequence[p1..=p2].reverse();
653 result.objective = f64::NEG_INFINITY;
654 result
655 }
656
657 pub fn apply_rotation<R: Rng>(&self, rng: &mut R) -> Self {
659 let mut result = self.clone();
660 if result.rotations.is_empty() || result.rotation_options <= 1 {
661 return result;
662 }
663
664 let idx = rng.random_range(0..result.rotations.len());
665 result.rotations[idx] = rng.random_range(0..result.rotation_options);
666 result.objective = f64::NEG_INFINITY;
667 result
668 }
669
670 pub fn apply_mirror_flip<R: Rng>(&self, rng: &mut R) -> Self {
675 let mut result = self.clone();
676 if result.mirrors.is_empty() {
677 return result;
678 }
679
680 let idx = rng.random_range(0..result.mirrors.len());
681 result.mirrors[idx] = !result.mirrors[idx];
682 result.objective = f64::NEG_INFINITY;
683 result
684 }
685
686 pub fn apply_chain<R: Rng>(&self, rng: &mut R) -> Self {
688 let mut result = self.clone();
689 let n = result.sequence.len();
690 if n < 4 {
691 return self.apply_swap(rng);
693 }
694
695 let mut positions: Vec<usize> = (0..n).collect();
697 positions.shuffle(rng);
698 let mut selected: Vec<usize> = positions.into_iter().take(3).collect();
699 selected.sort();
700
701 let (p1, p2, p3) = (selected[0], selected[1], selected[2]);
702
703 let seg1: Vec<usize> = result.sequence[..p1].to_vec();
706 let seg2: Vec<usize> = result.sequence[p1..p2].to_vec();
707 let seg3: Vec<usize> = result.sequence[p2..p3].to_vec();
708 let seg4: Vec<usize> = result.sequence[p3..].to_vec();
709
710 result.sequence = [seg1, seg3, seg2, seg4].concat();
711 result.objective = f64::NEG_INFINITY;
712 result
713 }
714}
715
716impl SaSolution for PermutationSolution {
717 fn objective(&self) -> f64 {
718 self.objective
719 }
720
721 fn set_objective(&mut self, value: f64) {
722 self.objective = value;
723 }
724}
725
726#[cfg(test)]
727mod tests {
728 use super::*;
729
730 struct SimpleMaxProblem {
731 size: usize,
732 }
733
734 impl SaProblem for SimpleMaxProblem {
735 type Solution = PermutationSolution;
736
737 fn initial_solution<R: Rng>(&self, rng: &mut R) -> Self::Solution {
738 PermutationSolution::random(self.size, 1, rng)
739 }
740
741 fn neighbor<R: Rng>(
742 &self,
743 solution: &Self::Solution,
744 operator: NeighborhoodOperator,
745 rng: &mut R,
746 ) -> Self::Solution {
747 match operator {
748 NeighborhoodOperator::Swap => solution.apply_swap(rng),
749 NeighborhoodOperator::Relocate => solution.apply_relocate(rng),
750 NeighborhoodOperator::Inversion => solution.apply_inversion(rng),
751 NeighborhoodOperator::Rotation => solution.apply_rotation(rng),
752 NeighborhoodOperator::Chain => solution.apply_chain(rng),
753 NeighborhoodOperator::MirrorFlip => solution.apply_mirror_flip(rng),
754 }
755 }
756
757 fn evaluate(&self, solution: &mut Self::Solution) {
758 let mut inversions = 0i64;
761 for i in 0..solution.sequence.len() {
762 for j in (i + 1)..solution.sequence.len() {
763 if solution.sequence[i] > solution.sequence[j] {
764 inversions += 1;
765 }
766 }
767 }
768 solution.set_objective(-inversions as f64);
769 }
770 }
771
772 #[test]
773 fn test_sa_basic() {
774 let config = SaConfig::default()
775 .with_initial_temp(100.0)
776 .with_final_temp(0.1)
777 .with_cooling_rate(0.9)
778 .with_iterations_per_temp(50)
779 .with_max_iterations(5000);
780
781 let problem = SimpleMaxProblem { size: 10 };
782 let runner = SaRunner::new(config, problem).expect("valid config");
783 let result = runner.run();
784
785 assert!(result.best.objective() > -20.0);
787 assert!(result.iterations > 0);
788 }
789
790 #[test]
791 fn test_cooling_schedules() {
792 let problem = SimpleMaxProblem { size: 5 };
793
794 for schedule in [
795 CoolingSchedule::Geometric,
796 CoolingSchedule::Linear,
797 CoolingSchedule::Adaptive,
798 CoolingSchedule::LundyMees,
799 ] {
800 let config = SaConfig::default()
801 .with_cooling_schedule(schedule)
802 .with_max_iterations(1000);
803
804 let runner = SaRunner::new(config, problem.clone()).expect("valid config");
805 let result = runner.run();
806
807 assert!(result.iterations > 0);
809 }
810 }
811
812 #[test]
813 fn test_neighborhood_operators() {
814 let mut rng = rand::rng();
815 let solution = PermutationSolution::random(10, 4, &mut rng);
816
817 let swap = solution.apply_swap(&mut rng);
819 let relocate = solution.apply_relocate(&mut rng);
820 let inversion = solution.apply_inversion(&mut rng);
821 let rotation = solution.apply_rotation(&mut rng);
822 let chain = solution.apply_chain(&mut rng);
823
824 for sol in [&swap, &relocate, &inversion, &rotation, &chain] {
825 let mut sorted = sol.sequence.clone();
826 sorted.sort();
827 assert_eq!(sorted, (0..10).collect::<Vec<_>>());
828 }
829 }
830
831 #[test]
832 fn test_reheating() {
833 let config = SaConfig::default()
834 .with_initial_temp(10.0)
835 .with_final_temp(0.1)
836 .with_max_iterations(500)
837 .with_reheating(50, 1.5);
838
839 let problem = SimpleMaxProblem { size: 8 };
840 let runner = SaRunner::new(config, problem).expect("valid config");
841 let result = runner.run();
842
843 assert!(result.iterations > 0);
845 }
846
847 impl Clone for SimpleMaxProblem {
848 fn clone(&self) -> Self {
849 Self { size: self.size }
850 }
851 }
852}