/*
  Exhaustive search for scheduling n jobs on a single machine,
  minimizing sum wj*f(Cj) for some penalty function f.

  Uses A* with basic lower bound.

  C.Durr, O.Vasquez, 2013

  log:
  xtof: Jeu 18 oct 2012 17:45:26 CEST -- first version of this program
  xtof: Mar 23 oct 2012 18:11:16 CEST -- implemented rule of WibkeHöhn in order to compare
  xtof: traduire en C++
  xtof: change Dijkstra to A*
  xtof: implemented forward solve
  xtof: extend experiments, new rules local-implies-global
  xtof: [6sep13] allow negative beta
  xtof: [15sep]  simplify program, check conjecture
  xtof: [18sep]  made globalPrec strict
  xtof: [26sep]  constraint can be specified from command line
  xtof: [25oct]  local precedence improved by table with separation times
*/

#include <iostream>
#include <fstream>
#include <map>
#include <set>
#include <queue>
#include <math.h>
#include <cstdlib>
#include <sstream>
#include <cassert>
#include <cstring>
#include <algorithm>

using namespace std;



// toggle flag to produce the exploration graph
#define DUMPGRAPH   false
#define TIMEOUT     1000000

const int BASIC       = 0;
const int HOHNJACOBS  = 1;
const int DURRVASQUEZ = 2;
int rule        = DURRVASQUEZ; // which rule should be applied

const int FORWARD     = 0;
const int BACKWARD    = 1;
const int BEST        = 2;
int direction         = BEST;        // in which direction should the schedules be constructed
bool      forward;                   // the direction that has been choosen according to direction

// ------------------------------------- problem type

// 45 bits are enough for the test cases with 35 jobs and beta=3
typedef double COST;

double beta = 2, sgnBeta = 1;
inline COST penalty(int t) { return sgnBeta*pow(t,beta); }

// tolerant comparison
bool isEqual(COST a, COST b) {
  double EPS = 1e-6;
  return ((a-EPS) / b <= 1+EPS &&
	  (b-EPS) / a <= 1+EPS);
}


// ------------------------------------- instance

int nbJobs =-1;
int p[64];
double w[64];

int optOrder[64];

int makespan=-1; // sum of all processing times

void initMakespan() {
  makespan = 0;
  for (int i=0; i<nbJobs; i++) 
    makespan += p[i];
}

/* ------------------------------------- job sets
   ...are represented as bit vectors.
   That's ok since the number of jobs will be small (<=40)  
   0 = the empty set
   allJobs = {0,1,...,nbJobs-1}
*/

typedef long long SET;

SET allJobs;

// input 2^i output i
int log2(SET I) {
  int i = 0;
  while (I>1) {
    i++;
    I/=2;
  }
  return i;
}


// iterator: use with brackets
#define forall(i,S) for(int i=0; i<nbJobs; i++) if (setContains(S,i))

// remove or add an element to a set
SET setToggle(SET S, int i) {return S ^ (1LL<<i);}

// test membership
bool setContains(SET S, int i) {return S & (1LL<<i);}
    
// produce a printable representation
string set2str(SET S) {
  ostringstream out;
  if (S==0) return "{}";
  char sep = '{';
  forall(i,S) {
    out << sep << i;
    sep = ',';
  }
  out <<  "}";
  return out.str();
}

void testSet() {
  if (nbJobs>62) { 
    cout << "ERROR: too many jobs for internal set datastructure" << endl;
    exit(1);
  }   
}

// ------------------------------------- prunning rules

COST tolerance = 1e-6;

// does i have to be before j if scheduled next to each other and completing at time t
bool localPrec(int i, int j, int t) {
  return 
    w[i]*penalty(t-p[j])+w[j]*penalty(t)  + tolerance
    < w[j]*penalty(t-p[i])+w[i]*penalty(t);
}

double timeChange(int i, int j) {
  double t0 = p[i]+p[j], t1=makespan;
  if (localPrec(i,j,t0)==localPrec(i,j,t1))
    return t1;
  while (t1-t0 > 0.01) { // limited precision
    // invariant localPrec(i,j,t0) != localPrec(i,j,t1)
    double m = (t0+t1)/2.;
    if (localPrec(i,j,t0)==localPrec(i,j,m))
      t0=m;
    else
      t1=m;
  }
  return t0;
}

