/* ======================================================================
           SCKNAP.c, David Pisinger   jan 1996, revised apr. 1999
   ====================================================================== */

/* This code solves strongly correlated knapsack problems, which
 * may be formulated as
 *
 *           maximize   \sum_{j=1}^{n} (w_{j}+k) x_{j}
 *           subject to \sum_{j=1}^{n} w_{j} x_{j} \leq c
 *                      x_{j} \in \{0,1\}, j = 1,\ldots,n
 *
 * A description of this code is found in the following paper:
 *
 *   D.Pisinger (1998)
 *   "A fast algorithm for strongly correlated knapsack problems"
 *   Discrete Applied Mathematics, 89, 197-212 
 *
 * The present code is written in ANSI-C, and has been compiled with
 * the GNU-C compiler using option "-ansi -pedantic" as well as the
 * HP-UX C compiler using option "-Aa" (ansi standard).
 *
 * This file contains the callable routine scknap with prototype
 *
 *   int scknap(int n, int *w, int k, int c, int *x);
 *
 * the meaning of the parameters is the following:
 *   n         Size of problem, i.e. number of items to pack.
 *   w         Integer array of length n, where w[j] is
 *             the weight of box j for j=0,..,n-1. 
 *   k         An integer stating the "fixed profit", thus the
 *             profit of an item is p[j] = w[j] + k.
 *   c         Capacity of knapsack.
 *   x         Integer array of length n, where x[j] upon return
 *             is set to 1 if items should be included in the knapsack
 *             and 0 otherwise.
 * the return value of scknap is the optimal objective value z.
 *
 * (c) Copyright 1999,
 *
 *   David Pisinger                      
 *   DIKU, University of Copenhagen      
 *   Universitetsparken 1              
 *   Copenhagen, Denmark              
 *
 * This code can be used free of charge for research and academic purposes
 * only.
 */

#include <stdlib.h>
#include <stdio.h>
#include <time.h>
#include <stdarg.h>
#include <values.h>
#include <string.h>
#include <math.h>
#include <malloc.h>


/* ======================================================================
				   macros
   ====================================================================== */

#define LOGSIZE     50000      /* size of log */

#define SORTSTACK     200      /* depth of stack used for sorting */
#define MINMED        100      /* find exact median if larger size */

#define TRUE  1
#define FALSE 0

#define LEFT  1
#define RIGHT 2

#define PARTIATE  1
#define LEFTREC   2
#define RIGHTREC  3
#define SORTALL   4

#define SWAP(a,b)    {register item t;t=*(a);*(a)=*(b);*(b)=t;}
#define NO(a,p)      ((int) ((p) - (a)->fitem + 1))
#define PTR(a,w)     ((a)->h + ((w) + (a)->res))
#define WGT(a,i)     ((int)((i) - (a)->h - (a)->res))


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

typedef short        boolean;
typedef int          ntype;   /* number of stages         */
typedef short        itype;   /* item weights             */
typedef long         stype;   /* sum of profit or weight  */

/* item record */
typedef struct irec {
  itype     w;       /* weight */
  int       *x;      /* solution variable */
} item;

/* s-rec */
typedef struct {
  item      *s;      /* solution with sum x_{j} = \beta defined on [s,t] */
  item      *s1;     /* previous value of s */
  item      *p;      /* solution with sum x_{j} = \beta+1 defined on [s,t] */
} srec;

/* log record */
typedef struct {
  stype     w;       /* weight sum */
  item      *i;      /* changed item */
} lrec;

/* interval-stack */
typedef struct {
  item     *f;       /* first item in interval */
  item     *l;       /* last item in interval */
} interval;

