/* Tabu search algorithm for constructing  nonresolvable PBIBD(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>
#include <sys/time.h>
#define maxv 62
#define maxb 620

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


int V, B, R, K, fminimo=0, MaxIt=0, ntabu, ind_lis;
int int1, int2, n1,n2, l1, l2, borra1, num, n;
int aux[2][maxv],Int[maxv][maxv];

int num_iter, pn = 0, Np = 1;
int best,mejor_g;
char vb[maxb][maxv];


struct ELEM {
  int p,a,b,c,d;
  int m;
} Ltabu[25], el_li;


int cuenta(C)
char C[][maxv];
{
  int i, j, a, cont=0;

  for (i=1; i<V; i++)
    for (j=i+1; j<=V; j++) {
      if ((int)C[i][j] == l2)
	a = Int[i][j] - int2;
      else
	a = Int[i][j] - int1;
      cont += abs(a);
    }
  return cont;
}

int cuenta1(C)
char C[][maxv];
{
  int i, j, 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 = (int)C[i][j];
      }
      else
	if (i == aux[1][0]) {
	  if (aux[1][j] != -1)
	    a = aux[1][j];
	  else a = (int)C[i][j];
	}
	else
	  if (j == aux[0][0]) {
	    if (aux[0][i] != -1)
	      a = aux[0][i];
	    else a = (int)C[i][j];
	  }
	  else
	    if (j == aux[1][0]) {
	      if (aux[1][i] != -1)
		a = aux[1][i];
	      else a = (int)C[i][j];
	    }
	    else a = (int)C[i][j];

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


/* cambia la configuracion vb[][] */
void Cambia(co1,co2,ele1,ele2)
int co1,co2,ele1,ele2;
{
   vb[co1][ele1] = vb[co2][ele2] = (char)0;
   vb[co1][ele2] = vb[co2][ele1] = (char)1;
}

/* Despliega conjunto y solucion actual */
void despsol(C)
char C[][maxv];
{
  int i, j, l;

  fprintf (stdout,"\n");
  fprintf (stdout,"     TS found an optima solution\n");
  fprintf (stdout,"     The blocks are \n");
   for (i=1; i<=B; i++) {
     for (j=1; j<=V; j++)
       if (vb[i][j] ==(char)1) {
	 fprintf (stdout,"%3d",j);
       }
     if (i % 2 == 0) fprintf(stdout,";\n");
     else fprintf(stdout,";  ");
   }
  fprintf (stdout,"\n");
  fprintf (stdout,"     Enter for Concurrence Matrix \n");
  getchar ();
  getchar ();
   for (l=1; l<=V; l++) {
     for (i=1; i<=V; i++)
       fprintf(stdout,"%3d",(int)C[l][i]);
     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++)
        if (i!=j)
       	  fprintf(stdout,"%3d",Int[i][j]);
        else
          fprintf(stdout,"  *");
     fprintf(stdout,"\n");
   }
}

/*   genera nueva solucion inicial vb[][] */
void solution_initial()
{
  int i,j,l,ni,j0,cambio;
  int ra[maxv],sa[maxv];
  int ordered;

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

  for(j=1; j<=V; j++) {
    sa[j] = j;
    ra[j] = 0;
  }
  for (i=1; i<=B; i++) {
    ni = 0;
    for(j=V; j>0; j--) {
      j0 = ((int)lrand48() % V) + 1;
      cambio = sa[j0];
      sa[j0] = sa[j];
      sa[j] = cambio;
    }
    do {
      ordered = 1;
      for (j=1; j<V; j++)
	if (ra[sa[j]] > ra[sa[j+1]]) {
	  cambio = sa[j];
	  sa[j] = sa[j+1];
	  sa[j+1] = cambio;
	  ordered = 0;
	}
    } while(ordered != 1);
    l = 1;
    while(ni < K) {
      if (ra[sa[l]] < R) {
	vb[i][sa[l]] = (char)1;
	ra[sa[l]]++;
	ni++;
      }
      l++;
    }
  }
}

