spin_sim/simulation/

mod.rs

pub(crate) mod driver;
pub mod realization;
pub mod xy;

pub use realization::Realization;

use std::sync::atomic::AtomicBool;

use crate::config::{ClusterAction, OverlapClusterBuildMode, SimConfig, SweepMode};
use crate::geometry::Lattice;
use crate::statistics::graph::{fold_histograms, histograms, ClusterObservationAccums};
use crate::statistics::physics::{CachedTotals, PhysicsCollector, PhysicsOptions};
use crate::statistics::results::{
    add_move_counts, add_rmc_counts, fold_counts, OverlapMoveCounts, RmcEdgeCounts,
};
use crate::statistics::{
    AutocorrAccum, ClusterSnapshot, ClusterStats, Diagnostics, DisorderSample,
    EquilDiagnosticAccum, OverlapAccum, Statistics, SweepResult,
};
use crate::{clusters, mcmc, spins};

/// Energies and magnetization sums of every system, kept in step with the spins.
///
/// With unit couplings a sweep reports exact integer deltas, so the cache survives
/// single-spin sweeps; every other mutation invalidates it until the next refresh.
struct ObservableCache<'a> {
    lattice: &'a Lattice,
    lookup: Option<&'a mcmc::sweep::UnitCouplingLookup>,
    sequential: bool,
    /// Exact -H per system, maintained only with a unit-coupling lookup.
    unit_totals: Vec<i64>,
    deltas: Vec<mcmc::sweep::UnitSweepDelta>,
    magnetization_sums: Vec<i64>,
    valid: bool,
}

impl<'a> ObservableCache<'a> {
    fn new(
        lattice: &'a Lattice,
        lookup: Option<&'a mcmc::sweep::UnitCouplingLookup>,
        n_systems: usize,
        config: &SimConfig,
    ) -> Self {
        let n_tracked = if lookup.is_some() { n_systems } else { 0 };
        Self {
            lattice,
            lookup,
            sequential: config.sequential,
            unit_totals: vec![0; n_tracked],
            deltas: vec![mcmc::sweep::UnitSweepDelta::default(); n_tracked],
            magnetization_sums: vec![0; n_systems],
            valid: false,
        }
    }

    /// Where the coming sweep reports its deltas. They are tracked only while the
    /// cache is valid and `consumed` before a cluster move invalidates it anyway.
    fn begin_sweep(&mut self, consumed: bool) -> Option<&mut [mcmc::sweep::UnitSweepDelta]> {
        self.valid &= self.lookup.is_some() && consumed;
        self.valid.then_some(self.deltas.as_mut_slice())
    }

    /// Applies the deltas of a tracked sweep.
    fn end_sweep(&mut self, energies: &mut [f32]) {
        if !self.valid {
            return;
        }
        for (system, delta) in self.deltas.iter().enumerate() {
            self.unit_totals[system] += delta.interaction;
            self.magnetization_sums[system] += delta.magnetization;
            energies[system] =
                spins::energy::unit_energy(self.unit_totals[system], self.lattice.n_spins);
        }
    }

    fn invalidate(&mut self) {
        self.valid = false;
    }

    /// Recomputes energies and magnetization sums unless they are current.
    fn ensure(&mut self, real: &mut Realization) {
        if self.valid {
            return;
        }
        self.valid = true;
        let Some(lookup) = self.lookup else {
            spins::energy::refresh_energies_and_magnetizations(
                self.lattice,
                &real.spins,
                &real.couplings,
                &mut real.energies,
                &mut self.magnetization_sums,
                self.sequential,
            );
            return;
        };
        spins::energy::refresh_unit_totals(
            self.lattice,
            &real.spins,
            lookup.couplings(),
            &mut self.unit_totals,
            &mut self.magnetization_sums,
            self.sequential,
        );
        for (energy, &total) in real.energies.iter_mut().zip(&self.unit_totals) {
            *energy = spins::energy::unit_energy(total, self.lattice.n_spins);
        }
    }
}

/// Run the full Monte Carlo loop (warmup + measurement) for one [`Realization`].
///
/// Each sweep consists of:
/// 1. A full single-spin pass (`sweep_mode`: Metropolis or Gibbs)
/// 2. An optional cluster update (every `cluster_update.interval` sweeps)
/// 3. Measurement (after `warmup_sweeps`)
/// 4. Optional overlap cluster move (every `overlap_cluster.interval` sweeps; `n_replicas` must reach the largest mode's group size)
/// 5. Optional parallel tempering (every `pt_interval` sweeps)
///
/// `on_sweep` is called once per sweep (useful for progress bars).
#[allow(clippy::too_many_arguments)]
pub fn run_sweep_loop(
    lattice: &Lattice,
    real: &mut Realization,
    n_replicas: usize,
    n_temps: usize,
    config: &SimConfig,
    interrupted: &AtomicBool,
    on_sweep: &(dyn Fn() + Sync),
    realization_idx: usize,
) -> Result<SweepResult, String> {
    driver::validate_layout(config, n_replicas, n_temps)?;
    real.validate(lattice, n_replicas, n_temps)?;
    run_sweep_loop_impl(
        lattice,
        real,
        n_replicas,
        n_temps,
        config,
        interrupted,
        on_sweep,
        realization_idx,
        true,
        None,
    )
}

