discretizationsolver.cpp

  model.add(ctaryLrefs);
  model.add(ctaryFullRef);
  //
  // set as not in clique, all the vertices that do not belong to one of the possible clique
  //        std::map < Vertex*, bool> okforclique;
  //        for (Vertex* v : this->inst_->vertices_) okforclique[v] = false;
  //        for (Clique* c : this->cliques_) {
  //            for (Vertex* v : c->vertices_) {
  //                okforclique[v] = true;
  //            }
  //        }
  //        for (Vertex* v : this->inst_->vertices_) {
  //            if (!okforclique[v]) {
  //                for (int k = 0; k <= L; k++) {
  //
  //                    model.add(hasrank_[v][k] == 0);
  //                }
  //            }
  //        }
}

//
// define the max fully-referenced cplex model

void DiscretizationSolver::defineminpartialip(IloModel & model) {

  IloEnv env = model.getEnv();
  char name[50];
  int nbvertices = this->inst_->nbvertices_;
  int L = this->inst_->L();
  int U = this->inst_->U();
  int nbcliques = this->cliques_.size();
  int nbcons = 0;

  //////////////////////////////////////////////////////////////////////////////
  // Declare the variables
  //////////////////////////////////////////////////////////////////////////////

  this->declareipvariables(model);


  //////////////////////////////////////////////////////////////////////////////
  // Set the cplex model that will be solved
  //////////////////////////////////////////////////////////////////////////////
  //
  // objective function
  IloExpr obj(env);
  obj += 1;
  for (Vertex* v : this->inst_->vertices_) {
    obj += 1 - isfullref_[v];
  }
  model.add(IloMinimize(env, obj));

  //////////////////////////////////////////////////////////////////////////////
  // Constraints that set the vertices belonging to the initial clique
  // Enumerate all the cliques and choose one
  //////////////////////////////////////////////////////////////////////////////
  //
  // choose one clique among those that have been enumerated
  if (this->modeltype_ != WITNESS) {
    IloRange ctOneClique(env, 1, 1, "ctOneClique");
    IloExpr sumcliques(env);
    for (int c = 0; c < nbcliques; c++) sumcliques += isclique_[c];
    ctOneClique = (sumcliques == 1);
    model.add(ctOneClique);

    if (this->modeltype_ == CCG) {
      for (int c = 0 ; c < nbcliques ; c++) {
        Clique* initialclique = this->cliques_[c];
        if (initialclique->greedyobjvalue_ >= 2) {
          IloExpr sumpartref(env);
          for (Vertex* v: initialclique->initialpartialrefs_) sumpartref += (1-isfullref_[v]);
          model.add(sumpartref >= isclique_[c]);
          for (Vertex* v: initialclique->initialfullrefs_) model.add(isfullref_[v] >= isclique_[c]);
        }
      }
    }
  }
  //
  // clique constraints in WITNESS: potential cliques are not enumerated
  if (this->modeltype_ == WITNESS) {
    //
    // vertices in the clique are not counted as part. ref
    IloRangeArray ctaryCliquesAreFull(env);
    nbcons = 0;
    for (Vertex* v : this->inst_->vertices_) {
      ctaryCliquesAreFull.add(isfullref_[v] - isvertexinclique_[v] >= 0);
      sprintf(name, "ctaryCliquesAreFull%i", v->id_);
      ctaryCliquesAreFull[nbcons++].setName(name);
    }
    model.add(ctaryCliquesAreFull);
    //
    // there must exactly L+1 vertices in the initial clique
    IloRange ctLPlusOneInClique(env, 1, 1, "ctLPlusOneInClique");
    IloExpr sumvertex(env);
    for (Vertex* v : this->inst_->vertices_) sumvertex += isvertexinclique_[v];
    ctLPlusOneInClique = (sumvertex == this->inst_->L() + 1);
    model.add(ctLPlusOneInClique);
    sumvertex.end();
    //
    // two non-adjacent vertices are not in the initial clique
    IloRangeArray ctaryCliquesAreAdjacent(env);
    nbcons = 0;
    for (Vertex* u : this->inst_->vertices_) {
      for (Vertex* v : this->inst_->vertices_) {
        if ((u == v) || (u->isneighbor(v))) continue;
        else {
          ctaryCliquesAreAdjacent.add(isvertexinclique_[u] + isvertexinclique_[v] <= 1);
          sprintf(name, "ctaryCliquesAreAdjacent%i_%i", u->id_, v->id_);
          ctaryCliquesAreAdjacent[nbcons++].setName(name);
        }
      }
    }
    model.add(ctaryCliquesAreAdjacent);
    //
    // vertices in the clique are witness to all their neighbors
    IloRangeArray ctaryCliqueIsWitness(env);
    nbcons = 0;
    for (Vertex* u : this->inst_->vertices_) {
      for (Vertex* v : u->neighbors_) {
        ctaryCliqueIsWitness.add(isbefore_[u][v] - isvertexinclique_[u] >= 0);
        sprintf(name, "ctaryCliqueIsWitness_%i_%i", u->id_, v->id_);
        ctaryCliqueIsWitness[nbcons++].setName(name);
      }
    }
    model.add(ctaryCliqueIsWitness);
    //
    // set as not in clique, all the vertices that do not belong to one of the possible clique
    //            std::map < Vertex*, bool> okforclique;
    //            for (Vertex* v : this->inst_->vertices_) okforclique[v] = false;
    //            for (Clique* c : this->cliques_) {
    //                for (Vertex* v : c->vertices_) {
    //                    okforclique[v] = true;
    //                }
    //            }
    //            for (Vertex* v : this->inst_->vertices_) {
    //                if (!okforclique[v]) {
    //                    model.add(isvertexinclique_[v] == 0);
    //                }
    //            }
  }

  //////////////////////////////////////////////////////////////////////////////
  // Constraints that guarantee that the solution digraph is acyclic
  // this is where the different models mostly differ:
  // - CYCLES forbid all the two-cycles and most three-cycles
  // - CCG, WITNESS: forbid only some cycles of the distance graph,  and
  // generate the other cycles through a callback
  // - RANKS: includes MTZ kinds of constraints instead of cycle
  // constraints, but add cycle cuts as valid inequalities anyways
  //////////////////////////////////////////////////////////////////////////////
  IloRangeArray ctaryNoTwoCycle(env);
  IloRangeArray ctaryNoThreeCycle(env);
  IloRangeArray ctaryNoFourCycle(env);
  //
  // in model CYCLES, we need to make sure that the solution is acyclic by
  // guaranteeing total order and transitivity
  if (this->modeltype_ == CYCLES) {
    nbcons = 0;
    for (Vertex* u : this->inst_->vertices_) {
      for (Vertex* v : this->inst_->vertices_) {
        if (v->id_ <= u->id_) continue;
        ctaryNoTwoCycle.add(isbefore_[u][v] + isbefore_[v][u] == 1);
        sprintf(name, "ctaryNoTwoCycle_%i_%i", u->id_, v->id_);
        ctaryNoTwoCycle[nbcons++].setName(name);
      }
    }
    model.add(ctaryNoTwoCycle);
    //
    // forbid every three-cycle but those that have no common edge with
    // the distance graph (we proved that this is sufficient)
    nbcons = 0;
    for (Vertex* u : this->inst_->vertices_) {
      for (Vertex* v : u->neighbors_) {
        if (v->id_ <= u->id_) continue;
        for (Vertex* w : this->inst_->vertices_) {
          if (w->id_ <= v->id_) continue;
          ctaryNoThreeCycle.add(
              isbefore_[u][v] + isbefore_[v][w] + isbefore_[w][u]
                  <= 2);
          sprintf(name, "ctaryNoThreeCycle_%i_%i_%i", u->id_, v->id_, w->id_);
          ctaryNoThreeCycle[nbcons++].setName(name);
          ctaryNoThreeCycle.add(
              isbefore_[w][v] + isbefore_[v][u] + isbefore_[u][w]
                  <= 2);
          sprintf(name, "ctaryNoThreeCycle_%i_%i_%i",
                  w->id_, v->id_, u->id_);
          ctaryNoThreeCycle[nbcons++].setName(name);
        }
      }
    }
    model.add(ctaryNoThreeCycle);
  }
  //
  // initially forbid cycles of the distance graphs with size <=
  // initialcyclesize_ in models CCG, WITNESS and RANKS;
  // in CCG and WITNESS; the remaining cycles cuts are added as lazy
  // constraints
  if ((this->modeltype_ == CCG || this->modeltype_ == WITNESS || this->modeltype_ == RANKS) && (this->initialcyclesize_ >= 2)) {
    // - forbid every two-cycle
    // - only the edges of the distance graph need to be considered in
    // the two models CCG and RANKS
    nbcons = 0;
    for (Vertex* u : this->inst_->vertices_) {
      for (Vertex* v : u->neighbors_) {
        if (v->id_ <= u->id_) continue;

        if (this->modeltype_ == WITNESS) {
          ctaryNoTwoCycle.add(
              isbefore_[u][v] + isbefore_[v][u] - isvertexinclique_[u]
                  <= 1);
        } else {
          ctaryNoTwoCycle.add(isbefore_[u][v] + isbefore_[v][u] == 1);
        }
        sprintf(name, "ctaryNoTwoCycle_%i_%i", u->id_, v->id_);
        ctaryNoTwoCycle[nbcons++].setName(name);
      }
    }
    model.add(ctaryNoTwoCycle);

    // forbid the three-cycles of the distance graph
    if (this->initialcyclesize_ >= 3) {
      nbcons = 0;
      for (Vertex* u : this->inst_->vertices_) {
        for (Vertex* v : u->neighbors_) {
          if (v->id_ <= u->id_) continue;
          for (Vertex* w : v->neighbors_) {
            if (!w->isneighbor(u)) continue;
            if (w->id_ <= v->id_) continue;
            if (this->modeltype_ == WITNESS) {
              ctaryNoThreeCycle.add(
                  isbefore_[u][v] + isbefore_[v][w]
                      + isbefore_[w][u] - isvertexinclique_[u] <= 2);
              ctaryNoThreeCycle.add(
                  isbefore_[w][v] + isbefore_[v][u]
                      + isbefore_[u][w] - isvertexinclique_[u] <= 2);
            } else {
              ctaryNoThreeCycle.add(
                  isbefore_[u][v] + isbefore_[v][w] + isbefore_[w][u] <= 2);
              ctaryNoThreeCycle.add(
                  isbefore_[w][v] + isbefore_[v][u] + isbefore_[u][w] <= 2);

            }
            sprintf(name, "ctaryNoThreeCycle_%i_%i_%i",
                    u->id_, v->id_, w->id_);
            ctaryNoThreeCycle[nbcons++].setName(name);
            sprintf(name, "ctaryNoThreeCycle_%i_%i_%i",
                    w->id_, v->id_, u->id_);
            ctaryNoThreeCycle[nbcons++].setName(name);
          }
        }
      }
      model.add(ctaryNoThreeCycle);
    }
    //
    // forbid every four-cycle: the complex verifications aim at adding only
    // one constraint per four-cycle that is not implied by two three-cycles
    if (this->initialcyclesize_ >= 4) {
      nbcons = 0;
      for (Vertex* u : this->inst_->vertices_) {
        for (Vertex* v : u->neighbors_) {
          if (v->id_ <= u->id_) continue;
          for (Vertex* w : v->neighbors_) {
            for (Vertex* x : w->neighbors_) {
              if (!x->isneighbor(u)) continue;
              if (w->id_ <= x->id_ && w->id_ <= v->id_) continue;
              if (x->isneighbor(v) && w->isneighbor(u)) continue;

              if (this->modeltype_ == WITNESS) {
                ctaryNoFourCycle.add(
                    isbefore_[u][v] + isbefore_[v][w] + isbefore_[w][x]
                        + isbefore_[x][u] - isvertexinclique_[u] <= 3);
                ctaryNoFourCycle.add(
                    isbefore_[x][w] + isbefore_[w][v] + isbefore_[v][u]
                        + isbefore_[u][x] - isvertexinclique_[u] <= 3);
              } else {
                ctaryNoFourCycle.add(
                    isbefore_[u][v] + isbefore_[v][w] + isbefore_[w][x]
                        + isbefore_[x][u] <= 3);
                ctaryNoFourCycle.add(
                    isbefore_[x][w] + isbefore_[w][v] + isbefore_[v][u]
                        + isbefore_[u][x] <= 3);

              }
              sprintf(name, "ctaryNoFourCycle_%i_%i_%i_%i",
                      u->id_, v->id_, w->id_, x->id_);
              ctaryNoFourCycle[nbcons++].setName(name);
              sprintf(name, "ctaryNoFourCycle_%i_%i_%i_%i",
                      x->id_, w->id_, v->id_, u->id_);
              ctaryNoFourCycle[nbcons++].setName(name);
            }
          }
        }
      }
      model.add(ctaryNoFourCycle);
    }
  }

  // - RANKS: first, make sure that each rank is taken by exaclty one vertex
  //	then, use MTZ constraints instead of cycle constraints
  if (this->modeltype_ == RANKS) {
    //
    // each rank of the order is taken by exactly one vertex
    IloRangeArray ctaryOneVertexPerRank(env);
    for (int k = 0; k < nbvertices; k++) {
      IloExpr sumvertices(env);
      for (Vertex* v : this->inst_->vertices_) sumvertices += hasrank_[v][k];
      ctaryOneVertexPerRank.add(sumvertices == 1);
      sprintf(name, "ctaryOneVertexPerRank_%i", k);
      ctaryOneVertexPerRank[k].setName(name);
    }
    model.add(ctaryOneVertexPerRank);
    //
    // each vertex has exactly one rank and his rank is computed
    IloRangeArray ctaryOneRankPerVertex(env);
    IloRangeArray ctaryComputeRank(env);
    nbcons = 0;
    for (Vertex* u : this->inst_->vertices_) {
      IloExpr sumhasrank(env);
      IloExpr sumranks(env);

      for (int k = 0; k < nbvertices; k++) sumhasrank += hasrank_[u][k];
      for (int k = 0; k < nbvertices; k++) sumranks += k * hasrank_[u][k];
      ctaryOneRankPerVertex.add(sumhasrank == 1);
      ctaryComputeRank.add(sumranks - rank_[u] == 0);
      sprintf(name, "ctaryOneRankPerVertex_%i", u->id_);
      ctaryOneRankPerVertex[nbcons].setName(name);
      sprintf(name, "ctaryComputeRank_%i", u->id_);
      ctaryComputeRank[nbcons++].setName(name);
    }
    model.add(ctaryOneRankPerVertex);
    model.add(ctaryComputeRank);
    //
    // Miller-Tucker-Zemlin constraints to ensure that the direction of the
    // distance graph edges respect the order of the ranks
    IloRangeArray ctaryRankTransitivity(env);
    nbcons = 0;
    for (Vertex* u : this->inst_->vertices_) {
      for (Vertex* v : u->neighbors_) {
        if (v->id_ <= u->id_) continue;
        ctaryRankTransitivity.add(
            1 <= nbvertices * isbefore_[u][v] - rank_[v] + rank_[u]
                <= nbvertices - 1);
        sprintf(name, "ctaryRankTransitivity_%i_%i", u->id_, v->id_);
        ctaryRankTransitivity[nbcons++].setName(name);
      }
    }
    model.add(ctaryRankTransitivity);
  }


  //////////////////////////////////////////////////////////////////////////////
  // Constraints that ensure that the ordering is a discretization ordering
  // every vertex has at least L references except if it is a member of the initial clique
  //////////////////////////////////////////////////////////////////////////////

  nbcons = 0;
  int nbrefcons = 0;
  IloRangeArray ctaryLrefs(env);
  //
  // in WITNESS, the constraint reflects the difference between witnesses and references
  if (this->modeltype_ == WITNESS) {
    for (Vertex* u : this->inst_->vertices_) {
      IloExpr sumwitness(env);
      for (Vertex* v : u->neighbors_) sumwitness += isbefore_[v][u];
      ctaryLrefs.add(
          sumwitness + (U - L) * (isvertexinclique_[u] - isfullref_[u]) == L);
      sprintf(name, "ctaryLRefs_%i", u->id_);
      ctaryLrefs[nbrefcons++].setName(name);
      sumwitness.end();
    }
  } else {
    for (Vertex* u : this->inst_->vertices_) {
      IloExpr sumrefs(env);
      IloExpr sumcliques(env);
      for (int c : this->cliqueidslist_[u]) {
        sumcliques += isclique_[c];
      }
      for (Vertex* v : u->neighbors_) sumrefs += isbefore_[v][u];
      ctaryLrefs.add(sumrefs + U * sumcliques - (U - L) * isfullref_[u] >= L);
      sprintf(name, "ctaryLRefs_%i", u->id_);
      ctaryLrefs[nbrefcons++].setName(name);
      sumrefs.end();
      sumcliques.end();
    }
  }
  model.add(ctaryLrefs);


  if (this->modeltype_ == RANKS) {
    IloRangeArray ctaryRankOfClique(env);
    IloExprArray sumrank(env, L + 1);
    for (int k = 0; k < L + 1; k++) {
      sumrank[k] = IloExpr(env);
    }
    for (Vertex* u : this->inst_->vertices_) {
      for (int c : this->cliqueidslist_[u]) {
        for (int k = 0; k < L + 1; k++) {
          if (this->cliques_[c]->vertices_[k] == u) sumrank[k] += isclique_[c];
        }
      }
      for (int k = 0; k < L + 1; k++) {
        ctaryRankOfClique.add(hasrank_[u][k] - sumrank[k] == 0);
        sprintf(name, "ctaryRankOfClique_%i_%i", u->id_, k);
        ctaryRankOfClique[nbcons++].setName(name);
      }
    }
  }

  //////////////////////////////////////////////////////////////////////////////
  // Constraints that strengthen the formulation
  //////////////////////////////////////////////////////////////////////////////
  //
  // if we could identify cycles composed of low degree vertices add a cut
  // specifying that the vertices included in such cycles must include at least
  // as many partially referenced vertices as there are cycles
  if (this->modeltype_ == CCG) {
    IloRangeArray ctaryLowDegreeCycles(env);
    for (int i = 0; i < this->lowdegreecycles_.size(); i++) {
      std::vector<Vertex*> cycle = this->lowdegreecycles_[i];
      int cyclesize = cycle.size();
      IloExpr sumfullrefincycle(env);
      IloExpr sumcliques(env);
      for (Vertex* v : cycle) {
        sumfullrefincycle += isfullref_[v];
        for (int c : this->cliqueidslist_[v]) {
          sumcliques += isclique_[c];
        }
      }
      ctaryLowDegreeCycles.add(
          sumfullrefincycle - sumcliques
              <= cyclesize - this->nbpartialsincycle_[i]);
      sprintf(name, "ctaryLowDegreeCycles_%i", i);
      ctaryLowDegreeCycles[i].setName(name);
    }
    model.add(ctaryLowDegreeCycles);
    //
    // if we could identify cliques that contain at least one
    // partially-referenced vertex, add a constraint to enforce this;  the
    // two-cliques need not be checked
    IloRangeArray ctaryCliquesWithPartialRef(env);
    nbcons = 0;
    for (Clique* clique : this->cliqueswithpartialref_) {
      IloExpr sumfullref(env);
      IloExpr sumcliques(env);
      for (Vertex* v : clique->vertices_) {
        sumfullref += isfullref_[v];
        for (int c : this->cliqueidslist_[v]) {
          sumcliques += isclique_[c];
        }
      }
      ctaryCliquesWithPartialRef.add(
          sumfullref - sumcliques <=
          clique->nbvertices() - this->nbpartialinclique_[clique]);
      sprintf(name, "ctaryCliquesWithPartialRef_%i", nbcons);
      ctaryCliquesWithPartialRef[nbcons++].setName(name);
    }
    model.add(ctaryCliquesWithPartialRef);
  }
}
//
// Integer programming model for the search of a discretization order that
// maximizes the number of vertices with more than free references
// modeltype:
// = 0 :	variables that indicate whether a vertex is before another one for
// 	every pair of vertices
// = 1 : variables only for each edge, but no enumeration of the cliques
// = 2 :	variables only for each edge, with enumeration of the cliques

