Changeset 965


Ignore:
Timestamp:
Jun 5, 2008 11:13:18 AM (11 years ago)
Author:
jpgoncal
Message:

Copied code from OsiClpSolverInterface?.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/dynamicbranching/dynamicbranching.cpp

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