#[allow(clippy::too_many_arguments)]
fn run_sweep_loop_impl(
    lattice: &Lattice,
    real: &mut Realization,
    n_replicas: usize,
    n_temps: usize,
    config: &SimConfig,
    interrupted: &AtomicBool,
    on_sweep: &(dyn Fn() + Sync),
    realization_idx: usize,
    materialize_disabled_stats: bool,
    physics_options: Option<&PhysicsOptions>,
) -> Result<SweepResult, String> {
    let acceptance = match config.sweep_mode {
        SweepMode::Metropolis => Some(mcmc::sweep::Acceptance::Metropolis),
        SweepMode::Gibbs => Some(mcmc::sweep::Acceptance::Gibbs),
        SweepMode::None => None,
    };
    let unit_lookup = acceptance.and_then(|acceptance| {
        mcmc::sweep::UnitCouplingLookup::new(
            &real.couplings,
            &real.temperatures,
            lattice.n_neighbors,
            acceptance,
        )
    });

    let n_spins = lattice.n_spins;
    let n_systems = n_replicas * n_temps;
    let n_sweeps = config.n_sweeps;

    let n_modes = config.overlap_cluster.as_ref().map_or(0, |h| h.modes.len());

    let n_pairs = n_replicas / 2;
    let overlap_active: Vec<bool> = real.temperatures[..n_temps]
        .iter()
        .map(|&t| {
            config
                .overlap_cluster
                .as_ref()
                .and_then(|oc| oc.max_temperature)
                .is_none_or(|max| t <= max)
        })
        .collect();
    let observe_fk = config
        .cluster_update
        .as_ref()
        .is_some_and(|cluster| cluster.action == ClusterAction::Observe);
    let observe_overlap = config
        .overlap_cluster
        .as_ref()
        .is_some_and(|cluster| cluster.action == ClusterAction::Observe);
    let collect_fk = config
        .cluster_update
        .as_ref()
        .is_some_and(|c| c.collect_stats || c.action == ClusterAction::Observe);
    let collect_overlap = config
        .overlap_cluster
        .as_ref()
        .is_some_and(|h| h.collect_stats || h.action == ClusterAction::Observe)
        && n_pairs > 0;

    let mut observation_accums = ClusterObservationAccums::default();
    if observe_fk {
        observation_accums.enable_fk(n_temps);
    }
    if let Some(overlap) = config.overlap_cluster.as_ref().filter(|_| observe_overlap) {
        for mode in &overlap.modes {
            observation_accums.enable_overlap(mode, n_temps);
        }
    }

    // The moves add to these histograms in place, one per output slot; they are
    // folded per temperature once the run ends.
    let mut fk_csd_slots = histograms(if collect_fk { n_systems } else { 0 }, n_spins);
    // Spins visited by recorded FK updates, per temperature slot.
    let fk_updates_spins = config
        .cluster_update
        .as_ref()
        .is_some_and(|cluster| cluster.action == ClusterAction::Update);
    let mut fk_visits = vec![0u64; if fk_updates_spins { n_systems } else { 0 }];
    let mut fk_recorded_updates = 0u64;
    let mut fk_observation_buf = if observe_fk {
        vec![clusters::GraphObservationSlot::default(); n_systems]
    } else {
        vec![]
    };

    let mut overlap_csd_slots: Vec<Vec<Vec<u64>>> = (0..n_modes)
        .map(|_| {
            histograms(
                if collect_overlap {
                    n_temps * n_pairs
                } else {
                    0
                },
                n_spins,
            )
        })
        .collect();

    let collect_top = collect_overlap;
    // Move counters are kept with the other overlap statistics, per output slot.
    let collect_moves = config
        .overlap_cluster
        .as_ref()
        .is_some_and(|overlap| overlap.collect_stats && overlap.action == ClusterAction::Update);
    let mut move_slots: Vec<Vec<OverlapMoveCounts>> = (0..n_modes)
        .map(|_| {
            vec![OverlapMoveCounts::default(); if collect_moves { n_temps * n_pairs } else { 0 }]
        })
        .collect();
    let mut rmc_slots = vec![RmcEdgeCounts::default(); if collect_moves { n_systems } else { 0 }];

    let mut top4_accum: Vec<Vec<[f64; 4]>> = if collect_top {
        (0..n_modes).map(|_| vec![[0.0; 4]; n_temps]).collect()
    } else {
        vec![]
    };
    let mut top4_n: Vec<usize> = if collect_top {
        vec![0; n_modes]
    } else {
        vec![]
    };
    let mut top4_buf: Vec<[u32; 4]> = if collect_top {
        vec![[0u32; 4]; n_temps * n_pairs]
    } else {
        vec![]
    };
    let mut overlap_observation_buf = if observe_overlap {
        vec![clusters::GraphObservationSlot::default(); n_temps * n_pairs]
    } else {
        vec![]
    };

    // A snapshot waits for the next overlap call that publishes a graph.
    let mut snapshot_due = false;

    let snapshot_interval = if realization_idx == 0 {
        config
            .overlap_cluster
            .as_ref()
            .and_then(|oc| oc.snapshot_interval)
    } else {
        None
    };
    let n_pair_slots = n_temps * n_pairs;
    let mut snap_buf: Vec<Vec<u32>> = if snapshot_interval.is_some() {
        (0..n_pair_slots)
            .map(|_| Vec::with_capacity(n_spins))
            .collect()
    } else {
        vec![]
    };
    let mut blue_snap_buf: Vec<Vec<u32>> = if snapshot_interval.is_some() {
        (0..n_pair_slots)
            .map(|_| Vec::with_capacity(n_spins))
            .collect()
    } else {
        vec![]
    };
    let mut spin_snap_buf: Vec<Vec<[Vec<i8>; 2]>> = if snapshot_interval.is_some() {
        (0..n_pair_slots).map(|_| Vec::new()).collect()
    } else {
        vec![]
    };
    let mut sid_snap_buf: Vec<Vec<[usize; 2]>> = if snapshot_interval.is_some() {
        (0..n_pair_slots).map(|_| Vec::new()).collect()
    } else {
        vec![]
    };
    let mut cluster_snapshots: Vec<ClusterSnapshot> = Vec::new();

    let mut mags_stat = Statistics::new(n_temps, 1);
    let mut mags2_stat = Statistics::new(n_temps, 1);
    let mut mags4_stat = Statistics::new(n_temps, 1);
    let mut energies_stat = Statistics::new(n_temps, 1);
    let mut energies2_stat = Statistics::new(n_temps, 2);
    let n_measurement_sweeps = config.n_measurement_sweeps();
    let ac_max_lag = config.autocorrelation_lag();
    let autocorrelation = || {
        ac_max_lag.map(|k| {
            AutocorrAccum::with_backend(
                k,
                n_temps,
                config.autocorrelation_backend,
                n_measurement_sweeps,
            )
        })
    };
    let mut m2_accum = autocorrelation();
    let mut energy_accum = autocorrelation();
    let mut q2_accum = autocorrelation().filter(|_| n_pairs > 0);
    let collect_ac = ac_max_lag.is_some();
    let collect_q2_ac = q2_accum.is_some();
    // Replica-averaged m² and energy of the sweep being measured.
    let mut m2_ac_buf = vec![0.0f64; if collect_ac { n_temps } else { 0 }];
    let mut energy_ac_buf = m2_ac_buf.clone();

    let equil_diag = config.equilibration_diagnostic;
    let mut equil_accum = if equil_diag {
        Some(EquilDiagnosticAccum::new(n_temps, n_sweeps))
    } else {
        None
    };
    let mut diag_e_buf = if equil_diag {
        vec![0.0f32; n_temps]
    } else {
        vec![]
    };

    let mut ov_accum = OverlapAccum::new(
        n_temps,
        n_spins,
        n_pairs,
        lattice.n_neighbors,
        equil_diag,
        collect_q2_ac,
    )
    .with_replicas(n_replicas);
    if !config.conditional_link_overlap {
        ov_accum = ov_accum.without_conditional_link_overlap();
    }

    let mut mags_buf = vec![0.0f32; n_temps];
    let mut mags2_buf = vec![0.0f32; n_temps];
    let mut mags4_buf = vec![0.0f32; n_temps];
    let mut energies_buf = vec![0.0f32; n_temps];
    let mut cache = ObservableCache::new(lattice, unit_lookup.as_ref(), n_systems, config);

    let mut physics = physics_options
        .map(|options| PhysicsCollector::new(lattice, n_temps, options, 1))
        .transpose()?
        .map(|collector| {
            if n_pairs > 0 {
                collector.with_overlap(lattice)
            } else {
                collector
            }
        });
    let occupied = if physics.is_some() {
        vec![true; n_spins]
    } else {
        vec![]
    };

    let first_sweep = real.schedule.sweeps;
    driver::run_sweeps(config, interrupted, on_sweep, first_sweep, |step| {
        let sweep_id = step.index;
        let record = step.record;
        snapshot_due |= record && snapshot_interval.is_some_and(|si| step.clock.is_multiple_of(si));
        if let Some(acceptance) = acceptance {
            // Track only if the deltas can be consumed before a cluster move
            // invalidates them anyway.
            let cluster_mutates = step.cluster
                && config
                    .cluster_update
                    .as_ref()
                    .is_some_and(|cluster| cluster.action == ClusterAction::Update);
            let overlap_mutates_first = !(record || equil_diag)
                && config
                    .overlap_cluster
                    .as_ref()
                    .is_some_and(|overlap| overlap.action == ClusterAction::Update && step.overlap);
            let deltas = cache.begin_sweep(!cluster_mutates && !overlap_mutates_first);
            mcmc::sweep::single_spin_sweep(
                lattice,
                &mut real.spins,
                &real.couplings,
                &real.temperatures,
                &real.system_ids,
                &mut real.rngs,
                config.sequential,
                acceptance,
                unit_lookup.as_ref(),
                deltas,
            );
            cache.end_sweep(&mut real.energies);
        }

        // An observation changes nothing, so one that is not recorded is skipped.
        let cluster_cfg = config
            .cluster_update
            .as_ref()
            .filter(|cluster| step.cluster && (cluster.action == ClusterAction::Update || record));
        if let Some(cluster_cfg) = cluster_cfg {
            let wolff = cluster_cfg.mode == crate::config::ClusterMode::Wolff;
            let csd_out = (collect_fk && record).then_some(fk_csd_slots.as_mut_slice());
            let observation_out = if cluster_cfg.action == ClusterAction::Observe && record {
                fk_observation_buf.fill(clusters::GraphObservationSlot::default());
                Some(fk_observation_buf.as_mut_slice())
            } else {
                None
            };

            clusters::fk_update(
                lattice,
                &mut real.spins,
                &real.couplings,
                &real.temperatures,
                &real.system_ids,
                &mut real.rngs,
                wolff,
                cluster_cfg.action,
                csd_out,
                observation_out,
                config.sequential,
                (fk_updates_spins && record).then_some(fk_visits.as_mut_slice()),
            );
            if cluster_cfg.action == ClusterAction::Update {
                cache.invalidate();
                fk_recorded_updates += u64::from(record);
            }

            if let Some(accum) = observation_accums.fk.as_mut().filter(|_| record) {
                for (slot, observation) in fk_observation_buf.iter().copied().enumerate() {
                    accum.push(slot % n_temps, observation);
                }
            }
        }

        // A single temperature has no edge to swap, so skip the energy refresh too.
        let pt_this_sweep = step.temper && n_temps > 1;

        // Tempering refreshes for itself below, after any overlap move.
        if record || equil_diag {
            cache.ensure(real);
        }

        if equil_diag {
            diag_e_buf.fill(0.0);
            #[allow(clippy::needless_range_loop)]
            for r in 0..n_replicas {
                let offset = r * n_temps;
                for t in 0..n_temps {
                    let system_id = real.system_ids[offset + t];
                    diag_e_buf[t] += real.energies[system_id];
                }
            }
            let inv = 1.0 / n_replicas as f32;
            for v in diag_e_buf.iter_mut() {
                *v *= inv;
            }
        }

        if (equil_diag || record) && n_pairs > 0 {
            ov_accum.collect(lattice, &real.spins, &real.system_ids, record);
        }

        if equil_diag {
            equil_accum.as_mut().unwrap().push(
                &diag_e_buf,
                &ov_accum.diag_ql_buf,
                &ov_accum.diag_q2_buf,
            );
        }

        if record {
            for t in 0..n_temps {
                mags_buf[t] = 0.0;
                mags2_buf[t] = 0.0;
                mags4_buf[t] = 0.0;
                energies_buf[t] = 0.0;
            }

            if collect_ac {
                m2_ac_buf.fill(0.0);
                energy_ac_buf.fill(0.0);
            }

            for r in 0..n_replicas {
                let offset = r * n_temps;
                for t in 0..n_temps {
                    let system_id = real.system_ids[offset + t];
                    let mag = cache.magnetization_sums[system_id] as f32 / n_spins as f32;
                    let m2 = mag * mag;
                    mags_buf[t] = mag;
                    mags2_buf[t] = m2;
                    mags4_buf[t] = m2 * m2;
                    energies_buf[t] = real.energies[system_id];
                }

                if collect_ac {
                    for t in 0..n_temps {
                        m2_ac_buf[t] += mags2_buf[t] as f64;
                        energy_ac_buf[t] += energies_buf[t] as f64;
                    }
                }

                mags_stat.update(&mags_buf);
                mags2_stat.update(&mags2_buf);
                mags4_stat.update(&mags4_buf);
                energies_stat.update(&energies_buf);
                energies2_stat.update(&energies_buf);
            }

            if let (Some(m2), Some(energy)) = (m2_accum.as_mut(), energy_accum.as_mut()) {
                let inv = 1.0 / n_replicas as f64;
                for v in m2_ac_buf.iter_mut().chain(energy_ac_buf.iter_mut()) {
                    *v *= inv;
                }
                m2.push(&m2_ac_buf);
                energy.push(&energy_ac_buf);
            }

            if let Some(ref mut acc) = q2_accum {
                acc.push(&ov_accum.q2_ac_buf);
            }
        }

        if let Some(collector) = physics.as_mut().filter(|_| record) {
            // Reuse this sweep's energy/magnetization pass; `energies` hold -H/N.
            for (slot, &system) in real.system_ids.iter().enumerate() {
                collector.measure_cached(
                    lattice,
                    &real.spins[system * n_spins..(system + 1) * n_spins],
                    &real.couplings,
                    &occupied,
                    slot % n_temps,
                    CachedTotals {
                        energy: -(real.energies[system] as f64),
                        magnetization: [cache.magnetization_sums[system] as f64, 0.0],
                    },
                );
            }
            // Same pairs (2p, 2p + 1) as the legacy overlap statistics.
            let weight = n_replicas as f64 / n_pairs.max(1) as f64;
            for pair in 0..n_pairs {
                for t in 0..n_temps {
                    let a = real.system_ids[2 * pair * n_temps + t];
                    let b = real.system_ids[(2 * pair + 1) * n_temps + t];
                    collector.measure_overlap(
                        lattice,
                        &real.spins[a * n_spins..(a + 1) * n_spins],
                        &real.spins[b * n_spins..(b + 1) * n_spins],
                        &occupied,
                        t,
                        weight,
                    );
                }
            }
            if n_pairs > 0 {
                for t in 0..n_temps {
                    let replicas: Vec<&[i8]> = (0..n_replicas)
                        .map(|r| {
                            let system = real.system_ids[r * n_temps + t];
                            &real.spins[system * n_spins..(system + 1) * n_spins]
                        })
                        .collect();
                    collector.measure_overlap_all_pairs(&replicas, t, n_replicas as f64);
                }
            }
            collector.end_sweep(n_replicas);
        }

        if let Some(oc_cfg) = config.overlap_cluster.as_ref().filter(|_| step.overlap) {
            let mode_idx = real.schedule.overlap_calls % n_modes;
            let mode = &oc_cfg.modes[mode_idx];
            // A skipped observation still counts, which keeps the mode rotation fixed.
            real.schedule.overlap_calls += 1;
            if oc_cfg.action == ClusterAction::Update || record {
                if oc_cfg.action == ClusterAction::Update {
                    cache.invalidate();
                }

                let ov_csd_out = (collect_overlap && record)
                    .then_some(overlap_csd_slots[mode_idx].as_mut_slice());
                let observation_out = if oc_cfg.action == ClusterAction::Observe && record {
                    overlap_observation_buf.fill(clusters::GraphObservationSlot::default());
                    Some(overlap_observation_buf.as_mut_slice())
                } else {
                    None
                };

                let top4_out = if collect_top && record {
                    for slot in top4_buf.iter_mut() {
                        *slot = [0u32; 4];
                    }
                    Some(top4_buf.as_mut_slice())
                } else {
                    None
                };

                let is_rmc = *mode == OverlapClusterBuildMode::Rmc;
                // Replica Monte Carlo builds no overlap graph to publish.
                let take_snapshot = snapshot_due && !is_rmc;
                snapshot_due &= !take_snapshot;

                let is_cmr = matches!(mode, crate::config::OverlapClusterBuildMode::Cmr);

                let snap = if take_snapshot {
                    for buf in spin_snap_buf.iter_mut() {
                        buf.clear();
                    }
                    for buf in sid_snap_buf.iter_mut() {
                        buf.clear();
                    }
                    Some(snap_buf.as_mut_slice())
                } else {
                    None
                };
                let blue_snap = if take_snapshot && is_cmr {
                    Some(blue_snap_buf.as_mut_slice())
                } else {
                    None
                };
                let spin_snap = if take_snapshot {
                    Some(spin_snap_buf.as_mut_slice())
                } else {
                    None
                };
                let sid_snap = if take_snapshot {
                    Some(sid_snap_buf.as_mut_slice())
                } else {
                    None
                };

                if is_rmc {
                    clusters::rmc_update(
                        lattice,
                        &mut real.spins,
                        &real.couplings,
                        &real.temperatures,
                        &real.system_ids,
                        n_replicas,
                        n_temps,
                        &overlap_active,
                        &mut real.rngs,
                        oc_cfg.cluster_mode,
                        config.sequential,
                        (collect_moves && record).then_some(rmc_slots.as_mut_slice()),
                    );
                } else {
                    clusters::overlap_update(
                        lattice,
                        &mut real.spins,
                        &real.couplings,
                        &real.temperatures,
                        &real.system_ids,
                        n_replicas,
                        n_temps,
                        &overlap_active,
                        &mut real.pair_rngs,
                        mode,
                        oc_cfg.cluster_mode,
                        oc_cfg.action,
                        ov_csd_out,
                        top4_out,
                        observation_out,
                        config.sequential,
                        snap,
                        blue_snap,
                        spin_snap,
                        sid_snap,
                        (collect_moves && record && n_pairs > 0)
                            .then_some(move_slots[mode_idx].as_mut_slice()),
                    );
                }

                if take_snapshot {
                    let ids: Vec<Vec<u32>> = (0..n_temps)
                        .map(|t| snap_buf[t * n_pairs].clone())
                        .collect();
                    let blue = if is_cmr {
                        Some(
                            (0..n_temps)
                                .map(|t| blue_snap_buf[t * n_pairs].clone())
                                .collect(),
                        )
                    } else {
                        None
                    };
                    let spins: Vec<[Vec<i8>; 2]> = (0..n_temps)
                        .map(|t| {
                            spin_snap_buf[t * n_pairs]
                                .first()
                                .cloned()
                                .unwrap_or_else(|| [vec![], vec![]])
                        })
                        .collect();
                    let sids: Vec<[usize; 2]> = (0..n_temps)
                        .map(|t| sid_snap_buf[t * n_pairs].first().copied().unwrap_or([0, 0]))
                        .collect();
                    cluster_snapshots.push(ClusterSnapshot {
                        sweep_id,
                        mode_idx,
                        cluster_ids: ids,
                        blue_ids: blue,
                        spins,
                        system_ids: sids,
                    });
                }

                if collect_top && record {
                    for t in 0..n_temps {
                        for p in 0..n_pairs {
                            let raw = top4_buf[t * n_pairs + p];
                            for (k, &v) in raw.iter().enumerate() {
                                top4_accum[mode_idx][t][k] += v as f64 / n_spins as f64;
                            }
                        }
                    }
                    top4_n[mode_idx] += 1;
                }

                if record && oc_cfg.action == ClusterAction::Observe {
                    let accum = observation_accums.overlap_mut(mode);
                    for (slot, observation) in overlap_observation_buf.iter().copied().enumerate() {
                        accum.push(slot / n_pairs, observation);
                    }
                }
            }
        }

        if pt_this_sweep {
            cache.ensure(real);
            real.temper(n_spins, n_replicas, n_temps, config.pt_schedule);
        }
        real.schedule.sweeps += 1;
    })?;

    let top_cluster_sizes: Vec<Vec<[f64; 4]>> = if collect_top {
        let modes = &config.overlap_cluster.as_ref().unwrap().modes;
        top4_accum
            .iter()
            .zip(top4_n.iter())
            .zip(modes)
            .map(|((mode_accum, &count), mode)| {
                // Each temperature publishes one top-4 record per replica group.
                let n_groups = match mode {
                    OverlapClusterBuildMode::Rmc => 0,
                    _ => n_replicas / mode.group_size(),
                };
                if count == 0 || n_groups == 0 {
                    return vec![[0.0; 4]; n_temps];
                }
                let denom = (count * n_groups) as f64;
                mode_accum
                    .iter()
                    .map(|arr| {
                        [
                            arr[0] / denom,
                            arr[1] / denom,
                            arr[2] / denom,
                            arr[3] / denom,
                        ]
                    })
                    .collect()
            })
            .collect()
    } else {
        vec![]
    };

    let (mags2_tau, mags2_tau_converged) = driver::tau_estimates(m2_accum.as_ref());
    let (overlap2_tau, overlap2_tau_converged) = driver::tau_estimates(q2_accum.as_ref());
    let (energy_tau, energy_tau_converged) = driver::tau_estimates(energy_accum.as_ref());

    let equil_checkpoints = equil_accum.map(|acc| acc.finish()).unwrap_or_default();
    // Wrapping is tracked along every axis of canonical square and cubic lattices.
    let winding_dims = if lattice.square_shape().is_some() {
        2
    } else if lattice.cubic_shape().is_some() {
        3
    } else {
        0
    };

    // Disabled histograms keep their shape only where the caller asked for it.
    let keep_shape = |collected: bool| collected || materialize_disabled_stats;
    let fk_csd = if keep_shape(collect_fk) {
        fold_histograms(&fk_csd_slots, n_temps, n_spins, |slot| slot % n_temps)
    } else {
        vec![]
    };
    let overlap_csd: Vec<Vec<Vec<u64>>> = if keep_shape(collect_overlap) {
        overlap_csd_slots
            .iter()
            .map(|slots| fold_histograms(slots, n_temps, n_spins, |slot| slot / n_pairs))
            .collect()
    } else {
        vec![]
    };
    let cluster_observations = observation_accums.finish(
        n_spins,
        lattice.n_neighbors,
        winding_dims,
        &fk_csd,
        config.overlap_cluster.as_ref().map_or(&[], |o| &o.modes),
        &overlap_csd,
    );

    let overlap_stats = ov_accum.finish();
    let mags2 = mags2_stat.average();
    let mags4 = mags4_stat.average();
    let energies = energies_stat.average();
    let energy_variance = Statistics::variance(&energies_stat, &energies2_stat);
    let sample = DisorderSample {
        mags2: mags2.clone(),
        mags4: mags4.clone(),
        energies: energies.clone(),
        energy_variance: energy_variance.clone(),
        overlap2: overlap_stats.overlap2.clone(),
        overlap4: overlap_stats.overlap4.clone(),
        link_overlap: overlap_stats.link_overlap.clone(),
        link_overlap2: overlap_stats.link_overlap2.clone(),
        link_overlap4: overlap_stats.link_overlap4.clone(),
        mags2_tau: mags2_tau.clone(),
        overlap2_tau: overlap2_tau.clone(),
        energy_tau: energy_tau.clone(),
        mags2_tau_converged: mags2_tau_converged.clone(),
        overlap2_tau_converged: overlap2_tau_converged.clone(),
        energy_tau_converged: energy_tau_converged.clone(),
    };

    Ok(SweepResult {
        mags: mags_stat.average(),
        mags2,
        mags4,
        energies,
        energies2: energies2_stat.average(),
        energy_variance,
        overlap_stats,
        cluster_stats: ClusterStats {
            fk_csd,
            overlap_csd,
            top_cluster_sizes,
            fk_updates: if fk_updates_spins {
                vec![fk_recorded_updates * n_replicas as u64; n_temps]
            } else {
                vec![]
            },
            overlap_moves: if collect_moves {
                config
                    .overlap_cluster
                    .as_ref()
                    .map_or(&[][..], |overlap| &overlap.modes)
                    .iter()
                    .zip(&move_slots)
                    .map(|(mode, slots)| match mode {
                        OverlapClusterBuildMode::Rmc => vec![],
                        _ => fold_counts(
                            slots,
                            n_temps,
                            |slot| slot / n_pairs.max(1),
                            add_move_counts,
                        ),
                    })
                    .collect()
            } else {
                vec![]
            },
            rmc_edges: if collect_moves
                && config
                    .overlap_cluster
                    .as_ref()
                    .is_some_and(|overlap| overlap.modes.contains(&OverlapClusterBuildMode::Rmc))
            {
                let mut edges =
                    fold_counts(&rmc_slots, n_temps, |slot| slot % n_temps, add_rmc_counts);
                edges.pop();
                edges
            } else {
                vec![]
            },
            fk_visited_spins: if fk_updates_spins {
                (0..n_temps)
                    .map(|t| fk_visits.iter().skip(t).step_by(n_temps).sum())
                    .collect()
            } else {
                vec![]
            },
        },
        per_disorder_cluster_observations: vec![cluster_observations],
        per_disorder_physics: physics.map(|p| vec![p.finish()]).unwrap_or_default(),
        per_disorder: vec![sample],
        diagnostics: Diagnostics {
            mags2_tau,
            overlap2_tau,
            energy_tau,
            mags2_tau_converged,
            overlap2_tau_converged,
            energy_tau_converged,
            autocorrelation_lag: ac_max_lag,
            equil_checkpoints,
        },
        cluster_snapshots,
    })
}

