source: branches/devel/Cbc/src/CbcCutGenerator.cpp @ 491

Last change on this file since 491 was 482, checked in by forrest, 14 years ago

for nonlinear and sprint

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 14.8 KB
Line 
1// Copyright (C) 2003, 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
11#include "OsiSolverInterface.hpp"
12#include "CbcModel.hpp"
13#include "CbcMessage.hpp"
14#include "CbcCutGenerator.hpp"
15#include "CbcBranchDynamic.hpp"
16#include "CglProbing.hpp"
17#include "CoinTime.hpp"
18
19// Default Constructor
20CbcCutGenerator::CbcCutGenerator ()
21  : model_(NULL),
22    generator_(NULL),
23    whenCutGenerator_(-1),
24    whenCutGeneratorInSub_(-100),
25    switchOffIfLessThan_(0),
26    depthCutGenerator_(-1),
27    depthCutGeneratorInSub_(-1),
28    generatorName_(NULL),
29    normal_(true),
30    atSolution_(false),
31    whenInfeasible_(false),
32    mustCallAgain_(false),
33    switchedOff_(false),
34    timing_(false),
35    timeInCutGenerator_(0.0),
36    numberTimes_(0),
37    numberCuts_(0),
38    numberColumnCuts_(0),
39    numberCutsActive_(0),
40    numberCutsAtRoot_(0),
41    numberActiveCutsAtRoot_(0)
42{
43}
44// Normal constructor
45CbcCutGenerator::CbcCutGenerator(CbcModel * model,CglCutGenerator * generator,
46                                 int howOften, const char * name,
47                                 bool normal, bool atSolution, 
48                                 bool infeasible, int howOftenInSub,
49                                 int whatDepth, int whatDepthInSub,
50                                 int switchOffIfLessThan)
51  : 
52    depthCutGenerator_(whatDepth),
53    depthCutGeneratorInSub_(whatDepthInSub),
54    mustCallAgain_(false),
55    switchedOff_(false),
56    timing_(false),
57    timeInCutGenerator_(0.0),
58    numberTimes_(0),
59    numberCuts_(0),
60    numberColumnCuts_(0),
61    numberCutsActive_(0),
62    numberCutsAtRoot_(0),
63    numberActiveCutsAtRoot_(0)
64{
65  model_ = model;
66  generator_=generator->clone();
67  generator_->refreshSolver(model_->solver());
68  whenCutGenerator_=howOften;
69  whenCutGeneratorInSub_ = howOftenInSub;
70  switchOffIfLessThan_=switchOffIfLessThan;
71  if (name)
72    generatorName_=strdup(name);
73  else
74    generatorName_ = strdup("Unknown");
75  normal_=normal;
76  atSolution_=atSolution;
77  whenInfeasible_=infeasible;
78}
79
80// Copy constructor
81CbcCutGenerator::CbcCutGenerator ( const CbcCutGenerator & rhs)
82{
83  model_ = rhs.model_;
84  generator_=rhs.generator_->clone();
85  generator_->refreshSolver(model_->solver());
86  whenCutGenerator_=rhs.whenCutGenerator_;
87  whenCutGeneratorInSub_ = rhs.whenCutGeneratorInSub_;
88  switchOffIfLessThan_ = rhs.switchOffIfLessThan_;
89  depthCutGenerator_=rhs.depthCutGenerator_;
90  depthCutGeneratorInSub_ = rhs.depthCutGeneratorInSub_;
91  generatorName_=strdup(rhs.generatorName_);
92  normal_=rhs.normal_;
93  atSolution_=rhs.atSolution_;
94  whenInfeasible_=rhs.whenInfeasible_;
95  mustCallAgain_ = rhs.mustCallAgain_;
96  switchedOff_ = rhs.switchedOff_;
97  timing_ = rhs.timing_;
98  timeInCutGenerator_ = rhs.timeInCutGenerator_;
99  numberTimes_ = rhs.numberTimes_;
100  numberCuts_ = rhs.numberCuts_;
101  numberColumnCuts_ = rhs.numberColumnCuts_;
102  numberCutsActive_ = rhs.numberCutsActive_;
103  numberCutsAtRoot_  = rhs.numberCutsAtRoot_;
104  numberActiveCutsAtRoot_ = rhs.numberActiveCutsAtRoot_;
105}
106
107// Assignment operator
108CbcCutGenerator & 
109CbcCutGenerator::operator=( const CbcCutGenerator& rhs)
110{
111  if (this!=&rhs) {
112    delete generator_;
113    free(generatorName_);
114    model_ = rhs.model_;
115    generator_=rhs.generator_->clone();
116    generator_->refreshSolver(model_->solver());
117    whenCutGenerator_=rhs.whenCutGenerator_;
118    whenCutGeneratorInSub_ = rhs.whenCutGeneratorInSub_;
119    switchOffIfLessThan_ = rhs.switchOffIfLessThan_;
120    depthCutGenerator_=rhs.depthCutGenerator_;
121    depthCutGeneratorInSub_ = rhs.depthCutGeneratorInSub_;
122    generatorName_=strdup(rhs.generatorName_);
123    normal_=rhs.normal_;
124    atSolution_=rhs.atSolution_;
125    whenInfeasible_=rhs.whenInfeasible_;
126    mustCallAgain_ = rhs.mustCallAgain_;
127    switchedOff_ = rhs.switchedOff_;
128    timing_ = rhs.timing_;
129    timeInCutGenerator_ = rhs.timeInCutGenerator_;
130    numberTimes_ = rhs.numberTimes_;
131    numberCuts_ = rhs.numberCuts_;
132    numberColumnCuts_ = rhs.numberColumnCuts_;
133    numberCutsActive_ = rhs.numberCutsActive_;
134    numberCutsAtRoot_  = rhs.numberCutsAtRoot_;
135    numberActiveCutsAtRoot_ = rhs.numberActiveCutsAtRoot_;
136  }
137  return *this;
138}
139
140// Destructor
141CbcCutGenerator::~CbcCutGenerator ()
142{
143  free(generatorName_);
144  delete generator_;
145}
146
147/* This is used to refresh any inforamtion.
148   It also refreshes the solver in the cut generator
149   in case generator wants to do some work
150*/
151void 
152CbcCutGenerator::refreshModel(CbcModel * model)
153{
154  model_=model;
155  generator_->refreshSolver(model_->solver());
156}
157/* Generate cuts for the model data contained in si.
158   The generated cuts are inserted into and returned in the
159   collection of cuts cs.
160*/
161bool
162CbcCutGenerator::generateCuts( OsiCuts & cs , bool fullScan, CbcNode * node)
163{
164  int howOften = whenCutGenerator_;
165  if (howOften==-100)
166    return false;
167  if (howOften>0)
168    howOften = howOften % 1000000;
169  else 
170    howOften=1;
171  if (!howOften)
172    howOften=1;
173  bool returnCode=false;
174  OsiSolverInterface * solver = model_->solver();
175  int depth;
176  if (node)
177    depth=node->depth();
178  else
179    depth=0;
180  int pass=model_->getCurrentPassNumber()-1;
181  bool doThis=(model_->getNodeCount()%howOften)==0;
182  if (depthCutGenerator_>0) {
183    doThis = (depth % depthCutGenerator_) ==0;
184    if (depth<depthCutGenerator_)
185      doThis=true; // and also at top of tree
186  }
187  // But turn off if 100
188  if (howOften==100)
189    doThis=false;
190  // Switch off if special setting
191  if (whenCutGeneratorInSub_==-200) {
192    fullScan=false;
193    doThis=false;
194  }
195  if (fullScan||doThis) {
196    double time1=0.0;
197    if (timing_)
198      time1 = CoinCpuTime();
199    //#define CBC_DEBUG
200    int numberRowCutsBefore = cs.sizeRowCuts() ;
201    int cutsBefore = cs.sizeCuts();
202    CglTreeInfo info;
203    info.level = depth;
204    info.pass = pass;
205    info.formulation_rows = model_->numberRowsAtContinuous();
206    info.inTree = node!=NULL;
207    incrementNumberTimesEntered();
208    CglProbing* generator =
209      dynamic_cast<CglProbing*>(generator_);
210    if (!generator) {
211      // Pass across model information in case it could be useful
212      //OsiSolverInterface * solver = model_->solver();
213      //void * saveData = solver->getApplicationData();
214      //solver->setApplicationData(model_);
215      generator_->generateCuts(*solver,cs,info);
216      //solver->setApplicationData(saveData);
217    } else {
218      // Probing - return tight column bounds
219      CglTreeProbingInfo * info2 = model_->probingInfo();
220      if (info2) {
221        info2->level = depth;
222        info2->pass = pass;
223        info2->formulation_rows = model_->numberRowsAtContinuous();
224        info2->inTree = node!=NULL;
225        generator->generateCutsAndModify(*solver,cs,info2);
226      } else {
227        generator->generateCutsAndModify(*solver,cs,&info);
228      }
229      const double * tightLower = generator->tightLower();
230      const double * lower = solver->getColLower();
231      const double * tightUpper = generator->tightUpper();
232      const double * upper = solver->getColUpper();
233      const double * solution = solver->getColSolution();
234      int j;
235      int numberColumns = solver->getNumCols();
236      double primalTolerance = 1.0e-8;
237#if 0
238      int numberChanged=0,ifCut=0;
239      CoinPackedVector lbs;
240      CoinPackedVector ubs;
241      for (j=0;j<numberColumns;j++) {
242        if (solver->isInteger(j)) {
243          if (tightUpper[j]<upper[j]) {
244            numberChanged++;
245            assert (tightUpper[j]==floor(tightUpper[j]+0.5));
246            ubs.insert(j,tightUpper[j]);
247            if (tightUpper[j]<solution[j]-primalTolerance)
248              ifCut=1;
249          }
250          if (tightLower[j]>lower[j]) {
251            numberChanged++;
252            assert (tightLower[j]==floor(tightLower[j]+0.5));
253            lbs.insert(j,tightLower[j]);
254            if (tightLower[j]>solution[j]+primalTolerance)
255              ifCut=1;
256          }
257        } else {
258          if (tightUpper[j]==tightLower[j]&&
259              upper[j]>lower[j]) {
260            // fix
261            //solver->setColLower(j,tightLower[j]);
262            //solver->setColUpper(j,tightUpper[j]);
263            double value = tightUpper[j];
264            numberChanged++;
265            if (value<upper[j])
266              ubs.insert(j,value);
267            if (value>lower[j])
268              lbs.insert(j,value);
269          }
270        }
271      }
272      if (numberChanged) {
273        OsiColCut cc;
274        cc.setUbs(ubs);
275        cc.setLbs(lbs);
276        if (ifCut) {
277          cc.setEffectiveness(100.0);
278        } else {
279          cc.setEffectiveness(1.0e-5);
280        }
281        cs.insert(cc);
282      }
283      // need to resolve if some bounds changed
284      returnCode = !solver->basisIsAvailable();
285      assert (!returnCode);
286#else
287      const char * tightenBounds = generator->tightenBounds();
288      for (j=0;j<numberColumns;j++) {
289        if (solver->isInteger(j)) {
290          if (tightUpper[j]<upper[j]) {
291            double nearest = floor(tightUpper[j]+0.5);
292            //assert (fabs(tightUpper[j]-nearest)<1.0e-5); may be infeasible
293            solver->setColUpper(j,nearest);
294            if (nearest<solution[j]-primalTolerance)
295              returnCode=true;
296          }
297          if (tightLower[j]>lower[j]) {
298            double nearest = floor(tightLower[j]+0.5);
299            //assert (fabs(tightLower[j]-nearest)<1.0e-5); may be infeasible
300            solver->setColLower(j,nearest);
301            if (nearest>solution[j]+primalTolerance)
302              returnCode=true;
303          }
304        } else {
305          if (upper[j]>lower[j]) {
306            if (tightUpper[j]==tightLower[j]) {
307              // fix
308              solver->setColLower(j,tightLower[j]);
309              solver->setColUpper(j,tightUpper[j]);
310              if (tightLower[j]>solution[j]+primalTolerance||
311                  tightUpper[j]<solution[j]-primalTolerance)
312                returnCode=true;
313            } else if (tightenBounds&&tightenBounds[j]) {
314              solver->setColLower(j,CoinMax(tightLower[j],lower[j]));
315              solver->setColUpper(j,CoinMin(tightUpper[j],upper[j]));
316              if (tightLower[j]>solution[j]+primalTolerance||
317                  tightUpper[j]<solution[j]-primalTolerance)
318                returnCode=true;
319            }
320          }
321        }
322      }
323      //if (!solver->basisIsAvailable())
324      //returnCode=true;
325#endif
326#if 0
327      // Pass across info to pseudocosts
328      char * mark = new char[numberColumns];
329      memset(mark,0,numberColumns);
330      int nLook = generator->numberThisTime();
331      const int * lookedAt = generator->lookedAt();
332      const int * fixedDown = generator->fixedDown();
333      const int * fixedUp = generator->fixedUp();
334      for (j=0;j<nLook;j++)
335        mark[lookedAt[j]]=1;
336      int numberObjects = model_->numberObjects();
337      for (int i=0;i<numberObjects;i++) {
338        CbcSimpleIntegerDynamicPseudoCost * obj1 =
339          dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(model_->modifiableObject(i)) ;
340        if (obj1) {
341          int iColumn = obj1->columnNumber();
342          if (mark[iColumn])
343            obj1->setProbingInformation(fixedDown[iColumn],fixedUp[iColumn]);
344        }
345      }
346      delete [] mark;
347#endif
348    }
349    CbcCutModifier * modifier = model_->cutModifier();
350    if (modifier) {
351      int numberRowCutsAfter = cs.sizeRowCuts() ;
352      int k ;
353      int nOdd=0;
354      const OsiSolverInterface * solver = model_->solver();
355      for (k = numberRowCutsAfter-1;k>=numberRowCutsBefore;k--) {
356        OsiRowCut & thisCut = cs.rowCut(k) ;
357        int returnCode = modifier->modify(solver,thisCut);
358        if (returnCode) {
359          nOdd++;
360          if (returnCode==3)
361            cs.eraseRowCut(k);
362        }
363      }
364      if (nOdd) 
365        printf("Cut generator %s produced %d cuts of which %d were modified\n",
366                 generatorName_,numberRowCutsAfter-numberRowCutsBefore,nOdd);
367    }
368#ifdef CBC_DEBUG
369    {
370      int numberRowCutsAfter = cs.sizeRowCuts() ;
371      int k ;
372      int nBad=0;
373      for (k = numberRowCutsBefore;k<numberRowCutsAfter;k++) {
374        OsiRowCut thisCut = cs.rowCut(k) ;
375        if (thisCut.lb()<=thisCut.ub()) {
376          /* check size of elements.
377             We can allow smaller but this helps debug generators as it
378             is unsafe to have small elements */
379          int n=thisCut.row().getNumElements();
380          const int * column = thisCut.row().getIndices();
381          const double * element = thisCut.row().getElements();
382          assert (n);
383          for (int i=0;i<n;i++) {
384            double value = element[i];
385            if (fabs(value)<=1.0e-12||fabs(value)>=1.0e20)
386              nBad++;
387          }
388        }
389        if (nBad) 
390          printf("Cut generator %s produced %d cuts of which %d had tiny or large elements\n",
391                 generatorName_,numberRowCutsAfter-numberRowCutsBefore,nBad);
392      }
393    }
394#endif
395    if (timing_)
396      timeInCutGenerator_ += CoinCpuTime()-time1;
397    // switch off if first time and no good
398    if (node==NULL&&!pass) {
399      if (cs.sizeCuts()-cutsBefore<CoinAbs(switchOffIfLessThan_)) {
400        whenCutGenerator_=-99;
401        whenCutGeneratorInSub_ = -200;
402      }
403    }
404  }
405  return returnCode;
406}
407void 
408CbcCutGenerator::setHowOften(int howOften) 
409{
410 
411  if (howOften>=1000000) {
412    // leave Probing every 10
413    howOften = howOften % 1000000;
414    CglProbing* generator =
415      dynamic_cast<CglProbing*>(generator_);
416   
417    if (generator&&howOften>10) 
418      howOften=10+1000000;
419    else
420      howOften += 1000000;
421  }
422  whenCutGenerator_ = howOften;
423}
424void 
425CbcCutGenerator::setWhatDepth(int value) 
426{
427  depthCutGenerator_ = value;
428}
429void 
430CbcCutGenerator::setWhatDepthInSub(int value) 
431{
432  depthCutGeneratorInSub_ = value;
433}
434
435
436// Default Constructor
437CbcCutModifier::CbcCutModifier() 
438{
439}
440
441
442// Destructor
443CbcCutModifier::~CbcCutModifier ()
444{
445}
446
447// Copy constructor
448CbcCutModifier::CbcCutModifier ( const CbcCutModifier & rhs)
449{
450}
451
452// Assignment operator
453CbcCutModifier & 
454CbcCutModifier::operator=( const CbcCutModifier& rhs)
455{
456  if (this!=&rhs) {
457  }
458  return *this;
459}
460
461// Default Constructor
462CbcCutSubsetModifier::CbcCutSubsetModifier ()
463  : CbcCutModifier(),
464    firstOdd_(INT_MAX)
465{
466}
467
468// Useful constructor
469CbcCutSubsetModifier::CbcCutSubsetModifier (int firstOdd)
470  : CbcCutModifier()
471{
472  firstOdd_=firstOdd;
473}
474
475// Copy constructor
476CbcCutSubsetModifier::CbcCutSubsetModifier ( const CbcCutSubsetModifier & rhs)
477  :CbcCutModifier(rhs)
478{
479  firstOdd_ = rhs.firstOdd_;
480}
481
482// Clone
483CbcCutModifier *
484CbcCutSubsetModifier::clone() const
485{
486  return new CbcCutSubsetModifier(*this);
487}
488
489// Assignment operator
490CbcCutSubsetModifier & 
491CbcCutSubsetModifier::operator=( const CbcCutSubsetModifier& rhs)
492{
493  if (this!=&rhs) {
494    CbcCutModifier::operator=(rhs);
495    firstOdd_ = rhs.firstOdd_;
496  }
497  return *this;
498}
499
500// Destructor
501CbcCutSubsetModifier::~CbcCutSubsetModifier ()
502{
503}
504/* Returns
505   0 unchanged
506   1 strengthened
507   2 weakened
508   3 deleted
509*/
510int 
511CbcCutSubsetModifier::modify(const OsiSolverInterface * solver, OsiRowCut & cut) 
512{
513  int n=cut.row().getNumElements();
514  if (!n)
515    return 0;
516  const int * column = cut.row().getIndices();
517  //const double * element = cut.row().getElements();
518  int returnCode=0;
519  for (int i=0;i<n;i++) {
520    if (column[i]>=firstOdd_) {
521      returnCode=3;
522      break;
523    }
524  }
525  if (!returnCode) {
526    const double * element = cut.row().getElements();
527    printf("%g <= ",cut.lb());
528    for (int i=0;i<n;i++) {
529      printf("%g*x%d ",element[i],column[i]);
530    }
531    printf("<= %g\n",cut.ub());
532  }
533  //return 3;
534  return returnCode;
535}
536
Note: See TracBrowser for help on using the repository browser.