/*
	Author:	James Pate Williams (c) 2000
	
	Simulated annealing of a combinatorial objective function.
	n-queens problem
*/

#include <fstream.h>
#include <iomanip.h>
#include <iostream.h>
#include <limits.h>
#include <math.h>
#include <stdlib.h>
#include <time.h>

extern "C" void sgenrand(unsigned long seed);
extern "C" unsigned long genrand(void);

const int MaxN = 10000;	// maximum number of variables

double doubleRand(void) {
	// random number in the closed interval [0, 1]
	return (double) genrand() / ULONG_MAX;
}

int intRand(int modulus) {
	// returns a random integer in the closed interval [0, modulus - 1]
	int r = genrand() % modulus;

	if (r < 0)
		r += modulus;
	return r;
}

inline
bool violation(int i, int j, int queen[MaxN]) {
	int a = queen[i], b = queen[j];

	if (a == - 1 || b == - 1)
		return false;
	if (a == b)
		return true;
	if (abs(i - j) == abs(a - b)) 
		return true;
	return false;
}

inline
int f(int N, int cv[MaxN], int queen[MaxN]) {
	// calculate the total number of constraint violations
	int c, i, j, v = 0;

	for (i = 0; i < N; i++) {
		for (c = 0, j = 0; j < N; j++)
			if (i != j && violation(i, j, queen))
				c++;
		cv[i] = c;
		v += c;
	}
	return v;
}

