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

#ifndef CLASSICAL_OU_MODEL
#define CLASSICAL_OU_MODEL

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



#pragma mark -
#pragma mark Classical (1D) Ornstein-Uhlenbeck process
#pragma mark 



class OUDynamics{
private:
	double r;				//relaxation rate
	double K;				//mean value ('equilibrium')
	double SD;				//long-term standard deviation around equilibrium
	double KPeriod; 		//Oscillation period of equilibrium.
	double KAmplitude; 		//Oscillation amplitude of equilibrium.
	double cycleAmplitude;	//Actual deterministic oscillation amplitude (smaller than Kamplitude due to finite relaxation rate)
	double cyclePhaseLag;	//Cycle phase lage w.r.t. forcing

public:
	
	bool setup(	double 	_r,
				double	_K,
				double 	_SD,
				double	_cyclePeriod,
				double	_cycleAmplitude){
		r = _r;
		K = _K;
		SD = abs(_SD);
		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))));
		}
		return true;
	}
	
	
	double getDeterministicTrajectory(double time) const{
		return K + cycleAmplitude*sin(time*2*PI/KPeriod - cyclePhaseLag);
	}
			
	
	bool operator()(double time, const XPoint &N, XPoint &velocity, XPoint &finalN) const{
		velocity.X = -r*(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;
	}
	
	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 - r*Ktemp));
		ratesDown[0] 	= max(0.0, +r*(N.X - r*Ktemp));
		return true;
	}
	
	void writeDetailedLogDescription(ostream &streamOut, const string &linePrefix) const{
		streamOut 	<< linePrefix << "Relaxation rate, r = " << r << "\n"
					<< linePrefix << "Equilibrium, K = " << K << (abs(KAmplitude)>EPSILON ? " + " + makeString(KAmplitude)+"*sin(t 2PI * "+makeString(1/KPeriod)+")" : ")") << "\n"
					<< linePrefix << "Cycle amplitude, alpha = " << cycleAmplitude << "\n"
					<< linePrefix << "Stationary SD around equilibrium, sigma = " << SD << "\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);
		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; }
	double getPowerAtZero() const{ return SQ(SD)*2.0/r; } // power spectrum at zero-frequency
	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 = getDeterministicTrajectory(time) + standardNormal()*SD; return p; } // returns random starting point in stationary case
};



			

#endif