Hex work!

.drone.yml

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

.gitignore

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

AIM.md

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

Cargo.toml

@@ -0,0 +1,71 @@
+[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
+
+[lib]
+name = "dla"
+crate-type = ["staticlib"]
+path = "src/clib.rs"
+
+[[bin]]
+name = "model"
+path = "src/main.rs"
+
+[[bin]]
+name = "ui"
+path = "src/ui.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"] }
+bevy = { version = "0.10.0" }
+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 = "0.32.2"
+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

@@ -1,319 +0,0 @@
-//
-//  DLASystem.cpp
-//
-
-#include "DLASystem.h"
-
-// colors
-namespace colours {
-	GLfloat blue[] = { 0.1, 0.3, 0.9, 1.0 };   // blue
-	GLfloat red[] = { 1.0, 0.2, 0.1, 0.2 };   // red
-	GLfloat green[] = { 0.3, 0.6, 0.3, 1.0 };     // green
-	GLfloat paleGrey[] = { 0.7, 0.7, 0.7, 1.0 };     // green
-	GLfloat darkGrey[] = { 0.2, 0.2, 0.2, 1.0 };     // green
-}
-
-
-// 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();
-	}
-	if (lastParticleIsActive == 0 || slowNotFast == 1)
-		glutPostRedisplay(); //Tell GLUT that the display has changed
-}
-
-
-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
-	running = 0;
-
-	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 view
-	int InitialViewSize = 40;
-	setViewSize(InitialViewSize);
-
-}
-
-// 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) {
-		pauseRunning();
-		cout << "STOP" << endl;
-		glutPostRedisplay(); // update display
-		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
-		cout << "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;
-	}
-}
-
-// if the view is smaller than the kill circle then increase the view area (zoom out)
-void DLASystem::updateViewSize() {
-	double mult = 1.2;
-	if (viewSize < 2.0*killCircle) {
-		setViewSize(viewSize * mult);
-	}
-}
-
-// set the view to be the size of the add circle (ie zoom in on the cluster)
-void DLASystem::viewAddCircle() {
-	setViewSize(2.0*addCircle);  // factor of 2 is to go from radius to diameter
-}
-
-// 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;
-			updateViewSize();
-		}
-		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) {
-		//cout << "#deleting particle" << endl;
-		setGrid(lastP->pos, 0);
-		particleList.pop_back();  // remove particle from particleList
-		numParticles--;
-		setParticleInactive();
-	}
-	// check if destination is empty
-	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()) {
-			//cout << "stick" << endl;
-			setParticleInactive();  // make the particle inactive (stuck)
-			updateClusterRadius(lastP->pos);  // update the cluster radius, addCircle, etc.
-
-			if (numParticles % 100 == 0 && logfile.is_open()) {
-				logfile << numParticles << " " << clusterRadius << endl;
-			}
-		}
-	}
-	else {
-		// if we get to here then we are trying to move to an occupied site
-		// (this should never happen as long as the sticking probability is 1.0)
-		cout << "reject " << rr << endl;
-		cout << lastP->pos[0] << " " << lastP->pos[1] << endl;
-		//cout << newpos[0] << " " << newpos[1] << " " << (int)newpos[0] << endl;
-		//printOccupied();
-	}
-}
-
-// check if the last particle should stick (to a neighbour)
-int DLASystem::checkStick() {
-	Particle *lastP = particleList[numParticles - 1];
-	int result = 0;
-	// loop over neighbours
-	for (int i = 0; i < 4; i++) {
-		double checkpos[2];
-		setPosNeighbour(checkpos, lastP->pos, i);
-		// if the neighbour is occupied...
-		if (readGrid(checkpos) == 1)
-			result = 1;
-	}
-	return result;
-}
-
-
-// constructor
-DLASystem::DLASystem(Window *set_win) {
-	cout << "creating system, gridSize " << gridSize << endl;
-	win = set_win;
-	numParticles = 0;
-	endNum = 1000;
-
-	// 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];
-	}
-	slowNotFast = 1;
-	// 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
-
-	// this opens a logfile, if we want to...
-	//logfile.open("opfile.txt");
-}
-
-// destructor
-DLASystem::~DLASystem() {
-	// strictly we should not print inside the destructor but never mind...
-	cout << "deleting system" << endl;
-	// delete the particles
-	clearParticles();
-	// delete the grid
-	for (int i = 0; i < gridSize; i++)
-		delete[] grid[i];
-	delete[] grid;
-
-	if (logfile.is_open())
-		logfile.close();
-}
-
-
-
-// this draws the system
-void DLASystem::DrawSquares() {
-
-	// draw the particles
-	double halfSize = 0.5;
-	for (int p = 0; p < numParticles; p++) {
-		double *vec = particleList[p]->pos;
-		glPushMatrix();
-		if (p == numParticles - 1 && lastParticleIsActive == 1)
-			glColor4fv(colours::red);
-		else if (p == 0)
-			glColor4fv(colours::green);
-		else
-			glColor4fv(colours::blue);
-		glRectd(drawScale*(vec[0] - halfSize),
-			drawScale*(vec[1] - halfSize),
-			drawScale*(vec[0] + halfSize),
-			drawScale*(vec[1] + halfSize));
-		glPopMatrix();
-	}
-
-	// print some information (at top left)
-	// this ostringstream is a way to create a string with numbers and words (similar to cout << ... )
-	ostringstream str;
-	str << "num " << numParticles << " size " << clusterRadius;
-
-	// print the string
-	win->displayString(str, -0.9, 0.9, colours::red);
-
-	// if we are paused then print this (at bottom left)
-	if (running == 0) {
-		ostringstream pauseStr;
-		pauseStr << "paused";
-		win->displayString(pauseStr, -0.9, -0.9, colours::red);
-	}
-
-}

