/* ======================================================================
                          PROJEKT UDVAELGELSE  
             projektopgave 2 paa Videregaaende Algoritmik 2006/07
  			    David Pisinger
   ====================================================================== */

/* Dette C program loeser projektudvaelgelsesproblemet
 *
 *    maximer      sum_{i=1,2,..,n} p_i x_i
 *    saaledes at  sum_{i=1,2,..,n} r_i x_i <= d
 *                 x_i <= x_j          , (i,j) tilhoerer E
 *                 x_i tilhoerer {0,1} , i tilhoerer V
 * hvor 
 *   E = (V,E) er en kravgraf.
 *   p_1,...,p_n er profitter af projekterne (positive/negative heltal)
 *   w_1,...,w_n er ressourceforbrug af projekterne (positive heltal)
 *   d er en oevre graense paa ressourceforbruget
 *
 * programmet oversaettes paa DIKU's linux pc'er med 
 *
 *   gcc -o projekt -O3 projekt.c
 *
 * Den udfoerbare kode kaldes med
 * 
 *   projekt instans
 *
 * hvor "instans" er et filnavn
 */

#define MSIZE 200             /* maximalt antal projekter */
#define PATH  "."             /* sti til instanser */
#define INFTY 9999            /* uendelig */
#define SCALE 10              /* skalering af flydende tal til heltal */

 
/* ======================================================================
		         inkludering af biblioteksfiler
   ====================================================================== */

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


/* ======================================================================
				 type erklaeringer
   ====================================================================== */

#define MS    (MSIZE+2)  /* foroeg med 2 saa plads til s og t  */
#define TRUE  1          /* logisk sand                        */
#define FALSE 0          /* logisk falsk                       */
#define FREE  -1         /* variabel er fri i branch-and-bound */

typedef int     boolean; /* boolsk variabel                    */

typedef int edge[MS][MS];
typedef int profit[MS];
typedef int resource[MS];
typedef int capacity[MS][MS];
typedef int flow[MS][MS];
typedef int list[MS];


/* ======================================================================
 		               globale variable
   ====================================================================== */

char    text[100];       /* tekst svarende til instans                   */
list    xstar;           /* (hidtil) bedste loesning                     */
int     zstar;           /* hidtil bedste loesningsvaerdi                */
int     bbnodes;         /* antal knuder besoegt i branch-and-bound      */
int     ubroot;          /* oevre graensevaerdi i rodknuden              */
double  lambdaroot;      /* bedste vaerdi af lambda i rodknuden          */
double  tottime;         /* cpu-tid brugt i algoritmen                   */


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

/* Den foelgende rutine bruges til maaling af brugt cpu-tid.
 *   cpu_time(NULL)    - starter stopuret
 *   cpu_time(&t)      - returnerer brugte cpu-tid i variablen 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);
  }
}


/* =======================================================================
	              	  udskriv fejl og afslut
   ======================================================================= */

/* Udskriv fejlmaelding og terminer hvis en fejl opdages */

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

  va_start(args, str);
  vprintf(str, args);
  printf("\n");
  va_end(args);
  exit(-1);
} 


/* ======================================================================
				show_instance
   ====================================================================== */

/* Vis nuvaerende instans paa skaermen */
void show_instance(int n, int d, profit p, resource r, edge e)
{
  int i, j;

  printf("n %d, d %d\n\n", n, d);

  printf("(i,p,r)\n");
  for (i = 1; i <= n; i++) { 
    printf("(%3d,%3d,%3d)\n", i, p[i], r[i]); 
  }
  printf("\n");

  printf("(i,j)\n");
  for (i = 1; i <= n; i++) {
    for (j = 1; j <= n; j++) {
      if (e[i][j]) printf("(%3d,%3d)\n", i, j);
    }
  }
  printf("\n");
}


/* ======================================================================
				show_flow
   ====================================================================== */

/* Vis kant-kapaciteter og flow paa skaermen */
void show_flow(int n, capacity c, flow f)
{
  int i, j;

  printf("(i,j) capacity flow\n");
  for (i = 0; i <= n+1; i++) {
    for (j = 0; j <= n+1; j++) {
      if (c[i][j] > 0) {
        printf("(%3d,%3d) %4d %4d\n", i, j, c[i][j], f[i][j]);
      }
    }
  }
}


