#ifndef AUXILIARY_FUNCTIONS
#define AUXILIARY_FUNCTIONS

#include <vector>
#include <cmath>
#include <iostream>
#include <sstream>
#include <string>
#include <string.h>
#include <sys/stat.h>
#include <stdarg.h>
#include <algorithm>
#include <stdlib.h>
#include "InternalDefs.h"

using namespace std;


#define NORMAL_CYCLES_COUNT 30 //Random-Variables used to create Standard-Normal-Distribution



#pragma mark -
#pragma mark Math
#pragma mark



template<typename TYPE>
inline TYPE average(const vector<TYPE> &values, long start, long end){
	TYPE s(0);
	start = max(0l, start); end=min(end,(long)values.size()-1);
	for(long i=start; i<=end; ++i){ s += values[i]; }
	return (end>=start ? s/(end-start+1.0) : s);
}


template<class TYPE>
inline TYPE SQ(TYPE value){
	return value * value;
}

template<class TYPE>
inline TYPE Qube(TYPE value){
	return value * value * value;
}

template<class TYPE>
inline TYPE QTR(TYPE value){
	return SQ(SQ(value));
}




//multiply matrices A*B
//A has size rA*cA, B has size cA*cB
//both A & B are given in row-major format
//result AB is also given in row-major format and has size rA*cB (must be pre-allocated)
void matmul(long rA, long cA, long cB, const double A[], const double B[], double AB[]){
	double sum;
	long r,c,k;
	for(r=0; r<rA; ++r){
		for(c=0; c<cB; ++c){
			for(sum=0, k=0; k<cA; ++k){
				sum += A[r*cA+k] * B[k*cB+c];
			}
			AB[r*cB+c] = sum;
		}
	}
}


double varianceOfSeries(const vector<double> &signal){
	double sum=0, sum2=0;
	for(long n=0; n<signal.size(); ++n){
		sum += signal[n];
		sum2 += SQ(signal[n]);
	}
	sum /= signal.size();
	sum2 /= signal.size();
	return sum2 - SQ(sum);
}



double linearExtrapolate(double x1, double x2, double y1, double y2, double x){
	return y1 + (y2-y1)*(x-x1)/(x2-x1);
}




#pragma mark -
#pragma mark Random numbers
#pragma mark



inline double uniformWithinInclusive(double minimum, double maximum){
	return minimum + ((maximum-minimum)*rand())/RAND_MAX;
}

inline double uniformWithinExclusive(double minimum, double maximum){
	return minimum + (maximum-minimum)*EPSILON + ((maximum - minimum)*(1-2*EPSILON)*rand())/RAND_MAX;
}

inline double uniformWithinInclusiveRight(double minimum, double maximum){
	return minimum + EPSILON + ((maximum - minimum - EPSILON)*rand())/RAND_MAX;
}

double standardNormal(){
	return sqrt(-2.0*log(uniformWithinInclusiveRight(0, 1)))*cos(2.0*PI*uniformWithinInclusiveRight(0,1));
}

inline double exponential(double mean){
	return - mean * log(uniformWithinExclusive(0.0, 1.0));
}

inline double correlateOneStandardNormalToAnother(double rho, double x, double y){
	double alpha = rho/sqrt(1-SQ(rho));
	return (y + alpha*x)/sqrt(1+SQ(alpha));
}

void correlatedStandardNormalPair(double rho, double &x, double &y){
	x = standardNormal();
	y = standardNormal();
	y = correlateOneStandardNormalToAnother(rho, x, y);
}


inline double chiSQ(long k){
	double S=0;
	for(long i=0; i<k; ++i){ S += SQ(standardNormal()); }
	return S;
}


#pragma mark -
#pragma mark Vectorized operations
#pragma mark



template<class CONTAINER>
inline typename CONTAINER::value_type maximum(const CONTAINER &values){
	if(values.empty()) return 0;
	typename CONTAINER::value_type maxValue = *values.begin();
	for(typename CONTAINER::const_iterator i=values.begin(); i!=values.end(); ++i){
		maxValue = max(maxValue, *i);
	}
	return maxValue;
}