DLASystem.h

@@ -1,154 +0,0 @@
-#pragma once
-
-#include <GLUT/glut.h>
-#include <iostream>
-#include <fstream>
-#include <stdio.h>
-#include <vector>
-#define _USE_MATH_DEFINES
-#include <math.h>
-#include <random>
-#include <string>
-#include <sstream>
-
-#include "Window.h"
-#include "Particle.h"
-#include "rnd.h"
-
-using namespace std;
-
-
-class DLASystem {
-  private:
-  // these are private variables and functions that the user will not see
-  
-    Window *win;  // window in which the system is running
-  
-    // list of particles
-    vector<Particle*> particleList;
-    int numParticles;
-
-    // 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
-  
-    // the window draws only part of the grid, viewSize controls how much...
-    double viewSize;
-    double drawScale;
-  
-    // random number generator, class name is rnd, instance is rgen
-    rnd rgen;
-  
-    // output file (not used at the moment)
-    ofstream logfile;
-  
-    // 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();
-
-    // draw particles as squares
-    void DrawSquares();
-  
-    // is the simulation running (1) or paused (0) ?
-    int running;
-  
-    // slowNotFast is +1 for slow running, 0 for fast
-    int slowNotFast;
-
-    // lastParticleIsActive is +1 if there is an active particle in the system, otherwise 0
-    int lastParticleIsActive;
-  
-    // constructor
-    DLASystem(Window *set_win);
-    // 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();
-  
-    // stop/start the algorithm
-    void setRunning() { if ( checkStop()==0 ) running = 1; }
-    void pauseRunning() { running = 0; }
-    // set whether it runs fast or slow
-    void setSlow() { slowNotFast = 1; }
-    void setFast() { slowNotFast = 0; }
-    void setSuperFast() { slowNotFast = -1; }
-
-    // set which part of the grid is visible on the screen
-    // basically the screen shows co-ordinates -vv < x < vv
-    // where vv is the input value
-    void setViewSize(double vv) { viewSize = vv; drawScale = 2.0/viewSize; }
-  
-    // if the killcircle is almost as big as the view then increase the view
-    void updateViewSize();
-
-    // set the view to be the approx size of the addCircle
-    void viewAddCircle();
-
-    // 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();
-
-    // set the background colour for the window
-    // it would be better for an OOP philosophy to make these member functions for the Window class
-    // but we are being a bit lazy here
-    void setWinBackgroundWhite() { glClearColor(1.0, 1.0, 1.0, 1.0); }
-    void setWinBackgroundBlack() { glClearColor(0.0, 0.0, 0.0, 0.0); }
-};

Makefile

@@ -1,72 +0,0 @@
-# ====================================================================================== #
-#                                   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 
-
-IFLAGS = -I/usr/local/include -I/usr/include
-
-LFLAGS = -L/usr/local/lib -lm -framework OpenGL -framework GLUT	
-
-# ------------------------------------------
-# 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

@@ -1,20 +0,0 @@
-#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; }
-};

Window.cpp

