source: branches/heur/Cbc/src/CbcCutGenerator.cpp @ 885

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

for deterministic parallel

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