source: branches/sandbox/Cbc/src/CbcHeuristic.cpp @ 1393

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

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

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 109.8 KB
RevLine 
[1271]1/* $Id: CbcHeuristic.cpp 1393 2009-12-11 00:29:38Z lou $ */
[2]2// Copyright (C) 2002, 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
[311]8
[325]9#include "CbcConfig.h"
[311]10
[2]11#include <cassert>
[904]12#include <cstdlib>
[2]13#include <cmath>
14#include <cfloat>
15
[912]16//#define PRINT_DEBUG
[311]17#ifdef COIN_HAS_CLP
[197]18#include "OsiClpSolverInterface.hpp"
[246]19#endif
[2]20#include "CbcModel.hpp"
21#include "CbcMessage.hpp"
22#include "CbcHeuristic.hpp"
[833]23#include "CbcHeuristicFPump.hpp"
[197]24#include "CbcStrategy.hpp"
25#include "CglPreProcess.hpp"
[1059]26#include "CglProbing.hpp"
[264]27#include "OsiAuxInfo.hpp"
[833]28#include "OsiPresolve.hpp"
[912]29#include "CbcBranchActual.hpp"
[940]30#include "CbcCutGenerator.hpp"
[912]31//==============================================================================
[2]32
[912]33CbcHeuristicNode::CbcHeuristicNode(const CbcHeuristicNode& rhs)
34{
[1286]35    numObjects_ = rhs.numObjects_;
36    brObj_ = new CbcBranchingObject*[numObjects_];
37    for (int i = 0; i < numObjects_; ++i) {
38        brObj_[i] = rhs.brObj_[i]->clone();
39    }
[912]40}
41
42void
43CbcHeuristicNodeList::gutsOfDelete()
44{
[1286]45    for (int i = nodes_.size() - 1; i >= 0; --i) {
46        delete nodes_[i];
47    }
[912]48}
49
50void
51CbcHeuristicNodeList::gutsOfCopy(const CbcHeuristicNodeList& rhs)
52{
[1286]53    append(rhs);
[912]54}
55
56CbcHeuristicNodeList::CbcHeuristicNodeList(const CbcHeuristicNodeList& rhs)
57{
[1286]58    gutsOfCopy(rhs);
[912]59}
60
61CbcHeuristicNodeList& CbcHeuristicNodeList::operator=
[1286]62(const CbcHeuristicNodeList & rhs)
[912]63{
[1286]64    if (this != &rhs) {
65        gutsOfDelete();
66        gutsOfCopy(rhs);
67    }
68    return *this;
[912]69}
70
71CbcHeuristicNodeList::~CbcHeuristicNodeList()
72{
[1286]73    gutsOfDelete();
[912]74}
75
76void
77CbcHeuristicNodeList::append(CbcHeuristicNode*& node)
78{
[1286]79    nodes_.push_back(node);
80    node = NULL;
[912]81}
82
83void
84CbcHeuristicNodeList::append(const CbcHeuristicNodeList& nodes)
85{
[1286]86    nodes_.reserve(nodes_.size() + nodes.size());
87    for (int i = 0; i < nodes.size(); ++i) {
88        CbcHeuristicNode* node = new CbcHeuristicNode(*nodes.node(i));
89        append(node);
90    }
[912]91}
92
93//==============================================================================
[1271]94#define DEFAULT_WHERE ((255-2-16)*(1+256))
[2]95// Default Constructor
[912]96CbcHeuristic::CbcHeuristic() :
[1286]97        model_(NULL),
98        when_(2),
99        numberNodes_(200),
100        feasibilityPumpOptions_(-1),
101        fractionSmall_(1.0),
102        heuristicName_("Unknown"),
103        howOften_(1),
104        decayFactor_(0.0),
105        switches_(0),
106        whereFrom_(DEFAULT_WHERE),
107        shallowDepth_(1),
108        howOftenShallow_(1),
109        numInvocationsInShallow_(0),
110        numInvocationsInDeep_(0),
111        lastRunDeep_(0),
112        numRuns_(0),
113        minDistanceToRun_(1),
114        runNodes_(),
115        numCouldRun_(0),
116        numberSolutionsFound_(0),
117        inputSolution_(NULL)
[2]118{
[1286]119    // As CbcHeuristic virtual need to modify .cpp if above change
[2]120}
121
122// Constructor from model
[912]123CbcHeuristic::CbcHeuristic(CbcModel & model) :
[1286]124        model_(&model),
125        when_(2),
126        numberNodes_(200),
127        feasibilityPumpOptions_(-1),
128        fractionSmall_(1.0),
129        heuristicName_("Unknown"),
130        howOften_(1),
131        decayFactor_(0.0),
132        switches_(0),
133        whereFrom_(DEFAULT_WHERE),
134        shallowDepth_(1),
135        howOftenShallow_(1),
136        numInvocationsInShallow_(0),
137        numInvocationsInDeep_(0),
138        lastRunDeep_(0),
139        numRuns_(0),
140        minDistanceToRun_(1),
141        runNodes_(),
142        numCouldRun_(0),
143        numberSolutionsFound_(0),
144        inputSolution_(NULL)
[912]145{}
146
147void
148CbcHeuristic::gutsOfCopy(const CbcHeuristic & rhs)
[2]149{
[1286]150    model_ = rhs.model_;
151    when_ = rhs.when_;
152    numberNodes_ = rhs.numberNodes_;
153    feasibilityPumpOptions_ = rhs.feasibilityPumpOptions_;
154    fractionSmall_ = rhs.fractionSmall_;
155    randomNumberGenerator_ = rhs.randomNumberGenerator_;
156    heuristicName_ = rhs.heuristicName_;
157    howOften_ = rhs.howOften_;
158    decayFactor_ = rhs.decayFactor_;
159    switches_ = rhs.switches_;
160    whereFrom_ = rhs.whereFrom_;
161    shallowDepth_ = rhs.shallowDepth_;
162    howOftenShallow_ = rhs.howOftenShallow_;
163    numInvocationsInShallow_ = rhs.numInvocationsInShallow_;
164    numInvocationsInDeep_ = rhs.numInvocationsInDeep_;
165    lastRunDeep_ = rhs.lastRunDeep_;
166    numRuns_ = rhs.numRuns_;
167    numCouldRun_ = rhs.numCouldRun_;
168    minDistanceToRun_ = rhs.minDistanceToRun_;
169    runNodes_ = rhs.runNodes_;
170    numberSolutionsFound_ = rhs.numberSolutionsFound_;
171    if (rhs.inputSolution_) {
172        int numberColumns = model_->getNumCols();
173        setInputSolution(rhs.inputSolution_, rhs.inputSolution_[numberColumns]);
174    }
[2]175}
[1286]176// Copy constructor
[640]177CbcHeuristic::CbcHeuristic(const CbcHeuristic & rhs)
178{
[1286]179    inputSolution_ = NULL;
180    gutsOfCopy(rhs);
[640]181}
[912]182
[1286]183// Assignment operator
184CbcHeuristic &
185CbcHeuristic::operator=( const CbcHeuristic & rhs)
[640]186{
[1286]187    if (this != &rhs) {
188        gutsOfDelete();
189        gutsOfCopy(rhs);
190    }
191    return *this;
[640]192}
193
[912]194void CbcHeurDebugNodes(CbcModel* model_)
195{
[1286]196    CbcNode* node = model_->currentNode();
197    CbcNodeInfo* nodeInfo = node->nodeInfo();
198    std::cout << "===============================================================\n";
199    while (nodeInfo) {
200        const CbcNode* node = nodeInfo->owner();
201        printf("nodeinfo: node %i\n", nodeInfo->nodeNumber());
202        {
203            const CbcIntegerBranchingObject* brPrint =
204                dynamic_cast<const CbcIntegerBranchingObject*>(nodeInfo->parentBranch());
205            if (!brPrint) {
206                printf("    parentBranch: NULL\n");
207            } else {
208                const double* downBounds = brPrint->downBounds();
209                const double* upBounds = brPrint->upBounds();
210                int variable = brPrint->variable();
211                int way = brPrint->way();
212                printf("   parentBranch: var %i downBd [%i,%i] upBd [%i,%i] way %i\n",
213                       variable, static_cast<int>(downBounds[0]), static_cast<int>(downBounds[1]),
214                       static_cast<int>(upBounds[0]), static_cast<int>(upBounds[1]), way);
215            }
216        }
217        if (! node) {
218            printf("    owner: NULL\n");
219        } else {
220            printf("    owner: node %i depth %i onTree %i active %i",
221                   node->nodeNumber(), node->depth(), node->onTree(), node->active());
222            const OsiBranchingObject* osibr =
223                nodeInfo->owner()->branchingObject();
224            const CbcBranchingObject* cbcbr =
225                dynamic_cast<const CbcBranchingObject*>(osibr);
226            const CbcIntegerBranchingObject* brPrint =
227                dynamic_cast<const CbcIntegerBranchingObject*>(cbcbr);
228            if (!brPrint) {
229                printf("        ownerBranch: NULL\n");
230            } else {
231                const double* downBounds = brPrint->downBounds();
232                const double* upBounds = brPrint->upBounds();
233                int variable = brPrint->variable();
234                int way = brPrint->way();
235                printf("        ownerbranch: var %i downBd [%i,%i] upBd [%i,%i] way %i\n",
236                       variable, static_cast<int>(downBounds[0]), static_cast<int>(downBounds[1]),
237                       static_cast<int>(upBounds[0]), static_cast<int>(upBounds[1]), way);
238            }
239        }
240        nodeInfo = nodeInfo->parent();
[912]241    }
242}
243
244void
245CbcHeuristic::debugNodes()
246{
[1286]247    CbcHeurDebugNodes(model_);
[912]248}
249
250void
251CbcHeuristic::printDistanceToNodes()
252{
[1286]253    const CbcNode* currentNode = model_->currentNode();
254    if (currentNode != NULL) {
255        CbcHeuristicNode* nodeDesc = new CbcHeuristicNode(*model_);
256        for (int i = runNodes_.size() - 1; i >= 0; --i) {
257            nodeDesc->distance(runNodes_.node(i));
258        }
259        runNodes_.append(nodeDesc);
[912]260    }
261}
262
263bool
[1271]264CbcHeuristic::shouldHeurRun(int whereFrom)
[912]265{
[1286]266    assert (whereFrom >= 0 && whereFrom < 16);
267    // take off 8 (code - likes new solution)
268    whereFrom &= 7;
269    if ((whereFrom_&(1 << whereFrom)) == 0)
270        return false;
271    // No longer used for original purpose - so use for ever run at all JJF
[1393]272#ifndef JJF_ONE
[1286]273    // Don't run if hot start
274    if (model_ && model_->hotstartSolution())
275        return false;
276    else
277        return true;
[1013]278#else
[1393]279#ifdef JJF_ZERO
[1286]280    const CbcNode* currentNode = model_->currentNode();
281    if (currentNode == NULL) {
282        return false;
283    }
[912]284
[1286]285    debugNodes();
[912]286//   return false;
287
[1286]288    const int depth = currentNode->depth();
[912]289#else
[1286]290    int depth = model_->currentDepth();
[912]291#endif
292
[1286]293    const int nodeCount = model_->getNodeCount();  // FIXME: check that this is
294    // correct in parallel
[912]295
[1286]296    if (nodeCount == 0 || depth <= shallowDepth_) {
297        // what to do when we are in the shallow part of the tree
298        if (model_->getCurrentPassNumber() == 1) {
299            // first time in the node...
300            numInvocationsInShallow_ = 0;
301        }
302        ++numInvocationsInShallow_;
303        // Very large howOftenShallow_ will give the original test:
304        // (model_->getCurrentPassNumber() != 1)
305        //    if ((numInvocationsInShallow_ % howOftenShallow_) != 1) {
306        if ((numInvocationsInShallow_ % howOftenShallow_) != 0) {
307            return false;
308        }
309        // LL: should we save these nodes in the list of nodes where the heur was
310        // LL: run?
[1393]311#ifndef JJF_ONE
[1286]312        if (currentNode != NULL) {
313            // Get where we are and create the appropriate CbcHeuristicNode object
314            CbcHeuristicNode* nodeDesc = new CbcHeuristicNode(*model_);
315            runNodes_.append(nodeDesc);
316        }
[912]317#endif
[1286]318    } else {
319        // deeper in the tree
320        if (model_->getCurrentPassNumber() == 1) {
321            // first time in the node...
322            ++numInvocationsInDeep_;
323        }
324        if (numInvocationsInDeep_ - lastRunDeep_ < howOften_) {
325            return false;
326        }
327        if (model_->getCurrentPassNumber() != 1) {
328            // Run the heuristic only when first entering the node.
329            // LL: I don't think this is right. It should run just before strong
330            // LL: branching, I believe.
331            return false;
332        }
333        // Get where we are and create the appropriate CbcHeuristicNode object
334        CbcHeuristicNode* nodeDesc = new CbcHeuristicNode(*model_);
335        //#ifdef PRINT_DEBUG
[1393]336#ifndef JJF_ONE
[1286]337        const double minDistanceToRun = 1.5 * log((double)depth) / log((double)2);
[915]338#else
339    const double minDistanceToRun = minDistanceToRun_;
340#endif
341#ifdef PRINT_DEBUG
[1286]342        double minDistance = nodeDesc->minDistance(runNodes_);
343        std::cout << "minDistance = " << minDistance
344                  << ", minDistanceToRun = " << minDistanceToRun << std::endl;
[912]345#endif
[1286]346        if (nodeDesc->minDistanceIsSmall(runNodes_, minDistanceToRun)) {
347            delete nodeDesc;
348            return false;
349        }
350        runNodes_.append(nodeDesc);
351        lastRunDeep_ = numInvocationsInDeep_;
352        //    ++lastRunDeep_;
[912]353    }
[1286]354    ++numRuns_;
355    return true;
[1013]356#endif
[912]357}
358
359bool
360CbcHeuristic::shouldHeurRun_randomChoice()
361{
[1286]362    if (!when_)
363        return false;
364    int depth = model_->currentDepth();
365    // when_ -999 is special marker to force to run
366    if (depth != 0 && when_ != -999) {
367        const double numerator = depth * depth;
368        const double denominator = exp(depth * log(2.0));
369        double probability = numerator / denominator;
370        double randomNumber = randomNumberGenerator_.randomDouble();
371        int when = when_ % 100;
372        if (when > 2 && when < 8) {
373            /* JJF adjustments
374            3 only at root and if no solution
375            4 only at root and if this heuristic has not got solution
[1315]376            5 as 3 but decay more
[1286]377            6 decay
378            7 run up to 2 times if solution found 4 otherwise
379            */
380            switch (when) {
381            case 3:
382            default:
383                if (model_->bestSolution())
384                    probability = -1.0;
385                break;
386            case 4:
387                if (numberSolutionsFound_)
388                    probability = -1.0;
389                break;
390            case 5:
[1315]391                assert (decayFactor_);
392                if (model_->bestSolution()) {
[1286]393                    probability = -1.0;
[1315]394                } else if (numCouldRun_ > 1000) {
395                    decayFactor_ *= 0.99;
396                    probability *= decayFactor_;
397                }
[1286]398                break;
399            case 6:
400                if (depth >= 3) {
401                    if ((numCouldRun_ % howOften_) == 0 &&
402                            numberSolutionsFound_*howOften_ < numCouldRun_) {
[1013]403#ifdef COIN_DEVELOP
[1286]404                        int old = howOften_;
[1013]405#endif
[1315]406                        howOften_ = CoinMin(CoinMax(static_cast<int> (howOften_ * 1.1), howOften_ + 1), 1000000);
[961]407#ifdef COIN_DEVELOP
[1286]408                        printf("Howoften changed from %d to %d for %s\n",
409                               old, howOften_, heuristicName_.c_str());
[961]410#endif
[1286]411                    }
412                    probability = 1.0 / howOften_;
413                    if (model_->bestSolution())
414                        probability *= 0.5;
415                }
416                break;
417            case 7:
418                if ((model_->bestSolution() && numRuns_ >= 2) || numRuns_ >= 4)
419                    probability = -1.0;
420                break;
421            }
422        }
423        if (randomNumber > probability)
424            return false;
425
426        if (model_->getCurrentPassNumber() > 1)
427            return false;
[961]428#ifdef COIN_DEVELOP
[1286]429        printf("Running %s, random %g probability %g\n",
430               heuristicName_.c_str(), randomNumber, probability);
[961]431#endif
[1286]432    } else {
[961]433#ifdef COIN_DEVELOP
[1286]434        printf("Running %s, depth %d when %d\n",
435               heuristicName_.c_str(), depth, when_);
[961]436#endif
[1286]437    }
438    ++numRuns_;
439    return true;
[912]440}
441
[6]442// Resets stuff if model changes
[1286]443void
[6]444CbcHeuristic::resetModel(CbcModel * model)
445{
[1286]446    model_ = model;
447}
[838]448// Set seed
449void
450CbcHeuristic::setSeed(int value)
451{
[1286]452    randomNumberGenerator_.setSeed(value);
[838]453}
[2]454
[640]455// Create C++ lines to get to current state
[1286]456void
457CbcHeuristic::generateCpp( FILE * fp, const char * heuristic)
[640]458{
[1286]459    // hard coded as CbcHeuristic virtual
460    if (when_ != 2)
461        fprintf(fp, "3  %s.setWhen(%d);\n", heuristic, when_);
462    else
463        fprintf(fp, "4  %s.setWhen(%d);\n", heuristic, when_);
464    if (numberNodes_ != 200)
465        fprintf(fp, "3  %s.setNumberNodes(%d);\n", heuristic, numberNodes_);
466    else
467        fprintf(fp, "4  %s.setNumberNodes(%d);\n", heuristic, numberNodes_);
468    if (feasibilityPumpOptions_ != -1)
469        fprintf(fp, "3  %s.setFeasibilityPumpOptions(%d);\n", heuristic, feasibilityPumpOptions_);
470    else
471        fprintf(fp, "4  %s.setFeasibilityPumpOptions(%d);\n", heuristic, feasibilityPumpOptions_);
472    if (fractionSmall_ != 1.0)
473        fprintf(fp, "3  %s.setFractionSmall(%g);\n", heuristic, fractionSmall_);
474    else
475        fprintf(fp, "4  %s.setFractionSmall(%g);\n", heuristic, fractionSmall_);
476    if (heuristicName_ != "Unknown")
477        fprintf(fp, "3  %s.setHeuristicName(\"%s\");\n",
478                heuristic, heuristicName_.c_str()) ;
479    else
480        fprintf(fp, "4  %s.setHeuristicName(\"%s\");\n",
481                heuristic, heuristicName_.c_str()) ;
482    if (decayFactor_ != 0.0)
483        fprintf(fp, "3  %s.setDecayFactor(%g);\n", heuristic, decayFactor_);
484    else
485        fprintf(fp, "4  %s.setDecayFactor(%g);\n", heuristic, decayFactor_);
486    if (switches_ != 0)
487        fprintf(fp, "3  %s.setSwitches(%d);\n", heuristic, switches_);
488    else
489        fprintf(fp, "4  %s.setSwitches(%d);\n", heuristic, switches_);
490    if (whereFrom_ != DEFAULT_WHERE)
491        fprintf(fp, "3  %s.setWhereFrom(%d);\n", heuristic, whereFrom_);
492    else
493        fprintf(fp, "4  %s.setWhereFrom(%d);\n", heuristic, whereFrom_);
494    if (shallowDepth_ != 1)
495        fprintf(fp, "3  %s.setShallowDepth(%d);\n", heuristic, shallowDepth_);
496    else
497        fprintf(fp, "4  %s.setShallowDepth(%d);\n", heuristic, shallowDepth_);
498    if (howOftenShallow_ != 1)
499        fprintf(fp, "3  %s.setHowOftenShallow(%d);\n", heuristic, howOftenShallow_);
500    else
501        fprintf(fp, "4  %s.setHowOftenShallow(%d);\n", heuristic, howOftenShallow_);
502    if (minDistanceToRun_ != 1)
503        fprintf(fp, "3  %s.setMinDistanceToRun(%d);\n", heuristic, minDistanceToRun_);
504    else
505        fprintf(fp, "4  %s.setMinDistanceToRun(%d);\n", heuristic, minDistanceToRun_);
[640]506}
[1286]507// Destructor
[2]508CbcHeuristic::~CbcHeuristic ()
509{
[1286]510    delete [] inputSolution_;
[2]511}
512
513// update model
514void CbcHeuristic::setModel(CbcModel * model)
515{
[1286]516    model_ = model;
[2]517}
[1271]518/* Clone but ..
519   type 0 clone solver, 1 clone continuous solver
520   Add 2 to say without integer variables which are at low priority
521   Add 4 to say quite likely infeasible so give up easily.*/
[1286]522OsiSolverInterface *
[1271]523CbcHeuristic::cloneBut(int type)
524{
[1286]525    OsiSolverInterface * solver;
526    if ((type&1) == 0 || !model_->continuousSolver())
527        solver = model_->solver()->clone();
528    else
529        solver = model_->continuousSolver()->clone();
[1271]530#ifdef COIN_HAS_CLP
[1286]531    OsiClpSolverInterface * clpSolver
[1271]532    = dynamic_cast<OsiClpSolverInterface *> (solver);
533#endif
[1286]534    if ((type&2) != 0) {
535        int n = model_->numberObjects();
536        int priority = model_->continuousPriority();
537        if (priority < COIN_INT_MAX) {
538            for (int i = 0; i < n; i++) {
539                const OsiObject * obj = model_->object(i);
540                const CbcSimpleInteger * thisOne =
541                    dynamic_cast <const CbcSimpleInteger *> (obj);
542                if (thisOne) {
543                    int iColumn = thisOne->columnNumber();
544                    if (thisOne->priority() >= priority)
545                        solver->setContinuous(iColumn);
546                }
547            }
548        }
549#ifdef COIN_HAS_CLP
550        if (clpSolver) {
551            for (int i = 0; i < n; i++) {
552                const OsiObject * obj = model_->object(i);
553                const CbcSimpleInteger * thisOne =
554                    dynamic_cast <const CbcSimpleInteger *> (obj);
555                if (thisOne) {
556                    int iColumn = thisOne->columnNumber();
557                    if (clpSolver->isOptionalInteger(iColumn))
558                        clpSolver->setContinuous(iColumn);
559                }
560            }
561        }
562#endif
[1271]563    }
564#ifdef COIN_HAS_CLP
[1286]565    if ((type&4) != 0 && clpSolver) {
566        int options = clpSolver->getModelPtr()->moreSpecialOptions();
567        clpSolver->getModelPtr()->setMoreSpecialOptions(options | 64);
[1271]568    }
569#endif
[1286]570    return solver;
[1271]571}
[1121]572// Whether to exit at once on gap
[1286]573bool
[1121]574CbcHeuristic::exitNow(double bestObjective) const
575{
[1286]576    if ((switches_&2048) != 0) {
577        // exit may be forced - but unset for next time
578        switches_ &= ~2048;
579        if ((switches_&1024) != 0)
580            return true;
581    } else if ((switches_&1) == 0) {
582        return false;
583    }
584    // See if can stop on gap
585    OsiSolverInterface * solver = model_->solver();
586    double bestPossibleObjective = solver->getObjValue() * solver->getObjSense();
587    double absGap = CoinMax(model_->getAllowableGap(),
588                            model_->getHeuristicGap());
589    double fracGap = CoinMax(model_->getAllowableFractionGap(),
590                             model_->getHeuristicFractionGap());
591    double testGap = CoinMax(absGap, fracGap *
592                             CoinMax(fabs(bestObjective),
593                                     fabs(bestPossibleObjective)));
[1271]594
[1286]595    if (bestObjective - bestPossibleObjective < testGap
596            && model_->getCutoffIncrement() >= 0.0) {
597        return true;
598    } else {
599        return false;
600    }
[1121]601}
[1013]602#ifdef HISTORY_STATISTICS
[833]603extern bool getHistoryStatistics_;
604#endif
[1286]605static double sizeRatio(int numberRowsNow, int numberColumnsNow,
606                        int numberRowsStart, int numberColumnsStart)
[1132]607{
[1286]608    double valueNow;
609    if (numberRowsNow*10 > numberColumnsNow || numberColumnsNow < 200) {
610        valueNow = 2 * numberRowsNow + numberColumnsNow;
611    } else {
612        // long and thin - rows are more important
613        if (numberRowsNow*40 > numberColumnsNow)
614            valueNow = 10 * numberRowsNow + numberColumnsNow;
615        else
616            valueNow = 200 * numberRowsNow + numberColumnsNow;
617    }
618    double valueStart;
619    if (numberRowsStart*10 > numberColumnsStart || numberColumnsStart < 200) {
620        valueStart = 2 * numberRowsStart + numberColumnsStart;
621    } else {
622        // long and thin - rows are more important
623        if (numberRowsStart*40 > numberColumnsStart)
624            valueStart = 10 * numberRowsStart + numberColumnsStart;
625        else
626            valueStart = 200 * numberRowsStart + numberColumnsStart;
627    }
628    //printf("sizeProblem Now %g, %d rows, %d columns\nsizeProblem Start %g, %d rows, %d columns\n",
629    // valueNow,numberRowsNow,numberColumnsNow,
630    // valueStart,numberRowsStart,numberColumnsStart);
631    if (10*numberRowsNow < 8*numberRowsStart)
632        return valueNow / valueStart;
633    else if (10*numberRowsNow < 9*numberRowsStart)
634        return 1.1*(valueNow / valueStart);
635    else if (numberRowsNow < numberRowsStart)
636        return 1.5*(valueNow / valueStart);
[1132]637    else
[1286]638        return 2.0*(valueNow / valueStart);
[1132]639}
640
[1286]641
[197]642// Do mini branch and bound (return 1 if solution)
[1286]643int
644CbcHeuristic::smallBranchAndBound(OsiSolverInterface * solver, int numberNodes,
[197]645                                  double * newSolution, double & newSolutionValue,
646                                  double cutoff, std::string name) const
647{
[1286]648    // size before
649    int shiftRows = 0;
650    if (numberNodes < 0)
651        shiftRows = solver->getNumRows() - numberNodes_;
652    int numberRowsStart = solver->getNumRows() - shiftRows;
653    int numberColumnsStart = solver->getNumCols();
[1020]654#ifdef CLP_INVESTIGATE
[1286]655    printf("%s has %d rows, %d columns\n",
656           name.c_str(), solver->getNumRows(), solver->getNumCols());
[1020]657#endif
[1286]658    // Use this fraction
659    double fractionSmall = fractionSmall_;
660    double before = 2 * numberRowsStart + numberColumnsStart;
661    if (before > 40000.0) {
662        // fairly large - be more conservative
663        double multiplier = 1.0 - 0.3 * CoinMin(100000.0, before - 40000.0) / 100000.0;
664        if (multiplier < 1.0) {
665            fractionSmall *= multiplier;
[1020]666#ifdef CLP_INVESTIGATE
[1286]667            printf("changing fractionSmall from %g to %g for %s\n",
668                   fractionSmall_, fractionSmall, name.c_str());
[931]669#endif
[1286]670        }
[1020]671    }
[311]672#ifdef COIN_HAS_CLP
[1286]673    OsiClpSolverInterface * osiclp = dynamic_cast< OsiClpSolverInterface*> (solver);
674    if (osiclp && (osiclp->specialOptions()&65536) == 0) {
675        // go faster stripes
676        if (osiclp->getNumRows() < 300 && osiclp->getNumCols() < 500) {
677            osiclp->setupForRepeatedUse(2, 0);
678        } else {
679            osiclp->setupForRepeatedUse(0, 0);
680        }
681        // Turn this off if you get problems
682        // Used to be automatically set
683        osiclp->setSpecialOptions(osiclp->specialOptions() | (128 + 64 - 128));
684        ClpSimplex * lpSolver = osiclp->getModelPtr();
685        lpSolver->setSpecialOptions(lpSolver->specialOptions() | 0x01000000); // say is Cbc (and in branch and bound)
686        lpSolver->setSpecialOptions(lpSolver->specialOptions() |
687                                    (/*16384+*/4096 + 512 + 128));
[197]688    }
[246]689#endif
[1013]690#ifdef HISTORY_STATISTICS
[1286]691    getHistoryStatistics_ = false;
[833]692#endif
[1286]693    int status = 0;
694    int logLevel = model_->logLevel();
[863]695#define LEN_PRINT 250
[1286]696    char generalPrint[LEN_PRINT];
697    // Do presolve to see if possible
698    int numberColumns = solver->getNumCols();
699    char * reset = NULL;
700    int returnCode = 1;
701    int saveModelOptions = model_->specialOptions();
702    assert ((saveModelOptions&2048) == 0);
703    model_->setSpecialOptions(saveModelOptions | 2048);
704    {
705        int saveLogLevel = solver->messageHandler()->logLevel();
706        if (saveLogLevel == 1)
707            solver->messageHandler()->setLogLevel(0);
708        OsiPresolve * pinfo = new OsiPresolve();
709        int presolveActions = 0;
710        // Allow dual stuff on integers
711        presolveActions = 1;
712        // Do not allow all +1 to be tampered with
713        //if (allPlusOnes)
714        //presolveActions |= 2;
715        // allow transfer of costs
716        // presolveActions |= 4;
717        pinfo->setPresolveActions(presolveActions);
718        OsiSolverInterface * presolvedModel = pinfo->presolvedModel(*solver, 1.0e-8, true, 2);
719        delete pinfo;
720        // see if too big
721
722        if (presolvedModel) {
723            int afterRows = presolvedModel->getNumRows();
724            int afterCols = presolvedModel->getNumCols();
725            //#define COIN_DEVELOP
[1121]726#ifdef COIN_DEVELOP_z
[1286]727            if (numberNodes < 0) {
728                solver->writeMpsNative("before.mps", NULL, NULL, 2, 1);
729                presolvedModel->writeMpsNative("after1.mps", NULL, NULL, 2, 1);
730            }
[1121]731#endif
[1286]732            delete presolvedModel;
733            double ratio = sizeRatio(afterRows - shiftRows, afterCols,
734                                     numberRowsStart, numberColumnsStart);
735            double after = 2 * afterRows + afterCols;
736            if (ratio > fractionSmall && after > 300 && numberNodes >= 0) {
737                // Need code to try again to compress further using used
738                const int * used =  model_->usedInSolution();
739                int maxUsed = 0;
740                int iColumn;
741                const double * lower = solver->getColLower();
742                const double * upper = solver->getColUpper();
743                for (iColumn = 0; iColumn < numberColumns; iColumn++) {
744                    if (upper[iColumn] > lower[iColumn]) {
745                        if (solver->isBinary(iColumn))
746                            maxUsed = CoinMax(maxUsed, used[iColumn]);
747                    }
748                }
749                if (maxUsed) {
750                    reset = new char [numberColumns];
751                    int nFix = 0;
752                    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
753                        reset[iColumn] = 0;
754                        if (upper[iColumn] > lower[iColumn]) {
755                            if (solver->isBinary(iColumn) && used[iColumn] == maxUsed) {
756                                bool setValue = true;
757                                if (maxUsed == 1) {
758                                    double randomNumber = randomNumberGenerator_.randomDouble();
759                                    if (randomNumber > 0.3)
760                                        setValue = false;
761                                }
762                                if (setValue) {
763                                    reset[iColumn] = 1;
764                                    solver->setColLower(iColumn, 1.0);
765                                    nFix++;
766                                }
767                            }
768                        }
769                    }
770                    pinfo = new OsiPresolve();
771                    presolveActions = 0;
772                    // Allow dual stuff on integers
773                    presolveActions = 1;
774                    // Do not allow all +1 to be tampered with
775                    //if (allPlusOnes)
776                    //presolveActions |= 2;
777                    // allow transfer of costs
778                    // presolveActions |= 4;
779                    pinfo->setPresolveActions(presolveActions);
780                    presolvedModel = pinfo->presolvedModel(*solver, 1.0e-8, true, 2);
781                    delete pinfo;
782                    if (presolvedModel) {
783                        // see if too big
784                        int afterRows2 = presolvedModel->getNumRows();
785                        int afterCols2 = presolvedModel->getNumCols();
786                        delete presolvedModel;
787                        double ratio = sizeRatio(afterRows2 - shiftRows, afterCols2,
788                                                 numberRowsStart, numberColumnsStart);
789                        double after = 2 * afterRows2 + afterCols2;
790                        if (ratio > fractionSmall && (after > 300 || numberNodes < 0)) {
791                            sprintf(generalPrint, "Full problem %d rows %d columns, reduced to %d rows %d columns - %d fixed gives %d, %d - still too large",
792                                    solver->getNumRows(), solver->getNumCols(),
793                                    afterRows, afterCols, nFix, afterRows2, afterCols2);
794                            // If much too big - give up
795                            if (ratio > 0.75)
796                                returnCode = -1;
797                        } else {
798                            sprintf(generalPrint, "Full problem %d rows %d columns, reduced to %d rows %d columns - %d fixed gives %d, %d - ok now",
799                                    solver->getNumRows(), solver->getNumCols(),
800                                    afterRows, afterCols, nFix, afterRows2, afterCols2);
801                        }
802                        model_->messageHandler()->message(CBC_GENERAL, model_->messages())
803                        << generalPrint
804                        << CoinMessageEol;
805                    } else {
806                        returnCode = 2; // infeasible
807                    }
808                }
809            } else if (ratio > fractionSmall && after > 300) {
810                returnCode = -1;
811            }
812        } else {
813            returnCode = 2; // infeasible
814        }
815        solver->messageHandler()->setLogLevel(saveLogLevel);
[833]816    }
[1286]817    if (returnCode == 2 || returnCode == -1) {
818        model_->setSpecialOptions(saveModelOptions);
819        delete [] reset;
[1013]820#ifdef HISTORY_STATISTICS
[1286]821        getHistoryStatistics_ = true;
[833]822#endif
[1286]823        //printf("small no good\n");
824        return returnCode;
[961]825    }
[1286]826    // Reduce printout
827    bool takeHint;
828    OsiHintStrength strength;
829    solver->getHintParam(OsiDoReducePrint, takeHint, strength);
830    solver->setHintParam(OsiDoReducePrint, true, OsiHintTry);
831    solver->setHintParam(OsiDoPresolveInInitial, false, OsiHintTry);
832    solver->setDblParam(OsiDualObjectiveLimit, cutoff*solver->getObjSense());
833    solver->initialSolve();
834    if (solver->isProvenOptimal()) {
835        CglPreProcess process;
836        /* Do not try and produce equality cliques and
837           do up to 2 passes (normally) 5 if restart */
838        int numberPasses = 2;
839        if (numberNodes < 0) {
840            numberPasses = 5;
841            // Say some rows cuts
842            int numberRows = solver->getNumRows();
843            if (numberNodes_ < numberRows && true /* think */) {
844                char * type = new char[numberRows];
845                memset(type, 0, numberNodes_);
846                memset(type + numberNodes_, 1, numberRows - numberNodes_);
847                process.passInRowTypes(type, numberRows);
848                delete [] type;
849            }
850        }
851        if (logLevel <= 1)
852            process.messageHandler()->setLogLevel(0);
853        OsiSolverInterface * solver2 = process.preProcessNonDefault(*solver, false,
854                                       numberPasses);
855        if (!solver2) {
856            if (logLevel > 1)
857                printf("Pre-processing says infeasible\n");
858            returnCode = 2; // so will be infeasible
859        } else {
[1121]860#ifdef COIN_DEVELOP_z
[1286]861            if (numberNodes < 0) {
862                solver2->writeMpsNative("after2.mps", NULL, NULL, 2, 1);
863            }
[1121]864#endif
[1286]865            // see if too big
866            double ratio = sizeRatio(solver2->getNumRows() - shiftRows, solver2->getNumCols(),
867                                     numberRowsStart, numberColumnsStart);
868            double after = 2 * solver2->getNumRows() + solver2->getNumCols();
869            if (ratio > fractionSmall && (after > 300 || numberNodes < 0)) {
870                sprintf(generalPrint, "Full problem %d rows %d columns, reduced to %d rows %d columns - too large",
871                        solver->getNumRows(), solver->getNumCols(),
872                        solver2->getNumRows(), solver2->getNumCols());
873                model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
874                << generalPrint
875                << CoinMessageEol;
876                returnCode = -1;
877                //printf("small no good2\n");
878            } else {
879                sprintf(generalPrint, "Full problem %d rows %d columns, reduced to %d rows %d columns",
880                        solver->getNumRows(), solver->getNumCols(),
881                        solver2->getNumRows(), solver2->getNumCols());
882                model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
883                << generalPrint
884                << CoinMessageEol;
885            }
886            if (returnCode == 1) {
887                solver2->resolve();
888                CbcModel model(*solver2);
889                if (numberNodes >= 0) {
890                    // normal
891                    model.setSpecialOptions(saveModelOptions | 2048);
892                    if (logLevel <= 1)
893                        model.setLogLevel(0);
894                    else
895                        model.setLogLevel(logLevel);
896                    // No small fathoming
897                    model.setFastNodeDepth(-1);
898                    model.setCutoff(cutoff);
899                    model.setMaximumNodes(numberNodes);
900                    model.solver()->setHintParam(OsiDoReducePrint, true, OsiHintTry);
901                    // Lightweight
902                    CbcStrategyDefaultSubTree strategy(model_, 1, 5, 1, 0);
903                    model.setStrategy(strategy);
904                    model.solver()->setIntParam(OsiMaxNumIterationHotStart, 10);
905                    model.setMaximumCutPassesAtRoot(CoinMin(20, CoinAbs(model_->getMaximumCutPassesAtRoot())));
906                    model.setMaximumCutPasses(CoinMin(10, model_->getMaximumCutPasses()));
907                } else {
908                    model.setSpecialOptions(saveModelOptions);
909                    model_->messageHandler()->message(CBC_RESTART, model_->messages())
910                    << solver2->getNumRows() << solver2->getNumCols()
911                    << CoinMessageEol;
912                    // going for full search and copy across more stuff
913                    model.gutsOfCopy(*model_, 2);
914                    for (int i = 0; i < model.numberCutGenerators(); i++) {
915                        model.cutGenerator(i)->setTiming(true);
916                        // Turn on if was turned on
917                        int iOften = model_->cutGenerator(i)->howOften();
[1053]918#ifdef CLP_INVESTIGATE
[1286]919                        printf("Gen %d often %d %d\n",
920                               i, model.cutGenerator(i)->howOften(),
921                               iOften);
[1053]922#endif
[1286]923                        if (iOften > 0)
924                            model.cutGenerator(i)->setHowOften(iOften % 1000000);
925                        if (model_->cutGenerator(i)->howOftenInSub() == -200)
926                            model.cutGenerator(i)->setHowOften(-100);
927                    }
928                    model.setCutoff(cutoff);
929                    // make sure can't do nested search! but allow heuristics
930                    model.setSpecialOptions((model.specialOptions()&(~(512 + 2048))) | 1024);
931                    bool takeHint;
932                    OsiHintStrength strength;
933                    // Switch off printing if asked to
934                    model_->solver()->getHintParam(OsiDoReducePrint, takeHint, strength);
935                    model.solver()->setHintParam(OsiDoReducePrint, takeHint, strength);
936                    CbcStrategyDefault strategy(1, model_->numberStrong(),
937                                                model_->numberBeforeTrust());
938                    // Set up pre-processing - no
939                    strategy.setupPreProcessing(0); // was (4);
940                    model.setStrategy(strategy);
941                    //model.solver()->writeMps("crunched");
942                    int numberCuts = process.cuts().sizeRowCuts();
943                    if (numberCuts) {
944                        // add in cuts
945                        CglStored cuts = process.cuts();
946                        model.addCutGenerator(&cuts, 1, "Stored from first");
947                    }
948                }
949                // Do search
950                if (logLevel > 1)
951                    model_->messageHandler()->message(CBC_START_SUB, model_->messages())
952                    << name
953                    << model.getMaximumNodes()
954                    << CoinMessageEol;
955                // probably faster to use a basis to get integer solutions
956                model.setSpecialOptions(model.specialOptions() | 2);
[687]957#ifdef CBC_THREAD
[1286]958                if (model_->getNumberThreads() > 0 && (model_->getThreadMode()&4) != 0) {
959                    // See if at root node
960                    bool atRoot = model_->getNodeCount() == 0;
961                    int passNumber = model_->getCurrentPassNumber();
962                    if (atRoot && passNumber == 1)
963                        model.setNumberThreads(model_->getNumberThreads());
964                }
[687]965#endif
[1286]966                model.setParentModel(*model_);
967                model.setOriginalColumns(process.originalColumns());
968                model.setSearchStrategy(-1);
969                // If no feasibility pump then insert a lightweight one
970                if (feasibilityPumpOptions_ >= 0) {
971                    bool gotPump = false;
972                    for (int i = 0; i < model.numberHeuristics(); i++) {
973                        const CbcHeuristicFPump* pump =
974                            dynamic_cast<const CbcHeuristicFPump*>(model.heuristic(i));
975                        if (pump)
976                            gotPump = true;
977                    }
978                    if (!gotPump) {
979                        CbcHeuristicFPump heuristic4;
980                        heuristic4.setMaximumPasses(10);
981                        int pumpTune = feasibilityPumpOptions_;
982                        if (pumpTune > 0) {
983                            /*
984                            >=10000000 for using obj
985                            >=1000000 use as accumulate switch
986                            >=1000 use index+1 as number of large loops
987                            >=100 use 0.05 objvalue as increment
988                            %100 == 10,20 etc for experimentation
989                            1 == fix ints at bounds, 2 fix all integral ints, 3 and continuous at bounds
990                            4 and static continuous, 5 as 3 but no internal integers
991                            6 as 3 but all slack basis!
992                            */
993                            double value = solver2->getObjSense() * solver2->getObjValue();
994                            int w = pumpTune / 10;
995                            int ix = w % 10;
996                            w /= 10;
997                            int c = w % 10;
998                            w /= 10;
999                            int r = w;
1000                            int accumulate = r / 1000;
1001                            r -= 1000 * accumulate;
1002                            if (accumulate >= 10) {
1003                                int which = accumulate / 10;
1004                                accumulate -= 10 * which;
1005                                which--;
1006                                // weights and factors
1007                                double weight[] = {0.1, 0.1, 0.5, 0.5, 1.0, 1.0, 5.0, 5.0};
1008                                double factor[] = {0.1, 0.5, 0.1, 0.5, 0.1, 0.5, 0.1, 0.5};
1009                                heuristic4.setInitialWeight(weight[which]);
1010                                heuristic4.setWeightFactor(factor[which]);
1011                            }
1012                            // fake cutoff
1013                            if (c) {
1014                                double cutoff;
1015                                solver2->getDblParam(OsiDualObjectiveLimit, cutoff);
1016                                cutoff = CoinMin(cutoff, value + 0.1 * fabs(value) * c);
1017                                heuristic4.setFakeCutoff(cutoff);
1018                            }
1019                            if (r) {
1020                                // also set increment
1021                                //double increment = (0.01*i+0.005)*(fabs(value)+1.0e-12);
1022                                double increment = 0.0;
1023                                heuristic4.setAbsoluteIncrement(increment);
1024                                heuristic4.setAccumulate(accumulate);
1025                                heuristic4.setMaximumRetries(r + 1);
1026                            }
1027                            pumpTune = pumpTune % 100;
1028                            if (pumpTune == 6)
1029                                pumpTune = 13;
1030                            if (pumpTune != 13)
1031                                pumpTune = pumpTune % 10;
1032                            heuristic4.setWhen(pumpTune);
1033                            if (ix) {
1034                                heuristic4.setFeasibilityPumpOptions(ix*10);
1035                            }
1036                        }
1037                        model.addHeuristic(&heuristic4, "feasibility pump", 0);
1038                    }
1039                }
1040                //printf("sol %x\n",inputSolution_);
1041                if (inputSolution_) {
1042                    // translate and add a serendipity heuristic
1043                    int numberColumns = solver2->getNumCols();
1044                    const int * which = process.originalColumns();
1045                    OsiSolverInterface * solver3 = solver2->clone();
1046                    for (int i = 0; i < numberColumns; i++) {
1047                        if (solver3->isInteger(i)) {
1048                            int k = which[i];
1049                            double value = inputSolution_[k];
1050                            //if (value)
1051                            //printf("orig col %d now %d val %g\n",
1052                            //       k,i,value);
1053                            solver3->setColLower(i, value);
1054                            solver3->setColUpper(i, value);
1055                        }
1056                    }
1057                    solver3->setDblParam(OsiDualObjectiveLimit, COIN_DBL_MAX);
1058                    solver3->resolve();
1059                    if (!solver3->isProvenOptimal()) {
1060                        // Try just setting nonzeros
1061                        OsiSolverInterface * solver4 = solver2->clone();
1062                        for (int i = 0; i < numberColumns; i++) {
1063                            if (solver4->isInteger(i)) {
1064                                int k = which[i];
1065                                double value = floor(inputSolution_[k] + 0.5);
1066                                if (value) {
1067                                    solver3->setColLower(i, value);
1068                                    solver3->setColUpper(i, value);
1069                                }
1070                            }
1071                        }
1072                        solver4->setDblParam(OsiDualObjectiveLimit, COIN_DBL_MAX);
1073                        solver4->resolve();
1074                        int nBad = -1;
1075                        if (solver4->isProvenOptimal()) {
1076                            nBad = 0;
1077                            const double * solution = solver4->getColSolution();
1078                            for (int i = 0; i < numberColumns; i++) {
1079                                if (solver4->isInteger(i)) {
1080                                    double value = floor(solution[i] + 0.5);
1081                                    if (fabs(value - solution[i]) > 1.0e-6)
1082                                        nBad++;
1083                                }
1084                            }
1085                        }
1086                        if (nBad) {
1087                            delete solver4;
1088                        } else {
1089                            delete solver3;
1090                            solver3 = solver4;
1091                        }
1092                    }
1093                    if (solver3->isProvenOptimal()) {
1094                        // good
1095                        CbcSerendipity heuristic(model);
1096                        double value = solver3->getObjSense() * solver3->getObjValue();
1097                        heuristic.setInputSolution(solver3->getColSolution(), value);
1098                        model.setCutoff(value + 1.0e-7*(1.0 + fabs(value)));
1099                        model.addHeuristic(&heuristic, "Previous solution", 0);
1100                        //printf("added seren\n");
1101                    } else {
1102                        double value = model_->getMinimizationObjValue();
1103                        model.setCutoff(value + 1.0e-7*(1.0 + fabs(value)));
[1132]1104#ifdef CLP_INVESTIGATE
[1286]1105                        printf("NOT added seren\n");
1106                        solver3->writeMps("bad_seren");
1107                        solver->writeMps("orig_seren");
[1150]1108#endif
[1286]1109                    }
1110                    delete solver3;
1111                }
1112                if (model_->searchStrategy() == 2) {
1113                    model.setNumberStrong(5);
1114                    model.setNumberBeforeTrust(5);
1115                }
1116                if (model.getNumCols()) {
1117                    if (numberNodes >= 0) {
1118                        setCutAndHeuristicOptions(model);
1119                        // not too many iterations
1120                        model.setMaximumNumberIterations(100*(numberNodes + 10));
1121                        // Not fast stuff
1122                        model.setFastNodeDepth(-1);
1123                    } else if (model.fastNodeDepth() >= 1000000) {
1124                        // already set
1125                        model.setFastNodeDepth(model.fastNodeDepth() - 1000000);
1126                    }
1127                    model.setWhenCuts(999998);
1128                    model.branchAndBound();
[1013]1129#ifdef COIN_DEVELOP
[1286]1130                    printf("sub branch %d nodes, %d iterations - max %d\n",
1131                           model.getNodeCount(), model.getIterationCount(),
1132                           100*(numberNodes + 10));
[1013]1133#endif
[1286]1134                    if (numberNodes < 0) {
1135                        model_->incrementIterationCount(model.getIterationCount());
1136                        model_->incrementNodeCount(model.getNodeCount());
1137                        for (int iGenerator = 0; iGenerator < model.numberCutGenerators(); iGenerator++) {
1138                            CbcCutGenerator * generator = model.cutGenerator(iGenerator);
1139                            sprintf(generalPrint,
1140                                    "%s was tried %d times and created %d cuts of which %d were active after adding rounds of cuts (%.3f seconds)",
1141                                    generator->cutGeneratorName(),
1142                                    generator->numberTimesEntered(),
1143                                    generator->numberCutsInTotal() +
1144                                    generator->numberColumnCuts(),
1145                                    generator->numberCutsActive(),
1146                                    generator->timeInCutGenerator());
1147                            CglStored * stored = dynamic_cast<CglStored*>(generator->generator());
1148                            if (stored && !generator->numberCutsInTotal())
1149                                continue;
[1059]1150#ifndef CLP_INVESTIGATE
[1286]1151                            CglImplication * implication = dynamic_cast<CglImplication*>(generator->generator());
1152                            if (implication)
1153                                continue;
[1059]1154#endif
[1286]1155                            model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
1156                            << generalPrint
1157                            << CoinMessageEol;
1158                        }
1159                    }
1160                } else {
1161                    // empty model
1162                    model.setMinimizationObjValue(model.solver()->getObjSense()*model.solver()->getObjValue());
1163                }
1164                if (logLevel > 1)
1165                    model_->messageHandler()->message(CBC_END_SUB, model_->messages())
1166                    << name
1167                    << CoinMessageEol;
1168                if (model.getMinimizationObjValue() < CoinMin(cutoff, 1.0e30)) {
1169                    // solution
1170                    if (model.getNumCols())
1171                        returnCode = model.isProvenOptimal() ? 3 : 1;
1172                    else
1173                        returnCode = 3;
1174                    // post process
[700]1175#ifdef COIN_HAS_CLP
[1286]1176                    OsiClpSolverInterface * clpSolver = dynamic_cast< OsiClpSolverInterface*> (model.solver());
1177                    if (clpSolver) {
1178                        ClpSimplex * lpSolver = clpSolver->getModelPtr();
1179                        lpSolver->setSpecialOptions(lpSolver->specialOptions() | 0x01000000); // say is Cbc (and in branch and bound)
1180                    }
[700]1181#endif
[1286]1182                    process.postProcess(*model.solver());
1183                    if (solver->isProvenOptimal() && solver->getObjValue() < cutoff) {
1184                        // Solution now back in solver
1185                        int numberColumns = solver->getNumCols();
1186                        memcpy(newSolution, solver->getColSolution(),
1187                               numberColumns*sizeof(double));
1188                        newSolutionValue = model.getMinimizationObjValue();
1189                    } else {
1190                        // odd - but no good
1191                        returnCode = 0; // so will be infeasible
1192                    }
1193                } else {
1194                    // no good
1195                    returnCode = model.isProvenInfeasible() ? 2 : 0; // so will be infeasible
1196                }
1197                int totalNumberIterations = model.getIterationCount() +
1198                                            process.numberIterationsPre() +
1199                                            process.numberIterationsPost();
1200                if (totalNumberIterations > 100*(numberNodes + 10)) {
1201                    // only allow smaller problems
1202                    fractionSmall = fractionSmall_;
1203                    fractionSmall_ *= 0.9;
[1013]1204#ifdef CLP_INVESTIGATE
[1286]1205                    printf("changing fractionSmall from %g to %g for %s as %d iterations\n",
1206                           fractionSmall, fractionSmall_, name.c_str(), totalNumberIterations);
[1013]1207#endif
[1286]1208                }
1209                if (model.status() == 5)
1210                    returnCode = -2; // stop
1211                if (model.isProvenInfeasible())
1212                    status = 1;
1213                else if (model.isProvenOptimal())
1214                    status = 2;
1215            }
1216        }
1217    } else {
1218        returnCode = 2; // infeasible finished
[197]1219    }
[1286]1220    model_->setSpecialOptions(saveModelOptions);
1221    model_->setLogLevel(logLevel);
1222    if (returnCode == 1 || returnCode == 2) {
1223        OsiSolverInterface * solverC = model_->continuousSolver();
1224        if (false && solverC) {
1225            const double * lower = solver->getColLower();
1226            const double * upper = solver->getColUpper();
1227            const double * lowerC = solverC->getColLower();
1228            const double * upperC = solverC->getColUpper();
1229            bool good = true;
1230            for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
1231                if (solverC->isInteger(iColumn)) {
1232                    if (lower[iColumn] > lowerC[iColumn] &&
1233                            upper[iColumn] < upperC[iColumn]) {
1234                        good = false;
1235                        printf("CUT - can't add\n");
1236                        break;
1237                    }
1238                }
1239            }
1240            if (good) {
1241                double * cut = new double [numberColumns];
1242                int * which = new int [numberColumns];
1243                double rhs = -1.0;
1244                int n = 0;
1245                for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
1246                    if (solverC->isInteger(iColumn)) {
1247                        if (lower[iColumn] == upperC[iColumn]) {
1248                            rhs += lower[iColumn];
1249                            cut[n] = 1.0;
1250                            which[n++] = iColumn;
1251                        } else if (upper[iColumn] == lowerC[iColumn]) {
1252                            rhs -= upper[iColumn];
1253                            cut[n] = -1.0;
1254                            which[n++] = iColumn;
1255                        }
1256                    }
1257                }
1258                printf("CUT has %d entries\n", n);
1259                OsiRowCut newCut;
1260                newCut.setLb(-COIN_DBL_MAX);
1261                newCut.setUb(rhs);
1262                newCut.setRow(n, which, cut, false);
1263                model_->makeGlobalCut(newCut);
1264                delete [] cut;
1265                delete [] which;
1266            }
1267        }
[1271]1268#ifdef COIN_DEVELOP
[1286]1269        if (status == 1)
1270            printf("heuristic could add cut because infeasible (%s)\n", heuristicName_.c_str());
1271        else if (status == 2)
1272            printf("heuristic could add cut because optimal (%s)\n", heuristicName_.c_str());
[1271]1273#endif
[833]1274    }
[1286]1275    if (reset) {
1276        for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
1277            if (reset[iColumn])
1278                solver->setColLower(iColumn, 0.0);
1279        }
1280        delete [] reset;
1281    }
[1013]1282#ifdef HISTORY_STATISTICS
[1286]1283    getHistoryStatistics_ = true;
[1013]1284#endif
[1286]1285    solver->setHintParam(OsiDoReducePrint, takeHint, strength);
1286    return returnCode;
[197]1287}
[940]1288// Set input solution
[1286]1289void
[940]1290CbcHeuristic::setInputSolution(const double * solution, double objValue)
1291{
[1286]1292    delete [] inputSolution_;
1293    inputSolution_ = NULL;
1294    if (model_ && solution) {
1295        int numberColumns = model_->getNumCols();
1296        inputSolution_ = new double [numberColumns+1];
1297        memcpy(inputSolution_, solution, numberColumns*sizeof(double));
1298        inputSolution_[numberColumns] = objValue;
1299    }
[940]1300}
[2]1301
[912]1302//##############################################################################
1303
1304inline int compare3BranchingObjects(const CbcBranchingObject* br0,
[1286]1305                                    const CbcBranchingObject* br1)
[912]1306{
[1286]1307    const int t0 = br0->type();
1308    const int t1 = br1->type();
1309    if (t0 < t1) {
1310        return -1;
1311    }
1312    if (t0 > t1) {
1313        return 1;
1314    }
1315    return br0->compareOriginalObject(br1);
[912]1316}
1317
1318//==============================================================================
1319
1320inline bool compareBranchingObjects(const CbcBranchingObject* br0,
[1286]1321                                    const CbcBranchingObject* br1)
[912]1322{
[1286]1323    return compare3BranchingObjects(br0, br1) < 0;
[912]1324}
1325
1326//==============================================================================
1327
1328void
1329CbcHeuristicNode::gutsOfConstructor(CbcModel& model)
1330{
[1286]1331    //  CbcHeurDebugNodes(&model);
1332    CbcNode* node = model.currentNode();
1333    brObj_ = new CbcBranchingObject*[node->depth()];
1334    CbcNodeInfo* nodeInfo = node->nodeInfo();
1335    int cnt = 0;
1336    while (nodeInfo->parentBranch() != NULL) {
1337        const OsiBranchingObject* br = nodeInfo->parentBranch();
1338        const CbcBranchingObject* cbcbr = dynamic_cast<const CbcBranchingObject*>(br);
1339        if (! cbcbr) {
1340            throw CoinError("CbcHeuristicNode can be used only with CbcBranchingObjects.\n",
1341                            "gutsOfConstructor",
1342                            "CbcHeuristicNode",
1343                            __FILE__, __LINE__);
1344        }
1345        brObj_[cnt] = cbcbr->clone();
1346        brObj_[cnt]->previousBranch();
1347        ++cnt;
1348        nodeInfo = nodeInfo->parent();
[915]1349    }
[1286]1350    std::sort(brObj_, brObj_ + cnt, compareBranchingObjects);
1351    if (cnt <= 1) {
1352        numObjects_ = cnt;
1353    } else {
1354        numObjects_ = 0;
1355        CbcBranchingObject* br = NULL; // What should this be?
1356        for (int i = 1; i < cnt; ++i) {
1357            if (compare3BranchingObjects(brObj_[numObjects_], brObj_[i]) == 0) {
[1357]1358                int comp = brObj_[numObjects_]->compareBranchingObject(brObj_[i], br != 0);
[1286]1359                switch (comp) {
1360                case CbcRangeSame: // the same range
1361                case CbcRangeDisjoint: // disjoint decisions
1362                    // should not happen! we are on a chain!
1363                    abort();
1364                case CbcRangeSubset: // brObj_[numObjects_] is a subset of brObj_[i]
1365                    delete brObj_[i];
1366                    break;
1367                case CbcRangeSuperset: // brObj_[i] is a subset of brObj_[numObjects_]
1368                    delete brObj_[numObjects_];
1369                    brObj_[numObjects_] = brObj_[i];
1370                    break;
1371                case CbcRangeOverlap: // overlap
1372                    delete brObj_[i];
1373                    delete brObj_[numObjects_];
1374                    brObj_[numObjects_] = br;
1375                    break;
1376                }
1377                continue;
1378            } else {
1379                brObj_[++numObjects_] = brObj_[i];
1380            }
1381        }
1382        ++numObjects_;
[912]1383    }
1384}
1385
1386//==============================================================================
1387
1388CbcHeuristicNode::CbcHeuristicNode(CbcModel& model)
1389{
[1286]1390    gutsOfConstructor(model);
[912]1391}
1392
1393//==============================================================================
1394
1395double
[1286]1396CbcHeuristicNode::distance(const CbcHeuristicNode* node) const
[912]1397{
1398
[1286]1399    const double disjointWeight = 1;
1400    const double overlapWeight = 0.4;
1401    const double subsetWeight = 0.2;
1402    int countDisjointWeight = 0;
1403    int countOverlapWeight = 0;
1404    int countSubsetWeight = 0;
1405    int i = 0;
1406    int j = 0;
1407    double dist = 0.0;
[912]1408#ifdef PRINT_DEBUG
[1286]1409    printf(" numObjects_ = %i, node->numObjects_ = %i\n",
1410           numObjects_, node->numObjects_);
[912]1411#endif
[1286]1412    while ( i < numObjects_ && j < node->numObjects_) {
1413        CbcBranchingObject* br0 = brObj_[i];
1414        const CbcBranchingObject* br1 = node->brObj_[j];
[912]1415#ifdef PRINT_DEBUG
[1286]1416        const CbcIntegerBranchingObject* brPrint0 =
1417            dynamic_cast<const CbcIntegerBranchingObject*>(br0);
1418        const double* downBounds = brPrint0->downBounds();
1419        const double* upBounds = brPrint0->upBounds();
1420        int variable = brPrint0->variable();
1421        int way = brPrint0->way();
1422        printf("   br0: var %i downBd [%i,%i] upBd [%i,%i] way %i\n",
1423               variable, static_cast<int>(downBounds[0]), static_cast<int>(downBounds[1]),
1424               static_cast<int>(upBounds[0]), static_cast<int>(upBounds[1]), way);
1425        const CbcIntegerBranchingObject* brPrint1 =
1426            dynamic_cast<const CbcIntegerBranchingObject*>(br1);
1427        downBounds = brPrint1->downBounds();
1428        upBounds = brPrint1->upBounds();
1429        variable = brPrint1->variable();
1430        way = brPrint1->way();
1431        printf("   br1: var %i downBd [%i,%i] upBd [%i,%i] way %i\n",
1432               variable, static_cast<int>(downBounds[0]), static_cast<int>(downBounds[1]),
1433               static_cast<int>(upBounds[0]), static_cast<int>(upBounds[1]), way);
[912]1434#endif
[1286]1435        const int brComp = compare3BranchingObjects(br0, br1);
1436        if (brComp < 0) {
1437            dist += subsetWeight;
1438            countSubsetWeight++;
1439            ++i;
1440        } else if (brComp > 0) {
1441            dist += subsetWeight;
1442            countSubsetWeight++;
1443            ++j;
1444        } else {
1445            const int comp = br0->compareBranchingObject(br1, false);
1446            switch (comp) {
1447            case CbcRangeSame:
1448                // do nothing
1449                break;
1450            case CbcRangeDisjoint: // disjoint decisions
1451                dist += disjointWeight;
1452                countDisjointWeight++;
1453                break;
1454            case CbcRangeSubset: // subset one way or another
1455            case CbcRangeSuperset:
1456                dist += subsetWeight;
1457                countSubsetWeight++;
1458                break;
1459            case CbcRangeOverlap: // overlap
1460                dist += overlapWeight;
1461                countOverlapWeight++;
1462                break;
1463            }
1464            ++i;
1465            ++j;
1466        }
[912]1467    }
[1286]1468    dist += subsetWeight * (numObjects_ - i + node->numObjects_ - j);
1469    countSubsetWeight += (numObjects_ - i + node->numObjects_ - j);
1470    printf("subset = %i, overlap = %i, disjoint = %i\n", countSubsetWeight,
1471           countOverlapWeight, countDisjointWeight);
1472    return dist;
[912]1473}
1474
1475//==============================================================================
1476
1477CbcHeuristicNode::~CbcHeuristicNode()
1478{
[1286]1479    for (int i = 0; i < numObjects_; ++i) {
1480        delete brObj_[i];
1481    }
1482    delete [] brObj_;
[912]1483}
1484
1485//==============================================================================
1486
1487double
[915]1488CbcHeuristicNode::minDistance(const CbcHeuristicNodeList& nodeList) const
[912]1489{
[1286]1490    double minDist = COIN_DBL_MAX;
1491    for (int i = nodeList.size() - 1; i >= 0; --i) {
1492        minDist = CoinMin(minDist, distance(nodeList.node(i)));
1493    }
1494    return minDist;
[912]1495}
1496
1497//==============================================================================
1498
[915]1499bool
1500CbcHeuristicNode::minDistanceIsSmall(const CbcHeuristicNodeList& nodeList,
[1286]1501                                     const double threshold) const
[915]1502{
[1286]1503    for (int i = nodeList.size() - 1; i >= 0; --i) {
1504        if (distance(nodeList.node(i)) >= threshold) {
1505            continue;
1506        } else {
1507            return true;
1508        }
[915]1509    }
[1286]1510    return false;
[915]1511}
1512
1513//==============================================================================
1514
[912]1515double
[915]1516CbcHeuristicNode::avgDistance(const CbcHeuristicNodeList& nodeList) const
[912]1517{
[1286]1518    if (nodeList.size() == 0) {
1519        return COIN_DBL_MAX;
1520    }
1521    double sumDist = 0;
1522    for (int i = nodeList.size() - 1; i >= 0; --i) {
1523        sumDist += distance(nodeList.node(i));
1524    }
1525    return sumDist / nodeList.size();
[912]1526}
1527
1528//##############################################################################
1529
[2]1530// Default Constructor
[1286]1531CbcRounding::CbcRounding()
1532        : CbcHeuristic()
[2]1533{
[1286]1534    // matrix and row copy will automatically be empty
1535    seed_ = 7654321;
1536    down_ = NULL;
1537    up_ = NULL;
1538    equal_ = NULL;
[1315]1539    //whereFrom_ |= 16; // allow more often
[2]1540}
1541
1542// Constructor from model
1543CbcRounding::CbcRounding(CbcModel & model)
[1286]1544        : CbcHeuristic(model)
[2]1545{
[1286]1546    // Get a copy of original matrix (and by row for rounding);
1547    assert(model.solver());
1548    if (model.solver()->getNumRows()) {
1549        matrix_ = *model.solver()->getMatrixByCol();
1550        matrixByRow_ = *model.solver()->getMatrixByRow();
1551        validate();
1552    }
1553    down_ = NULL;
1554    up_ = NULL;
1555    equal_ = NULL;
1556    seed_ = 7654321;
[1315]1557    //whereFrom_ |= 16; // allow more often
[2]1558}
1559
[1286]1560// Destructor
[2]1561CbcRounding::~CbcRounding ()
1562{
[1286]1563    delete [] down_;
1564    delete [] up_;
1565    delete [] equal_;
[2]1566}
1567
1568// Clone
1569CbcHeuristic *
1570CbcRounding::clone() const
1571{
[1286]1572    return new CbcRounding(*this);
[2]1573}
[356]1574// Create C++ lines to get to current state
[1286]1575void
1576CbcRounding::generateCpp( FILE * fp)
[356]1577{
[1286]1578    CbcRounding other;
1579    fprintf(fp, "0#include \"CbcHeuristic.hpp\"\n");
1580    fprintf(fp, "3  CbcRounding rounding(*cbcModel);\n");
1581    CbcHeuristic::generateCpp(fp, "rounding");
1582    if (seed_ != other.seed_)
1583        fprintf(fp, "3  rounding.setSeed(%d);\n", seed_);
1584    else
1585        fprintf(fp, "4  rounding.setSeed(%d);\n", seed_);
1586    fprintf(fp, "3  cbcModel->addHeuristic(&rounding);\n");
[356]1587}
[838]1588//#define NEW_ROUNDING
[1286]1589// Copy constructor
[2]1590CbcRounding::CbcRounding(const CbcRounding & rhs)
[1286]1591        :
1592        CbcHeuristic(rhs),
1593        matrix_(rhs.matrix_),
1594        matrixByRow_(rhs.matrixByRow_),
1595        seed_(rhs.seed_)
[2]1596{
[838]1597#ifdef NEW_ROUNDING
[1286]1598    int numberColumns = matrix_.getNumCols();
1599    down_ = CoinCopyOfArray(rhs.down_, numberColumns);
1600    up_ = CoinCopyOfArray(rhs.up_, numberColumns);
1601    equal_ = CoinCopyOfArray(rhs.equal_, numberColumns);
[838]1602#else
[1286]1603    down_ = NULL;
1604    up_ = NULL;
1605    equal_ = NULL;
1606#endif
[2]1607}
[6]1608
[1286]1609// Assignment operator
1610CbcRounding &
1611CbcRounding::operator=( const CbcRounding & rhs)
[640]1612{
[1286]1613    if (this != &rhs) {
1614        CbcHeuristic::operator=(rhs);
1615        matrix_ = rhs.matrix_;
1616        matrixByRow_ = rhs.matrixByRow_;
[838]1617#ifdef NEW_ROUNDING
[1286]1618        delete [] down_;
1619        delete [] up_;
1620        delete [] equal_;
1621        int numberColumns = matrix_.getNumCols();
1622        down_ = CoinCopyOfArray(rhs.down_, numberColumns);
1623        up_ = CoinCopyOfArray(rhs.up_, numberColumns);
1624        equal_ = CoinCopyOfArray(rhs.equal_, numberColumns);
[838]1625#else
[1286]1626        down_ = NULL;
1627        up_ = NULL;
1628        equal_ = NULL;
1629#endif
1630        seed_ = rhs.seed_;
1631    }
1632    return *this;
[640]1633}
1634
[6]1635// Resets stuff if model changes
[1286]1636void
[6]1637CbcRounding::resetModel(CbcModel * model)
1638{
[1286]1639    model_ = model;
1640    // Get a copy of original matrix (and by row for rounding);
1641    assert(model_->solver());
1642    matrix_ = *model_->solver()->getMatrixByCol();
1643    matrixByRow_ = *model_->solver()->getMatrixByRow();
1644    validate();
[6]1645}
[2]1646// See if rounding will give solution
1647// Sets value of solution
1648// Assumes rhs for original matrix still okay
[1286]1649// At present only works with integers
[2]1650// Fix values if asked for
1651// Returns 1 if solution, 0 if not
1652int
1653CbcRounding::solution(double & solutionValue,
[1286]1654                      double * betterSolution)
[2]1655{
1656
[1286]1657    numCouldRun_++;
1658    // See if to do
1659    if (!when() || (when() % 10 == 1 && model_->phase() != 1) ||
1660            (when() % 10 == 2 && (model_->phase() != 2 && model_->phase() != 3)))
1661        return 0; // switched off
1662    numRuns_++;
1663    OsiSolverInterface * solver = model_->solver();
1664    double direction = solver->getObjSense();
1665    double newSolutionValue = direction * solver->getObjValue();
1666    return solution(solutionValue, betterSolution, newSolutionValue);
[838]1667}
1668// See if rounding will give solution
1669// Sets value of solution
1670// Assumes rhs for original matrix still okay
[1286]1671// At present only works with integers
[838]1672// Fix values if asked for
1673// Returns 1 if solution, 0 if not
1674int
1675CbcRounding::solution(double & solutionValue,
[1286]1676                      double * betterSolution,
1677                      double newSolutionValue)
[838]1678{
1679
[1286]1680    // See if to do
1681    if (!when() || (when() % 10 == 1 && model_->phase() != 1) ||
1682            (when() % 10 == 2 && (model_->phase() != 2 && model_->phase() != 3)))
1683        return 0; // switched off
1684    OsiSolverInterface * solver = model_->solver();
1685    const double * lower = solver->getColLower();
1686    const double * upper = solver->getColUpper();
1687    const double * rowLower = solver->getRowLower();
1688    const double * rowUpper = solver->getRowUpper();
1689    const double * solution = solver->getColSolution();
1690    const double * objective = solver->getObjCoefficients();
1691    double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
1692    double primalTolerance;
1693    solver->getDblParam(OsiPrimalTolerance, primalTolerance);
[2]1694
[1286]1695    int numberRows = matrix_.getNumRows();
1696    assert (numberRows <= solver->getNumRows());
1697    int numberIntegers = model_->numberIntegers();
1698    const int * integerVariable = model_->integerVariable();
1699    int i;
1700    double direction = solver->getObjSense();
1701    //double newSolutionValue = direction*solver->getObjValue();
1702    int returnCode = 0;
1703    // Column copy
1704    const double * element = matrix_.getElements();
1705    const int * row = matrix_.getIndices();
1706    const CoinBigIndex * columnStart = matrix_.getVectorStarts();
1707    const int * columnLength = matrix_.getVectorLengths();
1708    // Row copy
1709    const double * elementByRow = matrixByRow_.getElements();
1710    const int * column = matrixByRow_.getIndices();
1711    const CoinBigIndex * rowStart = matrixByRow_.getVectorStarts();
1712    const int * rowLength = matrixByRow_.getVectorLengths();
[2]1713
[1286]1714    // Get solution array for heuristic solution
1715    int numberColumns = solver->getNumCols();
1716    double * newSolution = new double [numberColumns];
1717    memcpy(newSolution, solution, numberColumns*sizeof(double));
[2]1718
[1286]1719    double * rowActivity = new double[numberRows];
1720    memset(rowActivity, 0, numberRows*sizeof(double));
1721    for (i = 0; i < numberColumns; i++) {
1722        int j;
1723        double value = newSolution[i];
1724        if (value < lower[i]) {
1725            value = lower[i];
1726            newSolution[i] = value;
1727        } else if (value > upper[i]) {
1728            value = upper[i];
1729            newSolution[i] = value;
1730        }
1731        if (value) {
1732            for (j = columnStart[i];
1733                    j < columnStart[i] + columnLength[i]; j++) {
1734                int iRow = row[j];
1735                rowActivity[iRow] += value * element[j];
1736            }
1737        }
[64]1738    }
[1286]1739    // check was feasible - if not adjust (cleaning may move)
1740    for (i = 0; i < numberRows; i++) {
1741        if (rowActivity[i] < rowLower[i]) {
1742            //assert (rowActivity[i]>rowLower[i]-1000.0*primalTolerance);
1743            rowActivity[i] = rowLower[i];
1744        } else if (rowActivity[i] > rowUpper[i]) {
1745            //assert (rowActivity[i]<rowUpper[i]+1000.0*primalTolerance);
1746            rowActivity[i] = rowUpper[i];
1747        }
[2]1748    }
[1286]1749    for (i = 0; i < numberIntegers; i++) {
1750        int iColumn = integerVariable[i];
1751        double value = newSolution[iColumn];
1752        if (fabs(floor(value + 0.5) - value) > integerTolerance) {
1753            double below = floor(value);
1754            double newValue = newSolution[iColumn];
1755            double cost = direction * objective[iColumn];
1756            double move;
1757            if (cost > 0.0) {
1758                // try up
1759                move = 1.0 - (value - below);
1760            } else if (cost < 0.0) {
1761                // try down
1762                move = below - value;
[91]1763            } else {
[1286]1764                // won't be able to move unless we can grab another variable
1765                double randomNumber = randomNumberGenerator_.randomDouble();
1766                // which way?
1767                if (randomNumber < 0.5)
1768                    move = below - value;
1769                else
1770                    move = 1.0 - (value - below);
[91]1771            }
[1286]1772            newValue += move;
1773            newSolution[iColumn] = newValue;
1774            newSolutionValue += move * cost;
1775            int j;
1776            for (j = columnStart[iColumn];
1777                    j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1778                int iRow = row[j];
1779                rowActivity[iRow] += move * element[j];
[91]1780            }
1781        }
[1286]1782    }
1783
1784    double penalty = 0.0;
1785    const char * integerType = model_->integerType();
1786    // see if feasible - just using singletons
1787    for (i = 0; i < numberRows; i++) {
1788        double value = rowActivity[i];
1789        double thisInfeasibility = 0.0;
1790        if (value < rowLower[i] - primalTolerance)
1791            thisInfeasibility = value - rowLower[i];
1792        else if (value > rowUpper[i] + primalTolerance)
1793            thisInfeasibility = value - rowUpper[i];
1794        if (thisInfeasibility) {
1795            // See if there are any slacks I can use to fix up
1796            // maybe put in coding for multiple slacks?
1797            double bestCost = 1.0e50;
1798            int k;
1799            int iBest = -1;
1800            double addCost = 0.0;
1801            double newValue = 0.0;
1802            double changeRowActivity = 0.0;
1803            double absInfeasibility = fabs(thisInfeasibility);
1804            for (k = rowStart[i]; k < rowStart[i] + rowLength[i]; k++) {
1805                int iColumn = column[k];
1806                // See if all elements help
1807                if (columnLength[iColumn] == 1) {
1808                    double currentValue = newSolution[iColumn];
1809                    double elementValue = elementByRow[k];
1810                    double lowerValue = lower[iColumn];
1811                    double upperValue = upper[iColumn];
1812                    double gap = rowUpper[i] - rowLower[i];
1813                    double absElement = fabs(elementValue);
1814                    if (thisInfeasibility*elementValue > 0.0) {
1815                        // we want to reduce
1816                        if ((currentValue - lowerValue)*absElement >= absInfeasibility) {
1817                            // possible - check if integer
1818                            double distance = absInfeasibility / absElement;
1819                            double thisCost = -direction * objective[iColumn] * distance;
1820                            if (integerType[iColumn]) {
1821                                distance = ceil(distance - primalTolerance);
1822                                if (currentValue - distance >= lowerValue - primalTolerance) {
1823                                    if (absInfeasibility - distance*absElement < -gap - primalTolerance)
1824                                        thisCost = 1.0e100; // no good
1825                                    else
1826                                        thisCost = -direction * objective[iColumn] * distance;
1827                                } else {
1828                                    thisCost = 1.0e100; // no good
1829                                }
1830                            }
1831                            if (thisCost < bestCost) {
1832                                bestCost = thisCost;
1833                                iBest = iColumn;
1834                                addCost = thisCost;
1835                                newValue = currentValue - distance;
1836                                changeRowActivity = -distance * elementValue;
1837                            }
1838                        }
1839                    } else {
1840                        // we want to increase
1841                        if ((upperValue - currentValue)*absElement >= absInfeasibility) {
1842                            // possible - check if integer
1843                            double distance = absInfeasibility / absElement;
1844                            double thisCost = direction * objective[iColumn] * distance;
1845                            if (integerType[iColumn]) {
1846                                distance = ceil(distance - 1.0e-7);
1847                                assert (currentValue - distance <= upperValue + primalTolerance);
1848                                if (absInfeasibility - distance*absElement < -gap - primalTolerance)
1849                                    thisCost = 1.0e100; // no good
1850                                else
1851                                    thisCost = direction * objective[iColumn] * distance;
1852                            }
1853                            if (thisCost < bestCost) {
1854                                bestCost = thisCost;
1855                                iBest = iColumn;
1856                                addCost = thisCost;
1857                                newValue = currentValue + distance;
1858                                changeRowActivity = distance * elementValue;
1859                            }
1860                        }
1861                    }
1862                }
[91]1863            }
[1286]1864            if (iBest >= 0) {
1865                /*printf("Infeasibility of %g on row %d cost %g\n",
1866                  thisInfeasibility,i,addCost);*/
1867                newSolution[iBest] = newValue;
1868                thisInfeasibility = 0.0;
1869                newSolutionValue += addCost;
1870                rowActivity[i] += changeRowActivity;
1871            }
1872            penalty += fabs(thisInfeasibility);
[91]1873        }
1874    }
[1286]1875    if (penalty) {
1876        // see if feasible using any
1877        // first continuous
1878        double penaltyChange = 0.0;
1879        int iColumn;
1880        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
1881            if (integerType[iColumn])
1882                continue;
1883            double currentValue = newSolution[iColumn];
1884            double lowerValue = lower[iColumn];
1885            double upperValue = upper[iColumn];
1886            int j;
1887            int anyBadDown = 0;
1888            int anyBadUp = 0;
1889            double upImprovement = 0.0;
1890            double downImprovement = 0.0;
1891            for (j = columnStart[iColumn];
1892                    j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1893                int iRow = row[j];
1894                if (rowUpper[iRow] > rowLower[iRow]) {
1895                    double value = element[j];
1896                    if (rowActivity[iRow] > rowUpper[iRow] + primalTolerance) {
1897                        // infeasible above
1898                        downImprovement += value;
1899                        upImprovement -= value;
1900                        if (value > 0.0)
1901                            anyBadUp++;
1902                        else
1903                            anyBadDown++;
1904                    } else if (rowActivity[iRow] > rowUpper[iRow] - primalTolerance) {
1905                        // feasible at ub
1906                        if (value > 0.0) {
1907                            upImprovement -= value;
1908                            anyBadUp++;
1909                        } else {
1910                            downImprovement += value;
1911                            anyBadDown++;
1912                        }
1913                    } else if (rowActivity[iRow] > rowLower[iRow] + primalTolerance) {
1914                        // feasible in interior
1915                    } else if (rowActivity[iRow] > rowLower[iRow] - primalTolerance) {
1916                        // feasible at lb
1917                        if (value < 0.0) {
1918                            upImprovement += value;
1919                            anyBadUp++;
1920                        } else {
1921                            downImprovement -= value;
1922                            anyBadDown++;
1923                        }
1924                    } else {
1925                        // infeasible below
1926                        downImprovement -= value;
1927                        upImprovement += value;
1928                        if (value < 0.0)
1929                            anyBadUp++;
1930                        else
1931                            anyBadDown++;
1932                    }
1933                } else {
1934                    // equality row
1935                    double value = element[j];
1936                    if (rowActivity[iRow] > rowUpper[iRow] + primalTolerance) {
1937                        // infeasible above
1938                        downImprovement += value;
1939                        upImprovement -= value;
1940                        if (value > 0.0)
1941                            anyBadUp++;
1942                        else
1943                            anyBadDown++;
1944                    } else if (rowActivity[iRow] < rowLower[iRow] - primalTolerance) {
1945                        // infeasible below
1946                        downImprovement -= value;
1947                        upImprovement += value;
1948                        if (value < 0.0)
1949                            anyBadUp++;
1950                        else
1951                            anyBadDown++;
1952                    } else {
1953                        // feasible - no good
1954                        anyBadUp = -1;
1955                        anyBadDown = -1;
1956                        break;
1957                    }
1958                }
[91]1959            }
[1286]1960            // could change tests for anyBad
1961            if (anyBadUp)
1962                upImprovement = 0.0;
1963            if (anyBadDown)
1964                downImprovement = 0.0;
1965            double way = 0.0;
1966            double improvement = 0.0;
1967            if (downImprovement > 0.0 && currentValue > lowerValue) {
1968                way = -1.0;
1969                improvement = downImprovement;
1970            } else if (upImprovement > 0.0 && currentValue < upperValue) {
1971                way = 1.0;
1972                improvement = upImprovement;
[91]1973            }
[1286]1974            if (way) {
1975                // can improve
1976                double distance;
1977                if (way > 0.0)
1978                    distance = upperValue - currentValue;
1979                else
1980                    distance = currentValue - lowerValue;
1981                for (j = columnStart[iColumn];
1982                        j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1983                    int iRow = row[j];
1984                    double value = element[j] * way;
1985                    if (rowActivity[iRow] > rowUpper[iRow] + primalTolerance) {
1986                        // infeasible above
1987                        assert (value < 0.0);
1988                        double gap = rowActivity[iRow] - rowUpper[iRow];
1989                        if (gap + value*distance < 0.0)
1990                            distance = -gap / value;
1991                    } else if (rowActivity[iRow] < rowLower[iRow] - primalTolerance) {
1992                        // infeasible below
1993                        assert (value > 0.0);
1994                        double gap = rowActivity[iRow] - rowLower[iRow];
1995                        if (gap + value*distance > 0.0)
1996                            distance = -gap / value;
1997                    } else {
1998                        // feasible
1999                        if (value > 0) {
2000                            double gap = rowActivity[iRow] - rowUpper[iRow];
2001                            if (gap + value*distance > 0.0)
2002                                distance = -gap / value;
2003                        } else {
2004                            double gap = rowActivity[iRow] - rowLower[iRow];
2005                            if (gap + value*distance < 0.0)
2006                                distance = -gap / value;
2007                        }
2008                    }
2009                }
2010                //move
2011                penaltyChange += improvement * distance;
2012                distance *= way;
2013                newSolution[iColumn] += distance;
2014                newSolutionValue += direction * objective[iColumn] * distance;
2015                for (j = columnStart[iColumn];
2016                        j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2017                    int iRow = row[j];
2018                    double value = element[j];
2019                    rowActivity[iRow] += distance * value;
2020                }
[91]2021            }
2022        }
[1286]2023        // and now all if improving
2024        double lastChange = penaltyChange ? 1.0 : 0.0;
2025        while (lastChange > 1.0e-2) {
2026            lastChange = 0;
2027            for (iColumn = 0; iColumn < numberColumns; iColumn++) {
2028                bool isInteger = (integerType[iColumn] != 0);
2029                double currentValue = newSolution[iColumn];
2030                double lowerValue = lower[iColumn];
2031                double upperValue = upper[iColumn];
2032                int j;
2033                int anyBadDown = 0;
2034                int anyBadUp = 0;
2035                double upImprovement = 0.0;
2036                double downImprovement = 0.0;
2037                for (j = columnStart[iColumn];
2038                        j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2039                    int iRow = row[j];
2040                    double value = element[j];
2041                    if (isInteger) {
2042                        if (value > 0.0) {
2043                            if (rowActivity[iRow] + value > rowUpper[iRow] + primalTolerance)
2044                                anyBadUp++;
2045                            if (rowActivity[iRow] - value < rowLower[iRow] - primalTolerance)
2046                                anyBadDown++;
2047                        } else {
2048                            if (rowActivity[iRow] - value > rowUpper[iRow] + primalTolerance)
2049                                anyBadDown++;
2050                            if (rowActivity[iRow] + value < rowLower[iRow] - primalTolerance)
2051                                anyBadUp++;
2052                        }
2053                    }
2054                    if (rowUpper[iRow] > rowLower[iRow]) {
2055                        if (rowActivity[iRow] > rowUpper[iRow] + primalTolerance) {
2056                            // infeasible above
2057                            downImprovement += value;
2058                            upImprovement -= value;
2059                            if (value > 0.0)
2060                                anyBadUp++;
2061                            else
2062                                anyBadDown++;
2063                        } else if (rowActivity[iRow] > rowUpper[iRow] - primalTolerance) {
2064                            // feasible at ub
2065                            if (value > 0.0) {
2066                                upImprovement -= value;
2067                                anyBadUp++;
2068                            } else {
2069                                downImprovement += value;
2070                                anyBadDown++;
2071                            }
2072                        } else if (rowActivity[iRow] > rowLower[iRow] + primalTolerance) {
2073                            // feasible in interior
2074                        } else if (rowActivity[iRow] > rowLower[iRow] - primalTolerance) {
2075                            // feasible at lb
2076                            if (value < 0.0) {
2077                                upImprovement += value;
2078                                anyBadUp++;
2079                            } else {
2080                                downImprovement -= value;
2081                                anyBadDown++;
2082                            }
2083                        } else {
2084                            // infeasible below
2085                            downImprovement -= value;
2086                            upImprovement += value;
2087                            if (value < 0.0)
2088                                anyBadUp++;
2089                            else
2090                                anyBadDown++;
2091                        }
2092                    } else {
2093                        // equality row
2094                        if (rowActivity[iRow] > rowUpper[iRow] + primalTolerance) {
2095                            // infeasible above
2096                            downImprovement += value;
2097                            upImprovement -= value;
2098                            if (value > 0.0)
2099                                anyBadUp++;
2100                            else
2101                                anyBadDown++;
2102                        } else if (rowActivity[iRow] < rowLower[iRow] - primalTolerance) {
2103                            // infeasible below
2104                            downImprovement -= value;
2105                            upImprovement += value;
2106                            if (value < 0.0)
2107                                anyBadUp++;
2108                            else
2109                                anyBadDown++;
2110                        } else {
2111                            // feasible - no good
2112                            anyBadUp = -1;
2113                            anyBadDown = -1;
2114                            break;
2115                        }
2116                    }
2117                }
2118                // could change tests for anyBad
2119                if (anyBadUp)
2120                    upImprovement = 0.0;
2121                if (anyBadDown)
2122                    downImprovement = 0.0;
2123                double way = 0.0;
2124                double improvement = 0.0;
2125                if (downImprovement > 0.0 && currentValue > lowerValue) {
2126                    way = -1.0;
2127                    improvement = downImprovement;
2128                } else if (upImprovement > 0.0 && currentValue < upperValue) {
2129                    way = 1.0;
2130                    improvement = upImprovement;
2131                }
2132                if (way) {
2133                    // can improve
2134                    double distance = COIN_DBL_MAX;
2135                    for (j = columnStart[iColumn];
2136                            j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2137                        int iRow = row[j];
2138                        double value = element[j] * way;
2139                        if (rowActivity[iRow] > rowUpper[iRow] + primalTolerance) {
2140                            // infeasible above
2141                            assert (value < 0.0);
2142                            double gap = rowActivity[iRow] - rowUpper[iRow];
2143                            if (gap + value*distance < 0.0) {
2144                                // If integer then has to move by 1
2145                                if (!isInteger)
2146                                    distance = -gap / value;
2147                                else
2148                                    distance = CoinMax(-gap / value, 1.0);
2149                            }
2150                        } else if (rowActivity[iRow] < rowLower[iRow] - primalTolerance) {
2151                            // infeasible below
2152                            assert (value > 0.0);
2153                            double gap = rowActivity[iRow] - rowLower[iRow];
2154                            if (gap + value*distance > 0.0) {
2155                                // If integer then has to move by 1
2156                                if (!isInteger)
2157                                    distance = -gap / value;
2158                                else
2159                                    distance = CoinMax(-gap / value, 1.0);
2160                            }
2161                        } else {
2162                            // feasible
2163                            if (value > 0) {
2164                                double gap = rowActivity[iRow] - rowUpper[iRow];
2165                                if (gap + value*distance > 0.0)
2166                                    distance = -gap / value;
2167                            } else {
2168                                double gap = rowActivity[iRow] - rowLower[iRow];
2169                                if (gap + value*distance < 0.0)
2170                                    distance = -gap / value;
2171                            }
2172                        }
2173                    }
2174                    if (isInteger)
2175                        distance = floor(distance + 1.05e-8);
2176                    if (!distance) {
2177                        // should never happen
2178                        //printf("zero distance in CbcRounding - debug\n");
2179                    }
2180                    //move
2181                    lastChange += improvement * distance;
2182                    distance *= way;
2183                    newSolution[iColumn] += distance;
2184                    newSolutionValue += direction * objective[iColumn] * distance;
2185                    for (j = columnStart[iColumn];
2186                            j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2187                        int iRow = row[j];
2188                        double value = element[j];
2189                        rowActivity[iRow] += distance * value;
2190                    }
2191                }
2192            }
2193            penaltyChange += lastChange;
[91]2194        }
[1286]2195        penalty -= penaltyChange;
2196        if (penalty < 1.0e-5*fabs(penaltyChange)) {
2197            // recompute
2198            penalty = 0.0;
2199            for (i = 0; i < numberRows; i++) {
2200                double value = rowActivity[i];
2201                if (value < rowLower[i] - primalTolerance)
2202                    penalty += rowLower[i] - value;
2203                else if (value > rowUpper[i] + primalTolerance)
2204                    penalty += value - rowUpper[i];
[91]2205            }
2206        }
2207    }
[2]2208
[1286]2209    // Could also set SOS (using random) and repeat
2210    if (!penalty) {
2211        // See if we can do better
2212        //seed_++;
2213        //CoinSeedRandom(seed_);
2214        // Random number between 0 and 1.
2215        double randomNumber = randomNumberGenerator_.randomDouble();
2216        int iPass;
2217        int start[2];
2218        int end[2];
2219        int iRandom = static_cast<int> (randomNumber * (static_cast<double> (numberIntegers)));
2220        start[0] = iRandom;
2221        end[0] = numberIntegers;
2222        start[1] = 0;
2223        end[1] = iRandom;
2224        for (iPass = 0; iPass < 2; iPass++) {
2225            int i;
2226            for (i = start[iPass]; i < end[iPass]; i++) {
2227                int iColumn = integerVariable[i];
[201]2228#ifndef NDEBUG
[1286]2229                double value = newSolution[iColumn];
2230                assert (fabs(floor(value + 0.5) - value) < integerTolerance);
[201]2231#endif
[1286]2232                double cost = direction * objective[iColumn];
2233                double move = 0.0;
2234                if (cost > 0.0)
2235                    move = -1.0;
2236                else if (cost < 0.0)
2237                    move = 1.0;
2238                while (move) {
2239                    bool good = true;
2240                    double newValue = newSolution[iColumn] + move;
2241                    if (newValue < lower[iColumn] - primalTolerance ||
2242                            newValue > upper[iColumn] + primalTolerance) {
2243                        move = 0.0;
2244                    } else {
2245                        // see if we can move
2246                        int j;
2247                        for (j = columnStart[iColumn];
2248                                j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2249                            int iRow = row[j];
2250                            double newActivity = rowActivity[iRow] + move * element[j];
2251                            if (newActivity < rowLower[iRow] - primalTolerance ||
2252                                    newActivity > rowUpper[iRow] + primalTolerance) {
2253                                good = false;
2254                                break;
2255                            }
2256                        }
2257                        if (good) {
2258                            newSolution[iColumn] = newValue;
2259                            newSolutionValue += move * cost;
2260                            int j;
2261                            for (j = columnStart[iColumn];
2262                                    j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2263                                int iRow = row[j];
2264                                rowActivity[iRow] += move * element[j];
2265                            }
2266                        } else {
2267                            move = 0.0;
2268                        }
2269                    }
2270                }
2271            }
2272        }
2273        // Just in case of some stupidity
2274        double objOffset = 0.0;
2275        solver->getDblParam(OsiObjOffset, objOffset);
2276        newSolutionValue = -objOffset;
2277        for ( i = 0 ; i < numberColumns ; i++ )
2278            newSolutionValue += objective[i] * newSolution[i];
2279        newSolutionValue *= direction;
2280        //printf("new solution value %g %g\n",newSolutionValue,solutionValue);
2281        if (newSolutionValue < solutionValue) {
2282            // paranoid check
2283            memset(rowActivity, 0, numberRows*sizeof(double));
2284            for (i = 0; i < numberColumns; i++) {
2285                int j;
2286                double value = newSolution[i];
2287                if (value) {
2288                    for (j = columnStart[i];
2289                            j < columnStart[i] + columnLength[i]; j++) {
2290                        int iRow = row[j];
2291                        rowActivity[iRow] += value * element[j];
2292                    }
2293                }
2294            }
2295            // check was approximately feasible
2296            bool feasible = true;
2297            for (i = 0; i < numberRows; i++) {
2298                if (rowActivity[i] < rowLower[i]) {
2299                    if (rowActivity[i] < rowLower[i] - 1000.0*primalTolerance)
2300                        feasible = false;
2301                } else if (rowActivity[i] > rowUpper[i]) {
2302                    if (rowActivity[i] > rowUpper[i] + 1000.0*primalTolerance)
2303                        feasible = false;
2304                }
2305            }
2306            if (feasible) {
2307                // new solution
2308                memcpy(betterSolution, newSolution, numberColumns*sizeof(double));
2309                solutionValue = newSolutionValue;
2310                //printf("** Solution of %g found by rounding\n",newSolutionValue);
2311                returnCode = 1;
2312            } else {
2313                // Can easily happen
2314                //printf("Debug CbcRounding giving bad solution\n");
2315            }
2316        }
[2]2317    }
[838]2318#ifdef NEW_ROUNDING
[1286]2319    if (!returnCode) {
[1393]2320#ifdef JJF_ZERO
[1286]2321        // back to starting point
2322        memcpy(newSolution, solution, numberColumns*sizeof(double));
2323        memset(rowActivity, 0, numberRows*sizeof(double));
2324        for (i = 0; i < numberColumns; i++) {
2325            int j;
2326            double value = newSolution[i];
2327            if (value < lower[i]) {
2328                value = lower[i];
2329                newSolution[i] = value;
2330            } else if (value > upper[i]) {
2331                value = upper[i];
2332                newSolution[i] = value;
2333            }
2334            if (value) {
2335                for (j = columnStart[i];
2336                        j < columnStart[i] + columnLength[i]; j++) {
2337                    int iRow = row[j];
2338                    rowActivity[iRow] += value * element[j];
2339                }
2340            }
2341        }
2342        // check was feasible - if not adjust (cleaning may move)
2343        for (i = 0; i < numberRows; i++) {
2344            if (rowActivity[i] < rowLower[i]) {
2345                //assert (rowActivity[i]>rowLower[i]-1000.0*primalTolerance);
2346                rowActivity[i] = rowLower[i];
2347            } else if (rowActivity[i] > rowUpper[i]) {
2348                //assert (rowActivity[i]<rowUpper[i]+1000.0*primalTolerance);
2349                rowActivity[i] = rowUpper[i];
2350            }
2351        }
[838]2352#endif
[1286]2353        int * candidate = new int [numberColumns];
2354        int nCandidate = 0;
2355        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
2356            bool isInteger = (integerType[iColumn] != 0);
2357            if (isInteger) {
2358                double currentValue = newSolution[iColumn];
2359                if (fabs(currentValue - floor(currentValue + 0.5)) > 1.0e-8)
2360                    candidate[nCandidate++] = iColumn;
2361            }
2362        }
2363        if (true) {
2364            // Rounding as in Berthold
2365            while (nCandidate) {
2366                double infeasibility = 1.0e-7;
2367                int iRow = -1;
2368                for (i = 0; i < numberRows; i++) {
2369                    double value = 0.0;
2370                    if (rowActivity[i] < rowLower[i]) {
2371                        value = rowLower[i] - rowActivity[i];
2372                    } else if (rowActivity[i] > rowUpper[i]) {
2373                        value = rowActivity[i] - rowUpper[i];
2374                    }
2375                    if (value > infeasibility) {
2376                        infeasibility = value;
2377                        iRow = i;
2378                    }
2379                }
2380                if (iRow >= 0) {
2381                    // infeasible
2382                } else {
2383                    // feasible
2384                }
2385            }
2386        } else {
2387            // Shifting as in Berthold
2388        }
2389        delete [] candidate;
[838]2390    }
2391#endif
[1286]2392    delete [] newSolution;
2393    delete [] rowActivity;
2394    return returnCode;
[2]2395}
2396// update model
2397void CbcRounding::setModel(CbcModel * model)
2398{
[1286]2399    model_ = model;
2400    // Get a copy of original matrix (and by row for rounding);
2401    assert(model_->solver());
2402    if (model_->solver()->getNumRows()) {
2403        matrix_ = *model_->solver()->getMatrixByCol();
2404        matrixByRow_ = *model_->solver()->getMatrixByRow();
2405        // make sure model okay for heuristic
2406        validate();
2407    }
[2]2408}
[72]2409// Validate model i.e. sets when_ to 0 if necessary (may be NULL)
[1286]2410void
2411CbcRounding::validate()
[72]2412{
[1286]2413    if (model_ && (when() % 100) < 10) {
2414        if (model_->numberIntegers() !=
2415                model_->numberObjects() && (model_->numberObjects() ||
2416                                            (model_->specialOptions()&1024) == 0)) {
2417            int numberOdd = 0;
2418            for (int i = 0; i < model_->numberObjects(); i++) {
2419                if (!model_->object(i)->canDoHeuristics())
2420                    numberOdd++;
2421            }
2422            if (numberOdd)
2423                setWhen(0);
2424        }
[1271]2425    }
[838]2426#ifdef NEW_ROUNDING
[1286]2427    int numberColumns = matrix_.getNumCols();
2428    down_ = new unsigned short [numberColumns];
2429    up_ = new unsigned short [numberColumns];
2430    equal_ = new unsigned short [numberColumns];
2431    // Column copy
2432    const double * element = matrix_.getElements();
2433    const int * row = matrix_.getIndices();
2434    const CoinBigIndex * columnStart = matrix_.getVectorStarts();
2435    const int * columnLength = matrix_.getVectorLengths();
2436    const double * rowLower = model.solver()->getRowLower();
2437    const double * rowUpper = model.solver()->getRowUpper();
2438    for (int i = 0; i < numberColumns; i++) {
2439        int down = 0;
2440        int up = 0;
2441        int equal = 0;
2442        if (columnLength[i] > 65535) {
2443            equal[0] = 65535;
2444            break; // unlikely to work
2445        }
2446        for (CoinBigIndex j = columnStart[i];
2447                j < columnStart[i] + columnLength[i]; j++) {
2448            int iRow = row[j];
2449            if (rowLower[iRow] > -1.0e20 && rowUpper[iRow] < 1.0e20) {
2450                equal++;
2451            } else if (element[j] > 0.0) {
2452                if (rowUpper[iRow] < 1.0e20)
2453                    up++;
2454                else
2455                    down--;
2456            } else {
2457                if (rowLower[iRow] > -1.0e20)
2458                    up++;
2459                else
2460                    down--;
2461            }
2462        }
2463        down_[i] = (unsigned short) down;
2464        up_[i] = (unsigned short) up;
2465        equal_[i] = (unsigned short) equal;
[838]2466    }
2467#else
[1286]2468    down_ = NULL;
2469    up_ = NULL;
2470    equal_ = NULL;
2471#endif
[72]2472}
[2]2473
[264]2474// Default Constructor
[1286]2475CbcHeuristicPartial::CbcHeuristicPartial()
2476        : CbcHeuristic()
[854]2477{
[1286]2478    fixPriority_ = 10000;
[854]2479}
2480
2481// Constructor from model
2482CbcHeuristicPartial::CbcHeuristicPartial(CbcModel & model, int fixPriority, int numberNodes)
[1286]2483        : CbcHeuristic(model)
[854]2484{
[1286]2485    fixPriority_ = fixPriority;
2486    setNumberNodes(numberNodes);
2487    validate();
[854]2488}
2489
[1286]2490// Destructor
[854]2491CbcHeuristicPartial::~CbcHeuristicPartial ()
2492{
2493}
2494
2495// Clone
2496CbcHeuristic *
2497CbcHeuristicPartial::clone() const
2498{
[1286]2499    return new CbcHeuristicPartial(*this);
[854]2500}
2501// Create C++ lines to get to current state
[1286]2502void
2503CbcHeuristicPartial::generateCpp( FILE * fp)
[854]2504{
[1286]2505    CbcHeuristicPartial other;
2506    fprintf(fp, "0#include \"CbcHeuristic.hpp\"\n");
2507    fprintf(fp, "3  CbcHeuristicPartial partial(*cbcModel);\n");
2508    CbcHeuristic::generateCpp(fp, "partial");
2509    if (fixPriority_ != other.fixPriority_)
2510        fprintf(fp, "3  partial.setFixPriority(%d);\n", fixPriority_);
2511    else
2512        fprintf(fp, "4  partial.setFixPriority(%d);\n", fixPriority_);
2513    fprintf(fp, "3  cbcModel->addHeuristic(&partial);\n");
[854]2514}
2515//#define NEW_PARTIAL
[1286]2516// Copy constructor
[854]2517CbcHeuristicPartial::CbcHeuristicPartial(const CbcHeuristicPartial & rhs)
[1286]2518        :
2519        CbcHeuristic(rhs),
2520        fixPriority_(rhs.fixPriority_)
[854]2521{
2522}
2523
[1286]2524// Assignment operator
2525CbcHeuristicPartial &
2526CbcHeuristicPartial::operator=( const CbcHeuristicPartial & rhs)
[854]2527{
[1286]2528    if (this != &rhs) {
2529        CbcHeuristic::operator=(rhs);
2530        fixPriority_ = rhs.fixPriority_;
2531    }
2532    return *this;
[854]2533}
2534
2535// Resets stuff if model changes
[1286]2536void
[854]2537CbcHeuristicPartial::resetModel(CbcModel * model)
2538{
[1286]2539    model_ = model;
2540    // Get a copy of original matrix (and by row for partial);
2541    assert(model_->solver());
2542    validate();
[854]2543}
2544// See if partial will give solution
2545// Sets value of solution
2546// Assumes rhs for original matrix still okay
[1286]2547// At present only works with integers
[854]2548// Fix values if asked for
2549// Returns 1 if solution, 0 if not
2550int
2551CbcHeuristicPartial::solution(double & solutionValue,
[1286]2552                              double * betterSolution)
[854]2553{
[1286]2554    // Return if already done
2555    if (fixPriority_ < 0)
2556        return 0; // switched off
2557    const double * hotstartSolution = model_->hotstartSolution();
2558    const int * hotstartPriorities = model_->hotstartPriorities();
2559    if (!hotstartSolution)
2560        return 0;
2561    OsiSolverInterface * solver = model_->solver();
[854]2562
[1286]2563    int numberIntegers = model_->numberIntegers();
2564    const int * integerVariable = model_->integerVariable();
[854]2565
[1286]2566    OsiSolverInterface * newSolver = model_->continuousSolver()->clone();
2567    const double * colLower = newSolver->getColLower();
2568    const double * colUpper = newSolver->getColUpper();
2569
2570    double primalTolerance;
2571    solver->getDblParam(OsiPrimalTolerance, primalTolerance);
2572
2573    int i;
2574    int numberFixed = 0;
2575    int returnCode = 0;
2576
2577    for (i = 0; i < numberIntegers; i++) {
2578        int iColumn = integerVariable[i];
2579        if (abs(hotstartPriorities[iColumn]) <= fixPriority_) {
2580            double value = hotstartSolution[iColumn];
2581            double lower = colLower[iColumn];
2582            double upper = colUpper[iColumn];
2583            value = CoinMax(value, lower);
2584            value = CoinMin(value, upper);
2585            if (fabs(value - floor(value + 0.5)) < 1.0e-8) {
2586                value = floor(value + 0.5);
2587                newSolver->setColLower(iColumn, value);
2588                newSolver->setColUpper(iColumn, value);
2589                numberFixed++;
2590            }
2591        }
[854]2592    }
[1286]2593    if (numberFixed > numberIntegers / 5 - 100000000) {
[854]2594#ifdef COIN_DEVELOP
[1286]2595        printf("%d integers fixed\n", numberFixed);
[854]2596#endif
[1286]2597        returnCode = smallBranchAndBound(newSolver, numberNodes_, betterSolution, solutionValue,
2598                                         model_->getCutoff(), "CbcHeuristicPartial");
2599        if (returnCode < 0)
2600            returnCode = 0; // returned on size
2601        //printf("return code %d",returnCode);
2602        if ((returnCode&2) != 0) {
2603            // could add cut
2604            returnCode &= ~2;
2605            //printf("could add cut with %d elements (if all 0-1)\n",nFix);
2606        } else {
2607            //printf("\n");
2608        }
[854]2609    }
[1286]2610    fixPriority_ = -1; // switch off
2611
2612    delete newSolver;
2613    return returnCode;
[854]2614}
2615// update model
2616void CbcHeuristicPartial::setModel(CbcModel * model)
2617{
[1286]2618    model_ = model;
2619    assert(model_->solver());
2620    // make sure model okay for heuristic
2621    validate();
[854]2622}
2623// Validate model i.e. sets when_ to 0 if necessary (may be NULL)
[1286]2624void
2625CbcHeuristicPartial::validate()
[854]2626{
[1286]2627    if (model_ && (when() % 100) < 10) {
2628        if (model_->numberIntegers() !=
2629                model_->numberObjects())
2630            setWhen(0);
2631    }
[854]2632}
[1013]2633bool
[1271]2634CbcHeuristicPartial::shouldHeurRun(int /*whereFrom*/)
[1013]2635{
[1286]2636    return true;
[1013]2637}
[854]2638
2639// Default Constructor
[1286]2640CbcSerendipity::CbcSerendipity()
2641        : CbcHeuristic()
[264]2642{
2643}
2644
2645// Constructor from model
2646CbcSerendipity::CbcSerendipity(CbcModel & model)
[1286]2647        : CbcHeuristic(model)
[264]2648{
2649}
2650
[1286]2651// Destructor
[264]2652CbcSerendipity::~CbcSerendipity ()
2653{
2654}
2655
2656// Clone
2657CbcHeuristic *
2658CbcSerendipity::clone() const
2659{
[1286]2660    return new CbcSerendipity(*this);
[264]2661}
[356]2662// Create C++ lines to get to current state
[1286]2663void
2664CbcSerendipity::generateCpp( FILE * fp)
[356]2665{
[1286]2666    fprintf(fp, "0#include \"CbcHeuristic.hpp\"\n");
2667    fprintf(fp, "3  CbcSerendipity serendipity(*cbcModel);\n");
2668    CbcHeuristic::generateCpp(fp, "serendipity");
2669    fprintf(fp, "3  cbcModel->addHeuristic(&serendipity);\n");
[356]2670}
[264]2671
[1286]2672// Copy constructor
[264]2673CbcSerendipity::CbcSerendipity(const CbcSerendipity & rhs)
[1286]2674        :
2675        CbcHeuristic(rhs)
[264]2676{
2677}
2678
[1286]2679// Assignment operator
2680CbcSerendipity &
2681CbcSerendipity::operator=( const CbcSerendipity & rhs)
[640]2682{
[1286]2683    if (this != &rhs) {
2684        CbcHeuristic::operator=(rhs);
2685    }
2686    return *this;
[640]2687}
2688
[264]2689// Returns 1 if solution, 0 if not
2690int
2691CbcSerendipity::solution(double & solutionValue,
[1286]2692                         double * betterSolution)
[264]2693{
[1286]2694    if (!model_)
2695        return 0;
2696    if (!inputSolution_) {
2697        // get information on solver type
2698        OsiAuxInfo * auxInfo = model_->solver()->getAuxiliaryInfo();
2699        OsiBabSolver * auxiliaryInfo = dynamic_cast< OsiBabSolver *> (auxInfo);
2700        if (auxiliaryInfo) {
2701            return auxiliaryInfo->solution(solutionValue, betterSolution, model_->solver()->getNumCols());
2702        } else {
2703            return 0;
2704        }
[940]2705    } else {
[1286]2706        int numberColumns = model_->getNumCols();
2707        double value = inputSolution_[numberColumns];
2708        int returnCode = 0;
2709        if (value < solutionValue) {
2710            solutionValue = value;
2711            memcpy(betterSolution, inputSolution_, numberColumns*sizeof(double));
2712            returnCode = 1;
2713        }
2714        delete [] inputSolution_;
2715        inputSolution_ = NULL;
2716        model_ = NULL; // switch off
2717        return returnCode;
[940]2718    }
[264]2719}
2720// update model
2721void CbcSerendipity::setModel(CbcModel * model)
2722{
[1286]2723    model_ = model;
[264]2724}
2725// Resets stuff if model changes
[1286]2726void
[264]2727CbcSerendipity::resetModel(CbcModel * model)
2728{
[1286]2729    model_ = model;
[264]2730}
[961]2731
[1286]2732
[961]2733// Default Constructor
[1286]2734CbcHeuristicJustOne::CbcHeuristicJustOne()
2735        : CbcHeuristic(),
2736        probabilities_(NULL),
2737        heuristic_(NULL),
2738        numberHeuristics_(0)
[961]2739{
2740}
2741
2742// Constructor from model
2743CbcHeuristicJustOne::CbcHeuristicJustOne(CbcModel & model)
[1286]2744        : CbcHeuristic(model),
2745        probabilities_(NULL),
2746        heuristic_(NULL),
2747        numberHeuristics_(0)
[961]2748{
2749}
2750
[1286]2751// Destructor
[961]2752CbcHeuristicJustOne::~CbcHeuristicJustOne ()
2753{
[1286]2754    for (int i = 0; i < numberHeuristics_; i++)
2755        delete heuristic_[i];
2756    delete [] heuristic_;
2757    delete [] probabilities_;
[961]2758}
2759
2760// Clone
2761CbcHeuristicJustOne *
2762CbcHeuristicJustOne::clone() const
2763{
[1286]2764    return new CbcHeuristicJustOne(*this);
[961]2765}
2766
2767// Create C++ lines to get to current state
[1286]2768void
2769CbcHeuristicJustOne::generateCpp( FILE * fp)
[961]2770{
[1286]2771    CbcHeuristicJustOne other;
2772    fprintf(fp, "0#include \"CbcHeuristicJustOne.hpp\"\n");
2773    fprintf(fp, "3  CbcHeuristicJustOne heuristicJustOne(*cbcModel);\n");
2774    CbcHeuristic::generateCpp(fp, "heuristicJustOne");
2775    fprintf(fp, "3  cbcModel->addHeuristic(&heuristicJustOne);\n");
[961]2776}
2777
[1286]2778// Copy constructor
[961]2779CbcHeuristicJustOne::CbcHeuristicJustOne(const CbcHeuristicJustOne & rhs)
[1286]2780        :
2781        CbcHeuristic(rhs),
2782        probabilities_(NULL),
2783        heuristic_(NULL),
2784        numberHeuristics_(rhs.numberHeuristics_)
[961]2785{
[1286]2786    if (numberHeuristics_) {
2787        probabilities_ = CoinCopyOfArray(rhs.probabilities_, numberHeuristics_);
2788        heuristic_ = new CbcHeuristic * [numberHeuristics_];
2789        for (int i = 0; i < numberHeuristics_; i++)
2790            heuristic_[i] = rhs.heuristic_[i]->clone();
2791    }
[961]2792}
2793
[1286]2794// Assignment operator
2795CbcHeuristicJustOne &
2796CbcHeuristicJustOne::operator=( const CbcHeuristicJustOne & rhs)
[961]2797{
[1286]2798    if (this != &rhs) {
2799        CbcHeuristic::operator=(rhs);
2800        for (int i = 0; i < numberHeuristics_; i++)
2801            delete heuristic_[i];
2802        delete [] heuristic_;
2803        delete [] probabilities_;
2804        probabilities_ = NULL;
2805        heuristic_ = NULL;
2806        numberHeuristics_ = rhs.numberHeuristics_;
2807        if (numberHeuristics_) {
2808            probabilities_ = CoinCopyOfArray(rhs.probabilities_, numberHeuristics_);
2809            heuristic_ = new CbcHeuristic * [numberHeuristics_];
2810            for (int i = 0; i < numberHeuristics_; i++)
2811                heuristic_[i] = rhs.heuristic_[i]->clone();
2812        }
[961]2813    }
[1286]2814    return *this;
[961]2815}
2816// Sets value of solution
2817// Returns 1 if solution, 0 if not
2818int
2819CbcHeuristicJustOne::solution(double & solutionValue,
[1286]2820                              double * betterSolution)
[961]2821{
2822#ifdef DIVE_DEBUG
[1286]2823    std::cout << "solutionValue = " << solutionValue << std::endl;
[961]2824#endif
[1286]2825    ++numCouldRun_;
[961]2826
[1286]2827    // test if the heuristic can run
2828    if (!shouldHeurRun_randomChoice() || !numberHeuristics_)
2829        return 0;
2830    double randomNumber = randomNumberGenerator_.randomDouble();
2831    int i;
2832    for (i = 0; i < numberHeuristics_; i++) {
2833        if (randomNumber < probabilities_[i])
2834            break;
2835    }
2836    assert (i < numberHeuristics_);
2837    int returnCode;
2838    //model_->unsetDivingHasRun();
[961]2839#ifdef COIN_DEVELOP
[1286]2840    printf("JustOne running %s\n",
2841           heuristic_[i]->heuristicName());
[961]2842#endif
[1286]2843    returnCode = heuristic_[i]->solution(solutionValue, betterSolution);
[961]2844#ifdef COIN_DEVELOP
[1286]2845    if (returnCode)
2846        printf("JustOne running %s found solution\n",
2847               heuristic_[i]->heuristicName());
[961]2848#endif
[1286]2849    return returnCode;
[961]2850}
2851// Resets stuff if model changes
[1286]2852void
[961]2853CbcHeuristicJustOne::resetModel(CbcModel * model)
2854{
[1286]2855    CbcHeuristic::resetModel(model);
2856    for (int i = 0; i < numberHeuristics_; i++)
2857        heuristic_[i]->resetModel(model);
[961]2858}
2859// update model (This is needed if cliques update matrix etc)
[1286]2860void
[961]2861CbcHeuristicJustOne::setModel(CbcModel * model)
2862{
[1286]2863    CbcHeuristic::setModel(model);
2864    for (int i = 0; i < numberHeuristics_; i++)
2865        heuristic_[i]->setModel(model);
[961]2866}
2867// Validate model i.e. sets when_ to 0 if necessary (may be NULL)
[1286]2868void
[961]2869CbcHeuristicJustOne::validate()
2870{
[1286]2871    CbcHeuristic::validate();
2872    for (int i = 0; i < numberHeuristics_; i++)
2873        heuristic_[i]->validate();
[961]2874}
2875// Adds an heuristic with probability
[1286]2876void
[961]2877CbcHeuristicJustOne::addHeuristic(const CbcHeuristic * heuristic, double probability)
2878{
[1286]2879    CbcHeuristic * thisOne = heuristic->clone();
2880    thisOne->setWhen(-999);
2881    CbcHeuristic ** tempH = CoinCopyOfArrayPartial(heuristic_, numberHeuristics_ + 1,
2882                            numberHeuristics_);
2883    delete [] heuristic_;
2884    heuristic_ = tempH;
2885    heuristic_[numberHeuristics_] = thisOne;
2886    double * tempP = CoinCopyOfArrayPartial(probabilities_, numberHeuristics_ + 1,
2887                                            numberHeuristics_);
2888    delete [] probabilities_;
2889    probabilities_ = tempP;
2890    probabilities_[numberHeuristics_] = probability;
2891    numberHeuristics_++;
[961]2892}
2893// Normalize probabilities
[1286]2894void
[961]2895CbcHeuristicJustOne::normalizeProbabilities()
2896{
[1286]2897    double sum = 0.0;
2898    for (int i = 0; i < numberHeuristics_; i++)
2899        sum += probabilities_[i];
2900    double multiplier = 1.0 / sum;
2901    sum = 0.0;
2902    for (int i = 0; i < numberHeuristics_; i++) {
2903        sum += probabilities_[i];
2904        probabilities_[i] = sum * multiplier;
2905    }
2906    assert (fabs(probabilities_[numberHeuristics_-1] - 1.0) < 1.0e-5);
2907    probabilities_[numberHeuristics_-1] = 1.000001;
[961]2908}
Note: See TracBrowser for help on using the repository browser.