//===---------------- PBQP.cpp --------- PBQP Solver ------------*- C++ -*-===//
//
//                     The LLVM Compiler Infrastructure
//
// This file is distributed under the University of Illinois Open Source
// License. See LICENSE.TXT for details.
//
//===----------------------------------------------------------------------===//
//
// Developed by:                   Bernhard Scholz
//                             The University of Sydney
//                         http://www.it.usyd.edu.au/~scholz
//===----------------------------------------------------------------------===//

#include "PBQP.h"
#include "llvm/Config/alloca.h"
#include <limits>
#include <cassert>
#include <cstring>

namespace llvm {

/**************************************************************************
 * Data Structures 
 **************************************************************************/

/* edge of PBQP graph */
typedef struct adjnode {
  struct adjnode *prev,      /* doubly chained list */ 
                 *succ, 
                 *reverse;   /* reverse edge */
  int adj;                   /* adj. node */
  PBQPMatrix *costs;         /* cost matrix of edge */

  bool tc_valid;              /* flag whether following fields are valid */
  int *tc_safe_regs;          /* safe registers */
  int tc_impact;              /* impact */ 
} adjnode;

/* bucket node */
typedef struct bucketnode {
  struct bucketnode *prev;   /* doubly chained list */
  struct bucketnode *succ;   
  int u;                     /* node */
} bucketnode;

/* data structure of partitioned boolean quadratic problem */
struct pbqp {
  int num_nodes;             /* number of nodes */
  int max_deg;               /* maximal degree of a node */
  bool solved;               /* flag that indicates whether PBQP has been solved yet */
  bool optimal;              /* flag that indicates whether PBQP is optimal */
  PBQPNum min;
  bool changed;              /* flag whether graph has changed in simplification */

                             /* node fields */
  PBQPVector **node_costs;   /* cost vectors of nodes */
  int *node_deg;             /* node degree of nodes */
  int *solution;             /* solution for node */
  adjnode **adj_list;        /* adj. list */
  bucketnode **bucket_ptr;   /* bucket pointer of a node */

                             /* node stack */
  int *stack;                /* stack of nodes */
  int stack_ptr;             /* stack pointer */

                             /* bucket fields */
  bucketnode **bucket_list;  /* bucket list */

