/* 	This is the main source code for the simulations described in the article by Louca et al. - Detecting cyclicity in ecological time series
	The code is a Monte Carlo simulation of the models. 
	For each run the periodogram is calculated and evaluated as described in the article. 
	Summary statistics of this process are printed out at the end.
	Model parameters and simulation settings are set in the source file Config.h.

	Written by Stilianos Louca, December 21, 2013
*/




#include <iostream>
#include <fstream>
#include <vector>
#include <cmath>
#include <limits>
#include "Config.h"
#include "InternalDefs.h"
#include "Points.h"
#include "LombScargleSpectrum.h"
#include "StochasticRungeKutta.h"
#include "LogisticGrowthModel.h"
#include "GompertzGrowthModel.h"
#include "OrnsteinUhlenbeckModel.h"
#include "Auxiliary.cpp"
#include "PeakSignificance.cpp"
//#include "Grid3DInterpolator.h"
#include "GridMultilinearInterpolator.h"
#include "STPlot/sources/STPlot.h"
using namespace std;
using namespace ST;



#pragma mark -
#pragma mark Main
#pragma mark -






int main(int argc, char* argv[]){
	long						n,m,k,mode,series_size;
	bool						model_is_cyclic, shortcutted_simulation;
	ofstream 					fout, sub_log;
	string						fileBaseName, model_description, category;
	vector<double>				times,spectrum,signal;
	vector<XPoint>				signalTemp;
	double						error_SD, alpha, expected_power_o, expected_lambda, frequency, df, dummyD, dummyD1, dummyD2;
	double						significanceN, significanceOU, significanceWN,  fitted_power_o_OU, fitted_theta_OU, fitted_power_e_OU, fitted_power_WN;
	long						failedSpectrumRejects, lowModeRejects, averagedPeriodograms, failedSignificanceRejects, peakNotRealFreqRejects, failedSetup;
	long						alpha_epsilon_bin, countExamined, countCyclic[2][2], bin_size; //countCyclic[a][b] for a=1 iff non-OU, b=1 iff non-WN
	int							cyclicByOU, cyclicByWN;
	LGDynamics					LGmodel;
	OUDynamics					OUmodel;
	vector<double> 				eFAP_CDF_corrected, eFAP_CDF_uncorrected, eFAP_CDF_uncorrected_binnedByRho_and_NSR;
	GridMultilinearInterpolator	CDFInterpolator, FAPcorrector;
	string						CDFgrid_output_dir, all_plots_script;
	
	const vector<double> 		FAP_correction_grid_eFAPs = string2vector<double>(FAP_CORRECTION_GRID_eFAPs);
	const vector<double> 		FAP_correction_grid_rhos = string2vector<double>(FAP_CORRECTION_GRID_RHOs);
	const vector<double> 		FAP_correction_grid_NSRs = string2vector<double>(FAP_CORRECTION_GRID_NSRs);
	const long gridN_rhos = FAP_correction_grid_rhos.size();
	const long gridN_NSRs = FAP_correction_grid_NSRs.size();
	const long gridN_eFAPs = FAP_correction_grid_eFAPs.size();
	vector<double>				estimatedPvaluesOU, estimatedPvaluesOUuncorrected;
	vector<vector<vector<double> > >		estimatedPvaluesOUuncorrected_binnedByRho_and_NSR(gridN_rhos, vector<vector<double> >(gridN_NSRs));
	
	if((FAP_CDF_CORRECTION!="") && USE_STATE_SPACE_MODEL){
		cout << "  Found FAP corrector for OUSS model: Will be correcting eFAPs\n";
	}
	if(FAP_CDF_CORRECTION!=""){
		if(!FAPcorrector.load(4,FAP_CDF_CORRECTION)){
			cout << "  ERROR: Could not load FAP_CDF_CORRECTION string\n";
			return 0;
		}
	}
	
	
	
	
	
	//debug
	/*
	vector<double> dtimes = string2vector<double>("0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 6 6.5 7 7.5 8 8.5 9 9.5 10 10.5 11 11.5 12 12.5 13 13.5 14 14.5 15 15.5 16 16.5 17 17.5 18 18.5 19 19.5 20 20.5 21 21.5 22 22.5 23 23.5 24 24.5 25 25.5 26 26.5 27 27.5 28 28.5 29 29.5 30 30.5 31 31.5 32 32.5 33 33.5 34 34.5 35 35.5 36 36.5 37 37.5 38 38.5 39 39.5 40 40.5 41 41.5 42 42.5 43 43.5 44 44.5 45 45.5 46 46.5 47 47.5 48 48.5 49 49.5 50");
	vector<double> dvalues = string2vector<double>("0.941695 1.10381 0.927612 1.06133 0.724401 0.818952 0.833921 0.789225 1.06215 1.16206 0.862879 0.864727 1.35531 1.50136 1.09196 1.1722 1.20303 0.945381 0.828455 0.755115 1.06531 0.811853 0.445822 0.748415 0.749771 0.780996 1.0017 1.12162 1.2175 0.780049 0.910242 1.16985 1.35972 1.03736 1.13655 1.1875 0.810759 0.863924 0.96702 0.96132 0.814535 1.04363 0.743377 0.876182 0.968381 1.09869 1.12217 0.95409 0.716318 0.888672 1.07044 1.11015 0.927606 0.886407 0.93287 1.2959 1.47829 0.997813 1.25022 0.76487 1.06667 1.62357 1.04338 0.859184 1.12778 0.939522 0.717104 0.748877 0.704865 0.853298 1.26531 1.2543 1.23894 0.953348 1.06061 0.935198 1.37771 1.10347 1.05327 1.16348 1.03826 1.08092 0.950312 1.16824 1.27124 0.889717 0.659243 0.74293 1.15556 0.952268 1.18929 1.08938 0.800359 0.69744 0.918542 0.703612 1.05833 1.29011 1.20261 1.02755 1.05131");
	//calculate periodogram of sample path
	long N 		= dtimes.size();
	double T 	= (dtimes.back()-dtimes.front()), meandt = T/(N-1);
	if(!LombScargleSpectrum<double>(dvalues, 
									dtimes,
									0, 
									dtimes.size()-1,
									0,
									NULL,
									NULL,
									&spectrum,
									df)){
		cout << "error LSS\n"; return 0;
	}
	//get significance
	bool dummyB;
	if(!significanceOfSpectralPeak_OUSS(df, 
										dtimes, 
										dvalues, 
										spectrum, 
										min_mode, 
										dummyD1, 
										dummyD2, 
										1,
										1,
										1,
										significanceOU, 
										fitted_theta_OU, 
										fitted_power_o_OU, 
										fitted_power_e_OU, 
										significanceWN, 
										fitted_power_WN,
										dummyB)){
		cout << "error signif\n"; return 0;
	}
	const double fitted_rho = exp(-fitted_theta_OU*meandt);
	const double fitted_NSR = fitted_power_e_OU/fitted_power_o_OU;
	double fitted_sigma = sqrt((fitted_power_o_OU/meandt)*(1-exp(-fitted_theta_OU*meandt))/(1+exp(-fitted_theta_OU*meandt)));;
	double fitted_epsilon = sqrt(fitted_power_e_OU/meandt);

	double signifCorrected = FAPcorrector.interpolate(0, (double) dtimes.size(), (double) fitted_rho, (double) fitted_NSR, (double) significanceOU);
	cout 	<< " debug: fitted_rho=" << fitted_rho << ", fitted_NSR=" << fitted_NSR << ", fitted_lambda=" << fitted_theta_OU << ", fitted_power_o=" << fitted_power_o_OU << ", fitted_power_e=" << fitted_power_e_OU << "\n"
			<< "        df=" << df << ", dt=" << meandt << ", significanceWN=" << significanceWN << ", fitted_power_WN=" << fitted_power_WN << "\n"
			<< "        fitted sigma=" << fitted_sigma << ", fitted_epsilon=" << fitted_epsilon << "\n"
			<< "        significanceOU=" << significanceOU << ", corrected=" << signifCorrected << endl;
	return 0;
	*/
	
	
	
	
	//seed random number generator
	srand(time(0));

	//figure out time series sizes to consider
	vector<long> series_sizes;
	if(BUILD_FAP_CDF_GRID){
		if(simul_alpha_max>EPSILON*simul_K){
			cout << "  Warning: alpha_max>0 (i.e. model will be cyclic), but CDF grid calculation requested for FAP estimator. This makes little sense..\n";
		}
		string2vector<long>(FAP_CORRECTION_GRID_TSS, series_sizes);
		CDFInterpolator.setLayerPoints(4, FAP_CORRECTION_GRID_RHOs, FAP_CORRECTION_GRID_NSRs, FAP_CORRECTION_GRID_eFAPs);
		CDFgrid_output_dir = getNextNonExistingChild(OUTPUT_DIR, "FAP_CDF_calculation_", "", 2);
	}else{
		series_sizes.push_back(default_series_size);
	}
	
	for(long s=0; s<series_sizes.size(); ++s){
		// reset counts for this time series size
		series_size = series_sizes[s];
		estimatedPvaluesOUuncorrected_binnedByRho_and_NSR.assign(gridN_rhos,vector<vector<double> >(gridN_NSRs));
		estimatedPvaluesOUuncorrected.clear();
		estimatedPvaluesOU.clear();
		countExamined = 0;
		countCyclic[0][0]=countCyclic[0][1]=countCyclic[1][0]=countCyclic[1][1]=0;
		failedSpectrumRejects=lowModeRejects=failedSignificanceRejects=peakNotRealFreqRejects=failedSetup=averagedPeriodograms=0;
	
		//figure out output directory and files
		string output_subdir = (BUILD_FAP_CDF_GRID ? CDFgrid_output_dir+"/TS_"+makeString(series_size) : getNextNonExistingChild(OUTPUT_DIR, "run_", "", 2));
		createDirectory(output_subdir);
		all_plots_script = (BUILD_FAP_CDF_GRID ? CDFgrid_output_dir : output_subdir)+"/all_plots_script.sh";
		if(s==0){
			//make a copy of the config file (if available)
			copySilently("source/Config.h", (BUILD_FAP_CDF_GRID ? CDFgrid_output_dir : output_subdir)+"/Config.h");
		}
		
		if(BUILD_FAP_CDF_GRID){
			cout << "  Performing statistics for time series size " << series_size << "\n";
		}
		cout << "  Output will be written to: " << output_subdir << "\n";
		
		
		// write some feedback
		const string logFileName = output_subdir + string("/log.txt");
		openOutputFile(logFileName, sub_log, cout);
		sub_log	<< comment_prefix << " " << MODEL << " model: " << model_description << " (with variable parameters, see below)\n"
				<< comment_prefix << "   r = " << simul_r << (USE_OUSS_RHO_INSTEAD_OF_R && (MODEL=="OU") ? " (not relevant)" : "") << "\n"
				<< comment_prefix << "   OUSS rho = [" << simul_OUSS_rho_min << ", " << simul_OUSS_rho_max << "]" << (USE_OUSS_RHO_INSTEAD_OF_R && (MODEL=="OU") ? "" : " (not relevant)") << "\n"
				<< comment_prefix << "   K = " << simul_K << "\n"
				<< comment_prefix << "   SD, sigma = [" << simul_sigma_min << ", " << simul_sigma_max << "]" << (SCALE_LG_NOISE && (MODEL=="LG") ? " at equilibrium, but proportional to N/K in general" : "") << "\n"
				<< comment_prefix << "   cycle amplitude, alpha = [" << simul_alpha_min << ", " << simul_alpha_max << "]" << (simul_alpha_relative2SD ? " relative to sigma" : "") << "\n"
				<< comment_prefix << "   period, T = [" << simul_T_min << ", " << simul_T_max << "]\n"
				<< comment_prefix << "   epsilon = " << (USE_STATE_SPACE_MODEL ? stringprintf("[%.3g, %.3g]%s (state space model)", double(simul_epsilon_min), double(simul_epsilon_max), (simul_epsilon_relative2SD ? " relative to sigma" : "")) : "0") << "\n"
				<< comment_prefix << "   simulation time = " << simulation_time << "\n"
				<< comment_prefix << "   time series size = " << series_size << "\n"
				<< comment_prefix << "   frequency tolerance = " << frequency_tolerance << "\n"
				<< comment_prefix << "   min mode = " << min_mode << "\n"
				<< comment_prefix << "   OU significance threshold = " << significance_level_OU << "\n"
				<< comment_prefix << "   WN significance threshold = " << significance_level_WN << "\n"
				<< comment_prefix << "   FAP estimation using " << SIGNIF_BERNOULLI_TRIALS << " Bernoulli trials\n"
				<< comment_prefix << "   Monte Carlo count = " << Monte_Carlo_count << "\n";
		if(MODEL=="OU" && SHORTCUT_OUSS_SIMULATIONS){
			sub_log << comment_prefix << "   Time series generated using correlated draws\n";
		}else{
			sub_log << comment_prefix << "   Runge-Kutta integration time step = " << time_step << "\n";
		}
		sub_log << flush;


		long count, totalCount 	= Monte_Carlo_count;
		Scanner2D alpha_epsilon_scanner;
		alpha_epsilon_scanner.setup(2, totalCount, max(1,simul_alpha_grid_size), max(1,simul_epsilon_grid_size), simul_alpha_min, simul_alpha_max, simul_epsilon_min, simul_epsilon_max);
		long feedbackPeriod = max((long)1, (long) ((1.0*totalCount) / feedback_count));
		
		//Generate time series
		for(count=1, n=0; n<Monte_Carlo_count; ++n, ++count){
			//figure out some parameters
			double sigma = uniformWithinInclusive(simul_sigma_min, simul_sigma_max);
			double period = uniformWithinInclusive(simul_T_min, simul_T_max);
			
			// figure out cycle amplitude
			alpha_epsilon_bin = alpha_epsilon_scanner.get_bin_for_iteration(count-1);
			alpha = (simul_alpha_grid_size<=1 	? uniformWithinInclusive(simul_alpha_min, simul_alpha_max) // pick at random
												: alpha_epsilon_scanner.get_x_for_bin(alpha_epsilon_bin)); // pick sequentially on a regular grid
			if(simul_alpha_relative2SD){ alpha *= sigma; }
			
			//figure out how much 'measurement error' to add to the time series
			if(USE_STATE_SPACE_MODEL){
				error_SD = (simul_epsilon_grid_size<=1 	? uniformWithinInclusive((double)simul_epsilon_min, (double)simul_epsilon_max) 
														: alpha_epsilon_scanner.get_y_for_bin(alpha_epsilon_bin));
				if(simul_epsilon_relative2SD){ error_SD *= sigma; }
			}else{
				error_SD = 0;
			}
			double delta = simulation_time/(series_size-1.0); // first estimate of what the time step will be
			double simul_OUSS_rho, OUSS_lambda;
			if(USE_OUSS_RHO_INSTEAD_OF_R){
				simul_OUSS_rho = uniformWithinExclusive((double)simul_OUSS_rho_min, (double)simul_OUSS_rho_max);
				OUSS_lambda = -log(simul_OUSS_rho)/delta;
			}else{
				OUSS_lambda = simul_r;
				simul_OUSS_rho = exp(-OUSS_lambda*delta);
			}
			
			
			//setup & run population model
			bool dummyB;
			if(MODEL=="OU"){
				
				// Ornstein-Uhlenbeck model
				if(!OUmodel.setup(	OUSS_lambda,
									simul_K,
									sigma,
									period,
									alpha)){
					++failedSetup;
					continue;
				}
				if(SHORTCUT_OUSS_SIMULATIONS){
					//generate centralized OUSS time series using correlated draws
					generateOUSStimeSeries(	OUSS_lambda,
											sigma,
											0, // uncorrelated error terms will be added later
											delta,
											series_size,
											signal);
					times.resize(signal.size()); signalTemp.resize(signal.size());
					for(long i=0; i<times.size(); ++i){ times[i] = i*delta; }
					//superimpose the OUSS time series to the deterministic trajectory
					for(long i=0; i<times.size(); ++i){ signalTemp[i].X = signal[i] + OUmodel.getDeterministicTrajectory(times[i]); }
				
				}else{
					//use numerical SDE solver instead
					if(!StochasticRungeKutta(	OUmodel.getRandomValue(0), 
												simulation_time, 
												time_step, 
												OUmodel, 
												1, 
												false, 
												series_size,
												times, 
												signalTemp, 
												dummyB)){
						cout << "An error has occured during OU run " << (n+1) << "!\n";
						return 1;
					}
				}
				model_description 	= OUmodel.getShortFootnoteDescription() + ", rho="+makeString(exp(-OUSS_lambda*delta));
				expected_power_o 	= OUmodel.getPowerAtZero(time_step);
				expected_lambda 	= OUmodel.getRelaxationRate();
				model_is_cyclic 	= OUmodel.isCyclic();
			}else{
				// Logistic Growth model
				if(!LGmodel.setup(	simul_r,
									simul_K,
									sigma,
									period,
									alpha,
									SCALE_LG_NOISE)){
					++failedSetup;
					continue;
				}
				if(!StochasticRungeKutta_VN(LGmodel.getRandomValue(0), 
											simulation_time, 
											time_step, 
											LGmodel, 
											1, 
											false, 
											series_size,
											times, 
											signalTemp, 
											dummyB)){
					cout << "An error has occured during LG run " << (n+1) << "!\n";
					return 1;
				}
				model_description 	= LGmodel.getShortFootnoteDescription();
				expected_power_o 	= LGmodel.getPowerAtZero(time_step);
				expected_lambda 	= LGmodel.getRelaxationRate();
				model_is_cyclic 	= LGmodel.isCyclic();
			}
			model_description += ", error_SD =" + makeString(error_SD);
			delta = (times.back()-times.front())/(times.size()-1.0);
			
			//extract time series from simulation & add uncorrelated error terms
			signal.resize(signalTemp.size());
			for(m=0; m<signalTemp.size(); ++m){ 
				signal[m] = signalTemp[m].X + error_SD*standardNormal(); 
			}

			
			//calculate periodogram of sample path
			long N 		= times.size();
			double T 	= (times.back()-times.front()), meandt = T/(N-1);
			if(!LombScargleSpectrum<double>(signal, 
											times,
											0, 
											times.size()-1,
											0,
											NULL,
											NULL,
											&spectrum,
											df)){
				++failedSpectrumRejects;
				continue;
			}
		
			
			//analyze periodogram
			mode = maximumIndex(spectrum, 1, spectrum.size()-1);
			frequency = df*mode;
			bool reject = false;
			if((simul_alpha_max>EPSILON) && (frequency_tolerance*frequency<abs(frequency-1/period)) && (frequency_tolerance>0)){
				++peakNotRealFreqRejects;
				reject = true;
				category = "_false_frequency";
			}else if(mode < min_mode){
				++lowModeRejects;
				reject = true;
				category = "_low_mode";
			}else{
				if(USE_STATE_SPACE_MODEL){
					bool averagedPeriodogram;
					if(!significanceOfSpectralPeak_OUSS(df, 
														times, 
														signal, 
														spectrum, 
														min_mode, 
														dummyD1, 
														dummyD2, 
														OUSS_lambda,	// provide true params if needed for debugging
														sigma,			// provide true params if needed for debugging
														error_SD,		// provide true params if needed for debugging
														significanceOU, 
														fitted_theta_OU, 
														fitted_power_o_OU, 
														fitted_power_e_OU, 
														significanceWN, 
														fitted_power_WN,
														averagedPeriodogram)){
						++failedSignificanceRejects; 
						if(averagedPeriodogram) ++averagedPeriodograms;
						reject=true; 
						category = "_failed_fit";
					}
				}else{
					if(!significanceOfSpectralPeak_OU(df, T, spectrum, min_mode, dummyD1, dummyD2, significanceOU, fitted_power_o_OU, fitted_theta_OU, significanceWN, fitted_power_WN)){
						++failedSignificanceRejects;
						reject=true; 
						category = "_failed_fit";
					}
					fitted_power_e_OU = 0;
				}
			}
			const double fitted_rho = exp(-fitted_theta_OU*meandt);
			const double fitted_NSR = fitted_power_e_OU/fitted_power_o_OU;
			
			// count towards statistics
			if(!reject){
				const long rho_bin = GridMultilinearInterpolator::findOn1DGrid(FAP_correction_grid_rhos, fitted_rho);
				const long NSR_bin = GridMultilinearInterpolator::findOn1DGrid(FAP_correction_grid_NSRs, fitted_NSR);
				estimatedPvaluesOUuncorrected_binnedByRho_and_NSR[max(0l,min(gridN_rhos-1,rho_bin))][max(0l,min(gridN_NSRs-1,NSR_bin))].push_back(significanceOU);
				estimatedPvaluesOUuncorrected.push_back(significanceOU);
				if((FAP_CDF_CORRECTION!="") && USE_STATE_SPACE_MODEL){
					//correct OUSS FAP using pre-calculated correction grid
					significanceOU = FAPcorrector.interpolate(0, (double) series_size, (double) fitted_rho, (double) fitted_NSR, (double) significanceOU);
				}
				cyclicByOU = (significanceOU < significance_level_OU);
				cyclicByWN = (significanceWN < significance_level_WN);	
				category = string(cyclicByOU==1 ? "_OU" : "_") + (cyclicByWN==1 ? "_WN" : "_");			
				++countExamined;
				++countCyclic[(cyclicByOU ? 1 : 0)][(cyclicByWN ? 1 : 0)];
				estimatedPvaluesOU.push_back(significanceOU);
				alpha_epsilon_scanner.add_values_to_bin(alpha_epsilon_bin, (cyclicByOU ? 1.0 : 0.0), (cyclicByWN ? 1.0 : 0.0));
			}
			
			
			//also calculate averaged periodogram (if applicable)
			bool averaged_spectrum = false;
			long N_averaged = N;
			double dt_averaged = delta, df_averaged=df;
			vector<double> spectrum_averaged;
			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_theta_OU)*(N-PS_MAX_OMITTED_DATA_POINTS_IN_AVERAGING)/T));
				if(segmentsAim>=PS_MIN_NUMBER_OF_SEGMENTS_FOR_AVERAGING){
					long segmentsAchieved;
					if(averagedLombScargleSpectrum(	signal,
													times,
													0,
													times.size()-1,
													segmentsAim,
													PS_MIN_NUMBER_OF_SEGMENTS_FOR_AVERAGING,
													2,
													-1,
													segmentsAchieved,
													spectrum_averaged,
													df_averaged,
													dt_averaged,
													N_averaged)){
						averaged_spectrum = true;
						// correct for different effective delta
						for(long m=0; m<spectrum_averaged.size(); ++m){ spectrum_averaged[m] *= delta/dt_averaged; }
					}
				}
			}
			
			//save this time series if desired
			if(SAVE_TIME_SERIES){
				const string runFileName = output_subdir + string("/time_series") + category + string("/") + makeString(n+1);
				if(openOutputFile(runFileName, fout, cout)){
					//header
					fout 	<< comment_prefix << "Time series generated by " << (model_is_cyclic ? "cyclic" : "non-cyclic") << " " << MODEL << " model:\n" 
							<< comment_prefix << "  " << model_description << "\n"
							<< comment_prefix << "  Expected OU lambda = " << expected_lambda << ",  s_o=" << expected_power_o << (USE_STATE_SPACE_MODEL ? ",  s_e=" + makeString(LombScargleSpectrum_of_normal_errors(error_SD,T,N)) : "") << ", delta=" << delta << "\n"
							<< comment_prefix << "  Sample size = " << times.size() << ", spanning " << T << " time units\n"
							<< comment_prefix << "  Periodogram peak frequency = " << frequency << "\n";
					if(!reject){
						fout	<< comment_prefix << "  Periodogram peak significance " << significanceOU << " vs OU (fitted s_o=" << fitted_power_o_OU << ", lambda=" << fitted_theta_OU << (USE_STATE_SPACE_MODEL ? ", s_e=" + makeString(fitted_power_e_OU) : "") << "), " << significanceWN << " vs WN (fitted s=" << fitted_power_WN << ")\n"
								<< comment_prefix << "  --> Classified as not " << (cyclicByOU==1 ? "OU" : "..") << (cyclicByWN==1 ? ",WN" : ",..") << "\n";
					}
					
					//time series
					fout		<< comment_prefix << "Data index 0\n"
								<< comment_prefix << "Data columns: time  population_size\n";
					for(m=0; m<signal.size(); ++m){
						fout << times[m] << " " << signal[m] << "\n";
					}
					
					//periodogram & fitted spectra
					fout 	<< "\n\n"
							<< comment_prefix << "Data index 1\n"
							<< comment_prefix << "Periodogram\n"
							<< comment_prefix << "Data columns: frequency  periodogram" << (reject ? "" : "  fitted_OU_power_spectrum") << "\n";
					for(m=1; m<spectrum.size(); ++m){
						fout << m*df << " " << spectrum[m];
						if(!reject) fout << " " << (CORRECT_FOR_DISCRETE_FINITE_TIME_SERIES ? PS_DFTS_OUSS(fitted_theta_OU, fitted_power_o_OU, fitted_power_e_OU, delta, N, m*df) : OUSS_power_spectrum(fitted_theta_OU, fitted_power_o_OU, fitted_power_e_OU, m*df));
						fout <<  "\n";
					}

					// averaged periodogram
					if(averaged_spectrum){
						fout 	<< "\n\n"
								<< comment_prefix << "Data index 2\n"
								<< comment_prefix << "Averaged Periodogram\n"
								<< comment_prefix << "Data columns: frequency  averaged_periodogram\n";
						for(m=1; m<spectrum_averaged.size(); ++m){
							fout << m*df_averaged << " " << spectrum_averaged[m];
							fout <<  "\n";
						}
					}
					fout.close();
					
					//plot
					if(PLOT_TIME_SERIES){
						MultiPlot multiPlot;
						Plot plot;
						//time series
						plot.title = "Time series";
						plot.axisLabels = MultiLabel("time", "X");
						plot.sources.push_back(PlotSource(PlotType2DCurve, "OU", runFileName, 0, 1, 2));
						for(long i=0; i<plot.sources.size(); ++i){ plot.sources[i].setLineWidth(2); }
						multiPlot.plots.push_back(plot);
						//periodogram and fitted power spectra
						plot.showKeys = true;
						plot.title = "Periodogram & fitted power spectra";
						plot.axisLabels = MultiLabel("frequency", "power");
						plot.range.xRange.setRangeMin(0);
						plot.yLogarithmic=false;
						plot.sources.clear();
						plot.sources.push_back(PlotSource(PlotType2DCurve, "periodogram", runFileName, 1, 1, 2));
						if(averaged_spectrum) plot.sources.push_back(PlotSource(PlotType2DCurve, "averaged", runFileName, 2, 1, 2));
						if(!reject){
							plot.showKeys = true;
							plot.sources.push_back(PlotSource(PlotType2DCurve, "OU PS", runFileName, 1, 1, 3));
							plot.sources.push_back(PlotSource(PlotType2DAnalytic, "WN PS", makeString(fitted_power_WN), 0));
						}
						for(long i=0; i<plot.sources.size(); ++i){ plot.sources[i].setLineWidth(2); }
						multiPlot.plots.push_back(plot);
						multiPlot.footnote = stringprintf("%s\nP_{OU}=%.3g, P_{WN}=%.3g", model_description.c_str(), significanceOU, significanceWN);
						//plot
						string errorMessage;
						plotWithLog(multiPlot, 
									runFileName+"_plot_script.gp",
									runFileName+"_plot.pdf",
									runFileName+"_plot.svg",
									runFileName, 
									NULL, 
									false, 
									all_plots_script,
									errorMessage);
					}
					
				}
			}
			
			//feedback on progress
			if((count % feedbackPeriod) == 0){
				cout << "\n    " << int((100.0 * count)/totalCount) << "% finished.." << (count == totalCount? "\n" : "") << flush;
				sub_log << "\n#    " << int((100.0 * count)/totalCount) << "% finished.." << (count == totalCount? "\n" : "") << flush;
			}
			
		}	
		cout << endl;
		
		
		//estimate cumulative distribution function of raw (uncorrected) eFAP
		getEmpiricalCumulativeDF(estimatedPvaluesOUuncorrected, FAP_correction_grid_eFAPs, eFAP_CDF_uncorrected);
		//do the same for the rho-binned ones, for each rho-bin separately
		//store in row-major (=rho-major) format
		eFAP_CDF_uncorrected_binnedByRho_and_NSR.resize(gridN_rhos*gridN_NSRs*gridN_eFAPs);
		vector<double> tempV;
		for(long r=0; r<gridN_rhos; ++r){
			for(long n=0; n<gridN_NSRs; ++n){
				// get empirical CDF for rho x NSR bin
				getEmpiricalCumulativeDF(estimatedPvaluesOUuncorrected_binnedByRho_and_NSR[r][n], FAP_correction_grid_eFAPs, tempV);
				for(long i=0; i<gridN_eFAPs; ++i){ 
					eFAP_CDF_uncorrected_binnedByRho_and_NSR[r*gridN_NSRs*gridN_eFAPs + n*gridN_eFAPs + i] = tempV[i]; 
				}
			}
		}
		if(BUILD_FAP_CDF_GRID){
			CDFInterpolator.appendLayer(series_size, eFAP_CDF_uncorrected_binnedByRho_and_NSR);
		}
		//estimate CDF for corrected eFAP. This might be used for comparing the 'efficiency' of the correction algorithms
		getEmpiricalCumulativeDF(estimatedPvaluesOU, FAP_correction_grid_eFAPs, eFAP_CDF_corrected);
		
		

		//feedback on results
		long cOU = countCyclic[1][0]+countCyclic[1][1], cWN = countCyclic[0][1]+countCyclic[1][1], cOUWN = countCyclic[1][1];
		long countRejected=Monte_Carlo_count - countExamined;
		sub_log	<< comment_prefix << "\n"
				<< comment_prefix << " Dismissed (i.e. did not count) " << countRejected << " out of " << Monte_Carlo_count << " simulations: " << failedSetup << " failed setup, " << failedSpectrumRejects << " failed spectrum, " << lowModeRejects << " low-mode(<" << min_mode << "), " << failedSignificanceRejects << " failed fitting and P-value estimation, " << (frequency_tolerance>0 ? makeString(peakNotRealFreqRejects) + " peak not real freq" : "no frequency filtering") << ", " << averagedPeriodograms << " periodograms were averaged for fitting\n"
				<< comment_prefix << " Of " << countExamined << " examined, cyclic by:\n"
				<< comment_prefix << "   OU: " << cOU << " (" << (100.0*cOU)/countExamined << "%)\n"
				<< comment_prefix << "   WN: " << cWN << " (" << (100.0*cWN)/countExamined << "%)\n"
				<< comment_prefix << "   OU & WN: " << cOUWN << " (" << (100.0*cOUWN)/countExamined << "%)\n";
		//cycle detection rate over varying alpha
		if(simul_alpha_grid_size>1){
			sub_log 	<< "\n\n" << comment_prefix << " Detected cycles by OU & WN for varying alpha" << (simul_alpha_relative2SD ? " (relative to OU SD)" : "") << ", averaged over all epsilon\n"
					<< comment_prefix << " Data columns: alpha\t\tOU-cycles(fraction)\t\tWN-cycles(fraction)\t\tN\n";
			for(long i=0; i<alpha_epsilon_scanner.get_N_bins_x(); ++i){
				dummyD = alpha_epsilon_scanner.get_average_value_for_bin_x(i,0,bin_size);
				sub_log << alpha_epsilon_scanner.get_x_for_bin_x(i) << "\t" << dummyD << "\t" << alpha_epsilon_scanner.get_average_value_for_bin_x(i,1,bin_size) << " " << bin_size << "\n";
			}
		}
		//cycle detection rate over varying epsilon
		if(USE_STATE_SPACE_MODEL && (simul_epsilon_grid_size>1)){
			sub_log 	<< "\n\n" << comment_prefix << " Detected cycles by OU & WN for varying epsilon" << (simul_epsilon_relative2SD ? " (relative to OU SD)" : "") << ", averaged over all alpha\n"
					<< comment_prefix << " Data columns: epsilon\t\tOU-cycles(fraction)\t\tWN-cycles(fraction)\t\tN\n";
			for(long i=0; i<alpha_epsilon_scanner.get_N_bins_y(); ++i){
				dummyD = alpha_epsilon_scanner.get_average_value_for_bin_y(i,0,bin_size);
				sub_log << alpha_epsilon_scanner.get_y_for_bin_y(i) << "\t" << dummyD << "\t" << alpha_epsilon_scanner.get_average_value_for_bin_y(i,1,bin_size) << " " << bin_size << "\n";
			}
		}
		//cycle detection rate over varying alpha & epsilon (2D grid)
		if(simul_alpha_grid_size && (USE_STATE_SPACE_MODEL && (simul_epsilon_grid_size>1))){
			sub_log 	<< "\n\n" << comment_prefix << " Detected cycles by OU & WN for varying alpha" << (simul_alpha_relative2SD ? " (relative to OU SD)" : "") << " and/or epsilon" << (simul_epsilon_relative2SD ? " (relative to OU SD)" : "") << "\n"
					<< comment_prefix << " Data columns: alpha\t\tepsilon\t\tOU-cycles(fraction)\t\tWN-cycles(fraction)\t\tN\n";
			for(long r=0; r<alpha_epsilon_scanner.get_N_bins_x(); ++r){
				for(long c=0; c<alpha_epsilon_scanner.get_N_bins_y(); ++c){
					sub_log << alpha_epsilon_scanner.get_x_for_bin_x(r) << "\t"  << alpha_epsilon_scanner.get_y_for_bin_y(c) << "\t" << alpha_epsilon_scanner.get_average_value_for_bin_xy(r,c,0) << "\t" << alpha_epsilon_scanner.get_average_value_for_bin_xy(r,c,1) << " " << alpha_epsilon_scanner.get_bin_xy_size(r,c) << "\n";
				}
				sub_log << "\n";
			}
		}
		//empirical CDF
		sub_log 	<< "\n\n" << comment_prefix << " Empirical cumulative distribution function of raw (uncorrected) eFAPs\n"
				<< comment_prefix << " Data columns: eFAP   CDF\n";
		for(long i=0; i<gridN_eFAPs; ++i){
			sub_log << FAP_correction_grid_eFAPs[i] << "\t" << eFAP_CDF_uncorrected[i] << "\n";
		}
		//empirical CDF in one line
		sub_log 	<< "\n\n" << comment_prefix << " CDF values as above but in a single line\n";
		for(long i=0; i<gridN_eFAPs; ++i){
			sub_log << (i==0 ? "" : " ") << eFAP_CDF_uncorrected[i];
		}
		sub_log << "\n";
		//empirical CDF of corrected eFAPs
		sub_log 	<< "\n\n" << comment_prefix << " Empirical CDF of (potentially) corrected eFAPs\n"
				<< comment_prefix << " Data columns: corrected_eFAP   CDF\n";
		for(long i=0; i<gridN_eFAPs; ++i){
			sub_log << FAP_correction_grid_eFAPs[i] << "\t" << eFAP_CDF_corrected[i] << "\n";
		}
		//empirical CDF in one line
		sub_log 	<< "\n\n" << comment_prefix << " CDF values as above but in a single line\n";
		for(long i=0; i<gridN_eFAPs; ++i){
			sub_log << (i==0 ? "" : " ") << eFAP_CDF_corrected[i];
		}
		sub_log << "\n";
		if((FAP_CDF_CORRECTION!="") && USE_STATE_SPACE_MODEL){
			// compare the empirical CDF of the uncorrected eFAPs and the FAPcorrector's predictions
			double SSR = 0, SS=0, S=0;
			for(long r=0; r<gridN_rhos; ++r){
				for(long n=0; n<gridN_NSRs; ++n){
					for(long i=0; i<gridN_eFAPs; ++i){
						SSR += SQ(eFAP_CDF_uncorrected_binnedByRho_and_NSR[r*gridN_eFAPs*gridN_NSRs + n*gridN_eFAPs + i] - FAPcorrector.interpolate(0, (double) series_size, (double) FAP_correction_grid_rhos[r], (double) FAP_correction_grid_NSRs[n], (double) FAP_correction_grid_eFAPs[i]));
						SS 	+= SQ(eFAP_CDF_uncorrected_binnedByRho_and_NSR[r*gridN_eFAPs*gridN_NSRs + n*gridN_eFAPs + i]);
						S	+= eFAP_CDF_uncorrected_binnedByRho_and_NSR[r*gridN_eFAPs*gridN_NSRs + n*gridN_eFAPs + i];
					}
				}
			}
			const double SST = SS - SQ(S)/(gridN_rhos*gridN_NSRs*gridN_eFAPs);
			const double R2  = 1 - SSR/SST;
			sub_log 	<< "\n\n" << comment_prefix << " R2 (coefficient of determination) of FAP corrector, tested against empirical CDF of uncorrected eFAPs: R2=" << R2 << "  (SSR=" << SSR << ", SST=" << SST << ")\n";
		}
		//empirical CDF of uncorrected eFAPs, binned by rho
		sub_log << "\n\n" << comment_prefix << " Empirical CDF of raw (uncorrected) eFAPs, binned by fitted_rho and fitted_NSR\n"
				<< comment_prefix << " rho bin counts:\n";
		for(long r=0; r<gridN_rhos; ++r){ 
			for(long n=0; n<gridN_NSRs; ++n){ 
				sub_log << comment_prefix << "   rho>=" << FAP_correction_grid_rhos[r] << ",  NSR>=" << FAP_correction_grid_NSRs[n] << ": " << estimatedPvaluesOUuncorrected_binnedByRho_and_NSR[r][n].size() << "\n"; 
			}
		}
		sub_log	<< comment_prefix << " Data columns: fitted_rho  fitted_NSR   eFAP   CDF\n";
		for(long r=0; r<gridN_rhos; ++r){
			for(long n=0; n<gridN_NSRs; ++n){ 
				for(long i=0; i<gridN_eFAPs; ++i){
					sub_log << FAP_correction_grid_rhos[r] << "\t" << FAP_correction_grid_NSRs[n] << "\t" << FAP_correction_grid_eFAPs[i] << "\t" << eFAP_CDF_uncorrected_binnedByRho_and_NSR[r*gridN_NSRs*gridN_eFAPs + n*gridN_eFAPs + i] << "\n";
				}
				sub_log << "\n";
			}
		}
		//empirical CDF, binned by rho, in one line
		sub_log << "\n\n" << comment_prefix << " CDF of uncorrected eFAPs same as above, but in a single line (" << gridN_rhos << " rhos  x  " << gridN_NSRs << " NSRs  x  " << gridN_eFAPs << " eFAPs, sheet in row-major format)\n";
		for(long r=0; r<gridN_rhos; ++r){
			for(long n=0; n<gridN_NSRs; ++n){ 
				for(long i=0; i<gridN_eFAPs; ++i){
					sub_log << (i+r+n==0 ? "" : " ") << eFAP_CDF_uncorrected_binnedByRho_and_NSR[r*gridN_NSRs*gridN_eFAPs + n*gridN_eFAPs + i];
				}
			}
		}
		sub_log << "\n";
		sub_log.close();

		
		
		if(PLOT_SUMMARY){
			//generate summary plots
			MultiPlot multiPlot;
			Plot plot;
			plot.showKeys = true;
			long index=0;
			//cycle detection over alpha
			if(simul_alpha_grid_size>1){
				plot.title = "Cycle detection rate over $alpha$";
				plot.axisLabels = MultiLabel("$alpha$", "fraction");
				plot.range = MultiRange(Range(), Range(0,1));
				plot.sources.push_back(PlotSource(PlotType2DCurve, "OU", logFileName, index, 1, 2));
				plot.sources.push_back(PlotSource(PlotType2DCurve, "WN", logFileName, index, 1, 3));
				multiPlot.plots.push_back(plot);
				++index;
			}
			//cycle detection over epsilon
			if(USE_STATE_SPACE_MODEL && (simul_epsilon_grid_size>1)){
				plot.title = "Cycle detection rate over $epsilon$";
				plot.axisLabels = MultiLabel("$epsilon$", "fraction");
				plot.range = MultiRange(Range(), Range(0,1));
				plot.sources.clear();
				plot.sources.push_back(PlotSource(PlotType2DCurve, "OU", logFileName, index, 1, 2));
				plot.sources.push_back(PlotSource(PlotType2DCurve, "WN", logFileName, index, 1, 3));
				multiPlot.plots.push_back(plot);
				++index;
			}
			if(simul_alpha_grid_size && (USE_STATE_SPACE_MODEL && (simul_epsilon_grid_size>1))) ++index;
			//empirical CDF for uncorrected eFAPs
			plot.title = "Empirical CDF for estimated (uncorrected) eFAPs";
			plot.axisLabels = MultiLabel("eFAP", "CDF");
			plot.range = MultiRange(Range(), Range());
			plot.sources.clear();
			plot.xLogarithmic=plot.yLogarithmic=true;
			plot.sources.push_back(PlotSource(PlotType2DAnalytic, "ideal", "x", 0)); plot.sources.back().setColor(Color(0.5,0.5,0.5));
			plot.sources.push_back(PlotSource(PlotType2DCurve, "FAP CDF", logFileName, index, 1, 2)); plot.sources.back().setColor(Color(0.9,0,0));
			multiPlot.plots.push_back(plot);
			++index; ++index;
			
			//empirical CDF for corrected eFAPs
			plot.title = "Empirical CDF for estimated (corrected) eFAPs";
			plot.axisLabels = MultiLabel("eFAP", "CDF");
			plot.range = MultiRange(Range(), Range());
			plot.sources.clear();
			plot.xLogarithmic=plot.yLogarithmic=true;
			plot.sources.push_back(PlotSource(PlotType2DAnalytic, "ideal", "x", 0)); plot.sources.back().setColor(Color(0.5,0.5,0.5));
			plot.sources.push_back(PlotSource(PlotType2DCurve, "FAP CDF", logFileName, index, 1, 2)); plot.sources.back().setColor(Color(0,0,0.8));
			multiPlot.plots.push_back(plot);
			++index; ++index;
			//plot
			string errorMessage;
			plotWithLog(multiPlot, 
						output_subdir+"/summary_plot_script.gp",
						output_subdir+"/summary_plot.pdf",	
						output_subdir+"/summary_plot.svg",	
						"summary", 
						NULL, 
						false, 
						all_plots_script,
						errorMessage);
		}
	}
	
	//summarize CDF estimation and write to log
	if(BUILD_FAP_CDF_GRID){
		const string superlogFileName = CDFgrid_output_dir+"/log.txt";
		if(openOutputFile(superlogFileName, fout, cout)){
			fout 	<< comment_prefix << " CDF estimation for estimated, raw (uncorrected) FAPs (eFAP)\n"
					<< comment_prefix << " CDF specified on a TSS x fitted_rho x fitted_NSR x eFAP grid\n"
					<< comment_prefix << " Time series sizes (TSS): " << FAP_CORRECTION_GRID_TSS << "\n"
					<< comment_prefix << " rhos (correlation between time steps): " << FAP_CORRECTION_GRID_RHOs << "\n"
					<< comment_prefix << " NSRs (noise-to-signal ratio): " << FAP_CORRECTION_GRID_NSRs << "\n"
					<< comment_prefix << " eFAPs: " << FAP_CORRECTION_GRID_eFAPs << "\n";
			fout << "\n\n" << comment_prefix << " Index 0\n" << comment_prefix << " As single line string:\n" << CDFInterpolator.getGridAsString() << "\n";
			fout << "\n\n" << comment_prefix << " Index 1\n" << comment_prefix << " As scan\nData columns: TSS  fitted_rho  fitted_NSR  eFAP  CDF\n";
			CDFInterpolator.printGridScann(fout);
			fout.close();
		}
	}
	
	return 0;	
}