.drone.yml

@@ -1,9 +0,0 @@
-kind: pipeline
-name: default
-
-steps:
-  - name: test
-    image: rust:1.67
-    commands:
-      - cargo build --verbose --all
-      - cargo test --verbose --all

.gitignore

@@ -1,6 +1,9 @@
-/target
+# Final executable
+/run
+
+# IDE files
 .idea
+
+# Data files
 *.csv
-/*.json
-*.svg
-/data
+*.o

AIM.md

@@ -1,4 +0,0 @@
-Make a gridded structure which can support:
-
-- 3D
-- Hexagons

Cargo.toml

@@ -1,61 +0,0 @@
-[package]
-name = "rust-codebase"
-version = "0.1.0"
-edition = "2021"
-
-# See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html
-
-[[bin]]
-name = "model"
-path = "src/main.rs"
-
-[[bin]]
-name = "tools"
-path = "src/tools_cli.rs"
-
-# Set the default for crate.
-[profile.dev]
-opt-level = 1
-
-# Set the default for crate.
-[profile.release]
-debug = true
-
-# Set the default for dependencies.
-[profile.dev.package."*"]
-opt-level = 3
-
-[dependencies]
-clap = { version = "4.1.8", features = ["derive"] }
-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"

DLASystem.cpp

@@ -0,0 +1,244 @@
+//
+//  DLASystem.cpp
+//
+
+#include "DLASystem.h"
+
+// this function gets called every step,
+//   if there is an active particle then it gets moved,
+//   if not then add a particle
+void DLASystem::update() {
+    if (lastParticleIsActive == 1) {
+        moveLastParticle();
+    } else if (numParticles < endNum) {
+        addParticleOnAddCircle();
+        setParticleActive();
+    } else {
+        this->running = false;
+    }
+}
+
+
+void DLASystem::clearParticles() {
+    // delete particles and the particle list
+    for (int i = 0; i < numParticles; i++) {
+        delete particleList[i];
+    }
+    particleList.clear();
+    numParticles = 0;
+}
+
+// remove any existing particles and setup initial condition
+void DLASystem::Reset() {
+    // stop running
+    this->running = false;
+
+    clearParticles();
+
+    lastParticleIsActive = 0;
+
+    // set the grid to zero
+    for (int i = 0; i < gridSize; i++) {
+        for (int j = 0; j < gridSize; j++) {
+            grid[i][j] = 0;
+        }
+    }
+
+    // setup initial condition and parameters
+    addCircle = 10;
+    killCircle = 2.0 * addCircle;
+    clusterRadius = 0.0;
+    // add a single particle at the origin
+    double pos[] = {0.0, 0.0};
+    addParticle(pos);
+}
+
+// set the value of a grid cell for a particular position
+// note the position has the initial particle at (0,0)
+// but this corresponds to the middle of the grid array ie grid[ halfGrid ][ halfGrid ]
+void DLASystem::setGrid(double pos[], int val) {
+    int halfGrid = gridSize / 2;
+    grid[(int) (pos[0] + halfGrid)][(int) (pos[1] + halfGrid)] = val;
+}
+
+// read the grid cell for a given position
+int DLASystem::readGrid(double pos[]) {
+    int halfGrid = gridSize / 2;
+    return grid[(int) (pos[0] + halfGrid)][(int) (pos[1] + halfGrid)];
+}
+
+// check if the cluster is big enough and we should stop:
+// to be safe, we need the killCircle to be at least 2 less than the edge of the grid
+int DLASystem::checkStop() {
+    if (killCircle + 2 >= gridSize / 2) {
+        this->running = false;
+        cerr << "STOP" << endl;
+        return 1;
+    } else return 0;
+}
+
+// add a particle to the system at a specific position
+void DLASystem::addParticle(double pos[]) {
+    // create a new particle
+    Particle *p = new Particle(pos);
+    // push_back means "add this to the end of the list"
+    particleList.push_back(p);
+    numParticles++;
+
+    // pos coordinates should be -gridSize/2 < x < gridSize/2
+    setGrid(pos, 1);
+}
+
+// add a particle to the system at a random position on the addCircle
+// if we hit an occupied site then we do nothing except print a message
+// (this should never happen)
+void DLASystem::addParticleOnAddCircle() {
+    double pos[2];
+    double theta = rgen.random01() * 2 * M_PI;
+    pos[0] = ceil(addCircle * cos(theta));
+    pos[1] = ceil(addCircle * sin(theta));
+    if (readGrid(pos) == 0)
+        addParticle(pos);
+    else
+        cerr << "FAIL " << pos[0] << " " << pos[1] << endl;
+}
+
+// send back the position of a neighbour of a given grid cell
+// NOTE: there is no check that the neighbour is inside the grid,
+// this has to be done separately...
+void DLASystem::setPosNeighbour(double setpos[], double pos[], int val) {
+    switch (val) {
+        case 0:
+            setpos[0] = pos[0] + 1.0;
+            setpos[1] = pos[1];
+            break;
+        case 1:
+            setpos[0] = pos[0] - 1.0;
+            setpos[1] = pos[1];
+            break;
+        case 2:
+            setpos[0] = pos[0];
+            setpos[1] = pos[1] + 1.0;
+            break;
+        case 3:
+            setpos[0] = pos[0];
+            setpos[1] = pos[1] - 1.0;
+            break;
+    }
+}
+
+// when we add a particle to the cluster, we should update the cluster radius
+// and the sizes of the addCircle and the killCircle
+void DLASystem::updateClusterRadius(double pos[]) {
+
+    double rr = distanceFromOrigin(pos);
+    if (rr > clusterRadius) {
+        clusterRadius = rr;
+        // this is how big addCircle is supposed to be:
+        //   either 20% more than cluster radius, or at least 5 bigger.
+        double check = clusterRadius * addRatio;
+        if (check < clusterRadius + 5)
+            check = clusterRadius + 5;
+        // if it is smaller then update everything...
+        if (addCircle < check) {
+            addCircle = check;
+            killCircle = killRatio * addCircle;
+        }
+        checkStop();
+    }
+}
+
+// make a random move of the last particle in the particleList
+void DLASystem::moveLastParticle() {
+    int rr = rgen.randomInt(4);  // pick a random number in the range 0-3, which direction do we hop?
+    double newpos[2];
+
+    Particle *lastP = particleList[numParticles - 1];
+
+    setPosNeighbour(newpos, lastP->pos, rr);
+
+    if (distanceFromOrigin(newpos) > killCircle) {
+        //cerr << "#deleting particle" << endl;
+        setGrid(lastP->pos, 0);
+        particleList.pop_back();  // remove particle from particleList
+        numParticles--;
+        setParticleInactive();
+    } else if (readGrid(newpos) == 0) {
+        setGrid(lastP->pos, 0);  // set the old grid site to empty
+        // update the position
+        particleList[numParticles - 1]->pos[0] = newpos[0];
+        particleList[numParticles - 1]->pos[1] = newpos[1];
+        setGrid(lastP->pos, 1);  // set the new grid site to be occupied
+
+        // check if we stick
+        if (checkStick()) {
+            //cerr << "stick" << endl;
+            setParticleInactive();  // make the particle inactive (stuck)
+            updateClusterRadius(lastP->pos);  // update the cluster radius, addCircle, etc.
+        }
+    }
+}
+
+// check if the last particle should stick (to a neighbour)
+int DLASystem::checkStick() {
+    Particle *lastP = particleList[numParticles - 1];
+    // loop over neighbours
+    for (int i = 0; i < 4; i++) {
+        double checkpos[2];
+        setPosNeighbour(checkpos, lastP->pos, i);
+        // if the neighbour is occupied and the particle will stick probabilistically.
+        if (readGrid(checkpos) == 1 && rgen.random01() < stickProbability) {
+            return 1;
+        }
+    }
+
+    return 0;
+}
+
+
+// constructor
+DLASystem::DLASystem(const int maxParticles, const string &csvPath, const double stickProbability) {
+    cerr << "creating system, gridSize " << gridSize << endl;
+    numParticles = 0;
+    endNum = maxParticles;
+    this->stickProbability = stickProbability;
+
+    // allocate memory for the grid, remember to free the memory in destructor
+    grid = new int *[gridSize];
+    for (int i = 0; i < gridSize; i++) {
+        grid[i] = new int[gridSize];
+    }
+
+    // reset initial parameters
+    Reset();
+
+    addRatio = 1.2;   // how much bigger the addCircle should be, compared to cluster radius
+    killRatio = 1.7;   // how much bigger is the killCircle, compared to the addCircle
+
+    csv.open(csvPath);
+}
+
+// destructor
+DLASystem::~DLASystem() {
+    // strictly we should not print inside the destructor but never mind...
+    cerr << "deleting system" << endl;
+    // delete the particles
+    clearParticles();
+    // delete the grid
+    for (int i = 0; i < gridSize; i++)
+        delete[] grid[i];
+    delete[] grid;
+
+    if (csv.is_open()) {
+        csv.flush();
+        csv.close();
+    }
+}
+
+void DLASystem::exportData() {
+    csv << "x,y" << endl;
+
+    for (auto particle: particleList) {
+        csv << particle->pos[0] << "," << particle->pos[1] << endl;
+    }
+}

DLASystem.h

@@ -0,0 +1,122 @@
+#pragma once
+
+#include <iostream>
+#include <fstream>
+#include <stdio.h>
+#include <vector>
+#define _USE_MATH_DEFINES
+#include <math.h>
+#include <random>
+#include <string>
+#include <sstream>
+
+#include "Particle.h"
+#include "rnd.h"
+
+using namespace std;
+
+
+class DLASystem {
+  private:
+  // these are private variables and functions that the user will not see
+
+    // list of particles
+    vector<Particle*> particleList;
+    int numParticles;
+
+    /*
+     * The probability that a particle will stick at a given site when adjacent.
+     * */
+    double stickProbability;
+
+    // delete particles and clear the particle list
+    void clearParticles();
+
+    // size of cluster
+    double clusterRadius;
+    // these are related to the DLA algorithm
+    double addCircle;
+    double killCircle;
+  
+    // size of grid
+    static const int gridSize = 1600;
+    int **grid;  // this will be a 2d array that stores whether each site is occupied
+  
+    // random number generator, class name is rnd, instance is rgen
+    rnd rgen;
+  
+    // CSV ouput
+    ofstream csv;
+  
+    // number of particles at which the simulation will stop
+    // (the value is set in constructor)
+    int endNum;
+  
+    // the values of these variables are set in the constructor
+    double addRatio;    // how much bigger the addCircle should be, compared to cluster radius
+    double killRatio;   // how much bigger is the killCircle, compared to the addCircle
+
+  
+  public:
+  // these are public variables and functions
+
+    // update the system: if there is an active particle then move it,
+    // else create a new particle (on the adding circle)
+    void update();
+
+    // is the simulation running
+    bool running;
+  
+    // lastParticleIsActive is +1 if there is an active particle in the system, otherwise 0
+    int lastParticleIsActive;
+  
+    // constructor
+    explicit DLASystem(int maxParticles, const string &csvPath, double stickProbability);
+    // destructor
+    ~DLASystem();
+  
+    // delete all particles and reset
+    void Reset();
+
+    // this sets the seed for the random numbers
+    void setSeed(int s) { rgen.setSeed(s); }
+
+    // check whether we should stop (eg the cluster has reached the edge of the grid)
+    int checkStop();
+
+    // if pos is outside the cluster radius then set clusterRadius to be the distance to pos.
+    void updateClusterRadius( double pos[] );
+
+    // set and read grid entries associated with a given position
+    void setGrid(double pos[], int val);
+    int readGrid(double pos[]);
+
+    // return the distance of a given point from the origin
+    double distanceFromOrigin(double pos[]) {
+      return sqrt( pos[0]*pos[0] + pos[1]*pos[1] );
+    }
+
+    // set whether there is an active particle in the system or not
+    void setParticleActive()   { lastParticleIsActive = 1; }
+    void setParticleInactive() { lastParticleIsActive = 0; }
+
+    // add a particle at pos
+    void addParticle(double pos[]);
+    // add a particle at a random point on the addCircle
+    void addParticleOnAddCircle();
+
+    // assign setpos to the position of a neighbour of pos
+    // which neighbour we look at is determined by val (=0,1,2,3)
+    void setPosNeighbour(double setpos[], double pos[], int val);
+
+    // this attempts to move the last particle in the List to a random neighbour
+    // if the neighbour is occupied then nothing happens
+    // the function also checks if the moving particle should stick.
+    void moveLastParticle();
+
+    // check whether the last particle should stick
+    // currently it sticks whenever it touches another particle
+    int checkStick();
+
+    void exportData();
+};

