/********************************************************************************************
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);
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);
}