Add a README

.gitignore

@@ -1,3 +1,6 @@
 /target
 .idea
 *.csv
+/*.json
+*.svg
+/data

Cargo.toml

@@ -5,18 +5,13 @@ edition = "2021"
 
 # See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html
 
-[lib]
-name = "dla"
-crate-type = ["staticlib"]
-path = "src/clib.rs"
-
 [[bin]]
 name = "model"
 path = "src/main.rs"
 
 [[bin]]
-name = "ui"
-path = "src/ui.rs"
+name = "tools"
+path = "src/tools_cli.rs"
 
 # Set the default for crate.
 [profile.dev]
@@ -32,12 +27,35 @@ opt-level = 3
 
 [dependencies]
 clap = { version = "4.1.8", features = ["derive"] }
-bevy = { version = "0.9.1" }
 nd_array = "0.1.0"
 num-integer = "0.1.45"
 rand = { version = "0.8.5", features = ["default", "small_rng"] }
 csv = "1.1"
 serde = { version = "1.0.152", features = ["derive"] }
+serde_json = "1.0.93"
+kd-tree = { version = "0.5.1", features = ["nalgebra"] }
+nalgebra = { version = "0.32.2", features = ["serde-serialize"] }
+kiddo = "0.2.5"
+anyhow = "1.0.69"
+itertools = "0.10.5"
+svg = "0.13.0"
+polars = { version = "0.27.2", features = [
+    "docs",
+    "zip_with",
+    "csv-file",
+    "temporal",
+    "fmt",
+    "dtype-slim",
+    "csv-file",
+    "lazy",
+    "nightly",
+    "performant"
+] }
+rayon = "1.7.0"
+walkdir = "2.3.3"
+parquet = { version = "35.0.0", features = ["serde"] }
+colorous = "1.0.10"
+num-traits = "0.2.15"
 
 [build-dependencies]
 cbindgen = "0.24.3"

README.md

@@ -0,0 +1,8 @@
+# DLA Generic Model
+
+A generic pluggable model for diffusion limited aggregation. Produces two executables,
+
+- `./target/release/model`, built from `./src/main.rs`
+- `./target/release/tools`, built from `./src/tools_cli.rs`
+
+Build with `cargo build --release`. Requires a working rust installation.

build.rs

@@ -1,16 +0,0 @@
-extern crate cbindgen;
-
-use std::env;
-use std::path::Path;
-use cbindgen::{Config, Builder};
-
-fn main() {
-    let crate_env = env::var("CARGO_MANIFEST_DIR").unwrap();
-    let crate_path = Path::new(&crate_env);
-    let config = Config::from_root_or_default(crate_path);
-    Builder::new().with_crate(crate_path.to_str().unwrap())
-        .with_config(config)
-        .generate()
-        .expect("Cannot generate header file!")
-        .write_to_file("libdla.h");
-}

libdla.h

@@ -1,24 +0,0 @@
-#include <cstdarg>
-#include <cstdint>
-#include <cstdlib>
-#include <ostream>
-#include <new>
-
-struct CStorage;
-
-struct CPosition {
-  int32_t _0;
-  int32_t _1;
-};
-
-extern "C" {
-
-CStorage *storage_new(uint32_t grid_size);
-
-bool storage_at(const CStorage *storage, int32_t i, int32_t j);
-
-void storage_deposit(CStorage *storage, int32_t i, int32_t j, uint8_t val);
-
-CPosition walk(uint32_t d, int32_t i, int32_t j);
-
-} // extern "C"

src/cli/cli.rs

@@ -0,0 +1,88 @@
+use std::path::PathBuf;
+use clap::{Parser, Args, Subcommand, ValueEnum};
+
+#[derive(ValueEnum, Clone, Debug, Copy)]
+pub enum OutputFormat {
+    FullDataJson,
+    Positions,
+}
+
+#[derive(Args, Debug)]
+pub struct InitialCli {
+    pub grid_size: u32,
+}
+
+#[derive(Args, Debug)]
+pub struct StickProbabilityCli {
+    pub grid_size: u32,
+    pub stick_probability: f32,
+}
+
+#[derive(Args, Debug)]
+pub struct KDStickProbabilityCli {
+    pub stick_probability: f32,
+}
+
+#[derive(Args, Debug)]
+pub struct BallsCli {
+    pub ball_radius: f32,
+    pub stick_distance: f32,
+    pub walk_step: f32,
+}
+
+#[derive(Args, Debug)]
+pub struct SurfaceProbabilityMeasureCli {
+    pub grid_size: u32,
+    pub stick_probability: f32,
+    pub initial_data: PathBuf,
+    pub particles: u32,
+}
+
+#[derive(Subcommand, Debug)]
+pub enum PCM {
+    /// Traditional DLA model in a fixed-size 2D grid
+    Initial(InitialCli),
+
+    /// Traditional DLA model in a fixed-size 2D grid with a probabilistic sticking component
+    StickProbability(StickProbabilityCli),
+
+    /// Traditional DLA model in unbounded 3D space with probabilistic sticking
+    Grid3d(KDStickProbabilityCli),
+
+    /// Traditional DLA model in unbounded 3D space with probabilistic sticking and off-axis walks
+    Grid3dOffAxis(KDStickProbabilityCli),
+
+    Hex(StickProbabilityCli),
+
+    Balls2d(BallsCli),
+    Balls(BallsCli),
+
+    SurfaceProbabilityMeasure(SurfaceProbabilityMeasureCli),
+}
+
+#[derive(Parser, Debug)]
+pub struct ModelCli {
+    #[command(subcommand)]
+    pub preconfigured_model: PCM,
+
+    #[arg(long, help = "Set the maximum number of ticks the system will run before exiting, will still produce data after this has been reached.")]
+    pub max_frames: Option<usize>,
+
+    #[arg(short = 'N', long)]
+    pub max_particles: usize,
+
+    #[arg(long, short, help = "Specifying a seed allows for deterministic runs (when combined with other parameters such as -N)")]
+    pub seed: u64,
+
+    #[arg(value_enum, short, long, default_value_t = OutputFormat::Positions)]
+    pub format: OutputFormat,
+
+    #[arg(value_enum, short, long)]
+    pub output: PathBuf,
+
+    #[arg(long, default_value_t = false, help = "If the output file already exists, skipping running. Useful when running deterministically.")]
+    pub preserve_existing: bool,
+
+    #[arg(short, long)]
+    pub notify_every: Option<usize>,
+}

src/cli/mod.rs

@@ -0,0 +1,45 @@
+use rand::Rng;
+use crate::system::model::DLASystem;
+use crate::system::{Position, Storage};
+use crate::system::spawner::Spawner;
+use crate::system::sticker::Sticker;
+use crate::system::walker::Walker;
+use std::time::SystemTime;
+
+pub mod cli;
+pub mod output;
+
+pub fn drive_system<R: Rng, P: Position, S: Storage<P>, W: Walker<P>, Sp: Spawner<P>, St: Sticker<P, S>>(
+    sys: &mut DLASystem<R, P, S, W, Sp, St>,
+    max_frames: Option<usize>,
+    notify_every: Option<usize>,
+) {
+    let start = SystemTime::now();
+    let mut prev = start.clone();
+    let mut previous_n: usize = 0;
+    while sys.running {
+        sys.update();
+
+        if let Some(notify_every) = notify_every && (sys.history.len() % notify_every) == 0 && previous_n != sys.history.len() {
+            let now = SystemTime::now();
+
+            println!("[{}ms, d = {}ms] On frame {}, deposited {} particles",
+                     now.duration_since(start).unwrap().as_millis(),
+                     now.duration_since(prev).unwrap().as_millis(),
+                     sys.frame, sys.history.len()
+            );
+
+            prev = now;
+            previous_n = sys.history.len();
+        }
+
+        match max_frames {
+            Some(max_frames) if max_frames <= sys.frame => {
+                sys.running = false;
+                eprintln!("System halted as it ran to {frame} frames (max_frames = {max_frames}) and did not complete", frame = sys.frame)
+            }
+            _ => {}
+        }
+    }
+}
+

src/cli/output.rs

@@ -0,0 +1,55 @@
+use std::fs::{create_dir_all, File};
+use std::path::Path;
+use rand::Rng;
+use crate::cli::cli::OutputFormat;
+use crate::system::model::DLASystem;
+use crate::system::{Position, Storage};
+use crate::system::spawner::Spawner;
+use crate::system::sticker::Sticker;
+use crate::system::walker::Walker;
+
+pub fn write<R: Rng, P: Position, S: Storage<P>, W: Walker<P>, Sp: Spawner<P>, St: Sticker<P, S>>(
+    sys: &DLASystem<R, P, S, W, Sp, St>,
+    format: OutputFormat,
+    output: &Path,
+) {
+    // If the parent does not exist (and we are not in the root or operating a relative 1-component-length path))
+    // create it. This greatly simplifies the harness code.
+    if let Some(parent) = output.parent() &&
+        parent.to_str().map(|x| x != "").unwrap_or(true) &&
+        !parent.exists() {
+        create_dir_all(parent).expect("Failed to create path to output");
+    }
+
+    match format {
+        OutputFormat::FullDataJson => write_json_full_data(sys, output),
+        OutputFormat::Positions => write_csv_positions(sys, output),
+    }
+}
+
+fn write_csv_positions<R: Rng, P: Position, S: Storage<P>, W: Walker<P>, Sp: Spawner<P>, St: Sticker<P, S>>(sys: &DLASystem<R, P, S, W, Sp, St>, csv_path: &Path) {
+    let mut wtr = csv::Writer::from_path(csv_path)
+        .expect("Failed to open file");
+
+    // CSVs can only store the raw positions
+    let positions: Vec<&P> = sys.history
+        .iter()
+        .map(|line| &line.position)
+        .collect();
+
+    positions
+        .iter()
+        .for_each(|position|
+            wtr.serialize(position).expect("Failed to write row")
+        );
+
+    wtr.flush()
+        .unwrap();
+}
+
+fn write_json_full_data<R: Rng, P: Position, S: Storage<P>, W: Walker<P>, Sp: Spawner<P>, St: Sticker<P, S>>(sys: &DLASystem<R, P, S, W, Sp, St>, output_path: &Path) {
+    let file = File::create(output_path).expect("Failed to open file");
+
+    serde_json::to_writer_pretty(file, &sys.history)
+        .expect("Failed to write json");
+}

src/clib.rs

@@ -1,75 +0,0 @@
-#![feature(array_zip)]
-
-use system::Storage;
-use crate::system::grid::{Position, VectorStorage};
-
-mod system;
-
-#[derive(Eq, PartialEq, Debug)]
-#[repr(C)]
-pub struct CPosition(i32, i32);
-
-pub struct CStorage(VectorStorage);
-
-#[no_mangle]
-pub extern "C" fn storage_new(grid_size: u32) -> &'static mut CStorage {
-    let mut pntr = Box::new(CStorage(VectorStorage::new(grid_size)));
-    Box::leak(pntr)
-}
-
-#[no_mangle]
-pub extern "C" fn storage_at(storage: &CStorage, i: i32, j: i32) -> bool {
-    storage.0.at(&Position { x: i, y: j })
-}
-
-#[no_mangle]
-pub extern "C" fn storage_deposit(storage: &mut CStorage, i: i32, j: i32, val: u8) {
-    storage.0.write(&Position { x: i, y: j }, val == 1);
-}
-
-#[no_mangle]
-pub extern "C" fn walk(d: u32, i: i32, j: i32) -> CPosition {
-    return test::b(d, i, j);
-}
-
-mod test {
-    use num_integer::Integer;
-    use crate::CPosition;
-    use crate::system::grid::Position;
-
-    pub(crate) fn a(d: u32, i: i32, j: i32) -> CPosition {
-        match d {
-            0 => CPosition(i + 1, j),
-            1 => CPosition(i - 1, j),
-            2 => CPosition(i, j + 1),
-            3 => CPosition(i, j - 1),
-            _ => panic!("Ahh"),
-        }
-    }
-
-    pub(crate) fn b(d: u32, i: i32, j: i32) -> CPosition {
-        let (dim, sign) = d.div_rem(&2);
-        let sign = if sign == 0 { 1 } else { -1 };
-        // HACK: Our conventin and the MVA are different, since we are trying to strangle fig this, quick hack.
-        let offset = Position::in_direction(dim, sign);
-        let next = Position { x: i, y: j } + offset;
-
-        CPosition(next.x, next.y)
-    }
-
-    #[test]
-    fn test() {
-        let d = [0, 1, 2, 3];
-        d.iter()
-            .map(|d| d.div_rem(&2))
-            .for_each(|p| println!("{p:?}"));
-    }
-
-    #[test]
-    fn alignment() {
-        let d = [0, 1, 2, 3];
-        d.iter()
-            .map(|d| (a(*d, 0, 0), b(*d, 0, 0)))
-            .for_each(|p| assert_eq!(p.0, p.1));
-    }
-}

src/example_systems.rs

@@ -1,36 +0,0 @@
-use rand::rngs::SmallRng;
-use rand::SeedableRng;
-use crate::system::grid::{Position, VectorStorage};
-use crate::system::model::DLASystem;
-use crate::system::nd::{NDPosition, NDVectorStorage};
-use crate::system::walker::LocalRandomWalker;
-
-pub fn initial_config(seed: u64, max_particles: usize) -> DLASystem<SmallRng, Position, VectorStorage, LocalRandomWalker> {
-    DLASystem::new_g(
-        SmallRng::seed_from_u64(seed),
-        VectorStorage::new(1600),
-        LocalRandomWalker,
-        1.0,
-        max_particles,
-    )
-}
-
-pub fn stick_probability(seed: u64, max_particles: usize, stick_probability: f32) -> DLASystem<SmallRng, Position, VectorStorage, LocalRandomWalker> {
-    DLASystem::new_g(
-        SmallRng::seed_from_u64(seed),
-        VectorStorage::new(1600),
-        LocalRandomWalker,
-        stick_probability,
-        max_particles,
-    )
-}
-
-pub fn three_dimensional(seed: u64, max_particles: usize, stick_probability: f32) -> DLASystem<SmallRng, NDPosition<3>, NDVectorStorage<3>, LocalRandomWalker> {
-    DLASystem::new_g(
-        SmallRng::seed_from_u64(seed),
-        NDVectorStorage::new(1600),
-        LocalRandomWalker,
-        stick_probability,
-        max_particles,
-    )
-}

src/main.rs

@@ -1,36 +1,171 @@
 #![feature(array_zip)]
-
-use std::path::PathBuf;
-
-mod system;
-mod example_systems;
+#![feature(generic_const_exprs)]
+#![feature(let_chains)]
 
 use clap::Parser;
-use crate::example_systems::stick_probability;
+use rand::prelude::*;
+use crate::cli::{drive_system};
+use crate::cli::cli::{StickProbabilityCli, InitialCli, BallsCli, PCM, ModelCli, SurfaceProbabilityMeasureCli, KDStickProbabilityCli};
+use crate::cli::cli::PCM::Balls2d;
+use crate::cli::output::write;
+use crate::surface_probability_measure::{LoggerSticker, ReadOnlyVectorStorage};
+use crate::system::model::DLASystem;
+use crate::system::spaces::continuous_3d::{ContinuousSticker, ContinuousStorage, ContinuousWalker};
+use crate::system::spaces::hexagonal::HexPosition;
+use crate::system::spaces::kd_grid::{KDSpace};
+use crate::system::spaces::square_grid::{Grid2D, Grid3D};
+use crate::system::spaces::VectorStorage;
+use crate::system::spawner::UniformSpawner;
+use crate::system::sticker::{ProbabilisticSticking, SimpleSticking};
+use crate::system::walker::{DiagonalRandomWalker, LocalRandomWalker, Walker};
 
-#[derive(Parser, Debug)]
-struct Cli {
-    seed: u64,
-    max_particles: usize,
-    stick_probability: f32,
-    csv_path: PathBuf,
-}
+mod system;
+mod surface_probability_measure;
+mod cli;
 
 fn main() {
-    let cli = Cli::parse();
+    let cli = ModelCli::parse();
+
+    if cli.preserve_existing && cli.output.exists() {
+        println!("Skipping model execution as output file '{}' already exists", cli.output.to_str().unwrap());
+        return;
+    }
 
     println!("Running: {:?}", cli);
 
-    let mut sys = stick_probability(
-        cli.seed,
-        cli.max_particles,
-        cli.stick_probability
-    );
+    match cli.preconfigured_model {
+        PCM::Initial(InitialCli { grid_size }) => {
+            let mut sys = DLASystem::<_, Grid2D, _, _, _, _>::new(
+                SmallRng::seed_from_u64(cli.seed),
+                VectorStorage::new(grid_size, 2),
+                LocalRandomWalker,
+                UniformSpawner,
+                SimpleSticking,
+                cli.max_particles,
+            );
 
-    while sys.running {
-        sys.update();
+            drive_system(&mut sys, cli.max_frames, cli.notify_every);
+            write(&sys, cli.format, &cli.output);
+        }
+
+        PCM::StickProbability(StickProbabilityCli { grid_size, stick_probability }) => {
+            let mut sys = DLASystem::<_, Grid2D, _, _, _, _>::new(
+                SmallRng::seed_from_u64(cli.seed),
+                VectorStorage::new(grid_size, 2),
+                LocalRandomWalker,
+                UniformSpawner,
+                ProbabilisticSticking::new(stick_probability).unwrap(),
+                cli.max_particles,
+            );
+
+            drive_system(&mut sys, cli.max_frames, cli.notify_every);
+            write(&sys, cli.format, &cli.output);
+        }
+
+        PCM::Grid3d(KDStickProbabilityCli { stick_probability }) => {
+            let mut sys = DLASystem::<_, Grid3D, _, _, _, _>::new(
+                SmallRng::seed_from_u64(cli.seed),
+                KDSpace::new(),
+                LocalRandomWalker,
+                UniformSpawner,
+                ProbabilisticSticking::new(stick_probability).unwrap(),
+                cli.max_particles,
+            );
+
+            drive_system(&mut sys, cli.max_frames, cli.notify_every);
+            write(&sys, cli.format, &cli.output);
+        }
+
+        PCM::Grid3dOffAxis(KDStickProbabilityCli { stick_probability }) => {
+            let mut sys = DLASystem::<_, Grid3D, _, _, _, _>::new(
+                SmallRng::seed_from_u64(cli.seed),
+                KDSpace::new(),
+                DiagonalRandomWalker,
+                UniformSpawner,
+                ProbabilisticSticking::new(stick_probability).unwrap(),
+                cli.max_particles,
+            );
+
+            drive_system(&mut sys, cli.max_frames, cli.notify_every);
+            write(&sys, cli.format, &cli.output);
+        }
+
+        PCM::Hex(StickProbabilityCli { grid_size, stick_probability }) => {
+            let mut sys = DLASystem::<_, HexPosition, _, _, _, _>::new(
+                SmallRng::seed_from_u64(cli.seed),
+                VectorStorage::new(grid_size, 2),
+                LocalRandomWalker,
+                UniformSpawner,
+                ProbabilisticSticking::new(stick_probability).unwrap(),
+                cli.max_particles,
+            );
+
+            drive_system(&mut sys, cli.max_frames, cli.notify_every);
+            write(&sys, cli.format, &cli.output);
+        }
+
+        PCM::Balls2d(BallsCli { ball_radius, stick_distance, walk_step }) => {
+            use system::spaces::continuous_2d;
+
+            let mut sys = DLASystem::new(
+                SmallRng::seed_from_u64(cli.seed),
+                continuous_2d::ContinuousStorage { inner: kiddo::KdTree::new(), ball_radius_sq: ball_radius * ball_radius },
+                continuous_2d::ContinuousWalker { walk_step },
+                UniformSpawner,
+                continuous_2d::ContinuousSticker { range_sq: stick_distance * stick_distance },
+                cli.max_particles,
+            );
+
+            drive_system(&mut sys, cli.max_frames, cli.notify_every);
+            write(&sys, cli.format, &cli.output);
+        }
+
+        PCM::Balls(BallsCli { ball_radius, stick_distance, walk_step }) => {
+            let mut sys = DLASystem::new(
+                SmallRng::seed_from_u64(cli.seed),
+                ContinuousStorage { inner: kiddo::KdTree::new(), ball_radius_sq: ball_radius * ball_radius },
+                ContinuousWalker { walk_step },
+                UniformSpawner,
+                ContinuousSticker { range_sq: stick_distance * stick_distance },
+                cli.max_particles,
+            );
+
+            drive_system(&mut sys, cli.max_frames, cli.notify_every);
+            write(&sys, cli.format, &cli.output);
+        }
+
+        PCM::SurfaceProbabilityMeasure(SurfaceProbabilityMeasureCli { grid_size, stick_probability, particles, initial_data }) => {
+            let logger_sticker = LoggerSticker::new(stick_probability);
+
+            let mut sys = DLASystem::new(
+                SmallRng::seed_from_u64(cli.seed),
+                ReadOnlyVectorStorage::new(&initial_data, grid_size),
+                LocalRandomWalker,
+                UniformSpawner,
+                &logger_sticker,
+                cli.max_particles,
+            );
+
+            let particles = particles as usize;
+
+            while sys.running && logger_sticker.stick_positions.borrow().len() < particles {
+                sys.update();
+            }
+
+            let mut writer = csv::Writer::from_path(cli.output)
+                .expect("Failed to open output path");
+
+            let stick_positions = logger_sticker.stick_positions.borrow();
+
+            stick_positions
+                .iter()
+                .for_each(|pos|
+                    writer.serialize(pos)
+                        .expect("Failed to write position")
+                );
+
+            writer.flush()
+                .unwrap();
+        }
     }
-
-    sys.export_data(&cli.csv_path)
-        .expect("Failed to write");
 }

src/surface_probability_measure.rs

@@ -0,0 +1,70 @@
+use std::cell::RefCell;
+use std::path::{Path, PathBuf};
+use rand::rngs::SmallRng;
+use rand::{Rng, SeedableRng};
+use crate::system::model::DLASystem;
+use crate::system::{Storage};
+use crate::system::spaces::square_grid::Grid2D;
+use crate::system::spaces::VectorStorage;
+use crate::system::spawner::UniformSpawner;
+use crate::system::sticker::{ProbabilisticSticking, SimpleSticking, Sticker};
+use crate::system::walker::LocalRandomWalker;
+
+pub struct LoggerSticker {
+    inner: ProbabilisticSticking,
+    pub stick_positions: RefCell<Vec<Grid2D>>,
+}
+
+impl LoggerSticker {
+    pub fn new(stick_probability: f32) -> LoggerSticker {
+        LoggerSticker {
+            inner: ProbabilisticSticking { stick_probability },
+            stick_positions: RefCell::new(Vec::new()),
+        }
+    }
+}
+
+impl<S: Storage<Grid2D>> Sticker<Grid2D, S> for &LoggerSticker {
+    fn should_stick<R: Rng>(&self, rng: &mut R, space: &S, position: &Grid2D) -> bool {
+        let should_stick = self.inner.should_stick(rng, space, position);
+
+        if should_stick {
+            self.stick_positions.borrow_mut()
+                .push(position.clone());
+        }
+
+        // Writes are ignored so we deposit to delete the particle but it is not written
+        should_stick
+    }
+}
+
+pub struct ReadOnlyVectorStorage {
+    inner: VectorStorage,
+}
+
+impl ReadOnlyVectorStorage {
+    pub fn new(path: &Path, grid_size: u32) -> ReadOnlyVectorStorage {
+        let mut inner = VectorStorage::new(grid_size, 2);
+        let positions: Vec<Grid2D> = csv::Reader::from_path(path)
+            .expect("Failed to read initial data")
+            .deserialize::<Grid2D>()
+            .map(|row| row.expect("Failed to read row"))
+            .collect();
+
+        for pos in positions {
+            inner.deposit(&pos);
+        }
+
+        ReadOnlyVectorStorage { inner }
+    }
+}
+
+impl Storage<Grid2D> for ReadOnlyVectorStorage {
+    fn is_occupied(&self, position: &Grid2D) -> bool {
+        self.inner.is_occupied(position)
+    }
+
+    fn deposit(&mut self, position: &Grid2D) {
+        eprintln!("Write ignored for space at {position:?}");
+    }
+}

src/system/grid.rs

@@ -1,92 +0,0 @@
-use std::f32::consts::PI;
-use std::ops::Add;
-use num_integer::Integer;
-use rand::Rng;
-use crate::system::{GriddedPosition, Storage};
-use serde::{Serialize, Deserialize};
-
-pub struct VectorStorage {
-    backing: Vec<bool>,
-    grid_size: u32,
-}
-
-#[derive(Clone, Debug, PartialEq, Eq, Serialize, Deserialize)]
-pub struct Position {
-    pub x: i32,
-    pub y: i32,
-}
-
-impl Position {
-    pub fn in_direction(direction: u32, value: i32) -> Self {
-        if direction == 0 { Position { x: value, y: 0 } } else { Position { x: 0, y: value } }
-    }
-}
-
-impl Add for Position {
-    type Output = Position;
-
-    fn add(self, rhs: Self) -> Self::Output {
-        Position { x: self.x + rhs.x, y: self.y + rhs.y }
-    }
-}
-
-impl GriddedPosition for Position {
-    const NEIGHBOURS: u32 = 4;
-
-    fn zero() -> Position {
-        Position { x: 0, y: 0 }
-    }
-
-    fn spawn<R: Rng>(rng: &mut R, radius: f32) -> Self {
-        let theta = rng.gen_range(0f32..1.0) * 2.0 * PI;
-        let (x, y) = (radius * theta.cos(), radius * theta.sin());
-        Position { x: x as i32, y: y as i32 }
-    }
-
-    fn abs(&self) -> f32 {
-        ((self.x.pow(2) + self.y.pow(2)) as f32).powf(0.5)
-    }
-
-    fn neighbour(&self, neighbour_index: u32) -> Self {
-        let (dim, sign) = neighbour_index.div_rem(&2);
-        let sign = if sign == 0 { 1 } else { -1 };
-        let offset = Position::in_direction(dim, sign);
-
-        self.clone() + offset
-    }
-
-    fn linear_index(&self, grid_size: u32) -> usize {
-        let grid_size = grid_size as i32;
-
-        assert!(self.x <= grid_size && -(grid_size) <= self.x);
-        assert!(self.y <= grid_size && -(grid_size) <= self.y);
-
-        let x = (self.x + (grid_size) / 2) as usize;
-        let y = (self.y + (grid_size) / 2) as usize;
-
-        return grid_size as usize * y + x;
-    }
-}
-
-impl VectorStorage {
-    pub fn new(grid_size: u32) -> VectorStorage {
-        VectorStorage { grid_size, backing: vec![false; grid_size.pow(2) as usize] }
-    }
-
-    /*
-     * Convenience function for c-binding
-     * */
-    pub fn write(&mut self, position: &Position, val: bool) {
-        self.backing[position.linear_index(self.grid_size)] = val;
-    }
-}
-
-impl Storage<Position> for VectorStorage {
-    fn at(&self, position: &Position) -> bool {
-        return self.backing[position.linear_index(self.grid_size)]
-    }
-
-    fn deposit(&mut self, position: &Position) {
-        self.backing[position.linear_index(self.grid_size)] = true;
-    }
-}

