terminals_core/substrate/

thermal.rs

//! ThermalProjection -- Ising-model spin state for substrate atoms.
//!
//! Each atom carries spin (+1/-1), local bias, inverse temperature,
//! local energy, and running magnetization. Together these enable
//! Gibbs sampling over the substrate lattice, complementing the
//! Kuramoto phase oscillator with a thermodynamic ensemble view.

use super::projection::{Projection, ProjectionId};

/// 5 floats: spin, bias, temperature, energy, magnetization = 20 bytes.
const THERMAL_BYTES: usize = 5 * 4;

/// Ising-model thermal projection for a single ComputeAtom.
#[derive(Debug, Clone, Copy)]
pub struct ThermalProjection {
    /// Spin state: +1.0 or -1.0.
    pub spin: f32,
    /// Local field bias b_i.
    pub bias: f32,
    /// Inverse temperature beta (higher = more ordered).
    pub temperature: f32,
    /// Local energy contribution.
    pub energy: f32,
    /// Running magnetization average.
    pub magnetization: f32,
}

impl Default for ThermalProjection {
    fn default() -> Self {
        Self {
            spin: 1.0,
            bias: 0.0,
            temperature: 1.0,
            energy: 0.0,
            magnetization: 0.0,
        }
    }
}

impl Projection for ThermalProjection {
    fn byte_size() -> usize {
        THERMAL_BYTES
    }

    fn id() -> ProjectionId {
        ProjectionId::Thermal
    }

    fn read(buf: &[u8]) -> Self {
        assert!(buf.len() >= THERMAL_BYTES);
        Self {
            spin: f32::from_le_bytes([buf[0], buf[1], buf[2], buf[3]]),
            bias: f32::from_le_bytes([buf[4], buf[5], buf[6], buf[7]]),
            temperature: f32::from_le_bytes([buf[8], buf[9], buf[10], buf[11]]),
            energy: f32::from_le_bytes([buf[12], buf[13], buf[14], buf[15]]),
            magnetization: f32::from_le_bytes([buf[16], buf[17], buf[18], buf[19]]),
        }
    }

    fn write(&self, buf: &mut [u8]) {
        assert!(buf.len() >= THERMAL_BYTES);
        buf[0..4].copy_from_slice(&self.spin.to_le_bytes());
        buf[4..8].copy_from_slice(&self.bias.to_le_bytes());
        buf[8..12].copy_from_slice(&self.temperature.to_le_bytes());
        buf[12..16].copy_from_slice(&self.energy.to_le_bytes());
        buf[16..20].copy_from_slice(&self.magnetization.to_le_bytes());
    }

    fn shape_hash_contribution(&self) -> u32 {
        let mut hash = 0x811c_9dc5u32;
        for byte in self.spin.to_bits().to_le_bytes() {
            hash ^= byte as u32;
            hash = hash.wrapping_mul(0x0100_0193);
        }
        for byte in self.magnetization.to_bits().to_le_bytes() {
            hash ^= byte as u32;
            hash = hash.wrapping_mul(0x0100_0193);
        }
        hash
    }
}

