source: trunk/Cbc/src/CbcHeuristic.cpp @ 2110

Last change on this file since 2110 was 2110, checked in by tkr, 3 years ago

Make sure number of rows is not zero in rounding heuristic

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