# Physical simulation of diffusion-limited aggregates

The following code generates Diffusion Limited Aggregates on a two-dimensional square lattice. Some of the code has been omitted (e.g. support for differing lattice types and dimensions) for code-length reasons and due to the fact that the key thing here is optimisation improvement suggestions. The full code is on GitHub.

#pragma once
#include <algorithm>
#include <cmath>
#include <unordered_map>
#include <queue>
#include <random>
#include <stdexcept>
#include <utility>
#include <vector>

// hash function-object used for std::unordered_map of std::pair
struct pair_hash {
template<typename _Ty1, typename _Ty2> std::size_t operator()(const
std::pair<_Ty1, _Ty2>& _p) const {
return 51 + std::hash<_Ty1>()(_p.first)*51 + std::hash<_Ty2>()(_p.second);
}
}

// comparator function-object used for comparing distances from origin
struct distance_comparator {
bool operator()(const std::pair<int,int>& _lhs, const std::pair<int,int>& _rhs) const {
return (_lhs.first*_lhs.first + _lhs.second*_lhs.second)
< (_rhs.first*_rhs.first + _rhs.second*_rhs.second);
}
}

class DLAContainer {
public:
DLAContainer(const double& _coeff_stick = 1.0);
DLAContainer(const DLAContainer& _other);
DLAContainer(DLAContainer&& _other);

size_t size() const noexcept;
double get_coeff_stick() const noexcept;
void set_coeff_stick(const double& _coeff_stick);

// generates aggregate of _n particles
void generate(size_t _n);
void clear();
private:
// map to store aggregate point co-ordinates as keys and
// order of adding to container as values
std::unordered_map<std::pair<int,int>, size_t, pair_hash> aggregate_map;
// priority queue for retrieving co-ordinates of aggregate
// particle furthest from origin in constant time
std::priority_queue<std::pair<int,int>, std::vector<std::pair<int,int>>, distance_comparator> aggregate_pq;
double coeff_stick;
std::mt19937 mt_eng;

void spawn_particle(int& _x, int& _y, int& _spawn_diam,
std::uniform_real_distribution<>& _dist) noexcept;
void update_particle_position(int& _x, int& _y, const double& _movement_choice) const noexcept;
bool lattice_boundary_collision(int& _x, int& _y, const int& _prev_x,
const int& _prev_y, const int& _spawn_diam) const noexcept;
bool aggregate_collision(const int& _x, const int& _y, const int& _prev_x,
const int& _prev_y, const double& _sticky_pr, size_t& _count);
};


# DLAContainer Class Implementation File

#include "DLAContainer.h"

DLAContainer::DLAContainer(const double& _coeff_stick) {
set_coeff_stick(_coeff_stick);
std::random_device rd;
mt_eng = std::mt19937(rd());
}

DLAContainer::DLAContainer(const DLAContainer& _other) {
coeff_stick = _other.coeff_stick;
mt_eng = _other.mt_eng;
aggregate_map = _other.aggregate_map;
aggregate_pq = _other.aggregate_pq;
}

DLAContainer::DLAContainer(DLAContainer&& _other) {
coeff_stick = _other.coeff_stick;
mt_eng = _other.mt_eng;
aggregate_map = _other.aggregate_map;
aggregate_pq = _other.aggregate_pq;
_other.coeff_stick = 1.0;
_other.mt_eng = std::mt19937();
_other.aggregate_map = std::unordered_map<std::pair<int, int>, size_t, pair_hash>();
_other.aggregate_pq = std::priority_queue<std::pair<int, int>, std::vector<std::pair<int, int>>, distance_comparator>();
}

size_t DLAContainer::size() const noexcept {
return aggregate_map.size();
}

double DLAContainer::get_coeff_stick() const noexcept {
return coeff_stick;
}

void DLAContainer::set_coeff_stick(const double& _coeff_stick) {
if (_coeff_stick <= 0.0 || _coeff_stick > 1.0)
throw std::invalid_argument("_coeff_stick must be in (0,1]");
coeff_stick = _coeff_stick;
}

void DLAContainer::clear() {
aggregate_map.clear();
aggregate_pq = std::priority_queue<std::pair<int,int>, std::vector<std::pair<int,int>>, distance_comparator>();
}

