Skip to main content

bench_xy/
bench_xy.rs

1//! XY sweep timing through the public `XySimulation` API; prints ms/sweep and a state
2//! checksum that behavior-neutral changes must leave unchanged.
3//!
4//! `cargo run --release -p spin-sim --example bench_xy`, configured by
5//! environment variables (defaults first):
6//! `PEAPODS_L` (64), `PEAPODS_DIMS` (2), `PEAPODS_N_TEMPS` (16),
7//! `PEAPODS_N_REPLICAS` (2), `PEAPODS_N_SWEEPS` (50), `PEAPODS_N_REALIZATIONS` (128),
8//! `PEAPODS_SWEEP_MODE` (metropolis | gibbs | none), `PEAPODS_CLUSTER_MODE`
9//! (none | sw | wolff), `PEAPODS_OVERRELAXATION` (0), `PEAPODS_PT_INTERVAL` (unset),
10//! `PEAPODS_WARMUP_RATIO` (0.25), `PEAPODS_COUPLINGS` (ferro | bimodal | gaussian),
11//! `PEAPODS_DILUTION` (vacancy probability, 0), `PEAPODS_SEQUENTIAL` (flag).
12//! Temperatures are geometric on [0.5, 2], spanning the 2D KT point.
13use std::env;
14use std::hash::{DefaultHasher, Hash, Hasher};
15use std::str::FromStr;
16use std::sync::atomic::AtomicBool;
17use std::time::Instant;
18
19use rand::{Rng, RngCore, SeedableRng};
20use rand_xoshiro::Xoshiro256StarStar;
21use spin_sim::config::*;
22use spin_sim::statistics::physics::PhysicsOptions;
23use spin_sim::{XyConfig, XySimulation};
24
25fn var<T: FromStr>(name: &str, default: T) -> T {
26    let Ok(value) = env::var(name) else {
27        return default;
28    };
29    value
30        .parse()
31        .unwrap_or_else(|_| panic!("{name}={value} is not a valid value"))
32}
33
34fn main() {
35    let l = var("PEAPODS_L", 64usize);
36    let dims = var("PEAPODS_DIMS", 2usize);
37    let n_temps = var("PEAPODS_N_TEMPS", 16usize);
38    let n_replicas = var("PEAPODS_N_REPLICAS", 2usize);
39    let n_sweeps = var("PEAPODS_N_SWEEPS", 50usize);
40    let n_realizations = var("PEAPODS_N_REALIZATIONS", 128usize);
41    let sweep_mode = var("PEAPODS_SWEEP_MODE", "metropolis".to_string());
42    let cluster_mode = var("PEAPODS_CLUSTER_MODE", "none".to_string());
43    let overrelaxation_sweeps = var("PEAPODS_OVERRELAXATION", 0usize);
44    let pt_interval = env::var("PEAPODS_PT_INTERVAL")
45        .ok()
46        .map(|v| v.parse::<usize>().expect("PEAPODS_PT_INTERVAL"));
47    let warmup_ratio = var("PEAPODS_WARMUP_RATIO", 0.25f64);
48    let couplings_mode = var("PEAPODS_COUPLINGS", "ferro".to_string());
49    let dilution = var("PEAPODS_DILUTION", 0.0f64);
50    let sequential = env::var_os("PEAPODS_SEQUENTIAL").is_some();
51
52    let shape = vec![l; dims];
53    let n_spins = l.pow(dims as u32);
54    let temperatures: Vec<f64> = (0..n_temps)
55        .map(|i| 0.5 * 4.0f64.powf(i as f64 / (n_temps.max(2) - 1) as f64))
56        .collect();
57    let mut rng = Xoshiro256StarStar::seed_from_u64(0x5eed);
58    let couplings: Vec<Vec<f64>> = (0..n_realizations)
59        .map(|_| {
60            (0..n_spins * dims)
61                .map(|_| match couplings_mode.as_str() {
62                    "ferro" => 1.0,
63                    "bimodal" => {
64                        if rng.gen::<bool>() {
65                            1.0
66                        } else {
67                            -1.0
68                        }
69                    }
70                    "gaussian" => {
71                        // Box-Muller; only the distribution shape matters here.
72                        let (u, v) = (1.0 - rng.gen::<f64>(), rng.gen::<f64>());
73                        (-2.0 * u.ln()).sqrt() * (std::f64::consts::TAU * v).cos()
74                    }
75                    other => panic!("unknown PEAPODS_COUPLINGS '{other}'"),
76                })
77                .collect()
78        })
79        .collect();
80    let occupation = (dilution > 0.0).then(|| {
81        (0..n_realizations)
82            .map(|_| (0..n_spins).map(|_| rng.gen::<f64>() >= dilution).collect())
83            .collect()
84    });
85    let mut simulation =
86        XySimulation::new(shape, couplings, &temperatures, n_replicas, occupation, 42)
87            .expect("invalid XY benchmark configuration");
88
89    let cluster_update = match cluster_mode.as_str() {
90        "none" => None,
91        mode => Some(ClusterConfig {
92            interval: 1,
93            mode: ClusterMode::try_from(mode).expect("PEAPODS_CLUSTER_MODE"),
94            action: ClusterAction::Update,
95            collect_stats: false,
96        }),
97    };
98    let warmup_sweeps = ((n_sweeps as f64 * warmup_ratio).round() as usize).min(n_sweeps - 1);
99    let config = XyConfig {
100        simulation: SimConfig {
101            n_sweeps,
102            warmup_sweeps,
103            sweep_mode: SweepMode::try_from(sweep_mode.as_str()).expect("PEAPODS_SWEEP_MODE"),
104            cluster_update,
105            pt_interval,
106            pt_schedule: PtSchedule::SingleRandomEdge,
107            overlap_cluster: None,
108            autocorrelation_max_lag: None,
109            autocorrelation_backend: AutocorrelationBackend::Ring,
110            sequential,
111            equilibration_diagnostic: false,
112            conditional_link_overlap: true,
113        },
114        cluster_updates: 1,
115        overrelaxation_sweeps,
116        physics: PhysicsOptions::default(),
117        measure_interval: 1,
118    };
119
120    let interrupted = AtomicBool::new(false);
121    let start = Instant::now();
122    let result = match simulation.sample(&config, &interrupted, &|| {}) {
123        Ok(result) => result,
124        Err(error) => {
125            eprintln!("unsupported configuration: {error}");
126            std::process::exit(2);
127        }
128    };
129    let elapsed = start.elapsed();
130    let energy: f64 = result.values["energies"].iter().map(|row| row[0]).sum();
131
132    let mut state_hash = DefaultHasher::new();
133    for realization in &simulation.realizations {
134        for spin in &realization.spins {
135            spin.map(f64::to_bits).hash(&mut state_hash);
136        }
137        realization.system_ids.hash(&mut state_hash);
138        for rng in &realization.rngs {
139            rng.clone().next_u64().hash(&mut state_hash);
140        }
141    }
142    for (name, rows) in &result.values {
143        name.hash(&mut state_hash);
144        for value in rows.iter().flatten() {
145            value.to_bits().hash(&mut state_hash);
146        }
147    }
148    for disorder in &result.per_disorder {
149        for value in disorder.physics.moments.iter().flatten() {
150            value.to_bits().hash(&mut state_hash);
151        }
152        disorder.visited_spins.hash(&mut state_hash);
153    }
154    println!(
155        "L={l} dims={dims} temps={n_temps} replicas={n_replicas} realizations={n_realizations} \
156         sweeps={n_sweeps} warmup={warmup_sweeps} sweep={sweep_mode} cluster={cluster_mode} \
157         overrelaxation={overrelaxation_sweeps} pt={pt_interval:?} couplings={couplings_mode} \
158         dilution={dilution} sequential={sequential} threads={}",
159        rayon::current_num_threads()
160    );
161    println!("mean energy: {:.6}", energy / n_temps as f64);
162    println!("State checksum: {:016x}", state_hash.finish());
163    println!(
164        "{:.3} ms/sweep ({:.3} s total)",
165        elapsed.as_secs_f64() * 1e3 / n_sweeps as f64,
166        elapsed.as_secs_f64()
167    );
168}