/*        Tabu search algorithm for constructing RPBIBD(2)s       */
/*           This algortithm was described in the paper           */
/*       Constructing some PBIBD(2)s by Tabu Search Algortihm     */
/*                        Luis B. Morales  2000                   */
/*                                                                */
/* In order to use this program you will need a Unix-based sytem  */

#include <stdlib.h>
#include <string.h>
#include <stdio.h>
#include <time.h>

#define funcion(a) ((a)==(l2) ? (1) : (a-l1))

#define mini(a,b) (a < b ? a : b)
#define Maxv 51 
#define Maxr 21
#define Maxg 21 

int v,b,r,kv,g,n1,n2,l1,l2;
int fmin,ntabu,int1,int2;
int Xt[Maxr][Maxv];

int costot,mejor,nits,borra1,itera;
int aux[2][Maxv], Ct[Maxv][Maxv], Int[Maxv][Maxv];

int pn=0, Np=1;

struct TABU
{
  int p,m,i,j,k;
} tabu[100],el_li;


int cuenta() {
  int i=2, j=0, a, cont=0;

  for (i=1; i<v; i++)
    for (j=i+1; j<=v; j++) {
      if (Ct[i][j] == l1) {
	if (Int[i][j] != int1) {
	  a = Int[i][j] - int1;
	  cont += abs(a);
	}
      }
      else {
	if (Ct[i][j] == l2) {
	  if (Int[i][j] != int2) {
	    a = Int[i][j] - int2;
	    cont += abs(a);
	  }
	}
	else {
	  a = Int[i][j] - int1;
	  cont += abs(a);
	}
      }
    }
  return (cont);
}

int cuenta1() {
  int i=1, j=0, cont=0, a, b;

  for (i=1; i<v; i++)
    for (j=i+1; j<=v; j++) {
      if (i == aux[0][0]) {
	if (aux[0][j] != -1)
	  a = aux[0][j];
	else a = Ct[i][j];
      }
      else
	if (i == aux[1][0]) {
	  if (aux[1][j] != -1)
	    a = aux[1][j];
	  else a = Ct[i][j];
	}
	else
	  if (j == aux[0][0]) {
	    if (aux[0][i] != -1)
	      a = aux[0][i];
	    else a = Ct[i][j];
	  }
	  else
	    if (j == aux[1][0]) {
	      if (aux[1][i] != -1)
		a = aux[1][i];
	      else a = Ct[i][j];
	    }
	    else a = Ct[i][j];

      if (a == l1) {
	if (Int[i][j] != int1) {
	  b = Int[i][j] - int1;
	  cont += abs(b);
	}
      }
      else {
	if (a == l2) {
	  if (Int[i][j] != int2) {
	    b = Int[i][j] - int2;
	    cont += abs(b);
	  }
	}
	else {
	  b = Int[i][j] - int1;
	  cont += abs(b);
	}
      }
    }
  return (cont);
}

int Costo()
{
  int i, j, costo=0, p=0;

  for (i=1; i<v; i++)
    for (j=i+1; j<=v; j++) {
      p = funcion(Ct[i][j]);
      costo += p*p;
    }
  borra1 = cuenta();
  costo += borra1;
  return costo;
}

void cambiomatriz()
{
  int i,j,k,a,b;

  for(i=1; i<v; i++) {
    Ct[i][i] = -1;
    for(j=i+1; j<=v; j++)
      Ct[i][j] = 0;
  }
  Ct[v][v] =-1;

  for (k=1; k<=r; k++)
    for (i=1; i<v; i++)
      for (j=i+1; j<=v; j++)
	if (Xt[k][i] == Xt[k][j])
	 { Ct[i][j]++;  Ct[j][i]++;}

    for (i=1; i<v; i++)
      for (j=i+1; j<=v; j++)
	for (k=1; k<=v; k++)
	  if ((k!=i) && (k!=j)) {
	    if (k > i) a = Ct[i][k];
	    else a = Ct[k][i];
	    if (k > j) b = Ct[j][k];
	    else b = Ct[k][j];
	    if ((a == l1) && (a == b)) {
	      Int[i][j]++;
	      Int[j][i]++;
	    }
	}
}