void DLAContainer::generate(size_t _n) {
// initial seed set up
size_t count = 0;
std::pair<int,int> origin_sticky = std::make_pair(0,0);
aggregate_map.insert(std::make_pair(origin_sticky, count));
aggregate_pq.push(origin_sticky);
// initialise particle position vars
int x = 0;
int y = 0;
int prev_x = x;
int prev_y = y;
bool has_next_spawned = false;
// var to store current allowed size of bounding box spawn zone
int spawn_diameter = 0;
// uniform_distribution in [0,1] for probability generation
std::uniform_real_distribution<> dist(0.0, 1.0);

while (size() < _n) {
if (!has_next_spawned) {
spawn_particle(x,y,spawn_diameter,dist);
has_next_spawned = true;
}
double movement_choice = dist(mt_eng);
prev_x = x;
prev_y = y;
// update position of particle via unbiased random walk
update_particle_positions(x,y,movement_choice);
// check for collision with bounding walls
lattice_boundary_collision(x,y,prev_x,prev_y,spawn_diameter);
// generate stickiness probability
double sticky_pr = dist(mt_eng);
// check for collision with aggregate structure and add particle
// the aggregate if true
if (aggregate_collision(x,y,prev_x,prev_y,sticky_pr,count))
has_next_spawned = false;
}
}

void DLAContainer::spawn_particle(int& _x, int& _y, int& _spawn_diam, std::uniform_real_distribution<>& _dist) noexcept {
const int boundary_offset = 16;
// set diameter of spawn zone to double the maximum of the largest distance co-ordinate
// pair currently in the aggregate structure plus an offset to avoid direct sticking spawns
_spawn_diam = 2 * static_cast<int>(std::hypot(aggregate_pq.top().first, aggregate_pq.top().second)) + boundary_offset;

double placement_pr = _dist(mt_eng);
// spawn on upper line of lattice boundary
if (placement_pr < 0.25) {
_x = static_cast<int>(_spawn_diam*(_dist(mt_eng) - 0.5));
_y = _spawn_diam / 2;
}
// spawn on lower line of lattice boundary
else if (placement_pr >= 0.25 && placement_pr < 0.5) {
_x = static_cast<int>(_spawn_diam*(_dist(mt_eng) - 0.5));
_y = -_spawn_diam / 2;
}
// spawn on right line of lattice boundary
else if (placement_pr >= 0.5 && placement_pr < 0.75) {
_x = _spawn_diam / 2;
_y = static_cast<int>(_spawn_diam*(_dist(mt_eng) - 0.5));
}
// spawn on left line of lattice boundary
else {
_x = -_spawn_diam / 2;
_y = static_cast<int>(_spawn_diam*(_dist(mt_eng) - 0.5));
}
}

void DLAContainer::update_particle_position(int& _x, int& _y, const double& _movement_choice) const noexcept {
if (_movement_choice < 0.25)
++_x;
else if (_movement_choice >= 0.25 && _movement_choice < 0.5)
--_x;
else if (_movement_choice >= 0.5 && _movement_choice < 0.75)
++_y;
else if (_movement_choice >= 0.75 && _movement_choice < 1.0)
--_y;
}

bool DLAContainer::lattice_boundary_collision(int& _x, int& _y, const int& _prev_x, const int& _prev_y,
const int& _spawn_diam) const noexcept {
// small offset for boundary correction
const int epsilon = 2;
if (std::abs(_x) > ((_spawn_diam/2) + epsilon) ||
std::abs(_y) > ((_spawn_diam/2) + epsilon)) {
_x = _prev_x;
_y = _prev_y;
return true;
}
return false;
}

bool DLAContainer::aggregate_collision(const int& _x, const int& _y, const int& _prev_x,
const int& _prev_y, const double& _sticky_pr, size_t& _count) {
auto search = aggregate_map.find(std::make_pair(_x,_y));
// co-ordinates (_x,_y) occur in aggregate, collision occurred
if (search != aggregate_map.end() && sticky_pr <= coeff_stick) {

Currently I'm considering altering the private method signatures to take std::pair instances instead of passing lots of arguments to them which I'm doing at the moment, passing pairs containing int references may help to clean it up a little.