// This header describes an object used to do Gillespie
// simulations.
//
# ifndef _GillespieSimulation_ // There can be only one
# define _GillespieSimulation_
// Standard headers
# include <vector>
# include <map>
# include <iostream> // cout, etc.
// My headers
# include "Reaction.h"
# include <chrono>
#include <random>
# include "ziggurat_inline.hpp"
# include "TwoDsearch.h"
typedef TwoDsearch Search_obj;//use 2D search (rates are put in a matrix of size sqrt(n_react))
class GillespieSimulation
{
public:
GillespieSimulation(
const ReactionsArray &RX,Species &P) ;//need Species to pass to the chosen Reaction->do_reaction. Other use should only be via Species member functions
// void get_populations( Polymer &P) const
// {
// P = this->P ;
// }
double advance_time( double time_step ) ;
double advance_step( size_t n_steps ) ;
private:
ReactionsArray reactions ; // list of reactions
Species& P ; // populations
double SimTime;//simulation time//unused
size_t n_react;//instead of reactions.size()
std::vector<double> rate ; // the rates, only for temporary storage. rates are stored in the seacrh obj
double rate_sum ; // sum of rates
std::shared_ptr<Search_obj> search_obj; //binary tree or 2D search object
//random number generation
unsigned long int seed;
std::mt19937_64 gen;//for uniform random numbers to choose which reaction occurs. 64bit to ensure enough significant digits to select reactions with very small rates see Mauch & Stalzer
std::uniform_real_distribution<double> unidist;
uint32_t jsr_value;// for the ziggurat rng. used for exponentially distributed random numbers
uint32_t jcong_value;
uint32_t w_value;
uint32_t z_value;
size_t n_comp;
// Private functions for the dynamics
double time_to_next_event( void ) ;
std::tuple<size_t,double,double> choose_next_reaction( void ) ;//return the rx number, scaled random number and the rate (propensity)