bbsolver.cpp

/********************************************************************************************
 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