/*
	Author:	Pate Williams (c) 2000

	MGIDA applied to static n-queens problem
	See "Solving Small and Large Scale Constraint
	Satisfaction Problems Using a Heuristic-Based
	Microgenetic Algorithm" by Gerry Dozier et al
*/

#include <assert.h>
#include <iomanip.h>
#include <iostream.h>
#include <limits.h>
#include <stdlib.h>
#include <time.h>

//#define DEBUG 1

const int Size = 4;	// size of arrays
const int MaxN = 1000;	// maximum number of queens
const int MaxBreakouts = 10 * MaxN;	// maximum number of breakouts
const int MaxCandSol = 100000;	// maximum candidate solutions
const int MaxPopulation = 1000;	// maximum population

int candidateSolutions, constraintChecks;

class Breakout {
private:
	int ltTuple, ltValue, rtTuple, rtValue, weight;
public:
	Breakout(void) { ltTuple = rtTuple = weight = - 1; };
	Breakout(int ltTup, int ltVal, int rtTup, int rtVal) {
		ltTuple = ltTup;
		ltValue = ltVal;
		rtTuple = rtTup;
		rtValue = rtVal;
	};
	Breakout(int ltTup, int ltVal, int rtTup, int rtVal, int w) {
		ltTuple = ltTup;
		ltValue = ltVal;
		rtTuple = rtTup;
		rtValue = rtVal;
		weight = w;
	};
	Breakout(Breakout &breakout) {
		ltTuple = breakout.ltTuple;
		ltValue = breakout.ltValue;
		rtTuple = breakout.rtTuple;
		rtValue = breakout.rtValue;
		weight = breakout.weight;
	};
	bool operator == (Breakout &breakout) {
		return ltTuple == breakout.ltTuple &&
			   ltValue == breakout.ltValue &&
			   rtTuple == breakout.rtTuple &&
			   rtValue == breakout.rtValue;
	};
	bool operator < (Breakout &breakout) {
		if (ltTuple < breakout.ltTuple)
			return true;
		if (ltTuple > breakout.ltTuple)
			return false;
		if (ltValue < breakout.ltValue)
			return true;
		if (ltValue > breakout.ltValue)
			return false;
		if (rtTuple < breakout.rtTuple)
			return true;
		if (rtTuple > breakout.rtTuple)
			return false;
		if (rtValue < breakout.rtValue)
			return true;
		return false;
	};
	int getLtTuple(void) { return ltTuple; };
	int getLtValue(void) { return ltValue; };
	int getRtTuple(void) { return rtTuple; };
	int getRtValue(void) { return rtValue; };
	int getWeight(void) { return weight; };
	void setLtTuple(int lt) { ltTuple = lt; };
	void setRtTuple(int rt) { rtTuple = rt; };
	void setWeight(int w) { weight = w; };
	void incrementWeight(void) { weight++; };
	friend ostream &operator << (ostream &output, const Breakout &breakout) {
		cout << breakout.ltTuple << ' ' << breakout.ltValue << ' '
			 << breakout.rtTuple << ' ' << breakout.rtValue << ' '
			 << breakout.weight;
		return output;
	};
};

