source: branches/dynamicbranching/dynamicbranching.cpp @ 971

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

change how way is set. its value is now the currently processed child (instead of pointing to what will be processed next)

File size: 31.7 KB
Line 
1#include "CoinTime.hpp"
2#include "OsiClpSolverInterface.hpp"
3
4
5// below needed for pathetic branch and bound code
6#include <vector>
7#include <map>
8
9// Trivial class for Branch and Bound
10
11class DBNodeSimple  {
12public:
13  enum DBNodeWay {
14    DOWN_UP__NOTHING_DONE,
15    DOWN_UP__DOWN_DONE,
16    DOWN_UP__BOTH_DONE,
17    UP_DOWN__NOTHING_DONE,
18    UP_DOWN__UP_DONE,
19    UP_DOWN__BOTH_DONE,
20    WAY_UNSET
21  };
22 
23public:
24   
25  // Default Constructor
26  DBNodeSimple ();
27
28  // Constructor from current state (and list of integers)
29  // Also chooses branching variable (if none set to -1)
30  DBNodeSimple (OsiSolverInterface &model,
31                 int numberIntegers, int * integer,
32                 CoinWarmStart * basis);
33  void gutsOfConstructor (OsiSolverInterface &model,
34                 int numberIntegers, int * integer,
35                 CoinWarmStart * basis);
36  // Copy constructor
37  DBNodeSimple ( const DBNodeSimple &);
38   
39  // Assignment operator
40  DBNodeSimple & operator=( const DBNodeSimple& rhs);
41
42  // Destructor
43  ~DBNodeSimple ();
44  // Work of destructor
45  void gutsOfDestructor();
46  // Extension - true if other extension of this
47  bool extension(const DBNodeSimple & other,
48                 const double * originalLower,
49                 const double * originalUpper) const;
50  inline void incrementDescendants()
51  { descendants_++;}
52  // Public data
53  // Basis (should use tree, but not as wasteful as bounds!)
54  CoinWarmStart * basis_;
55  // Objective value (COIN_DBL_MAX) if spare node
56  double objectiveValue_;
57  // Branching variable (0 is first integer)
58  int variable_;
59  // Way to branch: see enum DBNodeWay
60  int way_;
61  // Number of integers (for length of arrays)
62  int numberIntegers_;
63  // Current value
64  double value_;
65  // Number of descendant nodes (so 2 is in interior)
66  int descendants_;
67  // Parent
68  int parent_;
69  // Left child
70  int child0_;
71  // Right child
72  int child1_;
73  // Previous in chain
74  int previous_;
75  // Next in chain
76  int next_;
77  // Now I must use tree
78  // Bounds stored in full (for integers)
79  int * lower_;
80  int * upper_;
81};
82
83
84DBNodeSimple::DBNodeSimple() :
85  basis_(NULL),
86  objectiveValue_(COIN_DBL_MAX),
87  variable_(-100),
88  way_(WAY_UNSET),
89  numberIntegers_(0),
90  value_(0.5),
91  descendants_(-1),
92  parent_(-1),
93  child0_(-1),
94  child1_(-1),
95  previous_(-1),
96  next_(-1),
97  lower_(NULL),
98  upper_(NULL)
99{
100}
101DBNodeSimple::DBNodeSimple(OsiSolverInterface & model,
102                 int numberIntegers, int * integer,CoinWarmStart * basis)
103{
104  gutsOfConstructor(model,numberIntegers,integer,basis);
105}
106void
107DBNodeSimple::gutsOfConstructor(OsiSolverInterface & model,
108                 int numberIntegers, int * integer,CoinWarmStart * basis)
109{
110  basis_ = basis;
111  variable_=-1;
112  way_=WAY_UNSET;
113  numberIntegers_=numberIntegers;
114  value_=0.0;
115  descendants_ = 0;
116  parent_ = -1;
117  child0_ = -1;
118  child1_ = -1;
119  previous_ = -1;
120  next_ = -1;
121  if (model.isProvenOptimal()&&!model.isDualObjectiveLimitReached()) {
122    objectiveValue_ = model.getObjSense()*model.getObjValue();
123  } else {
124    objectiveValue_ = 1.0e100;
125    lower_ = NULL;
126    upper_ = NULL;
127    return; // node cutoff
128  }
129  lower_ = new int [numberIntegers_];
130  upper_ = new int [numberIntegers_];
131  assert (upper_!=NULL);
132  const double * lower = model.getColLower();
133  const double * upper = model.getColUpper();
134  const double * solution = model.getColSolution();
135  int i;
136  // Hard coded integer tolerance
137#define INTEGER_TOLERANCE 1.0e-6
138  ///////// Start of Strong branching code - can be ignored
139  // Number of strong branching candidates
140#define STRONG_BRANCHING 5
141#ifdef STRONG_BRANCHING
142  double upMovement[STRONG_BRANCHING];
143  double downMovement[STRONG_BRANCHING];
144  double solutionValue[STRONG_BRANCHING];
145  int chosen[STRONG_BRANCHING];
146  int iSmallest=0;
147  // initialize distance from integer
148  for (i=0;i<STRONG_BRANCHING;i++) {
149    upMovement[i]=0.0;
150    chosen[i]=-1;
151  }
152  variable_=-1;
153  // This has hard coded integer tolerance
154  double mostAway=INTEGER_TOLERANCE;
155  int numberAway=0;
156  for (i=0;i<numberIntegers;i++) {
157    int iColumn = integer[i];
158    lower_[i]=(int)lower[iColumn];
159    upper_[i]=(int)upper[iColumn];
160    double value = solution[iColumn];
161    value = max(value,(double) lower_[i]);
162    value = min(value,(double) upper_[i]);
163    double nearest = floor(value+0.5);
164    if (fabs(value-nearest)>INTEGER_TOLERANCE)
165      numberAway++;
166    if (fabs(value-nearest)>mostAway) {
167      double away = fabs(value-nearest);
168      if (away>upMovement[iSmallest]) {
169        //add to list
170        upMovement[iSmallest]=away;
171        solutionValue[iSmallest]=value;
172        chosen[iSmallest]=i;
173        int j;
174        iSmallest=-1;
175        double smallest = 1.0;
176        for (j=0;j<STRONG_BRANCHING;j++) {
177          if (upMovement[j]<smallest) {
178            smallest=upMovement[j];
179            iSmallest=j;
180          }
181        }
182      }
183    }
184  }
185  int numberStrong=0;
186  for (i=0;i<STRONG_BRANCHING;i++) {
187    if (chosen[i]>=0) { 
188      numberStrong ++;
189      variable_ = chosen[i];
190    }
191  }
192  // out strong branching if bit set
193  OsiClpSolverInterface* clp =
194    dynamic_cast<OsiClpSolverInterface*>(&model);
195  if (clp&&(clp->specialOptions()&16)!=0&&numberStrong>1) {
196    int j;
197    int iBest=-1;
198    double best = 0.0;
199    for (j=0;j<STRONG_BRANCHING;j++) {
200      if (upMovement[j]>best) {
201        best=upMovement[j];
202        iBest=j;
203      }
204    }
205    numberStrong=1;
206    variable_=chosen[iBest];
207  }
208  if (numberStrong==1) {
209    // just one - makes it easy
210    int iColumn = integer[variable_];
211    double value = solution[iColumn];
212    value = max(value,(double) lower_[variable_]);
213    value = min(value,(double) upper_[variable_]);
214    double nearest = floor(value+0.5);
215    value_=value;
216    if (value<=nearest)
217      way_=UP_DOWN__NOTHING_DONE; // up
218    else
219      way_=DOWN_UP__NOTHING_DONE; // down
220  } else if (numberStrong) {
221    // more than one - choose
222    bool chooseOne=true;
223    model.markHotStart();
224    for (i=0;i<STRONG_BRANCHING;i++) {
225      int iInt = chosen[i];
226      if (iInt>=0) {
227        int iColumn = integer[iInt];
228        double value = solutionValue[i]; // value of variable in original
229        double objectiveChange;
230        value = max(value,(double) lower_[iInt]);
231        value = min(value,(double) upper_[iInt]);
232
233        // try down
234
235        model.setColUpper(iColumn,floor(value));
236        model.solveFromHotStart();
237        model.setColUpper(iColumn,upper_[iInt]);
238        if (model.isProvenOptimal()&&!model.isDualObjectiveLimitReached()) {
239          objectiveChange = model.getObjSense()*model.getObjValue()
240            - objectiveValue_;
241        } else {
242          objectiveChange = 1.0e100;
243        }
244        assert (objectiveChange>-1.0e-5);
245        objectiveChange = CoinMax(objectiveChange,0.0);
246        downMovement[i]=objectiveChange;
247
248        // try up
249
250        model.setColLower(iColumn,ceil(value));
251        model.solveFromHotStart();
252        model.setColLower(iColumn,lower_[iInt]);
253        if (model.isProvenOptimal()&&!model.isDualObjectiveLimitReached()) {
254          objectiveChange = model.getObjSense()*model.getObjValue()
255            - objectiveValue_;
256        } else {
257          objectiveChange = 1.0e100;
258        }
259        assert (objectiveChange>-1.0e-5);
260        objectiveChange = CoinMax(objectiveChange,0.0);
261        upMovement[i]=objectiveChange;
262       
263        /* Possibilities are:
264           Both sides feasible - store
265           Neither side feasible - set objective high and exit
266           One side feasible - change bounds and resolve
267        */
268        bool solveAgain=false;
269        if (upMovement[i]<1.0e100) {
270          if(downMovement[i]<1.0e100) {
271            // feasible - no action
272          } else {
273            // up feasible, down infeasible
274            solveAgain = true;
275            model.setColLower(iColumn,ceil(value));
276          }
277        } else {
278          if(downMovement[i]<1.0e100) {
279            // down feasible, up infeasible
280            solveAgain = true;
281            model.setColUpper(iColumn,floor(value));
282          } else {
283            // neither side feasible
284            objectiveValue_=1.0e100;
285            chooseOne=false;
286            break;
287          }
288        }
289        if (solveAgain) {
290          // need to solve problem again - signal this
291          variable_ = numberIntegers;
292          chooseOne=false;
293          break;
294        }
295      }
296    }
297    if (chooseOne) {
298      // choose the one that makes most difference both ways
299      double best = -1.0;
300      double best2 = -1.0;
301      for (i=0;i<STRONG_BRANCHING;i++) {
302        int iInt = chosen[i];
303        if (iInt>=0) {
304          //std::cout<<"Strong branching on "
305          //   <<i<<""<<iInt<<" down "<<downMovement[i]
306          //   <<" up "<<upMovement[i]
307          //   <<" value "<<solutionValue[i]
308          //   <<std::endl;
309          bool better = false;
310          if (min(upMovement[i],downMovement[i])>best) {
311            // smaller is better
312            better=true;
313          } else if (min(upMovement[i],downMovement[i])>best-1.0e-5) {
314            if (max(upMovement[i],downMovement[i])>best2+1.0e-5) {
315              // smaller is about same, but larger is better
316              better=true;
317            }
318          }
319          if (better) {
320            best = min(upMovement[i],downMovement[i]);
321            best2 = max(upMovement[i],downMovement[i]);
322            variable_ = iInt;
323            double value = solutionValue[i];
324            value = max(value,(double) lower_[variable_]);
325            value = min(value,(double) upper_[variable_]);
326            value_=value;
327            if (upMovement[i]<=downMovement[i])
328              way_=UP_DOWN__NOTHING_DONE; // up
329            else
330              way_=DOWN_UP__NOTHING_DONE; // down
331          }
332        }
333      }
334    }
335    // Delete the snapshot
336    model.unmarkHotStart();
337  }
338  ////// End of Strong branching
339#else
340  variable_=-1;
341  // This has hard coded integer tolerance
342  double mostAway=INTEGER_TOLERANCE;
343  int numberAway=0;
344  for (i=0;i<numberIntegers;i++) {
345    int iColumn = integer[i];
346    lower_[i]=(int)lower[iColumn];
347    upper_[i]=(int)upper[iColumn];
348    double value = solution[iColumn];
349    value = max(value,(double) lower_[i]);
350    value = min(value,(double) upper_[i]);
351    double nearest = floor(value+0.5);
352    if (fabs(value-nearest)>INTEGER_TOLERANCE)
353      numberAway++;
354    if (fabs(value-nearest)>mostAway) {
355      mostAway=fabs(value-nearest);
356      variable_=i;
357      value_=value;
358      if (value<=nearest)
359        way_=UP_DOWN__NOTHING_DONE; // up
360      else
361        way_=DOWN_UP__NOTHING_DONE; // down
362    }
363  }
364#endif
365}
366
367DBNodeSimple::DBNodeSimple(const DBNodeSimple & rhs) 
368{
369  if (rhs.basis_)
370    basis_=rhs.basis_->clone();
371  else
372    basis_ = NULL;
373  objectiveValue_=rhs.objectiveValue_;
374  variable_=rhs.variable_;
375  way_=rhs.way_;
376  numberIntegers_=rhs.numberIntegers_;
377  value_=rhs.value_;
378  descendants_ = rhs.descendants_;
379  parent_ = rhs.parent_;
380  child0_ = rhs.child0_;
381  child1_ = rhs.child1_;
382  previous_ = rhs.previous_;
383  next_ = rhs.next_;
384  lower_=NULL;
385  upper_=NULL;
386  if (rhs.lower_!=NULL) {
387    lower_ = new int [numberIntegers_];
388    upper_ = new int [numberIntegers_];
389    assert (upper_!=NULL);
390    memcpy(lower_,rhs.lower_,numberIntegers_*sizeof(int));
391    memcpy(upper_,rhs.upper_,numberIntegers_*sizeof(int));
392  }
393}
394
395DBNodeSimple &
396DBNodeSimple::operator=(const DBNodeSimple & rhs)
397{
398  if (this != &rhs) {
399    gutsOfDestructor();
400    if (rhs.basis_)
401      basis_=rhs.basis_->clone();
402    objectiveValue_=rhs.objectiveValue_;
403    variable_=rhs.variable_;
404    way_=rhs.way_;
405    numberIntegers_=rhs.numberIntegers_;
406    value_=rhs.value_;
407    descendants_ = rhs.descendants_;
408    parent_ = rhs.parent_;
409    child0_ = rhs.child0_;
410    child1_ = rhs.child1_;
411    previous_ = rhs.previous_;
412    next_ = rhs.next_;
413    if (rhs.lower_!=NULL) {
414      lower_ = new int [numberIntegers_];
415      upper_ = new int [numberIntegers_];
416      assert (upper_!=NULL);
417      memcpy(lower_,rhs.lower_,numberIntegers_*sizeof(int));
418      memcpy(upper_,rhs.upper_,numberIntegers_*sizeof(int));
419    }
420  }
421  return *this;
422}
423
424
425DBNodeSimple::~DBNodeSimple ()
426{
427  gutsOfDestructor();
428}
429// Work of destructor
430void 
431DBNodeSimple::gutsOfDestructor()
432{
433  delete [] lower_;
434  delete [] upper_;
435  delete basis_;
436  lower_ = NULL;
437  upper_ = NULL;
438  basis_ = NULL;
439  objectiveValue_ = COIN_DBL_MAX;
440}
441// Extension - true if other extension of this
442bool 
443DBNodeSimple::extension(const DBNodeSimple & other,
444                         const double * originalLower,
445                         const double * originalUpper) const
446{
447  bool ok=true;
448  for (int i=0;i<numberIntegers_;i++) {
449    if (upper_[i]<originalUpper[i]||
450        lower_[i]>originalLower[i]) {
451      if (other.upper_[i]>upper_[i]||
452          other.lower_[i]<lower_[i]) {
453        ok=false;
454        break;
455      }
456    }
457  }
458  return ok;
459}
460
461#include <vector>
462#define FUNNY_BRANCHING 1
463#define FUNNY_TREE
464#ifndef FUNNY_TREE
465// Vector of DBNodeSimples
466typedef std::vector<DBNodeSimple>    DBVectorNode;
467#else
468// Must code up by hand
469class DBVectorNode  {
470 
471public:
472   
473  // Default Constructor
474  DBVectorNode ();
475
476  // Copy constructor
477  DBVectorNode ( const DBVectorNode &);
478   
479  // Assignment operator
480  DBVectorNode & operator=( const DBVectorNode& rhs);
481
482  // Destructor
483  ~DBVectorNode ();
484  // Size
485  inline int size() const
486  { return size_-sizeDeferred_;}
487  // Push
488  void push_back(const DBNodeSimple & node);
489  // Last one in (or other criterion)
490  DBNodeSimple back() const;
491  // Get rid of last one
492  void pop_back();
493  // Works out best one
494  int best() const;
495 
496  // Public data
497  // Maximum size
498  int maximumSize_;
499  // Current size
500  int size_;
501  // Number still hanging around
502  int sizeDeferred_;
503  // First spare
504  int firstSpare_;
505  // First
506  int first_;
507  // Last
508  int last_;
509  // Chosen one
510  mutable int chosen_;
511  // Nodes
512  DBNodeSimple * nodes_;
513};
514
515
516DBVectorNode::DBVectorNode() :
517  maximumSize_(10),
518  size_(0),
519  sizeDeferred_(0),
520  firstSpare_(0),
521  first_(-1),
522  last_(-1)
523{
524  nodes_ = new DBNodeSimple[maximumSize_];
525  for (int i=0;i<maximumSize_;i++) {
526    nodes_[i].previous_=i-1;
527    nodes_[i].next_=i+1;
528  }
529}
530
531DBVectorNode::DBVectorNode(const DBVectorNode & rhs) 
532{ 
533  maximumSize_ = rhs.maximumSize_;
534  size_ = rhs.size_;
535  sizeDeferred_ = rhs.sizeDeferred_;
536  firstSpare_ = rhs.firstSpare_;
537  first_ = rhs.first_;
538  last_ = rhs.last_;
539  nodes_ = new DBNodeSimple[maximumSize_];
540  for (int i=0;i<maximumSize_;i++) {
541    nodes_[i] = rhs.nodes_[i];
542  }
543}
544
545DBVectorNode &
546DBVectorNode::operator=(const DBVectorNode & rhs)
547{
548  if (this != &rhs) {
549    delete [] nodes_;
550    maximumSize_ = rhs.maximumSize_;
551    size_ = rhs.size_;
552    sizeDeferred_ = rhs.sizeDeferred_;
553    firstSpare_ = rhs.firstSpare_;
554    first_ = rhs.first_;
555    last_ = rhs.last_;
556    nodes_ = new DBNodeSimple[maximumSize_];
557    for (int i=0;i<maximumSize_;i++) {
558      nodes_[i] = rhs.nodes_[i];
559    }
560  }
561  return *this;
562}
563
564
565DBVectorNode::~DBVectorNode ()
566{
567  delete [] nodes_;
568}
569// Push
570void 
571DBVectorNode::push_back(const DBNodeSimple & node)
572{
573  if (size_==maximumSize_) {
574    assert (firstSpare_==size_);
575    maximumSize_ = (maximumSize_*3)+10;
576    DBNodeSimple * temp = new DBNodeSimple[maximumSize_];
577    int i;
578    for (i=0;i<size_;i++) {
579      temp[i]=nodes_[i];
580    }
581    delete [] nodes_;
582    nodes_ = temp;
583    //firstSpare_=size_;
584    int last = -1;
585    for ( i=size_;i<maximumSize_;i++) {
586      nodes_[i].previous_=last;
587      nodes_[i].next_=i+1;
588      last = i;
589    }
590  }
591  assert (firstSpare_<maximumSize_);
592  assert (nodes_[firstSpare_].previous_<0);
593  int next = nodes_[firstSpare_].next_;
594  nodes_[firstSpare_]=node;
595  if (last_>=0) {
596    assert (nodes_[last_].next_==-1);
597    nodes_[last_].next_=firstSpare_;
598  }
599  nodes_[firstSpare_].previous_=last_;
600  nodes_[firstSpare_].next_=-1;
601  if (last_==-1) {
602    assert (first_==-1);
603    first_ = firstSpare_;
604  }
605  last_=firstSpare_;
606  if (next>=0&&next<maximumSize_) {
607    firstSpare_ = next;
608    nodes_[firstSpare_].previous_=-1;
609  } else {
610    firstSpare_=maximumSize_;
611  }
612  chosen_ = -1;
613  //best();
614  size_++;
615  assert (node.descendants_<=2);
616  if (node.descendants_==2)
617    sizeDeferred_++;
618}
619// Works out best one
620int 
621DBVectorNode::best() const
622{
623  // can modify
624  chosen_=-1;
625  if (chosen_<0) {
626    chosen_=last_;
627#if FUNNY_BRANCHING
628    while (nodes_[chosen_].descendants_==2) {
629      chosen_ = nodes_[chosen_].previous_;
630      assert (chosen_>=0);
631    }
632#endif
633  }
634  return chosen_;
635}
636// Last one in (or other criterion)
637DBNodeSimple
638DBVectorNode::back() const
639{
640  assert (last_>=0);
641  return nodes_[best()];
642}
643// Get rid of last one
644void 
645DBVectorNode::pop_back()
646{
647  // Temporary until more sophisticated
648  //assert (last_==chosen_);
649  if (nodes_[chosen_].descendants_==2)
650    sizeDeferred_--;
651  int previous = nodes_[chosen_].previous_;
652  int next = nodes_[chosen_].next_;
653  nodes_[chosen_].gutsOfDestructor();
654  if (previous>=0) {
655    nodes_[previous].next_=next;
656  } else {
657    first_ = next;
658  }
659  if (next>=0) {
660    nodes_[next].previous_ = previous;
661  } else {
662    last_ = previous;
663  }
664  nodes_[chosen_].previous_=-1;
665  if (firstSpare_>=0) {
666    nodes_[chosen_].next_ = firstSpare_;
667  } else {
668    nodes_[chosen_].next_ = -1;
669  }
670  firstSpare_ = chosen_;
671  chosen_ = -1;
672  assert (size_>0);
673  size_--;
674}
675#endif
676
677bool moveNodes(OsiSolverInterface & model,
678               DBVectorNode & branchingTree,
679               int kNode)
680{
681
682  DBNodeSimple & node = branchingTree[kNode];
683  DBNodeSimple & grandParent = branchingTree[node.parent_];
684  int grandParentVariable = grandParent.variable_;
685  // check if branching constraint of grandparent is tight
686  bool canMoveNodes = checkGrandparent();
687
688  if(!canMoveNodes)
689    return false;
690
691  node.parent_ = grandParent.parent_;
692  grandParent.parent_ = kNode;
693  // change bounds of grandParent
694  grandParent.lower_[
695
696
697}
698
699// Invoke solver's built-in enumeration algorithm
700void 
701branchAndBound(OsiSolverInterface & model) {
702  double time1 = CoinCpuTime();
703  // solve LP
704  model.initialSolve();
705  int funnyBranching=FUNNY_BRANCHING;
706
707  if (model.isProvenOptimal()&&!model.isDualObjectiveLimitReached()) {
708    // Continuous is feasible - find integers
709    int numberIntegers=0;
710    int numberColumns = model.getNumCols();
711    int iColumn;
712    int i;
713    for (iColumn=0;iColumn<numberColumns;iColumn++) {
714      if( model.isInteger(iColumn))
715        numberIntegers++;
716    }
717    if (!numberIntegers) {
718      std::cout<<"No integer variables"
719               <<std::endl;
720      return;
721    }
722    int * which = new int[numberIntegers]; // which variables are integer
723    // original bounds
724    int * originalLower = new int[numberIntegers];
725    int * originalUpper = new int[numberIntegers];
726    int * relaxedLower = new int[numberIntegers];
727    int * relaxedUpper = new int[numberIntegers];
728    {
729      const double * lower = model.getColLower();
730      const double * upper = model.getColUpper();
731      numberIntegers=0;
732      for (iColumn=0;iColumn<numberColumns;iColumn++) {
733        if( model.isInteger(iColumn)) {
734          originalLower[numberIntegers]=(int) lower[iColumn];
735          originalUpper[numberIntegers]=(int) upper[iColumn];
736          which[numberIntegers++]=iColumn;
737        }
738      }
739    }
740    double direction = model.getObjSense();
741    // empty tree
742    DBVectorNode branchingTree;
743   
744    // Add continuous to it;
745    DBNodeSimple rootNode(model,numberIntegers,which,model.getWarmStart());
746    // something extra may have been fixed by strong branching
747    // if so go round again
748    while (rootNode.variable_==numberIntegers) {
749      model.resolve();
750      rootNode = DBNodeSimple(model,numberIntegers,which,model.getWarmStart());
751    }
752    if (rootNode.objectiveValue_<1.0e100) {
753      // push on stack
754      branchingTree.push_back(rootNode);
755    }
756   
757    // For printing totals
758    int numberIterations=0;
759    int numberNodes =0;
760    int nRedundantUp=0;
761    int nRedundantDown=0;
762    int nRedundantUp2=0;
763    int nRedundantDown2=0;
764    DBNodeSimple bestNode;
765    ////// Start main while of branch and bound
766    // while until nothing on stack
767    while (branchingTree.size()) {
768      // last node
769      DBNodeSimple node = branchingTree.back();
770      int kNode = branchingTree.chosen_;
771      branchingTree.pop_back();
772      assert (node.descendants_<2);
773      numberNodes++;
774      if (node.variable_>=0) {
775        // branch - do bounds
776        for (i=0;i<numberIntegers;i++) {
777          iColumn=which[i];
778          model.setColBounds( iColumn,node.lower_[i],node.upper_[i]);
779        }
780        // move basis
781        model.setWarmStart(node.basis_);
782        // do branching variable
783        node.incrementDescendants();
784        bool down_branch = true;
785        switch (node.way_) {
786        case WAY_UNSET:
787        case DOWN_UP__BOTH_DONE:
788        case UP_DOWN__BOTH_DONE:
789          abort();
790        case DOWN_UP__NOTHING_DONE:
791          node.way_ = DOWN_UP__DOWN_DONE;
792          break;
793        case DOWN_UP__DOWN_DONE:
794          node.way_ = DOWN_UP__BOTH_DONE:
795          down_branch = false;
796          break;
797        case UP_DOWN__NOTHING_DONE:
798          node.way_ = UP_DOWN__UP_DONE:
799          down_branch = false;
800          break;
801        case UP_DOWN__UP_DONE:
802          node.way_ = UP_DOWN__BOTH_DONE;
803          break;
804        }
805        if (down_branch) {
806          model.setColUpper(which[node.variable_],floor(node.value_));
807        } else {
808          model.setColLower(which[node.variable_],ceil(node.value_));
809        }
810        // put back on tree anyway regardless whether any processing is left
811        // to be done. We want the whole tree all the time.
812        branchingTree.push_back(node);
813       
814        // solve
815        model.resolve();
816        CoinWarmStart * ws = model.getWarmStart();
817        const CoinWarmStartBasis* wsb =
818          dynamic_cast<const CoinWarmStartBasis*>(ws);
819        assert (wsb!=NULL); // make sure not volume
820        numberIterations += model.getIterationCount();
821        // fix on reduced costs
822        int nFixed0=0,nFixed1=0;
823        double cutoff;
824        model.getDblParam(OsiDualObjectiveLimit,cutoff);
825        double gap=(cutoff-model.getObjValue())*direction+1.0e-4;
826        //        double gap=(cutoff-modelPtr_->objectiveValue())*direction+1.0e-4;
827        if (gap<1.0e10&&model.isProvenOptimal()&&!model.isDualObjectiveLimitReached()) {
828          const double * dj = model.getReducedCost();
829          const double * lower = model.getColLower();
830          const double * upper = model.getColUpper();
831          for (i=0;i<numberIntegers;i++) {
832            iColumn=which[i];
833            if (upper[iColumn]>lower[iColumn]) {
834              double djValue = dj[iColumn]*direction;
835              if (wsb->getStructStatus(iColumn)==CoinWarmStartBasis::atLowerBound&&
836                  djValue>gap) {
837                nFixed0++;
838                model.setColUpper(iColumn,lower[iColumn]);
839              } else if (wsb->getStructStatus(iColumn)==CoinWarmStartBasis::atUpperBound&&
840                         -djValue>gap) {
841                nFixed1++;
842                model.setColLower(iColumn,upper[iColumn]);
843              }
844            }
845          }
846          //if (nFixed0+nFixed1)
847          //printf("%d fixed to lower, %d fixed to upper\n",nFixed0,nFixed1);
848        }
849        if (!model.isIterationLimitReached()) {
850          if (model.isProvenOptimal()&&!model.isDualObjectiveLimitReached()) {
851#if FUNNY_BRANCHING
852            // See if branched variable off bounds
853            const double * dj = model.getReducedCost();
854            const double * lower = model.getColLower();
855            const double * upper = model.getColUpper();
856            const double * solution = model.getColSolution();
857            // Better to use "natural" value - need flag to say fixed
858            for (i=0;i<numberIntegers;i++) {
859              iColumn=which[i];
860              relaxedLower[i]=originalLower[i];
861              relaxedUpper[i]=originalUpper[i];
862              double djValue = dj[iColumn]*direction;
863              if (djValue>1.0e-6) {
864                // wants to go down
865                if (lower[iColumn]>originalLower[i]) {
866                  // Lower bound active
867                  relaxedLower[i]=(int) lower[iColumn];
868                }
869                if (upper[iColumn]<originalUpper[i]) {
870                  // Upper bound NOT active
871                }
872              } else if (djValue<-1.0e-6) {
873                // wants to go up
874                if (lower[iColumn]>originalLower[i]) {
875                  // Lower bound NOT active
876                }
877                if (upper[iColumn]<originalUpper[i]) {
878                  // Upper bound active
879                  relaxedUpper[i]=(int) upper[iColumn];
880                }
881              }
882            }
883            // See if can do anything
884            {
885              /*
886                If kNode is on second branch then
887                a) If other feasible could free up as well
888                b) If other infeasible could do something clever.
889                For now - we have to give up
890              */
891              int jNode=branchingTree.nodes_[kNode].parent_;
892              bool canDelete = (branchingTree.nodes_[kNode].descendants_<2);
893              while (jNode>=0) {
894                DBNodeSimple & node = branchingTree.nodes_[jNode];
895                int next = node.parent_;
896                if (node.descendants_<2) {
897                  int variable = node.variable_;
898                  iColumn=which[variable];
899                  double value = node.value_;
900                  double djValue = dj[iColumn]*direction;
901                  assert (node.way_==2||node.way_==-2);
902                  // we don't know which branch it was - look at current bounds
903                  if (upper[iColumn]<value&&node.lower_[variable]<upper[iColumn]) {
904                    // must have been down branch
905                    if (djValue>1.0e-3||solution[iColumn]<upper[iColumn]-1.0e-5) {
906                      if (canDelete) {
907                        nRedundantDown++;
908#if 1
909                        printf("%d redundant branch down with value %g current upper %g solution %g dj %g\n",
910                               variable,node.value_,upper[iColumn],solution[iColumn],djValue);
911#endif
912                        node.descendants_=2; // ignore
913                        branchingTree.sizeDeferred_++;
914                        int newUpper = originalUpper[variable];
915                        if (next>=0) {
916                          DBNodeSimple & node2 = branchingTree.nodes_[next];
917                          newUpper = node2.upper_[variable];
918                        }
919                        if (branchingTree.nodes_[jNode].parent_!=next)
920                          assert (newUpper>upper[iColumn]);
921                        model.setColUpper(iColumn,newUpper);
922                        int kNode2=next;
923                        int jNode2=branchingTree.nodes_[kNode].parent_;
924                        assert (newUpper>branchingTree.nodes_[kNode].upper_[variable]);
925                        branchingTree.nodes_[kNode].upper_[variable]= newUpper;
926                        while (jNode2!=kNode2) {
927                          DBNodeSimple & node2 = branchingTree.nodes_[jNode2];
928                          int next = node2.parent_;
929                          if (next!=kNode2)
930                            assert (newUpper>node2.upper_[variable]);
931                          node2.upper_[variable]= newUpper;
932                          jNode2=next;
933                        }
934                      } else {
935                        // can't delete but can add other way to jNode
936                        nRedundantDown2++;
937                        DBNodeSimple & node2 = branchingTree.nodes_[kNode];
938                        assert (node2.way_==2||node2.way_==-2);
939                        double value2 = node2.value_;
940                        int variable2 = node2.variable_;
941                        int iColumn2 = which[variable2];
942                        if (variable != variable2) {
943                          if (node2.way_==2&&upper[iColumn2]<value2) {
944                            // must have been down branch which was done - carry over
945                            int newUpper = (int) floor(value2);
946                            assert (newUpper<node.upper_[variable2]);
947                            node.upper_[variable2]=newUpper;
948                          } else if (node2.way_==-2&&lower[iColumn2]>value2) {
949                            // must have been up branch which was done - carry over
950                            int newLower = (int) ceil(value2);
951                            assert (newLower>node.lower_[variable2]);
952                            node.lower_[variable2]=newLower;
953                          }
954                          if (node.lower_[variable2]>node.upper_[variable2]) {
955                            // infeasible
956                            node.descendants_=2; // ignore
957                            branchingTree.sizeDeferred_++;
958                          }
959                        }
960                      }
961                      break;
962                    }           
963                    // we don't know which branch it was - look at current bounds
964                  } else if (lower[iColumn]>value&&node.upper_[variable]>lower[iColumn]) {
965                    // must have been up branch
966                    if (djValue<-1.0e-3||solution[iColumn]>lower[iColumn]+1.0e-5) {
967                      if (canDelete) {
968                        nRedundantUp++;
969#if 1
970                        printf("%d redundant branch up with value %g current lower %g solution %g dj %g\n",
971                               variable,node.value_,lower[iColumn],solution[iColumn],djValue);
972#endif
973                        node.descendants_=2; // ignore
974                        branchingTree.sizeDeferred_++;
975                        int newLower = originalLower[variable];
976                        if (next>=0) {
977                          DBNodeSimple & node2 = branchingTree.nodes_[next];
978                          newLower = node2.lower_[variable];
979                        }
980                        if (branchingTree.nodes_[jNode].parent_!=next)
981                          assert (newLower<lower[iColumn]);
982                        model.setColLower(iColumn,newLower);
983                        int kNode2=next;
984                        int jNode2=branchingTree.nodes_[kNode].parent_;
985                        assert (newLower<branchingTree.nodes_[kNode].lower_[variable]);
986                        branchingTree.nodes_[kNode].lower_[variable]= newLower;
987                        while (jNode2!=kNode2) {
988                          DBNodeSimple & node2 = branchingTree.nodes_[jNode2];
989                          int next = node2.parent_;
990                          if (next!=kNode2)
991                            assert (newLower<node2.lower_[variable]);
992                          node2.lower_[variable]=newLower;
993                          jNode2=next;
994                        }
995                      } else {
996                        // can't delete but can add other way to jNode
997                        nRedundantUp2++;
998                        DBNodeSimple & node2 = branchingTree.nodes_[kNode];
999                        assert (node2.way_==2||node2.way_==-2);
1000                        double value2 = node2.value_;
1001                        int variable2 = node2.variable_;
1002                        int iColumn2 = which[variable2];
1003                        if (variable != variable2) {
1004                          if (node2.way_==2&&upper[iColumn2]<value2) {
1005                            // must have been down branch which was done - carry over
1006                            int newUpper = (int) floor(value2);
1007                            assert (newUpper<node.upper_[variable2]);
1008                            node.upper_[variable2]=newUpper;
1009                          } else if (node2.way_==-2&&lower[iColumn2]>value2) {
1010                            // must have been up branch which was done - carry over
1011                            int newLower = (int) ceil(value2);
1012                            assert (newLower>node.lower_[variable2]);
1013                            node.lower_[variable2]=newLower;
1014                          }
1015                          if (node.lower_[variable2]>node.upper_[variable2]) {
1016                            // infeasible
1017                            node.descendants_=2; // ignore
1018                            branchingTree.sizeDeferred_++;
1019                          }
1020                        }
1021                      }
1022                      break;
1023                    }           
1024                  }
1025                } else {
1026                  break;
1027                }
1028                jNode=next;
1029              }
1030            }
1031            // solve
1032            //resolve();
1033            //assert(!getIterationCount());
1034            if ((numberNodes%1000)==0) 
1035              printf("%d nodes, redundant down %d (%d) up %d (%d) tree size %d\n",
1036                     numberNodes,nRedundantDown,nRedundantDown2,nRedundantUp,nRedundantUp2,branchingTree.size());
1037#else
1038            if ((numberNodes%1000)==0) 
1039              printf("%d nodes, tree size %d\n",
1040                     numberNodes,branchingTree.size());
1041#endif
1042            if (CoinCpuTime()-time1>3600.0) {
1043              printf("stopping after 3600 seconds\n");
1044              exit(77);
1045            }
1046            DBNodeSimple newNode(model,numberIntegers,which,ws);
1047            // something extra may have been fixed by strong branching
1048            // if so go round again
1049            while (newNode.variable_==numberIntegers) {
1050              model.resolve();
1051              newNode = DBNodeSimple(model,numberIntegers,which,model.getWarmStart());
1052            }
1053            if (newNode.objectiveValue_<1.0e100) {
1054              if (newNode.variable_>=0) 
1055                assert (fabs(newNode.value_-floor(newNode.value_+0.5))>1.0e-6);
1056              newNode.parent_ = kNode;
1057              // push on stack
1058              branchingTree.push_back(newNode);
1059              if(branchingTree.nodes_[kNode].child0_ < 0)
1060                branchingTree.nodes_[kNode].child0_ = branchingTree.last_;
1061              else
1062                branchingTree.nodes_[kNode].child1_ = branchingTree.last_;
1063#if 0
1064              } else {
1065                // integer solution - save
1066                bestNode = node;
1067                // set cutoff (hard coded tolerance)
1068                model.setDblParam(OsiDualObjectiveLimit,(bestNode.objectiveValue_-1.0e-5)*direction);
1069                std::cout<<"Integer solution of "
1070                         <<bestNode.objectiveValue_
1071                         <<" found after "<<numberIterations
1072                         <<" iterations and "<<numberNodes<<" nodes"
1073                 <<std::endl;
1074              }
1075#endif
1076            }
1077          }
1078        } else {
1079          // maximum iterations - exit
1080          std::cout<<"Exiting on maximum iterations"
1081                   <<std::endl;
1082          break;
1083        }
1084      } else {
1085        // integer solution - save
1086        bestNode = node;
1087        // set cutoff (hard coded tolerance)
1088        model.setDblParam(OsiDualObjectiveLimit,(bestNode.objectiveValue_-1.0e-5)*direction);
1089        std::cout<<"Integer solution of "
1090                 <<bestNode.objectiveValue_
1091                 <<" found after "<<numberIterations
1092                 <<" iterations and "<<numberNodes<<" nodes"
1093                 <<std::endl;
1094      }
1095    }
1096    ////// End main while of branch and bound
1097    std::cout<<"Search took "
1098             <<numberIterations
1099             <<" iterations and "<<numberNodes<<" nodes"
1100             <<std::endl;
1101    if (bestNode.numberIntegers_) {
1102      // we have a solution restore
1103      // do bounds
1104      for (i=0;i<numberIntegers;i++) {
1105        iColumn=which[i];
1106        model.setColBounds( iColumn,bestNode.lower_[i],bestNode.upper_[i]);
1107      }
1108      // move basis
1109      model.setWarmStart(bestNode.basis_);
1110      // set cutoff so will be good (hard coded tolerance)
1111      model.setDblParam(OsiDualObjectiveLimit,(bestNode.objectiveValue_+1.0e-5)*direction);
1112      model.resolve();
1113    } else {
1114      OsiClpSolverInterface* clp =
1115        dynamic_cast<OsiClpSolverInterface*>(&model);
1116      if (clp) {
1117        ClpSimplex* modelPtr_ = clp->getModelPtr();
1118        modelPtr_->setProblemStatus(1);
1119      }
1120    }
1121    delete [] which;
1122    delete [] originalLower;
1123    delete [] originalUpper;
1124    delete [] relaxedLower;
1125    delete [] relaxedUpper;
1126  } else {
1127    std::cout<<"The LP relaxation is infeasible"
1128             <<std::endl;
1129    OsiClpSolverInterface* clp =
1130      dynamic_cast<OsiClpSolverInterface*>(&model);
1131    if (clp) {
1132      ClpSimplex* modelPtr_ = clp->getModelPtr();
1133      modelPtr_->setProblemStatus(1);
1134    }
1135    //throw CoinError("The LP relaxation is infeasible or too expensive",
1136    //"branchAndBound", "OsiClpSolverInterface");
1137  }
1138}
1139
1140
1141
1142int main(int argc, char* argv[])
1143{
1144  return 0;
1145}
Note: See TracBrowser for help on using the repository browser.