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

Last change on this file since 648 was 648, checked in by forrest, 12 years ago

for parallel

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 15.2 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, OsiSolverInterface * solver, 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      //void * saveData = solver->getApplicationData();
213      //solver->setApplicationData(model_);
214      generator_->generateCuts(*solver,cs,info);
215      //solver->setApplicationData(saveData);
216    } else {
217      // Probing - return tight column bounds
218      CglTreeProbingInfo * info2 = model_->probingInfo();
219      if (info2) {
220        info2->level = depth;
221        info2->pass = pass;
222        info2->formulation_rows = model_->numberRowsAtContinuous();
223        info2->inTree = node!=NULL;
224        generator->generateCutsAndModify(*solver,cs,info2);
225      } else {
226        generator->generateCutsAndModify(*solver,cs,&info);
227      }
228      const double * tightLower = generator->tightLower();
229      const double * lower = solver->getColLower();
230      const double * tightUpper = generator->tightUpper();
231      const double * upper = solver->getColUpper();
232      const double * solution = solver->getColSolution();
233      int j;
234      int numberColumns = solver->getNumCols();
235      double primalTolerance = 1.0e-8;
236      const char * tightenBounds = generator->tightenBounds();
237      if ((model_->getThreadMode()&2)==0) {
238        for (j=0;j<numberColumns;j++) {
239          if (solver->isInteger(j)) {
240            if (tightUpper[j]<upper[j]) {
241              double nearest = floor(tightUpper[j]+0.5);
242              //assert (fabs(tightUpper[j]-nearest)<1.0e-5); may be infeasible
243              solver->setColUpper(j,nearest);
244              if (nearest<solution[j]-primalTolerance)
245                returnCode=true;
246            }
247            if (tightLower[j]>lower[j]) {
248              double nearest = floor(tightLower[j]+0.5);
249              //assert (fabs(tightLower[j]-nearest)<1.0e-5); may be infeasible
250              solver->setColLower(j,nearest);
251              if (nearest>solution[j]+primalTolerance)
252                returnCode=true;
253            }
254          } else {
255            if (upper[j]>lower[j]) {
256              if (tightUpper[j]==tightLower[j]) {
257                // fix
258                solver->setColLower(j,tightLower[j]);
259                solver->setColUpper(j,tightUpper[j]);
260                if (tightLower[j]>solution[j]+primalTolerance||
261                    tightUpper[j]<solution[j]-primalTolerance)
262                  returnCode=true;
263              } else if (tightenBounds&&tightenBounds[j]) {
264                solver->setColLower(j,CoinMax(tightLower[j],lower[j]));
265                solver->setColUpper(j,CoinMin(tightUpper[j],upper[j]));
266                if (tightLower[j]>solution[j]+primalTolerance||
267                    tightUpper[j]<solution[j]-primalTolerance)
268                  returnCode=true;
269              }
270            }
271          }
272        }
273      } else {
274        CoinPackedVector lbs;
275        CoinPackedVector ubs;
276        int numberChanged=0;
277        bool ifCut=false;
278        for (j=0;j<numberColumns;j++) {
279          if (solver->isInteger(j)) {
280            if (tightUpper[j]<upper[j]) {
281              double nearest = floor(tightUpper[j]+0.5);
282              //assert (fabs(tightUpper[j]-nearest)<1.0e-5); may be infeasible
283              ubs.insert(j,nearest);
284              numberChanged++;
285              if (nearest<solution[j]-primalTolerance)
286                ifCut=true;
287            }
288            if (tightLower[j]>lower[j]) {
289              double nearest = floor(tightLower[j]+0.5);
290              //assert (fabs(tightLower[j]-nearest)<1.0e-5); may be infeasible
291              lbs.insert(j,nearest);
292              numberChanged++;
293              if (nearest>solution[j]+primalTolerance)
294                ifCut=true;
295            }
296          } else {
297            if (upper[j]>lower[j]) {
298              if (tightUpper[j]==tightLower[j]) {
299                // fix
300                lbs.insert(j,tightLower[j]);
301                ubs.insert(j,tightUpper[j]);
302                if (tightLower[j]>solution[j]+primalTolerance||
303                    tightUpper[j]<solution[j]-primalTolerance)
304                  ifCut=true;
305              } else if (tightenBounds&&tightenBounds[j]) {
306                lbs.insert(j,CoinMax(tightLower[j],lower[j]));
307                ubs.insert(j,CoinMin(tightUpper[j],upper[j]));
308                if (tightLower[j]>solution[j]+primalTolerance||
309                    tightUpper[j]<solution[j]-primalTolerance)
310                  ifCut=true;
311              }
312            }
313          }
314        }
315        if (numberChanged) {
316          OsiColCut cc;
317          cc.setUbs(ubs);
318          cc.setLbs(lbs);
319          if (ifCut) {
320            cc.setEffectiveness(100.0);
321          } else {
322            cc.setEffectiveness(1.0e-5);
323          }
324          cs.insert(cc);
325        }
326      }
327      //if (!solver->basisIsAvailable())
328      //returnCode=true;
329#if 0
330      // Pass across info to pseudocosts
331      char * mark = new char[numberColumns];
332      memset(mark,0,numberColumns);
333      int nLook = generator->numberThisTime();
334      const int * lookedAt = generator->lookedAt();
335      const int * fixedDown = generator->fixedDown();
336      const int * fixedUp = generator->fixedUp();
337      for (j=0;j<nLook;j++)
338        mark[lookedAt[j]]=1;
339      int numberObjects = model_->numberObjects();
340      for (int i=0;i<numberObjects;i++) {
341        CbcSimpleIntegerDynamicPseudoCost * obj1 =
342          dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(model_->modifiableObject(i)) ;
343        if (obj1) {
344          int iColumn = obj1->columnNumber();
345          if (mark[iColumn])
346            obj1->setProbingInformation(fixedDown[iColumn],fixedUp[iColumn]);
347        }
348      }
349      delete [] mark;
350#endif
351    }
352    CbcCutModifier * modifier = model_->cutModifier();
353    if (modifier) {
354      int numberRowCutsAfter = cs.sizeRowCuts() ;
355      int k ;
356      int nOdd=0;
357      //const OsiSolverInterface * solver = model_->solver();
358      for (k = numberRowCutsAfter-1;k>=numberRowCutsBefore;k--) {
359        OsiRowCut & thisCut = cs.rowCut(k) ;
360        int returnCode = modifier->modify(solver,thisCut);
361        if (returnCode) {
362          nOdd++;
363          if (returnCode==3)
364            cs.eraseRowCut(k);
365        }
366      }
367      if (nOdd) 
368        printf("Cut generator %s produced %d cuts of which %d were modified\n",
369                 generatorName_,numberRowCutsAfter-numberRowCutsBefore,nOdd);
370    }
371    { 
372      // make all row cuts without test for duplicate
373      int numberRowCutsAfter = cs.sizeRowCuts() ;
374      int k ;
375      for (k = numberRowCutsBefore;k<numberRowCutsAfter;k++) {
376        OsiRowCut * thisCut = cs.rowCutPtr(k) ;
377        thisCut->mutableRow().setTestForDuplicateIndex(false);
378      }
379    }
380#ifdef CBC_DEBUG
381    {
382      int numberRowCutsAfter = cs.sizeRowCuts() ;
383      int k ;
384      int nBad=0;
385      for (k = numberRowCutsBefore;k<numberRowCutsAfter;k++) {
386        OsiRowCut thisCut = cs.rowCut(k) ;
387        if (thisCut.lb()>thisCut.ub()||
388            thisCut.lb()>1.0e8||
389            thisCut.ub()<-1.0e8)
390          printf("cut from %s has bounds %g and %g!\n",
391                 generatorName_,thisCut.lb(),thisCut.ub());
392        if (thisCut.lb()<=thisCut.ub()) {
393          /* check size of elements.
394             We can allow smaller but this helps debug generators as it
395             is unsafe to have small elements */
396          int n=thisCut.row().getNumElements();
397          const int * column = thisCut.row().getIndices();
398          const double * element = thisCut.row().getElements();
399          assert (n);
400          for (int i=0;i<n;i++) {
401            double value = element[i];
402            if (fabs(value)<=1.0e-12||fabs(value)>=1.0e20)
403              nBad++;
404          }
405        }
406        if (nBad) 
407          printf("Cut generator %s produced %d cuts of which %d had tiny or large elements\n",
408                 generatorName_,numberRowCutsAfter-numberRowCutsBefore,nBad);
409      }
410    }
411#endif
412    if (timing_)
413      timeInCutGenerator_ += CoinCpuTime()-time1;
414    // switch off if first time and no good
415    if (node==NULL&&!pass) {
416      if (cs.sizeCuts()-cutsBefore<CoinAbs(switchOffIfLessThan_)) {
417        whenCutGenerator_=-99;
418        whenCutGeneratorInSub_ = -200;
419      }
420    }
421  }
422  return returnCode;
423}
424void 
425CbcCutGenerator::setHowOften(int howOften) 
426{
427 
428  if (howOften>=1000000) {
429    // leave Probing every 10
430    howOften = howOften % 1000000;
431    CglProbing* generator =
432      dynamic_cast<CglProbing*>(generator_);
433   
434    if (generator&&howOften>10) 
435      howOften=10+1000000;
436    else
437      howOften += 1000000;
438  }
439  whenCutGenerator_ = howOften;
440}
441void 
442CbcCutGenerator::setWhatDepth(int value) 
443{
444  depthCutGenerator_ = value;
445}
446void 
447CbcCutGenerator::setWhatDepthInSub(int value) 
448{
449  depthCutGeneratorInSub_ = value;
450}
451
452
453// Default Constructor
454CbcCutModifier::CbcCutModifier() 
455{
456}
457
458
459// Destructor
460CbcCutModifier::~CbcCutModifier ()
461{
462}
463
464// Copy constructor
465CbcCutModifier::CbcCutModifier ( const CbcCutModifier & rhs)
466{
467}
468
469// Assignment operator
470CbcCutModifier & 
471CbcCutModifier::operator=( const CbcCutModifier& rhs)
472{
473  if (this!=&rhs) {
474  }
475  return *this;
476}
477
478// Default Constructor
479CbcCutSubsetModifier::CbcCutSubsetModifier ()
480  : CbcCutModifier(),
481    firstOdd_(INT_MAX)
482{
483}
484
485// Useful constructor
486CbcCutSubsetModifier::CbcCutSubsetModifier (int firstOdd)
487  : CbcCutModifier()
488{
489  firstOdd_=firstOdd;
490}
491
492// Copy constructor
493CbcCutSubsetModifier::CbcCutSubsetModifier ( const CbcCutSubsetModifier & rhs)
494  :CbcCutModifier(rhs)
495{
496  firstOdd_ = rhs.firstOdd_;
497}
498
499// Clone
500CbcCutModifier *
501CbcCutSubsetModifier::clone() const
502{
503  return new CbcCutSubsetModifier(*this);
504}
505
506// Assignment operator
507CbcCutSubsetModifier & 
508CbcCutSubsetModifier::operator=( const CbcCutSubsetModifier& rhs)
509{
510  if (this!=&rhs) {
511    CbcCutModifier::operator=(rhs);
512    firstOdd_ = rhs.firstOdd_;
513  }
514  return *this;
515}
516
517// Destructor
518CbcCutSubsetModifier::~CbcCutSubsetModifier ()
519{
520}
521/* Returns
522   0 unchanged
523   1 strengthened
524   2 weakened
525   3 deleted
526*/
527int 
528CbcCutSubsetModifier::modify(const OsiSolverInterface * solver, OsiRowCut & cut) 
529{
530  int n=cut.row().getNumElements();
531  if (!n)
532    return 0;
533  const int * column = cut.row().getIndices();
534  //const double * element = cut.row().getElements();
535  int returnCode=0;
536  for (int i=0;i<n;i++) {
537    if (column[i]>=firstOdd_) {
538      returnCode=3;
539      break;
540    }
541  }
542  if (!returnCode) {
543    const double * element = cut.row().getElements();
544    printf("%g <= ",cut.lb());
545    for (int i=0;i<n;i++) {
546      printf("%g*x%d ",element[i],column[i]);
547    }
548    printf("<= %g\n",cut.ub());
549  }
550  //return 3;
551  return returnCode;
552}
553
Note: See TracBrowser for help on using the repository browser.