source: trunk/Cbc/src/CbcCutGenerator.cpp @ 1015

Last change on this file since 1015 was 1015, checked in by forrest, 11 years ago

add ifdefs for future exploration

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 24.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 <cstdlib>
10#include <cmath>
11#include <cfloat>
12
13#ifdef COIN_HAS_CLP
14#include "OsiClpSolverInterface.hpp"
15#else
16#include "OsiSolverInterface.hpp"
17#endif
18#include "CbcModel.hpp"
19#include "CbcMessage.hpp"
20#include "CbcCutGenerator.hpp"
21#include "CbcBranchDynamic.hpp"
22#include "CglProbing.hpp"
23#include "CoinTime.hpp"
24
25// Default Constructor
26CbcCutGenerator::CbcCutGenerator ()
27  : model_(NULL),
28    generator_(NULL),
29    whenCutGenerator_(-1),
30    whenCutGeneratorInSub_(-100),
31    switchOffIfLessThan_(0),
32    depthCutGenerator_(-1),
33    depthCutGeneratorInSub_(-1),
34    generatorName_(NULL),
35    normal_(true),
36    atSolution_(false),
37    whenInfeasible_(false),
38    mustCallAgain_(false),
39    switchedOff_(false),
40    timing_(false),
41    timeInCutGenerator_(0.0),
42    numberTimes_(0),
43    numberCuts_(0),
44    numberColumnCuts_(0),
45    numberCutsActive_(0),
46    numberCutsAtRoot_(0),
47    numberActiveCutsAtRoot_(0)
48{
49}
50// Normal constructor
51CbcCutGenerator::CbcCutGenerator(CbcModel * model,CglCutGenerator * generator,
52                                 int howOften, const char * name,
53                                 bool normal, bool atSolution, 
54                                 bool infeasible, int howOftenInSub,
55                                 int whatDepth, int whatDepthInSub,
56                                 int switchOffIfLessThan)
57  : 
58    depthCutGenerator_(whatDepth),
59    depthCutGeneratorInSub_(whatDepthInSub),
60    mustCallAgain_(false),
61    switchedOff_(false),
62    timing_(false),
63    timeInCutGenerator_(0.0),
64    numberTimes_(0),
65    numberCuts_(0),
66    numberColumnCuts_(0),
67    numberCutsActive_(0),
68    numberCutsAtRoot_(0),
69    numberActiveCutsAtRoot_(0)
70{
71  model_ = model;
72  generator_=generator->clone();
73  generator_->refreshSolver(model_->solver());
74  whenCutGenerator_=howOften;
75  whenCutGeneratorInSub_ = howOftenInSub;
76  switchOffIfLessThan_=switchOffIfLessThan;
77  if (name)
78    generatorName_=strdup(name);
79  else
80    generatorName_ = strdup("Unknown");
81  normal_=normal;
82  atSolution_=atSolution;
83  whenInfeasible_=infeasible;
84}
85
86// Copy constructor
87CbcCutGenerator::CbcCutGenerator ( const CbcCutGenerator & rhs)
88{
89  model_ = rhs.model_;
90  generator_=rhs.generator_->clone();
91  //generator_->refreshSolver(model_->solver());
92  whenCutGenerator_=rhs.whenCutGenerator_;
93  whenCutGeneratorInSub_ = rhs.whenCutGeneratorInSub_;
94  switchOffIfLessThan_ = rhs.switchOffIfLessThan_;
95  depthCutGenerator_=rhs.depthCutGenerator_;
96  depthCutGeneratorInSub_ = rhs.depthCutGeneratorInSub_;
97  generatorName_=strdup(rhs.generatorName_);
98  normal_=rhs.normal_;
99  atSolution_=rhs.atSolution_;
100  whenInfeasible_=rhs.whenInfeasible_;
101  mustCallAgain_ = rhs.mustCallAgain_;
102  switchedOff_ = rhs.switchedOff_;
103  timing_ = rhs.timing_;
104  timeInCutGenerator_ = rhs.timeInCutGenerator_;
105  numberTimes_ = rhs.numberTimes_;
106  numberCuts_ = rhs.numberCuts_;
107  numberColumnCuts_ = rhs.numberColumnCuts_;
108  numberCutsActive_ = rhs.numberCutsActive_;
109  numberCutsAtRoot_  = rhs.numberCutsAtRoot_;
110  numberActiveCutsAtRoot_ = rhs.numberActiveCutsAtRoot_;
111}
112
113// Assignment operator
114CbcCutGenerator & 
115CbcCutGenerator::operator=( const CbcCutGenerator& rhs)
116{
117  if (this!=&rhs) {
118    delete generator_;
119    free(generatorName_);
120    model_ = rhs.model_;
121    generator_=rhs.generator_->clone();
122    generator_->refreshSolver(model_->solver());
123    whenCutGenerator_=rhs.whenCutGenerator_;
124    whenCutGeneratorInSub_ = rhs.whenCutGeneratorInSub_;
125    switchOffIfLessThan_ = rhs.switchOffIfLessThan_;
126    depthCutGenerator_=rhs.depthCutGenerator_;
127    depthCutGeneratorInSub_ = rhs.depthCutGeneratorInSub_;
128    generatorName_=strdup(rhs.generatorName_);
129    normal_=rhs.normal_;
130    atSolution_=rhs.atSolution_;
131    whenInfeasible_=rhs.whenInfeasible_;
132    mustCallAgain_ = rhs.mustCallAgain_;
133    switchedOff_ = rhs.switchedOff_;
134    timing_ = rhs.timing_;
135    timeInCutGenerator_ = rhs.timeInCutGenerator_;
136    numberTimes_ = rhs.numberTimes_;
137    numberCuts_ = rhs.numberCuts_;
138    numberColumnCuts_ = rhs.numberColumnCuts_;
139    numberCutsActive_ = rhs.numberCutsActive_;
140    numberCutsAtRoot_  = rhs.numberCutsAtRoot_;
141    numberActiveCutsAtRoot_ = rhs.numberActiveCutsAtRoot_;
142  }
143  return *this;
144}
145
146// Destructor
147CbcCutGenerator::~CbcCutGenerator ()
148{
149  free(generatorName_);
150  delete generator_;
151}
152
153/* This is used to refresh any inforamtion.
154   It also refreshes the solver in the cut generator
155   in case generator wants to do some work
156*/
157void 
158CbcCutGenerator::refreshModel(CbcModel * model)
159{
160  model_=model;
161  generator_->refreshSolver(model_->solver());
162}
163/* Generate cuts for the model data contained in si.
164   The generated cuts are inserted into and returned in the
165   collection of cuts cs.
166*/
167bool
168CbcCutGenerator::generateCuts( OsiCuts & cs , int fullScan, OsiSolverInterface * solver, CbcNode * node)
169{
170#define PROBE1 0
171#define PROBE2 0
172#define PROBE3 0
173  int depth;
174  if (node)
175    depth=node->depth();
176  else
177    depth=0;
178  int howOften = whenCutGenerator_;
179  if (dynamic_cast<CglProbing*>(generator_)&&PROBE1) {
180    if (howOften==-100&&model_->doCutsNow(3)) {
181      howOften=1; // do anyway
182    }
183  }
184  if (howOften==-100)
185    return false;
186  if (howOften>0)
187    howOften = howOften % 1000000;
188  else 
189    howOften=1;
190  if (!howOften)
191    howOften=1;
192  bool returnCode=false;
193  //OsiSolverInterface * solver = model_->solver();
194  int pass=model_->getCurrentPassNumber()-1;
195  bool doThis=(model_->getNodeCount()%howOften)==0;
196  CoinThreadRandom * randomNumberGenerator=NULL;
197#ifdef COIN_HAS_CLP
198  {
199    OsiClpSolverInterface * clpSolver
200      = dynamic_cast<OsiClpSolverInterface *> (solver);
201    if (clpSolver) 
202      randomNumberGenerator = clpSolver->getModelPtr()->randomNumberGenerator();
203  }
204#endif
205  if (depthCutGenerator_>0) {
206    doThis = (depth % depthCutGenerator_) ==0;
207    if (depth<depthCutGenerator_)
208      doThis=true; // and also at top of tree
209  }
210  // But turn off if 100
211  if (howOften==100)
212    doThis=false;
213  // Switch off if special setting
214  if (whenCutGeneratorInSub_==-200) {
215    fullScan=0;
216    doThis=false;
217  }
218  if (fullScan||doThis) {
219    double time1=0.0;
220    if (timing_)
221      time1 = CoinCpuTime();
222    //#define CBC_DEBUG
223    int numberRowCutsBefore = cs.sizeRowCuts() ;
224    int numberColumnCutsBefore = cs.sizeColCuts() ;
225#if 0
226    int cutsBefore = cs.sizeCuts();
227#endif
228    CglTreeInfo info;
229    info.level = depth;
230    info.pass = pass;
231    info.formulation_rows = model_->numberRowsAtContinuous();
232    info.inTree = node!=NULL;
233    info.randomNumberGenerator=randomNumberGenerator;
234    incrementNumberTimesEntered();
235    CglProbing* generator =
236      dynamic_cast<CglProbing*>(generator_);
237    if (!generator) {
238      // Pass across model information in case it could be useful
239      //void * saveData = solver->getApplicationData();
240      //solver->setApplicationData(model_);
241      generator_->generateCuts(*solver,cs,info);
242      //solver->setApplicationData(saveData);
243    } else {
244      // Probing - return tight column bounds
245      CglTreeProbingInfo * info2 = model_->probingInfo();
246      bool doCuts=false;
247      if (info2&&!depth) {
248        info2->level = depth;
249        info2->pass = pass;
250        info2->formulation_rows = model_->numberRowsAtContinuous();
251        info2->inTree = node!=NULL;
252        info2->randomNumberGenerator=randomNumberGenerator;
253        generator->generateCutsAndModify(*solver,cs,info2);
254        doCuts=true;
255      } else if (depth||PROBE2) {
256#if PROBE3
257        if ((numberTimes_==200||(numberTimes_>200&&(numberTimes_%2000)==0))
258             &&!model_->parentModel()) {
259          // in tree, maxStack, maxProbe
260          int test[]= {
261            100123,
262            199999,
263            200123,
264            299999,
265            500123,
266            599999,
267            1000123,
268            1099999,
269            2000123,
270            2099999};
271          int n = (int) (sizeof(test)/sizeof(int));
272          int saveStack = generator->getMaxLook();
273          int saveNumber = generator->getMaxProbe();
274#ifdef CLP_INVESTIGATE
275          int kr1=0;
276          int kc1=0;
277#endif
278          int bestStackTree=-1;
279          int bestNumberTree=-1;
280          for (int i=0;i<n;i++) {
281            OsiCuts cs2 = cs;
282            int stack = test[i]/100000;
283            int number = test[i] - 100000*stack;
284            generator->setMaxLook(stack);
285            generator->setMaxProbe(number);
286            generator_->generateCuts(*solver,cs2,info);
287#ifdef CLP_INVESTIGATE
288            int numberRowCuts = cs2.sizeRowCuts()-numberRowCutsBefore ;
289            int numberColumnCuts= cs2.sizeColCuts()-numberColumnCutsBefore ;
290            if (numberRowCuts<kr1||numberColumnCuts<kc1)
291              printf("Odd ");
292            if (numberRowCuts>kr1||numberColumnCuts>kc1) {
293              printf("*** ");
294              kr1=numberRowCuts;
295              kc1=numberColumnCuts;
296              bestStackTree=stack;
297              bestNumberTree=number;
298            }
299            printf("maxStack %d number %d gives %d row cuts and %d column cuts\n",
300                   stack,number,numberRowCuts,numberColumnCuts);
301#endif
302          }
303          generator->setMaxLook(saveStack);
304          generator->setMaxProbe(saveNumber);
305          if (bestStackTree>0) {
306            generator->setMaxLook(bestStackTree);
307            generator->setMaxProbe(bestNumberTree);
308#ifdef CLP_INVESTIGATE
309            printf("RRNumber %d -> %d, stack %d -> %d\n",
310                   saveNumber,bestNumberTree,saveStack,bestStackTree);
311#endif
312          } else {
313            // no good
314            generator->setMaxLook(0);
315#ifdef CLP_INVESTIGATE
316            printf("RRSwitching off number %d -> %d, stack %d -> %d\n",
317                   saveNumber,saveNumber,saveStack,1);
318#endif
319          }
320        }
321#endif
322        if (generator->getMaxLook()>0) {
323          generator->generateCutsAndModify(*solver,cs,&info);
324          doCuts=true;
325        }
326      } else {
327        // at root - don't always do
328        if (!PROBE3||pass<15||(pass&1)==0) {
329          generator->generateCutsAndModify(*solver,cs,&info);
330          doCuts=true;
331        }
332      }
333      if (doCuts) {
334        const double * tightLower = generator->tightLower();
335        const double * lower = solver->getColLower();
336        const double * tightUpper = generator->tightUpper();
337        const double * upper = solver->getColUpper();
338        const double * solution = solver->getColSolution();
339        int j;
340        int numberColumns = solver->getNumCols();
341        double primalTolerance = 1.0e-8;
342        const char * tightenBounds = generator->tightenBounds();
343        if ((model_->getThreadMode()&2)==0) {
344          for (j=0;j<numberColumns;j++) {
345            if (solver->isInteger(j)) {
346              if (tightUpper[j]<upper[j]) {
347                double nearest = floor(tightUpper[j]+0.5);
348                //assert (fabs(tightUpper[j]-nearest)<1.0e-5); may be infeasible
349                solver->setColUpper(j,nearest);
350                if (nearest<solution[j]-primalTolerance)
351                  returnCode=true;
352              }
353              if (tightLower[j]>lower[j]) {
354                double nearest = floor(tightLower[j]+0.5);
355                //assert (fabs(tightLower[j]-nearest)<1.0e-5); may be infeasible
356                solver->setColLower(j,nearest);
357                if (nearest>solution[j]+primalTolerance)
358                  returnCode=true;
359              }
360            } else {
361              if (upper[j]>lower[j]) {
362                if (tightUpper[j]==tightLower[j]) {
363                  // fix
364                  solver->setColLower(j,tightLower[j]);
365                  solver->setColUpper(j,tightUpper[j]);
366                  if (tightLower[j]>solution[j]+primalTolerance||
367                      tightUpper[j]<solution[j]-primalTolerance)
368                    returnCode=true;
369                } else if (tightenBounds&&tightenBounds[j]) {
370                  solver->setColLower(j,CoinMax(tightLower[j],lower[j]));
371                  solver->setColUpper(j,CoinMin(tightUpper[j],upper[j]));
372                  if (tightLower[j]>solution[j]+primalTolerance||
373                      tightUpper[j]<solution[j]-primalTolerance)
374                    returnCode=true;
375                }
376              }
377            }
378          }
379        } else {
380          CoinPackedVector lbs;
381          CoinPackedVector ubs;
382          int numberChanged=0;
383          bool ifCut=false;
384          for (j=0;j<numberColumns;j++) {
385            if (solver->isInteger(j)) {
386              if (tightUpper[j]<upper[j]) {
387                double nearest = floor(tightUpper[j]+0.5);
388                //assert (fabs(tightUpper[j]-nearest)<1.0e-5); may be infeasible
389                ubs.insert(j,nearest);
390                numberChanged++;
391                if (nearest<solution[j]-primalTolerance)
392                  ifCut=true;
393              }
394              if (tightLower[j]>lower[j]) {
395                double nearest = floor(tightLower[j]+0.5);
396                //assert (fabs(tightLower[j]-nearest)<1.0e-5); may be infeasible
397                lbs.insert(j,nearest);
398                numberChanged++;
399                if (nearest>solution[j]+primalTolerance)
400                  ifCut=true;
401              }
402            } else {
403              if (upper[j]>lower[j]) {
404                if (tightUpper[j]==tightLower[j]) {
405                  // fix
406                  lbs.insert(j,tightLower[j]);
407                  ubs.insert(j,tightUpper[j]);
408                  if (tightLower[j]>solution[j]+primalTolerance||
409                      tightUpper[j]<solution[j]-primalTolerance)
410                    ifCut=true;
411                } else if (tightenBounds&&tightenBounds[j]) {
412                  lbs.insert(j,CoinMax(tightLower[j],lower[j]));
413                  ubs.insert(j,CoinMin(tightUpper[j],upper[j]));
414                  if (tightLower[j]>solution[j]+primalTolerance||
415                      tightUpper[j]<solution[j]-primalTolerance)
416                    ifCut=true;
417                }
418              }
419            }
420          }
421          if (numberChanged) {
422            OsiColCut cc;
423            cc.setUbs(ubs);
424            cc.setLbs(lbs);
425            if (ifCut) {
426              cc.setEffectiveness(100.0);
427            } else {
428              cc.setEffectiveness(1.0e-5);
429            }
430            cs.insert(cc);
431          }
432        }
433      }
434      //if (!solver->basisIsAvailable())
435      //returnCode=true;
436#if 0
437      // Pass across info to pseudocosts
438      char * mark = new char[numberColumns];
439      memset(mark,0,numberColumns);
440      int nLook = generator->numberThisTime();
441      const int * lookedAt = generator->lookedAt();
442      const int * fixedDown = generator->fixedDown();
443      const int * fixedUp = generator->fixedUp();
444      for (j=0;j<nLook;j++)
445        mark[lookedAt[j]]=1;
446      int numberObjects = model_->numberObjects();
447      for (int i=0;i<numberObjects;i++) {
448        CbcSimpleIntegerDynamicPseudoCost * obj1 =
449          dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(model_->modifiableObject(i)) ;
450        if (obj1) {
451          int iColumn = obj1->columnNumber();
452          if (mark[iColumn])
453            obj1->setProbingInformation(fixedDown[iColumn],fixedUp[iColumn]);
454        }
455      }
456      delete [] mark;
457#endif
458    }
459    CbcCutModifier * modifier = model_->cutModifier();
460    if (modifier) {
461      int numberRowCutsAfter = cs.sizeRowCuts() ;
462      int k ;
463      int nOdd=0;
464      //const OsiSolverInterface * solver = model_->solver();
465      for (k = numberRowCutsAfter-1;k>=numberRowCutsBefore;k--) {
466        OsiRowCut & thisCut = cs.rowCut(k) ;
467        int returnCode = modifier->modify(solver,thisCut);
468        if (returnCode) {
469          nOdd++;
470          if (returnCode==3)
471            cs.eraseRowCut(k);
472        }
473      }
474      if (nOdd) 
475        printf("Cut generator %s produced %d cuts of which %d were modified\n",
476                 generatorName_,numberRowCutsAfter-numberRowCutsBefore,nOdd);
477    }
478    { 
479      // make all row cuts without test for duplicate
480      int numberRowCutsAfter = cs.sizeRowCuts() ;
481      int k ;
482      for (k = numberRowCutsBefore;k<numberRowCutsAfter;k++) {
483        OsiRowCut * thisCut = cs.rowCutPtr(k) ;
484        thisCut->mutableRow().setTestForDuplicateIndex(false);
485      }
486    }
487    {
488      int numberRowCutsAfter = cs.sizeRowCuts() ;
489      if (numberRowCutsBefore<numberRowCutsAfter) {
490#if 0
491        printf("generator %s generated %d row cuts\n",
492               generatorName_,numberRowCutsAfter-numberRowCutsBefore);
493#endif
494        numberCuts_ += numberRowCutsAfter-numberRowCutsBefore;
495      }
496      int numberColumnCutsAfter = cs.sizeColCuts() ;
497      if (numberColumnCutsBefore<numberColumnCutsAfter) {
498#if 0
499        printf("generator %s generated %d column cuts\n",
500               generatorName_,numberColumnCutsAfter-numberColumnCutsBefore);
501#endif
502        numberColumnCuts_ += numberColumnCutsAfter-numberColumnCutsBefore;
503      }
504    }
505    {
506      int numberRowCutsAfter = cs.sizeRowCuts() ;
507      int k ;
508      int nEls=0;
509      int nCuts= numberRowCutsAfter-numberRowCutsBefore;
510      // Remove NULL cuts!
511      int nNull=0;
512      const double * solution = solver->getColSolution();
513      bool feasible=true;
514      for (k = numberRowCutsAfter-1;k>=numberRowCutsBefore;k--) {
515        const OsiRowCut * thisCut = cs.rowCutPtr(k) ;
516        double sum=0.0;
517        if (thisCut->lb()<=thisCut->ub()) {
518          int n=thisCut->row().getNumElements();
519          const int * column = thisCut->row().getIndices();
520          const double * element = thisCut->row().getElements();
521          if (n<=0) {
522            // infeasible cut - give up
523            feasible=false;
524            break;
525          }
526          nEls+= n;
527          for (int i=0;i<n;i++) {
528            double value = element[i];
529            sum += value*solution[column[i]];
530          }
531          if (sum>thisCut->ub()) {
532            sum= sum-thisCut->ub();
533          } else if (sum<thisCut->lb()) {
534            sum= thisCut->lb()-sum;
535          } else {
536            sum=0.0;
537            cs.eraseRowCut(k);
538            nNull++;
539          }
540        }
541      }
542      //if (nNull)
543      //printf("%s has %d cuts and %d elements - %d null!\n",generatorName_,
544      //       nCuts,nEls,nNull);
545      numberRowCutsAfter = cs.sizeRowCuts() ;
546      nCuts= numberRowCutsAfter-numberRowCutsBefore;
547      nEls=0;
548      for (k = numberRowCutsBefore;k<numberRowCutsAfter;k++) {
549        const OsiRowCut * thisCut = cs.rowCutPtr(k) ;
550        int n=thisCut->row().getNumElements();
551        nEls+= n;
552      }
553      //printf("%s has %d cuts and %d elements\n",generatorName_,
554      //     nCuts,nEls);
555      int nElsNow = solver->getMatrixByCol()->getNumElements();
556      int nAdd = model_->parentModel() ? 200 : 10000;
557      int numberColumns = solver->getNumCols();
558      int nAdd2 = model_->parentModel() ? 2*numberColumns : 5*numberColumns;
559      if (/*nEls>CoinMax(nAdd2,nElsNow/8+nAdd)*/nCuts&&feasible) {
560        //printf("need to remove cuts\n");
561        // just add most effective
562        int nReasonable = CoinMax(nAdd2,nElsNow/8+nAdd);
563        int nDelete = nEls - nReasonable;
564       
565        nElsNow = nEls;
566        double * sort = new double [nCuts];
567        int * which = new int [nCuts];
568        // For parallel cuts
569        double * element2 = new double [numberColumns];
570        CoinZeroN(element2,numberColumns);
571        for (k = numberRowCutsBefore;k<numberRowCutsAfter;k++) {
572          const OsiRowCut * thisCut = cs.rowCutPtr(k) ;
573          double sum=0.0;
574          if (thisCut->lb()<=thisCut->ub()) {
575            int n=thisCut->row().getNumElements();
576            const int * column = thisCut->row().getIndices();
577            const double * element = thisCut->row().getElements();
578            assert (n);
579            double norm=0.0;
580            for (int i=0;i<n;i++) {
581              double value = element[i];
582              sum += value*solution[column[i]];
583              norm += value*value;
584            }
585            if (sum>thisCut->ub()) {
586              sum= sum-thisCut->ub();
587            } else if (sum<thisCut->lb()) {
588              sum= thisCut->lb()-sum;
589            } else {
590              sum=0.0;
591            }
592            // normalize
593            sum /= sqrt(norm);
594            // adjust for length
595            //sum /= sqrt((double) n);
596            // randomize
597            //double randomNumber =
598            //model_->randomNumberGenerator()->randomDouble();
599            //sum *= (0.5+randomNumber);
600          } else {
601            // keep
602            sum=COIN_DBL_MAX;
603          }
604          sort[k-numberRowCutsBefore]=sum;
605          which[k-numberRowCutsBefore]=k;
606        }
607        CoinSort_2(sort,sort+nCuts,which);
608        // Now see which ones are too similar
609        int nParallel=0;
610        for (k = 0;k<nCuts;k++) {
611          int j=which[k];
612          const OsiRowCut * thisCut = cs.rowCutPtr(j) ;
613          if (thisCut->lb()>thisCut->ub()) 
614            break; // cut is infeasible
615          int n=thisCut->row().getNumElements();
616          const int * column = thisCut->row().getIndices();
617          const double * element = thisCut->row().getElements();
618          assert (n);
619          double norm=0.0;
620          double lb = thisCut->lb();
621          double ub = thisCut->ub();
622          for (int i=0;i<n;i++) {
623            double value = element[i];
624            element2[column[i]]=value;
625            norm += value*value;
626          }
627          int kkk = CoinMin(nCuts,k+5);
628          double testValue = (depth>1) ? 0.9 : 0.99999;
629          for (int kk=k+1;kk<kkk;kk++) { 
630            int jj=which[kk];
631            const OsiRowCut * thisCut2 = cs.rowCutPtr(jj) ;
632            if (thisCut2->lb()>thisCut2->ub()) 
633              break; // cut is infeasible
634            int nB=thisCut2->row().getNumElements();
635            const int * columnB = thisCut2->row().getIndices();
636            const double * elementB = thisCut2->row().getElements();
637            assert (nB);
638            double normB=0.0;
639            double product=0.0;
640            for (int i=0;i<nB;i++) {
641              double value = elementB[i];
642              normB += value*value;
643              product += value*element2[columnB[i]];
644            }
645            if (product>0.0&&product*product>testValue*norm*normB) {
646              bool parallel=true;
647              double lbB = thisCut2->lb();
648              double ubB = thisCut2->ub();
649              if ((lb<-1.0e20&&lbB>-1.0e20)||
650                  (lbB<-1.0e20&&lb>-1.0e20))
651                parallel = false;
652              double tolerance;
653              tolerance = CoinMax(fabs(lb),fabs(lbB))+1.0e-6;
654              if (fabs(lb-lbB)>tolerance)
655                parallel=false;
656              if ((ub>1.0e20&&ubB<1.0e20)||
657                  (ubB>1.0e20&&ub<1.0e20))
658                parallel = false;
659              tolerance = CoinMax(fabs(ub),fabs(ubB))+1.0e-6;
660              if (fabs(ub-ubB)>tolerance)
661                parallel=false;
662              if (parallel) {
663                nParallel++;
664                sort[k]=0.0;
665                break;
666              }
667            }
668          }
669          for (int i=0;i<n;i++) {
670            element2[column[i]]=0.0;
671          }
672        }
673        delete [] element2;
674        CoinSort_2(sort,sort+nCuts,which);
675        k=0;
676        while (nDelete>0||!sort[k]) {
677          int iCut=which[k];
678          const OsiRowCut * thisCut = cs.rowCutPtr(iCut) ;
679          int n=thisCut->row().getNumElements();
680          nDelete-=n; 
681          k++;
682          if (k>=nCuts)
683            break;
684        }
685        std::sort(which,which+k);
686        k--;
687        for (;k>=0;k--) {
688          cs.eraseRowCut(which[k]);
689        }
690        delete [] sort;
691        delete [] which;
692        numberRowCutsAfter = cs.sizeRowCuts() ;
693#if 0 //def CLP_INVESTIGATE
694        nEls=0;
695        int nCuts2= numberRowCutsAfter-numberRowCutsBefore;
696        for (k = numberRowCutsBefore;k<numberRowCutsAfter;k++) {
697          const OsiRowCut * thisCut = cs.rowCutPtr(k) ;
698          int n=thisCut->row().getNumElements();
699          nEls+= n;
700        }
701        if (!model_->parentModel()&&nCuts!=nCuts2)
702          printf("%s NOW has %d cuts and %d elements( down from %d cuts and %d els) - %d parallel\n",
703                 generatorName_,
704                 nCuts2,nEls,nCuts,nElsNow,nParallel);
705#endif
706      }
707    }
708#ifdef CBC_DEBUG
709    {
710      int numberRowCutsAfter = cs.sizeRowCuts() ;
711      int k ;
712      int nBad=0;
713      for (k = numberRowCutsBefore;k<numberRowCutsAfter;k++) {
714        OsiRowCut thisCut = cs.rowCut(k) ;
715        if (thisCut.lb()>thisCut.ub()||
716            thisCut.lb()>1.0e8||
717            thisCut.ub()<-1.0e8)
718          printf("cut from %s has bounds %g and %g!\n",
719                 generatorName_,thisCut.lb(),thisCut.ub());
720        if (thisCut.lb()<=thisCut.ub()) {
721          /* check size of elements.
722             We can allow smaller but this helps debug generators as it
723             is unsafe to have small elements */
724          int n=thisCut.row().getNumElements();
725          const int * column = thisCut.row().getIndices();
726          const double * element = thisCut.row().getElements();
727          assert (n);
728          for (int i=0;i<n;i++) {
729            double value = element[i];
730            if (fabs(value)<=1.0e-12||fabs(value)>=1.0e20)
731              nBad++;
732          }
733        }
734        if (nBad) 
735          printf("Cut generator %s produced %d cuts of which %d had tiny or large elements\n",
736                 generatorName_,numberRowCutsAfter-numberRowCutsBefore,nBad);
737      }
738    }
739#endif
740    if (timing_)
741      timeInCutGenerator_ += CoinCpuTime()-time1;
742#if 0
743    // switch off if first time and no good
744    if (node==NULL&&!pass) {
745      if (cs.sizeCuts()-cutsBefore<CoinAbs(switchOffIfLessThan_)) {
746        whenCutGenerator_=-99;
747        whenCutGeneratorInSub_ = -200;
748      }
749    }
750#endif
751  }
752  return returnCode;
753}
754void 
755CbcCutGenerator::setHowOften(int howOften) 
756{
757 
758  if (howOften>=1000000) {
759    // leave Probing every SCANCUTS_PROBING
760    howOften = howOften % 1000000;
761    CglProbing* generator =
762      dynamic_cast<CglProbing*>(generator_);
763   
764    if (generator&&howOften>SCANCUTS_PROBING) 
765      howOften=SCANCUTS_PROBING+1000000;
766    else
767      howOften += 1000000;
768  }
769  whenCutGenerator_ = howOften;
770}
771void 
772CbcCutGenerator::setWhatDepth(int value) 
773{
774  depthCutGenerator_ = value;
775}
776void 
777CbcCutGenerator::setWhatDepthInSub(int value) 
778{
779  depthCutGeneratorInSub_ = value;
780}
781
782
783// Default Constructor
784CbcCutModifier::CbcCutModifier() 
785{
786}
787
788
789// Destructor
790CbcCutModifier::~CbcCutModifier ()
791{
792}
793
794// Copy constructor
795CbcCutModifier::CbcCutModifier ( const CbcCutModifier & rhs)
796{
797}
798
799// Assignment operator
800CbcCutModifier & 
801CbcCutModifier::operator=( const CbcCutModifier& rhs)
802{
803  if (this!=&rhs) {
804  }
805  return *this;
806}
807
808// Default Constructor
809CbcCutSubsetModifier::CbcCutSubsetModifier ()
810  : CbcCutModifier(),
811    firstOdd_(COIN_INT_MAX)
812{
813}
814
815// Useful constructor
816CbcCutSubsetModifier::CbcCutSubsetModifier (int firstOdd)
817  : CbcCutModifier()
818{
819  firstOdd_=firstOdd;
820}
821
822// Copy constructor
823CbcCutSubsetModifier::CbcCutSubsetModifier ( const CbcCutSubsetModifier & rhs)
824  :CbcCutModifier(rhs)
825{
826  firstOdd_ = rhs.firstOdd_;
827}
828
829// Clone
830CbcCutModifier *
831CbcCutSubsetModifier::clone() const
832{
833  return new CbcCutSubsetModifier(*this);
834}
835
836// Assignment operator
837CbcCutSubsetModifier & 
838CbcCutSubsetModifier::operator=( const CbcCutSubsetModifier& rhs)
839{
840  if (this!=&rhs) {
841    CbcCutModifier::operator=(rhs);
842    firstOdd_ = rhs.firstOdd_;
843  }
844  return *this;
845}
846
847// Destructor
848CbcCutSubsetModifier::~CbcCutSubsetModifier ()
849{
850}
851/* Returns
852   0 unchanged
853   1 strengthened
854   2 weakened
855   3 deleted
856*/
857int 
858CbcCutSubsetModifier::modify(const OsiSolverInterface * solver, OsiRowCut & cut) 
859{
860  int n=cut.row().getNumElements();
861  if (!n)
862    return 0;
863  const int * column = cut.row().getIndices();
864  //const double * element = cut.row().getElements();
865  int returnCode=0;
866  for (int i=0;i<n;i++) {
867    if (column[i]>=firstOdd_) {
868      returnCode=3;
869      break;
870    }
871  }
872  if (!returnCode) {
873    const double * element = cut.row().getElements();
874    printf("%g <= ",cut.lb());
875    for (int i=0;i<n;i++) {
876      printf("%g*x%d ",element[i],column[i]);
877    }
878    printf("<= %g\n",cut.ub());
879  }
880  //return 3;
881  return returnCode;
882}
883
Note: See TracBrowser for help on using the repository browser.