discretizationsolver.cpp

/********************************************************************************************
 Name:       DiscretizationSolver
 Discrete optimization for graph discretization
 Author:     J.Omer
 Sources:    C++
 License:    GNU General Public License v.2
 History:
 *********************************************************************************************/

#include "discretizationsolver.hpp"
#include <ctime>



//

ILOLAZYCONSTRAINTCALLBACK2(CyclesLazyConstraints, Instance &, inst, DiscretizationSolver&, solver) {

#ifdef VERBOSE
std::cout << "revorder: check lazy dicycle constraints" << std::endl;
#endif

// Search for cycles in the current solution digraph
//
// generate the adjacency lists of the vertices corresponding to the current
// cplex solution
float tolerance = 1e-06;
std::map<Vertex*, std::vector<Vertex*> > adjlist;
std::map<Vertex*, int > nbrefs;
for (Vertex* v : inst.vertices_) {
adjlist[v] = std::vector<Vertex*>();
nbrefs[v] = 0;
}
for (Vertex* u : inst.vertices_) {
if (solver.modeltype_ == WITNESS) {
if (getValue(solver.isvertexinclique(u) >= 1 - tolerance)) continue;
}
for (Vertex* v : u->neighbors_) {
if (getValue(solver.isbefore(u, v)) >= 1 - tolerance) {
adjlist[u].push_back(v);
nbrefs[v]++;
}
}
}
//
// search for the root of the digraph
Vertex* root = nullptr;
int minnbrefs = inst.nbvertices_;
if (solver.modeltype_ == WITNESS) {
for (Vertex* v : inst.vertices_) {
if (getValue(solver.isvertexinclique(v) >= 1 - tolerance)) {
root = v;
break;
}
}
} else {
for (Vertex* v : inst.vertices_) {
if (nbrefs[v] < minnbrefs) {
minnbrefs = nbrefs[v];
root = v;
if (minnbrefs < inst.L()) break;
}
}
}
//
// call the enumeration of cycles, returns false if there are none
std::vector<std::vector<Vertex*> > cycles;
bool iscyclic = solver.enumeratecycles(adjlist, root, cycles);

// Look for dicycle inequalities : we restrict the search to the dicycles
// that contain at least one edge in the distance graph
// besides that, the search is performed using brute force
int nbaddedcuts = 0;
IloEnv env = getEnv();
if (iscyclic) {
for (std::vector<Vertex*> path : cycles) {
IloExpr sumedges(env);
int cyclelength = path.size();
Vertex* u = path.back();
for (int i = 0; i < cyclelength; i++) {
Vertex* v = path[i];
sumedges += solver.isbefore(u, v);
u = v;
}
if ((solver.modeltype_ == WITNESS) && (cyclelength <= inst.L() + 1)) {
add(sumedges - solver.isvertexinclique(path.front()) <= cyclelength - 1).end();
nbaddedcuts++;
} else {
add(sumedges <= cyclelength - 1).end();
nbaddedcuts++;
}
}
#ifdef VERBOSE
std::cout << "revorder: number of violated dicycle constraints : " << nbaddedcuts << std::endl;
#endif
} else {
#ifdef VERBOSE
std::cout << "revorder: the graph is acyclic " << std::endl;
#endif
}
return;
}

// User branch callback that gives priority to branching on clique variables