double sepTime[64][64];
bool   sepOrd[64][64];

void initSepTime() {
  for (int i=0; i<nbJobs; i++)
    for (int j=0; j<nbJobs; j++) {
      sepTime[i][j] = timeChange(i,j);
      sepOrd[i][j]  = localPrec(i,j, sepTime[i][j]);
    } 
}

bool localPrecQuick(int i, int j, int t) {
  return sepOrd[i][j] ^ (t>sepTime[i][j]);
}



// in the interval [t0,t1] is i globaly before j?
bool globalPrec(int i, int j, int t0, int t1) {
  switch (rule) {
  case BASIC: 
      return (p[i]<p[j] && w[i]>=w[j]);
  case HOHNJACOBS:
      return (p[i]<p[j] && w[i]>=w[j])   ||
             (beta==2  && 
	            (p[i]<p[j] && w[i]*p[j] > tolerance + beta * w[j]*p[i] ||
	             p[i]>p[j] && w[i]*p[j] > tolerance + w[j]*p[i] ));
  case DURRVASQUEZ:
      return (p[i]<p[j] && w[i]>=w[j])   ||
             (beta==2 && ((p[i]<p[j] && localPrecQuick(i,j,t0)) || 
                          (p[i]>p[j] && w[i]*p[j] > tolerance + w[j]*p[i]) )) ||	
             (beta>=1 && ((p[i]<p[j] && localPrecQuick(i,j,t0) ) || 
                          (p[i]>p[j] && w[i]*pow(p[j],beta) > tolerance + w[j]*pow(p[i],beta)) )) ||
             (0<beta && beta<=1 && 
                         ((p[i]<p[j] && localPrecQuick(i,j,t1) ) || 
                          (p[i]>p[j] && w[i]*pow(p[j],2) > tolerance + w[j]*pow(p[i],2)) )) ||
	     (beta==-1 &&((p[i]<p[j] && localPrecQuick(i,j,t1) ) ||
                          (p[i]>p[j] && localPrecQuick(i,j,t0)) ));
  default:
      cerr << "ERROR: bad rule" << endl;
      return false;
  }
}

// ------------------------------------- conjecture

// is j placed before i, for a job pair for which i prec_local j?
void checkConjecture() {
  if (beta==2 || beta<=0) return; // these cases are already proven
  for (int a=0; a<nbJobs; a++)
    for (int b=a+2; b<nbJobs; b++) {   // try all non-adjacent pairs
      int j = optOrder[a];
      int i = optOrder[b];
      if (p[i]>p[j] && 
        ((beta>1 && w[i]*p[j]>=tolerance+w[j]*p[i] && w[i]*pow(p[j],beta)+tolerance <= w[j]*pow(p[i],beta)) ||
         (beta<1 && localPrec(i,j,p[i]+p[j]) && w[i]*pow(p[j],2) +tolerance<= w[j]*pow(p[i],2) ))) {
         cerr << "COUNTEREXAMPLE: i="<<i<<" j="<<j;
         for (int k=0; k<nbJobs; k++)
             cerr << " " << p[k] << "/" << w[k];
         cerr  << endl;
      }
    }
}

// ------------------------------------- produce explored DAG

ofstream graphFile;

// produce a DOT file readable with graphviz
void graphOpen() {
  char graphFileName[80];
  static int graphNumber=0;
  sprintf(graphFileName, "graph_%d.dot", graphNumber++);
  graphFile.open(graphFileName);
  graphFile << "digraph G {\n";
}

void graphVertex(SET v) {
  graphFile << v << " [label=\"" << set2str(v) << "\"]\n";
}
    
void graphArc(SET u, SET v, COST c) {
  graphFile << u << " -> " << v << "  [label=\"" << c << "\"]\n";
}

void graphClose() {
  graphFile << "}\n";
  graphFile.close();
}

// ------------------------------------- for the solvers

