ising/

lattice.rs

use rand::{
    distributions::{Bernoulli, Distribution, Uniform},
    Rng,
};

/// A lattice spin.
pub struct Spin {
    pub x: usize,
    pub y: usize,
}

/// The lattice holding the spin values.
pub struct Lattice {
    side_length: usize,
    /// Quadratic lattice -> Number of spins is the square of the side length.
    n_spins: i64,
    /// Inverse stored to prevent repeated casting to float and division.
    inverse_n_spins: f64,
    /// Number of positive spins used to calculate the magnetization.
    n_positive_spins: i64,
    /// Vector of spin values. True is a positive spin +1, false is a negative spin -1.
    vec: Vec<bool>,
    /// Random spin distribution to get random spins.
    random_spin_dist: Uniform<usize>,
}

impl Lattice {
    /// Creates a new lattice initialized randomly with the given
    /// random number generator `rng`.
    ///
    /// The lattice is quadratic with `side_length` as the side length.
    /// It is initialized with the same probability for positive and negative spins.
    ///
    /// # Example
    /// ```
    /// use rand::{rngs::SmallRng, SeedableRng};
    /// use ising::lattice::Lattice;
    ///
    /// let side_length = 8;
    /// let mut rng = SmallRng::seed_from_u64(42);
    /// let lattice = Lattice::new(side_length, &mut rng);
    ///```
    pub fn new(side_length: usize, rng: &mut impl Rng) -> Self {
        // Quadratic lattice.
        let n_spins = side_length * side_length;
        let mut vec = Vec::with_capacity(n_spins);

        // Initialize lattice with random spins.
        // with equal probability (0.5) for spin up +1 or spin down -1.
        // unwrap: 0.5 is obviously between 0 and 1.
        let dist = Bernoulli::new(0.5).unwrap();
        let mut n_positive_spins = 0;
        vec.resize_with(n_spins, || {
            let spin_value = dist.sample(rng);

            // Update the number of positive spins.
            if spin_value {
                n_positive_spins += 1;
            }

            spin_value
        });

        Self {
            side_length,
            n_spins: n_spins as i64,
            inverse_n_spins: 1.0 / n_spins as f64,
            n_positive_spins,
            vec,
            random_spin_dist: Uniform::new(0, side_length),
        }
    }

    /// The side length of the lattice.
    ///
    /// # Example
    /// ```
    /// # use rand::{rngs::SmallRng, SeedableRng};
    /// # use ising::lattice::Lattice;
    /// #
    /// let side_length = 8;
    /// # let mut rng = SmallRng::seed_from_u64(42);
    /// let lattice = Lattice::new(side_length, &mut rng);
    /// assert_eq!(side_length, lattice.side_length());
    ///```
    #[must_use]
    pub const fn side_length(&self) -> usize {
        self.side_length
    }

    /// Vector holding the spin values
    /// (`true` -> `+1`, `false` -> `-1`).
    ///
    /// Can be used to plot the lattice.
    ///
    /// # Example
    /// ```
    /// # use rand::{rngs::SmallRng, SeedableRng};
    /// # use ising::lattice::Lattice;
    /// #
    /// let side_length = 8;
    /// # let mut rng = SmallRng::seed_from_u64(42);
    /// let lattice = Lattice::new(side_length, &mut rng);
    /// let magnetization_sum =
    ///     lattice.vec().iter().fold(
    ///         0,
    ///         |mag_sum, value| if *value { mag_sum + 1 } else { mag_sum - 1 },
    ///     );
    /// assert_eq!(magnetization_sum, lattice.magnetization());
    ///```
    #[must_use]
    pub const fn vec(&self) -> &Vec<bool> {
        &self.vec
    }

    /// Returns a random spin.
    ///
    /// Since the returned spin is granted to be inside the lattice,
    /// [`value_unchecked`](Lattice::value_unchecked) and
    /// [`set_value_unchecked`](Lattice::set_value_unchecked) are safe
    /// to be used with the returned spin on the same lattice.
    ///
    /// # Example
    /// ```
    /// # use rand::{rngs::SmallRng, SeedableRng};
    /// # use ising::lattice::Lattice;
    ///
    /// let side_length = 8;
    /// # let mut rng = SmallRng::seed_from_u64(42);
    /// let mut lattice = Lattice::new(side_length, &mut rng);
    ///
    /// let (random_spin, random_spin_value) = lattice.rand_spin_with_value(&mut rng);
    ///
    /// // Safe to use the unchecked methods on the same lattice.
    /// unsafe { lattice.set_value_unchecked(&random_spin, !random_spin_value) };
    /// ```
    pub fn rand_spin_with_value(&self, rng: &mut impl Rng) -> (Spin, bool) {
        let spin = Spin {
            x: self.random_spin_dist.sample(rng),
            y: self.random_spin_dist.sample(rng),
        };

        // Safety: The random spin is granted to be inside the lattice.
        let value = unsafe { self.value_unchecked(&spin) };

        (spin, value)
    }

    /// Index of the given spin in the vector holding the spin values.
    const fn index(&self, spin: &Spin) -> usize {
        // Matrix indexing.
        spin.y * self.side_length + spin.x
    }

    /// The value of a spin (`true` -> `+1`, `false` -> `-1`).
    #[must_use]
    pub fn value(&self, spin: &Spin) -> bool {
        let ind = self.index(spin);

        self.vec[ind]
    }

