/********************************************************************************************
Name: DiscretizationSolver
Discrete optimization for graph discretization
Author: J.Omer
Sources: C++
License: GNU General Public License v.2
History:
*********************************************************************************************/
#include "bbsolver.hpp"
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 ((float) n1->cost_ / n1->nbordered() > (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);
}
//
// add partially referenced vertices to the list of potential choices for branching
else if (this->nbrefs_[v] >= inst.L()) {
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
}
}
//
// 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 three
// references or exactly one unordered vertex with exactly three 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;
//
// add a vertex with more than U references
if (!this->fullrefs_.empty()) {
u = this->fullrefs_.back();
this->fullrefs_.pop_back();
this->isfullref_[u] = false;
ispropagate = true;
}
//
// if there is a partially referenced vertex with less than U neighbors, it can be ordered right away: it will never be fully referenced
else if (!this->readytoorder_.empty()) {
u = this->readytoorder_.back();
this->readytoorder_.pop_back();
ispropagate = true;
}
//
// if there is exactly one partially referenced vertex: add it to the order
else if (this->potentialchildren_.size() == 1) {
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);
}
//
// otherwise, compute dual and primal bound and prepare for next node
else {
//
// 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);
}
//
// otherwise, create new nodes by branching
else {
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
// make_heap pour bien ranger les noeuds de branch-and-bound au fur et a mesure
}
//
// propagate an initial clique by iteratively adding the fully-referenced vertices until
// one is necessarily free
void BBSolver::branch(BBNode& node) {
//
// create the children of the input node
// propagate the order starting from each incomplete order to detect dominated
// orders
std::vector<BBNode*> childrennodes;
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;
childrennodes.push_back(child);
}
//
// check for classical dominance
for (BBNode* n1 : childrennodes) {
if (!this->checkdominance_) continue; // do not check dominance if option is deactivated
if (isdominated[n1]) continue;
for (BBNode* n2 : childrennodes) {
if (n2 == n1) continue;
else if (isdominated[n2]) continue;
//
// classical dominance
else if (this->dominates(*n2,*n1)) {
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 : childrennodes) {
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 : childrennodes) {
if (isdominated[n]) delete n;
//
// make sure to maintain the heap while pushing the new node in the queue
else {
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);
}