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

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

chnages to try and make faster

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