lattice_qcd_rs/simulation/monte_carlo/
mod.rs

1//! Module for Monte-Carlo algrorithme, see the trait [`MonteCarlo`].
2//!
3//! This is one of the way to carry out simulation. This work by taking a state and progressively changing it (most of the time randomly).
4//!
5//! # Examples
6//! see [`MetropolisHastingsSweep`], [`HeatBathSweep`], [`overrelaxation`] etc...
7
8use std::marker::PhantomData;
9
10use na::ComplexField;
11use rand_distr::Distribution;
12#[cfg(feature = "serde-serialize")]
13use serde::{Deserialize, Serialize};
14
15use super::state::{LatticeState, LatticeStateDefault};
16use crate::{
17    field::LinkMatrix,
18    lattice::{Direction, LatticeCyclic, LatticeLink, LatticeLinkCanonical},
19    Complex, Real,
20};
21
22pub mod heat_bath;
23pub mod hybrid;
24pub mod hybrid_monte_carlo;
25pub mod metropolis_hastings;
26pub mod metropolis_hastings_sweep;
27pub mod overrelaxation;
28
29pub use heat_bath::*;
30pub use hybrid::*;
31pub use hybrid_monte_carlo::*;
32pub use metropolis_hastings::*;
33pub use metropolis_hastings_sweep::*;
34pub use overrelaxation::*;
35
36/// Monte-Carlo algorithm, giving the next element in the simulation.
37/// It is also a Markov chain.
38///
39/// # Example
40/// ```
41/// # use std::error::Error;
42/// #
43/// # fn main() -> Result<(), Box<dyn Error>> {
44/// use lattice_qcd_rs::error::ImplementationError;
45/// use lattice_qcd_rs::simulation::{
46///     LatticeState, LatticeStateDefault, MetropolisHastingsSweep, MonteCarlo,
47/// };
48/// use rand::SeedableRng;
49///
50/// let rng = rand::rngs::StdRng::seed_from_u64(0); // change with your seed
51/// let mut mh = MetropolisHastingsSweep::new(1, 0.1_f64, rng)
52///     .ok_or(ImplementationError::OptionWithUnexpectedNone)?;
53/// // Realistically you want more steps than 10
54///
55/// let mut state = LatticeStateDefault::<3>::new_cold(1_f64, 6_f64, 4)?;
56/// for _ in 0..10 {
57///     state = mh.next_element(state)?;
58///     // or state.monte_carlo_step(&mut hmc)?;
59///     // operation to track the progress or the evolution
60/// }
61/// // operation at the end of the simulation
62/// #     Ok(())
63/// # }
64/// ```
65pub trait MonteCarlo<State, const D: usize>
66where
67    State: LatticeState<D>,
68{
69    /// Error returned while getting the next element.
70    type Error;
71
72    /// Do one Monte Carlo simulation step.
73    ///
74    /// # Errors
75    /// Return an error if the simulation failed
76    fn next_element(&mut self, state: State) -> Result<State, Self::Error>;
77}
78
79/// Some times it is easier to just implement a potential next element, the rest is done automatically using an [`McWrapper`].
80///
81/// To get an [`MonteCarlo`] use the wrapper [`McWrapper`].
82/// # Example
83/// ```
84/// # use std::error::Error;
85/// #
86/// # fn main() -> Result<(), Box<dyn Error>> {
87/// use lattice_qcd_rs::error::ImplementationError;
88/// use lattice_qcd_rs::simulation::{
89///     LatticeState, LatticeStateDefault, McWrapper, MetropolisHastingsDiagnostic, MonteCarlo,
90/// };
91/// use rand::SeedableRng;
92///
93/// let rng = rand::rngs::StdRng::seed_from_u64(0); // change with your seed
94/// let mh = MetropolisHastingsDiagnostic::new(1, 0.1_f64)
95///     .ok_or(ImplementationError::OptionWithUnexpectedNone)?;
96/// let mut wrapper = McWrapper::new(mh, rng);
97///
98/// // Realistically you want more steps than 10
99///
100/// let mut state = LatticeStateDefault::<3>::new_cold(1_f64, 6_f64, 4)?;
101/// for _ in 0..100 {
102///     state = state.monte_carlo_step(&mut wrapper)?;
103///     println!(
104///         "probability of acceptance during last step {}, does it accepted the change ? {}",
105///         mh.prob_replace_last(),
106///         mh.has_replace_last()
107///     );
108///     // or state.monte_carlo_step(&mut wrapper)?;
109///     // operation to track the progress or the evolution
110/// }
111/// // operation at the end of the simulation
112/// let (_, rng) = wrapper.deconstruct(); // get the rng back
113///                                       // continue with further operation using the same rng ...
114/// #     Ok(())
115/// # }
116/// ```
117pub trait MonteCarloDefault<State, const D: usize>
118where
119    State: LatticeState<D>,
120{
121    /// Error returned while getting the next element.
122    type Error;
123
124    /// Generate a radom element from the previous element (like a Markov chain).
125    ///
126    /// # Errors
127    /// Gives an error if a potential next element cannot be generated.
128    fn potential_next_element<Rng>(
129        &mut self,
130        state: &State,
131        rng: &mut Rng,
132    ) -> Result<State, Self::Error>
133    where
134        Rng: rand::Rng + ?Sized;
135
136    /// probability of the next element to replace the current one.
137    ///
138    /// by default it is Exp(-H_old) / Exp(-H_new).
139    fn probability_of_replacement(old_state: &State, new_state: &State) -> Real {
140        (old_state.hamiltonian_links() - new_state.hamiltonian_links())
141            .exp()
142            .min(1_f64)
143            .max(0_f64)
144    }
145
146    /// Get the next element in the chain either the old state or a new one replacing it.
147    ///
148    /// # Errors
149    /// Gives an error if a potential next element cannot be generated.
150    fn next_element_default<Rng>(
151        &mut self,
152        state: State,
153        rng: &mut Rng,
154    ) -> Result<State, Self::Error>
155    where
156        Rng: rand::Rng + ?Sized,
157    {
158        let potential_next = self.potential_next_element(&state, rng)?;
159        let proba = Self::probability_of_replacement(&state, &potential_next)
160            .min(1_f64)
161            .max(0_f64);
162        let d = rand::distributions::Bernoulli::new(proba).unwrap();
163        if d.sample(rng) {
164            Ok(potential_next)
165        }
166        else {
167            Ok(state)
168        }
169    }
170}
171
172/// A wrapper used to implement [`MonteCarlo`] from a [`MonteCarloDefault`]
173///
174/// # Example
175/// ```
176/// # use std::error::Error;
177/// #
178/// # fn main() -> Result<(), Box<dyn Error>> {
179/// use lattice_qcd_rs::error::ImplementationError;
180/// use lattice_qcd_rs::simulation::{
181///     LatticeState, LatticeStateDefault, McWrapper, MetropolisHastingsDiagnostic, MonteCarlo,
182/// };
183/// use rand::SeedableRng;
184///
185/// let rng = rand::rngs::StdRng::seed_from_u64(0); // change with your seed
186/// let mh = MetropolisHastingsDiagnostic::new(1, 0.1_f64)
187///     .ok_or(ImplementationError::OptionWithUnexpectedNone)?;
188/// let mut wrapper = McWrapper::new(mh, rng);
189///
190/// // Realistically you want more steps than 10
191///
192/// let mut state = LatticeStateDefault::<3>::new_cold(1_f64, 6_f64, 4)?;
193/// for _ in 0..100 {
194///     state = state.monte_carlo_step(&mut wrapper)?;
195///     println!(
196///         "probability of acceptance during last step {}, does it accepted the change ? {}",
197///         mh.prob_replace_last(),
198///         mh.has_replace_last()
199///     );
200///     // or state.monte_carlo_step(&mut wrapper)?;
201///     // operation to track the progress or the evolution
202/// }
203/// // operation at the end of the simulation
204/// let (_, rng) = wrapper.deconstruct(); // get the rng back
205/// #     Ok(())
206/// # }
207/// ```
208#[derive(Clone, Debug, PartialEq, Eq, Hash)]
209#[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))]
210pub struct McWrapper<MCD, State, Rng, const D: usize>
211where
212    MCD: MonteCarloDefault<State, D>,
213    State: LatticeState<D>,
214    Rng: rand::Rng,
215{
216    mcd: MCD,
217    rng: Rng,
218    _phantom: PhantomData<State>,
219}
220
221impl<MCD, State, Rng, const D: usize> McWrapper<MCD, State, Rng, D>
222where
223    MCD: MonteCarloDefault<State, D>,
224    State: LatticeState<D>,
225    Rng: rand::Rng,
226{
227    getter!(
228        /// Get a ref to the rng.
229        pub const,
230        rng,
231        Rng
232    );
233
234    /// Create the wrapper.
235    pub const fn new(mcd: MCD, rng: Rng) -> Self {
236        Self {
237            mcd,
238            rng,
239            _phantom: PhantomData,
240        }
241    }
242
243    /// deconstruct the structure to get back the rng if necessary
244    #[allow(clippy::missing_const_for_fn)] // false positive
245    pub fn deconstruct(self) -> (MCD, Rng) {
246        (self.mcd, self.rng)
247    }
248
249    /// Get a reference to the [`MonteCarloDefault`] inside the wrapper.
250    pub const fn mcd(&self) -> &MCD {
251        &self.mcd
252    }
253
254    /// Get a mutable reference to the rng
255    pub fn rng_mut(&mut self) -> &mut Rng {
256        &mut self.rng
257    }
258}
259
260impl<MCD, State, Rng, const D: usize> AsRef<Rng> for McWrapper<MCD, State, Rng, D>
261where
262    MCD: MonteCarloDefault<State, D>,
263    State: LatticeState<D>,
264    Rng: rand::Rng,
265{
266    fn as_ref(&self) -> &Rng {
267        self.rng()
268    }
269}
270
271impl<MCD, State, Rng, const D: usize> AsMut<Rng> for McWrapper<MCD, State, Rng, D>
272where
273    MCD: MonteCarloDefault<State, D>,
274    State: LatticeState<D>,
275    Rng: rand::Rng,
276{
277    fn as_mut(&mut self) -> &mut Rng {
278        self.rng_mut()
279    }
280}
281
282impl<T, State, Rng, const D: usize> MonteCarlo<State, D> for McWrapper<T, State, Rng, D>
283where
284    T: MonteCarloDefault<State, D>,
285    State: LatticeState<D>,
286    Rng: rand::Rng,
287{
288    type Error = T::Error;
289
290    fn next_element(&mut self, state: State) -> Result<State, Self::Error> {
291        self.mcd.next_element_default(state, &mut self.rng)
292    }
293}
294
295impl<MCD, State, Rng, const D: usize> Default for McWrapper<MCD, State, Rng, D>
296where
297    MCD: MonteCarloDefault<State, D> + Default,
298    State: LatticeState<D>,
299    Rng: rand::Rng + Default,
300{
301    fn default() -> Self {
302        Self::new(MCD::default(), Rng::default())
303    }
304}
305
306impl<MCD, State, Rng, const D: usize> std::fmt::Display for McWrapper<MCD, State, Rng, D>
307where
308    MCD: MonteCarloDefault<State, D> + std::fmt::Display,
309    State: LatticeState<D>,
310    Rng: rand::Rng + std::fmt::Display,
311{
312    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
313        write!(
314            f,
315            "Monte Carlo wrapper method {} with rng {}",
316            self.mcd(),
317            self.rng(),
318        )
319    }
320}
321
322/// Get the delta of energy by changing a link.
323#[inline]
324fn delta_s_old_new_cmp<const D: usize>(
325    link_matrix: &LinkMatrix,
326    lattice: &LatticeCyclic<D>,
327    link: &LatticeLinkCanonical<D>,
328    new_link: &na::Matrix3<Complex>,
329    beta: Real,
330    old_matrix: &na::Matrix3<Complex>,
331) -> Real {
332    let a = staple(link_matrix, lattice, link);
333    -((new_link - old_matrix) * a).trace().real() * beta / LatticeStateDefault::<D>::CA
334}
335
336// TODO move in state
337/// return the staple
338#[inline]
339fn staple<const D: usize>(
340    link_matrix: &LinkMatrix,
341    lattice: &LatticeCyclic<D>,
342    link: &LatticeLinkCanonical<D>,
343) -> na::Matrix3<Complex> {
344    let dir_j = link.dir();
345    Direction::<D>::positive_directions()
346        .iter()
347        .filter(|dir_i| *dir_i != dir_j)
348        .map(|dir_i| {
349            let el_1 = link_matrix
350                .sij(link.pos(), dir_j, dir_i, lattice)
351                .unwrap()
352                .adjoint();
353            let l_1 = LatticeLink::new(lattice.add_point_direction(*link.pos(), dir_j), -dir_i);
354            let u1 = link_matrix.matrix(&l_1, lattice).unwrap();
355            let l_2 = LatticeLink::new(lattice.add_point_direction(*link.pos(), &-dir_i), *dir_j);
356            let u2 = link_matrix.matrix(&l_2, lattice).unwrap().adjoint();
357            let l_3 = LatticeLink::new(lattice.add_point_direction(*link.pos(), &-dir_i), *dir_i);
358            let u3 = link_matrix.matrix(&l_3, lattice).unwrap();
359            el_1 + u1 * u2 * u3
360        })
361        .sum()
362}
lattice_qcd_rs/simulation/monte_carlo/mod.rs

lattice_qcd_rs/simulation/monte_carlo/
mod.rs
