/*========================================================================

  File: REA.c

  Program REA - Version 1.1 - July 1999

  ========================================================================

  This module contains an implementation of the Recursive Enumeration
  Algorithm (REA) that enumerates (by increasing weight) the N shortest
  paths in weighted graphs. The algorithm is described in: 

    "Computing the K Shortest Paths: a New Algorithm and an
    Experimental Comparison" by Victor Jimenez and Andres Marzal,
    3rd Workshop on Algorithm Engineering, London, July 1999.
    To be published by Springer-Verlag in the LNCS series.

  The sets of candidate paths are implemented using binary heaps.

  ========================================================================

    Copyright (C) 1999 Victor Jimenez and Andres Marzal

    This program is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 2 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program (file COPYING); if not, write to the Free
    Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.

    You can contact the authors at:

           Victor Jimenez and Andres Marzal
           Dept. de Informatica, Universitat Jaume I, Castellon, Spain
           {vjimenez,amarzal}@inf.uji.es

    Updated information about this software will be available at the
    following www address:

           http://terra.act.uji.es/REA
 
  ========================================================================= */

// MODIFIED by Mark Ziegelmann  <mark@mpi-sb.mpg.de>

#define VERSION "1.1"

#include <stdlib.h>
#include <stdio.h>
#include <time.h>  
#include <unistd.h>  

#include <LEDA/array.h>
#include <LEDA/tuple.h>
#include <LEDA/graph.h>

#include <assert.h>

#include <CNOP/REA_k.h>
#include <CNOP/dijkstra_k.h>


/*============================================================================
  GLOBAL VARIABLES
  ============================================================================*/

/* A buffer for heap elements is dinamically allocated by the main function  */
/* with size the number of arcs in the graph (which is the worst case bound) */
/* New heaps are initialized using this buffer by the function CreateHeap_k    */
Path **hepsBuffer_k;
Path **firstFreeHeap_k;

/* For experimental purposes we avoid to use the malloc system for new paths */
/* New paths are taken from this buffer by the function CreatePath           */
#define MAX_PATHS 1000000
Path pathsBuffer_k[MAX_PATHS];
int  numberUsedPaths_k = 0;
Path *firstFreePath_k = pathsBuffer_k;

typedef array<double> Point; 
extern double get_plane_coeff(const array<Point>& current_points, int d, int coeff_num, int add=0);
extern bool point_in_simplex(const array<Point>& current_points, const Point& p, int k);

/*===========================================================================*/
 Path *CreatePath (Node *node, Path *backPath, COST_TYPE cost, COST_TYPE length, const array<COST_TYPE>& resource) {
   int k=resource.size();
  Path   *newPath;
  
#ifdef DEBUG
  assert (node != NULL);
#endif
  
  if (++numberUsedPaths_k > MAX_PATHS) {
    fprintf (stderr, "Redefine MAX_PATHS");
    exit(1);
  }
  newPath = firstFreePath_k++;
  
  newPath->backPath = backPath;
  newPath->cost     = cost;
  newPath->length     = length;
  newPath->resource = new array<COST_TYPE>(k);
  for (int i=0;i<k;i++)
    (*(newPath->resource))[i]     = resource[i];

  newPath->lastNode = node;
  
  return newPath;
  
}

//(*(aa.A))[1]=2;
/*===========================================================================*/
void PrintPath_k (Path *path) {
  /* Prints the path in reverse order (with the final node first). */
  Path *backPath;
  
#ifdef DEBUG
  assert (path != NULL); 
#endif
 if (path->cost < INFINITY_COST) {
  for (backPath = path; backPath != NULL; backPath = backPath->backPath)
    printf ("%i-", backPath->lastNode->name);
  printf (" \t(Cost: %f)", path->cost);
 }
  
}



#define COST(_i_) (H->elt[_i_]->cost)