bool DiscretizationSolver::minpartialip(float timelimit) {

  IloEnv env;
  IloTimer cpuClockTotal(env);
  IloTimer cpuClockInit(env);
  cpuClockTotal.start();
  cpuClockInit.start();

  this->inst_->computeneighbors();
  //    this->preallocatelowdegreevertices();
  //
  // enumerate the potential initial cliques
  if (!this->enumeratecliques(this->inst_->L() + 1)) {
    this->isfeasible_ = false;
    return false;
  }
  this->eliminateredundantcliques();
  //
  // solve with greedy and remove every clique that cannot be completed
  if (!this->greedysolve()) {
    std::cout << "Greedy solve found the problem to be infeasible at the"
                 " time of mip warm start." << std::endl;
    return false;
  }
  if (this->modeltype_ == CCG) {
    //
    // identify clique of vertices that will contain at least one
    // partially-referenced vertex
    this->computecliquecuts(this->cliquecutsmaxsize_);
    //
    // identify cycle of vertices with degree smaller than U+1: in each
    // of these cycles, at least one vertex is partially referenced
    this->enumeratelowdegreecycles();
  }

  IloModel model(env);
  if (this->modeltype_ == VERTEXRANK) {
    this->defineminpartial_vertexrank(model);
  } else {
    this->defineminpartialip(model);
  }

  //
  // create CPLEX object, load the model and provide initial solution
  IloCplex cplex(env);
  cplex.extract(model);
  //
  if (!this->greedymipstart(cplex)) {
    this->isfeasible_ = false;
    return false;
  }
  if ((this->modeltype_ == CCG) || (this->modeltype_ == WITNESS)) {
    cplex.use(CyclesLazyConstraints(env, *(this->inst_), *this));
  }
  if (this->modeltype_ == CCG) {
    IloCplex::MIPCallbackI::NodeId nodeid;
    cplex.use(BranchOnCliqueVariables(env, *this, nodeid));
  }

  cpuClockInit.stop();

  //////////////////////////////////////////////////////////////////////////////
  // Solve the instance
  //////////////////////////////////////////////////////////////////////////////
  //
  // set CPLEX parameters
  cplex.setParam(IloCplex::MIPDisplay, 4);
  cplex.setParam(IloCplex::TiLim, timelimit - cpuClockInit.getTime());
  cplex.setParam(IloCplex::Threads, 1);


#ifdef VERBOSE
  cplex.exportModel("model.lp");
#endif
  //
  // solve the integer program
  std::cout << "\nrevorder: ----------------------------------------------------------" << std::endl;
  std::cout << "revorder: solve an integer program to find an optimal discretization order" << std::endl;
  std::cout << "revorder: time limit has been set to " << timelimit << " seconds" << std::endl;
  std::cout << "revorder: cplex log coming below\n" << std::endl;
  try {
    cplex.solve();
  } catch (IloException& e) {
    std::cerr << e << std::endl;
  }
  cpuClockTotal.stop();


  IloAlgorithm::Status status = cplex.getStatus();
  std::cout << "revorder: ----------------------------------------------------------\n" << std::endl;
  std::cout << "revorder: CPLEX solution status = " << status << std::endl;
  isfeasible_ = (status == IloAlgorithm::Feasible) || (status == IloAlgorithm::Optimal);
  isoptimal_ = (status == IloAlgorithm::Optimal);

  if (isfeasible_) {
    //
    // retrieve the information on the solution process
    objvalue_ = std::min((int) cplex.getObjValue(), objvalue_);
    relativegap_ = cplex.getMIPRelativeGap();
    totaltime_ = cpuClockTotal.getTime();
    treatednodes_ = cplex.getNnodes();
#ifdef VERBOSE
    cplex.writeSolution("solution.sol");
#endif

    std::cout << "revorder: total cpu time                    = " << cpuClockTotal.getTime() << " s" << std::endl;
    std::cout << "revorder: initialization cpu time           = " << cpuClockInit.getTime() << " s" << std::endl;
    std::cout << "revorder: number of branching nodes         = " << treatednodes_ << std::endl;
    std::cout << "revorder: value of the objective function   = " << objvalue_ << std::endl;
    std::cout << std::setprecision(1) << "revorder: value of the relative duality gap = " << 100.0 * relativegap_ << "%" << std::endl;
  } else if (!isfeasible_) {
    std::string fileInfeasible("InfeasibleModel.lp");
    cplex.exportModel(fileInfeasible.c_str());
    std::cout << "revorder: no feasible solution was found during the optimization!" << std::endl;
    std::cout << "revorder: check the file " << fileInfeasible << " to see the corresponding model" << std::endl;
    cplex.end();
    env.end();
    return false;
  }
  //
  // build a discretization order that has the same value as the solution of
  // the integer program
  objvalue_ += this->inst_->preallocatedvertices_.size();
  this->buildoptimalorder(cplex);
  std::cout << "revorder: verification: number of part. ref. vertices = " << this->inst_->nbvertices_ - this->inst_->L() - this->bestnbfullref_ << std::endl;
  std::cout << "revorder: with preallocated vertices, we get " << objvalue_ << " part. ref. vertices" << std::endl;
  // terminate the cplex objects
  cplex.end();
  env.end();

  return isfeasible_;
}

