Changeset 988


Ignore:
Timestamp:
Jun 24, 2008 6:17:58 PM (11 years ago)
Author:
ladanyi
Message:

still not there...

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/dynamicbranching/dynamicbranching.cpp

    r987 r988  
    4343#include "OsiClpSolverInterface.hpp"
    4444
    45 #define DEBUG_DYNAMIC_BRANCHING 100
     45#define DEBUG_DYNAMIC_BRANCHING
     46
     47#ifdef DEBUG_DYNAMIC_BRANCHING
     48int dyn_debug = 1;
     49#endif
    4650
    4751// below needed for pathetic branch and bound code
     
    4953#include <map>
    5054
     55using namespace std;
     56
    5157class DBVectorNode;
     58
     59class LPresult {
     60private:
     61  LPresult(const LPresult& rhs);
     62  LPresult& operator=(const LPresult& rhs);
     63  void gutsOfConstructor(const OsiSolverInterface& model);
     64public:
     65  LPresult(const OsiSolverInterface& model);
     66  ~LPresult();
     67public:
     68  bool isAbandoned;
     69  bool isProvenDualInfeasible;
     70  bool isPrimalObjectiveLimitReached;
     71  bool isProvenOptimal;
     72  bool isDualObjectiveLimitReached;
     73  bool isIterationLimitReached;
     74  bool isProvenPrimalInfeasible;
     75  double getObjSense;
     76  double* getReducedCost;
     77  double* getColLower;
     78  double* getColUpper;
     79  double* getObjCoefficients;
     80  double yb_plus_rl_minus_su;
     81};
     82
     83void
     84LPresult::gutsOfConstructor(const OsiSolverInterface& model)
     85{
     86  isAbandoned = model.isAbandoned();
     87  isProvenDualInfeasible = model.isProvenDualInfeasible();
     88  isPrimalObjectiveLimitReached = model.isPrimalObjectiveLimitReached();
     89  isProvenOptimal = model.isProvenOptimal();
     90  isDualObjectiveLimitReached = model.isDualObjectiveLimitReached();
     91  isIterationLimitReached = model.isIterationLimitReached();
     92  isProvenPrimalInfeasible = model.isProvenPrimalInfeasible();
     93  getObjSense = model.getObjSense();
     94
     95  getReducedCost = new double[model.getNumCols()];
     96  CoinDisjointCopyN(model.getReducedCost(), model.getNumCols(), getReducedCost);
     97  getColLower = new double[model.getNumCols()];
     98  CoinDisjointCopyN(model.getColLower(), model.getNumCols(), getColLower);
     99  getColUpper = new double[model.getNumCols()];
     100  CoinDisjointCopyN(model.getColUpper(), model.getNumCols(), getColUpper);
     101
     102  getObjCoefficients = NULL;
     103  yb_plus_rl_minus_su = 0;
     104 
     105  if (!isProvenOptimal &&
     106      !isDualObjectiveLimitReached &&
     107      !isIterationLimitReached &&
     108      !isAbandoned &&
     109      !isProvenDualInfeasible &&
     110      !isPrimalObjectiveLimitReached) {
     111    assert(isProvenPrimalInfeasible);
     112    getObjCoefficients = new double[model.getNumCols()];
     113    CoinDisjointCopyN(model.getObjCoefficients(), model.getNumCols(), getObjCoefficients);
     114    const std::vector<double*> dual_rays = model.getDualRays(1);
     115    if (dual_rays.size() == 0) {
     116      printf("WARNING: LP is infeas, but no dual ray is returned!\n");
     117      return;
     118    }
     119    const double* dual_ray = dual_rays[0];
     120    const double direction = model.getObjSense();
     121    const double* rub = model.getRowUpper();
     122    const double* rlb = model.getRowLower();
     123    const double* cub = model.getColUpper();
     124    const double* clb = model.getColLower();
     125    const double* dj  = model.getReducedCost();
     126    const double* obj = model.getObjCoefficients();
     127    const int numRows = model.getNumRows();
     128    const int numCols = model.getNumCols();
     129    for (int i = 0; i < numRows; ++i) {
     130      const double ray_i = dual_ray[i];
     131      if (ray_i > 1e-6) {
     132        yb_plus_rl_minus_su += ray_i*rlb[i];
     133      } else if (ray_i < 1e-6) {
     134        yb_plus_rl_minus_su += ray_i*rub[i];
     135      }
     136    }
     137    for (int j = 0; j < numCols; ++j) {
     138      const double yA_j = dj[j] - obj[j];
     139      if (direction * yA_j > 1e-6) {
     140        yb_plus_rl_minus_su -= yA_j*cub[j];
     141      } else if (direction * yA_j < -1e-6) {
     142        yb_plus_rl_minus_su -= yA_j*clb[j];
     143      }
     144    }
     145    for (int i = dual_rays.size()-1; i >= 0; --i) {
     146      delete[] dual_rays[i];
     147    }
     148  }
     149}
     150 
     151LPresult::LPresult(const OsiSolverInterface& model)
     152{
     153  gutsOfConstructor(model);
     154}
     155
     156LPresult::~LPresult()
     157{
     158  delete[] getReducedCost;
     159  delete[] getColLower;
     160  delete[] getColUpper;
     161  delete[] getObjCoefficients;
     162}
    52163
    53164// Trivial class for Branch and Bound
     
    70181   
    71182  };
     183
     184private:
     185  void gutsOfCopy(const DBNodeSimple& rhs);
     186  void gutsOfConstructor (OsiSolverInterface &model,
     187                          int numberIntegers, int * integer,
     188                          CoinWarmStart * basis);
     189  void gutsOfDestructor(); 
    72190 
    73191public:
     
    81199                int numberIntegers, int * integer,
    82200                CoinWarmStart * basis);
    83   void gutsOfConstructor (OsiSolverInterface &model,
    84                           int numberIntegers, int * integer,
    85                           CoinWarmStart * basis);
    86201  // Copy constructor
    87202  DBNodeSimple ( const DBNodeSimple &);
     
    92207  // Destructor
    93208  ~DBNodeSimple ();
    94   // Work of destructor
    95   void gutsOfDestructor();
     209
    96210  // Extension - true if other extension of this
    97211  bool extension(const DBNodeSimple & other,
     
    100214  // Tests if we can switch this node (this is the parent) with its parent
    101215  bool canSwitchParentWithGrandparent(const int* which,
    102                                       OsiSolverInterface & model,
     216                                      const LPresult& lpres,
    103217                                      const int * original_lower,
    104218                                      const int * original_upper,
     
    458572}
    459573
    460 DBNodeSimple::DBNodeSimple(const DBNodeSimple & rhs)
     574void
     575DBNodeSimple::gutsOfCopy(const DBNodeSimple & rhs)
    461576{
    462577  node_id_=rhs.node_id_;
    463   if (rhs.basis_)
    464     basis_=rhs.basis_->clone();
    465   else
    466     basis_ = NULL;
     578  basis_ = rhs.basis_ ? rhs.basis_->clone() : NULL;
    467579  objectiveValue_=rhs.objectiveValue_;
    468580  variable_=rhs.variable_;
     
    488600}
    489601
     602DBNodeSimple::DBNodeSimple(const DBNodeSimple & rhs)
     603{
     604  gutsOfCopy(rhs);
     605}
     606
    490607DBNodeSimple &
    491608DBNodeSimple::operator=(const DBNodeSimple & rhs)
    492609{
    493610  if (this != &rhs) {
     611    // LL: in original code basis/lower/upper was left alone if rhs did not
     612    // have them. Was that intentional?
    494613    gutsOfDestructor();
    495     node_id_=rhs.node_id_;
    496     if (rhs.basis_)
    497       basis_=rhs.basis_->clone();
    498     objectiveValue_=rhs.objectiveValue_;
    499     variable_=rhs.variable_;
    500     way_=rhs.way_;
    501     numberIntegers_=rhs.numberIntegers_;
    502     value_=rhs.value_;
    503     parent_ = rhs.parent_;
    504     child_down_ = rhs.child_down_;
    505     child_up_ = rhs.child_up_;
    506     strong_branching_fixed_vars_ = rhs.strong_branching_fixed_vars_;
    507     reduced_cost_fixed_vars_ = rhs.reduced_cost_fixed_vars_;
    508     previous_ = rhs.previous_;
    509     next_ = rhs.next_;
    510     if (rhs.lower_!=NULL) {
    511       lower_ = new int [numberIntegers_];
    512       upper_ = new int [numberIntegers_];
    513       assert (upper_!=NULL);
    514       memcpy(lower_,rhs.lower_,numberIntegers_*sizeof(int));
    515       memcpy(upper_,rhs.upper_,numberIntegers_*sizeof(int));
    516     }
     614    gutsOfCopy(rhs);
    517615  }
    518616  return *this;
     
    580678  // Push
    581679  void push_back(DBNodeSimple & node); // the node_id_ of node will change
     680  /* Remove a single node from the tree and adjust the previos_/next_/first_
     681     etc fields. Does NOT update child/parent relationship */
     682  void removeNode(int n);
     683  /* Remove the subtree rooted at node n. properly adjusts previos_/next_ etc
     684     fields. Does NOT update child/parent relationships. */
     685  void removeSubTree(int n);
    582686  // Last one in (or other criterion)
    583687  DBNodeSimple back() const;
     
    589693  void moveNodeUp(const int* which,
    590694                  OsiSolverInterface& model, DBNodeSimple & node);
    591   // It changes the bounds of the descendants of node with node_id
    592   void changeDescendantBounds(const int node_id, const bool lower_bd,
    593                               const int brvar, const int new_bd);
     695  // Fix the bounds in the descendants of subroot
     696  void adjustBounds(int subroot, int brvar, int brlb, int brub);
     697
     698  // Check that the bounds correspond to the branching decisions...
     699  void checkTree() const;
     700  void checkNode(int node) const;
    594701
    595702  // Public data
     
    738845  return nodes_[best()];
    739846}
     847
     848void
     849DBVectorNode::removeNode(int n)
     850{
     851  DBNodeSimple& node = nodes_[n];
     852  if (node.bothChildDone())
     853    sizeDeferred_--;
     854  int previous = node.previous_;
     855  int next = node.next_;
     856  node.~DBNodeSimple();
     857  if (previous >= 0) {
     858    nodes_[previous].next_=next;
     859  } else {
     860    first_ = next;
     861  }
     862  if (next >= 0) {
     863    nodes_[next].previous_ = previous;
     864  } else {
     865    last_ = previous;
     866  }
     867  node.previous_ = -1;
     868  if (firstSpare_ >= 0) {
     869    node.next_ = firstSpare_;
     870  } else {
     871    node.next_ = -1;
     872  }
     873  firstSpare_ = n;
     874  assert (size_>0);
     875  size_--;
     876}
     877
     878void
     879DBVectorNode::removeSubTree(int n)
     880{
     881  if (nodes_[n].child_down_ >= 0) {
     882    removeSubTree(nodes_[n].child_down_);
     883  }
     884  if (nodes_[n].child_up_ >= 0) {
     885    removeSubTree(nodes_[n].child_up_);
     886  }
     887  removeNode(n);
     888}
     889
    740890// Get rid of last one
    741891void
     
    744894  // Temporary until more sophisticated
    745895  //assert (last_==chosen_);
    746   if (nodes_[chosen_].bothChildDone())
    747     sizeDeferred_--;
    748   int previous = nodes_[chosen_].previous_;
    749   int next = nodes_[chosen_].next_;
    750   nodes_[chosen_].gutsOfDestructor();
    751   if (previous>=0) {
    752     nodes_[previous].next_=next;
    753   } else {
    754     first_ = next;
    755   }
    756   if (next>=0) {
    757     nodes_[next].previous_ = previous;
    758   } else {
    759     last_ = previous;
    760   }
    761   nodes_[chosen_].previous_=-1;
    762   if (firstSpare_>=0) {
    763     nodes_[chosen_].next_ = firstSpare_;
    764   } else {
    765     nodes_[chosen_].next_ = -1;
    766   }
    767   firstSpare_ = chosen_;
     896  removeNode(chosen_);
    768897  chosen_ = -1;
    769   assert (size_>0);
    770   size_--;
    771 }
    772 
    773 static double
    774 compute_val_for_ray(const OsiSolverInterface& model)
    775 {
    776   const std::vector<double*> dual_rays = model.getDualRays(1);
    777   if (dual_rays.size() == 0) {
    778     printf("WARNING: LP is infeas, but no dual ray is returned!\n");
    779     return 0;
    780   }
    781   const double* dual_ray = dual_rays[0];
    782   const double direction = model.getObjSense();
    783   const double* rub = model.getRowUpper();
    784   const double* rlb = model.getRowLower();
    785   const double* cub = model.getColUpper();
    786   const double* clb = model.getColLower();
    787   const double* dj  = model.getReducedCost();
    788   const double* obj = model.getObjCoefficients();
    789   const int numRows = model.getNumRows();
    790   const int numCols = model.getNumCols();
    791   double yb_plus_rl_minus_su = 0;
    792   for (int i = 0; i < numRows; ++i) {
    793     const double ray_i = dual_ray[i];
    794     if (ray_i > 1e-6) {
    795       yb_plus_rl_minus_su += ray_i*rlb[i];
    796     } else if (ray_i < 1e-6) {
    797       yb_plus_rl_minus_su += ray_i*rub[i];
    798     }
    799   }
    800   for (int j = 0; j < numCols; ++j) {
    801     const double yA_j = dj[j] - obj[j];
    802     if (direction * yA_j > 1e-6) {
    803       yb_plus_rl_minus_su -= yA_j*cub[j];
    804     } else if (direction * yA_j < -1e-6) {
    805       yb_plus_rl_minus_su -= yA_j*clb[j];
    806     }
    807   }
    808   for (int i = dual_rays.size()-1; i >= 0; --i) {
    809     delete[] dual_rays[i];
    810   }
    811   return yb_plus_rl_minus_su;
    812898}
    813899
    814900bool
    815901DBNodeSimple::canSwitchParentWithGrandparent(const int* which,
    816                                              OsiSolverInterface & model,
     902                                             const LPresult& lpres,
    817903                                             const int * original_lower,
    818904                                             const int * original_upper,
     
    821907  /*
    822908    The current node ('this') is really the parent (P) and the solution in the
    823     model represents the child. The grandparent (GP) is this.parent_. Let's have
     909    lpres represents the child. The grandparent (GP) is this.parent_. Let's have
    824910    variable names respecting the truth.
    825911  */
     
    840926  // THINK: these are not going to happen (hopefully), but if they do, what
    841927  // should we do???
    842   if (model.isAbandoned()) {
    843     printf("WARNING: model.isAbandoned() true!\n");
     928  if (lpres.isAbandoned) {
     929    printf("WARNING: lpres.isAbandoned true!\n");
    844930    return false;
    845931  }
    846   if (model.isProvenDualInfeasible()) {
    847     printf("WARNING: model.isProvenDualInfeasible() true!\n");
     932  if (lpres.isProvenDualInfeasible) {
     933    printf("WARNING: lpres.isProvenDualInfeasible true!\n");
    848934    return false;
    849935  }
    850   if (model.isPrimalObjectiveLimitReached()) {
    851     printf("WARNING: model.isPrimalObjectiveLimitReached() true!\n");
     936  if (lpres.isPrimalObjectiveLimitReached) {
     937    printf("WARNING: lpres.isPrimalObjectiveLimitReached true!\n");
    852938    return false;
    853939  }
     
    858944  }
    859945 
    860   const double direction = model.getObjSense();
     946  const double direction = lpres.getObjSense;
    861947
    862948  const int GP_brvar = grandparent.variable_;
     
    864950  const bool parent_is_down_child = parent_id == grandparent.child_down_;
    865951
    866   if (model.isProvenOptimal() ||
    867       model.isDualObjectiveLimitReached() ||
    868       model.isIterationLimitReached()) {
     952  if (lpres.isProvenOptimal ||
     953      lpres.isDualObjectiveLimitReached ||
     954      lpres.isIterationLimitReached) {
    869955    // Dual feasible, and in this case we don't care how we have
    870956    // stopped (iteration limit, obj val limit, time limit, optimal solution,
    871957    // etc.), we can just look at the reduced costs to figure out if the
    872958    // grandparent is irrelevant. Remember, we are using dual simplex!
    873     double djValue = model.getReducedCost()[GP_brvar_fullid]*direction;
     959    double djValue = lpres.getReducedCost[GP_brvar_fullid]*direction;
    874960    if (djValue > 1.0e-6) {
    875961      // wants to go down
     
    877963        return true;
    878964      }
    879       if (model.getColLower()[GP_brvar_fullid] > std::ceil(grandparent.value_)) {
     965      if (lpres.getColLower[GP_brvar_fullid] > std::ceil(grandparent.value_)) {
    880966        return true;
    881967      }
     
    885971        return true;
    886972      }
    887       if (model.getColUpper()[GP_brvar_fullid] < std::floor(grandparent.value_)) {
     973      if (lpres.getColUpper[GP_brvar_fullid] < std::floor(grandparent.value_)) {
    888974        return true;
    889975      }
     
    891977    return false;
    892978  } else {
    893     assert(model.isProvenPrimalInfeasible());
     979    assert(lpres.isProvenPrimalInfeasible);
    894980    return false;
    895981    const int greatgrandparent_id = grandparent.parent_;
     
    903989    */
    904990   
    905     const double* dj = model.getReducedCost();
    906     const double* obj = model.getObjCoefficients();
     991    const double* dj = lpres.getReducedCost;
     992    const double* obj = lpres.getObjCoefficients;
    907993    const double yA_x = dj[x] - obj[x];
    908994
     
    913999          return true;
    9141000        }
    915         const double yb_plus_rl_minus_su = compute_val_for_ray(model);
    916         if (yb_plus_rl_minus_su < 1e-8) {
    917           printf("WARNING: yb_plus_rl_minus_su is not positive!\n");
     1001        if (lpres.yb_plus_rl_minus_su < 1e-8) {
     1002          printf("WARNING: lpres.yb_plus_rl_minus_su is not positive!\n");
    9181003          return false;
    9191004        }
    920         const double max_u_x = yb_plus_rl_minus_su / s_x + model.getColUpper()[x];
     1005        const double max_u_x = lpres.yb_plus_rl_minus_su / s_x + lpres.getColUpper[x];
    9211006        const double u_x_without_GP = greatgrandparent_id >= 0 ?
    9221007          branchingTree.nodes_[greatgrandparent_id].upper_[GP_brvar] :
     
    9281013          return true;
    9291014        }
    930         const double yb_plus_rl_minus_su = compute_val_for_ray(model);
    931         if (yb_plus_rl_minus_su < 1e-8) {
    932           printf("WARNING: yb_plus_rl_minus_su is not positive!\n");
     1015        if (lpres.yb_plus_rl_minus_su < 1e-8) {
     1016          printf("WARNING: lpres.yb_plus_rl_minus_su is not positive!\n");
    9331017          return false;
    9341018        }
    935         const double min_l_x = - yb_plus_rl_minus_su / r_x + model.getColLower()[x];
     1019        const double min_l_x = - lpres.yb_plus_rl_minus_su / r_x + lpres.getColLower[x];
    9361020        const double l_x_without_GP = greatgrandparent_id >= 0 ?
    9371021          branchingTree.nodes_[greatgrandparent_id].lower_[GP_brvar] :
     
    9451029          return true;
    9461030        }
    947         const double yb_plus_rl_minus_su = compute_val_for_ray(model);
    948         if (yb_plus_rl_minus_su > -1e-8) {
    949           printf("WARNING: yb_plus_rl_minus_su is not negative!\n");
     1031        if (lpres.yb_plus_rl_minus_su > -1e-8) {
     1032          printf("WARNING: lpres.yb_plus_rl_minus_su is not negative!\n");
    9501033          return false;
    9511034        }
    952         const double max_u_x = yb_plus_rl_minus_su / s_x + model.getColUpper()[x];
     1035        const double max_u_x = lpres.yb_plus_rl_minus_su / s_x + lpres.getColUpper[x];
    9531036        const double u_x_without_GP = greatgrandparent_id >= 0 ?
    9541037          branchingTree.nodes_[greatgrandparent_id].upper_[GP_brvar] :
     
    9601043          return true;
    9611044        }
    962         const double yb_plus_rl_minus_su = compute_val_for_ray(model);
    963         if (yb_plus_rl_minus_su > -1e-8) {
    964           printf("WARNING: yb_plus_rl_minus_su is not negative!\n");
     1045        if (lpres.yb_plus_rl_minus_su > -1e-8) {
     1046          printf("WARNING: lpres.yb_plus_rl_minus_su is not negative!\n");
    9651047          return false;
    9661048        }
    967         const double min_l_x = - yb_plus_rl_minus_su / r_x + model.getColLower()[x];
     1049        const double min_l_x = - lpres.yb_plus_rl_minus_su / r_x + lpres.getColLower[x];
    9681050        const double l_x_without_GP = greatgrandparent_id >= 0 ?
    9691051          branchingTree.nodes_[greatgrandparent_id].lower_[GP_brvar] :
     
    9781060
    9791061void
     1062DBVectorNode::adjustBounds(int subroot, int brvar, int brlb, int brub)
     1063{
     1064  assert(subroot != -1);
     1065  DBNodeSimple& node = nodes_[subroot];
     1066  // Take the intersection of brvar's bounds in node and those passed in
     1067  brub = CoinMin(brub, node.upper_[brvar]);
     1068  brlb = CoinMax(brlb, node.lower_[brvar]);
     1069  if (brub < brlb) {
     1070    // This node became infeasible. Get rid of it and of its descendants
     1071    removeSubTree(subroot);
     1072    return;
     1073  }
     1074  if (node.variable_ == brvar) {
     1075    if (node.value_ < brlb) {
     1076      // child_down_ just became infeasible. Just cut out the current node
     1077      // together with its child_down_ from the tree and hang child_up_ on the
     1078      // parent of this node.
     1079      if (node.child_down_ >= 0) {
     1080        removeSubTree(node.child_down_);
     1081      }
     1082      const int parent_id = node.parent_;
     1083      const int child_remains = node.child_up_;
     1084      if (nodes_[parent_id].child_down_ == subroot) {
     1085        nodes_[parent_id].child_down_ = child_remains;
     1086      } else {
     1087        nodes_[parent_id].child_up_ = child_remains;
     1088      }
     1089      removeNode(subroot);
     1090      if (child_remains >= 0) {
     1091        nodes_[child_remains].parent_ = parent_id;
     1092        adjustBounds(child_remains, brvar, brlb, brub);
     1093      }
     1094      return;
     1095    }
     1096    if (node.value_ > brub) {
     1097      // child_up_ just became infeasible. Just cut out the current node
     1098      // together with its child_down_ from the tree and hang child_down_ on
     1099      // the parent of this node.
     1100      if (node.child_up_ >= 0) {
     1101        removeSubTree(node.child_up_);
     1102      }
     1103      const int parent_id = node.parent_;
     1104      const int child_remains = node.child_down_;
     1105      if (nodes_[parent_id].child_down_ == subroot) {
     1106        nodes_[parent_id].child_down_ = child_remains;
     1107      } else {
     1108        nodes_[parent_id].child_up_ = child_remains;
     1109      }
     1110      removeNode(subroot);
     1111      if (child_remains >= 0) {
     1112        nodes_[child_remains].parent_ = parent_id;
     1113        adjustBounds(child_remains, brvar, brlb, brub);
     1114      }
     1115      return;
     1116    }
     1117    // Now brlb < node.value_ < brub (value_ is fractional)
     1118  }
     1119  node.lower_[brvar] = brlb;
     1120  node.upper_[brvar] = brub;
     1121  if (node.child_down_ >= 0) {
     1122    adjustBounds(node.child_down_, brvar, brlb, brub);
     1123  }
     1124  if (node.child_up_ >= 0) {
     1125    adjustBounds(node.child_up_, brvar, brlb, brub);
     1126  }
     1127}
     1128
     1129void
    9801130DBVectorNode::moveNodeUp(const int* which,
    9811131                         OsiSolverInterface& model, DBNodeSimple & node)
    9821132{
    9831133  /*
    984     The current node ('this') is really the parent (P) and the solution in the
    985     model represents the child. The grandparent (GP) is this.parent_. Let's have
    986     variable names respecting the truth.
     1134    The current node ('this') is really the parent (P). The grandparent (GP)
     1135    is this.parent_. Let's have variable names respecting the truth.
    9871136  */
    988   const bool childWasInfeasible = model.isProvenPrimalInfeasible();
    9891137  const int parent_id = node.node_id_;
    9901138  DBNodeSimple& parent = nodes_[parent_id];
     
    9961144  const bool parent_is_down_child = parent_id == grandparent.child_down_;
    9971145
    998 #if defined(DEBUG_DYNAMIC_BRANCHING) && DEBUG_DYNAMIC_BRANCHING >= 1000
    999   printf("entered moveNodeUp\n");
    1000   printf("parent_id %d grandparent_id %d greatgrandparent_id %d\n",
    1001          parent_id, grandparent_id, greatgrandparent_id);
    1002   printf("parent.way_ %d\n", parent.way_);
     1146#if defined(DEBUG_DYNAMIC_BRANCHING)
     1147  if (dyn_debug >= 100) {
     1148    printf("entered moveNodeUp\n");
     1149    printf("parent_id %d grandparent_id %d greatgrandparent_id %d\n",
     1150           parent_id, grandparent_id, greatgrandparent_id);
     1151    printf("parent.way_ %d\n", parent.way_);
     1152  }
    10031153#endif
    10041154
     
    10121162  const bool child_to_move_is_processed = parent.bothChildDone();
    10131163
    1014 #if defined(DEBUG_DYNAMIC_BRANCHING) && DEBUG_DYNAMIC_BRANCHING >= 1000
    1015   printf("parent_is_down_child %d down_child_stays_with_parent %d child_to_move %d\n", parent_is_down_child, down_child_stays_with_parent, child_to_move);
     1164#if defined(DEBUG_DYNAMIC_BRANCHING)
     1165  if (dyn_debug >= 1000) {
     1166    printf("parent_is_down_child %d down_child_stays_with_parent %d child_to_move %d\n", parent_is_down_child, down_child_stays_with_parent, child_to_move);
     1167  }
    10161168#endif
    10171169
     
    10341186  }
    10351187
    1036 #if defined(DEBUG_DYNAMIC_BRANCHING) && DEBUG_DYNAMIC_BRANCHING >= 1000
    1037   printf("after exchange\n");
    1038   printf("parent.parent_ %d parent.child_down_ %d parent.child_up_ %d\n",
    1039          parent.parent_, parent.child_down_, parent.child_up_);
    1040   printf("grandparent.parent_ %d grandparent.child_down_ %d grandparent.child_up_ %d\n",
    1041          grandparent.parent_, grandparent.child_down_, grandparent.child_up_);
    1042   if (greatgrandparent_id >= 0) {
    1043     DBNodeSimple& greatgrandparent = nodes_[greatgrandparent_id];
    1044     printf("greatgrandparent.parent_ %d greatgrandparent.child_down_ %d greatgrandparent.child_up_ %d\n",
    1045          greatgrandparent.parent_, greatgrandparent.child_down_, greatgrandparent.child_up_);
    1046   }
    1047   printf("exiting moveNodeUp\n");
     1188#if defined(DEBUG_DYNAMIC_BRANCHING)
     1189  if (dyn_debug >= 1000) {
     1190    printf("after exchange\n");
     1191    printf("parent.parent_ %d parent.child_down_ %d parent.child_up_ %d\n",
     1192           parent.parent_, parent.child_down_, parent.child_up_);
     1193    printf("grandparent.parent_ %d grandparent.child_down_ %d grandparent.child_up_ %d\n",
     1194           grandparent.parent_, grandparent.child_down_, grandparent.child_up_);
     1195    if (greatgrandparent_id >= 0) {
     1196      DBNodeSimple& greatgrandparent = nodes_[greatgrandparent_id];
     1197      printf("greatgrandparent.parent_ %d greatgrandparent.child_down_ %d greatgrandparent.child_up_ %d\n",
     1198             greatgrandparent.parent_, greatgrandparent.child_down_, greatgrandparent.child_up_);
     1199    }
     1200    printf("exiting moveNodeUp\n");
     1201  }
    10481202#endif
    10491203
     
    11071261    const int old_upper = grandparent.upper_[GP_brvar];
    11081262    parent.upper_[GP_brvar] = old_upper;
    1109     model.setColUpper(which[GP_brvar], old_upper);
     1263    if ((GP_brvar != parent.variable_) ||
     1264        (GP_brvar == parent.variable_ && !down_child_stays_with_parent)) {
     1265      model.setColUpper(which[GP_brvar], old_upper);
     1266    }
    11101267  } else {
    11111268    const int old_lower = grandparent.lower_[GP_brvar];
    11121269    parent.lower_[GP_brvar] = old_lower;
    1113     model.setColLower(which[GP_brvar], old_lower);
    1114   }
    1115 #if 0
    1116   // THINK: this might not be necessary
    1117   model.resolve();
    1118 #endif
     1270    if ((GP_brvar != parent.variable_) ||
     1271        (GP_brvar == parent.variable_ && down_child_stays_with_parent)) {
     1272      model.setColLower(which[GP_brvar], old_lower);
     1273    }
     1274  }
    11191275
    11201276  // Now add the branching var bound change of P to GP and all of its
    11211277  // descendant
    1122   const int P_brvar = parent.variable_;
    1123   const double P_value = parent.value_;
    1124   int new_bd;
    1125 #if 0
    11261278  if (down_child_stays_with_parent) {
    1127     // Former up child of P is now the down child of GP, so we need to change
    1128     // bounds of GP, its up child and descendants of that one.
    1129     new_bd = (int)std::ceil(P_value);
    1130     grandparent.lower_[P_brvar] = new_bd;
    1131     changeDescendantBounds(grandparent.child_up_,
    1132                            true /*lower bd*/, P_brvar, new_bd);
     1279    adjustBounds(grandparent_id, parent.variable_,
     1280                 (int)ceil(parent.value_), parent.upper_[parent.variable_]);
    11331281  } else {
    1134     // Former down child of P is now the up child of GP, so we need to change
    1135     // bounds of GP, its down child and descendants of that one.
    1136     new_bd = (int)floor(P_value);
    1137     grandparent.upper_[P_brvar] = new_bd;
    1138     changeDescendantBounds(grandparent.child_down_,
    1139                            false /*lower bd*/, P_brvar, new_bd);
    1140   }
    1141 #else
    1142   if (down_child_stays_with_parent) {
    1143     // Former up child of P is now the down (or up) child of GP,
    1144     // so we need to change
    1145     // bounds of GP, its up (or down) child and descendants of that one.
    1146     new_bd = (int)std::ceil(P_value);
    1147     grandparent.lower_[P_brvar] = new_bd;
    1148     if(parent_is_down_child)
    1149       changeDescendantBounds(grandparent.child_up_,
    1150                              true /*lower bd*/, P_brvar, new_bd);
    1151     else
    1152       changeDescendantBounds(grandparent.child_down_,
    1153                              true /*lower bd*/, P_brvar, new_bd);
    1154   } else {
    1155     // Former down child of P is now the up (or down) child of GP,
    1156     // so we need to change
    1157     // bounds of GP, its down (or up) child and descendants of that one.
    1158     new_bd = (int)floor(P_value);
    1159     grandparent.upper_[P_brvar] = new_bd;
    1160     if(parent_is_down_child)
    1161       changeDescendantBounds(grandparent.child_up_,
    1162                              false /*lower bd*/, P_brvar, new_bd);
    1163     else
    1164       changeDescendantBounds(grandparent.child_down_,
    1165                              false /*lower bd*/, P_brvar, new_bd);
    1166   }
    1167 #endif
    1168 #if defined(DEBUG_DYNAMIC_BRANCHING) && DEBUG_DYNAMIC_BRANCHING >= 1000
    1169   if (childWasInfeasible) {
    1170     model.resolve();
    1171     assert(model.isProvenPrimalInfeasible());
    1172   }
    1173 #endif
     1282    adjustBounds(grandparent_id, parent.variable_,
     1283                 parent.lower_[parent.variable_], (int)floor(parent.value_));
     1284  }
    11741285}
    11751286
    11761287void
    1177 DBVectorNode::changeDescendantBounds(const int node_id, const bool lower_bd,
    1178                                      const int brvar, const int new_bd)
    1179 {
    1180   if (node_id == -1) {
     1288DBVectorNode::checkNode(int node) const
     1289{
     1290  if (node == -1) {
    11811291    return;
    11821292  }
    1183   changeDescendantBounds(nodes_[node_id].child_down_, lower_bd, brvar, new_bd);
    1184   changeDescendantBounds(nodes_[node_id].child_up_, lower_bd, brvar, new_bd);
    1185   if (lower_bd) {
    1186     nodes_[node_id].lower_[brvar] = new_bd;
    1187   } else {
    1188     nodes_[node_id].upper_[brvar] = new_bd;
    1189   }
    1190 }
    1191 
     1293  const DBNodeSimple* n = nodes_ + node;
     1294  const DBNodeSimple* p = nodes_ + n->parent_;
     1295  for (int i = n->numberIntegers_-1; i >= 0; --i) {
     1296    if (i == p->variable_) {
     1297      if (node == p->child_down_) {
     1298        assert(p->lower_[i] <= n->lower_[i]);
     1299        assert((int)(floor(p->value_)) == n->upper_[i]);
     1300      } else {
     1301        assert((int)(ceil(p->value_)) == n->lower_[i]);
     1302        assert(p->upper_[i] >= n->upper_[i]);
     1303      }
     1304    } else {
     1305      assert(p->lower_[i] <= n->lower_[i]);
     1306      assert(p->upper_[i] >= n->upper_[i]);
     1307    }
     1308  }
     1309  checkNode(n->child_down_);
     1310  checkNode(n->child_up_);
     1311}
     1312
     1313void
     1314DBVectorNode::checkTree() const
     1315{
     1316  // find the root
     1317  int root = first_;
     1318  while (true) {
     1319    if (nodes_[root].parent_ == -1) {
     1320      break;
     1321    }
     1322  }
     1323  checkNode(nodes_[root].child_down_);
     1324  checkNode(nodes_[root].child_up_);
     1325}
     1326
     1327std::string getTree(DBVectorNode& branchingTree)
     1328{
     1329  std::string tree;
     1330  char line[1000];
     1331  for(int k=0; k<branchingTree.size_; k++) {
     1332    DBNodeSimple& node = branchingTree.nodes_[k];
     1333    sprintf(line, "%d %d %d %d %f %d %d 0x%x %d %d\n",
     1334            k, node.node_id_, node.parent_, node.variable_,
     1335            node.value_, node.lower_[node.variable_],
     1336            node.upper_[node.variable_], node.way_,
     1337            node.child_down_, node.child_up_);
     1338    tree += line;
     1339  }
     1340  return tree;
     1341}
     1342
     1343void printTree(const std::string& tree, int levels)
     1344{
     1345  size_t pos = tree.size();
     1346  for (int i = levels-1; i >= 0; --i) {
     1347    pos = tree.rfind('\n', pos-1);
     1348    if (pos == std::string::npos) {
     1349      pos = 0;
     1350      break;
     1351    }
     1352  }
     1353  printf("%s", tree.c_str() + (pos+1));
     1354}
     1355
     1356void printChain(DBVectorNode& branchingTree, int k)
     1357{
     1358  while (k != -1) {
     1359    DBNodeSimple& node = branchingTree.nodes_[k];
     1360    printf("   %d %d %d %d %f %d %d 0x%x %d %d %c %c\n",
     1361           k, node.node_id_, node.parent_, node.variable_,
     1362           node.value_, node.lower_[node.variable_],
     1363           node.upper_[node.variable_], node.way_,
     1364           node.child_down_, node.child_up_,
     1365           node.strong_branching_fixed_vars_ ? 'T' : 'F',
     1366           node.reduced_cost_fixed_vars_ ? 'T' : 'F');
     1367    k = node.parent_;
     1368  }
     1369}
    11921370
    11931371// Invoke solver's built-in enumeration algorithm
     
    12551433    // while until nothing on stack
    12561434    while (branchingTree.size()) {
    1257 #if defined(DEBUG_DYNAMIC_BRANCHING) && DEBUG_DYNAMIC_BRANCHING >= 1000
    1258       printf("branchingTree.size = %d %d\n",branchingTree.size(),branchingTree.size_);
    1259       printf("i node_id parent child_down child_up\n");
    1260       for(int k=0; k<branchingTree.size_; k++) {
    1261         DBNodeSimple& node = branchingTree.nodes_[k];
    1262         printf("%d %d %d %d %d\n",k, node.node_id_, node.parent_,
    1263                node.child_down_, node.child_up_);
     1435#if defined(DEBUG_DYNAMIC_BRANCHING)
     1436      if (dyn_debug >= 1000) {
     1437        printf("branchingTree.size = %d %d\n",branchingTree.size(),branchingTree.size_);
     1438        printf("i node_id parent child_down child_up\n");
     1439        for(int k=0; k<branchingTree.size_; k++) {
     1440          DBNodeSimple& node = branchingTree.nodes_[k];
     1441          printf("%d %d %d %d %d\n",k, node.node_id_, node.parent_,
     1442                 node.child_down_, node.child_up_);
     1443        }
    12641444      }
    12651445#endif
     
    12681448      int kNode = branchingTree.chosen_;
    12691449      branchingTree.pop_back();
    1270 #if defined(DEBUG_DYNAMIC_BRANCHING) && DEBUG_DYNAMIC_BRANCHING >= 1000
    1271       printf("Deleted current parent %d %d\n",branchingTree.size(),branchingTree.size_);
     1450#if defined(DEBUG_DYNAMIC_BRANCHING)
     1451      if (dyn_debug >= 1000) {
     1452        printf("Deleted current parent %d %d\n",branchingTree.size(),branchingTree.size_);
     1453      }
    12721454#endif
    12731455      assert (! node.bothChildDone());
     
    12881470                 <<" iterations and "<<numberNodes<<" nodes"
    12891471                 <<std::endl;
     1472#if defined(DEBUG_DYNAMIC_BRANCHING)
     1473        if (dyn_debug >= 1) {
     1474          branchingTree.checkTree();
     1475        }
     1476#endif
     1477
    12901478      } else {
    12911479        // branch - do bounds
     
    13061494        // to be done. We want the whole tree all the time.
    13071495        branchingTree.push_back(node);
    1308 #if defined(DEBUG_DYNAMIC_BRANCHING) && DEBUG_DYNAMIC_BRANCHING >= 1000
    1309       printf("Added current parent %d %d\n",branchingTree.size(),branchingTree.size_);
     1496#if defined(DEBUG_DYNAMIC_BRANCHING)
     1497        if (dyn_debug >= 1000) {
     1498          printf("Added current parent %d %d\n",branchingTree.size(),branchingTree.size_);
     1499        }
    13101500#endif
    13111501       
     
    13601550        }
    13611551
    1362         bool canSwitch = node.canSwitchParentWithGrandparent(which, model,
     1552        LPresult lpres(model);
     1553
     1554        bool canSwitch = node.canSwitchParentWithGrandparent(which, lpres,
    13631555                                                             originalLower,
    13641556                                                             originalUpper,
    13651557                                                             branchingTree);
    13661558        int cnt = 0;
    1367 #if defined(DEBUG_DYNAMIC_BRANCHING) && DEBUG_DYNAMIC_BRANCHING >= 100
     1559
     1560#if defined(DEBUG_DYNAMIC_BRANCHING)
     1561        if (dyn_debug >= 1) {
     1562          branchingTree.checkTree();
     1563        }
     1564#endif
     1565
     1566#if defined(DEBUG_DYNAMIC_BRANCHING)
    13681567        std::string tree0;
    1369         char line[1000];
    1370         if (canSwitch) {
    1371           for(int k=0; k<branchingTree.size_; k++) {
    1372             DBNodeSimple& node = branchingTree.nodes_[k];
    1373             sprintf(line, "%d %d %d %d %d %d 0x%x %d %d\n",
    1374                     k, node.node_id_, node.parent_, node.variable_,
    1375                     node.lower_[node.variable_],
    1376                     node.upper_[node.variable_], node.way_,
    1377                     node.child_down_, node.child_up_);
    1378             tree0 += line;
     1568        if (dyn_debug >= 10) {
     1569          if (canSwitch) {
     1570            tree0 = getTree(branchingTree);
    13791571          }
    13801572        }
     
    13831575        while (canSwitch) {
    13841576          branchingTree.moveNodeUp(which, model, node);
    1385           canSwitch = node.canSwitchParentWithGrandparent(which, model,
     1577          canSwitch = node.canSwitchParentWithGrandparent(which, lpres,
    13861578                                                          originalLower,
    13871579                                                          originalUpper,
    13881580                                                          branchingTree);
     1581#if defined(DEBUG_DYNAMIC_BRANCHING)
     1582          if (dyn_debug >= 1) {
     1583            branchingTree.checkTree();
     1584          }
     1585#endif
    13891586          ++cnt;
    13901587        }
    1391 #if defined(DEBUG_DYNAMIC_BRANCHING) && DEBUG_DYNAMIC_BRANCHING >= 100
    13921588        if (cnt > 0) {
    1393           std::string tree1;
    1394           for(int k=0; k<branchingTree.size_; k++) {
    1395             DBNodeSimple& node = branchingTree.nodes_[k];
    1396             sprintf(line, "%d %d %d %d %d %d 0x%x %d %d\n",
    1397                     k, node.node_id_, node.parent_, node.variable_,
    1398                     node.lower_[node.variable_],
    1399                     node.upper_[node.variable_], node.way_,
    1400                     node.child_down_, node.child_up_);
    1401             tree1 += line;
     1589          model.resolve();
     1590          // This is horribly looking... Get rid of it when properly debugged...
     1591          assert(lpres.isAbandoned == model.isAbandoned());
     1592          assert(lpres.isDualObjectiveLimitReached == model.isDualObjectiveLimitReached());
     1593          assert(lpres.isDualObjectiveLimitReached ||
     1594                 (lpres.isProvenOptimal == model.isProvenOptimal()));
     1595          assert(lpres.isDualObjectiveLimitReached ||
     1596                 (lpres.isProvenPrimalInfeasible == model.isProvenPrimalInfeasible()));
     1597        }
     1598           
     1599#if defined(DEBUG_DYNAMIC_BRANCHING)
     1600        if (dyn_debug >= 10) {
     1601          if (cnt > 0) {
     1602            std::string tree1 = getTree(branchingTree);
     1603            printf("=====================================\n");
     1604            printf("It can move node %i up. way_: 0x%x   brvar: %i\n",
     1605                   node.node_id_, node.way_, node.variable_);
     1606            printTree(tree0, cnt+10);
     1607            printf("=====================================\n");
     1608            printf("Finished moving the node up by %i levels.\n", cnt);
     1609            printTree(tree1, cnt+10);
     1610            printf("=====================================\n");
    14021611          }
    1403           printf("=====================================\n");
    1404           printf("It can move node %i up. way_: 0x%x   brvar: %i\n",
    1405                  node.node_id_, node.way_, node.variable_);
    1406           printf("i node_id parent brvar lb ub way child_down child_up\n");
    1407           size_t pos = tree0.size();
    1408           for (int ii = cnt + 10; ii >= 0; --ii) {
    1409             pos = tree0.rfind('\n', pos-1);
    1410             if (pos == std::string::npos) {
    1411               pos = 0;
    1412               break;
    1413             }
    1414           }
    1415           printf("%s", tree0.c_str() + (pos+1));
    1416           printf("=====================================\n");
    1417           printf("Finished moving the node up by %i levels.\n", cnt);
    1418           pos = tree1.size();
    1419           for (int ii = cnt + 10; ii >= 0; --ii) {
    1420             pos = tree1.rfind('\n', pos-1);
    1421             if (pos == std::string::npos) {
    1422               pos = 0;
    1423               break;
    1424             }
    1425           }
    1426           printf("%s", tree1.c_str() + (pos+1));
    1427           printf("=====================================\n");
    14281612        }
    14291613#endif
     
    14481632          // push on stack
    14491633          branchingTree.push_back(newNode);
    1450 #if defined(DEBUG_DYNAMIC_BRANCHING) && DEBUG_DYNAMIC_BRANCHING >= 1000
    1451           printf("Added current child %d %d\n",branchingTree.size(),branchingTree.size_);
     1634#if defined(DEBUG_DYNAMIC_BRANCHING)
     1635          if (dyn_debug >= 1000) {
     1636            printf("Added current child %d %d\n",branchingTree.size(),branchingTree.size_);
     1637          }
    14521638#endif
    14531639          if (branchingTree.nodes_[kNode].workingOnDownChild()) {
     
    14731659#endif
    14741660        }
     1661#if defined(DEBUG_DYNAMIC_BRANCHING)
     1662        if (dyn_debug >= 1) {
     1663          branchingTree.checkTree();
     1664        }
     1665#endif
    14751666      }
    14761667    }
Note: See TracChangeset for help on using the changeset viewer.