1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
use rand::{
distributions::{Bernoulli, Distribution},
rngs::SmallRng,
};
use crate::{
algorithm::Algorithm,
lattice::{Lattice, Spin},
};
/// The Wolff cluster Ising algorithm.
pub struct Wolff {
/// Probability of adding a spin to the cluster.
addition_random_dist: Bernoulli,
/// Vector to store inverted spins in the while loop.
inverted_spins: Vec<Spin>,
}
impl Wolff {
#[must_use]
pub fn new(lattice: &Lattice, temperature: f64) -> Self {
let addition_probability = Self::addition_p_from_temperature(temperature);
let side_length = lattice.side_length();
let n_spins = side_length * side_length;
Self {
// unwrap: Addition probability is between 0 and 1 for positive temperatures.
addition_random_dist: Bernoulli::new(addition_probability).unwrap(),
// One spin is subtracted because the seed spin is popped directly and never
// added afterwards in the while loop. One can surely allocate less, but it
// is not an easy problem to think about. The maximum will stay O(L^2), only
// the constant of L^2 can be minimized. But the risk is not worth it, since
// the capacity is duplicated after a push beyond the capacity.
inverted_spins: Vec::with_capacity(n_spins - 1),
}
}
// Probability to add a spin to the cluster that satisfies detailed balance.
//
// Between 0 and 1 for positive temperatures.
fn addition_p_from_temperature(temperature: f64) -> f64 {
1.0 - (-2.0 / temperature).exp()
}
}
impl Algorithm for Wolff {
// Updates the temperature dependent addition probability.
fn set_temperature(&mut self, val: f64) {
let p = Self::addition_p_from_temperature(val);
// unwrap: Addition probability is between 0 and 1 for positive temperatures.
self.addition_random_dist = Bernoulli::new(p).unwrap();
}
// Builds and flips a cluster.
fn step(&mut self, lattice: &mut Lattice, rng: &mut SmallRng) {
// Get a random seed spin.
let (seed_spin, old_seed_spin_value) = lattice.rand_spin_with_value(rng);
let new_seed_spin_value = !old_seed_spin_value;
// Invert the seed spin.
// Safety: Random spin is always in bounds for the same lattice.
unsafe {
lattice.set_value_unchecked(&seed_spin, new_seed_spin_value);
}
self.inverted_spins.push(seed_spin);
// Use a stack to prevent recursion.
while let Some(center_spin) = self.inverted_spins.pop() {
// Get neighbouring spins.
let neighbouring_spins = lattice.neighbouring_spins(¢er_spin);
for neighbouring_spin in neighbouring_spins {
// Safety: Neighbouring spins of a valid spin are granted to be inside the lattice.
let spin_value = unsafe { lattice.value_unchecked(&neighbouring_spin) };
// Invert neighbouring spins with the old seed spin value with the addition probability.
if spin_value == old_seed_spin_value && self.addition_random_dist.sample(rng) {
// Safety: Neighbouring spin of a valid spin.
unsafe { lattice.set_value_unchecked(&neighbouring_spin, new_seed_spin_value) };
self.inverted_spins.push(neighbouring_spin);
}
}
}
}
}