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

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

compiler warnings, deterministic parallel and stability

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