// allanvar.cpp
// C++ Code to calculate the overlapping (modified) Allan Variance of a time series
// Assumes tab-delimited format
//Michael Nickerson 2012/08/08

/*-
 * Copyright (c) 2013 Michael Nickerson
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
 * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
 * PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS
 * BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 * POSSIBILITY OF SUCH DAMAGE.
 */


//Tested functional on data from "The Calculation of Time Domain Frequency Stability" by W. J. Riley
//Also located in NIST Special Publication 1065, "Handbook of Frequency Stability Analysis"
//9-point data from NBS Monograph 140
//Data available at http://www.wriley.com/
//Caveat: This does not interpolate between unevenly spaced bins, so there may be some jitter


#include <deque>
#include <fstream>
#include <math.h>
#include <float.h>

using namespace std;

//Minimum number of samples required
#define MIN_SAMPS 4
#define _max(a,b)  (((a) > (b)) ? (a) : (b))
#define _min(a,b)  (((a) < (b)) ? (a) : (b))

struct timedata {
	double time;
	double data;
};

//Local function declarations
void showHelp();
void doAVAR(char *infile, int tcol, int dcol, double tmult, int numtaus, double mintau, double maxtau, bool logtau, bool overlap, double jitter);	//Main loop
deque<timedata> movingMean(deque<timedata> *input, double tau, bool overlap);	//Moving mean of a dataset, window size tau
double avar(deque<timedata> *data, double tau, double jitter);	//Get the Allan variance of a prebinned dataset
double getMean(deque<double> *input);	//Average a dataset
timedata getMean(deque<timedata> *input);	//Average a dataset


int main(int argc, char* argv[]) {
	//Abort if invalid inputs given or help requested
	if (argc < 2) {
		showHelp();
		return 0;
	} else if (!strcmp(argv[1], "-h") || !strcmp(argv[1], "/?") || !strcmp(argv[1], "-?")) {
		showHelp();
		return 0;
	}
	
	int tcol = 1, dcol = 2, numtaus = 100;
	double tmult = 1.0, mintau = 0, maxtau = 0, jitter = 0;
	bool logtau = false, setdcol = false, overlap = true;
	//Process each argument
	for (int i=1; i<argc; i++) {
		if (!strcmp(argv[i], "-time")) {
			tcol = atoi(argv[++i]);
			printf("Using column %i as time.\n", tcol);
		} else if (!strcmp(argv[i], "-data")) {
			dcol = atoi(argv[++i]); setdcol = true;
			printf("Using column %i as data.\n", dcol);
		} else if (!strcmp(argv[i], "-tmult")) {
			tmult = atof(argv[++i]);
			printf("Multiplying time by %f before use.\n", tmult);
		} else if (!strcmp(argv[i], "-logtau")) {
			logtau = true;
			printf("Using logarithmic tau spacing.\n");
		} else if (!strcmp(argv[i], "-numtaus") || !strcmp(argv[i], "-n")) {
			numtaus = atoi(argv[++i]);
			printf("Calculating for %i tau values\n", numtaus);
		} else if (!strcmp(argv[i], "-mintau")) {
			mintau = atof(argv[++i]);
			printf("Using %f for minimum tau size.\n", mintau);
		} else if (!strcmp(argv[i], "-maxtau")) {
			maxtau = atof(argv[++i]);
			printf("Using %f for maximum tau size.\n", maxtau);
		} else if (!strcmp(argv[i], "-nooverlap")) {
			overlap = false;
			printf("Disabling overlap - not using moving average.\n");
		} else if (!strcmp(argv[i], "-jitter")) {
			jitter = atof(argv[++i]);
			printf("Allowing %f fractional jitter in bin spacing.\n", jitter);
		} else {
			if (!setdcol) dcol = tcol + 1;	//If not explicitly specified, then time+1
			doAVAR(argv[i], tcol - 1, dcol - 1, tmult, numtaus, mintau, maxtau, logtau, overlap, jitter);	//basecol input is 1-based
			break;
		}
	}
	
	return 0;
}

