source: branches/dynamicbranching/dynamicbranching.cpp @ 975

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