oxicuda_seq/mrf/
gibbs.rs

1//! Gibbs sampling on the Ising MRF with optional simulated annealing.
2
3use super::mrf::IsingModel;
4use crate::error::{SeqError, SeqResult};
5use crate::handle::LcgRng;
6
7/// Gibbs sampling configuration.
8#[derive(Debug, Clone, Copy)]
9pub struct GibbsConfig {
10    /// Number of full sweeps (each site visited once per sweep).
11    pub n_sweeps: usize,
12    /// Burn-in sweeps to discard before recording statistics.
13    pub burn_in: usize,
14    /// Optional inverse-temperature schedule: linearly anneal from
15    /// `(beta_start, beta_end)` over the sweep budget.
16    pub anneal: Option<(f64, f64)>,
17}
18
19impl Default for GibbsConfig {
20    fn default() -> Self {
21        Self {
22            n_sweeps: 100,
23            burn_in: 20,
24            anneal: None,
25        }
26    }
27}
28
29/// Run Gibbs sampling on an Ising model.  Returns the final spin configuration
30/// and the post-burn-in mean magnetisation.
31pub fn ising_gibbs(
32    model: &IsingModel,
33    init: &[i32],
34    cfg: &GibbsConfig,
35    rng: &mut LcgRng,
36) -> SeqResult<(Vec<i32>, f64)> {
37    let n_sites = model.n_rows * model.n_cols;
38    if init.len() != n_sites {
39        return Err(SeqError::ShapeMismatch {
40            expected: n_sites,
41            got: init.len(),
42        });
43    }
44    for &s in init {
45        if s != 1 && s != -1 {
46            return Err(SeqError::InvalidParameter {
47                name: "spin".to_string(),
48                value: s as f64,
49            });
50        }
51    }
52    let mut spins = init.to_vec();
53    let mut mag_sum = 0.0;
54    let mut sample_count = 0usize;
55
56    for sweep in 0..cfg.n_sweeps {
57        let beta = match cfg.anneal {
58            Some((b0, b1)) => {
59                let frac = sweep as f64 / cfg.n_sweeps.max(1) as f64;
60                b0 + (b1 - b0) * frac
61            }
62            None => model.beta,
63        };
64        for r in 0..model.n_rows {
65            for c in 0..model.n_cols {
66                let mut nb_sum = 0i32;
67                if r > 0 {
68                    nb_sum += spins[(r - 1) * model.n_cols + c];
69                }
70                if r + 1 < model.n_rows {
71                    nb_sum += spins[(r + 1) * model.n_cols + c];
72                }
73                if c > 0 {
74                    nb_sum += spins[r * model.n_cols + (c - 1)];
75                }
76                if c + 1 < model.n_cols {
77                    nb_sum += spins[r * model.n_cols + (c + 1)];
78                }
79                let local_field = model.coupling * nb_sum as f64 + model.field;
80                let p_up = 1.0 / (1.0 + (-2.0 * beta * local_field).exp());
81                let u = rng.next_f64();
82                spins[r * model.n_cols + c] = if u < p_up { 1 } else { -1 };
83            }
84        }
85        if sweep >= cfg.burn_in {
86            mag_sum += model.magnetisation(&spins)?;
87            sample_count += 1;
88        }
89    }
90    let mean_mag = if sample_count == 0 {
91        0.0
92    } else {
93        mag_sum / sample_count as f64
94    };
95    Ok((spins, mean_mag))
96}
97
98#[cfg(test)]
99mod tests {
100    use super::*;
101
102    #[test]
103    fn gibbs_low_temp_polarises() {
104        // High β + strong coupling + positive field — magnetisation should be near +1.
105        let m = IsingModel::new(5, 5, 0.1, 1.0, 2.0).expect("ok");
106        let init = vec![1i32; 25];
107        let cfg = GibbsConfig {
108            n_sweeps: 200,
109            burn_in: 50,
110            anneal: None,
111        };
112        let mut rng = LcgRng::new(42);
113        let (_, mag) = ising_gibbs(&m, &init, &cfg, &mut rng).expect("ok");
114        assert!(mag > 0.5, "magnetisation {mag} should be > 0.5");
115    }
116}
oxicuda_seq/mrf/gibbs.rs

oxicuda_seq/mrf/
gibbs.rs
