source: trunk/CbcBranchCut.cpp @ 129

Last change on this file since 129 was 129, checked in by forrest, 16 years ago

make createBranch non const

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 16.9 KB
Line 
1// Copyright (C) 2004, International Business Machines
2// Corporation and others.  All Rights Reserved.
3#if defined(_MSC_VER)
4// Turn off compiler warning about long names
5#  pragma warning(disable:4786)
6#endif
7#include <cassert>
8#include <cmath>
9#include <cfloat>
10//#define CBC_DEBUG
11
12#include "OsiSolverInterface.hpp"
13#include "CbcModel.hpp"
14#include "CbcMessage.hpp"
15#include "CbcBranchCut.hpp"
16#include "CoinSort.hpp"
17#include "CoinError.hpp"
18
19
20/** Default Constructor
21
22*/
23CbcBranchCut::CbcBranchCut ()
24  : CbcObject()
25{
26}
27
28/* Constructor so model can be passed up
29*/ 
30CbcBranchCut::CbcBranchCut (CbcModel * model)
31  : CbcObject(model)
32{
33}
34// Copy constructor
35CbcBranchCut::CbcBranchCut ( const CbcBranchCut & rhs)
36  :CbcObject(rhs)
37
38{
39}
40
41// Clone
42CbcObject *
43CbcBranchCut::clone() const
44{
45  return new CbcBranchCut(*this);
46}
47
48// Assignment operator
49CbcBranchCut & 
50CbcBranchCut::operator=( const CbcBranchCut& rhs)
51{
52  return *this;
53}
54
55// Destructor
56CbcBranchCut::~CbcBranchCut ()
57{
58}
59
60// Infeasibility - large is 0.5
61double 
62CbcBranchCut::infeasibility(int & preferredWay) const
63{
64  throw CoinError("Use of base class","infeasibility","CbcBranchCut");
65  preferredWay=-1;
66  return 0.0;
67}
68
69// This looks at solution and sets bounds to contain solution
70/** More precisely: it first forces the variable within the existing
71    bounds, and then tightens the bounds to fix the variable at the
72    nearest integer value.
73*/
74void 
75CbcBranchCut::feasibleRegion()
76{
77}
78/* Return true if branch created by object should fix variables
79 */
80bool 
81CbcBranchCut::boundBranch() const 
82{return false;}
83
84// Creates a branching object
85CbcBranchingObject * 
86CbcBranchCut::createBranch(int way) 
87{
88  throw CoinError("Use of base class","createBranch","CbcBranchCut");
89  return new CbcCutBranchingObject();
90}
91
92
93/* Given valid solution (i.e. satisfied) and reduced costs etc
94   returns a branching object which would give a new feasible
95   point in direction reduced cost says would be cheaper.
96   If no feasible point returns null
97*/
98CbcBranchingObject * 
99CbcBranchCut::preferredNewFeasible() const
100{
101  throw CoinError("Use of base class","preferredNewFeasible","CbcBranchCut");
102  return new CbcCutBranchingObject();
103}
104 
105/* Given valid solution (i.e. satisfied) and reduced costs etc
106   returns a branching object which would give a new feasible
107   point in direction opposite to one reduced cost says would be cheaper.
108   If no feasible point returns null
109*/
110CbcBranchingObject * 
111CbcBranchCut::notPreferredNewFeasible() const 
112{
113  throw CoinError("Use of base class","notPreferredNewFeasible","CbcBranchCut");
114  return new CbcCutBranchingObject();
115}
116 
117/*
118  Bounds may be tightened, so it may be good to be able to refresh the local
119  copy of the original bounds.
120 */
121void 
122CbcBranchCut::resetBounds()
123{
124}
125
126
127// Default Constructor
128CbcCutBranchingObject::CbcCutBranchingObject()
129  :CbcBranchingObject()
130{
131  down_=OsiRowCut();
132  up_=OsiRowCut();
133}
134
135// Useful constructor
136CbcCutBranchingObject::CbcCutBranchingObject (CbcModel * model, 
137                                              OsiRowCut & down,
138                                              OsiRowCut &up)
139  :CbcBranchingObject(model,0,-1,0.0)
140{
141  down_ = down;
142  up_ = up;
143}
144
145// Copy constructor
146CbcCutBranchingObject::CbcCutBranchingObject ( const CbcCutBranchingObject & rhs) :CbcBranchingObject(rhs)
147{
148  down_ = rhs.down_;
149  up_ = rhs.up_;
150}
151
152// Assignment operator
153CbcCutBranchingObject & 
154CbcCutBranchingObject::operator=( const CbcCutBranchingObject& rhs)
155{
156  if (this != &rhs) {
157    CbcBranchingObject::operator=(rhs);
158    down_ = rhs.down_;
159    up_ = rhs.up_;
160  }
161  return *this;
162}
163CbcBranchingObject * 
164CbcCutBranchingObject::clone() const
165{ 
166  return (new CbcCutBranchingObject(*this));
167}
168
169
170// Destructor
171CbcCutBranchingObject::~CbcCutBranchingObject ()
172{
173}
174
175/*
176  Perform a branch by adjusting bounds and/or adding a cut. Note
177  that each arm of the branch advances the object to the next arm by
178  advancing the value of way_.
179
180  Returns change in guessed objective on next branch
181*/
182double
183CbcCutBranchingObject::branch(bool normalBranch)
184{
185  if (model_->messageHandler()->logLevel()>2&&normalBranch)
186    print(normalBranch);
187  numberBranchesLeft_--;
188  OsiRowCut * cut;
189  if (way_<0) {
190    cut = &down_;
191    way_=1;
192  } else {
193    cut = &up_;
194    way_=-1;      // Swap direction
195  }
196  // See if cut just fixes variables
197  double lb = cut->lb();
198  double ub = cut->ub();
199  int n=cut->row().getNumElements();
200  const int * column = cut->row().getIndices();
201  const double * element = cut->row().getElements();
202  OsiSolverInterface * solver = model_->solver();
203  const double * upper = solver->getColUpper();
204  const double * lower = solver->getColLower();
205  double low = 0.0;
206  double high=0.0;
207  for (int i=0;i<n;i++) {
208    int iColumn = column[i];
209    double value = element[i];
210    if (value>0.0) {
211      high += upper[iColumn]*value;
212      low += lower[iColumn]*value;
213    } else {
214      high += lower[iColumn]*value;
215      low += upper[iColumn]*value;
216    }
217  }
218  // assume cut was cunningly constructed so we need not worry too much about tolerances
219  if (low+1.0e-8>=ub) {
220    // fix
221    for (int i=0;i<n;i++) {
222      int iColumn = column[i];
223      double value = element[i];
224      if (value>0.0) {
225        solver->setColUpper(iColumn,lower[iColumn]);
226      } else {
227        solver->setColLower(iColumn,upper[iColumn]);
228      }
229    }
230  } else if (high-1.0e-8<=lb) {
231    // fix
232    for (int i=0;i<n;i++) {
233      int iColumn = column[i];
234      double value = element[i];
235      if (value>0.0) {
236        solver->setColLower(iColumn,upper[iColumn]);
237      } else {
238        solver->setColUpper(iColumn,lower[iColumn]);
239      }
240    }
241  } else {
242    // leave as cut
243    model_->setNextRowCut(*cut);
244  }
245  return 0.0;
246}
247// Print what would happen 
248void
249CbcCutBranchingObject::print(bool normalBranch)
250{
251  OsiRowCut * cut;
252  if (way_<0) {
253    cut = &down_;
254    printf("CbcCut would branch down");
255    way_=1;
256  } else {
257    cut = &up_;
258    printf("CbcCut would branch up");
259    way_=-1;      // Swap direction
260  }
261  double lb = cut->lb();
262  double ub = cut->ub();
263  int n=cut->row().getNumElements();
264  const int * column = cut->row().getIndices();
265  const double * element = cut->row().getElements();
266  if (n>5) {
267    printf(" - %d elements, lo=%g, up=%g\n",n,lb,ub);
268  } else {
269    printf(" - %g <=",lb);
270    for (int i=0;i<n;i++) {
271      int iColumn = column[i];
272      double value = element[i];
273      printf(" (%d,%g)",iColumn,value);
274    }
275    printf(" <= %g\n",ub);
276  }
277}
278
279// Return true if branch should fix variables
280bool 
281CbcCutBranchingObject::boundBranch() const
282{
283  return false;
284}
285
286/** Default Constructor
287
288  Equivalent to an unspecified binary variable.
289*/
290CbcBranchToFixLots::CbcBranchToFixLots ()
291  : CbcBranchCut(),
292    djTolerance_(COIN_DBL_MAX),
293    fractionFixed_(1.0),
294    mark_(NULL),
295    depth_(-1),
296    numberClean_(0),
297    alwaysCreate_(false)
298{
299}
300
301/* Useful constructor - passed reduced cost tolerance and fraction we would like fixed.
302   Also depth level to do at.
303   Also passed number of 1 rows which when clean triggers fix
304   Always does if all 1 rows cleaned up and number>0 or if fraction columns reached
305   Also whether to create branch if can't reach fraction.
306*/ 
307CbcBranchToFixLots::CbcBranchToFixLots (CbcModel * model, double djTolerance,
308                                        double fractionFixed, int depth,
309                                        int numberClean,
310                                        const char * mark, bool alwaysCreate)
311  : CbcBranchCut(model)
312{
313  djTolerance_ = djTolerance;
314  fractionFixed_ = fractionFixed;
315  if (mark) {
316    int numberColumns = model->getNumCols();
317    mark_ = new char[numberColumns];
318    memcpy(mark_,mark,numberColumns);
319  }
320  depth_ = depth;
321  assert (model);
322  OsiSolverInterface * solver = model_->solver();
323  matrixByRow_ = *solver->getMatrixByRow();
324  numberClean_ = numberClean;
325  alwaysCreate_ = alwaysCreate;
326}
327// Copy constructor
328CbcBranchToFixLots::CbcBranchToFixLots ( const CbcBranchToFixLots & rhs)
329  :CbcBranchCut(rhs)
330{
331  djTolerance_ = rhs.djTolerance_;
332  fractionFixed_ = rhs.fractionFixed_;
333  int numberColumns = model_->getNumCols();
334  mark_ = CoinCopyOfArray(rhs.mark_,numberColumns);
335  matrixByRow_=rhs.matrixByRow_;
336  depth_ = rhs.depth_;
337  numberClean_ = rhs.numberClean_;
338  alwaysCreate_ = rhs.alwaysCreate_;
339}
340
341// Clone
342CbcObject *
343CbcBranchToFixLots::clone() const
344{
345  return new CbcBranchToFixLots(*this);
346}
347
348// Assignment operator
349CbcBranchToFixLots & 
350CbcBranchToFixLots::operator=( const CbcBranchToFixLots& rhs)
351{
352  if (this!=&rhs) {
353    CbcBranchCut::operator=(rhs);
354    djTolerance_ = rhs.djTolerance_;
355    fractionFixed_ = rhs.fractionFixed_;
356    int numberColumns = model_->getNumCols();
357    mark_ = CoinCopyOfArray(rhs.mark_,numberColumns);
358    matrixByRow_=rhs.matrixByRow_;
359    depth_ = rhs.depth_;
360    numberClean_ = rhs.numberClean_;
361    alwaysCreate_ = rhs.alwaysCreate_;
362  }
363  return *this;
364}
365
366// Destructor
367CbcBranchToFixLots::~CbcBranchToFixLots ()
368{
369  delete [] mark_;
370}
371// Creates a branching object
372CbcBranchingObject * 
373CbcBranchToFixLots::createBranch(int way) const
374{
375  // by default way must be -1
376  assert (way==-1);
377  OsiSolverInterface * solver = model_->solver();
378  const double * solution = model_->currentSolution();
379  const double * lower = solver->getColLower();
380  const double * upper = solver->getColUpper();
381  const double * dj = solver->getReducedCost();
382  int i;
383  int numberIntegers = model_->numberIntegers();
384  const int * integerVariable = model_->integerVariable();
385  double integerTolerance = 
386    model_->getDblParam(CbcModel::CbcIntegerTolerance);
387  // make smaller ?
388  double tolerance = CoinMin(1.0e-8,integerTolerance);
389  // How many fixed are we aiming at
390  int wantedFixed = (int) ((double)numberIntegers*fractionFixed_);
391  int nSort=0;
392  int numberFixed=0;
393  int numberColumns = solver->getNumCols();
394  int * sort = new int[numberColumns];
395  double * dsort = new double[numberColumns];
396  int type = shallWe();
397  assert (type);
398  // Take clean first
399  if (type==1) {
400    for (i=0;i<numberIntegers;i++) {
401      int iColumn = integerVariable[i];
402      if (upper[iColumn]>lower[iColumn]) {
403        if (!mark_||!mark_[iColumn]) {
404          if(solution[iColumn]<lower[iColumn]+tolerance) {
405            if (dj[iColumn]>djTolerance_) {
406              dsort[nSort]=-dj[iColumn];
407              sort[nSort++]=iColumn;
408            }
409          } else if (solution[iColumn]>upper[iColumn]-tolerance) {
410            if (dj[iColumn]<-djTolerance_) {
411              dsort[nSort]=dj[iColumn];
412              sort[nSort++]=iColumn;
413            }
414          }
415        }
416      } else {
417        numberFixed++;
418      }
419    }
420    // sort
421    CoinSort_2(dsort,dsort+nSort,sort);
422    nSort= CoinMin(nSort,wantedFixed-numberFixed);
423  } else {
424    int i;
425    //const double * rowLower = solver->getRowLower();
426    const double * rowUpper = solver->getRowUpper();
427    // Row copy
428    const double * elementByRow = matrixByRow_.getElements();
429    const int * column = matrixByRow_.getIndices();
430    const CoinBigIndex * rowStart = matrixByRow_.getVectorStarts();
431    const int * rowLength = matrixByRow_.getVectorLengths();
432    const double * columnLower = solver->getColLower();
433    const double * columnUpper = solver->getColUpper();
434    const double * solution = solver->getColSolution();
435    int numberColumns = solver->getNumCols();
436    int numberRows = solver->getNumRows();
437    for (i=0;i<numberColumns;i++) {
438      sort[i]=i;
439      if (columnLower[i]!=columnUpper[i]){
440        dsort[i]=1.0e100;
441      } else {
442        dsort[i]=1.0e50;
443        numberFixed++;
444      }
445    }
446    for (i=0;i<numberRows;i++) {
447      double rhsValue = rowUpper[i];
448      bool oneRow=true;
449      // check elements
450      int numberUnsatisfied=0;
451      for (int j=rowStart[i];j<rowStart[i]+rowLength[i];j++) {
452        int iColumn = column[j];
453        double value = elementByRow[j];
454        double solValue = solution[iColumn];
455        if (columnLower[iColumn]!=columnUpper[iColumn]) {
456          if (solValue<1.0-integerTolerance&&solValue>integerTolerance)
457            numberUnsatisfied++;
458          if (value!=1.0) {
459            oneRow=false;
460            break;
461          }
462        } else {
463          rhsValue -= value*floor(solValue+0.5);
464        }
465      }
466      if (oneRow&&rhsValue<=1.0+tolerance) {
467        if (!numberUnsatisfied) {
468          for (int j=rowStart[i];j<rowStart[i]+rowLength[i];j++) {
469            int iColumn = column[j];
470            if (dsort[iColumn]>1.0e50){
471              dsort[iColumn]=0;
472              nSort++;
473            }
474          }
475        }
476      }
477    }
478    // sort
479    CoinSort_2(dsort,dsort+numberColumns,sort);
480  }
481  OsiRowCut down;
482  down.setLb(-COIN_DBL_MAX);
483  double rhs=0.0;
484  for (i=0;i<nSort;i++) {
485    int iColumn = sort[i];
486    if(solution[iColumn]<lower[iColumn]+tolerance) {
487      rhs += lower[iColumn];
488      dsort[i]=1.0;
489      assert (!lower[iColumn]);
490    } else {
491      assert (solution[iColumn]>upper[iColumn]-tolerance);
492      rhs -= upper[iColumn];
493      dsort[i]=-1.0;
494      //printf("%d at ub of %g\n",iColumn,upper[iColumn]);
495    }
496  }
497  down.setUb(rhs);
498  down.setRow(nSort,sort,dsort);
499  delete [] sort;
500  delete [] dsort;
501  // up is same - just with rhs changed
502  OsiRowCut up = down;
503  up.setLb(rhs +1.0);
504  up.setUb(COIN_DBL_MAX);
505  CbcCutBranchingObject * newObject = 
506    new CbcCutBranchingObject(model_,down,up);
507  if (model_->messageHandler()->logLevel()>1)
508    printf("creating cut in CbcBranchCut\n");
509  return newObject;
510}
511/* Does a lot of the work,
512   Returns 0 if no good, 1 if dj, 2 if clean, 3 if both
513*/
514int 
515CbcBranchToFixLots::shallWe() const
516{
517  int returnCode=0;
518  OsiSolverInterface * solver = model_->solver();
519  int numberRows = matrixByRow_.getNumRows();
520  //if (numberRows!=solver->getNumRows())
521  //return 0;
522  const double * solution = model_->currentSolution();
523  const double * lower = solver->getColLower();
524  const double * upper = solver->getColUpper();
525  const double * dj = solver->getReducedCost();
526  int i;
527  int numberIntegers = model_->numberIntegers();
528  const int * integerVariable = model_->integerVariable();
529  double integerTolerance = 
530    model_->getDblParam(CbcModel::CbcIntegerTolerance);
531  // make smaller ?
532  double tolerance = CoinMin(1.0e-8,integerTolerance);
533  // How many fixed are we aiming at
534  int wantedFixed = (int) ((double)numberIntegers*fractionFixed_);
535  if (djTolerance_<1.0e10) {
536    int nSort=0;
537    int numberFixed=0;
538    for (i=0;i<numberIntegers;i++) {
539      int iColumn = integerVariable[i];
540      if (upper[iColumn]>lower[iColumn]) {
541        if (!mark_||!mark_[iColumn]) {
542          if(solution[iColumn]<lower[iColumn]+tolerance) {
543            if (dj[iColumn]>djTolerance_) {
544              nSort++;
545            }
546          } else if (solution[iColumn]>upper[iColumn]-tolerance) {
547            if (dj[iColumn]<-djTolerance_) {
548              nSort++;
549            }
550          }
551        }
552      } else {
553        numberFixed++;
554      }
555    }
556    if (numberFixed+nSort<wantedFixed&&!alwaysCreate_) {
557      returnCode = 0;
558    } else if (numberFixed<wantedFixed) {
559      returnCode = 1;
560    } else {
561      returnCode = 0;
562    }
563  }
564  if (numberClean_) {
565    // see how many rows clean
566    int i;
567    //const double * rowLower = solver->getRowLower();
568    const double * rowUpper = solver->getRowUpper();
569    // Row copy
570    const double * elementByRow = matrixByRow_.getElements();
571    const int * column = matrixByRow_.getIndices();
572    const CoinBigIndex * rowStart = matrixByRow_.getVectorStarts();
573    const int * rowLength = matrixByRow_.getVectorLengths();
574    const double * columnLower = solver->getColLower();
575    const double * columnUpper = solver->getColUpper();
576    const double * solution = solver->getColSolution();
577    int numberClean=0;
578    bool someToDoYet=false;
579    int numberColumns = solver->getNumCols();
580    char * mark = new char[numberColumns];
581    int numberFixed=0;
582    for (i=0;i<numberColumns;i++) {
583      if (columnLower[i]!=columnUpper[i]){
584        mark[i]=0;
585      } else {
586        mark[i]=1;
587        numberFixed++;
588      }
589    }
590    int numberNewFixed=0;
591    for (i=0;i<numberRows;i++) {
592      double rhsValue = rowUpper[i];
593      bool oneRow=true;
594      // check elements
595      int numberUnsatisfied=0;
596      for (int j=rowStart[i];j<rowStart[i]+rowLength[i];j++) {
597        int iColumn = column[j];
598        double value = elementByRow[j];
599        double solValue = solution[iColumn];
600        if (columnLower[iColumn]!=columnUpper[iColumn]) {
601          if (solValue<1.0-integerTolerance&&solValue>integerTolerance)
602            numberUnsatisfied++;
603          if (value!=1.0) {
604            oneRow=false;
605            break;
606          }
607        } else {
608          rhsValue -= value*floor(solValue+0.5);
609        }
610      }
611      if (oneRow&&rhsValue<=1.0+tolerance) {
612        if (numberUnsatisfied) {
613          someToDoYet=true;
614        } else {
615          numberClean++;
616          for (int j=rowStart[i];j<rowStart[i]+rowLength[i];j++) {
617            int iColumn = column[j];
618            if (columnLower[iColumn]!=columnUpper[iColumn]&&!mark[iColumn]){
619              mark[iColumn]=1;
620              numberNewFixed++;
621            }
622          }
623        }
624      }
625    }
626    delete [] mark;
627    //printf("%d clean, %d old fixed, %d new fixed\n",
628    //   numberClean,numberFixed,numberNewFixed);
629    if (someToDoYet&&numberClean<numberClean_
630        &&numberNewFixed+numberFixed<wantedFixed) {
631    } else if (numberFixed<wantedFixed) {
632      returnCode |= 2;
633    } else {
634    }
635  }
636  return returnCode;
637}
638// Infeasibility - large is 0.5
639double 
640CbcBranchToFixLots::infeasibility(int & preferredWay) const
641{
642  preferredWay=-1;
643  CbcNode * node = model_->currentNode();
644  int depth;
645  if (node)
646    depth=CoinMax(node->depth(),0);
647  else
648    return 0.0;
649  if (depth_<0) {
650    return 0.0;
651  } else if (depth_>0) {
652    if ((depth%depth_)!=0)
653      return 0.0;
654  }
655  if (!shallWe())
656    return 0.0;
657  else
658    return 1.0e20;
659}
Note: See TracBrowser for help on using the repository browser.