/// Run the sweep loop in parallel over multiple disorder realizations.
///
/// The batch is validated before any realization changes, each realization then
/// runs the same loop as [`run_sweep_loop`], and the results are averaged over
/// disorder with the per-sample overlap histograms kept. Replicas run on the
/// calling thread once the realizations alone fill the thread pool. A single
/// realization is run directly, skipping rayon thread-pool overhead.
pub fn run_sweep_parallel(
    lattice: &Lattice,
    realizations: &mut [Realization],
    n_replicas: usize,
    n_temps: usize,
    config: &SimConfig,
    interrupted: &AtomicBool,
    on_sweep: &(dyn Fn() + Sync),
) -> Result<SweepResult, String> {
    run_sweep_parallel_with_physics(
        lattice,
        realizations,
        n_replicas,
        n_temps,
        config,
        interrupted,
        on_sweep,
        None,
    )
}

/// Opt-in physical moments; the legacy entrypoint retains its measurement cost.
#[allow(clippy::too_many_arguments)]
pub fn run_sweep_parallel_with_physics(
    lattice: &Lattice,
    realizations: &mut [Realization],
    n_replicas: usize,
    n_temps: usize,
    config: &SimConfig,
    interrupted: &AtomicBool,
    on_sweep: &(dyn Fn() + Sync),
    physics: Option<&PhysicsOptions>,
) -> Result<SweepResult, String> {
    driver::validate_batch(lattice, realizations, n_replicas, n_temps, config)?;
    if let Some(options) = physics {
        options.validate(lattice, 1)?;
    }
    let single = realizations.len() == 1;
    let mut effective = config.clone();
    effective.sequential = driver::effective_sequential(config, realizations.len());
    let config = &effective;
    let mut results = driver::map_realizations(realizations, |idx, real| {
        run_sweep_loop_impl(
            lattice,
            real,
            n_replicas,
            n_temps,
            config,
            interrupted,
            on_sweep,
            idx,
            single,
            physics,
        )
    })?;
    if single {
        return Ok(results.remove(0));
    }
    let snapshots = std::mem::take(&mut results[0].cluster_snapshots);
    let physics = results
        .iter_mut()
        .flat_map(|r| std::mem::take(&mut r.per_disorder_physics))
        .collect();
    let mut agg = SweepResult::aggregate_without_overlap_samples(&results);
    agg.per_disorder_physics = physics;
    if !agg.overlap_stats.overlap.is_empty() {
        agg.overlap_stats.per_sample_histogram = results
            .iter_mut()
            .map(|r| std::mem::take(&mut r.overlap_stats.histogram))
            .collect();
    }
    if !agg.overlap_stats.ql_at_q_sum.is_empty() {
        agg.overlap_stats.per_sample_ql_at_q_sum = results
            .iter_mut()
            .map(|r| std::mem::take(&mut r.overlap_stats.ql_at_q_sum))
            .collect();
        agg.overlap_stats.per_sample_ql2_at_q_sum = results
            .iter_mut()
            .map(|r| std::mem::take(&mut r.overlap_stats.ql2_at_q_sum))
            .collect();
    }
    if agg.cluster_stats.fk_csd.is_empty() {
        agg.cluster_stats.fk_csd = (0..n_temps)
            .map(|_| vec![0u64; lattice.n_spins + 1])
            .collect();
    }
    let n_modes = config.overlap_cluster.as_ref().map_or(0, |h| h.modes.len());
    if n_modes > 0 && agg.cluster_stats.overlap_csd.is_empty() {
        agg.cluster_stats.overlap_csd = (0..n_modes)
            .map(|_| {
                (0..n_temps)
                    .map(|_| vec![0u64; lattice.n_spins + 1])
                    .collect()
            })
            .collect();
    }
    agg.cluster_snapshots = snapshots;
    Ok(agg)
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::config::{
        AutocorrelationBackend, ClusterConfig, ClusterMode, OverlapClusterBuildMode,
        OverlapClusterConfig, PtSchedule,
    };

    fn run_with_cluster_stats(collect_stats: bool) -> SweepResult {
        let lattice = Lattice::new(vec![4, 4]);
        let couplings = vec![1.0; lattice.n_spins * lattice.n_neighbors];
        let mut realization = Realization::new(&lattice, couplings, &[2.0], 1, 42);
        let config = SimConfig {
            n_sweeps: 1,
            warmup_sweeps: 0,
            sweep_mode: SweepMode::Metropolis,
            cluster_update: Some(ClusterConfig {
                interval: 1,
                mode: ClusterMode::Sw,
                action: ClusterAction::Update,
                collect_stats,
            }),
            pt_interval: None,
            pt_schedule: PtSchedule::SingleRandomEdge,
            overlap_cluster: None,
            autocorrelation_max_lag: None,
            autocorrelation_backend: AutocorrelationBackend::Ring,
            sequential: true,
            equilibration_diagnostic: false,
            conditional_link_overlap: true,
        };

        run_sweep_loop(
            &lattice,
            &mut realization,
            1,
            1,
            &config,
            &AtomicBool::new(false),
            &|| {},
            0,
        )
        .unwrap()
    }

    fn plain_config(n_sweeps: usize) -> SimConfig {
        SimConfig {
            n_sweeps,
            warmup_sweeps: 0,
            sweep_mode: SweepMode::Metropolis,
            cluster_update: None,
            pt_interval: None,
            pt_schedule: PtSchedule::SingleRandomEdge,
            overlap_cluster: None,
            autocorrelation_max_lag: None,
            autocorrelation_backend: AutocorrelationBackend::Ring,
            sequential: true,
            equilibration_diagnostic: false,
            conditional_link_overlap: true,
        }
    }

    fn run(
        lattice: &Lattice,
        real: &mut Realization,
        n_replicas: usize,
        config: &SimConfig,
    ) -> SweepResult {
        let n_temps = real.temperatures.len() / n_replicas;
        run_sweep_loop(
            lattice,
            real,
            n_replicas,
            n_temps,
            config,
            &AtomicBool::new(false),
            &|| {},
            0,
        )
        .unwrap()
    }

    /// An event due every two sweeps must stay every two sweeps when the run is
    /// split into one-sweep calls.
    #[test]
    fn interval_phase_spans_calls() {
        let lattice = Lattice::new(vec![4, 4]);
        let couplings = vec![1.0; lattice.n_spins * lattice.n_neighbors];
        let mut real = Realization::new(&lattice, couplings, &[1.5, 2.5], 1, 9);
        let mut config = plain_config(1);
        config.pt_interval = Some(2);
        for _ in 0..6 {
            run(&lattice, &mut real, 1, &config);
        }
        assert_eq!(real.pt_edge_attempts(), &[3]);

        real.reset(&lattice, 1, 2, 9);
        run(&lattice, &mut real, 1, &config);
        assert_eq!(real.pt_edge_attempts(), &[1], "reset restarts the phase");
    }

    #[test]
    fn conditional_link_overlap_can_be_turned_off() {
        let lattice = Lattice::new(vec![4, 4]);
        let couplings = vec![1.0; lattice.n_spins * lattice.n_neighbors];
        let mut config = plain_config(4);
        for enabled in [true, false] {
            config.conditional_link_overlap = enabled;
            let mut real = Realization::new(&lattice, couplings.clone(), &[2.0], 2, 3);
            let stats = run(&lattice, &mut real, 2, &config).overlap_stats;
            assert_eq!(stats.histogram[0].iter().sum::<u64>(), 4);
            assert_eq!(stats.ql_at_q_sum.is_empty(), !enabled);
            assert_eq!(stats.ql2_at_q_sum.is_empty(), !enabled);
        }
    }

    /// Move counters must be consistent: no update moves more sites than are active,
    /// and replica Monte Carlo accepts at most what it proposes.
    #[test]
    fn overlap_move_counters_are_consistent() {
        let lattice = Lattice::new(vec![6, 6]);
        let couplings: Vec<f32> = (0..lattice.n_spins * lattice.n_neighbors)
            .map(|i| if i % 3 == 0 { -1.0 } else { 1.0 })
            .collect();
        let temps = [0.8, 1.2, 2.0];
        let mut config = plain_config(60);
        config.warmup_sweeps = 10;
        config.pt_interval = Some(1);
        config.overlap_cluster = Some(OverlapClusterConfig {
            interval: 1,
            modes: vec![
                OverlapClusterBuildMode::Houdayer(2),
                OverlapClusterBuildMode::Cmr,
                OverlapClusterBuildMode::Rmc,
            ],
            cluster_mode: ClusterMode::Sw,
            action: ClusterAction::Update,
            collect_stats: true,
            snapshot_interval: None,
            max_temperature: None,
        });
        let mut real = Realization::new(&lattice, couplings, &temps, 2, 11);
        let stats = run(&lattice, &mut real, 2, &config).cluster_stats;
        assert_eq!(stats.overlap_moves.len(), 3);
        assert!(
            stats.overlap_moves[2].is_empty(),
            "rmc publishes no group moves"
        );
        // Recorded sweeps 10..60 rotate through the modes by sweep: Houdayer runs on
        // the 16 multiples of 3, CMR on 17 sweeps; one group per temperature.
        for (counts, calls) in stats.overlap_moves[..2].iter().zip([16, 17]) {
            for counts in counts {
                assert_eq!((counts.groups, counts.counted_groups), (calls, calls));
                assert!(counts.moved_sites <= counts.active_sites);
            }
        }
        let n = lattice.n_spins as u64;
        assert!(stats.overlap_moves[1]
            .iter()
            .all(|c| c.active_sites == 17 * n));
        assert_eq!(stats.rmc_edges.len(), temps.len() - 1);
        for edge in &stats.rmc_edges {
            assert!(edge.proposed > 0 && edge.accepted <= edge.proposed);
            assert!(edge.flipped_sites <= 2 * 16 * n);
        }
    }

    /// A call that only thermalizes measures nothing, so its moments are undefined.
    #[test]
    fn warmup_only_calls_report_nan_moments() {
        let lattice = Lattice::new(vec![4, 4]);
        let couplings = vec![1.0; lattice.n_spins * lattice.n_neighbors];
        let mut real = Realization::new(&lattice, couplings, &[2.0], 2, 3);
        let mut config = plain_config(5);
        config.warmup_sweeps = 5;
        config.autocorrelation_max_lag = Some(2);
        let result = run(&lattice, &mut real, 2, &config);
        assert!(result.mags2[0].is_nan() && result.energies[0].is_nan());
        assert!(result.overlap_stats.overlap2[0].is_nan());
        assert!(result.diagnostics.mags2_tau[0].is_nan());
        assert_eq!(result.diagnostics.mags2_tau_converged, vec![false]);
        assert_eq!(result.per_disorder.len(), 1);
    }

    /// Unit-coupling sweeps track energies from exact per-flip deltas; after a long run
    /// with tempering, overlap and cluster moves they must equal a fresh recomputation.
    #[test]
    fn tracked_energies_match_recomputation() {
        let lattice = Lattice::new(vec![6, 5]);
        let couplings: Vec<f32> = (0..lattice.n_spins * lattice.n_neighbors)
            .map(|i| [1.0, -1.0, 1.0, 0.0][i % 4])
            .collect();
        let temps = [0.6, 1.1, 1.9, 3.0];
        for (sweep_mode, with_clusters) in [
            (SweepMode::Metropolis, false),
            (SweepMode::Gibbs, false),
            (SweepMode::Metropolis, true),
        ] {
            let mut realization = Realization::new(&lattice, couplings.clone(), &temps, 2, 5);
            let config = SimConfig {
                n_sweeps: 300,
                warmup_sweeps: 100,
                sweep_mode,
                cluster_update: with_clusters.then_some(ClusterConfig {
                    interval: 7,
                    mode: ClusterMode::Sw,
                    action: ClusterAction::Update,
                    collect_stats: false,
                }),
                pt_interval: Some(1),
                pt_schedule: PtSchedule::FullLadder,
                overlap_cluster: with_clusters.then_some(OverlapClusterConfig {
                    interval: 5,
                    modes: vec![OverlapClusterBuildMode::Cmr, OverlapClusterBuildMode::Rmc],
                    cluster_mode: ClusterMode::Sw,
                    action: ClusterAction::Update,
                    collect_stats: false,
                    snapshot_interval: None,
                    max_temperature: None,
                }),
                autocorrelation_max_lag: None,
                autocorrelation_backend: AutocorrelationBackend::Ring,
                sequential: false,
                equilibration_diagnostic: false,
                conditional_link_overlap: true,
            };
            run_sweep_loop(
                &lattice,
                &mut realization,
                2,
                temps.len(),
                &config,
                &AtomicBool::new(false),
                &|| {},
                0,
            )
            .unwrap();
            let mut expected = vec![0.0; realization.energies.len()];
            spins::energy::compute_energies_into(
                &lattice,
                &realization.spins,
                &realization.couplings,
                &mut expected,
            );
            // Cluster intervals 7 and 5 skip the final sweep (index 299), whose tempering
            // step refreshes observables after the last mutating move.
            assert_eq!(
                realization.energies, expected,
                "{sweep_mode:?} clusters={with_clusters}"
            );
        }
    }

    /// Overlap moves above `max_temperature` must leave those systems untouched, while
    /// the colder temperatures keep moving.
    #[test]
    fn overlap_moves_respect_max_temperature() {
        let lattice = Lattice::new(vec![6, 6]);
        let couplings: Vec<f32> = (0..lattice.n_spins * lattice.n_neighbors)
            .map(|i| if i % 3 == 0 { -1.0 } else { 1.0 })
            .collect();
        let temps = [0.5, 1.0, 3.0];
        for mode in [
            OverlapClusterBuildMode::Houdayer(2),
            OverlapClusterBuildMode::Rmc,
        ] {
            let mut realization = Realization::new(&lattice, couplings.clone(), &temps, 2, 31);
            let initial = realization.spins.clone();
            let config = SimConfig {
                n_sweeps: 20,
                warmup_sweeps: 0,
                sweep_mode: SweepMode::None,
                cluster_update: None,
                pt_interval: None,
                pt_schedule: PtSchedule::SingleRandomEdge,
                overlap_cluster: Some(OverlapClusterConfig {
                    interval: 1,
                    modes: vec![mode],
                    cluster_mode: ClusterMode::Sw,
                    action: ClusterAction::Update,
                    collect_stats: false,
                    snapshot_interval: None,
                    max_temperature: Some(1.0),
                }),
                autocorrelation_max_lag: None,
                autocorrelation_backend: AutocorrelationBackend::Ring,
                sequential: true,
                equilibration_diagnostic: false,
                conditional_link_overlap: true,
            };
            run_sweep_loop(
                &lattice,
                &mut realization,
                2,
                temps.len(),
                &config,
                &AtomicBool::new(false),
                &|| {},
                0,
            )
            .unwrap();
            let n = lattice.n_spins;
            let changed = |t: usize| {
                (0..2).any(|r| {
                    let system = realization.system_ids[r * temps.len() + t];
                    realization.spins[system * n..(system + 1) * n]
                        != initial[system * n..(system + 1) * n]
                })
            };
            assert!(!changed(2), "{mode:?} moved the T = 3.0 systems");
            assert!(changed(0), "{mode:?} never moved the T = 0.5 systems");
        }
    }

    #[test]
    fn disabled_cluster_histograms_preserve_public_shape() {
        let without_stats = run_with_cluster_stats(false);
        assert_eq!(without_stats.cluster_stats.fk_csd.len(), 1);
        assert_eq!(without_stats.cluster_stats.fk_csd[0].iter().sum::<u64>(), 0);

        let with_stats = run_with_cluster_stats(true);
        assert_eq!(with_stats.cluster_stats.fk_csd.len(), 1);
        assert!(with_stats.cluster_stats.fk_csd[0].iter().sum::<u64>() > 0);
    }

    #[test]
    fn disabled_overlap_histograms_preserve_public_shape() {
        let lattice = Lattice::new(vec![4, 4]);
        let couplings = vec![1.0; lattice.n_spins * lattice.n_neighbors];

        for collect_stats in [false, true] {
            let mut realization = Realization::new(&lattice, couplings.clone(), &[2.0], 2, 42);
            let config = SimConfig {
                n_sweeps: 1,
                warmup_sweeps: 0,
                sweep_mode: SweepMode::Metropolis,
                cluster_update: None,
                pt_interval: None,
                pt_schedule: PtSchedule::SingleRandomEdge,
                overlap_cluster: Some(OverlapClusterConfig {
                    interval: 1,
                    modes: vec![OverlapClusterBuildMode::Houdayer(2)],
                    cluster_mode: ClusterMode::Sw,
                    action: ClusterAction::Update,
                    collect_stats,
                    snapshot_interval: None,
                    max_temperature: None,
                }),
                autocorrelation_max_lag: None,
                autocorrelation_backend: AutocorrelationBackend::Ring,
                sequential: true,
                equilibration_diagnostic: false,
                conditional_link_overlap: true,
            };

            let result = run_sweep_loop(
                &lattice,
                &mut realization,
                2,
                1,
                &config,
                &AtomicBool::new(false),
                &|| {},
                0,
            )
            .unwrap();

            if collect_stats {
                assert_eq!(result.cluster_stats.overlap_csd.len(), 1);
                assert!(result.cluster_stats.overlap_csd[0][0].iter().sum::<u64>() > 0);
            } else {
                assert_eq!(result.cluster_stats.overlap_csd.len(), 1);
                assert_eq!(
                    result.cluster_stats.overlap_csd[0][0].iter().sum::<u64>(),
                    0
                );
            }
        }
    }

    /// Top-4 fractions average over the replica groups that publish them: one group of
    /// four replicas yields the fractions of its own graph, not half of them.
    #[test]
    fn top_cluster_sizes_average_over_replica_groups() {
        let lattice = Lattice::new(vec![6, 6]);
        let n = lattice.n_spins;
        let couplings = vec![1.0; n * lattice.n_neighbors];
        let mut realization = Realization::new(&lattice, couplings, &[2.0], 4, 3);
        // The balanced-site graph does not depend on how the group is shuffled.
        let balanced: Vec<bool> = (0..n)
            .map(|i| {
                (0..4)
                    .map(|r| realization.spins[r * n + i] as i32)
                    .sum::<i32>()
                    == 0
            })
            .collect();
        let mut component = vec![usize::MAX; n];
        let mut sizes = Vec::new();
        for start in 0..n {
            if component[start] != usize::MAX {
                continue;
            }
            component[start] = sizes.len();
            let (mut stack, mut size) = (vec![start], 0);
            while let Some(i) = stack.pop() {
                size += 1;
                for d in 0..lattice.n_neighbors {
                    for j in [lattice.neighbor_fwd(i, d), lattice.neighbor_bwd(i, d)] {
                        if balanced[i] && balanced[j] && component[j] == usize::MAX {
                            component[j] = sizes.len();
                            stack.push(j);
                        }
                    }
                }
            }
            sizes.push(size);
        }
        sizes.sort_unstable_by(|a, b| b.cmp(a));
        let expected: Vec<f64> = sizes[..4].iter().map(|&s| s as f64 / n as f64).collect();

        let config = SimConfig {
            n_sweeps: 1,
            warmup_sweeps: 0,
            sweep_mode: SweepMode::None,
            cluster_update: None,
            pt_interval: None,
            pt_schedule: PtSchedule::SingleRandomEdge,
            overlap_cluster: Some(OverlapClusterConfig {
                interval: 1,
                modes: vec![OverlapClusterBuildMode::Houdayer(4)],
                cluster_mode: ClusterMode::Sw,
                action: ClusterAction::Update,
                collect_stats: true,
                snapshot_interval: None,
                max_temperature: None,
            }),
            autocorrelation_max_lag: None,
            autocorrelation_backend: AutocorrelationBackend::Ring,
            sequential: true,
            equilibration_diagnostic: false,
            conditional_link_overlap: true,
        };
        let result = run_sweep_loop(
            &lattice,
            &mut realization,
            4,
            1,
            &config,
            &AtomicBool::new(false),
            &|| {},
            0,
        )
        .unwrap();
        assert_eq!(
            result.cluster_stats.top_cluster_sizes[0][0].to_vec(),
            expected
        );
    }

    #[test]
    fn fk_observe_records_without_changing_spins() {
        let lattice = Lattice::new(vec![4, 4]);
        let couplings = vec![1.0; lattice.n_spins * lattice.n_neighbors];
        let mut observed = Realization::new(&lattice, couplings.clone(), &[2.0], 1, 19);
        let mut baseline = Realization::new(&lattice, couplings, &[2.0], 1, 19);
        let mut config = SimConfig {
            n_sweeps: 1,
            warmup_sweeps: 0,
            sweep_mode: SweepMode::Metropolis,
            cluster_update: None,
            pt_interval: None,
            pt_schedule: PtSchedule::SingleRandomEdge,
            overlap_cluster: None,
            autocorrelation_max_lag: None,
            autocorrelation_backend: AutocorrelationBackend::Ring,
            sequential: true,
            equilibration_diagnostic: false,
            conditional_link_overlap: true,
        };
        run_sweep_loop(
            &lattice,
            &mut baseline,
            1,
            1,
            &config,
            &AtomicBool::new(false),
            &|| {},
            0,
        )
        .unwrap();

        config.cluster_update = Some(ClusterConfig {
            interval: 1,
            mode: ClusterMode::Sw,
            action: ClusterAction::Observe,
            collect_stats: true,
        });
        let result = run_sweep_loop(
            &lattice,
            &mut observed,
            1,
            1,
            &config,
            &AtomicBool::new(false),
            &|| {},
            0,
        )
        .unwrap();

        assert_eq!(observed.spins, baseline.spins);
        assert_eq!(observed.system_ids, baseline.system_ids);
        let summary = result.per_disorder_cluster_observations[0]
            .fk
            .as_ref()
            .unwrap();
        assert_eq!(summary.observation_count, vec![1]);
        assert_eq!(
            summary.cluster_size_counts[0]
                .iter()
                .enumerate()
                .map(|(size, &count)| size as u64 * count)
                .sum::<u64>(),
            lattice.n_spins as u64
        );
    }

    #[test]
    fn cmr_observe_stops_before_cluster_actions() {
        let lattice = Lattice::new(vec![4, 4]);
        let couplings = vec![1.0; lattice.n_spins * lattice.n_neighbors];
        let mut observed = Realization::new(&lattice, couplings.clone(), &[1.5], 2, 23);
        let mut baseline = Realization::new(&lattice, couplings, &[1.5], 2, 23);
        let mut config = SimConfig {
            n_sweeps: 1,
            warmup_sweeps: 0,
            sweep_mode: SweepMode::Metropolis,
            cluster_update: None,
            pt_interval: None,
            pt_schedule: PtSchedule::SingleRandomEdge,
            overlap_cluster: None,
            autocorrelation_max_lag: None,
            autocorrelation_backend: AutocorrelationBackend::Ring,
            sequential: true,
            equilibration_diagnostic: false,
            conditional_link_overlap: true,
        };
        run_sweep_loop(
            &lattice,
            &mut baseline,
            2,
            1,
            &config,
            &AtomicBool::new(false),
            &|| {},
            0,
        )
        .unwrap();

        config.overlap_cluster = Some(OverlapClusterConfig {
            interval: 1,
            modes: vec![OverlapClusterBuildMode::Cmr],
            cluster_mode: ClusterMode::Sw,
            action: ClusterAction::Observe,
            collect_stats: true,
            snapshot_interval: None,
            max_temperature: None,
        });
        let result = run_sweep_loop(
            &lattice,
            &mut observed,
            2,
            1,
            &config,
            &AtomicBool::new(false),
            &|| {},
            0,
        )
        .unwrap();

        assert_eq!(observed.spins, baseline.spins);
        assert_eq!(observed.system_ids, baseline.system_ids);
        assert_eq!(
            result.per_disorder_cluster_observations[0]
                .cmr_blue
                .as_ref()
                .unwrap()
                .observation_count,
            vec![1]
        );
    }
}