#ifndef SIGNIFICANCE_SPECTRAL_PEAK
#define SIGNIFICANCE_SPECTRAL_PEAK



#include <cmath>
#include "Config.h"
#include "Auxiliary.cpp"
#include "InternalDefs.h"
#include "FitOU.cpp"
#include "FitOUSS.cpp"
#include "VectorArithmetics.h"
#include "LombScargleSpectrum.h"
#ifdef COMPILE_WITH_OMP
	#include <omp.h>
#endif
using namespace std;






// fit OUSS model to periodogram
// include correction for discrete & finite time series
bool fit_DFTS_OUSS_to_periodogram(	double 					df,
									double					dt,
									long					N,
									const vector<double> 	&spectrum,	// assumes that normalization used is the conventional for continuous-time stochastic-processes
									const double 			true_lambda,	// for debugging
									const double			true_power_o,	// for debugging
									const double			true_power_e,	// for debugging
									bool					spectrumIsAveraged,
									double					&fitted_lambda,
									double					&fitted_power_o,
									double					&fitted_power_e){
	
	
	
	// debug
//	fitted_lambda = true_lambda;
//	fitted_power_o = true_power_o;
//	fitted_power_e = true_power_e;
//	return true;
	
	
	//guess reasonable parameter values (use as starting values for fitting)
	const double suspected_power_e	= average(spectrum,max(1l,(long)spectrum.size()-5),spectrum.size()-1);
	const double powerInt 			= sum(spectrum,1,spectrum.size()-1)*df;
	const double suspected_power_o  = average(spectrum,1,min((long)spectrum.size()-1,5l));
	const double suspected_lambda 	= 4*powerInt/suspected_power_o;
				
	if(PS_OUSS_FIT_LOGLIKELIHOOD && (!spectrumIsAveraged)){
		// fit OUSS power spectrum to periodogram using maximum-likelihood
		if(!fit_PS_ML_OUSS(	df,
							dt,
							N,
							spectrum,
							suspected_lambda, // debug true_lambda,
							suspected_power_o, // debug true_power_o,
							suspected_power_e, // debug true_power_e,
							fitted_lambda,
							fitted_power_o,
							fitted_power_e)){
			return false;
		}
		

	}else if(!fit_PS_LS_OUSS(	df,
								dt,
								N,
								spectrum,
								suspected_lambda, // debug true_lambda,
								suspected_power_o,// debug true_power_o,
								suspected_power_e, // debug true_power_e,
								fitted_lambda,
								fitted_power_o,
								fitted_power_e)){
		return false;
	}
	return true;
}





//Monte Carlo estimation of significance of spectral peak, assuming an Ornstein-Uhlenbeck process + normally distributed errors (= OU state-space model)
//Only accurate for periodograms obtained from long, evenly sampled time series 
double significancePeakOUSS(double 	lambda,
							double 	power_o, 
							double	power_e, // power generated by errors. can be zero if no errors (i.e. classical OU).
							double 	df, 
							long 	maxMode, 
							long	minMode, // min mode below which periodograms are not even considered. If you dismiss periodograms with a peak below this mode, you need to evaluate P-value conditional upon a non-dismissal.
							double 	peakFreq,
							double 	peakPower){
	long n, m, countPos=0, countEvaluations=0;
	double SQlambda = SQ(lambda);
	double EpeakPower = power_e + power_o*SQlambda/(SQ(2*PI*peakFreq)+SQlambda);
	double r_peakPower, r_EpeakPower, r_peakMode, power, Epower;
	vector<double> Epowers(maxMode+1,0);
	for(m=0; m<=maxMode; ++m){ Epowers[m] = power_e + power_o*SQlambda/(SQ(2*PI*m*df)+SQlambda); }
	for(n=0; n<SIGNIF_BERNOULLI_TRIALS; ++n){
		for(m=1; m<=maxMode; ++m){
			power 	= exponential(Epowers[m]);
			if((r_peakPower<power) || (m==1)){
				r_peakPower 	= power;
				r_EpeakPower	= Epowers[m];
				r_peakMode		= m;
			}
		}
		if(r_peakMode<minMode) continue; // dismiss this trial
		++countEvaluations;
//		if((r_peakPower>=peakPower) && (r_peakPower/r_EpeakPower >= peakPower/EpeakPower)){
		if(SQ(r_peakPower)/r_EpeakPower >= SQ(peakPower)/EpeakPower){
			++countPos;
		}
	}
	return (1.0*countPos)/max(1l,countEvaluations);
}