@@ -1,26 +0,0 @@
-#include "Window.h"
-
-// constructor
-Window::Window(int set_size[],  string &set_title) {
-  size[0]=set_size[0]; size[1]=set_size[1];
-  title=set_title;
-
-  locateOnScreen();
-  glutInitWindowSize(size[0],size[1]);
-  glutInitWindowPosition(pos[0],pos[1]);
-  glutInitDisplayMode(GLUT_RGBA | GLUT_DOUBLE);
-  glutCreateWindow(title.c_str());
-
-  // sets the background to black
-  glClearColor(0.0, 0.0, 0.0, 0.0);
-}
-
-// print a string at a given position, don't worry about how this works...
-void Window::displayString( ostringstream &str, double x, double y, GLfloat col[]) {
-    string localString = str.str();
-    glColor4fv(col);
-    glRasterPos2d(x,y);
-    for (int i = 0; i < localString.length(); i++) {
-      glutBitmapCharacter(GLUT_BITMAP_8_BY_13, localString[i]);
-    }
-}

Window.h

@@ -1,34 +0,0 @@
-#pragma once
-
-#include <GLUT/glut.h>
-#include <iostream>
-#include <string>
-#include <sstream>
-
-using namespace std;
-
-class Window {
-  public:
-    string title;
-    int size[2];
-    int pos[2];
-  
-    void locateOnScreen() {
-      // the fx sets where on the screen the window will appear
-      // (values should be between 0 and 1)
-      double fx[] = { 0.7,0.5 };
-      pos[0] = (glutGet(GLUT_SCREEN_WIDTH) - size[0]) * fx[0];
-      pos[1] = (glutGet(GLUT_SCREEN_HEIGHT) - size[1]) * fx[1];
-    }
-  
-    // constructor, size is in pixels
-    Window(int set_size[],  string &set_title);
-
-    // function which prints a string to the screen, at a given position, with a given color
-    // note position is "absolute", not easy to get two strings spaced one above each other like this
-    void displayString( ostringstream &str, double x, double y, GLfloat col[]);
-
-
-};
-
-

build.rs

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

@@ -0,0 +1,11 @@
+#include <cstdarg>
+#include <cstdint>
+#include <cstdlib>
+#include <ostream>
+#include <new>
+
+extern "C" {
+
+bool dla_rust_disabled();
+
+} // extern "C"

mainDLA.cpp

