diff --git a/TODO.md b/TODO.md index 59b1f74..69f2c92 100644 --- a/TODO.md +++ b/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? \ No newline at end of file + - 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` \ No newline at end of file diff --git a/src/blur.rs b/src/blur.rs index 288d00d..2fa9ca9 100644 --- a/src/blur.rs +++ b/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 - /// 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], @@ -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 - /// 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(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 - /// 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], @@ -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 - // conditions, so we need to push the left index to the opposite side of the 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]; 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 - /// 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], @@ -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 - // access patterns of sequential iteration over individual columns. Instead, we iterate over - // rows via loop interchange. + // 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. let height_sub_radius = height - radius; let offset = (height_sub_radius - 1) * width; self.row_buffer diff --git a/src/grid.rs b/src/grid.rs index dd350b8..bebb91a 100644 --- a/src/grid.rs +++ b/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(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(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, diff --git a/src/main.rs b/src/main.rs index fd90602..3407e67 100644 --- a/src/main.rs +++ b/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; diff --git a/src/math.rs b/src/math.rs index b1fb4e3..1f364ba 100644 --- a/src/math.rs +++ b/src/math.rs @@ -1,18 +1,13 @@ -#[inline(always)] -fn to_radians(x: f32) -> f32 { - x * (std::f32::consts::PI / 180.0) -} - -/// Previously from trig.rs -/// From https://bits.stephan-brumme.com/absFloat.html +// 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 -/// Branchless floor implementation +// Previously from trig.rs +// 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 -/// Approximates `cos(x)` in radians with the maximum error of `0.002` -/// https://stackoverflow.com/posts/28050328/revisions +// Previously from trig.rs +// Approximates `cos(x)` in radians with the maximum error of `0.002` +// 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 -/// Approximates `sin(x)` in radians with the maximum error of `0.002` +// Previously from trig.rs +// Approximates `sin(x)` in radians with the maximum error of `0.002` #[allow(dead_code)] #[inline(always)] pub fn sin(x: f32) -> f32 { diff --git a/src/model.rs b/src/model.rs index 5d9b518..5281d63 100644 --- a/src/model.rs +++ b/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 -/// 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)] 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(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, @@ -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;