/* construye la matriz de incidencia */
void Crea_mat(C)
char C[][maxv];
{
  int l,i,j,k,a,b;

  for (l=0; l<=V; l++)
    for (i=0; i<=V; i++)
      C[l][i]=(char)0;

  for (l=1; l<=B; l++)
    for (i=1; i<V; i++)
      if (vb[l][i]== (char)1)
	for (j=i+1; j<=V; j++)
	  if (vb[l][j] == (char)1)
	    C[i][j] = ++C[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)) {
	    a = (int)C[i][k];
	    b = (int)C[j][k];
	    if ((a == l1) && (a == b))
	      Int[i][j] = ++Int[j][i];
	  }
}

/* evalua la funcion objetivo de acuerdo co la configuracion vb[][] */
int Eval_F_ob(C)
char C[][maxv];
{
  int i,j, dif ;
  int suma = 0;

   for (i=1; i<V; i++)
     for (j=i+1; j<=V; j++) {
       dif = funcion((int)C[i][j]);
       suma += dif * dif;
      }
  borra1 = cuenta(C);
  suma += borra1;
  return (suma);
}

/* pone un movimiento a la lista tabu */
void mete_ele(ele)
struct ELEM ele;
{
   ind_lis %= ntabu;
   Ltabu[ind_lis].a = ele.a;
   Ltabu[ind_lis].b = ele.b;
   Ltabu[ind_lis].c = ele.d;
   Ltabu[ind_lis].d = ele.c;
   ind_lis ++;
 }

/* busca un movimiento en la lista tabu */
int busca_ele(ele)
struct ELEM ele;
{
  int i;
  for (i=0; i<ntabu; i++)
    if ((((Ltabu[i].a==ele.a && Ltabu[i].b==ele.b) && 
	  (Ltabu[i].c==ele.c && Ltabu[i].d==ele.d))) ||
        (((Ltabu[i].a==ele.c && Ltabu[i].b==ele.d) && 
	  (Ltabu[i].c==ele.a && Ltabu[i].d==ele.b)))) 
      return 1;
  return 0;
}

