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"

#include <algorithm>


ILOSTLBEGIN
//

ILOLAZYCONSTRAINTCALLBACK2(CyclesLazyConstraintsBB,
                           Instance &, inst,
                           DiscretizationSolver&, solver) {
#ifdef VERBOSE
  std::cout << "revorder: check lazy dicycle constraints" << std::endl;
#endif
  // Search for cycles in the current solution digraph
  //
  // generate the adjacency lists of the vertices corresponding to the current
  // cplex solution
  float tolerance = 1e-06;
  std::map<Vertex*, std::vector<Vertex*> > adjlist;
  std::map<Vertex*, int > nbrefs;
  for (Vertex* v : inst.vertices_) {
    adjlist[v] = std::vector<Vertex*>();
    nbrefs[v] = 0;
  }
  for (Vertex* u : inst.vertices_) {
    if (solver.modeltype_ == WITNESS) {
      if (getValue(solver.isvertexinclique(u) >= 1 - tolerance)) continue;
    }
    for (Vertex* v : u->neighbors_) {
      if (getValue(solver.isbefore(u, v)) >= 1 - tolerance) {
        adjlist[u].push_back(v);
        nbrefs[v]++;
      }
    }
  }
  //
  // search for the root of the digraph
  Vertex* root = nullptr;
  int minnbrefs = inst.U();
  if (solver.modeltype_ == WITNESS) {
    for (Vertex* v : inst.vertices_) {
      if (getValue(solver.isvertexinclique(v) >= 1 - tolerance)) {
        root = v;
        break;
      }
    }
  } else {
    for (Vertex* v : inst.vertices_) {
      if (nbrefs[v] < minnbrefs) {
        minnbrefs = nbrefs[v];
        root = v;
        break;
      }
    }
  }
  //
  // call the enumeration of cycles, returns false if there are none
  std::vector<std::vector<Vertex*> > cycles;
  bool iscyclic = solver.enumeratecycles(adjlist, root, cycles);

  // Look for dicycle inequalities : we restrict the search to the dicycles
  // that contain at least one edge in the distance graph
  // besides that, the search is performed using brute force
  int nbaddedcuts = 0;
  IloEnv env = getEnv();
  if (iscyclic) {
    for (std::vector<Vertex*> path : cycles) {
      IloExpr sumedges(env);
      int cyclelength = path.size();
      Vertex* u = path.back();
      for (int i = 0; i < cyclelength; i++) {
        Vertex* v = path[i];
        sumedges += solver.isbefore(u, v);
        u = v;
      }
      if ((solver.modeltype_ == WITNESS) && (cyclelength <= inst.L() + 1)) {
        add(sumedges - solver.isvertexinclique(path.front())
                <= cyclelength - 1).end();
        nbaddedcuts++;
      } else {
        add(sumedges <= cyclelength - 1).end();
        nbaddedcuts++;
      }
    }
#ifdef VERBOSE
    std::cout << "revorder: number of violated dicycle constraints : " << nbaddedcuts << std::endl;
#endif
  } else {
#ifdef VERBOSE
    std::cout << "revorder: the graph is acyclic " << std::endl;
#endif
  }
  return;
}
/*************************************************************************************************************************************************************/
//
// sort the nodes by ascending number of ordered vertices
bool bbnodesnborderedcompare(BBNode* n1, BBNode* n2) {
  return (n1->nbordered() < n2->nbordered());
}
bool bbnodesbestcompare(BBNode*n1, BBNode* n2) {
  return (static_cast<float>(n1->cost_) / n1->nbordered()
      > static_cast<float>(n2->cost_) / n2->nbordered());
}

// Public methods
//
// constructor and destructor
//
// constructor used only for the root
BBNode::BBNode(Instance& inst):
    depth_(0),
    inst_(&inst),
    nbordered_(0),
    nbfullref_(0),
    istreated_(true) {}
