1use 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 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}