void simulatedAnnealing(double alpha, double T0, int N,
						int m, int n, int aF[MaxN],	int A0[MaxN],
						int &fF, int &count1, int &count2,
						int &count3, int &count4) {
	// algorithm by Professor and Chair Alice E. Smith
	// alpha = exponential annealing schedule parameter in (0, 1)
	// standard = standard deviation for randomGaussian function
	// T0 = initial temperature
	// N = number of facilities and sites
	// m = number of temperatures
	// n = number of moves at a given temperature
	// aF = final value of assignment vector
	// A0 = initial value of assignment vector
	// fF = final value of function f
	// count1 = number of descending moves
	// count2 = number of Metropolis accepted uphill moves
	// count3 = number of rejected Metropolis moves
	// count4 = number of changes in final function values
	// count1 + count2 + count3 = m * n
	bool found;
	double T;
	int ai, c, f0, f1, i, j, t, u, uTemp, v, vTemp; 
	int a0[MaxN], aTemp[MaxN], cv0[MaxN], tempCV[MaxN];

	// select random starting point (random permutation of 0 to N - 1)
	for (i = 0; i < N; i++) {
		do {
			ai = intRand(N);
			for (found = false, j = 0; !found && j < i; j++)
				found = ai == a0[j];
		} while (found);
		A0[i] = ai;
		a0[i] = ai;
		aF[i] = ai;
	}
	T = T0;
	count1 = count2 = count3 = count4 = 0;
	f0 = fF = f(N, cv0, a0);
	if (fF == 0)
		return;
	for (t = 1; t <= m; t++) {
		for (i = 1; i <= n; i++) {
			// calculate move (swap)
			do {
				u = intRand(N);
				v = intRand(N);
			} while (u == v);
			// copy a to aTemp
			for (j = 0; j < N; j++) {
				aTemp[j] = a0[j];
				tempCV[j] = cv0[j];
			}
			uTemp = a0[u];
			vTemp = a0[v];
			// propagate out-ness
			for (j = 0; j < N; j++) {
				if (j != u && violation(j, u, a0)) {
					tempCV[j]--;
					tempCV[u]--;
				}
			}
			a0[u] = - 1;
			// propagate out-ness
			for (j = 0; j < N; j++) {
				if (j != v && violation(j, v, a0)) {
					tempCV[j]--;
					tempCV[v]--;
				}
			}
			a0[u] = uTemp;
			// propagate in-ness
			aTemp[u] = vTemp;
			aTemp[v] = - 1;
			for (c = j = 0; j < N; j++) {
				if (j != u && violation(j, u, aTemp)) {
					c++;
					tempCV[j]++;
				}
			}
			tempCV[u] = c;
			// propagate in-ness
			aTemp[v] = uTemp;
			for (c = j = 0; j < N; j++) {
				if (j != v && violation(j, v, aTemp)) {
					c++;
					tempCV[j]++;
				}
			}
			tempCV[v] = c;
			// calculate f1
			for (f1 = j = 0; j < N; j++)
				f1 += tempCV[j];
			if (f1 <= f0) {
				// update a0 and f0
				for (j = 0; j < N; j++) {
					a0[j] = aTemp[j];
					cv0[j] = tempCV[j];
				}
				f0 = f1;
				count1++;
			}
			else {
				// Metropolis part of algorithm
				if (doubleRand() <= exp(- (f1 - f0) / T)) {
					// accept uphill move
					// update a0 and f0
					for (j = 0; j < N; j++) {
						a0[j] = aTemp[j];
						cv0[j] = tempCV[j];
					}
					f0 = f1;
					count2++;
				}
				else {
					// reject uphill move
					count3++;
				}
			}
			if (f1 <= fF) {
				// update final assignment vector and final f
				for (j = 0; j < N; j++)
					aF[j] = aTemp[j];
				fF = f1;
				count4++;
				if (fF == 0)
					return;
			}
		}
		// exponential annealing schedule
		T = alpha * T;
	}
}

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

	if (n <= 12) {
		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) {
	double alpha, T0;
	int count1, count2, count3, count4, fF;
	int N, m = 10000, n = 100;
	int aF[MaxN], a0[MaxN];
	unsigned long seed;
	time_t time0 = time(NULL);

	if (argc != 7) {
		cout << "usage: " << argv[0] << " alpha T0 m n N seed" << endl;
		cout << "where alpha is exponential annealing parameter" << endl;
		cout << "T0 is the initial high temperature" << endl;
		cout << "m is the number of temperatures" << endl;
		cout << "n is the number of moves at a given temperature" << endl;
		cout << "N is the number of queens" << endl;
		cout << "seed is a nonzero counting number" << endl;
		exit(0);
	}
	alpha = atof(argv[1]);
	T0 = atof(argv[2]);
	m = atoi(argv[3]);
	n = atoi(argv[4]);
	N = atoi(argv[5]);
	seed = atol(argv[6]);
	if (alpha <= 0.0 || alpha >= 1.0) {
		cout << "alpha must be in the open interval (0, 1)" << endl;
		exit(0);
	}
	if (T0 <= 0) {
		cout << "error: T0 must be nonzero positive number" << endl;
		exit(0);
	}
	if (m <= 0) {
		cout << "error: m must be a nonzero counting number" << endl;
		exit(0);
	}
	if (n <= 0) {
		cout << "error: n must be a nonzero counting number" << endl;
		exit(0);
	}
	if (N < 4 || N > MaxN) {
		cout << "error: 4 <= N <= " << MaxN << endl;
		exit(0);
	}
	if (seed <= 0) {
		cout << "error: seed must be nonzero counting number" << endl;
		exit(0);
	}
	sgenrand(seed);
	simulatedAnnealing(alpha, T0, N, m, n, aF, a0, fF,
		count1, count2, count3, count4);
	cout << "alpha = " << alpha << endl;
	cout << "T0 = " << T0 << endl;
	cout << "N = " << N << endl;
	cout << "m = " << m << endl;
	cout << "n = " << n << endl;
	cout << "seed = " << seed << endl;
	cout << "f = " << fF << endl;
	cout << "count1 = " << count1 << endl;
	cout << "count2 = " << count2 << endl;
	cout << "count3 = " << count3 << endl;
	cout << "count4 = " << count4 << endl;
	if (N <= 12) {
		cout << endl << "the initial board is:" << endl << endl;
		print(N, a0);
		cout << endl;
		cout << "the final board is:" << endl << endl;
		print(N, aF);
		cout << endl;
	}
	cout << "total time in seconds = "
		 << time(NULL) - time0 << endl;
	return 0;
}