/*  This is -*- C -*- code ! */
/* ************************************************************
   file 'mip.c'
   utility functions for MIP related routines
   Copyright 2004/01
   Joao Pedro Pedroso, DCC-FC and LIACC
   Universidade do Porto, Portugal

   http://www.ncc.up.pt/~jpp/mipts/mip.c
   ************************************************************ */
#include "mip.h"

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

//
// evaluate infeasibilities
// (stolen from glpk internal code)
//
double eval_lp_sum(LPX *lp) {
  int m, n, i, j, k, type;
  double lb, ub, x, sum;
  m = lpx_get_num_rows(lp);
  n = lpx_get_num_cols(lp);
  sum = 0.0;
  for (k = 1; k <= m+n; k++) {
    if (k <= m) {
      i = k;
      lpx_get_row_bnds(lp, i, &type, &lb, &ub);
      x = lpx_get_row_prim(lp, i);
    }
    else {
      j = k - m;
      lpx_get_col_bnds(lp, j, &type, &lb, &ub);
      x = lpx_get_col_prim(lp, j);
    }
    switch (type) {
    case LPX_FR:
      break;
    case LPX_LO:
      if (x < lb) sum += lb - x;
      break;
    case LPX_UP:
      if (x > ub) sum += x - ub;
      break;
    case LPX_DB:
    case LPX_FX:
      if (x < lb) sum += lb - x;
      if (x > ub) sum += x - ub;
      break;
    default:
      assert(type != type);
    }
  }
  return sum;
}

void solvelp(LPPROB *prob, double *z, double *zeta) {

  // this LP evaluation version is nice, as it does not require
  // artificial variables.  
  // however, it has some problems: the same (infeasible)
  // fixed variables may stop in different values of the 
  // violation/objective, depending on the starting solution

  lpx_simplex(prob->lp);
  *z=lpx_get_obj_val(prob->lp);
  int status = lpx_get_status(prob->lp);
  if (status == LPX_OPT)	// to avoid avoid precision problems
    *zeta = 0.;
  else if (status == LPX_NOFEAS) {
    *zeta = eval_lp_sum(prob->lp);
  }
  else {
    // the status is LPX_UNDEF
    // likely, there are no free variables
    printf("warning:  lp solution with undefined status %d\n", status);
    assert(status == LPX_UNDEF);
    *zeta = eval_lp_sum(prob->lp);
  }
  if (0) { // (LOG)
    int j;
    for (j=0;j<prob->nints;j++) {
      printf("%d: %f <= %f <= %f\n", 
	     j, lpx_get_col_lb(prob->lp, prob->lpi[j]), 
	     lpx_get_col_prim(prob->lp, prob->lpi[j]), 
	     lpx_get_col_ub(prob->lp, prob->lpi[j]));
    }
    printf("solvelp: z=%lf, zeta=%lf\n", *z, *zeta);
  }
  return;
}


//
// initialize problem's data
//
void initprob(char *modfile, char *datfile, LPPROB *prob) {
  // read problem from files
  if (strstr(modfile, ".mps") != NULL)
    prob->lp = lpx_read_mps(modfile);
  else if  (strstr(modfile, ".mod") != NULL)
    prob->lp = lpx_read_model(modfile,datfile,NULL);
  else if  (strstr(modfile, ".lp") != NULL)
    /* for version 4.7 use prob->lp = lpx_read_cpxlp(modfile); */
    prob->lp = lpx_read_cpxlp(modfile);
  else {
    fprintf(stderr, "unkown input format, file %s", modfile);
    exit(-1);
  }
    

  lpx_set_int_parm(prob->lp, LPX_K_MSGLEV, 0);  // silent /verbose modes
  lpx_set_int_parm(prob->lp, LPX_K_DUAL, 1);    // use dual simplex if sol is dual feasible
  prob->nints = lpx_get_num_int(prob->lp);
  prob->nvars = lpx_get_num_cols(prob->lp);
  prob->ncons = lpx_get_num_rows(prob->lp);

  int j, n=0;
  prob->idx = (int *) malloc(sizeof(int) * (lpx_get_num_cols(prob->lp)+1));
  prob->lpi = (int *) malloc(sizeof(int) * prob->nints);	
  prob->lb  = (int *) malloc(sizeof(int) * prob->nints);
  prob->ub  = (int *) malloc(sizeof(int) * prob->nints);

  for (j=1; j<=prob->nvars; j++) {
    if (lpx_get_col_kind(prob->lp, j) == LPX_IV) {
      if (lpx_get_col_type(prob->lp, j) != LPX_DB) {
	printf("cannot handle integer variable without double bounds\n");
	exit(-1);
      }
      
      prob->idx[j] = n;
      prob->lpi[n] = j;
      prob->lb[n] = lround(lpx_get_col_lb(prob->lp, j));
      prob->ub[n] = lround(lpx_get_col_ub(prob->lp, j));
      n++;
    }
  }

  assert(n == prob->nints);
}

