bbsolver.cpp

  }
  //
  // scan all the bb nodes in the queue to delete the dominated ones if
  // option is activated
  if ((this->checkdominancealltree_) && (this->checkdominance_)) {
    for (BBNode* n : childnodes) {
      if (isdominated[n]) continue;
      for (int deltacost = -this->maxdeltadominance_ ;
           deltacost <= this->maxdeltadominance_; deltacost++) {
        if (n->cost_ <= deltacost) continue;
        if (!isdominated[n]) {
          isdominated[n] =
              checkdominancewithlist(n,
                                     this->activenodes_[n->cost_ - deltacost]);
        }
      }
    }
  }
  //
  // delete the redundant nodes and add the others to the node queue
  for (BBNode* n : childnodes) {
    if (isdominated[n]) {
      delete n;
    } else {
      // make sure to maintain the heap while pushing the new node in the queue
      this->bbnodes_++;
      n->setid(this->bbnodes_);
      node.addtochildren(n);
      n->postpropagation();
      this->bbnodesqueue_.push_back(n);
      this->activenodes_[n->cost_].push_back(n);
      this->nbactivenodes_++;
      if (this->exploredepth_) {
        std::push_heap(this->bbnodesqueue_.begin(),
                       this->bbnodesqueue_.end(), bbnodesnborderedcompare);
      } else if (this->explorebest_) {
        std::push_heap(this->bbnodesqueue_.begin(),
                       this->bbnodesqueue_.end(), bbnodesbestcompare);
      } else if (this->exploredepthbeforebest_) {
        if (!this->isfeasible_) {
          std::push_heap(this->bbnodesqueue_.begin(),
                         this->bbnodesqueue_.end(), bbnodesnborderedcompare);
        } else {
          std::push_heap(this->bbnodesqueue_.begin(),
                         this->bbnodesqueue_.end(), bbnodesbestcompare);
        }
      }
    }
  }
#ifdef VERBOSE
  std::cout << "revorder: bb: " << this->bbnodesqueue_.size() << " nodes in the queue" << std::endl;
#endif
}
//
// check the dominance of input node with every node in the input vector
bool BBSolver::checkdominancewithlist(BBNode* n1,
                                      std::vector<BBNode*>& nodelist) {
  std::vector<BBNode*>::iterator itnode = nodelist.begin();

  while (itnode < nodelist.end()) {
    BBNode* n2 = (*itnode);
    if (this->dominates(*n2, *n1)) {
      //
#ifdef VERBOSE
      std::cout << "revorder: bb: prune node: dominance of potential child ";
            if (n2->istreated_) std::cout << " with treated node" <<  std::endl;
            else std::cout << "with queue " << std::endl;
            std::cout << "revorder: bb: depth: dominant = " << n2->depth() << ", dominated = " << n1->depth() << std::endl;
            std::cout << "revorder: bb: costs: dominant = " << n2->cost_ << ", dominated = " << n1->cost_ << std::endl;
#endif
      //
      return true;
    } else if (this->dominates(*n1, *n2)) {
      //
#ifdef VERBOSE
      std::cout << "revorder: bb: prune node: dominance of created node ";
            if (n2->istreated_) std::cout << " already treated" <<  std::endl;
            else std::cout << "still in queue" << std::endl;
            std::cout << "revorder: bb: depth: dominant = " << n1->depth() << ", dominated = " << n2->depth() << std::endl;
            std::cout << "revorder: bb: costs: dominant = " << n1->cost_ << ", dominated = " << n2->cost_ << std::endl;
#endif
      //
      erasenodefromall(*itnode);
    } else {
      itnode++;
    }
  }
  return false;
}

//
// erase a node from every list where it appears
void BBSolver::erasenodefromall(BBNode* node) {
  // first erase it from the list of children of its father
  for (auto itn = node->getfather()->children_.begin();
       itn < node->getfather()->children_.end(); itn++) {
    if (*itn == node) {
      node->getfather()->children_.erase(itn);
      break;
    }
  }

  // then erase its descendants if any
  for (BBNode* child: node->children_) {
    this->erasealldescendants(child);
  }

  // erase from active nodes
  for (auto itn = activenodes_[node->cost_].begin();
       itn < activenodes_[node->cost_].end(); itn++) {
    if (*itn == node) {
      activenodes_[node->cost_].erase(itn);
      this->nbactivenodes_--;
      break;
    }
  }
  if (!node->istreated_) {
    for (auto itn = bbnodesqueue_.begin(); itn < bbnodesqueue_.end(); itn++) {
      if (*itn == node) {
        bbnodesqueue_.erase(itn);
        break;
      }
    }
  }
  delete node;
}