/* procedimiento que contiene el metodo busqueda tabu */
void catabu(solucion)
int *solucion;
{
  int ii,i,j,l,m,ia,k,a,b,h,ni;
  int tabu, biter,iMinimo;
  int Minimo,MinimoT,Costo,Aspira;
  int aa;
  int Mfre[maxv][maxv];
  char Cij[maxv][maxv];
  struct ELEM el;
  int fparcial, cparcial;
  int dif1,df;
  int bandera,caux1,caux2;
  int indiceb[maxb],indicev[maxv];
  int i1,j1,ll,mm,cambio;


  ind_lis = 0;

  for (i=1; i<=B+1; i++)
    for (j=1; j<=V+1; j++)
      vb[i][j] = 0;

  for (i=1; i<=B; i++)
    for (j=1; j<=B; j++)
      Mfre[i][j] = 0;

  memset(Ltabu,0,sizeof(Ltabu));

  solution_initial();
  Crea_mat(Cij);
  cparcial = Eval_F_ob(Cij);
  cparcial -= borra1;
  biter = num_iter = 0;
  Aspira = 31000;
  best   = 31000;

  for (i=1; i<=B; i++)
   indiceb[i] = i;
  for (i=1; i<=V; i++)
   indicev[i] = i;


  while (((num_iter - biter) < MaxIt) && (best > fminimo)) {

  for(i=B; i>0; i--) {
      j = ((int)lrand48() % B) + 1;
      cambio = indiceb[j];
      indiceb[j] = indiceb[i];
      indiceb[i] = cambio;
    }
  for(i=V; i>0; i--) {
      j = ((int)lrand48() % V) + 1;
      cambio = indicev[j];
      indicev[j] = indicev[i];
      indicev[i] = cambio;
    }

    Minimo = 31000 ;

    for (i1=1; i1<B; i1++)
      for (ll=1; ll<=V; ll++){
        i = indiceb[i1]; l = indicev[ll];
        if (vb[i][l] == (char)1)
          for (j1=i1+1; j1<=B ; j1++)
            for(mm=1 ; mm<=V; mm++){
               j = indiceb[j1]; m = indicev[mm];
	      if ((vb[j][m]==(char)1)&&
                  (vb[i][m]==(char)0)&&(vb[j][l]==(char)0)) {
	       if ((int)lrand48() % 100 < 80) {
		 fparcial = 0;
		 memset(aux,-1,2*maxv*sizeof(int));
		 aux[0][0] = l;
		 aux[1][0] = m;

		 for (ia=1; ia<=V ;ia++) {
		   bandera = 0;
		   if ((ia !=l) && (ia !=m)) {
		     if (vb[i][ia]==(char)1) {
		       bandera = 1;
		       dif1 = (int)Cij[ia][l];
		       df = funcion(dif1);
		       aux[0][ia] = (int)Cij[ia][l] - 1;
		       fparcial -= df*df;
		       dif1--;
		       df = funcion(dif1);
		       fparcial += df*df;
		       dif1 = (int)Cij[ia][m];
		       df = funcion(dif1);
		       aux[1][ia] = (int)Cij[ia][m] + 1;
		       fparcial -= df*df;
		       dif1++;
		       df = funcion(dif1);
		       fparcial += df*df;
		     }
		     if (vb[j][ia]==(char)1) {
		       if (bandera == 1) {
			 caux1 = (int)Cij[ia][m]+1;
			 caux2 = (int)Cij[ia][l]-1;
		       }
		       else {
			 caux1 = (int)Cij[ia][m];
			 caux2 = (int)Cij[ia][l];
		       }
		       dif1 = caux1;
		       df = funcion(dif1);
		       aux[1][ia] = caux1 - 1;
		       fparcial -= df*df;
		       dif1--;
		       df = funcion(dif1);
		       fparcial += df*df;
		       dif1 = caux2;
		       df = funcion(dif1);
		       aux[0][ia] = caux2 + 1;
		       fparcial -= df*df;
		       dif1++;
		       df = funcion(dif1);
		       fparcial += df*df;
		     }
		   }
		 }

/* update Int[a][b] */
	      for (ii=1; ii<=V; ii++)
	       for (k=ii+1; k<=V; k++)
		for (h=0; h<=1; h++){
		  if (h==0 ) {a = 0; b = l;}
		  else     {a = 1; b = m;}
		  if (ii!=b && k!=b){

		   if ((aux[a][ii]==l1)&&(Cij[ii][b]==l1)){
		      if ((aux[a][k]!=l1)&&(aux[a][k]!=-1)&&(Cij[b][k]==l1)){
			 Int[ii][k]--;  }
		      if ((aux[a][k]==l1)&&(Cij[b][k]!=l1)){
			 Int[ii][k]++;   }
		   }
		   if ((aux[a][ii]!=l1)&&(aux[a][ii]!=-1)&&(Cij[b][ii]==l1)){
		      if ((aux[a][k]==l1)&&(Cij[k][b]==l1)){
			 Int[ii][k]--;  }
		      if ((aux[a][k]!=l1)&&(aux[a][k]!=-1)&&(Cij[b][k]==l1)){
			 Int[ii][k]--;  }
		      if ((aux[a][k]==-1)&&(Cij[b][k]==l1)){
			 Int[ii][k]--;  }
		   }
		   if ((aux[a][ii]==l1)&&(Cij[b][ii]!=l1)){
		      if ((aux[a][k]==l1)&&(Cij[b][k]==l1)){
			 Int[ii][k]++;  }
		      if ((aux[a][k]==l1)&&(Cij[b][k]!=l1)){
			 Int[ii][k]++;  }
		      if ((aux[a][k]==-1)&&(Cij[b][k]==l1)){
			 Int[ii][k]++;  }
		   }
		   if ((aux[a][ii]==-1)&&(Cij[b][ii]==l1)){
		      if ((aux[a][k]==l1)&&(Cij[b][k]!=l1)){
			 Int[ii][k]++;  }
		      if ((aux[a][k]!=l1)&&(aux[a][k]!=-1)&&(Cij[b][k]==l1)){
			 Int[ii][k]--;  }
		   }
		  }
		}
	     for (ii=1; ii<=V; ii++){
		if (ii != l && ii!=m)
		  for (k=1; k<=V; k++)
		     if ((k != l) && (k !=m) && (k!=ii)){
		       if (Cij[ii][k] ==l1)
			 for (h=0; h<=1; h++){
			  if (h==0 ) {a = 0; b = l;}
			  else     {a = 1; b = m;}
			  if (aux[a][k]==l1 && Cij[b][k]!=l1){
			      if (b<ii) Int[b][ii]++;
			      else Int[ii][b]++;
			      }
			  if ((aux[a][k]!=l1)&&(aux[a][k]!=-1)&&(Cij[b][k])==l1){
			      if (b<ii) Int[b][ii]--;
			      else Int[ii][b]--;
			      }
			 }
		       }
	     }
		for (k=1; k<=V; k++)
		 if (aux[0][k]!=-1 && aux[1][k]!=-1){
		      if (aux[0][k]==l1 && aux[1][k]==l1)
			 if (l<m) Int[l][m]++;
			 else Int[m][l]++;
		      if (Cij[l][k]==l1 && Cij[m][k]==l1)
			 if (l<m) Int[l][m]--;
			 else Int[m][l]--;
		 }
		 aa = cuenta1(Cij);
		 Costo = cparcial + fparcial + aa;

		 if (Costo >= Aspira)
		    Costo +=  Mfre[i][j];

		 if (Costo <= Minimo) {
		   tabu = 0;
		   if (Costo >= Aspira) {
		     el.a = i;
		     el.b = j;
		     el.c = l;
		     el.d = m;
		     tabu = busca_ele(el);
		   }
		   if ((Costo < Minimo) && (tabu==0)) {
		     Minimo = Costo ;
		     el_li.p = aa;
		     el_li.m = Minimo;
		     el_li.a = i;
		     el_li.b = j;
		     el_li.c = l;
		     el_li.d = m;
		   }
		 }
		 for (ii=1; ii<V; ii++)
		  for (k=ii+1; k<=V; k++)
		     Int[ii][k] = Int[k][ii];
	       }
	      }
          }
          }

    if (el_li.m >= Aspira)
	el_li.m -= (int)Mfre[el_li.a][el_li.b];

    cparcial = el_li.m;

     Mfre[el_li.a][el_li.b] = Mfre[el_li.a][el_li.b] + 1 ;
     Mfre[el_li.b][el_li.a] = Mfre[el_li.b][el_li.a] + 1 ;

    for (ia=1; ia<=V ;ia++) {
      bandera = 0;
      if ((ia !=el_li.c) && (ia !=el_li.d)) {
	if (vb[el_li.a][ia]==(char)1) {
	  Cij[el_li.c][ia] = --Cij[ia][el_li.c];
	  Cij[el_li.d][ia] = ++Cij[ia][el_li.d];
	}
	if (vb[el_li.b][ia]==(char)1)
	  if (bandera != 1) {
	    Cij[el_li.c][ia] = ++Cij[ia][el_li.c];
	    Cij[el_li.d][ia] = --Cij[ia][el_li.d];
	  }
      }
    }

     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 = (int)Cij[i][k];
	     else a = (int)Cij[k][i];
	     if (k > j) b = (int)Cij[j][k];
	     else b = (int)Cij[k][j];
	     if ((a == l1) && (a == b))
	       Int[i][j] = ++Int[j][i];
	   }

    mete_ele(el_li);
    Cambia(el_li.a,el_li.b,el_li.c,el_li.d);

    if (cparcial < best) {
	 best = cparcial;
	 biter = num_iter;
    }
    cparcial -= el_li.p;


    Aspira = best;
    num_iter++;
  }

   if (best == fminimo){
     despsol(Cij);
     *solucion = 1;
   }
}

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 PBIBD(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",&K);
  if (V>50 || R>30 || K>10) {
     fprintf(stdout," Error v>50, r>30 or k>10  \n");
     exit (10);
  }
  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",&MaxIt);
  fprintf(stdout,"    Enter size of tabu list ");
  fscanf (stdin, "%d",&ntabu);
  fprintf(stdout,"\n");
  fminimo = n2*V/2;

  if ((V*R!=B*K) || ((V-1)!=n1+n2) || (R*(K-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,K,l1,l2,n1,n2,int1,int2);

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