ILOBRANCHCALLBACK2(BranchOnCliqueVariables, DiscretizationSolver&, solver, IloCplex::MIPCallbackI::NodeId&, nodeId) {
if (getBranchType() != BranchOnVariable)
return;

if (nodeId == getNodeId())
return;
nodeId = getNodeId();

// Branch on the fractionnary clique variable closest to 1
IntegerFeasibilityArray feas;
IloNumArray val;
try {
val = IloNumArray(getEnv());
feas = IntegerFeasibilityArray(getEnv());
getValues(val, solver.isclique_);
getFeasibilities(feas, solver.isclique_);
//
// get the clique variable with maximum integer infeasibility
IloInt bestvar = -1;
IloNum maxval = 0.3;
IloInt cols = solver.isclique_.getSize();
for (IloInt j = 0; j < cols; j++) {
if (feas[j] == Infeasible) {
if (val[j] >= 1 - 1.0e-6) {
bestvar = -1;
break;
}
if (std::min(val[j], 1 - val[j]) > maxval) {
bestvar = j;
maxval = std::min(val[j], 1 - val[j]);
}
} else if (val[j] >= 1 - 1.0e-6) {
bestvar = -1;
break;
}
}
//
// if there is at least one fractionnary clique variable create two branches
if (bestvar >= 0) {
//            HeuristicNode* data = new HeuristicNode();
//            data->fixcliqueindex(bestvar);
makeBranch(solver.isclique_[bestvar], val[bestvar], IloCplex::BranchUp, getObjValue());
makeBranch(solver.isclique_[bestvar], val[bestvar], IloCplex::BranchDown, getObjValue());
#ifdef VERBOSE
std::cout << "revorder: clique " << bestvar << " bounds " << getLB(solver.isclique_[bestvar]) << " " << getUB(solver.isclique_[bestvar]) << ", value " << val[bestvar] << ", feasibility " << feas[bestvar] << std::endl;
            std::cout << "revorder: branch callback: branch on the variable isclique_" << bestvar << std::endl;
#endif
}
//        else {
//            Vertex* bestu = NULL;
//            Vertex* bestv = NULL;
//            for (Vertex* u : solver.inst_->vertices_) {
//                std::map<Vertex*, IloNum> isbeforeval;
//                float sumrefs = 0;
//                for (Vertex* v : u->neighbors_) {
//                    isbeforeval[v] = getValue(solver.isbefore_[v][u]);
//                    sumrefs += isbeforeval[v];
//                }
//                if (sumrefs == solver.inst_->U()) {
//                    for (Vertex* v : u->neighbors_) {
//                        if (getFeasibility(solver.isbefore_[v][u]) == Infeasible) {
//                            if (std::min(isbeforeval[v], 1 - isbeforeval[v]) > maxval) {
//                                bestu = u;
//                                bestv = v;
//                                maxval = std::min(isbeforeval[v], 1 - isbeforeval[v]);
//                            }
//                        }
//                    }
//                }
//            }
//            if (bestu) {
//                // HeuristicNode* data = new HeuristicNode();
//                // data->fixcliqueindex(bestvar);
//                std::cout << "revorder: branch on isbefore_" << bestv->id_ << "_" << bestu->id_ << std::endl;
//                makeBranch(solver.isbefore_[bestv][bestu], getValue(solver.isbefore_[bestv][bestu]), IloCplex::BranchUp, getObjValue());
//                makeBranch(solver.isbefore_[bestv][bestu], getValue(solver.isbefore_[bestv][bestu]), IloCplex::BranchDown, getObjValue());
//            }
//        }
} catch (...) {

val.end();
feas.end();
throw;
}
val.end();
feas.end();
}


//
// constructors and destructor

Clique::Clique(std::vector<Vertex *> vert) : nbvertices_(vert.size()), vertices_(vert) {
}

Clique::~Clique() {
}

DiscretizationSolver::DiscretizationSolver(Instance* inst, Algorithm algo) : algo_(algo), inst_(inst) {
  this->bestnbfullref_ = -1;
  this->bestcliqueid_ = -1;
  for (Vertex* v : this->inst_->vertices_) {

    this->bestrank_[v] = -1;
    this->bestnbrefs_[v] = -1;
  }
  this->objvalue_ = inst->nbvertices_;
  this->cliquecutsmaxsize_ = inst->nbvertices_;
}
//