/// Perform a deterministic Gibbs step on a single spin site.
///
/// Computes the local energy delta for flipping spin_i, then accepts
/// with probability sigmoid(2 * beta * delta_E). Uses a deterministic
/// entropy seed instead of random sampling for reproducibility.
///
/// # Arguments
/// * `spin_i` - current spin (+1 or -1)
/// * `bias_i` - local field at site i
/// * `neighbor_coupling_sum` - sum over j: J_ij * s_j
/// * `beta` - inverse temperature
/// * `entropy_seed` - deterministic seed in [0, 1) for acceptance
///
/// # Returns
/// New spin value (+1.0 or -1.0) and local energy contribution.
pub fn gibbs_step(
    _spin_i: f32,
    bias_i: f32,
    neighbor_coupling_sum: f32,
    beta: f32,
    entropy_seed: f32,
) -> (f32, f32) {
    // Local field: h_i = b_i + sum_j J_ij * s_j
    let local_field = bias_i + neighbor_coupling_sum;

    // Energy contribution for spin_i = +1: E_up = -h_i
    // Energy contribution for spin_i = -1: E_down = +h_i
    // P(s_i = +1) = sigmoid(2 * beta * h_i) = 1 / (1 + exp(-2 * beta * h_i))
    let arg = 2.0 * beta * local_field;
    let prob_up = 1.0 / (1.0 + (-arg).exp());

    let seed = entropy_seed.clamp(0.0, 0.9999);
    let new_spin = if seed < prob_up { 1.0 } else { -1.0 };

    // Local energy: -s_i * h_i
    let energy = -new_spin * local_field;

    (new_spin, energy)
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn test_thermal_byte_size() {
        assert_eq!(ThermalProjection::byte_size(), 20);
    }

    #[test]
    fn test_thermal_roundtrip() {
        let proj = ThermalProjection {
            spin: -1.0,
            bias: 0.5,
            temperature: 2.0,
            energy: -0.3,
            magnetization: 0.7,
        };
        let mut buf = vec![0u8; ThermalProjection::byte_size()];
        proj.write(&mut buf);
        let restored = ThermalProjection::read(&buf);
        assert!((restored.spin - (-1.0)).abs() < 1e-6);
        assert!((restored.bias - 0.5).abs() < 1e-6);
        assert!((restored.temperature - 2.0).abs() < 1e-6);
        assert!((restored.energy - (-0.3)).abs() < 1e-6);
        assert!((restored.magnetization - 0.7).abs() < 1e-6);
    }

    #[test]
    fn test_thermal_default() {
        let proj = ThermalProjection::default();
        assert!((proj.spin - 1.0).abs() < 1e-6);
        assert!((proj.bias).abs() < 1e-6);
        assert!((proj.temperature - 1.0).abs() < 1e-6);
    }

    #[test]
    fn test_gibbs_step_strong_field_up() {
        // Strong positive field at low temperature -> spin should go to +1
        let (spin, energy) = gibbs_step(1.0, 10.0, 0.0, 5.0, 0.5);
        assert!((spin - 1.0).abs() < 1e-6, "spin = {}", spin);
        assert!(energy < 0.0, "energy should be negative for aligned spin");
    }

    #[test]
    fn test_gibbs_step_strong_field_down() {
        // Strong negative field -> spin should go to -1
        let (spin, _energy) = gibbs_step(1.0, -10.0, 0.0, 5.0, 0.5);
        assert!((spin - (-1.0)).abs() < 1e-6, "spin = {}", spin);
    }

    #[test]
    fn test_gibbs_step_deterministic() {
        // Same inputs -> same outputs
        let a = gibbs_step(1.0, 0.5, 0.3, 2.0, 0.4);
        let b = gibbs_step(1.0, 0.5, 0.3, 2.0, 0.4);
        assert!((a.0 - b.0).abs() < 1e-6);
        assert!((a.1 - b.1).abs() < 1e-6);
    }

    #[test]
    fn test_gibbs_step_entropy_clamp() {
        // Seed values at boundaries should not panic
        let _ = gibbs_step(1.0, 0.0, 0.0, 1.0, 0.0);
        let _ = gibbs_step(1.0, 0.0, 0.0, 1.0, 1.0);
        let _ = gibbs_step(1.0, 0.0, 0.0, 1.0, -0.5);
        let _ = gibbs_step(1.0, 0.0, 0.0, 1.0, 1.5);
    }

    #[test]
    fn test_shape_hash_changes_with_state() {
        let a = ThermalProjection {
            spin: 1.0,
            ..Default::default()
        };
        let b = ThermalProjection {
            spin: -1.0,
            ..Default::default()
        };
        assert_ne!(a.shape_hash_contribution(), b.shape_hash_contribution());
    }
}