joshuacoles/compb-dla-model
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Compare commits

merge into: joshuacoles:main
Branches Tags
joshuacoles:main
joshuacoles:rust
joshuacoles:3d-off-axis
joshuacoles:3d
joshuacoles:3d-off-axis-neighbours
joshuacoles:3d-direct-neighbours
joshuacoles:rust-sticking-probability
joshuacoles:stick-probability
joshuacoles:minimal-viable-alteration
joshuacoles:initially-provided-code
...
pull from: joshuacoles:rust
Branches Tags
joshuacoles:rust
joshuacoles:3d-off-axis
joshuacoles:3d
joshuacoles:main
joshuacoles:3d-off-axis-neighbours
joshuacoles:3d-direct-neighbours
joshuacoles:rust-sticking-probability
joshuacoles:stick-probability
joshuacoles:minimal-viable-alteration
joshuacoles:initially-provided-code

No commits in common. "main" and "rust" have entirely different histories.
main ... rust

36 changed files with 4671 additions and 530 deletions
Show Stats Expand all files Collapse all files

9
.drone.yml Normal file
View File

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

11
.gitignore vendored
View File

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

4
AIM.md Normal file
View File

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

2589
Cargo.lock generated Normal file
View File

File diff suppressed because it is too large Load Diff

61
Cargo.toml Normal file
View File

@ -0,0 +1,61 @@
[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"

244
DLASystem.cpp
View File

@ -1,244 +0,0 @@
//
// 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;
}
}

122
DLASystem.h
View File

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

72
Makefile
View File

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

20
Particle.h
View File

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

10
README.md
View File

@ -1,4 +1,8 @@
# DLA Model # DLA Generic Model
This currently contains The Initially Provided Code for the DLA model which will be used as the source of truth for A generic pluggable model for diffusion limited aggregation. Produces two executables,
further alterations once verified.
- `./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.

29
mainDLA.cpp
View File

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

33
rnd.h
View File

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

2
rust-toolchain.toml Normal file
View File

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

88
src/cli/cli.rs Normal file
View File

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

45
src/cli/mod.rs Normal file
View File

@ -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)
}
_ => {}
}
}
}

55
src/cli/output.rs Normal file
View File

@ -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");
}

171
src/main.rs Normal file
View File

@ -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();
}
}
}

70
src/surface_probability_measure.rs Normal file
View File

@ -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:?}");
}
}

30
src/system/mod.rs Normal file
View File

@ -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);
}

140
src/system/model.rs Normal file
View File

@ -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() });
}
}

102
src/system/spaces/continuous_2d.rs Normal file
View File

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

106
src/system/spaces/continuous_3d.rs Normal file
View File

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

70
src/system/spaces/hexagonal.rs Normal file
View File

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

87
src/system/spaces/kd_grid.rs Normal file
View File

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

10
src/system/spaces/mod.rs Normal file
View File

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

148
src/system/spaces/nalg.rs Normal file
View File

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

199
src/system/spaces/square_grid.rs Normal file
View File

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

24
src/system/spaces/vector_storage.rs Normal file
View File

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

58
src/system/spawner.rs Normal file
View File

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

39
src/system/sticker.rs Normal file
View File

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

147
src/system/walker.rs Normal file
View File

@ -0,0 +1,147 @@
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);
}

33
src/tools/analysis.rs Normal file
View File

@ -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);
}

134
src/tools/boxcount.rs Normal file
View File

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

31
src/tools/mod.rs Normal file
View File

@ -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()
}

125
src/tools/render.rs Normal file
View File

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

83
src/tools_cli.rs Normal file
View File

@ -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(())
}