source: branches/dynamicbranching/dynamicbranching.cpp @ 989

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

still not there...

File size: 49.8 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 DEBUG_DYNAMIC_BRANCHING
46
47#ifdef DEBUG_DYNAMIC_BRANCHING
48int dyn_debug = 1;
49#endif
50
51// below needed for pathetic branch and bound code
52#include <vector>
53#include <map>
54
55using namespace std;
56
57class DBVectorNode;
58
59class LPresult {
60private:
61  LPresult(const LPresult& rhs);
62  LPresult& operator=(const LPresult& rhs);
63  void gutsOfConstructor(const OsiSolverInterface& model);
64public:
65  LPresult(const OsiSolverInterface& model);
66  ~LPresult();
67public:
68  bool isAbandoned;
69  bool isProvenDualInfeasible;
70  bool isPrimalObjectiveLimitReached;
71  bool isProvenOptimal;
72  bool isDualObjectiveLimitReached;
73  bool isIterationLimitReached;
74  bool isProvenPrimalInfeasible;
75  double getObjSense;
76  double* getReducedCost;
77  double* getColLower;
78  double* getColUpper;
79  double* getObjCoefficients;
80  double yb_plus_rl_minus_su;
81};
82
83void
84LPresult::gutsOfConstructor(const OsiSolverInterface& model)
85{
86  isAbandoned = model.isAbandoned();
87  isProvenDualInfeasible = model.isProvenDualInfeasible();
88  isPrimalObjectiveLimitReached = model.isPrimalObjectiveLimitReached();
89  isProvenOptimal = model.isProvenOptimal();
90  isDualObjectiveLimitReached = model.isDualObjectiveLimitReached();
91  isIterationLimitReached = model.isIterationLimitReached();
92  isProvenPrimalInfeasible = model.isProvenPrimalInfeasible();
93  getObjSense = model.getObjSense();
94
95  getReducedCost = new double[model.getNumCols()];
96  CoinDisjointCopyN(model.getReducedCost(), model.getNumCols(), getReducedCost);
97  getColLower = new double[model.getNumCols()];
98  CoinDisjointCopyN(model.getColLower(), model.getNumCols(), getColLower);
99  getColUpper = new double[model.getNumCols()];
100  CoinDisjointCopyN(model.getColUpper(), model.getNumCols(), getColUpper);
101
102  getObjCoefficients = NULL;
103  yb_plus_rl_minus_su = 0;
104 
105  if (!isProvenOptimal &&
106      !isDualObjectiveLimitReached &&
107      !isIterationLimitReached &&
108      !isAbandoned &&
109      !isProvenDualInfeasible &&
110      !isPrimalObjectiveLimitReached) {
111    assert(isProvenPrimalInfeasible);
112    getObjCoefficients = new double[model.getNumCols()];
113    CoinDisjointCopyN(model.getObjCoefficients(), model.getNumCols(), getObjCoefficients);
114    const std::vector<double*> dual_rays = model.getDualRays(1);
115    if (dual_rays.size() == 0) {
116      printf("WARNING: LP is infeas, but no dual ray is returned!\n");
117      return;
118    }
119    const double* dual_ray = dual_rays[0];
120    const double direction = model.getObjSense();
121    const double* rub = model.getRowUpper();
122    const double* rlb = model.getRowLower();
123    const double* cub = model.getColUpper();
124    const double* clb = model.getColLower();
125    const double* dj  = model.getReducedCost();
126    const double* obj = model.getObjCoefficients();
127    const int numRows = model.getNumRows();
128    const int numCols = model.getNumCols();
129    for (int i = 0; i < numRows; ++i) {
130      const double ray_i = dual_ray[i];
131      if (ray_i > 1e-6) {
132        yb_plus_rl_minus_su += ray_i*rlb[i];
133      } else if (ray_i < 1e-6) {
134        yb_plus_rl_minus_su += ray_i*rub[i];
135      }
136    }
137    for (int j = 0; j < numCols; ++j) {
138      const double yA_j = dj[j] - obj[j];
139      if (direction * yA_j > 1e-6) {
140        yb_plus_rl_minus_su -= yA_j*cub[j];
141      } else if (direction * yA_j < -1e-6) {
142        yb_plus_rl_minus_su -= yA_j*clb[j];
143      }
144    }
145    for (int i = dual_rays.size()-1; i >= 0; --i) {
146      delete[] dual_rays[i];
147    }
148  }
149}
150 
151LPresult::LPresult(const OsiSolverInterface& model)
152{
153  gutsOfConstructor(model);
154}
155
156LPresult::~LPresult()
157{
158  delete[] getReducedCost;
159  delete[] getColLower;
160  delete[] getColUpper;
161  delete[] getObjCoefficients;
162}
163
164// Trivial class for Branch and Bound
165
166class DBNodeSimple  {
167public:
168  enum DBNodeWay {
169    WAY_UNSET=0x00,
170    WAY_DOWN_UP__NOTHING_DONE=0x10,
171    WAY_DOWN_UP__DOWN_DONE=0x11,
172    WAY_DOWN_UP__BOTH_DONE=0x13,
173    WAY_DOWN_UP__=0x10,
174    WAY_UP_DOWN__NOTHING_DONE=0x20,
175    WAY_UP_DOWN__UP_DONE=0x22,
176    WAY_UP_DOWN__BOTH_DONE=0x23,
177    WAY_UP_DOWN__=0x20,
178    WAY_DOWN_CURRENT=0x01,
179    WAY_UP_CURRENT=0x02,
180    WAY_BOTH_DONE=0x03,
181   
182  };
183
184private:
185  void gutsOfCopy(const DBNodeSimple& rhs);
186  void gutsOfConstructor (OsiSolverInterface &model,
187                          int numberIntegers, int * integer,
188                          CoinWarmStart * basis);
189  void gutsOfDestructor(); 
190 
191public:
192   
193  // Default Constructor
194  DBNodeSimple ();
195
196  // Constructor from current state (and list of integers)
197  // Also chooses branching variable (if none set to -1)
198  DBNodeSimple (OsiSolverInterface &model,
199                int numberIntegers, int * integer,
200                CoinWarmStart * basis);
201  // Copy constructor
202  DBNodeSimple ( const DBNodeSimple &);
203   
204  // Assignment operator
205  DBNodeSimple & operator=( const DBNodeSimple& rhs);
206
207  // Destructor
208  ~DBNodeSimple ();
209
210  // Extension - true if other extension of this
211  bool extension(const DBNodeSimple & other,
212                 const double * originalLower,
213                 const double * originalUpper) const;
214  // Tests if we can switch this node (this is the parent) with its parent
215  bool canSwitchParentWithGrandparent(const int* which,
216                                      const LPresult& lpres,
217                                      const int * original_lower,
218                                      const int * original_upper,
219                                      DBVectorNode & branchingTree);
220  inline bool bothChildDone() const {
221    return (way_ & WAY_BOTH_DONE) == WAY_BOTH_DONE;
222  }
223  inline bool workingOnDownChild() const {
224    return (way_ == WAY_DOWN_UP__DOWN_DONE || way_ == WAY_UP_DOWN__BOTH_DONE);
225  }
226  /* Return whether the child that will be now processed is the down branch or
227     not. */
228  inline bool advanceWay() { 
229    switch (way_) {
230    case WAY_UNSET:
231    case WAY_DOWN_UP__BOTH_DONE:
232    case WAY_UP_DOWN__BOTH_DONE:
233      abort();
234    case WAY_DOWN_UP__NOTHING_DONE:
235      way_ = WAY_DOWN_UP__DOWN_DONE;
236      return true;
237    case WAY_DOWN_UP__DOWN_DONE:
238      way_ = WAY_DOWN_UP__BOTH_DONE;
239      return false;
240    case WAY_UP_DOWN__NOTHING_DONE:
241      way_ = WAY_UP_DOWN__UP_DONE;
242      return false;
243    case WAY_UP_DOWN__UP_DONE:
244      way_ = WAY_UP_DOWN__BOTH_DONE;
245      return true;
246    }
247    return true; // to placate some compilers...
248  }
249   
250 
251  // Public data
252  // The id of the node
253  int node_id_;
254  // Basis (should use tree, but not as wasteful as bounds!)
255  CoinWarmStart * basis_;
256  // Objective value (COIN_DBL_MAX) if spare node
257  double objectiveValue_;
258  // Branching variable (0 is first integer)
259  int variable_;
260  // Way to branch: see enum DBNodeWay
261  int way_;
262  // Number of integers (for length of arrays)
263  int numberIntegers_;
264  // Current value
265  double value_;
266  // Parent
267  int parent_;
268  // Left child
269  int child_down_;
270  // Right child
271  int child_up_;
272  // Whether strong branching fixed variables when we branched on this node
273  bool strong_branching_fixed_vars_;
274  // Whether reduced cost fixing fixed anything in this node
275  bool reduced_cost_fixed_vars_;
276  // Previous in chain
277  int previous_;
278  // Next in chain
279  int next_;
280  // Now I must use tree
281  // Bounds stored in full (for integers)
282  int * lower_;
283  int * upper_;
284};
285
286
287DBNodeSimple::DBNodeSimple() :
288  node_id_(-1),
289  basis_(NULL),
290  objectiveValue_(COIN_DBL_MAX),
291  variable_(-100),
292  way_(WAY_UNSET),
293  numberIntegers_(0),
294  value_(0.5),
295  parent_(-1),
296  child_down_(-1),
297  child_up_(-1),
298  strong_branching_fixed_vars_(false),
299  reduced_cost_fixed_vars_(false),
300  previous_(-1),
301  next_(-1),
302  lower_(NULL),
303  upper_(NULL)
304{
305}
306DBNodeSimple::DBNodeSimple(OsiSolverInterface & model,
307                           int numberIntegers, int * integer,CoinWarmStart * basis)
308{
309  gutsOfConstructor(model,numberIntegers,integer,basis);
310}
311void
312DBNodeSimple::gutsOfConstructor(OsiSolverInterface & model,
313                                int numberIntegers, int * integer,CoinWarmStart * basis)
314{
315  node_id_ = -1;
316  basis_ = basis;
317  variable_=-1;
318  way_ = WAY_UNSET;
319  numberIntegers_=numberIntegers;
320  value_=0.0;
321  parent_ = -1;
322  child_down_ = -1;
323  child_up_ = -1;
324  strong_branching_fixed_vars_ = false;
325  reduced_cost_fixed_vars_ = false;
326  previous_ = -1;
327  next_ = -1;
328  if (model.isProvenOptimal()&&!model.isDualObjectiveLimitReached()) {
329    objectiveValue_ = model.getObjSense()*model.getObjValue();
330  } else {
331    objectiveValue_ = 1.0e100;
332    lower_ = NULL;
333    upper_ = NULL;
334    return; // node cutoff
335  }
336  lower_ = new int [numberIntegers_];
337  upper_ = new int [numberIntegers_];
338  assert (upper_!=NULL);
339  const double * lower = model.getColLower();
340  const double * upper = model.getColUpper();
341  const double * solution = model.getColSolution();
342  int i;
343  // Hard coded integer tolerance
344#define INTEGER_TOLERANCE 1.0e-6
345  ///////// Start of Strong branching code - can be ignored
346  // Number of strong branching candidates
347#define STRONG_BRANCHING 5
348#ifdef STRONG_BRANCHING
349  double upMovement[STRONG_BRANCHING];
350  double downMovement[STRONG_BRANCHING];
351  double solutionValue[STRONG_BRANCHING];
352  int chosen[STRONG_BRANCHING];
353  int iSmallest=0;
354  // initialize distance from integer
355  for (i=0;i<STRONG_BRANCHING;i++) {
356    upMovement[i]=0.0;
357    chosen[i]=-1;
358  }
359  variable_=-1;
360  // This has hard coded integer tolerance
361  double mostAway=INTEGER_TOLERANCE;
362  int numberAway=0;
363  for (i=0;i<numberIntegers;i++) {
364    int iColumn = integer[i];
365    lower_[i]=(int)lower[iColumn];
366    upper_[i]=(int)upper[iColumn];
367    double value = solution[iColumn];
368    value = max(value,(double) lower_[i]);
369    value = min(value,(double) upper_[i]);
370    double nearest = floor(value+0.5);
371    if (fabs(value-nearest)>INTEGER_TOLERANCE)
372      numberAway++;
373    if (fabs(value-nearest)>mostAway) {
374      double away = fabs(value-nearest);
375      if (away>upMovement[iSmallest]) {
376        //add to list
377        upMovement[iSmallest]=away;
378        solutionValue[iSmallest]=value;
379        chosen[iSmallest]=i;
380        int j;
381        iSmallest=-1;
382        double smallest = 1.0;
383        for (j=0;j<STRONG_BRANCHING;j++) {
384          if (upMovement[j]<smallest) {
385            smallest=upMovement[j];
386            iSmallest=j;
387          }
388        }
389      }
390    }
391  }
392  int numberStrong=0;
393  for (i=0;i<STRONG_BRANCHING;i++) {
394    if (chosen[i]>=0) { 
395      numberStrong ++;
396      variable_ = chosen[i];
397    }
398  }
399  // out strong branching if bit set
400  OsiClpSolverInterface* clp =
401    dynamic_cast<OsiClpSolverInterface*>(&model);
402  if (clp&&(clp->specialOptions()&16)!=0&&numberStrong>1) {
403    int j;
404    int iBest=-1;
405    double best = 0.0;
406    for (j=0;j<STRONG_BRANCHING;j++) {
407      if (upMovement[j]>best) {
408        best=upMovement[j];
409        iBest=j;
410      }
411    }
412    numberStrong=1;
413    variable_=chosen[iBest];
414  }
415  if (numberStrong==1) {
416    // just one - makes it easy
417    int iColumn = integer[variable_];
418    double value = solution[iColumn];
419    value = max(value,(double) lower_[variable_]);
420    value = min(value,(double) upper_[variable_]);
421    double nearest = floor(value+0.5);
422    value_=value;
423    if (value<=nearest)
424      way_ = WAY_UP_DOWN__NOTHING_DONE; // up
425    else
426      way_ = WAY_DOWN_UP__NOTHING_DONE; // down
427  } else if (numberStrong) {
428    // more than one - choose
429    bool chooseOne=true;
430    model.markHotStart();
431    for (i=0;i<STRONG_BRANCHING;i++) {
432      int iInt = chosen[i];
433      if (iInt>=0) {
434        int iColumn = integer[iInt];
435        double value = solutionValue[i]; // value of variable in original
436        double objectiveChange;
437        value = max(value,(double) lower_[iInt]);
438        value = min(value,(double) upper_[iInt]);
439
440        // try down
441
442        model.setColUpper(iColumn,floor(value));
443        model.solveFromHotStart();
444        model.setColUpper(iColumn,upper_[iInt]);
445        if (model.isProvenOptimal()&&!model.isDualObjectiveLimitReached()) {
446          objectiveChange = model.getObjSense()*model.getObjValue()
447            - objectiveValue_;
448        } else {
449          objectiveChange = 1.0e100;
450        }
451        assert (objectiveChange>-1.0e-5);
452        objectiveChange = CoinMax(objectiveChange,0.0);
453        downMovement[i]=objectiveChange;
454
455        // try up
456
457        model.setColLower(iColumn,ceil(value));
458        model.solveFromHotStart();
459        model.setColLower(iColumn,lower_[iInt]);
460        if (model.isProvenOptimal()&&!model.isDualObjectiveLimitReached()) {
461          objectiveChange = model.getObjSense()*model.getObjValue()
462            - objectiveValue_;
463        } else {
464          objectiveChange = 1.0e100;
465        }
466        assert (objectiveChange>-1.0e-5);
467        objectiveChange = CoinMax(objectiveChange,0.0);
468        upMovement[i]=objectiveChange;
469       
470        /* Possibilities are:
471           Both sides feasible - store
472           Neither side feasible - set objective high and exit
473           One side feasible - change bounds and resolve
474        */
475        bool solveAgain=false;
476        if (upMovement[i]<1.0e100) {
477          if(downMovement[i]<1.0e100) {
478            // feasible - no action
479          } else {
480            // up feasible, down infeasible
481            solveAgain = true;
482            model.setColLower(iColumn,ceil(value));
483          }
484        } else {
485          if(downMovement[i]<1.0e100) {
486            // down feasible, up infeasible
487            solveAgain = true;
488            model.setColUpper(iColumn,floor(value));
489          } else {
490            // neither side feasible
491            objectiveValue_=1.0e100;
492            chooseOne=false;
493            break;
494          }
495        }
496        if (solveAgain) {
497          // need to solve problem again - signal this
498          variable_ = numberIntegers;
499          chooseOne=false;
500          break;
501        }
502      }
503    }
504    if (chooseOne) {
505      // choose the one that makes most difference both ways
506      double best = -1.0;
507      double best2 = -1.0;
508      for (i=0;i<STRONG_BRANCHING;i++) {
509        int iInt = chosen[i];
510        if (iInt>=0) {
511          //std::cout<<"Strong branching on "
512          //   <<i<<""<<iInt<<" down "<<downMovement[i]
513          //   <<" up "<<upMovement[i]
514          //   <<" value "<<solutionValue[i]
515          //   <<std::endl;
516          bool better = false;
517          if (min(upMovement[i],downMovement[i])>best) {
518            // smaller is better
519            better=true;
520          } else if (min(upMovement[i],downMovement[i])>best-1.0e-5) {
521            if (max(upMovement[i],downMovement[i])>best2+1.0e-5) {
522              // smaller is about same, but larger is better
523              better=true;
524            }
525          }
526          if (better) {
527            best = min(upMovement[i],downMovement[i]);
528            best2 = max(upMovement[i],downMovement[i]);
529            variable_ = iInt;
530            double value = solutionValue[i];
531            value = max(value,(double) lower_[variable_]);
532            value = min(value,(double) upper_[variable_]);
533            value_=value;
534            if (upMovement[i]<=downMovement[i])
535              way_ = WAY_UP_DOWN__NOTHING_DONE; // up
536            else
537              way_ = WAY_DOWN_UP__NOTHING_DONE; // down
538          }
539        }
540      }
541    }
542    // Delete the snapshot
543    model.unmarkHotStart();
544  }
545  ////// End of Strong branching
546#else
547  variable_=-1;
548  // This has hard coded integer tolerance
549  double mostAway=INTEGER_TOLERANCE;
550  int numberAway=0;
551  for (i=0;i<numberIntegers;i++) {
552    int iColumn = integer[i];
553    lower_[i]=(int)lower[iColumn];
554    upper_[i]=(int)upper[iColumn];
555    double value = solution[iColumn];
556    value = max(value,(double) lower_[i]);
557    value = min(value,(double) upper_[i]);
558    double nearest = floor(value+0.5);
559    if (fabs(value-nearest)>INTEGER_TOLERANCE)
560      numberAway++;
561    if (fabs(value-nearest)>mostAway) {
562      mostAway=fabs(value-nearest);
563      variable_=i;
564      value_=value;
565      if (value<=nearest)
566        way_ = WAY_UP_DOWN__NOTHING_DONE; // up
567      else
568        way_ = WAY_DOWN_UP__NOTHING_DONE; // down
569    }
570  }
571#endif
572}
573
574void
575DBNodeSimple::gutsOfCopy(const DBNodeSimple & rhs)
576{
577  node_id_=rhs.node_id_;
578  basis_ = rhs.basis_ ? rhs.basis_->clone() : NULL;
579  objectiveValue_=rhs.objectiveValue_;
580  variable_=rhs.variable_;
581  way_=rhs.way_;
582  numberIntegers_=rhs.numberIntegers_;
583  value_=rhs.value_;
584  parent_ = rhs.parent_;
585  child_down_ = rhs.child_down_;
586  child_up_ = rhs.child_up_;
587  strong_branching_fixed_vars_ = rhs.strong_branching_fixed_vars_;
588  reduced_cost_fixed_vars_ = rhs.reduced_cost_fixed_vars_;
589  previous_ = rhs.previous_;
590  next_ = rhs.next_;
591  lower_=NULL;
592  upper_=NULL;
593  if (rhs.lower_!=NULL) {
594    lower_ = new int [numberIntegers_];
595    upper_ = new int [numberIntegers_];
596    assert (upper_!=NULL);
597    memcpy(lower_,rhs.lower_,numberIntegers_*sizeof(int));
598    memcpy(upper_,rhs.upper_,numberIntegers_*sizeof(int));
599  }
600}
601
602DBNodeSimple::DBNodeSimple(const DBNodeSimple & rhs)
603{
604  gutsOfCopy(rhs);
605}
606
607DBNodeSimple &
608DBNodeSimple::operator=(const DBNodeSimple & rhs)
609{
610  if (this != &rhs) {
611    // LL: in original code basis/lower/upper was left alone if rhs did not
612    // have them. Was that intentional?
613    gutsOfDestructor();
614    gutsOfCopy(rhs);
615  }
616  return *this;
617}
618
619
620DBNodeSimple::~DBNodeSimple ()
621{
622  gutsOfDestructor();
623}
624// Work of destructor
625void 
626DBNodeSimple::gutsOfDestructor()
627{
628  delete [] lower_;
629  delete [] upper_;
630  delete basis_;
631  lower_ = NULL;
632  upper_ = NULL;
633  basis_ = NULL;
634  objectiveValue_ = COIN_DBL_MAX;
635}
636// Extension - true if other extension of this
637bool 
638DBNodeSimple::extension(const DBNodeSimple & other,
639                        const double * originalLower,
640                        const double * originalUpper) const
641{
642  bool ok=true;
643  for (int i=0;i<numberIntegers_;i++) {
644    if (upper_[i]<originalUpper[i]||
645        lower_[i]>originalLower[i]) {
646      if (other.upper_[i]>upper_[i]||
647          other.lower_[i]<lower_[i]) {
648        ok=false;
649        break;
650      }
651    }
652  }
653  return ok;
654}
655
656#include <vector>
657#define FUNNY_BRANCHING 1
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    removeSubTree(subroot);
1073    return;
1074  }
1075  if (node.variable_ == brvar) {
1076    if (node.value_ < brlb) {
1077      // child_down_ just became infeasible. Just cut out the current node
1078      // together with its child_down_ from the tree and hang child_up_ on the
1079      // parent of this node.
1080      if (node.child_down_ >= 0) {
1081        removeSubTree(node.child_down_);
1082      }
1083      const int parent_id = node.parent_;
1084      const int child_remains = node.child_up_;
1085      if (nodes_[parent_id].child_down_ == subroot) {
1086        nodes_[parent_id].child_down_ = child_remains;
1087      } else {
1088        nodes_[parent_id].child_up_ = child_remains;
1089      }
1090      removeNode(subroot);
1091      if (child_remains >= 0) {
1092        nodes_[child_remains].parent_ = parent_id;
1093        adjustBounds(child_remains, brvar, brlb, brub);
1094      }
1095      return;
1096    }
1097    if (node.value_ > brub) {
1098      // child_up_ just became infeasible. Just cut out the current node
1099      // together with its child_down_ from the tree and hang child_down_ on
1100      // the parent of this node.
1101      if (node.child_up_ >= 0) {
1102        removeSubTree(node.child_up_);
1103      }
1104      const int parent_id = node.parent_;
1105      const int child_remains = node.child_down_;
1106      if (nodes_[parent_id].child_down_ == subroot) {
1107        nodes_[parent_id].child_down_ = child_remains;
1108      } else {
1109        nodes_[parent_id].child_up_ = child_remains;
1110      }
1111      removeNode(subroot);
1112      if (child_remains >= 0) {
1113        nodes_[child_remains].parent_ = parent_id;
1114        adjustBounds(child_remains, brvar, brlb, brub);
1115      }
1116      return;
1117    }
1118    // Now brlb < node.value_ < brub (value_ is fractional)
1119  }
1120  node.lower_[brvar] = brlb;
1121  node.upper_[brvar] = brub;
1122  if (node.child_down_ >= 0) {
1123    adjustBounds(node.child_down_, brvar, brlb, brub);
1124  }
1125  if (node.child_up_ >= 0) {
1126    adjustBounds(node.child_up_, brvar, brlb, brub);
1127  }
1128}
1129
1130void
1131DBVectorNode::moveNodeUp(const int* which,
1132                         OsiSolverInterface& model, DBNodeSimple & node)
1133{
1134  /*
1135    The current node ('this') is really the parent (P). The grandparent (GP)
1136    is this.parent_. Let's have variable names respecting the truth.
1137  */
1138  const int parent_id = node.node_id_;
1139  DBNodeSimple& parent = nodes_[parent_id];
1140  const int grandparent_id = parent.parent_;
1141  assert(grandparent_id != -1);
1142  DBNodeSimple& grandparent = nodes_[grandparent_id];
1143  const int greatgrandparent_id = grandparent.parent_;
1144 
1145  const bool parent_is_down_child = parent_id == grandparent.child_down_;
1146
1147#if defined(DEBUG_DYNAMIC_BRANCHING)
1148  if (dyn_debug >= 100) {
1149    printf("entered moveNodeUp\n");
1150    printf("parent_id %d grandparent_id %d greatgrandparent_id %d\n",
1151           parent_id, grandparent_id, greatgrandparent_id);
1152    printf("parent.way_ %d\n", parent.way_);
1153  }
1154#endif
1155
1156
1157  // First hang the nodes where they belong.
1158  parent.parent_ = greatgrandparent_id;
1159  grandparent.parent_ = parent_id;
1160  const bool down_child_stays_with_parent = parent.workingOnDownChild();
1161  int& child_to_move =
1162    down_child_stays_with_parent ? parent.child_up_ : parent.child_down_;
1163  const bool child_to_move_is_processed = parent.bothChildDone();
1164
1165#if defined(DEBUG_DYNAMIC_BRANCHING)
1166  if (dyn_debug >= 1000) {
1167    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);
1168  }
1169#endif
1170
1171  if (parent_is_down_child) {
1172    grandparent.child_down_ = child_to_move;
1173  } else {
1174    grandparent.child_up_ = child_to_move;
1175  }
1176  if (child_to_move >= 0) {
1177    nodes_[child_to_move].parent_ = grandparent_id;
1178  }
1179  child_to_move = grandparent_id;
1180  if (greatgrandparent_id >= 0) {
1181    DBNodeSimple& greatgrandparent = nodes_[greatgrandparent_id];
1182    if (grandparent_id == greatgrandparent.child_down_) {
1183      greatgrandparent.child_down_ = parent_id;
1184    } else {
1185      greatgrandparent.child_up_ = parent_id;
1186    }
1187  }
1188
1189#if defined(DEBUG_DYNAMIC_BRANCHING)
1190  if (dyn_debug >= 1000) {
1191    printf("after exchange\n");
1192    printf("parent.parent_ %d parent.child_down_ %d parent.child_up_ %d\n",
1193           parent.parent_, parent.child_down_, parent.child_up_);
1194    printf("grandparent.parent_ %d grandparent.child_down_ %d grandparent.child_up_ %d\n",
1195           grandparent.parent_, grandparent.child_down_, grandparent.child_up_);
1196    if (greatgrandparent_id >= 0) {
1197      DBNodeSimple& greatgrandparent = nodes_[greatgrandparent_id];
1198      printf("greatgrandparent.parent_ %d greatgrandparent.child_down_ %d greatgrandparent.child_up_ %d\n",
1199             greatgrandparent.parent_, greatgrandparent.child_down_, greatgrandparent.child_up_);
1200    }
1201    printf("exiting moveNodeUp\n");
1202  }
1203#endif
1204
1205
1206  // Now make sure way_ is set properly
1207  if (down_child_stays_with_parent) {
1208    if (!parent.bothChildDone()) {
1209      parent.way_ = DBNodeSimple::WAY_UP_DOWN__BOTH_DONE;
1210      sizeDeferred_++;
1211    }
1212    if (grandparent.bothChildDone()) {
1213      if (! child_to_move_is_processed) {
1214        grandparent.way_ = parent_is_down_child ?
1215          DBNodeSimple::WAY_UP_DOWN__UP_DONE :
1216          DBNodeSimple::WAY_DOWN_UP__DOWN_DONE;
1217        sizeDeferred_--;
1218      }
1219    } else { // only parent is processed from the two children of grandparent
1220      if (! child_to_move_is_processed) {
1221        grandparent.way_ = parent_is_down_child ?
1222          DBNodeSimple::WAY_DOWN_UP__NOTHING_DONE :
1223          DBNodeSimple::WAY_UP_DOWN__NOTHING_DONE;
1224      } else {
1225        grandparent.way_ = parent_is_down_child ?
1226          DBNodeSimple::WAY_DOWN_UP__DOWN_DONE :
1227          DBNodeSimple::WAY_UP_DOWN__UP_DONE;
1228      }
1229    }
1230  } else {
1231    if (!parent.bothChildDone()) {
1232      parent.way_ = DBNodeSimple::WAY_DOWN_UP__BOTH_DONE;
1233      sizeDeferred_++;
1234    }
1235    if (grandparent.bothChildDone()) {
1236      if (! child_to_move_is_processed) {
1237        grandparent.way_ = parent_is_down_child ?
1238          DBNodeSimple::WAY_UP_DOWN__UP_DONE :
1239          DBNodeSimple::WAY_DOWN_UP__DOWN_DONE;
1240        sizeDeferred_--;
1241      }
1242    } else { // only parent is processed from the two children of grandparent
1243      if (! child_to_move_is_processed) {
1244        grandparent.way_ = parent_is_down_child ?
1245          DBNodeSimple::WAY_DOWN_UP__NOTHING_DONE :
1246          DBNodeSimple::WAY_UP_DOWN__NOTHING_DONE;
1247      } else {
1248        grandparent.way_ = parent_is_down_child ?
1249          DBNodeSimple::WAY_DOWN_UP__DOWN_DONE :
1250          DBNodeSimple::WAY_UP_DOWN__UP_DONE;
1251      }
1252    }
1253  }
1254 
1255  // Now modify bounds
1256
1257  // First, get rid of GP's bound change of its branching variable in the
1258  // bound list of P. Note: If greatgrandparent_id is < 0 then GP is the root
1259  // so its bounds are the original bounds.
1260  const int GP_brvar = grandparent.variable_;
1261  if (parent_is_down_child) {
1262    const int old_upper = grandparent.upper_[GP_brvar];
1263    parent.upper_[GP_brvar] = old_upper;
1264    if ((GP_brvar != parent.variable_) ||
1265        (GP_brvar == parent.variable_ && !down_child_stays_with_parent)) {
1266      model.setColUpper(which[GP_brvar], old_upper);
1267    }
1268  } else {
1269    const int old_lower = grandparent.lower_[GP_brvar];
1270    parent.lower_[GP_brvar] = old_lower;
1271    if ((GP_brvar != parent.variable_) ||
1272        (GP_brvar == parent.variable_ && down_child_stays_with_parent)) {
1273      model.setColLower(which[GP_brvar], old_lower);
1274    }
1275  }
1276
1277  // Now add the branching var bound change of P to GP and all of its
1278  // descendant
1279  if (down_child_stays_with_parent) {
1280    adjustBounds(grandparent_id, parent.variable_,
1281                 (int)ceil(parent.value_), parent.upper_[parent.variable_]);
1282  } else {
1283    adjustBounds(grandparent_id, parent.variable_,
1284                 parent.lower_[parent.variable_], (int)floor(parent.value_));
1285  }
1286}
1287
1288void
1289DBVectorNode::checkNode(int node) const
1290{
1291  if (node == -1) {
1292    return;
1293  }
1294  const DBNodeSimple* n = nodes_ + node;
1295  const DBNodeSimple* p = nodes_ + n->parent_;
1296  for (int i = n->numberIntegers_-1; i >= 0; --i) {
1297    if (i == p->variable_) {
1298      if (node == p->child_down_) {
1299        assert(p->lower_[i] <= n->lower_[i]);
1300        assert((int)(floor(p->value_)) == n->upper_[i]);
1301      } else {
1302        assert((int)(ceil(p->value_)) == n->lower_[i]);
1303        assert(p->upper_[i] >= n->upper_[i]);
1304      }
1305    } else {
1306      assert(p->lower_[i] <= n->lower_[i]);
1307      assert(p->upper_[i] >= n->upper_[i]);
1308    }
1309  }
1310  checkNode(n->child_down_);
1311  checkNode(n->child_up_);
1312}
1313
1314void
1315DBVectorNode::checkTree() const
1316{
1317  // find the root
1318  int root = first_;
1319  while (true) {
1320    if (nodes_[root].parent_ == -1) {
1321      break;
1322    } else {
1323      root = nodes_[root].parent_;
1324    }
1325  }
1326  checkNode(nodes_[root].child_down_);
1327  checkNode(nodes_[root].child_up_);
1328}
1329
1330std::string getTree(DBVectorNode& branchingTree)
1331{
1332  std::string tree;
1333  char line[1000];
1334  for(int k=0; k<branchingTree.size_; k++) {
1335    DBNodeSimple& node = branchingTree.nodes_[k];
1336    sprintf(line, "%d %d %d %d %f %d %d 0x%x %d %d\n",
1337            k, node.node_id_, node.parent_, node.variable_,
1338            node.value_, node.lower_[node.variable_],
1339            node.upper_[node.variable_], node.way_,
1340            node.child_down_, node.child_up_);
1341    tree += line;
1342  }
1343  return tree;
1344}
1345
1346void printTree(const std::string& tree, int levels)
1347{
1348  size_t pos = tree.size();
1349  for (int i = levels-1; i >= 0; --i) {
1350    pos = tree.rfind('\n', pos-1);
1351    if (pos == std::string::npos) {
1352      pos = 0;
1353      break;
1354    }
1355  }
1356  printf("%s", tree.c_str() + (pos+1));
1357}
1358
1359void printChain(DBVectorNode& branchingTree, int k)
1360{
1361  while (k != -1) {
1362    DBNodeSimple& node = branchingTree.nodes_[k];
1363    printf("   %d %d %d %d %f %d %d 0x%x %d %d %c %c\n",
1364           k, node.node_id_, node.parent_, node.variable_,
1365           node.value_, node.lower_[node.variable_],
1366           node.upper_[node.variable_], node.way_,
1367           node.child_down_, node.child_up_,
1368           node.strong_branching_fixed_vars_ ? 'T' : 'F',
1369           node.reduced_cost_fixed_vars_ ? 'T' : 'F');
1370    k = node.parent_;
1371  }
1372}
1373
1374// Invoke solver's built-in enumeration algorithm
1375void 
1376branchAndBound(OsiSolverInterface & model) {
1377  double time1 = CoinCpuTime();
1378  // solve LP
1379  model.initialSolve();
1380
1381  if (model.isProvenOptimal()&&!model.isDualObjectiveLimitReached()) {
1382    // Continuous is feasible - find integers
1383    int numberIntegers=0;
1384    int numberColumns = model.getNumCols();
1385    int iColumn;
1386    int i;
1387    for (iColumn=0;iColumn<numberColumns;iColumn++) {
1388      if( model.isInteger(iColumn))
1389        numberIntegers++;
1390    }
1391    if (!numberIntegers) {
1392      std::cout<<"No integer variables"
1393               <<std::endl;
1394      return;
1395    }
1396    int * which = new int[numberIntegers]; // which variables are integer
1397    // original bounds
1398    int * originalLower = new int[numberIntegers];
1399    int * originalUpper = new int[numberIntegers];
1400    int * relaxedLower = new int[numberIntegers];
1401    int * relaxedUpper = new int[numberIntegers];
1402    {
1403      const double * lower = model.getColLower();
1404      const double * upper = model.getColUpper();
1405      numberIntegers=0;
1406      for (iColumn=0;iColumn<numberColumns;iColumn++) {
1407        if( model.isInteger(iColumn)) {
1408          originalLower[numberIntegers]=(int) lower[iColumn];
1409          originalUpper[numberIntegers]=(int) upper[iColumn];
1410          which[numberIntegers++]=iColumn;
1411        }
1412      }
1413    }
1414    double direction = model.getObjSense();
1415    // empty tree
1416    DBVectorNode branchingTree;
1417   
1418    // Add continuous to it;
1419    DBNodeSimple rootNode(model,numberIntegers,which,model.getWarmStart());
1420    // something extra may have been fixed by strong branching
1421    // if so go round again
1422    while (rootNode.variable_==numberIntegers) {
1423      model.resolve();
1424      rootNode = DBNodeSimple(model,numberIntegers,which,model.getWarmStart());
1425    }
1426    if (rootNode.objectiveValue_<1.0e100) {
1427      // push on stack
1428      branchingTree.push_back(rootNode);
1429    }
1430   
1431    // For printing totals
1432    int numberIterations=0;
1433    int numberNodes =0;
1434    DBNodeSimple bestNode;
1435    ////// Start main while of branch and bound
1436    // while until nothing on stack
1437    while (branchingTree.size()) {
1438#if defined(DEBUG_DYNAMIC_BRANCHING)
1439      if (dyn_debug >= 1000) {
1440        printf("branchingTree.size = %d %d\n",branchingTree.size(),branchingTree.size_);
1441        printf("i node_id parent child_down child_up\n");
1442        for(int k=0; k<branchingTree.size_; k++) {
1443          DBNodeSimple& node = branchingTree.nodes_[k];
1444          printf("%d %d %d %d %d\n",k, node.node_id_, node.parent_,
1445                 node.child_down_, node.child_up_);
1446        }
1447      }
1448#endif
1449      // last node
1450      DBNodeSimple node = branchingTree.back();
1451      int kNode = branchingTree.chosen_;
1452      branchingTree.pop_back();
1453#if defined(DEBUG_DYNAMIC_BRANCHING)
1454      if (dyn_debug >= 1000) {
1455        printf("Deleted current parent %d %d\n",branchingTree.size(),branchingTree.size_);
1456      }
1457#endif
1458      assert (! node.bothChildDone());
1459      numberNodes++;
1460      if (node.variable_ < 0) {
1461        // put back on tree and pretend both children are done. We want the
1462        // whole tree all the time.
1463        node.way_ = DBNodeSimple::WAY_UP_DOWN__BOTH_DONE;
1464        branchingTree.push_back(node);
1465        // Integer solution - save
1466        bestNode = node;
1467        // set cutoff (hard coded tolerance)
1468        const double limit = (bestNode.objectiveValue_-1.0e-5)*direction;
1469        model.setDblParam(OsiDualObjectiveLimit, limit);
1470        std::cout<<"Integer solution of "
1471                 <<bestNode.objectiveValue_
1472                 <<" found after "<<numberIterations
1473                 <<" iterations and "<<numberNodes<<" nodes"
1474                 <<std::endl;
1475#if defined(DEBUG_DYNAMIC_BRANCHING)
1476        if (dyn_debug >= 1) {
1477          branchingTree.checkTree();
1478        }
1479#endif
1480
1481      } else {
1482        // branch - do bounds
1483        for (i=0;i<numberIntegers;i++) {
1484          iColumn=which[i];
1485          model.setColBounds( iColumn,node.lower_[i],node.upper_[i]);
1486        }
1487        // move basis
1488        model.setWarmStart(node.basis_);
1489        // do branching variable
1490        const bool down_branch = node.advanceWay();
1491        if (down_branch) {
1492          model.setColUpper(which[node.variable_],floor(node.value_));
1493        } else {
1494          model.setColLower(which[node.variable_],ceil(node.value_));
1495        }
1496        // put back on tree anyway regardless whether any processing is left
1497        // to be done. We want the whole tree all the time.
1498        branchingTree.push_back(node);
1499#if defined(DEBUG_DYNAMIC_BRANCHING)
1500        if (dyn_debug >= 1000) {
1501          printf("Added current parent %d %d\n",branchingTree.size(),branchingTree.size_);
1502        }
1503#endif
1504       
1505        // solve
1506        model.resolve();
1507        CoinWarmStart * ws = model.getWarmStart();
1508        const CoinWarmStartBasis* wsb =
1509          dynamic_cast<const CoinWarmStartBasis*>(ws);
1510        assert (wsb!=NULL); // make sure not volume
1511        numberIterations += model.getIterationCount();
1512        // fix on reduced costs
1513        int nFixed0=0,nFixed1=0;
1514        double cutoff;
1515        model.getDblParam(OsiDualObjectiveLimit,cutoff);
1516        double gap=(cutoff-model.getObjValue())*direction+1.0e-4;
1517        //        double
1518        //        gap=(cutoff-modelPtr_->objectiveValue())*direction+1.0e-4;
1519        bool did_reduced_cost_fixing_for_child = false;
1520        if (gap<1.0e10&&model.isProvenOptimal()&&!model.isDualObjectiveLimitReached()) {
1521          const double * dj = model.getReducedCost();
1522          const double * lower = model.getColLower();
1523          const double * upper = model.getColUpper();
1524          for (i=0;i<numberIntegers;i++) {
1525            iColumn=which[i];
1526            if (upper[iColumn]>lower[iColumn]) {
1527              double djValue = dj[iColumn]*direction;
1528              if (wsb->getStructStatus(iColumn)==CoinWarmStartBasis::atLowerBound&&
1529                  djValue>gap) {
1530                nFixed0++;
1531                model.setColUpper(iColumn,lower[iColumn]);
1532              } else if (wsb->getStructStatus(iColumn)==CoinWarmStartBasis::atUpperBound&&
1533                         -djValue>gap) {
1534                nFixed1++;
1535                model.setColLower(iColumn,upper[iColumn]);
1536              }
1537            }
1538          }
1539          if (nFixed0+nFixed1) {
1540            //      printf("%d fixed to lower, %d fixed to upper\n",
1541            //             nFixed0,nFixed1);
1542            did_reduced_cost_fixing_for_child = true;
1543          }
1544        }
1545        if (model.isAbandoned()) {
1546          // THINK: What the heck should we do???
1547          abort();
1548        }
1549        if (model.isIterationLimitReached()) {
1550          // maximum iterations - exit
1551          std::cout<<"Exiting on maximum iterations\n";
1552          break;
1553        }
1554
1555        LPresult lpres(model);
1556
1557        bool canSwitch = node.canSwitchParentWithGrandparent(which, lpres,
1558                                                             originalLower,
1559                                                             originalUpper,
1560                                                             branchingTree);
1561        int cnt = 0;
1562
1563#if defined(DEBUG_DYNAMIC_BRANCHING)
1564        if (dyn_debug >= 1) {
1565          branchingTree.checkTree();
1566        }
1567#endif
1568
1569#if defined(DEBUG_DYNAMIC_BRANCHING)
1570        std::string tree0;
1571        if (dyn_debug >= 10) {
1572          if (canSwitch) {
1573            tree0 = getTree(branchingTree);
1574          }
1575        }
1576#endif
1577
1578        while (canSwitch) {
1579          branchingTree.moveNodeUp(which, model, node);
1580          canSwitch = node.canSwitchParentWithGrandparent(which, lpres,
1581                                                          originalLower,
1582                                                          originalUpper,
1583                                                          branchingTree);
1584#if defined(DEBUG_DYNAMIC_BRANCHING)
1585          if (dyn_debug >= 1) {
1586            branchingTree.checkTree();
1587          }
1588#endif
1589          ++cnt;
1590        }
1591        if (cnt > 0) {
1592          model.resolve();
1593          // This is horribly looking... Get rid of it when properly debugged...
1594          assert(lpres.isAbandoned == model.isAbandoned());
1595          assert(lpres.isDualObjectiveLimitReached == model.isDualObjectiveLimitReached());
1596          assert(lpres.isDualObjectiveLimitReached ||
1597                 (lpres.isProvenOptimal == model.isProvenOptimal()));
1598          assert(lpres.isDualObjectiveLimitReached ||
1599                 (lpres.isProvenPrimalInfeasible == model.isProvenPrimalInfeasible()));
1600        }
1601           
1602#if defined(DEBUG_DYNAMIC_BRANCHING)
1603        if (dyn_debug >= 10) {
1604          if (cnt > 0) {
1605            std::string tree1 = getTree(branchingTree);
1606            printf("=====================================\n");
1607            printf("It can move node %i up. way_: 0x%x   brvar: %i\n",
1608                   node.node_id_, node.way_, node.variable_);
1609            printTree(tree0, cnt+10);
1610            printf("=====================================\n");
1611            printf("Finished moving the node up by %i levels.\n", cnt);
1612            printTree(tree1, cnt+10);
1613            printf("=====================================\n");
1614          }
1615        }
1616#endif
1617        if ((numberNodes%1000)==0) 
1618          printf("%d nodes, tree size %d\n",
1619                 numberNodes,branchingTree.size());
1620        if (CoinCpuTime()-time1>3600.0) {
1621          printf("stopping after 3600 seconds\n");
1622          exit(77);
1623        }
1624        DBNodeSimple newNode(model,numberIntegers,which,ws);
1625        // something extra may have been fixed by strong branching
1626        // if so go round again
1627        while (newNode.variable_==numberIntegers) {
1628          model.resolve();
1629          newNode = DBNodeSimple(model,numberIntegers,which,model.getWarmStart());
1630          newNode.strong_branching_fixed_vars_ = true;
1631        }
1632        newNode.reduced_cost_fixed_vars_ = did_reduced_cost_fixing_for_child;
1633        if (newNode.objectiveValue_<1.0e100) {
1634          newNode.parent_ = kNode;
1635          // push on stack
1636          branchingTree.push_back(newNode);
1637#if defined(DEBUG_DYNAMIC_BRANCHING)
1638          if (dyn_debug >= 1000) {
1639            printf("Added current child %d %d\n",branchingTree.size(),branchingTree.size_);
1640          }
1641#endif
1642          if (branchingTree.nodes_[kNode].workingOnDownChild()) {
1643            branchingTree.nodes_[kNode].child_down_ = branchingTree.last_;
1644          } else {
1645            branchingTree.nodes_[kNode].child_up_ = branchingTree.last_;
1646          }
1647          if (newNode.variable_>=0) {
1648            assert (fabs(newNode.value_-floor(newNode.value_+0.5))>1.0e-6);
1649          }
1650#if 0
1651          else {
1652            // integer solution - save
1653            bestNode = node;
1654            // set cutoff (hard coded tolerance)
1655            model.setDblParam(OsiDualObjectiveLimit,(bestNode.objectiveValue_-1.0e-5)*direction);
1656            std::cout<<"Integer solution of "
1657                     <<bestNode.objectiveValue_
1658                     <<" found after "<<numberIterations
1659                     <<" iterations and "<<numberNodes<<" nodes"
1660                     <<std::endl;
1661          }
1662#endif
1663        }
1664#if defined(DEBUG_DYNAMIC_BRANCHING)
1665        if (dyn_debug >= 1) {
1666          branchingTree.checkTree();
1667        }
1668#endif
1669      }
1670    }
1671    ////// End main while of branch and bound
1672    std::cout<<"Search took "
1673             <<numberIterations
1674             <<" iterations and "<<numberNodes<<" nodes"
1675             <<std::endl;
1676    if (bestNode.numberIntegers_) {
1677      // we have a solution restore
1678      // do bounds
1679      for (i=0;i<numberIntegers;i++) {
1680        iColumn=which[i];
1681        model.setColBounds( iColumn,bestNode.lower_[i],bestNode.upper_[i]);
1682      }
1683      // move basis
1684      model.setWarmStart(bestNode.basis_);
1685      // set cutoff so will be good (hard coded tolerance)
1686      model.setDblParam(OsiDualObjectiveLimit,(bestNode.objectiveValue_+1.0e-5)*direction);
1687      model.resolve();
1688    } else {
1689      OsiClpSolverInterface* clp =
1690        dynamic_cast<OsiClpSolverInterface*>(&model);
1691      if (clp) {
1692        ClpSimplex* modelPtr_ = clp->getModelPtr();
1693        modelPtr_->setProblemStatus(1);
1694      }
1695    }
1696    delete [] which;
1697    delete [] originalLower;
1698    delete [] originalUpper;
1699    delete [] relaxedLower;
1700    delete [] relaxedUpper;
1701  } else {
1702    std::cout<<"The LP relaxation is infeasible"
1703             <<std::endl;
1704    OsiClpSolverInterface* clp =
1705      dynamic_cast<OsiClpSolverInterface*>(&model);
1706    if (clp) {
1707      ClpSimplex* modelPtr_ = clp->getModelPtr();
1708      modelPtr_->setProblemStatus(1);
1709    }
1710    //throw CoinError("The LP relaxation is infeasible or too expensive",
1711    //"branchAndBound", "OsiClpSolverInterface");
1712  }
1713}
1714
1715
1716
1717int main(int argc, char* argv[])
1718{
1719  OsiClpSolverInterface model;
1720  model.readMps(argv[1]);
1721  branchAndBound(model);
1722  return 0;
1723}
Note: See TracBrowser for help on using the repository browser.