#include <random>
#include <stdio.h>
#include <omp.h>

using namespace std;

class latticeinfo
{
public:
  latticeinfo(vector<int> dim, double new_s, char new_obs, char new_process, vector<double> rates, mt19937* new_twister);
        ~latticeinfo();
        int obs1;
	int obs2;
	int obs3;
        int N; /*Holds number of particles*/
        int L; /*Holds number of sites*/
	mt19937 *twister; /*Pseudo-random number generator*/
        double a;
        int q;
        std::uniform_real_distribution<double> cdist{0, 1};
        std::uniform_int_distribution<int> bdist{0, 1};
        std::uniform_int_distribution<int> pdist{0, 1};
        std::uniform_int_distribution<int> sdist{0, 1};
        double intesc; /*Integral of Escape Rate*/
	double intscesc; /*Integral of Scaled Escape Rate*/
	double exps; /*Holds exp(-s)*/
	char proc; /*Process*/
	double rate1; /*Process Rates*/
	double rate2;
	char obs; /*Observable we are biasing*/
	double tscale; /*Timescale to use when scaling*/
	double r1scale; /*Scaled rates*/
	double r2scale;
        double randc(); /*Returns a continuous random number on the range 0 to 1*/
        int randb(); /*Returns a random binary number*/
        int randp(); /*Returns a random particle index*/
        int rands(); /*Returns a random site index*/
        void init(lattice *Lat);
        int boundconds(int val);
        void sep(double tmax, lattice *lat);
        void asep(double tmax, lattice *lat);
        void cp(double tmax, lattice *lat);
        void upsh(double *upsilon); //Calculates Cloning Factor when Observing Hops
        void upsc(double *upsilon); //Calculates Cloning Factor when Observing Current
        void (latticeinfo::*dynamics)(double, lattice*); //Pointer to process to run
        void (latticeinfo::*calcups)(double*); //Pointer to function that calculates cloning factor Y
        int size; //Size in bytes of one lattice
};

latticeinfo::latticeinfo(vector<int> dim, double new_s, char new_obs, char new_process, vector<double> rates, mt19937* new_twister) {
	obs1 = 0;
	obs2 = 0;
	obs3 = 0;
	N = dim[1];
	L = dim[0];
	rate1 = rates[0];
	rate2 = rates[1];
	twister = new_twister;
	decltype(pdist.param()) p_range(0, N - 1);
	pdist.param(p_range); /*Set range of random particles*/
	decltype(sdist.param()) s_range(0, L - 1);
	sdist.param(s_range); /*Set range of random sites*/
	intesc = 0;
	intscesc = 0;
	exps = exp(-1 * new_s);
	obs = new_obs;
	proc = new_process;
	size = 2*N*sizeof(int) + L*sizeof(int) + sizeof(double); //Size of lattice components in bytes
	if (proc == 'S')
	{
	  dynamics = &latticeinfo::sep;
	}
	if (proc == 'A')
	{
	  dynamics = &latticeinfo::asep;
	}
	if (proc == 'C')
	{
	  dynamics = &latticeinfo::cp;
	}
	if (obs == 'h')
	{
		tscale = 1 / ((rate1 + rate2)*exp(-1 * new_s));
		r1scale = rate1*exp(-1 * new_s);
		r2scale = rate2*exp(-1 * new_s);
		calcups = &latticeinfo::upsh;
	}
	else if (obs == 'c')
	{
		tscale = 1 / (rate1*exp(new_s) + rate2*exp(-1 * new_s));
		r1scale = rate1*exp(new_s);
		r2scale = rate2*exp(-1 * new_s);
		calcups = &latticeinfo::upsc;
	}
	else
	{
		tscale = 1;
		r1scale = rate1;
		r2scale = rate2;
	}
}

double latticeinfo::randc()
{
  a = cdist(*twister);
  return a;
}

int latticeinfo::randb()
{
  q = bdist(*twister);
  return q;
}

int latticeinfo::randp()
{
  q = pdist(*twister);
  return q;
}

int latticeinfo::rands()
{
  q = sdist(*twister);
  return q;
}

