cleanup

TODO.md

@@ -5,4 +5,7 @@
     - Tried [ArrayFire-rust](https://github.com/arrayfire/arrayfire-rust) didn't work well, looking for another library
     - Try using [emu](https://github.com/calebwin/emu) (seems to be a very good option)
 - sin and cos optimizations
-    - sin/cos table?
+    - sin/cos table?
+- Make colisions for walls of grid
+- Add config and cmd arguments when running the binary to adjust simulation settings
+- Rewrite `grid.rs`

src/blur.rs

@@ -21,8 +21,7 @@ impl Blur {
         }
     }
 
-    /// Blur an image with 2 box filter passes. The result will be written to the src slice, while
+    // 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.
-    /// the buf slice is used as a scratch space.
     pub fn run(
         &mut self,
         src: &mut [f32],
@@ -37,8 +36,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
+    // 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.
-    /// 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;
@@ -54,8 +52,7 @@ impl Blur {
         result
     }
 
-    /// Perform one pass of the 2D box filter of the given radius. The result will be written to the
+    // 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.
-    /// src slice, while the buf slice is used as a scratch space.
     fn box_blur(
         &mut self,
         src: &mut [f32],
@@ -69,15 +66,14 @@ 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;
 
         src.par_chunks_exact(width)
             .zip(dst.par_chunks_exact_mut(width))
             .for_each(|(src_row, dst_row)| {
-                // First we build a value for the beginning of each row. We assume periodic boundary
+                // 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.
-                // 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];
                 for j in 0..radius {
@@ -93,8 +89,7 @@ impl Blur {
             })
     }
 
-    /// Perform one pass of the 1D box filter of the given radius along y axis. Applies the decay
+    // Perform one pass of the 1D box filter of the given radius along y axis. Applies the decay factor to the destination buffer.
-    /// factor to the destination buffer.
     fn box_blur_v(
         &mut self,
         src: &[f32],
@@ -106,9 +101,7 @@ impl Blur {
     ) {
         let weight = decay / (2 * radius + 1) as f32;
 
-        // We don't replicate the horizontal filter logic because of the cache-unfriendly memory
+        // We don't replicate the horizontal filter logic because of the cache-unfriendly memory  access patterns of sequential iteration over individual columns. Instead, we iterate over rows via loop interchange.
-        // access patterns of sequential iteration over individual columns. Instead, we iterate over
-        // rows via loop interchange.
         let height_sub_radius = height - radius;
         let offset = (height_sub_radius - 1) * width;
         self.row_buffer

src/grid.rs

@@ -4,7 +4,7 @@ use rand::{distributions::Uniform, Rng};
 
 use std::fmt::{Display, Formatter};
 
-/// A population configuration.
+// A population configuration.
 #[derive(Debug)]
 pub struct PopulationConfig {
     pub sensor_distance: f32,
@@ -58,7 +58,7 @@ impl PopulationConfig {
     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),
@@ -76,7 +76,7 @@ impl PopulationConfig {
     }
 }
 
-/// 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.
+// 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)]
 pub struct Grid {
     pub config: PopulationConfig,
@@ -104,7 +104,7 @@ impl Clone for Grid {
 }
 
 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) -> Self {
         if !width.is_power_of_two() || !height.is_power_of_two() {
             panic!("Grid dimensions must be a power of two.");
@@ -122,7 +122,7 @@ impl Grid {
         }
     }
 
-    /// 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 {
         // 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);
@@ -130,18 +130,18 @@ impl Grid {
         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)]
     }
 
-    /// 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,

src/main.rs

@@ -5,7 +5,8 @@ fn main() {
     let n_iterations = 2048;
 
     // Size of grid and pictures
-    let (width, height) = (256, 256);
+    // let (width, height) = (256, 256);
+    let (width, height) = (1024, 1024);
 
     // # of agents
     let n_particles = 1 << 20;

src/math.rs

@@ -1,18 +1,13 @@
-#[inline(always)]
+// Previously from trig.rs
-fn to_radians(x: f32) -> f32 {
+// From https://bits.stephan-brumme.com/absFloat.html
-    x * (std::f32::consts::PI / 180.0)
-}
-
-/// Previously from trig.rs
-/// From https://bits.stephan-brumme.com/absFloat.html
 #[allow(dead_code)]
 #[inline(always)]
 fn abs(x: f32) -> f32 {
     return f32::from_bits(x.to_bits() & 0x7FFF_FFFF);
 }
 
-/// Previously from trig.rs
+// Previously from trig.rs
-/// Branchless floor implementation
+// Branchless floor implementation
 #[allow(dead_code)]
 #[inline(always)]
 fn floor(x: f32) -> f32 {
@@ -21,9 +16,9 @@ fn floor(x: f32) -> f32 {
     return x_trunc;
 }
 
-/// Previously from trig.rs
+// Previously from trig.rs
-/// Approximates `cos(x)` in radians with the maximum error of `0.002`
+// Approximates `cos(x)` in radians with the maximum error of `0.002`
-/// https://stackoverflow.com/posts/28050328/revisions
+// https://stackoverflow.com/posts/28050328/revisions
 #[allow(dead_code)]
 #[inline(always)]
 pub fn cos(mut x: f32) -> f32 {
@@ -35,8 +30,8 @@ pub fn cos(mut x: f32) -> f32 {
     return x;
 }
 
-/// Previously from trig.rs
+// Previously from trig.rs
-/// Approximates `sin(x)` in radians with the maximum error of `0.002`
+// Approximates `sin(x)` in radians with the maximum error of `0.002`
 #[allow(dead_code)]
 #[inline(always)]
 pub fn sin(x: f32) -> f32 {

src/model.rs

@@ -11,12 +11,11 @@ use rayon::prelude::*;
 use itertools::multizip;
 use std::f32::consts::TAU;
 use std::time::{Instant};
-use rayon::iter::{ParallelIterator,};
+use rayon::iter::ParallelIterator;
 use indicatif::{ParallelProgressIterator, ProgressBar, ProgressStyle};
 use std::path::Path;
 
-/// A single Physarum agent. The x and y positions are continuous, hence we use floating point
+// A single Physarum agent. The x and y positions are continuous, hence we use floating point numbers instead of integers.
-/// numbers instead of integers.
 #[derive(Debug)]
 struct Agent {
     x: f32,
@@ -27,7 +26,7 @@ struct Agent {
 }
 
 impl Agent {
-    /// Construct a new agent with random parameters.
+    // Construct a new agent with random parameters.
     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)>();
         Agent {
@@ -107,7 +106,7 @@ impl PartialEq for Agent {
 }
 
 
-/// Top-level simulation class.
+// Top-level simulation class.
 pub struct Model {
     // Physarum agents.
     agents: Vec<Agent>,
@@ -144,7 +143,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,
@@ -190,7 +189,7 @@ impl Model {
     }
 
 
-    /// Simulates `steps` # of steps
+    // Simulates `steps` # of steps
     #[inline]
     pub fn run(&mut self, steps: usize) {
         let debug: bool = false;