//
// constructor used only for the initial cliques
BBNode::BBNode(Instance& inst,
               BBNode* root,
               Clique& c,
               bool branchonsmallerindex):
    depth_(1),
    father_(root),
    inst_(&inst),
    nbordered_(0),
    nbfullref_(0),
    istreated_(false) {
  this->orderedvertices_.clear();
  for (Vertex* v : inst.vertices_) {
    this->rank_[v] = -1;
    this->nbrefs_[v] = 0;
    this->isfullref_[v] = false;
    this->isordered_[v] = false;
    this->ispotentialchild_[v] = false;
    this->mustbefullref_[v] = false;
    this->maxrefs_[v] = v->degree_;
  }
  int currentrank = 1;
  for (Vertex* u : c.vertices_) {
    this->rank_[u] = currentrank++;
    this->orderedvertices_.push_back(u);
    this->switchtoordered(u);
    this->nbordered_++;
    //
    // update the number of references of the neighbors of the vertex
    for (Vertex* v : u->neighbors_) {
      if (!this->isordered_[v]) {
        this->nbrefs_[v]++;
      }
    }
  }

  // record the list of partially and fully referenced (unordered) vertices
  for (Vertex* v : inst.vertices_) {
    if (!this->isordered_[v]) {
      // record the list of unordered fully-referenced vertices
      if (this->nbrefs_[v] >= inst.U()) {
        this->isfullref_[v] = true;
        this->fullrefs_.push_back(v);
      } else if (this->nbrefs_[v] >= inst.L()) {
        // add partially referenced vertices to the list of potential choices
        // for branching
        if (branchonsmallerindex) {
          bool largerindex = true;
          Vertex* u = c.vertices_.back();
          if (!u->isneighbor(v)) {
            if (u->id_ >= v->id_) largerindex = false;
          }
          //
          // to break symmetries, make sure that the chosen potential child
          // is always that with biggest index
          if (largerindex) {
            this->ispotentialchild_[v] = true;
            this->potentialchildren_.push_back(v);
            //
            // partially referenced vertices with less than U neighbors can
            // be ordered right away since they will never be fully referenced
            if (this->maxrefs_[v] <= inst.U() - 1) {
              this->readytoorder_.push_back(v);
#ifdef VERBOSE
              std::cout << "revorder: bb: order vertex " << v->id_ << " without branching, maximum nb of references = " << this->maxrefs_[v] << std::endl;
#endif
            }
          } else {
            this->mustbefullref_[v] = true;
            // if a vertex must be fully-referenced, but it has exactly U
            // neighbors, it will for sure come after all its neighbors in
            // the order
            if (this->maxrefs_[v] <= inst.U()) {
              for (Vertex* u : v->neighbors_) {
                if (this->mustbefullref_[u]) {
                  if (this->maxrefs_[u] >= inst.U()+1) {
                    this->maxrefs_[u]--;
                  }
                } else {
                  this->maxrefs_[u]--;
                }
              }
            }
          }
        } else {
          this->ispotentialchild_[v] = true;
          this->potentialchildren_.push_back(v);

          // partially referenced vertices with less than U neighbors can be
          // ordered right away since they will never be fully referenced
          if (this->maxrefs_[v] <= inst.U() - 1) {
            this->readytoorder_.push_back(v);
#ifdef VERBOSE
            std::cout << "revorder: bb: order vertex " << v->id_ << " without branching, maximum nb of references = " << this->maxrefs_[v] << std::endl;
#endif
          }
        }
      }
    }
  }
  this->dualbound_ = -1;
  root->addtochildren(this);
}