DiscretizationSolver::DiscretizationSolver(Instance* inst, Algorithm algo, Problem pb, Model mod, std::string optionfile, float timelimit) :
    problem_(pb), algo_(algo), modeltype_(mod), timelimit_(timelimit), inst_(inst) {
  this->bestnbfullref_ = -1;
  this->bestcliqueid_ = -1;
  for (Vertex* v : this->inst_->vertices_) {

    this->bestrank_[v] = -1;
    this->bestnbrefs_[v] = -1;
  }
  this->objvalue_ = inst->nbvertices_;
  if (!optionfile.empty())
  {
    std::ifstream infile(optionfile.c_str());
    if (!infile.is_open()) {
      throwError("revorder: error: the option file could not be opened");
    }
    // read the file line by line
    // one line contains the data relative to one option
    std::string line;
    while (std::getline(infile, line))
    {
      std::istringstream iss(line);
      char buf[100];
      std::string optionname;
      int optionvalue;

      //read the line and detect format errors
      if (!(iss >> buf)) {
        std::cerr << "revorder: error: there is a mistake with the  "
                     "format of the option file" <<std::endl;
        throw;
      } // error
      optionname = std::string(buf);
      //
      if (optionname.find("cliquecutsmaxsize") !=std::string::npos) {
        if (!(iss >> optionvalue)) {
          std::cerr << "revorder: error: there is a mistake with "
                       "the format of the option file" <<std::endl;
          throw;
        } // error
        this->cliquecutsmaxsize_ = std::min(optionvalue, inst->nbvertices_);
        continue;
      } else if (optionname.find("initialcyclesize") !=
          std::string::npos) {
        if (!(iss >> optionvalue)) {
          std::cerr << "revorder: error: there is a mistake with "
                       "the format of the option file" << std::endl;
          throw;
        } // error
        this->initialcyclesize_ =
            std::min(optionvalue, inst->nbvertices_);
        continue;
      } else if (optionname.find("decomposecomponents") !=
          std::string::npos) {
        if (!(iss >> optionvalue)) {
          std::cerr << "revorder: error: there is a mistake with "
                       "the format of the option file" <<std::endl;
          throw;
        } // error
        this->decomposecomponents_ = (bool) optionvalue;
        continue;
      }
    }
  } else {
    this->decomposecomponents_ = false;
    this->cliquecutsmaxsize_ = inst->nbvertices_;
    this->initialcyclesize_ = 3;
  }
}
//

DiscretizationSolver::~DiscretizationSolver() {

  for (Clique* clique : this->cliques_) delete clique;
  cliques_.clear();
}
//
// abstract solve method

int DiscretizationSolver::solve() {
  return 0;
}
//
// get the objective value of a given vertex order

int DiscretizationSolver::getobjvalue(std::map<Vertex*, int> rank, std::map<Vertex*, int> nbrefs) {

  int nbfullref = 0;
  for (Vertex* u : this->inst_->vertices_) {

    if (nbrefs[u] >= this->inst_->U()) nbfullref += 1;
  }
  return this->inst_->nbvertices_ - this->inst_->L() - nbfullref;

}
//
// get number of fully referenced vertices in an order

int DiscretizationSolver::getnbfullref(std::map<Vertex*, int> rank, std::map<Vertex*, int> nbrefs) {

  int nbfullref = 0;
  for (Vertex* u : this->inst_->vertices_) {

    if (nbrefs[u] >= this->inst_->U()) nbfullref += 1;
  }
  return nbfullref;

}
//
// return true if objval1 is better than objval2: this method is necessary, because the problem can be a minimization or a maximization

bool DiscretizationSolver::isabetterobjvalue(int objval1, int objval2) {

  return objval1 > objval2;
}
//
// compute the best greedy solution among those starting from the potential initial cliques and return true if the instance is feasible
// if allcliques is set to true, check every clique, otherwise, stop with the first feasible revorder

