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

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

final changes before cleaning

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