/*  This is -*- C -*- code ! */
/* ************************************************************
   file 'tabusearch.c'
   tabu-search based metaheuristics for MIP
   Copyright 2004/01
   Joao Pedro Pedroso, DCC-FC and LIACC
   Universidade do Porto, Portugal

   http://www.ncc.up.pt/~jpp/mipts/tabusearch.c
   ************************************************************ */

#include "mip.h"

#define LOG 0	// whether to display verbose output or not

void construct(LPPROB *prob, int x[], double *z, double *zeta) {
  if (LOG)
    printf("constructing...\n");

  solvelp(prob, z, zeta);

  if (LOG)
    printf("%lf \nfixing...\n",*z);
  
  /* fix variables in a random order*/
  int p[prob->nints];
  permutation(p, prob->nints);
  roundfix(prob, x, z, zeta, p, prob->nints);

}

void diversify(LPPROB *prob, int x[], double *z, double *zeta, int tabu[]) {
  int p[prob->nints];
  permutation(p, prob->nints);
  int size = 1 + rand()%prob->nints;
  
  //if (LOG)
    printf("diversifying... releasing %d variable\n", size);

  /* release 'size' variables */
  int j;
  for (j=0;j<size;j++) {
    int k = p[j];
    lpx_set_col_bnds(prob->lp, prob->lpi[k], LPX_DB, prob->lb[k], prob->ub[k]);
  }

  /* fix variables again */
  roundfix(prob, x, z, zeta, p, size);

  /* reset tabu list */
/*   for (j=0; j<prob->nints; j++) */
/*     tabu[j] = -prob->nints; */

}

void intensify(LPPROB *prob, int x[], double *zin, double *zetain, int tabu[], int iter) {
  double z = *zin;
  double zeta= *zetain;

  //  if (LOG)
  printf("starting intensifying (CPU=%g)\n", get_cpu());
  // printivectz(x, prob->nints, z, zeta);
  // assert(checksolution(prob, x, z, zeta));

  int p[prob->nints];	  // indices for tabu variables
  int q[prob->nints];	  // indices for non-tabu variables

  int j, h = 0, n = 0;
  for (j=0; j< prob->nints; j++)
    if ((iter-tabu[j]) > prob->nints)	// !!! prob->nints -> tabu tenure
      p[h++] = j;
    else
      q[n++] = j;
  random_shuffle(p, h);
  random_shuffle(q, n);

  int size = h;		// h -> number of tabu variables
  printf("freeing %d variables\n", size);

  /* release 'size' non-tabu variables */
  for (j=0;j<size;j++) {
    int k = p[j];
    lpx_set_col_bnds(prob->lp, prob->lpi[k], LPX_DB, prob->lb[k], prob->ub[k]);
  }

  /* check if LP became feasible */
  lpx_set_int_parm(prob->lp, LPX_K_MSGLEV, 2);    	// silent/verbose mode
  solvelp(prob, &z, &zeta);

  /* if LP is not feasible yet, release more variables (until it becames feasible) */
  int super=0;	// # added variables which ARE tabu
  while (zeta != 0) {
    if (super == n)	// reached the last variable; problem is infeasible
      assert("problem is infeasible");
    int k = q[super++];
    lpx_set_col_bnds(prob->lp, prob->lpi[k], LPX_DB, prob->lb[k], prob->ub[k]);
    solvelp(prob, &z, &zeta);
  }
  printf("freeing %d extra variables\n", super);
  
  /* do a branch-and-bound on the released variables */
  if (zeta == 0.) {
    lpx_set_real_parm(prob->lp, LPX_K_TMLIM, prob->nints);	// time limit = nints
    int ret = lpx_integer(prob->lp);
    int stat = lpx_mip_status(prob->lp);
    lpx_set_real_parm(prob->lp, LPX_K_TMLIM, -1);
    printf("\t(mip stat: %d;  return value: %d )\n", stat, ret );
    if (ret == LPX_E_OK && (stat == LPX_I_OPT || stat == LPX_I_FEAS )) {
      // subproblem was solved to optimality
      z = lpx_mip_obj_val(prob->lp);
      printf("intensifying: exiting (found mip sol) (CPU=%g)\n", get_cpu());
      *zin = z;
      *zetain = zeta;
      printivectz(x, prob->nints, *zin, *zetain);
      lpx_set_int_parm(prob->lp, LPX_K_MSGLEV, 0);    	// verbose mode

      // fix the released vars at their new values
      for (j=0;j<size;j++) {
	int k = p[j];
	x[k] = lround(lpx_mip_col_val(prob->lp, prob->lpi[k]));
	lpx_set_col_bnds(prob->lp, prob->lpi[k], LPX_FX, x[k], x[k]);
      }
      for (j=0;j<super;j++) {
	int k = q[j];
	x[k] = lround(lpx_mip_col_val(prob->lp, prob->lpi[k]));
	lpx_set_col_bnds(prob->lp, prob->lpi[k], LPX_FX, x[k], x[k]);
      }
      return;
    }
  }
  
  // mip solution did not work;
  // fix released variables again....
  for (j=0;j<size;j++) {
    int k = p[j];
    lpx_set_col_bnds(prob->lp, prob->lpi[k], LPX_FX, x[k], x[k]);
  }
  for (j=0;j<super;j++) {
    int k = q[j];
    lpx_set_col_bnds(prob->lp, prob->lpi[k], LPX_FX, x[k], x[k]);
  }

  printf("intensifying: exiting... (CPU=%g)\n", get_cpu());
  // printivectz(x, prob->nints, z, zeta);
  lpx_set_int_parm(prob->lp, LPX_K_MSGLEV, 0);    	// silent mode 

}