/* all problem information */
typedef struct {
  ntype    n;                 /* number of items in problem */
  itype    k;                 /* p[i] = w[i] + k */

  item     *item0;            /* start of item array */
  item     *fitem;            /* pointer to first item */
  item     *litem;            /* pointer to last item */
  item     *t;                /* pointer to item "t" */
  item     *b;                /* pointer to break item */
  stype    res;               /* residual capacity */
  item     *beta;             /* pointer to beta item */
  itype    wbeta;             /* weight of break item */
  stype    rm;                /* upper bound on weights = wbeta+1 */
  itype    mw;                /* min weight */

  stype    c;                 /* capacity of knapsack */
  stype    z;                 /* incumbent solution */
  stype    wsumb;             /* weight sum up to break item b */
  stype    psumb;             /* profit sum up to break item b */
  stype    dantzig;           /* dantzig upper bound */
  stype    ub;                /* current upper bound */
  stype    maxwsum;           /* max weight sum of beta items */

  lrec     *fset;             /* start of log */
  lrec     *lset;             /* end of log */
  lrec     *curr;             /* current position of log */

  srec     *h;                /* table of s-values,  s[h ] = s[v] */
  srec     *h0;               /* first pos. weight,  s[h0] = s[c+1] */
  srec     *h1;               /* end of all weights, s[h1] = s[c+w{beta}+1] */

  interval *intv1, *intv2;    /* interval stacks */
  interval *intv1b, *intv2b;

  item     *fpart;            /* first item, partially sorted */
  item     *lpart;            /* last item, partially sorted */
  stype    wfpart;            /* weight sum up to first partial sorted */

  /* debug */
  long     tablesz;           /* size of table s[t], equal to res + wb */
  long     totspace;          /* size of log */
  long     coresize;          /* size of core */
  boolean  greedyok;          /* was greedy solution ok */
} allinfo;


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

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

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


/* ======================================================================
				  palloc
   ====================================================================== */

void pfree(void *p)
{
  if (p == NULL) error("freeing null");
  free(p);
}


void * palloc(long sz, long no)
{
  long *p, size;

  size = sz * no;
  if (size == 0) size = 1;
  if (size != (size_t) size) error("Alloc too big %ld", size);
  p = malloc(size);
  if (p == NULL) error("no memory size %ld", size);
  return p;
}



/* ======================================================================
                                push/pop
   ====================================================================== */

void push(allinfo *a, int side, item *f, item *l)
{
  interval *pos;
  if (l+1 == f) return;
  switch (side) {
    case LEFT : pos = a->intv1; (a->intv1)++; break;
    case RIGHT: pos = a->intv2; (a->intv2)--; break;
  }
  if (a->intv1 == a->intv2) error("interval stack full");
  pos->f = f; pos->l = l; 
}

void pop(allinfo *a, int side, item **f, item **l)
{
  interval *pos;
  switch (side) {
    case LEFT : if (a->intv1 == a->intv1b) error("pop left");
                (a->intv1)--; pos = a->intv1; break;
    case RIGHT: if (a->intv2 == a->intv2b) error("pop right");
                (a->intv2)++; pos = a->intv2; break;
  }
  *f = pos->f; *l = pos->l; 
}


/* ======================================================================
                                median
   ====================================================================== */

item *median(item *f1, item *l1, ntype s)
{
  /* Find median r of items [f1, f1+s, f1+2s, ... l1], */
  /* and ensure the ordering f1 >= r >= l1.            */
  register itype mw;
  register item *i, *j;
  item *f, *l, *k, *m, *q;
  ntype n, d;
  static item r;

  n = (l1 - f1) / s;              /* number of values      */
  f = f1;                         /* calculated first item */
  l = f1 + s * n;                 /* calculated last item  */
  k = l;                          /* saved last item       */
  q = f + s * (n / 2);            /* middle value          */

  for (;;) {
    d = (l - f + s) / s;
    m = f + s * (d / 2);

    if (d > 1) {
      if (f->w > m->w) SWAP(f, m);
      if (d > 2) {
        if (m->w > l->w) {
          SWAP(m, l);
          if (f->w > m->w) SWAP(f, m);
        }
      }
    }
    if (d <= 3) { r = *q; break; }

    r.w = mw = m->w; i = f; j = l;
    for (;;) {
      do { i += s; } while (i->w < mw);
      do { j -= s; } while (j->w > mw);
      if (i > j) break;
      SWAP(i, j);
    }

    if ((j <= q) && (q <= i)) break;
    if (i > q) l = j; else f = i;
  }
  SWAP(k, l1);
  return &r;
}


/* ======================================================================
                                partsort
   ====================================================================== */

