1use crate::lattice::Lattice;
2use crate::stats::Statistics;
3use crate::{clusters, energy, sweep, tempering};
4use rand::{Rng, SeedableRng};
5use rand_xoshiro::Xoshiro256StarStar;
6
7pub struct Realization {
8 pub couplings: Vec<f32>,
9 pub spins: Vec<i8>,
10 pub temperatures: Vec<f32>,
11 pub system_ids: Vec<usize>,
12 pub rngs: Vec<Xoshiro256StarStar>,
13 pub energies: Vec<f32>,
14 pub interactions: Vec<f32>,
15}
16
17impl Realization {
18 pub fn new(
19 lattice: &Lattice,
20 couplings: Vec<f32>,
21 temps: &[f32],
22 n_replicas: usize,
23 base_seed: u64,
24 ) -> Self {
25 let n_spins = lattice.n_spins;
26 let n_temps = temps.len();
27 let n_systems = n_replicas * n_temps;
28
29 let temperatures = temps.repeat(n_replicas);
30
31 let mut rngs = Vec::with_capacity(n_systems);
32 for i in 0..n_systems {
33 rngs.push(Xoshiro256StarStar::seed_from_u64(base_seed + i as u64));
34 }
35
36 let mut spins = vec![0i8; n_systems * n_spins];
37 for (i, rng) in rngs.iter_mut().enumerate() {
38 for j in 0..n_spins {
39 spins[i * n_spins + j] = if rng.gen::<f32>() < 0.5 { -1 } else { 1 };
40 }
41 }
42
43 let system_ids: Vec<usize> = (0..n_systems).collect();
44
45 let (energies, interactions) =
46 energy::compute_energies(lattice, &spins, &couplings, n_systems, true);
47 let interactions = interactions.unwrap();
48
49 Self {
50 couplings,
51 spins,
52 temperatures,
53 system_ids,
54 rngs,
55 energies,
56 interactions,
57 }
58 }
59
60 pub fn reset(&mut self, lattice: &Lattice, n_replicas: usize, n_temps: usize, base_seed: u64) {
61 let n_spins = lattice.n_spins;
62 let n_systems = n_replicas * n_temps;
63
64 for i in 0..n_systems {
65 self.rngs[i] = Xoshiro256StarStar::seed_from_u64(base_seed + i as u64);
66 for j in 0..n_spins {
67 self.spins[i * n_spins + j] = if self.rngs[i].gen::<f32>() < 0.5 {
68 -1
69 } else {
70 1
71 };
72 }
73 }
74
75 self.system_ids = (0..n_systems).collect();
76
77 let (energies, interactions) =
78 energy::compute_energies(lattice, &self.spins, &self.couplings, n_systems, true);
79 self.energies = energies;
80 self.interactions = interactions.unwrap();
81 }
82}
83
84pub struct SweepResult {
85 pub mags: Vec<f64>,
86 pub mags2: Vec<f64>,
87 pub mags4: Vec<f64>,
88 pub energies: Vec<f64>,
89 pub energies2: Vec<f64>,
90 pub overlap: Vec<f64>,
91 pub overlap2: Vec<f64>,
92 pub overlap4: Vec<f64>,
93}
94
95#[allow(clippy::too_many_arguments)]
96pub fn run_sweep_loop(
97 lattice: &Lattice,
98 real: &mut Realization,
99 n_replicas: usize,
100 n_temps: usize,
101 n_sweeps: usize,
102 warmup_sweeps: usize,
103 sweep_mode: &str,
104 cluster_update_interval: Option<usize>,
105 cluster_mode: &str,
106 pt_interval: Option<usize>,
107 houdayer_interval: Option<usize>,
108 on_sweep: &(dyn Fn() + Sync),
109) -> SweepResult {
110 let n_spins = lattice.n_spins;
111 let n_systems = n_replicas * n_temps;
112
113 let mut mags_stat = Statistics::new(n_temps, 1);
114 let mut mags2_stat = Statistics::new(n_temps, 1);
115 let mut mags4_stat = Statistics::new(n_temps, 1);
116 let mut energies_stat = Statistics::new(n_temps, 1);
117 let mut energies2_stat = Statistics::new(n_temps, 2);
118
119 let n_pairs = n_replicas / 2;
120 let mut overlap_stat = Statistics::new(n_temps, 1);
121 let mut overlap2_stat = Statistics::new(n_temps, 1);
122 let mut overlap4_stat = Statistics::new(n_temps, 1);
123
124 for sweep_id in 0..n_sweeps {
125 on_sweep();
126 let record = sweep_id >= warmup_sweeps;
127
128 match sweep_mode {
129 "metropolis" => sweep::metropolis_sweep(
130 lattice,
131 &mut real.spins,
132 &real.couplings,
133 &real.temperatures,
134 &real.system_ids,
135 &mut real.rngs,
136 ),
137 "gibbs" => sweep::gibbs_sweep(
138 lattice,
139 &mut real.spins,
140 &real.couplings,
141 &real.temperatures,
142 &real.system_ids,
143 &mut real.rngs,
144 ),
145 _ => unreachable!(),
146 }
147
148 let do_cluster = cluster_update_interval.is_some_and(|interval| sweep_id % interval == 0);
149
150 if do_cluster {
151 match cluster_mode {
152 "wolff" => {
153 clusters::wolff_update(
154 lattice,
155 &mut real.spins,
156 &real.couplings,
157 &real.temperatures,
158 &real.system_ids,
159 &mut real.rngs,
160 );
161 (real.energies, _) = energy::compute_energies(
162 lattice,
163 &real.spins,
164 &real.couplings,
165 n_systems,
166 false,
167 );
168 }
169 "sw" => {
170 let (energies, interactions) = energy::compute_energies(
171 lattice,
172 &real.spins,
173 &real.couplings,
174 n_systems,
175 true,
176 );
177 real.energies = energies;
178 real.interactions = interactions.unwrap();
179
180 clusters::sw_update(
181 lattice,
182 &mut real.spins,
183 &real.interactions,
184 &real.temperatures,
185 &real.system_ids,
186 &mut real.rngs,
187 );
188
189 (real.energies, _) = energy::compute_energies(
190 lattice,
191 &real.spins,
192 &real.couplings,
193 n_systems,
194 false,
195 );
196 }
197 _ => unreachable!(),
198 }
199 } else {
200 (real.energies, _) =
201 energy::compute_energies(lattice, &real.spins, &real.couplings, n_systems, false);
202 }
203
204 if record {
205 let mut mags = vec![0.0f32; n_temps];
206 let mut mags2 = vec![0.0f32; n_temps];
207 let mut mags4 = vec![0.0f32; n_temps];
208 let mut energies_ordered = vec![0.0f32; n_temps];
209
210 for r in 0..n_replicas {
211 let offset = r * n_temps;
212 for t in 0..n_temps {
213 let system_id = real.system_ids[offset + t];
214 let spin_base = system_id * n_spins;
215 let mut sum = 0i64;
216 for j in 0..n_spins {
217 sum += real.spins[spin_base + j] as i64;
218 }
219 let mag = sum as f32 / n_spins as f32;
220 let m2 = mag * mag;
221 mags[t] = mag;
222 mags2[t] = m2;
223 mags4[t] = m2 * m2;
224 energies_ordered[t] = real.energies[system_id];
225 }
226
227 mags_stat.update(&mags);
228 mags2_stat.update(&mags2);
229 mags4_stat.update(&mags4);
230 energies_stat.update(&energies_ordered);
231 energies2_stat.update(&energies_ordered);
232 }
233
234 for pair_idx in 0..n_pairs {
235 let r_a = 2 * pair_idx;
236 let r_b = 2 * pair_idx + 1;
237 let mut overlaps = vec![0.0f32; n_temps];
238 let mut overlaps2 = vec![0.0f32; n_temps];
239 let mut overlaps4 = vec![0.0f32; n_temps];
240
241 for t in 0..n_temps {
242 let sys_a = real.system_ids[r_a * n_temps + t];
243 let sys_b = real.system_ids[r_b * n_temps + t];
244 let base_a = sys_a * n_spins;
245 let base_b = sys_b * n_spins;
246 let mut dot = 0i64;
247 for j in 0..n_spins {
248 dot += (real.spins[base_a + j] as i64) * (real.spins[base_b + j] as i64);
249 }
250 let q = dot as f32 / n_spins as f32;
251 let q2 = q * q;
252 overlaps[t] = q;
253 overlaps2[t] = q2;
254 overlaps4[t] = q2 * q2;
255 }
256
257 overlap_stat.update(&overlaps);
258 overlap2_stat.update(&overlaps2);
259 overlap4_stat.update(&overlaps4);
260 }
261 }
262
263 if let Some(interval) = pt_interval {
264 if sweep_id % interval == 0 {
265 for r in 0..n_replicas {
266 let offset = r * n_temps;
267 let sid_slice = &mut real.system_ids[offset..offset + n_temps];
268 let temp_slice = &real.temperatures[offset..offset + n_temps];
269 tempering::parallel_tempering(
270 &real.energies,
271 temp_slice,
272 sid_slice,
273 n_spins,
274 &mut real.rngs[offset],
275 );
276 }
277 }
278 }
279
280 if let Some(interval) = houdayer_interval {
281 if sweep_id % interval == 0 && n_replicas >= 2 {
282 clusters::houdayer_update(
283 lattice,
284 &mut real.spins,
285 &real.system_ids,
286 n_replicas,
287 n_temps,
288 &mut real.rngs[0],
289 );
290 (real.energies, _) = energy::compute_energies(
291 lattice,
292 &real.spins,
293 &real.couplings,
294 n_systems,
295 false,
296 );
297 }
298 }
299 }
300
301 SweepResult {
302 mags: mags_stat.average(),
303 mags2: mags2_stat.average(),
304 mags4: mags4_stat.average(),
305 energies: energies_stat.average(),
306 energies2: energies2_stat.average(),
307 overlap: if n_pairs > 0 {
308 overlap_stat.average()
309 } else {
310 vec![]
311 },
312 overlap2: if n_pairs > 0 {
313 overlap2_stat.average()
314 } else {
315 vec![]
316 },
317 overlap4: if n_pairs > 0 {
318 overlap4_stat.average()
319 } else {
320 vec![]
321 },
322 }
323}
324
325pub fn aggregate_results(results: &[SweepResult]) -> SweepResult {
326 let n = results.len() as f64;
327 let n_temps = results[0].mags.len();
328 let n_overlap = results[0].overlap.len();
329
330 let mut agg = SweepResult {
331 mags: vec![0.0; n_temps],
332 mags2: vec![0.0; n_temps],
333 mags4: vec![0.0; n_temps],
334 energies: vec![0.0; n_temps],
335 energies2: vec![0.0; n_temps],
336 overlap: vec![0.0; n_overlap],
337 overlap2: vec![0.0; n_overlap],
338 overlap4: vec![0.0; n_overlap],
339 };
340
341 for r in results {
342 for (a, &v) in agg.mags.iter_mut().zip(r.mags.iter()) {
343 *a += v;
344 }
345 for (a, &v) in agg.mags2.iter_mut().zip(r.mags2.iter()) {
346 *a += v;
347 }
348 for (a, &v) in agg.mags4.iter_mut().zip(r.mags4.iter()) {
349 *a += v;
350 }
351 for (a, &v) in agg.energies.iter_mut().zip(r.energies.iter()) {
352 *a += v;
353 }
354 for (a, &v) in agg.energies2.iter_mut().zip(r.energies2.iter()) {
355 *a += v;
356 }
357 for (a, &v) in agg.overlap.iter_mut().zip(r.overlap.iter()) {
358 *a += v;
359 }
360 for (a, &v) in agg.overlap2.iter_mut().zip(r.overlap2.iter()) {
361 *a += v;
362 }
363 for (a, &v) in agg.overlap4.iter_mut().zip(r.overlap4.iter()) {
364 *a += v;
365 }
366 }
367
368 for v in agg
369 .mags
370 .iter_mut()
371 .chain(agg.mags2.iter_mut())
372 .chain(agg.mags4.iter_mut())
373 .chain(agg.energies.iter_mut())
374 .chain(agg.energies2.iter_mut())
375 .chain(agg.overlap.iter_mut())
376 .chain(agg.overlap2.iter_mut())
377 .chain(agg.overlap4.iter_mut())
378 {
379 *v /= n;
380 }
381
382 agg
383}