//
// Declare the variables used in the ip models

void DiscretizationSolver::declareipvariables(IloModel & model) {

  IloEnv env = model.getEnv();
  char name[50];
  int nbvertices = this->inst_->nbvertices_;
  int L = this->inst_->L();
  int nbcliques = this->cliques_.size();

  //
  // variables that represent the order relations between the vertices
  // - isbefore_[i,j] = 1 if vertex i is before j in the order
  for (Vertex* u : this->inst_->vertices_) {
    for (Vertex* v : this->inst_->vertices_) {
      if (u == v) continue;
      isbefore_[u][v] = IloBoolVar(env);
      sprintf(name, "isbefore_%i_%i", u->id_, v->id_);
      isbefore_[u][v].setName(name);
    }
  }
  //
  // isfullref_[v] =1 if vertex v is fully-referenced and 0 otherwise
  for (Vertex* v : this->inst_->vertices_) {
    isfullref_[v] = IloBoolVar(env);
    sprintf(name, "isfullref_%i", v->id_);
    isfullref_[v].setName(name);
  }
  //
  // variables that set the initial clique
  // - variables that set the initial clique from the enumeration of the cliques
  isclique_ = IloBoolVarArray(env, nbcliques);
  for (int c = 0; c < nbcliques; c++) {
    sprintf(name, "isclique_%i", c);
    isclique_[c].setName(name);
  }
  // - when no enumeration, isvertexinclique_[v]=1 if v is in the initial clique, 0 otherwise
  for (Vertex* v : this->inst_->vertices_) {
    isvertexinclique_[v] = IloBoolVar(env);
    sprintf(name, "isvertexinclique_%i", v->id_);
    isvertexinclique_[v].setName(name);
  }

  //
  // variables that set the rank of each vertex in the order
  // - hasrank_[v,k] = 1 if vertex i has rank k in the discretization order
  //  in this model, we use these variables only for the first L vertices
  //  of the cliques (if modeltype <= 1)
  int nbranks = (this->modeltype_ == RANKS) ? nbvertices : L + 1;
  for (Vertex* v : this->inst_->vertices_) {
    hasrank_[v] = IloBoolVarArray(env, nbranks);
    for (int k = 0; k < nbranks; k++) {
      sprintf(name, "hasrank_%i_%i", v->id_, k);
      hasrank_[v][k].setName(name);
    }
  }
  //
  // - rank[v] = integer variable representing the rank of v in the order
  if (this->modeltype_ == RANKS) {
    for (Vertex* v : this->inst_->vertices_) {

      rank_[v] = IloIntVar(env, 0, nbvertices - 1);
      sprintf(name, "rank_%i", v->id_);
      rank_[v].setName(name);
    }
  }
}
//
// create the cplex model for the relaxation of the IP