bool DiscretizationSolver::greedysolve(bool allcliques) {

  int U = this->inst_->U();
  int L = this->inst_->L();
  int nbvertices = this->inst_->nbvertices_;
  this->bestnbfullref_ = -1;
  this->bestcliqueid_ = -1;
  for (Vertex* u : this->inst_->vertices_) {
    this->bestisfullref_[u] = false;
    this->bestrank_[u] = -1;
    this->bestnbrefs_[u] = 0;
  }
  this->bestfullref_.clear();
  this->objvalue_ = this->inst_->nbvertices_;
  std::vector<Clique*> feasiblecliques;

  //
  // for each clique apply the constructive greedy that iteratively adds the
  // vertex with largest number of references
  std::map<Vertex*, int> rank;
  std::map<Vertex*, int> nbrefs;
  this->greedynbfullrefperclique_.clear();
  int feasiblecliqueid = -1;
  for (int initialcliqueindex = 0; initialcliqueindex < this->cliques_.size(); initialcliqueindex++) {
    //
    // data initialization
    Clique* initialclique = this->cliques_[initialcliqueindex];
    int currentrank = 1;
    std::vector<Vertex*> verticesleft = this->inst_->vertices_;
    int nbfullref = 0;
    for (Vertex * v : this->inst_->vertices_) {
      nbrefs[v] = 0;
      rank[v] = -1;
    }
    //
    // treat the initial clique first
    //
    // set the ranks of the vertices in the clique and set them as references of their neighbors
    for (Vertex* u : initialclique->vertices_) {
      rank[u] = currentrank++;
      for (Vertex* v : u->neighbors_) {
        if (rank[v] == -1) nbrefs[v]++;
      }
    }
    //
    // record the partially and fully referenced vertices
    std::vector<Vertex*> partialref;
    std::vector<Vertex*> fullref;
    for (Vertex* v : this->inst_->vertices_) {
      if (rank[v] == -1) {
        if (nbrefs[v] >= U) {
          fullref.push_back(v);
        } else if (nbrefs[v] >= L) {
          partialref.push_back(v);
        }
      }
    }
    //
    // include the vertex with largest number of references in the order until every vertex is treated or infeasibility of the initial clique is proved
    while (currentrank <= nbvertices) {
      if (!fullref.empty()) {
        Vertex* u = fullref.back();
        if (initialclique->initialpartialrefs_.empty()) {
          initialclique->initialfullrefs_.push_back(u);
        }
        rank[u] = currentrank++;
        nbfullref++;
        fullref.pop_back();
        for (Vertex* v : u->neighbors_) {
          if (rank[v] >= 1) continue;
          nbrefs[v]++;
          if (nbrefs[v] == U) {
            fullref.push_back(v);
          } else if (nbrefs[v] == L) {
            partialref.push_back(v);
          }
        }
      } else if (!partialref.empty()) {
        while (rank[partialref.back()] >= 1) {
          partialref.pop_back();
          if (partialref.empty()) {
            break;
          }
        }
        if (partialref.empty()) break;

        if (initialclique->initialpartialrefs_.empty()) {
          for (Vertex* v: partialref) {
            initialclique->initialpartialrefs_.push_back(v);
          }
        }
        Vertex * u = partialref.back();
        rank[u] = currentrank++;
        partialref.pop_back();
        for (Vertex* v : u->neighbors_) {
          if (rank[v] >= 0) continue;
          nbrefs[v]++;
          if (nbrefs[v] == U) {
            fullref.push_back(v);
          } else if (nbrefs[v] == L) {
            partialref.push_back(v);
          }
        }
      } else break;
    }

    if (currentrank - 1 == nbvertices) {
      feasiblecliqueid++;
      feasiblecliques.push_back(initialclique);
      this->greedynbfullrefperclique_.push_back(nbfullref);
      initialclique->greedyobjvalue_ = this->getobjvalue(rank, nbrefs);
      //
#ifdef VERBOSE
      std::cout << "revorder: greedy: clique computed solution with value " << initialclique->greedyobjvalue_ << std::endl;
#endif
      //
      // record the solution if a new incumbent is found
      if ((currentrank - 1 == nbvertices) && (initialclique->greedyobjvalue_ < this->objvalue_)) {
        //
#ifdef VERBOSE
        std::cout << "revorder: greedy: one new clique with cost " << initialclique->greedyobjvalue_ << '\n';
#endif
        //
        this->objvalue_ = initialclique->greedyobjvalue_;
        this->bestcliqueid_ = feasiblecliqueid;
        this->bestfullref_.clear();
        this->bestnbfullref_ = 0;
        for (Vertex* v : this->inst_->vertices_) {
          this->bestrank_[v] = rank[v];
          this->bestnbrefs_[v] = nbrefs[v];
          if (this->bestnbrefs_[v] >= U) {
            this->bestisfullref_[v] = true;
            this->bestfullref_.push_back(v);
            this->bestnbfullref_++;
          } else {
            this->bestisfullref_[v] = false;
          }
        }
        //
        // if we do not need to check every clique, stop with the first feasible revorder
        if (!allcliques) break;
      }
    } else {
#ifdef VERBOSE
      std::cout << "revorder: greedy: this clique is not a feasible start" << std::endl;
#endif
    }
  }
  //
  // build the list of cliques that will be used in the optimization
  this->cliques_.clear();
  for (Clique* clique : feasiblecliques) {
    this->cliques_.push_back(clique);
  }
  //
  // update the list of cliques ids each vertex is a member of
  for (Vertex* v : this->inst_->vertices_) {
    this->cliqueslist_[v].clear();
    this->cliqueidslist_[v].clear();
  }
  int id = 0;
  for (Clique* clique : this->cliques_) {
    for (Vertex* v : clique->vertices_) {
      this->cliqueslist_[v].push_back(clique);
      this->cliqueidslist_[v].push_back(id);
    }
    id++;
  }
  std::cout << "revorder: greedy: number of cliques after greedy solution = " << this->cliques_.size() << std::endl;
  //
  // Stop here and return false if the problem is infeasible
  if (this->bestnbfullref_ == -1) {
    std::cout << "revorder: greedy: no solution was found with the greedy, the problem is infeasible" << std::endl;
    return false;
  }
  //
  // Record the solution found
  std::cout << "revorder: greedy: number of part. ref. vertices in the best solution = " << this->objvalue_ << std::endl;
  this->greedynbfullref_ = this->bestnbfullref_;

  return true;
}
//
//
// run the search for a discretization order with the input options

