Support multiple populations.

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"

Cargo.toml

@@ -5,11 +5,12 @@ authors = ["mindv0rtex <mindv0rtex@users.noreply.github.com>"]
 edition = "2018"
 
 [dependencies]
-image = "*"
+image = "0.23"
 indicatif = "0.15"
-itertools = "*"
+itertools = "0.10"
-rand = "*"
+rand = "0.8"
-rayon = "*"
+rand_distr = "0.4"
+rayon = "1.5"
 
 [dev-dependencies]
 criterion = "0.3"

src/blur.rs

@@ -3,27 +3,30 @@ use rayon::prelude::*;
 
 #[derive(Debug)]
 pub struct Blur {
-    width: usize,
-    height: usize,
     row_buffer: Vec<f32>,
 }
 
 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) {
+    pub fn run(
-        let boxes = Blur::boxes_for_gaussian::<3>(sigma);
+        &mut self,
-        self.box_blur(src, buf, boxes[0], 1.0);
+        src: &mut [f32],
-        self.box_blur(src, buf, boxes[1], 1.0);
+        buf: &mut [f32],
-        self.box_blur(src, buf, boxes[2], decay);
+        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) {
+    fn box_blur(
-        self.box_blur_h(src, buf, radius);
+        &mut self,
-        self.box_blur_v(buf, src, radius, decay);
+        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<f32> = 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,

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<R: Rng + ?Sized>(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,
+    pub config: PopulationConfig,
-    height: usize,
+    pub width: usize,
+    pub height: usize,
+
     data: Vec<f32>,
 
     // 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<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.");
         }
-        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) {
+    pub fn diffuse(&mut self, radius: usize) {
-        self.blur
+        self.blur.run(
-            .run(&mut self.data, &mut self.buf, radius as f32, decay_factor);
+            &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<T>(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<f32> as *mut Vec<f32>;
+        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);

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);
+    let n_populations = 1 + rand::thread_rng().gen_range(1..4);
-    println!("Model configuration: {:#?}", model.config);
+    let mut model = model::Model::new(width, height, n_particles, n_populations, diffusivity);
+    model.print_configurations();
 
     for i in 0..n_iterations {
         model.step();

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<R: Rng + ?Sized>(width: usize, height: usize, rng: &mut R) -> Self {
+    fn new<R: Rng + ?Sized>(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<R: Rng + ?Sized>(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<Agent>,
 
     // The grid they move on.
-    grid: Grid,
+    grids: Vec<Grid>,
 
-    // Simulation parameters.
+    // Attraction table governs interaction across populations
+    attraction_table: Vec<Vec<f32>>,
+
+    // 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.
+        // Combine grids
-        let PopulationConfig {
+        let grids = &mut self.grids;
-            sensor_distance,
+        let attraction_table = &self.attraction_table;
-            sensor_angle,
+        combine(grids, attraction_table);
-            rotation_angle,
-            step_distance,
-            ..
-        } = self.config;
-        let (width, height) = (self.width, self.height);
-        let grid = &self.grid;
 
         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_c = grid.get_buf(xc, yc);
-            let trail_l = grid.get(xl, yl);
+            let trail_l = grid.get_buf(xl, yl);
-            let trail_r = grid.get(xr, yr);
+            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
+            self.grids[agent.population_id].deposit(agent.x, agent.y);
-                .add(agent.x, agent.y, self.config.deposition_amount);
         }
+
         // Diffuse + Decay
-        self.grid
+        let diffusivity = self.diffusivity;
-            .diffuse(self.diffusivity, self.config.decay_factor);
+        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 mut img =
-        let max_value = self.grid.quantile(0.999);
+            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() {
+        for (i, value) in self.grids[0].data().iter().enumerate() {
-            let x = (i % self.width) as u32;
+            let x = (i % self.grids[0].width) as u32;
-            let y = (i / self.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]));
         }