template<class TYPE>
inline TYPE maximum(const vector<TYPE> &values, long start, long end){
	if((values.size()<=end) || (end<start) || (start<0) || (end<0)) return 0;
	TYPE maxValue = values[start];
	for(long i=start; i<=end; ++i){
		maxValue = max(maxValue, values[i]);
	}
	return maxValue;
}


template<typename TYPE>
inline long maximumIndex(const vector<TYPE> &values, long start, long end){
	if((values.size()<=end) || (end<0) || (start<0) || (end<start)) return -1;
	long i, maxIndex = start;
	for(i=start+1; i<=end; ++i){
		if(values[maxIndex] < values[i]) maxIndex = i;
	}
	return maxIndex;
}


template<typename TYPE>
inline long minimumIndex(const vector<TYPE> &values){
	if(values.empty()) return -1;
	unsigned long i, minIndex = 0;
	for(i=1; i<values.size(); ++i){
		if(values[minIndex] > values[i]) minIndex = i;
	}
	return minIndex;
}



template<class TYPE>
inline TYPE minimum(const vector<TYPE> &values, long start, long end){
	if((values.size()<=end) || (end<start) || (start<0) || (end<0)) return 0;
	TYPE minValue = values[start];
	for(long i=start; i<=end; ++i){
		minValue = min(minValue, values[i]);
	}
	return minValue;
}

template<typename CONTAINER>
inline typename CONTAINER::value_type sum(CONTAINER values){
	typename CONTAINER::const_iterator i;
	typename CONTAINER::value_type s(0);
	for(i=values.begin(); i!=values.end(); ++i){
		s += *i;
	}
	return s;
}




template<typename TYPE>
inline TYPE sum(const vector<TYPE> &values, long start, long end){
	TYPE s(0);
	for(long i=max(0l,start); i<=end; ++i) s += values[i];
	return s;
}




#pragma mark -
#pragma mark String creation
#pragma mark -


template<class TYPE> 
string makeString(TYPE data){
	ostringstream stream;
	stream << data;
	return stream.str();
}



string vstringprintf(const char *format, va_list args){
	va_list temp;
	va_copy(temp, args);
	char *buffer = new char[vsnprintf(NULL, 0, format, temp) + 1];
	va_end(temp);
	
	vsprintf(buffer, format, args);
	string s(buffer);
	delete [] buffer;
	return s;
}



string stringprintf(const char *format, ...){
	string s;
	va_list args;
	va_start(args, format);
	s = vstringprintf(format, args);
	va_end(args);
	return s;
}



template<class TYPE>
void string2vector(const string &s, vector<TYPE> &v){
	istringstream ss(s);
	v.clear();
	TYPE d;
	do{
		ss >> d;
		if(ss.fail()) break;
		v.push_back(d);
	}while(true);
}


template<class TYPE>
vector<TYPE> string2vector(const string &s){
	vector<TYPE> v;
	string2vector<TYPE>(s,v);
	return v;
}


string vector2string(const vector<double> &v, string separator){
	ostringstream ss;
	for(long n=0; n<v.size(); ++n){
		ss << (n==0 ? "" : separator) << v[n];
	}
	return ss.str();
}



#pragma mark -
#pragma mark File IO
#pragma mark


string parentOf(const string &path){
	if(path == "") return "..";
	//string s = path.substr(0, path.length()-1);
	int p = path.find_last_of(DIRECTORY_SEPARATOR);
	if(p == path.npos){
		return ".";
	}else if(p == 0){
		return "/";
	}else if(p == path.length()-1){
		return parentOf(path.substr(0, path.length()-1));
	}else{
		return path.substr(0, p);
	}
}




bool fileExists(const string &path) { 
	struct stat fileInfo; 
	if(stat(path.c_str(),&fileInfo) == 0){ 
		return true; 
	}else{ 
		return false; 
	}
}



