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

Last change on this file since 1309 was 1309, checked in by forrest, 11 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: 34.2 KB
Line 
1/* $Id: CbcCutGenerator.cpp 1309 2009-11-18 09:55:13Z 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        bool feasible=true;
380        if ((model_->getThreadMode()&2)==0) {
381          for (j=0;j<numberColumns;j++) {
382            if (solver->isInteger(j)) {
383              if (tightUpper[j]<upper[j]) {
384                double nearest = floor(tightUpper[j]+0.5);
385                //assert (fabs(tightUpper[j]-nearest)<1.0e-5); may be infeasible
386                solver->setColUpper(j,nearest);
387                if (nearest<solution[j]-primalTolerance)
388                  returnCode=true;
389              }
390              if (tightLower[j]>lower[j]) {
391                double nearest = floor(tightLower[j]+0.5);
392                //assert (fabs(tightLower[j]-nearest)<1.0e-5); may be infeasible
393                solver->setColLower(j,nearest);
394                if (nearest>solution[j]+primalTolerance)
395                  returnCode=true;
396              }
397            } else {
398              if (upper[j]>lower[j]) {
399                if (tightUpper[j]==tightLower[j]) {
400                  // fix
401                  solver->setColLower(j,tightLower[j]);
402                  solver->setColUpper(j,tightUpper[j]);
403                  if (tightLower[j]>solution[j]+primalTolerance||
404                      tightUpper[j]<solution[j]-primalTolerance)
405                    returnCode=true;
406                } else if (tightenBounds&&tightenBounds[j]) {
407                  solver->setColLower(j,CoinMax(tightLower[j],lower[j]));
408                  solver->setColUpper(j,CoinMin(tightUpper[j],upper[j]));
409                  if (tightLower[j]>solution[j]+primalTolerance||
410                      tightUpper[j]<solution[j]-primalTolerance)
411                    returnCode=true;
412                }
413              }
414            }
415            if(upper[j]<lower[j]-1.0e-3) {
416              feasible=false;
417              break;
418            }
419          }
420        } else {
421          CoinPackedVector lbs;
422          CoinPackedVector ubs;
423          int numberChanged=0;
424          bool ifCut=false;
425          for (j=0;j<numberColumns;j++) {
426            if (solver->isInteger(j)) {
427              if (tightUpper[j]<upper[j]) {
428                double nearest = floor(tightUpper[j]+0.5);
429                //assert (fabs(tightUpper[j]-nearest)<1.0e-5); may be infeasible
430                ubs.insert(j,nearest);
431                numberChanged++;
432                if (nearest<solution[j]-primalTolerance)
433                  ifCut=true;
434              }
435              if (tightLower[j]>lower[j]) {
436                double nearest = floor(tightLower[j]+0.5);
437                //assert (fabs(tightLower[j]-nearest)<1.0e-5); may be infeasible
438                lbs.insert(j,nearest);
439                numberChanged++;
440                if (nearest>solution[j]+primalTolerance)
441                  ifCut=true;
442              }
443            } else {
444              if (upper[j]>lower[j]) {
445                if (tightUpper[j]==tightLower[j]) {
446                  // fix
447                  lbs.insert(j,tightLower[j]);
448                  ubs.insert(j,tightUpper[j]);
449                  if (tightLower[j]>solution[j]+primalTolerance||
450                      tightUpper[j]<solution[j]-primalTolerance)
451                    ifCut=true;
452                } else if (tightenBounds&&tightenBounds[j]) {
453                  lbs.insert(j,CoinMax(tightLower[j],lower[j]));
454                  ubs.insert(j,CoinMin(tightUpper[j],upper[j]));
455                  if (tightLower[j]>solution[j]+primalTolerance||
456                      tightUpper[j]<solution[j]-primalTolerance)
457                    ifCut=true;
458                }
459              }
460            }
461            if(upper[j]<lower[j]-1.0e-3) {
462              feasible=false;
463              break;
464            }
465          }
466          if (numberChanged) {
467            OsiColCut cc;
468            cc.setUbs(ubs);
469            cc.setLbs(lbs);
470            if (ifCut) {
471              cc.setEffectiveness(100.0);
472            } else {
473              cc.setEffectiveness(1.0e-5);
474            }
475            cs.insert(cc);
476          }
477        }
478        if (!feasible) {
479          // not feasible -add infeasible cut
480          OsiRowCut rc;
481          rc.setLb(DBL_MAX);
482          rc.setUb(0.0);   
483          cs.insert(rc);
484        }
485      }
486      //if (!solver->basisIsAvailable())
487      //returnCode=true;
488#if 0
489      // Pass across info to pseudocosts
490      char * mark = new char[numberColumns];
491      memset(mark,0,numberColumns);
492      int nLook = generator->numberThisTime();
493      const int * lookedAt = generator->lookedAt();
494      const int * fixedDown = generator->fixedDown();
495      const int * fixedUp = generator->fixedUp();
496      for (j=0;j<nLook;j++)
497        mark[lookedAt[j]]=1;
498      int numberObjects = model_->numberObjects();
499      for (int i=0;i<numberObjects;i++) {
500        CbcSimpleIntegerDynamicPseudoCost * obj1 =
501          dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(model_->modifiableObject(i)) ;
502        if (obj1) {
503          int iColumn = obj1->columnNumber();
504          if (mark[iColumn])
505            obj1->setProbingInformation(fixedDown[iColumn],fixedUp[iColumn]);
506        }
507      }
508      delete [] mark;
509#endif
510    }
511    CbcCutModifier * modifier = model_->cutModifier();
512    if (modifier) {
513      int numberRowCutsAfter = cs.sizeRowCuts() ;
514      int k ;
515      int nOdd=0;
516      //const OsiSolverInterface * solver = model_->solver();
517      for (k = numberRowCutsAfter-1;k>=numberRowCutsBefore;k--) {
518        OsiRowCut & thisCut = cs.rowCut(k) ;
519        int returnCode = modifier->modify(solver,thisCut);
520        if (returnCode) {
521          nOdd++;
522          if (returnCode==3)
523            cs.eraseRowCut(k);
524        }
525      }
526      if (nOdd) 
527        printf("Cut generator %s produced %d cuts of which %d were modified\n",
528                 generatorName_,numberRowCutsAfter-numberRowCutsBefore,nOdd);
529    }
530    { 
531      // make all row cuts without test for duplicate
532      int numberRowCutsAfter = cs.sizeRowCuts() ;
533      int k ;
534      for (k = numberRowCutsBefore;k<numberRowCutsAfter;k++) {
535        OsiRowCut * thisCut = cs.rowCutPtr(k) ;
536        thisCut->mutableRow().setTestForDuplicateIndex(false);
537      }
538    }
539    // Add in saved cuts if violated
540    if (false&&!depth) {
541      const double * solution = solver->getColSolution();
542      double primalTolerance = 1.0e-7;
543      int numberCuts = savedCuts_.sizeRowCuts() ;
544      for (int k = numberCuts-1;k>=0;k--) {
545        const OsiRowCut * thisCut = savedCuts_.rowCutPtr(k) ;
546        double sum=0.0;
547        int n=thisCut->row().getNumElements();
548        const int * column = thisCut->row().getIndices();
549        const double * element = thisCut->row().getElements();
550        assert (n);
551        for (int i=0;i<n;i++) {
552          double value = element[i];
553          sum += value*solution[column[i]];
554        }
555        if (sum>thisCut->ub()+primalTolerance) {
556          sum= sum-thisCut->ub();
557        } else if (sum<thisCut->lb()-primalTolerance) {
558          sum= thisCut->lb()-sum;
559        } else {
560          sum=0.0;
561        }
562        if (sum) {
563          // add to candidates and take out here
564          cs.insert(*thisCut);
565          savedCuts_.eraseRowCut(k);
566        }
567      }
568    }
569    if (!atSolution()) {
570      int numberRowCutsAfter = cs.sizeRowCuts() ;
571      int k ;
572      int nEls=0;
573      int nCuts= numberRowCutsAfter-numberRowCutsBefore;
574      // Remove NULL cuts!
575      int nNull=0;
576      const double * solution = solver->getColSolution();
577      bool feasible=true;
578      double primalTolerance = 1.0e-7;
579      int shortCut = (depth) ? -1 : generator_->maximumLengthOfCutInTree();
580      for (k = numberRowCutsAfter-1;k>=numberRowCutsBefore;k--) {
581        const OsiRowCut * thisCut = cs.rowCutPtr(k) ;
582        double sum=0.0;
583        if (thisCut->lb()<=thisCut->ub()) {
584          int n=thisCut->row().getNumElements();
585          if (n<=shortCut)
586            numberShortCutsAtRoot_++;
587          const int * column = thisCut->row().getIndices();
588          const double * element = thisCut->row().getElements();
589          if (n<=0) {
590            // infeasible cut - give up
591            feasible=false;
592            break;
593          }
594          nEls+= n;
595          for (int i=0;i<n;i++) {
596            double value = element[i];
597            sum += value*solution[column[i]];
598          }
599          if (sum>thisCut->ub()+primalTolerance) {
600            sum= sum-thisCut->ub();
601          } else if (sum<thisCut->lb()-primalTolerance) {
602            sum= thisCut->lb()-sum;
603          } else {
604            sum=0.0;
605            cs.eraseRowCut(k);
606            nNull++;
607          }
608        }
609      }
610      //if (nNull)
611      //printf("%s has %d cuts and %d elements - %d null!\n",generatorName_,
612      //       nCuts,nEls,nNull);
613      numberRowCutsAfter = cs.sizeRowCuts() ;
614      nCuts= numberRowCutsAfter-numberRowCutsBefore;
615      nEls=0;
616      for (k = numberRowCutsBefore;k<numberRowCutsAfter;k++) {
617        const OsiRowCut * thisCut = cs.rowCutPtr(k) ;
618        int n=thisCut->row().getNumElements();
619        nEls+= n;
620      }
621      //printf("%s has %d cuts and %d elements\n",generatorName_,
622      //     nCuts,nEls);
623      int nElsNow = solver->getMatrixByCol()->getNumElements();
624      int numberColumns = solver->getNumCols();
625      int numberRows = solver->getNumRows();
626      //double averagePerRow = static_cast<double>(nElsNow)/
627      //static_cast<double>(numberRows);
628      int nAdd;
629      int nAdd2;
630      int nReasonable;
631      if (!model_->parentModel()&&depth<2) {
632        if (inaccuracy_<3) {
633          nAdd=10000;
634          if (pass>0&&numberColumns>-500)
635            nAdd = CoinMin(nAdd,nElsNow+2*numberRows);
636        } else {
637          nAdd=10000;
638          if (pass>0)
639            nAdd = CoinMin(nAdd,nElsNow+2*numberRows);
640        }
641        nAdd2 = 5*numberColumns;
642        nReasonable = CoinMax(nAdd2,nElsNow/8+nAdd);
643        if (!depth&&!pass) {
644          // allow more
645          nAdd += nElsNow/2;
646          nAdd2 += nElsNow/2;
647          nReasonable += nElsNow/2;
648        }
649        //if (!depth&&ineffectualCuts())
650        //nReasonable *= 2;
651      } else {
652        nAdd = 200;
653        nAdd2 = 2*numberColumns;
654        nReasonable = CoinMax(nAdd2,nElsNow/8+nAdd);
655      }
656      //#define UNS_WEIGHT 0.1
657#ifdef UNS_WEIGHT
658      const double * colLower = solver->getColLower();
659      const double * colUpper = solver->getColUpper();
660#endif
661      if (/*nEls>CoinMax(nAdd2,nElsNow/8+nAdd)*/nCuts&&feasible) {
662        //printf("need to remove cuts\n");
663        // just add most effective
664#if 1
665        int nDelete = nEls - nReasonable;
666       
667        nElsNow = nEls;
668        double * sort = new double [nCuts];
669        int * which = new int [nCuts];
670        // For parallel cuts
671        double * element2 = new double [numberColumns];
672        //#define USE_OBJECTIVE 2
673#ifdef USE_OBJECTIVE
674        const double *objective = solver->getObjCoefficients() ;
675#if USE_OBJECTIVE>1
676        double objNorm=0.0;
677        for (int i=0;i<numberColumns;i++)
678          objNorm += objective[i]*objective[i];
679        if (objNorm)
680          objNorm = 1.0/sqrt(objNorm);
681        else
682          objNorm=1.0;
683        objNorm *= 0.01; // downgrade
684#endif
685#endif
686        CoinZeroN(element2,numberColumns);
687        for (k = numberRowCutsBefore;k<numberRowCutsAfter;k++) {
688          const OsiRowCut * thisCut = cs.rowCutPtr(k) ;
689          double sum=0.0;
690          if (thisCut->lb()<=thisCut->ub()) {
691            int n=thisCut->row().getNumElements();
692            const int * column = thisCut->row().getIndices();
693            const double * element = thisCut->row().getElements();
694            assert (n);
695#ifdef UNS_WEIGHT
696            double normU=0.0;
697            double norm=1.0e-3;
698            int nU=0;
699            for (int i=0;i<n;i++) {
700              double value = element[i];
701              int iColumn = column[i];
702              double solValue = solution[iColumn];
703              sum += value*solValue;
704              value *= value;
705              norm += value;
706              if (solValue>colLower[iColumn]+1.0e-6&&
707                  solValue<colUpper[iColumn]-1.0e-6) {
708                normU += value;
709                nU++;
710              }
711            }
712#if 0
713            int nS=n-nU;
714            if (numberColumns>20000) {
715              if (nS>50) {
716                double ratio = 50.0/nS;
717                normU /= ratio;
718              }
719            }
720#endif
721            norm += UNS_WEIGHT*(normU-norm);
722#else
723            double norm=1.0e-3;
724#ifdef USE_OBJECTIVE
725            double obj=0.0;
726#endif
727            for (int i=0;i<n;i++) {
728              int iColumn = column[i];
729              double value = element[i];
730              sum += value*solution[iColumn];
731              norm += value*value;
732#ifdef USE_OBJECTIVE
733              obj += value*objective[iColumn];
734#endif
735            }
736#endif
737            if (sum>thisCut->ub()) {
738              sum= sum-thisCut->ub();
739            } else if (sum<thisCut->lb()) {
740              sum= thisCut->lb()-sum;
741            } else {
742              sum=0.0;
743            }
744#ifdef USE_OBJECTIVE
745            if (sum) {
746#if USE_OBJECTIVE==1
747              obj = CoinMax(1.0e-6,fabs(obj));
748              norm=sqrt(obj*norm);
749              //sum += fabs(obj)*invObjNorm;
750              //printf("sum %g norm %g normobj %g invNorm %g mod %g\n",
751              //     sum,norm,obj,invObjNorm,obj*invObjNorm);
752              // normalize
753              sum /= sqrt(norm);
754#else
755              // normalize
756              norm = 1.0/sqrt(norm);
757              sum = (sum+objNorm*obj)*norm;
758#endif
759            }
760#else
761            // normalize
762            sum /= sqrt(norm);
763#endif
764            //sum /= pow(norm,0.3);
765            // adjust for length
766            //sum /= pow(reinterpret_cast<double>(n),0.2);
767            //sum /= sqrt((double) n);
768            // randomize
769            //double randomNumber =
770            //model_->randomNumberGenerator()->randomDouble();
771            //sum *= (0.5+randomNumber);
772          } else {
773            // keep
774            sum=COIN_DBL_MAX;
775          }
776          sort[k-numberRowCutsBefore]=sum;
777          which[k-numberRowCutsBefore]=k;
778        }
779        CoinSort_2(sort,sort+nCuts,which);
780        // Now see which ones are too similar
781        int nParallel=0;
782        double testValue = (depth>1) ? 0.99 : 0.999999;
783        for (k = 0;k<nCuts;k++) {
784          int j=which[k];
785          const OsiRowCut * thisCut = cs.rowCutPtr(j) ;
786          if (thisCut->lb()>thisCut->ub()) 
787            break; // cut is infeasible
788          int n=thisCut->row().getNumElements();
789          const int * column = thisCut->row().getIndices();
790          const double * element = thisCut->row().getElements();
791          assert (n);
792          double norm=0.0;
793          double lb = thisCut->lb();
794          double ub = thisCut->ub();
795          for (int i=0;i<n;i++) {
796            double value = element[i];
797            element2[column[i]]=value;
798            norm += value*value;
799          }
800          int kkk = CoinMin(nCuts,k+5);
801          for (int kk=k+1;kk<kkk;kk++) { 
802            int jj=which[kk];
803            const OsiRowCut * thisCut2 = cs.rowCutPtr(jj) ;
804            if (thisCut2->lb()>thisCut2->ub()) 
805              break; // cut is infeasible
806            int nB=thisCut2->row().getNumElements();
807            const int * columnB = thisCut2->row().getIndices();
808            const double * elementB = thisCut2->row().getElements();
809            assert (nB);
810            double normB=0.0;
811            double product=0.0;
812            for (int i=0;i<nB;i++) {
813              double value = elementB[i];
814              normB += value*value;
815              product += value*element2[columnB[i]];
816            }
817            if (product>0.0&&product*product>testValue*norm*normB) {
818              bool parallel=true;
819              double lbB = thisCut2->lb();
820              double ubB = thisCut2->ub();
821              if ((lb<-1.0e20&&lbB>-1.0e20)||
822                  (lbB<-1.0e20&&lb>-1.0e20))
823                parallel = false;
824              double tolerance;
825              tolerance = CoinMax(fabs(lb),fabs(lbB))+1.0e-6;
826              if (fabs(lb-lbB)>tolerance)
827                parallel=false;
828              if ((ub>1.0e20&&ubB<1.0e20)||
829                  (ubB>1.0e20&&ub<1.0e20))
830                parallel = false;
831              tolerance = CoinMax(fabs(ub),fabs(ubB))+1.0e-6;
832              if (fabs(ub-ubB)>tolerance)
833                parallel=false;
834              if (parallel) {
835                nParallel++;
836                sort[k]=0.0;
837                break;
838              }
839            }
840          }
841          for (int i=0;i<n;i++) {
842            element2[column[i]]=0.0;
843          }
844        }
845        delete [] element2;
846        CoinSort_2(sort,sort+nCuts,which);
847        k=0;
848        while (nDelete>0||!sort[k]) {
849          int iCut=which[k];
850          const OsiRowCut * thisCut = cs.rowCutPtr(iCut) ;
851          int n=thisCut->row().getNumElements();
852          // may be best, just to save if short
853          if (false&&n&&sort[k]) {
854            // add to saved cuts
855            savedCuts_.insert(*thisCut);
856          }
857          nDelete-=n; 
858          k++;
859          if (k>=nCuts)
860            break;
861        }
862        std::sort(which,which+k);
863        k--;
864        for (;k>=0;k--) {
865          cs.eraseRowCut(which[k]);
866        }
867        delete [] sort;
868        delete [] which;
869        numberRowCutsAfter = cs.sizeRowCuts() ;
870#else
871        double * norm = new double [nCuts];
872        int * which = new int [2*nCuts];
873        double * score = new double [nCuts];
874        double * ortho = new double [nCuts];
875        int nIn=0;
876        int nOut=nCuts;
877        // For parallel cuts
878        double * element2 = new double [numberColumns];
879        const double *objective = solver->getObjCoefficients() ;
880        double objNorm=0.0;
881        for (int i=0;i<numberColumns;i++)
882          objNorm += objective[i]*objective[i];
883        if (objNorm)
884          objNorm = 1.0/sqrt(objNorm);
885        else
886          objNorm=1.0;
887        objNorm *= 0.1; // weight of 0.1
888        CoinZeroN(element2,numberColumns);
889        int numberRowCuts = numberRowCutsAfter-numberRowCutsBefore;
890        int iBest=-1;
891        double best=0.0;
892        int nPossible=0;
893        double testValue = (depth>1) ? 0.7 : 0.5;
894        for (k = 0;k<numberRowCuts;k++) {
895          const OsiRowCut * thisCut = cs.rowCutPtr(k+numberRowCutsBefore) ;
896          double sum=0.0;
897          if (thisCut->lb()<=thisCut->ub()) {
898            int n=thisCut->row().getNumElements();
899            const int * column = thisCut->row().getIndices();
900            const double * element = thisCut->row().getElements();
901            assert (n);
902            double normThis=1.0e-6;
903            double obj=0.0;
904            for (int i=0;i<n;i++) {
905              int iColumn = column[i];
906              double value = element[i];
907              sum += value*solution[iColumn];
908              normThis += value*value;
909              obj += value*objective[iColumn];
910            }
911            if (sum>thisCut->ub()) {
912              sum= sum-thisCut->ub();
913            } else if (sum<thisCut->lb()) {
914              sum= thisCut->lb()-sum;
915            } else {
916              sum=0.0;
917            }
918            if (sum) {
919              normThis =1.0/sqrt(normThis);
920              norm[k]=normThis;
921              sum *= normThis;
922              obj *= normThis;
923              score[k]=sum+obj*objNorm;
924              ortho[k]=1.0;
925            }
926          } else {
927            // keep and discard others
928            nIn=1;
929            which[0]=k;
930            for (int j = 0;j<numberRowCuts;j++) {
931              if (j!=k)
932                which[nOut++]=j;
933            }
934            iBest=-1;
935            break;
936          }
937          if (sum) {
938            if (score[k]>best) {
939              best=score[k];
940              iBest=nPossible;
941            }
942            which[nPossible++]=k;
943          } else {
944            which[nOut++]=k;
945          }
946        }
947        while (iBest>=0) {
948          int kBest=which[iBest];
949          int j=which[nIn];
950          which[iBest]=j;
951          which[nIn++]=kBest;
952          const OsiRowCut * thisCut = cs.rowCutPtr(kBest+numberRowCutsBefore) ;
953          int n=thisCut->row().getNumElements();
954          nReasonable -= n;
955          if (nReasonable<=0) {
956            for (k=nIn;k<nPossible;k++)
957              which[nOut++]=which[k];
958            break;
959          }
960          // Now see which ones are too similar and choose next
961          iBest=-1;
962          best=0.0;
963          int nOld=nPossible;
964          nPossible=nIn;
965          const int * column = thisCut->row().getIndices();
966          const double * element = thisCut->row().getElements();
967          assert (n);
968          double normNew = norm[kBest];
969          for (int i=0;i<n;i++) {
970            double value = element[i];
971            element2[column[i]]=value;
972          }
973          for (int j=nIn;j<nOld;j++) {
974            k = which[j];
975            const OsiRowCut * thisCut2 = cs.rowCutPtr(k+numberRowCutsBefore) ;
976            int nB=thisCut2->row().getNumElements();
977            const int * columnB = thisCut2->row().getIndices();
978            const double * elementB = thisCut2->row().getElements();
979            assert (nB);
980            double normB=norm[k];
981            double product=0.0;
982            for (int i=0;i<nB;i++) {
983              double value = elementB[i];
984              product += value*element2[columnB[i]];
985            }
986            double orthoScore = 1.0-product*normNew*normB;
987            if (orthoScore>=testValue) {
988              ortho[k]=CoinMin(orthoScore,ortho[k]);
989              double test=score[k]+ortho[k];
990              if (test>best) {
991                best=score[k];
992                iBest=nPossible;
993              }
994              which[nPossible++]=k;
995            } else {
996              which[nOut++]=k;
997            }
998          }
999          for (int i=0;i<n;i++) {
1000            element2[column[i]]=0.0;
1001          }
1002        }
1003        delete [] score;
1004        delete [] ortho;
1005        std::sort(which+nCuts,which+nOut);
1006        k=nOut-1;
1007        for (;k>=nCuts;k--) {
1008          cs.eraseRowCut(which[k]+numberRowCutsBefore);
1009        }
1010        delete [] norm;
1011        delete [] which;
1012        numberRowCutsAfter = cs.sizeRowCuts() ;
1013#endif
1014      }
1015    }
1016#ifdef CBC_DEBUG
1017    {
1018      int numberRowCutsAfter = cs.sizeRowCuts() ;
1019      int k ;
1020      int nBad=0;
1021      for (k = numberRowCutsBefore;k<numberRowCutsAfter;k++) {
1022        OsiRowCut thisCut = cs.rowCut(k) ;
1023        if (thisCut.lb()>thisCut.ub()||
1024            thisCut.lb()>1.0e8||
1025            thisCut.ub()<-1.0e8)
1026          printf("cut from %s has bounds %g and %g!\n",
1027                 generatorName_,thisCut.lb(),thisCut.ub());
1028        if (thisCut.lb()<=thisCut.ub()) {
1029          /* check size of elements.
1030             We can allow smaller but this helps debug generators as it
1031             is unsafe to have small elements */
1032          int n=thisCut.row().getNumElements();
1033          const int * column = thisCut.row().getIndices();
1034          const double * element = thisCut.row().getElements();
1035          assert (n);
1036          for (int i=0;i<n;i++) {
1037            double value = element[i];
1038            if (fabs(value)<=1.0e-12||fabs(value)>=1.0e20)
1039              nBad++;
1040          }
1041        }
1042        if (nBad) 
1043          printf("Cut generator %s produced %d cuts of which %d had tiny or large elements\n",
1044                 generatorName_,numberRowCutsAfter-numberRowCutsBefore,nBad);
1045      }
1046    }
1047#endif
1048    {
1049      int numberRowCutsAfter = cs.sizeRowCuts() ;
1050      if (numberRowCutsBefore<numberRowCutsAfter) {
1051        for (int k = numberRowCutsBefore;k<numberRowCutsAfter;k++) {
1052          OsiRowCut thisCut = cs.rowCut(k) ;
1053          int n=thisCut.row().getNumElements();
1054          numberElements_ += n;
1055        }
1056#if 0
1057        printf("generator %s generated %d row cuts\n",
1058               generatorName_,numberRowCutsAfter-numberRowCutsBefore);
1059#endif
1060        numberCuts_ += numberRowCutsAfter-numberRowCutsBefore;
1061      }
1062      int numberColumnCutsAfter = cs.sizeColCuts() ;
1063      if (numberColumnCutsBefore<numberColumnCutsAfter) {
1064#if 0
1065        printf("generator %s generated %d column cuts\n",
1066               generatorName_,numberColumnCutsAfter-numberColumnCutsBefore);
1067#endif
1068        numberColumnCuts_ += numberColumnCutsAfter-numberColumnCutsBefore;
1069      }
1070    }
1071    if (timing())
1072      timeInCutGenerator_ += CoinCpuTime()-time1;
1073#if 0
1074    // switch off if first time and no good
1075    if (node==NULL&&!pass) {
1076      if (cs.sizeCuts()-cutsBefore<CoinAbs(switchOffIfLessThan_)) {
1077        whenCutGenerator_=-99;
1078        whenCutGeneratorInSub_ = -200;
1079      }
1080    }
1081#endif
1082  }
1083  return returnCode;
1084}
1085void 
1086CbcCutGenerator::setHowOften(int howOften) 
1087{
1088 
1089  if (howOften>=1000000) {
1090    // leave Probing every SCANCUTS_PROBING
1091    howOften = howOften % 1000000;
1092    CglProbing* generator =
1093      dynamic_cast<CglProbing*>(generator_);
1094   
1095    if (generator&&howOften>SCANCUTS_PROBING) 
1096      howOften=SCANCUTS_PROBING+1000000;
1097    else
1098      howOften += 1000000;
1099  }
1100  whenCutGenerator_ = howOften;
1101}
1102void 
1103CbcCutGenerator::setWhatDepth(int value) 
1104{
1105  depthCutGenerator_ = value;
1106}
1107void 
1108CbcCutGenerator::setWhatDepthInSub(int value) 
1109{
1110  depthCutGeneratorInSub_ = value;
1111}
1112// Create C++ lines to get to current state
1113void 
1114CbcCutGenerator::generateTuning( FILE * fp) 
1115{
1116  fprintf(fp,"// Cbc tuning for generator %s\n",generatorName_);
1117  fprintf(fp,"   generator->setHowOften(%d);\n",whenCutGenerator_);
1118  fprintf(fp,"   generator->setSwitchOffIfLessThan(%d);\n",switchOffIfLessThan_);
1119  fprintf(fp,"   generator->setWhatDepth(%d);\n",depthCutGenerator_);
1120  fprintf(fp,"   generator->setInaccuracy(%d);\n",inaccuracy_);
1121  if (timing())
1122    fprintf(fp,"   generator->setTiming(true);\n");
1123  if (normal())
1124    fprintf(fp,"   generator->setNormal(true);\n");
1125  if (atSolution())
1126    fprintf(fp,"   generator->setAtSolution(true);\n");
1127  if (whenInfeasible())
1128    fprintf(fp,"   generator->setWhenInfeasible(true);\n");
1129  if (needsOptimalBasis())
1130    fprintf(fp,"   generator->setNeedsOptimalBasis(true);\n");
1131  if (mustCallAgain())
1132    fprintf(fp,"   generator->setMustCallAgain(true);\n");
1133  if (whetherToUse())
1134    fprintf(fp,"   generator->setWhetherToUse(true);\n");
1135}
1136
1137
1138// Default Constructor
1139CbcCutModifier::CbcCutModifier() 
1140{
1141}
1142
1143
1144// Destructor
1145CbcCutModifier::~CbcCutModifier ()
1146{
1147}
1148
1149// Copy constructor
1150CbcCutModifier::CbcCutModifier ( const CbcCutModifier & /*rhs*/)
1151{
1152}
1153
1154// Assignment operator
1155CbcCutModifier & 
1156CbcCutModifier::operator=( const CbcCutModifier& rhs)
1157{
1158  if (this!=&rhs) {
1159  }
1160  return *this;
1161}
1162
1163// Default Constructor
1164CbcCutSubsetModifier::CbcCutSubsetModifier ()
1165  : CbcCutModifier(),
1166    firstOdd_(COIN_INT_MAX)
1167{
1168}
1169
1170// Useful constructor
1171CbcCutSubsetModifier::CbcCutSubsetModifier (int firstOdd)
1172  : CbcCutModifier()
1173{
1174  firstOdd_=firstOdd;
1175}
1176
1177// Copy constructor
1178CbcCutSubsetModifier::CbcCutSubsetModifier ( const CbcCutSubsetModifier & rhs)
1179  :CbcCutModifier(rhs)
1180{
1181  firstOdd_ = rhs.firstOdd_;
1182}
1183
1184// Clone
1185CbcCutModifier *
1186CbcCutSubsetModifier::clone() const
1187{
1188  return new CbcCutSubsetModifier(*this);
1189}
1190
1191// Assignment operator
1192CbcCutSubsetModifier & 
1193CbcCutSubsetModifier::operator=( const CbcCutSubsetModifier& rhs)
1194{
1195  if (this!=&rhs) {
1196    CbcCutModifier::operator=(rhs);
1197    firstOdd_ = rhs.firstOdd_;
1198  }
1199  return *this;
1200}
1201
1202// Destructor
1203CbcCutSubsetModifier::~CbcCutSubsetModifier ()
1204{
1205}
1206/* Returns
1207   0 unchanged
1208   1 strengthened
1209   2 weakened
1210   3 deleted
1211*/
1212int 
1213CbcCutSubsetModifier::modify(const OsiSolverInterface * /*solver*/, 
1214                             OsiRowCut & cut) 
1215{
1216  int n=cut.row().getNumElements();
1217  if (!n)
1218    return 0;
1219  const int * column = cut.row().getIndices();
1220  //const double * element = cut.row().getElements();
1221  int returnCode=0;
1222  for (int i=0;i<n;i++) {
1223    if (column[i]>=firstOdd_) {
1224      returnCode=3;
1225      break;
1226    }
1227  }
1228  if (!returnCode) {
1229    const double * element = cut.row().getElements();
1230    printf("%g <= ",cut.lb());
1231    for (int i=0;i<n;i++) {
1232      printf("%g*x%d ",element[i],column[i]);
1233    }
1234    printf("<= %g\n",cut.ub());
1235  }
1236  //return 3;
1237  return returnCode;
1238}
1239
Note: See TracBrowser for help on using the repository browser.