Makefile

@@ -0,0 +1,72 @@
+# ====================================================================================== #
+#                                   From the Author                                      #
+# ====================================================================================== #
+# ! The purpose of this Makefile is to build the DLASystem project
+# ! This makefile was adapted to work with any cpp project on OSX
+
+# ====================================================================================== #
+#                                 Variables of the Makefile                              #
+# ====================================================================================== #
+
+CXX     = clang++
+
+CXXFLAGS = -Wall -Wextra -g -O0  --std=c++20
+
+IFLAGS = -I/usr/local/include -I/usr/include
+
+LFLAGS = -L/usr/local/lib -lm
+
+# ------------------------------------------
+# FOR GENERIC MAKEFILE:
+# 1 - Binary directory
+# 2 - Source directory
+# 3 - Executable name
+# 4 - Sources names
+# 5 - Dependencies names
+# ------------------------------------------
+BIN = .
+SOURCE = .
+EXEC = ./run
+SOURCES = $(wildcard $(SOURCE)/*.cpp)
+OBJECTS = $(SOURCES:.cpp=.o)
+
+# ====================================================================================== #
+#                                   Targets of the Makefile                              #
+# target_name : dependency                                                               #
+# <tabulation> command                                                                   #
+# ====================================================================================== #
+
+# ------------------------------------------
+# ! - all : Compiles everything
+# ! - help : Shows this help
+# ! - clean : erases all object files *.o
+# !           and all binary executables
+# ------------------------------------------
+all : $(BIN)/run
+
+test: $(BIN)/hllc_test
+
+help :
+	@grep -E "^# !" Makefile | sed -e 's/# !/ /g'
+
+clean:
+	rm -f $(EXEC) $(OBJECTS)
+
+# ------------------------------------------
+# Executable
+# ------------------------------------------
+$(EXEC): $(OBJECTS)
+	$(CXX) $(OBJECTS) -o $(EXEC) $(IFLAGS) $(LFLAGS)
+
+# ------------------------------------------
+# Temorary files (*.o) (IFLAGS should be added here)
+# ------------------------------------------
+$(SOURCE)/%.o: $(SOURCE)/%.cpp
+	$(CXX) $(CXXFLAGS) -c $< -o $@ $(IFLAGS) $(LFLAGS)
+
+
+
+
+
+
+

Particle.h

@@ -0,0 +1,20 @@
+#pragma once
+
+class Particle {
+  public:
+    static const int dim = 2;  // we are in two dimensions
+    double *pos;  // pointer to an array of size dim, to store the position
+  
+    // default constructor
+    Particle() {
+      pos = new double[dim];
+    }
+    // constructor, with a specified initial position
+    Particle(double set_pos[]) {
+      pos = new double[dim];
+      for (int d=0;d<dim;d++)
+        pos[d]=set_pos[d];
+    }
+    // destructor
+    ~Particle() { delete[] pos; }
+};

README.md

@@ -1,8 +1,4 @@
-# DLA Generic Model
+# DLA 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.
+This currently contains The Initially Provided Code for the DLA model which will be used as the source of truth for
+further alterations once verified.

mainDLA.cpp

@@ -0,0 +1,29 @@
+#include <iostream>
+
+#include "DLASystem.h"
+
+int main(int argc, char **argv) {
+    if (argc != 5) {
+        cerr << "Usage " << argv[0] << " <seed> <maxParticles> <stickProbability> <csvPath>" << endl;
+        return 1;
+    }
+
+    int seed = std::stoi(argv[1]);
+    int maxParticles = std::stoi(argv[2]);
+    double stickProbability = std::stod(argv[3]);
+    string csvPath = argv[4];
+    std::cerr << "Setting seed " << seed << endl;
+
+    // create the system
+    auto *sys = new DLASystem(maxParticles, csvPath, stickProbability);
+    sys->setSeed(seed);
+    sys->running = true;
+
+    while (sys->running) {
+        sys->update();
+    }
+
+    sys->exportData();
+
+    return 0;
+}

rnd.h

@@ -0,0 +1,33 @@
+#pragma once
+#include <random>
+
+// ... don't worry how this all works
+// ... member functions that you may want to use:
+//       random01() returns a random double between 0 and 1
+//       randomInt(max) returns a random int between 0 and max-1 (inclusive)
+
+class rnd {
+private:
+	// nuts and bolts.. should not need to touch this.
+	std::default_random_engine generator;
+	int genMax;
+	std::uniform_int_distribution<int> *intmax;
+	std::uniform_real_distribution<double> *real01;
+
+public:
+	// constructor
+	rnd() {
+		genMax = 0x7fffffff;
+		//cout << "genMax is " << generator.max() << endl;
+		intmax = new std::uniform_int_distribution<int>(0, genMax);
+		real01 = new std::uniform_real_distribution<double>(0.0, 1.0);
+	}
+	// destructor
+	~rnd() { delete intmax; delete real01; }
+
+	// set the random seed
+	void   setSeed(int seed) { generator.seed(seed); }
+	// member functions for generating random double in [0,1] and random integer in [0,max-1]
+	double random01() { return (*real01)(generator); }
+	int    randomInt(int max) { return (*intmax)(generator) % max; }
+};

rust-toolchain.toml

@@ -1,2 +0,0 @@
-[toolchain]
-channel = "nightly"

src/cli/cli.rs

@@ -1,88 +0,0 @@
-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

@@ -1,45 +0,0 @@
-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

@@ -1,55 +0,0 @@
-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/main.rs

@@ -1,171 +0,0 @@
-#![feature(array_zip)]
-#![feature(generic_const_exprs)]
-#![feature(let_chains)]
-
-use clap::Parser;
-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};
-
-mod system;
-mod surface_probability_measure;
-mod cli;
-
-fn main() {
-    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);
-
-    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,
-            );
-
-            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();
-        }
-    }
-}

src/surface_probability_measure.rs

@@ -1,70 +0,0 @@
-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/mod.rs

@@ -1,30 +0,0 @@
-use std::ops::Add;
-use serde::Serialize;
-
-pub mod walker;
-pub mod spawner;
-pub mod sticker;
-pub mod model;
-pub mod spaces;
-
-pub trait Position: Add<Output=Self> + Serialize + Clone {
-    const DIM: usize;
-    // type Cartesian;
-
-    fn zero() -> 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: Position> {
-    fn is_occupied(&self, position: &P) -> bool;
-    fn deposit(&mut self, position: &P);
-}

src/system/model.rs

@@ -1,140 +0,0 @@
-use rand::prelude::*;
-use serde::{Deserialize, Serialize};
-use crate::system::{Position, Storage};
-use crate::system::spawner::Spawner;
-use crate::system::sticker::Sticker;
-use crate::system::walker::Walker;
-
-#[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,
-
-    pub frame: usize,
-    pub running: bool,
-    pub history: Vec<HistoryLine<P>>,
-
-    max_particles: usize,
-    particles: Vec<P>,
-    active_particle: Option<P>,
-
-    add_ratio: f32,
-    add_circle: f32,
-    kill_ratio: f32,
-    kill_circle: f32,
-    cluster_radius: f32,
-}
-
-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,
-            max_particles,
-            running: true,
-            spawner,
-            sticker,
-
-            space,
-            walker,
-
-            frame: 0,
-            particles: vec![],
-            history: vec![],
-            active_particle: None,
-
-            add_ratio: 1.2,
-            kill_ratio: 1.7,
-
-            add_circle: 10.0,
-            kill_circle: 20.0,
-            cluster_radius: 0.0,
-        };
-
-        sys.deposit(&P::zero());
-
-        sys
-    }
-
-    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 {
-            self.spawn_particle();
-        } else {
-            self.running = false;
-        }
-    }
-
-    fn move_particle(&mut self) {
-        let current_position = &self
-            .active_particle
-            .clone()
-            .expect("No active particle");
-
-        let next_position = self.walker.walk(&mut self.rng, current_position);
-        let distance = next_position.abs();
-
-        if distance > self.kill_circle {
-            self.active_particle = None;
-        } 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;
-            } else {
-                self.active_particle.replace(next_position);
-            }
-        }
-    }
-
-    fn spawn_particle(&mut self) {
-        let position = self.spawner.spawn(&mut self.rng, self.add_circle);
-
-        if !self.space.is_occupied(&position) {
-            self.active_particle = Some(position);
-        }
-    }
-
-    fn deposit(&mut self, p0: &P) {
-        self.particles.push(p0.clone());
-        self.space.deposit(p0);
-
-        let distance = p0.abs();
-        if distance > self.cluster_radius {
-            self.cluster_radius = distance;
-
-            let new_add_circle = (self.cluster_radius * self.add_ratio).max(self.cluster_radius + 5.0);
-            if self.add_circle < new_add_circle {
-                self.add_circle = new_add_circle;
-                self.kill_circle = self.kill_ratio * self.add_circle;
-            }
-        }
-
-        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/spaces/continuous_2d.rs

@@ -1,102 +0,0 @@
-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

@@ -1,106 +0,0 @@
-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,70 +0,0 @@
-use std::ops::Add;
-use num_integer::Roots;
-use num_traits::Pow;
-use crate::system::{GriddedPosition, Position};
-use serde::{Serialize, Deserialize};
-
-#[derive(Clone, Debug, PartialEq, Eq, Serialize, Deserialize)]
-pub struct HexPosition {
-    pub q: i32,
-    pub r: i32,
-}
-
-impl Add for HexPosition {
-    type Output = Self;
-
-    fn add(self, rhs: Self) -> Self::Output {
-        HexPosition { q: self.q + rhs.q, r: self.r + rhs.r }
-    }
-}
-
-impl GriddedPosition for HexPosition {
-    const NEIGHBOURS: u32 = 6;
-
-    fn neighbour(&self, neighbour_index: u32) -> Self {
-        let neighbour_index = neighbour_index as usize;
-
-        const OFFSETS: [(i32, i32); 6] = [
-            (1, 0), (1, -1), (0, -1),
-            (-1, 0), (-1, 1), (0, 1),
-        ];
-
-        self.clone() + HexPosition { q: OFFSETS[neighbour_index].0, r: OFFSETS[neighbour_index].1 }
-    }
-
-    fn linear_index(&self, grid_size: u32) -> 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

@@ -1,87 +0,0 @@
-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

@@ -1,10 +0,0 @@
-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

@@ -1,148 +0,0 @@
-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

@@ -1,199 +0,0 @@
-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

@@ -1,24 +0,0 @@
-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

@@ -1,58 +0,0 @@
-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

@@ -1,39 +0,0 @@
-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,147 +0,0 @@
-use std::hash::Hash;
-use itertools::Itertools;
-use rand::distributions::Slice;
-use rand::prelude::Rng;
-use crate::system::{GriddedPosition, Position};
-use crate::system::spaces::square_grid::{Grid2D, Grid3D};
-
-pub trait Walker<P: Position> {
-    fn walk<R: Rng>(&self, rng: &mut R, position: &P) -> P;
-}
-
-pub struct LocalRandomWalker;
-
-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))
-    }
-}
-
-pub struct DiagonalRandomWalker;
-
-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 },
-        ];
-
-        position.clone() + OFFSETS[rng.gen_range(0..OFFSETS.len())].clone()
-    }
-}
-
-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 })
-            }
-        }
-    }
-
-    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

@@ -1,33 +0,0 @@
-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

@@ -1,134 +0,0 @@
-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

@@ -1,31 +0,0 @@
-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

@@ -1,125 +0,0 @@
-#![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

@@ -1,83 +0,0 @@
-#![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(())
-}