/*===========================================================================*/
 void CreateHeap_k (Heap *H, int dimension) {
  /* Allocates memory for the heap elements from the heaps buffer 
     and initializes the number of element inside the heap to 0 */
#ifdef DEBUG
  assert (dimension >= 0);
  H->dimension = dimension;
#endif
  H->size = 0;
  H->elt = firstFreeHeap_k;
  firstFreeHeap_k += dimension;
}


/*===========================================================================*/
 void PreInsertInHeap_k(Heap *H, Path *elt) {
  /* Inserts a new element in the heap in time O(1)
     without preserving the heap property */
#ifdef DEBUG
  assert (H->size < H->dimension);
#endif
  
  H->elt[H->size++] = elt;
}


/*===========================================================================*/
 void BuildHeap_k(Heap *H) {
  /* Obtains the heap property in time linear with the heap size */
  
  register int i;
  for (i=H->size/2; i>=0; i--) {
    register Path * item = H->elt[i];
    register int j = 2*i+1;
    while (j<H->size) {
      if (j<H->size-1 && COST(j+1) < COST(j)) j++;
      if (item->cost <= COST(j)) break;
      H->elt[(j-1)/2] = H->elt[j];
      j = 2*j+1;
    }
    H->elt[(j-1)/2] = item;
  }
}


/*===========================================================================*/
 Path * DeleteBestInHeap_k (Heap *H) {
  /* Returns the element with minimum cost in the heap and deletes
     it from the heap, restoring the heap property in worst case time
     logarithmic with the heap size */
  Path *best;
  
  if (H->size == 0) return (NULL);
  best = H->elt[0];
  H->elt[0] = H->elt[--H->size];
  {
    register Path * item = H->elt[0];
    register int j = 1;
    while (j<H->size) {
      if (j<H->size-1 && COST(j+1) < COST(j)) j++;
      if (item->cost <= COST(j)) break;
      H->elt[(j-1)/2] = H->elt[j];
      j = 2*j+1;
    }
    H->elt[(j-1)/2] = item;
  }
  return (best);
}

/*===========================================================================*/
 Path * BestInHeap_k (Heap *H) {
     /* Returns the element with minimum cost in the heap without deleting it*/

  if (H->size == 0) return (NULL);
  return H->elt[0];
}

/*===========================================================================*/
 void UpdateBestInHeap_k (Heap *H, Path *elt) {
     /* Deletes the element with minimum cost in the heap and inserts a new one
	preserving the heap property, in worst case time logarithmic with the
	heap size */

  register Path * item = elt;
  register int j = 1;
  while (j<H->size) {
    if (j<H->size-1 && COST(j+1) < COST(j)) j++;
    if (item->cost <= COST(j)) break;
    H->elt[(j-1)/2] = H->elt[j];
    j = 2*j+1;
  }
  H->elt[(j-1)/2] = item;
}


/*===========================================================================*/
 void InitializeCandidateSet (Node *node, int k) {
  /* The set of candidates is  initialized with the best path from
     each  predecessor node, except  the  one from which  the best
     path at the current node  comes. If the current node is the 
     initial node in  the graph, all its predecessor nodes (if some
     exists) provide a candidate */

  PtrArc *ptrArc;
  int    numberArcs;
  Path   *path,
         *newCand;
  
  CreateHeap_k (&node->heap, node->numberArcsIn);
  
  for (ptrArc = node->firstArcIn, numberArcs = node->numberArcsIn;
       numberArcs != 0; ptrArc++, numberArcs--)
    if ( (path = (*ptrArc)->source->bestPath) != node->bestPath->backPath) {
      /* It is important to compare pointers and not costs or node names,
         because  several candidates could came from the same node (if
         the graph has parallel arcs) and/or with the same cost. */
      array<COST_TYPE> new_resA(k);
      for(int i=0;i<k;i++)
        new_resA[i]=(*(path->resource))[i]+(*((*ptrArc)->resource))[i];
      newCand = CreatePath (node, path, path->cost + (*ptrArc)->cost, path->length + (*ptrArc)->length, new_resA);
      PreInsertInHeap_k (&node->heap, newCand);
    }
  BuildHeap_k(&node->heap);
}