void partsort(allinfo *a, item *f, item *l, stype ws, int what)
{
  register itype mw;
  register item *i, *j;
  register stype wi;
  item *m;
  int d;

  d = l - f + 1;
  if (d < 1) error("negative interval in partsort");
  if (d > MINMED) {
    m = median(f, l, MINMED);
  } else {
    if (d > 1) {
      m = f + d / 2;
      if (f->w > m->w) SWAP(f, m);
      if (d > 2) {
        if (m->w > l->w) {
          SWAP(m, l);
          if (f->w > m->w) SWAP(f, m);
        }
      }
    }
  }

  if (d > 3) {
    mw = m->w; i = f; j = l; wi = ws;
    for (;;) {
      do { wi += i->w; i++; } while (i->w < mw);
      do {             j--; } while (j->w > mw);
      if (i > j) break;
      SWAP(i, j);
    }

    if (wi <= a->c) {
      if (what ==  SORTALL) partsort(a, f, i-1, ws, what);
      if (what == PARTIATE) push(a, LEFT, f, i-1);
      partsort(a, i, l, wi, what);
    } else {
      if (what ==  SORTALL) partsort(a, i, l, wi, what);
      if (what == PARTIATE) push(a, RIGHT, i,  l);
      if (what ==  LEFTREC) push(a, LEFT , i,  l);
      if (what == RIGHTREC) push(a, RIGHT, i,  l);
      partsort(a, f, i-1, ws, what);
    }
  }

  if ((d <= 3) || (what == SORTALL)) {
    a->fpart  = f;
    a->lpart  = l;
    a->wfpart = ws;
  }
}


/* ======================================================================
                                  sortmore
   ====================================================================== */

void sortmore(allinfo *a, int side, item **last)
{
  item *f, *l;
  itype mw;

  switch (side) {
    case LEFT : pop(a, LEFT, &f, &l);
                partsort(a, f, l, a->c+1, LEFTREC);
                *last = a->lpart;
                if (*last >= a->b) error("bad last");
                break;
    case RIGHT: pop(a,RIGHT, &f, &l);
                partsort(a, f, l, a->c+1, RIGHTREC);
                *last = a->lpart;
                break;
  }
}


/* ======================================================================
                                  sortrest
   ====================================================================== */

void sortrest(allinfo *a)
{
  item *f, *l;
  itype mw;

  while (a->intv1 != a->intv1b) { 
    pop(a, LEFT, &f, &l); 
    partsort(a, f, l, 0, SORTALL);
  }
  while (a->intv2 != a->intv2b) {
    pop(a, RIGHT, &f, &l); 
    partsort(a, f, l, 0, SORTALL);
  }
}


/* ======================================================================
                                  maxweight
   ====================================================================== */

stype maxweight(allinfo *a, ntype beta)
{
  register item *i, *m;
  register stype wsum;

  wsum = 0;
  for (i = a->litem-beta+1, m = a->litem+1; i != m; i++) wsum += i->w;
  if (wsum > a->c) return a->c;
  if (wsum + beta * a->k > a->z) {
    a->z = wsum + beta * a->k;
    for (i = a->fitem, m = a->litem-beta+1; i != m; i++) *(i->x) = 0;
    for (i = a->litem-beta+1, m = a->litem+1; i != m; i++) *(i->x) = 1;
  }
  return wsum;
}



/* ======================================================================
                                  heuristic
   ====================================================================== */

boolean heuristic(allinfo *a)
{
  register item *i, *j, *m, *b;
  register itype r, d;
  item *m1, *m2, *f, *l;  /* m1: first sorted item, m2: last sorted item */

  r = a->c - a->wsumb;
  b = a->b;

  /* partition around residual gap */
  d = b->w - r;
  for (i = a->fitem, j = a->b-1; i <= j; ) {
    *(i->x) = 1;
    if (i->w < d) { i++; } else { SWAP(i,j); j--; }
  }
  if (r == 0) return TRUE;
  if (i == b) return FALSE;

  /* prepare for partitioning */
  a->intv1 = a->intv1b;
  push(a, LEFT, a->fitem, i-1);
  push(a, LEFT, i, b-1);

  /* now run through items */
  m2 = a->lpart; /* Important: no other sorting before this */
  sortmore(a, LEFT, &m1); 
  for (j = b, m = a->litem + 1; j != m; j++) {
    if (j > m2) sortmore(a, RIGHT, &m2);
    d = j->w - r;
    while (i->w < d) {
      i++;
      if (i > m1) { 
        if (i >= b) return FALSE;
        sortmore(a, LEFT, &m1); 
      }
    }
    if (i->w == d) {
      *(i->x) = 0; *(j->x) = 1; return TRUE;
    }
  }
  return FALSE;
}