// Make a copy of the input BB node: used for every other node
BBNode::BBNode(BBNode& node): id_(node.id_),
                              depth_(node.depth_+1),
                              father_(&node),
                              inst_(node.inst_),
                              nbordered_(node.nbordered_),
                              nbfullref_(node.nbfullref_),
                              istreated_(false) {
  cost_ = node.cost_;
  dualbound_ = node.dualbound_;
  primalbound_ = node.primalbound_;
  for (Vertex* v : inst_->vertices_) {
    rank_[v] = node.rank_.at(v);
    isordered_[v] = node.isordered_.at(v);
    nbrefs_[v] = node.nbrefs_.at(v);
    isfullref_[v] = node.isfullref_.at(v);
    ispotentialchild_[v] = node.ispotentialchild_.at(v);
    mustbefullref_[v] = node.mustbefullref_.at(v);
    maxrefs_[v] = node.maxrefs_.at(v);
  }
  fullrefs_.clear();
  orderedvertices_.clear();
  potentialchildren_.clear();
  children_.clear();
  for (Vertex* v : node.readytoorder_) readytoorder_.push_back(v);
  for (Vertex* v : node.fullrefs_) fullrefs_.push_back(v);
  for (Vertex* v : node.orderedvertices_) orderedvertices_.push_back(v);
  for (Vertex* v : node.potentialchildren_) potentialchildren_.push_back(v);
}
//
// Destroy the node
BBNode::~BBNode() {


}
//
// remove a vertex from the list of potential choice for branching
void BBNode::removefrompotentialchildren(Vertex* v) {
  //
  // no need to proceed if the vertex is not marked as potential child
  if (!this->ispotentialchild_[v]) return;

  this->ispotentialchild_[v] = false;
  std::vector<Vertex*>::iterator itv = this->potentialchildren_.begin();
  while (itv < this->potentialchildren_.end()) {
    if (*itv == v) {
      this->potentialchildren_.erase(itv);
      break;
    } else {
      itv++;

      // detect error if the vertex was not found in the list of potential
      // children
      if (itv == this->potentialchildren_.end()) {
        std::cout << "revorder: bb: did not find a potential child in the list"
                  << '\n';
        throw;
      }
    }
  }
  itv = this->readytoorder_.begin();
  while (itv < this->readytoorder_.end()) {
    if (*itv == v) {
      this->readytoorder_.erase(itv);
      break;
    }
    itv++;
  }
}
//
// assign the next available rank to the input vertex
void BBNode::assignnextrank(Vertex* u) {
  //
  // assign the next available rank
  for (Vertex* v : this->orderedvertices_) {
    if (u == v) {
      std::cout << "Les ennuis commencent ici." << std::endl;
    }
  }
  this->rank_[u] = this->nbordered_ + 1;
  this->orderedvertices_.push_back(u);
  this->switchtoordered(u);
  this->nbordered_++;
  this->removefrompotentialchildren(u);

  if (this->nbrefs_[u] <= this->inst_->L() - 1) {
    std::cout << "revorder: bb: trying to assign a rank to a vertex with less"
                 " than L references" << '\n';
    throw;
  } else if (this->nbrefs_[u] >= this->inst_->U()) {
    this->nbfullref_++;
  }

  // update the number of references of the neighbors of v
  for (Vertex* v : u->neighbors_) {
    if (this->rank_[v] < 0) {
      this->nbrefs_[v]++;

      // add partially-referenced vertices to the list of potential choices
      // for branching
      if (this->nbrefs_[v] == this->inst_->L()) {
        this->potentialchildren_.push_back(v);
        this->ispotentialchild_[v] = true;

        // partially referenced vertices with less than U neighbors can be
        // ordered right away since they will never be fully referenced
        if (this->maxrefs_[v] <= this->inst_->U() - 1) {
          this->readytoorder_.push_back(v);
#ifdef VERBOSE
          std::cout << "revorder: bb: order vertex " << v->id_ << " without branching, maximum nb of references = " << this->maxrefs_[v] << std::endl;
#endif
        }
      }
      // a vertex is fully-referenced if it has more than U references, and it
      // is not a potential choice for branching anymore
      if (this->nbrefs_[v] == this->inst_->U()) {
        this->isfullref_[v] = true;
        this->fullrefs_.push_back(v);
        if (this->ispotentialchild_[v]) this->removefrompotentialchildren(v);
      }
    }
  }
}

// get the current objective value at a node, depending on the ordered vertices
void BBNode::setnodecost(Problem pb) {
  int nbpartials = 0;
  for (Vertex* u : this->inst_->vertices_) {
    if (this->isordered_[u]) {
      if (this->nbrefs_[u] <= this->inst_->U() - 1) nbpartials += 1;
    }
  }
  this->cost_ =  nbpartials;
}

