source: branches/dynamicbranching/dynamicbranching.cpp @ 991

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

one more check...

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