bool DiscretizationSolver::discretizationorder(float timelimit) {

  bool feas_status = true;

  if (this->algo_ == GREEDY) {
    //
    // enumerate the potential initial cliques
    this->inst_->computeneighbors();
    feas_status = this->enumeratecliques(this->inst_->L() + 1);
    if (feas_status == false) {
      this->isfeasible_ = false;
      return false;
    }
    this->eliminateredundantcliques();
    //
    // run the greedy search
    feas_status = greedysolve();
    this->isfeasible_ = feas_status;
    return feas_status;
  }


  switch (this->algo_) {
    case IP:
      feas_status = minpartialip(timelimit);
      break;
    case CP:
      feas_status = cpvertex(timelimit);
      break;

    default:
      std::cout << "revorder: error with the problem and/or algorithm fed to the solver" << std::endl;
      throw;
  }


  return feas_status;
}
//
// greedy algorithm to get a reference primal bound
// try all the possible initial cliques as starting points of the greedy

bool DiscretizationSolver::greedymipstart(IloCplex & cplex) {
  //
  // set the best greedy solution as primal solution for the mip solution
  IloNumArray mipstartvalues(cplex.getEnv());
  IloNumVarArray mipstartvariables(cplex.getEnv());
  if (this->modeltype_ == VERTEXRANK) {
    for (Vertex* v : this->inst_->vertices_) {
      for (int k = 0; k < this->inst_->nbvertices_; k++) {
        mipstartvariables.add(this->hasrank_[v][k]);
        if (this->bestrank_[v] == k + 1) mipstartvalues.add(1);
        else mipstartvalues.add(0);
      }
    }
  } else {
    if (this->modeltype_ == WITNESS) {
      std::map<Vertex*, int> nbwitness;
      for (Vertex* u : this->inst_->vertices_) nbwitness[u] = 0;

      for (Vertex* u : this->inst_->vertices_) {
        if (this->bestrank_[u] <= this->inst_->L() + 1) {
          for (Vertex* v : u->neighbors_) {
            mipstartvariables.add(this->isbefore_[u][v]);
            mipstartvalues.add(1);
            nbwitness[v] += 1;
          }
        }
      }
      for (Vertex* u : this->inst_->vertices_) {
        if (this->bestrank_[u] >= this->inst_->L() + 2) {
          for (Vertex* v : u->neighbors_) {
            mipstartvariables.add(this->isbefore_[u][v]);
            if ((this->bestrank_[u] < this->bestrank_[v])) {
              if ((this->bestisfullref_[v]) && (nbwitness[v] < this->inst_->U())) {
                mipstartvalues.add(1);
                nbwitness[v] += 1;
              } else if ((!this->bestisfullref_[v]) && (nbwitness[v] < this->inst_->L())) {
                mipstartvalues.add(1);
                nbwitness[v] += 1;
              } else mipstartvalues.add(0);
            } else mipstartvalues.add(0);
          }
        }
      }
    } else {
      for (Vertex* u : this->inst_->vertices_) {
        for (Vertex* v : u->neighbors_) {
          if (u->id_ >= v->id_) continue;
          mipstartvariables.add(this->isbefore_[u][v]);
          mipstartvariables.add(this->isbefore_[v][u]);
          if (this->bestrank_[u] < this->bestrank_[v]) {
            mipstartvalues.add(1);
            mipstartvalues.add(0);
          } else {
            mipstartvalues.add(0);
            mipstartvalues.add(1);
          }
        }
      }
    }
    //
    // set the clique variables and corresponding clique variables
    if (this->modeltype_ == WITNESS) {
      for (Vertex* v : this->inst_->vertices_) {
        mipstartvariables.add(this->isvertexinclique_[v]);
        if (this->bestrank_[v] <= this->inst_->L() + 1) mipstartvalues.add(1);
        else mipstartvalues.add(0);
      }
    } else {
      for (int c = 0; c < this->cliques_.size(); c++) {
        mipstartvariables.add(this->isclique_[c]);
        if (c == this->bestcliqueid_) mipstartvalues.add(1);
        else mipstartvalues.add(0);
      }
    }
    //
    // full-ref variables
    for (Vertex* u : this->inst_->vertices_) {
      mipstartvariables.add(this->isfullref_[u]);
      if ((this->bestisfullref_[u]) || (this->bestrank_[u] <= this->inst_->L() + 1)) mipstartvalues.add(1);
      else mipstartvalues.add(0);
    }
    //
    // rank variables for RANKS model
    for (Vertex* v : this->inst_->vertices_) {
      if (this->modeltype_ == RANKS) {
        for (int k = 0; k < this->inst_->nbvertices_; k++) {
          mipstartvariables.add(this->hasrank_[v][k]);
          if (this->bestrank_[v] == k + 1) mipstartvalues.add(1);
          else mipstartvalues.add(0);
        }
        mipstartvariables.add(this->rank_[v]);
        mipstartvalues.add(this->bestrank_[v] - 1);
      }
    }
  }

  cplex.addMIPStart(mipstartvariables, mipstartvalues);
  mipstartvariables.end();
  mipstartvalues.end();

  return true;
}


