/* ======================================================================
             ipsolver.c   O2-opgave DAT2A  2002    David Pisinger
   ====================================================================== */

/* Dette C program loeser generel heltalsprogrammering 
 *
 *    maximize    sum_{j=1,2,..,n} c_j x_j
 *    subject to  sum_{j=1,2,..,n} a_ij x_j <= b_i   for i=1,..,m
 *                 x_j heltal >= 0
 * hvor
 *   c = (c_j) er en omkostningsvektor af laengde n
 *   b = (b_i) er hoejresiderne givet ved en vektor af laengde m
 *   a = (a_ij) er en matrix af format n x n
 *
 * programmet oversaettes paa DIKU's linux pc'er med
 *
 *   gcc -ansi -o ipsolver ipsolver.c -lm
 *
 * naar tests skal koeres kan option -O4 (kodeoptimering) tilfoejes.
 * Den udfoerbare kode kaldes med
 *
 *   ipsolver indfil
 *
 * hvor "indfil" er et filnavn som overholder foelgende syntax (vist ved
 * eksempel)
 *   
 *    maximize    2 x1 - 3 x2 + 3 x3        \ kommentar: Cormen s. 781
 *    subject to    x1 +   x2 -   x3 <= 7
 *                 -x1 -   x2 +   x3 <= -7
 *                  x1 - 2 x2 + 2 x3 <= 4
 *    end
 * 
 * Det underfostaas at alle beslutningsvariable er ikke-negative.
 */

#include <stdlib.h>
#include <stdio.h>
#include <stdarg.h>
#include <string.h>
#include <math.h>
#include <sys/times.h>
#include <unistd.h>


/* ======================================================================
				 type declarations
   ====================================================================== */

#define MAXM       500
#define MAXN       500
#define EPSILON    1E-10
#define INFTY      1E+10

#define FALSE      0
#define TRUE       1

#define FEASIBLE   0          /* LP problem is feasible */
#define INFEASIBLE 1          /* LP problem is infeasilbe */
#define UNBOUNDED  2          /* LP problem is unbounded */

#define MAXSYM     100        /* max length of a variable */
#define MAXHASH    1007       /* size of hash table for parser */
#define FEOF       127        /* end of file */

#define MAXIMIZE   1          /* symbols in parser */
#define SUBJECT    2
#define TO         3
#define END        4
#define LESSEQUAL  5
#define PLUS       6
#define MINUS      7
#define VARNAME    8
#define VALUE      9
#define UNKNOWN    10

#define MAXPRINT   10         /* max number of columns to be printed nice */

#define rint(x)    (floor((x)+0.5))


/* ======================================================================
				  structure
   ====================================================================== */

typedef struct hashr {
   char         *name;
   int          no;
   struct hashr *nxt;
} hashrec;

typedef hashrec *hashptr;
typedef int boolean;
typedef double matrix[MAXM][MAXN];
typedef double mvector[MAXM];
typedef double nvector[MAXN];
typedef int indices[MAXM];


/* ======================================================================
				  global variables
   ====================================================================== */

int n, m;         /* problem has m rows, n columns */
matrix  a;        /* matrix of coefficients */
mvector b;        /* vector of righ-hand sides */
nvector c;        /* vector of costs */
nvector x;        /* integer solution */
double  zstar;    /* currently best solution */
boolean intcost;  /* objective function coefficients are integers */
int     bbnodes;  /* number of branch and bound nodes searched */
  

/* ======================================================================
                           internal variables for parser
   ====================================================================== */

hashptr htab[MAXHASH];  
hashptr varptr[MAXN];
int val, sym, varcounter;
double valf;
boolean ended;
FILE *in, *trace;
char ch;


/* =======================================================================
                                  error
   ======================================================================= */

static void error(char *str, ...)
{
  va_list args;

  va_start(args, str);
  vprintf(str, args); vfprintf(trace, str, args);
  printf("\nsee trace.out for additional info\n"); fprintf(trace, "\n");
  va_end(args);
  exit(-1);
}


/* ======================================================================
                                   cpu_time
   ====================================================================== */

/* The following routine is used for measuring the cpu time used
 *   cpu_time(NULL)    - start timer
 *   cpu_time(&t)      - return the elapsed cpu-time in variable t (double).
 */