// erase all the descendants of a node from all list of nodes and delete them
void BBSolver::erasealldescendants(BBNode* node) {

  for (BBNode* child: node->children_) {
    erasealldescendants(child);
  }
  for (auto itn = activenodes_[node->cost_].begin(); itn < activenodes_[node->cost_].end(); itn++) {
    if (*itn == node) {
      activenodes_[node->cost_].erase(itn);
      this->nbactivenodes_--;
      break;
    }
  }
  if (!node->istreated_) {
    for (auto itn = bbnodesqueue_.begin(); itn < bbnodesqueue_.end(); itn++) {
      if (*itn == node) {
        bbnodesqueue_.erase(itn);
        break;
      }
    }
  }
  delete node;
}

// compute a dual bound from the partial order described in the input node
double BBSolver::computedualbound(BBNode& node) {
  //
  // compute a simple bound by just looking at the ordered vertices
  //
  // first get the number of partially referenced vertices in the order
  int nbpartials = node.nbordered() - node.nbfullref() - this->inst_->L();
  //
  // then aknowledge that vertices with less than U neighbors cannot be
  // fully-referenced; but beware that we know that one potential child will
  // be partially-referenced, so do not count it twice
  int nbpotentialpartials = 1;
  for (Vertex* v: this->inst_->vertices_) {
    if (!node.isordered(v)) {
      if (v->degree_ <= this->inst_->U() - 1) {
        if (node.ispotentialchild_[v]) {
          nbpartials++;
          nbpotentialpartials = 0;
        }
        else{
          nbpartials++;
        }
      }
    }
  }

  // finally for each edge linking two vertices with exactly U neigbors, one
  // will be partially referenced; still need to beware with the potential
  // children of the node
  std::map<Vertex*,bool> inpartialedge;
  for (Vertex* u: this->inst_->vertices_) inpartialedge[u] = false;
  for (Vertex* u: this->inst_->vertices_) {
    if (!node.isordered(u)
        && (!inpartialedge[u])
        && (u->degree_ == this->inst_->U())) {
      for (Vertex* v: u->neighbors_) {
        if (!node.isordered(v)
            && (!inpartialedge[v])
            && (v->degree_ == this->inst_->U())) {
          nbpartials++;
          inpartialedge[u] = true;
          inpartialedge[v] = true;
          if (node.ispotentialchild_[v]) nbpotentialpartials = 0;
          break;
        }
      }
    }
  }
  nbpartials += nbpotentialpartials;
  int trivialbound = nbpartials;
  //
  // update current dual bound if improved
  node.dualbound_ = std::max(node.dualbound_, trivialbound);


  if ((!this->improvedualbound_) && (!this->userelaxbound_) ) {
    return node.dualbound_;
  }

  /////////////////////////////////////////////////////////////////////////////
  // Improve the trivial bound by considering several cuts
  /////////////////////////////////////////////////////////////////////////////

  if (this->improvedualbound_) {
    for (Vertex* v : this->inst_->vertices_) {
      if (node.isordered(v)) {
        if (node.nbrefs(v) >= this->inst_->U()) {
          this->ispartial_[v].setBounds(0, 0);
        }
        else {
          this->ispartial_[v].setBounds(1,1);
        }
      }
      else if (node.mustbefullref_[v]) {
        this->ispartial_[v].setBounds(0, 0);
      }
      else if (v->degree_ <= this->inst_->U() - 1) {
        this->ispartial_[v].setBounds(1,1);
      }
      else {
        this->ispartial_[v].setBounds(0, 1);
      }
    }
    //
    // at least one among the potential children will be partially referenced
    IloExpr sumpartialinpotentials(this->envdual_);
    for (Vertex* v : node.potentialchildren_) {
      sumpartialinpotentials += this->ispartial_[v];
    }
    IloRange ctPartialsInPotentials(this->envdual_, -sumpartialinpotentials, -1);
    this->modeldual_.add(ctPartialsInPotentials);
    sumpartialinpotentials.end();
    //
    // solve the IP
    this->cplexdual_.solve();
    IloAlgorithm::Status statusimproved = this->cplexdual_.getStatus();
    int improvedbound =  this->inst_->nbvertices_;
    if (statusimproved==IloAlgorithm::Infeasible) {
      node.dualbound_ =  this->inst_->nbvertices_;
#ifdef VERBOSE
      std::cout << "revorder: bb: improved dual bound IP is infeasible: prune the node" << std::endl;
#endif
    }
    else {
      improvedbound = ceil(this->cplexdual_.getObjValue()) - this->inst_->L();
      if (improvedbound > node.dualbound_ ){
        node.dualbound_ =  improvedbound;
      }
    }

    this->modeldual_.remove(ctPartialsInPotentials);
    ctPartialsInPotentials.end();

#ifdef VERBOSE
    std::cout << "revorder: node dual bound = " << improvedbound << " ; trivialbound = " << trivialbound << std::endl;
#endif
  }

  /////////////////////////////////////////////////////////////////////////////
  // Compute a dual bound by solving the LP relaxation of a chosen model
  /////////////////////////////////////////////////////////////////////////////

  // Start searching for better bounds only if trivial bound is not too far from primal bound and a significant number of nodes have been treated to improve primal bound
  if ( (this->userelaxbound_) && (this->treatednodes_ >= 1000) && (trivialbound * 2.0 >= this->primalbound_)) {
    // fix all the variables related to the vertices that already in the current incomplete order
    for (Vertex* v1 : this->inst_->vertices_) {
      if (node.rank(v1) >= this->inst_->L() + 2) {
        if (node.nbrefs(v1) >= this->inst_->U()) {
          isfullref_[v1].setBounds(1, 1);
        }
        else {
          isfullref_[v1].setBounds(0, 0);
        }
      }
      else if ((node.rank(v1) >= 1) && (node.rank(v1) <= this->inst_->L()+1)) {
        isfullref_[v1].setBounds(1, 1);
      }
      else {
        if (node.mustbefullref_[v1] == true) {
          isfullref_[v1].setBounds(1,1);
        }
        else {
          isfullref_[v1].setBounds(0, 1);
        }
      }
      for (Vertex* v2 : v1->neighbors_) {
        if ((node.rank(v1) < 0) && (node.rank(v2) < 0)) isbefore_[v1][v2].setBounds(0, 1);
        else if ( (node.rank(v1) < 0) && (node.rank(v2) >= 1)) isbefore_[v1][v2].setBounds(0, 0);
        else if ((node.rank(v1) >= 1) && (node.rank(v2) < 0)) isbefore_[v1][v2].setBounds(1, 1);
        else if ((node.rank(v1) >= 1) && (node.rank(v2) >= 1) && (node.rank(v1) < node.rank(v2)) ) isbefore_[v1][v2].setBounds(1, 1);
        else if ((node.rank(v1) >= 1) && (node.rank(v2) >= 1) && (node.rank(v1) > node.rank(v2)) ) isbefore_[v1][v2].setBounds(0, 0);
        else {
          cout << node.rank(v1) << "; " << node.rank(v2) << endl;
          throwError("abnormal rank values");
        }
      }
    }

    IloExpr sumnotallfullref(this->env_);
    int maxnbfullref;
    for (Vertex* v : node.potentialchildren_) {
      sumnotallfullref += isfullref_[v];
    }
    maxnbfullref = node.potentialchildren_.size() - 1;
    IloRange cons(this->env_, sumnotallfullref, maxnbfullref);
    relax_.add(cons);
    cplexrelax_.setParam(IloCplex::SimDisplay, 0);
//        cplexrelax_.setParam(IloCplex::Param::ParamDisplay, 0);
    cplexrelax_.setParam(IloCplex::Threads, 1);
    cplexrelax_.setParam(IloCplex::MIPDisplay, 0);
    cplexrelax_.setParam(IloCplex::TuningDisplay, 0);
    cplexrelax_.setOut(cplexrelax_.getEnv().getNullStream());
    bool updateprimal_ = false;
    if ( updateprimal_ &&  (node.nbordered() >= 0.7 * this->inst_->nbvertices_) ) {
      cplexrelax_.setParam(IloCplex::NodeLim, 9223372036800000000);
      cplexrelax_.setParam(IloCplex::MIPSearch, 1);
      if ((this->modeltype_ == CCG) || (this->modeltype_ == WITNESS)) {
        cplexrelax_.use(CyclesLazyConstraintsBB(cplexrelax_.getEnv(), *(this->inst_), *this));
      }
    }
    else {
      cplexrelax_.setParam(IloCplex::NodeLim, 0);
    }

    cplexrelax_.solve();
    IloAlgorithm::Status status = cplexrelax_.getStatus();
    //
    // if the relaxation is infeasible, the node can be pruned : set dual bound to negative value to state this
    if (status != IloAlgorithm::Infeasible) {
      int relaxbound = ceil(cplexrelax_.getBestObjValue());
      //
#ifdef VERBOSE
      std::cout << "revorder: bb: trivial dual bound = " << trivialbound << " ; relaxation bound = " << relaxbound;
            std::cout << " (best primal bound = " << this->primalbound_ << ")" << std::endl;
#endif
      if ( updateprimal_ &&  (node.nbordered() >= 0.7 * this->inst_->nbvertices_) ) {
        if (relaxbound < this->primalbound_) {
          this->primalbound_ = relaxbound;
          std::cout << "best primal bound updated = " << this->primalbound_ << std::endl;
          this->bestnbfullref_ = this->inst_->nbvertices_ - this->primalbound_;
          float tolerance = 1e-06;
          std::map<Vertex*, std::vector<Vertex*> > adjlist;
          for (Vertex* v : this->inst_->vertices_) {
            adjlist[v] = std::vector<Vertex*>();
          }

          for (Vertex* u : this->inst_->vertices_) {
            for (Vertex* v : u->neighbors_) {
              if (cplexrelax_.getValue(this->isbefore_[u][v]) >= 1 - tolerance) {
                adjlist[u].push_back(v);
              }
            }
          }
          //
          // search for the root of the digraph
          Vertex* root = nullptr;
          for (Vertex* v : this->inst_->vertices_) {
            if (node.rank(v) == 1) {
              root = v;
            }
          }
          // Search for a topological order that will be valid only if the digraph is
          // acyclic
          std::vector<Vertex*> reverseorder;
          std::map<Vertex*, int> rank;
          this->topologicalorder(root, adjlist, reverseorder, rank);
          // determine whether the digraph is cyclic or not
          //
          std::vector<std::pair<Vertex*, Vertex*> > reverseedges;
          bool iscyclic = this->getreverseedges(adjlist, reverseorder, rank, reverseedges);
          if (iscyclic) {
            std::cout << "revorder: error: the final integer solution is cyclic" << std::endl;
            throw;
          }
          for (Vertex* v : this->inst_->vertices_) {
            this->bestrank_[v] = rank[v];
          }
        }
      }
      //
      // update current dual bound if improved
      if (relaxbound > node.dualbound_) {
        node.dualbound_ = relaxbound;
        //
#ifdef VERBOSE
        std::cout << "revorder: bb: new dual bound = " << relaxbound << std::endl;
                std::cout << "revorder: bb: number of ordered vertices = " << node.nbordered() << std::endl;
#endif
        //
      }
    }
    else{
      node.dualbound_ =  this->inst_->nbvertices_;
#ifdef VERBOSE
      std::cout << "revorder: bb: dual bound relaxation is infeasible: prune the node" << std::endl;
#endif
    }

    relax_.remove(cons);
    cons.end();
  }
  return node.dualbound_;
}

