source: branches/dynamicbranching/dynamicbranching.cpp @ 965

Last change on this file since 965 was 965, checked in by jpgoncal, 12 years ago

Copied code from OsiClpSolverInterface?.

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