class BreakoutList {
private:
	int length;
	Breakout breakout[MaxBreakouts];
public:
	BreakoutList(void) { length = 0; };
	int getLength(void) { return length; };
	void setLength(int len) { length = len; };
	void add(Breakout &bo) {
		add(bo.getLtTuple(), bo.getLtValue(),
			bo.getRtTuple(), bo.getRtValue(), bo.getWeight());
	};
	void add(int ltTuple, int ltValue, int rtTuple, int rtValue, int weight) {
		bool found;
		int i, index;
		Breakout bo(ltTuple, ltValue, rtTuple, rtValue, weight);

		assert(length < MaxBreakouts);
		for (found = false, i = index = 0; !found && i < length; i++) {
			if (bo < breakout[i]) {
				found = true;
				index = i;
			}
		}
		if (!found) {
			breakout[length++] = bo;
			return;
		}
		for (i = length; i >= index + 1; i--)
			breakout[i] = breakout[i - 1];
		breakout[index] = bo;
		length++;
	};
	int binarySearch(Breakout bo) {
		int lo = 0, hi = length - 1, mid;

		while (hi >= lo) {
			mid = (lo + hi) / 2;
			if (bo == breakout[mid])
				return mid;
			if (bo < breakout[mid])
				hi = mid - 1;
			else
				lo = mid + 1;
		}
		return - 1;
	};
	Breakout &operator [](int i) { return breakout[i]; };
	friend ostream &operator << (ostream &output, const BreakoutList &breakoutList) {
		for (int i = 0; i < breakoutList.length; i++)
			output << breakoutList.breakout[i] << endl; 
		return output;
	};
};

class Chromosome {
private:
	int breakouts, fitness, iterations, n, pivot, value;
	int cv[MaxN], hValue[MaxN], queen[MaxN];
public:
	Chromosome(void) {}
	bool violation(int i, int j) {
		int a = queen[i], b = queen[j];

		constraintChecks++;
		if (a == b)
			return true;
		if (abs(i - j) == abs(a - b)) 
			return true;
		return false;
	}
	int calculateFitness(void) {
		int c, i, j;
		
		for (i = fitness = 0; i < n; i++) {
			for (c = j = 0; j < n; j++)
				if (j != i && violation(i, j))
					c++;
			cv[i] = c;
			fitness += c;
		}
		fitness /= 2;
		reset();
		return fitness;
	}
	int calculateFitness(BreakoutList &list) {
		int c, i, j, wt;
		
		for (i = fitness = 0; i < n; i++) {
			for (c = j = 0; j < n; j++)
				if (j != i && violation(i, j))
					c++;
			cv[i] = c;
			fitness += c;
		}
		for (i = wt = 0; i < list.getLength(); i++) {
			Breakout breakout = list[i];
			if (queen[breakout.getLtTuple()] == breakout.getLtValue() &&
			    queen[breakout.getRtTuple()] == breakout.getRtValue() &&
				violation(breakout.getLtTuple(), breakout.getRtTuple()))
				wt += breakout.getWeight();
		}
		fitness = fitness / 2 + wt;
		reset();
		return fitness;
	}
	int init(int nn) {
		int i;
		
		breakouts = iterations = 0;
		n = nn;
		for (i = 0; i < n; i++)
			queen[i] = rand() % n;
		return calculateFitness();
	}
	void decHValue(void) {
		int array[MaxN], count, i, max, s;

		hValue[pivot]--;
		max = cv[0] + hValue[0];
		for (i = 1; i < n; i++) {
			s = cv[i] + hValue[i];
			if (s > max)
				max = s;
		}
		for (count = i = 0; i < n; i++)
			if ((cv[i] + hValue[i]) == max)
				array[count++] = i;
		pivot = array[rand() % count];
	}
	int evaluate(BreakoutList &list) {
		int i, c, t, wt;

		for (i = 0; i < n; i++)
			if (i != pivot && violation(pivot, i))
				cv[i]--;
		queen[pivot] = value;
		for (c = i = 0; i < n; i++) {
			if (i != pivot && violation(pivot, i)) {
				c++;
				cv[i]++;
			}
		}
		cv[pivot] = c;
		for (i = t = 0; i < n; i++)
			t += cv[i];
		for (i = wt = 0; i < list.getLength(); i++) {
			Breakout breakout = list[i];
			if (queen[breakout.getLtTuple()] == breakout.getLtValue() &&
			    queen[breakout.getRtTuple()] == breakout.getRtValue() &&
				violation(breakout.getLtTuple(), breakout.getRtTuple()))
				wt += breakout.getWeight();
		}
		fitness = t / 2 + wt;
		return fitness;
	}
	void mutate(void) {
		value = rand() % n;
	}
	int numberConflictingVariables(void) {
		int count, i;

		for (count = i = 0; i < n; i++)
			if (cv[i] != 0)
				count++;
		return count;
	}
	void reset(void) {
		int array[MaxN], count, i, max;
		
		for (i = 0; i < n; i++)
			hValue[i] = 0;
		max = cv[0];
		for (i = 1; i < n; i++)
			if (cv[i] > max)
				max = cv[i];
		for (count = i = 0; i < n; i++)
			if (cv[i] == max)
				array[count++] = i;
		pivot = array[rand() % count];
	}
	int getBreakouts(void) {
		return breakouts;
	}
	void setBreakouts(int b) {
		breakouts = b;
	}
	int &operator[](int i) {
		return queen[i];
	}
};

