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

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

changes for heuristic clone

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