quantrs2_anneal/
simulator.rs

1//! Simulated quantum annealing simulator
2//!
3//! This module provides a simulator for quantum annealing, which can be used
4//! to solve optimization problems formulated as Ising models or QUBO problems.
5
6use scirs2_core::random::prelude::*;
7use scirs2_core::random::ChaCha8Rng;
8use scirs2_core::random::{Rng, SeedableRng};
9use std::time::{Duration, Instant};
10use thiserror::Error;
11
12use crate::ising::{IsingError, IsingModel};
13
14/// Errors that can occur during simulated annealing
15#[derive(Error, Debug, Clone)]
16pub enum AnnealingError {
17    /// Error in the underlying Ising model
18    #[error("Ising error: {0}")]
19    IsingError(#[from] IsingError),
20
21    /// Error when the annealing schedule is invalid
22    #[error("Invalid annealing schedule: {0}")]
23    InvalidSchedule(String),
24
25    /// Error when the annealing parameters are invalid
26    #[error("Invalid annealing parameter: {0}")]
27    InvalidParameter(String),
28
29    /// Error when the annealing process times out
30    #[error("Annealing timeout after {0:?}")]
31    Timeout(Duration),
32}
33
34/// Result type for annealing operations
35pub type AnnealingResult<T> = Result<T, AnnealingError>;
36
37/// Transverse field strength schedule for quantum annealing
38///
39/// The transverse field represents the quantum tunneling term in the Hamiltonian,
40/// which decreases over time during the annealing process.
41#[derive(Debug, Clone)]
42pub enum TransverseFieldSchedule {
43    /// Linear schedule: A(t) = `A_0` * (1 - `t/t_f`)
44    Linear,
45
46    /// Exponential schedule: A(t) = `A_0` * exp(-alpha * `t/t_f`)
47    Exponential(f64), // alpha parameter
48
49    /// Custom schedule function: A(t, `t_f`) -> A
50    Custom(fn(f64, f64) -> f64),
51}
52
53impl TransverseFieldSchedule {
54    /// Calculate the transverse field strength at time t
55    #[must_use]
56    pub fn calculate(&self, t: f64, t_f: f64, a_0: f64) -> f64 {
57        match self {
58            Self::Linear => a_0 * (1.0 - t / t_f),
59            Self::Exponential(alpha) => a_0 * (-alpha * t / t_f).exp(),
60            Self::Custom(func) => func(t, t_f),
61        }
62    }
63}
64
65/// Temperature schedule for simulated quantum annealing
66///
67/// The temperature controls the probability of accepting non-improving moves,
68/// and typically decreases over time during the annealing process.
69#[derive(Debug, Clone)]
70pub enum TemperatureSchedule {
71    /// Linear schedule: T(t) = `T_0` * (1 - `t/t_f`)
72    Linear,
73
74    /// Exponential schedule: T(t) = `T_0` * exp(-alpha * `t/t_f`)
75    Exponential(f64), // alpha parameter
76
77    /// Geometric schedule: T(t) = `T_0` * `alpha^(t/delta_t)`
78    Geometric(f64, f64), // alpha and delta_t parameters
79
80    /// Custom schedule function: T(t, `t_f`) -> T
81    Custom(fn(f64, f64) -> f64),
82}
83
84impl TemperatureSchedule {
85    /// Calculate the temperature at time t
86    #[must_use]
87    pub fn calculate(&self, t: f64, t_f: f64, t_0: f64) -> f64 {
88        match self {
89            Self::Linear => t_0 * (1.0 - t / t_f),
90            Self::Exponential(alpha) => t_0 * (-alpha * t / t_f).exp(),
91            Self::Geometric(alpha, delta_t) => t_0 * alpha.powf(t / delta_t),
92            Self::Custom(func) => func(t, t_f),
93        }
94    }
95}
96
97/// Parameters for simulated quantum annealing
98#[derive(Debug, Clone)]
99pub struct AnnealingParams {
100    /// Initial transverse field strength
101    pub initial_transverse_field: f64,
102
103    /// Transverse field schedule
104    pub transverse_field_schedule: TransverseFieldSchedule,
105
106    /// Initial temperature
107    pub initial_temperature: f64,
108
109    /// Final temperature
110    pub final_temperature: f64,
111
112    /// Temperature schedule
113    pub temperature_schedule: TemperatureSchedule,
114
115    /// Number of Monte Carlo steps
116    pub num_sweeps: usize,
117
118    /// Number of spins to update per sweep
119    pub updates_per_sweep: Option<usize>,
120
121    /// Number of repetitions/restarts
122    pub num_repetitions: usize,
123
124    /// Random seed for reproducibility
125    pub seed: Option<u64>,
126
127    /// Maximum runtime in seconds
128    pub timeout: Option<f64>,
129
130    /// Number of Trotter slices for quantum annealing
131    pub trotter_slices: usize,
132}
133
134impl AnnealingParams {
135    /// Create new annealing parameters with default values
136    #[must_use]
137    pub const fn new() -> Self {
138        Self {
139            initial_transverse_field: 2.0,
140            transverse_field_schedule: TransverseFieldSchedule::Linear,
141            initial_temperature: 2.0,
142            final_temperature: 0.01,
143            temperature_schedule: TemperatureSchedule::Exponential(3.0),
144            num_sweeps: 1000,
145            updates_per_sweep: None,
146            num_repetitions: 10,
147            seed: None,
148            timeout: Some(60.0), // 60 seconds default timeout
149            trotter_slices: 20,
150        }
151    }
152
153    /// Validate the annealing parameters
154    pub fn validate(&self) -> AnnealingResult<()> {
155        // Check transverse field
156        if self.initial_transverse_field <= 0.0 || !self.initial_transverse_field.is_finite() {
157            return Err(AnnealingError::InvalidParameter(format!(
158                "Initial transverse field must be positive and finite, got {}",
159                self.initial_transverse_field
160            )));
161        }
162
163        // Check temperature
164        if self.initial_temperature <= 0.0 || !self.initial_temperature.is_finite() {
165            return Err(AnnealingError::InvalidParameter(format!(
166                "Initial temperature must be positive and finite, got {}",
167                self.initial_temperature
168            )));
169        }
170
171        // Check final temperature
172        if self.final_temperature <= 0.0 || !self.final_temperature.is_finite() {
173            return Err(AnnealingError::InvalidParameter(format!(
174                "Final temperature must be positive and finite, got {}",
175                self.final_temperature
176            )));
177        }
178
179        // Check sweeps
180        if self.num_sweeps == 0 {
181            return Err(AnnealingError::InvalidParameter(
182                "Number of sweeps must be positive".to_string(),
183            ));
184        }
185
186        // Check repetitions
187        if self.num_repetitions == 0 {
188            return Err(AnnealingError::InvalidParameter(
189                "Number of repetitions must be positive".to_string(),
190            ));
191        }
192
193        // Check timeout
194        if let Some(timeout) = self.timeout {
195            if timeout <= 0.0 || !timeout.is_finite() {
196                return Err(AnnealingError::InvalidParameter(format!(
197                    "Timeout must be positive and finite, got {timeout}"
198                )));
199            }
200        }
201
202        // Check Trotter slices
203        if self.trotter_slices == 0 {
204            return Err(AnnealingError::InvalidParameter(
205                "Number of Trotter slices must be positive".to_string(),
206            ));
207        }
208
209        Ok(())
210    }
211}
212
213impl Default for AnnealingParams {
214    fn default() -> Self {
215        Self::new()
216    }
217}
218
219/// Result of a simulated quantum annealing run
220#[derive(Debug, Clone)]
221pub struct AnnealingSolution {
222    /// Best spin configuration found
223    pub best_spins: Vec<i8>,
224
225    /// Energy of the best configuration
226    pub best_energy: f64,
227
228    /// Number of repetitions performed
229    pub repetitions: usize,
230
231    /// Total number of sweeps performed
232    pub total_sweeps: usize,
233
234    /// Time taken for the annealing process
235    pub runtime: Duration,
236
237    /// Information about the annealing process
238    pub info: String,
239}
240
241/// Simulated quantum annealing solver
242///
243/// This uses path integral Monte Carlo to simulate quantum annealing,
244/// which can be used to find low-energy states of Ising models.
245#[derive(Debug, Clone)]
246pub struct QuantumAnnealingSimulator {
247    /// Parameters for the annealing process
248    params: AnnealingParams,
249}
250
251impl QuantumAnnealingSimulator {
252    /// Create a new quantum annealing simulator with the given parameters
253    pub fn new(params: AnnealingParams) -> AnnealingResult<Self> {
254        // Validate parameters
255        params.validate()?;
256
257        Ok(Self { params })
258    }
259
260    /// Create a new quantum annealing simulator with default parameters
261    pub fn with_default_params() -> AnnealingResult<Self> {
262        Self::new(AnnealingParams::default())
263    }
264}
265
266impl Default for QuantumAnnealingSimulator {
267    fn default() -> Self {
268        Self::with_default_params().expect("Default parameters should be valid")
269    }
270}
271
272impl QuantumAnnealingSimulator {
273    /// Find the ground state of an Ising model using simulated quantum annealing
274    pub fn solve(&self, model: &IsingModel) -> AnnealingResult<AnnealingSolution> {
275        // Start timer
276        let start_time = Instant::now();
277
278        // Create random number generator
279        let mut rng = match self.params.seed {
280            Some(seed) => ChaCha8Rng::seed_from_u64(seed),
281            None => ChaCha8Rng::seed_from_u64(thread_rng().gen()),
282        };
283
284        // Initialize best result
285        let num_qubits = model.num_qubits;
286        let mut best_spins = vec![1; num_qubits]; // Start with all spins up
287        let mut best_energy = match model.energy(&best_spins) {
288            Ok(energy) => energy,
289            Err(err) => return Err(AnnealingError::IsingError(err)),
290        };
291
292        // Track sweeps and repetitions
293        let mut total_sweeps = 0;
294        let mut completed_repetitions = 0;
295
296        // Determine updates per sweep
297        let updates_per_sweep = self.params.updates_per_sweep.unwrap_or(num_qubits);
298
299        // Prepare for quantum annealing with path integral Monte Carlo
300        let trotter_slices = self.params.trotter_slices;
301
302        // Perform multiple repetitions
303        for _ in 0..self.params.num_repetitions {
304            // Initialize random spin configuration for each Trotter slice
305            let mut trotter_spins = vec![vec![0; num_qubits]; trotter_slices];
306            for slice in &mut trotter_spins {
307                for spin in slice.iter_mut() {
308                    *spin = if rng.random_bool(0.5) { 1 } else { -1 };
309                }
310            }
311
312            // Perform simulated quantum annealing
313            for sweep in 0..self.params.num_sweeps {
314                // Check timeout
315                if let Some(timeout) = self.params.timeout {
316                    let elapsed = start_time.elapsed().as_secs_f64();
317                    if elapsed > timeout {
318                        return Err(AnnealingError::Timeout(Duration::from_secs_f64(elapsed)));
319                    }
320                }
321
322                // Calculate current normalized time
323                let t = sweep as f64 / self.params.num_sweeps as f64;
324                let t_f = 1.0; // Final time (normalized)
325
326                // Calculate current transverse field and temperature
327                let transverse_field = self.params.transverse_field_schedule.calculate(
328                    t,
329                    t_f,
330                    self.params.initial_transverse_field,
331                );
332                let temperature = self.params.temperature_schedule.calculate(
333                    t,
334                    t_f,
335                    self.params.initial_temperature,
336                );
337
338                // Calculate coupling between Trotter slices (J_perp)
339                let j_perp = -0.5
340                    * temperature
341                    * (trotter_slices as f64)
342                    * (transverse_field / temperature).ln_1p().abs();
343
344                // Perform Monte Carlo updates
345                for _ in 0..updates_per_sweep {
346                    // Choose a random qubit and Trotter slice
347                    let qubit = rng.random_range(0..num_qubits);
348                    let slice = rng.random_range(0..trotter_slices);
349
350                    // Calculate energy change from flipping the spin
351                    let current_spin = trotter_spins[slice][qubit];
352                    let new_spin = -current_spin;
353
354                    // Temporary change to calculate energy difference
355                    trotter_spins[slice][qubit] = new_spin;
356
357                    // Calculate energy of this Trotter slice
358                    let mut delta_e = match model.energy(&trotter_spins[slice]) {
359                        Ok(energy_new) => {
360                            // Revert change to calculate original energy
361                            trotter_spins[slice][qubit] = current_spin;
362                            let energy_old = model.energy(&trotter_spins[slice])?;
363                            energy_new - energy_old
364                        }
365                        Err(err) => return Err(AnnealingError::IsingError(err)),
366                    };
367
368                    // Add contribution from neighboring Trotter slices
369                    let prev_slice = (slice + trotter_slices - 1) % trotter_slices;
370                    let next_slice = (slice + 1) % trotter_slices;
371
372                    // Convert spins to f64 for the calculations
373                    let new_spin_f64 = f64::from(new_spin);
374                    let current_spin_f64 = f64::from(current_spin);
375                    let neighbor_sum = f64::from(
376                        trotter_spins[prev_slice][qubit] + trotter_spins[next_slice][qubit],
377                    );
378
379                    delta_e += j_perp * new_spin_f64 * neighbor_sum;
380                    delta_e -= j_perp * current_spin_f64 * neighbor_sum;
381
382                    // Metropolis acceptance criterion
383                    let accept = delta_e <= 0.0 || {
384                        let p = (-delta_e / temperature).exp();
385                        rng.random_range(0.0..1.0) < p
386                    };
387
388                    // Apply the spin flip if accepted
389                    if accept {
390                        trotter_spins[slice][qubit] = new_spin;
391                    }
392                }
393
394                // Increment sweep counter
395                total_sweeps += 1;
396            }
397
398            // After annealing, compute the average spin configuration
399            let mut avg_spins = vec![0; num_qubits];
400            for qubit in 0..num_qubits {
401                let sum: i32 = trotter_spins
402                    .iter()
403                    .map(|slice| i32::from(slice[qubit]))
404                    .sum();
405                avg_spins[qubit] = if sum >= 0 { 1 } else { -1 };
406            }
407
408            // Check if this is a better solution
409            match model.energy(&avg_spins) {
410                Ok(energy) => {
411                    if energy < best_energy {
412                        best_energy = energy;
413                        best_spins = avg_spins;
414                    }
415                }
416                Err(err) => return Err(AnnealingError::IsingError(err)),
417            }
418
419            // Increment repetition counter
420            completed_repetitions += 1;
421        }
422
423        // Calculate runtime
424        let runtime = start_time.elapsed();
425
426        // Build result
427        Ok(AnnealingSolution {
428            best_spins,
429            best_energy,
430            repetitions: completed_repetitions,
431            total_sweeps,
432            runtime,
433            info: format!(
434                "Performed {completed_repetitions} repetitions with {total_sweeps} total sweeps in {runtime:?}"
435            ),
436        })
437    }
438}
439
440/// Classical simulated annealing solver
441///
442/// This uses Metropolis-Hastings algorithm for simulated annealing,
443/// which can be used to find low-energy states of Ising models.
444#[derive(Debug, Clone)]
445pub struct ClassicalAnnealingSimulator {
446    /// Parameters for the annealing process
447    params: AnnealingParams,
448}
449
450impl ClassicalAnnealingSimulator {
451    /// Create a new classical annealing simulator with the given parameters
452    pub fn new(params: AnnealingParams) -> AnnealingResult<Self> {
453        // Validate parameters
454        params.validate()?;
455
456        Ok(Self { params })
457    }
458
459    /// Create a new classical annealing simulator with default parameters
460    pub fn with_default_params() -> AnnealingResult<Self> {
461        Self::new(AnnealingParams::default())
462    }
463}
464
465impl Default for ClassicalAnnealingSimulator {
466    fn default() -> Self {
467        Self::with_default_params().expect("Default parameters should be valid")
468    }
469}
470
471impl ClassicalAnnealingSimulator {
472    /// Find the ground state of an Ising model using classical simulated annealing
473    pub fn solve(&self, model: &IsingModel) -> AnnealingResult<AnnealingSolution> {
474        // Start timer
475        let start_time = Instant::now();
476
477        // Create random number generator
478        let mut rng = match self.params.seed {
479            Some(seed) => ChaCha8Rng::seed_from_u64(seed),
480            None => ChaCha8Rng::seed_from_u64(thread_rng().gen()),
481        };
482
483        // Initialize best result
484        let num_qubits = model.num_qubits;
485        let mut best_spins = vec![1; num_qubits]; // Start with all spins up
486        let mut best_energy = match model.energy(&best_spins) {
487            Ok(energy) => energy,
488            Err(err) => return Err(AnnealingError::IsingError(err)),
489        };
490
491        // Track sweeps and repetitions
492        let mut total_sweeps = 0;
493        let mut completed_repetitions = 0;
494
495        // Determine updates per sweep
496        let updates_per_sweep = self.params.updates_per_sweep.unwrap_or(num_qubits);
497
498        // Perform multiple repetitions
499        for _ in 0..self.params.num_repetitions {
500            // Initialize random spin configuration
501            let mut spins = vec![0; num_qubits];
502            for spin in &mut spins {
503                *spin = if rng.random_bool(0.5) { 1 } else { -1 };
504            }
505
506            // Calculate initial energy
507            let mut current_energy = match model.energy(&spins) {
508                Ok(energy) => energy,
509                Err(err) => return Err(AnnealingError::IsingError(err)),
510            };
511
512            // Perform simulated annealing
513            for sweep in 0..self.params.num_sweeps {
514                // Check timeout
515                if let Some(timeout) = self.params.timeout {
516                    let elapsed = start_time.elapsed().as_secs_f64();
517                    if elapsed > timeout {
518                        return Err(AnnealingError::Timeout(Duration::from_secs_f64(elapsed)));
519                    }
520                }
521
522                // Calculate current normalized time
523                let t = sweep as f64 / self.params.num_sweeps as f64;
524                let t_f = 1.0; // Final time (normalized)
525
526                // Calculate current temperature
527                let temperature = self.params.temperature_schedule.calculate(
528                    t,
529                    t_f,
530                    self.params.initial_temperature,
531                );
532
533                // Perform Monte Carlo updates
534                for _ in 0..updates_per_sweep {
535                    // Choose a random qubit
536                    let qubit = rng.random_range(0..num_qubits);
537
538                    // Calculate energy change from flipping the spin
539                    let current_spin = spins[qubit];
540                    let new_spin = -current_spin;
541
542                    // Temporary change to calculate energy difference
543                    spins[qubit] = new_spin;
544
545                    // Calculate new energy
546                    let new_energy = match model.energy(&spins) {
547                        Ok(energy) => energy,
548                        Err(err) => return Err(AnnealingError::IsingError(err)),
549                    };
550
551                    let delta_e = new_energy - current_energy;
552
553                    // Metropolis acceptance criterion
554                    let accept = delta_e <= 0.0 || {
555                        let p = (-delta_e / temperature).exp();
556                        rng.random_range(0.0..1.0) < p
557                    };
558
559                    // Apply the spin flip if accepted
560                    if accept {
561                        current_energy = new_energy;
562                    } else {
563                        // Revert the change
564                        spins[qubit] = current_spin;
565                    }
566                }
567
568                // Increment sweep counter
569                total_sweeps += 1;
570            }
571
572            // Check if this is a better solution
573            if current_energy < best_energy {
574                best_energy = current_energy;
575                best_spins = spins.clone();
576            }
577
578            // Increment repetition counter
579            completed_repetitions += 1;
580        }
581
582        // Calculate runtime
583        let runtime = start_time.elapsed();
584
585        // Build result
586        Ok(AnnealingSolution {
587            best_spins,
588            best_energy,
589            repetitions: completed_repetitions,
590            total_sweeps,
591            runtime,
592            info: format!(
593                "Performed {completed_repetitions} repetitions with {total_sweeps} total sweeps in {runtime:?}"
594            ),
595        })
596    }
597}
598
599#[cfg(test)]
600mod tests {
601    use super::*;
602    #[allow(unused_imports)]
603    use crate::ising::QuboModel;
604
605    #[test]
606    fn test_annealing_params() {
607        // Create default parameters
608        let params = AnnealingParams::default();
609
610        // Validate parameters
611        assert!(params.validate().is_ok());
612
613        // Test invalid parameters
614        let mut invalid_params = params.clone();
615        invalid_params.initial_temperature = 0.0;
616        assert!(invalid_params.validate().is_err());
617
618        invalid_params = params.clone();
619        invalid_params.num_sweeps = 0;
620        assert!(invalid_params.validate().is_err());
621    }
622
623    #[test]
624    fn test_classical_annealing_simple() {
625        // Create a simple 2-qubit ferromagnetic Ising model
626        let mut model = IsingModel::new(2);
627        model
628            .set_coupling(0, 1, -1.0)
629            .expect("Failed to set coupling"); // Ferromagnetic coupling
630
631        // Create annealing simulator with fixed seed for reproducibility
632        let mut params = AnnealingParams::default();
633        params.seed = Some(42);
634        params.num_sweeps = 100;
635        params.num_repetitions = 5;
636
637        let simulator =
638            ClassicalAnnealingSimulator::new(params).expect("Failed to create simulator");
639
640        // Solve the model
641        let result = simulator.solve(&model).expect("Failed to solve model");
642
643        // Check that we found the ground state (all spins aligned)
644        assert_eq!(result.best_spins.len(), 2);
645        assert!(
646            (result.best_spins[0] == 1 && result.best_spins[1] == 1)
647                || (result.best_spins[0] == -1 && result.best_spins[1] == -1)
648        );
649
650        // Check energy
651        assert_eq!(result.best_energy, -1.0);
652    }
653
654    #[test]
655    fn test_quantum_annealing_simple() {
656        // Create a simple 2-qubit ferromagnetic Ising model
657        let mut model = IsingModel::new(2);
658        model
659            .set_coupling(0, 1, -1.0)
660            .expect("Failed to set coupling"); // Ferromagnetic coupling
661
662        // Create annealing simulator with fixed seed for reproducibility
663        let mut params = AnnealingParams::default();
664        params.seed = Some(42);
665        params.num_sweeps = 100;
666        params.num_repetitions = 5;
667        params.trotter_slices = 10;
668
669        let simulator =
670            QuantumAnnealingSimulator::new(params).expect("Failed to create quantum simulator");
671
672        // Solve the model
673        let result = simulator.solve(&model).expect("Failed to solve model");
674
675        // Check that we found the ground state (all spins aligned)
676        assert_eq!(result.best_spins.len(), 2);
677        assert!(
678            (result.best_spins[0] == 1 && result.best_spins[1] == 1)
679                || (result.best_spins[0] == -1 && result.best_spins[1] == -1)
680        );
681
682        // Check energy
683        assert_eq!(result.best_energy, -1.0);
684    }
685
686    #[test]
687    fn test_classical_annealing_frustrated() {
688        // Create a 3-qubit frustrated Ising model
689        let mut model = IsingModel::new(3);
690        model
691            .set_coupling(0, 1, -1.0)
692            .expect("Failed to set coupling"); // Ferromagnetic coupling
693        model
694            .set_coupling(1, 2, -1.0)
695            .expect("Failed to set coupling"); // Ferromagnetic coupling
696        model
697            .set_coupling(0, 2, 1.0)
698            .expect("Failed to set coupling"); // Antiferromagnetic coupling
699
700        // Create annealing simulator with fixed seed for reproducibility
701        let mut params = AnnealingParams::default();
702        params.seed = Some(42);
703        params.num_sweeps = 200;
704        params.num_repetitions = 10;
705
706        let simulator =
707            ClassicalAnnealingSimulator::new(params).expect("Failed to create simulator");
708
709        // Solve the model
710        let result = simulator.solve(&model).expect("Failed to solve model");
711
712        // Check energy (should be -1.0 for the ground state)
713        assert!(result.best_energy <= -1.0);
714    }
715}
quantrs2_anneal/simulator.rs

quantrs2_anneal/
simulator.rs
