/* See English below */
/*********************************************************************/
/* MBMUDs.c, Generador de MBMUDs a base de fixar el nombre d'uns per fila
   i columna i minimitzar les quadruples amb moviments de tipus 
   (10,01)->(01,10). Busca l'optim o, en el seu defecte, el millor minim
   local despres d'un numero d'iteracions donat. Mes detalls a l'article 
   mencionat a continuacio.
   
   Si us plau, aquest programa es codi lliure, per referenciar-lo citeu:
   Bofill, P. & Torras, C., "MBMUDs: a combinatorial extension of BIBDs
   showing good optimality behaviour", Journal of Statistical Planning
   and Inference 124 (2004) 185-204.

   Per compilar: cc -o MBMUDs -O MBMUDs.c -lm

   Pau, Maig 2002 (pau@ac.upc.es)    
   Ultima revisio, 29 de Juliol 2004 (la versio de l'1 de Juliol tenia
   un error: a cada iteracio no recorria exhaustivament totes les
   quadruples possibles). */
/*********************************************************************/
/* MBMUDs.c, MBMUD generator based on setting the right number of ones
   per row and column, and then minimizing the number of quadruples
   using unit rearrangements of the kind (10,01)->(01,10). Searches
   the optimal configuration or, in case of failure, the best local 
   minimum within a given number of iterations. See the reference below 
   for details.
   
   This program is free code. Please, use the following reference for
   citation:
   Bofill, P. & Torras, C., "MBMUDs: a combinatorial extension of BIBDs
   showing good optimality behaviour", Journal of Statistical Planning
   and Inference 124 (2004) 185-204.

   Compile with: cc -o MBMUDs -O MBMUDs.c -lm

   Pau, May 2002 (pau@ac.upc.es)  
   Last revision, July 29th 2004 (The version of July 1st had a bug:
   Not all possible quadruples were tested at every iteration). */
/*********************************************************************/

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

#define MAXB 100
#define MAXV 75
typedef struct {
  int queden; 
  int llista[MAXB]; 
}sacboles; /* A bag with balls */

typedef struct {
  int erre,v0; 
  int ka,b0;
  int lambda,f0;
  int Qstar;
}descrdis; /* Design descriptors */

long seed;

void SacNou(int nboles,sacboles *sac){
  /* Inicializes a bag with nboles balls */

  int i;
  for (i=0;i<nboles;i++) sac->llista[i]=i;
  sac->queden=nboles;
}

int Bola(sacboles *sac){
  /* Takes the next ball from the bag, at random */
  int i,quina;

  if (sac->queden==0) return(-1);
  quina = (rand()%sac->queden);
  for (i=quina;i<sac->queden-1;i++) sac->llista[i]=sac->llista[i+1];
  sac->queden--;
  return(quina);
}

LlistaSac(sacboles sac){
  /* List the balls in the bag */
  int i;
  printf("%d (",sac.queden);
  for (i=0;i<sac.queden;i++) printf("%d ",sac.llista[i]);
  printf(")\n");
}

Llistat(FILE *sort,int matriu[][MAXB],int vfils,int bcols){
  /* List the incidence matrix, ones per row, coocurrences, etc. */

  int i,j,kk;
  int ri,lamik;
  int kj[MAXB];
  int lams[MAXV][MAXV];
  int xij;
  int nuns;
  int quads;

    fprintf(sort,"\n");
    nuns=0; quads=0;
    for (i=0;i<vfils;i++) {
      for (kk=i+1;kk<vfils;kk++) {
	lamik=0;
	for (j=0;j<bcols;j++) {
	  lamik+=matriu[i][j]*matriu[kk][j];
	}
	lams[i][kk]=lamik;
	quads+=(int)(0.5*lamik*(lamik-1));
      }
    }
    for (j=0;j<bcols;j++) kj[j]=0;
    for (i=0;i<vfils;i++) {
      ri=0;
      for (j=0;j<bcols;j++) {
	xij=matriu[i][j];
	ri+=xij;
	kj[j]+=xij;
	nuns+=xij;
	fprintf(sort,"%2d",xij);
      }
      fprintf(sort,"%4d  ",ri);
      for (kk=i+1;kk<vfils;kk++) {
	fprintf(sort,"%2d",lams[i][kk]);
      }
      fprintf(sort,"\n");
    }
    fprintf(sort,"\n");
    for (j=0;j<bcols;j++) fprintf(sort,"%2d",kj[j]);
    fprintf(sort,"%8d",quads);
    fprintf(sort,"\n\n");
}