map<SET,int>  visitedFromJob;  // maps S to job j that lead to S in shortest path

/*
// the actual schedule
int order[64];
int tOrder;
bool compareByLocalPrec(int i, int j) {
  return localPrecEq(i,j,tOrder) && (!localPrecEq(j,i,tOrder) || i<j);
}

// t0 is the starting time from which on we need to schedule S
// does the local order induce a total order on the jobs S remaining to be scheduled?
bool isTrivial(SET S, int pS, int t0) {
  //  return false;
  forall(i,S)
    forall(j,S)
    if (i<j && localPrecEq(i,j,t0+p[i]+p[j]) != localPrecEq(i,j,t0+pS))
      return false;
  return true;
}

// solve the jobs according to this total order
void solveTrivial(SET S, int pS, int t0, double &opt, bool forward) {
  int n=0;
  forall(i, S)
    order[n++] = i;
  tOrder = t0+pS;
  sort(order, order+n, compareByLocalPrec);
  int t1 = tOrder;
  if (forward)
    for (int r=0; r<n; r++) {
      int j = order[r];
      t1 += p[j];
      opt += w[j] * penalty(t1);
      S = setToggle(S, j);
      visitedFromJob[S] = j;
  }
  else
    for (int r=n-1; r>=0; r--) {
      int j = order[r];
      opt += w[j] * penalty(t1);
      t1 -= p[j];
      S = setToggle(S, j);
      visitedFromJob[S] = j;
  }
}

string orderString(int t) {
  for (int i=0; i<nbJobs; i++)
    order[i] = i;
  tOrder = t;
  sort(order, order+nbJobs, compareByLocalPrec);
  ostringstream out;
  string sep = "[";
  for (int r=0; r<nbJobs; r++) {
    out << sep << order[r];
    sep = ",";
  }
  out << "]";
  return out.str();
}

void dumpVisitedFrom() {
  cout << "dumpVisisted:"<< endl;
  for (map<SET,int>::iterator i = visitedFromJob.begin();
       i!=visitedFromJob.end(); i++)
    cout << "visitedFromJob["<<i->first
	 << set2str(i->first)<<"]="<<i->second<<endl;
}

*/

// ------------------------------------- explore dag

struct Node {
  // c = cost of shortest path from root to this node 
  // cl = c + lower bound on distance to destination
  COST c, cl; 
  int j;  // vertex that lead to that node
  SET S;  // set of remaining jobs to be scheduled
  int pS; // total processing time of them
  Node(COST _c, COST _cl, int _j, SET _S, int _pS) {
    c=_c; cl=_cl; j=_j; S=_S; pS=_pS;
  }
  bool operator<(const Node &T) const {
    return (cl>T.cl) || (cl==T.cl) && 
      (S>T.S || S==T.S && j>=T.j);
  }
};

// returns optimal value of a schedule
COST solveBackward(int &nbNodes) {
  COST opt=-1;        // optimal cost, -1=not yet computed

  // queue contains (c,j,S) with c the cost of a path from source to S with last job j not in S
  priority_queue<Node> queue;
  //                                          -- prepare root node
  COST lowerBound   =0; 
  forall(i,allJobs) 
    lowerBound    += w[i]*penalty(p[i]);
  queue.push(Node(0, lowerBound, -1, allJobs, makespan));
  
#if DUMPGRAPH
  graphOpen();
#endif
  //                                          -- main loop
  while ( !queue.empty() ) {
    Node T = queue.top(); queue.pop();
    if (visitedFromJob.find(T.S) != visitedFromJob.end())
      continue;
    //                                        -- mark new visited vertex        
    visitedFromJob[T.S] = T.j;
    nbNodes += 1;
    if (nbNodes>=TIMEOUT)
      break;

#if DUMPGRAPH
    graphVertex(T.S);
#endif
    //                                        -- close to leaf ?
    if (T.S==0LL) {
      opt = T.c;
      break;
    }
    // if (isTrivial(T.S, T.pS, 0)) {
    //   solveTrivial(T.S, T.pS, 0, opt, false);
    //   break;
    // }
    //                                        -- prepare list of neighbors
    SET forbidden = 0LL;
    int t1;
    if (T.j!=-1) t1 = T.pS+p[T.j];
    forall(i, T.S) {
      if (T.j!=-1 && localPrec(T.j, i, t1))   // -- local rule
	forbidden = setToggle(forbidden, i);
      else {
	forall(k, T.S)                            // -- global rule
	  if (i!=k && globalPrec(i,k,p[i]+p[k], T.pS)) {
	    forbidden = setToggle(forbidden,i);
	    break;
	  }
      }
    }
    //                                        -- generate arcs to neighbors
    forall(j, T.S) {
      if (!setContains(forbidden, j)) {
	//          add real cost of j        remove lower bound
	COST c  = T.c  + w[j]*penalty(T.pS);
	COST cl = T.cl + w[j]*penalty(T.pS) - w[j]*penalty(p[j]);
	queue.push(Node(c, cl, j, setToggle(T.S, j), T.pS - p[j]));
#if DUMPGRAPH
	graphArc(T.S, setToggle(T.S, j), c);
#endif
      }
    }
  }
#if DUMPGRAPH
  graphClose();
#endif

  return opt;
}
// ------------------------------------- solve forward