void DiscretizationSolver::createrelaxationmodel(IloModel& model, IloModel & relax) {
  IloEnv env = model.getEnv();
  relax.add(model);
  for (Vertex* u : this->inst_->vertices_) {
    relax.add(IloConversion(env, isfullref_[u], ILOFLOAT));
    if (this->modeltype_ == CCG || this->modeltype_ == WITNESS || this->modeltype_ == RANKS) {
      for (Vertex* v : u->neighbors_) {
        if (v == u) continue;
        relax.add(IloConversion(env, isbefore_[u][v], ILOFLOAT));
      }
    } else {
      for (Vertex* v : this->inst_->vertices_) {
        if (v == u) continue;
        relax.add(IloConversion(env, isbefore_[u][v], ILOFLOAT));
      }
    }
    for (int k = 0; k < this->inst_->L()+1; k++) {
      relax.add(IloConversion(env, hasrank_[u][k], ILOFLOAT));
    }
    if (this->modeltype_ == RANKS) {
      for (int k = this->inst_->L()+1; k < this->inst_->nbvertices_; k++) {
        relax.add(IloConversion(env, hasrank_[u][k], ILOFLOAT));
      }
    }
  }
  for (int n = 0; n < this->cliques_.size(); n++) {

    relax.add(IloConversion(env, isclique_[n], ILOFLOAT));
  }
}
//
// Solve the linear relaxation of the input model

