extract τ per realization in Rust, average τ across disorder (literature-standard)

Signed-off-by: PeaBrane <yanrpei@gmail.com>

Author: PeaBrane

python/peapods/spin_models.py

@@ -226,58 +226,13 @@ def sample(
         if "top_cluster_sizes" in result:
             self.top_cluster_sizes = result["top_cluster_sizes"]
 
-        if "mags2_autocorr" in result:
-            self.mags2_autocorr = result["mags2_autocorr"]
-        if "overlap2_autocorr" in result:
-            self.overlap2_autocorr = result["overlap2_autocorr"]
+        if "mags2_tau" in result:
+            self.mags2_tau = result["mags2_tau"]
+        if "overlap2_tau" in result:
+            self.overlap2_tau = result["overlap2_tau"]
 
         return result
 
-    @staticmethod
-    def _sokal_tau(gamma):
-        """Integrated autocorrelation time via Sokal's self-consistent windowing.
-
-        Args:
-            gamma: 1D array, normalized autocorrelation Γ(δ) with Γ[0]=1.
-
-        Returns:
-            τ_int estimate (float).
-        """
-        tau = 0.5
-        for w in range(1, len(gamma)):
-            tau += gamma[w]
-            if w >= 5 * tau:
-                return tau
-        return tau
-
-    def mags2_autocorrelation_time(self):
-        """Integrated autocorrelation time of m² per temperature.
-
-        Returns:
-            Array of shape ``(n_temps,)``.
-        """
-        if not hasattr(self, "mags2_autocorr"):
-            raise RuntimeError(
-                "No autocorrelation data; call sample() with autocorrelation_max_lag first"
-            )
-        return np.array(
-            [self._sokal_tau(self.mags2_autocorr[t]) for t in range(self.n_temps)]
-        )
-
-    def overlap2_autocorrelation_time(self):
-        """Integrated autocorrelation time of q² per temperature. Requires n_replicas >= 2.
-
-        Returns:
-            Array of shape ``(n_temps,)``.
-        """
-        if not hasattr(self, "overlap2_autocorr"):
-            raise RuntimeError(
-                "No overlap autocorrelation data; call sample() with autocorrelation_max_lag and n_replicas >= 2"
-            )
-        return np.array(
-            [self._sokal_tau(self.overlap2_autocorr[t]) for t in range(self.n_temps)]
-        )
-
     def get_energies(self):
         """Return the mean energies per temperature from the last sample run."""
         return self.energies_avg

python/peapods/sweep.py

@@ -50,10 +50,10 @@ def _save_data(models, config_label, temperatures, output_dir):
             save_dict[f"{prefix}_mean_cluster_size"] = model.mean_cluster_size
         if hasattr(model, "top_cluster_sizes"):
             save_dict[f"{prefix}_top_cluster_sizes"] = model.top_cluster_sizes
-        if hasattr(model, "mags2_autocorr"):
-            save_dict[f"{prefix}_mags2_tau"] = model.mags2_autocorrelation_time()
-        if hasattr(model, "overlap2_autocorr"):
-            save_dict[f"{prefix}_overlap2_tau"] = model.overlap2_autocorrelation_time()
+        if hasattr(model, "mags2_tau"):
+            save_dict[f"{prefix}_mags2_tau"] = model.mags2_tau
+        if hasattr(model, "overlap2_tau"):
+            save_dict[f"{prefix}_overlap2_tau"] = model.overlap2_tau
 
     path = Path(output_dir) / f"sweep_{config_label}.npz"
     np.savez(path, **save_dict)
@@ -139,9 +139,9 @@ def _plot_autocorrelation_time(all_results, temperatures, plot_temp, output_dir)
     else:
         t_idx = None
 
-    for obs_name, attr, tau_method in [
-        ("m2", "mags2_autocorr", "mags2_autocorrelation_time"),
-        ("q2", "overlap2_autocorr", "overlap2_autocorrelation_time"),
+    for obs_name, attr in [
+        ("m2", "mags2_tau"),
+        ("q2", "overlap2_tau"),
     ]:
         has_any = any(
             hasattr(m, attr) for models in all_results.values() for m in models.values()
@@ -156,7 +156,7 @@ def _plot_autocorrelation_time(all_results, temperatures, plot_temp, output_dir)
             for size_label, model in models.items():
                 if not hasattr(model, attr):
                     continue
-                tau_arr = getattr(model, tau_method)()
+                tau_arr = getattr(model, attr)
                 L = max(model.lattice_shape)
                 sizes_L.append(L)
                 if t_idx is not None:

spin-sim/src/simulation.rs

@@ -88,6 +88,17 @@ impl AutocorrAccum {
     }
 }
 
+fn sokal_tau(gamma: &[f64]) -> f64 {
+    let mut tau = 0.5;
+    for (w, &g) in gamma.iter().enumerate().skip(1) {
+        tau += g;
+        if w as f64 >= 5.0 * tau {
+            return tau;
+        }
+    }
+    tau
+}
+
 /// Mutable state for one disorder realization.
 ///
 /// Holds the coupling array (fixed after construction), spin configurations for
@@ -603,13 +614,13 @@ pub fn run_sweep_loop(
         vec![]
     };
 
-    let mags2_autocorrelation = m2_accum
+    let mags2_tau = m2_accum
         .as_ref()
-        .map(|acc| acc.finish())
+        .map(|acc| acc.finish().iter().map(|g| sokal_tau(g)).collect())
         .unwrap_or_default();
-    let overlap2_autocorrelation = q2_accum
+    let overlap2_tau = q2_accum
         .as_ref()
-        .map(|acc| acc.finish())
+        .map(|acc| acc.finish().iter().map(|g| sokal_tau(g)).collect())
         .unwrap_or_default();
 
     Ok(SweepResult {
@@ -636,8 +647,8 @@ pub fn run_sweep_loop(
         fk_csd: fk_csd_accum,
         overlap_csd: overlap_csd_accum,
         top_cluster_sizes,
-        mags2_autocorrelation,
-        overlap2_autocorrelation,
+        mags2_tau,
+        overlap2_tau,
     })
 }
 

spin-sim/src/statistics/results.rs

@@ -26,12 +26,12 @@ pub struct SweepResult {
     /// Average relative size of k-th largest blue cluster per temperature.
     /// Shape: [n_temps][4]. Empty if collect_top_clusters=false.
     pub top_cluster_sizes: Vec<[f64; 4]>,
-    /// Normalized autocorrelation Γ(Δt) of m², shape [n_temps][max_lag+1].
+    /// Integrated autocorrelation time τ_int(m²) per temperature.
     /// Empty if autocorrelation_max_lag is None.
-    pub mags2_autocorrelation: Vec<Vec<f64>>,
-    /// Normalized autocorrelation Γ(Δt) of q², shape [n_temps][max_lag+1].
+    pub mags2_tau: Vec<f64>,
+    /// Integrated autocorrelation time τ_int(q²) per temperature.
     /// Empty if autocorrelation_max_lag is None or n_replicas < 2.
-    pub overlap2_autocorrelation: Vec<Vec<f64>>,
+    pub overlap2_tau: Vec<f64>,
 }
 
 impl SweepResult {
@@ -48,16 +48,8 @@ impl SweepResult {
 
         let n_top = results[0].top_cluster_sizes.len();
 
-        let m2_ac_ntemps = results[0].mags2_autocorrelation.len();
-        let m2_ac_len = results[0]
-            .mags2_autocorrelation
-            .first()
-            .map_or(0, |v| v.len());
-        let q2_ac_ntemps = results[0].overlap2_autocorrelation.len();
-        let q2_ac_len = results[0]
-            .overlap2_autocorrelation
-            .first()
-            .map_or(0, |v| v.len());
+        let m2_tau_len = results[0].mags2_tau.len();
+        let q2_tau_len = results[0].overlap2_tau.len();
 
         let mut agg = SweepResult {
             mags: vec![0.0; n_temps],
@@ -71,8 +63,8 @@ impl SweepResult {
             fk_csd: (0..n_fk_csd).map(|_| vec![0u64; fk_len]).collect(),
             overlap_csd: (0..n_ov_csd).map(|_| vec![0u64; ov_len]).collect(),
             top_cluster_sizes: vec![[0.0; 4]; n_top],
-            mags2_autocorrelation: (0..m2_ac_ntemps).map(|_| vec![0.0; m2_ac_len]).collect(),
-            overlap2_autocorrelation: (0..q2_ac_ntemps).map(|_| vec![0.0; q2_ac_len]).collect(),
+            mags2_tau: vec![0.0; m2_tau_len],
+            overlap2_tau: vec![0.0; q2_tau_len],
         };
 
         for r in results {
@@ -119,23 +111,11 @@ impl SweepResult {
                     a[k] += s[k];
                 }
             }
-            for (a_row, r_row) in agg
-                .mags2_autocorrelation
-                .iter_mut()
-                .zip(r.mags2_autocorrelation.iter())
-            {
-                for (a, &v) in a_row.iter_mut().zip(r_row.iter()) {
-                    *a += v;
-                }
+            for (a, &v) in agg.mags2_tau.iter_mut().zip(r.mags2_tau.iter()) {
+                *a += v;
             }
-            for (a_row, r_row) in agg
-                .overlap2_autocorrelation
-                .iter_mut()
-                .zip(r.overlap2_autocorrelation.iter())
-            {
-                for (a, &v) in a_row.iter_mut().zip(r_row.iter()) {
-                    *a += v;
-                }
+            for (a, &v) in agg.overlap2_tau.iter_mut().zip(r.overlap2_tau.iter()) {
+                *a += v;
             }
         }
 
@@ -159,15 +139,11 @@ impl SweepResult {
             }
         }
 
-        for row in agg.mags2_autocorrelation.iter_mut() {
-            for v in row.iter_mut() {
-                *v /= n;
-            }
+        for v in agg.mags2_tau.iter_mut() {
+            *v /= n;
         }
-        for row in agg.overlap2_autocorrelation.iter_mut() {
-            for v in row.iter_mut() {
-                *v /= n;
-            }
+        for v in agg.overlap2_tau.iter_mut() {
+            *v /= n;
         }
 
         agg

src/lib.rs

@@ -270,19 +270,15 @@ impl IsingSimulation {
             dict.set_item("top_cluster_sizes", arr)?;
         }
 
-        if !agg.mags2_autocorrelation.is_empty() {
-            let k = agg.mags2_autocorrelation[0].len();
-            let arr = Array2::from_shape_fn((n_temps, k), |(t, d)| agg.mags2_autocorrelation[t][d])
-                .into_pyarray(py);
-            dict.set_item("mags2_autocorr", arr)?;
+        if !agg.mags2_tau.is_empty() {
+            dict.set_item("mags2_tau", Array1::from(agg.mags2_tau).into_pyarray(py))?;
         }
 
-        if !agg.overlap2_autocorrelation.is_empty() {
-            let k = agg.overlap2_autocorrelation[0].len();
-            let arr =
-                Array2::from_shape_fn((n_temps, k), |(t, d)| agg.overlap2_autocorrelation[t][d])
-                    .into_pyarray(py);
-            dict.set_item("overlap2_autocorr", arr)?;
+        if !agg.overlap2_tau.is_empty() {
+            dict.set_item(
+                "overlap2_tau",
+                Array1::from(agg.overlap2_tau).into_pyarray(py),
+            )?;
         }
 
         Ok(dict)