void delprob(LPPROB *prob) {
  free(prob->idx);
  free(prob->lpi);
  free(prob->lb );
  free(prob->ub );
  lpx_delete_prob(prob->lp);
  free(prob);
}  

void update_sol(LPPROB *prob,
		int x[], double *z, double *zeta,
		int xnew[], double znew, double zetanew) {
  *z = znew;
  *zeta = zetanew;
  int j;
  for (j=0; j<prob->nints; j++) {
    x[j] = xnew[j];
    lpx_set_col_bnds(prob->lp, prob->lpi[j], LPX_FX, x[j], x[j]);
  }
}

void roundfix(LPPROB *prob, int x[], double *z, double *zeta,
	      int p[], int size) {
  int j,k;
  for (j=0;j<size;j++) {
    k = p[j];	// get indices in random order
    
    double sol = lpx_get_col_prim(prob->lp,prob->lpi[k]);
    int isol = (int) sol;
    double u = uniform();
    //      printf("sol: %lf, u: %lf, isol:%d, rounding ", sol, u, isol);
    if (u < sol-isol) {
      x[k] = isol + 1;
      //	printf("up\n");
    }
    else { 
      x[k] = isol;
      //	printf("down\n");
    }
    x[k] = min(x[k], prob->ub[k]);
    x[k] = max(x[k], prob->lb[k]);
    
    lpx_set_col_bnds(prob->lp, prob->lpi[k], LPX_FX, x[k], x[k]);
    if (LOG)
      printf(" %i (%d), %s -> \t %d <= %d <= %d\n", k, prob->lpi[k],
	     lpx_get_col_name(prob->lp, prob->lpi[k]),
	     prob->lb[k], x[k], prob->ub[k]);
    
    // progressive fixing: solve lp again...
    solvelp(prob, z, zeta);
    // lpx_print_sol(prob->lp, "lpsol.sol");
    if (LOG)
      printivectz(x, prob->nints, *z, *zeta);
  }
}

//
//	utility functions
//

// random number in [0,1] with uniform distribution:
double uniform() { return rand()/(double) RAND_MAX; }

// place the elements of a vector in random order
void random_shuffle(int vector[], int size) {
  if (size<=1) return;
  int i;
  for (i=0; i<size-1; i++) {
    int j = i + 1 + rand()%(size-i-1);
    int tmp = vector[i];
    vector[i] = vector[j];
    vector[j] = tmp;
  }
}

// obtain a random permutation of the numbers 0,1,2,3,...,size-1
void permutation(int vector[], int size) {
  int i;
  for (i=0; i<size; i++)
    vector[i] = i;
  random_shuffle(vector, size);
}



//
// CPU information
//

// cpu time since the begin of the program, in seconds:
double get_cpu() {
  struct tms buf;
  times(&buf);
  return (double) buf.tms_utime / sysconf(_SC_CLK_TCK);
}
  


//
// printing utility functions
//
void reportsol(double z, double zeta) {
  double cpu = get_cpu();
  printf("cpu,z,zeta:\t%f\t%f\t%f\n", cpu, z, zeta);
}

void printivectz(int x[], int len, double z, double zeta) {
  int i=-1;
  while (++i < len)
    printf("%d ", (int) *(x+i));
  printf("\t%12lg", z);
  printf("\t%12lg", zeta);
  printf("\n");
}