// a bb node n1 dominates another one, n2, if the corresponding partial
// solution has a smaller or equal cost but its list of ordered vertices
// contains that of n2
bool BBSolver::dominates(BBNode& n1, BBNode& n2) {
  if ( n2.cost_ < n1.cost_ )  {
    return false;
  }

  for (Vertex* v : n2.orderedvertices_) {
    if (!(n1.isordered(v))) {
      return false;
    }
  }
  return true;
}
//
// propagate the set of ordered vertices by iteratively ordering
// fully-referenced vertices
void BBNode::prepropagation() {
#ifdef VERBOSE
  cout << "prepropagation" << endl;
#endif
  // iteratively propagate while there are fully-referenced vertices
  while (!this->fullrefs_.empty()) {
    Vertex* u = this->fullrefs_.back();
    this->fullrefs_.pop_back();
    this->isfullref_[u] = false;
    this->assignnextrank(u);
  }
  this->setnodecost(MINPARTIAL);
}
//
// propagate the set of ordered vertices by iteratively ordering vertices
// as long as there is either at least one unordered vertex with more than U
// references or exactly one unordered vertex with exactly U references
void BBNode::postpropagation() {
#ifdef VERBOSE
  cout << "in postpropagation" << endl;
#endif
  //
  // no need of propagation if all the vertices are ordered
  if (this->nbordered_ == this->inst_->nbvertices_) return;
  //
  // iteratively propagate while there are fully-referenced vertices or
  // exactly one partially-referenced vertex
  while (true) {
    bool ispropagate = false;
    Vertex* u;
    if (!this->fullrefs_.empty()) {
      // add a vertex with more than U references
      u = this->fullrefs_.back();
      this->fullrefs_.pop_back();
      this->isfullref_[u] = false;
      ispropagate = true;
    } else if (!this->readytoorder_.empty()) {
      // if there is a partially referenced vertex with less than U
      // neighbors, it can be ordered right away: it will never be fully
      // referenced
      u = this->readytoorder_.back();
      this->readytoorder_.pop_back();
      ispropagate = true;
    } else if (this->potentialchildren_.size() == 1) {
      // if there is exactly one partially referenced vertex: add it to the
      // order
      u = this->potentialchildren_.back();
      this->potentialchildren_.pop_back();
      this->ispotentialchild_[u] = false;
      ispropagate = true;
    }

    if (!ispropagate) break;
    //
    // add the identified vertex at the end of the order
    this->assignnextrank(u);
  }
  //
  // update the cost of the node
  this->setnodecost(MINPARTIAL);
}

/*******************************************************************************
 Class BBSolver
 - an instance of BBSolver is a solver of the discretization problem using
 branch-and-bound
 *******************************************************************************/
