use av1 for encoding

This reverts commit 650edff95c.

.gitignore

@@ -1,7 +1,6 @@
 /target
-/tmp
-/test.mp4
+/output.mp4
 /.vscode
 
 # CLion
-**/.idea
+**/.idea

Cargo.toml

@@ -2,33 +2,26 @@
 name = "physarum"
 version = "0.1.0"
 authors = ["Simon Gardling <titaniumtown@gmail.com>", "mindv0rtex <mindv0rtex@users.noreply.github.com>"]
-edition = "2018"
+edition = "2021"
 
 [dependencies]
-image = "0.23.14"
-indicatif = {version = "0.15.0", features = ["rayon"]}
-itertools = "0.10"
-rand = "0.8.3"
-rand_distr = "0.4.0"
-rayon = {git = "https://github.com/rayon-rs/rayon.git"}
-fastapprox = "0.3.0"
+image = "0.25"
+indicatif = { version = "0.17", features = [ "rayon" ] }
+itertools = "0.14"
+rand = "0.9"
+rand_distr = "0.5"
+rayon = "1.10"
+fastapprox = "0.3"
 
 [dev-dependencies]
-criterion = "0.3.4"
+criterion = { version = "0.5", features = ["html_reports"] }
 
-[patch.crates-io]
-rayon = {git = "https://github.com/rayon-rs/rayon.git"} #fixes dependency issues
-
-[profile.dev]
-codegen-units = 1
-opt-level = 3
-target-cpu = "native"
-lto = "fat"
-debug = true
+[[bench]]
+name = "benchmark"
+harness = false
 
 [profile.release]
 codegen-units = 1
 opt-level = 3
-target-cpu = "native"
 lto = "fat"
-debug = false
+debug = false

Notes.md

@@ -14,4 +14,7 @@
 - fast_approx::fast::sin + fast_approx::fast::cos
     - 0.000018658ms (2.56% slower)
 - fast_approx::faster::sin + fast_approx::faster::cos
