// Test of closest pair algorithms
// David Eppstein, UC Irvine, 19 Apr 1997
//
// Main entry point.
//
// Determines
//		how many points to generate
//		which point set generator to use
//		which algorithm to call on that point set
//		which closest pair data structure to use in the algorithm
//
// Then puts them all together and times the results.
// The primitive user interface is hopefully more portable than a Mac dialog.

#define MAX_LINE_LENGTH 1024		// how big to make buf for tty input?

#include "AllClosestPairs.h"		// find all implementations of closest pairs
#include "AllPointSets.h"			// find all point set generators
#include "AllAlgorithms.h"			// find all applications of closest pairs
#include "Error.h"
#include "Random.h"
#include <time.h>
#include <iostream.h>

unsigned long gDistances = 0;
unsigned long gPairs = 0;
unsigned long gInsertions = 0;
unsigned long gDeletions = 0;
double gTotalTime = 0.0;
double gTimeSquare = 0.0;
int gDimension = 0;

int gNegating = 0;

clock_t clock_value;

void StartTiming()
{
	clock_value = clock();
}

void EndTiming()
{
	clock_value = clock() - clock_value;
	double DoubleTime = (double) clock_value / (double) CLOCKS_PER_SEC;
	gTotalTime += DoubleTime;
	gTimeSquare += DoubleTime*DoubleTime;
}

void ReportOperationCounts()
{
	cout << "\n" << gInsertions << " point insertions.\n"
				 << gDeletions << " point deletions.\n"
				 << gDistances << " distance computations.\n"
				 << gPairs << " closest pair computations.\n"
				 << "Total (wall clock) time: "
				 << (double) clock_value / (double) CLOCKS_PER_SEC << "s.\n";
}

char * NextWord()
{
	static char buffer[MAX_LINE_LENGTH];
	cin >> buffer;
	if (!cin.good()) error("unexpected input termination");
	return buffer;
}

unsigned long ReadNum(const char * prompt)
{
	cout << prompt;
	char * s = NextWord();
	long retval = 0;
	while (s != 0 && *s != 0) {
		if (*s < '0' || *s > '9') {
			cout << "Unrecognized number format.\n\n";
			return ReadNum(prompt);			
		}
		retval = (retval * 10) + *s++ - '0';
	}
	return retval;
}

unsigned long ReadDimension()
{
	while (gDimension == 0) gDimension = ReadNum("Dimension: ");
	return gDimension;
}

PointSet * FindPointSet(unsigned long np)
{
	static char c = 0;
	if (c == 0) {
		cout << "\nChoose a point set generation method.\n"
			 << "  1 = points in R^d, L1 metric\n"
			 << "  2 = points in R^d, (squared) L2 metric\n"
			 << "  d = points in R^d, dot product\n"
			 << "  i = points in R^d, L(inf) metric\n"
			 << "  r = random distance matrix\n"
			 << "  s = points in generalized Sierpinski tetrahedron\n"
			 << "  t = two-adic non-archimedean metric\n"
			 << "Your choice: ";
		char * s = NextWord();
		if (s != 0 && *s != 0 && s[1] == 0) c = s[0];
	}
	switch(c) {
		case '1':
			return new VectorL1(np, ReadDimension());

		case '2':
			return new VectorL2(np, ReadDimension());

		case 'd': case 'D':
			return new VectorDot(np, ReadDimension());

		case 'i': case 'I':
			return new VectorLinf(np, ReadDimension());

		case 'r': case 'R':
			return new RandomDistance(np);
			
		case 's': case 'S':
			return new SierpinskiTetrahedron(np, ReadDimension());

		case 't': case 'T':
			return new TwoAdic(np);

		default:
			c = 0;
			cout << "Unrecognized point set generation option.\n\n";
			return FindPointSet(np);
	}
}

