lattice_qcd_rs/simulation/monte_carlo/
overrelaxation.rs

1//! Overrelaxation method
2//!
3//! The goal of the overrelaxation is to move thought the phase space as much as possible but conserving the hamiltonian.
4//! It can be used to improve the speed of thermalisation by vising more states.
5//! Alone it can't advance the simulation as it preserved the hamiltonian. You need to use other method with this one.
6//! You can look at [`super::HybridMethodVec`] and [`super::HybridMethodCouple`].
7//!
8//! In my limited experience [`OverrelaxationSweepReverse`] moves a bit more though the phase space than [`OverrelaxationSweepRotation`].
9//! The difference is slight though.
10//!
11//! # Example
12//! ```
13//! # use std::error::Error;
14//! #
15//! # fn main() -> Result<(), Box<dyn Error>> {
16//! use lattice_qcd_rs::simulation::{HeatBathSweep, LatticeState, LatticeStateDefault, OverrelaxationSweepReverse};
17//! use rand::SeedableRng;
18//!
19//! let rng = rand::rngs::StdRng::seed_from_u64(0); // change with your seed
20//! let mut heat_bath = HeatBathSweep::new(rng);
21//! let mut overrelax = OverrelaxationSweepReverse::default();
22//!
23//! let mut state = LatticeStateDefault::<3>::new_cold(1_f64, 8_f64, 4)?; // 1_f64 : size, 8_f64: beta, 4 number of points.
24//! for _ in 0..2 {
25//!     state = state.monte_carlo_step(&mut heat_bath)?;
26//!     state = state.monte_carlo_step(&mut overrelax)?;
27//!     // operation to track the progress or the evolution
28//! }
29//! // operation at the end of the simulation
30//! #     Ok(())
31//! # }
32//! ```
33
34#[cfg(feature = "serde-serialize")]
35use serde::{Deserialize, Serialize};
36
37use super::{
38    super::{
39        super::{error::Never, lattice::LatticeLinkCanonical, su3, Complex},
40        state::{LatticeState, LatticeStateDefault},
41    },
42    staple, MonteCarlo,
43};
44
45/// Overrelaxation algorithm using rotation method.
46///
47/// Alone it can't advance the simulation as it preserved the hamiltonian.
48/// You need to use other method with this one.
49/// You can look at [`super::HybridMethodVec`] and [`super::HybridMethodCouple`].
50///
51/// The algorithm is described <https://arxiv.org/abs/hep-lat/0503041> in section 2.1 up to step 2 using `\hat X_{NN}`.
52///
53/// # Example
54/// ```
55/// # use std::error::Error;
56/// #
57/// # fn main() -> Result<(), Box<dyn Error>> {
58/// use lattice_qcd_rs::simulation::{HeatBathSweep, LatticeState, LatticeStateDefault, OverrelaxationSweepRotation};
59/// use rand::SeedableRng;
60///
61/// let rng = rand::rngs::StdRng::seed_from_u64(0); // change with your seed
62/// let mut heat_bath = HeatBathSweep::new(rng);
63/// let mut overrelax = OverrelaxationSweepRotation::default();
64///
65/// let mut state = LatticeStateDefault::<3>::new_cold(1_f64, 8_f64, 4)?; // 1_f64 : size, 8_f64: beta, 4 number of points.
66/// for _ in 0..2 {
67///     state = state.monte_carlo_step(&mut heat_bath)?;
68///     state = state.monte_carlo_step(&mut overrelax)?;
69///     // operation to track the progress or the evolution
70/// }
71/// // operation at the end of the simulation
72/// #     Ok(())
73/// # }
74/// ```
75#[derive(Clone, Copy, Debug, PartialEq, Eq, Hash, PartialOrd, Ord)]
76#[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))]
77pub struct OverrelaxationSweepRotation;
78
79impl OverrelaxationSweepRotation {
80    /// Create a new Self with an given RNG
81    pub const fn new() -> Self {
82        Self {}
83    }
84
85    #[inline]
86    fn get_modif<const D: usize>(
87        state: &LatticeStateDefault<D>,
88        link: &LatticeLinkCanonical<D>,
89    ) -> na::Matrix3<Complex> {
90        let link_matrix = state
91            .link_matrix()
92            .matrix(&link.into(), state.lattice())
93            .unwrap();
94        let a = staple(state.link_matrix(), state.lattice(), link).adjoint();
95        let svd = na::SVD::<Complex, na::U3, na::U3>::new(a, true, true);
96        let rot = svd.u.unwrap() * svd.v_t.unwrap();
97        rot * link_matrix.adjoint() * rot
98    }
99
100    #[inline]
101    fn next_element_default<const D: usize>(
102        mut state: LatticeStateDefault<D>,
103    ) -> LatticeStateDefault<D> {
104        let lattice = state.lattice().clone();
105        lattice.get_links().for_each(|link| {
106            let potential_modif = Self::get_modif(&state, &link);
107            *state.link_mut(&link).unwrap() = potential_modif;
108        });
109        state
110    }
111}
112
113impl std::fmt::Display for OverrelaxationSweepRotation {
114    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
115        write!(f, "overrelaxation method by rotation")
116    }
117}
118
119impl Default for OverrelaxationSweepRotation {
120    fn default() -> Self {
121        Self::new()
122    }
123}
124
125impl<const D: usize> MonteCarlo<LatticeStateDefault<D>, D> for OverrelaxationSweepRotation {
126    type Error = Never;
127
128    #[inline]
129    fn next_element(
130        &mut self,
131        state: LatticeStateDefault<D>,
132    ) -> Result<LatticeStateDefault<D>, Self::Error> {
133        Ok(Self::next_element_default(state))
134    }
135}
136
137/// Overrelaxation algorithm using the reverse method.
138///
139/// Alone it can't advance the simulation as it preserved the hamiltonian.
140/// You need to use other method with this one.
141/// You can look at [`super::HybridMethodVec`] and [`super::HybridMethodCouple`].
142///
143/// The algorithm is described in <https://doi.org/10.1016/0370-2693(90)90032-2>.
144///
145/// # Example
146/// see level module documentation [`super::overrelaxation`]
147#[derive(Clone, Copy, Debug, PartialEq, Eq, Hash, PartialOrd, Ord)]
148#[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))]
149pub struct OverrelaxationSweepReverse;
150
151impl OverrelaxationSweepReverse {
152    /// Create a new Self with an given RNG
153    pub const fn new() -> Self {
154        Self {}
155    }
156
157    #[inline]
158    fn get_modif<const D: usize>(
159        state: &LatticeStateDefault<D>,
160        link: &LatticeLinkCanonical<D>,
161    ) -> na::Matrix3<Complex> {
162        let link_matrix = state
163            .link_matrix()
164            .matrix(&link.into(), state.lattice())
165            .unwrap();
166        let a = staple(state.link_matrix(), state.lattice(), link).adjoint();
167        let svd = na::SVD::<Complex, na::U3, na::U3>::new(a, true, true);
168        svd.u.unwrap()
169            * su3::reverse(svd.u.unwrap().adjoint() * link_matrix * svd.v_t.unwrap().adjoint())
170            * svd.v_t.unwrap()
171    }
172
173    #[inline]
174    fn next_element_default<const D: usize>(
175        mut state: LatticeStateDefault<D>,
176    ) -> LatticeStateDefault<D> {
177        let lattice = state.lattice().clone();
178        lattice.get_links().for_each(|link| {
179            let potential_modif = Self::get_modif(&state, &link);
180            *state.link_mut(&link).unwrap() = potential_modif;
181        });
182        state
183    }
184}
185
186impl Default for OverrelaxationSweepReverse {
187    fn default() -> Self {
188        Self::new()
189    }
190}
191
192impl std::fmt::Display for OverrelaxationSweepReverse {
193    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
194        write!(f, "overrelaxation method by reverse")
195    }
196}
197
198impl<const D: usize> MonteCarlo<LatticeStateDefault<D>, D> for OverrelaxationSweepReverse {
199    type Error = Never;
200
201    #[inline]
202    fn next_element(
203        &mut self,
204        state: LatticeStateDefault<D>,
205    ) -> Result<LatticeStateDefault<D>, Self::Error> {
206        Ok(Self::next_element_default(state))
207    }
208}
209
210#[cfg(test)]
211mod test {
212    use rand::SeedableRng;
213
214    use super::super::super::state::{LatticeState, LatticeStateDefault};
215    use super::super::MonteCarlo;
216    use super::*;
217
218    const SEED_RNG: u64 = 0x45_78_93_f4_4a_b0_67_f0;
219
220    fn test_same_energy<MC>(mc: &mut MC, rng: &mut impl rand::Rng)
221    where
222        MC: MonteCarlo<LatticeStateDefault<3>, 3>,
223        MC::Error: core::fmt::Debug,
224    {
225        let state = LatticeStateDefault::<3>::new_determinist(1_f64, 1_f64, 4, rng).unwrap();
226        let h = state.hamiltonian_links();
227        let state2 = mc.next_element(state).unwrap();
228        let h2 = state2.hamiltonian_links();
229        println!("h1 {}, h2 {}", h, h2);
230        // Relative assert : we need to multi by the mean value of h
231        // TODO use crate approx ?
232        assert!((h - h2).abs() < f64::EPSILON * 100_f64 * 4_f64.powi(3) * (h + h2) * 0.5_f64);
233    }
234
235    /// Here we test that OverrelaxationSweepReverse conserve the energy.
236    #[test]
237    fn same_energy_reverse() {
238        let mut rng = rand::rngs::StdRng::seed_from_u64(SEED_RNG);
239        let mut overrelax = OverrelaxationSweepReverse::new();
240        for _ in 0_u32..10_u32 {
241            test_same_energy(&mut overrelax, &mut rng);
242        }
243    }
244
245    /// Here we test that OverrelaxationSweepRotation conserve the energy.
246    #[test]
247    fn same_energy_rotation() {
248        let mut rng = rand::rngs::StdRng::seed_from_u64(SEED_RNG);
249        let mut overrelax = OverrelaxationSweepRotation::new();
250        for _ in 0_u32..10_u32 {
251            test_same_energy(&mut overrelax, &mut rng);
252        }
253    }
254}
lattice_qcd_rs/simulation/monte_carlo/overrelaxation.rs

lattice_qcd_rs/simulation/monte_carlo/
overrelaxation.rs
