/* Copyright (C) 2000  Dale Schuurmans, Finnegan Southey                     */
/* This work is licensed under the Gnu General Public License (see gpl.txt). */

#include <stdlib.h>
#include <stdio.h>
#include "churn.h"

#define USE_CHURN_SANITY_CHECKS 1

int main(int argc, char **argv)
	{
	FILE *in, *out;
	int totalRepeats;
	int numRepeats;
	int numVars;
	int numAss;
	int maxLength;
	ATYPE *initialAss;
	ATYPE *chunk;
	WINTYPE *window = NULL;
	WINTYPE *windowAvg = NULL;
	int targetWindowLength = -1;
	int windowLength = -1;
	int gap, diameter;
	double assAvg;
	double covRate;
	double gapAvg, diameterAvg;
	double velocityAvg, velocityAvgAvg;
	double covRateAvg;
	int repeatCount;
	int width;
	char outName[2000];
	int doAutoCorr;

	if (argc < 7)
		{
		fprintf(stderr, "Usage: <ass file> <output window name> <max processing length> <do auto-corr (0/1)> <window length> <repeats>\n");
		exit(1);
		}

	in = fopen(argv[1], "r");
	if (in == NULL)
		{
		fprintf(stderr, "Error opening assignment file: %s\n", argv[1]);
		exit(1);
		}

	if (sscanf(argv[3], "%d", &maxLength) != 1)
		{
		fprintf(stderr, "Must specify maximum length.\n");
		exit(1);
		}
	
	if (sscanf(argv[4], "%d", &doAutoCorr) != 1)
		{
		fprintf(stderr, "Must specify 0 or 1 indicating whether auto-correlation should be computed.\n");
		exit(1);
		}

	if (sscanf(argv[5], "%d", &targetWindowLength) != 1)
		{
		fprintf(stderr, "Must specify the window length for the auto-correlation.\n");
		exit(1);
		}

	if (sscanf(argv[6], "%d", &numRepeats) != 1)
		{
		fprintf(stderr, "Must specify the maximum number of repeats to be processed (0 == all of them).\n");
		exit(1);
		}

	readAssFileHeader(in, &totalRepeats, &numVars);

	if (numRepeats == 0)
		numRepeats = totalRepeats;

	if (totalRepeats < numRepeats)
		numRepeats = totalRepeats;

	initialAss = malloc(numVars * sizeof(ATYPE));
	if (initialAss == NULL)
		{
		fprintf(stderr, "Failed to allocate buffer for storing intial assignments.\n");
		exit(1);
		}

	if (doAutoCorr)
		{
		if (numRepeats > 1)
			{
			windowAvg = calloc(maxLength, sizeof(WINTYPE));
			if (windowAvg == NULL)
				{
				fprintf(stderr, "Failed to allocate buffer for storing window average.\n");
				exit(1);
				}
			}
		}

	assAvg = 0;
	gapAvg = 0;
	diameterAvg = 0;
	velocityAvgAvg = 0;
	covRateAvg = 0;

	fprintf(stdout, "#Flips\tGap\tDiam\tVelAvg\tCovRate\n");
	for (repeatCount = 0; repeatCount < numRepeats; repeatCount++)
		{
		chunk = malloc(maxLength * sizeof(ATYPE));
		if (chunk == NULL)
			{
			fprintf(stderr, "Failed to allocate buffer for storing flips.\n");
			exit(1);
			}
	
		readAssFileChunk(in, numVars, maxLength, chunk, initialAss, &numAss);
	
		if (doAutoCorr)
			{
			if (targetWindowLength > numAss)
				windowLength = numAss;
			else
				windowLength = targetWindowLength;
	
			window = calloc(windowLength, sizeof(WINTYPE));
			if (window == NULL)
				{
				fprintf(stderr, "Failed to allocate buffer for storing window.\n");
				exit(1);
				}
			}
	
		mangle(numVars, numAss, chunk, initialAss, &gap, &diameter, &velocityAvg, &covRate, doAutoCorr, windowLength, window);
	
		fprintf(stdout, "%d\t%d\t%d\t%g\t%g\n", numAss, gap, diameter, velocityAvg, covRate);
		assAvg += numAss;
		gapAvg += gap;
		diameterAvg += diameter;
		velocityAvgAvg += velocityAvg;
		covRateAvg += covRate;

		if (doAutoCorr)
			{
			sprintf(outName, "%s.%d", argv[2], repeatCount);

			out = fopen(outName, "w");
			if (out == NULL)
				{
				fprintf(stderr, "Error opening output window file: %s\n", argv[2]);
				exit(1);
				}
	
			for (width = 0; width < windowLength; width++)
				fprintf(out, "%g\n", window[width]);

			fclose(out);
	
			if (numRepeats > 1)
				for (width = 0; width < windowLength; width++)
					windowAvg[width] += window[width];

			free(window);
			}
	
		free(chunk);
		}

	assAvg /= numRepeats;
	gapAvg /= numRepeats;
	diameterAvg /= numRepeats;
	velocityAvgAvg /= numRepeats;
	covRateAvg /= numRepeats;
	
	fprintf(stdout, "\n%g\t%g\t%g\t%g\t%g\n", assAvg, gapAvg, diameterAvg, velocityAvgAvg, covRateAvg);

	if (doAutoCorr && numRepeats > 1)
		{
		sprintf(outName, "%s.avg", argv[2]);

		out = fopen(outName, "w");
		if (out == NULL)
			{
			fprintf(stderr, "Error opening output window avg file: %s\n", outName);
			exit(1);
			}
	
		for (width = 0; width < windowLength; width++)
			{
			windowAvg[width] /= numRepeats;
			fprintf(out, "%g\n", windowAvg[width]);
			}

		fclose(out);

		free(windowAvg);
		}
	
	free(initialAss);

	fclose(in);

	return 0;
	}

