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

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

update branches/devel for threads

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