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

Last change on this file since 1433 was 1393, checked in by lou, 10 years ago

Mark #if 0 with JJF_ZERO and #if 1 with JJF_ONE as a historical reference
point.

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