//
// choose a candidate move using a tabu list
//
int tabumove(LPPROB *prob, PARAM *par, double zbest, double zetabest,
	     double *zcurr, double *zetacurr, int x[],
	     int tabu[], int iter) {
  double z = *zcurr;		// current (entering) solution's objective
  double zeta = *zetacurr;	// current (entering) solution's constraint violation
  int xcand[prob->nints];	// best non-tabu move	!!!!!!!!!!!!!! should be replaced by a scalar, x[move] !!!!!!!!
  double zcand = Infinity;	// best move's objective
  double zetacand = Infinity;	// best move's constraint violation
  int alltabu = 1;	// to know if we are blocked, i.e. all possible moves are tabu
  int move = -1;	// candidate move
  int j, k, len, step;

  // backup starting solution ??????????? required??????????
  memcpy(xcand, x, prob->nints*sizeof(int));

  for (j=0; j< prob->nints; j++)
    if ((iter-tabu[j]) > prob->nints) {	// !!! prob->nints -> tabu tenure
      alltabu = 0;
      break;
    }
  if (alltabu) {
    printf("all variables tabu, blocked\n");
    j = rand()%prob->nints;
    printf("randomly choosing variable %d to change\n", j);
    x[j] = rand() % (prob->ub[j] - prob->lb[j] + 1) + prob->lb[j];
    lpx_set_col_bnds(prob->lp,prob->lpi[j],LPX_FX,x[j],x[j]);
    
    solvelp(prob, &z, &zeta);
    tabu[j] = iter;
    *zcurr = z;
    *zetacurr = zeta;
    if (zeta < zetabest || (zeta==zetabest && z < zbest))
      return 1;
    return 0;
  }

  // permutation: useful for depth first
  // to avoid having always the same vars being explored
#ifdef DEPTH_FIRST
  int p[prob->nints];
  permutation(p, prob->nints);
#endif

  for (k=0; k < prob->nints; k++) {
    j = k;
#ifdef DEPTH_FIRST
    j = p[k];
#endif

    int tmp = x[j]; 

    // N1 neighborhood: add or subtract 1 to x[j]
    for (step=-1; step<=1; step+=2) {
      x[j] += step;
      if (prob->lb[j] <= x[j]  &&  x[j] <= prob->ub[j]) {
	lpx_set_col_bnds(prob->lp,prob->lpi[j],LPX_FX,x[j],x[j]);
	// printf("fixing x[%d] to %d;\n",j, x[j]);
	
	
	solvelp(prob, &z, &zeta);
	if (zeta < zetabest || (zeta==zetabest && z < zbest)) {	// aspiration criterion
	  // printf("TS^; x=");
	  // printivectz(x, prob->nints, z, zeta);
	  *zcurr = z;
	  *zetacurr = zeta;
	  return 1;
	}

	len = (par->tenure == 0 ?
	       (int) (prob->nints * (1.*rand()/RAND_MAX)) :  /* !!! prob->nints-->tabu tenure */
	       par->tenure);

	if ( (iter-tabu[j])>len /* not tabu */  &&
	     (zeta < zetacand || (zeta==zetacand && z < zcand)) /* solution is beter */) {
	  
	  move = j;
	  memcpy(xcand, x, prob->nints*sizeof(int)); 
	  zcand = z;
	  zetacand = zeta;
	  //	  if (LOG) {
	  //  printf("current best move: x[%d] -> %d\t%f\t%f\n", j, x[j], zcand, zetacand);
	  //  printivectz(x, prob->nints, z, zeta);
	  //	  }
	  
	  // depth first: return as soon as a non-tabu improving move is found
#ifdef DEPTH_FIRST
	  if (zeta < *zetacurr || (zeta==*zetacurr && z < *zcurr)) {	/* improved curr sol */
	    //	    printf("breaking\n");
	    break;
	  }
#endif
	}

	// printf("fixing x[%d] to %d;\n",j, tmp);
	lpx_set_col_bnds(prob->lp,prob->lpi[j],LPX_FX,tmp,tmp);

      }
      
      x[j] -= step;
    }
  }

  /*  update tabu list */
  assert(move >= 0  &&  move < prob->nints);
  tabu[move] = iter;

  // copy best candidate into tabu 
  memcpy(x, xcand, prob->nints*sizeof(int)); 
  *zcurr=zcand;
  *zetacurr=zetacand;
  lpx_set_col_bnds(prob->lp, prob->lpi[move], LPX_FX, x[move], x[move]);
/*   printf("move is %d, x[]=%d\n", move, x[move]); */
/*   printf("fixing x[%d] to %d;\n", move, x[move]); */

  if (LOG) {
    printf("iter %d, var %d made tabu\n", iter, move);
    printivectz(x, prob->nints, zcand, zetacand);
  }

    //    assert(checksolution(prob, x, *zcurr, *zetacurr));

  return 0;
}  