-    - 0.000015878ms (14.57% faster)
+    - 0.000015878ms (14.57% faster)
+## sin and cos lookup table:
+- With lookup tables: 0.000044201500713825226ms
+- Without lookup tables: 0.000043705105781555176ms`

README.md

@@ -1,14 +1,6 @@
-![Physarum distribution](assets/example.gif)
+# Physarum simulation
+
+based on: https://github.com/yan-zaretskiy/physarum
 
 
-I recently came across [this repo](https://github.com/fogleman/physarum) that uses Go to implement an extended version of the Physarum transport networks model, first described in:
-
-> Jones, J. (2010). Characteristics of pattern formation and evolution in approximations of physarum transport networks. Artificial Life, 16(2), 127-153. https://doi.org/10.1162/artl.2010.16.2.16202
-
-I was amazed by the images this model produces, so I decided to figure out how it works and reimplement it in Rust for no other reason than to have some fun and practice some Rust.
-
-# Physarum Simulation
-
-The simulation consists of a 2D grid that stores information about the spatial distribution of some chemotactic sensory stimuli. This is called the _trail_ layer in the original paper. The stimuli attract discrete agents that are meant to represent a particle of Physarum plasmodium gel-sol structure. Each agent knows its position and orientation on the said 2D grid. They are attracted to the stimuli. They have sensors to tell them which way they should turn and move to reach the location with the highest level of the thing that attracts them. As the agents move, they leave a trail of the same stimuli behind, effectively signaling other agents where they should follow. The trail layer is also subjected to a simple diffusion-decay operator after every agent got their change to move to a new location. This image (shamelessly taken from [here](https://sagejenson.com/physarum)) shows all the steps that happen during one iteration of the simulation:
-
-![Physarum simulation steps](assets/physarum-steps.jpg)
+![](./assets/test.gif)

benches/benchmark.rs

@@ -0,0 +1,107 @@
+use criterion::{criterion_group, criterion_main, BatchSize, BenchmarkId, Criterion};
+use physarum::{
+    agent::Agent,
+    grid::{combine, Grid},
+    model,
+};
+use rand::{rngs::StdRng, SeedableRng};
+
+// Benchmark agent movement and deposition
+fn agent_benchmark(c: &mut Criterion) {
+    let mut group = c.benchmark_group("Agent Tick");
+    let n_agents = [1_000, 10_000, 100_000];
+
+    for &n in &n_agents {
+        group.bench_with_input(BenchmarkId::from_parameter(n), &n, |b, &n| {
+            let mut rng = StdRng::seed_from_u64(42);
+            let agents = (0..n).map(|_| Agent::new(256, 256, &mut rng)).collect();
+            let mut grid = Grid::new(256, 256, &mut rng, agents);
+
+            b.iter(|| {
+                grid.tick();
+            });
+        });
+    }
+    group.finish();
+}
+
+// Benchmark grid diffusion (blur)
+fn diffusion_benchmark(c: &mut Criterion) {
+    let mut group = c.benchmark_group("Grid Diffusion");
+    let sizes = [(256, 256), (512, 512)];
+    let radii = [1, 3];
+
+    for &(w, h) in &sizes {
+        for &r in &radii {
+            group.bench_with_input(
+                BenchmarkId::new("diffuse", format!("{}x{}_r{}", w, h, r)),
+                &(w, h, r),
+                |b, &(w, h, r)| {
+                    b.iter_batched(
+                        || {
+                            let mut rng = StdRng::seed_from_u64(42);
+                            Grid::new(w, h, &mut rng, vec![])
+                        },
+                        |mut grid| grid.diffuse(r),
+                        BatchSize::SmallInput,
+                    );
+                },
+            );
+        }
+    }
+    group.finish();
+}
+
+// Benchmark grid combining
+fn combine_benchmark(c: &mut Criterion) {
+    let mut group = c.benchmark_group("Combine Grids");
+    let populations = [2, 4];
+
+    for &np in &populations {
+        group.bench_with_input(BenchmarkId::from_parameter(np), &np, |b, &np| {
+            b.iter_batched(
+                || {
+                    let mut rng = StdRng::seed_from_u64(42);
+                    let grids = (0..np)
+                        .map(|_| Grid::new(256, 256, &mut rng, vec![]))
+                        .collect::<Vec<_>>();
+                    let attraction_table = vec![vec![1.0; np]; np];
+                    (grids, attraction_table)
+                },
+                |(mut grids, table)| combine(&mut grids, &table),
+                BatchSize::SmallInput,
+            );
+        });
+    }
+    group.finish();
+}
+
+// Benchmark full model step
+fn model_step_benchmark(c: &mut Criterion) {
+    let mut group = c.benchmark_group("Model Step");
+    let params = [(256, 256, 2), (512, 512, 4)];
+
+    for &(w, h, np) in &params {
+        group.bench_with_input(
+            BenchmarkId::new("step", format!("{}x{}_p{}", w, h, np)),
+            &(w, h, np),
+            |b, &(w, h, np)| {
+                b.iter_batched(
+                    || model::Model::new(w, h, 1 << 16, np, 1),
+                    |mut model| model.step(),
+                    BatchSize::SmallInput,
+                );
+            },
+        );
+    }
+    group.finish();
+}
+
+criterion_group!(
+    benches,
+    agent_benchmark,
+    diffusion_benchmark,
+    combine_benchmark,
+    model_step_benchmark
+);
+criterion_main!(benches);

rust-toolchain

@@ -1 +0,0 @@
-nightly

rustfmt.toml

@@ -1,3 +0,0 @@
-unstable_features = true
-imports_granularity = "Crate"
-wrap_comments = false

src/agent.rs

@@ -1,84 +1,48 @@
-use crate::{
-    util::wrap,
-    buffer::Buf,
-};
-
-use rand::{seq::SliceRandom, Rng};
-use std::f32::consts::TAU;
+use crate::grid::PopulationConfig;
+use crate::{buffer::Buf, util::wrap};
 use fastapprox::faster::{cos, sin};
+use rand::prelude::IndexedRandom;
+use rand::Rng;
+use std::f32::consts::TAU;
 use std::fmt::{Display, Formatter};
 
-// A single Physarum agent. The x and y positions are continuous, hence we use floating point numbers instead of integers.
-#[derive(Debug)]
+/// A single Physarum agent. The x and y positions are continuous, hence we use floating point numbers instead of integers.
+#[derive(Debug, Clone, PartialEq)]
 pub struct Agent {
     pub x: f32,
     pub y: f32,
-    pub angle: f32,
-    pub population_id: usize,
-    pub i: usize,
+    heading: f32,
 }
 
 impl Display for Agent {
     fn fmt(&self, f: &mut Formatter<'_>) -> std::fmt::Result {
-        write!(
-            f,
-            "{{\n(x,y): ({},{})\nangle: {}\npopulation id: {}\ni: {}}}",
-            self.x, self.y, self.angle, self.population_id, self.i,
-        )
+        write!(f, "{:?}", self)
     }
 }
 
 impl Agent {
-    // Construct a new agent with random parameters.
-    pub fn new<R: Rng + ?Sized>(width: usize, height: usize, id: usize, rng: &mut R, i: usize) -> Self {
-        let (x, y, angle) = rng.gen::<(f32, f32, f32)>();
+    /// Construct a new agent with random parameters.
+    pub fn new<R: Rng + ?Sized>(width: usize, height: usize, rng: &mut R) -> Self {
+        let (x, y, angle) = rng.random::<(f32, f32, f32)>();
         Agent {
             x: x * width as f32,
             y: y * height as f32,
-            angle: angle * TAU,
-            population_id: id,
-            i,
+            heading: angle * TAU,
         }
     }
 
-    pub fn get_from_table(i: f32, table: &Vec<f32>) -> f32 {
-        let interval: f32 = 1.0/0.01;
+    /// Tick an agent
+    pub fn tick(&mut self, buf: &Buf, pop_config: PopulationConfig, width: usize, height: usize) {
+        let xc = self.x + cos(self.heading) * pop_config.sensor_distance;
+        let yc = self.y + sin(self.heading) * pop_config.sensor_distance;
 
-        let i_proc = (i * interval).round()+(360.0*interval);
-        let output = table[i_proc as usize];
+        let agent_add_sens = self.heading + pop_config.sensor_angle;
+        let agent_sub_sens = self.heading - pop_config.sensor_angle;
 
-        return output;
-    }
-
-    // Tick an agent
-    #[inline]
-    pub fn tick(&mut self, buf: &Buf, sensor_distance: f32, sensor_angle: f32, rotation_angle: f32, step_distance: f32, width: usize, height: usize, sin_table: &Vec<f32>, cos_table: &Vec<f32>) {
-        /*
-        let xc = self.x + cos(self.angle) * sensor_distance;
-        let yc = self.y + sin(self.angle) * sensor_distance;
-        */
-
-        // /*
-        let xc = self.x + Self::get_from_table(self.angle, &cos_table) * sensor_distance;
-        let yc = self.y + Self::get_from_table(self.angle, &sin_table) * sensor_distance;
-        // */
-
-        let agent_add_sens = self.angle + sensor_angle;
-        let agent_sub_sens = self.angle - sensor_angle;
-
-        // /*
-        let xl = self.x + Self::get_from_table(agent_sub_sens, &cos_table) * sensor_distance;
-        let yl = self.y + Self::get_from_table(agent_sub_sens, &sin_table) * sensor_distance;
-        let xr = self.x + Self::get_from_table(agent_add_sens, &cos_table) * sensor_distance;
-        let yr = self.y + Self::get_from_table(agent_add_sens, &sin_table) * sensor_distance;
-        // */
-
-        /*
-        let xl = self.x + cos(agent_sub_sens) * sensor_distance;
-        let yl = self.y + sin(agent_sub_sens) * sensor_distance;
-        let xr = self.x + cos(agent_add_sens) * sensor_distance;
-        let yr = self.y + sin(agent_add_sens) * sensor_distance;
-        */
+        let xl = self.x + cos(agent_sub_sens) * pop_config.sensor_distance;
+        let yl = self.y + sin(agent_sub_sens) * pop_config.sensor_distance;
+        let xr = self.x + cos(agent_add_sens) * pop_config.sensor_distance;
+        let yr = self.y + sin(agent_add_sens) * pop_config.sensor_distance;
 
         // We sense from the buffer because this is where we previously combined data from all the grid.
         let center = buf.get_buf(xc, yc);
@@ -86,53 +50,30 @@ impl Agent {
         let right = buf.get_buf(xr, yr);
 
         // Rotate and move logic
-        let mut rng = rand::thread_rng();
-        let mut direction: f32 = 0.0;
-
-        if (center > left) && (center > right) {
-            direction = 0.0;
+        let direction = if (center > left) && (center > right) {
+            0.0
         } else if (center < left) && (center < right) {
-            direction = *[-1.0, 1.0].choose(&mut rng).unwrap();
+            *[-1.0, 1.0]
+                .choose(&mut rand::rng())
+                .expect("unable to choose random direction")
         } else if left < right {
-            direction = 1.0;
+            1.0
         } else if right < left {
-            direction = -1.0;
-        }
+            -1.0
+        } else {
+            0.0
+        };
 
-        let delta_angle = rotation_angle * direction;
+        let delta_angle = pop_config.rotation_angle * direction;
 
-        self.angle = wrap(self.angle + delta_angle, TAU);
+        self.heading = wrap(self.heading + delta_angle, TAU);
         self.x = wrap(
-            // self.x + step_distance * cos(self.angle),
-            self.x + step_distance * Self::get_from_table(self.angle, &cos_table),
+            self.x + pop_config.step_distance * cos(self.heading),
             width as f32,
         );
         self.y = wrap(
-            // self.y + step_distance * sin(self.angle),
-            self.y + step_distance * Self::get_from_table(self.angle, &sin_table),
+            self.y + pop_config.step_distance * sin(self.heading),
             height as f32,
         );
     }
 }
-
-impl Clone for Agent {
-    fn clone(&self) -> Agent {
-        Agent {
-            x: self.x,
-            y: self.y,
-            angle: self.angle,
-            population_id: self.population_id,
-            i: self.i,
-        }
-    }
-}
-
-impl PartialEq for Agent {
-    fn eq(&self, other: &Self) -> bool {
-        self.x == other.x
-            && self.y == other.y
-            && self.angle == other.angle
-            && self.population_id == other.population_id
-            && self.i == other.i
-    }
-}

src/blur.rs

@@ -1,19 +1,11 @@
 use itertools::multizip;
 use rayon::prelude::*;
 
-#[derive(Debug)]
+#[derive(Debug, Clone)]
 pub struct Blur {
     row_buffer: Vec<f32>,
 }
 
-impl Clone for Blur {
-    fn clone(&self) -> Blur {
-        Blur {
-            row_buffer: self.row_buffer.clone(),
-        }
-    }
-}
-
 impl Blur {
     pub fn new(width: usize) -> Self {
         Blur {
@@ -21,7 +13,7 @@ impl Blur {
         }
     }
 
-    // Blur an image with 2 box filter passes. The result will be written to the src slice, while the buf slice is used as a scratch space.
+    /// Blur an image with 2 box filter passes. The result will be written to the src slice, while the buf slice is used as a scratch space.
     pub fn run(
         &mut self,
         src: &mut [f32],
@@ -36,23 +28,22 @@ impl Blur {
         self.box_blur(src, buf, width, height, boxes[1], decay);
     }
 
-    // Approximate 1D Gaussian filter of standard deviation sigma with N box filter passes. Each element in the output array contains the radius of the box filter for the corresponding pass.
-    fn boxes_for_gaussian<const N: usize>(sigma: f32) -> ([usize; N]) {
-        let w_ideal = (12.0 * sigma * sigma / N as f32 + 1.0).sqrt();
+    /// Approximate 1D Gaussian filter of standard deviation sigma with N box filter passes. Each element in the output array contains the radius of the box filter for the corresponding pass.
+    fn boxes_for_gaussian<const N: usize>(sigma: f32) -> [usize; N] {
+        let sigma_sq = sigma.powi(2);
+        let w_ideal = (12.0 * sigma_sq / N as f32 + 1.0).sqrt();
         let mut w = w_ideal as usize;
         w -= 1 - (w & 1);
-        let mut m = 0.25 * (N * (w + 3)) as f32;
-        m -= 3.0 * sigma * sigma / (w + 1) as f32;
-        let m = m.round() as usize;
+        let m = (0.25 * (N * (w + 3)) as f32 - 3.0 * sigma_sq / (w + 1) as f32).round() as usize;
 
-        let mut result = [0; N];
-        for (i, value) in result.iter_mut().enumerate() {
-            *value = (if i < m { w - 1 } else { w + 1 }) / 2;
-        }
-        result
+        (0..N)
+            .map(|i| (w + 1 - 2 * (i < m) as usize) / 2)
+            .collect::<Vec<_>>()
+            .try_into()
+            .unwrap()
     }
 
-    // Perform one pass of the 2D box filter of the given radius. The result will be written to the src slice, while the buf slice is used as a scratch space.
+    /// Perform one pass of the 2D box filter of the given radius. The result will be written to the src slice, while the buf slice is used as a scratch space.
     fn box_blur(
         &mut self,
         src: &mut [f32],
@@ -66,7 +57,7 @@ impl Blur {
         self.box_blur_v(buf, src, width, height, radius, decay);
     }
 
-    // Perform one pass of the 1D box filter of the given radius along x axis.
+    /// Perform one pass of the 1D box filter of the given radius along x axis.
     fn box_blur_h(&mut self, src: &[f32], dst: &mut [f32], width: usize, radius: usize) {
         let weight = 1.0 / (2 * radius + 1) as f32;
 
@@ -75,7 +66,7 @@ impl Blur {
             .for_each(|(src_row, dst_row)| {
                 // First we build a value for the beginning of each row. We assume periodic boundary conditions, so we need to push the left index to the opposite side of the row.
                 let width_sub_radius = width - radius;
-                let mut value = src_row[width - radius - 1];
+                let mut value = src_row[width_sub_radius - 1];
                 for j in 0..radius {
                     value += src_row[width_sub_radius + j] + src_row[j];
                 }
@@ -89,7 +80,7 @@ impl Blur {
             })
     }
 
-    // Perform one pass of the 1D box filter of the given radius along y axis. Applies the decay factor to the destination buffer.
+    /// Perform one pass of the 1D box filter of the given radius along y axis. Applies the decay factor to the destination buffer.
     fn box_blur_v(
         &mut self,
         src: &[f32],
@@ -148,15 +139,69 @@ mod tests {
         // ndimage.uniform_filter(a, size=3, mode='wrap') # 2D blur
 
         let mut src: Vec<f32> = vec![
-            0.32352856, 0.06571674, 0.01939427, 0.06352045, 0.708_527, 0.617_221_7, 0.16638431,
-            0.628_400_74, 0.554_893_9, 0.240_076_77, 0.325_009_94, 0.08515139, 0.679_840_9, 0.6975669,
-            0.736_234_25, 0.55053085, 0.692_227_66, 0.22727048, 0.13594262, 0.10002105, 0.16099514,
-            0.07719103, 0.23984282, 0.9083058, 0.642_227_4, 0.968_934_2, 0.74662715, 0.715_620_1,
-            0.736_546_5, 0.70610344, 0.221_011_18, 0.755_721_87, 0.691_958_84, 0.837_414, 0.27583158,
-            0.572_570_5, 0.681_606, 0.392_373_38, 0.33524343, 0.893_968_34, 0.602_969_35, 0.171_301_13,
-            0.1733834, 0.771_278_2, 0.99537134, 0.915_049_6, 0.493_121_1, 0.430_352_03, 0.70297265,
-            0.367_341_8, 0.4551964, 0.471_043_14, 0.603_747_8, 0.738_726_85, 0.5630592, 0.974_402_25,
-            0.633_682_85, 0.841_092_94, 0.24447136, 0.750384, 0.16893725, 0.542_256_65, 0.435_607_82,
+            0.32352856,
+            0.06571674,
+            0.01939427,
+            0.06352045,
+            0.708_527,
+            0.617_221_7,
+            0.16638431,
+            0.628_400_74,
+            0.554_893_9,
+            0.240_076_77,
+            0.325_009_94,
+            0.08515139,
+            0.679_840_9,
+            0.6975669,
+            0.736_234_25,
+            0.55053085,
+            0.692_227_66,
+            0.22727048,
+            0.13594262,
+            0.10002105,
+            0.16099514,
+            0.07719103,
+            0.23984282,
+            0.9083058,
+            0.642_227_4,
+            0.968_934_2,
+            0.74662715,
+            0.715_620_1,
+            0.736_546_5,
+            0.70610344,
+            0.221_011_18,
+            0.755_721_87,
+            0.691_958_84,
+            0.837_414,
+            0.27583158,
+            0.572_570_5,
+            0.681_606,
+            0.392_373_38,
+            0.33524343,
+            0.893_968_34,
+            0.602_969_35,
+            0.171_301_13,
+            0.1733834,
+            0.771_278_2,
+            0.99537134,
+            0.915_049_6,
+            0.493_121_1,
+            0.430_352_03,
+            0.70297265,
+            0.367_341_8,
+            0.4551964,
+            0.471_043_14,
+            0.603_747_8,
+            0.738_726_85,
+            0.5630592,
+            0.974_402_25,
+            0.633_682_85,
+            0.841_092_94,
+            0.24447136,
+            0.750384,
+            0.16893725,
+            0.542_256_65,
+            0.435_607_82,
             0.414_971_23,
         ];
         let (width, height) = (8, 8);
@@ -165,53 +210,215 @@ mod tests {
 
         blur.box_blur_h(&src, &mut dst, width, 1);
         let mut sol: Vec<f32> = vec![
-            0.339_215_37, 0.136_213_18, 0.04954382, 0.263_813_9, 0.46308973, 0.497_377_7, 0.470_668_94,
-            0.372_771_2, 0.448_500_5, 0.373_326_87, 0.21674603, 0.363_334_1, 0.48751974, 0.70454735,
-            0.661_444, 0.613_886_36, 0.609_268, 0.351_813_58, 0.15441138, 0.1323196, 0.11273574,
-            0.159343, 0.40844655, 0.613_458_75, 0.788_961_2, 0.785_929_56, 0.810_393_8, 0.732_931_26,
-            0.719_423_35, 0.554_553_7, 0.560_945_5, 0.539_653_5, 0.807_780_4, 0.601_734_8, 0.561_938_7,
-            0.510_002_7, 0.548_849_94, 0.46974093, 0.540_528_4, 0.640_390_2, 0.40154082, 0.315_884_62,
-            0.371_987_58, 0.646_677_6, 0.893_899_74, 0.801_180_66, 0.612_840_9, 0.508_814_16, 0.681_572_2,
-            0.508_503_6, 0.431_193_77, 0.509_995_76, 0.604_505_9, 0.635_177_9, 0.758_729_4, 0.746_811_33,
-            0.629_915_65, 0.573_082_4, 0.611_982_76, 0.38793087, 0.48719263, 0.38226724, 0.464_278_58,
+            0.339_215_37,
+            0.136_213_18,
+            0.04954382,
+            0.263_813_9,
+            0.46308973,
+            0.497_377_7,
+            0.470_668_94,
+            0.372_771_2,
+            0.448_500_5,
+            0.373_326_87,
+            0.21674603,
+            0.363_334_1,
+            0.48751974,
+            0.70454735,
+            0.661_444,
+            0.613_886_36,
+            0.609_268,
+            0.351_813_58,
+            0.15441138,
+            0.1323196,
+            0.11273574,
+            0.159343,
+            0.40844655,
+            0.613_458_75,
+            0.788_961_2,
+            0.785_929_56,
+            0.810_393_8,
+            0.732_931_26,
+            0.719_423_35,
+            0.554_553_7,
+            0.560_945_5,
+            0.539_653_5,
+            0.807_780_4,
+            0.601_734_8,
+            0.561_938_7,
+            0.510_002_7,
+            0.548_849_94,
+            0.46974093,
+            0.540_528_4,
+            0.640_390_2,
+            0.40154082,
+            0.315_884_62,
+            0.371_987_58,
+            0.646_677_6,
+            0.893_899_74,
+            0.801_180_66,
+            0.612_840_9,
+            0.508_814_16,
+            0.681_572_2,
+            0.508_503_6,
+            0.431_193_77,
+            0.509_995_76,
+            0.604_505_9,
+            0.635_177_9,
+            0.758_729_4,
+            0.746_811_33,
+            0.629_915_65,
+            0.573_082_4,
+            0.611_982_76,
+            0.38793087,
+            0.48719263,
+            0.38226724,
+            0.464_278_58,
             0.494_753_96,
         ];
         for (v1, v2) in dst.iter().zip(sol) {
-            assert!((v1 - v2).abs() < 1e-6);
+            assert!((v1 - v2).abs() < 1e-6, "box_blur_h failure");
         }
 
         blur.box_blur_v(&src, &mut dst, width, height, 1, 1.0);
         sol = vec![
-            0.504_035_1, 0.382_295_5, 0.19629186, 0.299_685_27, 0.519_101_74, 0.619_015_1, 0.446_075_47,
-            0.531_300_96, 0.523_550_03, 0.177688, 0.16011561, 0.08289763, 0.516_454_34, 0.46399322,
-            0.38082045, 0.695_745_8, 0.629_783_03, 0.47876048, 0.402_526_56, 0.300_264_18, 0.5257942,
-            0.49362046, 0.3990294, 0.738_186_2, 0.675_471_3, 0.677_872_9, 0.386_133_8, 0.46273723,
-            0.526_382_57, 0.391_889_3, 0.265_365_8, 0.852_665_36, 0.645_718_5, 0.659_216_46, 0.39861405,
-            0.686_489_6, 0.804_508, 0.671_175_5, 0.349_791_88, 0.693_347_4, 0.665_966_9, 0.458_685_64,
-            0.30147046, 0.604_963_96, 0.760_241_7, 0.682_05, 0.463_807_9, 0.766_240_9, 0.6465416,
-            0.459_911_97, 0.291_017_06, 0.664_235_1, 0.589_352_13, 0.732_011, 0.497_262_72, 0.606_575_13,
-            0.553_394_7, 0.42471716, 0.23968734, 0.428_315_88, 0.493_737_34, 0.632_735_1, 0.388_350_46,
+            0.504_035_1,
+            0.382_295_5,
+            0.19629186,
+            0.299_685_27,
+            0.519_101_74,
+            0.619_015_1,
+            0.446_075_47,
+            0.531_300_96,
+            0.523_550_03,
+            0.177688,
+            0.16011561,
+            0.08289763,
+            0.516_454_34,
+            0.46399322,
+            0.38082045,
+            0.695_745_8,
+            0.629_783_03,
+            0.47876048,
+            0.402_526_56,
+            0.300_264_18,
+            0.5257942,
+            0.49362046,
+            0.3990294,
+            0.738_186_2,
+            0.675_471_3,
+            0.677_872_9,
+            0.386_133_8,
+            0.46273723,
+            0.526_382_57,
+            0.391_889_3,
+            0.265_365_8,
+            0.852_665_36,
+            0.645_718_5,
+            0.659_216_46,
+            0.39861405,
+            0.686_489_6,
+            0.804_508,
+            0.671_175_5,
+            0.349_791_88,
+            0.693_347_4,
+            0.665_966_9,
+            0.458_685_64,
+            0.30147046,
+            0.604_963_96,
+            0.760_241_7,
+            0.682_05,
+            0.463_807_9,
+            0.766_240_9,
+            0.6465416,
+            0.459_911_97,
+            0.291_017_06,
+            0.664_235_1,
+            0.589_352_13,
+            0.732_011,
+            0.497_262_72,
+            0.606_575_13,
+            0.553_394_7,
+            0.42471716,
+            0.23968734,
+            0.428_315_88,
+            0.493_737_34,
+            0.632_735_1,
+            0.388_350_46,
             0.672_591_45,
         ];
         for (v1, v2) in dst.iter().zip(sol) {
-            assert!((v1 - v2).abs() < 1e-6);
+            assert!((v1 - v2).abs() < 1e-6, "box_blur_v failure");
         }
 
         blur.box_blur(&mut src, &mut dst, width, height, 1, 1.0);
         sol = vec![
-            0.472_543_84, 0.36087415, 0.29275754, 0.338_359_62, 0.47926736, 0.528_064_1, 0.5321305,
-            0.493_803_83, 0.465_661_3, 0.287_117_9, 0.140_233_76, 0.253_155_86, 0.3544484, 0.453_756,
-            0.513_519_8, 0.5333721, 0.615_576_57, 0.503_69, 0.393_850_42, 0.40952832, 0.43989295,
-            0.472_814_68, 0.543_612, 0.588_999_5, 0.735_336_54, 0.579826, 0.508_914_65, 0.458_417_86,
-            0.460_336_36, 0.39454588, 0.503_306_8, 0.597_834_17, 0.666_094_1, 0.567_849_7, 0.581_440_03,
-            0.62987053, 0.720_724_34, 0.608_491_8, 0.571_438_25, 0.562_952_64, 0.630_297_84, 0.475_374_34,
-            0.455_04, 0.5555587, 0.682_418_5, 0.635_366_5, 0.63736624, 0.632_005_2, 0.571_009_6,
-            0.465_823_53, 0.471_721_38, 0.5148681, 0.661_866_07, 0.606_208_6, 0.611_949_6, 0.583_459_8,
-            0.550_234_44, 0.405_933_05, 0.364_240_14, 0.38724685, 0.518_262_74, 0.504_940_9, 0.564_559,
+            0.472_543_84,
+            0.36087415,
+            0.29275754,
+            0.338_359_62,
+            0.47926736,
+            0.528_064_1,
+            0.5321305,
+            0.493_803_83,
+            0.465_661_3,
+            0.287_117_9,
+            0.140_233_76,
+            0.253_155_86,
+            0.3544484,
+            0.453_756,
+            0.513_519_8,
+            0.5333721,
+            0.615_576_57,
+            0.503_69,
+            0.393_850_42,
+            0.40952832,
+            0.43989295,
+            0.472_814_68,
+            0.543_612,
+            0.588_999_5,
+            0.735_336_54,
+            0.579826,
+            0.508_914_65,
+            0.458_417_86,
+            0.460_336_36,
+            0.39454588,
+            0.503_306_8,
+            0.597_834_17,
+            0.666_094_1,
+            0.567_849_7,
+            0.581_440_03,
+            0.62987053,
+            0.720_724_34,
+            0.608_491_8,
+            0.571_438_25,
+            0.562_952_64,
+            0.630_297_84,
+            0.475_374_34,
+            0.455_04,
+            0.5555587,
+            0.682_418_5,
+            0.635_366_5,
+            0.63736624,
+            0.632_005_2,
+            0.571_009_6,
+            0.465_823_53,
+            0.471_721_38,
+            0.5148681,
+            0.661_866_07,
+            0.606_208_6,
+            0.611_949_6,
+            0.583_459_8,
+            0.550_234_44,
+            0.405_933_05,
+            0.364_240_14,
+            0.38724685,
+            0.518_262_74,
+            0.504_940_9,
+            0.564_559,
             0.538_112_16,
         ];
         for (v1, v2) in src.iter().zip(sol) {
-            assert!((v1 - v2).abs() < 1e-6);
+            assert!((v1 - v2).abs() < 1e-6, "box_blur failure");
         }
     }
 
@@ -226,4 +433,129 @@ mod tests {
         let boxes = Blur::boxes_for_gaussian::<3>(2.5);
         assert_eq!(boxes, [2, 2, 2]);
     }
+
+    #[test]
+    fn total_blur_test() {
+        let height = 10;
+        let width = 10;
+        let mut src = (1..=(height * width))
+            .map(|i| (i as f32).recip())
+            .collect::<Vec<_>>();
+        let mut blur = Blur::new(width);
+
+        blur.run(
+            &mut src,
+            &mut vec![0.0; width * height],
+            width,
+            height,
+            2_f32,
+            0.1,
+        );
+
+        let sol = vec![
+            0.050528992,
+            0.044103604,
+            0.038919702,
+            0.032494307,
+            0.027310405,
+            0.020885015,
+            0.023104476,
+            0.020885015,
+            0.023104476,
+            0.028288381,
+            0.043704934,
+            0.038152207,
+            0.033674292,
+            0.028121557,
+            0.023643643,
+            0.018090911,
+            0.020009955,
+            0.018090911,
+            0.020009955,
+            0.024487872,
+            0.03968891,
+            0.03461781,
+            0.03053501,
+            0.025463907,
+            0.021381106,
+            0.016310005,
+            0.018066125,
+            0.016310005,
+            0.018066125,
+            0.022148928,
+            0.032864854,
+            0.028666414,
+            0.025289603,
+            0.021091158,
+            0.017714344,
+            0.013515903,
+            0.014971604,
+            0.013515901,
+            0.014971604,
+            0.01834842,
+            0.02884883,
+            0.025132021,
+            0.022150321,
+            0.018433508,
+            0.015451807,
+            0.011734996,
+            0.013027772,
+            0.011734993,
+            0.013027772,
+            0.016009476,
+            0.022024775,
+            0.019180624,
+            0.016904911,
+            0.014060758,
+            0.011785044,
+            0.008940893,
+            0.009933252,
+            0.00894089,
+            0.009933252,
+            0.012208968,
+            0.02513346,
+            0.021875666,
+            0.019268055,
+            0.016010256,
+            0.013402643,
+            0.010144847,
+            0.011281048,
+            0.010144845,
+            0.011281048,
+            0.013888664,
+            0.022024775,
+            0.019180622,
+            0.016904911,
+            0.014060758,
+            0.011785044,
+            0.008940893,
+            0.009933252,
+            0.00894089,
+            0.009933252,
+            0.012208967,
+            0.02513346,
+            0.021875666,
+            0.019268055,
+            0.016010256,
+            0.013402643,
+            0.010144847,
+            0.011281048,
+            0.010144845,
+            0.011281048,
+            0.013888664,
+            0.029149484,
+            0.02541006,
+            0.022407336,
+            0.018667907,
+            0.015665181,
+            0.011925754,
+            0.013224879,
+            0.011925753,
+            0.013224879,
+            0.016227607,
+        ];
+        for (v1, v2) in src.iter().zip(sol) {
+            assert!((v1 - v2).abs() < 1e-6, "run failure {} vs {}", v1, v2);
+        }
+    }
 }

src/buffer.rs

@@ -1,39 +1,26 @@
-#[derive(Debug)]
+#[derive(Debug, Clone)]
 pub struct Buf {
-    pub width: usize,
-    pub height: usize,
+    width: usize,
+    height: usize,
     pub buf: Vec<f32>,
 }
 
-impl Clone for Buf {
-    fn clone(&self) -> Buf {
-        Buf {
-            width: self.width,
-            height: self.height,
-            buf: self.buf.clone(),
-        }
-    }
-}
-
 impl Buf {
-    pub fn new(width: usize, height: usize, buf: Vec<f32>) -> Self {
+    pub fn new(width: usize, height: usize) -> Self {
         Buf {
             width,
             height,
-            buf,
+            buf: vec![0.0; height * width],
         }
     }
 
-    // Truncate x and y and return a corresponding index into the data slice.
-    fn index(&self, x: f32, y: f32) -> usize {
-        // x/y can come in negative, hence we shift them by width/height.
-        let i = (x + self.width as f32) as usize & (self.width - 1);
-        let j = (y + self.height as f32) as usize & (self.height - 1);
-        j * self.width + i
+    /// Truncate x and y and return a corresponding index into the data slice.
+    const fn index(&self, x: f32, y: f32) -> usize {
+        crate::util::index(self.width, self.height, x, y)
     }
 
-    // Get the buffer value at a given position. The implementation effectively treats data as periodic, hence any finite position will produce a value.
+    /// Get the buffer value at a given position. The implementation effectively treats data as periodic, hence any finite position will produce a value.
     pub fn get_buf(&self, x: f32, y: f32) -> f32 {
         self.buf[self.index(x, y)]
     }
-}
+}

src/grid.rs

@@ -1,87 +1,41 @@
-use crate::{
-    blur::Blur,
-    agent::Agent,
-    buffer::Buf,
-};
-
-use rand::{distributions::Uniform, Rng};
-use std::fmt::{Display, Formatter};
+use crate::{agent::Agent, blur::Blur, buffer::Buf};
+use rand::Rng;
+use rand_distr::Uniform;
 use rayon::{iter::ParallelIterator, prelude::*};
+use std::fmt::{Display, Formatter};
 
-// A population configuration.
-#[derive(Debug)]
+/// A population configuration.
+#[derive(Debug, Clone, Copy)]
 pub struct PopulationConfig {
     pub sensor_distance: f32,
     pub step_distance: f32,
     pub sensor_angle: f32,
     pub rotation_angle: f32,
 
-    decay_factor: f32,
     deposition_amount: f32,
 }
 
-impl Clone for PopulationConfig {
-    fn clone(&self) -> PopulationConfig {
-        PopulationConfig {
-            sensor_distance: self.sensor_distance,
-            step_distance: self.step_distance,
-            sensor_angle: self.sensor_angle,
-            rotation_angle: self.rotation_angle,
-            decay_factor: self.decay_factor,
-            deposition_amount: self.deposition_amount,
-        }
-    }
-}
-
 impl Display for PopulationConfig {
     fn fmt(&self, f: &mut Formatter<'_>) -> std::fmt::Result {
-        write!(
-            f,
-            "{{\nSensor Distance: {},\nStep Distance: {},\nSensor Angle: {},\nRotation Angle: {},\nDecay Factor: {},\nDeposition Amount: {},\n}}",
-            self.sensor_distance,
-            self.step_distance,
-            self.sensor_angle,
-            self.rotation_angle,
-            self.decay_factor,
-            self.deposition_amount
-        )
+        write!(f, "{:?}", self)
     }
 }
 
 impl PopulationConfig {
-    const SENSOR_ANGLE_MIN: f32 = 0.0;
-    const SENSOR_ANGLE_MAX: f32 = 120.0;
-    const SENSOR_DISTANCE_MIN: f32 = 0.0;
-    const SENSOR_DISTANCE_MAX: f32 = 64.0;
-    const ROTATION_ANGLE_MIN: f32 = 0.0;
-    const ROTATION_ANGLE_MAX: f32 = 120.0;
-    const STEP_DISTANCE_MIN: f32 = 0.2;
-    const STEP_DISTANCE_MAX: f32 = 2.0;
-    const DEPOSITION_AMOUNT_MIN: f32 = 5.0;
-    const DEPOSITION_AMOUNT_MAX: f32 = 5.0;
-    const DECAY_FACTOR_MIN: f32 = 0.1;
-    const DECAY_FACTOR_MAX: f32 = 0.1;
-
-    // Construct a random configuration.
+    /// Construct a random configuration.
     pub fn new<R: Rng + ?Sized>(rng: &mut R) -> Self {
         PopulationConfig {
-            sensor_distance: rng.gen_range(Self::SENSOR_DISTANCE_MIN..=Self::SENSOR_DISTANCE_MAX),
-            step_distance: rng.gen_range(Self::STEP_DISTANCE_MIN..=Self::STEP_DISTANCE_MAX),
-            decay_factor: rng.gen_range(Self::DECAY_FACTOR_MIN..=Self::DECAY_FACTOR_MAX),
-            sensor_angle: rng
-                .gen_range(Self::SENSOR_ANGLE_MIN..=Self::SENSOR_ANGLE_MAX)
-                .to_radians(),
-            rotation_angle: rng
-                .gen_range(Self::ROTATION_ANGLE_MIN..=Self::ROTATION_ANGLE_MAX)
-                .to_radians(),
-            deposition_amount: rng
-                .gen_range(Self::DEPOSITION_AMOUNT_MIN..=Self::DEPOSITION_AMOUNT_MAX),
+            sensor_distance: rng.random_range(0.0..=64.0),
+            step_distance: rng.random_range(0.2..=2.0),
+            sensor_angle: rng.random_range(0.0_f32..=120.0).to_radians(),
+            rotation_angle: rng.random_range(0.0_f32..=120.0).to_radians(),
+            deposition_amount: rng.random_range(5.0..=5.0),
         }
     }
 }
 
 // A 2D grid with a scalar value per each grid block. Each grid is occupied by a single population, hence we store the population config inside the grid.
-#[derive(Debug)]
+#[derive(Debug, Clone)]
 pub struct Grid {
     pub config: PopulationConfig,
     pub width: usize,
@@ -90,33 +44,21 @@ pub struct Grid {
     pub data: Vec<f32>,
 
     // Scratch space for the blur operation.
-    // pub buf: Vec<f32>,
-    pub buf: Buf,
-    pub blur: Blur,
-    pub agents: Vec<Agent>
-}
+    buf: Buf,
 
-impl Clone for Grid {
-    fn clone(&self) -> Grid {
-        Grid {
-            config: self.config.clone(),
-            width: self.width,
-            height: self.height,
-            data: self.data.clone(),
-            buf: self.buf.clone(),
-            blur: self.blur.clone(),
-            agents: self.agents.clone(),
-        }
-    }
+    blur: Blur,
+    pub agents: Vec<Agent>,
 }
 
 impl Grid {
-    // Create a new grid filled with random floats in the [0.0..1.0) range.
-    pub fn new<R: Rng + ?Sized>(width: usize, height: usize, rng: &mut R, agents: Vec<Agent>) -> Self {
-        if !width.is_power_of_two() || !height.is_power_of_two() {
-            panic!("Grid dimensions must be a power of two.");
-        }
-        let range = Uniform::from(0.0..1.0);
+    /// Create a new grid filled with random floats in the [0.0..1.0) range.
+    pub fn new<R: Rng + ?Sized>(
+        width: usize,
+        height: usize,
+        rng: &mut R,
+        agents: Vec<Agent>,
+    ) -> Self {
+        let range = Uniform::new(0.0, 1.0).expect("unable to create uniform distr");
         let data = rng.sample_iter(range).take(width * height).collect();
 
         Grid {
@@ -124,34 +66,18 @@ impl Grid {
             height,
             data,
             config: PopulationConfig::new(rng),
-            buf: Buf::new(width, height, vec![0.0; width * height]),
+            buf: Buf::new(width, height),
             blur: Blur::new(width),
             agents,
         }
     }
 
-    // Truncate x and y and return a corresponding index into the data slice.
-    fn index(&self, x: f32, y: f32) -> usize {
-        // x/y can come in negative, hence we shift them by width/height.
-        let i = (x + self.width as f32) as usize & (self.width - 1);
-        let j = (y + self.height as f32) as usize & (self.height - 1);
-        j * self.width + i
+    /// Truncate x and y and return a corresponding index into the data slice.
+    const fn index(&self, x: f32, y: f32) -> usize {
+        crate::util::index(self.width, self.height, x, y)
     }
 
-    /*
-    // Get the buffer value at a given position. The implementation effectively treats data as periodic, hence any finite position will produce a value.
-    pub fn get_buf(&self, x: f32, y: f32) -> f32 {
-        self.buf.buf[self.index(x, y)]
-    }
-    */
-
-    // Add a value to the grid data at a given position.
-    pub fn deposit(&mut self, x: f32, y: f32) {
-        let idx = self.index(x, y);
-        self.data[idx] += self.config.deposition_amount;
-    }
-
-    // Diffuse grid data and apply a decay multiplier.
+    /// Diffuse grid data and apply a decay multiplier.
     pub fn diffuse(&mut self, radius: usize) {
         self.blur.run(
             &mut self.data,
@@ -159,59 +85,23 @@ impl Grid {
             self.width,
             self.height,
             radius as f32,
-            self.config.decay_factor,
+            0.1, // decay is always 0.1
         );
     }
 
-    #[inline]
-    pub fn tick(&mut self, sin_table: Vec<f32>, cos_table: Vec<f32>) {
-        let (width, height) = (self.width, self.height);
-        let PopulationConfig {
-            sensor_distance,
-            sensor_angle,
-            rotation_angle,
-            step_distance,
-            ..
-        } = self.config;
-
-        let buf = self.buf.clone();
-
+    pub fn tick(&mut self) {
         self.agents.par_iter_mut().for_each(|agent| {
-            agent.tick(&buf, 
-                sensor_distance, sensor_angle, 
-                rotation_angle, step_distance,
-                width, height, &sin_table, &cos_table);
+            agent.tick(&self.buf, self.config, self.width, self.height);
         });
         self.deposit_all();
     }
 
-    #[inline]
     pub fn deposit_all(&mut self) {
-        let agent_list = self.agents.clone();
-        for agent in agent_list.iter() {
-            self.deposit(agent.x, agent.y);
+        for agent in self.agents.iter() {
+            let idx = self.index(agent.x, agent.y);
+            self.data[idx] += self.config.deposition_amount;
         }
     }
-
-    // No longer needed (moved to imgdata.rs)
-    /*
-    pub fn quantile(&self, fraction: f32) -> f32 {
-        let index = if (fraction - 1.0_f32).abs() < f32::EPSILON {
-            self.data.len() - 1
-        } else {
-            (self.data.len() as f32 * fraction) as usize
-        };
-        let mut sorted = self.data.clone();
-        sorted
-            .as_mut_slice()
-            .select_nth_unstable_by(index, |a, b| a.partial_cmp(b).unwrap());
-        sorted[index]
-    }
-
-    pub fn data(&self) -> &[f32] {
-        &self.data
-    }
-    */
 }
 
 pub fn combine<T>(grids: &mut [Grid], attraction_table: &[T])
@@ -222,14 +112,16 @@ where
     let bufs: Vec<_> = grids.iter().map(|grid| &grid.buf.buf).collect();
 
     // We mutate grid buffers and read grid data. We use unsafe because we need shared/unique borrows on different fields of the same Grid struct.
-    bufs.iter().enumerate().for_each(|(i, buf)| unsafe {
+    bufs.iter().enumerate().for_each(|(i, buf)| {
         let buf_ptr = *buf as *const Vec<f32> as *mut Vec<f32>;
-        buf_ptr.as_mut().unwrap().fill(0.0);
+        // SAFETY! we can take these are raw pointers because we are
+        // getting it from a `&mut [Grid]`
+        let buf_ptr_mut = unsafe { buf_ptr.as_mut().unwrap_unchecked() };
+
+        buf_ptr_mut.fill(0.0);
         datas.iter().enumerate().for_each(|(j, other)| {
             let multiplier = attraction_table[i].as_ref()[j];
-            buf_ptr
-                .as_mut()
-                .unwrap()
+            buf_ptr_mut
                 .iter_mut()
                 .zip(*other)
                 .for_each(|(to, from)| *to += from * multiplier)
@@ -241,16 +133,9 @@ where
 mod tests {
     use super::*;
 
-    #[test]
-    #[should_panic]
-    fn test_grid_new_panics() {
-        let mut rng = rand::thread_rng();
-        let _ = Grid::new(5, 5, &mut rng, vec![]);
-    }
-
     #[test]
     fn test_grid_new() {
-        let mut rng = rand::thread_rng();
+        let mut rng = rand::rng();
         let grid = Grid::new(8, 8, &mut rng, vec![]);
         assert_eq!(grid.index(0.5, 0.6), 0);
         assert_eq!(grid.index(1.5, 0.6), 1);
@@ -258,6 +143,5 @@ mod tests {
         assert_eq!(grid.index(2.5, 0.6), 2);
         assert_eq!(grid.index(2.5, 1.6), 10);
         assert_eq!(grid.index(7.9, 7.9), 63);
-        assert_eq!(grid.index(-0.5, -0.6), 0);
     }
 }

src/imgdata.rs

@@ -1,50 +1,30 @@
 use crate::{grid::Grid, palette::Palette};
-
+use image::RgbImage;
 use itertools::multizip;
 
-// Stores data that is located in grids that is used for image generation
+/// Stores data that is located in grids that is used for image generation
+#[derive(Clone)]
 pub struct ThinGridData {
-    pub width: usize,
-    pub height: usize,
-    pub data: Vec<f32>,
-}
-
-
-impl Clone for ThinGridData {
-    fn clone(&self) -> ThinGridData {
-        ThinGridData {
-            width: self.width,
-            height: self.height,
-            data: self.data.clone(),
-        }
-    }
+    width: usize,
+    height: usize,
+    data: Vec<f32>,
 }
 
 impl ThinGridData {
-    // Convert Grid to ThinGridData
-    pub fn new_from_grid(in_grid: &Grid) -> Self {
+    /// Convert Grid to ThinGridData
+    pub fn new_from_grid(in_grid: Grid) -> Self {
         ThinGridData {
             width: in_grid.width,
             height: in_grid.height,
-            data: in_grid.data.clone(),
+            data: in_grid.data,
         }
     }
 
-    #[allow(dead_code)]
-    pub fn new_from_grid_vec(in_grids: Vec<Grid>) -> Vec<Self> {
-        return in_grids.iter().map(|grid|{
-            Self::new_from_grid(grid)
-        }).collect();
+    pub fn new_from_grid_vec(in_grids: &[Grid]) -> Vec<Self> {
+        in_grids.iter().cloned().map(Self::new_from_grid).collect()
     }
 
-    // from grid.rs (needed in image gen)
-    #[allow(dead_code)]
-    pub fn data(&self) -> &[f32] {
-        &self.data
-    }
-
-    // from grid.rs (needed in image gen)
-    #[allow(dead_code)]
+    /// from grid.rs (needed in image gen)
     pub fn quantile(&self, fraction: f32) -> f32 {
         let index = if (fraction - 1.0_f32).abs() < f32::EPSILON {
             self.data.len() - 1
@@ -57,65 +37,33 @@ impl ThinGridData {
             .select_nth_unstable_by(index, |a, b| a.partial_cmp(b).unwrap());
         sorted[index]
     }
-
-    pub fn size_of(&self) -> usize {
-        let mut output: usize = 0;
-        output = output + std::mem::size_of_val(&self.width);
-        output = output + std::mem::size_of_val(&self.height);
-        for i in self.data.iter() {
-            output = output + std::mem::size_of_val(&i);
-        }
-        return output;
-    }
 }
 
-// Class for storing data that will be used to create images
+/// Class for storing data that will be used to create images
+#[derive(Clone)]
 pub struct ImgData {
-    pub grids: Vec<ThinGridData>,
-    pub palette: Palette,
-    pub iteration: i32,
-}
-
-impl Clone for ImgData {
-    fn clone(&self) -> ImgData {
-        ImgData {
-            grids: self.grids.clone(),
-            palette: self.palette,
-            iteration: self.iteration,
-        }
-    }
+    grids: Vec<ThinGridData>,
+    palette: Palette,
 }
 
 impl ImgData {
-    pub fn new(in_grids: Vec<ThinGridData>, in_palette: Palette, in_iteration: i32) -> Self {
+    pub const fn new(in_grids: Vec<ThinGridData>, in_palette: Palette) -> Self {
         ImgData {
             grids: in_grids,
             palette: in_palette,
-            iteration: in_iteration,
         }
     }
 
-    pub fn size_of(&self) -> usize {
-        let mut output: usize = 0;
-        output = output + std::mem::size_of_val(&self.iteration);
-        output = output + std::mem::size_of_val(&self.palette);
-        for grid in self.grids.iter() {
-            output = output + grid.size_of();
-        }
-        return output;
-    }
-
-    #[inline]
-    pub fn save_to_image(&self) {
+    pub fn to_image(&self) -> RgbImage {
         let (width, height) = (self.grids[0].width, self.grids[0].height);
         let mut img = image::RgbImage::new(width as u32, height as u32);
-    
+
         let max_values: Vec<_> = self
             .grids
             .iter()
             .map(|grid| grid.quantile(0.999) * 1.5)
             .collect();
-    
+
         for y in 0..height {
             for x in 0..width {
                 let i = y * width + x;
@@ -123,7 +71,7 @@ impl ImgData {
                 for (grid, max_value, color) in
                     multizip((&self.grids, &max_values, &self.palette.colors))
                 {
-                    let mut t = (grid.data()[i] / max_value).clamp(0.0, 1.0);
+                    let mut t = (grid.data[i] / max_value).clamp(0.0, 1.0);
                     t = t.powf(1.0 / 2.2); // gamma correction
                     r += color.0[0] as f32 * t;
                     g += color.0[1] as f32 * t;
@@ -135,8 +83,6 @@ impl ImgData {
                 img.put_pixel(x as u32, y as u32, image::Rgb([r as u8, g as u8, b as u8]));
             }
         }
-    
-        img.save(format!("./tmp/out_{}.png", self.iteration).as_str())
-            .unwrap();
+        img
     }
 }

src/lib.rs

@@ -1,8 +1,8 @@
+pub mod agent;
 mod blur;
-mod grid;
-mod imgdata; // for storing image data
+mod buffer;
+pub mod grid;
+pub mod imgdata; // for storing image data
 pub mod model;
 mod palette;
 mod util; // for math things
-mod agent;
-mod buffer;

src/main.rs

@@ -1,39 +1,57 @@
-use physarum::model;
+use physarum::{
+    imgdata::{ImgData, ThinGridData},
+    model,
+};
+use std::io::Write;
 
 fn main() {
-    // # of iterations to go through
-    let n_iterations = 128;
-    // let n_iterations = 2048;
-    // let n_iterations = 1 << 14;
-
-    // Size of grid and pictures
-    let (width, height) = (256, 256);
-    // let (width, height) = (512, 512);
-    // let (width, height) = (1024, 1024);
-
-    // # of agents
-    // let n_particles = 1 << 10;
-    // let n_particles = 1 << 16;
-    // let n_particles = 1 << 20;
-    let n_particles = 1 << 24;
-    println!("n_particles: {}", n_particles);
-
+    let n_iterations = 1024;
+    let (width, height) = (1024, 1024);
+    let n_particles = 1 << 22;
     let diffusivity = 1;
+    let n_populations = 3;
 
-    // `n_populations` is the # of types of agents
-    // let n_populations = 4;
-    let n_populations = 1;
-    // let n_populations = 1 + rng.gen_range(1..4); // make # of populations between 2 and 5
+    let mut model = model::Model::new(width, height, n_particles, n_populations, diffusivity);
+    model.print_configurations();
 
-    let mut model = model::Model::new(width, height, n_particles, n_populations, diffusivity); // Create the model
+    // Setup ffmpeg
+    let mut ffmpeg = std::process::Command::new("ffmpeg")
+        .args([
+            "-y",
+            "-f",
+            "rawvideo",
+            "-pix_fmt",
+            "rgb24",
+            "-s",
+            &format!("{}x{}", width, height),
+            "-r",
+            "30",
+            "-i",
+            "-",
+            "-c:v",
+            "libsvtav1",
+            "output.mp4",
+        ])
+        .stdin(std::process::Stdio::piped())
+        .spawn()
+        .expect("Failed to start ffmpeg");
+    let mut stdin = ffmpeg.stdin.take().unwrap();
 
-    model.print_configurations(); // Print config for model
+    for _ in 0..n_iterations {
+        model.step();
 
-    model.run(n_iterations); // Actually run the model
+        // Generate image
+        let grids = ThinGridData::new_from_grid_vec(model.population_grids());
+        let img_data = ImgData::new(grids, model.palette());
+        let img = img_data.to_image();
+        let raw_data = img.into_raw();
 
-    // export saved image data
-    println!("Rendering all saved image data....");
-    model.render_all_imgdata();
-    model.flush_image_data();
+        // Write to ffmpeg
+        stdin.write_all(&raw_data).unwrap();
+    }
+
+    // Cleanup
+    drop(stdin);
+    ffmpeg.wait().unwrap();
     println!("Done!");
 }

src/model.rs

@@ -1,36 +1,32 @@
 use crate::{
-    grid::{combine, Grid},
-    imgdata::{ImgData, ThinGridData},
-    palette::{random_palette, Palette},
     agent::Agent,
+    grid::{combine, Grid},
+    palette::{random_palette, Palette},
 };
-
-use indicatif::{ParallelProgressIterator, ProgressBar, ProgressStyle};
-// use rand::Rng;
+use indicatif::{ProgressBar, ProgressStyle};
 use rand_distr::{Distribution, Normal};
 use rayon::{iter::ParallelIterator, prelude::*};
-use std::{path::Path, time::Instant};
-use fastapprox::faster::{cos, sin};
+use std::time::Instant;
 
-// Top-level simulation class.
+/// Top-level simulation class.
 pub struct Model {
-    // The grid they move on.
-    grids: Vec<Grid>,
+    /// per-population grid (one for each population)
+    population_grids: Vec<Grid>,
 
-    // Attraction table governs interaction across populations
+    /// Attraction table governs interaction across populations
     attraction_table: Vec<Vec<f32>>,
 
-    // Global grid diffusivity.
+    /// Global grid diffusivity.
     diffusivity: usize,
 
-    // Current model iteration.
-    iteration: i32,
+    /// Current model iteration.
+    iteration: usize,
 
-    // Color palette
+    /// Color palette
     palette: Palette,
 
-    // List of ImgData to be processed post-simulation into images
-    img_data_vec: Vec<ImgData>,
+    time_per_agent_list: Vec<f64>,
+    time_per_step_list: Vec<f64>,
 }
 
 impl Model {
@@ -40,13 +36,13 @@ impl Model {
     const REPULSION_FACTOR_STD: f32 = 0.1;
 
     pub fn print_configurations(&self) {
-        for (i, grid) in self.grids.iter().enumerate() {
+        for (i, grid) in self.population_grids.iter().enumerate() {
             println!("Grid {}: {}", i, grid.config);
         }
         println!("Attraction table: {:#?}", self.attraction_table);
     }
 
-    // Construct a new model with random initial conditions and random configuration.
+    /// Construct a new model with random initial conditions and random configuration.
     pub fn new(
         width: usize,
         height: usize,
@@ -57,222 +53,92 @@ impl Model {
         let particles_per_grid = (n_particles as f64 / n_populations as f64).ceil() as usize;
         let _n_particles = particles_per_grid * n_populations;
 
-        let mut rng = rand::thread_rng();
+        let mut rng = rand::rng();
 
         let attraction_distr =
             Normal::new(Self::ATTRACTION_FACTOR_MEAN, Self::ATTRACTION_FACTOR_STD).unwrap();
-        let repulstion_distr =
+        let repulsion_distr =
             Normal::new(Self::REPULSION_FACTOR_MEAN, Self::REPULSION_FACTOR_STD).unwrap();
 
         let mut attraction_table = Vec::with_capacity(n_populations);
         for i in 0..n_populations {
             attraction_table.push(Vec::with_capacity(n_populations));
             for j in 0..n_populations {
-                attraction_table[i].push(if i == j {
-                    attraction_distr.sample(&mut rng)
-                } else {
-                    repulstion_distr.sample(&mut rng)
-                });
+                attraction_table[i].push(
+                    if i == j {
+                        &attraction_distr
+                    } else {
+                        &repulsion_distr
+                    }
+                    .sample(&mut rng),
+                );
             }
         }
 
-        let mut grids: Vec<Grid> = Vec::new();
-        for pop in 0..n_populations {
-            let agents = (0..particles_per_grid)
-                .map(|i| {
-                    Agent::new(width, height, pop, &mut rng, i)
-                }).collect();
-            grids.push(Grid::new(width, height, &mut rng, agents));
-        }
+        let grids = (0..n_populations)
+            .map(|_| {
+                let agents = (0..particles_per_grid)
+                    .map(|_| Agent::new(width, height, &mut rng))
+                    .collect();
+                Grid::new(width, height, &mut rng, agents)
+            })
+            .collect();
 
         Model {
-            grids,
+            population_grids: grids,
             attraction_table,
             diffusivity,
             iteration: 0,
             palette: random_palette(),
-            img_data_vec: Vec::new(),
+            time_per_agent_list: Vec::new(),
+            time_per_step_list: Vec::new(),
         }
     }
 
-    pub fn gen_table(opt: i32) -> Vec<f32> {
-        let mut output: Vec<f32> = Vec::new();
+    pub fn step(&mut self) {
+        combine(&mut self.population_grids, &self.attraction_table);
 
-        let max: f32 = 360.0;
-        let min: f32 = -360.0;
-        let interval: f32 = 0.01;
-        
+        let agents_tick_time = Instant::now();
+        self.population_grids.par_iter_mut().for_each(|grid| {
+            grid.tick();
+            grid.diffuse(self.diffusivity);
+        });
+        let agents_tick_elapsed = agents_tick_time.elapsed().as_millis() as f64;
+        let agents_num: usize = self.population_grids.iter().map(|g| g.agents.len()).sum();
+        let ms_per_agent = agents_tick_elapsed / agents_num as f64;
 
-        let mut i: f32;
-        let mut tmp: f32 = 0.0;
-        for i1 in ((min/interval)as i32)..((max/interval)as i32) {
-            i = (i1 as f32)*interval;
-            if opt == 0 {
-                tmp = sin(i);
-            } else if opt == 1 {
-                tmp = cos(i);
-            }
-            output.push(tmp);
-
-        }
-        println!("{}", output.len());
-
-        return output;
+        self.time_per_agent_list.push(ms_per_agent);
+        self.time_per_step_list.push(agents_tick_elapsed);
+        self.iteration += 1;
     }
 
-
-    // Simulates `steps` # of steps
-    #[inline]
     pub fn run(&mut self, steps: usize) {
-        let debug: bool = false;
-
         let pb = ProgressBar::new(steps as u64);
-        pb.set_style(
-            ProgressStyle::default_bar()
-                .template(
-                    "{spinner:.green} [{elapsed_precise}] [{bar:40.cyan/blue}] {pos}/{len} ({eta} {percent}%, {per_sec})",
-                )
-                .progress_chars("#>-"),
-        );
+        pb.set_style(ProgressStyle::default_bar()
+            .template("{spinner:.green} [{elapsed_precise}] [{bar:40.cyan/blue}] {pos}/{len} ({eta} {percent}%, {per_sec})").expect("invalid progresstyle template")
+            .progress_chars("#>-"));
 
-        let mut time_per_agent_list: Vec<f64> = Vec::new();
-        let mut time_per_step_list: Vec<f64> = Vec::new();
-
-        let agents_num: usize = self.grids.iter().map(|grid| grid.agents.len()).sum();
-        let sin_table1 = Self::gen_table(0);
-        let cos_table1 = Self::gen_table(1);
-        let sin_table = sin_table1.clone();
-        let cos_table = cos_table1.clone();
-        for i in 0..steps {
-            if debug {
-                println!("Starting tick for all agents...")
-            };
-
-            // Combine grids
-            let grids = &mut self.grids;
-            combine(grids, &self.attraction_table);
-            let agents_tick_time = Instant::now();
-
-            // Tick agents
-            let diffusivity = self.diffusivity;            
-            self.grids.par_iter_mut().for_each(|grid| {
-                grid.tick(sin_table.clone(), cos_table.clone());
-                grid.diffuse(diffusivity); // Diffuse + Decay
-            });
-
-            self.save_image_data();
-
-            let agents_tick_elapsed: f64 = agents_tick_time.elapsed().as_millis() as f64;
-            let ms_per_agent: f64 = (agents_tick_elapsed as f64) / (agents_num as f64);
-            time_per_agent_list.push(ms_per_agent);
-            time_per_step_list.push(agents_tick_elapsed);
-
-            if debug {
-                println!(
-                    "Finished tick for all agents. took {}ms\nTime per agent: {}ms\n",
-                    agents_tick_elapsed, ms_per_agent
-                )
-            };
-
-            self.iteration += 1;
-            pb.set_position(i as u64);
+        for _ in 0..steps {
+            self.step();
+            pb.inc(1);
         }
         pb.finish();
 
         let avg_per_step: f64 =
-            time_per_step_list.iter().sum::<f64>() as f64 / time_per_step_list.len() as f64;
+            self.time_per_step_list.iter().sum::<f64>() / self.time_per_step_list.len() as f64;
         let avg_per_agent: f64 =
-            time_per_agent_list.iter().sum::<f64>() as f64 / time_per_agent_list.len() as f64;
+            self.time_per_agent_list.iter().sum::<f64>() / self.time_per_agent_list.len() as f64;
         println!(
             "Average time per step: {}ms\nAverage time per agent: {}ms",
             avg_per_step, avg_per_agent
         );
     }
 
-    fn size_of_imgdata_vec(&self) -> usize {
-        return (self.img_data_vec[0].size_of() as usize) * (self.img_data_vec.len() as usize);
+    pub fn population_grids(&self) -> &[Grid] {
+        &self.population_grids
     }
 
-    fn save_image_data(&mut self) {
-        let grids = ThinGridData::new_from_grid_vec(self.grids.clone());
-        let img_data = ImgData::new(grids, self.palette, self.iteration);
-        self.img_data_vec.push(img_data);
-        let size: usize = self.size_of_imgdata_vec();
-        let mb = size/1024/1024;
-        // println!("{} B | {} KB | {} MB", size, size/1024, size/1024/1024);
-
-        let max_mb = 6000;
-        if mb >= max_mb {
-            println!("ram usage is over {} MB (and len of {}), flushing to disk\n", max_mb, self.img_data_vec.len());
-            self.render_all_imgdata();
-            self.flush_image_data();
-        }
+    pub fn palette(&self) -> Palette {
+        self.palette
     }
-
-    pub fn flush_image_data(&mut self) {
-        self.img_data_vec.clear();
-    }
-
-    pub fn render_all_imgdata(&self) {
-        if !Path::new("./tmp").exists() {
-            std::fs::create_dir("./tmp").expect("could create directory");
-        }
-
-        let pb = ProgressBar::new(self.img_data_vec.len() as u64);
-        pb.set_style(ProgressStyle::default_bar().template(
-            "{spinner:.green} [{elapsed_precise}] [{bar:40.cyan/blue}] ({pos}/{len}, {percent}%, {per_sec})",
-        ));
-
-        /*
-        for img in &self.img_data_vec {
-            Self::save_to_image(img.to_owned());
-            pb.inc(1);
-        }
-        pb.finish();
-        */
-
-        (&self.img_data_vec)
-            .par_iter()
-            .progress_with(pb)
-            .for_each(|img| {
-                // Self::save_to_image(img.to_owned());
-                img.save_to_image();
-            });
-    }
-
-    /*
-    pub fn save_to_image(imgdata: ImgData) {
-        let (width, height) = (imgdata.grids[0].width, imgdata.grids[0].height);
-        let mut img = image::RgbImage::new(width as u32, height as u32);
-
-        let max_values: Vec<_> = imgdata
-            .grids
-            .iter()
-            .map(|grid| grid.quantile(0.999) * 1.5)
-            .collect();
-
-        for y in 0..height {
-            for x in 0..width {
-                let i = y * width + x;
-                let (mut r, mut g, mut b) = (0.0_f32, 0.0_f32, 0.0_f32);
-                for (grid, max_value, color) in
-                    multizip((&imgdata.grids, &max_values, &imgdata.palette.colors))
-                {
-                    let mut t = (grid.data()[i] / max_value).clamp(0.0, 1.0);
-                    t = t.powf(1.0 / 2.2); // gamma correction
-                    r += color.0[0] as f32 * t;
-                    g += color.0[1] as f32 * t;
-                    b += color.0[2] as f32 * t;
-                }
-                r = r.clamp(0.0, 255.0);
-                g = g.clamp(0.0, 255.0);
-                b = b.clamp(0.0, 255.0);
-                img.put_pixel(x as u32, y as u32, image::Rgb([r as u8, g as u8, b as u8]));
-            }
-        }
-
-        img.save(format!("./tmp/out_{}.png", imgdata.iteration).as_str())
-            .unwrap();
-    }
-    */
 }

src/palette.rs

@@ -1,4 +1,4 @@
-use rand::{seq::SliceRandom, thread_rng, Rng};
+use rand::{seq::SliceRandom, Rng};
 
 #[derive(Clone, Copy)]
 pub struct Palette {
@@ -6,8 +6,8 @@ pub struct Palette {
 }
 
 pub fn random_palette() -> Palette {
-    let mut rng = thread_rng();
-    let mut palette = PALETTES[rng.gen_range(0..PALETTES.len())];
+    let mut rng = rand::rng();
+    let mut palette = PALETTES[rng.random_range(0..PALETTES.len())];
     palette.colors.shuffle(&mut rng);
     palette
 }
@@ -15,7 +15,7 @@ pub fn random_palette() -> Palette {
 const fn hex_to_color(c: usize) -> image::Rgb<u8> {
     let r = (c >> 16) & 0xff;
     let g = (c >> 8) & 0xff;
-    let b = (c >> 0) & 0xff;
+    let b = c & 0xff;
     image::Rgb::<u8>([r as u8, g as u8, b as u8])
 }
 

src/util.rs

@@ -1,4 +1,111 @@
 #[inline]
 pub fn wrap(x: f32, max: f32) -> f32 {
-    x - max * ((x > max) as i32 as f32 - (x < 0.0_f32) as i32 as f32)
+    // x - max * ((x > max) as i32 - x.is_sign_negative() as i32) as f32
+    x.rem_euclid(max)
+}
+
+/// Truncate x and y and return a corresponding index into the data slice.
+#[inline]
+pub const fn index(width: usize, height: usize, x: f32, y: f32) -> usize {
+    // x/y can come in negative, hence we shift them by width/height.
+    let i = (x + width as f32) as usize & (width - 1);
+    let j = (y + height as f32) as usize & (height - 1);
+    j * width + i
+}
+
+#[cfg(test)]
+mod test {
+    use super::*;
+
+    mod wrap {
+        use super::*;
+        #[test]
+        fn over() {
+            assert_eq!(wrap(1.1, 1.0), 0.100000024); // floating point weirdness
+        }
+
+        #[test]
+        fn middle() {
+            assert_eq!(wrap(0.5, 1.0), 0.5);
+        }
+
+        #[test]
+        fn under() {
+            assert_eq!(wrap(-1.0, 2.0), 1.0);
+        }
+    }
+
+    mod index {
+        use super::*;
+
+        #[test]
+        fn basic_positive_coordinates() {
+            let width = 4;
+            let height = 4;
+            assert_eq!(index(width, height, 1.5, 2.5), 9);
+        }
+
+        #[test]
+        fn negative_x_coordinate() {
+            let width = 8;
+            let height = 8;
+            assert_eq!(index(width, height, -3.2, 5.6), 44);
+        }
+
+        #[test]
+        fn exact_boundary_values() {
+            let width = 16;
+            let height = 16;
+            assert_eq!(index(width, height, 16.0, 0.0), 0);
+        }
+
+        #[test]
+        fn large_coordinates() {
+            let width = 2;
+            let height = 2;
+            assert_eq!(index(width, height, 1000.0, 2000.0), 0);
+        }
+
+        #[test]
+        fn negative_x_and_y() {
+            let width = 4;
+            let height = 4;
+            assert_eq!(index(width, height, -1.5, -0.5), 14);
+        }
+
+        #[test]
+        fn fractional_truncation() {
+            let width = 4;
+            let height = 4;
+            assert_eq!(index(width, height, 3.9, 3.999), 15);
+        }
+
+        #[test]
+        fn zero_coordinates() {
+            let width = 4;
+            let height = 4;
+            assert_eq!(index(width, height, 0.0, 0.0), 0);
+        }
+
+        #[test]
+        fn x_equals_width() {
+            let width = 8;
+            let height = 8;
+            assert_eq!(index(width, height, 8.0, 0.0), 0);
+        }
+
+        #[test]
+        fn y_negative_beyond_height() {
+            let width = 4;
+            let height = 4;
+            assert_eq!(index(width, height, 0.0, -4.5), 0);
+        }
+
+        #[test]
+        fn width_and_height_one() {
+            let width = 1;
+            let height = 1;
+            assert_eq!(index(width, height, 123.4, -56.7), 0);
+        }
+    }
 }