/*  Tabu search algorithm for constructing difference families   */
/*           This algortithm was described in the paper          */
/*  Constructing Difference Faliies through an Optimization      */
/*                     Approach. Six new BIBD's                  */
/*                        Luis B. Morales  1998                   */
/*                                                                */
/* In order to use this program you will need a Unix-based sytem  */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#include <sys/time.h>
#include <sys/resource.h>       
#define RANDMAX 0x7fffffff
#define maxm 160 
#define maxn 10
#define maxt 20
#define maxk 40 
#define max2n 80
#define maxsets 1500 

int kn;
int ss[maxsets][maxk];
int o;
int v,b,r,k,lam;
int m,n,t;
int n_tabu,max_iter,iter;

struct move {
   int i,j,l,jp;
};

void write_sol(x,c)
char x[maxt][maxm][maxn];
int c[maxn][maxn][maxm];
{ int i,j,l,l1;

   fprintf (stdout,"\n");
   fprintf (stdout,"     TS found an optima solution\n");
   fprintf (stdout,"     Base blocks are \n");
   for (i=1;i<=t; i++){
     for (l=1; l<=n; l++)
      for (j=0; j<m; j++)
	if (x[i][j][l] == (char)1)
	   if (n>1) fprintf (stdout,"%3d_%d",j,l);
	   else fprintf (stdout,"%3d",j);
	fprintf(stdout,"\n");
   }

   fprintf (stdout,"\n");
   fprintf (stdout,"     (n x n x m)-array B\n");

  for (l=1;l<=n;l++)
    for (l1=1;l1<=n;l1++){
      fprintf (stdout,"%3d %3d  ",l,l1);
      for (j=0;j<m;j++)
	fprintf (stdout,"%3d",c[l][l1][j]);

      fprintf (stdout,"\n");
    }
}

void set_s(y)
int *y;
{
 int j,p;
 int contar = 0;
 int a[maxk],b,m_1l;

 memset(a,0,maxk*sizeof(int));
 m_1l = (m-1)*lam;
 p = 1;
 for (j=1;j<=kn;j++){
   if (y[j]==1) p++;
   else{
      if (o<maxsets){
         a[p]++;
         b = a[p]*(a[p]-1);
         if (b<=m_1l){
	    contar++;
	    ss[o][contar] = p;
         }
      }
   }
 }
 if (contar==k) o++;
}

void sub(sum,h,y)
int sum,h;
int *y;
{
 int i,s;

 for (i=0;i<=sum;i++)
   if (i<=1){
     y[h] = i;
     s = sum - i;
     if (h==(kn-1)){
       if (s<=1){
	   y[kn] = s;
	   set_s(y);
       }
     }
     else{
	sub(s,(h+1),y);
    }
 }
}

int back_tracking(a,b,node,n_s)
int a[max2n][maxsets],b[max2n],node[maxt];
int n_s;
{
 int i,l,h;
 int true,i1;
 int n_node;
 int suma[max2n];


 n_node = 1;
 i = 1;
 while (i<=n_s){
   node[n_node] = i;
   for (h=1;h<=2*n;h++)
	  suma[h] = a[h][node[n_node]];
   i1 = i;
   while ((n_node<t)&&(n_node>0)){
     while ((i1<=n_s)){
       true = 1;
       while ((n_node<=t)&&(true !=0)){
	   n_node++;
	   node[n_node] = i1;

	   for (h=1;h<=2*n;h++)
	      suma[h] += a[h][node[n_node]];
	   l=1;
	   while (l<=2*n && true !=0){
	      if (suma[l]>b[l]){
		 true = 0;
		 for (h=1;h<=2*n;h++)
		   suma[h] -= a[h][node[n_node]];
		 n_node--;
	      }
	      l++;
	   }
       }
       i1++;
     }
     i1 = node[n_node]+1;
     if (n_node == t){
	   true = 1;
	   l=1;
	   while (l<=2*n && true !=0){
	      if (suma[l] != b[l])
		 true = 0;
	      l++;
	   }
	  if (true == 1){
	    return (1);
	  }
     }
     for (h=1;h<=2*n;h++)
	 suma[h] -= a[h][node[n_node]];
     n_node--;
   }
   i++;
   n_node = 1;
 }
 return (0);
}

