source: stable/2.3/Cbc/src/CbcCutGenerator.cpp @ 1310

Last change on this file since 1310 was 1310, checked in by forrest, 10 years ago

check for infeasible bound tightening

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