int fitness[MaxPopulation];
BreakoutList list;
Chromosome child, chromosome[MaxPopulation];

bool findSolution(int n, int population, int u) {
	bool equal;
	double rate;
	int childFitness, count, f, i, indx, j, max, min, p;
	int parent, parent0, parent1;
	int breakouts = 0, iterations = 0;
	int array[MaxPopulation];

	for (;;) {
		/*bool used[MaxPopulation] = {false};
		int tournamentSize = 3;
		int limit = tournamentSize < population ? tournamentSize : population;
		f = INT_MAX;
		for (i = 0; i < limit; i++) {
			do
				parent1 = rand() % population;
			while (used[parent1]);
			used[parent1] = true;
			if (fitness[parent1] < f) {
				f = fitness[parent1];
				parent = parent1;
			}
		}*/
		// binary tournament selection
		parent0 = rand() % population;
		parent1 = rand() % population;
		parent = fitness[parent0] < fitness[parent1] ? parent0 : parent1;
		// copy parent to child
		child = chromosome[parent];
		rate = 1.0 / (child.numberConflictingVariables() + u + 1);
		if ((double) rand() / RAND_MAX <= rate)
			child.reset();
		child.mutate();
		childFitness = child.evaluate(list);
		if (childFitness == 0) {
			child.setBreakouts(breakouts);
			return true;
		}
		if (fitness[parent] >= childFitness)
			child.decHValue();
		for (i = 0; i < population; i++) {
			equal = true;
			for (j = 0; equal && j < n; j++)
				equal = child[j] == chromosome[i][j];
			if (equal)
				break;
		}
		if (!equal) {
			max = fitness[0];
			for (i = 1; i < population; i++)
				if (fitness[i] > max)
					max = fitness[i];
			for (count = i = 0; i < population; i++)
				if (fitness[i] == max)
					array[count++] = i;
			i = array[rand() % count];
			chromosome[i] = child;
			fitness[i] = childFitness;
		}
		candidateSolutions++;
		if (candidateSolutions == MaxCandSol)
			break;
		// find minimum fitness
		min = fitness[0];
		for (i = 1; i < population; i++)
			if (fitness[i] < min)
				min = fitness[i];
		for (count = i = 0; i < population; i++)
			if (fitness[i] == min)
				array[count++] = i;
		min = chromosome[array[rand() % count]].numberConflictingVariables();
		iterations++;
		if (iterations >= min * n) {
			breakouts++;
			for (p = 0; p < population; p++) {
				for (f = i = 0; i < n - 1; i++) {
					for (j = i + 1; j < n; j++) {
						if (chromosome[p].violation(i, j)) {
							Breakout breakout(i, chromosome[p][i],
								j, chromosome[p][j], 1);
							indx = list.binarySearch(breakout);
							if (indx != - 1) {
								list[indx].incrementWeight();
								f += 1 + list[indx].getWeight();
							}
							else {
								list.add(breakout);
								f += 2;
							}
						}
					}
				}
				fitness[p] = f;
			}
			iterations = 0;
		}
	}
	return false;
}

