source: branches/dynamicbranching/dynamicbranching.cpp @ 981

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

Added print statements for debugging.

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