//Monte Carlo estimation of significance of spectral peak, assuming an Ornstein-Uhlenbeck process + normally distributed errors (= OU state-space model)
// Only accurate for periodograms obtained from long, evenly sampled time series 
// But tries to correct expected periodogram power for discrete finite time series
double significancePeakOUSS(const double 	lambda,
							const double 	power_o, 
							const double	power_e, 	// power generated by errors. can be zero if no errors (i.e. classical OU).
							const double 	df, 
							const long 		maxMode, 
							const long		minMode, 	// min mode below which periodograms are not even considered. If you dismiss periodograms with a peak below this mode, you need to evaluate P-value conditional upon a non-dismissal.
							const double	dt,			// time step used for the time series. This is needed for a correction of the expected periodogram.
							const long		N,			// time series length. This is needed for a correction the the expected periodogram.
							const double 	peakFreq,
							const double 	peakPower){
	long n, m, countPos=0, countEvaluations=0;
	double EpeakPower = PS_DFTS_OUSS(lambda, power_o, power_e, dt, N, peakFreq);
	double r_peakPower, r_EpeakPower, r_peakMode, power, Epower;
	vector<double> Epowers(maxMode+1,0);
	for(m=0; m<=maxMode; ++m){ Epowers[m] = PS_DFTS_OUSS(lambda, power_o, power_e, dt, N, m*df); }
	
	#pragma omp parallel for schedule(static) private(n,m,power,r_peakPower,r_EpeakPower,r_peakMode) reduction(+:countPos,countEvaluations) 
	for(long n=0; n<SIGNIF_BERNOULLI_TRIALS; ++n){
		for(m=1; m<=maxMode; ++m){
			power = exponential(Epowers[m]);
			if((r_peakPower<power) || (m==1)){
				r_peakPower 	= power;
				r_EpeakPower	= Epowers[m];
				r_peakMode		= m;
			}
		}
		if(r_peakMode<minMode) continue; // dismiss this trial
		++countEvaluations;
//		if((r_peakPower>=peakPower) && (r_peakPower/r_EpeakPower >= peakPower/EpeakPower)){
		if(SQ(r_peakPower)/r_EpeakPower >= SQ(peakPower)/EpeakPower){
			++countPos;
		}
	}
	return (1.0*countPos)/max(1l,countEvaluations);
}



#pragma mark -
#pragma mark Simulation-based P-value estimation
#pragma mark 




// generate time series at times 0, dt, 2dt...
// of stationary OUSS process Y, defined as:
//		dX = -lambda*X*dt  +  sqrt(2*sigma^2)*dW,   Y = X + epsilon*N,   N:standard-normal
void generateOUSStimeSeries(double lambda,
							double sigma,
							double epsilon,
							double dt,
							double count,
							vector<double> &signal){
	signal.resize(count);
	long n;
	double rho = exp(-lambda*dt);
	//random start point
	signal[0] = standardNormal()*sigma;
	//generate correlated time series
	for(n=1; n<count; ++n){
		signal[n] = rho*signal[n-1] + sqrt(1-SQ(rho))*standardNormal()*sigma;
	}
	//add normal error terms to OU signal
	for(n=0; n<count; ++n){
		signal[n] += epsilon * standardNormal();
	}
}



void generateOUSStimeSeries(double 					lambda,
							double 					sigma,
							double 					epsilon,
							const vector<double> 	&times,
							vector<double> 			&signal){
	const long N = times.size();
	signal.resize(N);
	long n;
	double rho;
	//random start point
	signal[0] = standardNormal()*sigma;
	//generate correlated time series
	for(n=1; n<N; ++n){
		rho = exp(-lambda*(times[n]-times[n-1]));
		signal[n] = rho*signal[n-1] + sqrt(1-SQ(rho))*standardNormal()*sigma;
	}
	//add normal error terms to OU signal
	for(n=0; n<N; ++n){
		signal[n] += epsilon * standardNormal();
	}
}


//Calculate periodogram using classical DFT formula
bool DFTPeriodogram(double dt,
					const vector<double> &signal,
					vector<double> &periodogram,
					double &df){
	long N = signal.size();
	if(N<2) return false;
	double T = dt*(N-1); //data time span
	long n,m;
	df = 1.0/(N*dt);
	long maxMode = N/2-1; //upper limit given by Nyquist frequency
	periodogram.resize(1+maxMode);
	double FRe, FIm;
	
	for(m=0; m<=maxMode; ++m){
		for(n=0, FRe=FIm=0; n<N; ++n){
			FRe += signal[n]*cos(-(2.0*PI*m*n)/N);
			FIm += signal[n]*sin(-(2.0*PI*m*n)/N);
		}
		FRe = FRe/N;
		FIm = FIm/N;
		
		periodogram[m] = T*(FRe*FRe + FIm*FIm);
	}
	return true;
}



