source: branches/dynamicbranching/dynamicbranching.cpp @ 985

Last change on this file since 985 was 985, checked in by ladanyi, 11 years ago

slowly debugging

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