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

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

trying to make go faster

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 24.7 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 1
171#define PROBE2 0
172#define PROBE3 1
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&&model_->parentModel()) {
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()&&info.formulation_rows>500) {
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#undef CLP_INVESTIGATE
275#ifdef CLP_INVESTIGATE
276          int kr1=0;
277          int kc1=0;
278#endif
279          int bestStackTree=-1;
280          int bestNumberTree=-1;
281          for (int i=0;i<n;i++) {
282            OsiCuts cs2 = cs;
283            int stack = test[i]/100000;
284            int number = test[i] - 100000*stack;
285            generator->setMaxLook(stack);
286            generator->setMaxProbe(number);
287            generator_->generateCuts(*solver,cs2,info);
288#ifdef CLP_INVESTIGATE
289            int numberRowCuts = cs2.sizeRowCuts()-numberRowCutsBefore ;
290            int numberColumnCuts= cs2.sizeColCuts()-numberColumnCutsBefore ;
291            if (numberRowCuts<kr1||numberColumnCuts<kc1)
292              printf("Odd ");
293            if (numberRowCuts>kr1||numberColumnCuts>kc1) {
294              printf("*** ");
295              kr1=numberRowCuts;
296              kc1=numberColumnCuts;
297              bestStackTree=stack;
298              bestNumberTree=number;
299            }
300            printf("maxStack %d number %d gives %d row cuts and %d column cuts\n",
301                   stack,number,numberRowCuts,numberColumnCuts);
302#endif
303          }
304          generator->setMaxLook(saveStack);
305          generator->setMaxProbe(saveNumber);
306          if (bestStackTree>0) {
307            generator->setMaxLook(bestStackTree);
308            generator->setMaxProbe(bestNumberTree);
309#ifdef CLP_INVESTIGATE
310            printf("RRNumber %d -> %d, stack %d -> %d\n",
311                   saveNumber,bestNumberTree,saveStack,bestStackTree);
312#endif
313          } else {
314            // no good
315            generator->setMaxLook(0);
316#ifdef CLP_INVESTIGATE
317            printf("RRSwitching off number %d -> %d, stack %d -> %d\n",
318                   saveNumber,saveNumber,saveStack,0);
319#endif
320          }
321        }
322#endif
323        if (generator->getMaxLook()>0) {
324          generator->generateCutsAndModify(*solver,cs,&info);
325          doCuts=true;
326        }
327      } else {
328        // at root - don't always do
329        if (!PROBE3||pass<15||(pass&1)==0) {
330          generator->generateCutsAndModify(*solver,cs,&info);
331          doCuts=true;
332        }
333      }
334      if (doCuts) {
335        const double * tightLower = generator->tightLower();
336        const double * lower = solver->getColLower();
337        const double * tightUpper = generator->tightUpper();
338        const double * upper = solver->getColUpper();
339        const double * solution = solver->getColSolution();
340        int j;
341        int numberColumns = solver->getNumCols();
342        double primalTolerance = 1.0e-8;
343        const char * tightenBounds = generator->tightenBounds();
344        if ((model_->getThreadMode()&2)==0) {
345          for (j=0;j<numberColumns;j++) {
346            if (solver->isInteger(j)) {
347              if (tightUpper[j]<upper[j]) {
348                double nearest = floor(tightUpper[j]+0.5);
349                //assert (fabs(tightUpper[j]-nearest)<1.0e-5); may be infeasible
350                solver->setColUpper(j,nearest);
351                if (nearest<solution[j]-primalTolerance)
352                  returnCode=true;
353              }
354              if (tightLower[j]>lower[j]) {
355                double nearest = floor(tightLower[j]+0.5);
356                //assert (fabs(tightLower[j]-nearest)<1.0e-5); may be infeasible
357                solver->setColLower(j,nearest);
358                if (nearest>solution[j]+primalTolerance)
359                  returnCode=true;
360              }
361            } else {
362              if (upper[j]>lower[j]) {
363                if (tightUpper[j]==tightLower[j]) {
364                  // fix
365                  solver->setColLower(j,tightLower[j]);
366                  solver->setColUpper(j,tightUpper[j]);
367                  if (tightLower[j]>solution[j]+primalTolerance||
368                      tightUpper[j]<solution[j]-primalTolerance)
369                    returnCode=true;
370                } else if (tightenBounds&&tightenBounds[j]) {
371                  solver->setColLower(j,CoinMax(tightLower[j],lower[j]));
372                  solver->setColUpper(j,CoinMin(tightUpper[j],upper[j]));
373                  if (tightLower[j]>solution[j]+primalTolerance||
374                      tightUpper[j]<solution[j]-primalTolerance)
375                    returnCode=true;
376                }
377              }
378            }
379          }
380        } else {
381          CoinPackedVector lbs;
382          CoinPackedVector ubs;
383          int numberChanged=0;
384          bool ifCut=false;
385          for (j=0;j<numberColumns;j++) {
386            if (solver->isInteger(j)) {
387              if (tightUpper[j]<upper[j]) {
388                double nearest = floor(tightUpper[j]+0.5);
389                //assert (fabs(tightUpper[j]-nearest)<1.0e-5); may be infeasible
390                ubs.insert(j,nearest);
391                numberChanged++;
392                if (nearest<solution[j]-primalTolerance)
393                  ifCut=true;
394              }
395              if (tightLower[j]>lower[j]) {
396                double nearest = floor(tightLower[j]+0.5);
397                //assert (fabs(tightLower[j]-nearest)<1.0e-5); may be infeasible
398                lbs.insert(j,nearest);
399                numberChanged++;
400                if (nearest>solution[j]+primalTolerance)
401                  ifCut=true;
402              }
403            } else {
404              if (upper[j]>lower[j]) {
405                if (tightUpper[j]==tightLower[j]) {
406                  // fix
407                  lbs.insert(j,tightLower[j]);
408                  ubs.insert(j,tightUpper[j]);
409                  if (tightLower[j]>solution[j]+primalTolerance||
410                      tightUpper[j]<solution[j]-primalTolerance)
411                    ifCut=true;
412                } else if (tightenBounds&&tightenBounds[j]) {
413                  lbs.insert(j,CoinMax(tightLower[j],lower[j]));
414                  ubs.insert(j,CoinMin(tightUpper[j],upper[j]));
415                  if (tightLower[j]>solution[j]+primalTolerance||
416                      tightUpper[j]<solution[j]-primalTolerance)
417                    ifCut=true;
418                }
419              }
420            }
421          }
422          if (numberChanged) {
423            OsiColCut cc;
424            cc.setUbs(ubs);
425            cc.setLbs(lbs);
426            if (ifCut) {
427              cc.setEffectiveness(100.0);
428            } else {
429              cc.setEffectiveness(1.0e-5);
430            }
431            cs.insert(cc);
432          }
433        }
434      }
435      //if (!solver->basisIsAvailable())
436      //returnCode=true;
437#if 0
438      // Pass across info to pseudocosts
439      char * mark = new char[numberColumns];
440      memset(mark,0,numberColumns);
441      int nLook = generator->numberThisTime();
442      const int * lookedAt = generator->lookedAt();
443      const int * fixedDown = generator->fixedDown();
444      const int * fixedUp = generator->fixedUp();
445      for (j=0;j<nLook;j++)
446        mark[lookedAt[j]]=1;
447      int numberObjects = model_->numberObjects();
448      for (int i=0;i<numberObjects;i++) {
449        CbcSimpleIntegerDynamicPseudoCost * obj1 =
450          dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(model_->modifiableObject(i)) ;
451        if (obj1) {
452          int iColumn = obj1->columnNumber();
453          if (mark[iColumn])
454            obj1->setProbingInformation(fixedDown[iColumn],fixedUp[iColumn]);
455        }
456      }
457      delete [] mark;
458#endif
459    }
460    CbcCutModifier * modifier = model_->cutModifier();
461    if (modifier) {
462      int numberRowCutsAfter = cs.sizeRowCuts() ;
463      int k ;
464      int nOdd=0;
465      //const OsiSolverInterface * solver = model_->solver();
466      for (k = numberRowCutsAfter-1;k>=numberRowCutsBefore;k--) {
467        OsiRowCut & thisCut = cs.rowCut(k) ;
468        int returnCode = modifier->modify(solver,thisCut);
469        if (returnCode) {
470          nOdd++;
471          if (returnCode==3)
472            cs.eraseRowCut(k);
473        }
474      }
475      if (nOdd) 
476        printf("Cut generator %s produced %d cuts of which %d were modified\n",
477                 generatorName_,numberRowCutsAfter-numberRowCutsBefore,nOdd);
478    }
479    { 
480      // make all row cuts without test for duplicate
481      int numberRowCutsAfter = cs.sizeRowCuts() ;
482      int k ;
483      for (k = numberRowCutsBefore;k<numberRowCutsAfter;k++) {
484        OsiRowCut * thisCut = cs.rowCutPtr(k) ;
485        thisCut->mutableRow().setTestForDuplicateIndex(false);
486      }
487    }
488    {
489      int numberRowCutsAfter = cs.sizeRowCuts() ;
490      if (numberRowCutsBefore<numberRowCutsAfter) {
491#if 0
492        printf("generator %s generated %d row cuts\n",
493               generatorName_,numberRowCutsAfter-numberRowCutsBefore);
494#endif
495        numberCuts_ += numberRowCutsAfter-numberRowCutsBefore;
496      }
497      int numberColumnCutsAfter = cs.sizeColCuts() ;
498      if (numberColumnCutsBefore<numberColumnCutsAfter) {
499#if 0
500        printf("generator %s generated %d column cuts\n",
501               generatorName_,numberColumnCutsAfter-numberColumnCutsBefore);
502#endif
503        numberColumnCuts_ += numberColumnCutsAfter-numberColumnCutsBefore;
504      }
505    }
506    {
507      int numberRowCutsAfter = cs.sizeRowCuts() ;
508      int k ;
509      int nEls=0;
510      int nCuts= numberRowCutsAfter-numberRowCutsBefore;
511      // Remove NULL cuts!
512      int nNull=0;
513      const double * solution = solver->getColSolution();
514      bool feasible=true;
515      double primalTolerance = 1.0e-7;
516      for (k = numberRowCutsAfter-1;k>=numberRowCutsBefore;k--) {
517        const OsiRowCut * thisCut = cs.rowCutPtr(k) ;
518        double sum=0.0;
519        if (thisCut->lb()<=thisCut->ub()) {
520          int n=thisCut->row().getNumElements();
521          const int * column = thisCut->row().getIndices();
522          const double * element = thisCut->row().getElements();
523          if (n<=0) {
524            // infeasible cut - give up
525            feasible=false;
526            break;
527          }
528          nEls+= n;
529          for (int i=0;i<n;i++) {
530            double value = element[i];
531            sum += value*solution[column[i]];
532          }
533          if (sum>thisCut->ub()+primalTolerance) {
534            sum= sum-thisCut->ub();
535          } else if (sum<thisCut->lb()-primalTolerance) {
536            sum= thisCut->lb()-sum;
537          } else {
538            sum=0.0;
539            cs.eraseRowCut(k);
540            nNull++;
541          }
542        }
543      }
544      //if (nNull)
545      //printf("%s has %d cuts and %d elements - %d null!\n",generatorName_,
546      //       nCuts,nEls,nNull);
547      numberRowCutsAfter = cs.sizeRowCuts() ;
548      nCuts= numberRowCutsAfter-numberRowCutsBefore;
549      nEls=0;
550      for (k = numberRowCutsBefore;k<numberRowCutsAfter;k++) {
551        const OsiRowCut * thisCut = cs.rowCutPtr(k) ;
552        int n=thisCut->row().getNumElements();
553        nEls+= n;
554      }
555      //printf("%s has %d cuts and %d elements\n",generatorName_,
556      //     nCuts,nEls);
557      int nElsNow = solver->getMatrixByCol()->getNumElements();
558      int nAdd = model_->parentModel() ? 200 : 10000;
559      int numberColumns = solver->getNumCols();
560      int nAdd2 = model_->parentModel() ? 2*numberColumns : 5*numberColumns;
561      if (/*nEls>CoinMax(nAdd2,nElsNow/8+nAdd)*/nCuts&&feasible) {
562        //printf("need to remove cuts\n");
563        // just add most effective
564        int nReasonable = CoinMax(nAdd2,nElsNow/8+nAdd);
565        int nDelete = nEls - nReasonable;
566       
567        nElsNow = nEls;
568        double * sort = new double [nCuts];
569        int * which = new int [nCuts];
570        // For parallel cuts
571        double * element2 = new double [numberColumns];
572        CoinZeroN(element2,numberColumns);
573        for (k = numberRowCutsBefore;k<numberRowCutsAfter;k++) {
574          const OsiRowCut * thisCut = cs.rowCutPtr(k) ;
575          double sum=0.0;
576          if (thisCut->lb()<=thisCut->ub()) {
577            int n=thisCut->row().getNumElements();
578            const int * column = thisCut->row().getIndices();
579            const double * element = thisCut->row().getElements();
580            assert (n);
581            double norm=0.0;
582            for (int i=0;i<n;i++) {
583              double value = element[i];
584              sum += value*solution[column[i]];
585              norm += value*value;
586            }
587            if (sum>thisCut->ub()) {
588              sum= sum-thisCut->ub();
589            } else if (sum<thisCut->lb()) {
590              sum= thisCut->lb()-sum;
591            } else {
592              sum=0.0;
593            }
594            // normalize
595            sum /= sqrt(norm);
596            // adjust for length
597            //sum /= sqrt((double) n);
598            // randomize
599            //double randomNumber =
600            //model_->randomNumberGenerator()->randomDouble();
601            //sum *= (0.5+randomNumber);
602          } else {
603            // keep
604            sum=COIN_DBL_MAX;
605          }
606          sort[k-numberRowCutsBefore]=sum;
607          which[k-numberRowCutsBefore]=k;
608        }
609        CoinSort_2(sort,sort+nCuts,which);
610        // Now see which ones are too similar
611        int nParallel=0;
612        double testValue = (depth>1) ? 0.99 : 0.999999;
613        for (k = 0;k<nCuts;k++) {
614          int j=which[k];
615          const OsiRowCut * thisCut = cs.rowCutPtr(j) ;
616          if (thisCut->lb()>thisCut->ub()) 
617            break; // cut is infeasible
618          int n=thisCut->row().getNumElements();
619          const int * column = thisCut->row().getIndices();
620          const double * element = thisCut->row().getElements();
621          assert (n);
622          double norm=0.0;
623          double lb = thisCut->lb();
624          double ub = thisCut->ub();
625          for (int i=0;i<n;i++) {
626            double value = element[i];
627            element2[column[i]]=value;
628            norm += value*value;
629          }
630          int kkk = CoinMin(nCuts,k+5);
631          for (int kk=k+1;kk<kkk;kk++) { 
632            int jj=which[kk];
633            const OsiRowCut * thisCut2 = cs.rowCutPtr(jj) ;
634            if (thisCut2->lb()>thisCut2->ub()) 
635              break; // cut is infeasible
636            int nB=thisCut2->row().getNumElements();
637            const int * columnB = thisCut2->row().getIndices();
638            const double * elementB = thisCut2->row().getElements();
639            assert (nB);
640            double normB=0.0;
641            double product=0.0;
642            for (int i=0;i<nB;i++) {
643              double value = elementB[i];
644              normB += value*value;
645              product += value*element2[columnB[i]];
646            }
647            if (product>0.0&&product*product>testValue*norm*normB) {
648              bool parallel=true;
649              double lbB = thisCut2->lb();
650              double ubB = thisCut2->ub();
651              if ((lb<-1.0e20&&lbB>-1.0e20)||
652                  (lbB<-1.0e20&&lb>-1.0e20))
653                parallel = false;
654              double tolerance;
655              tolerance = CoinMax(fabs(lb),fabs(lbB))+1.0e-6;
656              if (fabs(lb-lbB)>tolerance)
657                parallel=false;
658              if ((ub>1.0e20&&ubB<1.0e20)||
659                  (ubB>1.0e20&&ub<1.0e20))
660                parallel = false;
661              tolerance = CoinMax(fabs(ub),fabs(ubB))+1.0e-6;
662              if (fabs(ub-ubB)>tolerance)
663                parallel=false;
664              if (parallel) {
665                nParallel++;
666                sort[k]=0.0;
667                break;
668              }
669            }
670          }
671          for (int i=0;i<n;i++) {
672            element2[column[i]]=0.0;
673          }
674        }
675        delete [] element2;
676        CoinSort_2(sort,sort+nCuts,which);
677        k=0;
678        while (nDelete>0||!sort[k]) {
679          int iCut=which[k];
680          const OsiRowCut * thisCut = cs.rowCutPtr(iCut) ;
681          int n=thisCut->row().getNumElements();
682          nDelete-=n; 
683          k++;
684          if (k>=nCuts)
685            break;
686        }
687        std::sort(which,which+k);
688        k--;
689        for (;k>=0;k--) {
690          cs.eraseRowCut(which[k]);
691        }
692        delete [] sort;
693        delete [] which;
694        numberRowCutsAfter = cs.sizeRowCuts() ;
695#if 0 //def CLP_INVESTIGATE
696        nEls=0;
697        int nCuts2= numberRowCutsAfter-numberRowCutsBefore;
698        for (k = numberRowCutsBefore;k<numberRowCutsAfter;k++) {
699          const OsiRowCut * thisCut = cs.rowCutPtr(k) ;
700          int n=thisCut->row().getNumElements();
701          nEls+= n;
702        }
703        if (!model_->parentModel()&&nCuts!=nCuts2)
704          printf("%s NOW has %d cuts and %d elements( down from %d cuts and %d els) - %d parallel\n",
705                 generatorName_,
706                 nCuts2,nEls,nCuts,nElsNow,nParallel);
707#endif
708      }
709    }
710#ifdef CBC_DEBUG
711    {
712      int numberRowCutsAfter = cs.sizeRowCuts() ;
713      int k ;
714      int nBad=0;
715      for (k = numberRowCutsBefore;k<numberRowCutsAfter;k++) {
716        OsiRowCut thisCut = cs.rowCut(k) ;
717        if (thisCut.lb()>thisCut.ub()||
718            thisCut.lb()>1.0e8||
719            thisCut.ub()<-1.0e8)
720          printf("cut from %s has bounds %g and %g!\n",
721                 generatorName_,thisCut.lb(),thisCut.ub());
722        if (thisCut.lb()<=thisCut.ub()) {
723          /* check size of elements.
724             We can allow smaller but this helps debug generators as it
725             is unsafe to have small elements */
726          int n=thisCut.row().getNumElements();
727          const int * column = thisCut.row().getIndices();
728          const double * element = thisCut.row().getElements();
729          assert (n);
730          for (int i=0;i<n;i++) {
731            double value = element[i];
732            if (fabs(value)<=1.0e-12||fabs(value)>=1.0e20)
733              nBad++;
734          }
735        }
736        if (nBad) 
737          printf("Cut generator %s produced %d cuts of which %d had tiny or large elements\n",
738                 generatorName_,numberRowCutsAfter-numberRowCutsBefore,nBad);
739      }
740    }
741#endif
742    if (timing_)
743      timeInCutGenerator_ += CoinCpuTime()-time1;
744#if 0
745    // switch off if first time and no good
746    if (node==NULL&&!pass) {
747      if (cs.sizeCuts()-cutsBefore<CoinAbs(switchOffIfLessThan_)) {
748        whenCutGenerator_=-99;
749        whenCutGeneratorInSub_ = -200;
750      }
751    }
752#endif
753  }
754  return returnCode;
755}
756void 
757CbcCutGenerator::setHowOften(int howOften) 
758{
759 
760  if (howOften>=1000000) {
761    // leave Probing every SCANCUTS_PROBING
762    howOften = howOften % 1000000;
763    CglProbing* generator =
764      dynamic_cast<CglProbing*>(generator_);
765   
766    if (generator&&howOften>SCANCUTS_PROBING) 
767      howOften=SCANCUTS_PROBING+1000000;
768    else
769      howOften += 1000000;
770  }
771  whenCutGenerator_ = howOften;
772}
773void 
774CbcCutGenerator::setWhatDepth(int value) 
775{
776  depthCutGenerator_ = value;
777}
778void 
779CbcCutGenerator::setWhatDepthInSub(int value) 
780{
781  depthCutGeneratorInSub_ = value;
782}
783
784
785// Default Constructor
786CbcCutModifier::CbcCutModifier() 
787{
788}
789
790
791// Destructor
792CbcCutModifier::~CbcCutModifier ()
793{
794}
795
796// Copy constructor
797CbcCutModifier::CbcCutModifier ( const CbcCutModifier & rhs)
798{
799}
800
801// Assignment operator
802CbcCutModifier & 
803CbcCutModifier::operator=( const CbcCutModifier& rhs)
804{
805  if (this!=&rhs) {
806  }
807  return *this;
808}
809
810// Default Constructor
811CbcCutSubsetModifier::CbcCutSubsetModifier ()
812  : CbcCutModifier(),
813    firstOdd_(COIN_INT_MAX)
814{
815}
816
817// Useful constructor
818CbcCutSubsetModifier::CbcCutSubsetModifier (int firstOdd)
819  : CbcCutModifier()
820{
821  firstOdd_=firstOdd;
822}
823
824// Copy constructor
825CbcCutSubsetModifier::CbcCutSubsetModifier ( const CbcCutSubsetModifier & rhs)
826  :CbcCutModifier(rhs)
827{
828  firstOdd_ = rhs.firstOdd_;
829}
830
831// Clone
832CbcCutModifier *
833CbcCutSubsetModifier::clone() const
834{
835  return new CbcCutSubsetModifier(*this);
836}
837
838// Assignment operator
839CbcCutSubsetModifier & 
840CbcCutSubsetModifier::operator=( const CbcCutSubsetModifier& rhs)
841{
842  if (this!=&rhs) {
843    CbcCutModifier::operator=(rhs);
844    firstOdd_ = rhs.firstOdd_;
845  }
846  return *this;
847}
848
849// Destructor
850CbcCutSubsetModifier::~CbcCutSubsetModifier ()
851{
852}
853/* Returns
854   0 unchanged
855   1 strengthened
856   2 weakened
857   3 deleted
858*/
859int 
860CbcCutSubsetModifier::modify(const OsiSolverInterface * solver, OsiRowCut & cut) 
861{
862  int n=cut.row().getNumElements();
863  if (!n)
864    return 0;
865  const int * column = cut.row().getIndices();
866  //const double * element = cut.row().getElements();
867  int returnCode=0;
868  for (int i=0;i<n;i++) {
869    if (column[i]>=firstOdd_) {
870      returnCode=3;
871      break;
872    }
873  }
874  if (!returnCode) {
875    const double * element = cut.row().getElements();
876    printf("%g <= ",cut.lb());
877    for (int i=0;i<n;i++) {
878      printf("%g*x%d ",element[i],column[i]);
879    }
880    printf("<= %g\n",cut.ub());
881  }
882  //return 3;
883  return returnCode;
884}
885
Note: See TracBrowser for help on using the repository browser.