source: branches/sandbox/Cbc/src/CbcCutGenerator.cpp @ 1314

Last change on this file since 1314 was 1314, checked in by EdwinStraver, 10 years ago

Broke up CbcCompareActual?.cpp into CbcCompareDepth?, CbcCompareDefault?, CbcCompareObjective? and CbcCompareEstimate?.
Carved CbcCutModifier? and CbcCutSubsetModifier? out of CbcCutGenerator?.
Updated spreadsheets.

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