// initialize the IP model for improved dual bound
void BBSolver::initializedualboundIP() {
  this->modeldual_ = IloModel(this->envdual_);
  this->cplexdual_ = IloCplex(this->modeldual_);

  IloExpr obj(this->envdual_);
  char name[256];
  for (Vertex* v : this->inst_->vertices_) {
    this->ispartial_[v] = IloBoolVar(this->envdual_);
    sprintf(name, "ispartial%i", v->id_);
    this->ispartial_[v].setName(name);
    obj += this->ispartial_[v];
  }
  this->modeldual_.add(IloMinimize(this->envdual_, obj));
  obj.end();
  //
  // if we could identify cliques that contain at least one partially-referenced vertex, add a constraint to enforce this; the two-cliques need not be checked
  for (Clique* c : this->cliqueswithpartialref_) {
    IloExpr sumpartialsinclique(this->envdual_);
    for (Vertex* v : c->vertices_) sumpartialsinclique += this->ispartial_[v];
    this->modeldual_.add(sumpartialsinclique >= this->nbpartialinclique_[c]);
    sumpartialsinclique.end();
  }
  //
  // if we could identify cycles composed of low degree vertices add a cut specifying that the vertices included in such cycles must include at least as many partially referenced vertices as there are cycles
  IloRangeArray ctaryLowDegreeCycles(this->envdual_);
  int nbcons = 0;
  for (int i = 0; i < this->lowdegreecycles_.size(); i++) {
    std::vector<Vertex*> cycle = this->lowdegreecycles_[i];
    IloExpr sumpartialsincycle(this->envdual_);
    for (Vertex* v : cycle) {
      sumpartialsincycle += this->ispartial_[v];
    }
    ctaryLowDegreeCycles.add(sumpartialsincycle >= this->nbpartialsincycle_[i]);
    sprintf(name, "ctaryLowDegreeCycles_%i", nbcons);
    ctaryLowDegreeCycles[nbcons++].setName(name);
    sumpartialsincycle.end();
  }
  this->modeldual_.add(ctaryLowDegreeCycles);
  ctaryLowDegreeCycles.end();
  //
  // load the integer problem and set display parameters
  this->cplexdual_.setParam(IloCplex::MIPDisplay, 0);
  this->cplexdual_.setParam(IloCplex::SimDisplay, 0);
  this->cplexdual_.setParam(IloCplex::TiLim, 10);
  this->cplexdual_.setParam(IloCplex::Threads, 1);
//    this->cplexdual_.setParam(IloCplex::Param::ParamDisplay, 0);
}