source: branches/dynamicbranching/dynamicbranching.cpp @ 996

Last change on this file since 996 was 996, checked in by ladanyi, 13 years ago
File size: 53.0 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  // solve LP
1407  model.initialSolve();
1408  model.setLogLevel(0);
1409
1410  if (model.isProvenOptimal()&&!model.isDualObjectiveLimitReached()) {
1411    // Continuous is feasible - find integers
1412    int numberIntegers=0;
1413    int numberColumns = model.getNumCols();
1414    int iColumn;
1415    int i;
1416    for (iColumn=0;iColumn<numberColumns;iColumn++) {
1417      if( model.isInteger(iColumn))
1418        numberIntegers++;
1419    }
1420    if (!numberIntegers) {
1421      std::cout<<"No integer variables"
1422               <<std::endl;
1423      return;
1424    }
1425    int * which = new int[numberIntegers]; // which variables are integer
1426    // original bounds
1427    int * originalLower = new int[numberIntegers];
1428    int * originalUpper = new int[numberIntegers];
1429    int * relaxedLower = new int[numberIntegers];
1430    int * relaxedUpper = new int[numberIntegers];
1431    {
1432      const double * lower = model.getColLower();
1433      const double * upper = model.getColUpper();
1434      numberIntegers=0;
1435      for (iColumn=0;iColumn<numberColumns;iColumn++) {
1436        if( model.isInteger(iColumn)) {
1437          originalLower[numberIntegers]=(int) lower[iColumn];
1438          originalUpper[numberIntegers]=(int) upper[iColumn];
1439          which[numberIntegers++]=iColumn;
1440        }
1441      }
1442    }
1443    double direction = model.getObjSense();
1444    // empty tree
1445    DBVectorNode branchingTree;
1446   
1447    // Add continuous to it;
1448    DBNodeSimple rootNode(0,model,numberIntegers,which,model.getWarmStart());
1449    // something extra may have been fixed by strong branching
1450    // if so go round again
1451    while (rootNode.variable_==numberIntegers) {
1452      model.resolve();
1453      rootNode = DBNodeSimple(0,model,numberIntegers,which,model.getWarmStart());
1454    }
1455    if (rootNode.objectiveValue_<1.0e100) {
1456      // push on stack
1457      branchingTree.push_back(rootNode);
1458    }
1459   
1460    // For printing totals
1461    int numberIterations=0;
1462    int numberNodes =0;
1463    DBNodeSimple bestNode;
1464    ////// Start main while of branch and bound
1465    // while until nothing on stack
1466    while (branchingTree.size()) {
1467#if defined(DEBUG_DYNAMIC_BRANCHING)
1468      if (dyn_debug >= 1000) {
1469        printf("branchingTree.size = %d %d\n",branchingTree.size(),branchingTree.size_);
1470        printf("i node_id parent child_down child_up\n");
1471        for(int k=0; k<branchingTree.size_; k++) {
1472          DBNodeSimple& node = branchingTree.nodes_[k];
1473          printf("%d %d %d %d %d\n",k, node.node_id_, node.parent_,
1474                 node.child_down_, node.child_up_);
1475        }
1476      }
1477#endif
1478      // last node
1479      DBNodeSimple node = branchingTree.back();
1480      int kNode = branchingTree.chosen_;
1481      branchingTree.pop_back();
1482#if defined(DEBUG_DYNAMIC_BRANCHING)
1483      if (dyn_debug >= 1000) {
1484        printf("Deleted current parent %d %d\n",branchingTree.size(),branchingTree.size_);
1485      }
1486#endif
1487      assert (! node.bothChildDone());
1488      numberNodes++;
1489      if (node.variable_ < 0) {
1490        // put back on tree and pretend both children are done. We want the
1491        // whole tree all the time.
1492        node.way_ = DBNodeSimple::WAY_UP_DOWN__BOTH_DONE;
1493        branchingTree.push_back(node);
1494        // Integer solution - save
1495        bestNode = node;
1496        // set cutoff (hard coded tolerance)
1497        const double limit = (bestNode.objectiveValue_-1.0e-5)*direction;
1498        model.setDblParam(OsiDualObjectiveLimit, limit);
1499        std::cout<<"Integer solution of "
1500                 <<bestNode.objectiveValue_
1501                 <<" found after "<<numberIterations
1502                 <<" iterations and "<<numberNodes<<" nodes"
1503                 <<std::endl;
1504#if defined(DEBUG_DYNAMIC_BRANCHING)
1505        if (dyn_debug >= 1) {
1506          branchingTree.checkTree();
1507        }
1508#endif
1509
1510      } else {
1511        // branch - do bounds
1512        for (i=0;i<numberIntegers;i++) {
1513          iColumn=which[i];
1514          model.setColBounds( iColumn,node.lower_[i],node.upper_[i]);
1515        }
1516        // move basis
1517        model.setWarmStart(node.basis_);
1518        // do branching variable
1519        const bool down_branch = node.advanceWay();
1520        if (down_branch) {
1521          model.setColUpper(which[node.variable_],floor(node.value_));
1522        } else {
1523          model.setColLower(which[node.variable_],ceil(node.value_));
1524        }
1525        // put back on tree anyway regardless whether any processing is left
1526        // to be done. We want the whole tree all the time.
1527        branchingTree.push_back(node);
1528#if defined(DEBUG_DYNAMIC_BRANCHING)
1529        if (dyn_debug >= 1000) {
1530          printf("Added current parent %d %d\n",branchingTree.size(),branchingTree.size_);
1531        }
1532#endif
1533       
1534        // solve
1535        model.resolve();
1536        CoinWarmStart * ws = model.getWarmStart();
1537        const CoinWarmStartBasis* wsb =
1538          dynamic_cast<const CoinWarmStartBasis*>(ws);
1539        assert (wsb!=NULL); // make sure not volume
1540        numberIterations += model.getIterationCount();
1541        // fix on reduced costs
1542        int nFixed0=0,nFixed1=0;
1543        double cutoff;
1544        model.getDblParam(OsiDualObjectiveLimit,cutoff);
1545        double gap=(cutoff-model.getObjValue())*direction+1.0e-4;
1546        //        double
1547        //        gap=(cutoff-modelPtr_->objectiveValue())*direction+1.0e-4;
1548        bool did_reduced_cost_fixing_for_child = false;
1549        if (gap<1.0e10&&model.isProvenOptimal()&&!model.isDualObjectiveLimitReached()) {
1550          const double * dj = model.getReducedCost();
1551          const double * lower = model.getColLower();
1552          const double * upper = model.getColUpper();
1553          for (i=0;i<numberIntegers;i++) {
1554            iColumn=which[i];
1555            if (upper[iColumn]>lower[iColumn]) {
1556              double djValue = dj[iColumn]*direction;
1557              if (wsb->getStructStatus(iColumn)==CoinWarmStartBasis::atLowerBound&&
1558                  djValue>gap) {
1559                nFixed0++;
1560                model.setColUpper(iColumn,lower[iColumn]);
1561              } else if (wsb->getStructStatus(iColumn)==CoinWarmStartBasis::atUpperBound&&
1562                         -djValue>gap) {
1563                nFixed1++;
1564                model.setColLower(iColumn,upper[iColumn]);
1565              }
1566            }
1567          }
1568          if (nFixed0+nFixed1) {
1569            //      printf("%d fixed to lower, %d fixed to upper\n",
1570            //             nFixed0,nFixed1);
1571            did_reduced_cost_fixing_for_child = true;
1572          }
1573        }
1574        if (model.isAbandoned()) {
1575          // THINK: What the heck should we do???
1576          abort();
1577        }
1578        if (model.isIterationLimitReached()) {
1579          // maximum iterations - exit
1580          std::cout<<"Exiting on maximum iterations\n";
1581          break;
1582        }
1583
1584        const double parentGap = (cutoff-node.objectiveValue_)*direction + 1.0e-4;
1585        assert (parentGap >= 0);
1586        const bool smallGap = false; // parentGap / fabs(cutoff) < 0.05;
1587
1588        // We are not going to do any switching unless the gap is small
1589        if (smallGap) {
1590          // THINK: If goes over the dual obj limit then clp sets dual obj limit
1591          // reached *AND* poven primal infeas. So we need to run it to
1592          // completition to know whether it's really infeas or just over the
1593          // bound. Something faster would be better...
1594          const bool mustResolve =
1595            model.isDualObjectiveLimitReached() && !model.isProvenOptimal();
1596          double oldlimit = 0;
1597          if (mustResolve) {
1598            model.getDblParam(OsiDualObjectiveLimit, oldlimit);
1599            model.setDblParam(OsiDualObjectiveLimit, 1e100);
1600            model.resolve();
1601          }
1602          LPresult lpres(model);
1603          if (mustResolve) {
1604            lpres.isDualObjectiveLimitReached = true;
1605            model.setDblParam(OsiDualObjectiveLimit, oldlimit);
1606          }
1607          bool canSwitch =
1608            node.canSwitchParentWithGrandparent(which, lpres, originalLower,
1609                                                originalUpper, branchingTree);
1610          int cnt = 0;
1611
1612#if defined(DEBUG_DYNAMIC_BRANCHING)
1613          if (dyn_debug >= 1) {
1614            branchingTree.checkTree();
1615          }
1616#endif
1617
1618#if defined(DEBUG_DYNAMIC_BRANCHING)
1619          std::string tree0;
1620          if (dyn_debug >= 10) {
1621            if (canSwitch) {
1622              tree0 = getTree(branchingTree);
1623            }
1624          }
1625#endif
1626
1627          while (canSwitch) {
1628            branchingTree.moveNodeUp(which, model, node);
1629            canSwitch = node.canSwitchParentWithGrandparent(which, lpres,
1630                                                            originalLower,
1631                                                            originalUpper,
1632                                                            branchingTree);
1633#if defined(DEBUG_DYNAMIC_BRANCHING)
1634            if (dyn_debug >= 1) {
1635              branchingTree.checkTree();
1636            }
1637#endif
1638            ++cnt;
1639          }
1640          if (cnt > 0) {
1641            printf("Before switch: opt: %c   inf: %c   dual_bd: %c\n",
1642                   lpres.isProvenOptimal ? 'T' : 'F',
1643                   lpres.isProvenPrimalInfeasible ? 'T' : 'F',
1644                   lpres.isDualObjectiveLimitReached ? 'T' : 'F');
1645            model.resolve();
1646            // This is horribly looking... Get rid of it when properly
1647            // debugged...
1648#if 1
1649            const bool mustResolve =
1650              model.isDualObjectiveLimitReached() && !model.isProvenOptimal();
1651            double oldlimit = 0;
1652           
1653            if (mustResolve) {
1654              // THINK: Something faster would be better...
1655              model.getDblParam(OsiDualObjectiveLimit, oldlimit);
1656              model.setDblParam(OsiDualObjectiveLimit, 1e100);
1657              model.resolve();
1658            }
1659            LPresult lpres1(model);
1660            if (mustResolve) {
1661              lpres.isDualObjectiveLimitReached = true;
1662              model.setDblParam(OsiDualObjectiveLimit, oldlimit);
1663            }
1664            printf("After resolve: opt: %c   inf: %c   dual_bd: %c\n",
1665                   lpres1.isProvenOptimal ? 'T' : 'F',
1666                   lpres1.isProvenPrimalInfeasible ? 'T' : 'F',
1667                   lpres1.isDualObjectiveLimitReached ? 'T' : 'F');
1668            assert(lpres.isAbandoned == model.isAbandoned());
1669            assert(lpres.isDualObjectiveLimitReached == model.isDualObjectiveLimitReached());
1670            assert(lpres.isDualObjectiveLimitReached ||
1671                   (lpres.isProvenOptimal == model.isProvenOptimal()));
1672            assert(lpres.isDualObjectiveLimitReached ||
1673                   (lpres.isProvenPrimalInfeasible == model.isProvenPrimalInfeasible()));
1674            assert(!lpres.isProvenOptimal || ! model.isProvenOptimal() ||
1675                   (lpres.isProvenOptimal && model.isProvenOptimal() &&
1676                    lpres.getObjValue == model.getObjValue()));
1677#endif
1678            printf("Finished moving node %d up by %i levels.\n", node.node_id_, cnt);
1679          }
1680         
1681#if defined(DEBUG_DYNAMIC_BRANCHING)
1682          if (dyn_debug >= 10) {
1683            if (cnt > 0) {
1684              std::string tree1 = getTree(branchingTree);
1685              printf("=====================================\n");
1686              printf("It can move node %i up. way_: 0x%x   brvar: %i\n",
1687                     node.node_id_, node.way_, node.variable_);
1688              printTree(tree0, cnt+10);
1689              printf("=====================================\n");
1690              printf("Finished moving the node up by %i levels.\n", cnt);
1691              printTree(tree1, cnt+10);
1692              printf("=====================================\n");
1693            }
1694          }
1695#endif
1696        }
1697        if ((numberNodes%10)==0) 
1698          printf("%d nodes, tree size %d\n",
1699                 numberNodes,branchingTree.size());
1700        if (CoinCpuTime()-time1>3600.0) {
1701          printf("stopping after 3600 seconds\n");
1702          exit(77);
1703        }
1704        DBNodeSimple newNode(numberNodes, model,numberIntegers,which,ws);
1705        // something extra may have been fixed by strong branching
1706        // if so go round again
1707        while (newNode.variable_==numberIntegers) {
1708          model.resolve();
1709          newNode = DBNodeSimple(numberNodes, model,numberIntegers,which,model.getWarmStart());
1710          newNode.strong_branching_fixed_vars_ = true;
1711        }
1712        newNode.reduced_cost_fixed_vars_ = did_reduced_cost_fixing_for_child;
1713        if (newNode.objectiveValue_<1.0e100) {
1714          newNode.parent_ = kNode;
1715          // push on stack
1716          branchingTree.push_back(newNode);
1717#if defined(DEBUG_DYNAMIC_BRANCHING)
1718          if (dyn_debug >= 1000) {
1719            printf("Added current child %d %d\n",branchingTree.size(),branchingTree.size_);
1720          }
1721#endif
1722          if (branchingTree.nodes_[kNode].workingOnDownChild()) {
1723            branchingTree.nodes_[kNode].child_down_ = branchingTree.last_;
1724          } else {
1725            branchingTree.nodes_[kNode].child_up_ = branchingTree.last_;
1726          }
1727          if (newNode.variable_>=0) {
1728            assert (fabs(newNode.value_-floor(newNode.value_+0.5))>1.0e-6);
1729          }
1730#if 0
1731          else {
1732            // integer solution - save
1733            bestNode = node;
1734            // set cutoff (hard coded tolerance)
1735            model.setDblParam(OsiDualObjectiveLimit,(bestNode.objectiveValue_-1.0e-5)*direction);
1736            std::cout<<"Integer solution of "
1737                     <<bestNode.objectiveValue_
1738                     <<" found after "<<numberIterations
1739                     <<" iterations and "<<numberNodes<<" nodes"
1740                     <<std::endl;
1741          }
1742#endif
1743        }
1744#if defined(DEBUG_DYNAMIC_BRANCHING)
1745        if (dyn_debug >= 1) {
1746          branchingTree.checkTree();
1747        }
1748#endif
1749      }
1750    }
1751    ////// End main while of branch and bound
1752    std::cout<<"Search took "
1753             <<numberIterations
1754             <<" iterations and "<<numberNodes<<" nodes"
1755             <<std::endl;
1756    if (bestNode.numberIntegers_) {
1757      // we have a solution restore
1758      // do bounds
1759      for (i=0;i<numberIntegers;i++) {
1760        iColumn=which[i];
1761        model.setColBounds( iColumn,bestNode.lower_[i],bestNode.upper_[i]);
1762      }
1763      // move basis
1764      model.setWarmStart(bestNode.basis_);
1765      // set cutoff so will be good (hard coded tolerance)
1766      model.setDblParam(OsiDualObjectiveLimit,(bestNode.objectiveValue_+1.0e-5)*direction);
1767      model.resolve();
1768    } else {
1769      OsiClpSolverInterface* clp =
1770        dynamic_cast<OsiClpSolverInterface*>(&model);
1771      if (clp) {
1772        ClpSimplex* modelPtr_ = clp->getModelPtr();
1773        modelPtr_->setProblemStatus(1);
1774      }
1775    }
1776    delete [] which;
1777    delete [] originalLower;
1778    delete [] originalUpper;
1779    delete [] relaxedLower;
1780    delete [] relaxedUpper;
1781  } else {
1782    std::cout<<"The LP relaxation is infeasible"
1783             <<std::endl;
1784    OsiClpSolverInterface* clp =
1785      dynamic_cast<OsiClpSolverInterface*>(&model);
1786    if (clp) {
1787      ClpSimplex* modelPtr_ = clp->getModelPtr();
1788      modelPtr_->setProblemStatus(1);
1789    }
1790    //throw CoinError("The LP relaxation is infeasible or too expensive",
1791    //"branchAndBound", "OsiClpSolverInterface");
1792  }
1793}
1794
1795
1796
1797int main(int argc, char* argv[])
1798{
1799  OsiClpSolverInterface model;
1800  model.readMps(argv[1]);
1801  branchAndBound(model);
1802  return 0;
1803}
Note: See TracBrowser for help on using the repository browser.