// schedule all jobs from S right shifted. 
// basic lower bound over these schedules
COST lowerBound(SET S, int pS) {
  assert(makespan!=-1);
  COST lb=0;
  forall(i,S)
    lb += w[i] * penalty(makespan-pS+p[i]);
  return lb;
}

// returns optimal value of a schedule
COST solveForward(int &nbNodes) {
  COST opt=-1;        // optimal cost, -1=not yet computed

  // queue contains (c,j,S) with c the cost of a path from source to S with last job j not in S
  priority_queue<Node> queue;
  //                                          -- prepare root node
  queue.push(Node(0, lowerBound(allJobs, makespan), -1, allJobs, makespan));
  
#if DUMPGRAPH
  graphOpen();
#endif
  //                                          -- main loop
  while ( !queue.empty() ) {
    Node T = queue.top(); queue.pop();
    if (visitedFromJob.find(T.S) != visitedFromJob.end())
      //    if (visitedFromJob.count(T.S) > 0)
      continue;
    //                                        -- mark new visited vertex        
    visitedFromJob[T.S] = T.j;
    nbNodes += 1;
    if (nbNodes>=TIMEOUT)
      break;

#if DUMPGRAPH
    graphVertex(T.S);
#endif
    //                                        -- close to leaf ?
    if (T.S==0LL) {
      opt = T.c;
      break;
    }
    // if (isTrivial(T.S, T.pS, makespan-T.pS)) {
    //   solveTrivial(T.S, T.pS, makespan-T.pS, opt, true);
    //   break;
    // }
    //                                        -- prepare list of neighbors
    SET forbidden = 0LL;
    forall(i, T.S) {
      int t0 = makespan-T.pS+p[i];
      if (T.j!=-1 && localPrec(i, T.j, t0))   // -- local rule
	forbidden = setToggle(forbidden, i);
      else {
	forall(k, T.S)                            // -- global rule
	  if (i!=k && globalPrec(k,i,t0+p[k], makespan)) {
	    forbidden = setToggle(forbidden,i);
	    break;
	  }
      }
    }
    //                                        -- generate arcs to neighbors
    forall(j, T.S) {
      if (!setContains(forbidden, j)) {
	//          add real cost of j        remove lower bound
	int new_pS  = T.pS-p[j];
	SET new_S   = setToggle(T.S, j);
	COST c  = T.c + w[j]*penalty(makespan-new_pS);
	COST cl = c + lowerBound(new_S, new_pS);
	queue.push(Node(c, cl, j, new_S, new_pS));
#if DUMPGRAPH
	graphArc(T.S, setToggle(T.S, j), c);
#endif
      }
    }
  }
#if DUMPGRAPH
  graphClose();
#endif

  return opt;
}

// ------------------------------------- pretty print

