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

Last change on this file since 1750 was 1656, checked in by forrest, 8 years ago

allow end cuts and lagomory

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