source: branches/dynamicbranching/dynamicbranching.cpp @ 970

Last change on this file since 970 was 970, checked in by jpgoncal, 14 years ago

Added function to do move.

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