//UNTESTED!!
//Similar to significancePeakOUSS(..) above, but here OU time series is generated and from that a periodogram is calculated
//This estimator is much slower than the one above, but expected to be more accurate, especially for short time series
double significancePeakOUSS_simul(	const double 			lambda,
									const double 			power_o, 
									const double 			power_e, 	// power generated by errors. can be zero if no errors (i.e. classical OU).
									const double 			dt,			// time step used for the time series from which periodogram peak was calculated
									const long 	 			N,			// size of time series
									const vector<double> 	&times,
									const long 				maxMode, 	// maximum mode below which periodogram peak was calculated
									const long				minMode, 	// min mode below which periodograms are not even considered.
									const double 			peakFreq,
									const double 			peakPower){
	vector<double> OUsignal, periodogram;
	const double T = (N-1)*dt;
	const double EpeakPower = PS_DFTS_OUSS(lambda, power_o, power_e, dt, N, peakFreq);
	double r_peakPower, r_EpeakPower, power, Epower;
	long r_peakMode;
	// figure out OU process parameters
	const double epsilon = sqrt(power_e/dt);
	const double rho = exp(-lambda*dt);
	const double sigma = sqrt((power_o/dt)*(1-rho)/(1+rho)); // standard deviation of stationary OU process
	double df;
	long countPos=0, countEvaluations=0;
	
	// debug
//	double average = 0;
//	double last_df;
//	cout << "  minMode=" << minMode << ", maxMode=" << maxMode << ", peakFreq=" << peakFreq << endl;
	
	#pragma omp parallel for private(r_peakPower,r_EpeakPower,r_peakMode,OUsignal,periodogram,df) schedule(static) reduction(+:countPos,countEvaluations)
	// run Monte Carlo
	for(long n=0; n<SIGNIF_BERNOULLI_TRIALS; ++n){
		// generate OU time series
		generateOUSStimeSeries(	lambda,
								sigma,
								epsilon,
								dt,
								N,
								OUsignal);
		
		// calculate periodogram
		DFTPeriodogram(dt, OUsignal, periodogram, df);
		
		// debug: use LS periodogram
/*		LombScargleSpectrum<double>(OUsignal, 
									times,
									0, 
									N-1,
									0,
									NULL,
									NULL,
									&periodogram,
									df); 
		*/
		bool OK = true;
		
		// DEPRECATED
/*		if(P_VALUE_USING_SIMULATION_REFIT_EVERY_TRIALS){
			#pragma omp critical
			// the fitting routines are not thread safe because they use global variables
			{
				OK = OK && fit_DFTS_OUSS_to_periodogram(df,
														dt,
														N,
														periodogram,	
														trial_lambda,
														trial_power_o,
														trial_power_e);
			}
		}
*/
		if(!OK) continue;
//		average = (n==0 ? 0 : average) + periodogram[10]; last_df = df; // debug
//		if(n==0) cout << "  power[" << df << "] = " << periodogram[1] << endl;
		
		// find mode (maximum) of generated periodogram
		r_peakMode = -1;
		for(long m=minMode; m<periodogram.size(); ++m){
			if((r_peakPower<periodogram[m]) || (r_peakMode<0)){
				r_peakPower 	= periodogram[m];
				r_EpeakPower	= PS_DFTS_OUSS(lambda, power_o, power_e, dt, N, m*df);
				r_peakMode		= m;
			}
		}
		++countEvaluations;
//		if((r_peakPower>=peakPower) && (r_peakPower/r_EpeakPower >= peakPower/EpeakPower)){
		if(SQ(r_peakPower)/r_EpeakPower >= SQ(peakPower)/EpeakPower){
			++countPos;
		}
	}
	
	// debug
//	average /= SIGNIF_BERNOULLI_TRIALS;
//	cout << "    average[50] = " << average << ", vs expected = " << PS_DFTS_OUSS(lambda, power_o, power_e, dt, N, 10*last_df) << ", vs asymptotic = " << PS_DTS_OUSS(lambda, power_o, power_e, dt, 50*last_df) << endl;
	return (1.0*countPos)/max(1l,countEvaluations);
}







#pragma mark -
#pragma mark Wrappers for significance testing
#pragma mark