/* ======================================================================
				show_path
   ====================================================================== */

/* Vis fundne augmenting path */
void show_path(int n, list path)
{
  int i, t;

  printf("path =");
  t = n+1;
  for (i = 0; i <= n; i++) {
    printf(" %d", path[i]);
    if (path[i] == t) break;
  }
  printf("\n");
}


/* ======================================================================
				show_cut
   ====================================================================== */

/* Vis maengden S i snittet (S,T) */
void show_cut(int n, list visited)
{
  int i;

  printf("cut S = {");
  for (i = 0; i <= n+1; i++) {
    if (visited[i]) printf(" %d", i);
  }
  printf("}\n");
}


/* ======================================================================
				read_instance
   ====================================================================== */

/* Laes en instans til tabellerne p, r, og e.
 * Instansen laeses fra filen "navn" med soegesti svarende til 
 * makroen PATH defineret foerst i programmet. Sidst i proceduren
 * kontrolleres at inddata overholder de givne antagelser. 
 */

void read_instance(char *name, int *n, int *d, 
		   profit p, resource r, edge e)
{
  int i, j, t, u, v, m;
  char s[100];
  FILE *in;
  char c;

  /* konstruer filnavn og aaben fil */
  sprintf(s, "%s/%s", PATH, name);
  in = fopen(s, "r"); if (in == NULL) error("%s no such file", s);

  /* skip en eventuel kommentar i foerste linie */
  c = getc(in); strcpy(text, "unknown instance id\n");
  if (c == '#') { fgets(text,100,in); printf("#%s",text); } else ungetc(c,in);

  fscanf(in, "%d", &t); *n = t; 
  fscanf(in, "%d", &t); m = t; 
  fscanf(in, "%d", &t); *d = t; 

  /* initialiser kant matrix */
  for (i = 1; i <= *n; i++) {
    for (j = 1; j <= *n; j++) e[i][j] = 0;
  }

  /* laes data */
  for (i = 1; i <= *n; i++) { fscanf(in, "%d", &t); p[i] = t; }
  for (i = 1; i <= *n; i++) { fscanf(in, "%d", &t); r[i] = t; }
  for (i = 1; i <= m; i++) { fscanf(in, "%d,%d", &u, &v); e[u][v] = 1; }

  /* kontroller at instansen overholder antagne egenskaber */
  if (*n <= 0) error("empty instance");
  for (i = 1; i <= *n; i++) {
    if (r[i] < 0) error("bad resource\n");
  }
  fclose(in);
}


/* ======================================================================
				find_cut
   ====================================================================== */

/* returner kapaciteten af snittet (S,T), hvor S er givet ved tabel visited */
int find_cut(int n, list visited, capacity c)
{
  int i, j;
  int cap;

  cap = 0;
  for (i = 0; i <= n+1; i++) {
    for (j = 0; j <= n+1; j++) {
      if (visited[i] && !visited[j]) cap += c[i][j];
    }
  }
  return cap;
}


/* ======================================================================
                                increase_flow
   ====================================================================== */

/* Foroeg stroemning f med vaerdien cr langs stien path */
void increase_flow(int n, list path, flow f, capacity c, int cr)
{
  int i, u, v;
  int e;
 
  for (i = 0; path[i] != n+1; i++) {
    u = path[i]; v = path[i+1]; e = cr;
    if (f[v][u] > e) { f[v][u] -= e; continue; }
    if (f[v][u] > 0) { e -= f[v][u]; f[v][u] = 0; }
    f[u][v] += e;    
    if (f[u][v] > c[u][v]) error("flow > cap");
    if (f[v][u] > c[v][u]) error("flow > cap");
  }
}


/* ======================================================================
                                current_flow
   ====================================================================== */

/* Find vaerdien af nuvaerende stroemning f */
int current_flow(int n, flow f)
{
  int i;
  int fval;

  fval = 0;
  for (i = 1; i <= n+1; i++) fval += f[0][i];
  return fval;
}


