/************************************************************************
 *                                                                      *
 *      AN EXPERIMENTAL TABU SEARCH CODE FOR THE N-QUEENS PROBLEM       *
 *                                                                      *
 *               Copyright (C) 1992  by  Manuel Laguna                  *
 *                                                                      *
 *  Last revision: May 4, 1992.                          Version: 1.1   *
 *                                                                      *
 *  TSQUEEN.C :   <n>          number of queens                         *
 *                <tabu_size>  size of the tabu list                    *
 *                <max_iter>   maximum number of iterations             *
 *                <seed>       seed for the random number generator     *
 *                                                                      * 
 ************************************************************************/

#include <stdio.h>

#define TRACE 1     /* set = 1 to display information every iteration   */
#define ZERO  0     /* set = 1 to disallow zero moves                   */
#define FIS   0     /* set = 1 to implement first-improving strategy    */
#define MAX(X,Y)    ( (X) > (Y) ? (X) : (Y) )

struct move_info {
	int i;
	int j;
	int value;
};

/************************************************************************
 *                                                                      *
 *                       FUNCTION PROTOTYPES                            *
 *                                                                      *
 ************************************************************************/

void tabu_search(int n, int tabu_size, int max_iter);
int initialization(int n, int *pi, int *d1, int *d2, int **tabu_time);
void indexx(int n, int *arrin, int *indx);
void update_best(int n, int *pi_best, int *pi_curr, int *c_best, int c_curr);
void best_move(int n, int iter, int c_best, int c_curr, int *pi, int *d1,
int *d2, int **tabu_time, struct move_info *move);
int evaluate(int i, int j, int *pi, int *d1, int *d2);
void execute_move(struct move_info move, int *pi, int *d1, int *d2);
void print_solution(int flag, int n, int *pi, int c, int iter);

/***************** Memory Allocation Prototypes *************************/

int *ivector(int nl, int nh);
void free_ivector(int *v, int nl);
int **imatrix(int nrl, int nrh, int ncl, int nch);
void free_imatrix(int **m, int nrl, int nrh, int ncl);
void nrerror(char *error_text);

/************************************************************************
 *                                                                      *
 *                         MAIN FUNCTION                                *
 *                                                                      *
 ************************************************************************/

void main(int argc, char **argv)

{
	int n;                  /* number of queens                         */
	int tabu_size;          /* size of the short term memory function   */
	int max_iter;           /* maximum number of iterations             */
	int seed;               /* seed for random number generator         */

	if (argc != 5) nrerror("TSQUEEN: <n> <tabu_size> <max_iter> <seed>");
	n = atoi(argv[1]);
	tabu_size = atoi(argv[2]);
	max_iter = atoi(argv[3]);
	seed = atoi(argv[4]);
	srandom(seed*314);
	tabu_search(n, tabu_size, max_iter);
}

/************************************************************************
 *                                                                      *
 *                           TABU SEARCH                                *
 *                                                                      *
 ************************************************************************/

void tabu_search(int n, int tabu_size, int max_iter)

{
	int iter;                  /* iteration counter                     */
	int iter_best;             /* iteration where best was found        */
	int *d1;                   /* queens in negative diagonals          */
	int *d2;                   /* queens in positive diagonals          */
	int *pi_curr;              /* current queen configuration           */
	int *pi_best;              /* best queen configuration              */
	int c_curr;                /* number of collisions in pi_curr       */
	int c_best;                /* number of collisions in pi_best       */
	int **tabu_time;           /* tabu structure for swap moves         */
	
	struct move_info move;     /* information related to best move      */
	
	d1 = ivector(2, 2*n);
	d2 = ivector(-n+1, n-1);
	pi_curr = ivector(1, n);
	pi_best = ivector(1, n);
	tabu_time = imatrix(1, n, 1, n);

	do {	
		c_curr = initialization(n, pi_curr, d1, d2, tabu_time);
	} while (c_curr == 0);
	update_best(n, pi_best, pi_curr, &c_best, c_curr);
	iter_best = 0;
	iter = 0;
	print_solution(0, n, pi_best, c_best, iter_best);
	while (iter < max_iter && c_best > 0) {
		++iter;
		best_move(n, iter, c_best, c_curr, pi_curr, d1, d2,
        tabu_time, &move);
		execute_move(move, pi_curr, d1, d2);
		tabu_time[move.i][move.j] = iter + tabu_size;
		c_curr += move.value;
#if TRACE
		printf("SWAP(%3d, %3d) = %3d\n", move.i, move.j, move.value);
#endif
		if (c_curr < c_best) {
			update_best(n, pi_best, pi_curr, &c_best, c_curr);
			iter_best = iter;
		}
#if TRACE
	print_solution(1, n, pi_curr, c_curr, iter); 
#endif
	}
	print_solution(2, n, pi_best, c_best, iter_best); 
	free_ivector(d1, 2);
	free_ivector(d2, -n+1);
	free_ivector(pi_curr, 1);
	free_ivector(pi_best, 1);
	free_imatrix(tabu_time, 1, n, 1);
}