int CheckRowCol(int matriu[][MAXB],int vfils,
		int bcols,descrdis descr){
  /* Checks that the number of ones per row and column are correct */

  int i,j;
  int ri;
  int kj[MAXB];
  int xij;
  int v0,b0;

    v0=0;
    for (j=0;j<bcols;j++) kj[j]=0;
    for (i=0;i<vfils;i++) {
      ri=0;
      for (j=0;j<bcols;j++) {
	xij=matriu[i][j];
	ri+=xij;
	kj[j]+=xij;
      }
      if (ri==descr.erre+1) v0++;
      else if (ri!=descr.erre) return(0);
    }
    if (v0!=descr.v0) return(0);
    b0=0;
    for (j=0;j<bcols;j++) 
      if (kj[j]==descr.ka+1) b0++;
      else if (kj[j]!=descr.ka) return(0);
    if (b0!=descr.b0) return(0);
    return(1);
}

void CopiaMatriu(int origen[][MAXB],int desti[][MAXB],int vfils,int bcols){
  /* Copies the incidence matrix origen into the incidence matrix desti,
     both with vfils rows and bcols columns */
  int i,j;
  for (i=0;i<vfils;i++)
    for (j=0;j<bcols;j++)
      desti[i][j]=origen[i][j];
}

int InitMatriu(int matriu[][MAXB],int vfils,int bcols,int uuns){
  /* Inicializes and incidence matrix with the right number of ones
     per row and column. Returns the total number of quadruples. */

  int u,i,j,kk;
  int quads;
  int lamik;
  /* int lams[MAXV][MAXV];*/

  for (i=0;i<vfils;i++)
    for (j=0;j<bcols;j++) matriu[i][j]=0;

  i=0;j=0;
  for (u=0;u<uuns;u++){
    matriu[i][j]++;
    i=(i+1)%vfils;
    j=(j+1)%bcols; 
    if (matriu[i][j]==1)j++;
  }

  quads=0;
  for (i=0;i<vfils-1;i++) {
    for (kk=i+1;kk<vfils;kk++) {
	lamik=0;
	for (j=0;j<bcols;j++) {
	  lamik+=matriu[i][j]*matriu[kk][j];
	}
	/* lams[i][kk]=lamik; */
	quads+=(int)(0.5*lamik*(lamik-1));
    }
  }
  return(quads);
}

void UpdateMatr(int matriu[][MAXB],
     int ii1,int ii2,int jj1,int jj2){ 
     /* Updates the matrix with a transformation of the kind 
	(10,01)->(01,10) */

  matriu[ii1][jj1]=0;
  matriu[ii1][jj2]=1;
  matriu[ii2][jj1]=1;
  matriu[ii2][jj2]=0;
}

int DeltaQuad(int matriu[][MAXB],int vfils,int bcols,
     int ii1,int ii2,int jj1,int jj2){
  /* Computes the increment in quadruples that *would* 
     take place if the proposed transformation was taken */

  int i,j;
  int dquads=0;

  for (i=0;i<vfils;i++)
    if ((i!=ii1)&&(i!=ii2))
      for (j=0;j<bcols;j++)
        if ((j!=jj1)&&(j!=jj2))
          if (matriu[i][j]){
            if (matriu[ii1][j]&&!matriu[ii2][j]&&
	      matriu[i][jj1]&&!matriu[i][jj2])dquads--;
            if (!matriu[ii1][j]&&matriu[ii2][j]&&
	      matriu[i][jj1]&&!matriu[i][jj2])dquads++;
            if (matriu[ii1][j]&&!matriu[ii2][j]&&
	      !matriu[i][jj1]&&matriu[i][jj2])dquads++;
            if (!matriu[ii1][j]&&matriu[ii2][j]&&
	      !matriu[i][jj1]&&matriu[i][jj2])dquads--;
      }
  return(dquads);
}

void Descriptors(int vfils,int bcols,int uuns,descrdis *des){
  /* Given the parameters, computes the descriptors of the design */

  int efe,pvstar;

  des->erre=uuns/vfils;
  des->v0=uuns%vfils;
  des->ka=uuns/bcols;
  des->b0=uuns%bcols;
  efe=(int)(0.5*vfils*(vfils-1));
  pvstar=(int)(0.5*bcols*des->ka*(des->ka-1)+des->b0*des->ka);
  des->lambda=pvstar/efe;
  des->f0=pvstar%efe;
  des->Qstar=(int)(0.5*efe*des->lambda*(des->lambda-1)+des->f0*des->lambda);
}

