use av1 for encoding

Cargo.toml

@@ -5,14 +5,21 @@ authors = ["Simon Gardling <titaniumtown@gmail.com>", "mindv0rtex <mindv0rtex@us
 edition = "2021"
 
 [dependencies]
-image = "0.23"
-indicatif = { version = "0.15", features = [ "rayon" ] }
-itertools = "0.10"
-rand = "0.8"
-rand_distr = "0.4"
+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 = { 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,6 +1,8 @@
+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};
 
@@ -9,9 +11,7 @@ use std::fmt::{Display, Formatter};
 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 {
@@ -22,44 +22,27 @@ impl Display for Agent {
 
 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)>();
+    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,
         }
     }
 
     /// Tick an agent
-    pub fn tick(
-        &mut self,
-        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;
+    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 agent_add_sens = self.angle + sensor_angle;
-        let agent_sub_sens = self.angle - sensor_angle;
+        let agent_add_sens = self.heading + pop_config.sensor_angle;
+        let agent_sub_sens = self.heading - pop_config.sensor_angle;
 
-        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);
@@ -67,25 +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)
-                .expect("unable to choose random direction");
+            *[-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.x = wrap(self.x + step_distance * cos(self.angle), width as f32);
-        self.y = wrap(self.y + step_distance * sin(self.angle), height as f32);
+        self.heading = wrap(self.heading + delta_angle, TAU);
+        self.x = wrap(
+            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

@@ -30,18 +30,17 @@ impl Blur {
 
     /// 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;
-        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.
@@ -67,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];
                 }
@@ -449,7 +448,7 @@ mod tests {
             &mut vec![0.0; width * height],
             width,
             height,
-            2 as f32,
+            2_f32,
             0.1,
         );
 

src/grid.rs

@@ -1,6 +1,6 @@
 use crate::{agent::Agent, blur::Blur, buffer::Buf};
-
-use rand::{distributions::Uniform, Rng};
+use rand::Rng;
+use rand_distr::Uniform;
 use rayon::{iter::ParallelIterator, prelude::*};
 use std::fmt::{Display, Formatter};
 
@@ -12,7 +12,6 @@ pub struct PopulationConfig {
     pub sensor_angle: f32,
     pub rotation_angle: f32,
 
-    decay_factor: f32,
     deposition_amount: f32,
 }
 
@@ -26,12 +25,11 @@ impl PopulationConfig {
     /// Construct a random configuration.
     pub fn new<R: Rng + ?Sized>(rng: &mut R) -> Self {
         PopulationConfig {
-            sensor_distance: rng.gen_range(0.0..=64.0),
-            step_distance: rng.gen_range(0.2..=2.0),
-            decay_factor: rng.gen_range(0.1..=0.1),
-            sensor_angle: rng.gen_range(0.0_f32..=120.0).to_radians(),
-            rotation_angle: rng.gen_range(0.0_f32..=120.0).to_radians(),
-            deposition_amount: rng.gen_range(5.0..=5.0),
+            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),
         }
     }
 }
@@ -60,7 +58,7 @@ impl Grid {
         rng: &mut R,
         agents: Vec<Agent>,
     ) -> Self {
-        let range = Uniform::from(0.0..1.0);
+        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 {
@@ -79,12 +77,6 @@ impl Grid {
         crate::util::index(self.width, self.height, 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(
@@ -93,38 +85,21 @@ impl Grid {
             self.width,
             self.height,
             radius as f32,
-            self.config.decay_factor,
+            0.1, // decay is always 0.1
         );
     }
 
     pub fn tick(&mut self) {
-        let (width, height) = (self.width, self.height);
-        let PopulationConfig {
-            sensor_distance,
-            sensor_angle,
-            rotation_angle,
-            step_distance,
-            ..
-        } = self.config;
-
         self.agents.par_iter_mut().for_each(|agent| {
-            agent.tick(
-                &self.buf,
-                sensor_distance,
-                sensor_angle,
-                rotation_angle,
-                step_distance,
-                width,
-                height,
-            );
+            agent.tick(&self.buf, self.config, self.width, self.height);
         });
         self.deposit_all();
     }
 
     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;
         }
     }
 }
@@ -139,11 +114,14 @@ where
     // 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)| {
         let buf_ptr = *buf as *const Vec<f32> as *mut Vec<f32>;
-        unsafe { 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];
-            unsafe { buf_ptr.as_mut() }
-                .unwrap()
+            buf_ptr_mut
                 .iter_mut()
                 .zip(*other)
                 .for_each(|(to, from)| *to += from * multiplier)
@@ -157,7 +135,7 @@ mod tests {
 
     #[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,28 +1,27 @@
 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,
-    pub height: usize,
-    pub data: Vec<f32>,
+    width: usize,
+    height: usize,
+    data: Vec<f32>,
 }
 
 impl 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.iter().map(Self::new_from_grid).collect()
+        in_grids.iter().cloned().map(Self::new_from_grid).collect()
     }
 
     /// from grid.rs (needed in image gen)
@@ -43,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>,
-    pub palette: Palette,
+    grids: Vec<ThinGridData>,
+    palette: Palette,
 }
 
 impl ImgData {

src/lib.rs

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

src/main.rs

@@ -16,7 +16,7 @@ fn main() {
 
     // Setup ffmpeg
     let mut ffmpeg = std::process::Command::new("ffmpeg")
-        .args(&[
+        .args([
             "-y",
             "-f",
             "rawvideo",
@@ -29,9 +29,7 @@ fn main() {
             "-i",
             "-",
             "-c:v",
-            "libx264",
-            "-preset",
-            "fast",
+            "libsvtav1",
             "output.mp4",
         ])
         .stdin(std::process::Stdio::piped())

src/model.rs

@@ -3,9 +3,7 @@ use crate::{
     grid::{combine, Grid},
     palette::{random_palette, Palette},
 };
-
 use indicatif::{ProgressBar, ProgressStyle};
-// use rand::Rng;
 use rand_distr::{Distribution, Normal};
 use rayon::{iter::ParallelIterator, prelude::*};
 use std::time::Instant;
@@ -55,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_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 {
             population_grids: grids,
@@ -113,7 +115,7 @@ impl Model {
     pub fn run(&mut self, steps: usize) {
         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})")
+            .template("{spinner:.green} [{elapsed_precise}] [{bar:40.cyan/blue}] {pos}/{len} ({eta} {percent}%, {per_sec})").expect("invalid progresstyle template")
             .progress_chars("#>-"));
 
         for _ in 0..steps {

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
 }

src/util.rs

@@ -1,6 +1,7 @@
 #[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.
@@ -11,3 +12,100 @@ pub const fn index(width: usize, height: usize, x: f32, y: f32) -> usize {
     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);
+        }
+    }
+}