/************************************************************************
 *                                                                      *
 *                         INITIALIZATION                               *
 *                                                                      *
 ************************************************************************/
	
int initialization(int n, int *pi, int *d1, int *d2, int **tabu_time)

{
	int i;                  /* queen index                              */
	int *r;                 /* uniform random number                    */
	int collisions;         /* total number of collisions               */
	
	r = ivector(1, n);
	for (i = 1; i <= n; ++i) r[i] = random();
	indexx(n, r, pi); 
	for (i = 2; i <= 2*n; ++i) {
		d1[i] = 0;
		d2[i-n-1] = 0;
	}
	for (i = 1; i <= n; ++i) {
		++d1[i+pi[i]];
		++d2[i-pi[i]];
	}
	collisions = 0;
	for (i = 2; i <= 2*n; ++i) {
		collisions += MAX(0, d1[i]-1);
		collisions += MAX(0, d2[i-n-1]-1);
	}
	free_ivector(r, 1);
	return collisions;
}
		
/************************************************************************
 *                                                                      *
 *                  CONSTRUCTION OF AN INDEX TABLE                      *
 *                                                                      *
 ************************************************************************/

void indexx(int n, int *arrin, int *indx)

{
	int l, j, ir, indxt, i;
	int q;

	for (j = 1; j <= n; j++) indx[j] = j;
	if (n == 1) return;
	l = (n >> 1) + 1;
	ir = n;
	for (;;) {
		if (l > 1)
			q = arrin[(indxt = indx[--l])];
		else {
			q = arrin[(indxt = indx[ir])];
			indx[ir] = indx[1];
			if (--ir == 1) {
				indx[1] = indxt;
				return;
			}
		}
		i = l;
		j = l << 1;
		while (j <= ir) {
			if ( j < ir && arrin[indx[j]] < arrin[indx[j+1]]) j++;
			if (q < arrin[indx[j]]) {
				indx[i] = indx[j];
				j += (i = j);
			}
			else j = ir + 1;
		}
		indx[i] = indxt;
	}
}
		
/************************************************************************
 *                                                                      *
 *                      UPDATE BEST SOLUTION                            *
 *                                                                      *
 ************************************************************************/

void update_best(int n, int *pi_best, int *pi_curr, int *c_best, int c_curr)

{
	int i;                  /* queen index                              */

	for (i = 1; i <= n; ++i)
		pi_best[i] = pi_curr[i];
	*c_best = c_curr;
}

/************************************************************************
 *                                                                      *
 *                           BEST MOVE                                  *
 *                                                                      *
 ************************************************************************/

void best_move(int n, int iter, int c_best, int c_curr, int *pi, int *d1,
int *d2, int **tabu_time, struct move_info *move)

{
	int i;                  /* queen i index                            */
	int j;                  /* queen j index                            */
	int value;              /* value of move under consideration        */

	move->value = n + 1;
	for (i = 1; i <= n-1; ++i) {
		for (j = i+1; j <= n; ++j) {
			value = evaluate(i, j, pi, d1, d2);
#if ZERO
			if (value) {
#endif
			if (tabu_time[i][j] < iter || c_curr + value < c_best) {
				if (value < move->value) {
					move->value = value;
					move->i = i;
					move->j = j;
				}
			}
#if ZERO
			}
#endif
#if FIS
			if (move->value < 0) break; 
#endif
		}
#if FIS
		if (move->value < 0) break; 
#endif
	}
}

/************************************************************************
 *                                                                      *
 *                       MOVE EVALUATION                                *
 *                                                                      *
 ************************************************************************/

int evaluate(int i, int j, int *pi, int *d1, int *d2)

