use av1 for encoding

.gitignore

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

Cargo.toml

@@ -5,14 +5,21 @@ authors = ["Simon Gardling <titaniumtown@gmail.com>", "mindv0rtex <mindv0rtex@us
 edition = "2021"
 
 [dependencies]
-image = "0.23"
+image = "0.25"
-indicatif = { version = "0.15", features = [ "rayon" ] }
+indicatif = { version = "0.17", features = [ "rayon" ] }
-itertools = "0.10"
+itertools = "0.14"
-rand = "0.8"
+rand = "0.9"
-rand_distr = "0.4"
+rand_distr = "0.5"
 rayon = "1.10"
 fastapprox = "0.3"
 
+[dev-dependencies]
+criterion = { version = "0.5", features = ["html_reports"] }
+
+[[bench]]
+name = "benchmark"
+harness = false
+
 [profile.release]
 codegen-units = 1
 opt-level = 3

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);

src/agent.rs

@@ -1,18 +1,17 @@
+use crate::grid::PopulationConfig;
 use crate::{buffer::Buf, util::wrap};
-
 use fastapprox::faster::{cos, sin};
-use rand::{seq::SliceRandom, Rng};
+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.
+/// 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,
+    heading: f32,
-    pub population_id: usize,
-    pub i: usize,
 }
 
 impl Display for Agent {
@@ -22,46 +21,28 @@ impl Display for Agent {
 }
 
 impl Agent {
-    // Construct a new agent with random parameters.
+    /// Construct a new agent with random parameters.
-    pub fn new<R: Rng + ?Sized>(
+    pub fn new<R: Rng + ?Sized>(width: usize, height: usize, rng: &mut R) -> Self {
-        width: usize,
+        let (x, y, angle) = rng.random::<(f32, f32, f32)>();
-        height: usize,
-        id: usize,
-        rng: &mut R,
-        i: usize,
-    ) -> Self {
-        let (x, y, angle) = rng.gen::<(f32, f32, f32)>();
         Agent {
             x: x * width as f32,
             y: y * height as f32,
-            angle: angle * TAU,
+            heading: angle * TAU,
-            population_id: id,
-            i,
         }
     }
 
-    // Tick an agent
+    /// Tick an agent
-    #[inline]
+    pub fn tick(&mut self, buf: &Buf, pop_config: PopulationConfig, width: usize, height: usize) {
-    pub fn tick(
+        let xc = self.x + cos(self.heading) * pop_config.sensor_distance;
-        &mut self,
+        let yc = self.y + sin(self.heading) * pop_config.sensor_distance;
-        buf: &Buf,
-        sensor_distance: f32,
-        sensor_angle: f32,
-        rotation_angle: f32,
-        step_distance: f32,
-        width: usize,
-        height: usize,
-    ) {
-        let xc = self.x + cos(self.angle) * sensor_distance;
-        let yc = self.y + sin(self.angle) * sensor_distance;
 
-        let agent_add_sens = self.angle + sensor_angle;
+        let agent_add_sens = self.heading + pop_config.sensor_angle;
-        let agent_sub_sens = self.angle - sensor_angle;
+        let agent_sub_sens = self.heading - pop_config.sensor_angle;
 
-        let xl = self.x + cos(agent_sub_sens) * sensor_distance;
+        let xl = self.x + cos(agent_sub_sens) * pop_config.sensor_distance;
-        let yl = self.y + sin(agent_sub_sens) * sensor_distance;
+        let yl = self.y + sin(agent_sub_sens) * pop_config.sensor_distance;
-        let xr = self.x + cos(agent_add_sens) * sensor_distance;
+        let xr = self.x + cos(agent_add_sens) * pop_config.sensor_distance;
-        let yr = self.y + sin(agent_add_sens) * 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);
@@ -69,23 +50,30 @@ impl Agent {
         let right = buf.get_buf(xr, yr);
 
         // Rotate and move logic
-        let mut rng = rand::thread_rng();
+        let direction = if (center > left) && (center > right) {
-        let mut direction: f32 = 0.0;
+            0.0
-
-        if (center > left) && (center > right) {
-            direction = 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), width as f32);
+        self.x = wrap(
-        self.y = wrap(self.y + step_distance * sin(self.angle), height as f32);
+            self.x + pop_config.step_distance * cos(self.heading),
+            width as f32,
+        );
+        self.y = wrap(
+            self.y + pop_config.step_distance * sin(self.heading),
+            height as f32,
+        );
     }
 }

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.
+    /// 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();
+        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;
+        let m = (0.25 * (N * (w + 3)) as f32 - 3.0 * sigma_sq / (w + 1) as f32).round() as usize;
-        m -= 3.0 * sigma * sigma / (w + 1) as f32;
-        let m = m.round() as usize;
 
-        let mut result = [0; N];
+        (0..N)
-        for (i, value) in result.iter_mut().enumerate() {
+            .map(|i| (w + 1 - 2 * (i < m) as usize) / 2)
-            *value = (if i < m { w - 1 } else { w + 1 }) / 2;
+            .collect::<Vec<_>>()
-        }
+            .try_into()
-        result
+            .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],
@@ -285,7 +276,7 @@ mod tests {
             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);
@@ -356,7 +347,7 @@ mod tests {
             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);
@@ -427,7 +418,7 @@ mod tests {
             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");
         }
     }
 
@@ -442,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,24 +1,25 @@
 #[derive(Debug, Clone)]
 pub struct Buf {
-    pub width: usize,
+    width: usize,
-    pub height: usize,
+    height: usize,
     pub buf: Vec<f32>,
 }
 
 impl Buf {
-    pub const fn new(width: usize, height: usize, buf: Vec<f32>) -> Self {
+    pub fn new(width: usize, height: usize) -> Self {
-        Buf { width, height, buf }
+        Buf {
+            width,
+            height,
+            buf: vec![0.0; height * width],
+        }
     }
 
-    // Truncate x and y and return a corresponding index into the data slice.
+    /// Truncate x and y and return a corresponding index into the data slice.
     const fn index(&self, x: f32, y: f32) -> usize {
-        // x/y can come in negative, hence we shift them by width/height.
+        crate::util::index(self.width, self.height, x, y)
-        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
     }
 
-    // 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,10 +1,10 @@
 use crate::{agent::Agent, blur::Blur, buffer::Buf};
-
+use rand::Rng;
-use rand::{distributions::Uniform, Rng};
+use rand_distr::Uniform;
 use rayon::{iter::ParallelIterator, prelude::*};
 use std::fmt::{Display, Formatter};
 
-// A population configuration.
+/// A population configuration.
 #[derive(Debug, Clone, Copy)]
 pub struct PopulationConfig {
     pub sensor_distance: f32,
@@ -12,7 +12,6 @@ pub struct PopulationConfig {
     pub sensor_angle: f32,
     pub rotation_angle: f32,
 
-    decay_factor: f32,
     deposition_amount: f32,
 }
 
@@ -23,33 +22,14 @@ impl Display for PopulationConfig {
 }
 
 impl PopulationConfig {
-    const SENSOR_ANGLE_MIN: f32 = 0.0;
+    /// Construct a random configuration.
-    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.
     pub fn new<R: Rng + ?Sized>(rng: &mut R) -> Self {
         PopulationConfig {
-            sensor_distance: rng.gen_range(Self::SENSOR_DISTANCE_MIN..=Self::SENSOR_DISTANCE_MAX),
+            sensor_distance: rng.random_range(0.0..=64.0),
-            step_distance: rng.gen_range(Self::STEP_DISTANCE_MIN..=Self::STEP_DISTANCE_MAX),
+            step_distance: rng.random_range(0.2..=2.0),
-            decay_factor: rng.gen_range(Self::DECAY_FACTOR_MIN..=Self::DECAY_FACTOR_MAX),
+            sensor_angle: rng.random_range(0.0_f32..=120.0).to_radians(),
-            sensor_angle: rng
+            rotation_angle: rng.random_range(0.0_f32..=120.0).to_radians(),
-                .gen_range(Self::SENSOR_ANGLE_MIN..=Self::SENSOR_ANGLE_MAX)
+            deposition_amount: rng.random_range(5.0..=5.0),
-                .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),
         }
     }
 }
@@ -64,24 +44,21 @@ pub struct Grid {
     pub data: Vec<f32>,
 
     // Scratch space for the blur operation.
-    // pub buf: Vec<f32>,
+    buf: Buf,
-    pub buf: Buf,
+
-    pub blur: Blur,
+    blur: Blur,
     pub agents: Vec<Agent>,
 }
 
 impl Grid {
-    // Create a new grid filled with random floats in the [0.0..1.0) range.
+    /// 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() {
+        let range = Uniform::new(0.0, 1.0).expect("unable to create uniform distr");
-            panic!("Grid dimensions must be a power of two.");
-        }
-        let range = Uniform::from(0.0..1.0);
         let data = rng.sample_iter(range).take(width * height).collect();
 
         Grid {
@@ -89,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.
+    /// Truncate x and y and return a corresponding index into the data slice.
-    fn index(&self, x: f32, y: f32) -> usize {
+    const fn index(&self, x: f32, y: f32) -> usize {
-        // x/y can come in negative, hence we shift them by width/height.
+        crate::util::index(self.width, self.height, x, y)
-        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
     }
 
-    /*
+    /// Diffuse grid data and apply a decay multiplier.
-    // 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.
     pub fn diffuse(&mut self, radius: usize) {
         self.blur.run(
             &mut self.data,
@@ -124,64 +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) {
-        let (width, height) = (self.width, self.height);
-        let PopulationConfig {
-            sensor_distance,
-            sensor_angle,
-            rotation_angle,
-            step_distance,
-            ..
-        } = self.config;
-
-        let buf = self.buf.clone();
-
         self.agents.par_iter_mut().for_each(|agent| {
-            agent.tick(
+            agent.tick(&self.buf, self.config, self.width, self.height);
-                &buf,
-                sensor_distance,
-                sensor_angle,
-                rotation_angle,
-                step_distance,
-                width,
-                height,
-            );
         });
         self.deposit_all();
     }
 
-    #[inline]
     pub fn deposit_all(&mut self) {
-        let agent_list = self.agents.clone();
+        for agent in self.agents.iter() {
-        for agent in agent_list.iter() {
+            let idx = self.index(agent.x, agent.y);
-            self.deposit(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])
@@ -192,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
+            buf_ptr_mut
-                .as_mut()
-                .unwrap()
                 .iter_mut()
                 .zip(*other)
                 .for_each(|(to, from)| *to += from * multiplier)
@@ -211,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);

src/imgdata.rs

@@ -1,34 +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
 #[derive(Clone)]
 pub struct ThinGridData {
-    pub width: usize,
+    width: usize,
-    pub height: usize,
+    height: usize,
-    pub data: Vec<f32>,
+    data: Vec<f32>,
 }
 
 impl ThinGridData {
-    // Convert Grid to ThinGridData
+    /// Convert Grid to ThinGridData
-    pub fn new_from_grid(in_grid: &Grid) -> Self {
+    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,
         }
     }
 
     pub fn new_from_grid_vec(in_grids: &[Grid]) -> Vec<Self> {
-        in_grids
+        in_grids.iter().cloned().map(Self::new_from_grid).collect()
-            .iter()
-            .map(|grid| Self::new_from_grid(grid))
-            .collect()
     }
 
-    // from grid.rs (needed in image gen)
+    /// 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
@@ -46,8 +42,8 @@ impl ThinGridData {
 /// Class for storing data that will be used to create images
 #[derive(Clone)]
 pub struct ImgData {
-    pub grids: Vec<ThinGridData>,
+    grids: Vec<ThinGridData>,
-    pub palette: Palette,
+    palette: Palette,
 }
 
 impl ImgData {

src/lib.rs

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

src/main.rs

@@ -1,38 +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 = 1024;
-    // 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);
-
+    let n_particles = 1 << 22;
-    // # 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 diffusivity = 1;
+    let n_populations = 3;
 
-    // `n_populations` is the # of types of agents
+    let mut model = model::Model::new(width, height, n_particles, n_populations, diffusivity);
-    // let n_populations = 4;
+    model.print_configurations();
-    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); // 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
+        // Write to ffmpeg
-    println!("Rendering all saved image data....");
+        stdin.write_all(&raw_data).unwrap();
-    model.render_all_imgdata();
+    }
+
+    // Cleanup
+    drop(stdin);
+    ffmpeg.wait().unwrap();
     println!("Done!");
 }

src/model.rs

@@ -1,35 +1,32 @@
 use crate::{
     agent::Agent,
     grid::{combine, Grid},
-    imgdata::{ImgData, ThinGridData},
     palette::{random_palette, Palette},
 };
-
+use indicatif::{ProgressBar, ProgressStyle};
-use indicatif::{ParallelProgressIterator, ProgressBar, ProgressStyle};
-// use rand::Rng;
 use rand_distr::{Distribution, Normal};
 use rayon::{iter::ParallelIterator, prelude::*};
-use std::{path::Path, time::Instant};
+use std::time::Instant;
 
-// Top-level simulation class.
+/// Top-level simulation class.
 pub struct Model {
-    // per-population grid (one for each population)
+    /// 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.
+    /// Current model iteration.
     iteration: usize,
 
-    // Color palette
+    /// Color palette
     palette: Palette,
 
-    // List of ImgData to be processed post-simulation into images
+    time_per_agent_list: Vec<f64>,
-    img_data_vec: Vec<(usize, ImgData)>,
+    time_per_step_list: Vec<f64>,
 }
 
 impl Model {
@@ -45,7 +42,7 @@ impl Model {
         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,
@@ -56,32 +53,36 @@ 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_table[i].push(
-                    attraction_distr.sample(&mut rng)
+                    if i == j {
-                } else {
+                        &attraction_distr
-                    repulstion_distr.sample(&mut rng)
+                    } else {
-                });
+                        &repulsion_distr
+                    }
+                    .sample(&mut rng),
+                );
             }
         }
 
-        let mut grids: Vec<Grid> = Vec::new();
+        let grids = (0..n_populations)
-        for pop in 0..n_populations {
+            .map(|_| {
-            let agents = (0..particles_per_grid)
+                let agents = (0..particles_per_grid)
-                .map(|i| Agent::new(width, height, pop, &mut rng, i))
+                    .map(|_| Agent::new(width, height, &mut rng))
-                .collect();
+                    .collect();
-            grids.push(Grid::new(width, height, &mut rng, agents));
+                Grid::new(width, height, &mut rng, agents)
-        }
+            })
+            .collect();
 
         Model {
             population_grids: grids,
@@ -89,102 +90,55 @@ impl Model {
             diffusivity,
             iteration: 0,
             palette: random_palette(),
-            img_data_vec: Vec::new(),
+            time_per_agent_list: Vec::new(),
+            time_per_step_list: Vec::new(),
         }
     }
 
-    // Simulates `steps` # of steps
+    pub fn step(&mut self) {
-    #[inline]
+        combine(&mut self.population_grids, &self.attraction_table);
+
+        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;
+
+        self.time_per_agent_list.push(ms_per_agent);
+        self.time_per_step_list.push(agents_tick_elapsed);
+        self.iteration += 1;
+    }
+
     pub fn run(&mut self, steps: usize) {
         let pb = ProgressBar::new(steps as u64);
-        pb.set_style(
+        pb.set_style(ProgressStyle::default_bar()
-            ProgressStyle::default_bar()
+            .template("{spinner:.green} [{elapsed_precise}] [{bar:40.cyan/blue}] {pos}/{len} ({eta} {percent}%, {per_sec})").expect("invalid progresstyle template")
-                .template(
+            .progress_chars("#>-"));
-                    "{spinner:.green} [{elapsed_precise}] [{bar:40.cyan/blue}] {pos}/{len} ({eta} {percent}%, {per_sec})",
-                )
-                .progress_chars("#>-"),
-        );
 
-        let mut time_per_agent_list: Vec<f64> = Vec::new();
+        for _ in 0..steps {
-        let mut time_per_step_list: Vec<f64> = Vec::new();
+            self.step();
-
-        let agents_num: usize = self
-            .population_grids
-            .iter()
-            .map(|grid| grid.agents.len())
-            .sum();
-        (0..steps).for_each(|_| {
-            // Combine grids
-            combine(&mut self.population_grids, &self.attraction_table);
-            let agents_tick_time = Instant::now();
-
-            // Tick agents
-            self.population_grids.par_iter_mut().for_each(|grid| {
-                grid.tick();
-                grid.diffuse(self.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 / (agents_num as f64);
-            time_per_agent_list.push(ms_per_agent);
-            time_per_step_list.push(agents_tick_elapsed);
-
-            self.iteration += 1;
             pb.inc(1);
-        });
+        }
         pb.finish();
 
         let avg_per_step: f64 =
-            time_per_step_list.iter().sum::<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>() / 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 save_image_data(&mut self) {
+    pub fn population_grids(&self) -> &[Grid] {
-        let grids = ThinGridData::new_from_grid_vec(&self.population_grids);
+        &self.population_grids
-        let img_data = ImgData::new(grids, self.palette);
-        self.img_data_vec.push((self.iteration + 1, img_data));
-        let size: usize = std::mem::size_of_val(&self.img_data_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();
-        }
     }
 
-    pub fn render_all_imgdata(&mut self) {
+    pub fn palette(&self) -> Palette {
-        if !Path::new("./tmp").exists() {
+        self.palette
-            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})",
-        ));
-
-        self.img_data_vec
-            .drain(..)
-            .collect::<Vec<_>>()
-            .par_iter()
-            .progress_with(pb)
-            .for_each(|(i, img)| {
-                img.to_image()
-                    .save(format!("./tmp/out_{}.png", i).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 rng = rand::rng();
-    let mut palette = PALETTES[rng.gen_range(0..PALETTES.len())];
+    let mut palette = PALETTES[rng.random_range(0..PALETTES.len())];
     palette.colors.shuffle(&mut rng);
     palette
 }

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);
+        }
+    }
 }