source: branches/dynamicbranching/dynamicbranching.cpp @ 976

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

seems to work if FUNNY_BRANCHING is not defined

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    while (nodes_[chosen_].descendants_==2) {
703      chosen_ = nodes_[chosen_].previous_;
704      assert (chosen_>=0);
705    }
706  }
707  return chosen_;
708}
709// Last one in (or other criterion)
710DBNodeSimple
711DBVectorNode::back() const
712{
713  assert (last_>=0);
714  return nodes_[best()];
715}
716// Get rid of last one
717void 
718DBVectorNode::pop_back()
719{
720  // Temporary until more sophisticated
721  //assert (last_==chosen_);
722  if (nodes_[chosen_].descendants_==2)
723    sizeDeferred_--;
724  int previous = nodes_[chosen_].previous_;
725  int next = nodes_[chosen_].next_;
726  nodes_[chosen_].gutsOfDestructor();
727  if (previous>=0) {
728    nodes_[previous].next_=next;
729  } else {
730    first_ = next;
731  }
732  if (next>=0) {
733    nodes_[next].previous_ = previous;
734  } else {
735    last_ = previous;
736  }
737  nodes_[chosen_].previous_=-1;
738  if (firstSpare_>=0) {
739    nodes_[chosen_].next_ = firstSpare_;
740  } else {
741    nodes_[chosen_].next_ = -1;
742  }
743  firstSpare_ = chosen_;
744  chosen_ = -1;
745  assert (size_>0);
746  size_--;
747}
748
749static double
750compute_val_for_ray(const OsiSolverInterface& model)
751{
752  const std::vector<double*> dual_rays = model.getDualRays(1);
753  if (dual_rays.size() == 0) {
754    printf("WARNING: LP is infeas, but no dual ray is returned!\n");
755    return 0;
756  }
757  const double* dual_ray = dual_rays[0];
758  const double direction = model.getObjSense();
759  const double* rub = model.getRowUpper();
760  const double* rlb = model.getRowLower();
761  const double* cub = model.getColUpper();
762  const double* clb = model.getColLower();
763  const double* dj  = model.getReducedCost();
764  const double* obj = model.getObjCoefficients();
765  const int numRows = model.getNumRows();
766  const int numCols = model.getNumCols();
767  double yb_plus_rl_minus_su = 0;
768  for (int i = 0; i < numRows; ++i) {
769    const double ray_i = dual_ray[i];
770    if (ray_i > 1e-6) {
771      yb_plus_rl_minus_su += ray_i*rlb[i];
772    } else if (ray_i < 1e-6) {
773      yb_plus_rl_minus_su += ray_i*rub[i];
774    }
775  }
776  for (int j = 0; j < numCols; ++j) {
777    const double yA_j = dj[j] - obj[j];
778    if (direction * yA_j > 1e-6) {
779      yb_plus_rl_minus_su -= yA_j*cub[j];
780    } else if (direction * yA_j < -1e-6) {
781      yb_plus_rl_minus_su -= yA_j*clb[j];
782    }
783  }
784  for (int i = dual_rays.size()-1; i >= 0; --i) {
785    delete[] dual_rays[i];
786  }
787  return yb_plus_rl_minus_su;
788}
789
790bool
791DBNodeSimple::canSwitchParentWithGrandparent(const int* which,
792                                             OsiSolverInterface & model,
793                                             const int * original_lower,
794                                             const int * original_upper,
795                                             DBVectorNode & branchingTree)
796{
797  /*
798    The current node ('this') is really the parent (P) and the solution in the
799    model represents the child. The grandparent (GP) is this.parent_. Let's have
800    variable names respecting the truth.
801  */
802#if !defined(FUNNY_BRANCHING)
803  return false;
804#endif
805
806  const int parent_id = node_id_;
807  const DBNodeSimple& parent = branchingTree.nodes_[parent_id];
808  const int grandparent_id = parent.parent_;
809
810  if (grandparent_id == -1) {
811    // can't flip root higher up...
812    return false;
813  }
814  const DBNodeSimple& grandparent = branchingTree.nodes_[grandparent_id];
815 
816  if (model.isProvenDualInfeasible()) {
817    // THINK: this is not going to happen, but if it does, what should we do???
818    return false;
819  }
820
821  if (parent.strong_branching_fixed_vars_ || parent.reduced_cost_fixed_vars_ ||
822      grandparent.strong_branching_fixed_vars_) {
823    return false;
824  }
825 
826  const double direction = model.getObjSense();
827
828  const int GP_brvar = grandparent.variable_;
829  const int GP_brvar_fullid = which[GP_brvar];
830  const bool parent_is_down_child = parent_id == grandparent.child_down_;
831
832  if (model.isProvenPrimalInfeasible()) {
833    const int greatgrandparent_id = grandparent.parent_;
834    const int x = GP_brvar_fullid; // for easier reference... we'll use s_x
835
836    /*
837      Now we are going to check that if we relax the GP's branching decision
838      then the child's LP relaxation is still infeasible. The test is done by
839      checking whether the column (or its negative) of the GP's branching
840      variable cuts off the dual ray proving the infeasibility.
841    */
842   
843    const double* dj = model.getReducedCost();
844    const double* obj = model.getObjCoefficients();
845    const double yA_x = dj[x] - obj[x];
846    bool canSwitch = true;
847
848    if (direction > 0) { // minimization problem
849      if (parent_is_down_child) {
850        const double s_x = CoinMax(yA_x, -1e-8);
851        if (s_x > 1e-6) { // if 0 then canSwitch can stay true
852          const double yb_plus_rl_minus_su = compute_val_for_ray(model);
853          if (yb_plus_rl_minus_su < 1e-8) {
854            printf("WARNING: yb_plus_rl_minus_su is not positive!\n");
855            canSwitch = false;
856          } else {
857            const double max_u_x = yb_plus_rl_minus_su / s_x + model.getColUpper()[x];
858            const double u_x_without_GP = greatgrandparent_id >= 0 ?
859              branchingTree.nodes_[greatgrandparent_id].upper_[GP_brvar] :
860              original_upper[GP_brvar];
861            canSwitch = max_u_x > u_x_without_GP - 1e-8;
862          }
863        }
864      } else {
865        const double r_x = CoinMax(yA_x, -1e-8);
866        if (r_x > 1e-6) { // if 0 then canSwitch can stay true
867          const double yb_plus_rl_minus_su = compute_val_for_ray(model);
868          if (yb_plus_rl_minus_su < 1e-8) {
869            printf("WARNING: yb_plus_rl_minus_su is not positive!\n");
870            canSwitch = false;
871          } else {
872            const double min_l_x = - yb_plus_rl_minus_su / r_x + model.getColLower()[x];
873            const double l_x_without_GP = greatgrandparent_id >= 0 ?
874              branchingTree.nodes_[greatgrandparent_id].lower_[GP_brvar] :
875              original_lower[GP_brvar];
876            canSwitch = min_l_x < l_x_without_GP + 1e-8;
877          }
878        }
879      }
880    } else { // maximization problem
881      if (parent_is_down_child) {
882        const double s_x = CoinMin(yA_x, 1e-8);
883        if (s_x < -1e-6) { // if 0 then canSwitch can stay true
884          const double yb_plus_rl_minus_su = compute_val_for_ray(model);
885          if (yb_plus_rl_minus_su > -1e-8) {
886            printf("WARNING: yb_plus_rl_minus_su is not negative!\n");
887            canSwitch = false;
888          } else {
889            const double max_u_x = yb_plus_rl_minus_su / s_x + model.getColUpper()[x];
890            const double u_x_without_GP = greatgrandparent_id >= 0 ?
891              branchingTree.nodes_[greatgrandparent_id].upper_[GP_brvar] :
892              original_upper[GP_brvar];
893            canSwitch = max_u_x > u_x_without_GP - 1e-8;
894          }
895        }
896      } else {
897        const double r_x = CoinMin(yA_x, 1e-8);
898        if (r_x > -1e-6) { // if 0 then canSwitch can stay true
899          const double yb_plus_rl_minus_su = compute_val_for_ray(model);
900          if (yb_plus_rl_minus_su > -1e-8) {
901            printf("WARNING: yb_plus_rl_minus_su is not negative!\n");
902            canSwitch = false;
903          } else {
904            const double min_l_x = - yb_plus_rl_minus_su / r_x + model.getColLower()[x];
905            const double l_x_without_GP = greatgrandparent_id >= 0 ?
906              branchingTree.nodes_[greatgrandparent_id].lower_[GP_brvar] :
907              original_lower[GP_brvar];
908            canSwitch = min_l_x < l_x_without_GP + 1e-8;
909          }
910        }
911      }
912    }
913
914    return canSwitch;
915  }
916  // Both primal and dual feasible, and in this case we don't care how we have
917  // stopped (iteration limit, obj val limit, time limit, optimal solution,
918  // etc.), we can just look at the reduced costs to figure out if the
919  // grandparent is irrelevant. Remember, we are using dual simplex!
920  double djValue = model.getReducedCost()[GP_brvar_fullid]*direction;
921  if (djValue > 1.0e-6) {
922    // wants to go down
923    if (parent_is_down_child) {
924      return true;
925    }
926    if (model.getColLower()[GP_brvar_fullid] > std::ceil(grandparent.value_)) {
927      return true;
928    }
929    return false;
930  } else if (djValue < -1.0e-6) {
931    // wants to go up
932    if (! parent_is_down_child) {
933      return true;
934    }
935    if (model.getColUpper()[GP_brvar_fullid] < std::floor(grandparent.value_)) {
936      return true;
937    }
938    return false;
939  }
940  return false;
941}
942
943void
944DBVectorNode::moveNodeUp(const int* which,
945                         OsiSolverInterface& model, DBNodeSimple & node)
946{
947  /*
948    The current node ('this') is really the parent (P) and the solution in the
949    model represents the child. The grandparent (GP) is this.parent_. Let's have
950    variable names respecting the truth.
951  */
952  const int parent_id = node.node_id_;
953  DBNodeSimple& parent = nodes_[parent_id];
954  const int grandparent_id = parent.parent_;
955  assert(grandparent_id != -1);
956  DBNodeSimple& grandparent = nodes_[grandparent_id];
957  const int greatgrandparent_id = grandparent.parent_;
958 
959  const bool parent_is_down_child = parent_id == grandparent.child_down_;
960
961  // First hang the nodes where they belong.
962  parent.parent_ = greatgrandparent_id;
963  grandparent.parent_ = parent_id;
964  const bool down_child_stays_with_parent = parent.way_ & DBNodeSimple::DOWN_CURRENT;
965  int& child_to_move =
966    down_child_stays_with_parent ? parent.child_up_ : parent.child_down_;
967  if (parent_is_down_child) {
968    grandparent.child_down_ = child_to_move;
969  } else {
970    grandparent.child_up_ = child_to_move;
971  }
972  if (child_to_move >= 0) {
973    nodes_[child_to_move].parent_ = grandparent_id;
974  }
975  child_to_move = grandparent_id;
976  if (greatgrandparent_id >= 0) {
977    DBNodeSimple& greatgrandparent = nodes_[greatgrandparent_id];
978    if (grandparent_id == greatgrandparent.child_down_) {
979      greatgrandparent.child_down_ = parent_id;
980    } else {
981      greatgrandparent.child_up_ = parent_id;
982    }
983  }
984
985  // Now modify bounds
986
987  // First, get rid of GP's bound change of its branching variable in the
988  // bound list of P. Note: If greatgrandparent_id is < 0 then GP is the root
989  // so its bounds are the original bounds.
990  const int GP_brvar = grandparent.variable_;
991  if (parent_is_down_child) {
992    const int old_upper = greatgrandparent_id >= 0 ?
993      nodes_[greatgrandparent_id].upper_[GP_brvar] :
994      grandparent.upper_[GP_brvar];
995    parent.upper_[GP_brvar] = old_upper;
996    model.setColUpper(which[GP_brvar], old_upper);
997  } else {
998    const int old_lower = greatgrandparent_id >= 0 ?
999      nodes_[greatgrandparent_id].lower_[GP_brvar] :
1000      grandparent.lower_[GP_brvar];
1001    parent.lower_[GP_brvar] = old_lower;
1002    model.setColLower(which[GP_brvar], old_lower);
1003  }
1004  // THINK: this might not be necessary
1005  model.resolve();
1006
1007  // Now add the branching var bound change of P to GP and all of its
1008  // descendant
1009  const int P_brvar = parent.variable_;
1010  const double P_value = parent.value_;
1011  int new_bd;
1012  if (down_child_stays_with_parent) {
1013    // Former up child of P is now the down child of GP, so we need to change
1014    // bounds of GP, its up child and descendants of that one.
1015    new_bd = (int)std::ceil(P_value);
1016    grandparent.lower_[P_brvar] = new_bd;
1017    changeDescendantBounds(grandparent.child_up_,
1018                           true /*lower bd*/, P_brvar, new_bd);
1019  } else {
1020    // Former down child of P is now the up child of GP, so we need to change
1021    // bounds of GP, its down child and descendants of that one.
1022    new_bd = (int)floor(P_value);
1023    grandparent.upper_[P_brvar] = new_bd;
1024    changeDescendantBounds(grandparent.child_down_,
1025                           false /*lower bd*/, P_brvar, new_bd);
1026  }
1027}
1028
1029void
1030DBVectorNode::changeDescendantBounds(const int node_id, const bool lower_bd,
1031                                     const int brvar, const int new_bd)
1032{
1033  if (node_id == -1) {
1034    return;
1035  }
1036  changeDescendantBounds(nodes_[node_id].child_down_, lower_bd, brvar, new_bd);
1037  changeDescendantBounds(nodes_[node_id].child_up_, lower_bd, brvar, new_bd);
1038  if (lower_bd) {
1039    nodes_[node_id].lower_[brvar] = new_bd;
1040  } else {
1041    nodes_[node_id].upper_[brvar] = new_bd;
1042  }
1043}
1044
1045
1046// Invoke solver's built-in enumeration algorithm
1047void 
1048branchAndBound(OsiSolverInterface & model) {
1049  double time1 = CoinCpuTime();
1050  // solve LP
1051  model.initialSolve();
1052
1053  if (model.isProvenOptimal()&&!model.isDualObjectiveLimitReached()) {
1054    // Continuous is feasible - find integers
1055    int numberIntegers=0;
1056    int numberColumns = model.getNumCols();
1057    int iColumn;
1058    int i;
1059    for (iColumn=0;iColumn<numberColumns;iColumn++) {
1060      if( model.isInteger(iColumn))
1061        numberIntegers++;
1062    }
1063    if (!numberIntegers) {
1064      std::cout<<"No integer variables"
1065               <<std::endl;
1066      return;
1067    }
1068    int * which = new int[numberIntegers]; // which variables are integer
1069    // original bounds
1070    int * originalLower = new int[numberIntegers];
1071    int * originalUpper = new int[numberIntegers];
1072    int * relaxedLower = new int[numberIntegers];
1073    int * relaxedUpper = new int[numberIntegers];
1074    {
1075      const double * lower = model.getColLower();
1076      const double * upper = model.getColUpper();
1077      numberIntegers=0;
1078      for (iColumn=0;iColumn<numberColumns;iColumn++) {
1079        if( model.isInteger(iColumn)) {
1080          originalLower[numberIntegers]=(int) lower[iColumn];
1081          originalUpper[numberIntegers]=(int) upper[iColumn];
1082          which[numberIntegers++]=iColumn;
1083        }
1084      }
1085    }
1086    double direction = model.getObjSense();
1087    // empty tree
1088    DBVectorNode branchingTree;
1089   
1090    // Add continuous to it;
1091    DBNodeSimple rootNode(model,numberIntegers,which,model.getWarmStart());
1092    // something extra may have been fixed by strong branching
1093    // if so go round again
1094    while (rootNode.variable_==numberIntegers) {
1095      model.resolve();
1096      rootNode = DBNodeSimple(model,numberIntegers,which,model.getWarmStart());
1097    }
1098    if (rootNode.objectiveValue_<1.0e100) {
1099      // push on stack
1100      branchingTree.push_back(rootNode);
1101    }
1102   
1103    // For printing totals
1104    int numberIterations=0;
1105    int numberNodes =0;
1106    DBNodeSimple bestNode;
1107    ////// Start main while of branch and bound
1108    // while until nothing on stack
1109    while (branchingTree.size()) {
1110      // last node
1111      DBNodeSimple node = branchingTree.back();
1112      int kNode = branchingTree.chosen_;
1113      branchingTree.pop_back();
1114      assert (node.descendants_<2);
1115      numberNodes++;
1116      if (node.variable_>=0) {
1117        // branch - do bounds
1118        for (i=0;i<numberIntegers;i++) {
1119          iColumn=which[i];
1120          model.setColBounds( iColumn,node.lower_[i],node.upper_[i]);
1121        }
1122        // move basis
1123        model.setWarmStart(node.basis_);
1124        // do branching variable
1125        node.incrementDescendants();
1126        bool down_branch = true;
1127        switch (node.way_) {
1128        case DBNodeSimple::WAY_UNSET:
1129        case DBNodeSimple::DOWN_UP__BOTH_DONE:
1130        case DBNodeSimple::UP_DOWN__BOTH_DONE:
1131          abort();
1132        case DBNodeSimple::DOWN_UP__NOTHING_DONE:
1133          node.way_ = DBNodeSimple::DOWN_UP__DOWN_DONE;
1134          break;
1135        case DBNodeSimple::DOWN_UP__DOWN_DONE:
1136          node.way_ = DBNodeSimple::DOWN_UP__BOTH_DONE;
1137          down_branch = false;
1138          break;
1139        case DBNodeSimple::UP_DOWN__NOTHING_DONE:
1140          node.way_ = DBNodeSimple::UP_DOWN__UP_DONE;
1141          down_branch = false;
1142          break;
1143        case DBNodeSimple::UP_DOWN__UP_DONE:
1144          node.way_ = DBNodeSimple::UP_DOWN__BOTH_DONE;
1145          break;
1146        }
1147        if (down_branch) {
1148          model.setColUpper(which[node.variable_],floor(node.value_));
1149        } else {
1150          model.setColLower(which[node.variable_],ceil(node.value_));
1151        }
1152        // put back on tree anyway regardless whether any processing is left
1153        // to be done. We want the whole tree all the time.
1154        branchingTree.push_back(node);
1155       
1156        // solve
1157        model.resolve();
1158        CoinWarmStart * ws = model.getWarmStart();
1159        const CoinWarmStartBasis* wsb =
1160          dynamic_cast<const CoinWarmStartBasis*>(ws);
1161        assert (wsb!=NULL); // make sure not volume
1162        numberIterations += model.getIterationCount();
1163        // fix on reduced costs
1164        int nFixed0=0,nFixed1=0;
1165        double cutoff;
1166        model.getDblParam(OsiDualObjectiveLimit,cutoff);
1167        double gap=(cutoff-model.getObjValue())*direction+1.0e-4;
1168        //        double
1169        //        gap=(cutoff-modelPtr_->objectiveValue())*direction+1.0e-4;
1170        bool did_reduced_cost_fixing_for_child = false;
1171        if (gap<1.0e10&&model.isProvenOptimal()&&!model.isDualObjectiveLimitReached()) {
1172          const double * dj = model.getReducedCost();
1173          const double * lower = model.getColLower();
1174          const double * upper = model.getColUpper();
1175          for (i=0;i<numberIntegers;i++) {
1176            iColumn=which[i];
1177            if (upper[iColumn]>lower[iColumn]) {
1178              double djValue = dj[iColumn]*direction;
1179              if (wsb->getStructStatus(iColumn)==CoinWarmStartBasis::atLowerBound&&
1180                  djValue>gap) {
1181                nFixed0++;
1182                model.setColUpper(iColumn,lower[iColumn]);
1183              } else if (wsb->getStructStatus(iColumn)==CoinWarmStartBasis::atUpperBound&&
1184                         -djValue>gap) {
1185                nFixed1++;
1186                model.setColLower(iColumn,upper[iColumn]);
1187              }
1188            }
1189          }
1190          if (nFixed0+nFixed1) {
1191            //      printf("%d fixed to lower, %d fixed to upper\n",
1192            //             nFixed0,nFixed1);
1193            did_reduced_cost_fixing_for_child = true;
1194          }
1195        }
1196        if (model.isAbandoned()) {
1197          // THINK: What the heck should we do???
1198          abort();
1199        }
1200        if (model.isIterationLimitReached()) {
1201          // maximum iterations - exit
1202          std::cout<<"Exiting on maximum iterations\n";
1203          break;
1204        }
1205
1206        while (node.canSwitchParentWithGrandparent(which, model,
1207                                                   originalLower,
1208                                                   originalUpper,
1209                                                   branchingTree)) {
1210          branchingTree.moveNodeUp(which, model, node);
1211        }
1212        if ((numberNodes%1000)==0) 
1213          printf("%d nodes, tree size %d\n",
1214                 numberNodes,branchingTree.size());
1215        if (CoinCpuTime()-time1>3600.0) {
1216          printf("stopping after 3600 seconds\n");
1217          exit(77);
1218        }
1219        DBNodeSimple newNode(model,numberIntegers,which,ws);
1220        // something extra may have been fixed by strong branching
1221        // if so go round again
1222        while (newNode.variable_==numberIntegers) {
1223          model.resolve();
1224          newNode = DBNodeSimple(model,numberIntegers,which,model.getWarmStart());
1225          newNode.strong_branching_fixed_vars_ = true;
1226        }
1227        newNode.reduced_cost_fixed_vars_ = did_reduced_cost_fixing_for_child;
1228        if (newNode.objectiveValue_<1.0e100) {
1229          newNode.parent_ = kNode;
1230          // push on stack
1231          branchingTree.push_back(newNode);
1232          if(branchingTree.nodes_[kNode].child_down_ < 0)
1233            branchingTree.nodes_[kNode].child_down_ = branchingTree.last_;
1234          else
1235            branchingTree.nodes_[kNode].child_up_ = branchingTree.last_;
1236          if (newNode.variable_>=0) {
1237            assert (fabs(newNode.value_-floor(newNode.value_+0.5))>1.0e-6);
1238          }
1239#if 0
1240          else {
1241            // integer solution - save
1242            bestNode = node;
1243            // set cutoff (hard coded tolerance)
1244            model.setDblParam(OsiDualObjectiveLimit,(bestNode.objectiveValue_-1.0e-5)*direction);
1245            std::cout<<"Integer solution of "
1246                     <<bestNode.objectiveValue_
1247                     <<" found after "<<numberIterations
1248                     <<" iterations and "<<numberNodes<<" nodes"
1249                     <<std::endl;
1250          }
1251#endif
1252        }
1253      } else {
1254        // Integer solution - save
1255        bestNode = node;
1256        // set cutoff (hard coded tolerance)
1257        model.setDblParam(OsiDualObjectiveLimit,(bestNode.objectiveValue_-1.0e-5)*direction);
1258        std::cout<<"Integer solution of "
1259                 <<bestNode.objectiveValue_
1260                 <<" found after "<<numberIterations
1261                 <<" iterations and "<<numberNodes<<" nodes"
1262                 <<std::endl;
1263      }
1264    }
1265    ////// End main while of branch and bound
1266    std::cout<<"Search took "
1267             <<numberIterations
1268             <<" iterations and "<<numberNodes<<" nodes"
1269             <<std::endl;
1270    if (bestNode.numberIntegers_) {
1271      // we have a solution restore
1272      // do bounds
1273      for (i=0;i<numberIntegers;i++) {
1274        iColumn=which[i];
1275        model.setColBounds( iColumn,bestNode.lower_[i],bestNode.upper_[i]);
1276      }
1277      // move basis
1278      model.setWarmStart(bestNode.basis_);
1279      // set cutoff so will be good (hard coded tolerance)
1280      model.setDblParam(OsiDualObjectiveLimit,(bestNode.objectiveValue_+1.0e-5)*direction);
1281      model.resolve();
1282    } else {
1283      OsiClpSolverInterface* clp =
1284        dynamic_cast<OsiClpSolverInterface*>(&model);
1285      if (clp) {
1286        ClpSimplex* modelPtr_ = clp->getModelPtr();
1287        modelPtr_->setProblemStatus(1);
1288      }
1289    }
1290    delete [] which;
1291    delete [] originalLower;
1292    delete [] originalUpper;
1293    delete [] relaxedLower;
1294    delete [] relaxedUpper;
1295  } else {
1296    std::cout<<"The LP relaxation is infeasible"
1297             <<std::endl;
1298    OsiClpSolverInterface* clp =
1299      dynamic_cast<OsiClpSolverInterface*>(&model);
1300    if (clp) {
1301      ClpSimplex* modelPtr_ = clp->getModelPtr();
1302      modelPtr_->setProblemStatus(1);
1303    }
1304    //throw CoinError("The LP relaxation is infeasible or too expensive",
1305    //"branchAndBound", "OsiClpSolverInterface");
1306  }
1307}
1308
1309
1310
1311int main(int argc, char* argv[])
1312{
1313  OsiClpSolverInterface model;
1314  model.readMps(argv[1]);
1315  branchAndBound(model);
1316  return 0;
1317}
Note: See TracBrowser for help on using the repository browser.