source: stable/2.4/Cbc/src/CbcCutGenerator.cpp @ 1271

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

Creating new stable branch 2.4 from trunk (rev 1270)

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