/* Despliega conjunto y solucion actual */
void despsol()
{
  int i,j,l,ch=0,cont = 0;

  fprintf (stdout,"\n");
  fprintf (stdout,"     TS found an optima solution\n");
  fprintf (stdout,"     Blocks are \n");            
  for (l=1; l<=r; l++) {
    for (i=1; i<=g; i++){
      ch++;
      for (j=1; j<=v; j++)
	if (Xt[l][j] == i) {
	  fprintf (stdout,"%3d",j);
	}
	fprintf (stdout,";");
    }
  fprintf (stdout,"\n");
  }
  fprintf (stdout,"\n");
  fprintf (stdout,"     Enter for Concurrence Matrix \n");   
  getchar ();
  getchar ();
   for (i=1; i<=v; i++)
     for (j=i+1; j<=v; j++)
       Ct[j][i] = Ct[i][j];

   for (i=1; i<=v; i++) {
     for (j=1; j<=v; j++)
       fprintf (stdout,"%3d",Ct[i][j]);
     fprintf (stdout,"\n");
   }
  fprintf (stdout,"\n");
  fprintf (stdout,"     Enter for Intersection Matrix \n");   
  getchar ();
    for (i=1; i<=v; i++) {
      for (j=1; j<=v; j++)
	fprintf (stdout,"%3d",Int[i][j]);
     fprintf (stdout,"\n");
   }
}

void init()
{
  int i, j, j0, grupo0;

/*  INICIALIZA LA LISTA TABU                */
  memset(tabu,0,sizeof(tabu));
  memset(Ct,0,Maxv*Maxv*sizeof(int));
  memset(Xt,0,Maxr*Maxv*sizeof(int));
  memset(Int,0,Maxv*Maxv*sizeof(int));

/*  Encuentra la solucion inicial            */
  for (i=1; i<=r; i++)
    for (j=1; j<=v; j++)
      Xt[i][j] = j;

  for (i=1; i<=r; i++)
    for (j=v; j>=1; j--) {
      j0 = (int)lrand48() % j;
      j0++;
      grupo0 = Xt[i][j0];
      Xt[i][j0] = Xt[i][j];
      Xt[i][j] = grupo0;
    }

  for (i=1; i<=r; i++)
    for (j=1; j<=v; j++) {
      Xt[i][j] %= g;
      Xt[i][j]++;
    }

}

int estabu(i,j,k)
int i, j, k;
{
  int l;

  for (l=0; l<ntabu; l++)
    if ((tabu[l].i == i) && 
       ((tabu[l].j==j||tabu[l].k==k)||(tabu[l].j==k||tabu[l].k==j)))
      return 1;
    return 0;
}

