source: branches/dynamicbranching/dynamicbranching.cpp @ 990

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

one step closer...

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