void setOptOrder() {
  int n= ::forward ? nbJobs-1 : 0;
  SET S = 0;  //                                  -- extract optimal solution
  while (S!=allJobs) {
    if (visitedFromJob.find(S)==visitedFromJob.end()) {
       cout << "ERROR in setOptOrder: S=" << S << " n=" << n << endl;
       exit(1);
    }
    int j = visitedFromJob[S];
    assert (0<=n && n<nbJobs);
    if (::forward)
      optOrder[n--] = j;
    else
      optOrder[n++] = j;
    S = setToggle(S, j);
  }
}

// prints the optimal schedule in SVG format
void printSolution(double opt) {
  int horiz=max(2, 1400/makespan), vert=20;
  cout << "<svg xmlns=\"http://www.w3.org/2000/svg\" version=\"1.1\">"<< endl;
  //dumpVisitedFrom();

  int t=0;
  for (int r=0; r<nbJobs; r++) {
    int j = optOrder[r];
    cout << "<rect x='" << t*horiz << "' y='0'"  
	 << " width='" << int(p[j]*horiz) << "' height='"<<vert   
	 << "' style='fill:rgb(255,196,"<< j*255/nbJobs 
	 << ");stroke-width:1;stroke:black'/>" << endl 
	 << "<text x='" << t*horiz +2<< "' y='"<<vert-2<<"'>"<< j 
	 << "</text>" << endl;
    t += p[j];
  }
  int row=1;
  for (int j=0; j<nbJobs; j++)
    for (int i=0; i<j; i++) {
      int t = timeChange(i,j) - p[i] - p[j]; 
      // we shift because in the paper and the webpage, we use a different definition than for the prog.
      if (t>0) {	
	cout << "<line x1='"<< t*horiz << "' y1='" << vert 
	     << "' x2='"<< t*horiz << "' y2='"<< (row+1)*vert 
	     << "' style='stroke-width:1;stroke:grey'/>" << endl 
	     << "<text x='" << t*horiz +2 << "' y='" << (row+1)*vert-2
	     << "'>" << i << "," << j << ":" << t << "</text>" << endl;
	row++;
      }
    }
/*
  cout << "<text x='0'  y='" << (row+1)*vert << "'>" 
       << orderString(0) << "  " << orderString(makespan) 
       << "</text>" << endl;
  row++;
*/
  cout << "<text x='0'  y='" << (row+1)*vert << "'>" 
       << "objfct=sum wjCj^"<<beta<< " objvalue="<< opt 
       << "</text>" << endl;
  cout << "</svg>"<< endl;
}

// ------------------------------------- solver wrapper

COST solve(int &nbNodes) {
  double opt;
  allJobs = (1LL<<(long long)nbJobs)-1LL;

  initMakespan();
  initSepTime();
  visitedFromJob.clear();
  nbNodes = 0;

  if (direction==FORWARD || 
      direction==BEST && (rule==DURRVASQUEZ &&  beta>1 || 
                          rule!=DURRVASQUEZ && (beta<1 || beta==2))) {
    ::forward = true;
    opt = solveForward(nbNodes);
  }
  else {
    ::forward = false;
    opt = solveBackward(nbNodes);
  }
  if (nbNodes<TIMEOUT) {
    setOptOrder();
    checkConjecture();
  }
  return opt; 
}

/*
bool canSolve() {
  int nbNodes;
  return solve(nbNodes) != -1;
}
*/

// ------------------------------------- main program

bool compareSolutions = false;
bool statGlobalPrec   = false;

// command line given instance in pj/wj space separated format
void parseInstance(char *inst) {
    char *s = strtok(inst, " \t/");
    nbJobs = 0;
    while (s!=NULL) {
        p[nbJobs] = atoi(s);
        s = strtok(NULL, " \t/");
        if (p!=NULL) {
            w[nbJobs++] = atof(s);
        }
        s = strtok(NULL, " \t/");
    }
}