// Public methods
//
// constructor and destructor
BBSolver::BBSolver(Instance* inst,
                   Problem pb,
                   std::string  optionfile,
                   float timelimit):
    DiscretizationSolver(inst,
                         BRANCHANDBOUND,
                         pb,
                         NONE,
                         optionfile,
                         timelimit) {
  this->isfeasible_ = false;
  this->isoptimal_ = false;
  this->bbnodes_ = 1;
  this->treatednodes_ = 0;
  this->nbactivenodes_ = 0;
  this->primalbound_ = inst->nbvertices_;

  if (!optionfile.empty())
  {
    std::ifstream infile(optionfile.c_str());
    if (!infile.is_open()) {
      throwError("revorder: error: the option file could not be opened");
    }
    // read the file line by line
    // one line contains the data relative to one option
    std::string line;
    while (std::getline(infile, line))
    {
      std::istringstream iss(line);
      std::string optionname;
      char buf[100];
      int optionvalue;
      std::string optionstring;

      //read the line and detect format errors
      if (!(iss >> buf)){
        std::cerr << "revorder: error: there is a mistake with the format of "
                     "the option file" <<std::endl;
        throw;
      } // error
      optionname  = std::string(buf);
      if (optionname.find("relaxmodel") != std::string::npos) {
        if (!(iss >> buf)){
          std::cerr << "revorder: error: there is a mistake with the format "
                       "of the option file" <<std::endl;
          throw;
        }
        optionstring = std::string(buf);
        if (!strcmp(optionstring.c_str(),"ccg"))
        {
          this->modeltype_ = CCG;
        }
        else if (!strcmp(optionstring.c_str(),"cycles"))
        {
          this->modeltype_ = CYCLES;
        }
        else if (!strcmp(optionstring.c_str(),"ranks"))
        {
          this->modeltype_ = RANKS;
        }
        else if (!strcmp(optionstring.c_str(),"vertexrank"))
        {
          this->modeltype_ = VERTEXRANK;
        }
        else if (!strcmp(optionstring.c_str(),"witness"))
        {
          this->modeltype_ = WITNESS;
        }
        continue;
      }
      if (!(iss >> optionvalue)){
        std::cerr << "revorder: error: there is a mistake with the format of "
                     "the option file" <<std::endl;
        throw;
      } // error
      //
      if (optionname.find("branchonsmallerindex") !=std::string::npos) {
        branchonsmallerindex_ = (bool) optionvalue;
        continue;
      }
      else if (optionname.find("prunesameneighbors") !=std::string::npos) {
        prunesameneighbors_ = (bool) optionvalue;
        continue;
      }
      else if (optionname.find("checkdominance") !=std::string::npos) {
        if (optionname.find("checkdominancealltree") !=std::string::npos) {
          checkdominancealltree_ = (bool) optionvalue;
          continue;
        }
        checkdominance_ = (bool) optionvalue;
        continue;
      }
      else if (optionname.find("maxdeltadominance") !=std::string::npos) {
        maxdeltadominance_ = optionvalue;
        continue;
      }
      else if (optionname.find("improvedualbound") !=std::string::npos) {
        improvedualbound_ = (bool) optionvalue;
        continue;
      }
      else if (optionname.find("userelaxbound") !=std::string::npos) {
        userelaxbound_ = (bool) optionvalue;
        continue;
      }
      else if (optionname.find("explorebest") !=std::string::npos) {
        explorebest_ = (bool) optionvalue;
        continue;
      }
      else if (optionname.find("exploredepth") !=std::string::npos) {
        if (optionname.find("exploredepthbeforebest") !=std::string::npos) {
          exploredepthbeforebest_ = (bool) optionvalue;
          continue;
        }
        exploredepth_ = (bool) optionvalue;
        continue;
      }
    }
  }
  else
  {
    //
    // set some parameters for the solution of the relaxation in dual bounds
    // computations
    this->modeltype_ = CCG;
    this->branchonsmallerindex_ = true;
    this->prunesameneighbors_ = true;
    this->checkdominance_ = true;
    this->checkdominancealltree_ = true;
    this->maxdeltadominance_ = 1;
    this->improvedualbound_ = false;
    this->userelaxbound_ = false;
    this->explorebest_ = false;
    this->exploredepth_ = false;
    this->exploredepthbeforebest_ = true;
  }
  //
  // build a root node
  this->rootnode_ = new BBNode(*inst);
}
//
BBSolver::~BBSolver() {
  cplexrelax_.end();
  relax_.end();
  env_.end();
  cplexdual_.end();
  modeldual_.end();
  envdual_.end();
  delete rootnode_;
}
//
// main method : called to run the branch-and-bound algorithm
int BBSolver::solve() {
  //
  // run timer
  IloTimer cpuClockTotal(this->env_);
  IloTimer cpuClockInit(this->env_);
  cpuClockTotal.start();
  cpuClockInit.start();

  std::cout << std::endl;
  std::cout << "revorder: "
               "----------------------------------------------------------"
            << std::endl;
  std::cout << "\nrevorder: bb: Run preprocessing procedures "
            <<  std::endl << std::endl;
  //
  // try and reduce the size of the instance right from the beginnning
  this->inst_->computeneighbors();
  this->preallocatelowdegreevertices();
  //
  // enumerate the potential initial cliques
  if (!this->enumeratecliques(this->inst_->L() + 1)) {
    this->isfeasible_ = false;
    return false;
  }
  //
  // eliminate the redundant and dominated cliques
  this->eliminateredundantcliques();
  //
  // solve with greedy and remove cliques that cannot be completed
  this->isfeasible_ = this->greedysolve();
  if (!this->isfeasible_) return false;
  //
  // initialize the clique nodes
  for (Clique* c : this->cliques_) {
    BBNode* cliquenode =new BBNode(*inst_, rootnode_, *c,
                                   this->branchonsmallerindex_);
    this->bbnodesqueue_.push_back(cliquenode);
    cliquenode->primalbound_ = c->greedyobjvalue_;
    cliquenode->setid(this->bbnodes_);
    cliquenode->setnodecost(this->problem_);
    cliquenode->prepropagation();
    cliquenode->postpropagation();
    this->activenodes_[cliquenode->cost_].push_back(cliquenode);
    this->nbactivenodes_++;
    //
    // first propagate the partial order
    this->bbnodes_++;
  }
  std::cout << "revorder: bb: preprocessing created "
            <<  this->bbnodesqueue_.size()
            << " initial clique nodes" << std::endl;
  //
  // sort the bb nodes according to the choice of exploration method
  if (this->exploredepth_ || this->exploredepthbeforebest_) {
    std::make_heap(this->bbnodesqueue_.begin(), this->bbnodesqueue_.end(),
                   bbnodesnborderedcompare);
  } else if (this->explorebest_) {
    std::make_heap(this->bbnodesqueue_.begin(), this->bbnodesqueue_.end(),
                   bbnodesbestcompare);
  } else {
    std::cout << "revorder: bb: need to choose an exploration method in "
                 "branch-and-bound" << std::endl;
    throw;
  }
  //
#ifdef VERBOSE
  std::cout << "revorder: bb: greatest number of ordered vertices in a clique: ";
    std::cout << ", " << this->bbnodesqueue_.front()->nbordered() << " vertices" << std::endl;
#endif
  //
  // compute sets of vertices that must contain at least one
  // partially-referenced vertex
  this->computecliquecuts(this->cliquecutsmaxsize_);
  //
  // identify cycle of vertices with degree smaller than U+1: in each of
  // these cycles, at least one vertex is partially referenced
  this->enumeratelowdegreecycles();
  //
  // initialize the IP used for improved dual bound
  this->initializedualboundIP();
  //
  // set greedy solution as initial primal bound
  this->primalbound_ = this->objvalue_;
  //
  cpuClockInit.stop();

  // Run the branch-and-bound algorithm
  //
  // initialize the ip model to find dual bounds using the linear relaxation
//    IloModel model(this->env_);
  this->relax_ = IloModel(this->env_);
//    this->defineminpartialip(model);
  if (userelaxbound_) {
//        this->createrelaxationmodel(model, this->relax_);
    this->defineminpartialip(this->relax_);
    this->cplexrelax_ = IloCplex(this->env_);
    this->cplexrelax_.extract(this->relax_);
  }
  //
  // treat

  std::cout << std::endl;
  std::cout << "revorder: "
               "----------------------------------------------------------"
            << std::endl;
  std::cout << "\nrevorder: bb: Start the B&B algorithm " <<  std::endl
            << std::endl;

  while (!this->bbnodesqueue_.empty()) {
    this->treatnode(*(this->bbnodesqueue_.front()));
    if (cpuClockTotal.getTime() > this->timelimit_) {
      std::cout << "revorder: bb: time limit has been reached, stop the "
                   "solution algorithm" << std::endl;
      break;
    }
  }
  cpuClockTotal.stop();

  // Display the results
  this->totaltime_ = cpuClockTotal.getTime();
  this->relaxtime_ = 0.0;

  std::cout << std::endl;
  std::cout << "revorder: "
               "----------------------------------------------------------"
            << std::endl;
  std::cout << "revorder: "
               "----------------------------------------------------------"
            << std::endl;
  std::cout << "\nmdjeep: Branch-and-bound solution report: " << std::endl;

  if (isfeasible_)	{
    this->objvalue_ = this->primalbound_;

    // the solution is optimal only if all the bb nodes have been treated
    if (this->bbnodesqueue_.empty()) {
      isoptimal_ = true;
      std::cout << "revorder: the instance was solved to optimality"
                << std::endl;
    } else {
      std::cout << "revorder: the branch-and-bound was stopped before proving"
                   " optimality" << std::endl;
    }

    // reconstruct the solution by including the preallocated vertices
    objvalue_ += this->inst_->preallocatedvertices_.size();
    this->reconstructsolution();
    //
    // display the solution
    std::cout << "revorder: total cpu time                  = " << cpuClockTotal.getTime() << " s" << std::endl;
    std::cout << "revorder: initialization cpu time         = " << cpuClockInit.getTime() << " s" << std::endl;
    std::cout << "revorder: number of explored nodes        = " << this->treatednodes_ << std::endl;
    std::cout << "revorder: number of created nodes        = " << this->bbnodes_ << std::endl;
    std::cout << "revorder: value of the objective function = " << this->objvalue_ << " (after addition of preallocated vertices)" << std::endl;
    //
    // verify the resulting order is indeed a revorder
    this->verifyorder(this->bestrank_);
  }
  else {
    std::cout << "revorder: the instance is not discretizable!" << std::endl;
  }
  return isfeasible_;
}
//
// treat a branch-and-bound node (compute bounds and branch)
void BBSolver::treatnode(BBNode& node) {
  this->treatednodes_++;
  if (std::remainder(this->treatednodes_, 100) == 0) {
    std::cout << "revorder: bb: treated nodes = " << this->treatednodes_
              << ", nodes in the queue = " << this->bbnodesqueue_.size()
              << ", nodes in memory = " << this->nbactivenodes_
              << ", primal bound = " << this->primalbound_ <<  std::endl;
  }
  //
  // remove the bb node from the list
  //
  // make sure to maintain the heap first
  if (this->isfeasible_) {
    if (this->exploredepth_) {
      std::pop_heap(this->bbnodesqueue_.begin(),
                    this->bbnodesqueue_.end(), bbnodesnborderedcompare);
    } else if (this->explorebest_ || this->exploredepthbeforebest_) {
      std::pop_heap(this->bbnodesqueue_.begin(),
                    this->bbnodesqueue_.end(), bbnodesbestcompare);
    }
  } else {
    if (this->explorebest_) {
      std::pop_heap(this->bbnodesqueue_.begin(),
                    this->bbnodesqueue_.end(), bbnodesbestcompare);

    } else {
      std::pop_heap(this->bbnodesqueue_.begin(),
                    this->bbnodesqueue_.end(), bbnodesnborderedcompare);
    }
  }

  // actually remove the node from the queue
  this->bbnodesqueue_.pop_back();

  // update the primal bound and prune the node if it is a leaf of the
  // enumeration tree
  if (node.nbordered() == this->inst_->nbvertices_) {
#ifdef VERBOSE
    std::cout << "revorder: bb: reached a leaf at node " << node.id() << ": depth = " << node.depth();
#endif
    node.primalbound_ = this->getobjvalue(node.rank_, node.nbrefs_);
    //
    //  update the primal bound if improved
    if (node.primalbound_ < this->primalbound_) {
      //
      // swap the search method if this is the first feasible solution
      std::cout << "revorder: bb: the primal bound has been improved, value = "
                << node.primalbound_ << std::endl;
      if (!this->isfeasible_) {
        if (this->exploredepthbeforebest_) {
          std::make_heap(this->bbnodesqueue_.begin(),
                         this->bbnodesqueue_.end(), bbnodesbestcompare);
        }
        this->isfeasible_ = true;
      }
      this->primalbound_ = node.primalbound_;
      this->bestnbfullref_ = this->getnbfullref(node.rank_, node.nbrefs_);
      for (Vertex* v : inst_->vertices_) {
        this->bestrank_[v] = node.rank(v);
      }
    }
    this->erasenodefromall(&node);
  } else {
    // otherwise, compute dual and primal bound and prepare for next node

    // prune the node if no potential child node
    if (node.potentialchildren_.empty()) {
      //
#ifdef VERBOSE
      std::cout << "revorder: bb: prune node " << node.id() << ": ";
            std::cout << ": no feasible solution on this branch" << std::endl;
#endif
      //
      this->erasenodefromall(&node);
    } else {
      // compute a dual bound based on the ordered vertices
      this->computedualbound(node);
      //
      // prune the node by using bounds
      if (node.dualbound_ >= this->primalbound_) {
#ifdef VERBOSE
        std::cout << "revorder: bb: prune node " << node.id() << ": ";
                std::cout << "using bounds: primal = " << this->primalbound_;
                std::cout << " ; dual = " << node.dualbound_ << std::endl;
#endif
        this->erasenodefromall(&node);
      } else {
        // otherwise, create new nodes by branching
        this->branch(node);
        node.istreated_ = true;
//                if (!this->checkdominancealltree_) {
//                    this->erasenodefromall(&node);
//                }
      }
    }
  }
  //
  // treat next node in the queue
  // since it is a heap, the front node is always the greatest in the
  // implemented order
#ifdef VERBOSE
  if (!this->bbnodesqueue_.empty()) {
        BBNode* nextnode = this->bbnodesqueue_.front();
        std::cout << '\n' << "revorder: bb: treat node " << nextnode->id() << " ; depth = " << nextnode->depth() << " ; ";
        std::cout << nextnode->nbordered() << " ordered vertices ; " << nextnode->nbpartialref() << " partially referenced vertices" << std::endl;
    }
#endif
}