//
// tabu search ======================================
// do a tabu search on integer variables using the N1 neighborhood
//
void tabusearch(LPPROB *prob, PARAM *par, double *z, double *zeta, int x[]) {
  double zbest;
  double zetabest;
  int best[prob->nints];
  int tabu[prob->nints];
  int j;
  for (j=0; j<prob->nints; j++)
    tabu[j] = -prob->nints;
  
  // // fix all integer variables
  // printf("\nstarting tabu search; fixing variables:\n  ");
  // for (j=1; j<= prob->nints; j++) {
  //   // x[j]=lpx_get_mip_col(prob->lp,j);
  //   lpx_set_col_bnds(prob->lp,j,LPX_FX,x[j],x[j]);
  //   // lpx_set_col_kind(prob->lp,j,LPX_IV);     // integer variable
  //   if (LOG)
  //     printf("%d\t%d\n", j, (int) x[j]);
  // }
  
  zbest = *z;
  zetabest = *zeta;
  if (LOG)
    printf("z->%f, zeta->%f\n", *z, *zeta);
  memcpy(best,x,(prob->nints)*sizeof(int));
  
  int i=0;
  
  // measures on a search after a construction/diversification
  int consec = 0;	// consecutive iterations with decreasing quality
  double zold = *z;	// previous objective
  double zetaold = *zeta;	// previous violation
  int xold[prob->nints];
  memcpy(xold,x,(prob->nints)*sizeof(int));

  // printf("beginning tabu search...\n");
  // printivectz(x, prob->nints, *z, *zeta);
  while (i++ < par->niter) {
    // printf("%d [consec:%d]: ", i, consec);
    if (tabumove(prob, par, zbest, zetabest, z, zeta, x, tabu, i)) {
      // if improved, update best
      memcpy(best,x,(prob->nints)*sizeof(int));
      zbest=*z;
      zetabest=*zeta;
      reportsol(*z, *zeta);
    }
    
    // measures for diversification/intensification
    if (*zeta < zetaold || (*zeta==zetaold && *z < zold)) {
      // improved best solution ON THIS SEARCH thread
      consec = 0;
      zold = *z;
      zetaold = *zeta;
      memcpy(xold,x,(prob->nints)*sizeof(int));
    }
    else 
      consec++;
    
    if (consec == prob->nints) { // !!!par->ndivers) {
      // try an intensification on the best solution of this thread
      update_sol(prob, x, z, zeta, xold, zold, zetaold);
      
      printf("%d\t(%dV)\tx=", i, consec);
      printivectz(x, prob->nints, *z, *zeta);

      intensify(prob, x, z, zeta, tabu, i);
    
      printf("%d\t(%dV)\tx=", i, consec);
      printivectz(x, prob->nints, *z, *zeta);
    }
    else if (consec > prob->nints) { // !!!par->ndivers) {
      printf("%d\t(%dV)\tx=", i, consec);
      printivectz(x, prob->nints, *z, *zeta);

      consec = 0;
      /* par->ndivers++; not used in this version */
      // printf ("*** current diversification: %d\n", par->ndivers);
      // restorebounds(prob);
      // construct(prob, x, z, zeta);
      diversify(prob, x, z, zeta, tabu);
      // not necessary: update_sol(prob, x, z, zeta, xold, zold, zetaold);

      printf("%d\t(%dV)\tx=", i, consec);
      printivectz(x, prob->nints, *z, *zeta);

    }
    
    // check solution improvements on the intens/divers
    if (*zeta < zetabest || (*zeta==zetabest && *z < zbest)) {	/* solution is beter */
      // improved best solution, update best
      memcpy(best,x,(prob->nints)*sizeof(int));
      zbest=*z;
      zetabest=*zeta;
      reportsol(*z, *zeta);

      // update also thread's solution
      consec = 0;
      zold = *z;
      zetaold = *zeta;
      memcpy(xold,x,(prob->nints)*sizeof(int));
    }

    // printf("%d\t(%dV)\tx=", i, consec);
    // printivectz(x, prob->nints, *z, *zeta);

  }
  
  // printf("final TS solution:");
  // printivectz(best, prob->nints, zbest, zetabest);

  // update x[] and z, and return
  memcpy(x, best, prob->nints*sizeof(int)); 
  *z=zbest;
  *zeta=zetabest;
}