void doAVAR(char *infile, int tcol, int dcol, double tmult, int numtaus, double mintau, double maxtau, bool logtau, bool overlap, double jitter) {
	deque<timedata> dataArray, avgData;
	deque<timedata>::iterator datIt;
	double curVal, curtau = 0, aval;
	
	//Split filename into base and ext. for output
	char outbase[511], *outext, outfile[511];
	strcpy(outbase, infile);
	outext = strrchr(outbase, '.');		//".ext"
	outext[0] = '\0';	//Change the "." into a terminator
	outext++;			//"ext"
	sprintf(outfile, "%s_avar.%s", outbase, outext);
	
	//Open input file
	printf("Input: %s\n", infile);
	ifstream inStream(infile, ios::in);
	
	//Read input
	timedata samp;
	char line[511], *column;
	unsigned long i=1;
	printf("Reading data...");
	while (!inStream.eof()) {
		//Line by line
		inStream.getline(line, 511);
		//Skip comment
		if (line[0] == '#')
			continue;
		
		//Parse line into columns
		int j=0;
		column = strtok(line, "\t");
		samp.data = -DBL_MAX; samp.time = -DBL_MAX;
		while (column != NULL) {
			if (column[0] > '9') {	//Bad data
				j++;
				column = strtok(NULL, "\t");
				continue;
			} else
				curVal = atof(column);
			
			if (j == tcol)
				samp.time = curVal * tmult;
			if (j == dcol)
				samp.data = curVal;

			j++;
			column = strtok(NULL, "\t");
		}
		
		if (tcol < 0)
			samp.time = (i-1) * tmult;

		if ((samp.data > -DBL_MAX) && (samp.time > -DBL_MAX))
			dataArray.push_back(samp);
		
		i++;
	}
	inStream.close();
	printf(" found %i datapoints\n", dataArray.size());
	
	//Check if we should even continue
	if (dataArray.size() <= MIN_SAMPS) {
		printf("Insufficient valid datapoints to calculate Allan variance for %s\n", infile);
		printf("Only found %i, need at least %i\n", dataArray.size(), MIN_SAMPS);
		return;
	}
	

	//Calculate all tau values we need
	//printf("Calculating tau values...\n");
	//Find the minimum and maximum time values, to create appropriate tau sizes
	double min, max;
	if (tcol > -1) {
		min = dataArray[0].time;
		max = dataArray.back().time;
		for (datIt = dataArray.begin(); datIt != dataArray.end(); datIt++) {
			min = _min(min, datIt->time);
			max = _max(max, datIt->time);
		}
	} else {	//Simple: just use index
		min = 0;
		max = (double)(dataArray.size() - 1)*tmult;
	}
	
	//Calculate tau extents
	if (!mintau)	//Skip if specified
		mintau = dataArray[2].time - min;	//Time extent of first 3 samples, i.e. time(4th) - time(1st)
	if (!maxtau)	//Skip if specified
		maxtau = (max-min)/MIN_SAMPS;
	printf("Tau extent: %f to %f\n", mintau, maxtau);
	

	//Open output file
	ofstream outFile(outfile, ios::trunc);
	outFile.precision(15);	//Increase precision to double max
	printf("Output: %s\n\n", outfile);
	
	//Write header
	if (overlap)
		outFile << "#Overlapped";
	else
		outFile << "#Non-overlapped";
	outFile << " Allan variance of " << infile << "\n";
	if (tcol > -1)
		outFile << "#Using column " << tcol+1;
	else
		outFile << "#Using datapoint number";
	if (tmult != 1)
		outFile << "*" << tmult;
	outFile << " as time, column " << dcol+1 << " as data\n";
	if (logtau)
		outFile << "#Using logarithmically spaced taus\n";
	else
		outFile << "#Using linearly spaced taus\n";
	outFile << "#tau\tavar\tadev\n";
	printf("\ntau\tAVAR\tADEV\n");

	//Now loop through all taus, calculating the average bins and Allan variance
	double curjitter = jitter*mintau;	//Allowed jitter
	for (double taunum = 1; taunum <= numtaus; taunum++) {
		//Determine tau
		if (numtaus < 2)	//Only one point
			curtau = mintau;
		else if (logtau)	//Use logarithmic spacing between minimum and maximum
			curtau = maxtau * pow( mintau/maxtau, (numtaus - taunum)/((double)numtaus - 1.0) );
		else	//Linear distribution across the available values
			curtau = ((taunum-1)/((double)numtaus-1.0))*(maxtau - mintau) + mintau;
		if (!overlap)
			curjitter = jitter*curtau;	//Jitter should be relative to bin size if we're not overlapping
		
		//Find the allan variance
		avgData = movingMean(&dataArray, curtau, overlap);
		if (avgData.size() < 2)
			continue;
		aval = avar(&avgData, curtau, jitter);
		avgData.clear();

		//Output
		printf("%f\t%f\t%f\n", curtau, aval, pow(aval, 0.5));
		outFile << curtau << '\t' << aval << '\t' << pow(aval, 0.5) << '\n';
	}

	outFile.close();

	return;
}