  int num_r0;                /* counters for number statistics */
  int num_ri;
  int num_rii;
  int num_rn; 
  int num_rn_special;      
};

bool isInf(PBQPNum n) { return n == std::numeric_limits<PBQPNum>::infinity(); } 

/*****************************************************************************
 * allocation/de-allocation of pbqp problem 
 ****************************************************************************/

/* allocate new partitioned boolean quadratic program problem */
pbqp *alloc_pbqp(int num_nodes)
{
  pbqp *this_;
  int u;
  
  assert(num_nodes > 0);
  
  /* allocate memory for pbqp data structure */   
  this_ = (pbqp *)malloc(sizeof(pbqp));

  /* Initialize pbqp fields */
  this_->num_nodes = num_nodes;
  this_->solved = false;
  this_->optimal = true;
  this_->min = 0.0;
  this_->max_deg = 0;
  this_->changed = false;
  this_->num_r0 = 0;
  this_->num_ri = 0;
  this_->num_rii = 0;
  this_->num_rn = 0;
  this_->num_rn_special = 0;
  
  /* initialize/allocate stack fields of pbqp */ 
  this_->stack = (int *) malloc(sizeof(int)*num_nodes);
  this_->stack_ptr = 0;
  
  /* initialize/allocate node fields of pbqp */
  this_->adj_list = (adjnode **) malloc(sizeof(adjnode *)*num_nodes);
  this_->node_deg = (int *) malloc(sizeof(int)*num_nodes);
  this_->solution = (int *) malloc(sizeof(int)*num_nodes);
  this_->bucket_ptr = (bucketnode **) malloc(sizeof(bucketnode **)*num_nodes);
  this_->node_costs = (PBQPVector**) malloc(sizeof(PBQPVector*) * num_nodes);
  for(u=0;u<num_nodes;u++) {
    this_->solution[u]=-1;
    this_->adj_list[u]=NULL;
    this_->node_deg[u]=0;
    this_->bucket_ptr[u]=NULL;
    this_->node_costs[u]=NULL;
  }
  
  /* initialize bucket list */
  this_->bucket_list = NULL;
  
  return this_;
}

/* free pbqp problem */
void free_pbqp(pbqp *this_)
{
  int u;
  int deg;
  adjnode *adj_ptr,*adj_next;
  bucketnode *bucket,*bucket_next;
  
  assert(this_ != NULL);
  
  /* free node cost fields */
  for(u=0;u < this_->num_nodes;u++) {
    delete this_->node_costs[u];
  }
  free(this_->node_costs);
  
  /* free bucket list */
  for(deg=0;deg<=this_->max_deg;deg++) {
    for(bucket=this_->bucket_list[deg];bucket!=NULL;bucket=bucket_next) {
      this_->bucket_ptr[bucket->u] = NULL;
      bucket_next = bucket-> succ;
      free(bucket);
    }
  }
  free(this_->bucket_list);
  
  /* free adj. list */
  assert(this_->adj_list != NULL);
  for(u=0;u < this_->num_nodes; u++) {
    for(adj_ptr = this_->adj_list[u]; adj_ptr != NULL; adj_ptr = adj_next) {
      adj_next = adj_ptr -> succ;
      if (u < adj_ptr->adj) {
        assert(adj_ptr != NULL);
        delete adj_ptr->costs;
      }
      if (adj_ptr -> tc_safe_regs != NULL) {
           free(adj_ptr -> tc_safe_regs);
      }
      free(adj_ptr);
    }
  }
  free(this_->adj_list);
  
  /* free other node fields */
  free(this_->node_deg);
  free(this_->solution);
  free(this_->bucket_ptr);

  /* free stack */
  free(this_->stack);

  /* free pbqp data structure itself */
  free(this_);
}


/****************************************************************************
 * adj. node routines 
 ****************************************************************************/

/* find data structure of adj. node of a given node */
static
adjnode *find_adjnode(pbqp *this_,int u,int v)
{
  adjnode *adj_ptr;
  
  assert (this_ != NULL);
  assert (u >= 0 && u < this_->num_nodes);
  assert (v >= 0 && v < this_->num_nodes);
  assert(this_->adj_list != NULL);

  for(adj_ptr = this_ -> adj_list[u];adj_ptr != NULL; adj_ptr = adj_ptr -> succ) {
    if (adj_ptr->adj == v) {
      return adj_ptr;
    }
  }
  return NULL;
}

/* allocate a new data structure for adj. node */
static
adjnode *alloc_adjnode(pbqp *this_,int u, PBQPMatrix *costs)
{
  adjnode *p;

  assert(this_ != NULL);
  assert(costs != NULL);
  assert(u >= 0 && u < this_->num_nodes);

  p = (adjnode *)malloc(sizeof(adjnode));
  assert(p != NULL);
  
  p->adj = u;
  p->costs = costs;  

  p->tc_valid= false;
  p->tc_safe_regs = NULL;
  p->tc_impact = 0;

  return p;
}

/* insert adjacence node to adj. list */
static
void insert_adjnode(pbqp *this_, int u, adjnode *adj_ptr)
{

  assert(this_ != NULL);
  assert(adj_ptr != NULL);
  assert(u >= 0 && u < this_->num_nodes);

  /* if adjacency list of node is not empty -> update
     first node of the list */
  if (this_ -> adj_list[u] != NULL) {
    assert(this_->adj_list[u]->prev == NULL);
    this_->adj_list[u] -> prev = adj_ptr;
  }

  /* update doubly chained list pointers of pointers */
  adj_ptr -> succ = this_->adj_list[u];
  adj_ptr -> prev = NULL;

  /* update adjacency list pointer of node u */
  this_->adj_list[u] = adj_ptr;
}

/* remove entry in an adj. list */
static
void remove_adjnode(pbqp *this_, int u, adjnode *adj_ptr)
{
  assert(this_!= NULL);
  assert(u >= 0 && u <= this_->num_nodes);
  assert(this_->adj_list != NULL);
  assert(adj_ptr != NULL);
  
  if (adj_ptr -> prev == NULL) {
    this_->adj_list[u] = adj_ptr -> succ;
  } else {
    adj_ptr -> prev -> succ = adj_ptr -> succ;
  } 

  if (adj_ptr -> succ != NULL) {
    adj_ptr -> succ -> prev = adj_ptr -> prev;
  }

  if(adj_ptr->reverse != NULL) {
    adjnode *rev = adj_ptr->reverse;
    rev->reverse = NULL;
  }

  if (adj_ptr -> tc_safe_regs != NULL) {
     free(adj_ptr -> tc_safe_regs);
  }

  free(adj_ptr);
}

/*****************************************************************************
 * node functions 
 ****************************************************************************/

/* get degree of a node */
static
int get_deg(pbqp *this_,int u)
{
  adjnode *adj_ptr;
  int deg = 0;
  
  assert(this_ != NULL);
  assert(u >= 0 && u < this_->num_nodes);
  assert(this_->adj_list != NULL);

  for(adj_ptr = this_ -> adj_list[u];adj_ptr != NULL; adj_ptr = adj_ptr -> succ) {
    deg ++;
  }
  return deg;
}

/* reinsert node */
static
void reinsert_node(pbqp *this_,int u)
{
  adjnode *adj_u,
          *adj_v;

  assert(this_!= NULL);
  assert(u >= 0 && u <= this_->num_nodes);
  assert(this_->adj_list != NULL);

  for(adj_u = this_ -> adj_list[u]; adj_u != NULL; adj_u = adj_u -> succ) {
    int v = adj_u -> adj;
    adj_v = alloc_adjnode(this_,u,adj_u->costs);
    insert_adjnode(this_,v,adj_v);
  }
}

/* remove node */
static
void remove_node(pbqp *this_,int u)
{
  adjnode *adj_ptr;

  assert(this_!= NULL);
  assert(u >= 0 && u <= this_->num_nodes);
  assert(this_->adj_list != NULL);

  for(adj_ptr = this_ -> adj_list[u]; adj_ptr != NULL; adj_ptr = adj_ptr -> succ) {
    remove_adjnode(this_,adj_ptr->adj,adj_ptr -> reverse);
  }
}

/*****************************************************************************
 * edge functions
 ****************************************************************************/

/* insert edge to graph */
/* (does not check whether edge exists in graph */
static
void insert_edge(pbqp *this_, int u, int v, PBQPMatrix *costs)
{
  adjnode *adj_u,
          *adj_v;
  
  /* create adjanceny entry for u */
  adj_u = alloc_adjnode(this_,v,costs);
  insert_adjnode(this_,u,adj_u);


  /* create adjanceny entry for v */
  adj_v = alloc_adjnode(this_,u,costs);
  insert_adjnode(this_,v,adj_v);
  
  /* create link for reverse edge */
  adj_u -> reverse = adj_v;
  adj_v -> reverse = adj_u;
}

/* delete edge */
static
void delete_edge(pbqp *this_,int u,int v)
{
  adjnode *adj_ptr;
  adjnode *rev;
  
  assert(this_ != NULL);
  assert( u >= 0 && u < this_->num_nodes);
  assert( v >= 0 && v < this_->num_nodes);

  adj_ptr=find_adjnode(this_,u,v);
  assert(adj_ptr != NULL);
  assert(adj_ptr->reverse != NULL);

  delete adj_ptr -> costs;
 
  rev = adj_ptr->reverse; 
  remove_adjnode(this_,u,adj_ptr);
  remove_adjnode(this_,v,rev);
} 

/*****************************************************************************
 * cost functions 
 ****************************************************************************/

/* Note: Since cost(u,v) = transpose(cost(v,u)), it would be necessary to store 
   two matrices for both edges (u,v) and (v,u). However, we only store the 
   matrix for the case u < v. For the other case we transpose the stored matrix
   if required. 
*/

/* add costs to cost vector of a node */
void add_pbqp_nodecosts(pbqp *this_,int u, PBQPVector *costs)
{
  assert(this_ != NULL);
  assert(costs != NULL);
  assert(u >= 0 && u <= this_->num_nodes);
  
  if (!this_->node_costs[u]) {
    this_->node_costs[u] = new PBQPVector(*costs);
  } else {
    *this_->node_costs[u] += *costs;
  }
}

/* get cost matrix ptr */
static
PBQPMatrix *get_costmatrix_ptr(pbqp *this_, int u, int v)
{
  adjnode *adj_ptr;
  PBQPMatrix *m = NULL;

  assert (this_ != NULL);
  assert (u >= 0 && u < this_->num_nodes);
  assert (v >= 0 && v < this_->num_nodes); 

  adj_ptr = find_adjnode(this_,u,v);

  if (adj_ptr != NULL) {
    m = adj_ptr -> costs;
  } 

  return m;
}

/* get cost matrix ptr */
/* Note: only the pointer is returned for 
   cost(u,v), if u < v.
*/ 
static
PBQPMatrix *pbqp_get_costmatrix(pbqp *this_, int u, int v)
{
  adjnode *adj_ptr = find_adjnode(this_,u,v);
  
  if (adj_ptr != NULL) {
    if ( u < v) {
      return new PBQPMatrix(*adj_ptr->costs);
    } else {
      return new PBQPMatrix(adj_ptr->costs->transpose());
    }
  } else {
    return NULL;
  }  
}

/* add costs to cost matrix of an edge */
void add_pbqp_edgecosts(pbqp *this_,int u,int v, PBQPMatrix *costs)
{
  PBQPMatrix *adj_costs;

  assert(this_!= NULL);
  assert(costs != NULL);
  assert(u >= 0 && u <= this_->num_nodes);
  assert(v >= 0 && v <= this_->num_nodes);
  
  /* does the edge u-v exists ? */
  if (u == v) {
    PBQPVector *diag = new PBQPVector(costs->diagonalize());
    add_pbqp_nodecosts(this_,v,diag);
    delete diag;
  } else if ((adj_costs = get_costmatrix_ptr(this_,u,v))!=NULL) {
    if ( u < v) {
      *adj_costs += *costs;
    } else {
      *adj_costs += costs->transpose();
    }
  } else {
    adj_costs = new PBQPMatrix((u < v) ? *costs : costs->transpose());
    insert_edge(this_,u,v,adj_costs);
  } 
}

/* remove bucket from bucket list */
static
void pbqp_remove_bucket(pbqp *this_, bucketnode *bucket)
{
  int u = bucket->u;
  
  assert(this_ != NULL);
  assert(u >= 0 && u < this_->num_nodes);
  assert(this_->bucket_list != NULL);
  assert(this_->bucket_ptr[u] != NULL);
  
  /* update predecessor node in bucket list 
     (if no preceeding bucket exists, then
     the bucket_list pointer needs to be 
     updated.)
  */    
  if (bucket->prev != NULL) {
    bucket->prev-> succ = bucket->succ; 
  } else {
    this_->bucket_list[this_->node_deg[u]] = bucket -> succ;
  }
  
  /* update successor node in bucket list */ 
  if (bucket->succ != NULL) { 
    bucket->succ-> prev = bucket->prev;
  }
}

/**********************************************************************************
 * pop functions
 **********************************************************************************/

/* pop node of given degree */
static
int pop_node(pbqp *this_,int deg)
{
  bucketnode *bucket;
  int u;

  assert(this_ != NULL);
  assert(deg >= 0 && deg <= this_->max_deg);
  assert(this_->bucket_list != NULL);
   
  /* get first bucket of bucket list */
  bucket = this_->bucket_list[deg];
  assert(bucket != NULL);

  /* remove bucket */
  pbqp_remove_bucket(this_,bucket);
  u = bucket->u;
  free(bucket);
  return u;
}

/**********************************************************************************
 * reorder functions
 **********************************************************************************/

/* add bucket to bucketlist */
static
void add_to_bucketlist(pbqp *this_,bucketnode *bucket, int deg)
{
  bucketnode *old_head;
  
  assert(bucket != NULL);
  assert(this_ != NULL);
  assert(deg >= 0 && deg <= this_->max_deg);
  assert(this_->bucket_list != NULL);

  /* store node degree (for re-ordering purposes)*/
  this_->node_deg[bucket->u] = deg;
  
  /* put bucket to front of doubly chained list */
  old_head = this_->bucket_list[deg];
  bucket -> prev = NULL;
  bucket -> succ = old_head;
  this_ -> bucket_list[deg] = bucket;
  if (bucket -> succ != NULL ) {
    assert ( old_head -> prev == NULL);
    old_head -> prev = bucket;
  }
}


/* reorder node in bucket list according to 
   current node degree */
static
void reorder_node(pbqp *this_, int u)
{
  int deg; 
  
  assert(this_ != NULL);
  assert(u>= 0 && u < this_->num_nodes);
  assert(this_->bucket_list != NULL);
  assert(this_->bucket_ptr[u] != NULL);

  /* get current node degree */
  deg = get_deg(this_,u);
  
  /* remove bucket from old bucket list only
     if degree of node has changed. */
  if (deg != this_->node_deg[u]) {
    pbqp_remove_bucket(this_,this_->bucket_ptr[u]);
    add_to_bucketlist(this_,this_->bucket_ptr[u],deg);
  } 
}

/* reorder adj. nodes of a node */
static
void reorder_adjnodes(pbqp *this_,int u)
{
  adjnode *adj_ptr;
  
  assert(this_!= NULL);
  assert(u >= 0 && u <= this_->num_nodes);
  assert(this_->adj_list != NULL);

  for(adj_ptr = this_ -> adj_list[u]; adj_ptr != NULL; adj_ptr = adj_ptr -> succ) {
    reorder_node(this_,adj_ptr->adj);
  }
}

/**********************************************************************************
 * creation functions
 **********************************************************************************/

/* create new bucket entry */
/* consistency of the bucket list is not checked! */
static
void create_bucket(pbqp *this_,int u,int deg)
{
  bucketnode *bucket;
  
  assert(this_ != NULL);
  assert(u >= 0 && u < this_->num_nodes);
  assert(this_->bucket_list != NULL);
  
  bucket = (bucketnode *)malloc(sizeof(bucketnode));
  assert(bucket != NULL);

  bucket -> u = u;
  this_->bucket_ptr[u] = bucket;

  add_to_bucketlist(this_,bucket,deg);
}

/* create bucket list */
static
void create_bucketlist(pbqp *this_)
{
  int u;
  int max_deg;
  int deg;

  assert(this_ != NULL);
  assert(this_->bucket_list == NULL);

  /* determine max. degree of the nodes */
  max_deg = 2;  /* at least of degree two! */
  for(u=0;u<this_->num_nodes;u++) {
    deg = this_->node_deg[u] = get_deg(this_,u);
    if (deg > max_deg) {
      max_deg = deg;
    }
  }
  this_->max_deg = max_deg;
  
  /* allocate bucket list */
  this_ -> bucket_list = (bucketnode **)malloc(sizeof(bucketnode *)*(max_deg + 1));
  memset(this_->bucket_list,0,sizeof(bucketnode *)*(max_deg + 1));
  assert(this_->bucket_list != NULL);
  
  /* insert nodes to the list */
  for(u=0;u<this_->num_nodes;u++) {
    create_bucket(this_,u,this_->node_deg[u]);  
  }
}

/*****************************************************************************
 * PBQP simplification for trivial nodes
 ****************************************************************************/

/* remove trivial node with cost vector length of one */
static
void disconnect_trivialnode(pbqp *this_,int u)
{
  int v;
  adjnode *adj_ptr, 
          *next;
  PBQPMatrix *c_uv;
  PBQPVector *c_v;
  
  assert(this_ != NULL);
  assert(this_->node_costs != NULL);
  assert(u >= 0 && u < this_ -> num_nodes);
  assert(this_->node_costs[u]->getLength() == 1);
  
  /* add edge costs to node costs of adj. nodes */
  for(adj_ptr = this_->adj_list[u]; adj_ptr != NULL; adj_ptr = next){
    next = adj_ptr -> succ;
    v = adj_ptr -> adj;
    assert(v >= 0 && v < this_ -> num_nodes);
    
    /* convert matrix to cost vector offset for adj. node */
    c_uv = pbqp_get_costmatrix(this_,u,v);
    c_v = new PBQPVector(c_uv->getRowAsVector(0));
    *this_->node_costs[v] += *c_v;
    
    /* delete edge & free vec/mat */
    delete c_v;
    delete c_uv;
    delete_edge(this_,u,v);
  }   
}

/* find all trivial nodes and disconnect them */
static   
void eliminate_trivial_nodes(pbqp *this_)
{
   int u;
   
   assert(this_ != NULL);
   assert(this_ -> node_costs != NULL);
   
   for(u=0;u < this_ -> num_nodes; u++) {
     if (this_->node_costs[u]->getLength() == 1) {
       disconnect_trivialnode(this_,u); 
     }
   }
}

/*****************************************************************************
 * Normal form for PBQP 
 ****************************************************************************/

/* simplify a cost matrix. If the matrix
   is independent, then simplify_matrix
   returns true - otherwise false. In
   vectors u and v the offset values of
   the decomposition are stored. 
*/

static
bool normalize_matrix(PBQPMatrix *m, PBQPVector *u, PBQPVector *v)
{
  assert( m != NULL);
  assert( u != NULL);
  assert( v != NULL);
  assert( u->getLength() > 0);
  assert( v->getLength() > 0);
  
  assert(m->getRows() == u->getLength());
  assert(m->getCols() == v->getLength());

  /* determine u vector */
  for(unsigned r = 0; r < m->getRows(); ++r) {
    PBQPNum min = m->getRowMin(r);
    (*u)[r] += min;
    if (!isInf(min)) {
      m->subFromRow(r, min);
    } else {
      m->setRow(r, 0);
    }
  }
  
  /* determine v vector */
  for(unsigned c = 0; c < m->getCols(); ++c) {
    PBQPNum min = m->getColMin(c);
    (*v)[c] += min;
    if (!isInf(min)) {
      m->subFromCol(c, min);
    } else {
      m->setCol(c, 0);
    }
  }
  
  /* determine whether matrix is 
     independent or not. 
    */
  return m->isZero();
}

/* simplify single edge */
static
void simplify_edge(pbqp *this_,int u,int v)
{
  PBQPMatrix *costs;
  bool is_zero; 
  
  assert (this_ != NULL);
  assert (u >= 0 && u <this_->num_nodes);
  assert (v >= 0 && v <this_->num_nodes);
  assert (u != v);

  /* swap u and v  if u > v in order to avoid un-necessary
     tranpositions of the cost matrix */
  
  if (u > v) {
    int swap = u;
    u = v;
    v = swap;
  }
  
  /* get cost matrix and simplify it */  
  costs = get_costmatrix_ptr(this_,u,v);
  is_zero=normalize_matrix(costs,this_->node_costs[u],this_->node_costs[v]);

  /* delete edge */
  if(is_zero){
    delete_edge(this_,u,v);
    this_->changed = true;
  }
}

/* normalize cost matrices and remove 
   edges in PBQP if they ary independent, 
   i.e. can be decomposed into two 
   cost vectors.
*/
static
void eliminate_independent_edges(pbqp *this_)
{
  int u,v;
  adjnode *adj_ptr,*next;
  
  assert(this_ != NULL);
  assert(this_ -> adj_list != NULL);

  this_->changed = false;
  for(u=0;u < this_->num_nodes;u++) {
    for (adj_ptr = this_ -> adj_list[u]; adj_ptr != NULL; adj_ptr = next) {
      next = adj_ptr -> succ;
      v = adj_ptr -> adj;
      assert(v >= 0 && v < this_->num_nodes);
      if (u < v) {
        simplify_edge(this_,u,v);
      } 
    }
  }
}


/*****************************************************************************
 * PBQP reduction rules 
 ****************************************************************************/

/* RI reduction
   This reduction rule is applied for nodes 
   of degree one. */

static
void apply_RI(pbqp *this_,int x)
{
  int y;
  unsigned xlen,
           ylen;
  PBQPMatrix *c_yx;
  PBQPVector *c_x, *delta;
  
  assert(this_ != NULL);
  assert(x >= 0 && x < this_->num_nodes);
  assert(this_ -> adj_list[x] != NULL);
  assert(this_ -> adj_list[x] -> succ == NULL);

  /* get adjacence matrix */
  y = this_ -> adj_list[x] -> adj;
  assert(y >= 0 && y < this_->num_nodes);
  
  /* determine length of cost vectors for node x and y */
  xlen = this_ -> node_costs[x]->getLength();
  ylen = this_ -> node_costs[y]->getLength();

  /* get cost vector c_x and matrix c_yx */
  c_x = this_ -> node_costs[x];
  c_yx = pbqp_get_costmatrix(this_,y,x); 
  assert (c_yx != NULL);

  
  /* allocate delta vector */
  delta = new PBQPVector(ylen);

  /* compute delta vector */
  for(unsigned i = 0; i < ylen; ++i) {
    PBQPNum min =  (*c_yx)[i][0] + (*c_x)[0];
    for(unsigned j = 1; j < xlen; ++j) {
      PBQPNum c =  (*c_yx)[i][j] + (*c_x)[j];
      if ( c < min )  
         min = c;
    }
    (*delta)[i] = min; 
  } 

  /* add delta vector */
  *this_ -> node_costs[y] += *delta;

  /* delete node x */
  remove_node(this_,x);

  /* reorder adj. nodes of node x */
  reorder_adjnodes(this_,x);

  /* push node x on stack */
  assert(this_ -> stack_ptr < this_ -> num_nodes);
  this_->stack[this_ -> stack_ptr++] = x;

  /* free vec/mat */
  delete c_yx;
  delete delta;

  /* increment counter for number statistic */
  this_->num_ri++;
}

/* RII reduction
   This reduction rule is applied for nodes 
   of degree two. */

static
void apply_RII(pbqp *this_,int x)
{
  int y,z; 
  unsigned xlen,ylen,zlen;
  adjnode *adj_yz;

  PBQPMatrix *c_yx, *c_zx;
  PBQPVector *cx;
  PBQPMatrix *delta;
 
  assert(this_ != NULL);
  assert(x >= 0 && x < this_->num_nodes);
  assert(this_ -> adj_list[x] != NULL);
  assert(this_ -> adj_list[x] -> succ != NULL);
  assert(this_ -> adj_list[x] -> succ -> succ == NULL);

  /* get adjacence matrix */
  y = this_ -> adj_list[x] -> adj;
  z = this_ -> adj_list[x] -> succ -> adj;
  assert(y >= 0 && y < this_->num_nodes);
  assert(z >= 0 && z < this_->num_nodes);
  
  /* determine length of cost vectors for node x and y */
  xlen = this_ -> node_costs[x]->getLength();
  ylen = this_ -> node_costs[y]->getLength();
  zlen = this_ -> node_costs[z]->getLength();

  /* get cost vector c_x and matrix c_yx */
  cx = this_ -> node_costs[x];
  c_yx = pbqp_get_costmatrix(this_,y,x); 
  c_zx = pbqp_get_costmatrix(this_,z,x); 
  assert(c_yx != NULL);
  assert(c_zx != NULL);

  /* Colour Heuristic */
  if ( (adj_yz = find_adjnode(this_,y,z)) != NULL) {
    adj_yz->tc_valid = false;
    adj_yz->reverse->tc_valid = false; 
  }

  /* allocate delta matrix */
  delta = new PBQPMatrix(ylen, zlen);

  /* compute delta matrix */
  for(unsigned i=0;i<ylen;i++) {
    for(unsigned j=0;j<zlen;j++) {
      PBQPNum min = (*c_yx)[i][0] + (*c_zx)[j][0] + (*cx)[0];
      for(unsigned k=1;k<xlen;k++) {
        PBQPNum c = (*c_yx)[i][k] + (*c_zx)[j][k] + (*cx)[k];
        if ( c < min ) {
          min = c;
        }
      }
      (*delta)[i][j] = min;
    }
  }

  /* add delta matrix */
  add_pbqp_edgecosts(this_,y,z,delta);

  /* delete node x */
  remove_node(this_,x);

  /* simplify cost matrix c_yz */
  simplify_edge(this_,y,z);

  /* reorder adj. nodes */
  reorder_adjnodes(this_,x);

  /* push node x on stack */
  assert(this_ -> stack_ptr < this_ -> num_nodes);
  this_->stack[this_ -> stack_ptr++] = x;

  /* free vec/mat */
  delete c_yx;
  delete c_zx;
  delete delta;

  /* increment counter for number statistic */
  this_->num_rii++;

}

/* RN reduction */
static
void apply_RN(pbqp *this_,int x)
{
  unsigned xlen;

  assert(this_ != NULL);
  assert(x >= 0 && x < this_->num_nodes);
  assert(this_ -> node_costs[x] != NULL);

  xlen = this_ -> node_costs[x] -> getLength();

  /* after application of RN rule no optimality
     can be guaranteed! */
  this_ -> optimal = false;
  
  /* push node x on stack */
  assert(this_ -> stack_ptr < this_ -> num_nodes);
  this_->stack[this_ -> stack_ptr++] = x;

  /* delete node x */ 
  remove_node(this_,x);

  /* reorder adj. nodes of node x */
  reorder_adjnodes(this_,x);

  /* increment counter for number statistic */
  this_->num_rn++;
}


static
void compute_tc_info(pbqp *this_, adjnode *p)
{
   adjnode *r;
   PBQPMatrix *m;
   int x,y;
   PBQPVector *c_x, *c_y;
   int *row_inf_counts;

   assert(p->reverse != NULL);

   /* set flags */ 
   r = p->reverse;
   p->tc_valid = true;
   r->tc_valid = true;

   /* get edge */
   x = r->adj;
   y = p->adj;

   /* get cost vectors */
   c_x = this_ -> node_costs[x];
   c_y = this_ -> node_costs[y];

   /* get cost matrix */
   m = pbqp_get_costmatrix(this_, x, y);

  
   /* allocate allowed set for edge (x,y) and (y,x) */
   if (p->tc_safe_regs == NULL) {
     p->tc_safe_regs = (int *) malloc(sizeof(int) * c_x->getLength());
   } 

   if (r->tc_safe_regs == NULL ) {
     r->tc_safe_regs = (int *) malloc(sizeof(int) * c_y->getLength());
   }

   p->tc_impact = r->tc_impact = 0;

   row_inf_counts = (int *) alloca(sizeof(int) * c_x->getLength());

   /* init arrays */
   p->tc_safe_regs[0] = 0;
   row_inf_counts[0] = 0;
   for(unsigned i = 1; i < c_x->getLength(); ++i){
     p->tc_safe_regs[i] = 1;
     row_inf_counts[i] = 0;
   }

   r->tc_safe_regs[0] = 0;
   for(unsigned j = 1; j < c_y->getLength(); ++j){
     r->tc_safe_regs[j] = 1;
   }

   for(unsigned j = 0; j < c_y->getLength(); ++j) {
      int col_inf_counts = 0;
      for (unsigned i = 0; i < c_x->getLength(); ++i) {
         if (isInf((*m)[i][j])) {
              ++col_inf_counts;
              ++row_inf_counts[i];
         
              p->tc_safe_regs[i] = 0;
              r->tc_safe_regs[j] = 0;
         }
      }
      if (col_inf_counts > p->tc_impact) {
           p->tc_impact = col_inf_counts;
      }
   }

   for(unsigned i = 0; i < c_x->getLength(); ++i){
     if (row_inf_counts[i] > r->tc_impact)
     {
           r->tc_impact = row_inf_counts[i];
     }
   }
           
   delete m;
}

/* 
 * Checks whether node x can be locally coloured. 
 */
static 
int is_colorable(pbqp *this_,int x)
{
  adjnode *adj_ptr;
  PBQPVector *c_x;
  int result = 1;
  int *allowed;
  int num_allowed = 0;
  unsigned total_impact = 0;

  assert(this_ != NULL);
  assert(x >= 0 && x < this_->num_nodes);
  assert(this_ -> node_costs[x] != NULL);

  c_x = this_ -> node_costs[x];

  /* allocate allowed set */
  allowed = (int *)malloc(sizeof(int) * c_x->getLength());
  for(unsigned i = 0; i < c_x->getLength(); ++i){
    if (!isInf((*c_x)[i]) && i > 0) {
      allowed[i] = 1;
      ++num_allowed;
    } else { 
      allowed[i] = 0;
    }
  }

  /* determine local minimum */
  for(adj_ptr=this_->adj_list[x] ;adj_ptr != NULL; adj_ptr = adj_ptr -> succ) {
      if (!adj_ptr -> tc_valid) { 
          compute_tc_info(this_, adj_ptr);
      }

      total_impact += adj_ptr->tc_impact;

      if (num_allowed > 0) {
          for (unsigned i = 1; i < c_x->getLength(); ++i){
            if (allowed[i]){
              if (!adj_ptr->tc_safe_regs[i]){
                allowed[i] = 0;
                --num_allowed;
                if (num_allowed == 0)
                    break;
              }
            }
          }
      }
      
      if ( total_impact >= c_x->getLength() - 1 && num_allowed == 0 ) {
         result = 0;
         break;
      }
  }
  free(allowed);

  return result;
}

/* use briggs heuristic 
  note: this_ is not a general heuristic. it only is useful for 
  interference graphs.
 */
int pop_colorablenode(pbqp *this_)
{
  int deg;
  bucketnode *min_bucket=NULL;
  PBQPNum min = std::numeric_limits<PBQPNum>::infinity();
 
  /* select node where the number of colors is less than the node degree */
  for(deg=this_->max_deg;deg > 2;deg--) {
    bucketnode *bucket;
    for(bucket=this_->bucket_list[deg];bucket!= NULL;bucket = bucket -> succ) {
      int u = bucket->u;
      if (is_colorable(this_,u)) {
        pbqp_remove_bucket(this_,bucket);
        this_->num_rn_special++;
        free(bucket);
        return u; 
      } 
    }
  }

  /* select node with minimal ratio between average node costs and degree of node */
  for(deg=this_->max_deg;deg >2; deg--) {
    bucketnode *bucket;
    for(bucket=this_->bucket_list[deg];bucket!= NULL;bucket = bucket -> succ) {
      PBQPNum h;
      int u;
 
      u = bucket->u;
      assert(u>=0 && u < this_->num_nodes);
      h = (*this_->node_costs[u])[0] / (PBQPNum) deg;
      if (h < min) {
        min_bucket = bucket;
        min = h;
      } 
    }
  }

  /* return node and free bucket */
  if (min_bucket != NULL) {
    int u;

    pbqp_remove_bucket(this_,min_bucket);
    u = min_bucket->u;
    free(min_bucket);
    return u;
  } else {
    return -1;
  }
}


/*****************************************************************************
 * PBQP graph parsing
 ****************************************************************************/
 
/* reduce pbqp problem (first phase) */
static
void reduce_pbqp(pbqp *this_)
{
  int u;

  assert(this_ != NULL);
  assert(this_->bucket_list != NULL);

  for(;;){

    if (this_->bucket_list[1] != NULL) {
      u = pop_node(this_,1);
      apply_RI(this_,u); 
    } else if (this_->bucket_list[2] != NULL) {
      u = pop_node(this_,2);
      apply_RII(this_,u);
    } else if ((u = pop_colorablenode(this_)) != -1) {
      apply_RN(this_,u);
    } else {
      break;
    }
  } 
}

/*****************************************************************************
 * PBQP back propagation
 ****************************************************************************/

/* determine solution of a reduced node. Either
   RI or RII was applied for this_ node. */
static
void determine_solution(pbqp *this_,int x)
{
  PBQPVector *v = new PBQPVector(*this_ -> node_costs[x]);
  adjnode *adj_ptr;

  assert(this_ != NULL);
  assert(x >= 0 && x < this_->num_nodes);
  assert(this_ -> adj_list != NULL);
  assert(this_ -> solution != NULL);

  for(adj_ptr=this_->adj_list[x] ;adj_ptr != NULL; adj_ptr = adj_ptr -> succ) {
    int y = adj_ptr -> adj;
    int y_sol = this_ -> solution[y];

    PBQPMatrix *c_yx = pbqp_get_costmatrix(this_,y,x);
    assert(y_sol >= 0 && y_sol < (int)this_->node_costs[y]->getLength());
    (*v) += c_yx->getRowAsVector(y_sol);
    delete c_yx;
  }
  this_ -> solution[x] = v->minIndex();

  delete v;
}

/* back popagation phase of PBQP */
static
void back_propagate(pbqp *this_)
{
   int i;

   assert(this_ != NULL);
   assert(this_->stack != NULL);
   assert(this_->stack_ptr < this_->num_nodes);

   for(i=this_ -> stack_ptr-1;i>=0;i--) {
      int x = this_ -> stack[i];
      assert( x >= 0 && x < this_ -> num_nodes);
      reinsert_node(this_,x);
      determine_solution(this_,x);
   }
}

/* solve trivial nodes of degree zero */
static
void determine_trivialsolution(pbqp *this_)
{
  int u;
  PBQPNum delta;

  assert( this_ != NULL);
  assert( this_ -> bucket_list != NULL);

  /* determine trivial solution */
  while (this_->bucket_list[0] != NULL) {
    u = pop_node(this_,0);

    assert( u >= 0 && u < this_ -> num_nodes);

    this_->solution[u] = this_->node_costs[u]->minIndex();
    delta = (*this_->node_costs[u])[this_->solution[u]];
    this_->min = this_->min + delta;

    /* increment counter for number statistic */
    this_->num_r0++;
  }
}

/*****************************************************************************
 * debug facilities
 ****************************************************************************/
static
void check_pbqp(pbqp *this_)
{
  int u,v;
  PBQPMatrix *costs;
  adjnode *adj_ptr;
  
  assert( this_ != NULL);
  
  for(u=0;u< this_->num_nodes; u++) {
    assert (this_ -> node_costs[u] != NULL);
    for(adj_ptr = this_ -> adj_list[u];adj_ptr != NULL; adj_ptr = adj_ptr -> succ) {
      v = adj_ptr -> adj;
      assert( v>= 0 && v < this_->num_nodes);
      if (u < v ) {
        costs = adj_ptr -> costs;
        assert( costs->getRows() == this_->node_costs[u]->getLength() &&
                costs->getCols() == this_->node_costs[v]->getLength());
      }           
    }
  }
}

/*****************************************************************************
 * PBQP solve routines 
 ****************************************************************************/

/* solve PBQP problem */
void solve_pbqp(pbqp *this_)
{
  assert(this_ != NULL);
  assert(!this_->solved); 
  
  /* check vector & matrix dimensions */
  check_pbqp(this_);

  /* simplify PBQP problem */  
  
  /* eliminate trivial nodes, i.e.
     nodes with cost vectors of length one.  */
  eliminate_trivial_nodes(this_); 

  /* eliminate edges with independent 
     cost matrices and normalize matrices */
  eliminate_independent_edges(this_);
  
  /* create bucket list for graph parsing */
  create_bucketlist(this_);
  
  /* reduce phase */
  reduce_pbqp(this_);
  
  /* solve trivial nodes */
  determine_trivialsolution(this_);

  /* back propagation phase */
  back_propagate(this_); 
  
  this_->solved = true;
}

/* get solution of a node */
int get_pbqp_solution(pbqp *this_,int x)
{
  assert(this_ != NULL);
  assert(this_->solution != NULL);
  assert(this_ -> solved);
  
  return this_->solution[x];
}

/* is solution optimal? */
bool is_pbqp_optimal(pbqp *this_)
{
  assert(this_ -> solved);
  return this_->optimal;
}

} 

/* end of pbqp.c */