{
    int c;

    c = MAX(0, d1[i+pi[i]]-2) + MAX(0, d2[i-pi[i]]-2) - d1[i+pi[i]] -
        d2[i-pi[i]] + 2;
    if (i+pi[i] != j+pi[j]) c += MAX(0, d1[j+pi[j]]-2) - d1[j+pi[j]] + 1;
    if (i-pi[i] != j-pi[j]) c += MAX(0, d2[j-pi[j]]-2) - d2[j-pi[j]] + 1;
    if (i+pi[i] == j+pi[j] && d1[i+pi[i]] > 2) --c;
    if (i-pi[i] == j-pi[j] && d2[i-pi[i]] > 2) --c;
    c += d1[i+pi[j]] + d2[i-pi[j]] - MAX(0, d1[i+pi[j]]-1) -
         MAX(0, d2[i-pi[j]]-1);
    if (i+pi[j] != j+pi[i]) c += d1[j+pi[i]] - MAX(0, d1[j+pi[i]]-1);
    if (i-pi[j] != j-pi[i]) c += d2[j-pi[i]] - MAX(0, d2[j-pi[i]]-1);
    if (i+pi[j] == j+pi[i] || i-pi[j] == j-pi[i]) ++c;
    return c;
}

/************************************************************************
 *                                                                      *
 *                        EXECUTE MOVE                                  *
 *                                                                      *
 ************************************************************************/

void execute_move(struct move_info move, int *pi, int *d1, int *d2)

{
	int j;                  /* column index                             */
	
	--d1[move.i+pi[move.i]];
	--d2[move.i-pi[move.i]];
	--d1[move.j+pi[move.j]];
	--d2[move.j-pi[move.j]];
	++d1[move.i+pi[move.j]];
	++d2[move.i-pi[move.j]];
	++d1[move.j+pi[move.i]];
	++d2[move.j-pi[move.i]];
	j = pi[move.j];
	pi[move.j] = pi[move.i];
	pi[move.i] = j;
}
		
/************************************************************************
 *                                                                      *
 *                      PRINT BEST SOLUTION                             *
 *                                                                      *
 ************************************************************************/

void print_solution(int flag, int n, int *pi, int c, int iter)

{
	int i;                  /* queen index                              */
	int j;                  /* column index                             */

	if (flag == 0) printf("\nStarting solution:\n\n");
	if (flag == 1) printf("\nSolution found at iteration %4d:\n\n", iter);
	if (flag == 2) printf("\nBest solution found:\n\n");
	printf("   Number of collisions = %6d\n", c);
	if (flag == 2) printf("   Number of iterations = %6d\n", iter);
	printf("\n   ");
	for (j = 1; j <= n; ++j) printf("%3d", j);
	for (i = 1; i <= n; ++i) {
		printf("\n%3d", i);
		for (j = 1; j <= n; ++j)
			if (pi[i] == j) printf(" **"); else printf(" __");
	}
	printf("\n\n");
}

/************************************************************************
 *                                                                      *
 *                  MEMORY ALLOCATION FUNCTIONS                         *
 *                                                                      *
 ************************************************************************/

int *ivector(int nl, int nh)

{
	int *v;

	v = (int *)malloc((unsigned)(nh-nl+1)*sizeof(int));
	if (!v) nrerror("Memory allocation failure in ivector()");
	return v - nl;
}

void free_ivector(int *v, int nl)

{
	free((char*) (v + nl));
}

int **imatrix(int nrl, int nrh, int ncl, int nch)

{
	int i;
	int **m;

	m = (int **)malloc((unsigned)(nrh - nrl +1)*sizeof(int*));
	if (!m) nrerror("allocation failure 1 in imatrix()");
	m -= nrl;
	for(i = nrl; i <= nrh ; i++) {
		m[i]= (int*)malloc((unsigned)(nch - ncl +1)*sizeof(int));
		if (!m[i]) nrerror("allocation failure 2 in imatrix()");
		m[i] -= ncl;
	}
	return m;
}

void free_imatrix(int **m, int nrl, int nrh, int ncl)

{
	int i;

	for (i = nrh; i >= nrl; i--)
		free((char*) (m[i] + ncl));
	free((char *)(m + nrl));
}

void nrerror(char *error_text)

{
	fprintf(stderr,"\n%s\n",error_text);
	fprintf(stderr,"... now exiting to system ...\n");
	exit(1);
}