// create new nodes by branching
void BBSolver::branch(BBNode& node) {
  // create the children of the input node
  // propagate the order starting from each incomplete order to detect
  // dominatedorders
  std::vector<BBNode*> childnodes;
  std::map<BBNode*, bool> isdominated;
  std::map<Vertex*, bool> issymmetric;
  //
  // we start by checking diverse pruning rules that need to be applied
  // before dominance
  for (int i = 0; i < node.potentialchildren_.size() ; i++) {
    Vertex* u = node.potentialchildren_[i];
    issymmetric[u] = false;

    if (!this->prunesameneighbors_) continue;

    for (int j = 0; j < i; j++) {
      Vertex* v = node.potentialchildren_[j];
      if (issymmetric[v]) continue;

      // if two vertices that can be chosen for branching have the exact same
      // unordered neighbors, only that with smaller number of references
      // needs to be considered for branching at this stage
      bool sameneighbors = true;
      for (Vertex* neighbor : u->neighbors_) {
        if (!node.isordered(neighbor)) {
          if ( !v->isneighbor_[neighbor] ) {
            sameneighbors = false;
            break;
          }
        }
      }
      if (sameneighbors) {
#ifdef VERBOSE
        std::cout << "revorder: bb: prune node: two potential children have same neighbors: vertices " << u->id_ << " and " << v->id_ <<  std::endl;
#endif
        if (node.nbrefs_[u] < node.nbrefs_[v]) {
          issymmetric[v] = true;
        } else if (node.nbrefs_[v] < node.nbrefs_[u]) {
          issymmetric[u] = true;
        } else {
          if (u->id_ < v->id_) issymmetric[v] = true;
          else if (v->id_ < u->id_) issymmetric[u] =true;
        }
      }
    }
  }

  for (Vertex* v : node.potentialchildren_) {
    if (issymmetric[v]) continue;  // do not consider the vertices that are
    // symmetric

    BBNode* child = new BBNode(node);
    //
    // in MINPARTIAL, we can always make the arbitrary decision that if a
    // vertex must be ordered when it is still partially referenced, then it
    // is the vertex with smallest index among all the potential children; as
    // a consequence the vertices with smaller index must be removed from the
    // list of potential children
    if ( this->branchonsmallerindex_ ) {
      std::vector<Vertex*>::iterator itv = child->potentialchildren_.begin();
      while (itv < child->potentialchildren_.end()) {
        if ( ((*itv)->id_ < v->id_) ) {
          child->ispotentialchild_[*itv] = false;
          child->potentialchildren_.erase(itv);
          child->mustbefullref_[*itv] = true;
          //
          // if a vertex must be fully-referenced, but it has exactly U
          // neighbors, it will for sure come after all its neighbors in the
          // order
          if (child->maxrefs_[*itv] <= this->inst_->U()) {
            for (Vertex* u : v->neighbors_) {
              if (child->mustbefullref_[u]) {
                if (child->maxrefs_[u] >= this->inst_->U()+1) {
                  child->maxrefs_[u]--;
                }
              } else {
                child->maxrefs_[u]--;
              }
            }
          }
        } else {
          itv++;
        }
      }
    }

    child->assignnextrank(v);
    child->prepropagation();
    isdominated[child] = false;
    childnodes.push_back(child);
  }
  //
  // check for classical dominance
  for (BBNode* n1 : childnodes) {
    if (!this->checkdominance_) continue;  // do not check dominance if option
    // is deactivated

    if (isdominated[n1]) continue;
    for (BBNode* n2 : childnodes) {
      if (n2 == n1) {
        continue;
      } else if (isdominated[n2]) {
        continue;
      } else if (this->dominates(*n2, *n1)) {
        // classical domination
        isdominated[n1] = true;
        //
#ifdef VERBOSE
        std::cout << "revorder: bb: prune node: dominance among children" <<  std::endl;
#endif
        break;
      } else if (this->dominates(*n1, *n2)) {
        isdominated[n2] = true;
        //
#ifdef VERBOSE
        std::cout << "revorder: bb: prune node: dominance among children" <<  std::endl;
#endif
      }
    }