int xtabusearch(LPPROB *prob, PARAM *par, double *z, double *zeta, int x[]) {
  // read problem and initialize data
  /* not used in this version: 
     par -> ndivers = 1;	
     par -> nintens = 1;
  */

  // find initial solution
  construct(prob, x, z, zeta);
  reportsol(*z, *zeta);
  
  
  // // do a local+tabu search on the integer variables using N1 neighborhood
  tabusearch(prob, par, z, zeta, x);

  reportsol(*z, *zeta);
  
  printf("solution: "); printivectz(x, prob->nints, *z, *zeta);

  return 0;
}

int main(int argc, char *argv[]) {
  if (argc != 5  &&  argc != 6) {
    printf("usage: '%s seed niter len filename [filename.dat]'\n", argv[0]);
    printf("where:\n");
    printf("  seed -> initializer for random number generation (0 -> system default's seed)\n");
    printf("  niter -> number of tabu search iterations \n");
    printf("  len -> tabu tenure (0 for a random number in [1..nintegers]\n");
    printf("  filename -> name of the instance file:\n");
    printf("       *.mps -> mps format\n");
    printf("       *.mod -> ampl/gnu mathprog model)\n");
    printf("                (=> filename.dat -> NULL, or corresponting data)\n");
    printf("       *.lp -> cplex lp format\n");
    printf("\n");
    exit(-1);
  }

  int randseed = atoi(argv[1]);
  printf("random seed: %d\n", randseed);
  if (randseed > 0)
    srand(randseed);

  // read problem and initialize data
  LPPROB *prob = (LPPROB *) malloc(sizeof(LPPROB));
  PARAM *par = (PARAM *) malloc(sizeof(PARAM));
  char *modfile = argv[4];
  char *datfile = NULL;
  if (argc == 6)
    datfile = argv[5];
  initprob(modfile, datfile, prob);
  par->niter = atoi(argv[2]);	// tabu search iterations
  par->tenure = atoi(argv[3]);	// tabu tenure

  // construction + tabu search
  int x[prob->nints];
  double z, zeta;
  xtabusearch(prob, par, &z, &zeta, x);
  printivectz(x, prob->nints, z, zeta); 
  delprob(prob);
  
  return 0;
}