source: branches/dynamicbranching/dynamicbranching.cpp @ 997

Last change on this file since 997 was 997, checked in by jpgoncal, 11 years ago

removed printout from clp.

File size: 53.2 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* 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    return false;
997    if (lpres.getObjValue > grandparent.objectiveValue_ + 1e-8) {
998      double djValue = lpres.getReducedCost[GP_brvar_fullid]*direction;
999      if (djValue > 1.0e-6) {
1000        // wants to go down
1001        if (parent_is_down_child) {
1002          return true;
1003        }
1004        if (lpres.getColLower[GP_brvar_fullid] > std::ceil(grandparent.value_)) {
1005          return true;
1006        }
1007      } else if (djValue < -1.0e-6) {
1008        // wants to go up
1009        if (! parent_is_down_child) {
1010          return true;
1011        }
1012        if (lpres.getColUpper[GP_brvar_fullid] < std::floor(grandparent.value_)) {
1013          return true;
1014        }
1015      }
1016    }
1017    return false;
1018  }
1019
1020  // Problem is really infeasible and has a dual ray.
1021  assert(lpres.isProvenPrimalInfeasible);
1022  const int greatgrandparent_id = grandparent.parent_;
1023  const int x = GP_brvar_fullid; // for easier reference... we'll use s_x
1024
1025  /*
1026    Now we are going to check that if we relax the GP's branching decision
1027    then the child's LP relaxation is still infeasible. The test is done by
1028    checking whether the column (or its negative) of the GP's branching
1029    variable cuts off the dual ray proving the infeasibility.
1030  */
1031   
1032  const double* colDualRay = lpres.colDualRay;
1033  const double ub_without_GP = greatgrandparent_id >= 0 ?
1034    branchingTree.nodes_[greatgrandparent_id].upper_[GP_brvar] :
1035    original_upper[GP_brvar];
1036  const double lb_without_GP = greatgrandparent_id >= 0 ?
1037    branchingTree.nodes_[greatgrandparent_id].lower_[GP_brvar] :
1038    original_lower[GP_brvar];
1039
1040  if (parent_is_down_child) { // upper bound on x is forced
1041    if (colDualRay[x] < 1e-8) { // if yA_x] is non-positive then the var is
1042                                // at its lower bound (or has 0 red cost)
1043                                // so same dual ray still proves infeas
1044      return true;
1045    }
1046    const double max_ub_increase = lpres.yb_plus_rl_minus_su / colDualRay[x];
1047    assert(max_ub_increase >= 1e-8);
1048    const bool willSwitch =
1049      max_ub_increase > ub_without_GP - lpres.getColUpper[x] + 1e-8;
1050    if (willSwitch) {
1051      // The switch will happen so we might as well update the "infeasibility"
1052      lpres.yb_plus_rl_minus_su -= colDualRay[x] * (ub_without_GP - lpres.getColUpper[x]);
1053    }
1054    return willSwitch;
1055  } else { // lower bound on x is forced
1056    if (colDualRay[x] > -1e-8) { // if yA_x is non-negative then the var is
1057                                 // at its upper bound (or has 0 red cost)
1058                                 // so same dual ray still proves infeas
1059      return true;
1060    }
1061    const double max_lb_decrease = - lpres.yb_plus_rl_minus_su / colDualRay[x];
1062    assert(max_lb_decrease >= 1e-8);
1063    const bool willSwitch =
1064      max_lb_decrease > lpres.getColLower[x] - lb_without_GP + 1e-8;
1065    if (willSwitch) {
1066      // The switch will happen so we might as well update the "infeasibility"
1067      lpres.yb_plus_rl_minus_su += colDualRay[x] * (lpres.getColLower[x] - lb_without_GP);
1068    }
1069    return willSwitch;
1070  }
1071  // THINK: the above is definitely good for minimization problems. Is it good
1072  // for max?
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        return; // No bound changes are needed on the GO side as GP is
1249                // deleted...
1250      } else {
1251        grandparent.way_ = parent_is_down_child ?
1252          DBNodeSimple::WAY_DOWN_UP__DOWN_DONE :
1253          DBNodeSimple::WAY_UP_DOWN__UP_DONE;
1254      }
1255    }
1256  } else {
1257    if (!parent.bothChildDone()) {
1258      parent.way_ = DBNodeSimple::WAY_DOWN_UP__BOTH_DONE;
1259      sizeDeferred_++;
1260    }
1261    if (grandparent.bothChildDone()) {
1262      if (! child_to_move_is_processed) {
1263        grandparent.way_ = parent_is_down_child ?
1264          DBNodeSimple::WAY_UP_DOWN__UP_DONE :
1265          DBNodeSimple::WAY_DOWN_UP__DOWN_DONE;
1266        sizeDeferred_--;
1267      }
1268    } else { // only parent is processed from the two children of grandparent
1269      if (! child_to_move_is_processed) {
1270        removeNode(grandparent_id);
1271        parent.child_down_ = -1;
1272        parent.way_ = DBNodeSimple::WAY_UP_DOWN__UP_DONE;
1273        sizeDeferred_--;
1274        return;
1275      } else {
1276        grandparent.way_ = parent_is_down_child ?
1277          DBNodeSimple::WAY_DOWN_UP__DOWN_DONE :
1278          DBNodeSimple::WAY_UP_DOWN__UP_DONE;
1279      }
1280    }
1281  }
1282 
1283  // Now modify bounds
1284
1285  // First, get rid of GP's bound change of its branching variable in the
1286  // bound list of P. Note: If greatgrandparent_id is < 0 then GP is the root
1287  // so its bounds are the original bounds.
1288  const int GP_brvar = grandparent.variable_;
1289  if (parent_is_down_child) {
1290    const int old_upper = grandparent.upper_[GP_brvar];
1291    parent.upper_[GP_brvar] = old_upper;
1292    if ((GP_brvar != parent.variable_) ||
1293        (GP_brvar == parent.variable_ && !down_child_stays_with_parent)) {
1294      model.setColUpper(which[GP_brvar], old_upper);
1295    }
1296  } else {
1297    const int old_lower = grandparent.lower_[GP_brvar];
1298    parent.lower_[GP_brvar] = old_lower;
1299    if ((GP_brvar != parent.variable_) ||
1300        (GP_brvar == parent.variable_ && down_child_stays_with_parent)) {
1301      model.setColLower(which[GP_brvar], old_lower);
1302    }
1303  }
1304
1305  // Now add the branching var bound change of P to GP and all of its
1306  // descendant
1307  if (down_child_stays_with_parent) {
1308    adjustBounds(grandparent_id, parent.variable_,
1309                 (int)ceil(parent.value_), parent.upper_[parent.variable_]);
1310  } else {
1311    adjustBounds(grandparent_id, parent.variable_,
1312                 parent.lower_[parent.variable_], (int)floor(parent.value_));
1313  }
1314}
1315
1316void
1317DBVectorNode::checkNode(int node) const
1318{
1319  if (node == -1) {
1320    return;
1321  }
1322  const DBNodeSimple* n = nodes_ + node;
1323  const DBNodeSimple* p = nodes_ + n->parent_;
1324  for (int i = n->numberIntegers_-1; i >= 0; --i) {
1325    if (i == p->variable_) {
1326      if (node == p->child_down_) {
1327        assert(p->lower_[i] <= n->lower_[i]);
1328        assert((int)(floor(p->value_)) == n->upper_[i]);
1329      } else {
1330        assert((int)(ceil(p->value_)) == n->lower_[i]);
1331        assert(p->upper_[i] >= n->upper_[i]);
1332      }
1333    } else {
1334      assert(p->lower_[i] <= n->lower_[i]);
1335      assert(p->upper_[i] >= n->upper_[i]);
1336    }
1337  }
1338  checkNode(n->child_down_);
1339  checkNode(n->child_up_);
1340}
1341
1342void
1343DBVectorNode::checkTree() const
1344{
1345  // find the root
1346  int root = first_;
1347  while (true) {
1348    if (nodes_[root].parent_ == -1) {
1349      break;
1350    } else {
1351      root = nodes_[root].parent_;
1352    }
1353  }
1354  checkNode(nodes_[root].child_down_);
1355  checkNode(nodes_[root].child_up_);
1356}
1357
1358std::string getTree(DBVectorNode& branchingTree)
1359{
1360  std::string tree;
1361  char line[1000];
1362  for(int k=0; k<branchingTree.size_; k++) {
1363    DBNodeSimple& node = branchingTree.nodes_[k];
1364    sprintf(line, "%d %d %d %d %f %d %d 0x%x %d %d\n",
1365            k, node.node_id_, node.parent_, node.variable_,
1366            node.value_, node.lower_[node.variable_],
1367            node.upper_[node.variable_], node.way_,
1368            node.child_down_, node.child_up_);
1369    tree += line;
1370  }
1371  return tree;
1372}
1373
1374void printTree(const std::string& tree, int levels)
1375{
1376  size_t pos = tree.size();
1377  for (int i = levels-1; i >= 0; --i) {
1378    pos = tree.rfind('\n', pos-1);
1379    if (pos == std::string::npos) {
1380      pos = 0;
1381      break;
1382    }
1383  }
1384  printf("%s", tree.c_str() + (pos+1));
1385}
1386
1387void printChain(DBVectorNode& branchingTree, int k)
1388{
1389  while (k != -1) {
1390    DBNodeSimple& node = branchingTree.nodes_[k];
1391    printf("   %d %d %d %d %f %d %d 0x%x %d %d %c %c\n",
1392           k, node.node_id_, node.parent_, node.variable_,
1393           node.value_, node.lower_[node.variable_],
1394           node.upper_[node.variable_], node.way_,
1395           node.child_down_, node.child_up_,
1396           node.strong_branching_fixed_vars_ ? 'T' : 'F',
1397           node.reduced_cost_fixed_vars_ ? 'T' : 'F');
1398    k = node.parent_;
1399  }
1400}
1401
1402// Invoke solver's built-in enumeration algorithm
1403void 
1404branchAndBound(OsiSolverInterface & model) {
1405  double time1 = CoinCpuTime();
1406
1407  OsiClpSolverInterface* clp =
1408    dynamic_cast<OsiClpSolverInterface*>(&model);
1409  if(clp) {
1410    ClpSimplex* modelPtr = clp->getModelPtr();
1411    modelPtr->setLogLevel(0);
1412  }
1413
1414  // solve LP
1415  model.initialSolve();
1416
1417  if (model.isProvenOptimal()&&!model.isDualObjectiveLimitReached()) {
1418    // Continuous is feasible - find integers
1419    int numberIntegers=0;
1420    int numberColumns = model.getNumCols();
1421    int iColumn;
1422    int i;
1423    for (iColumn=0;iColumn<numberColumns;iColumn++) {
1424      if( model.isInteger(iColumn))
1425        numberIntegers++;
1426    }
1427    if (!numberIntegers) {
1428      std::cout<<"No integer variables"
1429               <<std::endl;
1430      return;
1431    }
1432    int * which = new int[numberIntegers]; // which variables are integer
1433    // original bounds
1434    int * originalLower = new int[numberIntegers];
1435    int * originalUpper = new int[numberIntegers];
1436    int * relaxedLower = new int[numberIntegers];
1437    int * relaxedUpper = new int[numberIntegers];
1438    {
1439      const double * lower = model.getColLower();
1440      const double * upper = model.getColUpper();
1441      numberIntegers=0;
1442      for (iColumn=0;iColumn<numberColumns;iColumn++) {
1443        if( model.isInteger(iColumn)) {
1444          originalLower[numberIntegers]=(int) lower[iColumn];
1445          originalUpper[numberIntegers]=(int) upper[iColumn];
1446          which[numberIntegers++]=iColumn;
1447        }
1448      }
1449    }
1450    double direction = model.getObjSense();
1451    // empty tree
1452    DBVectorNode branchingTree;
1453   
1454    // Add continuous to it;
1455    DBNodeSimple rootNode(0,model,numberIntegers,which,model.getWarmStart());
1456    // something extra may have been fixed by strong branching
1457    // if so go round again
1458    while (rootNode.variable_==numberIntegers) {
1459      model.resolve();
1460      rootNode = DBNodeSimple(0,model,numberIntegers,which,model.getWarmStart());
1461    }
1462    if (rootNode.objectiveValue_<1.0e100) {
1463      // push on stack
1464      branchingTree.push_back(rootNode);
1465    }
1466   
1467    // For printing totals
1468    int numberIterations=0;
1469    int numberNodes =0;
1470    DBNodeSimple bestNode;
1471    ////// Start main while of branch and bound
1472    // while until nothing on stack
1473    while (branchingTree.size()) {
1474#if defined(DEBUG_DYNAMIC_BRANCHING)
1475      if (dyn_debug >= 1000) {
1476        printf("branchingTree.size = %d %d\n",branchingTree.size(),branchingTree.size_);
1477        printf("i node_id parent child_down child_up\n");
1478        for(int k=0; k<branchingTree.size_; k++) {
1479          DBNodeSimple& node = branchingTree.nodes_[k];
1480          printf("%d %d %d %d %d\n",k, node.node_id_, node.parent_,
1481                 node.child_down_, node.child_up_);
1482        }
1483      }
1484#endif
1485      // last node
1486      DBNodeSimple node = branchingTree.back();
1487      int kNode = branchingTree.chosen_;
1488      branchingTree.pop_back();
1489#if defined(DEBUG_DYNAMIC_BRANCHING)
1490      if (dyn_debug >= 1000) {
1491        printf("Deleted current parent %d %d\n",branchingTree.size(),branchingTree.size_);
1492      }
1493#endif
1494      assert (! node.bothChildDone());
1495      numberNodes++;
1496      if (node.variable_ < 0) {
1497        // put back on tree and pretend both children are done. We want the
1498        // whole tree all the time.
1499        node.way_ = DBNodeSimple::WAY_UP_DOWN__BOTH_DONE;
1500        branchingTree.push_back(node);
1501        // Integer solution - save
1502        bestNode = node;
1503        // set cutoff (hard coded tolerance)
1504        const double limit = (bestNode.objectiveValue_-1.0e-5)*direction;
1505        model.setDblParam(OsiDualObjectiveLimit, limit);
1506        std::cout<<"Integer solution of "
1507                 <<bestNode.objectiveValue_
1508                 <<" found after "<<numberIterations
1509                 <<" iterations and "<<numberNodes<<" nodes"
1510                 <<std::endl;
1511#if defined(DEBUG_DYNAMIC_BRANCHING)
1512        if (dyn_debug >= 1) {
1513          branchingTree.checkTree();
1514        }
1515#endif
1516
1517      } else {
1518        // branch - do bounds
1519        for (i=0;i<numberIntegers;i++) {
1520          iColumn=which[i];
1521          model.setColBounds( iColumn,node.lower_[i],node.upper_[i]);
1522        }
1523        // move basis
1524        model.setWarmStart(node.basis_);
1525        // do branching variable
1526        const bool down_branch = node.advanceWay();
1527        if (down_branch) {
1528          model.setColUpper(which[node.variable_],floor(node.value_));
1529        } else {
1530          model.setColLower(which[node.variable_],ceil(node.value_));
1531        }
1532        // put back on tree anyway regardless whether any processing is left
1533        // to be done. We want the whole tree all the time.
1534        branchingTree.push_back(node);
1535#if defined(DEBUG_DYNAMIC_BRANCHING)
1536        if (dyn_debug >= 1000) {
1537          printf("Added current parent %d %d\n",branchingTree.size(),branchingTree.size_);
1538        }
1539#endif
1540       
1541        // solve
1542        model.resolve();
1543        CoinWarmStart * ws = model.getWarmStart();
1544        const CoinWarmStartBasis* wsb =
1545          dynamic_cast<const CoinWarmStartBasis*>(ws);
1546        assert (wsb!=NULL); // make sure not volume
1547        numberIterations += model.getIterationCount();
1548        // fix on reduced costs
1549        int nFixed0=0,nFixed1=0;
1550        double cutoff;
1551        model.getDblParam(OsiDualObjectiveLimit,cutoff);
1552        double gap=(cutoff-model.getObjValue())*direction+1.0e-4;
1553        //        double
1554        //        gap=(cutoff-modelPtr_->objectiveValue())*direction+1.0e-4;
1555        bool did_reduced_cost_fixing_for_child = false;
1556        if (gap<1.0e10&&model.isProvenOptimal()&&!model.isDualObjectiveLimitReached()) {
1557          const double * dj = model.getReducedCost();
1558          const double * lower = model.getColLower();
1559          const double * upper = model.getColUpper();
1560          for (i=0;i<numberIntegers;i++) {
1561            iColumn=which[i];
1562            if (upper[iColumn]>lower[iColumn]) {
1563              double djValue = dj[iColumn]*direction;
1564              if (wsb->getStructStatus(iColumn)==CoinWarmStartBasis::atLowerBound&&
1565                  djValue>gap) {
1566                nFixed0++;
1567                model.setColUpper(iColumn,lower[iColumn]);
1568              } else if (wsb->getStructStatus(iColumn)==CoinWarmStartBasis::atUpperBound&&
1569                         -djValue>gap) {
1570                nFixed1++;
1571                model.setColLower(iColumn,upper[iColumn]);
1572              }
1573            }
1574          }
1575          if (nFixed0+nFixed1) {
1576            //      printf("%d fixed to lower, %d fixed to upper\n",
1577            //             nFixed0,nFixed1);
1578            did_reduced_cost_fixing_for_child = true;
1579          }
1580        }
1581        if (model.isAbandoned()) {
1582          // THINK: What the heck should we do???
1583          abort();
1584        }
1585        if (model.isIterationLimitReached()) {
1586          // maximum iterations - exit
1587          std::cout<<"Exiting on maximum iterations\n";
1588          break;
1589        }
1590
1591        const double parentGap = (cutoff-node.objectiveValue_)*direction + 1.0e-4;
1592        assert (parentGap >= 0);
1593        const bool smallGap = false; // parentGap / fabs(cutoff) < 0.05;
1594
1595        // We are not going to do any switching unless the gap is small
1596        if (smallGap) {
1597          // THINK: If goes over the dual obj limit then clp sets dual obj limit
1598          // reached *AND* poven primal infeas. So we need to run it to
1599          // completition to know whether it's really infeas or just over the
1600          // bound. Something faster would be better...
1601          const bool mustResolve =
1602            model.isDualObjectiveLimitReached() && !model.isProvenOptimal();
1603          double oldlimit = 0;
1604          if (mustResolve) {
1605            model.getDblParam(OsiDualObjectiveLimit, oldlimit);
1606            model.setDblParam(OsiDualObjectiveLimit, 1e100);
1607            model.resolve();
1608          }
1609          LPresult lpres(model);
1610          if (mustResolve) {
1611            lpres.isDualObjectiveLimitReached = true;
1612            model.setDblParam(OsiDualObjectiveLimit, oldlimit);
1613          }
1614          bool canSwitch =
1615            node.canSwitchParentWithGrandparent(which, lpres, originalLower,
1616                                                originalUpper, branchingTree);
1617          int cnt = 0;
1618
1619#if defined(DEBUG_DYNAMIC_BRANCHING)
1620          if (dyn_debug >= 1) {
1621            branchingTree.checkTree();
1622          }
1623#endif
1624
1625#if defined(DEBUG_DYNAMIC_BRANCHING)
1626          std::string tree0;
1627          if (dyn_debug >= 10) {
1628            if (canSwitch) {
1629              tree0 = getTree(branchingTree);
1630            }
1631          }
1632#endif
1633
1634          while (canSwitch) {
1635            branchingTree.moveNodeUp(which, model, node);
1636            canSwitch = node.canSwitchParentWithGrandparent(which, lpres,
1637                                                            originalLower,
1638                                                            originalUpper,
1639                                                            branchingTree);
1640#if defined(DEBUG_DYNAMIC_BRANCHING)
1641            if (dyn_debug >= 1) {
1642              branchingTree.checkTree();
1643            }
1644#endif
1645            ++cnt;
1646          }
1647          if (cnt > 0) {
1648            printf("Before switch: opt: %c   inf: %c   dual_bd: %c\n",
1649                   lpres.isProvenOptimal ? 'T' : 'F',
1650                   lpres.isProvenPrimalInfeasible ? 'T' : 'F',
1651                   lpres.isDualObjectiveLimitReached ? 'T' : 'F');
1652            model.resolve();
1653            // This is horribly looking... Get rid of it when properly
1654            // debugged...
1655#if 1
1656            const bool mustResolve =
1657              model.isDualObjectiveLimitReached() && !model.isProvenOptimal();
1658            double oldlimit = 0;
1659           
1660            if (mustResolve) {
1661              // THINK: Something faster would be better...
1662              model.getDblParam(OsiDualObjectiveLimit, oldlimit);
1663              model.setDblParam(OsiDualObjectiveLimit, 1e100);
1664              model.resolve();
1665            }
1666            LPresult lpres1(model);
1667            if (mustResolve) {
1668              lpres.isDualObjectiveLimitReached = true;
1669              model.setDblParam(OsiDualObjectiveLimit, oldlimit);
1670            }
1671            printf("After resolve: opt: %c   inf: %c   dual_bd: %c\n",
1672                   lpres1.isProvenOptimal ? 'T' : 'F',
1673                   lpres1.isProvenPrimalInfeasible ? 'T' : 'F',
1674                   lpres1.isDualObjectiveLimitReached ? 'T' : 'F');
1675            assert(lpres.isAbandoned == model.isAbandoned());
1676            assert(lpres.isDualObjectiveLimitReached == model.isDualObjectiveLimitReached());
1677            assert(lpres.isDualObjectiveLimitReached ||
1678                   (lpres.isProvenOptimal == model.isProvenOptimal()));
1679            assert(lpres.isDualObjectiveLimitReached ||
1680                   (lpres.isProvenPrimalInfeasible == model.isProvenPrimalInfeasible()));
1681            assert(!lpres.isProvenOptimal || ! model.isProvenOptimal() ||
1682                   (lpres.isProvenOptimal && model.isProvenOptimal() &&
1683                    lpres.getObjValue == model.getObjValue()));
1684#endif
1685            printf("Finished moving node %d up by %i levels.\n", node.node_id_, cnt);
1686          }
1687         
1688#if defined(DEBUG_DYNAMIC_BRANCHING)
1689          if (dyn_debug >= 10) {
1690            if (cnt > 0) {
1691              std::string tree1 = getTree(branchingTree);
1692              printf("=====================================\n");
1693              printf("It can move node %i up. way_: 0x%x   brvar: %i\n",
1694                     node.node_id_, node.way_, node.variable_);
1695              printTree(tree0, cnt+10);
1696              printf("=====================================\n");
1697              printf("Finished moving the node up by %i levels.\n", cnt);
1698              printTree(tree1, cnt+10);
1699              printf("=====================================\n");
1700            }
1701          }
1702#endif
1703        }
1704        if ((numberNodes%10)==0) 
1705          printf("%d nodes, tree size %d\n",
1706                 numberNodes,branchingTree.size());
1707        if (CoinCpuTime()-time1>3600.0) {
1708          printf("stopping after 3600 seconds\n");
1709          exit(77);
1710        }
1711        DBNodeSimple newNode(numberNodes, model,numberIntegers,which,ws);
1712        // something extra may have been fixed by strong branching
1713        // if so go round again
1714        while (newNode.variable_==numberIntegers) {
1715          model.resolve();
1716          newNode = DBNodeSimple(numberNodes, model,numberIntegers,which,model.getWarmStart());
1717          newNode.strong_branching_fixed_vars_ = true;
1718        }
1719        newNode.reduced_cost_fixed_vars_ = did_reduced_cost_fixing_for_child;
1720        if (newNode.objectiveValue_<1.0e100) {
1721          newNode.parent_ = kNode;
1722          // push on stack
1723          branchingTree.push_back(newNode);
1724#if defined(DEBUG_DYNAMIC_BRANCHING)
1725          if (dyn_debug >= 1000) {
1726            printf("Added current child %d %d\n",branchingTree.size(),branchingTree.size_);
1727          }
1728#endif
1729          if (branchingTree.nodes_[kNode].workingOnDownChild()) {
1730            branchingTree.nodes_[kNode].child_down_ = branchingTree.last_;
1731          } else {
1732            branchingTree.nodes_[kNode].child_up_ = branchingTree.last_;
1733          }
1734          if (newNode.variable_>=0) {
1735            assert (fabs(newNode.value_-floor(newNode.value_+0.5))>1.0e-6);
1736          }
1737#if 0
1738          else {
1739            // integer solution - save
1740            bestNode = node;
1741            // set cutoff (hard coded tolerance)
1742            model.setDblParam(OsiDualObjectiveLimit,(bestNode.objectiveValue_-1.0e-5)*direction);
1743            std::cout<<"Integer solution of "
1744                     <<bestNode.objectiveValue_
1745                     <<" found after "<<numberIterations
1746                     <<" iterations and "<<numberNodes<<" nodes"
1747                     <<std::endl;
1748          }
1749#endif
1750        }
1751#if defined(DEBUG_DYNAMIC_BRANCHING)
1752        if (dyn_debug >= 1) {
1753          branchingTree.checkTree();
1754        }
1755#endif
1756      }
1757    }
1758    ////// End main while of branch and bound
1759    std::cout<<"Search took "
1760             <<numberIterations
1761             <<" iterations and "<<numberNodes<<" nodes"
1762             <<std::endl;
1763    if (bestNode.numberIntegers_) {
1764      // we have a solution restore
1765      // do bounds
1766      for (i=0;i<numberIntegers;i++) {
1767        iColumn=which[i];
1768        model.setColBounds( iColumn,bestNode.lower_[i],bestNode.upper_[i]);
1769      }
1770      // move basis
1771      model.setWarmStart(bestNode.basis_);
1772      // set cutoff so will be good (hard coded tolerance)
1773      model.setDblParam(OsiDualObjectiveLimit,(bestNode.objectiveValue_+1.0e-5)*direction);
1774      model.resolve();
1775    } else {
1776      OsiClpSolverInterface* clp =
1777        dynamic_cast<OsiClpSolverInterface*>(&model);
1778      if (clp) {
1779        ClpSimplex* modelPtr_ = clp->getModelPtr();
1780        modelPtr_->setProblemStatus(1);
1781      }
1782    }
1783    delete [] which;
1784    delete [] originalLower;
1785    delete [] originalUpper;
1786    delete [] relaxedLower;
1787    delete [] relaxedUpper;
1788  } else {
1789    std::cout<<"The LP relaxation is infeasible"
1790             <<std::endl;
1791    OsiClpSolverInterface* clp =
1792      dynamic_cast<OsiClpSolverInterface*>(&model);
1793    if (clp) {
1794      ClpSimplex* modelPtr_ = clp->getModelPtr();
1795      modelPtr_->setProblemStatus(1);
1796    }
1797    //throw CoinError("The LP relaxation is infeasible or too expensive",
1798    //"branchAndBound", "OsiClpSolverInterface");
1799  }
1800}
1801
1802
1803
1804int main(int argc, char* argv[])
1805{
1806  OsiClpSolverInterface model;
1807  model.readMps(argv[1]);
1808  branchAndBound(model);
1809  return 0;
1810}
Note: See TracBrowser for help on using the repository browser.