bool significanceOfSpectralPeak_OU( double 					df,
									double					T,
									const vector<double> 	&spectrum, // normalization used for periodogram is irrelevant
									long					minMode, //dismiss (i.e. return false) if periodogram peak is at mode below minMode. If you already filtered your samples by such a criterion, you still need to pass this value so that P-values are correctly estimated.
									double					&peakf,
									double					&peakPower,
									double					&significanceOU,
									double					&fitted_power_o_OU,
									double					&fitted_lambda_OU,
									double					&significanceWN,
									double					&fitted_power_WN){
	long p	 	= maximumIndex(spectrum,1,spectrum.size()-1);
	if(p<minMode) return false;
	peakPower 	= spectrum[p];
	peakf 		= p*df;

	//guess reasonable parameter values (use as starting values for fitting)
	double powerInt 			= sum(spectrum,1,spectrum.size()-1)*df;
	double suspected_power_o  	= average(spectrum,1,3);
	double suspected_lambda 	= 4*powerInt/suspected_power_o;
	
	//Fit Ornstein-Uhlenbeck
	if(!fit_PS_OU(	df,
					spectrum,
					suspected_power_o,
					suspected_lambda,
					fitted_power_o_OU,
					fitted_lambda_OU)){
		return false;
	}
	
	//Use fitted OU as null-hypothesis for P-value
	significanceOU = significancePeakOUSS(	fitted_power_o_OU, 
											fitted_lambda_OU,
											0, // power_e = 0 for classical OU
											df, 
											spectrum.size()-1,
											minMode,
											peakf, 
											peakPower);
	
	
	//linear least squares is equivalent to taking arithmetic mean
	fitted_power_WN 	= average(spectrum,1,spectrum.size()-1);
	significanceWN 		= 1 - pow(1 - exp(-peakPower/fitted_power_WN), (spectrum.size()-1.0));
	
	return true;
}







