source: branches/dynamicbranching/dynamicbranching.cpp @ 972

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