int sol_inicial(node)
int node[maxt];
{
 int i,j,jp,s;
 int w,n_s;
 int a[max2n][maxsets],b[max2n],y[max2n];
 int cambio[maxk];

  kn = k+n-1;
  w = n-1;
  o = 1;

  for (i=0;i<=1;i++){
	y[1]=i;
	s = w-i;
	sub(s,2,&y);
  }
  o--;

  n_s = o/2;

    for (i=o;i>0;i--){
      jp = ((int)lrand48()%o) + 1;
      for (j=1;j<=k;j++){
        cambio[j] = ss[jp][j];
        ss[jp][j] = ss[i][j];
        ss[i][j] = cambio[j];
      }
    }

/*
  for (i=1;i<=n_s;i++){
     for (j=1;j<=k;j++){
        cambio[j] = ss[i*2][j];
        ss[i*2][j] = ss[o-i+1][j];
        ss[o-i+1][j] = cambio[j];
     }
  }
*/

  memset(a,0,maxn*maxsets*sizeof(int));

  for (i=1;i<=o;i++){
     for (j=1;j<=k;j++){
        a[ss[i][j]][i]++;
     }
  }
  for (j=1;j<=o;j++){
    for (i=1;i<=n;i++){
	 a[i+n][j] =  a[i][j]*(a[i][j]-1);
    }
  }

  for (j=1;j<=o;j++)
	 a[(2*n+1)][j] = 1;

  for (i=1;i<=n;i++){
      b[n+i] = (m-1)*lam;
      b[i] = r;
  }

  b[2*n+1] = t;


  n_s = o;
/*
  for (i=1;i<=n_s;i++){
     for (j=1;j<=k;j++)
        fprintf (stdout,"%3d",ss[i][j]);
     fprintf (stdout,"\n");
  }

  fprintf (stdout,"A= \n");

  for (i=1;i<=(2*n+1);i++){
     for (j=1;j<=n_s;j++)
        fprintf (stdout,"%3d",a[i][j]);
     fprintf (stdout,"\n");
  }

*/
  if (n_s==1)
      for (i=1;i<=t;i++)
          node[i] = 1;
  if (n_s>1)
     if (back_tracking(a,b,node,n_s)==0){
        fprintf(stdout,"There is not subscript ocurrence matrix \n");
        return (0);
   }
/*
  for (i=1;i<=t;i++)
      fprintf (stdout,"node[%d] = %3d\n",i,node[i]);
  fprintf (stdout," \n");
*/

  return (1);
}

int ns(x,c)
char x[maxt][maxm][maxn];
int c[maxn][maxn][maxm];
{
  int i,j,j1,l,l1;
  int h;
  int suma,dif;

  memset(c,0,maxn*maxn*maxm*sizeof(int));

  for (i=1;i<=t;i++){
     for (j=0;j<m;j++)
	for (j1=0;j1<m;j1++)
	  for (l=1;l<=n;l++)
	    for (l1=1;l1<=n;l1++)

	    if ((x[i][j][l]==(char)1)&&(x[i][j1][l1]==(char)1)){
	       h = j-j1;
	       if (h<0)
		  h = m + h;
		c[l][l1][h]++;
		/*fprintf (stdout,"%3d%3d c[%d][%d][%d]=%3d\n"
				,j,j1,l,l1,h,c[l][l1][h]); */

	    }
  }
  suma = 0;
  for (l=1;l<=n;l++){
    c[l][l][0] = lam;
    for (l1=1;l1<=n;l1++){
      /*fprintf (stdout,"%3d %3d  ",l,l1);*/
      for (j=0;j<m;j++){
	 dif = c[l][l1][j] - lam;
	 dif = dif * dif;
	 suma += dif;
	/*fprintf (stdout,"%3d",c[l][l1][j]);*/
      }
      /*fprintf (stdout,"\n");*/
    }
  }
  return(suma);
}

int estatus_tabu(mov_tabu,i,j,l,jp,costo)
   struct move mov_tabu[100];
   int i,j,l,jp;
   int costo;
{
 int ii;
 for (ii=0;ii<n_tabu;ii++)
   if ((mov_tabu[ii].i==i && mov_tabu[ii].l==l) &&
       (mov_tabu[ii].jp==j))
       return(0);

  return(1);
}

