source: branches/dynamicbranching/dynamicbranching.cpp @ 992

Last change on this file since 992 was 992, checked in by ladanyi, 12 years ago

one more bug...

File size: 50.6 KB
Line 
1/**********************************************************************
2TODO:
3
4If we swap two nodes (parent and grandparent) then we should check if anything
5has been explored under GP after the swap, and if not, just get rid of GP and
6everything below.
7
8If strong branching fixed anything in the grandparent we may still swap it
9with the parent (we don't do it yet, see the test on
10strong_branching_fixed_vars_), but those fixings must be moved to the parent as
11well.
12
13Same thing for locally valid cuts, if GP has them and GP is switched with P
14then locally valid cuts must be treated as generated in P.
15
16Same for reduced cost fixing :-(. If P has done any then P and GP can't be
17switched.
18
19Maybe the best solution is to create local information (strong branching, rc
20fixing, local cuts, etc.) only every so often, say every 10th level. Then that
21would block switches, but everywhere else we would be safe. Good question
22what is worth more: the switches or the local info.
23
24Alternative solution: Do not add local info to the tree, but keep it in a set
25of excluded clauses ala Todias Achterberg: Conflict Analysis in Mixed Integer
26Programming.
27
28We may want to disable fixing by strong branching completely and if such
29situation occurs then simply do the branch and one side will be fathomed
30immediately and we can try to switch.
31
32Bound modifications when nodes are swapped could be avoided if always start
33from original bounds and go back to root to apply all branching decisions. On
34the other hand, reduced cost fixing would be lost. And so would fixings by
35strong branching. Although both could be stored in an array keeping track of
36changes implied by the branching decisions.
37
38**********************************************************************
39
40
41**********************************************************************/
42#include "CoinTime.hpp"
43#include "OsiClpSolverInterface.hpp"
44
45#define FUNNY_BRANCHING 1
46#define DEBUG_DYNAMIC_BRANCHING
47
48#ifdef DEBUG_DYNAMIC_BRANCHING
49int dyn_debug = 1;
50#endif
51
52// below needed for pathetic branch and bound code
53#include <vector>
54#include <map>
55
56using namespace std;
57
58class DBVectorNode;
59
60class LPresult {
61private:
62  LPresult(const LPresult& rhs);
63  LPresult& operator=(const LPresult& rhs);
64  void gutsOfConstructor(const OsiSolverInterface& model);
65public:
66  LPresult(const OsiSolverInterface& model);
67  ~LPresult();
68public:
69  bool isAbandoned;
70  bool isProvenDualInfeasible;
71  bool isPrimalObjectiveLimitReached;
72  bool isProvenOptimal;
73  bool isDualObjectiveLimitReached;
74  bool isIterationLimitReached;
75  bool isProvenPrimalInfeasible;
76  double getObjSense;
77  double getObjValue;
78  double* getReducedCost;
79  double* getColLower;
80  double* getColUpper;
81  double* getObjCoefficients;
82  double yb_plus_rl_minus_su;
83};
84
85void
86LPresult::gutsOfConstructor(const OsiSolverInterface& model)
87{
88  isAbandoned = model.isAbandoned();
89  isProvenDualInfeasible = model.isProvenDualInfeasible();
90  isPrimalObjectiveLimitReached = model.isPrimalObjectiveLimitReached();
91  isProvenOptimal = model.isProvenOptimal();
92  isDualObjectiveLimitReached = model.isDualObjectiveLimitReached();
93  isIterationLimitReached = model.isIterationLimitReached();
94  isProvenPrimalInfeasible = model.isProvenPrimalInfeasible();
95  getObjSense = model.getObjSense();
96  getObjValue = model.getObjValue();
97
98  getReducedCost = new double[model.getNumCols()];
99  CoinDisjointCopyN(model.getReducedCost(), model.getNumCols(), getReducedCost);
100  getColLower = new double[model.getNumCols()];
101  CoinDisjointCopyN(model.getColLower(), model.getNumCols(), getColLower);
102  getColUpper = new double[model.getNumCols()];
103  CoinDisjointCopyN(model.getColUpper(), model.getNumCols(), getColUpper);
104
105  getObjCoefficients = NULL;
106  yb_plus_rl_minus_su = 0;
107 
108  if (!isProvenOptimal &&
109      !isDualObjectiveLimitReached &&
110      !isIterationLimitReached &&
111      !isAbandoned &&
112      !isProvenDualInfeasible &&
113      !isPrimalObjectiveLimitReached) {
114    assert(isProvenPrimalInfeasible);
115    getObjCoefficients = new double[model.getNumCols()];
116    CoinDisjointCopyN(model.getObjCoefficients(), model.getNumCols(), getObjCoefficients);
117    const std::vector<double*> dual_rays = model.getDualRays(1);
118    if (dual_rays.size() == 0) {
119      printf("WARNING: LP is infeas, but no dual ray is returned!\n");
120      return;
121    }
122    const double* dual_ray = dual_rays[0];
123    const double direction = model.getObjSense();
124    const double* rub = model.getRowUpper();
125    const double* rlb = model.getRowLower();
126    const double* cub = model.getColUpper();
127    const double* clb = model.getColLower();
128    const double* dj  = model.getReducedCost();
129    const double* obj = model.getObjCoefficients();
130    const int numRows = model.getNumRows();
131    const int numCols = model.getNumCols();
132    for (int i = 0; i < numRows; ++i) {
133      const double ray_i = dual_ray[i];
134      if (ray_i > 1e-6) {
135        yb_plus_rl_minus_su += ray_i*rlb[i];
136      } else if (ray_i < 1e-6) {
137        yb_plus_rl_minus_su += ray_i*rub[i];
138      }
139    }
140    for (int j = 0; j < numCols; ++j) {
141      const double yA_j = dj[j] - obj[j];
142      if (direction * yA_j > 1e-6) {
143        yb_plus_rl_minus_su -= yA_j*cub[j];
144      } else if (direction * yA_j < -1e-6) {
145        yb_plus_rl_minus_su -= yA_j*clb[j];
146      }
147    }
148    for (int i = dual_rays.size()-1; i >= 0; --i) {
149      delete[] dual_rays[i];
150    }
151  }
152}
153 
154LPresult::LPresult(const OsiSolverInterface& model)
155{
156  gutsOfConstructor(model);
157}
158
159LPresult::~LPresult()
160{
161  delete[] getReducedCost;
162  delete[] getColLower;
163  delete[] getColUpper;
164  delete[] getObjCoefficients;
165}
166
167// Trivial class for Branch and Bound
168
169class DBNodeSimple  {
170public:
171  enum DBNodeWay {
172    WAY_UNSET=0x00,
173    WAY_DOWN_UP__NOTHING_DONE=0x10,
174    WAY_DOWN_UP__DOWN_DONE=0x11,
175    WAY_DOWN_UP__BOTH_DONE=0x13,
176    WAY_DOWN_UP__=0x10,
177    WAY_UP_DOWN__NOTHING_DONE=0x20,
178    WAY_UP_DOWN__UP_DONE=0x22,
179    WAY_UP_DOWN__BOTH_DONE=0x23,
180    WAY_UP_DOWN__=0x20,
181    WAY_DOWN_CURRENT=0x01,
182    WAY_UP_CURRENT=0x02,
183    WAY_BOTH_DONE=0x03,
184   
185  };
186
187private:
188  void gutsOfCopy(const DBNodeSimple& rhs);
189  void gutsOfConstructor (OsiSolverInterface &model,
190                          int numberIntegers, int * integer,
191                          CoinWarmStart * basis);
192  void gutsOfDestructor(); 
193 
194public:
195   
196  // Default Constructor
197  DBNodeSimple ();
198
199  // Constructor from current state (and list of integers)
200  // Also chooses branching variable (if none set to -1)
201  DBNodeSimple (OsiSolverInterface &model,
202                int numberIntegers, int * integer,
203                CoinWarmStart * basis);
204  // Copy constructor
205  DBNodeSimple ( const DBNodeSimple &);
206   
207  // Assignment operator
208  DBNodeSimple & operator=( const DBNodeSimple& rhs);
209
210  // Destructor
211  ~DBNodeSimple ();
212
213  // Extension - true if other extension of this
214  bool extension(const DBNodeSimple & other,
215                 const double * originalLower,
216                 const double * originalUpper) const;
217  // Tests if we can switch this node (this is the parent) with its parent
218  bool canSwitchParentWithGrandparent(const int* which,
219                                      const LPresult& lpres,
220                                      const int * original_lower,
221                                      const int * original_upper,
222                                      DBVectorNode & branchingTree);
223  inline bool bothChildDone() const {
224    return (way_ & WAY_BOTH_DONE) == WAY_BOTH_DONE;
225  }
226  inline bool workingOnDownChild() const {
227    return (way_ == WAY_DOWN_UP__DOWN_DONE || way_ == WAY_UP_DOWN__BOTH_DONE);
228  }
229  /* Return whether the child that will be now processed is the down branch or
230     not. */
231  inline bool advanceWay() { 
232    switch (way_) {
233    case WAY_UNSET:
234    case WAY_DOWN_UP__BOTH_DONE:
235    case WAY_UP_DOWN__BOTH_DONE:
236      abort();
237    case WAY_DOWN_UP__NOTHING_DONE:
238      way_ = WAY_DOWN_UP__DOWN_DONE;
239      return true;
240    case WAY_DOWN_UP__DOWN_DONE:
241      way_ = WAY_DOWN_UP__BOTH_DONE;
242      return false;
243    case WAY_UP_DOWN__NOTHING_DONE:
244      way_ = WAY_UP_DOWN__UP_DONE;
245      return false;
246    case WAY_UP_DOWN__UP_DONE:
247      way_ = WAY_UP_DOWN__BOTH_DONE;
248      return true;
249    }
250    return true; // to placate some compilers...
251  }
252   
253 
254  // Public data
255  // The id of the node
256  int node_id_;
257  // Basis (should use tree, but not as wasteful as bounds!)
258  CoinWarmStart * basis_;
259  // Objective value (COIN_DBL_MAX) if spare node
260  double objectiveValue_;
261  // Branching variable (0 is first integer)
262  int variable_;
263  // Way to branch: see enum DBNodeWay
264  int way_;
265  // Number of integers (for length of arrays)
266  int numberIntegers_;
267  // Current value
268  double value_;
269  // Parent
270  int parent_;
271  // Left child
272  int child_down_;
273  // Right child
274  int child_up_;
275  // Whether strong branching fixed variables when we branched on this node
276  bool strong_branching_fixed_vars_;
277  // Whether reduced cost fixing fixed anything in this node
278  bool reduced_cost_fixed_vars_;
279  // Previous in chain
280  int previous_;
281  // Next in chain
282  int next_;
283  // Now I must use tree
284  // Bounds stored in full (for integers)
285  int * lower_;
286  int * upper_;
287};
288
289
290DBNodeSimple::DBNodeSimple() :
291  node_id_(-1),
292  basis_(NULL),
293  objectiveValue_(COIN_DBL_MAX),
294  variable_(-100),
295  way_(WAY_UNSET),
296  numberIntegers_(0),
297  value_(0.5),
298  parent_(-1),
299  child_down_(-1),
300  child_up_(-1),
301  strong_branching_fixed_vars_(false),
302  reduced_cost_fixed_vars_(false),
303  previous_(-1),
304  next_(-1),
305  lower_(NULL),
306  upper_(NULL)
307{
308}
309DBNodeSimple::DBNodeSimple(OsiSolverInterface & model,
310                           int numberIntegers, int * integer,CoinWarmStart * basis)
311{
312  gutsOfConstructor(model,numberIntegers,integer,basis);
313}
314void
315DBNodeSimple::gutsOfConstructor(OsiSolverInterface & model,
316                                int numberIntegers, int * integer,CoinWarmStart * basis)
317{
318  node_id_ = -1;
319  basis_ = basis;
320  variable_=-1;
321  way_ = WAY_UNSET;
322  numberIntegers_=numberIntegers;
323  value_=0.0;
324  parent_ = -1;
325  child_down_ = -1;
326  child_up_ = -1;
327  strong_branching_fixed_vars_ = false;
328  reduced_cost_fixed_vars_ = false;
329  previous_ = -1;
330  next_ = -1;
331  if (model.isProvenOptimal()&&!model.isDualObjectiveLimitReached()) {
332    objectiveValue_ = model.getObjSense()*model.getObjValue();
333  } else {
334    objectiveValue_ = 1.0e100;
335    lower_ = NULL;
336    upper_ = NULL;
337    return; // node cutoff
338  }
339  lower_ = new int [numberIntegers_];
340  upper_ = new int [numberIntegers_];
341  assert (upper_!=NULL);
342  const double * lower = model.getColLower();
343  const double * upper = model.getColUpper();
344  const double * solution = model.getColSolution();
345  int i;
346  // Hard coded integer tolerance
347#define INTEGER_TOLERANCE 1.0e-6
348  ///////// Start of Strong branching code - can be ignored
349  // Number of strong branching candidates
350#define STRONG_BRANCHING 5
351#ifdef STRONG_BRANCHING
352  double upMovement[STRONG_BRANCHING];
353  double downMovement[STRONG_BRANCHING];
354  double solutionValue[STRONG_BRANCHING];
355  int chosen[STRONG_BRANCHING];
356  int iSmallest=0;
357  // initialize distance from integer
358  for (i=0;i<STRONG_BRANCHING;i++) {
359    upMovement[i]=0.0;
360    chosen[i]=-1;
361  }
362  variable_=-1;
363  // This has hard coded integer tolerance
364  double mostAway=INTEGER_TOLERANCE;
365  int numberAway=0;
366  for (i=0;i<numberIntegers;i++) {
367    int iColumn = integer[i];
368    lower_[i]=(int)lower[iColumn];
369    upper_[i]=(int)upper[iColumn];
370    double value = solution[iColumn];
371    value = max(value,(double) lower_[i]);
372    value = min(value,(double) upper_[i]);
373    double nearest = floor(value+0.5);
374    if (fabs(value-nearest)>INTEGER_TOLERANCE)
375      numberAway++;
376    if (fabs(value-nearest)>mostAway) {
377      double away = fabs(value-nearest);
378      if (away>upMovement[iSmallest]) {
379        //add to list
380        upMovement[iSmallest]=away;
381        solutionValue[iSmallest]=value;
382        chosen[iSmallest]=i;
383        int j;
384        iSmallest=-1;
385        double smallest = 1.0;
386        for (j=0;j<STRONG_BRANCHING;j++) {
387          if (upMovement[j]<smallest) {
388            smallest=upMovement[j];
389            iSmallest=j;
390          }
391        }
392      }
393    }
394  }
395  int numberStrong=0;
396  for (i=0;i<STRONG_BRANCHING;i++) {
397    if (chosen[i]>=0) { 
398      numberStrong ++;
399      variable_ = chosen[i];
400    }
401  }
402  // out strong branching if bit set
403  OsiClpSolverInterface* clp =
404    dynamic_cast<OsiClpSolverInterface*>(&model);
405  if (clp&&(clp->specialOptions()&16)!=0&&numberStrong>1) {
406    int j;
407    int iBest=-1;
408    double best = 0.0;
409    for (j=0;j<STRONG_BRANCHING;j++) {
410      if (upMovement[j]>best) {
411        best=upMovement[j];
412        iBest=j;
413      }
414    }
415    numberStrong=1;
416    variable_=chosen[iBest];
417  }
418  if (numberStrong==1) {
419    // just one - makes it easy
420    int iColumn = integer[variable_];
421    double value = solution[iColumn];
422    value = max(value,(double) lower_[variable_]);
423    value = min(value,(double) upper_[variable_]);
424    double nearest = floor(value+0.5);
425    value_=value;
426    if (value<=nearest)
427      way_ = WAY_UP_DOWN__NOTHING_DONE; // up
428    else
429      way_ = WAY_DOWN_UP__NOTHING_DONE; // down
430  } else if (numberStrong) {
431    // more than one - choose
432    bool chooseOne=true;
433    model.markHotStart();
434    for (i=0;i<STRONG_BRANCHING;i++) {
435      int iInt = chosen[i];
436      if (iInt>=0) {
437        int iColumn = integer[iInt];
438        double value = solutionValue[i]; // value of variable in original
439        double objectiveChange;
440        value = max(value,(double) lower_[iInt]);
441        value = min(value,(double) upper_[iInt]);
442
443        // try down
444
445        model.setColUpper(iColumn,floor(value));
446        model.solveFromHotStart();
447        model.setColUpper(iColumn,upper_[iInt]);
448        if (model.isProvenOptimal()&&!model.isDualObjectiveLimitReached()) {
449          objectiveChange = model.getObjSense()*model.getObjValue()
450            - objectiveValue_;
451        } else {
452          objectiveChange = 1.0e100;
453        }
454        assert (objectiveChange>-1.0e-5);
455        objectiveChange = CoinMax(objectiveChange,0.0);
456        downMovement[i]=objectiveChange;
457
458        // try up
459
460        model.setColLower(iColumn,ceil(value));
461        model.solveFromHotStart();
462        model.setColLower(iColumn,lower_[iInt]);
463        if (model.isProvenOptimal()&&!model.isDualObjectiveLimitReached()) {
464          objectiveChange = model.getObjSense()*model.getObjValue()
465            - objectiveValue_;
466        } else {
467          objectiveChange = 1.0e100;
468        }
469        assert (objectiveChange>-1.0e-5);
470        objectiveChange = CoinMax(objectiveChange,0.0);
471        upMovement[i]=objectiveChange;
472       
473        /* Possibilities are:
474           Both sides feasible - store
475           Neither side feasible - set objective high and exit
476           One side feasible - change bounds and resolve
477        */
478        bool solveAgain=false;
479        if (upMovement[i]<1.0e100) {
480          if(downMovement[i]<1.0e100) {
481            // feasible - no action
482          } else {
483            // up feasible, down infeasible
484            solveAgain = true;
485            model.setColLower(iColumn,ceil(value));
486          }
487        } else {
488          if(downMovement[i]<1.0e100) {
489            // down feasible, up infeasible
490            solveAgain = true;
491            model.setColUpper(iColumn,floor(value));
492          } else {
493            // neither side feasible
494            objectiveValue_=1.0e100;
495            chooseOne=false;
496            break;
497          }
498        }
499        if (solveAgain) {
500          // need to solve problem again - signal this
501          variable_ = numberIntegers;
502          chooseOne=false;
503          break;
504        }
505      }
506    }
507    if (chooseOne) {
508      // choose the one that makes most difference both ways
509      double best = -1.0;
510      double best2 = -1.0;
511      for (i=0;i<STRONG_BRANCHING;i++) {
512        int iInt = chosen[i];
513        if (iInt>=0) {
514          //std::cout<<"Strong branching on "
515          //   <<i<<""<<iInt<<" down "<<downMovement[i]
516          //   <<" up "<<upMovement[i]
517          //   <<" value "<<solutionValue[i]
518          //   <<std::endl;
519          bool better = false;
520          if (min(upMovement[i],downMovement[i])>best) {
521            // smaller is better
522            better=true;
523          } else if (min(upMovement[i],downMovement[i])>best-1.0e-5) {
524            if (max(upMovement[i],downMovement[i])>best2+1.0e-5) {
525              // smaller is about same, but larger is better
526              better=true;
527            }
528          }
529          if (better) {
530            best = min(upMovement[i],downMovement[i]);
531            best2 = max(upMovement[i],downMovement[i]);
532            variable_ = iInt;
533            double value = solutionValue[i];
534            value = max(value,(double) lower_[variable_]);
535            value = min(value,(double) upper_[variable_]);
536            value_=value;
537            if (upMovement[i]<=downMovement[i])
538              way_ = WAY_UP_DOWN__NOTHING_DONE; // up
539            else
540              way_ = WAY_DOWN_UP__NOTHING_DONE; // down
541          }
542        }
543      }
544    }
545    // Delete the snapshot
546    model.unmarkHotStart();
547  }
548  ////// End of Strong branching
549#else
550  variable_=-1;
551  // This has hard coded integer tolerance
552  double mostAway=INTEGER_TOLERANCE;
553  int numberAway=0;
554  for (i=0;i<numberIntegers;i++) {
555    int iColumn = integer[i];
556    lower_[i]=(int)lower[iColumn];
557    upper_[i]=(int)upper[iColumn];
558    double value = solution[iColumn];
559    value = max(value,(double) lower_[i]);
560    value = min(value,(double) upper_[i]);
561    double nearest = floor(value+0.5);
562    if (fabs(value-nearest)>INTEGER_TOLERANCE)
563      numberAway++;
564    if (fabs(value-nearest)>mostAway) {
565      mostAway=fabs(value-nearest);
566      variable_=i;
567      value_=value;
568      if (value<=nearest)
569        way_ = WAY_UP_DOWN__NOTHING_DONE; // up
570      else
571        way_ = WAY_DOWN_UP__NOTHING_DONE; // down
572    }
573  }
574#endif
575}
576
577void
578DBNodeSimple::gutsOfCopy(const DBNodeSimple & rhs)
579{
580  node_id_=rhs.node_id_;
581  basis_ = rhs.basis_ ? rhs.basis_->clone() : NULL;
582  objectiveValue_=rhs.objectiveValue_;
583  variable_=rhs.variable_;
584  way_=rhs.way_;
585  numberIntegers_=rhs.numberIntegers_;
586  value_=rhs.value_;
587  parent_ = rhs.parent_;
588  child_down_ = rhs.child_down_;
589  child_up_ = rhs.child_up_;
590  strong_branching_fixed_vars_ = rhs.strong_branching_fixed_vars_;
591  reduced_cost_fixed_vars_ = rhs.reduced_cost_fixed_vars_;
592  previous_ = rhs.previous_;
593  next_ = rhs.next_;
594  lower_=NULL;
595  upper_=NULL;
596  if (rhs.lower_!=NULL) {
597    lower_ = new int [numberIntegers_];
598    upper_ = new int [numberIntegers_];
599    assert (upper_!=NULL);
600    memcpy(lower_,rhs.lower_,numberIntegers_*sizeof(int));
601    memcpy(upper_,rhs.upper_,numberIntegers_*sizeof(int));
602  }
603}
604
605DBNodeSimple::DBNodeSimple(const DBNodeSimple & rhs)
606{
607  gutsOfCopy(rhs);
608}
609
610DBNodeSimple &
611DBNodeSimple::operator=(const DBNodeSimple & rhs)
612{
613  if (this != &rhs) {
614    // LL: in original code basis/lower/upper was left alone if rhs did not
615    // have them. Was that intentional?
616    gutsOfDestructor();
617    gutsOfCopy(rhs);
618  }
619  return *this;
620}
621
622
623DBNodeSimple::~DBNodeSimple ()
624{
625  gutsOfDestructor();
626}
627// Work of destructor
628void 
629DBNodeSimple::gutsOfDestructor()
630{
631  delete [] lower_;
632  delete [] upper_;
633  delete basis_;
634  lower_ = NULL;
635  upper_ = NULL;
636  basis_ = NULL;
637  objectiveValue_ = COIN_DBL_MAX;
638}
639// Extension - true if other extension of this
640bool 
641DBNodeSimple::extension(const DBNodeSimple & other,
642                        const double * originalLower,
643                        const double * originalUpper) const
644{
645  bool ok=true;
646  for (int i=0;i<numberIntegers_;i++) {
647    if (upper_[i]<originalUpper[i]||
648        lower_[i]>originalLower[i]) {
649      if (other.upper_[i]>upper_[i]||
650          other.lower_[i]<lower_[i]) {
651        ok=false;
652        break;
653      }
654    }
655  }
656  return ok;
657}
658
659#include <vector>
660
661// Must code up by hand
662class DBVectorNode  {
663 
664public:
665   
666  // Default Constructor
667  DBVectorNode ();
668
669  // Copy constructor
670  DBVectorNode ( const DBVectorNode &);
671   
672  // Assignment operator
673  DBVectorNode & operator=( const DBVectorNode& rhs);
674
675  // Destructor
676  ~DBVectorNode ();
677  // Size
678  inline int size() const
679  { return size_-sizeDeferred_;}
680  // Push
681  void push_back(DBNodeSimple & node); // the node_id_ of node will change
682  /* Remove a single node from the tree and adjust the previos_/next_/first_
683     etc fields. Does NOT update child/parent relationship */
684  void removeNode(int n);
685  /* Remove the subtree rooted at node n. properly adjusts previos_/next_ etc
686     fields. Does NOT update child/parent relationships. */
687  void removeSubTree(int n);
688  // Last one in (or other criterion)
689  DBNodeSimple back() const;
690  // Get rid of last one
691  void pop_back();
692  // Works out best one
693  int best() const;
694  // Rearranges the tree
695  void moveNodeUp(const int* which,
696                  OsiSolverInterface& model, DBNodeSimple & node);
697  // Fix the bounds in the descendants of subroot
698  void adjustBounds(int subroot, int brvar, int brlb, int brub);
699
700  // Check that the bounds correspond to the branching decisions...
701  void checkTree() const;
702  void checkNode(int node) const;
703
704  // Public data
705  // Maximum size
706  int maximumSize_;
707  // Current size
708  int size_;
709  // Number still hanging around
710  int sizeDeferred_;
711  // First spare
712  int firstSpare_;
713  // First
714  int first_;
715  // Last
716  int last_;
717  // Chosen one
718  mutable int chosen_;
719  // Nodes
720  DBNodeSimple * nodes_;
721};
722
723
724DBVectorNode::DBVectorNode() :
725  maximumSize_(10),
726  size_(0),
727  sizeDeferred_(0),
728  firstSpare_(0),
729  first_(-1),
730  last_(-1)
731{
732  nodes_ = new DBNodeSimple[maximumSize_];
733  for (int i=0;i<maximumSize_;i++) {
734    nodes_[i].previous_=i-1;
735    nodes_[i].next_=i+1;
736  }
737}
738
739DBVectorNode::DBVectorNode(const DBVectorNode & rhs) 
740{ 
741  maximumSize_ = rhs.maximumSize_;
742  size_ = rhs.size_;
743  sizeDeferred_ = rhs.sizeDeferred_;
744  firstSpare_ = rhs.firstSpare_;
745  first_ = rhs.first_;
746  last_ = rhs.last_;
747  nodes_ = new DBNodeSimple[maximumSize_];
748  for (int i=0;i<maximumSize_;i++) {
749    nodes_[i] = rhs.nodes_[i];
750  }
751}
752
753DBVectorNode &
754DBVectorNode::operator=(const DBVectorNode & rhs)
755{
756  if (this != &rhs) {
757    delete [] nodes_;
758    maximumSize_ = rhs.maximumSize_;
759    size_ = rhs.size_;
760    sizeDeferred_ = rhs.sizeDeferred_;
761    firstSpare_ = rhs.firstSpare_;
762    first_ = rhs.first_;
763    last_ = rhs.last_;
764    nodes_ = new DBNodeSimple[maximumSize_];
765    for (int i=0;i<maximumSize_;i++) {
766      nodes_[i] = rhs.nodes_[i];
767    }
768  }
769  return *this;
770}
771
772
773DBVectorNode::~DBVectorNode ()
774{
775  delete [] nodes_;
776}
777// Push
778void 
779DBVectorNode::push_back(DBNodeSimple & node)
780{
781  if (size_==maximumSize_) {
782    assert (firstSpare_==size_);
783    maximumSize_ = (maximumSize_*3)+10;
784    DBNodeSimple * temp = new DBNodeSimple[maximumSize_];
785    int i;
786    for (i=0;i<size_;i++) {
787      temp[i]=nodes_[i];
788    }
789    delete [] nodes_;
790    nodes_ = temp;
791    //firstSpare_=size_;
792    int last = -1;
793    for ( i=size_;i<maximumSize_;i++) {
794      nodes_[i].previous_=last;
795      nodes_[i].next_=i+1;
796      last = i;
797    }
798  }
799  assert (firstSpare_<maximumSize_);
800  assert (nodes_[firstSpare_].previous_<0);
801  int next = nodes_[firstSpare_].next_;
802  nodes_[firstSpare_]=node;
803  nodes_[firstSpare_].node_id_ = firstSpare_;
804  if (last_>=0) {
805    assert (nodes_[last_].next_==-1);
806    nodes_[last_].next_=firstSpare_;
807  }
808  nodes_[firstSpare_].previous_=last_;
809  nodes_[firstSpare_].next_=-1;
810  if (last_==-1) {
811    assert (first_==-1);
812    first_ = firstSpare_;
813  }
814  last_=firstSpare_;
815  if (next>=0&&next<maximumSize_) {
816    firstSpare_ = next;
817    nodes_[firstSpare_].previous_=-1;
818  } else {
819    firstSpare_=maximumSize_;
820  }
821  chosen_ = -1;
822  //best();
823  size_++;
824  if (node.bothChildDone())
825    sizeDeferred_++;
826}
827// Works out best one
828int 
829DBVectorNode::best() const
830{
831  // can modify
832  chosen_=-1;
833  if (chosen_<0) {
834    chosen_=last_;
835    while (nodes_[chosen_].bothChildDone()) {
836      chosen_ = nodes_[chosen_].previous_;
837      assert (chosen_>=0);
838    }
839  }
840  return chosen_;
841}
842// Last one in (or other criterion)
843DBNodeSimple
844DBVectorNode::back() const
845{
846  assert (last_>=0);
847  return nodes_[best()];
848}
849
850void
851DBVectorNode::removeNode(int n)
852{
853  DBNodeSimple& node = nodes_[n];
854  if (node.bothChildDone())
855    sizeDeferred_--;
856  int previous = node.previous_;
857  int next = node.next_;
858  node.~DBNodeSimple();
859  if (previous >= 0) {
860    nodes_[previous].next_=next;
861  } else {
862    first_ = next;
863  }
864  if (next >= 0) {
865    nodes_[next].previous_ = previous;
866  } else {
867    last_ = previous;
868  }
869  node.previous_ = -1;
870  if (firstSpare_ >= 0) {
871    node.next_ = firstSpare_;
872  } else {
873    node.next_ = -1;
874  }
875  firstSpare_ = n;
876  assert (size_>0);
877  size_--;
878}
879
880void
881DBVectorNode::removeSubTree(int n)
882{
883  if (nodes_[n].child_down_ >= 0) {
884    removeSubTree(nodes_[n].child_down_);
885  }
886  if (nodes_[n].child_up_ >= 0) {
887    removeSubTree(nodes_[n].child_up_);
888  }
889  removeNode(n);
890}
891
892// Get rid of last one
893void 
894DBVectorNode::pop_back()
895{
896  // Temporary until more sophisticated
897  //assert (last_==chosen_);
898  removeNode(chosen_);
899  chosen_ = -1;
900}
901
902bool
903DBNodeSimple::canSwitchParentWithGrandparent(const int* which,
904                                             const LPresult& lpres,
905                                             const int * original_lower,
906                                             const int * original_upper,
907                                             DBVectorNode & branchingTree)
908{
909  /*
910    The current node ('this') is really the parent (P) and the solution in the
911    lpres represents the child. The grandparent (GP) is this.parent_. Let's have
912    variable names respecting the truth.
913  */
914#if !defined(FUNNY_BRANCHING)
915  return false;
916#endif
917
918  const int parent_id = node_id_;
919  const DBNodeSimple& parent = branchingTree.nodes_[parent_id];
920  const int grandparent_id = parent.parent_;
921
922  if (grandparent_id == -1) {
923    // can't flip root higher up...
924    return false;
925  }
926  const DBNodeSimple& grandparent = branchingTree.nodes_[grandparent_id];
927 
928  // THINK: these are not going to happen (hopefully), but if they do, what
929  // should we do???
930  if (lpres.isAbandoned) {
931    printf("WARNING: lpres.isAbandoned true!\n");
932    return false;
933  }
934  if (lpres.isProvenDualInfeasible) {
935    printf("WARNING: lpres.isProvenDualInfeasible true!\n");
936    return false;
937  }
938  if (lpres.isPrimalObjectiveLimitReached) {
939    printf("WARNING: lpres.isPrimalObjectiveLimitReached true!\n");
940    return false;
941  }
942
943  if (parent.strong_branching_fixed_vars_ || parent.reduced_cost_fixed_vars_ ||
944      grandparent.strong_branching_fixed_vars_ ||
945      grandparent.reduced_cost_fixed_vars_) {
946    return false;
947  }
948 
949  const double direction = lpres.getObjSense;
950
951  const int GP_brvar = grandparent.variable_;
952  const int GP_brvar_fullid = which[GP_brvar];
953  const bool parent_is_down_child = parent_id == grandparent.child_down_;
954
955  if (lpres.isProvenOptimal ||
956      lpres.isIterationLimitReached) {
957    // THINK: should we do anything? like:
958#if 0
959    double djValue = lpres.getReducedCost[GP_brvar_fullid]*direction;
960    if (djValue > 1.0e-6) {
961      // wants to go down
962      return (parent_is_down_child);
963    } else if (djValue < -1.0e-6) {
964      return (! parent_is_down_child);
965    }
966#endif
967    return false;
968  }
969   
970  if (lpres.isDualObjectiveLimitReached) {
971    // Dual feasible, and in this case we don't care how we have
972    // stopped (iteration limit, obj val limit, time limit, optimal solution,
973    // etc.), we can just look at the reduced costs to figure out if the
974    // grandparent is irrelevant. Remember, we are using dual simplex!
975    double djValue = lpres.getReducedCost[GP_brvar_fullid]*direction;
976    if (djValue > 1.0e-6) {
977      // wants to go down
978      if (parent_is_down_child) {
979        return true;
980      }
981      if (lpres.getColLower[GP_brvar_fullid] > std::ceil(grandparent.value_)) {
982        return true;
983      }
984    } else if (djValue < -1.0e-6) {
985      // wants to go up
986      if (! parent_is_down_child) {
987        return true;
988      }
989      if (lpres.getColUpper[GP_brvar_fullid] < std::floor(grandparent.value_)) {
990        return true;
991      }
992    }
993    return false;
994  } else {
995    assert(lpres.isProvenPrimalInfeasible);
996    return false;
997    const int greatgrandparent_id = grandparent.parent_;
998    const int x = GP_brvar_fullid; // for easier reference... we'll use s_x
999
1000    /*
1001      Now we are going to check that if we relax the GP's branching decision
1002      then the child's LP relaxation is still infeasible. The test is done by
1003      checking whether the column (or its negative) of the GP's branching
1004      variable cuts off the dual ray proving the infeasibility.
1005    */
1006   
1007    const double* dj = lpres.getReducedCost;
1008    const double* obj = lpres.getObjCoefficients;
1009    const double yA_x = dj[x] - obj[x];
1010
1011    if (direction > 0) { // minimization problem
1012      if (parent_is_down_child) {
1013        const double s_x = CoinMax(yA_x, -1e-8);
1014        if (s_x < 1e-6) {
1015          return true;
1016        }
1017        if (lpres.yb_plus_rl_minus_su < 1e-8) {
1018          printf("WARNING: lpres.yb_plus_rl_minus_su is not positive!\n");
1019          return false;
1020        }
1021        const double max_u_x = lpres.yb_plus_rl_minus_su / s_x + lpres.getColUpper[x];
1022        const double u_x_without_GP = greatgrandparent_id >= 0 ?
1023          branchingTree.nodes_[greatgrandparent_id].upper_[GP_brvar] :
1024          original_upper[GP_brvar];
1025        return max_u_x > u_x_without_GP - 1e-8;
1026      } else {
1027        const double r_x = CoinMax(yA_x, -1e-8);
1028        if (r_x < 1e-6) {
1029          return true;
1030        }
1031        if (lpres.yb_plus_rl_minus_su < 1e-8) {
1032          printf("WARNING: lpres.yb_plus_rl_minus_su is not positive!\n");
1033          return false;
1034        }
1035        const double min_l_x = - lpres.yb_plus_rl_minus_su / r_x + lpres.getColLower[x];
1036        const double l_x_without_GP = greatgrandparent_id >= 0 ?
1037          branchingTree.nodes_[greatgrandparent_id].lower_[GP_brvar] :
1038          original_lower[GP_brvar];
1039        return min_l_x < l_x_without_GP + 1e-8;
1040      }
1041    } else { // maximization problem
1042      if (parent_is_down_child) {
1043        const double s_x = CoinMin(yA_x, 1e-8);
1044        if (s_x > -1e-6) {
1045          return true;
1046        }
1047        if (lpres.yb_plus_rl_minus_su > -1e-8) {
1048          printf("WARNING: lpres.yb_plus_rl_minus_su is not negative!\n");
1049          return false;
1050        }
1051        const double max_u_x = lpres.yb_plus_rl_minus_su / s_x + lpres.getColUpper[x];
1052        const double u_x_without_GP = greatgrandparent_id >= 0 ?
1053          branchingTree.nodes_[greatgrandparent_id].upper_[GP_brvar] :
1054          original_upper[GP_brvar];
1055        return max_u_x > u_x_without_GP - 1e-8;
1056      } else {
1057        const double r_x = CoinMin(yA_x, 1e-8);
1058        if (r_x < -1e-6) {
1059          return true;
1060        }
1061        if (lpres.yb_plus_rl_minus_su > -1e-8) {
1062          printf("WARNING: lpres.yb_plus_rl_minus_su is not negative!\n");
1063          return false;
1064        }
1065        const double min_l_x = - lpres.yb_plus_rl_minus_su / r_x + lpres.getColLower[x];
1066        const double l_x_without_GP = greatgrandparent_id >= 0 ?
1067          branchingTree.nodes_[greatgrandparent_id].lower_[GP_brvar] :
1068          original_lower[GP_brvar];
1069        return min_l_x < l_x_without_GP + 1e-8;
1070      }
1071    }
1072  }
1073
1074  return true; // to placate some compilers
1075}
1076
1077void
1078DBVectorNode::adjustBounds(int subroot, int brvar, int brlb, int brub)
1079{
1080  assert(subroot != -1);
1081  DBNodeSimple& node = nodes_[subroot];
1082  // Take the intersection of brvar's bounds in node and those passed in
1083  brub = CoinMin(brub, node.upper_[brvar]);
1084  brlb = CoinMax(brlb, node.lower_[brvar]);
1085  if (brub < brlb) {
1086    // This node became infeasible. Get rid of it and of its descendants
1087    const int parent_id = node.parent_;
1088    removeSubTree(subroot);
1089    if (nodes_[parent_id].child_down_ == subroot) {
1090      nodes_[parent_id].child_down_ = -1;
1091    } else {
1092      nodes_[parent_id].child_up_ = -1;
1093    }
1094    return;
1095  }
1096  if (node.variable_ == brvar) {
1097    if (node.value_ < brlb) {
1098      // child_down_ just became infeasible. Just cut out the current node
1099      // together with its child_down_ from the tree and hang child_up_ on the
1100      // parent of this node.
1101      if (node.child_down_ >= 0) {
1102        removeSubTree(node.child_down_);
1103      }
1104      const int parent_id = node.parent_;
1105      const int child_remains = node.child_up_;
1106      if (nodes_[parent_id].child_down_ == subroot) {
1107        nodes_[parent_id].child_down_ = child_remains;
1108      } else {
1109        nodes_[parent_id].child_up_ = child_remains;
1110      }
1111      removeNode(subroot);
1112      if (child_remains >= 0) {
1113        nodes_[child_remains].parent_ = parent_id;
1114        adjustBounds(child_remains, brvar, brlb, brub);
1115      }
1116      return;
1117    }
1118    if (node.value_ > brub) {
1119      // child_up_ just became infeasible. Just cut out the current node
1120      // together with its child_down_ from the tree and hang child_down_ on
1121      // the parent of this node.
1122      if (node.child_up_ >= 0) {
1123        removeSubTree(node.child_up_);
1124      }
1125      const int parent_id = node.parent_;
1126      const int child_remains = node.child_down_;
1127      if (nodes_[parent_id].child_down_ == subroot) {
1128        nodes_[parent_id].child_down_ = child_remains;
1129      } else {
1130        nodes_[parent_id].child_up_ = child_remains;
1131      }
1132      removeNode(subroot);
1133      if (child_remains >= 0) {
1134        nodes_[child_remains].parent_ = parent_id;
1135        adjustBounds(child_remains, brvar, brlb, brub);
1136      }
1137      return;
1138    }
1139    // Now brlb < node.value_ < brub (value_ is fractional)
1140  }
1141  node.lower_[brvar] = brlb;
1142  node.upper_[brvar] = brub;
1143  if (node.child_down_ >= 0) {
1144    adjustBounds(node.child_down_, brvar, brlb, brub);
1145  }
1146  if (node.child_up_ >= 0) {
1147    adjustBounds(node.child_up_, brvar, brlb, brub);
1148  }
1149}
1150
1151void
1152DBVectorNode::moveNodeUp(const int* which,
1153                         OsiSolverInterface& model, DBNodeSimple & node)
1154{
1155  /*
1156    The current node ('this') is really the parent (P). The grandparent (GP)
1157    is this.parent_. Let's have variable names respecting the truth.
1158  */
1159  const int parent_id = node.node_id_;
1160  DBNodeSimple& parent = nodes_[parent_id];
1161  const int grandparent_id = parent.parent_;
1162  assert(grandparent_id != -1);
1163  DBNodeSimple& grandparent = nodes_[grandparent_id];
1164  const int greatgrandparent_id = grandparent.parent_;
1165 
1166  const bool parent_is_down_child = parent_id == grandparent.child_down_;
1167
1168#if defined(DEBUG_DYNAMIC_BRANCHING)
1169  if (dyn_debug >= 100) {
1170    printf("entered moveNodeUp\n");
1171    printf("parent_id %d grandparent_id %d greatgrandparent_id %d\n",
1172           parent_id, grandparent_id, greatgrandparent_id);
1173    printf("parent.way_ %d\n", parent.way_);
1174  }
1175#endif
1176
1177
1178  // First hang the nodes where they belong.
1179  parent.parent_ = greatgrandparent_id;
1180  grandparent.parent_ = parent_id;
1181  const bool down_child_stays_with_parent = parent.workingOnDownChild();
1182  int& child_to_move =
1183    down_child_stays_with_parent ? parent.child_up_ : parent.child_down_;
1184  const bool child_to_move_is_processed = parent.bothChildDone();
1185
1186#if defined(DEBUG_DYNAMIC_BRANCHING)
1187  if (dyn_debug >= 1000) {
1188    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);
1189  }
1190#endif
1191
1192  if (parent_is_down_child) {
1193    grandparent.child_down_ = child_to_move;
1194  } else {
1195    grandparent.child_up_ = child_to_move;
1196  }
1197  if (child_to_move >= 0) {
1198    nodes_[child_to_move].parent_ = grandparent_id;
1199  }
1200  child_to_move = grandparent_id;
1201  if (greatgrandparent_id >= 0) {
1202    DBNodeSimple& greatgrandparent = nodes_[greatgrandparent_id];
1203    if (grandparent_id == greatgrandparent.child_down_) {
1204      greatgrandparent.child_down_ = parent_id;
1205    } else {
1206      greatgrandparent.child_up_ = parent_id;
1207    }
1208  }
1209
1210#if defined(DEBUG_DYNAMIC_BRANCHING)
1211  if (dyn_debug >= 1000) {
1212    printf("after exchange\n");
1213    printf("parent.parent_ %d parent.child_down_ %d parent.child_up_ %d\n",
1214           parent.parent_, parent.child_down_, parent.child_up_);
1215    printf("grandparent.parent_ %d grandparent.child_down_ %d grandparent.child_up_ %d\n",
1216           grandparent.parent_, grandparent.child_down_, grandparent.child_up_);
1217    if (greatgrandparent_id >= 0) {
1218      DBNodeSimple& greatgrandparent = nodes_[greatgrandparent_id];
1219      printf("greatgrandparent.parent_ %d greatgrandparent.child_down_ %d greatgrandparent.child_up_ %d\n",
1220             greatgrandparent.parent_, greatgrandparent.child_down_, greatgrandparent.child_up_);
1221    }
1222    printf("exiting moveNodeUp\n");
1223  }
1224#endif
1225
1226
1227  // Now make sure way_ is set properly
1228  if (down_child_stays_with_parent) {
1229    if (!parent.bothChildDone()) {
1230      parent.way_ = DBNodeSimple::WAY_UP_DOWN__BOTH_DONE;
1231      sizeDeferred_++;
1232    }
1233    if (grandparent.bothChildDone()) {
1234      if (! child_to_move_is_processed) {
1235        grandparent.way_ = parent_is_down_child ?
1236          DBNodeSimple::WAY_UP_DOWN__UP_DONE :
1237          DBNodeSimple::WAY_DOWN_UP__DOWN_DONE;
1238        sizeDeferred_--;
1239      }
1240    } else { // only parent is processed from the two children of grandparent
1241      if (! child_to_move_is_processed) {
1242        // remove grandparent, none of its children is processed now, why
1243        // force its branching decision?
1244        removeNode(grandparent_id);
1245        parent.child_up_ = -1;
1246        parent.way_ = DBNodeSimple::WAY_DOWN_UP__DOWN_DONE;
1247        sizeDeferred_--;
1248      } else {
1249        grandparent.way_ = parent_is_down_child ?
1250          DBNodeSimple::WAY_DOWN_UP__DOWN_DONE :
1251          DBNodeSimple::WAY_UP_DOWN__UP_DONE;
1252      }
1253    }
1254  } else {
1255    if (!parent.bothChildDone()) {
1256      parent.way_ = DBNodeSimple::WAY_DOWN_UP__BOTH_DONE;
1257      sizeDeferred_++;
1258    }
1259    if (grandparent.bothChildDone()) {
1260      if (! child_to_move_is_processed) {
1261        grandparent.way_ = parent_is_down_child ?
1262          DBNodeSimple::WAY_UP_DOWN__UP_DONE :
1263          DBNodeSimple::WAY_DOWN_UP__DOWN_DONE;
1264        sizeDeferred_--;
1265      }
1266    } else { // only parent is processed from the two children of grandparent
1267      if (! child_to_move_is_processed) {
1268        removeNode(grandparent_id);
1269        parent.child_down_ = -1;
1270        parent.way_ = DBNodeSimple::WAY_UP_DOWN__UP_DONE;
1271        sizeDeferred_--;
1272      } else {
1273        grandparent.way_ = parent_is_down_child ?
1274          DBNodeSimple::WAY_DOWN_UP__DOWN_DONE :
1275          DBNodeSimple::WAY_UP_DOWN__UP_DONE;
1276      }
1277    }
1278  }
1279 
1280  // Now modify bounds
1281
1282  // First, get rid of GP's bound change of its branching variable in the
1283  // bound list of P. Note: If greatgrandparent_id is < 0 then GP is the root
1284  // so its bounds are the original bounds.
1285  const int GP_brvar = grandparent.variable_;
1286  if (parent_is_down_child) {
1287    const int old_upper = grandparent.upper_[GP_brvar];
1288    parent.upper_[GP_brvar] = old_upper;
1289    if ((GP_brvar != parent.variable_) ||
1290        (GP_brvar == parent.variable_ && !down_child_stays_with_parent)) {
1291      model.setColUpper(which[GP_brvar], old_upper);
1292    }
1293  } else {
1294    const int old_lower = grandparent.lower_[GP_brvar];
1295    parent.lower_[GP_brvar] = old_lower;
1296    if ((GP_brvar != parent.variable_) ||
1297        (GP_brvar == parent.variable_ && down_child_stays_with_parent)) {
1298      model.setColLower(which[GP_brvar], old_lower);
1299    }
1300  }
1301
1302  // Now add the branching var bound change of P to GP and all of its
1303  // descendant
1304  if (down_child_stays_with_parent) {
1305    adjustBounds(grandparent_id, parent.variable_,
1306                 (int)ceil(parent.value_), parent.upper_[parent.variable_]);
1307  } else {
1308    adjustBounds(grandparent_id, parent.variable_,
1309                 parent.lower_[parent.variable_], (int)floor(parent.value_));
1310  }
1311}
1312
1313void
1314DBVectorNode::checkNode(int node) const
1315{
1316  if (node == -1) {
1317    return;
1318  }
1319  const DBNodeSimple* n = nodes_ + node;
1320  const DBNodeSimple* p = nodes_ + n->parent_;
1321  for (int i = n->numberIntegers_-1; i >= 0; --i) {
1322    if (i == p->variable_) {
1323      if (node == p->child_down_) {
1324        assert(p->lower_[i] <= n->lower_[i]);
1325        assert((int)(floor(p->value_)) == n->upper_[i]);
1326      } else {
1327        assert((int)(ceil(p->value_)) == n->lower_[i]);
1328        assert(p->upper_[i] >= n->upper_[i]);
1329      }
1330    } else {
1331      assert(p->lower_[i] <= n->lower_[i]);
1332      assert(p->upper_[i] >= n->upper_[i]);
1333    }
1334  }
1335  checkNode(n->child_down_);
1336  checkNode(n->child_up_);
1337}
1338
1339void
1340DBVectorNode::checkTree() const
1341{
1342  // find the root
1343  int root = first_;
1344  while (true) {
1345    if (nodes_[root].parent_ == -1) {
1346      break;
1347    } else {
1348      root = nodes_[root].parent_;
1349    }
1350  }
1351  checkNode(nodes_[root].child_down_);
1352  checkNode(nodes_[root].child_up_);
1353}
1354
1355std::string getTree(DBVectorNode& branchingTree)
1356{
1357  std::string tree;
1358  char line[1000];
1359  for(int k=0; k<branchingTree.size_; k++) {
1360    DBNodeSimple& node = branchingTree.nodes_[k];
1361    sprintf(line, "%d %d %d %d %f %d %d 0x%x %d %d\n",
1362            k, node.node_id_, node.parent_, node.variable_,
1363            node.value_, node.lower_[node.variable_],
1364            node.upper_[node.variable_], node.way_,
1365            node.child_down_, node.child_up_);
1366    tree += line;
1367  }
1368  return tree;
1369}
1370
1371void printTree(const std::string& tree, int levels)
1372{
1373  size_t pos = tree.size();
1374  for (int i = levels-1; i >= 0; --i) {
1375    pos = tree.rfind('\n', pos-1);
1376    if (pos == std::string::npos) {
1377      pos = 0;
1378      break;
1379    }
1380  }
1381  printf("%s", tree.c_str() + (pos+1));
1382}
1383
1384void printChain(DBVectorNode& branchingTree, int k)
1385{
1386  while (k != -1) {
1387    DBNodeSimple& node = branchingTree.nodes_[k];
1388    printf("   %d %d %d %d %f %d %d 0x%x %d %d %c %c\n",
1389           k, node.node_id_, node.parent_, node.variable_,
1390           node.value_, node.lower_[node.variable_],
1391           node.upper_[node.variable_], node.way_,
1392           node.child_down_, node.child_up_,
1393           node.strong_branching_fixed_vars_ ? 'T' : 'F',
1394           node.reduced_cost_fixed_vars_ ? 'T' : 'F');
1395    k = node.parent_;
1396  }
1397}
1398
1399// Invoke solver's built-in enumeration algorithm
1400void 
1401branchAndBound(OsiSolverInterface & model) {
1402  double time1 = CoinCpuTime();
1403  // solve LP
1404  model.initialSolve();
1405
1406  if (model.isProvenOptimal()&&!model.isDualObjectiveLimitReached()) {
1407    // Continuous is feasible - find integers
1408    int numberIntegers=0;
1409    int numberColumns = model.getNumCols();
1410    int iColumn;
1411    int i;
1412    for (iColumn=0;iColumn<numberColumns;iColumn++) {
1413      if( model.isInteger(iColumn))
1414        numberIntegers++;
1415    }
1416    if (!numberIntegers) {
1417      std::cout<<"No integer variables"
1418               <<std::endl;
1419      return;
1420    }
1421    int * which = new int[numberIntegers]; // which variables are integer
1422    // original bounds
1423    int * originalLower = new int[numberIntegers];
1424    int * originalUpper = new int[numberIntegers];
1425    int * relaxedLower = new int[numberIntegers];
1426    int * relaxedUpper = new int[numberIntegers];
1427    {
1428      const double * lower = model.getColLower();
1429      const double * upper = model.getColUpper();
1430      numberIntegers=0;
1431      for (iColumn=0;iColumn<numberColumns;iColumn++) {
1432        if( model.isInteger(iColumn)) {
1433          originalLower[numberIntegers]=(int) lower[iColumn];
1434          originalUpper[numberIntegers]=(int) upper[iColumn];
1435          which[numberIntegers++]=iColumn;
1436        }
1437      }
1438    }
1439    double direction = model.getObjSense();
1440    // empty tree
1441    DBVectorNode branchingTree;
1442   
1443    // Add continuous to it;
1444    DBNodeSimple rootNode(model,numberIntegers,which,model.getWarmStart());
1445    // something extra may have been fixed by strong branching
1446    // if so go round again
1447    while (rootNode.variable_==numberIntegers) {
1448      model.resolve();
1449      rootNode = DBNodeSimple(model,numberIntegers,which,model.getWarmStart());
1450    }
1451    if (rootNode.objectiveValue_<1.0e100) {
1452      // push on stack
1453      branchingTree.push_back(rootNode);
1454    }
1455   
1456    // For printing totals
1457    int numberIterations=0;
1458    int numberNodes =0;
1459    DBNodeSimple bestNode;
1460    ////// Start main while of branch and bound
1461    // while until nothing on stack
1462    while (branchingTree.size()) {
1463#if defined(DEBUG_DYNAMIC_BRANCHING)
1464      if (dyn_debug >= 1000) {
1465        printf("branchingTree.size = %d %d\n",branchingTree.size(),branchingTree.size_);
1466        printf("i node_id parent child_down child_up\n");
1467        for(int k=0; k<branchingTree.size_; k++) {
1468          DBNodeSimple& node = branchingTree.nodes_[k];
1469          printf("%d %d %d %d %d\n",k, node.node_id_, node.parent_,
1470                 node.child_down_, node.child_up_);
1471        }
1472      }
1473#endif
1474      // last node
1475      DBNodeSimple node = branchingTree.back();
1476      int kNode = branchingTree.chosen_;
1477      branchingTree.pop_back();
1478#if defined(DEBUG_DYNAMIC_BRANCHING)
1479      if (dyn_debug >= 1000) {
1480        printf("Deleted current parent %d %d\n",branchingTree.size(),branchingTree.size_);
1481      }
1482#endif
1483      assert (! node.bothChildDone());
1484      numberNodes++;
1485      if (node.variable_ < 0) {
1486        // put back on tree and pretend both children are done. We want the
1487        // whole tree all the time.
1488        node.way_ = DBNodeSimple::WAY_UP_DOWN__BOTH_DONE;
1489        branchingTree.push_back(node);
1490        // Integer solution - save
1491        bestNode = node;
1492        // set cutoff (hard coded tolerance)
1493        const double limit = (bestNode.objectiveValue_-1.0e-5)*direction;
1494        model.setDblParam(OsiDualObjectiveLimit, limit);
1495        std::cout<<"Integer solution of "
1496                 <<bestNode.objectiveValue_
1497                 <<" found after "<<numberIterations
1498                 <<" iterations and "<<numberNodes<<" nodes"
1499                 <<std::endl;
1500#if defined(DEBUG_DYNAMIC_BRANCHING)
1501        if (dyn_debug >= 1) {
1502          branchingTree.checkTree();
1503        }
1504#endif
1505
1506      } else {
1507        // branch - do bounds
1508        for (i=0;i<numberIntegers;i++) {
1509          iColumn=which[i];
1510          model.setColBounds( iColumn,node.lower_[i],node.upper_[i]);
1511        }
1512        // move basis
1513        model.setWarmStart(node.basis_);
1514        // do branching variable
1515        const bool down_branch = node.advanceWay();
1516        if (down_branch) {
1517          model.setColUpper(which[node.variable_],floor(node.value_));
1518        } else {
1519          model.setColLower(which[node.variable_],ceil(node.value_));
1520        }
1521        // put back on tree anyway regardless whether any processing is left
1522        // to be done. We want the whole tree all the time.
1523        branchingTree.push_back(node);
1524#if defined(DEBUG_DYNAMIC_BRANCHING)
1525        if (dyn_debug >= 1000) {
1526          printf("Added current parent %d %d\n",branchingTree.size(),branchingTree.size_);
1527        }
1528#endif
1529       
1530        // solve
1531        model.resolve();
1532        CoinWarmStart * ws = model.getWarmStart();
1533        const CoinWarmStartBasis* wsb =
1534          dynamic_cast<const CoinWarmStartBasis*>(ws);
1535        assert (wsb!=NULL); // make sure not volume
1536        numberIterations += model.getIterationCount();
1537        // fix on reduced costs
1538        int nFixed0=0,nFixed1=0;
1539        double cutoff;
1540        model.getDblParam(OsiDualObjectiveLimit,cutoff);
1541        double gap=(cutoff-model.getObjValue())*direction+1.0e-4;
1542        //        double
1543        //        gap=(cutoff-modelPtr_->objectiveValue())*direction+1.0e-4;
1544        bool did_reduced_cost_fixing_for_child = false;
1545        if (gap<1.0e10&&model.isProvenOptimal()&&!model.isDualObjectiveLimitReached()) {
1546          const double * dj = model.getReducedCost();
1547          const double * lower = model.getColLower();
1548          const double * upper = model.getColUpper();
1549          for (i=0;i<numberIntegers;i++) {
1550            iColumn=which[i];
1551            if (upper[iColumn]>lower[iColumn]) {
1552              double djValue = dj[iColumn]*direction;
1553              if (wsb->getStructStatus(iColumn)==CoinWarmStartBasis::atLowerBound&&
1554                  djValue>gap) {
1555                nFixed0++;
1556                model.setColUpper(iColumn,lower[iColumn]);
1557              } else if (wsb->getStructStatus(iColumn)==CoinWarmStartBasis::atUpperBound&&
1558                         -djValue>gap) {
1559                nFixed1++;
1560                model.setColLower(iColumn,upper[iColumn]);
1561              }
1562            }
1563          }
1564          if (nFixed0+nFixed1) {
1565            //      printf("%d fixed to lower, %d fixed to upper\n",
1566            //             nFixed0,nFixed1);
1567            did_reduced_cost_fixing_for_child = true;
1568          }
1569        }
1570        if (model.isAbandoned()) {
1571          // THINK: What the heck should we do???
1572          abort();
1573        }
1574        if (model.isIterationLimitReached()) {
1575          // maximum iterations - exit
1576          std::cout<<"Exiting on maximum iterations\n";
1577          break;
1578        }
1579
1580        LPresult lpres(model);
1581
1582        bool canSwitch = node.canSwitchParentWithGrandparent(which, lpres,
1583                                                             originalLower,
1584                                                             originalUpper,
1585                                                             branchingTree);
1586        int cnt = 0;
1587
1588#if defined(DEBUG_DYNAMIC_BRANCHING)
1589        if (dyn_debug >= 1) {
1590          branchingTree.checkTree();
1591        }
1592#endif
1593
1594#if defined(DEBUG_DYNAMIC_BRANCHING)
1595        std::string tree0;
1596        if (dyn_debug >= 10) {
1597          if (canSwitch) {
1598            tree0 = getTree(branchingTree);
1599          }
1600        }
1601#endif
1602
1603        while (canSwitch) {
1604          branchingTree.moveNodeUp(which, model, node);
1605          canSwitch = node.canSwitchParentWithGrandparent(which, lpres,
1606                                                          originalLower,
1607                                                          originalUpper,
1608                                                          branchingTree);
1609#if defined(DEBUG_DYNAMIC_BRANCHING)
1610          if (dyn_debug >= 1) {
1611            branchingTree.checkTree();
1612          }
1613#endif
1614          ++cnt;
1615        }
1616        if (cnt > 0) {
1617          model.resolve();
1618          // This is horribly looking... Get rid of it when properly debugged...
1619          assert(lpres.isAbandoned == model.isAbandoned());
1620          assert(lpres.isDualObjectiveLimitReached == model.isDualObjectiveLimitReached());
1621          assert(lpres.isDualObjectiveLimitReached ||
1622                 (lpres.isProvenOptimal == model.isProvenOptimal()));
1623          assert(lpres.isDualObjectiveLimitReached ||
1624                 (lpres.isProvenPrimalInfeasible == model.isProvenPrimalInfeasible()));
1625          assert(!lpres.isProvenOptimal || ! model.isProvenOptimal() ||
1626                 (lpres.isProvenOptimal && model.isProvenOptimal() &&
1627                  lpres.getObjValue == model.getObjValue()));
1628          printf("Finished moving node %d up by %i levels.\n", node.node_id_, cnt);
1629        }
1630           
1631#if defined(DEBUG_DYNAMIC_BRANCHING)
1632        if (dyn_debug >= 10) {
1633          if (cnt > 0) {
1634            std::string tree1 = getTree(branchingTree);
1635            printf("=====================================\n");
1636            printf("It can move node %i up. way_: 0x%x   brvar: %i\n",
1637                   node.node_id_, node.way_, node.variable_);
1638            printTree(tree0, cnt+10);
1639            printf("=====================================\n");
1640            printf("Finished moving the node up by %i levels.\n", cnt);
1641            printTree(tree1, cnt+10);
1642            printf("=====================================\n");
1643          }
1644        }
1645#endif
1646        if ((numberNodes%10)==0) 
1647          printf("%d nodes, tree size %d\n",
1648                 numberNodes,branchingTree.size());
1649        if (CoinCpuTime()-time1>3600.0) {
1650          printf("stopping after 3600 seconds\n");
1651          exit(77);
1652        }
1653        DBNodeSimple newNode(model,numberIntegers,which,ws);
1654        // something extra may have been fixed by strong branching
1655        // if so go round again
1656        while (newNode.variable_==numberIntegers) {
1657          model.resolve();
1658          newNode = DBNodeSimple(model,numberIntegers,which,model.getWarmStart());
1659          newNode.strong_branching_fixed_vars_ = true;
1660        }
1661        newNode.reduced_cost_fixed_vars_ = did_reduced_cost_fixing_for_child;
1662        if (newNode.objectiveValue_<1.0e100) {
1663          newNode.parent_ = kNode;
1664          // push on stack
1665          branchingTree.push_back(newNode);
1666#if defined(DEBUG_DYNAMIC_BRANCHING)
1667          if (dyn_debug >= 1000) {
1668            printf("Added current child %d %d\n",branchingTree.size(),branchingTree.size_);
1669          }
1670#endif
1671          if (branchingTree.nodes_[kNode].workingOnDownChild()) {
1672            branchingTree.nodes_[kNode].child_down_ = branchingTree.last_;
1673          } else {
1674            branchingTree.nodes_[kNode].child_up_ = branchingTree.last_;
1675          }
1676          if (newNode.variable_>=0) {
1677            assert (fabs(newNode.value_-floor(newNode.value_+0.5))>1.0e-6);
1678          }
1679#if 0
1680          else {
1681            // integer solution - save
1682            bestNode = node;
1683            // set cutoff (hard coded tolerance)
1684            model.setDblParam(OsiDualObjectiveLimit,(bestNode.objectiveValue_-1.0e-5)*direction);
1685            std::cout<<"Integer solution of "
1686                     <<bestNode.objectiveValue_
1687                     <<" found after "<<numberIterations
1688                     <<" iterations and "<<numberNodes<<" nodes"
1689                     <<std::endl;
1690          }
1691#endif
1692        }
1693#if defined(DEBUG_DYNAMIC_BRANCHING)
1694        if (dyn_debug >= 1) {
1695          branchingTree.checkTree();
1696        }
1697#endif
1698      }
1699    }
1700    ////// End main while of branch and bound
1701    std::cout<<"Search took "
1702             <<numberIterations
1703             <<" iterations and "<<numberNodes<<" nodes"
1704             <<std::endl;
1705    if (bestNode.numberIntegers_) {
1706      // we have a solution restore
1707      // do bounds
1708      for (i=0;i<numberIntegers;i++) {
1709        iColumn=which[i];
1710        model.setColBounds( iColumn,bestNode.lower_[i],bestNode.upper_[i]);
1711      }
1712      // move basis
1713      model.setWarmStart(bestNode.basis_);
1714      // set cutoff so will be good (hard coded tolerance)
1715      model.setDblParam(OsiDualObjectiveLimit,(bestNode.objectiveValue_+1.0e-5)*direction);
1716      model.resolve();
1717    } else {
1718      OsiClpSolverInterface* clp =
1719        dynamic_cast<OsiClpSolverInterface*>(&model);
1720      if (clp) {
1721        ClpSimplex* modelPtr_ = clp->getModelPtr();
1722        modelPtr_->setProblemStatus(1);
1723      }
1724    }
1725    delete [] which;
1726    delete [] originalLower;
1727    delete [] originalUpper;
1728    delete [] relaxedLower;
1729    delete [] relaxedUpper;
1730  } else {
1731    std::cout<<"The LP relaxation is infeasible"
1732             <<std::endl;
1733    OsiClpSolverInterface* clp =
1734      dynamic_cast<OsiClpSolverInterface*>(&model);
1735    if (clp) {
1736      ClpSimplex* modelPtr_ = clp->getModelPtr();
1737      modelPtr_->setProblemStatus(1);
1738    }
1739    //throw CoinError("The LP relaxation is infeasible or too expensive",
1740    //"branchAndBound", "OsiClpSolverInterface");
1741  }
1742}
1743
1744
1745
1746int main(int argc, char* argv[])
1747{
1748  OsiClpSolverInterface model;
1749  model.readMps(argv[1]);
1750  branchAndBound(model);
1751  return 0;
1752}
Note: See TracBrowser for help on using the repository browser.