source: branches/dynamicbranching/dynamicbranching.cpp @ 999

Last change on this file since 999 was 999, checked in by jpgoncal, 13 years ago

fixed getTree()

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