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_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.max(0.001);
102 self
103 }
104
105 pub fn with_final_temp(mut self, temp: f64) -> Self {
107 self.final_temp = temp.max(0.0001);
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<()> {
125 check_range("cooling_rate", self.cooling_rate, 0.001, 0.9999)?;
126 check_at_least("initial_temp", self.initial_temp, f64::MIN_POSITIVE)?;
127 check_at_least("final_temp", self.final_temp, f64::MIN_POSITIVE)?;
128 check_range("final_temp", self.final_temp, 0.0, self.initial_temp)?;
129 check_at_least("reheat_factor", self.reheat_factor, 1.0)?;
130 Ok(())
131 }
132
133 pub fn with_iterations_per_temp(mut self, iterations: usize) -> Self {
135 self.iterations_per_temp = iterations.max(1);
136 self
137 }
138
139 pub fn with_max_iterations(mut self, iterations: u64) -> Self {
141 self.max_iterations = Some(iterations);
142 self
143 }
144
145 pub fn with_cooling_schedule(mut self, schedule: CoolingSchedule) -> Self {
147 self.cooling_schedule = schedule;
148 self
149 }
150
151 pub fn with_time_limit(mut self, duration: Duration) -> Self {
153 self.time_limit = Some(duration);
154 self
155 }
156
157 pub fn with_target_fitness(mut self, fitness: f64) -> Self {
159 self.target_fitness = Some(fitness);
160 self
161 }
162
163 pub fn with_reheating(mut self, threshold: u64, factor: f64) -> Self {
165 self.enable_reheating = true;
166 self.reheat_threshold = threshold;
167 self.reheat_factor = factor.max(1.1);
168 self
169 }
170}
171
172pub trait SaSolution: Clone + Send + Sync {
174 fn objective(&self) -> f64;
177
178 fn set_objective(&mut self, value: f64);
180}
181
182#[derive(Debug, Clone, Copy, PartialEq, Eq)]
184pub enum NeighborhoodOperator {
185 Swap,
187 Relocate,
189 Inversion,
191 Rotation,
193 Chain,
195 MirrorFlip,
200}
201
202pub trait SaProblem: Send + Sync {
204 type Solution: SaSolution;
206
207 fn initial_solution<R: Rng>(&self, rng: &mut R) -> Self::Solution;
209
210 fn neighbor<R: Rng>(
212 &self,
213 solution: &Self::Solution,
214 operator: NeighborhoodOperator,
215 rng: &mut R,
216 ) -> Self::Solution;
217
218 fn evaluate(&self, solution: &mut Self::Solution);
220
221 fn available_operators(&self) -> Vec<NeighborhoodOperator> {
223 vec![
224 NeighborhoodOperator::Swap,
225 NeighborhoodOperator::Relocate,
226 NeighborhoodOperator::Inversion,
227 ]
228 }
229
230 fn on_temperature_change(
232 &self,
233 _temperature: f64,
234 _iteration: u64,
235 _best: &Self::Solution,
236 _current: &Self::Solution,
237 ) {
238 }
240}
241
242#[derive(Debug, Clone)]
244pub struct SaProgress {
245 pub temperature: f64,
247 pub iteration: u64,
249 pub best_fitness: f64,
251 pub current_fitness: f64,
253 pub acceptance_rate: f64,
255 pub elapsed: Duration,
257 pub running: bool,
259}
260
261#[derive(Debug, Clone)]
263pub struct SaResult<S: SaSolution> {
264 pub best: S,
266 pub final_temperature: f64,
268 pub iterations: u64,
270 pub elapsed: Duration,
272 pub target_reached: bool,
274 pub reheat_count: u32,
276 pub history: Vec<f64>,
278}
279
280pub struct SaRunner<P: SaProblem> {
282 config: SaConfig,
283 problem: P,
284 cancelled: Arc<AtomicBool>,
285}
286
287impl<P: SaProblem> SaRunner<P> {
288 pub fn new(config: SaConfig, problem: P) -> Result<Self> {
293 Self::with_cancellation(config, problem, Arc::new(AtomicBool::new(false)))
294 }
295
296 pub fn with_cancellation(
302 config: SaConfig,
303 problem: P,
304 cancelled: Arc<AtomicBool>,
305 ) -> Result<Self> {
306 config.validate()?;
307 Ok(Self {
308 config,
309 problem,
310 cancelled,
311 })
312 }
313
314 pub fn cancel_handle(&self) -> Arc<AtomicBool> {
316 self.cancelled.clone()
317 }
318
319 pub fn run(&self) -> SaResult<P::Solution> {
321 self.run_with_rng(&mut rand::rng())
322 }
323
324 pub fn run_with_rng<R: Rng>(&self, rng: &mut R) -> SaResult<P::Solution> {
326 let start = Timer::now();
327 let mut history = Vec::new();
328
329 let mut current = self.problem.initial_solution(rng);
331 self.problem.evaluate(&mut current);
332 let mut best = current.clone();
333 let mut best_fitness = best.objective();
334
335 let mut temperature = self.config.initial_temp;
336 let mut iteration = 0u64;
337 let mut target_reached = false;
338 let mut reheat_count = 0u32;
339 let mut stagnation_count = 0u64;
340
341 let operators = self.problem.available_operators();
342 let temp_delta = if matches!(self.config.cooling_schedule, CoolingSchedule::Linear) {
343 (self.config.initial_temp - self.config.final_temp)
344 / (self.config.max_iterations.unwrap_or(10000) as f64
345 / self.config.iterations_per_temp as f64)
346 } else {
347 0.0
348 };
349
350 let mut accepted_count = 0usize;
352 let mut total_count = 0usize;
353
354 while temperature > self.config.final_temp {
355 if self.cancelled.load(Ordering::Relaxed) {
357 break;
358 }
359
360 if expired(&start, self.config.time_limit) {
362 break;
363 }
364
365 if let Some(max) = self.config.max_iterations {
367 if iteration >= max {
368 break;
369 }
370 }
371
372 if let Some(target) = self.config.target_fitness {
374 if best_fitness >= target {
375 target_reached = true;
376 break;
377 }
378 }
379
380 for _ in 0..self.config.iterations_per_temp {
383 if expired(&start, self.config.time_limit) {
384 break;
385 }
386 iteration += 1;
387 total_count += 1;
388
389 let operator = operators[rng.random_range(0..operators.len())];
391
392 let mut neighbor = self.problem.neighbor(¤t, operator, rng);
394 self.problem.evaluate(&mut neighbor);
395
396 let current_obj = current.objective();
397 let neighbor_obj = neighbor.objective();
398 let delta = neighbor_obj - current_obj;
399
400 let accept = if delta >= 0.0 {
402 true
404 } else {
405 let probability = (delta / temperature).exp();
407 rng.random::<f64>() < probability
408 };
409
410 if accept {
411 accepted_count += 1;
412 current = neighbor;
413
414 if current.objective() > best_fitness {
416 best = current.clone();
417 best_fitness = best.objective();
418 stagnation_count = 0;
419 } else {
420 stagnation_count += 1;
421 }
422 } else {
423 stagnation_count += 1;
424 }
425
426 if let Some(max) = self.config.max_iterations {
428 if iteration >= max {
429 break;
430 }
431 }
432 }
433
434 history.push(best_fitness);
436
437 self.problem
439 .on_temperature_change(temperature, iteration, &best, ¤t);
440
441 if self.config.enable_reheating && stagnation_count >= self.config.reheat_threshold {
443 temperature *= self.config.reheat_factor;
444 temperature = temperature.min(self.config.initial_temp);
445 stagnation_count = 0;
446 reheat_count += 1;
447 }
448
449 temperature = self.cool_down(temperature, temp_delta, accepted_count, total_count);
451
452 accepted_count = 0;
454 total_count = 0;
455 }
456
457 history.push(best_fitness);
459
460 SaResult {
461 best,
462 final_temperature: temperature,
463 iterations: iteration,
464 elapsed: start.elapsed(),
465 target_reached,
466 reheat_count,
467 history,
468 }
469 }
470
471 fn cool_down(&self, current_temp: f64, delta: f64, accepted: usize, total: usize) -> f64 {
473 match self.config.cooling_schedule {
474 CoolingSchedule::Geometric => current_temp * self.config.cooling_rate,
475 CoolingSchedule::Linear => (current_temp - delta).max(self.config.final_temp),
476 CoolingSchedule::Adaptive => {
477 let acceptance_rate = if total > 0 {
479 accepted as f64 / total as f64
480 } else {
481 0.5
482 };
483
484 let adjusted_rate = if acceptance_rate > 0.5 {
486 self.config.cooling_rate * 0.95 } else if acceptance_rate < 0.1 {
488 self.config.cooling_rate.powf(0.5) } else {
490 self.config.cooling_rate
491 };
492
493 current_temp * adjusted_rate
494 }
495 CoolingSchedule::LundyMees => {
496 current_temp / (1.0 + self.config.cooling_rate * current_temp)
498 }
499 }
500 }
501
502 #[cfg(feature = "parallel")]
515 pub fn run_parallel(&self, num_restarts: usize) -> SaResult<P::Solution>
516 where
517 P: Clone,
518 {
519 let num_restarts = num_restarts.max(1);
520
521 let results: Vec<SaResult<P::Solution>> = (0..num_restarts)
523 .into_par_iter()
524 .map(|_| {
525 let mut rng = rand::rng();
526 self.run_with_rng(&mut rng)
527 })
528 .collect();
529
530 results
532 .into_iter()
533 .max_by(|a, b| {
534 a.best
535 .objective()
536 .partial_cmp(&b.best.objective())
537 .unwrap_or(std::cmp::Ordering::Equal)
538 })
539 .expect("At least one result should exist")
540 }
541}
542
543#[derive(Debug, Clone)]
545pub struct PermutationSolution {
546 pub sequence: Vec<usize>,
548 pub rotations: Vec<usize>,
550 pub rotation_options: usize,
552 pub mirrors: Vec<bool>,
555 objective: f64,
557}
558
559impl PermutationSolution {
560 pub fn new(size: usize, rotation_options: usize) -> Self {
562 Self {
563 sequence: (0..size).collect(),
564 rotations: vec![0; size],
565 rotation_options,
566 mirrors: vec![false; size],
567 objective: f64::NEG_INFINITY,
568 }
569 }
570
571 pub fn random<R: Rng>(size: usize, rotation_options: usize, rng: &mut R) -> Self {
573 let mut sequence: Vec<usize> = (0..size).collect();
574 sequence.shuffle(rng);
575
576 let rotations: Vec<usize> = (0..size)
577 .map(|_| rng.random_range(0..rotation_options.max(1)))
578 .collect();
579
580 let mirrors: Vec<bool> = (0..size).map(|_| rng.random()).collect();
581
582 Self {
583 sequence,
584 rotations,
585 rotation_options,
586 mirrors,
587 objective: f64::NEG_INFINITY,
588 }
589 }
590
591 pub fn len(&self) -> usize {
593 self.sequence.len()
594 }
595
596 pub fn is_empty(&self) -> bool {
598 self.sequence.is_empty()
599 }
600
601 pub fn apply_swap<R: Rng>(&self, rng: &mut R) -> Self {
603 let mut result = self.clone();
604 if result.sequence.len() < 2 {
605 return result;
606 }
607
608 let i = rng.random_range(0..result.sequence.len());
609 let j = rng.random_range(0..result.sequence.len());
610 result.sequence.swap(i, j);
611 result.objective = f64::NEG_INFINITY;
612 result
613 }
614
615 pub fn apply_relocate<R: Rng>(&self, rng: &mut R) -> Self {
617 let mut result = self.clone();
618 if result.sequence.len() < 2 {
619 return result;
620 }
621
622 let from = rng.random_range(0..result.sequence.len());
623 let to = rng.random_range(0..result.sequence.len());
624
625 if from != to {
626 let elem = result.sequence.remove(from);
627 let insert_pos = if to > from { to - 1 } else { to };
628 result
629 .sequence
630 .insert(insert_pos.min(result.sequence.len()), elem);
631 }
632
633 result.objective = f64::NEG_INFINITY;
634 result
635 }
636
637 pub fn apply_inversion<R: Rng>(&self, rng: &mut R) -> Self {
639 let mut result = self.clone();
640 let n = result.sequence.len();
641 if n < 2 {
642 return result;
643 }
644
645 let (mut p1, mut p2) = (rng.random_range(0..n), rng.random_range(0..n));
646 if p1 > p2 {
647 std::mem::swap(&mut p1, &mut p2);
648 }
649
650 result.sequence[p1..=p2].reverse();
651 result.objective = f64::NEG_INFINITY;
652 result
653 }
654
655 pub fn apply_rotation<R: Rng>(&self, rng: &mut R) -> Self {
657 let mut result = self.clone();
658 if result.rotations.is_empty() || result.rotation_options <= 1 {
659 return result;
660 }
661
662 let idx = rng.random_range(0..result.rotations.len());
663 result.rotations[idx] = rng.random_range(0..result.rotation_options);
664 result.objective = f64::NEG_INFINITY;
665 result
666 }
667
668 pub fn apply_mirror_flip<R: Rng>(&self, rng: &mut R) -> Self {
673 let mut result = self.clone();
674 if result.mirrors.is_empty() {
675 return result;
676 }
677
678 let idx = rng.random_range(0..result.mirrors.len());
679 result.mirrors[idx] = !result.mirrors[idx];
680 result.objective = f64::NEG_INFINITY;
681 result
682 }
683
684 pub fn apply_chain<R: Rng>(&self, rng: &mut R) -> Self {
686 let mut result = self.clone();
687 let n = result.sequence.len();
688 if n < 4 {
689 return self.apply_swap(rng);
691 }
692
693 let mut positions: Vec<usize> = (0..n).collect();
695 positions.shuffle(rng);
696 let mut selected: Vec<usize> = positions.into_iter().take(3).collect();
697 selected.sort();
698
699 let (p1, p2, p3) = (selected[0], selected[1], selected[2]);
700
701 let seg1: Vec<usize> = result.sequence[..p1].to_vec();
704 let seg2: Vec<usize> = result.sequence[p1..p2].to_vec();
705 let seg3: Vec<usize> = result.sequence[p2..p3].to_vec();
706 let seg4: Vec<usize> = result.sequence[p3..].to_vec();
707
708 result.sequence = [seg1, seg3, seg2, seg4].concat();
709 result.objective = f64::NEG_INFINITY;
710 result
711 }
712}
713
714impl SaSolution for PermutationSolution {
715 fn objective(&self) -> f64 {
716 self.objective
717 }
718
719 fn set_objective(&mut self, value: f64) {
720 self.objective = value;
721 }
722}
723
724#[cfg(test)]
725mod tests {
726 use super::*;
727
728 struct SimpleMaxProblem {
729 size: usize,
730 }
731
732 impl SaProblem for SimpleMaxProblem {
733 type Solution = PermutationSolution;
734
735 fn initial_solution<R: Rng>(&self, rng: &mut R) -> Self::Solution {
736 PermutationSolution::random(self.size, 1, rng)
737 }
738
739 fn neighbor<R: Rng>(
740 &self,
741 solution: &Self::Solution,
742 operator: NeighborhoodOperator,
743 rng: &mut R,
744 ) -> Self::Solution {
745 match operator {
746 NeighborhoodOperator::Swap => solution.apply_swap(rng),
747 NeighborhoodOperator::Relocate => solution.apply_relocate(rng),
748 NeighborhoodOperator::Inversion => solution.apply_inversion(rng),
749 NeighborhoodOperator::Rotation => solution.apply_rotation(rng),
750 NeighborhoodOperator::Chain => solution.apply_chain(rng),
751 NeighborhoodOperator::MirrorFlip => solution.apply_mirror_flip(rng),
752 }
753 }
754
755 fn evaluate(&self, solution: &mut Self::Solution) {
756 let mut inversions = 0i64;
759 for i in 0..solution.sequence.len() {
760 for j in (i + 1)..solution.sequence.len() {
761 if solution.sequence[i] > solution.sequence[j] {
762 inversions += 1;
763 }
764 }
765 }
766 solution.set_objective(-inversions as f64);
767 }
768 }
769
770 #[test]
771 fn test_sa_basic() {
772 let config = SaConfig::default()
773 .with_initial_temp(100.0)
774 .with_final_temp(0.1)
775 .with_cooling_rate(0.9)
776 .with_iterations_per_temp(50)
777 .with_max_iterations(5000);
778
779 let problem = SimpleMaxProblem { size: 10 };
780 let runner = SaRunner::new(config, problem).expect("valid config");
781 let result = runner.run();
782
783 assert!(result.best.objective() > -20.0);
785 assert!(result.iterations > 0);
786 }
787
788 #[test]
789 fn test_cooling_schedules() {
790 let problem = SimpleMaxProblem { size: 5 };
791
792 for schedule in [
793 CoolingSchedule::Geometric,
794 CoolingSchedule::Linear,
795 CoolingSchedule::Adaptive,
796 CoolingSchedule::LundyMees,
797 ] {
798 let config = SaConfig::default()
799 .with_cooling_schedule(schedule)
800 .with_max_iterations(1000);
801
802 let runner = SaRunner::new(config, problem.clone()).expect("valid config");
803 let result = runner.run();
804
805 assert!(result.iterations > 0);
807 }
808 }
809
810 #[test]
811 fn test_neighborhood_operators() {
812 let mut rng = rand::rng();
813 let solution = PermutationSolution::random(10, 4, &mut rng);
814
815 let swap = solution.apply_swap(&mut rng);
817 let relocate = solution.apply_relocate(&mut rng);
818 let inversion = solution.apply_inversion(&mut rng);
819 let rotation = solution.apply_rotation(&mut rng);
820 let chain = solution.apply_chain(&mut rng);
821
822 for sol in [&swap, &relocate, &inversion, &rotation, &chain] {
823 let mut sorted = sol.sequence.clone();
824 sorted.sort();
825 assert_eq!(sorted, (0..10).collect::<Vec<_>>());
826 }
827 }
828
829 #[test]
830 fn test_reheating() {
831 let config = SaConfig::default()
832 .with_initial_temp(10.0)
833 .with_final_temp(0.1)
834 .with_max_iterations(500)
835 .with_reheating(50, 1.5);
836
837 let problem = SimpleMaxProblem { size: 8 };
838 let runner = SaRunner::new(config, problem).expect("valid config");
839 let result = runner.run();
840
841 assert!(result.iterations > 0);
843 }
844
845 impl Clone for SimpleMaxProblem {
846 fn clone(&self) -> Self {
847 Self { size: self.size }
848 }
849 }
850}