void Tabu(solucion)
int *solucion;
{
  int ii,h,ni;
  int p,pena,iter = 0,minimo;
  int costop,sump,minimop;
  int i,j,i0,i1,k,k0,q,a,b;
  int Mfre[Maxr][Maxg][Maxg];
  int iMinimo, MinimoT;
  int indicer[Maxr],indicev[Maxv];
  int cambio,iii,kk,jj;

  cambiomatriz();
  mejor = costot = Costo();

  costop = costot - borra1;
  nits = 0;

  for (i=1; i<=r; i++)
    for (j=1; j<=20; j++)
      for (k=1; k<=20; k++)
	 Mfre[i][j][k] = 0;
  memset (Mfre,0,Maxr*Maxg*Maxg*sizeof(int));

  for (i=1; i<=r; i++)
   indicer[i] = i;
  for (i=1; i<=v; i++)
   indicev[i] = i;


/*  Inicio del proceso iterativo       */
  while (mejor > fmin && (nits - iter) < itera) {
     minimop = minimo = 30000;
     nits++;

    for (i=r; i>1; i--){
      j = ((int)lrand48() % r) + 1;
      cambio = indicer[j];
      indicer[j] = indicer[i];
      indicer[i] = cambio;
    }
    for (i=v; i>1; i--){
      j = ((int)lrand48() % v) + 1;
      cambio = indicev[j];
      indicev[j] = indicev[i];
      indicev[i] = cambio;
    }

/*      Busqueda del mejor movimiento en la iteracion repetir
	Calcula todos LOS MOVIMIENTOS        */

     for (kk=1; kk<=r; kk++){
       k0 = indicer[kk];
       for (iii=1; iii<v; iii++)
         for (jj=iii+1; jj<=v; jj++){
           i0 = indicev[iii]; i1 = indicev[jj];
           if (Xt[k0][i0] != Xt[k0][i1]) {
	    {
	     sump = 0;
	     memset(aux,-1,2*Maxv*sizeof(int));
	     aux[0][0] = i0;
	     aux[1][0] = i1;
	       for (i=1; i<=v; i++)
		 if((i!=i0) && (i!=i1)) {
		  if (Xt[k0][i0] == Xt[k0][i]) {
		   if (i > i0){
		      p = funcion(Ct[i0][i]);
		      aux[0][i] = Ct[i0][i]-1;
		   }
		   else{
		      p = funcion(Ct[i][i0]);
		      aux[0][i] = Ct[i][i0]-1;
		   }
		   sump -= p*p;
		   if (i > i0){
		      p = funcion(Ct[i0][i]-1);
		   }
		   else{
		      p = funcion(Ct[i][i0]-1);
		  }
		   sump += p*p;
		   if (i > i1){
		      p = funcion(Ct[i1][i]);
		      aux[1][i] = Ct[i1][i]+1;
		   }
		   else{
		      p = funcion(Ct[i][i1]);
		      aux[1][i] = Ct[i][i1]+1;
		  }
		   sump -= p*p;
		   if (i > i1){
		      p = funcion(Ct[i1][i]+1);
		   }
		   else{
		      p = funcion(Ct[i][i1]+1);
		   }
		   sump += p*p;
		 }
		 if (Xt[k0][i1] == Xt[k0][i]) {
		   if (i > i1){
		      p = funcion(Ct[i1][i]);
		      aux[1][i] = Ct[i1][i]-1;

		   }
		   else{
		      p = funcion(Ct[i][i1]);
		      aux[1][i] = Ct[i][i1]-1;

		  }
		   sump -= p*p;
		   if (i > i1){
		      p = funcion(Ct[i1][i]-1);
		   }
		   else{
		      p = funcion(Ct[i][i1]-1);

		   }
		   sump += p*p;
		   if (i > i0){
		      p = funcion(Ct[i0][i]);
		      aux[0][i] = Ct[i0][i]+1;
		   }
		   else{
		      p = funcion(Ct[i][i0]);
		      aux[0][i] = Ct[i][i0]+1;
		  }
		   sump -= p*p;
		   if (i > i0){
		      p = funcion(Ct[i0][i]+1);
		   }
		   else{
		      p = funcion(Ct[i][i0]+1);
		   }
		   sump += p*p;
		   }
		  }

	      /* update Int[a][b] */

	      for (ii=1; ii<=v; ii++)
	       for (q=ii+1; q<=v; q++)
		for (h=0; h<=1; h++){
		  if (h==0 ) {a = 0; b = i0;}
		  else     {a = 1; b = i1;}
		  if (ii!=b && q!=b){
		   if ((aux[a][ii]==l1)&&(Ct[ii][b]==l1)){
		      if ((aux[a][q]!=l1)&&(aux[a][q]!=-1)&&(Ct[b][q]==l1)){
			 Int[ii][q]--;  }
		      if ((aux[a][q]==l1)&&(Ct[b][q]!=l1)){
			 Int[ii][q]++;   }

		   }
		   if ((aux[a][ii]!=l1)&&(aux[a][ii]!=-1)&&(Ct[b][ii]==l1)){
		      if ((aux[a][q]==l1)&&(Ct[q][b]==l1)){
			 Int[ii][q]--;  }
		      if ((aux[a][q]!=l1)&&(aux[a][q]!=-1)&&(Ct[b][q]==l1)){
			 Int[ii][q]--;  }
		      if ((aux[a][q]==-1)&&(Ct[b][q]==l1)){
			 Int[ii][q]--;  }
		   }
		   if ((aux[a][ii]==l1)&&(Ct[b][ii]!=l1)){
		      if ((aux[a][q]==l1)&&(Ct[b][q]==l1)){
			 Int[ii][q]++;  }
		      if ((aux[a][q]==l1)&&(Ct[b][q]!=l1)){
			 Int[ii][q]++;  }
		      if ((aux[a][q]==-1)&&(Ct[b][q]==l1)){
			 Int[ii][q]++;  }
		   }
		   if ((aux[a][ii]==-1)&&(Ct[b][ii]==l1)){
		      if ((aux[a][q]==l1)&&(Ct[b][q]!=l1)){
			 Int[ii][q]++;  }
		      if ((aux[a][q]!=l1)&&(aux[a][q]!=-1)&&(Ct[b][q]==l1)){
			 Int[ii][q]--;  }
		   }
		  }
		}
	     for (ii=1; ii<=v; ii++){
		if (ii != i0 && ii!=i1)
		  for (q=1; q<=v; q++)
		     if ((q != i0) && (q !=i1) && (q!=ii)){
		       if (Ct[ii][q] ==l1)
			 for (h=0; h<=1; h++){
			  if (h==0 ) {a = 0; b = i0;}
			  else     {a = 1; b = i1;}
			  if (aux[a][q]==l1 && Ct[b][q]!=l1){
			      if (b<ii) Int[b][ii]++;
			      else Int[ii][b]++;
			      }
			  if ((aux[a][q]!=l1)&&(aux[a][q]!=-1)&&(Ct[b][q])==l1){
			      if (b<ii) Int[b][ii]--;
			      else Int[ii][b]--;
			      }
			 }
		       }
	     }
	     for (q=1; q<=v; q++)
		 if (aux[0][q]!=-1 && aux[1][q]!=-1){
		      if (aux[0][q]==l1 && aux[1][q]==l1)
			 if (i0<i1) Int[i0][i1]++;
			 else Int[i1][i0]++;
		      if (Ct[i0][q]==l1 && Ct[i1][q]==l1)
			 if (i0<i1) Int[i0][i1]--;
			 else Int[i1][i0]--;
		 }

	     a = cuenta1();
	     costot = costop + sump + a;

	     if (costot <= minimo) {
	      if (costot >= mejor){
		 pena = costot  + Mfre[k0][Xt[k0][i0]][Xt[k0][i1]];
		 if (!estabu(k0,i0,i1))
		   if (pena <= minimop) {
		     minimo = costot;
		     minimop = pena;
		     el_li.p = a;
		     el_li.m = pena;
		     el_li.i = k0;
		     el_li.j = i0;
		     el_li.k = i1;
		   }
	       }
	       else {
		 minimop = minimo = costot;
		 el_li.p = a;
		 el_li.m = costot;
		 el_li.i = k0;
		 el_li.j = i0;
		 el_li.k = i1;
	       }
	     }
	     for (i=1; i<v; i++)
	       for (q=i+1; q<=v; q++)
		 Int[i][q] = Int[q][i];
	   }
    }
   }
   }
/*      Fin del proceso de busqueda del mejor movimiento en una iteracion  */

     if (el_li.m >= mejor)
	 el_li.m = el_li.m  -
         Mfre[el_li.i]
             [Xt[el_li.i][el_li.j]][Xt[el_li.i][el_li.k]];

/*   ACTUALIZA LA MATRIZ DE FRECUENCIA           */
     Mfre[el_li.i]
         [Xt[el_li.i][el_li.j]][Xt[el_li.i][el_li.k]]++; 
      Mfre[el_li.i]
          [Xt[el_li.i][el_li.k]][Xt[el_li.i][el_li.j]]++; 

     costop = el_li.m;

     p = Xt[el_li.i][el_li.j];
     Xt[el_li.i][el_li.j] = Xt[el_li.i][el_li.k];
     Xt[el_li.i][el_li.k] = p;

     for (i=1; i<=v; i++)
       if((i!=el_li.j) && (i!=el_li.k)) {
	 if (Xt[el_li.i][el_li.j] == Xt[el_li.i][i]) {
	   if (i > el_li.j) Ct[el_li.j][i]++;
	   else Ct[i][el_li.j]++;
	   if (i > el_li.k) Ct[el_li.k][i]--;
	   else Ct[i][el_li.k]--;
	 }
	 if (Xt[el_li.i][el_li.k] == Xt[el_li.i][i]) {
	   if (i > el_li.k) Ct[el_li.k][i]++;
	   else Ct[i][el_li.k]++;
	   if (i > el_li.j) Ct[el_li.j][i]--;
	   else Ct[i][el_li.j]--;
	 }
       }

     for (i=1; i<v; i++)
      for(j=i+1; j<=v; j++)
	    Ct[j][i] = Ct[i][j];

     memset(Int,0,Maxv*Maxv*sizeof(int));

     for (i=1; i<v; i++)
       for (j=i+1; j<=v; j++)
	 for (k=1; k<=v; k++)
	   if ((k!=i) && (k!=j)) {
	     if (k > i) a = Ct[i][k];
	     else a = Ct[k][i];
	     if (k > j) b = Ct[j][k];
	     else b = Ct[k][j];
	     if ((a == l1) && (a == b)) {
	       Int[i][j]++;
	       Int[j][i]++;
	     }
	   }


/*   ACTUALIZA LA LISTA TABU           */
     tabu[nits%ntabu] = el_li;


/*   ACTUALIZA LA FUNCION ASPIRACION   */
     if (costop < mejor) {
       mejor = costop;
       iter = nits;
     }
/*
     fprintf(errorf,"%5d%5d%5d%5d ",nits,costop,el_li.p,mejor);
     fprintf(errorf,"%3d%3d%3d\n",el_li.i,el_li.j,el_li.k);
*/
     costop -= el_li.p;

  }
  if (mejor==fmin) {
     despsol(); 
     *solucion = 1;
  }
  fprintf (stdout,"The best value found is %d \n",mejor);
} 


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

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

  fprintf(stdout,"    Tabu search for constructing RPBIBD(2)s \n");
  fprintf(stdout,
          "    with parameters v,b,r,k,lambda1,lambda2,n1,n2, p11^1,p11^2 \n");

  fprintf(stdout,"    Number of TS runs to be performed ");
  fscanf (stdin, "%d",&nstatis);
  fprintf(stdout,"    Enter parameter v(<50)  ");
  fscanf (stdin, "%d",&v);
  fprintf(stdout,"    Enter parameter b  ");
  fscanf (stdin, "%d",&b);
  fprintf(stdout,"    Enter parameter r(<30)  ");
  fscanf (stdin, "%d",&r);
  fprintf(stdout,"    Enter parameter k(<10)  ");
  fscanf (stdin, "%d",&kv);
  if (v>50 || r>30 || kv>10) {
     fprintf(stdout," Error V>50, r>30 or k>10  \n");
     exit (10);
  }
  if (b%r != 0) {
      fprintf (stdout,"    The design is not resolvable  \n");
      exit(10);      
  }
  g = b/r;
  fprintf(stdout,"    Enter parameter lambda1 ");
  fscanf (stdin, "%d",&l1);
  fprintf(stdout,"    Enter parameter lambda2 ");
  fscanf (stdin, "%d",&l2);
  fprintf(stdout,"    Enter parameter n1 ");
  fscanf (stdin, "%d",&n1);
  fprintf(stdout,"    Enter parameter n2 ");
  fscanf (stdin, "%d",&n2);
  fprintf(stdout,"    Enter parameter p11^1 ");
  fscanf (stdin, "%d",&int1);
  fprintf(stdout,"    Enter parameter p11^2 ");
  fscanf (stdin, "%d",&int2);
  fprintf(stdout,"    Enter maximum of iterations ");
  fscanf (stdin, "%d",&itera);
  fprintf(stdout,"    Enter size of tabu list ");
  fscanf (stdin, "%d",&ntabu);
  fprintf(stdout,"\n");
  fmin = n2*v/2;

  if ((v*r!=b*kv) || ((v-1)!=n1+n2) || (r*(kv-1)!=l1*n1+l2*n2) ){
    fprintf(stdout,"Error in the parameters of first kind\n");
    exit (20);
  }
  fprintf (stdout, "\n");
  fprintf (stdout, "(%5d%5d%5d%5d%5d%5d%5d%5d%5d%5d)-RPBIBD(2)\n", 
                    v,b,r,kv,l1,l2,n1,n2,int1,int2);

  solucion = 0;
  for (j=1;j<=nstatis;j++){
      init();
      Tabu(&solucion);
      if (solucion==1) break;
  }
  if (solucion==0){ 
     fprintf (stdout,"     TS did not produce an optima solution\n");
     exit (30);
  }
 exit (40);
}