/* ======================================================================
                                augmenting_path
   ====================================================================== */

/* Find en augmenting path i grafen givet ved kapaciteter c og flow f */
/* Returner visited[i]=TRUE hvis en knude er blevet markeret */
/* Returner path som angiver en sti fra s til t hvis denne findes */

boolean augmenting_path_recurs(int k, int i, int n, 
		               list visited, list path, capacity c, flow f)
{
  int j, s, t;
 
  /* dybde foerst soegning efter augmenting path startende fra knude i */
  s = 0; t = n+1;
  visited[i] = TRUE;
  if (i == t) return TRUE;
  for (j = 1; j <= t; j++) { 
    if (!visited[j]) {
      if (c[i][j] - f[i][j] + f[j][i] > 0) {
        if (augmenting_path_recurs(k+1, j, n, visited, path, c, f)) { 
          path[k] = j; return TRUE; 
	}
      }
    }
  }
  return FALSE;
}

int augmenting_path(int n, list visited, list path, capacity c, flow f)
{
  int i, u, v, s, t;
  int cr, ce;
  boolean found;

  /* initialiser tabeller path, visited */
  s = 0; t = n+1;
  for (i = 0; i <= n+1; i++) visited[i] = FALSE;
  visited[s] = TRUE;
  path[0] = s; path[1] = t;

  /* kald dybde-foerst soegning */
  found = augmenting_path_recurs(1, 0, n, visited, path, c, f);
  if (!found) return 0;

  /* find kapacitet af agumenting path */
  cr = c[path[0]][path[1]];
  for (i = 0; path[i] != t; i++) {
    u = path[i]; v = path[i+1];
    ce = c[u][v] - f[u][v] + f[v][u];
    if (ce < cr) cr = ce;
  }
  return cr;
}


/* ======================================================================
                                maximum_flow
   ====================================================================== */

/* Find en maksimal stroemning i grafen givet ved kantkapaciteter c */
/* Det antages at s=0, mens t=n+1. Oevrige knuder er 1, 2, ..., n */
/* Vaerdien af stroemningen returneres */

int maximum_flow(int n, capacity c)
{
  int i, j;
  int cr, cap, fval;
  flow f;
  list visited;
  list path;

  /* check kapaciteter */
  for (i = 0; i <= n+1; i++) {
    for (j = 0; j <= n+1; j++) if (c[i][j] < 0) error("negative c");
  }

  /* initialiser stroemning f til nul */
  for (i = 0; i <= n+1; i++) {
    for (j = 0; j <= n+1; j++) f[i][j] = 0;
  }

  /* find augmenting path og foroeg stroemning */
  for (;;) {
    cr = augmenting_path(n, visited, path, c, f);
    if (cr == 0) break;
    increase_flow(n, path, f, c, cr);
  }

  /* tabel visited[] definerer minimalt snit (S,T) */
  cap = find_cut(n, visited, c);

  /* find vaerdi af stroemning */
  fval = current_flow(n, f); 
  if (cap != fval) error("max flow %d != min cut %d\n", fval, cap); 

  return fval;
}


/* ======================================================================
                                solve_relaxed
   ====================================================================== */

/* Relaxer project selection problem med lambda, og loes tilhoerende */
/* uncapacitated project selection problem. */
/* Da de fremkomne kapaciteter er reelle tal skaleres de med faktoren */
/* SCALE. Derfor er precisionen af lambda angivet ved 1/SCALE */

int solve_relaxed(int n, int d, double lambda, profit p, resource r, edge e)
{
  int i, j;
  int q, fval, qpos;
  int lambdai;
  capacity c;

  /* skaler og afrund lambda til heltallig vaerdi */
  lambdai = lambda * SCALE;     

  /* relaxer problemet med heltallig lambda */
  for (i = 0; i <= n+1; i++) { 
    c[0][i] = c[n+1][i] = c[i][0] = c[i][n+1] = 0; 
  }
  for (i = 1; i <= n; i++) {
    for (j = 1; j <= n; j++) {
      c[i][j] = INFTY * e[i][j];
    }
  }

  qpos = 0;
  for (j = 1; j <= n; j++) {
    q = SCALE * p[j] - lambdai * r[j];
    if (q >= 0) { c[0][j] = q; qpos += q; } else { c[j][n+1] = -q; }
  }

  /* loes maximum flow problemet */
  fval = maximum_flow(n, c);
  return (qpos-fval+lambdai*d)/SCALE;
}