void latticeinfo::init(lattice *Lat)
{
  int k, csite, lsite, rsite, neighbours;
  int ip = rands(); /*Random inital particle position*/
  Lat->esc=0; /*Set escape rate to 0*/
  for(k=0; k<N; k++)
  {
    while(Lat->s[ip] > -1) /*Until a valid position found*/
    {
      ip = rands(); /*Randomnly pick initial site*/
    }
    Lat->p[k]=ip; /*Initialise each particle position*/
    Lat->s[ip]=k;
  }
  for(k=0; k<N; k++) /*Calculate initial freedom of each particle*/
  {
    csite = Lat->p[k]; /*Site that this particle is in*/
    lsite = csite - 1;
    lsite = boundconds(lsite);
    rsite = csite + 1;
    rsite = boundconds(rsite);
    neighbours = (Lat->s[lsite]>-1) + (Lat->s[rsite]>-1); /*Calculate number of neighbours that the particle has*/
    Lat->freedom[k] = 2 - neighbours; /*Calculate particle freedom*/
    Lat->esc += 0.5*(Lat->freedom[k]);
  }
    if(proc == 'R')
    {
      Lat->esc = 2 * Lat->esc;
    }
    if(proc == 'A')
    {
      Lat->esc = (rate1+rate2)*Lat->esc;
    }
}

int latticeinfo::boundconds(int val){ /*Applies periodic boundary conditions*/
  if (val == -1) /*Left-hand Boundary Condition*/
  {
    val = L - 1;
  }
  if (val == L) /*Right-hand Boundary Condition*/
  {
    val = 0;
  }
  return val;
}

void latticeinfo::sep(double tmax, lattice* lat) /*Symmetric Exclusion Process*/
{
	int i, d, k;
	double t, r1;
	int r2, ns, np; /*Neighbouring site and neighbouring particle*/
	int as, ts; /*Particle being moved, moves away from as and towards ts*/
	int ap, tp; /*Particle being moved, moves away from ap and towards tp*/
	int hops = 0; /*Counts number of hops*/
	int current = 0; /*Records Current*/
	intesc = 0; /*Reset local time and integral of escape rates*/
	t = 0;
	double dt;
	r1 = randc(); /*Random number to determine size of time step*/
	dt = (-1. / (N))*log(r1)*tscale; /*Generate time step*/
	t += dt; /*Update time*/
	intesc += (lat->esc)*dt;/*Integrate escape rate*/
	if (t < tmax)
	{
		while (t < tmax) /*So as not to exceed time*/
		{
			r2 = randc()*(r1scale + r2scale); /*Random number to determine direction*/
			i = randp(); /*Random particle to move*/
			if (r2<r1scale)
			{
				d = -1; /*Direction to move particle*/
			}
			else
			{
				d = 1;
			}
			ns = lat->p[i] + d; /*Neighbouring site*/
			ns = boundconds(ns); /*Apply boundary conditions*/
			np = lat->s[ns]; /*Neighbouring particle*/

			ts = ns + d; /*Site moving towards*/
			ts = boundconds(ts);
			tp = lat->s[ts];

			as = lat->p[i] - d; /*Site moving away from*/
			as = boundconds(as);
			ap = lat->s[as];

			if (np == -1) /*If neighbouring site unoccupied*/
			{
				lat->s[lat->p[i]] = -1; /*Make current site empty*/
				if (ap>-1) /*If there is a particle in site being moved away from*/
				{
					lat->freedom[ap] += 1; /*Increase freedom of both particles*/
					lat->freedom[i] += 1;
					lat->esc += 1; /*Increase escape rate by 0.5 for each new possible hop*/
				}
				lat->s[ns] = i; /*Occupy new site*/
				lat->p[i] = ns; /*Update particle position*/

				if (tp>-1) /*If there is a particle in site being moved towards*/
				{
					lat->freedom[tp] -= 1; /*Decrease freedom of both*/
					lat->freedom[i] -= 1;
					lat->esc -= 1; /*Decrease escape rate by 0.5 for each hop now impossible*/
				}
				/*Apply boundary condition*/
				hops += 1; /*Update Number of Hops*/
				current += d; /*Update Current*/
			}
			r1 = randc(); /*Random number to determine size of time step*/
			dt = (-1. / N)*log(r1)*tscale; /*Generate time step*/
			intesc += dt*(lat->esc); /*Integrate Time Step*/
			t += dt; /*Update time*/
		}
	}
	intesc -= (t - tmax)*(lat->esc); /*Final time step overshoots max time*/
	/*For correct integral value remove the integral of the scape rate over this extra time*/
	obs1 = hops;
	obs2 = current;
};

