Program Listing for File GillespieCalculator.h¶
↰ Return to documentation for file (necsim/GillespieCalculator.h
)
//
// Created by sam on 07/09/19.
//
#ifndef NECSIM_GILLESPIECALCULATOR_H
#define NECSIM_GILLESPIECALCULATOR_H
#include <map>
#include <numeric>
#include <memory>
#include "Cell.h"
#include "RNGController.h"
#include "MapLocation.h"
namespace necsim
{
using namespace random_numbers;
enum class EventType
{
undefined, cell_event, map_event, sample_event
};
enum class CellEventType
{
undefined, dispersal_event, coalescence_event, speciation_event
};
class GillespieProbability
{
protected:
double dispersal_outside_cell_probability;
double coalescence_probability;
double speciation_probability;
double random_number;
MapLocation location;
public:
GillespieProbability() : GillespieProbability(MapLocation())
{ }
GillespieProbability(const MapLocation &c) : dispersal_outside_cell_probability(0.0),
coalescence_probability(0.0), speciation_probability(0.0),
random_number(0.0), location(c)
{
}
void setDispersalOutsideCellProbability(const double &d);
void setCoalescenceProbability(const double &c);
void setSpeciationProbability(const double &s);
void setRandomNumber(const double &r);
double getInCellProbability() const;
CellEventType generateRandomEvent(const shared_ptr<RNGController> &rng) const;
MapLocation &getMapLocation();
const MapLocation &getMapLocation() const;
double getLambda(const double &local_death_rate, const double &summed_death_rate, const unsigned long &n) const;
double calcTimeToNextEvent(const double &local_death_rate,
const double &mean_death_rate,
const unsigned long &n) const;
void reset();
friend ostream &operator<<(ostream &os, const GillespieProbability &gp);
friend std::istream &operator>>(std::istream &is, GillespieProbability &gp);
};
class GillespieHeapNode
{
private:
void updateHeapPosition()
{
// Enclosing heap vector not known or no locator to report to
if(heap == nullptr || locator == nullptr)
{
return;
}
// Node is outside its enclosing heap, i.e. it was probably moved to a temporary place
if(this < heap->data() || this >= (heap->data() + heap->size()))
{
return;
}
// Store the heap vector index with the locator
*locator = this - heap->data();
}
public:
Cell cell;
double time_of_event;
EventType event_type;
// Pointer to index in heap
vector<GillespieHeapNode>* heap;
unsigned long* locator;
GillespieHeapNode(const Cell cell,
const double time_of_event,
const EventType &e,
vector<GillespieHeapNode>* heap,
unsigned long* locator) : cell(cell), time_of_event(time_of_event), event_type(e), heap(heap),
locator(locator)
{
updateHeapPosition();
}
GillespieHeapNode(const double time_of_event, const EventType &e) : GillespieHeapNode(Cell(),
time_of_event,
e,
nullptr,
nullptr)
{ }
GillespieHeapNode() : GillespieHeapNode(Cell(), 0.0, EventType::undefined, nullptr, nullptr)
{ }
GillespieHeapNode(const GillespieHeapNode &other) noexcept
{
cell = other.cell;
time_of_event = other.time_of_event;
event_type = other.event_type;
heap = other.heap;
locator = other.locator;
updateHeapPosition();
}
GillespieHeapNode &operator=(const GillespieHeapNode &other)
{
cell = other.cell;
time_of_event = other.time_of_event;
event_type = other.event_type;
heap = other.heap;
locator = other.locator;
updateHeapPosition();
return *this;
}
GillespieHeapNode(GillespieHeapNode &&other) noexcept
{
cell = std::move(other.cell);
time_of_event = std::move(other.time_of_event);
event_type = std::move(other.event_type);
heap = std::move(other.heap);
locator = std::move(other.locator);
updateHeapPosition();
}
GillespieHeapNode &operator=(GillespieHeapNode &&other)
{
cell = std::move(other.cell);
time_of_event = std::move(other.time_of_event);
event_type = std::move(other.event_type);
heap = std::move(other.heap);
locator = std::move(other.locator);
updateHeapPosition();
return *this;
}
bool operator<(const GillespieHeapNode &n) const
{
return time_of_event > n.time_of_event;
}
};
}
#endif //NECSIM_GILLESPIECALCULATOR_H