//
// contraint programming CP^VERTEX from Bodur and MacNeil, 2019

bool DiscretizationSolver::cpvertex(float timelimit) {

  char name[50];
  int nbvertices = this->inst_->nbvertices_;
  int L = this->inst_->L();
  int U = this->inst_->U();
  IloEnv env;
  IloTimer cpuClockTotal(env);
  IloTimer cpuClockInit(env);
  cpuClockTotal.start();
  cpuClockInit.start();


  std::cout << std::endl;
  std::cout << "revorder: ----------------------------------------------------------" << std::endl;
  std::cout << "\nrevorder: cp: Run preprocessing procedures " << std::endl << std::endl;

  //
  // compute adjacency lists
  this->inst_->computeneighbors();
  //
  // identify the set of feasible initial cliques
  if (!this->enumeratecliques(this->inst_->L() + 1)) {
    this->isfeasible_ = false;
    return false;
  }
  this->eliminateredundantcliques();
  //
  // run the greedy algorithm for initial solution and reduction of initial cliques
  if (!this->greedysolve()) {
    std::cout << "Greedy solve found the problem to be infeasible at the time of mip warm start." << std::endl;
    return false;
  }

  //
  // initialize the model
  IloModel model(env);
  //
  // assign one unique index in {0,...,n-1} to each vertex
  std::map<Vertex*, int> indv;
  int ind = 0;
  for (Vertex* u : this->inst_->vertices_) {
    indv[u] = ind++;
  }

  //
  // initialize adjacency matrix
  IloIntArray adjmatrix(env, nbvertices * nbvertices);
  for (Vertex* u : this->inst_->vertices_) {
    for (Vertex* v : this->inst_->vertices_) {
      if (indv[v] < indv[u]) continue;
      if (u->isneighbor_[v]) {
        adjmatrix[nbvertices * indv[u] + indv[v]] = 1;
        adjmatrix[nbvertices * indv[v] + indv[u]] = 1;
      } else {
        adjmatrix[nbvertices * indv[u] + indv[v]] = 0;
        adjmatrix[nbvertices * indv[v] + indv[u]] = 0;
      }
    }
  }
  //
  // variables that set the rank of each vertex in the order
  //
  // isrankfullref[k] =1 if vertex v is fully-referenced and 0 otherwise
  IloBoolVarArray isrankfullref(env, nbvertices);
  for (int k = 0; k < nbvertices; k++) {
    isrankfullref[k] = IloBoolVar(env);
    sprintf(name, "isrankfullref_%i", k);
    isrankfullref[k].setName(name);
  }
  //
  // indatrank[r] = index of the vertex at rank r
  IloIntVarArray indatrank(env, nbvertices);
  for (int k = 0; k < nbvertices; k++) {
    indatrank[k] = IloIntVar(env, 0, nbvertices - 1);
    sprintf(name, "indatrank_%i", k);
    indatrank[k].setName(name);
  }
  //
  // objective function
  IloExpr obj(env);
  obj += 1;
  for (int k = L + 1; k < nbvertices; k++) {
    obj += (1 - isrankfullref[k]);
  }
  model.add(IloMinimize(env, obj));

  for (int i = 0; i <= L - 1; i++) {
    for (int j = i + 1; j <= L; j++) {
      model.add(IloElement(adjmatrix, nbvertices * indatrank[i] + indatrank[j]) == 1);
    }
  }

  for (int r = L + 1; r < nbvertices; r++) {
    IloExpr sumrefs(env);
    for (int i = 0; i <= r - 1; i++) {
      sumrefs += IloElement(adjmatrix, nbvertices * indatrank[i] + indatrank[r]);
    }
    model.add(sumrefs - (U - L) * isrankfullref[r] >= L);
    sumrefs.end();
  }

  model.add(IloAllDiff(env, indatrank));
  //
  // 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 r = 0; r <= L; r++) {
        model.add(indatrank[r] != indv[v]);
      }
    }
  }
  //
  // load model in the CP solver
  IloCP cp(model);
  //
  // provide starting point with greedy solution
  //
  IloSolution sol(env);
  for (Vertex* v : this->inst_->vertices_) {
    if (this->bestrank_[v] - 1 >= L + 1) {
      sol.setValue(isrankfullref[this->bestrank_[v] - 1], bestisfullref_[v]);
    }
    sol.setValue(indatrank[this->bestrank_[v] - 1], indv[v]);
  }
  cp.setStartingPoint(sol);

  cpuClockInit.stop();
  //
  // set CP parameter
  cp.setParameter(IloCP::LogVerbosity, IloCP::Normal);
  cp.setParameter(IloCP::Workers, 1);
  cp.setParameter(IloCP::TimeLimit, timelimit - cpuClockInit.getTime());
