Skip to content

Instantly share code, notes, and snippets.

@terasakisatoshi
Last active Jan 23, 2018
Embed
What would you like to do?
Ising Model Rust
[package]
name = "ising"
version = "0.1.0"
authors = ["SatoshiTeraaki <terasakisatoshi.math@gmail.com>"]
[dependencies]
ndarray = "0.11.0"
rand = "0.3"
ndarray-rand="0.7.0"
[profile.release]
opt-level = 3
debug = false
rpath = false
lto = true
debug-assertions = false
codegen-units = 1
panic = 'unwind'
extern crate rand;
extern crate ndarray_rand;
extern crate ndarray;
use std::f64;
use ndarray::prelude::*;
use rand::{random, Open01, XorShiftRng, Rng};
use std::time::Instant;
fn ising2d_sum_of_adjacent_spins(s: &Array2<i64>, m: usize, n: usize, i: usize, j: usize) -> i64 {
let i_bottom: usize;
let i_top: usize;
let j_right: usize;
let j_left: usize;
if i + 1 < m {
i_bottom = i + 1
} else {
i_bottom = 0
}
if i >= 1 { i_top = i - 1 } else { i_top = m - 1 }
if j + 1 < n {
j_right = j + 1
} else {
j_right = 0
}
if j >= 1 {
j_left = j - 1
} else {
j_left = n - 1
}
return s[(i_bottom, j)] + s[(i_top, j)] + s[(i, j_right)] + s[(i, j_left)];
}
fn ising2d_sweep(mut s: Array2<i64>, beta: f64, niter: i64) {
let m = s.shape()[0];
let n = s.shape()[1];
let prob = [
f64::exp(-2.0 * beta * (-4.0)),
f64::exp(-2.0 * beta * (-3.0)),
f64::exp(-2.0 * beta * (-2.0)),
f64::exp(-2.0 * beta * (-1.0)),
f64::exp(-2.0 * beta * (0.0)),
f64::exp(-2.0 * beta * (1.0)),
f64::exp(-2.0 * beta * (2.0)),
f64::exp(-2.0 * beta * (3.0)),
f64::exp(-2.0 * beta * (4.0)),
];
let mut rng = XorShiftRng::new_unseeded();
let iteration = niter / (m * n) as i64;
for _ in 0..iteration {
for i in 0..m {
for j in 0..n {
let s1 = s[(i, j)];
let k = s1 * ising2d_sum_of_adjacent_spins(&s, m, n, i, j);
//let Open01(val) = random::<Open01<f64>>();
let val = rng.gen_range(0., 1.);
if val < prob[(k + 4) as usize] {
s[(i, j)] = -(s1 as i64);
}
}
}
}
}
fn main() {
let n = 100;
let mut s = Array::<i64, _>::zeros(n * n).into_shape((n, n)).unwrap();
for i in 0..n {
for j in 0..n {
let Open01(val) = random::<Open01<f32>>();
if val < 0.5 {
s[(i, j)] = -1;
} else {
s[(i, j)] = 1;
}
}
}
let niter = 1000_000_000;
let beta = f64::ln(1.0 + f64::sqrt(2.0)) / 2.0;
let start = Instant::now();
ising2d_sweep(s, beta, niter);
let end = start.elapsed();
println!(
"{}.{:03}秒経過しました。",
end.as_secs(),
end.subsec_nanos() / 1_000_000
);
}
@terasakisatoshi
Copy link
Author

terasakisatoshi commented Jan 23, 2018

License is MIT

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment