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

Last change on this file since 1573 was 1573, checked in by lou, 8 years ago

Change to EPL license notice.

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