int parseArgs(char *argv[], int argc) {
  int argi=1;
  while (argi<argc && argv[argi][0]=='-') {
    if      (strncmp(argv[argi],"-beta=", 6)==0) {
      beta = atof(argv[argi]+6);
      sgnBeta = int(beta>0) - int(beta<0);
    }
    else if (strncmp(argv[argi],"-rule=", 6)==0) 
      rule = atoi(argv[argi]+6);
    else if (strcmp(argv[argi], "-globalPrec")==0)
      statGlobalPrec = true;
    else if (strcmp(argv[argi], "-compare")==0)
      compareSolutions = true;
    else if (strncmp(argv[argi], "-direction=",11)==0) 
      switch (argv[argi][11]) {
      case 'F':     direction = FORWARD;  break;
      case 'B':     direction = BACKWARD; break;
      default:      cout << "ERROR: wrong direction" << endl;
	exit(1);
      }
    else if (strncmp(argv[argi],"-instance=", 10)==0)
      parseInstance(argv[argi]+10);
    else {
      cout << "ERROR: wrong parameter "<<argv[argi] << endl;
      return -1;
    }
    argi++;
  } 
  return argi;
}

const char *optionsHelp = "Usage: AstarCompare <options> <instance files>\n"
	" options are \n"
	"     -beta=<beta>     [mandatory]\n"
	"     -rule=<c>        [which rule to use default:2=ours, 1=before]\n"
	"     -direction=<d>   [force resolution direction, F=forward, B=backward]\n"
	"     -compare         [to compare number of generated nodes with and without our rules]\n"
	"     -globalPrec      [just print fraction of job pairs with global precedence]\n"
	"     -instance=<inst> [instance is a space sep. list of pi/wi, and prog. prints the actual solution]\n";

int main(int argc, char *argv[]) {
  string s;
  double stdDeviation;
  int    nbInstances, k, argi;
  argi = parseArgs(argv,argc);
  if (argi==-1) {
    cout << optionsHelp;
    return 1;
  }
  testSet();
  if (nbJobs==-1) {
        cout << "beta=" << beta << " datafiles="<< argv[argi] << "..." << endl;
        cout << "beta\tsigma\tnbJobs\t#nodes\t#nodes-without-our-rules\n";
        for (; argi<argc; argi++) {
            ifstream in(argv[argi]);
            in >> s;
            if (s!="#") {
               cout << "ERROR: wrong file format or missing file " << argv[argi] << endl;
               continue;
            }
            // file format http://www.coga.tu-berlin.de/v-menue/projekte/complex_scheduling/generalized_min-sum_scheduling_instance_library/
            for (k=0; k<4; k++)  in >> s;
            in >> nbInstances;
            for (k=0; k<4; k++)  in >> s;
            in >> nbJobs;
            for (k=0; k<3; k++)  in >> s;
            in >> stdDeviation;
            for (k=0; k<11; k++) in >> s;
            for (int r=0; r<nbInstances; r++) {
                char c;
                for (int j=0; j<nbJobs; j++) 
                    in >> c >> c >> p[j] >> c >> w[j] >> c;
                in >> c;
		cout << beta << "\t" << stdDeviation << '\t' << nbJobs << '\t';
		if (compareSolutions) { 
                    int nb[2];
                    COST opt[2];
		    rule = DURRVASQUEZ;
		    opt[0] = solve(nb[0]);
		    rule = HOHNJACOBS;
		    opt[1] = solve(nb[1]);
                    cout << nb[0]  << '\t' <<  nb[1] << endl;
                    if (nb[0]<TIMEOUT && nb[1]<TIMEOUT && !isEqual(opt[0],opt[1]))
                        cout << "ERROR: Different objectives\t" 
                        << opt[0] << '\t' << opt[1] << endl;
                }
                else if (statGlobalPrec) {
                    int total=0, count=0, makespan=0;
		    for (int i=0; i<nbJobs; i++)
                        makespan += p[i];
		    for (int i=0; i<nbJobs; i++)
		       for (int j=i+1; j<nbJobs; j++) {
                           int t0 = p[i] + p[j];
                           total += 1;
                           if (globalPrec(i,j,t0,makespan) || globalPrec(j,i,t0,makespan))
                               count += 1;
                       }
                    cout << float(count) / total << endl;       
                }
                else {					
                    int nb;
                    solve(nb);
                    cout << nb << endl;
                }
            }
        }
  } else {
      int nb;
      printSolution(solve(nb));
  }
  return 0;
}