bool significanceOfSpectralPeak_OUSS( 	double 					df,
										const vector<double> 	&times, 	// used to fit the OUSS model. assumed to be in increasing order.
										const vector<double> 	&values, 	// used to fit the OUSS model
										const vector<double> 	&spectrum,	// assumes that normalization used is the conventional for continuous-time stochastic-processes
										long					minMode, 	//dismiss (i.e. return false) if periodogram peak is at mode below minMode. If you already filtered your samples by such a criterion, you still need to pass this value so that P-values are correctly estimated.
										double					&peakf,
										double					&peakPower,
										double					true_lambda,	// debugging
										double					true_sigma,		// debugging
										double					true_epsilon,	// debugging
										double					&significanceOU,
										double					&fitted_lambda,
										double					&fitted_power_o,
										double					&fitted_power_e,
										double					&significanceWN,
										double					&fitted_power_WN,
										bool					&averagedPeriodogram){
	const long N	= times.size();
	const double T 	= times.back() - times.front();
	const double dt = T/(N-1.0);
	long n, p 	= maximumIndex(spectrum,1,spectrum.size()-1);
	if(p<minMode) return false;
	peakPower 	= spectrum[p];
	peakf 		= p*df;	
	averagedPeriodogram = false;

//	cout << "  PS LS fit suggests s_o=" << 	fitted_power_o_OU << ", lambda=" << fitted_lambda_OU << ", s_e=" << fitted_power_e_OU << endl;
	
	// might be needed for debugging
	const double true_power_o = SQ(true_sigma) * dt *(1+exp(-dt*true_lambda))/(1-exp(-dt*true_lambda));
	const double true_power_e = SQ(true_epsilon) * dt;
	
	
	if(!fit_DFTS_OUSS_to_periodogram(	df,
										dt,
										N,
										spectrum,
										true_lambda, 	// for debugging
										true_power_o,	// for debugging
										true_power_e,	// for debugging
										false,
										fitted_lambda,
										fitted_power_o,
										fitted_power_e)) return false;
	
	if(PS_MIN_SERIES_SEGMENT_SIZE_FOR_AVERAGED_PERIODOGRAM>1){
		// starting from non-averaged periodogram fit, try to fit averaged periodogram
		const long segmentsAim = min(N/max(1,PS_MIN_SERIES_SEGMENT_SIZE_FOR_AVERAGED_PERIODOGRAM), (long)((- log(PS_MIN_RHO_FOR_AVERAGING)/fitted_lambda)*(N-PS_MAX_OMITTED_DATA_POINTS_IN_AVERAGING)/T));
		if(segmentsAim>=PS_MIN_NUMBER_OF_SEGMENTS_FOR_AVERAGING){
			double df_forFitting, dt_forFitting;
			long N_forFitting, segmentsAchieved;
			vector<double> spectrum_forFitting;
			if(averagedLombScargleSpectrum(	values,
											times,
											0,
											values.size()-1,
											segmentsAim,
											PS_MIN_NUMBER_OF_SEGMENTS_FOR_AVERAGING,
											2,
											-1,
											segmentsAchieved,
											spectrum_forFitting,
											df_forFitting,
											dt_forFitting,
											N_forFitting)){
				if(fit_DFTS_OUSS_to_periodogram(	df_forFitting,
													dt_forFitting,
													N_forFitting,
													spectrum_forFitting,
													true_lambda, 	// for debugging
													true_power_o,	// for debugging
													true_power_e,	// for debugging
													true,
													fitted_lambda,
													fitted_power_o,
													fitted_power_e)){
													
					// if fitting was done using a different (dt,N), adjust to original (dt,N)
					double fitted_sigma = sqrt((fitted_power_o/dt_forFitting)*(1-exp(-fitted_lambda*dt_forFitting))/(1+exp(-fitted_lambda*dt_forFitting)));;
					double fitted_epsilon = sqrt(fitted_power_e/dt_forFitting);
					fitted_power_o = SQ(fitted_sigma) * dt *(1+exp(-dt*fitted_lambda))/(1-exp(-dt*fitted_lambda));	
					fitted_power_e = SQ(fitted_epsilon) * dt;
					averagedPeriodogram = true;
				}
			}
		}
	}
	
	if(!OUSS_FIT_PS_INSTEAD_OF_TS){
		//Fit Ornstein-Uhlenbeck state-space model
		
		
		// fit OUSS model to time series (restricted maximum likelihood)
		// use ML fit as initial guess
		double suspected_lambda		= fitted_lambda;
		double suspected_sigma		= sqrt((fitted_power_o/dt)*(1-exp(-fitted_lambda*dt))/(1+exp(-fitted_lambda*dt)));
		double suspected_amplitude	= suspected_sigma * sqrt(2*suspected_lambda);
		double suspected_epsilon	= sqrt(fitted_power_e/dt);
		double fitted_mu, fitted_sigma, fitted_amplitude, fitted_epsilon;
//		cout << "    and suggests lambda=" << fitted_lambda_OU << ",  sigma=" << suspected_sigma << ",  epsilon=" << suspected_epsilon << endl; //debug
		if(!fit_TS_ML_OUSS(	times,
							values,
							suspected_lambda,
							suspected_amplitude,
							suspected_epsilon,
							fitted_mu,
							fitted_lambda,
							fitted_amplitude,
							fitted_epsilon)){
			return false;
		}
		fitted_sigma	= fitted_amplitude/sqrt(2*fitted_lambda);
		fitted_power_o 	= SQ(fitted_sigma) * dt *(1+exp(-dt*fitted_lambda))/(1-exp(-dt*fitted_lambda));
		fitted_power_e 	= SQ(fitted_epsilon) * dt;		
	}
	
	
	
	//Use fitted OU as null-hypothesis for P-value
	if(P_VALUE_USING_SIMULATIONS){
//		cout 	<< " debug: true lambda=" << true_lambda << ", true power_o = " << true_power_o << ", true_power_e=" << true_power_e << "\n"
//				<< "   fitted lambda=" << fitted_lambda << ", fitted power_o = " << fitted_power_o << ", fitted_power_e=" << fitted_power_e << "\n"
//				<< "   estimating P-value" << endl;
				
		significanceOU = significancePeakOUSS_simul(fitted_lambda,		// debug true_lambda,
													fitted_power_o, 	// debug true_power_o,
													fitted_power_e, 	// debug true_power_e,
													dt,
													N,
													times,
													spectrum.size()-1,
													minMode,
													peakf,
													peakPower);
	}else if(CORRECT_FOR_DISCRETE_FINITE_TIME_SERIES){
		significanceOU = significancePeakOUSS(	fitted_lambda,	// debug true_lambda,
												fitted_power_o, 	// debug true_power_o,
												fitted_power_e, 	// debug true_power_e,
												df,
												spectrum.size()-1,
												minMode,
												dt,
												N,
												peakf, 
												peakPower);
	}else{
		significanceOU = significancePeakOUSS(	fitted_lambda,
												fitted_power_o, 
												fitted_power_e,
												df, 
												spectrum.size()-1,
												minMode,
												peakf, 
												peakPower);
	}	
	
	//linear least squares is equivalent to taking arithmetic mean
	fitted_power_WN = average(spectrum,1,spectrum.size()-1);
	significanceWN 	= 1 - pow(1 - exp(-peakPower/fitted_power_WN), (spectrum.size()-1.0));
	return true;
}



// not correct for periodograms of discrete time series
double OUSS_power_spectrum(double lambda, double power_o, double power_e, double f){
	return power_o * SQ(lambda)/(SQ(2*PI*f) + SQ(lambda)) + power_e;
}




#endif