// The class defined in this file codes the mathematical models (non-cyclic and cyclic) described in the article.

#ifndef LOGISTIC_GROWTH_MODEL
#define LOGISTIC_GROWTH_MODEL

#include "Points.h"
#include <cmath>
#include <iostream>
using namespace std;



#pragma mark -
#pragma mark Logistic growth model (Verhulst)
#pragma mark 



class LGDynamics{
private:
	double 	r;	//intrinsic growth rate
	double 	K;	//carrying capacity
	double 	SD;	//long-term standard deviation around equilibrium resulting from additive white noise
	double 	KPeriod; //Oscillation period of carrying capacity.
	double 	KAmplitude; //Oscillation amplitude of carrying capacity
	double 	cycleAmplitude;
	bool 	scaleNoiseSDWithPopSize;
	double 	cyclePhaseLag;	//Cycle phase lage w.r.t. forcing (approximation for small-amplitude noise and cycles)
	
public:
	
	bool setup(	double 	_r,
				double	_K,
				double 	_SD,
				double	_cyclePeriod,
				double	_cycleAmplitude,
				bool	_scaleNoiseSDWithPopSize){
		r = _r;
		K = _K;
		SD = abs(_SD);
		scaleNoiseSDWithPopSize = _scaleNoiseSDWithPopSize;
		KPeriod = _cyclePeriod;
		if((_cycleAmplitude<K*EPSILON) || (_cyclePeriod*r<EPSILON)){ 
			KAmplitude = cycleAmplitude = 0;
		}else{
			cycleAmplitude = _cycleAmplitude;
			KAmplitude = sqrt((SQ(r)+SQ(2*PI/KPeriod))/SQ(r)) * cycleAmplitude; // the system's finite relaxation rate dampens the forcing and reduces the actual cycle amplitude
			cyclePhaseLag = atan2(sqrt(SQ(2*PI*r)/(SQ(2*PI*r)+SQ(r))), sqrt(SQ(r)/(SQ(2*PI*r)+SQ(r))));
		}
		if(KAmplitude>=K) return false;
		return true;
	}
			
	
	bool operator()(double time, const XPoint &N, XPoint &velocity, XPoint &finalN) const{
		velocity.X = r*N.X*(1.0 - N.X/(K + KAmplitude*sin(time*2*PI/KPeriod)));
		return true;
	}
	
	double D(unsigned int k) const{
		switch(k){
		case 0: return r*SQ(SD);
		}
		return 0;
	}
	
	double D(unsigned int k, const XPoint &N) const{
		switch(k){
		case 0: return r*SQ(SD*(scaleNoiseSDWithPopSize ? N.X/K : 1));
		}
		return 0;
	}
	
	bool rates(double time, const XPoint &N, double ratesUp[], double ratesDown[], XPoint &finalN) const{
		double Ktemp	= K + KAmplitude*sin(time*2*PI/KPeriod);
		ratesUp[0] 		= max(0.0, +r*N.X) + max(0.0, -r*SQ(N.X)/Ktemp);
		ratesDown[0] 	= max(0.0, -r*N.X) + max(0.0, +r*SQ(N.X)/Ktemp);
		return true;
	}
	
	void writeDetailedLogDescription(ostream &streamOut, const string &linePrefix) const{
		streamOut 	<< linePrefix << "Growth rate r = " << r << "\n"
					<< linePrefix << "Carrying capacity K = " << K << (abs(KAmplitude)>EPSILON ? " + " + makeString(KAmplitude)+"*sin(t 2PI * "+makeString(1/KPeriod)+")" : ")") << "\n"
					<< linePrefix << "Cycle amplitude (approx), alpha = " << cycleAmplitude << "\n"
					<< linePrefix << "Stationary SD around equilibrium sigma = " << SD << (scaleNoiseSDWithPopSize ? " (but noise amplitude scales with N)" : "") << "\n";
	}
	
	string getShortFootnoteDescription() const{
		string s;
		s = "r = " + makeString(r) 
			+ ", K = " + makeString(K) + (abs(KAmplitude)>EPSILON ? " + " + makeString(KAmplitude)+"*sin(t 2pi * "+makeString(1/KPeriod)+")" : "")
			+ ", alpha = " + makeString(cycleAmplitude)
			+ ", sigma = " + makeString(SD) + (scaleNoiseSDWithPopSize ? " (but noise amplitude scales with N)" : "");
		return s;
	}


	bool deterministic() const{ return (abs(SD)<EPSILON); }
	double getKAmplitude() const{ return KAmplitude; }
	double getCycleAmplitude() const{ return cycleAmplitude; }
	double getKPeriod() const{ return KPeriod; }
	bool isCyclic() const{ return abs(KAmplitude)>EPSILON; }
	double getRelaxationRate() const{ return r; } // relaxation rate near the carrying capacity
	double getPowerAtZero() const{ return SQ(SD)*2.0/r; } // power spectrum at zero-frequency (approximating the model by an OU process around the carrying capacity)
	double getPowerAtZero(double dt) const{ return getPowerAtZero() * dt*(r/2) * (1+exp(-r*dt))/(1-exp(-r*dt)); } // power spectrum at zero-frequency, corrected for discrete time series
	double getSD() const{ return SD; }
	XPoint getRandomValue(double time) const{ XPoint p; p.X = K + cycleAmplitude*sin(time*2*PI/KPeriod - cyclePhaseLag) + standardNormal()*SD; return p; } // returns random starting point in stationary case
};



			

#endif