src/system/mod.rs

@@ -1,24 +1,30 @@
 use std::ops::Add;
-use rand::Rng;
 use serde::Serialize;
 
 pub mod walker;
-pub mod grid;
+pub mod spawner;
+pub mod sticker;
 pub mod model;
-pub mod nd;
-mod hexagonal;
+pub mod spaces;
 
-pub trait GriddedPosition: Add<Output=Self> + Serialize + Clone {
-    const NEIGHBOURS: u32;
+pub trait Position: Add<Output=Self> + Serialize + Clone {
+    const DIM: usize;
+    // type Cartesian;
 
     fn zero() -> Self;
-    fn spawn<R: Rng>(rng: &mut R, radius: f32) -> Self;
     fn abs(&self) -> f32;
+    fn from_cartesian(cartesian: &[f32]) -> Self;
+    fn to_cartesian(&self) -> Vec<f32>;
+}
+
+pub trait GriddedPosition: Position {
+    const NEIGHBOURS: u32;
+
     fn neighbour(&self, neighbour_index: u32) -> Self;
     fn linear_index(&self, grid_size: u32) -> usize;
 }
 
-pub trait Storage<P: GriddedPosition> {
-    fn at(&self, position: &P) -> bool;
+pub trait Storage<P: Position> {
+    fn is_occupied(&self, position: &P) -> bool;
     fn deposit(&mut self, position: &P);
 }

src/system/model.rs

@@ -1,19 +1,44 @@
-use std::io;
-use std::path::Path;
 use rand::prelude::*;
-use crate::system::{GriddedPosition, Storage};
+use serde::{Deserialize, Serialize};
+use crate::system::{Position, Storage};
+use crate::system::spawner::Spawner;
+use crate::system::sticker::Sticker;
 use crate::system::walker::Walker;
 
-pub struct DLASystem<R: Rng, P: GriddedPosition, S: Storage<P>, W: Walker<P>> {
+#[derive(Serialize, Deserialize)]
+pub struct HistoryLine<P: Position> {
+    pub frame: usize,
+    pub cluster_radius: f32,
+    pub fd: f32,
+    pub position: P,
+}
+
+pub struct DLASystem<R: Rng, P: Position, S: Storage<P>, W: Walker<P>, Sp: Spawner<P>, St: Sticker<P, S>> {
     rng: R,
+
+    /*
+     * These object encapsulate the behaviour of our particular model.
+     */
+
+    /*
+     * The space, along with the chosen position choose the embedding space that the cluster grows
+     * in.
+     */
     space: S,
+
+    /*
+     * Walkers allow us to choose between different particle movement behaviours, eg random or
+     * biased walker
+     */
     walker: W,
+    spawner: Sp,
+    sticker: St,
 
-    stick_probability: f32,
-    max_particles: usize,
-
+    pub frame: usize,
     pub running: bool,
+    pub history: Vec<HistoryLine<P>>,
 
+    max_particles: usize,
     particles: Vec<P>,
     active_particle: Option<P>,
 
@@ -24,17 +49,21 @@ pub struct DLASystem<R: Rng, P: GriddedPosition, S: Storage<P>, W: Walker<P>> {
     cluster_radius: f32,
 }
 
-impl<R: Rng, P: GriddedPosition, S: Storage<P>, W: Walker<P>> DLASystem<R, P, S, W> {
-    pub fn new_g(rng: R, space: S, walker: W, stick_probability: f32, max_particles: usize) -> Self {
+impl<R: Rng, P: Position, S: Storage<P>, W: Walker<P>, Sp: Spawner<P>, St: Sticker<P, S>> DLASystem<R, P, S, W, Sp, St> {
+    pub fn new(rng: R, space: S, walker: W, spawner: Sp, sticker: St, max_particles: usize) -> Self {
         let mut sys = DLASystem {
             rng,
-            stick_probability,
             max_particles,
             running: true,
+            spawner,
+            sticker,
 
             space,
             walker,
+
+            frame: 0,
             particles: vec![],
+            history: vec![],
             active_particle: None,
 
             add_ratio: 1.2,
@@ -51,6 +80,7 @@ impl<R: Rng, P: GriddedPosition, S: Storage<P>, W: Walker<P>> DLASystem<R, P, S,
     }
 
     pub fn update(&mut self) {
+        self.frame += 1;
         if self.active_particle.is_some() {
             self.move_particle();
         } else if self.particles.len() < self.max_particles {
@@ -71,8 +101,8 @@ impl<R: Rng, P: GriddedPosition, S: Storage<P>, W: Walker<P>> DLASystem<R, P, S,
 
         if distance > self.kill_circle {
             self.active_particle = None;
-        } else if !self.space.at(&next_position) {
-            if self.check_stick(&next_position) {
+        } else if !self.space.is_occupied(&next_position) {
+            if self.sticker.should_stick(&mut self.rng, &self.space, &next_position) {
                 self.deposit(&next_position);
                 self.active_particle = None;
                 return;
@@ -82,19 +112,10 @@ impl<R: Rng, P: GriddedPosition, S: Storage<P>, W: Walker<P>> DLASystem<R, P, S,
         }
     }
 
-    fn check_stick(&mut self, position: &P) -> bool {
-        (0..P::NEIGHBOURS)
-            .map(|n| position.neighbour(n))
-            .any(|neighbour|
-                self.space.at(&neighbour)
-                    && self.rng.gen_range(0.0f32..=1.0) < self.stick_probability
-            )
-    }
-
     fn spawn_particle(&mut self) {
-        let position = P::spawn(&mut self.rng, self.add_circle);
+        let position = self.spawner.spawn(&mut self.rng, self.add_circle);
 
-        if !self.space.at(&position) {
+        if !self.space.is_occupied(&position) {
             self.active_particle = Some(position);
         }
     }
@@ -113,16 +134,7 @@ impl<R: Rng, P: GriddedPosition, S: Storage<P>, W: Walker<P>> DLASystem<R, P, S,
                 self.kill_circle = self.kill_ratio * self.add_circle;
             }
         }
-    }
 
-    pub fn export_data(&self, path: &Path) -> io::Result<()> {
-        let mut wtr = csv::Writer::from_path(path)?;
-
-        for particle in &self.particles {
-            wtr.serialize(particle)?;
-        }
-
-        wtr.flush()?;
-        Ok(())
+        self.history.push(HistoryLine { frame: self.frame, position: p0.clone(), cluster_radius: self.cluster_radius, fd: (self.particles.len() as f32).ln() / self.cluster_radius.ln() });
     }
 }

src/system/nd.rs

@@ -1,107 +0,0 @@
-use std::ops::Add;
-use num_integer::Integer;
-use rand::Rng;
-use serde::{Serialize, Serializer};
-use serde::ser::SerializeMap;
-use crate::system::GriddedPosition;
-use crate::system::Storage;
-
-pub struct NDVectorStorage<const DIM: usize> {
-    backing: Vec<bool>,
-    grid_size: u32,
-}
-
-impl<const DIM: usize> NDVectorStorage<DIM> {
-    pub fn new(grid_size: u32) -> Self {
-        Self { grid_size, backing: vec![false; grid_size.pow(DIM as u32) as usize] }
-    }
-}
-
-impl<const DIM: usize> Storage<NDPosition<DIM>> for NDVectorStorage<DIM> {
-    fn at(&self, position: &NDPosition<DIM>) -> bool {
-        return self.backing[position.linear_index(self.grid_size)];
-    }
-
-    fn deposit(&mut self, position: &NDPosition<DIM>) {
-        self.backing[position.linear_index(self.grid_size)] = true;
-    }
-}
-
-#[derive(Clone, Debug, PartialEq, Eq)]
-pub struct NDPosition<const DIM: usize>([i32; DIM]);
-
-impl<const DIM: usize> NDPosition<DIM> {
-    pub fn in_direction(direction: usize, value: i32) -> Self {
-        let mut arr = [0; DIM];
-        arr[direction] = value;
-        NDPosition(arr)
-    }
-}
-
-impl<const DIM: usize> Add for NDPosition<DIM> {
-    type Output = Self;
-
-    fn add(self, rhs: Self) -> Self::Output {
-        Self(self.0.zip(rhs.0).map(|(a, b)| a + b))
-    }
-}
-
-impl<const DIM: usize> Serialize for NDPosition<DIM> {
-    fn serialize<S>(&self, serializer: S) -> Result<S::Ok, S::Error> where S: Serializer {
-        let mut map = serializer.serialize_map(Some(DIM))?;
-        for (i, v) in self.0.iter().enumerate() {
-            map.serialize_entry(&format!("r{}", i), v)?;
-        }
-        map.end()
-    }
-}
-
-impl<const DIM: usize> GriddedPosition for NDPosition<DIM> {
-    const NEIGHBOURS: u32 = { 2u32.pow(DIM as u32) } as u32;
-
-    fn zero() -> Self {
-        NDPosition([0; DIM])
-    }
-
-    fn spawn<R: Rng>(rng: &mut R, radius: f32) -> Self {
-        let mut a: [f32; DIM] = [0f32; DIM];
-        let mut b: [i32; DIM] = [0i32; DIM];
-
-        for i in 0..DIM {
-            a[i] = rng.gen_range(0f32..1f32);
-        }
-
-        let norm = a.iter().sum::<f32>()
-            .sqrt();
-
-        for i in 0..DIM {
-            a[i] = a[i] * radius / norm;
-            b[i] = a[i] as i32;
-        }
-
-        return Self(b);
-    }
-
-    fn abs(&self) -> f32 {
-        let a: i32 = self.0.iter()
-            .map(|r| r.pow(2))
-            .sum();
-
-        (a as f32).powf(0.5)
-    }
-
-    fn neighbour(&self, neighbour_index: u32) -> Self {
-        let (dim, sign) = neighbour_index.div_rem(&(DIM as u32));
-        let sign = if sign == 0 { 1 } else { -1 };
-        let offset = Self::in_direction(dim as usize, sign);
-
-        self.clone() + offset
-    }
-
-    fn linear_index(&self, grid_size: u32) -> usize {
-        self.0.iter()
-            .enumerate()
-            .map(|(i, v)| (grid_size.pow(i as u32) as usize) * (v + ((grid_size / 2) as i32)) as usize)
-            .sum()
-    }
-}

src/system/spaces/continuous_2d.rs

@@ -0,0 +1,102 @@
+use std::f32::consts::PI;
+use std::ops::Add;
+use kiddo::distance::squared_euclidean;
+use rand::Rng;
+use serde::{Serialize, Deserialize};
+use crate::system::sticker::Sticker;
+use crate::system::{Position, Storage};
+use crate::system::walker::Walker;
+
+#[derive(Serialize, Deserialize, Debug, Clone)]
+pub struct P2 {
+    pub x: f32,
+    pub y: f32,
+}
+
+impl P2 {
+    fn as_arr(&self) -> [f32; 2] {
+        [self.x, self.y]
+    }
+
+    pub fn random_with_radius<R: Rng>(rng: &mut R, radius: f32) -> P2 {
+        let theta = rng.gen_range(0f32..1.0) * 2.0 * PI;
+
+        let (x, y) = (
+            radius * theta.sin(),
+            radius * theta.cos(),
+        );
+
+        P2 { x, y }
+    }
+}
+
+impl Add for P2 {
+    type Output = P2;
+
+    fn add(self, rhs: Self) -> Self::Output {
+        P2 {
+            x: self.x + rhs.x,
+            y: self.y + rhs.y,
+        }
+    }
+}
+
+pub struct ContinuousStorage {
+    pub inner: kiddo::KdTree<f32, (), 2>,
+    pub ball_radius_sq: f32,
+}
+
+impl Position for P2 {
+    const DIM: usize = 2;
+
+    fn zero() -> Self {
+        P2 { x: 0f32, y: 0f32 }
+    }
+
+    fn abs(&self) -> f32 {
+        (self.x.powi(2) + self.y.powi(2)).powf(0.5)
+    }
+
+    fn from_cartesian(cartesian: &[f32]) -> Self {
+        P2 { x: cartesian[0], y: cartesian[1] }
+    }
+
+    fn to_cartesian(&self) -> Vec<f32> {
+        vec![self.x, self.y]
+    }
+}
+
+impl Storage<P2> for ContinuousStorage {
+    fn is_occupied(&self, position: &P2) -> bool {
+        let (dist_sq, _) = self.inner.nearest_one(&position.as_arr(), &squared_euclidean).unwrap();
+
+        // Is the distance of this point to the next one less than twice the ball radius
+        dist_sq < 2.0 * self.ball_radius_sq
+    }
+
+    fn deposit(&mut self, position: &P2) {
+        self.inner.add(&position.as_arr(), ()).expect("Failed to write to space")
+    }
+}
+
+pub struct ContinuousSticker {
+    /// INVARIANT: THIS SHOULD BE GREATER THAN THE BALL_RADIUS_SQ value
+    pub range_sq: f32,
+}
+
+impl Sticker<P2, ContinuousStorage> for ContinuousSticker {
+    fn should_stick<R: Rng>(&self, _rng: &mut R, space: &ContinuousStorage, position: &P2) -> bool {
+        let (a, _) = space.inner.nearest_one(&position.as_arr(), &squared_euclidean).unwrap();
+        a < self.range_sq
+    }
+}
+
+pub struct ContinuousWalker {
+    pub walk_step: f32
+}
+
+impl Walker<P2> for ContinuousWalker {
+    fn walk<R: Rng>(&self, rng: &mut R, position: &P2) -> P2 {
+        position.clone() + P2::random_with_radius(rng, self.walk_step)
+    }
+}

src/system/spaces/continuous_3d.rs

@@ -0,0 +1,106 @@
+use std::f32::consts::PI;
+use std::ops::Add;
+use kiddo::distance::squared_euclidean;
+use rand::Rng;
+use serde::Serialize;
+use crate::system::sticker::Sticker;
+use crate::system::{Position, Storage};
+use crate::system::walker::Walker;
+
+#[derive(Serialize, Debug, Clone)]
+pub struct P3 {
+    x: f32,
+    y: f32,
+    z: f32,
+}
+
+impl P3 {
+    fn as_arr(&self) -> [f32; 3] {
+        [self.x, self.y, self.z]
+    }
+
+    pub fn random_with_radius<R: Rng>(rng: &mut R, radius: f32) -> P3 {
+        let theta = rng.gen_range(0f32..1.0) * 2.0 * PI;
+        let phi = rng.gen_range(0f32..1.0) * 2.0 * PI;
+
+        let (x, y, z) = (
+            radius * theta.sin() * phi.cos(),
+            radius * theta.sin() * phi.sin(),
+            radius * theta.cos()
+        );
+
+        P3 { x, y, z}
+    }
+}
+
+impl Add for P3 {
+    type Output = P3;
+
+    fn add(self, rhs: Self) -> Self::Output {
+        P3 {
+            x: self.x + rhs.x,
+            y: self.y + rhs.y,
+            z: self.z + rhs.z,
+        }
+    }
+}
+
+pub struct ContinuousStorage {
+    pub inner: kiddo::KdTree<f32, (), 3>,
+    pub ball_radius_sq: f32,
+}
+
+impl Position for P3 {
+    const DIM: usize = 3;
+
+    fn zero() -> Self {
+        P3 { x: 0f32, y: 0f32, z: 0f32 }
+    }
+
+    fn abs(&self) -> f32 {
+        (self.x.powi(2) + self.y.powi(2) + self.z.powi(2)).powf(0.5)
+    }
+
+    fn from_cartesian(cartesian: &[f32]) -> Self {
+        P3 { x: cartesian[0], y: cartesian[1], z: cartesian[3] }
+    }
+
+    fn to_cartesian(&self) -> Vec<f32> {
+        vec![self.x, self.y, self.z]
+    }
+}
+
+impl Storage<P3> for ContinuousStorage {
+    fn is_occupied(&self, position: &P3) -> bool {
+        let (dist_sq, _) = self.inner.nearest_one(&position.as_arr(), &squared_euclidean).unwrap();
+
+        // Is the distance of this point to the next one less than twice the ball radius
+        dist_sq < 2.0 * self.ball_radius_sq
+    }
+
+    fn deposit(&mut self, position: &P3) {
+        self.inner.add(&position.as_arr(), ()).expect("Failed to write to space")
+    }
+}
+
+pub struct ContinuousSticker {
+    /// INVARIANT: THIS SHOULD BE GREATER THAN THE BALL_RADIUS_SQ value
+    pub range_sq: f32,
+}
+
+impl Sticker<P3, ContinuousStorage> for ContinuousSticker {
+    fn should_stick<R: Rng>(&self, _rng: &mut R, space: &ContinuousStorage, position: &P3) -> bool {
+        let (a, _) = space.inner.nearest_one(&position.as_arr(), &squared_euclidean).unwrap();
+        a < self.range_sq
+    }
+}
+
+pub struct ContinuousWalker {
+    pub walk_step: f32
+}
+
+impl Walker<P3> for ContinuousWalker {
+    fn walk<R: Rng>(&self, rng: &mut R, position: &P3) -> P3 {
+        position.clone() + P3::random_with_radius(rng, self.walk_step)
+    }
+}

src/system/spaces/hexagonal.rs

@@ -1,7 +1,7 @@
 use std::ops::Add;
-use bevy::render::color::HexColorError::Hex;
-use rand::Rng;
-use crate::system::GriddedPosition;
+use num_integer::Roots;
+use num_traits::Pow;
+use crate::system::{GriddedPosition, Position};
 use serde::{Serialize, Deserialize};
 
 #[derive(Clone, Debug, PartialEq, Eq, Serialize, Deserialize)]
@@ -21,18 +21,6 @@ impl Add for HexPosition {
 impl GriddedPosition for HexPosition {
     const NEIGHBOURS: u32 = 6;
 
-    fn zero() -> Self {
-        HexPosition { q: 0, r: 0 }
-    }
-
-    fn spawn<R: Rng>(rng: &mut R, radius: f32) -> Self {
-        todo!()
-    }
-
-    fn abs(&self) -> f32 {
-        ((self.q.pow(2) + self.r.pow(2) + self.q * self.r) as f32).sqrt()
-    }
-
     fn neighbour(&self, neighbour_index: u32) -> Self {
         let neighbour_index = neighbour_index as usize;
 
@@ -41,11 +29,42 @@ impl GriddedPosition for HexPosition {
             (-1, 0), (-1, 1), (0, 1),
         ];
 
-        self.clone() + HexPosition { q: OFFSETS[neighbour_index].0, r: OFFSETS[neighbour_index].0 }
+        self.clone() + HexPosition { q: OFFSETS[neighbour_index].0, r: OFFSETS[neighbour_index].1 }
     }
 
     fn linear_index(&self, grid_size: u32) -> usize {
-        todo!()
-        // ((self.q + grid_size / 2) + grid_size * self.r) as usize
+        let q = (self.q + (grid_size as i32 / 2)) as usize;
+        let r = (self.r + (grid_size as i32 / 2)) as usize;
+
+        r * (grid_size as usize) + q
+    }
+}
+
+impl Position for HexPosition {
+    const DIM: usize = 2;
+
+    fn zero() -> Self {
+        HexPosition { q: 0, r: 0 }
+    }
+
+    fn abs(&self) -> f32 {
+        ((self.q.pow(2) + self.r.pow(2) + self.q * self.r) as f32).sqrt()
+    }
+
+    fn from_cartesian(cartesian: &[f32]) -> Self {
+        let q = (1.0f32 / 3.0f32).sqrt() * cartesian[0] - 1.0 / 3.0 * cartesian[1];
+        let r = 2.0 / 3.0 * cartesian[1];
+
+        Self { q: q as i32, r: r as i32 }
+    }
+
+    fn to_cartesian(&self) -> Vec<f32> {
+        let q = self.q as f32;
+        let r = self.r as f32;
+
+        vec![
+            3f32.sqrt() * q + 3f32.sqrt() / 2f32 * r,
+            (3. / 2.) * r
+        ]
     }
 }

src/system/spaces/kd_grid.rs

@@ -0,0 +1,87 @@
+use anyhow::anyhow;
+use rand::Rng;
+use crate::system::{GriddedPosition, Position, Storage};
+use crate::system::spaces::square_grid::{Grid2D, Grid3D};
+use crate::system::sticker::Sticker;
+
+fn taxicab_grid2(a: &[f32; 2], b: &[f32; 2]) -> f32 {
+    (a[0] - b[0]).abs() + (a[1] - b[1]).abs()
+}
+
+fn taxicab_grid3(a: &[f32; 3], b: &[f32; 3]) -> f32 {
+    (a[0] - b[0]).abs() + (a[1] - b[1]).abs() + (a[2] - b[2]).abs()
+}
+
+pub struct KDSpace<const N: usize> {
+    pub(crate) inner: kiddo::KdTree<f32, (), N>,
+}
+
+impl<const N: usize> KDSpace<N> {
+    pub fn new() -> KDSpace<N> {
+        KDSpace { inner: kiddo::KdTree::new() }
+    }
+}
+
+impl Storage<Grid2D> for KDSpace<2> {
+    fn is_occupied(&self, position: &Grid2D) -> bool {
+        let a = self.inner.best_n_within(&[position.x as f32, position.y as f32], 0f32, 1, &taxicab_grid2).unwrap();
+        !a.is_empty()
+    }
+
+    fn deposit(&mut self, position: &Grid2D) {
+        self.inner.add(&[position.x as f32, position.y as f32], ())
+            .expect("Failed to write to space")
+    }
+}
+
+impl Storage<Grid3D> for KDSpace<3> {
+    fn is_occupied(&self, position: &Grid3D) -> bool {
+        let a = self.inner.best_n_within(&[position.x as f32, position.y as f32, position.z as f32], 0f32, 1, &taxicab_grid3).unwrap();
+        !a.is_empty()
+    }
+
+    fn deposit(&mut self, position: &Grid3D) {
+        self.inner.add(&[position.x as f32, position.y as f32, position.z as f32], ())
+            .expect("Failed to write to space")
+    }
+}
+
+pub struct KDProbabilisticSticking {
+    pub(crate) stick_probability: f32,
+}
+
+impl KDProbabilisticSticking {
+    fn new(stick_probability: f32) -> anyhow::Result<KDProbabilisticSticking> {
+        return if 0f32 < stick_probability && stick_probability <= 1f32 {
+            Ok(KDProbabilisticSticking { stick_probability })
+        } else {
+            Err(anyhow!("Sticking probability outside of (0, 1] range."))
+        }
+    }
+}
+
+impl Sticker<Grid2D, KDSpace<2>> for KDProbabilisticSticking {
+    fn should_stick<R: Rng>(&self, rng: &mut R, space: &KDSpace<2>, position: &Grid2D) -> bool {
+        let a = space.inner.best_n_within(&[position.x as f32, position.y as f32], 1f32, Grid2D::NEIGHBOURS as usize, &taxicab_grid2)
+            .unwrap();
+
+        if a.len() == 0 {
+            return false;
+        }
+
+        let q = 1f32 - self.stick_probability;
+        let a = q.powi(a.len() as i32);
+        rng.gen_range(0f32..1f32) > a
+    }
+}
+
+impl Sticker<Grid3D, KDSpace<3>> for KDProbabilisticSticking {
+    fn should_stick<R: Rng>(&self, rng: &mut R, space: &KDSpace<3>, position: &Grid3D) -> bool {
+        let a = space.inner.best_n_within(&[position.x as f32, position.y as f32, position.z as f32], 1f32, Grid2D::NEIGHBOURS as usize, &taxicab_grid3)
+            .unwrap();
+
+        let q = (1f32 - self.stick_probability);
+        let a = q.powi(a.len() as i32);
+        rng.gen_range(0f32..1f32) > a
+    }
+}

src/system/spaces/mod.rs

@@ -0,0 +1,10 @@
+pub mod vector_storage;
+pub use vector_storage::VectorStorage;
+
+pub mod square_grid;
+pub mod kd_grid;
+pub mod hexagonal;
+
+pub mod continuous_3d;
+pub mod continuous_2d;
+pub mod nalg;

src/system/spaces/nalg.rs

@@ -0,0 +1,148 @@
+use std::ops::Add;
+use itertools::Itertools;
+use nalgebra::{EuclideanNorm, LpNorm, Matrix, Norm, OMatrix, SVector};
+use num_traits::Pow;
+use serde::{Serialize, Deserialize};
+use crate::system::{GriddedPosition, Position, Storage};
+
+#[derive(Clone, PartialEq, Eq, Serialize, Deserialize)]
+#[serde(transparent)]
+pub struct Gridded<const D: usize>(SVector<i32, D>);
+
+#[derive(Clone, PartialEq, Serialize, Deserialize)]
+#[serde(transparent)]
+pub struct Continuous<const D: usize>(SVector<f32, D>);
+
+impl<const D: usize> Add for Continuous<D> {
+    type Output = Continuous<D>;
+
+    fn add(self, rhs: Self) -> Self::Output {
+        Continuous(self.0 + rhs.0)
+    }
+}
+
+impl<const D: usize> Position for Continuous<D> {
+    const DIM: usize = 0;
+
+    fn zero() -> Self {
+        Continuous(SVector::<f32, D>::zeros())
+    }
+
+    fn abs(&self) -> f32 {
+        self.0.norm()
+    }
+
+    fn from_cartesian(cartesian: &[f32]) -> Self {
+        Continuous(SVector::<f32, D>::from_fn(|i, _| cartesian[i]))
+    }
+
+    fn to_cartesian(&self) -> Vec<f32> {
+        self.0.as_slice().to_vec()
+    }
+}
+
+impl<const D: usize> Add for Gridded<D> {
+    type Output = Gridded<D>;
+
+    fn add(self, rhs: Self) -> Self::Output {
+        Gridded(self.0 + rhs.0)
+    }
+}
+
+impl<const D: usize> Position for Gridded<D> {
+    const DIM: usize = 0;
+
+    fn zero() -> Self {
+        Gridded(SVector::<i32, D>::zeros())
+    }
+
+    fn abs(&self) -> f32 {
+        (self.0.fold(0, |r, c| r + c.pow(2)) as f32).sqrt()
+    }
+
+    fn from_cartesian(cartesian: &[f32]) -> Self {
+        Gridded(SVector::<i32, D>::from_fn(|i, _| cartesian[i] as i32))
+    }
+
+    fn to_cartesian(&self) -> Vec<f32> {
+        self.0.as_slice()
+            .iter()
+            .map(|a| *a as f32)
+            .collect_vec()
+    }
+}
+
+pub struct KDSpace<const N: usize>(pub(crate) kiddo::KdTree<f32, (), N>);
+
+impl<const D: usize> Storage<Gridded<D>> for KDSpace<D> {
+    fn is_occupied(&self, position: &Gridded<D>) -> bool {
+        let a = self.0.best_n_within(
+            &position.0.data.0[0].map(|i| i as f32),
+            0f32,
+            1,
+            &|a, b| {
+                LpNorm(1).metric_distance(
+                    &SVector::<f32, D>::from_row_slice(a),
+                    &SVector::<f32, D>::from_row_slice(b),
+                )
+            },
+        ).unwrap();
+
+        !a.is_empty()
+    }
+
+    fn deposit(&mut self, position: &Gridded<D>) {
+        self.0.add(&position.0.data.0[0].map(|i| i as f32), ())
+            .expect("Failed to write to space")
+    }
+}
+
+impl<const D: usize> Storage<Continuous<D>> for KDSpace<D> {
+    fn is_occupied(&self, position: &Continuous<D>) -> bool {
+        let a = self.0.best_n_within(
+            &position.0.data.0[0],
+            0f32,
+            1,
+            &|a, b| {
+                EuclideanNorm.metric_distance(
+                    &SVector::<f32, D>::from_row_slice(a),
+                    &SVector::<f32, D>::from_row_slice(b),
+                )
+            },
+        ).unwrap();
+
+        !a.is_empty()
+    }
+
+    fn deposit(&mut self, position: &Continuous<D>) {
+        self.0.add(&position.0.data.0[0], ())
+            .expect("Failed to write to space")
+    }
+}
+
+pub struct VectorStorage {
+    backing: Vec<bool>,
+    grid_size: usize,
+}
+
+impl<const D: usize> Storage<Gridded<D>> for VectorStorage {
+    fn is_occupied(&self, position: &Gridded<D>) -> bool {
+        let mut index: usize = 0;
+
+        for i in 0..D {
+            index += (position.0[i] + (self.grid_size as i32 / 2)) as usize * self.grid_size * i;
+        }
+
+        return self.backing[index];
+    }
+
+    fn deposit(&mut self, position: &Gridded<D>) {
+        let mut index: usize = 0;
+
+        for i in 0..D {
+            index += (position.0[i] + (self.grid_size as i32 / 2)) as usize * self.grid_size * i;
+        }
+
+        self.backing[index] = true;
+    }
+}

src/system/spaces/square_grid.rs

@@ -0,0 +1,199 @@
+use std::ops::Add;
+use num_integer::Integer;
+use crate::system::{GriddedPosition, Position};
+use serde::{Serialize, Deserialize};
+
+#[derive(Clone, Debug, Hash, PartialEq, Eq, Serialize, Deserialize)]
+pub struct Grid2D {
+    pub x: i32,
+    pub y: i32,
+}
+
+impl Grid2D {
+    pub fn in_direction(direction: u32, value: i32) -> Self {
+        match direction {
+            0 => Grid2D { x: value, y: 0 },
+            1 => Grid2D { x: 0, y: value },
+            _ => panic!("Invalid direction"),
+        }
+    }
+}
+
+impl Add for Grid2D {
+    type Output = Grid2D;
+
+    fn add(self, rhs: Self) -> Self::Output {
+        Grid2D { x: self.x + rhs.x, y: self.y + rhs.y }
+    }
+}
+
+impl Position for Grid2D {
+    const DIM: usize = 2;
+
+    fn zero() -> Self {
+        Grid2D { x: 0, y: 0 }
+    }
+
+    fn abs(&self) -> f32 {
+        (((self.x * self.x) + (self.y * self.y)) as f32).powf(0.5)
+    }
+
+    fn from_cartesian(cartesian: &[f32]) -> Self {
+        Grid2D {
+            x: cartesian[0] as i32,
+            y: cartesian[1] as i32,
+        }
+    }
+
+    fn to_cartesian(&self) -> Vec<f32> {
+        vec![self.x as f32, self.y as f32]
+    }
+}
+
+impl GriddedPosition for Grid2D {
+    const NEIGHBOURS: u32 = 4;
+
+    fn neighbour(&self, neighbour_index: u32) -> Self {
+        let (dim, sign) = neighbour_index.div_rem(&2);
+        let sign = if sign == 0 { 1 } else { -1 };
+        let offset = Self::in_direction(dim, sign);
+
+        self.clone() + offset
+    }
+
+    fn linear_index(&self, grid_size: u32) -> usize {
+        let grid_size = grid_size as i32;
+
+        assert!(self.x < grid_size && -(grid_size) < self.x);
+        assert!(self.y < grid_size && -(grid_size) < self.y);
+
+        let x = (self.x + (grid_size) / 2) as usize;
+        let y = (self.y + (grid_size) / 2) as usize;
+
+        let linear_index = grid_size as usize * y + x;
+
+        if linear_index >= (grid_size * grid_size) as usize {
+            eprintln!("AHHH SOMETHING WENT WRONG {:?} gives {}", self, linear_index);
+        }
+
+        return linear_index;
+    }
+}
+
+#[derive(Clone, Debug, PartialEq, Eq, Hash, Serialize, Deserialize)]
+pub struct Grid3D {
+    pub x: i32,
+    pub y: i32,
+    pub z: i32,
+}
+
+impl Grid3D {
+    pub fn in_direction(direction: u32, value: i32) -> Self {
+        match direction {
+            0 => Grid3D { x: value, y: 0, z: 0 },
+            1 => Grid3D { x: 0, y: value, z: 0 },
+            2 => Grid3D { x: 0, y: 0, z: value },
+            _ => panic!("Invalid direction"),
+        }
+    }
+}
+
+impl Add for Grid3D {
+    type Output = Grid3D;
+
+    fn add(self, rhs: Self) -> Self::Output {
+        Grid3D {
+            x: self.x + rhs.x,
+            y: self.y + rhs.y,
+            z: self.z + rhs.z,
+        }
+    }
+}
+
+#[test]
+fn grid3_add_test() {
+    assert_eq!(
+        Grid3D::from_cartesian(&[0f32, 0f32, 0f32]) + Grid3D::from_cartesian(&[0f32, 1f32, 0f32]),
+        Grid3D::from_cartesian(&[0f32, 1f32, 0f32])
+    );
+    assert_eq!(
+        Grid3D::from_cartesian(&[5.0, 3.0, 1.0]) + Grid3D::from_cartesian(&[-2.0, 5.0, 100.0]),
+        Grid3D::from_cartesian(&[3.0, 8.0, 101.0])
+    );
+}
+
+impl Position for Grid3D {
+    const DIM: usize = 3;
+
+    fn zero() -> Self {
+        Grid3D { x: 0, y: 0, z: 0 }
+    }
+
+    fn abs(&self) -> f32 {
+        (((self.x * self.x) + (self.y * self.y) + (self.z * self.z)) as f32).powf(0.5)
+    }
+
+    fn from_cartesian(cartesian: &[f32]) -> Self {
+        Self {
+            x: cartesian[0] as i32,
+            y: cartesian[1] as i32,
+            z: cartesian[2] as i32,
+        }
+    }
+
+    fn to_cartesian(&self) -> Vec<f32> {
+        vec![self.x as f32, self.y as f32, self.z as f32]
+    }
+}
+
+#[test]
+fn grid3_neighbours_test() {
+    let neighbours = (0..Grid3D::NEIGHBOURS)
+        .map(|n| Grid3D::neighbour(&Grid3D::zero(), n))
+        .collect::<Vec<_>>();
+
+    assert!(neighbours.contains(&Grid3D { x: 1, y: 0, z: 0 }));
+    assert!(neighbours.contains(&Grid3D { x: -1, y: 0, z: 0 }));
+    assert!(neighbours.contains(&Grid3D { x: 0, y: 1, z: 0 }));
+    assert!(neighbours.contains(&Grid3D { x: 0, y: -1, z: 0 }));
+    assert!(neighbours.contains(&Grid3D { x: 0, y: 0, z: 1 }));
+    assert!(neighbours.contains(&Grid3D { x: 0, y: 0, z: -1 }));
+}
+
+impl GriddedPosition for Grid3D {
+    const NEIGHBOURS: u32 = 6;
+
+    fn neighbour(&self, neighbour_index: u32) -> Self {
+        let (dim, sign) = neighbour_index.div_rem(&2);
+        let sign = if sign == 0 { 1 } else { -1 };
+        let offset = Self::in_direction(dim, sign);
+
+        self.clone() + offset
+    }
+
+    fn linear_index(&self, grid_size: u32) -> usize {
+        let grid_size = grid_size as i32;
+
+        assert!(self.x < grid_size && -(grid_size) < self.x);
+        assert!(self.y < grid_size && -(grid_size) < self.y);
+        assert!(self.z < grid_size && -(grid_size) < self.z);
+
+        let half_grid = (grid_size) / 2;
+        let x = (self.x + half_grid) as usize;
+        let y = (self.y + half_grid) as usize;
+        let z = (self.z + half_grid) as usize;
+
+        let grid_size_usize = grid_size as usize;
+
+        let linear_index =
+            (grid_size_usize * grid_size_usize) * z
+                + (grid_size_usize) * y
+                + x;
+
+        if linear_index >= (grid_size_usize * grid_size_usize * grid_size_usize) as usize {
+            eprintln!("AHHH SOMETHING WENT WRONG {:?} gives {}", self, linear_index);
+        }
+
+        return linear_index;
+    }
+}

src/system/spaces/vector_storage.rs

@@ -0,0 +1,24 @@
+use crate::system::{GriddedPosition, Storage};
+
+pub struct VectorStorage {
+    backing: Vec<bool>,
+    grid_size: u32,
+    dimension: u32,
+}
+
+impl VectorStorage {
+    pub fn new(grid_size: u32, dimension: u32) -> VectorStorage {
+        VectorStorage { grid_size, dimension, backing: vec![false; (grid_size as usize).pow(dimension) as usize] }
+    }
+}
+
+impl<P: GriddedPosition> Storage<P> for VectorStorage {
+    fn is_occupied(&self, position: &P) -> bool {
+        let i = position.linear_index(self.grid_size);
+        return self.backing[i];
+    }
+
+    fn deposit(&mut self, position: &P) {
+        self.backing[position.linear_index(self.grid_size)] = true;
+    }
+}

src/system/spawner.rs

@@ -0,0 +1,58 @@
+use std::f32::consts::PI;
+use rand::Rng;
+use crate::system::Position;
+use crate::system::spaces::continuous_2d::P2;
+use crate::system::spaces::continuous_3d::P3;
+use crate::system::spaces::hexagonal::HexPosition;
+use crate::system::spaces::square_grid::{Grid2D, Grid3D};
+
+pub trait Spawner<P: Position> {
+    fn spawn<R: Rng>(&self, rng: &mut R, radius: f32) -> P;
+}
+
+pub struct UniformSpawner;
+
+impl Spawner<Grid2D> for UniformSpawner {
+    fn spawn<R: Rng>(&self, rng: &mut R, radius: f32) -> Grid2D {
+        let theta = rng.gen_range(0f32..1.0) * 2.0 * PI;
+        let (x, y) = (radius * theta.cos(), radius * theta.sin());
+
+        Grid2D::from_cartesian(&[x, y])
+    }
+}
+
+impl Spawner<HexPosition> for UniformSpawner {
+    fn spawn<R: Rng>(&self, rng: &mut R, radius: f32) -> HexPosition {
+        let theta = rng.gen_range(0f32..1.0) * 2.0 * PI;
+        let (x, y) = (radius * theta.cos(), radius * theta.sin());
+
+        HexPosition::from_cartesian(&[x, y])
+    }
+}
+
+impl Spawner<Grid3D> for UniformSpawner {
+    fn spawn<R: Rng>(&self, rng: &mut R, radius: f32) -> Grid3D {
+        let theta = rng.gen_range(0f32..1.0) * 2.0 * PI;
+        let phi = rng.gen_range(0f32..1.0) * 2.0 * PI;
+
+        let (x, y, z) = (
+            radius * theta.sin() * phi.cos(),
+            radius * theta.sin() * phi.sin(),
+            radius * theta.cos()
+        );
+
+        Grid3D::from_cartesian(&[x, y, z])
+    }
+}
+
+impl Spawner<P3> for UniformSpawner {
+    fn spawn<R: Rng>(&self, rng: &mut R, radius: f32) -> P3 {
+        P3::random_with_radius(rng, radius)
+    }
+}
+
+impl Spawner<P2> for UniformSpawner {
+    fn spawn<R: Rng>(&self, rng: &mut R, radius: f32) -> P2 {
+        P2::random_with_radius(rng, radius)
+    }
+}

src/system/sticker.rs

@@ -0,0 +1,39 @@
+use anyhow::anyhow;
+use rand::Rng;
+use crate::system::{GriddedPosition, Position, Storage};
+
+pub trait Sticker<P: Position, S: Storage<P>> {
+    fn should_stick<R: Rng>(&self, rng: &mut R, space: &S, position: &P) -> bool;
+}
+
+pub struct SimpleSticking;
+
+pub struct ProbabilisticSticking {
+    pub(crate) stick_probability: f32
+}
+
+impl ProbabilisticSticking {
+    pub fn new(stick_probability: f32) -> anyhow::Result<ProbabilisticSticking> {
+        return if 0f32 < stick_probability && stick_probability <= 1f32 {
+            Ok(ProbabilisticSticking { stick_probability })
+        } else {
+            Err(anyhow!("Sticking probability outside of (0, 1] range."))
+        }
+    }
+}
+
+impl<P: GriddedPosition, S: Storage<P>> Sticker<P, S> for SimpleSticking {
+    fn should_stick<R: Rng>(&self, _rng: &mut R, space: &S, position: &P) -> bool {
+        (0..P::NEIGHBOURS)
+            .map(|n| position.neighbour(n))
+            .any(|neighbour| space.is_occupied(&neighbour))
+    }
+}
+
+impl<P: GriddedPosition, S: Storage<P>> Sticker<P, S> for ProbabilisticSticking {
+    fn should_stick<R: Rng>(&self, rng: &mut R, space: &S, position: &P) -> bool {
+        (0..P::NEIGHBOURS)
+            .map(|n| position.neighbour(n))
+            .any(|neighbour| space.is_occupied(&neighbour) && rng.gen_range(0.0f32..=1.0) < self.stick_probability)
+    }
+}

src/system/walker.rs

@@ -1,57 +1,147 @@
+use std::hash::Hash;
+use itertools::Itertools;
+use rand::distributions::Slice;
 use rand::prelude::Rng;
-use crate::system::GriddedPosition;
+use crate::system::{GriddedPosition, Position};
+use crate::system::spaces::square_grid::{Grid2D, Grid3D};
 
-pub trait Walker<P: GriddedPosition> {
+pub trait Walker<P: Position> {
     fn walk<R: Rng>(&self, rng: &mut R, position: &P) -> P;
 }
 
 pub struct LocalRandomWalker;
 
-impl<Position: GriddedPosition> Walker<Position> for LocalRandomWalker {
-    fn walk<R: Rng>(&self, rng: &mut R, position: &Position) -> Position {
-        position.neighbour(rng.gen_range(0u32..Position::NEIGHBOURS))
+impl<P: GriddedPosition> Walker<P> for LocalRandomWalker {
+    fn walk<R: Rng>(&self, rng: &mut R, position: &P) -> P {
+        position.neighbour(rng.gen_range(0u32..P::NEIGHBOURS))
     }
 }
 
-mod test {
-    use rand::rngs::SmallRng;
-    use rand::{SeedableRng, thread_rng};
-    use crate::system::{GriddedPosition, grid::Position};
-    use crate::system::walker::{LocalRandomWalker, Walker};
+pub struct DiagonalRandomWalker;
 
-    #[test]
-    fn uniformity() {
-        let walker = LocalRandomWalker;
-        let mut rng = SmallRng::from_rng(thread_rng()).unwrap();
-        let mut results: Vec<Position> = vec![];
+impl Walker<Grid3D> for DiagonalRandomWalker {
+    fn walk<R: Rng>(&self, rng: &mut R, position: &Grid3D) -> Grid3D {
+        static OFFSETS: [Grid3D; 26] = [
+            Grid3D { x: 1, y: 0, z: 0 },
+            Grid3D { x: 1, y: 1, z: 0 },
+            Grid3D { x: 1, y: -1, z: 0 },
+            Grid3D { x: -1, y: 0, z: 0 },
+            Grid3D { x: -1, y: 1, z: 0 },
+            Grid3D { x: -1, y: -1, z: 0 },
+            Grid3D { x: 0, y: 1, z: 0 },
+            Grid3D { x: 0, y: -1, z: 0 },
+            Grid3D { x: 1, y: 0, z: -1 },
+            Grid3D { x: 1, y: 1, z: -1 },
+            Grid3D { x: 1, y: -1, z: -1 },
+            Grid3D { x: -1, y: 0, z: -1 },
+            Grid3D { x: -1, y: 1, z: -1 },
+            Grid3D { x: -1, y: -1, z: -1 },
+            Grid3D { x: 0, y: 1, z: -1 },
+            Grid3D { x: 0, y: -1, z: -1 },
+            Grid3D { x: 0, y: 0, z: -1 },
+            Grid3D { x: 1, y: 0, z: 1 },
+            Grid3D { x: 1, y: 1, z: 1 },
+            Grid3D { x: 1, y: -1, z: 1 },
+            Grid3D { x: -1, y: 0, z: 1 },
+            Grid3D { x: -1, y: 1, z: 1 },
+            Grid3D { x: -1, y: -1, z: 1 },
+            Grid3D { x: 0, y: 1, z: 1 },
+            Grid3D { x: 0, y: -1, z: 1 },
+            Grid3D { x: 0, y: 0, z: 1 },
+        ];
 
-        let origin = &Position::zero();
+        position.clone() + OFFSETS[rng.gen_range(0..OFFSETS.len())].clone()
+    }
+}
 
-        let x: u32 = (1_000_000);
-        for i in 0..x {
-            results.push(walker.walk(&mut rng, origin));
+impl Walker<Grid2D> for DiagonalRandomWalker {
+    fn walk<R: Rng>(&self, rng: &mut R, position: &Grid2D) -> Grid2D {
+        static OFFSETS: [Grid2D; 8] = [
+            Grid2D { x: 1, y: 0 },
+            Grid2D { x: 1, y: 1 },
+            Grid2D { x: 1, y: -1 },
+            Grid2D { x: -1, y: 0 },
+            Grid2D { x: -1, y: 1 },
+            Grid2D { x: -1, y: -1 },
+            Grid2D { x: 0, y: 1 },
+            Grid2D { x: 0, y: -1 },
+        ];
+
+        position.clone() + OFFSETS[rng.gen_range(0..OFFSETS.len())].clone()
+    }
+}
+
+fn test_uniformity_and_range<W: Walker<P>, P>(walker: W, expected: &[P], n: usize, tolerance: f32) where P: GriddedPosition + Hash + Eq {
+    use rand::thread_rng;
+    let mut rng = thread_rng();
+    let origin = &P::zero();
+
+    let results: Vec<P> = (0..n)
+        .map(|_| walker.walk(&mut rng, origin))
+        .collect();
+
+    let groups = results.iter()
+        .into_group_map_by(|a| (*a).clone());
+
+    assert_eq!(groups.len(), expected.len(), "Wrong number of walk positions generated");
+    assert!(results.iter().unique().all(|a| expected.contains(a)), "Contains unexpected walk position");
+
+    for group in groups.values() {
+        let proportion = group.len() as f32 / n as f32;
+        assert!((proportion - (1.0 / expected.len() as f32)).abs() < tolerance, "Failed tolerance check");
+    }
+}
+
+#[test]
+fn uniformity_direct_grid2d() {
+    test_uniformity_and_range(LocalRandomWalker, &[
+        Grid2D { x: 1, y: 0 },
+        Grid2D { x: -1, y: 0 },
+        Grid2D { x: 0, y: 1 },
+        Grid2D { x: 0, y: -1 },
+    ], 1_000_000, 0.001);
+}
+
+#[test]
+fn uniformity_direct_grid3d() {
+    test_uniformity_and_range(LocalRandomWalker, &[
+        Grid3D { x: 1, y: 0, z: 0 },
+        Grid3D { x: -1, y: 0, z: 0 },
+        Grid3D { x: 0, y: 1, z: 0 },
+        Grid3D { x: 0, y: -1, z: 0 },
+        Grid3D { x: 0, y: 0, z: 1 },
+        Grid3D { x: 0, y: 0, z: -1 },
+    ], 1_000_000, 0.001);
+}
+
+#[test]
+fn diagonal_grid2d() {
+    let mut expected = Vec::new();
+
+    for x in -1..=1 {
+        for y in -1..=1 {
+            if !(x == 0 && y == 0) {
+                expected.push(Grid2D { x,  y })
+            }
         }
-
-        let a = results
-            .iter()
-            .filter(|a| **a == Position { x: 0, y: 1 })
-            .count();
-
-        let b = results
-            .iter()
-            .filter(|a| **a == Position { x: 0, y: -1 })
-            .count();
-
-        let c = results
-            .iter()
-            .filter(|a| **a == Position { x: 1, y: 0 })
-            .count();
-
-        let d = results
-            .iter()
-            .filter(|a| **a == Position { x: -1, y: 0 })
-            .count();
-
-        println!("{} {} {} {}", a as f32 / x as f32, b as f32 / x as f32, c as f32 / x as f32, d as f32 / x as f32);
     }
+
+    test_uniformity_and_range(DiagonalRandomWalker, &expected, 1_000_000, 0.001);
+}
+
+#[test]
+fn diagonal_grid3d() {
+    let mut expected = Vec::new();
+
+    for x in -1..=1 {
+        for y in -1..=1 {
+            for z in -1..=1 {
+                if !(x == 0 && y == 0 && z == 0) {
+                    expected.push(Grid3D { x,  y, z })
+                }
+            }
+        }
+    }
+
+    test_uniformity_and_range(DiagonalRandomWalker, &expected, 1_000_000, 0.001);
 }

src/tools/analysis.rs

@@ -0,0 +1,33 @@
+use crate::DataAnalysisCli;
+
+use polars::prelude::*;
+use polars::io::RowCount;
+use rayon::prelude::*;
+use walkdir::WalkDir;
+
+pub(crate) fn main(cli: &DataAnalysisCli) {
+    let a: Vec<_> = WalkDir::new(&cli.sp_dir)
+        .into_iter()
+        .par_bridge()
+        .filter(|run_file| run_file.as_ref().unwrap().path().is_file())
+        .map(|run_file| {
+            let run_file = run_file.unwrap().path().to_path_buf();
+            let probability = run_file.parent().unwrap().file_name().unwrap().to_str().unwrap().parse::<f32>().unwrap();
+            let run_seed = run_file.file_stem().unwrap().to_str().unwrap();
+
+            LazyCsvReader::new(&run_file)
+                .has_header(true)
+                // Add N, done at this stage to only count within single CSV file
+                .with_row_count(Some(RowCount { name: String::from("N"), offset: 0 }))
+                .finish().unwrap()
+                .with_columns([
+                    lit(probability).alias("probability"),
+                    lit(run_seed).alias("run")
+                ])
+        }).collect();
+
+    let p = concat(a, true, true).unwrap();
+    let p = p.collect().unwrap();
+
+    println!("{:?}", p);
+}

src/tools/boxcount.rs

@@ -0,0 +1,134 @@
+use std::fs::File;
+use std::os::unix::fs::symlink;
+use crate::system::spaces::square_grid::{Grid2D, Grid3D};
+use itertools::{Itertools, MinMaxResult};
+use clap::Parser;
+use serde::Serialize;
+use crate::BoxCountCli;
+use crate::cli::cli::OutputFormat;
+use crate::system::{GriddedPosition, Position};
+use crate::system::model::HistoryLine;
+use crate::tools::read;
+
+fn bb(data: &Vec<Grid2D>) -> ((i32, i32), (i32, i32)) {
+    let x = data
+        .iter().minmax_by(|a, b| a.x.cmp(&b.x));
+
+    let y = data
+        .iter().minmax_by(|a, b| a.y.cmp(&b.y));
+
+    match (x, y) {
+        (MinMaxResult::MinMax(min_x, max_x), MinMaxResult::MinMax(min_y, max_y)) => {
+            ((min_x.x, min_y.y), (max_x.x, max_y.y))
+        }
+
+        _ => panic!("Cannot determine bounding box")
+    }
+}
+
+fn box_count_2d(data: &Vec<Grid2D>, box_number: u32) -> (f64, usize) {
+    let ((x_min, y_min), (x_max, y_max)) = bb(data);
+    let x_range = (x_max - x_min) as f64;
+    let y_range = (y_max - y_min) as f64;
+
+    let w: f64 = x_range / (box_number as f64);
+
+    let boxes_occupied = data.iter()
+        .map(|Grid2D { x, y }| [((x - x_min) as f64 / w) as i32, ((y - y_min) as f64 / w) as i32])
+        .unique()
+        .count();
+
+    (w, boxes_occupied)
+}
+
+fn box_count_3d(data: &Vec<Grid3D>, size: u32) -> usize {
+    let n = data.len();
+
+    let x_min = data
+        .iter()
+        .min_by(|Grid3D { x: x1, .. }, Grid3D { x: x2, .. }| x1.cmp(x2))
+        .unwrap().x;
+
+    let x_max = data
+        .iter()
+        .max_by(|Grid3D { x: x1, .. }, Grid3D { x: x2, .. }| x1.cmp(x2))
+        .unwrap().x;
+
+    let y_min = data
+        .iter()
+        .min_by(|Grid3D { y: v1, .. }, Grid3D { y: v2, .. }| v1.cmp(v2))
+        .unwrap().y;
+
+    let y_max = data
+        .iter()
+        .max_by(|Grid3D { y: v1, .. }, Grid3D { y: v2, .. }| v1.cmp(v2))
+        .unwrap().y;
+
+    let z_min = data
+        .iter()
+        .min_by(|Grid3D { z: v1, .. }, Grid3D { z: v2, .. }| v1.cmp(v2))
+        .unwrap().y;
+
+    let z_max = data
+        .iter()
+        .max_by(|Grid3D { z: v1, .. }, Grid3D { z: v2, .. }| v1.cmp(v2))
+        .unwrap().y;
+
+    let x_range = (x_max - x_min) as f64;
+    let y_range = (y_max - y_min) as f64;
+    let z_range = (z_max - z_min) as f64;
+
+    let w: f64 = x_range / (size as f64);
+
+    let grid_points = data.iter()
+        .map(|Grid3D { x, y, z }| [
+            ((x - x_min) as f64 / w) as i32,
+            ((y - y_min) as f64 / w) as i32,
+            ((z - z_min) as f64 / w) as i32,
+        ])
+        .collect::<Vec<_>>();
+
+    return grid_points.iter()
+        .unique()
+        .count();
+}
+
+fn box_count_nd<const N: usize>(data: &Vec<[f32; N]>, size: u32) -> usize {
+    let ranges = (0..N).map(|n|
+        match data.iter()
+            .minmax_by(|a, b| a[n].total_cmp(&b[n])) {
+            MinMaxResult::NoElements => panic!("No data"),
+            MinMaxResult::OneElement(_) => panic!("Needs more than one point to compute boxcount"),
+            MinMaxResult::MinMax(min, max) => [min[n], min[n]],
+        }).collect::<Vec<_>>();
+
+    let w: f32 = (ranges[0][1] - ranges[0][0]) / (size as f32);
+
+    return data.iter()
+        .map(|point| -> [i32; N] { std::array::from_fn(|n| ((point[n] - ranges[n][0]) / w) as i32) })
+        .unique()
+        .count();
+}
+
+#[derive(Serialize)]
+struct FDRow {
+    w: f64,
+    n_occupied: usize,
+}
+
+pub(crate) fn main(cli: &BoxCountCli) {
+    let particles: Vec<Grid2D> = read(&cli.path, cli.format);
+    let n_particles = dbg!(particles.len());
+    let box_side_counts = 1..500;
+
+    let mut writer = csv::Writer::from_path(&cli.output)
+        .expect("Unable to create csv");
+
+    box_side_counts
+        .map(|box_side_count| box_count_2d(&particles, box_side_count))
+        // .filter(|(w, n_occupied)| *n_occupied < n_particles) // Remove saturated values
+        .map(|(w, n_occupied)| FDRow { w, n_occupied })
+        .for_each(|row| writer.serialize(row).expect("Failed to write row"));
+
+    writer.flush().unwrap();
+}

src/tools/mod.rs

@@ -0,0 +1,31 @@
+use std::fs::File;
+use std::path::Path;
+use serde::de::DeserializeOwned;
+use crate::cli::cli::OutputFormat;
+use crate::system::model::HistoryLine;
+use crate::system::Position;
+
+pub mod boxcount;
+pub mod render;
+pub mod analysis;
+
+pub fn read<T: Position>(path: &Path, format: OutputFormat) -> Vec<T> where T: DeserializeOwned {
+    match format {
+        OutputFormat::FullDataJson => read_json(path),
+        OutputFormat::Positions => read_csv(path)
+    }
+}
+
+pub fn read_json<T: Position>(path: &Path) -> Vec<T> where T: DeserializeOwned {
+    serde_json::from_reader::<_, Vec<HistoryLine<T>>>(File::open(path).expect("Failed to open file"))
+        .expect("Failed to read json")
+        .iter()
+        .map(|l| (l.position.clone()))
+        .collect::<Vec<_>>()
+}
+
+pub fn read_csv<T: Position>(path: &Path) -> Vec<T> where T: DeserializeOwned {
+    csv::Reader::from_path(path).expect("Failed to read positions csv").deserialize::<T>()
+        .collect::<Result<Vec<T>, _>>()
+        .unwrap()
+}

src/tools/render.rs

@@ -0,0 +1,125 @@
+#![feature(generic_const_exprs)]
+#![feature(let_chains)]
+
+use std::fs::File;
+use std::path::PathBuf;
+use anyhow::Context;
+use crate::cli::cli::OutputFormat;
+use crate::system::model::HistoryLine;
+use crate::system::spaces::square_grid::Grid2D;
+use clap::Parser;
+use colorous::Color;
+use itertools::Itertools;
+use num_traits::real::Real;
+use num_traits::Signed;
+use serde::de::DeserializeOwned;
+use serde::Deserialize;
+use svg::Node;
+use svg::node::element::{Circle, Polygon, Rectangle};
+use crate::{RenderCli, Space};
+use crate::system::{GriddedPosition, Position};
+use crate::system::spaces::continuous_2d::P2;
+use crate::system::spaces::hexagonal::HexPosition;
+use crate::tools::read;
+
+#[derive(Debug, Parser)]
+struct Args {
+    format: OutputFormat,
+    path: PathBuf,
+    output: PathBuf,
+}
+
+trait ToSvg {
+    fn to_svg(&self, colour: Color) -> Box<dyn Node>;
+}
+
+impl ToSvg for Grid2D {
+    fn to_svg(&self, colour: Color) -> Box<dyn Node> {
+        Box::new(Rectangle::new()
+            .set("fill", format!("rgb({}, {}, {})", colour.r, colour.g, colour.b))
+            .set("width", 1)
+            .set("height", 1)
+            .set("x", self.x)
+            .set("y", self.y))
+    }
+}
+
+impl ToSvg for P2 {
+    fn to_svg(&self, colour: Color) -> Box<dyn Node> {
+        Box::new(Circle::new()
+            .set("fill", format!("rgb({}, {}, {})", colour.r, colour.g, colour.b))
+            .set("r", 1)
+            .set("cx", self.x)
+            .set("cy", self.y))
+    }
+}
+
+impl ToSvg for HexPosition {
+    fn to_svg(&self, colour: Color) -> Box<dyn Node> {
+        let points = [
+            [25.045, 128.0], [256.0, 0.0], [486.955, 128.0], [486.955, 384.0], [256.0, 512.0], [25.045, 384.0]
+        ];
+
+        let a = 256.0;
+        // let a = 400.0;
+        let b = points.map(|x| [
+            (x[0] / a),
+            (x[1] / a)]
+        );
+
+        let c = b.map(|p| format!("{},{}", p[0], p[1])).join(" ");
+
+        let cartesian = self.to_cartesian();
+        Box::new(Polygon::new()
+            .set("fill", format!("rgb({}, {}, {})", colour.r, colour.g, colour.b))
+            .set("points", c)
+            .set("transform", format!("translate({}, {})", cartesian[0], cartesian[1])))
+    }
+}
+
+pub fn compute_max_size<P: Position>(positions: &[P]) -> f32 {
+    let mut positions = positions.iter().map(P::to_cartesian);
+    let mut maximums: Vec<f32> = positions.next().unwrap();
+
+    for cartesian in positions {
+        for i in 0..maximums.len() {
+            maximums[i] = maximums[i].abs().max(cartesian[i].abs());
+        }
+    }
+
+    maximums.into_iter()
+        .fold(0f32, |r, c| r.abs().max(c.abs()))
+}
+
+fn render<P: Position>(args: &RenderCli) where P: DeserializeOwned + ToSvg {
+    let positions = read::<P>(&args.path, args.format);
+    let max_size = compute_max_size(&positions) + 10.0;
+
+    let mut svg = svg::Document::new()
+        .set("width", args.image_size)
+        .set("height", args.image_size)
+        .set("viewBox", format!("{} {} {} {}", -max_size, -max_size, max_size * 2.0, max_size * 2.0));
+
+    svg.append(Rectangle::new()
+        .set("fill", "white")
+        .set("width", max_size * 2.0)
+        .set("height", max_size * 2.0)
+        .set("x", -max_size)
+        .set("y", -max_size)
+    );
+
+    for (n, position) in positions.iter().enumerate() {
+        let colour = if args.colour { colorous::VIRIDIS.eval_rational(n, positions.len()) } else { Color::default() };
+        svg.append(position.to_svg(colour));
+    }
+
+    svg::write(File::create(&args.output).unwrap(), &svg).unwrap();
+}
+
+pub(crate) fn main(args: &RenderCli) {
+    match args.space {
+        Space::Grid2D => render::<Grid2D>(args),
+        Space::Continuous2D => render::<P2>(args),
+        Space::Hex => render::<HexPosition>(args),
+    }
+}

src/tools_cli.rs

@@ -0,0 +1,83 @@
+#![feature(generic_const_exprs)]
+#![feature(let_chains)]
+#![feature(array_zip)]
+
+use std::fs::File;
+use std::ops::{Index, Mul};
+use std::path::PathBuf;
+use anyhow::Context;
+use crate::cli::cli::OutputFormat;
+use crate::system::model::HistoryLine;
+use crate::system::spaces::square_grid::Grid2D;
+use clap::{Parser, Command, Args, Subcommand};
+use itertools::{concat, Itertools};
+use polars::io::RowCount;
+use serde::de::DeserializeOwned;
+use serde::Deserialize;
+use svg::Node;
+use svg::node::element::Rectangle;
+use walkdir::WalkDir;
+use crate::system::Position;
+use crate::system::spaces::hexagonal::HexPosition;
+
+mod system;
+mod cli;
+mod tools;
+
+#[derive(Debug, Parser)]
+enum ToolsCli {
+    Render(RenderCli),
+    BoxCount(BoxCountCli),
+    DataAnalysis(DataAnalysisCli),
+}
+
+#[derive(clap::ValueEnum, Clone, Debug, Copy)]
+enum Space {
+    Grid2D,
+    Continuous2D,
+    Hex,
+}
+
+#[derive(Debug, Args)]
+struct RenderCli {
+    #[arg(value_enum, short, long, default_value_t = OutputFormat::Positions)]
+    format: OutputFormat,
+
+    #[arg(value_enum, short, long, default_value_t = Space::Grid2D)]
+    space: Space,
+
+    path: PathBuf,
+    output: PathBuf,
+
+    #[arg(short, long, default_value_t = 800)]
+    image_size: u32,
+
+    #[arg(short, long, default_value_t = false)]
+    colour: bool,
+}
+
+#[derive(Debug, Args)]
+struct BoxCountCli {
+    #[arg(value_enum, short, long, default_value_t = OutputFormat::Positions)]
+    format: OutputFormat,
+    path: PathBuf,
+    output: PathBuf,
+}
+
+#[derive(Debug, Args)]
+struct DataAnalysisCli {
+    sp_dir: PathBuf,
+}
+
+fn main() -> anyhow::Result<()> {
+    let args = ToolsCli::parse();
+    dbg!(&args);
+
+    match args {
+        ToolsCli::Render(cli) => tools::render::main(&cli),
+        ToolsCli::BoxCount(cli) => tools::boxcount::main(&cli),
+        ToolsCli::DataAnalysis(cli) => tools::analysis::main(&cli),
+    }
+
+    Ok(())
+}

src/ui.rs

@@ -1,94 +0,0 @@
-#![feature(array_zip)]
-
-use std::default::Default;
-use std::error::Error;
-use bevy::{prelude::*, sprite::MaterialMesh2dBundle};
-use csv::{Position, Reader, ReaderBuilder};
-use clap::Parser;
-
-#[derive(Parser, Resource)] // requires `derive` feature
-enum UICli {
-    CSV(CSVArgs)
-}
-
-#[derive(clap::Args)]
-struct CSVArgs {
-    path: std::path::PathBuf,
-}
-
-struct Position2D(i32, i32);
-
-fn main() -> Result<(), Box<dyn Error>> {
-    let cli = UICli::parse();
-
-    App::new()
-        .insert_resource(cli)
-        .add_plugins(DefaultPlugins.set(WindowPlugin {
-            window: WindowDescriptor {
-                title: "DLA Static 2D Render".to_string(),
-                width: 800.,
-                height: 800.,
-                ..default()
-            },
-            ..default()
-        }))
-        .add_startup_system(setup_ui)
-        .add_startup_system(read_csv)
-        .run();
-
-    Ok(())
-}
-
-fn read_csv(
-    cli: Res<UICli>,
-    mut commands: Commands
-) {
-    let csv_path = match &cli.into_inner() {
-        UICli::CSV(CSVArgs { path }) => path,
-        _ => panic!("Ahh"),
-    };
-
-    let mut reader = ReaderBuilder::new()
-        .has_headers(true)
-        .from_path(csv_path)
-        .expect("Failed to read csv");
-
-    let headers = reader.headers()
-        .expect("Failed to read headers");
-
-    let x_column = headers.iter().position(|name| name.trim() == "x")
-        .expect("Failed to find x column");
-
-    let y_column = headers.iter().position(|name| name.trim() == "y")
-        .expect("Failed to find x column");
-
-    let positions = reader
-        .records()
-        .map(|record| {
-            let record = record.expect("Failed to read position");
-            let x: i32 = record[x_column].trim().parse::<i32>().expect("Failed to read x");
-            let y: i32 = record[y_column].trim().parse::<i32>().expect("Failed to read y");
-
-            Position2D(x, y)
-        });
-
-    for Position2D(x, y) in positions {
-        let rect_size = 5.0;
-
-        commands.spawn(SpriteBundle {
-            sprite: Sprite {
-                color: Color::rgb(0.25, 0.25, 0.75),
-                custom_size: Some(Vec2::new(rect_size, rect_size)),
-                ..default()
-            },
-            transform: Transform::from_translation(Vec3::new((x as f32) * rect_size, (y as f32) * rect_size, 0.)),
-            ..default()
-        });
-    }
-}
-
-fn setup_ui(
-    mut commands: Commands,
-) {
-    commands.spawn(Camera2dBundle::default());
-}