int optimizador(y,solucion)
int y[max2n];
int *solucion;
{
 int i,j,l,i1,j1,l1;
 int jp,jp1;
 int i_min,j_min,l_min,jp_min;
 int costo,cparcial,fparcial;
 int minimo,b_iter;
 int best;
 char x[maxt][maxm][maxn];
 char best_x[maxt][maxm][maxn];
 int c[maxn][maxn][maxm];
 int ca[maxn][maxn][maxm];
 int d_t[maxt],d_m[maxm],d_n[maxm];
 int cambio;
 int na,ma,h1,h2,h3,h4,dif;
 int aspiracion,ind;
 struct move mov_tabu[100];

  memset(x,(char)0, maxt*maxm*maxn*sizeof(char));

  for (i=0;i<m;i++)
     d_m[i] = i;

  for (l=1;l<=t;l++){
    for (i=(m-1);i>=0;i--){
      j = (int)lrand48()%m;
      cambio = d_m[j];
      d_m[j] = d_m[i];
      d_m[i] = cambio;
    }
    for (j=1;j<=k;j++){
       /*fprintf (stdout,"%3d",ss[y[l]][j]);*/
       x[l][d_m[j-1]][ss[y[l]][j]] = (char)1;
    }
    /*fprintf (stdout,"\n");*/
  }


  costo = ns(x,c);
  best = 31000;
  aspiracion = 31000;

  for (i=0;i<n_tabu;i++){
    mov_tabu[i].i = -1;
    mov_tabu[i].j = -1;
    mov_tabu[i].l = -1;
    mov_tabu[i].jp = -1;
  }

  for (i=0;i<m;i++)
    d_m[i] = i;
  for (i=1;i<t;i++)
    d_t[i] = i;
  }
  for (i=1;i<t;i++)
    d_n[i] = i;

  iter = b_iter = 0;
  cparcial = costo;

  while ((iter-b_iter<max_iter)&&(best>0)){
    iter++;
    jp_min = minimo = 31000;

    for (i=t;i>0;i--){
      j = ((int)lrand48()%t) + 1;
      cambio = d_t[j];
      d_t[j] = d_t[i];
      d_t[i] = cambio;
    }
    for (i=(m-1);i>=0;i--){
      j = ((int)lrand48()%m);
      cambio = d_m[j];
      d_m[j] = d_m[i];
      d_m[i] = cambio;
    }
    for (i=n;i>0;i--){
      j = ((int)lrand48()%n) + 1;
      cambio = d_n[j];
      d_n[j] = d_n[i];
      d_n[i] = cambio;
    }

    for (i1=1;i1<=t;i1++)
      for (j1=0;j1<m;j1++)
	for (l1=1;l1<=n;l1++){

	   i = d_t[i1];
	   j = d_m[j1];
	   l = d_n[l1];

	   if (x[i][j][l]==(char)1){
	     for (jp1=0;jp1<m;jp1++){
	       jp = d_m[jp1];
	       if (x[i][jp][l]==(char)0){


		 for (na=1;na<=n;na++)
		   for (ma=0;ma<m;ma++)
		    if (ma!=j || na!=l)
		    if (x[i][ma][na]==(char)1){
		       h1 = j-ma;
		       if (h1<0) h1 = m+h1;
		       ca[l][na][h1] = c[l][na][h1];
		       h2 = jp - ma;
		       if (h2<0) h2 = m+h2;
		       ca[l][na][h2]  = c[l][na][h2];
		       h3 = ma-j;
		       if (h3<0) h3 = m+h3;
		       ca[na][l][h3] = c[na][l][h3];
		       h4 = ma - jp;
		       if (h4<0) h4 = m+h4;
		       ca[na][l][h4] = c[na][l][h4];
		    }


		 fparcial = 0;
		 for (na=1;na<=n;na++)
		   for (ma=0;ma<m;ma++)
		    if (ma!=j || na!=l)
		    if (x[i][ma][na]==(char)1){
		       h1 = j-ma;
		       h2 = jp - ma;
		       if (h1<0) h1 = m+h1;
		       dif = ca[l][na][h1] - lam;
		       fparcial -= dif*dif;
		       dif--;
		       ca[l][na][h1]--;
		       fparcial += dif*dif;

		       if (h2<0) h2 = m+h2;
		       dif = ca[l][na][h2] - lam;
		       fparcial -= dif*dif;
		       dif++;
		       ca[l][na][h2]++;
		       fparcial += dif*dif;

		       h3 = ma-j;
		       if (h3<0) h3 = m+h3;
		       dif = ca[na][l][h3] - lam;
		       fparcial -= dif*dif;
		       dif--;
		       ca[na][l][h3]--;
		       fparcial += dif*dif;

		       h4 = ma - jp;
		       if (h4<0) h4 = m+h4;
		       dif = ca[na][l][h4] - lam;
		       fparcial -= dif*dif;
		       dif++;
		       ca[na][l][h4]++;
		       fparcial += dif*dif;

		    }
		  costo = cparcial + fparcial;

		  if (costo >= aspiracion){
		    if (costo<minimo)
		     if (estatus_tabu(mov_tabu,i,j,l,jp,costo)){
			minimo = costo;
			i_min = i;
			j_min = j;
			jp_min = jp;
			l_min = l;
		     }
		  }
		  else{
		    if (costo<minimo){
			minimo = costo;
			i_min = i;
			j_min = j;
			jp_min = jp;
			l_min = l;
		     }
		  }
	       }
	     }
	   }
	}

	x[i_min][j_min][l_min]  = (char)0;
	x[i_min][jp_min][l_min] = (char)1;

	ind = iter % n_tabu;
	mov_tabu[ind].i = i_min;
	mov_tabu[ind].j = j_min;
	mov_tabu[ind].l = l_min;
	mov_tabu[ind].jp = jp_min;
     
	if (minimo<best){
	    best = minimo;
	    b_iter = iter;
	    for (i=1 ;i<=t;i++)
	      for (j=0;j<m;j++)
		 for (l=1;l<=n;l++)
		   best_x[i][j][l] = x[i][j][l];
	}
	aspiracion = best;
	cparcial = ns(x,c);
  }

  for (i=1;i<=t;i++)
     for (j=0;j<m;j++)
        for (l=1;l<=n;l++)
	   x[i][j][l] = best_x[i][j][l];

  costo = ns(x,c);

  if (costo==0) {
     write_sol(x,c);
     *solucion = 1;
  }
  fprintf (stdout,"The best value found is %d \n",best);
  return (1);
}