void mangle(int numVars, int numAss, ATYPE *chunk, ATYPE *initialAss, int
		*gapPtr, int *diameterPtr, double *velocityAvgPtr, double *covRatePtr,
		int doAutoCorr, int windowLength, WINTYPE *window)
	{
	int minCompDist;
	int currentDiameter = -1;
	int currentGap = 0;
	int outer, inner, var;
	int dist, compDist;
	ATYPE *outerPtr;
	ATYPE *innerPtr;
	ATYPE *outerAss;
	ATYPE *innerAss;
	int maxTravelDist;
	double vel, maxVel, velSum;
	int *windowAccum = NULL;
	int *windowObs = NULL;
	int lag;

	outerAss = malloc(numVars * sizeof(ATYPE));
	innerAss = malloc(numVars * sizeof(ATYPE));
	
	if (doAutoCorr)
		{
		if (windowLength > numAss)
			windowLength = numAss;

		windowAccum = calloc(windowLength, sizeof(int));
		if (windowAccum == NULL)
			{
			fprintf(stderr, "Failed to allocate buffer for storing accumulators.\n");
			exit(1);
			}
	
		windowObs = calloc(windowLength, sizeof(int));
		if (windowObs == NULL)
			{
			fprintf(stderr, "Failed to allocate buffer for storing observation counts.\n");
			exit(1);
			}
		}

	velSum = 0;

	for (var = 0; var < numVars; var++)
		outerAss[var] = initialAss[var];
	
	outerPtr = chunk;

		/* For each assignment, x. */
	for (outer = 0; outer < numAss; outer++)
		{
		minCompDist = numVars;
		maxTravelDist = 0;
		maxVel = 0;

			/* Initialize inner assignment and distance measure. */
		dist = 0;
		for (var = 0; var < numVars; var++)
			{
			innerAss[var] = initialAss[var];
			if (innerAss[var] != outerAss[var])
				dist++;
			}

		
		innerPtr = chunk;

			/* For each assignment, y. */
		for (inner = 0; inner < numAss; inner++)
			{
#if USE_CHURN_SANITY_CHECKS
			if (dist < 0)
				{
				fprintf(stdout, "BADDIST: Negative distance!\n");
				exit(1);
				}
			if (dist > numVars)
				{
				fprintf(stdout, "BADDIST: Excessive distance!\n");
				exit(1);
				}
#endif			
			lag = inner - outer;

			compDist = numVars - dist;
	
				/* Maintain minimum compdist (candidate for gap). */
			if (compDist < minCompDist)
				minCompDist = compDist;


				/* Everything under this condition only needs to be compared against
				 * all previous assignments, and not against all future assignments.
				 */
			if (inner < outer)
				{
					/* Maintain maximum dist (diameter). */
				if (dist > currentDiameter)
					currentDiameter = dist;
				}

				/* Compute velocity. */
			if (inner > outer && lag < numVars)
				{
				if (dist > maxTravelDist)
					maxTravelDist = dist;

				vel = (double)(maxTravelDist - dist) / lag;

				if (vel > maxVel)
					maxVel = vel;
				}

			if (doAutoCorr)
				{
					/* Accumulate for auto-correlation. */
				if (inner > outer && lag < windowLength)
					{
					windowAccum[lag] += dist;
					windowObs[lag]++;
					}
				}
			
				/* Make inner flip and update distance measure. */
			if (inner < numAss - 1)
				{
				innerAss[*innerPtr] *= -1;
				if (innerAss[*innerPtr] == outerAss[*innerPtr])
					dist--;
				else
					dist++;

				innerPtr++;
				}
			}

			/* If minCompDist > currentGap, update currentGap. */
		if (minCompDist > currentGap)
			currentGap = minCompDist;

			/* Accumulate velocity scores for averaging. */
		velSum += maxVel;

			/* Make outer flip. */
		if (outer < numAss - 1)
			{
			outerAss[*outerPtr] *= -1;
			outerPtr++;
			}
		}

	if (doAutoCorr)
		{
			/* Compute averages for window. */
		for (outer = 1; outer < windowLength; outer++)
			{
			window[outer] = (WINTYPE)windowAccum[outer] / windowObs[outer];
			}

		free(windowAccum);
		free(windowObs);
		}

	*gapPtr = currentGap;
	*diameterPtr = currentDiameter;
	*velocityAvgPtr = velSum / numAss;
	*covRatePtr = ((double)numVars - currentGap) / numVars / numAss;
	}

