diff --git a/Cargo.lock b/Cargo.lock index 7cd49cb..32a2731 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -1,5 +1,7 @@ # This file is automatically @generated by Cargo. # It is not intended for manual editing. +version = 3 + [[package]] name = "adler" version = "0.2.3" @@ -388,6 +390,12 @@ version = "0.2.86" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "b7282d924be3275cec7f6756ff4121987bc6481325397dde6ba3e7802b1a8b1c" +[[package]] +name = "libm" +version = "0.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c7d73b3f436185384286bd8098d17ec07c9a7d2388a6599f824d8502b529702a" + [[package]] name = "log" version = "0.4.14" @@ -481,6 +489,7 @@ source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "9a64b1ec5cda2586e284722486d802acf1f7dbdc623e2bfc57e65ca1cd099290" dependencies = [ "autocfg", + "libm", ] [[package]] @@ -514,6 +523,7 @@ dependencies = [ "indicatif", "itertools 0.10.0", "rand", + "rand_distr", "rayon", ] @@ -612,6 +622,16 @@ dependencies = [ "getrandom", ] +[[package]] +name = "rand_distr" +version = "0.4.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "da9e8f32ad24fb80d07d2323a9a2ce8b30d68a62b8cb4df88119ff49a698f038" +dependencies = [ + "num-traits", + "rand", +] + [[package]] name = "rand_hc" version = "0.3.0" diff --git a/Cargo.toml b/Cargo.toml index 49831a3..874a91c 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -5,11 +5,12 @@ authors = ["mindv0rtex "] edition = "2018" [dependencies] -image = "*" +image = "0.23" indicatif = "0.15" -itertools = "*" -rand = "*" -rayon = "*" +itertools = "0.10" +rand = "0.8" +rand_distr = "0.4" +rayon = "1.5" [dev-dependencies] criterion = "0.3" diff --git a/src/blur.rs b/src/blur.rs index 10c1eef..37ab8b7 100644 --- a/src/blur.rs +++ b/src/blur.rs @@ -3,27 +3,30 @@ use rayon::prelude::*; #[derive(Debug)] pub struct Blur { - width: usize, - height: usize, row_buffer: Vec, } impl Blur { - pub fn new(width: usize, height: usize) -> Self { + pub fn new(width: usize) -> Self { Blur { - width, - height, row_buffer: vec![0.0; width], } } /// Blur an image with 3 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], buf: &mut [f32], sigma: f32, decay: f32) { - let boxes = Blur::boxes_for_gaussian::<3>(sigma); - self.box_blur(src, buf, boxes[0], 1.0); - self.box_blur(src, buf, boxes[1], 1.0); - self.box_blur(src, buf, boxes[2], decay); + pub fn run( + &mut self, + src: &mut [f32], + buf: &mut [f32], + width: usize, + height: usize, + sigma: f32, + decay: f32, + ) { + let boxes = Blur::boxes_for_gaussian::<2>(sigma); + self.box_blur(src, buf, width, height, boxes[0], 1.0); + self.box_blur(src, buf, width, height, boxes[1], decay); } /// Approximate 1D Gaussian filter of standard deviation sigma with N box filter passes. Each @@ -46,15 +49,22 @@ impl Blur { /// 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], buf: &mut [f32], radius: usize, decay: f32) { - self.box_blur_h(src, buf, radius); - self.box_blur_v(buf, src, radius, decay); + fn box_blur( + &mut self, + src: &mut [f32], + buf: &mut [f32], + width: usize, + height: usize, + radius: usize, + decay: f32, + ) { + self.box_blur_h(src, buf, width, radius); + 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. - fn box_blur_h(&mut self, src: &[f32], dst: &mut [f32], radius: usize) { + fn box_blur_h(&mut self, src: &[f32], dst: &mut [f32], width: usize, radius: usize) { let weight = 1.0 / (2 * radius + 1) as f32; - let width = self.width; src.par_chunks_exact(width) .zip(dst.par_chunks_exact_mut(width)) @@ -66,20 +76,27 @@ impl Blur { value += src_row[width - radius + j] + src_row[j]; } - for i in 0..width { + for (i, dst_elem) in dst_row.iter_mut().enumerate() { let left = (i + width - radius - 1) & (width - 1); let right = (i + radius) & (width - 1); value += src_row[right] - src_row[left]; - dst_row[i] = value * weight; + *dst_elem = value * weight; } }) } /// 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], dst: &mut [f32], radius: usize, decay: f32) { + fn box_blur_v( + &mut self, + src: &[f32], + dst: &mut [f32], + width: usize, + height: usize, + radius: usize, + decay: f32, + ) { let weight = decay / (2 * radius + 1) as f32; - let (width, height) = (self.width, self.height); // 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 @@ -142,9 +159,9 @@ mod tests { ]; let (width, height) = (8, 8); let mut dst = vec![0.0; width * height]; - let mut blur = Blur::new(width, height); + let mut blur = Blur::new(width); - blur.box_blur_h(&src, &mut dst, 1); + blur.box_blur_h(&src, &mut dst, width, 1); let mut sol: Vec = vec![ 0.33921536, 0.13621319, 0.04954382, 0.26381392, 0.46308973, 0.49737768, 0.47066893, 0.37277121, 0.44850051, 0.37332688, 0.21674603, 0.36333409, 0.48751974, 0.70454735, @@ -161,7 +178,7 @@ mod tests { assert!((v1 - v2).abs() < 1e-6); } - blur.box_blur_v(&src, &mut dst, 1, 1.0); + blur.box_blur_v(&src, &mut dst, width, height, 1, 1.0); sol = vec![ 0.50403511, 0.38229549, 0.19629186, 0.29968528, 0.51910173, 0.61901508, 0.44607546, 0.53130095, 0.52355005, 0.177688, 0.16011561, 0.08289763, 0.51645436, 0.46399322, @@ -178,7 +195,7 @@ mod tests { assert!((v1 - v2).abs() < 1e-6); } - blur.box_blur(&mut src, &mut dst, 1, 1.0); + blur.box_blur(&mut src, &mut dst, width, height, 1, 1.0); sol = vec![ 0.47254385, 0.36087415, 0.29275754, 0.33835963, 0.47926736, 0.52806409, 0.5321305, 0.49380384, 0.46566129, 0.28711789, 0.14023375, 0.25315587, 0.3544484, 0.45375601, diff --git a/src/grid.rs b/src/grid.rs index 7902755..8b94d05 100644 --- a/src/grid.rs +++ b/src/grid.rs @@ -1,12 +1,76 @@ -use rand::{distributions::Uniform, Rng}; - use crate::blur::Blur; -/// A 2D grid with a scalar value per each grid block. +use rand::{distributions::Uniform, Rng}; + +use std::fmt::{Display, Formatter}; + +/// A population configuration. +#[derive(Debug)] +pub struct PopulationConfig { + pub sensor_distance: f32, + pub step_distance: f32, + pub sensor_angle: f32, + pub rotation_angle: f32, + + decay_factor: f32, + deposition_amount: f32, +} + +impl Display for PopulationConfig { + fn fmt(&self, f: &mut Formatter<'_>) -> std::fmt::Result { + write!( + f, + "{{\n Sensor Distance: {},\n Step Distance: {},\n Sensor Angle: {},\n Rotation Angle: {},\n Decay Factor: {},\n Deposition Amount: {},\n}}", + self.sensor_distance, + self.step_distance, + self.sensor_angle, + self.rotation_angle, + self.decay_factor, + self.deposition_amount + ) + } +} + +impl PopulationConfig { + const SENSOR_ANGLE_MIN: f32 = 0.0; + 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(rng: &mut R) -> Self { + PopulationConfig { + sensor_distance: rng.gen_range(Self::SENSOR_DISTANCE_MIN..=Self::SENSOR_DISTANCE_MAX), + step_distance: rng.gen_range(Self::STEP_DISTANCE_MIN..=Self::STEP_DISTANCE_MAX), + decay_factor: rng.gen_range(Self::DECAY_FACTOR_MIN..=Self::DECAY_FACTOR_MAX), + sensor_angle: rng + .gen_range(Self::SENSOR_ANGLE_MIN..=Self::SENSOR_ANGLE_MAX) + .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), + } + } +} + +/// 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 { - width: usize, - height: usize, + pub config: PopulationConfig, + pub width: usize, + pub height: usize, + data: Vec, // Scratch space for the blur operation. @@ -16,11 +80,10 @@ pub struct Grid { impl Grid { /// Create a new grid filled with random floats in the [0.0..1.0) range. - pub fn new(width: usize, height: usize) -> Self { + 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."); } - let rng = rand::thread_rng(); let range = Uniform::from(0.0..1.0); let data = rng.sample_iter(range).take(width * height).collect(); @@ -28,8 +91,9 @@ impl Grid { width, height, data, + config: PopulationConfig::new(rng), buf: vec![0.0; width * height], - blur: Blur::new(width, height), + blur: Blur::new(width), } } @@ -41,12 +105,6 @@ impl Grid { j * self.width + i } - /// Get the data value at a given position. The implementation effectively treats data as - /// periodic, hence any finite position will produce a value. - pub fn get(&self, x: f32, y: f32) -> f32 { - self.data[self.index(x, y)] - } - /// 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 { @@ -54,19 +112,25 @@ impl Grid { } /// Add a value to the grid data at a given position. - pub fn add(&mut self, x: f32, y: f32, value: f32) { + pub fn deposit(&mut self, x: f32, y: f32) { let idx = self.index(x, y); - self.data[idx] += value + self.data[idx] += self.config.deposition_amount; } /// Diffuse grid data and apply a decay multiplier. - pub fn diffuse(&mut self, radius: usize, decay_factor: f32) { - self.blur - .run(&mut self.data, &mut self.buf, radius as f32, decay_factor); + pub fn diffuse(&mut self, radius: usize) { + self.blur.run( + &mut self.data, + &mut self.buf, + self.width, + self.height, + radius as f32, + self.config.decay_factor, + ); } pub fn quantile(&self, fraction: f32) -> f32 { - let index = if fraction == 1.0 { + let index = if (fraction - 1.0_f32).abs() < f32::EPSILON { self.data.len() - 1 } else { (self.data.len() as f32 * fraction) as usize @@ -83,6 +147,30 @@ impl Grid { } } +pub fn combine(grids: &mut [Grid], attraction_table: &[T]) +where + T: AsRef<[f32]> + Sync, +{ + let datas: Vec<_> = grids.iter().map(|grid| &grid.data).collect(); + let bufs: Vec<_> = grids.iter().map(|grid| &grid.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 { + let buf_ptr = *buf as *const Vec as *mut Vec; + 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 + .as_mut() + .unwrap() + .iter_mut() + .zip(*other) + .for_each(|(to, from)| *to += from * multiplier) + }) + }); +} + #[cfg(test)] mod tests { use super::*; @@ -90,12 +178,14 @@ mod tests { #[test] #[should_panic] fn test_grid_new_panics() { - let _ = Grid::new(5, 5); + let mut rng = rand::thread_rng(); + let _ = Grid::new(5, 5, &mut rng); } #[test] fn test_grid_new() { - let grid = Grid::new(8, 8); + let mut rng = rand::thread_rng(); + let grid = Grid::new(8, 8, &mut rng); assert_eq!(grid.index(0.5, 0.6), 0); assert_eq!(grid.index(1.5, 0.6), 1); assert_eq!(grid.index(0.5, 1.6), 8); diff --git a/src/main.rs b/src/main.rs index 5a681d5..d645094 100644 --- a/src/main.rs +++ b/src/main.rs @@ -1,5 +1,6 @@ use indicatif::{ProgressBar, ProgressStyle}; use physarum::model; +use rand::Rng; fn main() { let n_iterations = 400; @@ -15,8 +16,9 @@ fn main() { let (width, height) = (1024, 1024); let n_particles = 1 << 22; let diffusivity = 1; - let mut model = model::Model::new(width, height, n_particles, diffusivity); - println!("Model configuration: {:#?}", model.config); + let n_populations = 1 + rand::thread_rng().gen_range(1..4); + let mut model = model::Model::new(width, height, n_particles, n_populations, diffusivity); + model.print_configurations(); for i in 0..n_iterations { model.step(); diff --git a/src/model.rs b/src/model.rs index 20e0f98..2df6762 100644 --- a/src/model.rs +++ b/src/model.rs @@ -1,7 +1,9 @@ -use crate::grid::Grid; +use crate::grid::{combine, Grid, PopulationConfig}; use rand::{seq::SliceRandom, Rng}; +use rand_distr::{Distribution, Normal}; use rayon::prelude::*; + use std::f32::consts::TAU; /// A single Physarum agent. The x and y positions are continuous, hence we use floating point @@ -11,16 +13,18 @@ struct Agent { x: f32, y: f32, angle: f32, + population_id: usize, } impl Agent { /// Construct a new agent with random parameters. - fn new(width: usize, height: usize, rng: &mut R) -> Self { + fn new(width: usize, height: usize, id: usize, rng: &mut R) -> Self { let (x, y, angle) = rng.gen::<(f32, f32, f32)>(); Agent { x: x * width as f32, y: y * height as f32, angle: angle * TAU, + population_id: id, } } @@ -41,50 +45,6 @@ impl Agent { } } -/// A model configuration. We make it into a separate type, because we will eventually have multiple -/// configurations in one model. -#[derive(Debug)] -pub struct PopulationConfig { - sensor_distance: f32, - step_distance: f32, - decay_factor: f32, - sensor_angle: f32, - rotation_angle: f32, - deposition_amount: f32, -} - -impl PopulationConfig { - const SENSOR_ANGLE_MIN: f32 = 0.0; - 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(rng: &mut R) -> Self { - PopulationConfig { - sensor_distance: rng.gen_range(Self::SENSOR_DISTANCE_MIN..=Self::SENSOR_DISTANCE_MAX), - step_distance: rng.gen_range(Self::STEP_DISTANCE_MIN..=Self::STEP_DISTANCE_MAX), - decay_factor: rng.gen_range(Self::DECAY_FACTOR_MIN..=Self::DECAY_FACTOR_MAX), - sensor_angle: rng - .gen_range(Self::SENSOR_ANGLE_MIN..=Self::SENSOR_ANGLE_MAX) - .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), - } - } -} - /// Top-level simulation class. #[derive(Debug)] pub struct Model { @@ -92,32 +52,71 @@ pub struct Model { agents: Vec, // The grid they move on. - grid: Grid, + grids: Vec, - // Simulation parameters. + // Attraction table governs interaction across populations + attraction_table: Vec>, + + // Global grid diffusivity. diffusivity: usize, - pub config: PopulationConfig, + // Current model iteration. iteration: i32, - width: usize, - height: usize, } impl Model { + const ATTRACTION_FACTOR_MEAN: f32 = 1.0; + const ATTRACTION_FACTOR_STD: f32 = 0.1; + const REPULSION_FACTOR_MEAN: f32 = -1.0; + const REPULSION_FACTOR_STD: f32 = 0.1; + + pub fn print_configurations(&self) { + for (i, grid) in self.grids.iter().enumerate() { + println!("Grid {}: {}", i, grid.config); + } + println!("Attraction table: {:#?}", self.attraction_table); + } + /// Construct a new model with random initial conditions and random configuration. - pub fn new(width: usize, height: usize, n_particles: usize, diffusivity: usize) -> Self { + pub fn new( + width: usize, + height: usize, + n_particles: usize, + n_populations: usize, + diffusivity: usize, + ) -> Self { + 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 attraction_distr = + Normal::new(Self::ATTRACTION_FACTOR_MEAN, Self::ATTRACTION_FACTOR_STD).unwrap(); + let repulstion_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) + }); + } + } + Model { agents: (0..n_particles) - .map(|_| Agent::new(width, height, &mut rng)) + .map(|i| Agent::new(width, height, i / particles_per_grid, &mut rng)) .collect(), - grid: Grid::new(width, height), + grids: (0..n_populations) + .map(|_| Grid::new(width, height, &mut rng)) + .collect(), + attraction_table, diffusivity, - config: PopulationConfig::new(&mut rng), iteration: 0, - width, - height, } } @@ -137,19 +136,22 @@ impl Model { /// Perform a single simulation step. pub fn step(&mut self) { - // To avoid borrow-checker errors inside the parallel loop. - let PopulationConfig { - sensor_distance, - sensor_angle, - rotation_angle, - step_distance, - .. - } = self.config; - let (width, height) = (self.width, self.height); - let grid = &self.grid; + // Combine grids + let grids = &mut self.grids; + let attraction_table = &self.attraction_table; + combine(grids, attraction_table); self.agents.par_iter_mut().for_each(|agent| { - let mut rng = rand::thread_rng(); + let grid = &grids[agent.population_id]; + let PopulationConfig { + sensor_distance, + sensor_angle, + rotation_angle, + step_distance, + .. + } = grid.config; + let (width, height) = (grid.width, grid.height); + let xc = agent.x + agent.angle.cos() * sensor_distance; let yc = agent.y + agent.angle.sin() * sensor_distance; let xl = agent.x + (agent.angle - sensor_angle).cos() * sensor_distance; @@ -158,34 +160,38 @@ impl Model { let yr = agent.y + (agent.angle + sensor_angle).sin() * sensor_distance; // Sense - let trail_c = grid.get(xc, yc); - let trail_l = grid.get(xl, yl); - let trail_r = grid.get(xr, yr); + let trail_c = grid.get_buf(xc, yc); + let trail_l = grid.get_buf(xl, yl); + let trail_r = grid.get_buf(xr, yr); // Rotate and move + let mut rng = rand::thread_rng(); let direction = Model::pick_direction(trail_c, trail_l, trail_r, &mut rng); agent.rotate_and_move(direction, rotation_angle, step_distance, width, height); }); // Deposit for agent in self.agents.iter() { - self.grid - .add(agent.x, agent.y, self.config.deposition_amount); + self.grids[agent.population_id].deposit(agent.x, agent.y); } + // Diffuse + Decay - self.grid - .diffuse(self.diffusivity, self.config.decay_factor); + let diffusivity = self.diffusivity; + self.grids.par_iter_mut().for_each(|grid| { + grid.diffuse(diffusivity); + }); self.iteration += 1; } /// Output the current trail layer as a grayscale image. pub fn save_to_image(&self, name: &str) { - let mut img = image::GrayImage::new(self.width as u32, self.height as u32); - let max_value = self.grid.quantile(0.999); + let mut img = + image::GrayImage::new(self.grids[0].width as u32, self.grids[0].height as u32); + let max_value = self.grids[0].quantile(0.999); - for (i, value) in self.grid.data().iter().enumerate() { - let x = (i % self.width) as u32; - let y = (i / self.width) as u32; + for (i, value) in self.grids[0].data().iter().enumerate() { + let x = (i % self.grids[0].width) as u32; + let y = (i / self.grids[0].width) as u32; let c = (value / max_value).clamp(0.0, 1.0) * 255.0; img.put_pixel(x, y, image::Luma([c as u8])); }