/* ======================================================================
                                upper_bound
   ====================================================================== */

/* find oevre graensevaerdi */
int upper_bound(int n, int d, profit p, resource r, edge e)
{
}



/* ======================================================================
                                improved
   ====================================================================== */

/* kontroller om loesning x er bedre end hidtil bedste loesning */
void improved(int n, list x, profit p)
{
  int i, psum;

  /* find profit sum */
  psum = 0;
  for (i = 1; i <= n; i++) { if (x[i]) psum += p[i]; }
  if (psum <= zstar) return;

  /* save improved solution */
  zstar = psum;
  for (i = 1; i <= n; i++) xstar[i] = x[i];
}


/* ======================================================================
                                infeasible
   ====================================================================== */

/* kontroller om vektoren x[] overskrider kapaciteten af ressourcen */
boolean infeasible(int n, int d, list x, resource r)
{
  int i;

  for (i = 1; i <= n; i++) { if (x[i] == TRUE) d -= r[i]; }
  return (d < 0);
}


/* ======================================================================
                                allfixed
   ====================================================================== */

/* kontroller om alle variable i x[] har faaet tildelt en vaerdi */
boolean allfixed(int n, list x)
{
  int i;

  for (i = 1; i <= n; i++) { if (x[i] == FREE) return FALSE; }
  return TRUE;
}


/* ======================================================================
                                branch-and-bound
   ====================================================================== */


void branch_and_bound(int v, int n, int d, list x, profit p, resource r, edge e)
{
  int i, j, m, bno;
  int pfix, rfix, u;
  list x1; /* det kan vaere hensigtsmaessigt at gemme en kopi af x */

  /* check om kapacitet overskrides */
  bbnodes++; 
  if (infeasible(n, d, x, r)) { 
    return; 
  }
 
  /* check om alle x[] er blevet tildelt en vaerdi */
  if (allfixed(n, x)) {
    improved(n, x, p); 
    return; 
  }

  /* udregn oevre graensevaerdi. Det kan vaere hensigtsmaessigt at */
  /* kopiere tabellerne til nye matricer svarende til de frie variable */

  /* upper bound test */

  /* branching */

  /* retablering af x-vektor */
}


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

/* Hovedprogram. Parametre til algoritmen indlaeses, instansen laeses fra
 * den angivne fil, og det oenskede problem loeses.
 * Tidsforbrug og antal knuder i branch-and-bound algoritmen udskrives.
 */

int main(int argc, char *argv[])
{
  char name[100];
  int n;
  int d;
  profit p;
  resource r;
  edge e;
  list x;
  int i, rsum;

  if (argc == 2) {
    strcpy(name, argv[1]); 
    printf("project selection '%s'\n", name);
  } else {
    printf("project selection\n");
    printf("filename: ");
    scanf("%s", name);
  }

  /* laes instansen fra fil med navn 'name' */
  read_instance(name, &n, &d, p, r, e);
  show_instance(n, d, p, r, e); 

  /* initialiser hidtil bedste loesning og kald branch-and-bound */
  cpu_time(NULL);
  zstar = 0; bbnodes = 0; 
  for (i = 1; i <= n; i++) x[i] = FREE;
  branch_and_bound(0, n, d, x, p, r, e);
  cpu_time(&tottime);

  /* find ressourceforbrug af optimale loesning */
  for (i = 1; i <= n; i++) if (xstar[i] == TRUE) rsum += r[i];

  printf("size: %d\n", n);
  printf("nodes in b&b tree: %d\n", bbnodes);
  printf("upper bound in root node: %d\n", ubroot);
  printf("lambda in root node: %.2lf\n", lambdaroot);
  printf("optimal solution: %d\n", zstar);
  printf("resources available: %d\n", d);
  printf("resources used: %d\n", rsum);
  printf("total cpu time (seconds): %.2lf\n\n", tottime);
  return 0;
}