main() 
{
     int vfils,bcols,uns; /* Number of rows, columns and ones */
     int MInc[MAXV][MAXB]; /* Incidence Matrix */
     int bestMInc[MAXV][MAXB]; /* Thes best Incidence Matrix so far */
     descrdis descr; /* Design descriptors */
     int Quads,Dquads; /* Absolut and incremetal number of quadruples */
     int bestQuads; /* The number of quadruples of the best local minimum
		       so far */
     int iter,niters; /* Iteration number, number of iterations */ 
     int i1,i2,j1,j2;
     int ii1,ii2,jj1,jj2;
     int updates,nmins; /* Updates per iteration, number of local minima */
     sacboles sac1,sac2,sac3,sac4;

     seed=time(0);
     srand(seed);
    
     printf("\nGENERATION OF MBMUDs. http://people.ac.upc.es/pau/mbmuds \n\n");

     printf("Number of rows? v=");
     scanf("%d%*c",&vfils);
     printf("Number of columns? b=");
     scanf("%d%*c",&bcols);
     printf("Number of ones? u=");
     scanf("%d%*c",&uns);
     printf("Maximum number of iterations? Niter=");
     scanf("%d%*c",&niters);

     /* Searches until a MBMUD is found or the maximum number
	of iterations (niters) is exceeded. When a local minimum
	is found the incidence matrix MInc is randomly reinicialized.
        At the end, the best configuration found is listed (either
        the optimal or the best local minimum). */
    
     Descriptors(vfils,bcols,uns,&descr);
     Quads=InitMatriu(MInc,vfils,bcols,uns);
     bestQuads=Quads;
     CopiaMatriu(MInc,bestMInc,vfils,bcols);
     printf("\ndiss(%d,%d,%d) %d\n",vfils,bcols,uns,descr.Qstar);
     nmins=0;
     for (iter=1;iter<=niters && Quads>descr.Qstar;iter++){
       updates=0;
       /* At every iteration tries all possible quadruples */
       for (i1=0,SacNou(vfils,&sac1);i1<vfils;i1++){
         ii1=Bola(&sac1);
	 for (j1=0,SacNou(bcols,&sac2);j1<bcols;j1++){
	   jj1=Bola(&sac2);
	   for (i2=0,SacNou(vfils,&sac3);i2<vfils;i2++){
	     ii2=Bola(&sac3);
             if (ii2!=ii1){
	       for (j2=0,SacNou(bcols,&sac4);j2<bcols;j2++){
	         jj2=Bola(&sac4);
                 if (jj2!=jj1){
		   /* If the quadruple is of the kind (10,01),
		     if the transformation reduces Quads, it is taken */
		   if (MInc[ii1][jj1]&&!MInc[ii1][jj2]&&
		        !MInc[ii2][jj1]&&MInc[ii2][jj2]){
		     Dquads=DeltaQuad(MInc,vfils,bcols,ii1,ii2,jj1,jj2);
		     if (Dquads<0) {
		       Quads+=Dquads;
		       UpdateMatr(MInc,ii1,ii2,jj1,jj2);
		       updates++;
		     }
		   }
	         }
	       }
	     }
	   }
	 }
       }
       if (updates==0 || Quads==descr.Qstar) {
	 /* if local minimum */
         nmins++;
	 if (Quads<bestQuads){
	   bestQuads=Quads;
	   CopiaMatriu(MInc,bestMInc,vfils,bcols);
	 }
       }
       printf("it=%3d, updts=%3d, nmins=%2d, Q=%3d, bestQ=%3d, Q*=%3d\n",
	      iter,updates,nmins,Quads,bestQuads,descr.Qstar);
       if (updates==0){
	 /* restart if not optimal */
         Quads=InitMatriu(MInc,vfils,bcols,uns);
       }
     }
     Llistat(stdout,bestMInc,vfils,bcols);
     if (bestQuads==descr.Qstar) 
       fprintf(stdout," OPTIMAL\n\n");
     printf("Press any key to exit\n");
     scanf("%*c"); 
}