// The class defined in this file codes the Gompertz growth models (non-cyclic and cyclic) described in the article.
// [Dennis and Ponciano (2014) - Density dependent state-space model for population abundance data with unequal time intervals, Eq. 3]

#ifndef GOMPERTZ_GROWTH_MODEL
#define GOMPERTZ_GROWTH_MODEL

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



#pragma mark -
#pragma mark Gompertz growth model
#pragma mark 



class GGDynamics{
private:
	double theta;		//equilibration speed (relaxation rate of log-transformed process)
	double K;			//carrying capacity
	double SD;			//long-term standard deviation around equilibrium resulting from additive white noise. Equal to beta^2K^2/(2theta) in [Dennis et al]
	double KPeriod; 	//Oscillation period of carrying capacity.
	double KAmplitude; 	//Oscillation amplitude of carrying capacity. Set to zero using setup() to obtain the non-cyclic model.
	
public:
	
	bool setup(	double 	_theta,
				double	_K,
				double 	_SD,
				double	_KPeriod,
				double	_KAmplitude){
		theta	= _theta;
		K 		= _K;
		SD 		= abs(_SD);
		KPeriod = _KPeriod;
		KAmplitude = _KAmplitude;
		return true;
	}
			
	
	bool operator()(double time, const XPoint &N, XPoint &velocity, XPoint &finalN) const{
		if(N.X<=0){
			velocity.X = 0;
		}else{
			velocity.X = theta*N.X*(log(K + KAmplitude*sin(time*2*PI/KPeriod)) - log(N.X));
		}
		return true;
	}
	
	double D(unsigned int k) const{
		switch(k){
		case 0: return theta*SQ(SD);
		}
		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);
		if(N.X<=0){
			ratesDown[0] = ratesUp[0] = 0;
		}else{
			ratesUp[0] 		= max(0.0, +theta*N.X*log(Ktemp)) + max(0.0, -theta*N.X*log(N.X));
			ratesDown[0] 	= max(0.0, -theta*N.X*log(Ktemp)) + max(0.0, +theta*N.X*log(N.X));
		}
		return true;
	}
	
	void writeDetailedLogDescription(ostream &streamOut, const string &linePrefix) const{
		streamOut 	<< linePrefix << "Equilibration rate theta = " << theta << "\n"
					<< linePrefix << "Carrying capacity K = " << K << (abs(KAmplitude)>EPSILON ? " + " + makeString(KAmplitude)+"*sin(t 2PI * "+makeString(1/KPeriod)+")" : ")") << "\n"
					<< linePrefix << "Stationary SD around equilibrium, sigma = " << SD << "\n";
	}
	
	string getShortFootnoteDescription() const{
		string s;
		s = "theta = " + makeString(theta) 
			+ ", K = " + makeString(K) + (abs(KAmplitude)>EPSILON ? " + " + makeString(KAmplitude)+"*sin(t 2pi * "+makeString(1/KPeriod)+")" : "")
			+ ", sigma = " + makeString(SD);
		return s;
	}


	bool deterministic() const{ return (abs(SD)<EPSILON); }
	double getKAmplitude() const{ return KAmplitude; }
	double getKPeriod() const{ return KPeriod; }
	bool isCyclic() const{ return abs(KAmplitude)>EPSILON; }
};



			

#endif