void main(argc, argv)
int argc;
char *argv[];
{
  int node[maxt];
  int j,nstatis,solucion;
  time_t t1;

  srand48((unsigned) time(&t1));

  fprintf(stdout,"    Tabu search for constructing difference family\n");
  fprintf(stdout,"           with parameters m,n,t,r,k,lambda\n");

  fprintf(stdout,"    Number of TS runs to be performed ");
  fscanf (stdin, "%d",&nstatis);
  fprintf(stdout,"    Enter parameter m  ");
  fscanf (stdin, "%d",&m);
  fprintf(stdout,"    Enter parameter 0<n<4  ");
  fscanf (stdin, "%d",&n);
  if (n>3) {
     fprintf(stdout," Error ");
     exit (10);
  }  
  fprintf(stdout,"    Enter parameter t  ");
  fscanf (stdin, "%d",&t);
  fprintf(stdout,"    Enter parameter r  ");
  fscanf (stdin, "%d",&r);
  fprintf(stdout,"    Enter parameter k  ");
  fscanf (stdin, "%d",&k);
  fprintf(stdout,"    Enter parameter lambda ");
  fscanf (stdin, "%d",&lam);
  fprintf(stdout,"    Enter maximum of iterations ");
  fscanf (stdin, "%d",&max_iter);
  fprintf(stdout,"    Enter size of tabu list ");
  fscanf (stdin, "%d",&n_tabu);
  fprintf(stdout,"\n");

  if ((k*t!=n*r) || (r*(k-1)!= (lam*(m*n-1)))){
    fprintf(stdout,"Error in the parameters \n");
    exit (20);
  }
  v = m*n;
  b = m*t;

  fprintf (stdout,"    Difference family with parameters\n");
  fprintf (stdout, "%5d%5d%5d%5d%5d%5d\n", m,n,t,r,k,lam);
  fprintf (stdout,"    (%5d%5d%5d%5d%5d)-BIBD\n",v,b,r,k,lam);

  solucion = 0;
  memset(node,0,maxt*sizeof(int));

  if (sol_inicial(node) != 0){ 
      for (j=1;j<=nstatis;j++){
          if (solucion==0){ 
             optimizador(node,&solucion);
          }
      }
  }
  else
    exit (50);
  fprintf (stdout,"     TS did not produce an optima solution\n");
  exit (30);
}