#ifdef VERBOSE
  cp.setParameter(IloCP::LogVerbosity, IloCP::Normal);
#endif
  IloTimer cpuClock(env);
  cpuClock.start();
  try {
    cp.solve();
  } catch (IloException& e) {
    std::cerr << e << std::endl;
  }
  cpuClockTotal.stop();

  IloAlgorithm::Status status = cp.getStatus();
  std::cout << "revorder: ----------------------------------------------------------\n" << std::endl;
  std::cout << "revorder: CP optimizer 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
    this->objvalue_ = cp.getObjValue();
    totaltime_ = cpuClock.getTime();
    treatednodes_ = cp.getInfo(IloCP::NumberOfBranches);

    //
    // get the optimal order
    for (int r = 0; r < this->inst_->nbvertices_; r++) {
      this->bestrank_[this->inst_->vertices_[cp.getValue(indatrank[r])]] = r;
    }
    //
    // get the number of references of each vertex
    for (Vertex* u : this->inst_->vertices_) {
      this->bestnbrefs_[u] = 0;
      for (Vertex* v : u->neighbors_) {
        if (this->bestrank_[v] < this->bestrank_[u]) this->bestnbrefs_[u]++;
      }
    }
    //
    // get fully referenced vertices
    this->bestfullref_.clear();
    this->bestnbfullref_ = 0;
    this->bestisfullref_.clear();
    for (Vertex* u : this->inst_->vertices_) {
      if (this->bestnbrefs_[u] >= this->inst_->U()) {
        this->bestisfullref_[u] = true;
        this->bestfullref_.push_back(u);
        this->bestnbfullref_++;
      } else this->bestisfullref_[u] = false;
    }
    //
    // summary of cp execution
    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 branches explored       = " << treatednodes_ << std::endl;
    std::cout << "revorder: value of the objective function   = " << objvalue_ << std::endl;
    std::cout << "revorder: verification: number of part. ref. vertices = " << this->inst_->nbvertices_ - this->inst_->L() - this->bestnbfullref_ << std::endl;
  } else {
    std::string fileInfeasible("InfeasibleModel.lp");
    cp.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;
  }
  //
  // delete CP objects
  cp.end();
  model.end();
  env.end();

  return isfeasible_;
}
//
// create the cplex model for the IP formulation described in Bodur an MacNeil, 2019