void latticeinfo::asep(double tmax, lattice *lat)
{
    
    double etascale = r1scale+r2scale;
    double eta = rate1 + rate2;
    int i, d, k;
    double t, r1, r2;
    int ns, np; /*Neighbouring site and neighbouring particle*/
    int as, ts; /*Particle being moved, moves away from site as and towards site ts*/
    int ap, tp; /*Particle being moved, moves away from particle ap and towards particle tp*/
    int hops = 0; /*Counts number of hops*/
    int current = 0; /*Records Current*/
    intesc = 0; /*Reset local time and integral of escape rates*/
    t = 0;
    double dt;
    r1 = randc(); /*Random number to determine size of time step*/
    dt = (-1. / (etascale*N))*log(r1); /*Generate time step*/
    t += dt; /*Update time*/
    intesc += (lat->esc)*dt;
    if (t < tmax)
    {
        while (t < tmax) /*So as not to exceed time*/
        {
            r2 = (r1scale+r2scale)*randc(); /*Random number to determine direction*/
            i = randp();/*Particle to move*/
            if (r2<r1scale) /*Direction to move particle*/
            {
                d = -1;
            }
            else
            {
                d = 1;
            }
            ns = lat->p[i] + d; /*Neighbouring site*/
            ns = boundconds(ns); /*Apply boundary conditions*/
            np = lat->s[ns]; /*Neighbouring particle*/
            
            ts = ns + d; /*Site moving towards*/
            ts = boundconds(ts);
            tp = lat->s[ts];
            
            as = lat->p[i] - d; /*Site moving away from*/
            as = boundconds(as);
            ap = lat->s[as];
            
            /*Apply periodic boundary conditions*/
            if (np == -1) /*If neighbouring site unoccupied*/
            {
                lat->s[lat->p[i]] = -1; /*Make current site empty*/
                if(ap>-1) /*If there is a particle in site being moved away from*/
                {
                    lat->freedom[ap] += 1; /*Increase freedom of both particles*/
                    lat->freedom[i] += 1;
                    lat->esc += eta; /*One new hop increases escape rate by p, the other by q*/
                }
                lat->s[ns] = i; /*Occupy new site*/
                lat->p[i] = ns; /*Update particle position*/
                if(tp>-1) /*If there is a particle in site being moved towards*/
                {
                    lat->freedom[tp] -= 1; /*Decrease freedom of both*/
                    lat->freedom[i] -= 1;
                    lat->esc -= eta; /*One hop inhibition decreases escape rate by p, the other by q*/
                }
                hops += 1; /*Update Number of Hops*/
                current += d; /*Update Current*/
            }
            r1 = randc(); /*Random number to determine size of time step*/
            dt = (-1. / (etascale*N))*log(r1); /*Generate time step*/
            intesc+=dt*(lat->esc); /*Integrate time step*/
            t += dt; /*Update time*/
        }
    }
    intesc-=(t-tmax)*(lat->esc); /*Final time step overshoots max time*/
    /*For correct integral value remove the integral of the escape rate over this extra time*/
    obs1 = hops;
    obs2 = current;
};

void latticeinfo::cp(double tmax, lattice *lat)
{
  double dt; /*Size of time step*/
  int i, nl ,nr;
  double r1, r2, r3;
  double t = 0;
  while (t < tmax) /*At each time step*/
  {
    i = rands(); /*Random number to determine site*/
    nl = i - 1; /*Nearest site to left*/
    nr = i + 1; /*Nearest site to right*/
    nl = boundconds(nl);
    nr = boundconds(nr);
    r2 = randc(); /*Random number to determine annihilation/not*/
    r3 = randc(); /*Random number to determine new occupation/not*/
    if (lat->s[i] > 0) /*If site occupied*/
    {
      if (r2 < 1. / (2 * rate1 + rate2 + 1)) /*With rate 1*/
      {
	lat->s[i] = -1; /*Annihilate*/
      }
    }
    else /*If site unoccupied*/
    {
      if (r3 < (rate1*((lat->s[nl]>0) + (lat->s[nr]>0)) + rate2) / (2 * rate1 + rate2 + 1))
      {
	lat->s[i] = 1; /*Occupy*/
      }
    }
    r1 = randc(); /*Random number to determine size of time step*/
    dt = (-1. / (L*(2 * rate1 + rate2 + 1)))*log(r1); /*Generate time step*/
    t += dt;
  }
};

void latticeinfo::upsh(double *upsilon){
  *upsilon = exp((exps - 1)*intesc);
}

void latticeinfo::upsc(double *upsilon){
  intscesc = (intesc / (rate1+rate2))*(rate1 / exps + rate2 * exps);
  *upsilon = exp(intscesc - intesc);
}

latticeinfo::~latticeinfo(){ /*Destructor*/
}