/* ======================================================================
				  definesolution
   ====================================================================== */

void definesolution(allinfo *a)
{
  register lrec *i, *m;
  register srec *j, *k;
  register item *s, *f, *t;
  register stype w, wm;
  stype z;

  /* find best solution and store info */
  for (j = a->h0-1, w = 0, f = a->item0; j->s == f; j--) w--;
  z = a->c + w + a->k * NO(a,a->beta);
  a->coresize = (a->t - a->b);
  if (a->res + a->wbeta > a->tablesz ) a->tablesz  = a->res + a->wbeta;
  if (a->curr - a->fset > a->totspace) a->totspace = a->curr - a->fset;

  /* check if optimal */
  if (z <= a->z) return;

  /* improved solution is found */
  for (s = a->fitem, t = a->beta+1; s != t; s++) *(s->x) = 1;
  for (s = a->beta+1, t = a->litem+1; s != t; s++) *(s->x) = 0;
  a->z = z;

  /* backtrack solution */
  wm = a->rm;
  for (i = a->curr-1, m = a->fset-1, t = NULL; i != m; i--) {
    if (i->w == w) {
      s = i->i;
      *(s->x) = 0; w += s->w; /* remove items s */
      *(t->x) = 1; w -= t->w; /* insert item t */
      /* avoid more choices from this group */
      do { i--; } while (i->w != wm);
    }
    if (i->w == wm) t = i->i;
  }
}


/* ======================================================================
				  iterate
   ====================================================================== */

void iterate(allinfo *a, item *t)
{
  register item *j, *n;
  register srec *i, *k, *l, *m;
  register stype w;
  srec *h, *h0;

  /* add weight t to states with beta items, remove from states with beta+1 */
  h0 = a->h + a->res;
  for (i = a->h, k = i + t->w, m = a->h0 + a->wbeta - t->w, w = WGT(a,k); 
       i != m; i++, k++, w++) {
    if (i->s1 > k->p) {
      n = k->p - 1; h = h0 + w; k->p = i->s1; 
      for (j = k->p-1; j != n; j--) {
        l = h - j->w;
        if (j > l->s) l->s = j; 
      }
    }
  }
}


/* ======================================================================
				  cleanup
   ====================================================================== */

void cleanup(allinfo *a)
{
  register lrec *k, *l;
  register srec *i, *m;
  register stype w;

  k = a->curr; l = a->lset;

  /* negative weights are stacked and rotated */
  for (i = a->h, m = a->h0, w = -a->res; i != m; i++, w++) {
    if (i->s != i->s1) {
      k->i  = i->s;
      k->w  = w;
      k++; if (k == l) break;
      i->s1 = i->s;
    }
  }
  if (k == l) error("overfull log");

  /* positive weights are kept unchanged */

  /* save current value of t */
  k->w = a->rm;  /* to not match any weights */
  k->i = a->t;
  k++; if (k == l) error("overfull log");
  a->curr = k;
}


/* ======================================================================
				  initfirst
   ====================================================================== */

void initfirst(allinfo *a)
{
  register srec *i, *m;
  register item *f, *f0, *j;
  register stype w;
  lrec *k;
  int size;

  /* initialize vector table */
  a->fset = palloc(LOGSIZE, sizeof(lrec));
  a->curr = a->fset;
  a->lset = a->fset + LOGSIZE;
  a->rm   = a->wbeta + 1;

  /* initialize hashtable */
  size = a->res + a->wbeta + 1;
  a->h = palloc(size, sizeof(srec));
  a->h0 = a->h + a->res + 1;
  a->h1 = a->h + size;
  for (i = a->h, m = a->h0, f = a->item0, f0 = a->fitem; i != m; i++) {
    i->s  = f;
    i->s1 = f;
    i->p  = f0;
  }
  for (i = a->h0, m = a->h1, w = 1; i != m; i++, w++) {
    while (f0->w < w) f0++; /* halts as w <= w{b} */
    i->s  = f;
    i->s1 = f;
    i->p  = f0;
  }

  /* set weight of item 0 to 0 */
  a->item0->w = 0;

  /* insert one vector */
  i = PTR(a, -a->res);
  i->s  = a->beta+1;
  i->s1 = a->beta+1;

  k = a->curr;
  k->w = a->rm;
  k->i = a->beta;
  a->curr++;
}