void DiscretizationSolver::solverelaxation(IloModel & model) {
  //
  // solve the relaxation alone first and disable every generation of cuts and
  // presolve of the integer programming problem

  IloEnv env = model.getEnv();
  IloModel relax(env);
  createrelaxationmodel(model, relax);
  IloCplex cplexrelax(env);

  cplexrelax.extract(relax);
  //
  // run cplex and retrieve the objective value
  std::cout << "\nrevorder: ----------------------------------------------------------" << std::endl;
  std::cout << "revorder: solve the relaxation of the discretization order problem" << std::endl;
  std::cout << "revorder: cplex log coming below\n" << std::endl;
  IloTimer cpuClock(env);
  cpuClock.start();
  cplexrelax.solve();
  cpuClock.stop();
  objvaluerelax_ = cplexrelax.getObjValue();
  std::cout << "revorder: value of the linear relaxation  = " << objvaluerelax_ << std::endl;
  std::cout << "revorder: total cplex cpu (relaxation) = " << cpuClock.getTime() << " s" << std::endl;
#ifdef VERBOSE
  cplexrelax.writeSolution("solutionrelax.sol");
#endif
  //
  cplexrelax.end();
}
//
// After solving the problem, build an order that produces the same objective
// value as the optimal solution

void DiscretizationSolver::buildoptimalorder(const IloCplex & cplex) {
  int nbvertices = this->inst_->nbvertices_;
  //
  // initialize the data relative to best solution
  this->bestfullref_.clear();
  for (Vertex* v : this->inst_->vertices_) {
    this->bestnbfullref_ = 0;
    this->bestnbrefs_[v] = 0;
    this->bestisfullref_[v] = false;
    this->bestrank_[v] = -1;
  }
  //
  // Compute a topological ordering of the graph
  //
  // generate the adjacency lists of the vertices corresponding to the current
  // cplex solution
  if (this->modeltype_ != VERTEXRANK) {
    float tolerance = 1e-06;
    std::map<Vertex*, std::vector<Vertex*> > adjlist;
    for (Vertex* v : this->inst_->vertices_) {
      adjlist[v] = std::vector<Vertex*>();
    }

    for (Vertex* u : this->inst_->vertices_) {
      for (Vertex* v : u->neighbors_) {
        if (cplex.getValue(this->isbefore_[u][v]) >= 1 - tolerance) {
          if ((this->modeltype_ == WITNESS) && (cplex.getValue(this->isvertexinclique(v)) >= 1 - tolerance)) {
            if (cplex.getValue(this->isvertexinclique(u)) >= 1 - tolerance) {
              if (u->id_ < v->id_) {
                adjlist[u].push_back(v);
                this->bestnbrefs_[v]++;
              }
            }
          } else {
            adjlist[u].push_back(v);
            this->bestnbrefs_[v]++;
          }
        }
      }
    }
    //
    // search for the root of the digraph
    Vertex* root = nullptr;
    int minnbrefs = this->inst_->nbvertices_;
    for (Vertex* v : this->inst_->vertices_) {
      if (this->bestnbrefs_[v] < minnbrefs) {
        minnbrefs = this->bestnbrefs_[v];
        root = v;
        if (minnbrefs == 0) break;
      }
    }
    //
    // call the enumeration of cycles, returns false if there are none
    std::vector<std::vector<Vertex*> > cycles;
    // Search for a topological order that will be valid only if the
    // digraph is acyclic
    std::vector<Vertex*> reverseorder;
    std::map<Vertex*, int> rank;
    this->topologicalorder(root, adjlist, reverseorder, rank);

    // determine whether the digraph is cyclic or not
    //
    std::vector<std::pair<Vertex*, Vertex*> > reverseedges;
    bool iscyclic = this->getreverseedges(adjlist, reverseorder, rank, reverseedges);
    if (iscyclic) {
      std::cout << "revorder: error: the final integer solution is cyclic" << std::endl;
      throw;
    }


    // modify the rank attributes of the vertices to represent their rank in the order
    for (Vertex* v : this->inst_->vertices_) {
      this->bestrank_[v] = rank[v];
      if (this->bestnbrefs_[v] >= this->inst_->U()) {
        this->bestisfullref_[v] = true;
        this->bestfullref_.push_back(v);
        this->bestnbfullref_++;
      }
    }
  } else {
    // get the rank of each vertex in the order
    for (Vertex* u : this->inst_->vertices_) {
      for (int k = 0; k < nbvertices; k++) {
        if (cplex.getValue(this->hasrank_[u][k]) >= 1 - 1.0e-6) {
          this->bestrank_[u] = k;
        }
      }
    }
    // get the number of references of each vertex
    for (Vertex* u : this->inst_->vertices_) {
      for (Vertex* v : u->neighbors_) {
        if (v->id_ <= u->id_) continue;
        if (this->bestrank_[v] < this->bestrank_[u]) this->bestnbrefs_[u]++;
        else this->bestnbrefs_[v]++;
      }
    }
    // get fully referenced vertices
    for (Vertex* u : this->inst_->vertices_) {
      if (this->bestnbrefs_[u] >= this->inst_->U()) {
        this->bestisfullref_[u] = true;
        this->bestfullref_.push_back(u);
        this->bestnbfullref_++;
      }
    }
  }
  //
  // verify that the order is indeed a referenced order
  this->reconstructsolution();
  // this->verifyorder(this->bestrank_);
  //
  // sort the input ranks map by ascening order of values
  std::set<std::pair<Vertex*, int>, Comparator> orderedranks(
      this->bestrank_.begin(), this->bestrank_.end(), valuecompare);
}