void DiscretizationSolver::defineminpartial_vertexrank(IloModel & model) {
  IloEnv env = model.getEnv();
  char name[50];
  int nbvertices = this->inst_->nbvertices_;
  int L = this->inst_->L();
  int U = this->inst_->U();

  //////////////////////////////////////////////////////////////////////////////
  // Declare the variables
  //////////////////////////////////////////////////////////////////////////////
  //
  // 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)
  for (Vertex* v : this->inst_->vertices_) {
    hasrank_[v] = IloBoolVarArray(env, nbvertices);
    for (int k = 0; k < nbvertices; k++) {
      sprintf(name, "hasrank_%i_%i", v->id_, k);
      hasrank_[v][k].setName(name);
    }
  }
  //
  // isrankfullref[k] =1 if vertex v is fully-referenced and 0 otherwise
  IloBoolVarArray isrankfullref(env, nbvertices);
  for (int k = 0; k < nbvertices; k++) {
    isrankfullref[k] = IloBoolVar(env);
    sprintf(name, "isrankfullref_%i", k);
    isrankfullref[k].setName(name);
  }
  //
  // isfullrefatrank[u,k] = 1 if vertex u is fully-referenced and is at rank k
  BoolVarMapArray isfullrefatrank;
  for (Vertex* v : this->inst_->vertices_) {
    isfullrefatrank[v] = IloBoolVarArray(env, nbvertices);
    for (int k = 0; k < this->inst_->nbvertices_; k++) {
      sprintf(name, "isfullrefatrank%i_%i", v->id_, k);
      isfullrefatrank[v][k].setName(name);
    }
  }

  //////////////////////////////////////////////////////////////////////////////
  // Set the cplex model that will be solved
  //////////////////////////////////////////////////////////////////////////////
  //
  // objective function
  IloExpr obj(env);
  obj += 1;
  for (int k = 0; k < nbvertices; k++) {
    obj += (1 - isrankfullref[k]);
  }
  model.add(IloMinimize(env, obj));

  //
  // 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);
  int nbcons = 0;
  for (Vertex* u : this->inst_->vertices_) {
    IloExpr sumhasrank(env);
    for (int k = 0; k < nbvertices; k++) sumhasrank += hasrank_[u][k];
    ctaryOneRankPerVertex.add(sumhasrank == 1);
    sprintf(name, "ctaryOneRankPerVertex_%i", u->id_);
    ctaryOneRankPerVertex[nbcons].setName(name);
  }
  model.add(ctaryOneRankPerVertex);
  //
  // Guarantee that isfullrefatrank[v][k] =1 only if v has rank k and has at least U references
  IloRangeArray ctaryIsFullRefAtLevel(env);
  for (int k = 0; k <= L; k++) {
    ctaryIsFullRefAtLevel.add(isrankfullref[k] == 1);
  }
  for (int k = 0; k < this->inst_->nbvertices_; k++) {
    for (Vertex* v : this->inst_->vertices_) {
      ctaryIsFullRefAtLevel.add(hasrank_[v][k] - (1 - isrankfullref[k]) - isfullrefatrank[v][k] <= 0);
    }
  }
  model.add(ctaryIsFullRefAtLevel);
  //
  // Discretization constraints: every vertex must be partially-referenced and those with more than U references are fully-referenced
  IloRangeArray ctaryLrefs(env);
  IloRangeArray ctaryClique(env);
  IloRangeArray ctaryFullRef(env);
  nbcons = 0;
  for (Vertex* u : this->inst_->vertices_) {
    //
    // first deal with the vertices of the initial clique
    for (int k = 0; k <= L; k++) {
      IloExpr sumrefs(env);
      for (Vertex* v : u->neighbors_) {
        for (int l = 0; l <= k - 1; l++) sumrefs += hasrank_[v][l];
      }
      ctaryClique.add(sumrefs - k * hasrank_[u][k] >= 0);
    }
    //
    // then deal with the rest of the order
    for (int k = L + 1; k < nbvertices; k++) {
      IloExpr sumrefs(env);
      for (Vertex* v : u->neighbors_) {
        for (int l = 0; l <= k - 1; l++) sumrefs += hasrank_[v][l];
      }
      ctaryLrefs.add(sumrefs - L * hasrank_[u][k] >= 0);
      ctaryFullRef.add(sumrefs - U * isfullrefatrank[u][k] >= 0);
    }
  }
  model.add(ctaryClique);