bbsolver.cpp

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