void readAssFileChunk(FILE *in, int numVars, int maxLength, ATYPE *chunk, ATYPE *initialAssignments, int *numAssPtr)
	{
	ATYPE *chunkPtr;
	ATYPE val;
	int varCount;
	int assCount;

	assCount = 1;
	chunkPtr = chunk;

	for (varCount = 0; varCount < numVars; varCount++)
		fscanf(in, " %d ", initialAssignments + varCount);

	fscanf(in, " %d ", &val);

	while (val != -1 && assCount < maxLength)
		{
		*chunkPtr = val;
		chunkPtr++;
		assCount++;
		fscanf(in, " %d ", &val);
		}

		/* Throw away remainder of chunk. */
	while (val != -1)
		fscanf(in, " %d ", &val);

	*numAssPtr = assCount;
	}

void readAssFileHeader(FILE *in, int *numRepeatsPtr, int *numVarsPtr)
	{
	int ch;

		/* Ditch comments. */
	ch = fgetc(in);
	while (ch == '#')
		{
		while (fgetc(in) != '\n');
		ch = fgetc(in);
		}

	if (ch != 'r')
		{
		fprintf(stderr, "Was expecting repeat count (e.g. 'r 5').\n");
		exit(1);
		}

	if (fscanf(in, " %d ", numRepeatsPtr) < 0)
		{
		fprintf(stderr, "Was expecting repeat count (e.g. 'r 5').\n");
		exit(1);
		}

	ch = fgetc(in);

	if (ch != 'v')
		{
		fprintf(stderr, "Was expecting var count (e.g. 'v 5').\n");
		exit(1);
		}

	if (fscanf(in, " %d ", numVarsPtr) < 0)
		{
		fprintf(stderr, "Was expecting var count (e.g. 'v 5').\n");
		exit(1);
		}
	}