source: branches/dynamicbranching/dynamicbranching.cpp @ 988

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

still not there...

File size: 49.7 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    return false;
944  }
945 
946  const double direction = lpres.getObjSense;
947
948  const int GP_brvar = grandparent.variable_;
949  const int GP_brvar_fullid = which[GP_brvar];
950  const bool parent_is_down_child = parent_id == grandparent.child_down_;
951
952  if (lpres.isProvenOptimal ||
953      lpres.isDualObjectiveLimitReached ||
954      lpres.isIterationLimitReached) {
955    // Dual feasible, and in this case we don't care how we have
956    // stopped (iteration limit, obj val limit, time limit, optimal solution,
957    // etc.), we can just look at the reduced costs to figure out if the
958    // grandparent is irrelevant. Remember, we are using dual simplex!
959    double djValue = lpres.getReducedCost[GP_brvar_fullid]*direction;
960    if (djValue > 1.0e-6) {
961      // wants to go down
962      if (parent_is_down_child) {
963        return true;
964      }
965      if (lpres.getColLower[GP_brvar_fullid] > std::ceil(grandparent.value_)) {
966        return true;
967      }
968    } else if (djValue < -1.0e-6) {
969      // wants to go up
970      if (! parent_is_down_child) {
971        return true;
972      }
973      if (lpres.getColUpper[GP_brvar_fullid] < std::floor(grandparent.value_)) {
974        return true;
975      }
976    }
977    return false;
978  } else {
979    assert(lpres.isProvenPrimalInfeasible);
980    return false;
981    const int greatgrandparent_id = grandparent.parent_;
982    const int x = GP_brvar_fullid; // for easier reference... we'll use s_x
983
984    /*
985      Now we are going to check that if we relax the GP's branching decision
986      then the child's LP relaxation is still infeasible. The test is done by
987      checking whether the column (or its negative) of the GP's branching
988      variable cuts off the dual ray proving the infeasibility.
989    */
990   
991    const double* dj = lpres.getReducedCost;
992    const double* obj = lpres.getObjCoefficients;
993    const double yA_x = dj[x] - obj[x];
994
995    if (direction > 0) { // minimization problem
996      if (parent_is_down_child) {
997        const double s_x = CoinMax(yA_x, -1e-8);
998        if (s_x < 1e-6) {
999          return true;
1000        }
1001        if (lpres.yb_plus_rl_minus_su < 1e-8) {
1002          printf("WARNING: lpres.yb_plus_rl_minus_su is not positive!\n");
1003          return false;
1004        }
1005        const double max_u_x = lpres.yb_plus_rl_minus_su / s_x + lpres.getColUpper[x];
1006        const double u_x_without_GP = greatgrandparent_id >= 0 ?
1007          branchingTree.nodes_[greatgrandparent_id].upper_[GP_brvar] :
1008          original_upper[GP_brvar];
1009        return max_u_x > u_x_without_GP - 1e-8;
1010      } else {
1011        const double r_x = CoinMax(yA_x, -1e-8);
1012        if (r_x < 1e-6) {
1013          return true;
1014        }
1015        if (lpres.yb_plus_rl_minus_su < 1e-8) {
1016          printf("WARNING: lpres.yb_plus_rl_minus_su is not positive!\n");
1017          return false;
1018        }
1019        const double min_l_x = - lpres.yb_plus_rl_minus_su / r_x + lpres.getColLower[x];
1020        const double l_x_without_GP = greatgrandparent_id >= 0 ?
1021          branchingTree.nodes_[greatgrandparent_id].lower_[GP_brvar] :
1022          original_lower[GP_brvar];
1023        return min_l_x < l_x_without_GP + 1e-8;
1024      }
1025    } else { // maximization problem
1026      if (parent_is_down_child) {
1027        const double s_x = CoinMin(yA_x, 1e-8);
1028        if (s_x > -1e-6) {
1029          return true;
1030        }
1031        if (lpres.yb_plus_rl_minus_su > -1e-8) {
1032          printf("WARNING: lpres.yb_plus_rl_minus_su is not negative!\n");
1033          return false;
1034        }
1035        const double max_u_x = lpres.yb_plus_rl_minus_su / s_x + lpres.getColUpper[x];
1036        const double u_x_without_GP = greatgrandparent_id >= 0 ?
1037          branchingTree.nodes_[greatgrandparent_id].upper_[GP_brvar] :
1038          original_upper[GP_brvar];
1039        return max_u_x > u_x_without_GP - 1e-8;
1040      } else {
1041        const double r_x = CoinMin(yA_x, 1e-8);
1042        if (r_x < -1e-6) {
1043          return true;
1044        }
1045        if (lpres.yb_plus_rl_minus_su > -1e-8) {
1046          printf("WARNING: lpres.yb_plus_rl_minus_su is not negative!\n");
1047          return false;
1048        }
1049        const double min_l_x = - lpres.yb_plus_rl_minus_su / r_x + lpres.getColLower[x];
1050        const double l_x_without_GP = greatgrandparent_id >= 0 ?
1051          branchingTree.nodes_[greatgrandparent_id].lower_[GP_brvar] :
1052          original_lower[GP_brvar];
1053        return min_l_x < l_x_without_GP + 1e-8;
1054      }
1055    }
1056  }
1057
1058  return true; // to placate some compilers
1059}
1060
1061void
1062DBVectorNode::adjustBounds(int subroot, int brvar, int brlb, int brub)
1063{
1064  assert(subroot != -1);
1065  DBNodeSimple& node = nodes_[subroot];
1066  // Take the intersection of brvar's bounds in node and those passed in
1067  brub = CoinMin(brub, node.upper_[brvar]);
1068  brlb = CoinMax(brlb, node.lower_[brvar]);
1069  if (brub < brlb) {
1070    // This node became infeasible. Get rid of it and of its descendants
1071    removeSubTree(subroot);
1072    return;
1073  }
1074  if (node.variable_ == brvar) {
1075    if (node.value_ < brlb) {
1076      // child_down_ just became infeasible. Just cut out the current node
1077      // together with its child_down_ from the tree and hang child_up_ on the
1078      // parent of this node.
1079      if (node.child_down_ >= 0) {
1080        removeSubTree(node.child_down_);
1081      }
1082      const int parent_id = node.parent_;
1083      const int child_remains = node.child_up_;
1084      if (nodes_[parent_id].child_down_ == subroot) {
1085        nodes_[parent_id].child_down_ = child_remains;
1086      } else {
1087        nodes_[parent_id].child_up_ = child_remains;
1088      }
1089      removeNode(subroot);
1090      if (child_remains >= 0) {
1091        nodes_[child_remains].parent_ = parent_id;
1092        adjustBounds(child_remains, brvar, brlb, brub);
1093      }
1094      return;
1095    }
1096    if (node.value_ > brub) {
1097      // child_up_ just became infeasible. Just cut out the current node
1098      // together with its child_down_ from the tree and hang child_down_ on
1099      // the parent of this node.
1100      if (node.child_up_ >= 0) {
1101        removeSubTree(node.child_up_);
1102      }
1103      const int parent_id = node.parent_;
1104      const int child_remains = node.child_down_;
1105      if (nodes_[parent_id].child_down_ == subroot) {
1106        nodes_[parent_id].child_down_ = child_remains;
1107      } else {
1108        nodes_[parent_id].child_up_ = child_remains;
1109      }
1110      removeNode(subroot);
1111      if (child_remains >= 0) {
1112        nodes_[child_remains].parent_ = parent_id;
1113        adjustBounds(child_remains, brvar, brlb, brub);
1114      }
1115      return;
1116    }
1117    // Now brlb < node.value_ < brub (value_ is fractional)
1118  }
1119  node.lower_[brvar] = brlb;
1120  node.upper_[brvar] = brub;
1121  if (node.child_down_ >= 0) {
1122    adjustBounds(node.child_down_, brvar, brlb, brub);
1123  }
1124  if (node.child_up_ >= 0) {
1125    adjustBounds(node.child_up_, brvar, brlb, brub);
1126  }
1127}
1128
1129void
1130DBVectorNode::moveNodeUp(const int* which,
1131                         OsiSolverInterface& model, DBNodeSimple & node)
1132{
1133  /*
1134    The current node ('this') is really the parent (P). The grandparent (GP)
1135    is this.parent_. Let's have variable names respecting the truth.
1136  */
1137  const int parent_id = node.node_id_;
1138  DBNodeSimple& parent = nodes_[parent_id];
1139  const int grandparent_id = parent.parent_;
1140  assert(grandparent_id != -1);
1141  DBNodeSimple& grandparent = nodes_[grandparent_id];
1142  const int greatgrandparent_id = grandparent.parent_;
1143 
1144  const bool parent_is_down_child = parent_id == grandparent.child_down_;
1145
1146#if defined(DEBUG_DYNAMIC_BRANCHING)
1147  if (dyn_debug >= 100) {
1148    printf("entered moveNodeUp\n");
1149    printf("parent_id %d grandparent_id %d greatgrandparent_id %d\n",
1150           parent_id, grandparent_id, greatgrandparent_id);
1151    printf("parent.way_ %d\n", parent.way_);
1152  }
1153#endif
1154
1155
1156  // First hang the nodes where they belong.
1157  parent.parent_ = greatgrandparent_id;
1158  grandparent.parent_ = parent_id;
1159  const bool down_child_stays_with_parent = parent.workingOnDownChild();
1160  int& child_to_move =
1161    down_child_stays_with_parent ? parent.child_up_ : parent.child_down_;
1162  const bool child_to_move_is_processed = parent.bothChildDone();
1163
1164#if defined(DEBUG_DYNAMIC_BRANCHING)
1165  if (dyn_debug >= 1000) {
1166    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);
1167  }
1168#endif
1169
1170  if (parent_is_down_child) {
1171    grandparent.child_down_ = child_to_move;
1172  } else {
1173    grandparent.child_up_ = child_to_move;
1174  }
1175  if (child_to_move >= 0) {
1176    nodes_[child_to_move].parent_ = grandparent_id;
1177  }
1178  child_to_move = grandparent_id;
1179  if (greatgrandparent_id >= 0) {
1180    DBNodeSimple& greatgrandparent = nodes_[greatgrandparent_id];
1181    if (grandparent_id == greatgrandparent.child_down_) {
1182      greatgrandparent.child_down_ = parent_id;
1183    } else {
1184      greatgrandparent.child_up_ = parent_id;
1185    }
1186  }
1187
1188#if defined(DEBUG_DYNAMIC_BRANCHING)
1189  if (dyn_debug >= 1000) {
1190    printf("after exchange\n");
1191    printf("parent.parent_ %d parent.child_down_ %d parent.child_up_ %d\n",
1192           parent.parent_, parent.child_down_, parent.child_up_);
1193    printf("grandparent.parent_ %d grandparent.child_down_ %d grandparent.child_up_ %d\n",
1194           grandparent.parent_, grandparent.child_down_, grandparent.child_up_);
1195    if (greatgrandparent_id >= 0) {
1196      DBNodeSimple& greatgrandparent = nodes_[greatgrandparent_id];
1197      printf("greatgrandparent.parent_ %d greatgrandparent.child_down_ %d greatgrandparent.child_up_ %d\n",
1198             greatgrandparent.parent_, greatgrandparent.child_down_, greatgrandparent.child_up_);
1199    }
1200    printf("exiting moveNodeUp\n");
1201  }
1202#endif
1203
1204
1205  // Now make sure way_ is set properly
1206  if (down_child_stays_with_parent) {
1207    if (!parent.bothChildDone()) {
1208      parent.way_ = DBNodeSimple::WAY_UP_DOWN__BOTH_DONE;
1209      sizeDeferred_++;
1210    }
1211    if (grandparent.bothChildDone()) {
1212      if (! child_to_move_is_processed) {
1213        grandparent.way_ = parent_is_down_child ?
1214          DBNodeSimple::WAY_UP_DOWN__UP_DONE :
1215          DBNodeSimple::WAY_DOWN_UP__DOWN_DONE;
1216        sizeDeferred_--;
1217      }
1218    } else { // only parent is processed from the two children of grandparent
1219      if (! child_to_move_is_processed) {
1220        grandparent.way_ = parent_is_down_child ?
1221          DBNodeSimple::WAY_DOWN_UP__NOTHING_DONE :
1222          DBNodeSimple::WAY_UP_DOWN__NOTHING_DONE;
1223      } else {
1224        grandparent.way_ = parent_is_down_child ?
1225          DBNodeSimple::WAY_DOWN_UP__DOWN_DONE :
1226          DBNodeSimple::WAY_UP_DOWN__UP_DONE;
1227      }
1228    }
1229  } else {
1230    if (!parent.bothChildDone()) {
1231      parent.way_ = DBNodeSimple::WAY_DOWN_UP__BOTH_DONE;
1232      sizeDeferred_++;
1233    }
1234    if (grandparent.bothChildDone()) {
1235      if (! child_to_move_is_processed) {
1236        grandparent.way_ = parent_is_down_child ?
1237          DBNodeSimple::WAY_UP_DOWN__UP_DONE :
1238          DBNodeSimple::WAY_DOWN_UP__DOWN_DONE;
1239        sizeDeferred_--;
1240      }
1241    } else { // only parent is processed from the two children of grandparent
1242      if (! child_to_move_is_processed) {
1243        grandparent.way_ = parent_is_down_child ?
1244          DBNodeSimple::WAY_DOWN_UP__NOTHING_DONE :
1245          DBNodeSimple::WAY_UP_DOWN__NOTHING_DONE;
1246      } else {
1247        grandparent.way_ = parent_is_down_child ?
1248          DBNodeSimple::WAY_DOWN_UP__DOWN_DONE :
1249          DBNodeSimple::WAY_UP_DOWN__UP_DONE;
1250      }
1251    }
1252  }
1253 
1254  // Now modify bounds
1255
1256  // First, get rid of GP's bound change of its branching variable in the
1257  // bound list of P. Note: If greatgrandparent_id is < 0 then GP is the root
1258  // so its bounds are the original bounds.
1259  const int GP_brvar = grandparent.variable_;
1260  if (parent_is_down_child) {
1261    const int old_upper = grandparent.upper_[GP_brvar];
1262    parent.upper_[GP_brvar] = old_upper;
1263    if ((GP_brvar != parent.variable_) ||
1264        (GP_brvar == parent.variable_ && !down_child_stays_with_parent)) {
1265      model.setColUpper(which[GP_brvar], old_upper);
1266    }
1267  } else {
1268    const int old_lower = grandparent.lower_[GP_brvar];
1269    parent.lower_[GP_brvar] = old_lower;
1270    if ((GP_brvar != parent.variable_) ||
1271        (GP_brvar == parent.variable_ && down_child_stays_with_parent)) {
1272      model.setColLower(which[GP_brvar], old_lower);
1273    }
1274  }
1275
1276  // Now add the branching var bound change of P to GP and all of its
1277  // descendant
1278  if (down_child_stays_with_parent) {
1279    adjustBounds(grandparent_id, parent.variable_,
1280                 (int)ceil(parent.value_), parent.upper_[parent.variable_]);
1281  } else {
1282    adjustBounds(grandparent_id, parent.variable_,
1283                 parent.lower_[parent.variable_], (int)floor(parent.value_));
1284  }
1285}
1286
1287void
1288DBVectorNode::checkNode(int node) const
1289{
1290  if (node == -1) {
1291    return;
1292  }
1293  const DBNodeSimple* n = nodes_ + node;
1294  const DBNodeSimple* p = nodes_ + n->parent_;
1295  for (int i = n->numberIntegers_-1; i >= 0; --i) {
1296    if (i == p->variable_) {
1297      if (node == p->child_down_) {
1298        assert(p->lower_[i] <= n->lower_[i]);
1299        assert((int)(floor(p->value_)) == n->upper_[i]);
1300      } else {
1301        assert((int)(ceil(p->value_)) == n->lower_[i]);
1302        assert(p->upper_[i] >= n->upper_[i]);
1303      }
1304    } else {
1305      assert(p->lower_[i] <= n->lower_[i]);
1306      assert(p->upper_[i] >= n->upper_[i]);
1307    }
1308  }
1309  checkNode(n->child_down_);
1310  checkNode(n->child_up_);
1311}
1312
1313void
1314DBVectorNode::checkTree() const
1315{
1316  // find the root
1317  int root = first_;
1318  while (true) {
1319    if (nodes_[root].parent_ == -1) {
1320      break;
1321    }
1322  }
1323  checkNode(nodes_[root].child_down_);
1324  checkNode(nodes_[root].child_up_);
1325}
1326
1327std::string getTree(DBVectorNode& branchingTree)
1328{
1329  std::string tree;
1330  char line[1000];
1331  for(int k=0; k<branchingTree.size_; k++) {
1332    DBNodeSimple& node = branchingTree.nodes_[k];
1333    sprintf(line, "%d %d %d %d %f %d %d 0x%x %d %d\n",
1334            k, node.node_id_, node.parent_, node.variable_,
1335            node.value_, node.lower_[node.variable_],
1336            node.upper_[node.variable_], node.way_,
1337            node.child_down_, node.child_up_);
1338    tree += line;
1339  }
1340  return tree;
1341}
1342
1343void printTree(const std::string& tree, int levels)
1344{
1345  size_t pos = tree.size();
1346  for (int i = levels-1; i >= 0; --i) {
1347    pos = tree.rfind('\n', pos-1);
1348    if (pos == std::string::npos) {
1349      pos = 0;
1350      break;
1351    }
1352  }
1353  printf("%s", tree.c_str() + (pos+1));
1354}
1355
1356void printChain(DBVectorNode& branchingTree, int k)
1357{
1358  while (k != -1) {
1359    DBNodeSimple& node = branchingTree.nodes_[k];
1360    printf("   %d %d %d %d %f %d %d 0x%x %d %d %c %c\n",
1361           k, node.node_id_, node.parent_, node.variable_,
1362           node.value_, node.lower_[node.variable_],
1363           node.upper_[node.variable_], node.way_,
1364           node.child_down_, node.child_up_,
1365           node.strong_branching_fixed_vars_ ? 'T' : 'F',
1366           node.reduced_cost_fixed_vars_ ? 'T' : 'F');
1367    k = node.parent_;
1368  }
1369}
1370
1371// Invoke solver's built-in enumeration algorithm
1372void 
1373branchAndBound(OsiSolverInterface & model) {
1374  double time1 = CoinCpuTime();
1375  // solve LP
1376  model.initialSolve();
1377
1378  if (model.isProvenOptimal()&&!model.isDualObjectiveLimitReached()) {
1379    // Continuous is feasible - find integers
1380    int numberIntegers=0;
1381    int numberColumns = model.getNumCols();
1382    int iColumn;
1383    int i;
1384    for (iColumn=0;iColumn<numberColumns;iColumn++) {
1385      if( model.isInteger(iColumn))
1386        numberIntegers++;
1387    }
1388    if (!numberIntegers) {
1389      std::cout<<"No integer variables"
1390               <<std::endl;
1391      return;
1392    }
1393    int * which = new int[numberIntegers]; // which variables are integer
1394    // original bounds
1395    int * originalLower = new int[numberIntegers];
1396    int * originalUpper = new int[numberIntegers];
1397    int * relaxedLower = new int[numberIntegers];
1398    int * relaxedUpper = new int[numberIntegers];
1399    {
1400      const double * lower = model.getColLower();
1401      const double * upper = model.getColUpper();
1402      numberIntegers=0;
1403      for (iColumn=0;iColumn<numberColumns;iColumn++) {
1404        if( model.isInteger(iColumn)) {
1405          originalLower[numberIntegers]=(int) lower[iColumn];
1406          originalUpper[numberIntegers]=(int) upper[iColumn];
1407          which[numberIntegers++]=iColumn;
1408        }
1409      }
1410    }
1411    double direction = model.getObjSense();
1412    // empty tree
1413    DBVectorNode branchingTree;
1414   
1415    // Add continuous to it;
1416    DBNodeSimple rootNode(model,numberIntegers,which,model.getWarmStart());
1417    // something extra may have been fixed by strong branching
1418    // if so go round again
1419    while (rootNode.variable_==numberIntegers) {
1420      model.resolve();
1421      rootNode = DBNodeSimple(model,numberIntegers,which,model.getWarmStart());
1422    }
1423    if (rootNode.objectiveValue_<1.0e100) {
1424      // push on stack
1425      branchingTree.push_back(rootNode);
1426    }
1427   
1428    // For printing totals
1429    int numberIterations=0;
1430    int numberNodes =0;
1431    DBNodeSimple bestNode;
1432    ////// Start main while of branch and bound
1433    // while until nothing on stack
1434    while (branchingTree.size()) {
1435#if defined(DEBUG_DYNAMIC_BRANCHING)
1436      if (dyn_debug >= 1000) {
1437        printf("branchingTree.size = %d %d\n",branchingTree.size(),branchingTree.size_);
1438        printf("i node_id parent child_down child_up\n");
1439        for(int k=0; k<branchingTree.size_; k++) {
1440          DBNodeSimple& node = branchingTree.nodes_[k];
1441          printf("%d %d %d %d %d\n",k, node.node_id_, node.parent_,
1442                 node.child_down_, node.child_up_);
1443        }
1444      }
1445#endif
1446      // last node
1447      DBNodeSimple node = branchingTree.back();
1448      int kNode = branchingTree.chosen_;
1449      branchingTree.pop_back();
1450#if defined(DEBUG_DYNAMIC_BRANCHING)
1451      if (dyn_debug >= 1000) {
1452        printf("Deleted current parent %d %d\n",branchingTree.size(),branchingTree.size_);
1453      }
1454#endif
1455      assert (! node.bothChildDone());
1456      numberNodes++;
1457      if (node.variable_ < 0) {
1458        // put back on tree and pretend both children are done. We want the
1459        // whole tree all the time.
1460        node.way_ = DBNodeSimple::WAY_UP_DOWN__BOTH_DONE;
1461        branchingTree.push_back(node);
1462        // Integer solution - save
1463        bestNode = node;
1464        // set cutoff (hard coded tolerance)
1465        const double limit = (bestNode.objectiveValue_-1.0e-5)*direction;
1466        model.setDblParam(OsiDualObjectiveLimit, limit);
1467        std::cout<<"Integer solution of "
1468                 <<bestNode.objectiveValue_
1469                 <<" found after "<<numberIterations
1470                 <<" iterations and "<<numberNodes<<" nodes"
1471                 <<std::endl;
1472#if defined(DEBUG_DYNAMIC_BRANCHING)
1473        if (dyn_debug >= 1) {
1474          branchingTree.checkTree();
1475        }
1476#endif
1477
1478      } else {
1479        // branch - do bounds
1480        for (i=0;i<numberIntegers;i++) {
1481          iColumn=which[i];
1482          model.setColBounds( iColumn,node.lower_[i],node.upper_[i]);
1483        }
1484        // move basis
1485        model.setWarmStart(node.basis_);
1486        // do branching variable
1487        const bool down_branch = node.advanceWay();
1488        if (down_branch) {
1489          model.setColUpper(which[node.variable_],floor(node.value_));
1490        } else {
1491          model.setColLower(which[node.variable_],ceil(node.value_));
1492        }
1493        // put back on tree anyway regardless whether any processing is left
1494        // to be done. We want the whole tree all the time.
1495        branchingTree.push_back(node);
1496#if defined(DEBUG_DYNAMIC_BRANCHING)
1497        if (dyn_debug >= 1000) {
1498          printf("Added current parent %d %d\n",branchingTree.size(),branchingTree.size_);
1499        }
1500#endif
1501       
1502        // solve
1503        model.resolve();
1504        CoinWarmStart * ws = model.getWarmStart();
1505        const CoinWarmStartBasis* wsb =
1506          dynamic_cast<const CoinWarmStartBasis*>(ws);
1507        assert (wsb!=NULL); // make sure not volume
1508        numberIterations += model.getIterationCount();
1509        // fix on reduced costs
1510        int nFixed0=0,nFixed1=0;
1511        double cutoff;
1512        model.getDblParam(OsiDualObjectiveLimit,cutoff);
1513        double gap=(cutoff-model.getObjValue())*direction+1.0e-4;
1514        //        double
1515        //        gap=(cutoff-modelPtr_->objectiveValue())*direction+1.0e-4;
1516        bool did_reduced_cost_fixing_for_child = false;
1517        if (gap<1.0e10&&model.isProvenOptimal()&&!model.isDualObjectiveLimitReached()) {
1518          const double * dj = model.getReducedCost();
1519          const double * lower = model.getColLower();
1520          const double * upper = model.getColUpper();
1521          for (i=0;i<numberIntegers;i++) {
1522            iColumn=which[i];
1523            if (upper[iColumn]>lower[iColumn]) {
1524              double djValue = dj[iColumn]*direction;
1525              if (wsb->getStructStatus(iColumn)==CoinWarmStartBasis::atLowerBound&&
1526                  djValue>gap) {
1527                nFixed0++;
1528                model.setColUpper(iColumn,lower[iColumn]);
1529              } else if (wsb->getStructStatus(iColumn)==CoinWarmStartBasis::atUpperBound&&
1530                         -djValue>gap) {
1531                nFixed1++;
1532                model.setColLower(iColumn,upper[iColumn]);
1533              }
1534            }
1535          }
1536          if (nFixed0+nFixed1) {
1537            //      printf("%d fixed to lower, %d fixed to upper\n",
1538            //             nFixed0,nFixed1);
1539            did_reduced_cost_fixing_for_child = true;
1540          }
1541        }
1542        if (model.isAbandoned()) {
1543          // THINK: What the heck should we do???
1544          abort();
1545        }
1546        if (model.isIterationLimitReached()) {
1547          // maximum iterations - exit
1548          std::cout<<"Exiting on maximum iterations\n";
1549          break;
1550        }
1551
1552        LPresult lpres(model);
1553
1554        bool canSwitch = node.canSwitchParentWithGrandparent(which, lpres,
1555                                                             originalLower,
1556                                                             originalUpper,
1557                                                             branchingTree);
1558        int cnt = 0;
1559
1560#if defined(DEBUG_DYNAMIC_BRANCHING)
1561        if (dyn_debug >= 1) {
1562          branchingTree.checkTree();
1563        }
1564#endif
1565
1566#if defined(DEBUG_DYNAMIC_BRANCHING)
1567        std::string tree0;
1568        if (dyn_debug >= 10) {
1569          if (canSwitch) {
1570            tree0 = getTree(branchingTree);
1571          }
1572        }
1573#endif
1574
1575        while (canSwitch) {
1576          branchingTree.moveNodeUp(which, model, node);
1577          canSwitch = node.canSwitchParentWithGrandparent(which, lpres,
1578                                                          originalLower,
1579                                                          originalUpper,
1580                                                          branchingTree);
1581#if defined(DEBUG_DYNAMIC_BRANCHING)
1582          if (dyn_debug >= 1) {
1583            branchingTree.checkTree();
1584          }
1585#endif
1586          ++cnt;
1587        }
1588        if (cnt > 0) {
1589          model.resolve();
1590          // This is horribly looking... Get rid of it when properly debugged...
1591          assert(lpres.isAbandoned == model.isAbandoned());
1592          assert(lpres.isDualObjectiveLimitReached == model.isDualObjectiveLimitReached());
1593          assert(lpres.isDualObjectiveLimitReached ||
1594                 (lpres.isProvenOptimal == model.isProvenOptimal()));
1595          assert(lpres.isDualObjectiveLimitReached ||
1596                 (lpres.isProvenPrimalInfeasible == model.isProvenPrimalInfeasible()));
1597        }
1598           
1599#if defined(DEBUG_DYNAMIC_BRANCHING)
1600        if (dyn_debug >= 10) {
1601          if (cnt > 0) {
1602            std::string tree1 = getTree(branchingTree);
1603            printf("=====================================\n");
1604            printf("It can move node %i up. way_: 0x%x   brvar: %i\n",
1605                   node.node_id_, node.way_, node.variable_);
1606            printTree(tree0, cnt+10);
1607            printf("=====================================\n");
1608            printf("Finished moving the node up by %i levels.\n", cnt);
1609            printTree(tree1, cnt+10);
1610            printf("=====================================\n");
1611          }
1612        }
1613#endif
1614        if ((numberNodes%1000)==0) 
1615          printf("%d nodes, tree size %d\n",
1616                 numberNodes,branchingTree.size());
1617        if (CoinCpuTime()-time1>3600.0) {
1618          printf("stopping after 3600 seconds\n");
1619          exit(77);
1620        }
1621        DBNodeSimple newNode(model,numberIntegers,which,ws);
1622        // something extra may have been fixed by strong branching
1623        // if so go round again
1624        while (newNode.variable_==numberIntegers) {
1625          model.resolve();
1626          newNode = DBNodeSimple(model,numberIntegers,which,model.getWarmStart());
1627          newNode.strong_branching_fixed_vars_ = true;
1628        }
1629        newNode.reduced_cost_fixed_vars_ = did_reduced_cost_fixing_for_child;
1630        if (newNode.objectiveValue_<1.0e100) {
1631          newNode.parent_ = kNode;
1632          // push on stack
1633          branchingTree.push_back(newNode);
1634#if defined(DEBUG_DYNAMIC_BRANCHING)
1635          if (dyn_debug >= 1000) {
1636            printf("Added current child %d %d\n",branchingTree.size(),branchingTree.size_);
1637          }
1638#endif
1639          if (branchingTree.nodes_[kNode].workingOnDownChild()) {
1640            branchingTree.nodes_[kNode].child_down_ = branchingTree.last_;
1641          } else {
1642            branchingTree.nodes_[kNode].child_up_ = branchingTree.last_;
1643          }
1644          if (newNode.variable_>=0) {
1645            assert (fabs(newNode.value_-floor(newNode.value_+0.5))>1.0e-6);
1646          }
1647#if 0
1648          else {
1649            // integer solution - save
1650            bestNode = node;
1651            // set cutoff (hard coded tolerance)
1652            model.setDblParam(OsiDualObjectiveLimit,(bestNode.objectiveValue_-1.0e-5)*direction);
1653            std::cout<<"Integer solution of "
1654                     <<bestNode.objectiveValue_
1655                     <<" found after "<<numberIterations
1656                     <<" iterations and "<<numberNodes<<" nodes"
1657                     <<std::endl;
1658          }
1659#endif
1660        }
1661#if defined(DEBUG_DYNAMIC_BRANCHING)
1662        if (dyn_debug >= 1) {
1663          branchingTree.checkTree();
1664        }
1665#endif
1666      }
1667    }
1668    ////// End main while of branch and bound
1669    std::cout<<"Search took "
1670             <<numberIterations
1671             <<" iterations and "<<numberNodes<<" nodes"
1672             <<std::endl;
1673    if (bestNode.numberIntegers_) {
1674      // we have a solution restore
1675      // do bounds
1676      for (i=0;i<numberIntegers;i++) {
1677        iColumn=which[i];
1678        model.setColBounds( iColumn,bestNode.lower_[i],bestNode.upper_[i]);
1679      }
1680      // move basis
1681      model.setWarmStart(bestNode.basis_);
1682      // set cutoff so will be good (hard coded tolerance)
1683      model.setDblParam(OsiDualObjectiveLimit,(bestNode.objectiveValue_+1.0e-5)*direction);
1684      model.resolve();
1685    } else {
1686      OsiClpSolverInterface* clp =
1687        dynamic_cast<OsiClpSolverInterface*>(&model);
1688      if (clp) {
1689        ClpSimplex* modelPtr_ = clp->getModelPtr();
1690        modelPtr_->setProblemStatus(1);
1691      }
1692    }
1693    delete [] which;
1694    delete [] originalLower;
1695    delete [] originalUpper;
1696    delete [] relaxedLower;
1697    delete [] relaxedUpper;
1698  } else {
1699    std::cout<<"The LP relaxation is infeasible"
1700             <<std::endl;
1701    OsiClpSolverInterface* clp =
1702      dynamic_cast<OsiClpSolverInterface*>(&model);
1703    if (clp) {
1704      ClpSimplex* modelPtr_ = clp->getModelPtr();
1705      modelPtr_->setProblemStatus(1);
1706    }
1707    //throw CoinError("The LP relaxation is infeasible or too expensive",
1708    //"branchAndBound", "OsiClpSolverInterface");
1709  }
1710}
1711
1712
1713
1714int main(int argc, char* argv[])
1715{
1716  OsiClpSolverInterface model;
1717  model.readMps(argv[1]);
1718  branchAndBound(model);
1719  return 0;
1720}
Note: See TracBrowser for help on using the repository browser.