    /// Same as [`value`](Self::value) but without bounds check.
    ///
    /// # Safety
    /// Calling this method with a spin out of the lattice is undefined behavior!
    ///
    /// Test your code with [`value`](Self::value) first and make sure
    /// that you access only valid lattice spins, then use this function
    /// if you need more performance.
    #[must_use]
    pub unsafe fn value_unchecked(&self, spin: &Spin) -> bool {
        let ind = self.index(spin);

        unsafe { *self.vec.get_unchecked(ind) }
    }

    // Update the number of positive spins.
    fn update_n_positive_spin(&mut self, new_val: bool) {
        if new_val {
            self.n_positive_spins += 1;
        } else {
            self.n_positive_spins -= 1;
        }
    }

    /// Sets the value of the spin.
    ///
    /// # Example
    /// ```
    /// # use rand::{rngs::SmallRng, SeedableRng};
    /// # use ising::lattice::Lattice;
    /// use ising::lattice::Spin;
    ///
    /// let side_length = 8;
    /// # let mut rng = SmallRng::seed_from_u64(42);
    /// let mut lattice = Lattice::new(side_length, &mut rng);
    /// let spin = Spin{x: 0, y: 0};
    /// let old_spin_value = lattice.value(&spin);
    ///
    /// lattice.set_value(&spin, !old_spin_value);
    /// assert_eq!(lattice.value(&spin), !old_spin_value);
    ///```
    pub fn set_value(&mut self, spin: &Spin, val: bool) {
        let ind = self.index(spin);
        self.vec[ind] = val;

        self.update_n_positive_spin(val);
    }

    /// Same as [`set_value`](Self::set_value) but without bounds check.
    ///
    /// # Safety
    /// Calling this method with a spin out of the lattice is undefined behavior!
    ///
    /// Like with [`value_unchecked`](Self::value_unchecked),
    /// test your code with [`set_value`](Self::set_value) first and make sure
    /// that you access only valid lattice spins, then use this function
    /// if you need more performance.
    pub unsafe fn set_value_unchecked(&mut self, spin: &Spin, val: bool) {
        let ind = self.index(spin);
        unsafe { *self.vec.get_unchecked_mut(ind) = val };

        self.update_n_positive_spin(val);
    }

    /// The previous and next coordinate to a given coordinate
    /// while applying the periodic boundary condition.
    const fn previous_next_coordinate(&self, coordinate: usize) -> (usize, usize) {
        let last_ind = self.side_length - 1;

        if coordinate == 0 {
            // First spin -> previous is last spin.
            (last_ind, 1)
        } else if coordinate == last_ind {
            // Last spin -> next is first spin.
            (last_ind - 1, 0)
        } else {
            (coordinate - 1, coordinate + 1)
        }
    }

    /// The four neighbouring spins of a given spin with the periodic boundary condition.
    ///
    /// If the center spin is inside the lattice, the neighbouring spins are also inside
    /// the lattice and safe to be used with the methods [`value_unchecked`](Self::value_unchecked)
    /// and [`set_value_unchecked`](Self::set_value_unchecked).
    #[must_use]
    pub const fn neighbouring_spins(&self, spin: &Spin) -> [Spin; 4] {
        let (previous_x, next_x) = self.previous_next_coordinate(spin.x);
        let (previous_y, next_y) = self.previous_next_coordinate(spin.y);

        // Respect order in memory for cache optimization.
        [
            Spin {
                x: spin.x,
                y: previous_y,
            },
            Spin {
                x: previous_x,
                y: spin.y,
            },
            Spin {
                x: next_x,
                y: spin.y,
            },
            Spin {
                x: spin.x,
                y: next_y,
            },
        ]
    }

    /// The current magnetization.
    ///
    /// The value is between `-N` and `N` (with step `2`) where `N` is the number of spins.
    ///
    /// # Example
    /// ```
    /// # use rand::{rngs::SmallRng, SeedableRng};
    /// # use ising::lattice::Lattice;
    /// # use ising::lattice::Spin;
    /// #
    /// let side_length = 8;
    /// # let mut rng = SmallRng::seed_from_u64(42);
    /// let mut lattice = Lattice::new(side_length, &mut rng);
    /// let spin = Spin{x: 0, y: 0};
    /// let old_value = lattice.value(&spin);
    /// let old_magnetization = lattice.magnetization();
    ///
    /// // Invert spin value
    /// lattice.set_value(&spin, !old_value);
    /// let new_value_sign = if old_value { -1 } else { 1 };
    /// assert_eq!(lattice.magnetization(), old_magnetization + 2 * new_value_sign);
    ///```
    ///
    /// See also [`average_magnetization`](Self::average_magnetization).
    #[must_use]
    pub const fn magnetization(&self) -> i64 {
        // magnetization = n_positive_spins - n_negative_spins
        //               = n_positive_spins - (n_spins - n_positive_spins)
        //               = 2 * n_positive_spins - n_spins
        2 * self.n_positive_spins - self.n_spins
    }

    /// The current magnetization averaged over the number of spins.
    ///
    /// The value is between `-1.0` and `1.0`.
    #[must_use]
    pub fn average_magnetization(&self) -> f64 {
        // average_magnetization = (n_positive_spins - n_negative_spins) / n_spins
        //                       = (n_positive_spins - (n_spins - n_positive_spins)) / n_spins
        //                       = (2 * n_positive_spins - n_spins) / n_spins
        //                       = 2 * n_positive_spins / n_spins - 1
        (2.0 * self.n_positive_spins as f64).mul_add(self.inverse_n_spins, -1.0)
    }
}