/* ======================================================================
                                  copyproblem
   ====================================================================== */

void copyproblem(item *f, item *l, int *w, int *x)
{
  register item *i, *m;
  register int *ww, *xx;

  for (i = f, m = l+1, ww = w, xx = x; i != m; i++, ww++, xx++) {
    i->w = *ww; i->x = xx;
  }
}


/* ======================================================================
				findbreak
   ====================================================================== */

void findbreak(allinfo *a)
{
  register item *i, *m;
  stype wsum, c;

  if (a->c < 0) error("negative capacity");
  wsum = a->wfpart; c = a->c; /* notice: x[i] is set to 1 in heuristic */
  for (i = a->fpart; wsum <= c; i++) wsum += i->w; 
  i--;  /* we went one item too far */
  wsum -= i->w;
  if (i > a->lpart) error("limits blown");

  a->item0   = a->fitem - 1;
  a->b       = i;
  a->beta    = i;
  a->wbeta   = i->w;
  a->wsumb   = wsum;
  a->psumb   = wsum + a->k * NO(a,i-1);
  a->z       = a->wsumb + a->k * NO(a,i-1);
  a->dantzig = a->psumb + ((i->w + a->k) * (double)(a->c - a->wsumb)) / i->w;
  a->ub      = a->c + a->k * NO(a,i-1);

  /* initialize remaining vector */
  for (m = a->litem + 1; i != m; i++) *(i->x) = 0;
}


/* ======================================================================
				subweights
   ====================================================================== */

void subweights(allinfo *a)
{
  register item *i, *m;
  register itype mw;

  mw = a->fitem->w - 1;
  for (i = a->fitem, m = a->litem+1; i != m; i++) i->w -= mw;
  a->mw = mw;
}


/* ======================================================================
				addweights
   ====================================================================== */

void addweights(allinfo *a)
{
  register item *i, *m;
  register itype mw;

  mw = a->mw;
  for (i = a->fitem, m = a->litem+1; i != m; i++) i->w += mw;
}


/* ======================================================================
				scknap
   ====================================================================== */

extern int scknap(int n, int *w, int k, int c, int *x)
{
  allinfo a;
  item *tab;
  interval *inttab;
  ntype beta;

  /* allocate space for test example and two border items */
  tab = (item *) palloc(sizeof(item), n+1);
  a.fitem = &tab[1]; a.litem = &tab[n];
  copyproblem(a.fitem, a.litem, w, x);
  a.n        = n;
  a.c        = c;
  a.k        = k;
  a.greedyok = FALSE;
  a.tablesz  = 0;
  a.totspace = 0;
  a.coresize = 0;

  inttab  = palloc(SORTSTACK, sizeof(interval));
  a.intv1 = a.intv1b = &inttab[0];
  a.intv2 = a.intv2b = &inttab[SORTSTACK - 1];
  partsort(&a, a.fitem, a.litem, 0, PARTIATE);
  findbreak(&a);

  /* solve greedy */
  a.greedyok = heuristic(&a);

  if (a.greedyok) {
    a.z = a.c + a.k * NO(&a,a.b-1);
  } else {
    sortrest(&a);
    a.res = a.c - a.wsumb;
    for (beta = NO(&a,a.b)-1; beta != 0; beta--, a.res += a.wbeta) {
      a.beta    = a.fitem + beta - 1;
      a.wbeta   = a.beta->w;
      a.maxwsum = maxweight(&a, beta);
      a.ub      = a.maxwsum + a.k * beta;
      if ((a.beta+1)->w-a.beta->w > a.res) continue; 
      if (a.ub <= a.z) break; 

      initfirst(&a);

      for (a.t = a.fitem + beta; a.t <= a.litem; a.t++) {
        if (a.t->w > a.res + a.wbeta) break; 
        iterate(&a, a.t);
        cleanup(&a);
        if ((a.h0-1)->s != a.item0) break; 
      }
      definesolution(&a);
      pfree(a.h);
      pfree(a.fset);
    }
  }
  pfree(tab);
  return a.z;
}