void cpu_time(double *t)
{
  static struct tms start, stop;

  if (t == NULL) {
    times(&start);
  } else {
    times(&stop);
    *t = (double) (stop.tms_utime - start.tms_utime) / sysconf(_SC_CLK_TCK);
  }
}


/* ======================================================================
   				    nextch
   ====================================================================== */

void nextch(void)
{
  ch = fgetc(in); if (ch == EOF) { ch = ' '; ended = TRUE; }
  fprintf(trace,"%c", ch);
  if ((ch == '\n') || (ended)) { ch = ' '; return; }
  if (ch == '\\') { 
    for (;;) {
      ch = fgetc(in); if (ch == EOF) { ch = ' '; ended = TRUE; }
      fprintf(trace,"%c", ch);
      if ((ch == '\n') || (ended)) { ch = ' '; return; }
    }
  }
}


/* ======================================================================
   				    lookup
   ====================================================================== */

int lookup(char *str)
{
  hashptr h, *hptr;
  char *i;
  int j, hash;

  /* find hash value */
  hash = 0;
  for (i = str, j = 0; *i != 0; i++, j++) hash = (hash + *i * j) % MAXHASH;
  /* search through hash table */
  h = htab[hash]; hptr = &(htab[hash]);
  while (h != NULL) {
    if (strcmp(str, h->name) == 0) return h->no;
    hptr = &(h->nxt); h = h->nxt;
  }
  /* if not found, insert at end */
  varcounter++; 
  *hptr = h = malloc(sizeof(hashrec));
  h->name = malloc(strlen(str)+1);
  h->nxt  = NULL;
  h->no = varcounter; 
  strcpy(h->name, str);
  varptr[varcounter] = h; 
  return (h->no);
}


/* ======================================================================
				    nextsym
   ====================================================================== */

void nextsym(void)
{
  char *i;
  char cur[MAXSYM];  /* current symbol read */
  boolean numeric, alphanum;

  while ((ch == ' ') && (!ended)) nextch();
  /* printf("skipped blank having %c\n", ch); */
  i = cur; numeric = FALSE; alphanum = FALSE;
  if ((ch == '<') || (ch == '+') || (ch == '-')) {
    *i = ch; 
    if (ch == '<') { nextch(); i++; *i = ch; }
    nextch(); i++;
  } else {
    if (((ch >= '0') && (ch <= '9')) || (ch == '.')) {
      for (numeric = TRUE; ; ) {
        *i = ch; nextch(); i++;
        if (i-cur > MAXSYM) error("too long symbol");
        if (((ch < '0') || (ch > '9')) && (ch != '.')) break;
      }
    } else {
      if (((ch >= 'a') && (ch <= 'z')) || ((ch >= 'A') && (ch <= 'Z'))) {
        for (alphanum = TRUE; ; ) {
          *i = ch; nextch(); i++;
          if (i-cur > MAXSYM) error("too long symbol");
          if (((ch < 'a') || (ch > 'z')) && ((ch < 'A') || (ch > 'Z')) &&
              ((ch < '0') || (ch > '9'))) break;
        }
      } 
    }
  } 
  *i = 0; sym = val = valf = 0;
  do {
    if (strcmp(cur,"<=")==0)      { sym = LESSEQUAL; break; }
    if (strcmp(cur,"+")==0)       { sym = PLUS; break; }
    if (strcmp(cur,"-")==0)       { sym = MINUS; break; }
    if (strcmp(cur,"maximize")==0){ sym = MAXIMIZE; break; }
    if (strcmp(cur,"subject")==0) { sym = SUBJECT; break; }
    if (strcmp(cur,"to")==0)      { sym = TO; break; }
    if (strcmp(cur,"end")==0)     { sym = END; break; }
    if (numeric)                  { sym = VALUE; valf = atof(cur); break; }
    if (alphanum)                 { sym = VARNAME; val = lookup(cur); break; }
    sym = UNKNOWN; 
  } while (FALSE);
}


/* ======================================================================
				   readsignconst
   ====================================================================== */

void readsignconst(double *fact)
{
  *fact = 1.0;
  if (sym == PLUS ) { nextsym(); }
  if (sym == MINUS) { *fact = -1.0; nextsym(); }
  if (sym != VALUE) error("expected +, - constant");
  *fact *= valf; nextsym(); 
}


/* ======================================================================
				   readsignconstvar
   ====================================================================== */