@@ -1,169 +0,0 @@
-#include <GLUT/glut.h>
-#include <iostream>
-#include <stdio.h>
-#include <vector>
-#include <math.h>
-#include <string>
-
-#include "DLASystem.h"
-#include "Window.h" 
-
-using namespace std;
-
-// functions which are needed for openGL go into a namespace so that we can identify them
-namespace drawFuncs {
-  void handleKeypress(unsigned char key, int x, int y);
-  void display(void);
-  void update(int val);
-  void introMessage();
-}
-
-// this is a global pointer, which is how we access the system itself
-DLASystem *sys;
-
-int main(int argc, char **argv) {
-  // turn on glut
-	glutInit(&argc, argv);
-
-  int window_size[] = { 480,480 };
-  string window_title("simple DLA simulation");
-
-  // create a window
-  Window *win = new Window(window_size,window_title);
-
-  // create the system
-  sys = new DLASystem(win);
-  
-  // this is the seed for the random numbers
-  int seed = 6;
-  cout << "setting seed " << seed << endl;
-  sys->setSeed(seed);
-  
-  // print the "help" message to the console
-  drawFuncs::introMessage();
-  
-  // tell openGL how to redraw the screen and respond to the keyboard
-	glutDisplayFunc(  drawFuncs::display );
-	glutKeyboardFunc( drawFuncs::handleKeypress );
-  
-  // tell openGL to do its first update after waiting 10ms
-  int wait = 10;
-  int val = 0;
-	glutTimerFunc(wait, drawFuncs::update, val);
-
-  // start the openGL stuff
- 	glutMainLoop();
-
-  return 0;
-}
-
-// this is just a help message
-void drawFuncs::introMessage() {
-        cout << "Keys (while in graphics window):" << endl << "  q or e to quit (or exit)" << endl;
-        cout << "  h to print this message (help)" << endl;
-        cout << "  u for a single update" << endl;
-        cout << "  g to start running (go)" << endl;
-        cout << "  p to pause running" << endl;
-        cout << "  s to run in slow-mode" << endl;
-        cout << "  f to run in fast-mode" << endl;
-        cout << "  r to clear everything (reset)" << endl;
-        cout << "  z to pause and zoom in" << endl;
-        cout << "  w or b to change background colour to white or black" << endl;
-}
-
-// openGL function deals with the keyboard
-void drawFuncs::handleKeypress(unsigned char key, int x, int y) {
-	switch (key) {
-  case 'h':
-    drawFuncs::introMessage();
-    break;
-	case 'q':
-	case 'e':
-		cout << "Exiting..." << endl;
-    // delete the system
-    delete sys;
-		exit(0);
-		break;
-  case 'p':
-    cout << "pause" << endl;
-    sys->pauseRunning();
-    break;
-  case 'g':
-    cout << "go" << endl;
-    sys->setRunning();
-    glutTimerFunc(0, drawFuncs::update, 0);
-    break;
-  case 's':
-    cout << "slow" << endl;
-    sys->setSlow();
-    break;
-  case 'w':
-    cout << "white" << endl;
-    sys->setWinBackgroundWhite();
-    break;
-  case 'b':
-    cout << "black" << endl;
-    sys->setWinBackgroundBlack();
-    break;
-  case 'f':
-    cout << "fast" << endl;
-    sys->setFast();
-    break;
-  case 'r':
-    cout << "reset" << endl;
-    sys->Reset();
-    break;
-  case 'z':
-    cout << "zoom" << endl;
-    sys->pauseRunning();
-    sys->viewAddCircle();
-    break;
-  case 'u':
-    cout << "upd" << endl;
-    sys->Update();
-    break;
-	}
-  // tell openGL to redraw the window
-	glutPostRedisplay();
-}
-
-// this function gets called whenever the algorithm should do its update
-void drawFuncs::update(int val) {
-  int wait;  // time to wait between updates (milliseconds)
-  
-  if ( sys->running ) {
-    if ( sys->slowNotFast == 1)
-      wait = 10;
-    else
-      wait = 0;
-    
-    sys->Update();
-  
-    // tell openGL to call this funtion again after "wait" milliseconds
-    glutTimerFunc(wait, drawFuncs::update, 0);
-  }
-  
-}
-
-// this function redraws the window when necessary
-void drawFuncs::display() {
-    	//  Clear the window or more specifically the frame buffer...
-      //  This happens by replacing all the contents of the frame
-      //  buffer by the clear color (black in our case)
-      glClear(GL_COLOR_BUFFER_BIT);
-
-      // this puts the camera at the origin (not sure why) with (I think) z axis out of page and y axis up
-      // there is also the question of the GL perspective which is not set up in any clear way at the moment
-      glMatrixMode(GL_MODELVIEW);
-      glLoadIdentity();
-      gluLookAt(0.0, 0.0, 1.0,   /* camera position */
-        0.0, 0.0, -1.0,        /* point to look at */
-        0.0, 1.0, 0.0);		   /* up direction */
-
-      //sys->DrawSpheres();
-      sys->DrawSquares();
-
-      //  Swap contents of backward and forward frame buffers
-      glutSwapBuffers();
-}
-

readme.txt

@@ -1,27 +0,0 @@
-To run the code on a Mac through the command line:
-
-1. Download Makefile from this section
-
-2. Put Makefile in same directory as DLASystem.cpp (! Remove the .txt extension if present !)
-
-3. Relative to MS Visual Studio code, these files include the following changes: DLASystem.h Window.h mainDLA.cpp
-
-change #include <gl/freeglut.h>
-
-to #include <GLUT/glut.h>
-
-4. Open the terminal in the same folder as these files
-
-5. Type "make" to build project
-
-6. Run the program by typing "./run"
-
-7. To edit your files, run gedit, save the code, and recompile with "make -B" before running the program again.
-
-
-
-Thanks to Eliot Ayache for figuring this out!
-
-
-
-PS For Coursework 2: In addition, there needs to be a change in Line 72 of MAIN_Source.cpp file from "glutLeaveMainLoop();" to "exit(1);"

rnd.h

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

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

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

@@ -0,0 +1,4 @@
+#[no_mangle]
+pub extern "C" fn dla_rust_disabled() -> bool {
+    true
+}

src/main.rs

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

@@ -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/mod.rs

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

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

@@ -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

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

@@ -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,9 @@
+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;

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

@@ -0,0 +1,148 @@
+use std::hash::Hash;
+use bevy::render::render_resource::encase::private::RuntimeSizedArray;
+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