int loopfree_k(Path* path, Node* node) {
  Path *backPath;
  for (backPath = path; backPath != NULL; backPath = backPath->backPath) {
    if (node==backPath->lastNode)
      return 0;
  }
  return 1;
}

   
/*===========================================================================*/
Path* NextPath (Path *path,int k, bool loopless_variant) {
  /* Central routine of the Recursive Enumeration Algorithm: computes the
     next path from the initial node to the same node in which the argument
     path ends, assuming that it has not been computed before.
  */
     
  Node      *node         = NULL;
  Path      *backPath     = NULL,
            *nextBackPath = NULL,
            *bestCand     = NULL,
            *newCand      = NULL;
  COST_TYPE arcCost,arcLength;
  
#ifdef DEBUG
  assert (path != NULL);
  assert (path->nextPath == NULL);
  assert (path->lastNode != NULL);
#endif
  
  node = path->lastNode;
  
#ifdef DEBUG
  printf ("\nComputing next path at node %i\n", node->name);
#endif
  
  if (node->heap.elt == NULL) 
    /* This is done here instead of in Dijkstra function, so that it is */
    /* only done for nodes in which alternative paths are required.     */
    InitializeCandidateSet (node,k);
  
  if ((backPath = path->backPath) != NULL) {
    nextBackPath = backPath->nextPath;
    if (nextBackPath == NULL)
      nextBackPath = NextPath (backPath,k, loopless_variant);
    if (nextBackPath != NULL) {
      arcCost = path->cost - backPath->cost;
      arcLength = path->length - backPath->length;
      array<COST_TYPE> new_resA(k);
      for(int i=0;i<k;i++)
        new_resA[i]=(*(nextBackPath->resource))[i]+(*(path->resource))[i] - (*(backPath->resource))[i];
      if (nextBackPath->nextPath!=NULL) {

	if (!loopless_variant || loopfree_k(nextBackPath,node))
	  newCand = CreatePath (node, nextBackPath, nextBackPath->cost + arcCost, nextBackPath->length + arcLength, new_resA);
	else
	  newCand = CreatePath (node, nextBackPath, MAXINT, nextBackPath->length + arcLength, new_resA);

	/*
	if (nextBackPath->nextPath->lastNode=node)
	  newCand = CreatePath (node, nextBackPath, nextBackPath->cost + arcCost+1000, nextBackPath->length + arcLength, new_resA);
	else
	  newCand = CreatePath (node, nextBackPath, nextBackPath->cost + arcCost, nextBackPath->length + arcLength, new_resA);
	*/
      }
      else
	newCand = CreatePath (node, nextBackPath, nextBackPath->cost + arcCost, nextBackPath->length + arcLength, new_resA);
    }
  }
  
  if (newCand == NULL)
    bestCand = DeleteBestInHeap_k (&node->heap);
  else {
    bestCand = BestInHeap_k (&node->heap);
    if (bestCand != NULL && bestCand->cost < newCand->cost)
      UpdateBestInHeap_k (&node->heap, newCand);
    else
      bestCand = newCand;
  }
  
  if (bestCand == NULL)
    return NULL;
  
  /* Adds the best candidate to the list of paths ending in the same node */
  bestCand->nextPath = NULL;
  path->nextPath     = bestCand;
  
  
  return bestCand;
}

/*===========================================================================*/
void Copyright_k () {
  printf ("\n=====================================================================\n");
  printf ("REA version %s, Copyright (C) 1999 Victor Jimenez and Andres Marzal\n",
	  VERSION);
  printf ("REA comes with ABSOLUTELY NO WARRANTY.\n");
  printf ("This is free software, and you are welcome to redistribute it under\n");
  printf ("certain conditions; see the README and COPYING files for more details.\n");
}





bool feasible(array<COST_TYPE> *res, const array<double>& res_limit) {
  for (int i=0;i<res_limit.size();i++)
    if ((*(res))[i]>res_limit[i])
      return false;
  return true;
}