void readsignconstvar(int *var, double *fact)
{
  *fact = 1.0;
  if (sym == PLUS) { nextsym(); }
  if (sym == MINUS) { *fact = -1.0; nextsym(); }
  if (sym == VALUE) { *fact *= valf; nextsym(); } 
  if (sym != VARNAME) error("expected +, -, constant or variable");
  *var = val; nextsym();
}


/* ======================================================================
				  readfromfile
   ====================================================================== */

void readfromfile(char *name)
{
  int i, j;
  double v;
 
  in = fopen(name, "r");
  if (in == NULL) error("infile not open");
   
  /* initialize parser */
  printf("reading instance from file: %s\n", name);
  for (i = 0; i < MAXHASH; i++) htab[i] = NULL; 
  varcounter = 0; nextch(); nextsym(); 
 
  if (sym != MAXIMIZE) error("expected 'maximize'"); 
  nextsym();
  for (;;) {
    readsignconstvar(&j, &v);
    c[j] = v;
    if (sym == SUBJECT) break;
    if ((sym == MINUS) || (sym == PLUS)) continue; 
    error("expected +, - or 'subject'"); 
  } 
  if (sym != SUBJECT) error("expected 'subject'"); 
  nextsym();
  if (sym != TO) error("expected 'to'"); 
  nextsym();
  for (i = 1; ; i++) {
    for (;;) {
      readsignconstvar(&j, &v);
      a[i][j] = v;
      if (sym == LESSEQUAL) break;
      if ((sym == MINUS) || (sym == PLUS)) continue; 
      error("expected +, - or <="); 
    }
    nextsym();
    readsignconst(&v);
    b[i] = v; 
    if (sym == END) break;
  }
  if (sym != END) error("expected 'end'"); 

  m = i; n = varcounter;
  printf("instance size m %d, n %d\n", m, n);
  fclose(in);
}   


/* ======================================================================
				  formatvar
   ====================================================================== */

void formatvar(double fac, char *var)
{
  if ((n > MAXPRINT) && (fac == 0)) return;
  if (fac == 0) { fprintf(trace, "           "); return; }
  if (fac > 0) fprintf(trace, " +");
  if (fac < 0) { fprintf(trace, " -"); fac = -fac; }
  fprintf(trace, "%4.2lf %-4s", fac, var);
}


/* ======================================================================
				     varname 
   ====================================================================== */

char *varname(int i)
{
  static char str[MAXSYM];
  if (i <= n) return varptr[i]->name;
  sprintf(str, "_s%d", i-n);
  return str;
}


/* ======================================================================
				  printinstance
   ====================================================================== */

void printinstance(int m, int n)
{
  int i, j;
  double v;

  fprintf(trace, "\nmaximize\n");
  for (j = 1; j <= n; j++) formatvar(c[j], varname(j));
  fprintf(trace, "\nsubject to\n");
  for (i = 1; i <= m; i++) {
    for (j = 1; j <= n; j++) formatvar(a[i][j], varname(j));
    fprintf(trace," <= ");
    fprintf(trace, "%4.2lf\n", b[i]);
  }
  fprintf(trace, "end\n");
  v = 0.0;
  for (j = 1; j <= n; j++) {
    fprintf(trace, "%-4s = %4.2lf\n", varname(j), x[j]);
    v += x[j] * c[j];
  }
  fprintf(trace, "v    = %4.2lf\n", v);
}


/* ======================================================================
				  printlp
   ====================================================================== */

void printlp(int m, int n, indices basis, matrix a, mvector b, nvector c, 
             double v)
{
  int i, j;

  if (n > MAXPRINT) return;
  fprintf(trace, "\n");
  for (j = 1; j <= n; j++) fprintf(trace,"%5d ", j);
  fprintf(trace,"\n");
  for (j = 1; j <= n; j++) fprintf(trace,"%5.2lf ", c[j]);
  fprintf(trace,"+ %5.2lf\n", v);
  for (i = 1; i <= m; i++) {
    for (j = 1; j <= n; j++) fprintf(trace,"%5.2lf ", a[i][j]);
    fprintf(trace,"= %5.2lf  basisvar %d\n", b[i], basis[i]);
    /* if (fabs(c[basis[i]]) > EPSILON) error("basis variable not zero"); */
  }
}


/* ======================================================================
				  pivot
   ====================================================================== */

void pivot(int m, int n, indices basis, matrix a, mvector b, nvector c, 
           double *v, int l, int e)
