Blur: update run test

src/agent.rs

@@ -4,7 +4,7 @@ use rand::{seq::SliceRandom, 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,
@@ -21,7 +21,7 @@ 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>(
         width: usize,
         height: usize,
@@ -39,7 +39,7 @@ impl Agent {
         }
     }
 
-    // Tick an agent
+    /// Tick an agent
     pub fn tick(
         &mut self,
         buf: &Buf,

src/blur.rs

@@ -13,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],
@@ -28,7 +28,7 @@ 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 mut w = w_ideal as usize;
@@ -44,7 +44,7 @@ impl Blur {
         result
     }
 
-    // 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],
@@ -58,7 +58,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;
 
@@ -81,7 +81,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],
@@ -277,7 +277,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);
@@ -348,7 +348,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);
@@ -419,7 +419,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");
         }
     }
 
@@ -434,4 +434,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 as 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

@@ -4,7 +4,7 @@ use rand::{distributions::Uniform, Rng};
 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,
@@ -23,7 +23,7 @@ impl Display for PopulationConfig {
 }
 
 impl PopulationConfig {
-    // Construct a random configuration.
+    /// Construct a random configuration.
     pub fn new<R: Rng + ?Sized>(rng: &mut R) -> Self {
         PopulationConfig {
             sensor_distance: rng.gen_range(0.0..=64.0),
@@ -46,23 +46,20 @@ 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() {
-            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();
 
@@ -71,27 +68,24 @@ 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
     }
 
-    // Add a value to the grid data at a given position.
+    /// 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,
@@ -143,13 +137,12 @@ 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);
+        unsafe { buf_ptr.as_mut() }.unwrap().fill(0.0);
         datas.iter().enumerate().for_each(|(j, other)| {
             let multiplier = attraction_table[i].as_ref()[j];
-            buf_ptr
+            unsafe { buf_ptr.as_mut() }
-                .as_mut()
                 .unwrap()
                 .iter_mut()
                 .zip(*other)
@@ -162,13 +155,6 @@ 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();

src/imgdata.rs

@@ -12,7 +12,7 @@ pub struct ThinGridData {
 }
 
 impl ThinGridData {
-    // Convert Grid to ThinGridData
+    /// Convert Grid to ThinGridData
     pub fn new_from_grid(in_grid: &Grid) -> Self {
         ThinGridData {
             width: in_grid.width,
@@ -22,13 +22,10 @@ impl ThinGridData {
     }
 
     pub fn new_from_grid_vec(in_grids: &[Grid]) -> Vec<Self> {
-        in_grids
+        in_grids.iter().map(Self::new_from_grid).collect()
-            .iter()
-            .map(Self::new_from_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

src/model.rs

@@ -10,21 +10,21 @@ use rand_distr::{Distribution, Normal};
 use rayon::{iter::ParallelIterator, prelude::*};
 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,
 
     time_per_agent_list: Vec<f64>,
@@ -44,7 +44,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,
@@ -132,7 +132,6 @@ impl Model {
         );
     }
 
-    // Accessors for rendering
     pub fn population_grids(&self) -> &[Grid] {
         &self.population_grids
     }

src/util.rs

@@ -2,3 +2,12 @@
 pub fn wrap(x: f32, max: f32) -> f32 {
     x - max * ((x > max) as i32 as f32 - (x < 0.0_f32) as i32 as f32)
 }
+
+/// 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
+}