ClosestPairs * FindClosestPairs(unsigned long np, unsigned long mp, Distance & d)
{
	static char c = 0;
	static unsigned long nconga;
	static int read_nconga = 0;
	if (c == 0) {
		cout << "\nChoose a closest pair data structure.\n"
			 << "  b = brute force\n"
			 << "  c = conga line\n"
			 << "  f = fastpair conga-neighbor hybrid\n"
			 << "  m = many-subset conga line\n"
			 << "  n = nearest neighbor heuristic\n"
			 << "  q = quadtree\n"
			 << "Your choice: ";
		char * s = NextWord();
		if (s != 0 && *s != 0 && s[1] == 0) c = s[0];
	}
	switch(c) {
		case 'b': case 'B':
			StartTiming();
			return new BruteForceCP(np, mp, d);

		case 'c': case 'C':
			if (!read_nconga) {
				nconga = ReadNum("Number of subsets (zero for log_2(n)): ");
				read_nconga = 1;
			}
			StartTiming();
			return new CongaLine(np, mp, d, nconga);

		case 'f': case 'F':
			StartTiming();
			return new FastPair(np, mp, d);

		case 'm': case 'M':
			StartTiming();
			return new MultiConga(np, mp, d);

		case 'n': case 'N':
			StartTiming();
			return new NeighborCP(np, mp, d);

		case 'q': case 'Q':
			StartTiming();
			return new QuadTreeCP(np, mp, d);
			
		default:
			cout << "Unrecognized data structure option.\n\n";
			c = 0;
			return FindClosestPairs(np, mp, d);
	}
}

CPApplication * FindApplication()
{
	static char c = 0;
	if (c == 0) {
		cout << "\nChoose a closest pair application.\n"
			 << "  g = greedy matching\n"
			 << "  h = hierarchical clustering\n"
			 << "  i = cheapest insertion tsp\n"
			 << "  m = multifragment tsp\n"
			 << "  r = ray-intersection diagram\n"
			 << "Use a minus sign (e.g. \"-g\") for maximization versions\n"
			 << "Your choice: ";
		char * s = NextWord();
		if (s != 0 && *s == '-') {
			gNegating = 1;
			s++;
		} else gNegating = 0;
		if (s != 0 && *s != 0 && s[1] == 0) c = s[0];
	}
	switch(c) {
		case 'g': case 'G':
			return GreedyMatching;

		case 'h': case 'H':
			return Cluster;

		case 'i': case 'I':
			return CheapestInsertion;

		case 'm': case 'M':
			return MultiFragment;
			
		case 'r': case 'R':
			return RayDiagram;

		default:
			cout << "Unrecognized application option.\n\n";
			c = 0;
			return FindApplication();
	}
}

void main()
{
	unsigned long seed = ReadNum("Random number seed: ");
	SeedRandom(seed);
	int num_iters = ReadNum("Number of iterations: ");
	unsigned long num_pts = ReadNum("Number of points: ");
	while (num_pts == 0) {
		cout << "Nonzero number of points expected.\n\n";
		num_pts = ReadNum("Number of points: ");
	}
	for (int i = 0; i < num_iters; i++) {
		gInsertions = gDeletions = gPairs = gDistances = 0;
		CPApplication * app = FindApplication();

		PointSet * ps;
		if (app == RayDiagram) ps = new Tron(num_pts);
		else ps = FindPointSet(num_pts);
		if (gNegating) ps = new Negation(num_pts, *ps);
		if (app == MultiFragment) ps = new MultiFragmentDistance(num_pts, *ps);
		else if (app == CheapestInsertion)
			ps = new CheapestInsDistance(num_pts, *ps);

		ClosestPairs * p = FindClosestPairs(num_pts, num_pts, *ps);
		
		// do it!
		double total = (*app)(num_pts, *ps, *p);
		if (gNegating) total = -total;
		EndTiming();
		ReportOperationCounts();
		cout << "Total weight = " << total << ".\n";
		
		delete p;
		delete ps;
	}
	if (num_iters > 1)
		cout << "\nAverage time = " << gTotalTime / num_iters << "s.\n"
			 << "Standard deviation = "
			 << sqrt((gTimeSquare - (gTotalTime*gTotalTime)/num_iters)
			 		 / ((double) (num_iters - 1))) << "s.\n";
}