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

Last change on this file since 539 was 529, checked in by forrest, 13 years ago

for nonlinear

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 15.1 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    { 
369      // make all row cuts without test for duplicate
370      int numberRowCutsAfter = cs.sizeRowCuts() ;
371      int k ;
372      for (k = numberRowCutsBefore;k<numberRowCutsAfter;k++) {
373        OsiRowCut * thisCut = cs.rowCutPtr(k) ;
374        thisCut->mutableRow().setTestForDuplicateIndex(false);
375      }
376    }
377#ifdef CBC_DEBUG
378    {
379      int numberRowCutsAfter = cs.sizeRowCuts() ;
380      int k ;
381      int nBad=0;
382      for (k = numberRowCutsBefore;k<numberRowCutsAfter;k++) {
383        OsiRowCut thisCut = cs.rowCut(k) ;
384        if (thisCut.lb()<=thisCut.ub()) {
385          /* check size of elements.
386             We can allow smaller but this helps debug generators as it
387             is unsafe to have small elements */
388          int n=thisCut.row().getNumElements();
389          const int * column = thisCut.row().getIndices();
390          const double * element = thisCut.row().getElements();
391          assert (n);
392          for (int i=0;i<n;i++) {
393            double value = element[i];
394            if (fabs(value)<=1.0e-12||fabs(value)>=1.0e20)
395              nBad++;
396          }
397        }
398        if (nBad) 
399          printf("Cut generator %s produced %d cuts of which %d had tiny or large elements\n",
400                 generatorName_,numberRowCutsAfter-numberRowCutsBefore,nBad);
401      }
402    }
403#endif
404    if (timing_)
405      timeInCutGenerator_ += CoinCpuTime()-time1;
406    // switch off if first time and no good
407    if (node==NULL&&!pass) {
408      if (cs.sizeCuts()-cutsBefore<CoinAbs(switchOffIfLessThan_)) {
409        whenCutGenerator_=-99;
410        whenCutGeneratorInSub_ = -200;
411      }
412    }
413  }
414  return returnCode;
415}
416void 
417CbcCutGenerator::setHowOften(int howOften) 
418{
419 
420  if (howOften>=1000000) {
421    // leave Probing every 10
422    howOften = howOften % 1000000;
423    CglProbing* generator =
424      dynamic_cast<CglProbing*>(generator_);
425   
426    if (generator&&howOften>10) 
427      howOften=10+1000000;
428    else
429      howOften += 1000000;
430  }
431  whenCutGenerator_ = howOften;
432}
433void 
434CbcCutGenerator::setWhatDepth(int value) 
435{
436  depthCutGenerator_ = value;
437}
438void 
439CbcCutGenerator::setWhatDepthInSub(int value) 
440{
441  depthCutGeneratorInSub_ = value;
442}
443
444
445// Default Constructor
446CbcCutModifier::CbcCutModifier() 
447{
448}
449
450
451// Destructor
452CbcCutModifier::~CbcCutModifier ()
453{
454}
455
456// Copy constructor
457CbcCutModifier::CbcCutModifier ( const CbcCutModifier & rhs)
458{
459}
460
461// Assignment operator
462CbcCutModifier & 
463CbcCutModifier::operator=( const CbcCutModifier& rhs)
464{
465  if (this!=&rhs) {
466  }
467  return *this;
468}
469
470// Default Constructor
471CbcCutSubsetModifier::CbcCutSubsetModifier ()
472  : CbcCutModifier(),
473    firstOdd_(INT_MAX)
474{
475}
476
477// Useful constructor
478CbcCutSubsetModifier::CbcCutSubsetModifier (int firstOdd)
479  : CbcCutModifier()
480{
481  firstOdd_=firstOdd;
482}
483
484// Copy constructor
485CbcCutSubsetModifier::CbcCutSubsetModifier ( const CbcCutSubsetModifier & rhs)
486  :CbcCutModifier(rhs)
487{
488  firstOdd_ = rhs.firstOdd_;
489}
490
491// Clone
492CbcCutModifier *
493CbcCutSubsetModifier::clone() const
494{
495  return new CbcCutSubsetModifier(*this);
496}
497
498// Assignment operator
499CbcCutSubsetModifier & 
500CbcCutSubsetModifier::operator=( const CbcCutSubsetModifier& rhs)
501{
502  if (this!=&rhs) {
503    CbcCutModifier::operator=(rhs);
504    firstOdd_ = rhs.firstOdd_;
505  }
506  return *this;
507}
508
509// Destructor
510CbcCutSubsetModifier::~CbcCutSubsetModifier ()
511{
512}
513/* Returns
514   0 unchanged
515   1 strengthened
516   2 weakened
517   3 deleted
518*/
519int 
520CbcCutSubsetModifier::modify(const OsiSolverInterface * solver, OsiRowCut & cut) 
521{
522  int n=cut.row().getNumElements();
523  if (!n)
524    return 0;
525  const int * column = cut.row().getIndices();
526  //const double * element = cut.row().getElements();
527  int returnCode=0;
528  for (int i=0;i<n;i++) {
529    if (column[i]>=firstOdd_) {
530      returnCode=3;
531      break;
532    }
533  }
534  if (!returnCode) {
535    const double * element = cut.row().getElements();
536    printf("%g <= ",cut.lb());
537    for (int i=0;i<n;i++) {
538      printf("%g*x%d ",element[i],column[i]);
539    }
540    printf("<= %g\n",cut.ub());
541  }
542  //return 3;
543  return returnCode;
544}
545
Note: See TracBrowser for help on using the repository browser.