//Moving mean of a dataset, window size tau
deque<timedata> movingMean(deque<timedata> *input, double tau, bool overlap) {
	deque<timedata>::const_iterator datIt;
	deque<timedata> movingBin, meanVals;
	
	//Move a bin along the queue, add the averages
	for (datIt = input->begin(); datIt != input->end(); datIt++) {
		//If the bin is empty, add the first point and move on
		if (movingBin.empty()) {
			movingBin.push_back(*datIt);
			continue;
		}

		//Check if we should add this point
		if ((datIt->time - movingBin.front().time) < tau)	//We're still in the bin, add the point
			movingBin.push_back(*datIt);
		else {
			//No longer in the bin: get the last bin's average, move the bin forward one point
			meanVals.push_back(getMean(&movingBin));
			movingBin.pop_front();
			if (!overlap)
				movingBin.clear();	//Not using moving average, empty the bin and start a new one
			movingBin.push_back(*datIt);
		}
	}

	//Last bin
	if ((movingBin.back().time - movingBin.front().time) > 0.9*tau)	//Only add if it's a mostly-complete bin
		meanVals.push_back(getMean(&movingBin));

	return meanVals;
}

//Get the 2-value Allan variance of a prebinned dataset
double avar(deque<timedata> *data, double tau, double jitter) {
	deque<timedata>::const_iterator dIt, lastIt;
	double retVal = 0, diff;
	long count = 0;

	//Allan variance = (1/2) * < (first difference_i)^2 >

	//Get the squared differences between each pair of vectors
	lastIt = data->begin();
	for (dIt = data->begin(); dIt != data->end(); dIt++) {
		if (dIt->data == -DBL_MAX)	//Skip bad data
			continue;

		//Wait until we're tau apart
		if ((dIt->time - lastIt->time) < tau-jitter)	//Allow some bin jitter
			continue;

		//First difference
		diff = (lastIt->data - dIt->data);
		retVal += pow(diff, 2);
		count++;
		lastIt++;
	}

	if (!count)	//Abort if no valid data at all
		return 0;

	//Average
	retVal /= count;
	retVal /= 2;	//Definition has a 1/2

	return retVal;
}

//Get the average of a data series
double getMean(deque<double> *input) {
	double retVal=0;
	deque<double>::const_iterator vecIt;
	long count=0;
	
	for (vecIt = input->begin(); vecIt != input->end(); vecIt++) {
		if (*vecIt == -DBL_MAX)	//Ignore bad data
			continue;
		count++;
		retVal += *vecIt;
	}
	
	if (count == 0)
		retVal = -DBL_MAX;	//Indicate bad data
	else
		//Do the averaging
		retVal /= count;

	return retVal;
}
//Overloaded for timedata
timedata getMean(deque<timedata> *input) {
	timedata retVal;
	deque<timedata>::const_iterator vecIt;
	long count=0;
	
	retVal.time = 0; retVal.data = 0;
	for (vecIt = input->begin(); vecIt != input->end(); vecIt++) {
		if ((vecIt->data == -DBL_MAX) || (vecIt->time == -DBL_MAX))	//Ignore bad data
			continue;
		count++;
		retVal.data += vecIt->data;
		retVal.time += vecIt->time;
	}
	
	if (count == 0) {
		retVal.data = -DBL_MAX;	//Indicate bad data
		retVal.time = -DBL_MAX;
	} else {
		//Do the averaging
		retVal.data /= count;
		retVal.time /= count;
	}

	return retVal;
}


void showHelp() {
	printf(" alanvar\n Calculate the overlapping Alan Variance of a time series\n\n");
	printf(" usage:\n   alanvar [-options] <datafile(s)>\n   Eg. 'alanvar run1234.dat -logtau run5678.dat'\n   Will produce run1234_avar.dat and run5678_avar.dat\n\n");
	printf(" Options:\n");
	printf("  -time <tcol=1>     Will use 'tcol' as time (1-based, 0 for row #)\n");
	printf("  -data <dcol=2>     Will use 'dcol' as data (1-based, 0 for row #)\n");
	printf("  -tmult <mult=1.0>  Will multiply timebase by 'mult' before using\n");
	printf("  -n[umtaus] <n=100> Will calculate for n taus spaced from minimum to maximum\n");
	printf("  -logtau            Will logarithmically spaced taus\n");
	printf("  -mintau <min>      Set the minimum tau to be different than T(3)-T(1)\n");
	printf("  -maxtau <min>      Set the maximum tau to be different than (max time)/4\n");
	printf("  -nooverlap         Do not use overlapped Alan variance (moving averages)\n");
	printf("  -jitter <tol>      Allow 'tol' fractional jitter in bin separation\n");
	return;
}