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

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

fixes

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