source: stable/2.0/Cbc/src/CbcCutGenerator.cpp @ 905

Last change on this file since 905 was 905, checked in by ladanyi, 13 years ago

include cstdlib before cmath to get things to compile on AIX with xlC (same as changeset 904 in trunk)

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