/********************************************************************************************
Name: DiscretizationSolver
Discrete optimization for graph discretization
Author: J.Omer
Sources: C++
License: GNU General Public License v.2
History:
*********************************************************************************************/
#include "bbsolver.hpp"
#include <algorithm>
ILOSTLBEGIN
//
ILOLAZYCONSTRAINTCALLBACK2(CyclesLazyConstraintsBB,
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.U();
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;
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;
}
/*************************************************************************************************************************************************************/
//
// sort the nodes by ascending number of ordered vertices
bool bbnodesnborderedcompare(BBNode* n1, BBNode* n2) {
return (n1->nbordered() < n2->nbordered());
}
bool bbnodesbestcompare(BBNode*n1, BBNode* n2) {
return (static_cast<float>(n1->cost_) / n1->nbordered()
> static_cast<float>(n2->cost_) / n2->nbordered());
}
// Public methods
//
// constructor and destructor
//
// constructor used only for the root
BBNode::BBNode(Instance& inst):
depth_(0),
inst_(&inst),
nbordered_(0),
nbfullref_(0),
istreated_(true) {}
//
// constructor used only for the initial cliques
BBNode::BBNode(Instance& inst,
BBNode* root,
Clique& c,
bool branchonsmallerindex):
depth_(1),
father_(root),
inst_(&inst),
nbordered_(0),
nbfullref_(0),
istreated_(false) {
this->orderedvertices_.clear();
for (Vertex* v : inst.vertices_) {
this->rank_[v] = -1;
this->nbrefs_[v] = 0;
this->isfullref_[v] = false;
this->isordered_[v] = false;
this->ispotentialchild_[v] = false;
this->mustbefullref_[v] = false;
this->maxrefs_[v] = v->degree_;
}
int currentrank = 1;
for (Vertex* u : c.vertices_) {
this->rank_[u] = currentrank++;
this->orderedvertices_.push_back(u);
this->switchtoordered(u);
this->nbordered_++;
//
// update the number of references of the neighbors of the vertex
for (Vertex* v : u->neighbors_) {
if (!this->isordered_[v]) {
this->nbrefs_[v]++;
}
}
}
// record the list of partially and fully referenced (unordered) vertices
for (Vertex* v : inst.vertices_) {
if (!this->isordered_[v]) {
// record the list of unordered fully-referenced vertices
if (this->nbrefs_[v] >= inst.U()) {
this->isfullref_[v] = true;
this->fullrefs_.push_back(v);
} else if (this->nbrefs_[v] >= inst.L()) {
// add partially referenced vertices to the list of potential choices
// for branching
if (branchonsmallerindex) {
bool largerindex = true;
Vertex* u = c.vertices_.back();
if (!u->isneighbor(v)) {
if (u->id_ >= v->id_) largerindex = false;
}
//
// to break symmetries, make sure that the chosen potential child
// is always that with biggest index
if (largerindex) {
this->ispotentialchild_[v] = true;
this->potentialchildren_.push_back(v);
//
// partially referenced vertices with less than U neighbors can
// be ordered right away since they will never be fully referenced
if (this->maxrefs_[v] <= inst.U() - 1) {
this->readytoorder_.push_back(v);
#ifdef VERBOSE
std::cout << "revorder: bb: order vertex " << v->id_ << " without branching, maximum nb of references = " << this->maxrefs_[v] << std::endl;
#endif
}
} else {
this->mustbefullref_[v] = true;
// if a vertex must be fully-referenced, but it has exactly U
// neighbors, it will for sure come after all its neighbors in
// the order
if (this->maxrefs_[v] <= inst.U()) {
for (Vertex* u : v->neighbors_) {
if (this->mustbefullref_[u]) {
if (this->maxrefs_[u] >= inst.U()+1) {
this->maxrefs_[u]--;
}
} else {
this->maxrefs_[u]--;
}
}
}
}
} else {
this->ispotentialchild_[v] = true;
this->potentialchildren_.push_back(v);
// partially referenced vertices with less than U neighbors can be
// ordered right away since they will never be fully referenced
if (this->maxrefs_[v] <= inst.U() - 1) {
this->readytoorder_.push_back(v);
#ifdef VERBOSE
std::cout << "revorder: bb: order vertex " << v->id_ << " without branching, maximum nb of references = " << this->maxrefs_[v] << std::endl;
#endif
}
}
}
}
}
this->dualbound_ = -1;
root->addtochildren(this);
}
// Make a copy of the input BB node: used for every other node
BBNode::BBNode(BBNode& node): id_(node.id_),
depth_(node.depth_+1),
father_(&node),
inst_(node.inst_),
nbordered_(node.nbordered_),
nbfullref_(node.nbfullref_),
istreated_(false) {
cost_ = node.cost_;
dualbound_ = node.dualbound_;
primalbound_ = node.primalbound_;
for (Vertex* v : inst_->vertices_) {
rank_[v] = node.rank_.at(v);
isordered_[v] = node.isordered_.at(v);
nbrefs_[v] = node.nbrefs_.at(v);
isfullref_[v] = node.isfullref_.at(v);
ispotentialchild_[v] = node.ispotentialchild_.at(v);
mustbefullref_[v] = node.mustbefullref_.at(v);
maxrefs_[v] = node.maxrefs_.at(v);
}
fullrefs_.clear();
orderedvertices_.clear();
potentialchildren_.clear();
children_.clear();
for (Vertex* v : node.readytoorder_) readytoorder_.push_back(v);
for (Vertex* v : node.fullrefs_) fullrefs_.push_back(v);
for (Vertex* v : node.orderedvertices_) orderedvertices_.push_back(v);
for (Vertex* v : node.potentialchildren_) potentialchildren_.push_back(v);
}
//
// Destroy the node
BBNode::~BBNode() {
}
//
// remove a vertex from the list of potential choice for branching
void BBNode::removefrompotentialchildren(Vertex* v) {
//
// no need to proceed if the vertex is not marked as potential child
if (!this->ispotentialchild_[v]) return;
this->ispotentialchild_[v] = false;
std::vector<Vertex*>::iterator itv = this->potentialchildren_.begin();
while (itv < this->potentialchildren_.end()) {
if (*itv == v) {
this->potentialchildren_.erase(itv);
break;
} else {
itv++;
// detect error if the vertex was not found in the list of potential
// children
if (itv == this->potentialchildren_.end()) {
std::cout << "revorder: bb: did not find a potential child in the list"
<< '\n';
throw;
}
}
}
itv = this->readytoorder_.begin();
while (itv < this->readytoorder_.end()) {
if (*itv == v) {
this->readytoorder_.erase(itv);
break;
}
itv++;
}
}
//
// assign the next available rank to the input vertex
void BBNode::assignnextrank(Vertex* u) {
//
// assign the next available rank
for (Vertex* v : this->orderedvertices_) {
if (u == v) {
std::cout << "Les ennuis commencent ici." << std::endl;
}
}
this->rank_[u] = this->nbordered_ + 1;
this->orderedvertices_.push_back(u);
this->switchtoordered(u);
this->nbordered_++;
this->removefrompotentialchildren(u);
if (this->nbrefs_[u] <= this->inst_->L() - 1) {
std::cout << "revorder: bb: trying to assign a rank to a vertex with less"
" than L references" << '\n';
throw;
} else if (this->nbrefs_[u] >= this->inst_->U()) {
this->nbfullref_++;
}
// update the number of references of the neighbors of v
for (Vertex* v : u->neighbors_) {
if (this->rank_[v] < 0) {
this->nbrefs_[v]++;
// add partially-referenced vertices to the list of potential choices
// for branching
if (this->nbrefs_[v] == this->inst_->L()) {
this->potentialchildren_.push_back(v);
this->ispotentialchild_[v] = true;
// partially referenced vertices with less than U neighbors can be
// ordered right away since they will never be fully referenced
if (this->maxrefs_[v] <= this->inst_->U() - 1) {
this->readytoorder_.push_back(v);
#ifdef VERBOSE
std::cout << "revorder: bb: order vertex " << v->id_ << " without branching, maximum nb of references = " << this->maxrefs_[v] << std::endl;
#endif
}
}
// a vertex is fully-referenced if it has more than U references, and it
// is not a potential choice for branching anymore
if (this->nbrefs_[v] == this->inst_->U()) {
this->isfullref_[v] = true;
this->fullrefs_.push_back(v);
if (this->ispotentialchild_[v]) this->removefrompotentialchildren(v);
}
}
}
}
// get the current objective value at a node, depending on the ordered vertices
void BBNode::setnodecost(Problem pb) {
int nbpartials = 0;
for (Vertex* u : this->inst_->vertices_) {
if (this->isordered_[u]) {
if (this->nbrefs_[u] <= this->inst_->U() - 1) nbpartials += 1;
}
}
this->cost_ = nbpartials;
}
// a bb node n1 dominates another one, n2, if the corresponding partial
// solution has a smaller or equal cost but its list of ordered vertices
// contains that of n2
bool BBSolver::dominates(BBNode& n1, BBNode& n2) {
if ( n2.cost_ < n1.cost_ ) {
return false;
}
for (Vertex* v : n2.orderedvertices_) {
if (!(n1.isordered(v))) {
return false;
}
}
return true;
}
//
// propagate the set of ordered vertices by iteratively ordering
// fully-referenced vertices
void BBNode::prepropagation() {
#ifdef VERBOSE
cout << "prepropagation" << endl;
#endif
// iteratively propagate while there are fully-referenced vertices
while (!this->fullrefs_.empty()) {
Vertex* u = this->fullrefs_.back();
this->fullrefs_.pop_back();
this->isfullref_[u] = false;
this->assignnextrank(u);
}
this->setnodecost(MINPARTIAL);
}
//
// propagate the set of ordered vertices by iteratively ordering vertices
// as long as there is either at least one unordered vertex with more than U
// references or exactly one unordered vertex with exactly U references
void BBNode::postpropagation() {
#ifdef VERBOSE
cout << "in postpropagation" << endl;
#endif
//
// no need of propagation if all the vertices are ordered
if (this->nbordered_ == this->inst_->nbvertices_) return;
//
// iteratively propagate while there are fully-referenced vertices or
// exactly one partially-referenced vertex
while (true) {
bool ispropagate = false;
Vertex* u;
if (!this->fullrefs_.empty()) {
// add a vertex with more than U references
u = this->fullrefs_.back();
this->fullrefs_.pop_back();
this->isfullref_[u] = false;
ispropagate = true;
} else if (!this->readytoorder_.empty()) {
// if there is a partially referenced vertex with less than U
// neighbors, it can be ordered right away: it will never be fully
// referenced
u = this->readytoorder_.back();
this->readytoorder_.pop_back();
ispropagate = true;
} else if (this->potentialchildren_.size() == 1) {
// if there is exactly one partially referenced vertex: add it to the
// order
u = this->potentialchildren_.back();
this->potentialchildren_.pop_back();
this->ispotentialchild_[u] = false;
ispropagate = true;
}
if (!ispropagate) break;
//
// add the identified vertex at the end of the order
this->assignnextrank(u);
}
//
// update the cost of the node
this->setnodecost(MINPARTIAL);
}
/*******************************************************************************
Class BBSolver
- an instance of BBSolver is a solver of the discretization problem using
branch-and-bound
*******************************************************************************/
// Public methods
//
// constructor and destructor
BBSolver::BBSolver(Instance* inst,
Problem pb,
std::string optionfile,
float timelimit):
DiscretizationSolver(inst,
BRANCHANDBOUND,
pb,
NONE,
optionfile,
timelimit) {
this->isfeasible_ = false;
this->isoptimal_ = false;
this->bbnodes_ = 1;
this->treatednodes_ = 0;
this->nbactivenodes_ = 0;
this->primalbound_ = 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);
std::string optionname;
char buf[100];
int optionvalue;
std::string optionstring;
//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("relaxmodel") != std::string::npos) {
if (!(iss >> buf)){
std::cerr << "revorder: error: there is a mistake with the format "
"of the option file" <<std::endl;
throw;
}
optionstring = std::string(buf);
if (!strcmp(optionstring.c_str(),"ccg"))
{
this->modeltype_ = CCG;
}
else if (!strcmp(optionstring.c_str(),"cycles"))
{
this->modeltype_ = CYCLES;
}
else if (!strcmp(optionstring.c_str(),"ranks"))
{
this->modeltype_ = RANKS;
}
else if (!strcmp(optionstring.c_str(),"vertexrank"))
{
this->modeltype_ = VERTEXRANK;
}
else if (!strcmp(optionstring.c_str(),"witness"))
{
this->modeltype_ = WITNESS;
}
continue;
}
if (!(iss >> optionvalue)){
std::cerr << "revorder: error: there is a mistake with the format of "
"the option file" <<std::endl;
throw;
} // error
//
if (optionname.find("branchonsmallerindex") !=std::string::npos) {
branchonsmallerindex_ = (bool) optionvalue;
continue;
}
else if (optionname.find("prunesameneighbors") !=std::string::npos) {
prunesameneighbors_ = (bool) optionvalue;
continue;
}
else if (optionname.find("checkdominance") !=std::string::npos) {
if (optionname.find("checkdominancealltree") !=std::string::npos) {
checkdominancealltree_ = (bool) optionvalue;
continue;
}
checkdominance_ = (bool) optionvalue;
continue;
}
else if (optionname.find("maxdeltadominance") !=std::string::npos) {
maxdeltadominance_ = optionvalue;
continue;
}
else if (optionname.find("improvedualbound") !=std::string::npos) {
improvedualbound_ = (bool) optionvalue;
continue;
}
else if (optionname.find("userelaxbound") !=std::string::npos) {
userelaxbound_ = (bool) optionvalue;
continue;
}
else if (optionname.find("explorebest") !=std::string::npos) {
explorebest_ = (bool) optionvalue;
continue;
}
else if (optionname.find("exploredepth") !=std::string::npos) {
if (optionname.find("exploredepthbeforebest") !=std::string::npos) {
exploredepthbeforebest_ = (bool) optionvalue;
continue;
}
exploredepth_ = (bool) optionvalue;
continue;
}
}
}
else
{
//
// set some parameters for the solution of the relaxation in dual bounds
// computations
this->modeltype_ = CCG;
this->branchonsmallerindex_ = true;
this->prunesameneighbors_ = true;
this->checkdominance_ = true;
this->checkdominancealltree_ = true;
this->maxdeltadominance_ = 1;
this->improvedualbound_ = false;
this->userelaxbound_ = false;
this->explorebest_ = false;
this->exploredepth_ = false;
this->exploredepthbeforebest_ = true;
}
//
// build a root node
this->rootnode_ = new BBNode(*inst);
}
//
BBSolver::~BBSolver() {
cplexrelax_.end();
relax_.end();
env_.end();
cplexdual_.end();
modeldual_.end();
envdual_.end();
delete rootnode_;
}
//
// main method : called to run the branch-and-bound algorithm
int BBSolver::solve() {
//
// run timer
IloTimer cpuClockTotal(this->env_);
IloTimer cpuClockInit(this->env_);
cpuClockTotal.start();
cpuClockInit.start();
std::cout << std::endl;
std::cout << "revorder: "
"----------------------------------------------------------"
<< std::endl;
std::cout << "\nrevorder: bb: Run preprocessing procedures "
<< std::endl << std::endl;
//
// try and reduce the size of the instance right from the beginnning
this->inst_->computeneighbors();
this->preallocatelowdegreevertices();
//
// enumerate the potential initial cliques
if (!this->enumeratecliques(this->inst_->L() + 1)) {
this->isfeasible_ = false;
return false;
}
//
// eliminate the redundant and dominated cliques
this->eliminateredundantcliques();
//
// solve with greedy and remove cliques that cannot be completed
this->isfeasible_ = this->greedysolve();
if (!this->isfeasible_) return false;
//
// initialize the clique nodes
for (Clique* c : this->cliques_) {
BBNode* cliquenode =new BBNode(*inst_, rootnode_, *c,
this->branchonsmallerindex_);
this->bbnodesqueue_.push_back(cliquenode);
cliquenode->primalbound_ = c->greedyobjvalue_;
cliquenode->setid(this->bbnodes_);
cliquenode->setnodecost(this->problem_);
cliquenode->prepropagation();
cliquenode->postpropagation();
this->activenodes_[cliquenode->cost_].push_back(cliquenode);
this->nbactivenodes_++;
//
// first propagate the partial order
this->bbnodes_++;
}
std::cout << "revorder: bb: preprocessing created "
<< this->bbnodesqueue_.size()
<< " initial clique nodes" << std::endl;
//
// sort the bb nodes according to the choice of exploration method
if (this->exploredepth_ || this->exploredepthbeforebest_) {
std::make_heap(this->bbnodesqueue_.begin(), this->bbnodesqueue_.end(),
bbnodesnborderedcompare);
} else if (this->explorebest_) {
std::make_heap(this->bbnodesqueue_.begin(), this->bbnodesqueue_.end(),
bbnodesbestcompare);
} else {
std::cout << "revorder: bb: need to choose an exploration method in "
"branch-and-bound" << std::endl;
throw;
}
//
#ifdef VERBOSE
std::cout << "revorder: bb: greatest number of ordered vertices in a clique: ";
std::cout << ", " << this->bbnodesqueue_.front()->nbordered() << " vertices" << std::endl;
#endif
//
// compute sets of vertices that must 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();
//
// initialize the IP used for improved dual bound
this->initializedualboundIP();
//
// set greedy solution as initial primal bound
this->primalbound_ = this->objvalue_;
//
cpuClockInit.stop();
// Run the branch-and-bound algorithm
//
// initialize the ip model to find dual bounds using the linear relaxation
// IloModel model(this->env_);
this->relax_ = IloModel(this->env_);
// this->defineminpartialip(model);
if (userelaxbound_) {
// this->createrelaxationmodel(model, this->relax_);
this->defineminpartialip(this->relax_);
this->cplexrelax_ = IloCplex(this->env_);
this->cplexrelax_.extract(this->relax_);
}
//
// treat
std::cout << std::endl;
std::cout << "revorder: "
"----------------------------------------------------------"
<< std::endl;
std::cout << "\nrevorder: bb: Start the B&B algorithm " << std::endl
<< std::endl;
while (!this->bbnodesqueue_.empty()) {
this->treatnode(*(this->bbnodesqueue_.front()));
if (cpuClockTotal.getTime() > this->timelimit_) {
std::cout << "revorder: bb: time limit has been reached, stop the "
"solution algorithm" << std::endl;
break;
}
}
cpuClockTotal.stop();
// Display the results
this->totaltime_ = cpuClockTotal.getTime();
this->relaxtime_ = 0.0;
std::cout << std::endl;
std::cout << "revorder: "
"----------------------------------------------------------"
<< std::endl;
std::cout << "revorder: "
"----------------------------------------------------------"
<< std::endl;
std::cout << "\nmdjeep: Branch-and-bound solution report: " << std::endl;
if (isfeasible_) {
this->objvalue_ = this->primalbound_;
// the solution is optimal only if all the bb nodes have been treated
if (this->bbnodesqueue_.empty()) {
isoptimal_ = true;
std::cout << "revorder: the instance was solved to optimality"
<< std::endl;
} else {
std::cout << "revorder: the branch-and-bound was stopped before proving"
" optimality" << std::endl;
}
// reconstruct the solution by including the preallocated vertices
objvalue_ += this->inst_->preallocatedvertices_.size();
this->reconstructsolution();
//
// display the solution
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 explored nodes = " << this->treatednodes_ << std::endl;
std::cout << "revorder: number of created nodes = " << this->bbnodes_ << std::endl;
std::cout << "revorder: value of the objective function = " << this->objvalue_ << " (after addition of preallocated vertices)" << std::endl;
//
// verify the resulting order is indeed a revorder
this->verifyorder(this->bestrank_);
}
else {
std::cout << "revorder: the instance is not discretizable!" << std::endl;
}
return isfeasible_;
}
//
// treat a branch-and-bound node (compute bounds and branch)
void BBSolver::treatnode(BBNode& node) {
this->treatednodes_++;
if (std::remainder(this->treatednodes_, 100) == 0) {
std::cout << "revorder: bb: treated nodes = " << this->treatednodes_
<< ", nodes in the queue = " << this->bbnodesqueue_.size()
<< ", nodes in memory = " << this->nbactivenodes_
<< ", primal bound = " << this->primalbound_ << std::endl;
}
//
// remove the bb node from the list
//
// make sure to maintain the heap first
if (this->isfeasible_) {
if (this->exploredepth_) {
std::pop_heap(this->bbnodesqueue_.begin(),
this->bbnodesqueue_.end(), bbnodesnborderedcompare);
} else if (this->explorebest_ || this->exploredepthbeforebest_) {
std::pop_heap(this->bbnodesqueue_.begin(),
this->bbnodesqueue_.end(), bbnodesbestcompare);
}
} else {
if (this->explorebest_) {
std::pop_heap(this->bbnodesqueue_.begin(),
this->bbnodesqueue_.end(), bbnodesbestcompare);
} else {
std::pop_heap(this->bbnodesqueue_.begin(),
this->bbnodesqueue_.end(), bbnodesnborderedcompare);
}
}
// actually remove the node from the queue
this->bbnodesqueue_.pop_back();
// update the primal bound and prune the node if it is a leaf of the
// enumeration tree
if (node.nbordered() == this->inst_->nbvertices_) {
#ifdef VERBOSE
std::cout << "revorder: bb: reached a leaf at node " << node.id() << ": depth = " << node.depth();
#endif
node.primalbound_ = this->getobjvalue(node.rank_, node.nbrefs_);
//
// update the primal bound if improved
if (node.primalbound_ < this->primalbound_) {
//
// swap the search method if this is the first feasible solution
std::cout << "revorder: bb: the primal bound has been improved, value = "
<< node.primalbound_ << std::endl;
if (!this->isfeasible_) {
if (this->exploredepthbeforebest_) {
std::make_heap(this->bbnodesqueue_.begin(),
this->bbnodesqueue_.end(), bbnodesbestcompare);
}
this->isfeasible_ = true;
}
this->primalbound_ = node.primalbound_;
this->bestnbfullref_ = this->getnbfullref(node.rank_, node.nbrefs_);
for (Vertex* v : inst_->vertices_) {
this->bestrank_[v] = node.rank(v);
}
}
this->erasenodefromall(&node);
} else {
// otherwise, compute dual and primal bound and prepare for next node
// prune the node if no potential child node
if (node.potentialchildren_.empty()) {
//
#ifdef VERBOSE
std::cout << "revorder: bb: prune node " << node.id() << ": ";
std::cout << ": no feasible solution on this branch" << std::endl;
#endif
//
this->erasenodefromall(&node);
} else {
// compute a dual bound based on the ordered vertices
this->computedualbound(node);
//
// prune the node by using bounds
if (node.dualbound_ >= this->primalbound_) {
#ifdef VERBOSE
std::cout << "revorder: bb: prune node " << node.id() << ": ";
std::cout << "using bounds: primal = " << this->primalbound_;
std::cout << " ; dual = " << node.dualbound_ << std::endl;
#endif
this->erasenodefromall(&node);
} else {
// otherwise, create new nodes by branching
this->branch(node);
node.istreated_ = true;
// if (!this->checkdominancealltree_) {
// this->erasenodefromall(&node);
// }
}
}
}
//
// treat next node in the queue
// since it is a heap, the front node is always the greatest in the
// implemented order
#ifdef VERBOSE
if (!this->bbnodesqueue_.empty()) {
BBNode* nextnode = this->bbnodesqueue_.front();
std::cout << '\n' << "revorder: bb: treat node " << nextnode->id() << " ; depth = " << nextnode->depth() << " ; ";
std::cout << nextnode->nbordered() << " ordered vertices ; " << nextnode->nbpartialref() << " partially referenced vertices" << std::endl;
}
#endif
}
// create new nodes by branching
void BBSolver::branch(BBNode& node) {
// create the children of the input node
// propagate the order starting from each incomplete order to detect
// dominatedorders
std::vector<BBNode*> childnodes;
std::map<BBNode*, bool> isdominated;
std::map<Vertex*, bool> issymmetric;
//
// we start by checking diverse pruning rules that need to be applied
// before dominance
for (int i = 0; i < node.potentialchildren_.size() ; i++) {
Vertex* u = node.potentialchildren_[i];
issymmetric[u] = false;
if (!this->prunesameneighbors_) continue;
for (int j = 0; j < i; j++) {
Vertex* v = node.potentialchildren_[j];
if (issymmetric[v]) continue;
// if two vertices that can be chosen for branching have the exact same
// unordered neighbors, only that with smaller number of references
// needs to be considered for branching at this stage
bool sameneighbors = true;
for (Vertex* neighbor : u->neighbors_) {
if (!node.isordered(neighbor)) {
if ( !v->isneighbor_[neighbor] ) {
sameneighbors = false;
break;
}
}
}
if (sameneighbors) {
#ifdef VERBOSE
std::cout << "revorder: bb: prune node: two potential children have same neighbors: vertices " << u->id_ << " and " << v->id_ << std::endl;
#endif
if (node.nbrefs_[u] < node.nbrefs_[v]) {
issymmetric[v] = true;
} else if (node.nbrefs_[v] < node.nbrefs_[u]) {
issymmetric[u] = true;
} else {
if (u->id_ < v->id_) issymmetric[v] = true;
else if (v->id_ < u->id_) issymmetric[u] =true;
}
}
}
}
for (Vertex* v : node.potentialchildren_) {
if (issymmetric[v]) continue; // do not consider the vertices that are
// symmetric
BBNode* child = new BBNode(node);
//
// in MINPARTIAL, we can always make the arbitrary decision that if a
// vertex must be ordered when it is still partially referenced, then it
// is the vertex with smallest index among all the potential children; as
// a consequence the vertices with smaller index must be removed from the
// list of potential children
if ( this->branchonsmallerindex_ ) {
std::vector<Vertex*>::iterator itv = child->potentialchildren_.begin();
while (itv < child->potentialchildren_.end()) {
if ( ((*itv)->id_ < v->id_) ) {
child->ispotentialchild_[*itv] = false;
child->potentialchildren_.erase(itv);
child->mustbefullref_[*itv] = true;
//
// if a vertex must be fully-referenced, but it has exactly U
// neighbors, it will for sure come after all its neighbors in the
// order
if (child->maxrefs_[*itv] <= this->inst_->U()) {
for (Vertex* u : v->neighbors_) {
if (child->mustbefullref_[u]) {
if (child->maxrefs_[u] >= this->inst_->U()+1) {
child->maxrefs_[u]--;
}
} else {
child->maxrefs_[u]--;
}
}
}
} else {
itv++;
}
}
}
child->assignnextrank(v);
child->prepropagation();
isdominated[child] = false;
childnodes.push_back(child);
}
//
// check for classical dominance
for (BBNode* n1 : childnodes) {
if (!this->checkdominance_) continue; // do not check dominance if option
// is deactivated
if (isdominated[n1]) continue;
for (BBNode* n2 : childnodes) {
if (n2 == n1) {
continue;
} else if (isdominated[n2]) {
continue;
} else if (this->dominates(*n2, *n1)) {
// classical domination
isdominated[n1] = true;
//
#ifdef VERBOSE
std::cout << "revorder: bb: prune node: dominance among children" << std::endl;
#endif
break;
} else if (this->dominates(*n1, *n2)) {
isdominated[n2] = true;
//
#ifdef VERBOSE
std::cout << "revorder: bb: prune node: dominance among children" << std::endl;
#endif
}
}
}
//
// scan all the bb nodes in the queue to delete the dominated ones if
// option is activated
if ((this->checkdominancealltree_) && (this->checkdominance_)) {
for (BBNode* n : childnodes) {
if (isdominated[n]) continue;
for (int deltacost = -this->maxdeltadominance_ ;
deltacost <= this->maxdeltadominance_; deltacost++) {
if (n->cost_ <= deltacost) continue;
if (!isdominated[n]) {
isdominated[n] =
checkdominancewithlist(n,
this->activenodes_[n->cost_ - deltacost]);
}
}
}
}
//
// delete the redundant nodes and add the others to the node queue
for (BBNode* n : childnodes) {
if (isdominated[n]) {
delete n;
} else {
// make sure to maintain the heap while pushing the new node in the queue
this->bbnodes_++;
n->setid(this->bbnodes_);
node.addtochildren(n);
n->postpropagation();
this->bbnodesqueue_.push_back(n);
this->activenodes_[n->cost_].push_back(n);
this->nbactivenodes_++;
if (this->exploredepth_) {
std::push_heap(this->bbnodesqueue_.begin(),
this->bbnodesqueue_.end(), bbnodesnborderedcompare);
} else if (this->explorebest_) {
std::push_heap(this->bbnodesqueue_.begin(),
this->bbnodesqueue_.end(), bbnodesbestcompare);
} else if (this->exploredepthbeforebest_) {
if (!this->isfeasible_) {
std::push_heap(this->bbnodesqueue_.begin(),
this->bbnodesqueue_.end(), bbnodesnborderedcompare);
} else {
std::push_heap(this->bbnodesqueue_.begin(),
this->bbnodesqueue_.end(), bbnodesbestcompare);
}
}
}
}
#ifdef VERBOSE
std::cout << "revorder: bb: " << this->bbnodesqueue_.size() << " nodes in the queue" << std::endl;
#endif
}
//
// check the dominance of input node with every node in the input vector
bool BBSolver::checkdominancewithlist(BBNode* n1,
std::vector<BBNode*>& nodelist) {
std::vector<BBNode*>::iterator itnode = nodelist.begin();
while (itnode < nodelist.end()) {
BBNode* n2 = (*itnode);
if (this->dominates(*n2, *n1)) {
//
#ifdef VERBOSE
std::cout << "revorder: bb: prune node: dominance of potential child ";
if (n2->istreated_) std::cout << " with treated node" << std::endl;
else std::cout << "with queue " << std::endl;
std::cout << "revorder: bb: depth: dominant = " << n2->depth() << ", dominated = " << n1->depth() << std::endl;
std::cout << "revorder: bb: costs: dominant = " << n2->cost_ << ", dominated = " << n1->cost_ << std::endl;
#endif
//
return true;
} else if (this->dominates(*n1, *n2)) {
//
#ifdef VERBOSE
std::cout << "revorder: bb: prune node: dominance of created node ";
if (n2->istreated_) std::cout << " already treated" << std::endl;
else std::cout << "still in queue" << std::endl;
std::cout << "revorder: bb: depth: dominant = " << n1->depth() << ", dominated = " << n2->depth() << std::endl;
std::cout << "revorder: bb: costs: dominant = " << n1->cost_ << ", dominated = " << n2->cost_ << std::endl;
#endif
//
erasenodefromall(*itnode);
} else {
itnode++;
}
}
return false;
}
//
// erase a node from every list where it appears
void BBSolver::erasenodefromall(BBNode* node) {
// first erase it from the list of children of its father
for (auto itn = node->getfather()->children_.begin();
itn < node->getfather()->children_.end(); itn++) {
if (*itn == node) {
node->getfather()->children_.erase(itn);
break;
}
}
// then erase its descendants if any
for (BBNode* child: node->children_) {
this->erasealldescendants(child);
}
// erase from active nodes
for (auto itn = activenodes_[node->cost_].begin();
itn < activenodes_[node->cost_].end(); itn++) {
if (*itn == node) {
activenodes_[node->cost_].erase(itn);
this->nbactivenodes_--;
break;
}
}
if (!node->istreated_) {
for (auto itn = bbnodesqueue_.begin(); itn < bbnodesqueue_.end(); itn++) {
if (*itn == node) {
bbnodesqueue_.erase(itn);
break;
}
}
}
delete node;
}
// erase all the descendants of a node from all list of nodes and delete them
void BBSolver::erasealldescendants(BBNode* node) {
for (BBNode* child: node->children_) {
erasealldescendants(child);
}
for (auto itn = activenodes_[node->cost_].begin(); itn < activenodes_[node->cost_].end(); itn++) {
if (*itn == node) {
activenodes_[node->cost_].erase(itn);
this->nbactivenodes_--;
break;
}
}
if (!node->istreated_) {
for (auto itn = bbnodesqueue_.begin(); itn < bbnodesqueue_.end(); itn++) {
if (*itn == node) {
bbnodesqueue_.erase(itn);
break;
}
}
}
delete node;
}
// compute a dual bound from the partial order described in the input node
double BBSolver::computedualbound(BBNode& node) {
//
// compute a simple bound by just looking at the ordered vertices
//
// first get the number of partially referenced vertices in the order
int nbpartials = node.nbordered() - node.nbfullref() - this->inst_->L();
//
// then aknowledge that vertices with less than U neighbors cannot be
// fully-referenced; but beware that we know that one potential child will
// be partially-referenced, so do not count it twice
int nbpotentialpartials = 1;
for (Vertex* v: this->inst_->vertices_) {
if (!node.isordered(v)) {
if (v->degree_ <= this->inst_->U() - 1) {
if (node.ispotentialchild_[v]) {
nbpartials++;
nbpotentialpartials = 0;
}
else{
nbpartials++;
}
}
}
}
// finally for each edge linking two vertices with exactly U neigbors, one
// will be partially referenced; still need to beware with the potential
// children of the node
std::map<Vertex*,bool> inpartialedge;
for (Vertex* u: this->inst_->vertices_) inpartialedge[u] = false;
for (Vertex* u: this->inst_->vertices_) {
if (!node.isordered(u)
&& (!inpartialedge[u])
&& (u->degree_ == this->inst_->U())) {
for (Vertex* v: u->neighbors_) {
if (!node.isordered(v)
&& (!inpartialedge[v])
&& (v->degree_ == this->inst_->U())) {
nbpartials++;
inpartialedge[u] = true;
inpartialedge[v] = true;
if (node.ispotentialchild_[v]) nbpotentialpartials = 0;
break;
}
}
}
}
nbpartials += nbpotentialpartials;
int trivialbound = nbpartials;
//
// update current dual bound if improved
node.dualbound_ = std::max(node.dualbound_, trivialbound);
if ((!this->improvedualbound_) && (!this->userelaxbound_) ) {
return node.dualbound_;
}
/////////////////////////////////////////////////////////////////////////////
// Improve the trivial bound by considering several cuts
/////////////////////////////////////////////////////////////////////////////
if (this->improvedualbound_) {
for (Vertex* v : this->inst_->vertices_) {
if (node.isordered(v)) {
if (node.nbrefs(v) >= this->inst_->U()) {
this->ispartial_[v].setBounds(0, 0);
}
else {
this->ispartial_[v].setBounds(1,1);
}
}
else if (node.mustbefullref_[v]) {
this->ispartial_[v].setBounds(0, 0);
}
else if (v->degree_ <= this->inst_->U() - 1) {
this->ispartial_[v].setBounds(1,1);
}
else {
this->ispartial_[v].setBounds(0, 1);
}
}
//
// at least one among the potential children will be partially referenced
IloExpr sumpartialinpotentials(this->envdual_);
for (Vertex* v : node.potentialchildren_) {
sumpartialinpotentials += this->ispartial_[v];
}
IloRange ctPartialsInPotentials(this->envdual_, -sumpartialinpotentials, -1);
this->modeldual_.add(ctPartialsInPotentials);
sumpartialinpotentials.end();
//
// solve the IP
this->cplexdual_.solve();
IloAlgorithm::Status statusimproved = this->cplexdual_.getStatus();
int improvedbound = this->inst_->nbvertices_;
if (statusimproved==IloAlgorithm::Infeasible) {
node.dualbound_ = this->inst_->nbvertices_;
#ifdef VERBOSE
std::cout << "revorder: bb: improved dual bound IP is infeasible: prune the node" << std::endl;
#endif
}
else {
improvedbound = ceil(this->cplexdual_.getObjValue()) - this->inst_->L();
if (improvedbound > node.dualbound_ ){
node.dualbound_ = improvedbound;
}
}
this->modeldual_.remove(ctPartialsInPotentials);
ctPartialsInPotentials.end();
#ifdef VERBOSE
std::cout << "revorder: node dual bound = " << improvedbound << " ; trivialbound = " << trivialbound << std::endl;
#endif
}
/////////////////////////////////////////////////////////////////////////////
// Compute a dual bound by solving the LP relaxation of a chosen model
/////////////////////////////////////////////////////////////////////////////
// Start searching for better bounds only if trivial bound is not too far from primal bound and a significant number of nodes have been treated to improve primal bound
if ( (this->userelaxbound_) && (this->treatednodes_ >= 1000) && (trivialbound * 2.0 >= this->primalbound_)) {
// fix all the variables related to the vertices that already in the current incomplete order
for (Vertex* v1 : this->inst_->vertices_) {
if (node.rank(v1) >= this->inst_->L() + 2) {
if (node.nbrefs(v1) >= this->inst_->U()) {
isfullref_[v1].setBounds(1, 1);
}
else {
isfullref_[v1].setBounds(0, 0);
}
}
else if ((node.rank(v1) >= 1) && (node.rank(v1) <= this->inst_->L()+1)) {
isfullref_[v1].setBounds(1, 1);
}
else {
if (node.mustbefullref_[v1] == true) {
isfullref_[v1].setBounds(1,1);
}
else {
isfullref_[v1].setBounds(0, 1);
}
}
for (Vertex* v2 : v1->neighbors_) {
if ((node.rank(v1) < 0) && (node.rank(v2) < 0)) isbefore_[v1][v2].setBounds(0, 1);
else if ( (node.rank(v1) < 0) && (node.rank(v2) >= 1)) isbefore_[v1][v2].setBounds(0, 0);
else if ((node.rank(v1) >= 1) && (node.rank(v2) < 0)) isbefore_[v1][v2].setBounds(1, 1);
else if ((node.rank(v1) >= 1) && (node.rank(v2) >= 1) && (node.rank(v1) < node.rank(v2)) ) isbefore_[v1][v2].setBounds(1, 1);
else if ((node.rank(v1) >= 1) && (node.rank(v2) >= 1) && (node.rank(v1) > node.rank(v2)) ) isbefore_[v1][v2].setBounds(0, 0);
else {
cout << node.rank(v1) << "; " << node.rank(v2) << endl;
throwError("abnormal rank values");
}
}
}
IloExpr sumnotallfullref(this->env_);
int maxnbfullref;
for (Vertex* v : node.potentialchildren_) {
sumnotallfullref += isfullref_[v];
}
maxnbfullref = node.potentialchildren_.size() - 1;
IloRange cons(this->env_, sumnotallfullref, maxnbfullref);
relax_.add(cons);
cplexrelax_.setParam(IloCplex::SimDisplay, 0);
// cplexrelax_.setParam(IloCplex::Param::ParamDisplay, 0);
cplexrelax_.setParam(IloCplex::Threads, 1);
cplexrelax_.setParam(IloCplex::MIPDisplay, 0);
cplexrelax_.setParam(IloCplex::TuningDisplay, 0);
cplexrelax_.setOut(cplexrelax_.getEnv().getNullStream());
bool updateprimal_ = false;
if ( updateprimal_ && (node.nbordered() >= 0.7 * this->inst_->nbvertices_) ) {
cplexrelax_.setParam(IloCplex::NodeLim, 9223372036800000000);
cplexrelax_.setParam(IloCplex::MIPSearch, 1);
if ((this->modeltype_ == CCG) || (this->modeltype_ == WITNESS)) {
cplexrelax_.use(CyclesLazyConstraintsBB(cplexrelax_.getEnv(), *(this->inst_), *this));
}
}
else {
cplexrelax_.setParam(IloCplex::NodeLim, 0);
}
cplexrelax_.solve();
IloAlgorithm::Status status = cplexrelax_.getStatus();
//
// if the relaxation is infeasible, the node can be pruned : set dual bound to negative value to state this
if (status != IloAlgorithm::Infeasible) {
int relaxbound = ceil(cplexrelax_.getBestObjValue());
//
#ifdef VERBOSE
std::cout << "revorder: bb: trivial dual bound = " << trivialbound << " ; relaxation bound = " << relaxbound;
std::cout << " (best primal bound = " << this->primalbound_ << ")" << std::endl;
#endif
if ( updateprimal_ && (node.nbordered() >= 0.7 * this->inst_->nbvertices_) ) {
if (relaxbound < this->primalbound_) {
this->primalbound_ = relaxbound;
std::cout << "best primal bound updated = " << this->primalbound_ << std::endl;
this->bestnbfullref_ = this->inst_->nbvertices_ - this->primalbound_;
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 (cplexrelax_.getValue(this->isbefore_[u][v]) >= 1 - tolerance) {
adjlist[u].push_back(v);
}
}
}
//
// search for the root of the digraph
Vertex* root = nullptr;
for (Vertex* v : this->inst_->vertices_) {
if (node.rank(v) == 1) {
root = v;
}
}
// 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;
}
for (Vertex* v : this->inst_->vertices_) {
this->bestrank_[v] = rank[v];
}
}
}
//
// update current dual bound if improved
if (relaxbound > node.dualbound_) {
node.dualbound_ = relaxbound;
//
#ifdef VERBOSE
std::cout << "revorder: bb: new dual bound = " << relaxbound << std::endl;
std::cout << "revorder: bb: number of ordered vertices = " << node.nbordered() << std::endl;
#endif
//
}
}
else{
node.dualbound_ = this->inst_->nbvertices_;
#ifdef VERBOSE
std::cout << "revorder: bb: dual bound relaxation is infeasible: prune the node" << std::endl;
#endif
}
relax_.remove(cons);
cons.end();
}
return node.dualbound_;
}
// initialize the IP model for improved dual bound
void BBSolver::initializedualboundIP() {
this->modeldual_ = IloModel(this->envdual_);
this->cplexdual_ = IloCplex(this->modeldual_);
IloExpr obj(this->envdual_);
char name[256];
for (Vertex* v : this->inst_->vertices_) {
this->ispartial_[v] = IloBoolVar(this->envdual_);
sprintf(name, "ispartial%i", v->id_);
this->ispartial_[v].setName(name);
obj += this->ispartial_[v];
}
this->modeldual_.add(IloMinimize(this->envdual_, obj));
obj.end();
//
// 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
for (Clique* c : this->cliqueswithpartialref_) {
IloExpr sumpartialsinclique(this->envdual_);
for (Vertex* v : c->vertices_) sumpartialsinclique += this->ispartial_[v];
this->modeldual_.add(sumpartialsinclique >= this->nbpartialinclique_[c]);
sumpartialsinclique.end();
}
//
// 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
IloRangeArray ctaryLowDegreeCycles(this->envdual_);
int nbcons = 0;
for (int i = 0; i < this->lowdegreecycles_.size(); i++) {
std::vector<Vertex*> cycle = this->lowdegreecycles_[i];
IloExpr sumpartialsincycle(this->envdual_);
for (Vertex* v : cycle) {
sumpartialsincycle += this->ispartial_[v];
}
ctaryLowDegreeCycles.add(sumpartialsincycle >= this->nbpartialsincycle_[i]);
sprintf(name, "ctaryLowDegreeCycles_%i", nbcons);
ctaryLowDegreeCycles[nbcons++].setName(name);
sumpartialsincycle.end();
}
this->modeldual_.add(ctaryLowDegreeCycles);
ctaryLowDegreeCycles.end();
//
// load the integer problem and set display parameters
this->cplexdual_.setParam(IloCplex::MIPDisplay, 0);
this->cplexdual_.setParam(IloCplex::SimDisplay, 0);
this->cplexdual_.setParam(IloCplex::TiLim, 10);
this->cplexdual_.setParam(IloCplex::Threads, 1);
// this->cplexdual_.setParam(IloCplex::Param::ParamDisplay, 0);
}