use ::ndarray::prelude::*;
use ::rand::prelude::*;
pub struct Lattice<R: RngCore> {
dims: (usize, usize),
n_of_spins: i32,
rng: R,
inner: Array2<i32>,
neighbors: Array2<[(usize, usize); 4]>,
}
impl<R: RngCore> Lattice<R> {
pub fn from_rng(size: (usize, usize), mut rng: R) -> Self {
let inner = Array2::from_shape_fn(size, |_| *[-1, 1].choose(&mut rng).unwrap());
Self::from_array_rng(inner, rng)
}
pub fn from_array_rng(array: Array2<i32>, rng: R) -> Self {
assert!(
array.iter().all(|spin| *spin == 1 || *spin == -1),
"Invalid spin value."
);
let roll_index = |ix: usize, amt: i32, max: usize| {
let max = max as i32;
((ix as i32 + amt + max) % max) as usize
};
let (width, height) = array.dim();
let neighbors = Array2::from_shape_fn((width, height), |ix| {
[
(roll_index(ix.0, 1, width), ix.1),
(ix.0, roll_index(ix.1, 1, height)),
(roll_index(ix.0, -1, width), ix.1),
(ix.0, roll_index(ix.1, -1, height)),
]
});
Lattice {
dims: (width, height),
inner: array,
n_of_spins: width as i32 * height as i32,
rng,
neighbors,
}
}
pub fn dims(&self) -> (usize, usize) {
self.dims
}
fn spin_times_all_neighbors(&self, ix: (usize, usize)) -> i32 {
self.inner[ix]
* self.neighbors[ix]
.iter()
.map(|n_ix| self.inner[*n_ix])
.sum::<i32>()
}
fn spin_times_two_neighbors(&self, ix: (usize, usize)) -> i32 {
self.inner[ix]
* self.neighbors[ix][0..2]
.iter()
.map(|n_ix| self.inner[*n_ix])
.sum::<i32>()
}
pub fn measure_E_diff(&self, (i, j): (usize, usize), J: f64) -> f64 {
2.0 * J * f64::from(self.spin_times_all_neighbors((i, j)))
}
pub fn measure_E(&self, J: f64) -> f64 {
-J * f64::from(
self.inner
.indexed_iter()
.map(|(ix, _)| self.spin_times_two_neighbors(ix))
.sum::<i32>(),
)
}
pub fn measure_I(&self) -> f64 {
f64::from(self.inner.sum().abs()) / f64::from(self.n_of_spins)
}
pub fn flip_spin(&mut self, ix: (usize, usize)) {
*self.inner.get_mut(ix).unwrap() *= -1;
}
pub fn gen_random_index(&mut self) -> (usize, usize) {
(
self.rng.gen_range(0, self.dims.0),
self.rng.gen_range(0, self.dims.1),
)
}
}
impl Lattice<SmallRng> {
pub fn new(dims: (usize, usize)) -> Self {
Self::from_rng(dims, SmallRng::from_entropy())
}
pub fn from_array(array: Array2<i32>) -> Self {
Self::from_array_rng(array, SmallRng::from_entropy())
}
}
#[cfg(test)]
mod test {
use ::pretty_assertions::assert_eq;
use super::*;
fn float_error(x: f64, t: f64) -> f64 {
(x - t).abs() / t
}
#[test]
fn test_lattice_from_rng() {
let rng = SmallRng::from_entropy();
let lattice = Lattice::from_rng((17, 10), rng);
assert_eq!(lattice.dims(), (17, 10));
}
#[test]
fn test_lattice_from_array_rng() {
let array = Array::from_shape_vec((2, 2), vec![1, -1, 1, -1]).unwrap();
let rng = SmallRng::from_entropy();
let lattice = Lattice::from_array_rng(array, rng);
assert_eq!(lattice.dims(), (2, 2));
}
#[test]
fn test_lattice_new() {
let lattice = Lattice::new((40, 20));
assert_eq!(lattice.dims(), (40, 20));
}
#[test]
fn test_lattice_from_array() {
let array = Array::from_shape_vec((2, 2), vec![1, -1, 1, -1]).unwrap();
let lattice = Lattice::from_array(array);
assert_eq!(lattice.dims(), (2, 2));
}
#[test]
fn test_spin_times_neighbors() {
let spins = [-1, -1, 1, 1, 1, 1, 1, 1, -1];
let array = Array::from_shape_vec((3, 3), spins.to_vec()).unwrap();
let lattice = Lattice::from_array(array);
let product = lattice.spin_times_all_neighbors((1, 1));
assert_eq!(product, 2);
}
#[test]
fn test_measure_E_difference() {
let array =
Array::from_shape_vec((3, 3), vec![-1, -1, 1, 1, 1, 1, -1, 1, 1])
.unwrap();
let lattice = Lattice::from_array(array);
let J = 1.0;
let E_diff = lattice.measure_E_diff((1, 1), J);
assert_eq!(E_diff, 4.0);
}
#[test]
fn test_measure_E() {
let array =
Array::from_shape_vec((3, 3), vec![-1, -1, -1, 1, 1, -1, 1, 1, -1])
.unwrap();
let lattice = Lattice::from_array(array);
let J = 1.0;
let E = lattice.measure_E(J);
assert_eq!(E, -2.0);
}
#[test]
fn test_measure_I() {
let array = Array::from_shape_vec((2, 2), vec![-1, -1, -1, 1]).unwrap();
let lattice = Lattice::from_array(array);
let I = lattice.measure_I();
assert_eq!(I, 0.5);
}
#[test]
fn test_flip_spin() {
let array = Array::from_shape_vec(
(3, 3),
vec![-1, -1, -1, -1, 1, 1, -1, -1, 1],
)
.unwrap();
let mut lattice = Lattice::from_array(array);
let J = 1.0;
let E_1 = lattice.measure_E(J);
lattice.flip_spin((1, 1));
let E_2 = lattice.measure_E(J);
assert!(float_error(E_2 - E_1, -4.0) < 0.01);
}
}