/* pivot procedure taken from Sedgewick page 618 */
{
  int i, j;

  /* compute objective value */
  *v += c[e] * b[l] / a[l][e];

  /* compute c coefficients */
  for (j = n; j >= 1; j--) {
    if (j != e) c[j] -= a[l][j] * c[e] / a[l][e];
  }
  c[e] = 0.0;

  /* compute a and b coefficients */
  for (i = 1; i <= m; i++) {
    if (i != l) {
      b[i] -= b[l] * a[i][e] / a[l][e];
      for (j = n; j >= 1; j--) {
        if (j != e) a[i][j] -= a[l][j] * a[i][e] / a[l][e];
      }
    }
  }
  b[l] /= a[l][e];

  /* update column e (entering variable) of matrix a */ 
  for (i = 1; i <= m; i++) {
    if (i != l) a[i][e] = 0.0;
  }

  /* update a row l (leaving row) of matrix a */
  for (j = 1; j <= n; j++) {
    if (j != e) a[l][j] /= a[l][e];
  }
  a[l][e] = 1.0;
  basis[l] = e;
}


/* ======================================================================
                                pivot_simplex
   ====================================================================== */

void pivot_simplex(int m, int n, indices basis, 
                   matrix a, mvector b, nvector c, double *v, int *status)
{
  double delta, mindelta;
  int i, j, e, l; 

  while (TRUE) {
    /* choose an entering column, i.e. an index e for which c[e] > 0 */
    for (j = 1; j <= n; j++) if (c[j] > EPSILON) break;
    if (j <= n) e = j; else break;  /* algorithm finished */

    /* choose a leaving variable l */
    l = 0;
    for (i = 1; i <= m; i++) {  
      if (a[i][e] > EPSILON) {
        delta = b[i] / a[i][e];
        if (l == 0) { 
          l = i; mindelta = delta; 
        } else {
          if (delta < mindelta) { l = i; mindelta = delta; }
          if (fabs(delta - mindelta) <= EPSILON) { 
            /* Blends rule for avoiding cycling */ 
            if (basis[l] < basis[e]) { l = i; mindelta = delta; }
          }
        }
      } 
    }
    if (l == 0) { *status = UNBOUNDED; *v = 0; return; }
    pivot(m, n, basis, a, b, c, v, l, e);
  }
}


/* ======================================================================
                               initialize_simplex
   ====================================================================== */

void initialize_simplex(int *m, int *n, matrix aa, mvector bb, nvector cc, 
                        indices basis, matrix a, mvector b, nvector c, 
                        double *v, int *status)
{
  int i, j, k, l;
  double bmin;

  /* initialize simplex algorithm from Cormen page 812 */
  /* convert to slack form */
  for (i = 1; i <= *m; i++) {
    for (j =    1; j <=    *n; j++) a[i][j] = aa[i][j];
    for (j = *n+1; j <= *n+*m; j++) a[i][j] = 0.0;
    a[i][*n+i] = 1.0;
  }
  for (i = 1; i <= *m; i++) b[i] = bb[i];
  for (i = 1; i <= *m; i++) basis[i] = *n + i;
  *v = 0.0;

  /* check if initial basis solution is feasible */
  for (l = 1, i = 2; i <= *m; i++) if (bb[i] < bb[l]) l = i;
  if (bb[l] > -EPSILON) {
    /* initial basis is feasible */
    for (j =    1; j <=    *n; j++) c[j] = cc[j];
    for (j = *n+1; j <= *n+*m; j++) c[j] = 0.0;
    *n = *n + *m;
  } else {
    /* define auxilary problem by adding extra variable */
    for (j =    1; j <= *n+*m; j++) c[j] = 0.0;
    c[*n+*m+1] = -1.0;
    for (i = 1; i <= *m; i++) a[i][*n+*m+1] = -1.0;
    pivot(*m, *n+*m+1, basis, a, b, c, v, l, *n+*m+1);
    pivot_simplex(*m, *n+*m+1, basis, a, b, c, v, status); 
    if (*status != FEASIBLE) return;
    /* check if extra variable has value 0 */
    for (i = 1; i <= *m; i++) {
      if ((basis[i] == *n+*m+1) && (fabs(b[i]) > EPSILON)) break;
    }
    if (i > *m) {
      /* auxilary problem is feasible with x_{m+n+1} = 0 */
      for (j =    1; j <=    *n; j++) c[j] = cc[j];
      for (j = *n+1; j <= *n+*m; j++) c[j] = 0.0;
      *n = *n + *m + 1;
      *v = 0.0;
      for (i = 1; i <= *m; i++) {
        k = basis[i]; if ((k > *n-*m) || (cc[k] == 0)) continue;
        for (j = 1; j <= *n; j++) if (j != k) c[j] -= a[i][j] * cc[k];
        *v += b[i] * cc[k];
        c[k] = 0; 
      }
      *n += 1; *m += 1;
      for (i = 1; i <= *m; i++) a[i][*n] = 0.0;
      for (j = 1; j <= *n; j++) a[*m][j] = 0.0;
      a[*m][*n] = 1.0; a[*m][*n-1] = 1.0; b[*m] = 0.0; basis[*m] = *n;
      for (i = 1; i <= *m; i++) if (basis[i] == *n-1) break;
      if (i <= *m) {
        for (j = 1; j <= *n; j++) a[*m][j] -= a[i][j];
        b[*m] -= b[i]; 
      }
    } else {
      *status = INFEASIBLE;
      *v = 0.0;
    }
  }
}


