2
\$\begingroup\$

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.

DLAContainer Class Header File

#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) {
        std::pair<int,int> added_particle = std::make_pair(_prev_x,_prev_y);
        // insert previous traversed position to aggregate structures
        aggregate_map.insert(std::make_pair(added_particle,++_count));
        aggregate_pq.push(added_particle);
        return true;
    }
    return false;
}

This code currently runs pretty damn fast but I'd like to get as much performance out of it as possible therefore any further optimisation suggestions are welcome — also any non-optimisation related suggestions would be useful if you spot something that should be changed.

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.

\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.