list<edge> make_edge_list_k(Path *UBpath, const array<node>& nA) {
  Path *backPath;
  int Target,Source;
  Target=UBpath->lastNode->name;
  node v1,v2;
  edge e;
  list<edge> oldUBpath;
  for (backPath = UBpath->backPath; backPath != NULL; backPath = backPath->backPath) {
    Source=backPath->lastNode->name;
    //no parallel arcs!
    //node=graph->node+target-1;
    v1=nA[Source];
    v2=nA[Target];
    forall_adj_edges(e,v1)
      if (target(e)==v2)
        break;
    Target=Source;
    oldUBpath.push(e);
  }
  return oldUBpath;
}

/*===========================================================================*/

bool REA_k(Graph gr, const graph& G, const array<node>& nA, int n_paths, const array<COST_TYPE>& res_limit, const double& n0, const array<double>& normal_vector, double& UB, double& LB, list<edge>& oldUBpath, bool (*is_feasible)(const list<edge>&, const array<double>&, const array<double>&), double delta_max, bool loopless_variant) { 
 
  Path   *path, *backPath, *UBpath;
  int    i, numberPaths = 1;

#ifdef DEBUG
  Node   *node;
  int    numberNodes;
#endif

  int k=normal_vector.size();

  numberUsedPaths_k=0;

  UBpath=NULL;

  numberPaths=n_paths;

  
  /********** Allocates memory for heaps ************************************/
  /* Memory is allocated for the worst case: one element per arc */
  hepsBuffer_k = (Path **) malloc(gr.numberArcs*sizeof(*hepsBuffer_k));
  if (hepsBuffer_k == NULL) {
    perror("Not enough memory for heaps.\n");
    exit(1);
  }
  firstFreeHeap_k = hepsBuffer_k;
  
  Dijkstra_k (&gr,k);


#ifdef DEBUG
  for (node = gr.node, numberNodes = gr.numberNodes; numberNodes != 0; 
       node++, numberNodes--) { 
    printf ("\nBest Path for FinalNode=%i: ", node->name);
    PrintPath_k (node->bestPath);
  }
#endif
  
  /******************** Computes the K shortest paths ***********************/
  i = 2;
  double L_u;
  double bound_sum=0;
  array<COST_TYPE> optRes(k); 
  for(int ii=0;ii<k;ii++)
    bound_sum+=normal_vector[ii]*res_limit[ii];
  path = gr.finalNode->bestPath;

  list<edge> eL;
  array<double> newpoint(k+1);
  
  while (i <= numberPaths && path != NULL) {
      if (LB>=UB || UB-LB<delta_max) {
        break;
      }
      else {
	L_u=(path->cost-bound_sum)/n0;
	LB=L_u; 
	if (feasible(path->resource,res_limit) && path->length<UB) {

	  eL=make_edge_list_k(path,nA);
	  newpoint[0]=path->length;
	  for(int j=0;j<k;j++)
	    newpoint[j]=(*(path->resource))[j];
	  if (is_feasible(eL,newpoint,res_limit)) {
	    UB=path->length;
	    UBpath=path;
	  }    
          
        } 

    }
    path = NextPath (path,k,loopless_variant);
    
    i++;
  }


  int Target;
  if (path!=NULL) {
    Target=path->lastNode->name;
  }

  if (UBpath==NULL) {
    if (i>numberPaths)
      return false;
    else
      return true;
  }
  else {
    Target=UBpath->lastNode->name;
    oldUBpath.clear();
  }
  int Source;
  edge e;
  node v1,v2;
  for (backPath = UBpath->backPath; backPath != NULL; backPath = backPath->backPath) {
    Source=backPath->lastNode->name;
    //no parallel arcs!
    //node=gr->node+target-1;
    v1=nA[Source];
    v2=nA[Target];
    forall_adj_edges(e,v1)
      if (G.target(e)==v2)
        break;
    Target=Source;
    oldUBpath.push(e);
  }
 
  numberUsedPaths_k=0;

  if (i>=numberPaths) 
    return false;
  else
    return true;  
}