@@ -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,135 @@
+use std::fs::File;
+use std::os::unix::fs::symlink;
+use bevy::tasks::ParallelSlice;
+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", max_size)
+        .set("height", max_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

@@ -0,0 +1,108 @@
+//! A simplified implementation of the classic game "Breakout".
+
+use bevy::{
+    prelude::*,
+};
+
+// These constants are defined in `Transform` units.
+// Using the default 2D camera they correspond 1:1 with screen pixels.
+const GAP_BETWEEN_PADDLE_AND_FLOOR: f32 = 60.0;
+
+const LEFT_WALL: f32 = -450.;
+const RIGHT_WALL: f32 = 450.;
+// y coordinates
+const BOTTOM_WALL: f32 = -300.;
+const TOP_WALL: f32 = 300.;
+
+const PARTICLE_SIZE: Vec2 = Vec2::new(10., 10.);
+
+// These values are exact
+const GAP_BETWEEN_PADDLE_AND_BRICKS: f32 = 270.0;
+const GAP_BETWEEN_BRICKS: f32 = 5.0;
+// These values are lower bounds, as the number of bricks is computed
+const GAP_BETWEEN_BRICKS_AND_CEILING: f32 = 20.0;
+const GAP_BETWEEN_BRICKS_AND_SIDES: f32 = 20.0;
+
+const BACKGROUND_COLOR: Color = Color::rgb(0.9, 0.9, 0.9);
+const BRICK_COLOR: Color = Color::rgb(0.5, 0.5, 1.0);
+
+fn main() {
+    App::new()
+        .add_plugins(DefaultPlugins)
+        .insert_resource(ClearColor(BACKGROUND_COLOR))
+        .add_startup_system(setup)
+        .add_system(bevy::window::close_on_esc)
+        .run();
+}
+
+#[derive(Component)]
+struct Brick;
+
+// Add the game's entities to our world
+fn setup(
+    mut commands: Commands,
+) {
+    // Camera
+    commands.spawn(Camera2dBundle::default());
+
+    // Paddle
+    let paddle_y = BOTTOM_WALL + GAP_BETWEEN_PADDLE_AND_FLOOR;
+
+    // Bricks
+    // Negative scales result in flipped sprites / meshes,
+    // which is definitely not what we want here
+    assert!(PARTICLE_SIZE.x > 0.0);
+    assert!(PARTICLE_SIZE.y > 0.0);
+
+    let total_width_of_bricks = (RIGHT_WALL - LEFT_WALL) - 2. * GAP_BETWEEN_BRICKS_AND_SIDES;
+    let bottom_edge_of_bricks = paddle_y + GAP_BETWEEN_PADDLE_AND_BRICKS;
+    let total_height_of_bricks = TOP_WALL - bottom_edge_of_bricks - GAP_BETWEEN_BRICKS_AND_CEILING;
+
+    assert!(total_width_of_bricks > 0.0);
+    assert!(total_height_of_bricks > 0.0);
+
+    // Given the space available, compute how many rows and columns of bricks we can fit
+    let n_columns = (total_width_of_bricks / (PARTICLE_SIZE.x + GAP_BETWEEN_BRICKS)).floor() as usize;
+    let n_rows = (total_height_of_bricks / (PARTICLE_SIZE.y + GAP_BETWEEN_BRICKS)).floor() as usize;
+    let n_vertical_gaps = n_columns - 1;
+
+    // Because we need to round the number of columns,
+    // the space on the top and sides of the bricks only captures a lower bound, not an exact value
+    let center_of_bricks = (LEFT_WALL + RIGHT_WALL) / 2.0;
+    let left_edge_of_bricks = center_of_bricks
+        // Space taken up by the bricks
+        - (n_columns as f32 / 2.0 * PARTICLE_SIZE.x)
+        // Space taken up by the gaps
+        - n_vertical_gaps as f32 / 2.0 * GAP_BETWEEN_BRICKS;
+
+    // In Bevy, the `translation` of an entity describes the center point,
+    // not its bottom-left corner
+    let offset_x = left_edge_of_bricks + PARTICLE_SIZE.x / 2.;
+    let offset_y = bottom_edge_of_bricks + PARTICLE_SIZE.y / 2.;
+
+    for row in 0..n_rows {
+        for column in 0..n_columns {
+            let brick_position = Vec2::new(
+                offset_x + column as f32 * (PARTICLE_SIZE.x + GAP_BETWEEN_BRICKS),
+                offset_y + row as f32 * (PARTICLE_SIZE.y + GAP_BETWEEN_BRICKS),
+            );
+
+            // brick
+            commands.spawn((
+                SpriteBundle {
+                    sprite: Sprite {
+                        color: BRICK_COLOR,
+                        ..default()
+                    },
+                    transform: Transform {
+                        translation: brick_position.extend(0.0),
+                        scale: Vec3::new(PARTICLE_SIZE.x, PARTICLE_SIZE.y, 1.0),
+                        ..default()
+                    },
+                    ..default()
+                },
+                Brick,
+            ));
+        }
+    }
+}