1use rand::prelude::*;
19use std::sync::atomic::{AtomicBool, Ordering};
20use std::sync::Arc;
21use std::time::Duration;
22
23use crate::error::{check_above, check_at_least, check_range, Result};
24use crate::result::StopReason;
25use crate::timing::{expired, Timer};
26
27#[cfg(feature = "serde")]
28use serde::{Deserialize, Serialize};
29
30#[derive(Debug, Clone, Copy, Default)]
32#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
33pub enum CoolingSchedule {
34 #[default]
36 Geometric,
37 Linear,
39 Adaptive,
41 LundyMees,
43}
44
45#[derive(Debug, Clone)]
47#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
48pub struct SaConfig {
49 pub initial_temp: f64,
51 pub final_temp: f64,
53 pub cooling_rate: f64,
55 pub iterations_per_temp: usize,
57 pub max_iterations: Option<u64>,
59 pub cooling_schedule: CoolingSchedule,
61 pub time_limit: Option<Duration>,
63 pub target_fitness: Option<f64>,
65 pub enable_reheating: bool,
67 pub reheat_threshold: u64,
69 pub reheat_factor: f64,
71}
72
73impl Default for SaConfig {
74 fn default() -> Self {
75 Self {
76 initial_temp: 1000.0,
77 final_temp: 0.001,
78 cooling_rate: 0.95,
79 iterations_per_temp: 100,
80 max_iterations: Some(100_000),
81 cooling_schedule: CoolingSchedule::Geometric,
82 time_limit: None,
83 target_fitness: None,
84 enable_reheating: false,
85 reheat_threshold: 1000,
86 reheat_factor: 2.0,
87 }
88 }
89}
90
91impl SaConfig {
92 pub fn new() -> Self {
94 Self::default()
95 }
96
97 pub fn with_initial_temp(mut self, temp: f64) -> Self {
99 self.initial_temp = temp;
100 self
101 }
102
103 pub fn with_final_temp(mut self, temp: f64) -> Self {
105 self.final_temp = temp;
106 self
107 }
108
109 pub fn with_cooling_rate(mut self, rate: f64) -> Self {
111 self.cooling_rate = rate;
112 self
113 }
114
115 pub fn validate(&self) -> Result<()> {
124 check_range("cooling_rate", self.cooling_rate, 0.001, 0.9999)?;
125 check_above("initial_temp", self.initial_temp, 0.0)?;
126 check_above("final_temp", self.final_temp, 0.0)?;
127 check_range("final_temp", self.final_temp, 0.0, self.initial_temp)?;
128 check_at_least("iterations_per_temp", self.iterations_per_temp as f64, 1.0)?;
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;
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;
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 stop_reason: StopReason,
275 pub reheat_count: u32,
277 pub history: Vec<f64>,
279}
280
281pub struct SaRunner<P: SaProblem> {
283 config: SaConfig,
284 problem: P,
285 cancelled: Arc<AtomicBool>,
286}
287
288impl<P: SaProblem> SaRunner<P> {
289 pub fn new(config: SaConfig, problem: P) -> Result<Self> {
294 Self::with_cancellation(config, problem, Arc::new(AtomicBool::new(false)))
295 }
296
297 pub fn with_cancellation(
303 config: SaConfig,
304 problem: P,
305 cancelled: Arc<AtomicBool>,
306 ) -> Result<Self> {
307 config.validate()?;
308 Ok(Self {
309 config,
310 problem,
311 cancelled,
312 })
313 }
314
315 pub fn cancel_handle(&self) -> Arc<AtomicBool> {
317 self.cancelled.clone()
318 }
319
320 pub fn run(&self) -> SaResult<P::Solution> {
322 self.run_with_rng(&mut rand::rng())
323 }
324
325 pub fn run_with_rng<R: Rng>(&self, rng: &mut R) -> SaResult<P::Solution> {
327 self.run_with_rng_and_progress::<R, fn(SaProgress)>(rng, None)
328 }
329
330 pub fn run_with_rng_and_progress<R: Rng, F>(
334 &self,
335 rng: &mut R,
336 progress_callback: Option<F>,
337 ) -> SaResult<P::Solution>
338 where
339 F: Fn(SaProgress),
340 {
341 let start = Timer::now();
342 let report = |progress: SaProgress| {
343 if let Some(ref callback) = progress_callback {
344 callback(progress);
345 }
346 };
347 let mut history = Vec::new();
348
349 let mut current = self.problem.initial_solution(rng);
351 self.problem.evaluate(&mut current);
352 let mut best = current.clone();
353 let mut best_fitness = best.objective();
354
355 let mut temperature = self.config.initial_temp;
356 let mut iteration = 0u64;
357 let reached = |fitness: f64| self.config.target_fitness.is_some_and(|t| fitness >= t);
358 let mut reheat_count = 0u32;
359 let mut stagnation_count = 0u64;
360
361 let operators = self.problem.available_operators();
362 let temp_delta = if matches!(self.config.cooling_schedule, CoolingSchedule::Linear) {
363 (self.config.initial_temp - self.config.final_temp)
364 / (self.config.max_iterations.unwrap_or(10000) as f64
365 / self.config.iterations_per_temp as f64)
366 } else {
367 0.0
368 };
369
370 let mut accepted_count = 0usize;
372 let mut total_count = 0usize;
373
374 let stop_reason = loop {
377 if reached(best_fitness) {
378 break StopReason::TargetReached;
379 }
380 if temperature <= self.config.final_temp {
381 break StopReason::Completed;
382 }
383 if self
384 .config
385 .max_iterations
386 .is_some_and(|max| iteration >= max)
387 {
388 break StopReason::IterationLimit;
389 }
390 if self.cancelled.load(Ordering::Relaxed) {
391 break StopReason::Cancelled;
392 }
393 if expired(&start, self.config.time_limit) {
394 break StopReason::TimeLimit;
395 }
396
397 for _ in 0..self.config.iterations_per_temp {
400 if expired(&start, self.config.time_limit) {
401 break;
402 }
403 iteration += 1;
404 total_count += 1;
405
406 let operator = operators[rng.random_range(0..operators.len())];
408
409 let mut neighbor = self.problem.neighbor(¤t, operator, rng);
411 self.problem.evaluate(&mut neighbor);
412
413 let current_obj = current.objective();
414 let neighbor_obj = neighbor.objective();
415 let delta = neighbor_obj - current_obj;
416
417 let accept = if delta >= 0.0 {
419 true
421 } else {
422 let probability = (delta / temperature).exp();
424 rng.random::<f64>() < probability
425 };
426
427 if accept {
428 accepted_count += 1;
429 current = neighbor;
430
431 if current.objective() > best_fitness {
433 best = current.clone();
434 best_fitness = best.objective();
435 stagnation_count = 0;
436 } else {
437 stagnation_count += 1;
438 }
439 } else {
440 stagnation_count += 1;
441 }
442
443 report(SaProgress {
447 temperature,
448 iteration,
449 best_fitness,
450 current_fitness: current.objective(),
451 acceptance_rate: accepted_count as f64 / total_count.max(1) as f64,
452 elapsed: start.elapsed(),
453 running: true,
454 });
455
456 if let Some(max) = self.config.max_iterations {
458 if iteration >= max {
459 break;
460 }
461 }
462 }
463
464 history.push(best_fitness);
466
467 self.problem
469 .on_temperature_change(temperature, iteration, &best, ¤t);
470
471 if self.config.enable_reheating && stagnation_count >= self.config.reheat_threshold {
473 temperature *= self.config.reheat_factor;
474 temperature = temperature.min(self.config.initial_temp);
475 stagnation_count = 0;
476 reheat_count += 1;
477 }
478
479 temperature = self.cool_down(temperature, temp_delta, accepted_count, total_count);
481
482 accepted_count = 0;
484 total_count = 0;
485 };
486
487 history.push(best_fitness);
489 report(SaProgress {
490 temperature,
491 iteration,
492 best_fitness,
493 current_fitness: current.objective(),
494 acceptance_rate: 0.0,
495 elapsed: start.elapsed(),
496 running: false,
497 });
498
499 SaResult {
500 best,
501 final_temperature: temperature,
502 iterations: iteration,
503 elapsed: start.elapsed(),
504 stop_reason,
505 reheat_count,
506 history,
507 }
508 }
509
510 fn cool_down(&self, current_temp: f64, delta: f64, accepted: usize, total: usize) -> f64 {
512 match self.config.cooling_schedule {
513 CoolingSchedule::Geometric => current_temp * self.config.cooling_rate,
514 CoolingSchedule::Linear => (current_temp - delta).max(self.config.final_temp),
515 CoolingSchedule::Adaptive => {
516 let acceptance_rate = if total > 0 {
518 accepted as f64 / total as f64
519 } else {
520 0.5
521 };
522
523 let adjusted_rate = if acceptance_rate > 0.5 {
525 self.config.cooling_rate * 0.95 } else if acceptance_rate < 0.1 {
527 self.config.cooling_rate.powf(0.5) } else {
529 self.config.cooling_rate
530 };
531
532 current_temp * adjusted_rate
533 }
534 CoolingSchedule::LundyMees => {
535 current_temp / (1.0 + self.config.cooling_rate * current_temp)
537 }
538 }
539 }
540}
541
542pub fn another_option<R: Rng>(current: usize, count: usize, rng: &mut R) -> usize {
544 let pick = rng.random_range(0..count - 1);
545 if pick >= current.min(count - 1) {
546 pick + 1
547 } else {
548 pick
549 }
550}
551
552#[derive(Debug, Clone)]
554pub struct PermutationSolution {
555 pub sequence: Vec<usize>,
557 pub rotations: Vec<usize>,
560 pub rotation_counts: std::sync::Arc<[usize]>,
563 pub mirrors: Vec<bool>,
566 objective: f64,
568}
569
570impl PermutationSolution {
571 pub fn new(rotation_counts: impl Into<std::sync::Arc<[usize]>>) -> Self {
574 let rotation_counts = rotation_counts.into();
575 let size = rotation_counts.len();
576 Self {
577 sequence: (0..size).collect(),
578 rotations: vec![0; size],
579 rotation_counts,
580 mirrors: vec![false; size],
581 objective: f64::NEG_INFINITY,
582 }
583 }
584
585 pub fn random<R: Rng>(
588 rotation_counts: impl Into<std::sync::Arc<[usize]>>,
589 rng: &mut R,
590 ) -> Self {
591 let rotation_counts = rotation_counts.into();
592 let size = rotation_counts.len();
593 let mut sequence: Vec<usize> = (0..size).collect();
594 sequence.shuffle(rng);
595
596 let rotations: Vec<usize> = rotation_counts
597 .iter()
598 .map(|&k| rng.random_range(0..k.max(1)))
599 .collect();
600
601 let mirrors: Vec<bool> = (0..size).map(|_| rng.random()).collect();
602
603 Self {
604 sequence,
605 rotations,
606 rotation_counts,
607 mirrors,
608 objective: f64::NEG_INFINITY,
609 }
610 }
611
612 pub fn len(&self) -> usize {
614 self.sequence.len()
615 }
616
617 pub fn is_empty(&self) -> bool {
619 self.sequence.is_empty()
620 }
621
622 pub fn apply_swap<R: Rng>(&self, rng: &mut R) -> Self {
624 let mut result = self.clone();
625 if result.sequence.len() < 2 {
626 return result;
627 }
628
629 let i = rng.random_range(0..result.sequence.len());
630 let j = rng.random_range(0..result.sequence.len());
631 result.sequence.swap(i, j);
632 result.objective = f64::NEG_INFINITY;
633 result
634 }
635
636 pub fn apply_relocate<R: Rng>(&self, rng: &mut R) -> Self {
638 let mut result = self.clone();
639 if result.sequence.len() < 2 {
640 return result;
641 }
642
643 let from = rng.random_range(0..result.sequence.len());
644 let to = rng.random_range(0..result.sequence.len());
645
646 if from != to {
647 let elem = result.sequence.remove(from);
648 let insert_pos = if to > from { to - 1 } else { to };
649 result
650 .sequence
651 .insert(insert_pos.min(result.sequence.len()), elem);
652 }
653
654 result.objective = f64::NEG_INFINITY;
655 result
656 }
657
658 pub fn apply_inversion<R: Rng>(&self, rng: &mut R) -> Self {
660 let mut result = self.clone();
661 let n = result.sequence.len();
662 if n < 2 {
663 return result;
664 }
665
666 let (mut p1, mut p2) = (rng.random_range(0..n), rng.random_range(0..n));
667 if p1 > p2 {
668 std::mem::swap(&mut p1, &mut p2);
669 }
670
671 result.sequence[p1..=p2].reverse();
672 result.objective = f64::NEG_INFINITY;
673 result
674 }
675
676 pub fn apply_rotation<R: Rng>(&self, rng: &mut R) -> Self {
680 let mut result = self.clone();
681 let turnable: Vec<usize> = (0..result.rotations.len())
682 .filter(|&i| result.rotation_counts.get(i).copied().unwrap_or(1) > 1)
683 .collect();
684 if turnable.is_empty() {
685 return result;
686 }
687
688 let idx = turnable[rng.random_range(0..turnable.len())];
689 result.rotations[idx] =
690 another_option(result.rotations[idx], result.rotation_counts[idx], rng);
691 result.objective = f64::NEG_INFINITY;
692 result
693 }
694
695 pub fn apply_mirror_flip<R: Rng>(&self, rng: &mut R) -> Self {
700 let mut result = self.clone();
701 if result.mirrors.is_empty() {
702 return result;
703 }
704
705 let idx = rng.random_range(0..result.mirrors.len());
706 result.mirrors[idx] = !result.mirrors[idx];
707 result.objective = f64::NEG_INFINITY;
708 result
709 }
710
711 pub fn apply_chain<R: Rng>(&self, rng: &mut R) -> Self {
713 let mut result = self.clone();
714 let n = result.sequence.len();
715 if n < 4 {
716 return self.apply_swap(rng);
718 }
719
720 let mut positions: Vec<usize> = (0..n).collect();
722 positions.shuffle(rng);
723 let mut selected: Vec<usize> = positions.into_iter().take(3).collect();
724 selected.sort();
725
726 let (p1, p2, p3) = (selected[0], selected[1], selected[2]);
727
728 let seg1: Vec<usize> = result.sequence[..p1].to_vec();
731 let seg2: Vec<usize> = result.sequence[p1..p2].to_vec();
732 let seg3: Vec<usize> = result.sequence[p2..p3].to_vec();
733 let seg4: Vec<usize> = result.sequence[p3..].to_vec();
734
735 result.sequence = [seg1, seg3, seg2, seg4].concat();
736 result.objective = f64::NEG_INFINITY;
737 result
738 }
739}
740
741impl SaSolution for PermutationSolution {
742 fn objective(&self) -> f64 {
743 self.objective
744 }
745
746 fn set_objective(&mut self, value: f64) {
747 self.objective = value;
748 }
749}
750
751#[cfg(test)]
752mod tests {
753 use super::*;
754
755 #[test]
756 fn a_rotation_move_turns_only_items_with_options_within_their_own() {
757 use rand::rngs::StdRng;
758 use rand::SeedableRng;
759 let mut rng = StdRng::seed_from_u64(7);
760 let counts: Vec<usize> = vec![1, 3, 1, 8];
761 let start = PermutationSolution::random(counts.clone(), &mut rng);
762 for (r, k) in start.rotations.iter().zip(&counts) {
763 assert!(r < k);
764 }
765 let mut reached = vec![std::collections::HashSet::new(); counts.len()];
766 let mut current = start;
767 for _ in 0..2000 {
768 let next = current.apply_rotation(&mut rng);
769 let changed: Vec<usize> = (0..counts.len())
770 .filter(|&i| next.rotations[i] != current.rotations[i])
771 .collect();
772 assert_eq!(changed.len(), 1, "a move changes exactly one item");
773 assert!(counts[changed[0]] > 1, "a fixed item was turned");
774 for (i, &r) in next.rotations.iter().enumerate() {
775 assert!(r < counts[i]);
776 reached[i].insert(r);
777 }
778 current = next;
779 }
780 assert_eq!(reached[1].len(), 3);
782 assert_eq!(reached[3].len(), 8);
783 assert_eq!(reached[0].len(), 1);
784 let fixed = PermutationSolution::new(vec![1; 4]);
786 assert_eq!(fixed.apply_rotation(&mut rng).rotations, vec![0; 4]);
787 }
788
789 struct SimpleMaxProblem {
790 size: usize,
791 }
792
793 impl SaProblem for SimpleMaxProblem {
794 type Solution = PermutationSolution;
795
796 fn initial_solution<R: Rng>(&self, rng: &mut R) -> Self::Solution {
797 PermutationSolution::random(vec![1; self.size], rng)
798 }
799
800 fn neighbor<R: Rng>(
801 &self,
802 solution: &Self::Solution,
803 operator: NeighborhoodOperator,
804 rng: &mut R,
805 ) -> Self::Solution {
806 match operator {
807 NeighborhoodOperator::Swap => solution.apply_swap(rng),
808 NeighborhoodOperator::Relocate => solution.apply_relocate(rng),
809 NeighborhoodOperator::Inversion => solution.apply_inversion(rng),
810 NeighborhoodOperator::Rotation => solution.apply_rotation(rng),
811 NeighborhoodOperator::Chain => solution.apply_chain(rng),
812 NeighborhoodOperator::MirrorFlip => solution.apply_mirror_flip(rng),
813 }
814 }
815
816 fn evaluate(&self, solution: &mut Self::Solution) {
817 let mut inversions = 0i64;
820 for i in 0..solution.sequence.len() {
821 for j in (i + 1)..solution.sequence.len() {
822 if solution.sequence[i] > solution.sequence[j] {
823 inversions += 1;
824 }
825 }
826 }
827 solution.set_objective(-inversions as f64);
828 }
829 }
830
831 #[test]
832 fn test_sa_basic() {
833 let config = SaConfig::default()
834 .with_initial_temp(100.0)
835 .with_final_temp(0.1)
836 .with_cooling_rate(0.9)
837 .with_iterations_per_temp(50)
838 .with_max_iterations(5000);
839
840 let problem = SimpleMaxProblem { size: 10 };
841 let runner = SaRunner::new(config, problem).expect("valid config");
842 let result = runner.run();
843
844 assert!(result.best.objective() > -20.0);
846 assert!(result.iterations > 0);
847 }
848
849 #[test]
850 fn test_cooling_schedules() {
851 let problem = SimpleMaxProblem { size: 5 };
852
853 for schedule in [
854 CoolingSchedule::Geometric,
855 CoolingSchedule::Linear,
856 CoolingSchedule::Adaptive,
857 CoolingSchedule::LundyMees,
858 ] {
859 let config = SaConfig::default()
860 .with_cooling_schedule(schedule)
861 .with_max_iterations(1000);
862
863 let runner = SaRunner::new(config, problem.clone()).expect("valid config");
864 let result = runner.run();
865
866 assert!(result.iterations > 0);
868 }
869 }
870
871 #[test]
872 fn test_neighborhood_operators() {
873 let mut rng = rand::rng();
874 let solution = PermutationSolution::random(vec![4; 10], &mut rng);
875
876 let swap = solution.apply_swap(&mut rng);
878 let relocate = solution.apply_relocate(&mut rng);
879 let inversion = solution.apply_inversion(&mut rng);
880 let rotation = solution.apply_rotation(&mut rng);
881 let chain = solution.apply_chain(&mut rng);
882
883 for sol in [&swap, &relocate, &inversion, &rotation, &chain] {
884 let mut sorted = sol.sequence.clone();
885 sorted.sort();
886 assert_eq!(sorted, (0..10).collect::<Vec<_>>());
887 }
888 }
889
890 #[test]
891 fn test_reheating() {
892 let config = SaConfig::default()
893 .with_initial_temp(10.0)
894 .with_final_temp(0.1)
895 .with_max_iterations(500)
896 .with_reheating(50, 1.5);
897
898 let problem = SimpleMaxProblem { size: 8 };
899 let runner = SaRunner::new(config, problem).expect("valid config");
900 let result = runner.run();
901
902 assert!(result.iterations > 0);
904 }
905
906 impl Clone for SimpleMaxProblem {
907 fn clone(&self) -> Self {
908 Self { size: self.size }
909 }
910 }
911}