bool simpleEA(int n, int population, int u, Chromosome &child) {
	int f, i;

	// initialize population of chromosomes to random domain values
	candidateSolutions = 0;
	constraintChecks = 0;
	list.setLength(0);	
	for (i = 0; i < population; i++) {
		f = fitness[i] = chromosome[i].init(n);
		candidateSolutions++;
		if (f == 0) {
			child.setBreakouts(0);
			child = chromosome[i];
			return true;
		}
	}
	return findSolution(n, population, u);
}

void print(int n, int *queen) {
	char hyphen[256];
	int column, i, i4, row;

	if (n <= 10) {
		for (i = 0; i < n; i++) {
			i4 = i * 4;
			hyphen[i4 + 0] = '-';
			hyphen[i4 + 1] = '-';
			hyphen[i4 + 2] = '-';
			hyphen[i4 + 3] = '-';
		}
		i4 = i * 4;
		hyphen[i4 + 0] = '-';
		hyphen[i4 + 1] = '\n';
		hyphen[i4 + 2] = '\0';
		for (row = 0; row < n; row++) {
			column = queen[row];
			cout << hyphen;
			for (i = 0; i < column; i++)
				cout << "|   ";
			cout << "| Q ";
			for (i = column + 1; i < n; i++)
				cout << "|   ";
			cout << '|' << endl;
		}
		cout << hyphen;
	}
	else
		for (row = 0; row < n; row++)
			cout << row << ' ' << queen[row] << endl;
}

int main(int argc, char *argv[]) {
	bool value;
	char name[8] = {"MGIDA"};
	double breakouts = 0, candSol = 0, cc = 0, runTime = 0, tm;
	int converged = 0, i, j, k, n, population, trials, u;
	clock_t clock0, clock1;
	time_t startTime = time(NULL), finishTime;

	srand(startTime);
#ifdef DEBUG
	if (argc != 3) {
		cout << "usage: " << argv[0] << " n population" << endl;
		exit(0);
	}
	int queen[MaxN];
	n = atoi(argv[1]);
	population = atoi(argv[2]);
	u = 1;
	simpleEA(n, population, u, child);
	for (i = 0; i < n; i++)
		queen[i] = child[i];
	print(n, queen);
	return 0;
#endif
	if (argc != 5) {
		cout << "usage: " << argv[0] << " n population trials u" << endl;
		cout << "where u is the heuristic reset rate parameter" << endl;
		exit(0);
	}
	n = atoi(argv[1]);
	if (n < 4 || n > MaxN) {
		cout << "error: 4 <= n <= " << MaxN << endl;
		exit(0);
	}
	population = atoi(argv[2]);
	trials = atoi(argv[3]);
	u = atoi(argv[4]);
	for (i = 0; i < trials; i++) {
		clock0 = clock();
		value = simpleEA(n, population, u, child);
		clock1 = clock();
		if (value) {
			converged++;
			breakouts += child.getBreakouts();
			candSol += candidateSolutions;
			cc += constraintChecks;
			tm = (double) (clock1 - clock0) / CLOCKS_PER_SEC;
			runTime += tm;
			for (j = 0; j < n - 1; j++) {
				for (k = j + 1; k < n; k++) {
					if (child.violation(j, k)) {
						cout << "constraint violation!" << endl;
						exit(0);
					}
				}
			}
		}
	}
	cout << "number of queens:             " << '\t' << n << endl;
	cout << "population:                   " << '\t' << population << endl;
	cout << "random number generator seed: " << '\t' << startTime << endl;
	cout << "number of trials:             " << '\t' << trials << endl;
	cout << "heuristic reset parameter:    " << '\t' << u << endl;
	cout << "algo\tcan sol\tbo\tcc\t%\ttime (s)" << endl;
	if (converged != 0)
		cout << name << '\t'
			 << candSol / converged << '\t'
			 << breakouts / converged << '\t'
			 << cc / converged << '\t'
			 << 100.0 * converged / trials << '\t'
			 << runTime / converged << endl;
	finishTime = time(NULL) - startTime;
	cout << "total time in seconds: " << finishTime << endl;
	return 0;
}