//make a copy of the passed file (if available)
void copySilently(const string &from, const string &to){
	if(fileExists(from)){
		ifstream src(from.c_str(), ios::binary);
		ofstream dst(to.c_str(), ios::binary);
		dst << src.rdbuf();
	}
}



bool createDirectorySupposingParentExists(const string &path){
	if(fileExists(path)) return true;
	#if defined(IS_UNIX) || defined(IS_APPLE)
		return (mkdir(path.c_str(), 0777) != -1);
	#elif defined(IS_WINDOWS)
		return (mkdir(path.c_str() != -1);
	#else
		return false;
	#endif
}
	
	
bool createDirectory(const string &path){
	int p;
	if(path.substr(0,strlen(DIRECTORY_SEPARATOR)) == DIRECTORY_SEPARATOR){
		p = strlen(DIRECTORY_SEPARATOR);
	}else{
		p=-1;
	}
	while((p = path.find_first_of(DIRECTORY_SEPARATOR,p+1)) != path.npos){
		if(! createDirectorySupposingParentExists(path.substr(0,p))){
			return false;
		}
	}
	return createDirectorySupposingParentExists(path);
}


bool openOutputFile(const string &path, ofstream &fout, ostream &feedbackStream){
	if(!createDirectory(parentOf(path))){
		feedbackStream << "Could not create directory " << parentOf(path) << "\n";
		return false;
	}
	fout.open(path.c_str());
	if(!fout.is_open()){
		feedbackStream << "Could not write to file " << path << "\n";
		return false;
	}
	return true;
}




string childOf(const string &parent, const string &childName){
	if((parent == "") || (parent == ".")){
		return childName;
	}else if(parent.substr(parent.length() - strlen(DIRECTORY_SEPARATOR)) == DIRECTORY_SEPARATOR){
		return parent + childName;
	}else{
		return parent + DIRECTORY_SEPARATOR + childName;
	}
}





string getNextNonExistingChild(	const string &parent, 
								const string &prefix, 
								const string &suffix, 
								unsigned int digits){
	unsigned int i=1;
	string path;
	do{
		path = childOf(parent, prefix + stringprintf("%0*d", digits, i) + suffix);
		++i;
	}while(fileExists(path));
	return path;
}



#pragma mark -
#pragma mark Statistics
#pragma mark 


double getCorrelation(const vector<double> &x, long lag){
	const long N = x.size()-abs(lag);
	lag = abs(lag);
	double sum1=0, sum2=0, sumSQ1=0, sumSQ2=0, sumP=0;
	for(long n=lag; n<x.size(); ++n){
		sum1 	+= x[n-lag];
		sum2 	+= x[n];
		sumSQ1 	+= SQ(x[n-lag]);
		sumSQ2 	+= SQ(x[n]);
		sumP	+= x[n-lag]*x[n];
	}
	
	double EP = sumP/N;
	double PE = (sum1/N)*(sum2/N);
	double COV = EP - PE;
	double V1 = sumSQ1/N - SQ(sum1/N);
	double V2 = sumSQ2/N - SQ(sum2/N);
	return COV/sqrt(V1*V2);
	
}



bool getEmpiricalQuantile(	vector<double> 			&samples, 	// (input) samples (x-values) from which to estimate quantile function. Samples might get re-ordered.
							const vector<double> 	&Pvalues, 	// (input) P-values for which to estimate x-thresholds. Must be in increasing order. Must be within [0,1].
							vector<double>			&Xvalues){	// (output) Estimated X-values at which the given P-values are realized
	const long NS = samples.size();
	const long NP = Pvalues.size();
	const double NSd = double(NS);
	long s,p;
	Xvalues.resize(NP);
	sort(samples.begin(), samples.end());
	for(s=p=0; s<NS; ++s){
		if(p>=NP) break;
		while((s+1)/NSd>Pvalues[p]){
			Xvalues[p] = (s==0 ? samples[0] : (samples[s]-samples[s-1]<EPSILON ? samples[s] : linearExtrapolate(s/NSd, (s+1)/NSd, samples[s-1], samples[s], Pvalues[p])));
			++p;
			if(p>=NP) break;
		}
	}
	for(; p<NP; ++p){
		Xvalues[p] = samples.back();
	}
	return true;
}


// Estimates cumulative distribution function (CDF) using a collection of samples drawn from a real probability distribution.
// All Xvalues on the left of all samples will be assigned a CDF=0, all Xvalues to the right of all samples will be assigned a CDF=1.
bool getEmpiricalCumulativeDF(	vector<double> 			&samples, 	 // (input) samples (x-values) from which to estimate cumulative distribution function. Samples might get re-ordered.
								const vector<double> 	&Xvalues, 	 // (input) X-values for which to estimate CDF. Must be in increasing order.
								vector<double>			&CDFvalues){ // (output) CDF values estimated for the givesn Xvalues.
	const long NS = samples.size();
	const long NX = Xvalues.size();
	const double NSd = double(NS);
	long s,x;
	CDFvalues.resize(NX);
	sort(samples.begin(), samples.end());
	for(s=x=0; s<NS; ++s){
		if(x>=NX) break;
		while(samples[s]>Xvalues[x]){
			CDFvalues[x] = (s==0 ? 0 : linearExtrapolate(samples[s-1], samples[s], s/NSd, (s+1)/NSd, Xvalues[x]));
			++x;
			if(x>=NX) break;
		}
	}
	for(; x<NX; ++x){
		CDFvalues[x] = 1;
	}
	return true;
}



#pragma mark -
#pragma mark Data structures
#pragma mark 




//Class for scanning a parameter on a regular 1D grid and performing statistics on each grid point
class Scanner1D{
private:
	vector<double> 	bin_values;
	vector<double>	bin_centres;
	vector<long> 	bin_sizes;
	double			minx;
	double			dx;
	long			ideal_bin_size;
public:
	double setup(long N_total, long N_bins, double _minx, double _maxx){
			minx = _minx; 
			ideal_bin_size = max(1l, N_total/N_bins);
			dx = (0 ? N_bins<2 : (_maxx - _minx)/(N_bins-1));
			bin_values.resize(N_bins,0.0);
			bin_sizes.resize(N_bins,0);
			bin_centres.resize(N_bins);
			for(long i=0; i<N_bins; ++i){ bin_centres[i] = _minx + i*dx; }
		}
		
	long get_bin_for_iteration(long i) const{ return min(i/ideal_bin_size, long(bin_centres.size()-1)); }
	double get_x_for_bin(long i) const{ return bin_centres[i]; }
	double add_y_to_bin(long i, double y){ ++bin_sizes[i]; bin_values[i] += y; }
	double get_overall_average() const{
			double sum=0;
			long count=0;
			for(long i=0; i<bin_centres.size(); ++i){ sum += bin_sizes[i] * bin_values[i]; count += bin_sizes[i]; }
			return (0 ? count==0 : sum/count);
		}
	double get_average_y_for_bin(long i) const{ return bin_values[i]/bin_sizes[i]; }
	long get_N_bins() const{ return bin_values.size(); }
};







//Class for scanning a 2 parameters on a regular 2D rectangular grid and performing statistics on each grid point
//Similar to Scanner1D
class Scanner2D{
private:
	vector<vector<double> > 	bin_values;
	vector<double>				bin_centres_x;
	vector<double>				bin_centres_y;
	vector<long> 				bin_sizes;
	double						minx, miny;
	double						dx, dy;
	long						ideal_bin_size;
	long						value_dim, N_bins_x, N_bins_y;
public:
	double setup(	long _value_dim, 
					long N_total, 
					long _N_bins_x, 
					long _N_bins_y, 
					double _minx, 
					double _maxx, 
					double _miny, 
					double _maxy){
			value_dim = _value_dim;
			N_bins_x = _N_bins_x;
			N_bins_y = _N_bins_y;
			minx = _minx;
			miny = _miny;
			long N_bins = N_bins_x*N_bins_y;
			ideal_bin_size = max(1l, N_total/N_bins);
			dx = (0 ? N_bins_x<2 : (_maxx - _minx)/(N_bins_x-1));
			dy = (0 ? N_bins_y<2 : (_maxy - _miny)/(N_bins_y-1));
			bin_values.resize(value_dim,vector<double>(N_bins,0.0));
			bin_sizes.resize(N_bins,0);
			bin_centres_x.resize(N_bins);
			bin_centres_y.resize(N_bins);
			
			long r,c;
			for(r=0; r<N_bins_x; ++r){
				for(c=0; c<N_bins_y; ++c){
					bin_centres_x[r*N_bins_y + c] = _minx + r*dx;
					bin_centres_y[r*N_bins_y + c] = _miny + c*dy;
				}
			}
		}
		
		
	long get_bin_for_iteration(long i) const{ return min(i/ideal_bin_size, long(bin_sizes.size()-1)); }
	double get_x_for_bin(long i) const{ return bin_centres_x[i]; }
	double get_x_for_bin_x(long i) const{ return bin_centres_x[i*N_bins_y + 0]; }
	double get_y_for_bin(long i) const{ return bin_centres_y[i]; }
	double get_y_for_bin_y(long i) const{ return bin_centres_y[0*N_bins_y + i]; }
	
	//expecting value_dim doubles as additional arguments
	double add_values_to_bin(long i, ...){ 
		++bin_sizes[i];
		va_list args;
		va_start(args, i);
		for(long j=0; j<value_dim; ++j){
			bin_values[j][i] += va_arg(args, double);
		}
		va_end(args);
	}
	
	double get_overall_average(long value_i, long &count) const{
			double sum=0;
			count=0;
			for(long i=0; i<bin_sizes.size(); ++i){ sum += bin_values[value_i][i]; count += bin_sizes[i]; }
			return (0 ? count==0 : sum/count);
		}
		
	long get_bin_size(long bin_i) const{ return bin_sizes[bin_i]; }
	long get_bin_xy_size(long bin_x_i, long bin_y_i) const{ return bin_sizes[bin_x_i*N_bins_y + bin_y_i]; }
	double get_average_value_for_bin(long bin_i, long value_i) const{ return bin_values[value_i][bin_i]/bin_sizes[bin_i]; }
	double get_average_value_for_bin_xy(long bin_x_i, long bin_y_i, long value_i) const{ long bin_i = bin_x_i*N_bins_y+bin_y_i; return bin_values[value_i][bin_i]/bin_sizes[bin_i]; }
	double get_average_value_for_bin_x(long bin_x_i, long value_i, long &count) const{
			double sum=0;
			long bin_i;
			count=0;
			for(long i=0; i<N_bins_y; ++i){ 
				bin_i = bin_x_i*N_bins_y + i;
				sum += bin_values[value_i][bin_i]; 
				count += bin_sizes[bin_i]; 
			}
			return (0 ? count==0 : sum/count);
		}
	double get_average_value_for_bin_y(long bin_y_i, long value_i, long &count) const{
			double sum=0;
			long bin_i;
			count=0;
			for(long i=0; i<N_bins_x; ++i){ 
				bin_i = i*N_bins_y + bin_y_i;
				sum += bin_values[value_i][bin_i]; 
				count += bin_sizes[bin_i]; 
			}
			return (0 ? count==0 : sum/count);
		}
	long get_N_bins() const{ return bin_sizes.size(); }
	long get_N_bins_x() const{ return N_bins_x; }
	long get_N_bins_y() const{ return N_bins_y; }
	
	
	void print_state() const{
			cout << "Scanner2D: " << N_bins_x << " x " << N_bins_y << " grid\n";
			for(long i=0; i<bin_sizes.size(); ++i){
				cout << " bin " << i << ": count " << bin_sizes[i] << ", x=" << bin_centres_x[i] << ", y=" << bin_centres_y[i] << "\n   values:";
				for(long j=0; j<value_dim; ++j){
					cout << "  " << bin_values[j][i];
				}
				cout << "\n";
			}
		}
};









#endif