source: branches/dynamicbranching/dynamicbranching.cpp @ 1422

Last change on this file since 1422 was 1021, checked in by jpgoncal, 11 years ago

Added methods to delete nodes above the current node.

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