/* ======================================================================
				  simplex
   ====================================================================== */

double simplex(int m, int n, matrix aa, mvector bb, nvector cc, nvector xx,
               int *status)
{
  matrix  a;
  mvector b;
  nvector c;
  indices basis;
  double v, delta, mindelta;
  int i, j, e, l; 

  /* simplex algorithm from Cormen page 797 */
  *status = FEASIBLE;
  initialize_simplex(&m, &n, aa, bb, cc, basis, a, b, c, &v, status);
  if ((*status) != FEASIBLE) return 0;

  /* make pivot iterations until no c[i] is positive */ 
  pivot_simplex(m, n, basis, a, b, c, &v, status);
  if ((*status) != FEASIBLE) return 0;
  /* copy solution back */
  for (j = 1; j <= n; j++) xx[j] = 0.0;
  for (i = 1; i <= m; i++) xx[basis[i]] = b[i];
  return v;
}


/* ======================================================================
                                integertest
   ====================================================================== */

boolean integertest(double d)
{
  /* test if d is an integer */
  return (fabs(rint(d)-d) < EPSILON);
}


/* ======================================================================
                               branchandbound
   ====================================================================== */

void branchandbound(int m, int n, matrix a, mvector b, nvector c, 
                    nvector xstar)
{
  int status; /* status returned by simplex algorithm */
  nvector x;  /* solution vector for linear problem */
  double z;   /* upper bound */

  /* derive upper bound from lp relaxation */
  bbnodes++;
  z = simplex(m, n, a, b, c, x, &status);

  /* status is FEASIBLE if a solution was found */
  /* status is INFEASIBLE if no solution exists */
  /* status is UNBOUNDED if the solution value can be arbitrarily large */

  /* ======================================================= */
  /* HERE THE BRANCH AND BOUND ALGORITHM SHOULD BE FILLED IN */
  /* ======================================================= */
}


/* ======================================================================
				   ipsolve
   ====================================================================== */

double ipsolve(int m, int n, matrix a, mvector b, nvector c, nvector xstar)
{
  int j;

  /* check if objective coefficients are integers */
  intcost = TRUE;
  for (j = 1; j <= n; j++) if (!integertest(c[j])) intcost = FALSE;
 
  /* initialize current best (incumbent) solution */
  bbnodes = 0;
  zstar = -INFTY; 
  for (j = 0; j <= n; j++) xstar[j] = 0.0;

  /* run branch and bound */
  branchandbound(m, n, a, b, c, xstar);
  return zstar;
}


/* ======================================================================
				    main
   ====================================================================== */

int main(int argc, char *argv[])
{
  double time, z;

  if (argc != 2) { printf("filename expected as argument\n"); exit(0); }
  trace = fopen("trace.out", "w");
  readfromfile(argv[1]);
  cpu_time(NULL);
  z = ipsolve(m, n, a, b, c, x);
  cpu_time(&time);
  printinstance(m, n);
  printf("solution value %lf\n", z);
  printf("time used %.2lf\n", time);
  printf("branch and bound nodes %d\n", bbnodes);
  fprintf(trace,"solution value %lf\n", z);
  fprintf(trace,"time used %.2lf\n", time);
  fprintf(trace,"branch and bound nodes %d\n", bbnodes);
  fclose(trace);
  return 0;
}