source: stable/2.9/Cbc/src/CbcHeuristic.cpp @ 2123

Last change on this file since 2123 was 2123, checked in by forrest, 4 years ago

fix elapsed subproblem maxtime

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 118.5 KB
Line 
1/* $Id: CbcHeuristic.cpp 2123 2015-01-29 17:58:10Z forrest $ */
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                double startTime=model_->getDblParam(CbcModel::CbcStartSeconds);
978                model.setDblParam(CbcModel::CbcStartSeconds,startTime);
979                // move seed across
980                model.randomNumberGenerator()->setSeed(model_->randomNumberGenerator()->getSeed());
981                if (numberNodes >= 0) {
982                    // normal
983                    model.setSpecialOptions(saveModelOptions | 2048);
984                    if (logLevel <= 1 && feasibilityPumpOptions_ != -3)
985                        model.setLogLevel(0);
986                    else
987                        model.setLogLevel(logLevel);
988                    // No small fathoming
989                    model.setFastNodeDepth(-1);
990                    model.setCutoff(signedCutoff);
991                    model.setStrongStrategy(0);
992                    // Don't do if original fraction > 1.0 and too large
993                    if (fractionSmall_>1.0 && fractionSmall_ < 1000000.0) {
994                      /* 1.4 means -1 nodes if >.4
995                         2.4 means -1 nodes if >.5 and 0 otherwise
996                         3.4 means -1 nodes if >.6 and 0 or 5
997                         4.4 means -1 nodes if >.7 and 0, 5 or 10
998                      */
999                      double fraction = fractionSmall_-floor(fractionSmall_);
1000                      if (ratio>fraction) {
1001                        int type = static_cast<int>(floor(fractionSmall_*0.1));
1002                        int over = static_cast<int>(ceil(ratio-fraction));
1003                        int maxNodes[]={-1,0,5,10};
1004                        if (type>over)
1005                          numberNodes=maxNodes[type-over];
1006                        else
1007                          numberNodes=-1;
1008                      }
1009                    }
1010                    model.setMaximumNodes(numberNodes);
1011                    model.solver()->setHintParam(OsiDoReducePrint, true, OsiHintTry);
1012                    if ((saveModelOptions&2048) == 0)
1013                      model.setMoreSpecialOptions(model_->moreSpecialOptions());
1014                      model.setMoreSpecialOptions2(model_->moreSpecialOptions2());
1015                    // off conflict analysis
1016                    model.setMoreSpecialOptions(model.moreSpecialOptions()&~4194304);
1017                   
1018                    // Lightweight
1019                    CbcStrategyDefaultSubTree strategy(model_, 1, 5, 1, 0);
1020                    model.setStrategy(strategy);
1021                    model.solver()->setIntParam(OsiMaxNumIterationHotStart, 10);
1022                    model.setMaximumCutPassesAtRoot(CoinMin(20, CoinAbs(model_->getMaximumCutPassesAtRoot())));
1023                    model.setMaximumCutPasses(CoinMin(10, model_->getMaximumCutPasses()));
1024                    // Set best solution (even if bad for this submodel)
1025                    if (model_->bestSolution()) {
1026                      const double * bestSolution = model_->bestSolution();
1027                      int numberColumns2 = model.solver()->getNumCols();
1028                      double * bestSolution2 = new double [numberColumns2];
1029                      const int * originalColumns = process.originalColumns();
1030                      for (int iColumn=0;iColumn<numberColumns2;iColumn++) {
1031                        int jColumn = originalColumns[iColumn];
1032                        bestSolution2[iColumn] = bestSolution[jColumn];
1033                      }
1034                      model.setBestSolution(bestSolution2,numberColumns2,
1035                                            1.0e50,
1036                                            false);
1037                      model.setSolutionCount(1);
1038                      maximumSolutions++; 
1039                      delete [] bestSolution2;
1040                    }
1041                } else {
1042                    model.setSpecialOptions(saveModelOptions);
1043                    model_->messageHandler()->message(CBC_RESTART, model_->messages())
1044                    << solver2->getNumRows() << solver2->getNumCols()
1045                    << CoinMessageEol;
1046                    // going for full search and copy across more stuff
1047                    model.gutsOfCopy(*model_, 2);
1048                    assert (!model_->heuristicModel());
1049                    model_->setHeuristicModel(&model);
1050                    for (int i = 0; i < model.numberCutGenerators(); i++) {
1051                        CbcCutGenerator * generator = model.cutGenerator(i);
1052                        CglGomory * gomory = dynamic_cast<CglGomory *>
1053                          (generator->generator());
1054                        if (gomory&&gomory->originalSolver()) 
1055                          gomory->passInOriginalSolver(model.solver());
1056                        generator->setTiming(true);
1057                        // Turn on if was turned on
1058                        int iOften = model_->cutGenerator(i)->howOften();
1059#ifdef CLP_INVESTIGATE
1060                        printf("Gen %d often %d %d\n",
1061                               i, generator->howOften(),
1062                               iOften);
1063#endif
1064                        if (iOften > 0)
1065                            generator->setHowOften(iOften % 1000000);
1066                        if (model_->cutGenerator(i)->howOftenInSub() == -200)
1067                            generator->setHowOften(-100);
1068                    }
1069                    model.setCutoff(signedCutoff);
1070                    // make sure can't do nested search! but allow heuristics
1071                    model.setSpecialOptions((model.specialOptions()&(~(512 + 2048))) | 1024);
1072                    bool takeHint;
1073                    OsiHintStrength strength;
1074                    // Switch off printing if asked to
1075                    model_->solver()->getHintParam(OsiDoReducePrint, takeHint, strength);
1076                    model.solver()->setHintParam(OsiDoReducePrint, takeHint, strength);
1077                    // no cut generators if none in parent
1078                    CbcStrategyDefault
1079                      strategy(model_->numberCutGenerators() ? 1 : -1, 
1080                               model_->numberStrong(),
1081                               model_->numberBeforeTrust());
1082                    // Set up pre-processing - no
1083                    strategy.setupPreProcessing(0); // was (4);
1084                    model.setStrategy(strategy);
1085                    //model.solver()->writeMps("crunched");
1086                    int numberCuts = process.cuts().sizeRowCuts();
1087                    if (numberCuts) {
1088                        // add in cuts
1089                        CglStored cuts = process.cuts();
1090                        model.addCutGenerator(&cuts, 1, "Stored from first");
1091                        model.cutGenerator(model.numberCutGenerators()-1)->setGlobalCuts(true);
1092                    }
1093                }
1094                // Do search
1095                if (logLevel > 1)
1096                    model_->messageHandler()->message(CBC_START_SUB, model_->messages())
1097                    << name
1098                    << model.getMaximumNodes()
1099                    << CoinMessageEol;
1100                // probably faster to use a basis to get integer solutions
1101                model.setSpecialOptions(model.specialOptions() | 2);
1102#ifdef CBC_THREAD
1103                if (model_->getNumberThreads() > 0 && (model_->getThreadMode()&4) != 0) {
1104                    // See if at root node
1105                    bool atRoot = model_->getNodeCount() == 0;
1106                    int passNumber = model_->getCurrentPassNumber();
1107                    if (atRoot && passNumber == 1)
1108                        model.setNumberThreads(model_->getNumberThreads());
1109                }
1110#endif
1111                model.setParentModel(*model_);
1112                model.setMaximumSolutions(maximumSolutions); 
1113                model.setOriginalColumns(process.originalColumns());
1114                model.setSearchStrategy(-1);
1115                // If no feasibility pump then insert a lightweight one
1116                if (feasibilityPumpOptions_ >= 0 || feasibilityPumpOptions_ == -2) {
1117                    CbcHeuristicFPump * fpump = NULL;
1118                    for (int i = 0; i < model.numberHeuristics(); i++) {
1119                        CbcHeuristicFPump* pump =
1120                            dynamic_cast<CbcHeuristicFPump*>(model.heuristic(i));
1121                        if (pump) {
1122                            fpump = pump;
1123                            break;
1124                        }
1125                    }
1126                    if (!fpump) {
1127                        CbcHeuristicFPump heuristic4;
1128                        // use any cutoff
1129                        heuristic4.setFakeCutoff(0.5*COIN_DBL_MAX);
1130                        if (fractionSmall_<=1.0) 
1131                          heuristic4.setMaximumPasses(10);
1132                        int pumpTune = feasibilityPumpOptions_;
1133                        if (pumpTune==-2)
1134                          pumpTune = 4; // proximity
1135                        if (pumpTune > 0) {
1136                            /*
1137                            >=10000000 for using obj
1138                            >=1000000 use as accumulate switch
1139                            >=1000 use index+1 as number of large loops
1140                            >=100 use 0.05 objvalue as increment
1141                            %100 == 10,20 etc for experimentation
1142                            1 == fix ints at bounds, 2 fix all integral ints, 3 and continuous at bounds
1143                            4 and static continuous, 5 as 3 but no internal integers
1144                            6 as 3 but all slack basis!
1145                            */
1146                            double value = solver2->getObjSense() * solver2->getObjValue();
1147                            int w = pumpTune / 10;
1148                            int ix = w % 10;
1149                            w /= 10;
1150                            int c = w % 10;
1151                            w /= 10;
1152                            int r = w;
1153                            int accumulate = r / 1000;
1154                            r -= 1000 * accumulate;
1155                            if (accumulate >= 10) {
1156                                int which = accumulate / 10;
1157                                accumulate -= 10 * which;
1158                                which--;
1159                                // weights and factors
1160                                double weight[] = {0.1, 0.1, 0.5, 0.5, 1.0, 1.0, 5.0, 5.0};
1161                                double factor[] = {0.1, 0.5, 0.1, 0.5, 0.1, 0.5, 0.1, 0.5};
1162                                heuristic4.setInitialWeight(weight[which]);
1163                                heuristic4.setWeightFactor(factor[which]);
1164                            }
1165                            // fake cutoff
1166                            if (c) {
1167                                double cutoff;
1168                                solver2->getDblParam(OsiDualObjectiveLimit, cutoff);
1169                                cutoff = CoinMin(cutoff, value + 0.1 * fabs(value) * c);
1170                                heuristic4.setFakeCutoff(cutoff);
1171                            }
1172                            if (r) {
1173                                // also set increment
1174                                //double increment = (0.01*i+0.005)*(fabs(value)+1.0e-12);
1175                                double increment = 0.0;
1176                                heuristic4.setAbsoluteIncrement(increment);
1177                                heuristic4.setAccumulate(accumulate);
1178                                heuristic4.setMaximumRetries(r + 1);
1179                            }
1180                            pumpTune = pumpTune % 100;
1181                            if (pumpTune == 6)
1182                                pumpTune = 13;
1183                            if (pumpTune != 13)
1184                                pumpTune = pumpTune % 10;
1185                            heuristic4.setWhen(pumpTune);
1186                            if (ix) {
1187                                heuristic4.setFeasibilityPumpOptions(ix*10);
1188                            }
1189                        }
1190                        model.addHeuristic(&heuristic4, "feasibility pump", 0);
1191       
1192                    }
1193                } else if (feasibilityPumpOptions_==-3) {
1194                  // add all (except this)
1195                  for (int i = 0; i < model_->numberHeuristics(); i++) {
1196                    if (strcmp(heuristicName(),model_->heuristic(i)->heuristicName()))
1197                      model.addHeuristic(model_->heuristic(i)); 
1198                  }
1199                }
1200                // modify heuristics
1201                for (int i = 0; i < model.numberHeuristics(); i++) {
1202                  // reset lastNode
1203                  CbcHeuristicRINS * rins =
1204                    dynamic_cast<CbcHeuristicRINS*>(model.heuristic(i));
1205                  if (rins) {
1206                    rins->setLastNode(-1000);
1207                    rins->setSolutionCount(0);
1208                  }
1209                }
1210                //printf("sol %x\n",inputSolution_);
1211                if (inputSolution_) {
1212                    // translate and add a serendipity heuristic
1213                    int numberColumns = solver2->getNumCols();
1214                    const int * which = process.originalColumns();
1215                    OsiSolverInterface * solver3 = solver2->clone();
1216                    for (int i = 0; i < numberColumns; i++) {
1217                        if (solver3->isInteger(i)) {
1218                            int k = which[i];
1219                            double value = inputSolution_[k];
1220                            //if (value)
1221                            //printf("orig col %d now %d val %g\n",
1222                            //       k,i,value);
1223                            solver3->setColLower(i, value);
1224                            solver3->setColUpper(i, value);
1225                        }
1226                    }
1227                    solver3->setDblParam(OsiDualObjectiveLimit, COIN_DBL_MAX);
1228                    solver3->resolve();
1229                    if (!solver3->isProvenOptimal()) {
1230                        // Try just setting nonzeros
1231                        OsiSolverInterface * solver4 = solver2->clone();
1232                        for (int i = 0; i < numberColumns; i++) {
1233                            if (solver4->isInteger(i)) {
1234                                int k = which[i];
1235                                double value = floor(inputSolution_[k] + 0.5);
1236                                if (value) {
1237                                    solver3->setColLower(i, value);
1238                                    solver3->setColUpper(i, value);
1239                                }
1240                            }
1241                        }
1242                        solver4->setDblParam(OsiDualObjectiveLimit, COIN_DBL_MAX);
1243                        solver4->resolve();
1244                        int nBad = -1;
1245                        if (solver4->isProvenOptimal()) {
1246                            nBad = 0;
1247                            const double * solution = solver4->getColSolution();
1248                            for (int i = 0; i < numberColumns; i++) {
1249                                if (solver4->isInteger(i)) {
1250                                    double value = floor(solution[i] + 0.5);
1251                                    if (fabs(value - solution[i]) > 1.0e-6)
1252                                        nBad++;
1253                                }
1254                            }
1255                        }
1256                        if (nBad) {
1257                            delete solver4;
1258                        } else {
1259                            delete solver3;
1260                            solver3 = solver4;
1261                        }
1262                    }
1263                    if (solver3->isProvenOptimal()) {
1264                        // good
1265                        CbcSerendipity heuristic(model);
1266                        double value = solver3->getObjSense() * solver3->getObjValue();
1267                        heuristic.setInputSolution(solver3->getColSolution(), value);
1268                        value = value + 1.0e-7*(1.0 + fabs(value));
1269                        value *= solver3->getObjSense();
1270                        model.setCutoff(value);
1271                        model.addHeuristic(&heuristic, "Previous solution", 0);
1272                        //printf("added seren\n");
1273                    } else {
1274                        double value = model_->getMinimizationObjValue();
1275                        value = value + 1.0e-7*(1.0 + fabs(value));
1276                        value *= solver3->getObjSense();
1277                        model.setCutoff(value);
1278                        sprintf(generalPrint, "Unable to insert previous solution - using cutoff of %g",
1279                                value);
1280                        model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
1281                        << generalPrint
1282                        << CoinMessageEol;
1283#ifdef CLP_INVESTIGATE
1284                        printf("NOT added seren\n");
1285                        solver3->writeMps("bad_seren");
1286                        solver->writeMps("orig_seren");
1287#endif
1288                    }
1289                    delete solver3;
1290                }
1291                if (model_->searchStrategy() == 2) {
1292                    model.setNumberStrong(5);
1293                    model.setNumberBeforeTrust(5);
1294                }
1295                if (model.getNumCols()) {
1296                    if (numberNodes >= 0) {
1297                        setCutAndHeuristicOptions(model);
1298                        // not too many iterations
1299                        model.setMaximumNumberIterations(iterationMultiplier*(numberNodes + 10));
1300                        // Not fast stuff
1301                        model.setFastNodeDepth(-1);
1302                        //model.solver()->writeMps("before");
1303                    } else if (model.fastNodeDepth() >= 1000000) {
1304                        // already set
1305                        model.setFastNodeDepth(model.fastNodeDepth() - 1000000);
1306                    }
1307                    model.setWhenCuts(999998);
1308#define ALWAYS_DUAL
1309#ifdef ALWAYS_DUAL
1310                    OsiSolverInterface * solverD = model.solver();
1311                    bool takeHint;
1312                    OsiHintStrength strength;
1313                    solverD->getHintParam(OsiDoDualInResolve, takeHint, strength);
1314                    solverD->setHintParam(OsiDoDualInResolve, true, OsiHintDo);
1315#endif
1316                    model.passInEventHandler(model_->getEventHandler());
1317                    // say model_ is sitting there
1318                    int saveOptions = model_->specialOptions();
1319                    model_->setSpecialOptions(saveOptions|1048576);
1320                    // and switch off debugger
1321                    model.setSpecialOptions(model.specialOptions()&(~1));
1322#if 0 //def COIN_HAS_CLP
1323                    OsiClpSolverInterface * clpSolver
1324                      = dynamic_cast<OsiClpSolverInterface *> (model.solver());
1325                    if (clpSolver)
1326                      clpSolver->zapDebugger();
1327#endif
1328#ifdef CONFLICT_CUTS
1329                    if ((model_->moreSpecialOptions()&4194304)!=0)
1330                      model.zapGlobalCuts();
1331#endif
1332                    model.branchAndBound();
1333                    model_->setHeuristicModel(NULL);
1334                    model_->setSpecialOptions(saveOptions);
1335#ifdef ALWAYS_DUAL
1336                    solverD = model.solver();
1337                    solverD->setHintParam(OsiDoDualInResolve, takeHint, strength);
1338#endif
1339                    numberNodesDone_ = model.getNodeCount();
1340#ifdef COIN_DEVELOP
1341                    printf("sub branch %d nodes, %d iterations - max %d\n",
1342                           model.getNodeCount(), model.getIterationCount(),
1343                           100*(numberNodes + 10));
1344#endif
1345                    if (numberNodes < 0) {
1346                        model_->incrementIterationCount(model.getIterationCount());
1347                        model_->incrementNodeCount(model.getNodeCount());
1348                        // update best solution (in case ctrl-c)
1349                        // !!! not a good idea - think a bit harder
1350                        //model_->setMinimizationObjValue(model.getMinimizationObjValue());
1351                        for (int iGenerator = 0; iGenerator < model.numberCutGenerators(); iGenerator++) {
1352                            CbcCutGenerator * generator = model.cutGenerator(iGenerator);
1353                            sprintf(generalPrint,
1354                                    "%s was tried %d times and created %d cuts of which %d were active after adding rounds of cuts (%.3f seconds)",
1355                                    generator->cutGeneratorName(),
1356                                    generator->numberTimesEntered(),
1357                                    generator->numberCutsInTotal() +
1358                                    generator->numberColumnCuts(),
1359                                    generator->numberCutsActive(),
1360                                    generator->timeInCutGenerator());
1361                            CglStored * stored = dynamic_cast<CglStored*>(generator->generator());
1362                            if (stored && !generator->numberCutsInTotal())
1363                                continue;
1364#ifndef CLP_INVESTIGATE
1365                            CglImplication * implication = dynamic_cast<CglImplication*>(generator->generator());
1366                            if (implication)
1367                                continue;
1368#endif
1369                            model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
1370                            << generalPrint
1371                            << CoinMessageEol;
1372                        }
1373                    }
1374                } else {
1375                    // empty model
1376                    model.setMinimizationObjValue(model.solver()->getObjSense()*model.solver()->getObjValue());
1377                }
1378                if (logLevel > 1)
1379                    model_->messageHandler()->message(CBC_END_SUB, model_->messages())
1380                    << name
1381                    << CoinMessageEol;
1382                if (model.getMinimizationObjValue() < CoinMin(cutoff, 1.0e30)) {
1383                    // solution
1384                    if (model.getNumCols())
1385                        returnCode = model.isProvenOptimal() ? 3 : 1;
1386                    else
1387                        returnCode = 3;
1388                    // post process
1389#ifdef COIN_HAS_CLP
1390                    OsiClpSolverInterface * clpSolver = dynamic_cast< OsiClpSolverInterface*> (model.solver());
1391                    if (clpSolver) {
1392                        ClpSimplex * lpSolver = clpSolver->getModelPtr();
1393                        lpSolver->setSpecialOptions(lpSolver->specialOptions() | 0x01000000); // say is Cbc (and in branch and bound)
1394                    }
1395#endif
1396                    //if (fractionSmall_ < 1000000.0)
1397                      process.postProcess(*model.solver());
1398                    if (solver->isProvenOptimal() && solver->getObjValue()*solver->getObjSense() < cutoff) {
1399                        // Solution now back in solver
1400                        int numberColumns = solver->getNumCols();
1401                        memcpy(newSolution, solver->getColSolution(),
1402                               numberColumns*sizeof(double));
1403                        newSolutionValue = model.getMinimizationObjValue();
1404                    } else {
1405                        // odd - but no good
1406                        returnCode = 0; // so will be infeasible
1407                    }
1408                } else {
1409                    // no good
1410                    returnCode = model.isProvenInfeasible() ? 2 : 0; // so will be infeasible
1411                }
1412                int totalNumberIterations = model.getIterationCount() +
1413                                            process.numberIterationsPre() +
1414                                            process.numberIterationsPost();
1415                if (totalNumberIterations > 100*(numberNodes + 10)
1416                    && fractionSmall_ < 1000000.0) {
1417                    // only allow smaller problems
1418                    fractionSmall = fractionSmall_;
1419                    fractionSmall_ *= 0.9;
1420#ifdef CLP_INVESTIGATE
1421                    printf("changing fractionSmall from %g to %g for %s as %d iterations\n",
1422                           fractionSmall, fractionSmall_, name.c_str(), totalNumberIterations);
1423#endif
1424                }
1425                if (model.status() == 5)
1426                   model_->sayEventHappened();
1427#ifdef COIN_DEVELOP
1428                if (model.isProvenInfeasible())
1429                    status = 1;
1430                else if (model.isProvenOptimal())
1431                    status = 2;
1432#endif
1433            }
1434        }
1435    } else {
1436        returnCode = 2; // infeasible finished
1437        if (logLevel > 1){
1438           printf("Infeasible on initial solve\n");
1439        }
1440    }
1441    model_->setSpecialOptions(saveModelOptions);
1442    model_->setLogLevel(logLevel);
1443    if (returnCode == 1 || returnCode == 2) {
1444        OsiSolverInterface * solverC = model_->continuousSolver();
1445        if (false && solverC) {
1446            const double * lower = solver->getColLower();
1447            const double * upper = solver->getColUpper();
1448            const double * lowerC = solverC->getColLower();
1449            const double * upperC = solverC->getColUpper();
1450            bool good = true;
1451            for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
1452                if (solverC->isInteger(iColumn)) {
1453                    if (lower[iColumn] > lowerC[iColumn] &&
1454                            upper[iColumn] < upperC[iColumn]) {
1455                        good = false;
1456                        printf("CUT - can't add\n");
1457                        break;
1458                    }
1459                }
1460            }
1461            if (good) {
1462                double * cut = new double [numberColumns];
1463                int * which = new int [numberColumns];
1464                double rhs = -1.0;
1465                int n = 0;
1466                for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
1467                    if (solverC->isInteger(iColumn)) {
1468                        if (lower[iColumn] == upperC[iColumn]) {
1469                            rhs += lower[iColumn];
1470                            cut[n] = 1.0;
1471                            which[n++] = iColumn;
1472                        } else if (upper[iColumn] == lowerC[iColumn]) {
1473                            rhs -= upper[iColumn];
1474                            cut[n] = -1.0;
1475                            which[n++] = iColumn;
1476                        }
1477                    }
1478                }
1479                printf("CUT has %d entries\n", n);
1480                OsiRowCut newCut;
1481                newCut.setLb(-COIN_DBL_MAX);
1482                newCut.setUb(rhs);
1483                newCut.setRow(n, which, cut, false);
1484                model_->makeGlobalCut(newCut);
1485                delete [] cut;
1486                delete [] which;
1487            }
1488        }
1489#ifdef COIN_DEVELOP
1490        if (status == 1)
1491            printf("heuristic could add cut because infeasible (%s)\n", heuristicName_.c_str());
1492        else if (status == 2)
1493            printf("heuristic could add cut because optimal (%s)\n", heuristicName_.c_str());
1494#endif
1495    }
1496    if (reset) {
1497        for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
1498            if (reset[iColumn])
1499                solver->setColLower(iColumn, 0.0);
1500        }
1501        delete [] reset;
1502    }
1503#ifdef HISTORY_STATISTICS
1504    getHistoryStatistics_ = true;
1505#endif
1506    solver->setHintParam(OsiDoReducePrint, takeHint, strength);
1507    return returnCode;
1508}
1509// Set input solution
1510void
1511CbcHeuristic::setInputSolution(const double * solution, double objValue)
1512{
1513    delete [] inputSolution_;
1514    inputSolution_ = NULL;
1515    if (model_ && solution) {
1516        int numberColumns = model_->getNumCols();
1517        inputSolution_ = new double [numberColumns+1];
1518        memcpy(inputSolution_, solution, numberColumns*sizeof(double));
1519        inputSolution_[numberColumns] = objValue;
1520    }
1521}
1522
1523//##############################################################################
1524
1525inline int compare3BranchingObjects(const CbcBranchingObject* br0,
1526                                    const CbcBranchingObject* br1)
1527{
1528    const int t0 = br0->type();
1529    const int t1 = br1->type();
1530    if (t0 < t1) {
1531        return -1;
1532    }
1533    if (t0 > t1) {
1534        return 1;
1535    }
1536    return br0->compareOriginalObject(br1);
1537}
1538
1539//==============================================================================
1540
1541inline bool compareBranchingObjects(const CbcBranchingObject* br0,
1542                                    const CbcBranchingObject* br1)
1543{
1544    return compare3BranchingObjects(br0, br1) < 0;
1545}
1546
1547//==============================================================================
1548
1549void
1550CbcHeuristicNode::gutsOfConstructor(CbcModel& model)
1551{
1552    //  CbcHeurDebugNodes(&model);
1553    CbcNode* node = model.currentNode();
1554    brObj_ = new CbcBranchingObject*[node->depth()];
1555    CbcNodeInfo* nodeInfo = node->nodeInfo();
1556    int cnt = 0;
1557    while (nodeInfo->parentBranch() != NULL) {
1558        const OsiBranchingObject* br = nodeInfo->parentBranch();
1559        const CbcBranchingObject* cbcbr = dynamic_cast<const CbcBranchingObject*>(br);
1560        if (! cbcbr) {
1561            throw CoinError("CbcHeuristicNode can be used only with CbcBranchingObjects.\n",
1562                            "gutsOfConstructor",
1563                            "CbcHeuristicNode",
1564                            __FILE__, __LINE__);
1565        }
1566        brObj_[cnt] = cbcbr->clone();
1567        brObj_[cnt]->previousBranch();
1568        ++cnt;
1569        nodeInfo = nodeInfo->parent();
1570    }
1571    std::sort(brObj_, brObj_ + cnt, compareBranchingObjects);
1572    if (cnt <= 1) {
1573        numObjects_ = cnt;
1574    } else {
1575        numObjects_ = 0;
1576        CbcBranchingObject* br = NULL; // What should this be?
1577        for (int i = 1; i < cnt; ++i) {
1578            if (compare3BranchingObjects(brObj_[numObjects_], brObj_[i]) == 0) {
1579                int comp = brObj_[numObjects_]->compareBranchingObject(brObj_[i], br != 0);
1580                switch (comp) {
1581                case CbcRangeSame: // the same range
1582                case CbcRangeDisjoint: // disjoint decisions
1583                    // should not happen! we are on a chain!
1584                    abort();
1585                case CbcRangeSubset: // brObj_[numObjects_] is a subset of brObj_[i]
1586                    delete brObj_[i];
1587                    break;
1588                case CbcRangeSuperset: // brObj_[i] is a subset of brObj_[numObjects_]
1589                    delete brObj_[numObjects_];
1590                    brObj_[numObjects_] = brObj_[i];
1591                    break;
1592                case CbcRangeOverlap: // overlap
1593                    delete brObj_[i];
1594                    delete brObj_[numObjects_];
1595                    brObj_[numObjects_] = br;
1596                    break;
1597                }
1598                continue;
1599            } else {
1600                brObj_[++numObjects_] = brObj_[i];
1601            }
1602        }
1603        ++numObjects_;
1604    }
1605}
1606
1607//==============================================================================
1608
1609CbcHeuristicNode::CbcHeuristicNode(CbcModel& model)
1610{
1611    gutsOfConstructor(model);
1612}
1613
1614//==============================================================================
1615
1616double
1617CbcHeuristicNode::distance(const CbcHeuristicNode* node) const
1618{
1619
1620    const double disjointWeight = 1;
1621    const double overlapWeight = 0.4;
1622    const double subsetWeight = 0.2;
1623    int countDisjointWeight = 0;
1624    int countOverlapWeight = 0;
1625    int countSubsetWeight = 0;
1626    int i = 0;
1627    int j = 0;
1628    double dist = 0.0;
1629#ifdef PRINT_DEBUG
1630    printf(" numObjects_ = %i, node->numObjects_ = %i\n",
1631           numObjects_, node->numObjects_);
1632#endif
1633    while ( i < numObjects_ && j < node->numObjects_) {
1634        CbcBranchingObject* br0 = brObj_[i];
1635        const CbcBranchingObject* br1 = node->brObj_[j];
1636#ifdef PRINT_DEBUG
1637        const CbcIntegerBranchingObject* brPrint0 =
1638            dynamic_cast<const CbcIntegerBranchingObject*>(br0);
1639        const double* downBounds = brPrint0->downBounds();
1640        const double* upBounds = brPrint0->upBounds();
1641        int variable = brPrint0->variable();
1642        int way = brPrint0->way();
1643        printf("   br0: var %i downBd [%i,%i] upBd [%i,%i] way %i\n",
1644               variable, static_cast<int>(downBounds[0]), static_cast<int>(downBounds[1]),
1645               static_cast<int>(upBounds[0]), static_cast<int>(upBounds[1]), way);
1646        const CbcIntegerBranchingObject* brPrint1 =
1647            dynamic_cast<const CbcIntegerBranchingObject*>(br1);
1648        downBounds = brPrint1->downBounds();
1649        upBounds = brPrint1->upBounds();
1650        variable = brPrint1->variable();
1651        way = brPrint1->way();
1652        printf("   br1: var %i downBd [%i,%i] upBd [%i,%i] way %i\n",
1653               variable, static_cast<int>(downBounds[0]), static_cast<int>(downBounds[1]),
1654               static_cast<int>(upBounds[0]), static_cast<int>(upBounds[1]), way);
1655#endif
1656        const int brComp = compare3BranchingObjects(br0, br1);
1657        if (brComp < 0) {
1658            dist += subsetWeight;
1659            countSubsetWeight++;
1660            ++i;
1661        } else if (brComp > 0) {
1662            dist += subsetWeight;
1663            countSubsetWeight++;
1664            ++j;
1665        } else {
1666            const int comp = br0->compareBranchingObject(br1, false);
1667            switch (comp) {
1668            case CbcRangeSame:
1669                // do nothing
1670                break;
1671            case CbcRangeDisjoint: // disjoint decisions
1672                dist += disjointWeight;
1673                countDisjointWeight++;
1674                break;
1675            case CbcRangeSubset: // subset one way or another
1676            case CbcRangeSuperset:
1677                dist += subsetWeight;
1678                countSubsetWeight++;
1679                break;
1680            case CbcRangeOverlap: // overlap
1681                dist += overlapWeight;
1682                countOverlapWeight++;
1683                break;
1684            }
1685            ++i;
1686            ++j;
1687        }
1688    }
1689    dist += subsetWeight * (numObjects_ - i + node->numObjects_ - j);
1690    countSubsetWeight += (numObjects_ - i + node->numObjects_ - j);
1691    COIN_DETAIL_PRINT(printf("subset = %i, overlap = %i, disjoint = %i\n", countSubsetWeight,
1692                             countOverlapWeight, countDisjointWeight));
1693    return dist;
1694}
1695
1696//==============================================================================
1697
1698CbcHeuristicNode::~CbcHeuristicNode()
1699{
1700    for (int i = 0; i < numObjects_; ++i) {
1701        delete brObj_[i];
1702    }
1703    delete [] brObj_;
1704}
1705
1706//==============================================================================
1707
1708double
1709CbcHeuristicNode::minDistance(const CbcHeuristicNodeList& nodeList) const
1710{
1711    double minDist = COIN_DBL_MAX;
1712    for (int i = nodeList.size() - 1; i >= 0; --i) {
1713        minDist = CoinMin(minDist, distance(nodeList.node(i)));
1714    }
1715    return minDist;
1716}
1717
1718//==============================================================================
1719
1720bool
1721CbcHeuristicNode::minDistanceIsSmall(const CbcHeuristicNodeList& nodeList,
1722                                     const double threshold) const
1723{
1724    for (int i = nodeList.size() - 1; i >= 0; --i) {
1725        if (distance(nodeList.node(i)) >= threshold) {
1726            continue;
1727        } else {
1728            return true;
1729        }
1730    }
1731    return false;
1732}
1733
1734//==============================================================================
1735
1736double
1737CbcHeuristicNode::avgDistance(const CbcHeuristicNodeList& nodeList) const
1738{
1739    if (nodeList.size() == 0) {
1740        return COIN_DBL_MAX;
1741    }
1742    double sumDist = 0;
1743    for (int i = nodeList.size() - 1; i >= 0; --i) {
1744        sumDist += distance(nodeList.node(i));
1745    }
1746    return sumDist / nodeList.size();
1747}
1748
1749//##############################################################################
1750
1751// Default Constructor
1752CbcRounding::CbcRounding()
1753        : CbcHeuristic()
1754{
1755    // matrix and row copy will automatically be empty
1756    seed_ = 7654321;
1757    down_ = NULL;
1758    up_ = NULL;
1759    equal_ = NULL;
1760    //whereFrom_ |= 16*(1+256); // allow more often
1761}
1762
1763// Constructor from model
1764CbcRounding::CbcRounding(CbcModel & model)
1765        : CbcHeuristic(model)
1766{
1767    // Get a copy of original matrix (and by row for rounding);
1768    assert(model.solver());
1769    if (model.solver()->getNumRows()) {
1770        matrix_ = *model.solver()->getMatrixByCol();
1771        matrixByRow_ = *model.solver()->getMatrixByRow();
1772        validate();
1773    }
1774    down_ = NULL;
1775    up_ = NULL;
1776    equal_ = NULL;
1777    seed_ = 7654321;
1778    //whereFrom_ |= 16*(1+256); // allow more often
1779}
1780
1781// Destructor
1782CbcRounding::~CbcRounding ()
1783{
1784    delete [] down_;
1785    delete [] up_;
1786    delete [] equal_;
1787}
1788
1789// Clone
1790CbcHeuristic *
1791CbcRounding::clone() const
1792{
1793    return new CbcRounding(*this);
1794}
1795// Create C++ lines to get to current state
1796void
1797CbcRounding::generateCpp( FILE * fp)
1798{
1799    CbcRounding other;
1800    fprintf(fp, "0#include \"CbcHeuristic.hpp\"\n");
1801    fprintf(fp, "3  CbcRounding rounding(*cbcModel);\n");
1802    CbcHeuristic::generateCpp(fp, "rounding");
1803    if (seed_ != other.seed_)
1804        fprintf(fp, "3  rounding.setSeed(%d);\n", seed_);
1805    else
1806        fprintf(fp, "4  rounding.setSeed(%d);\n", seed_);
1807    fprintf(fp, "3  cbcModel->addHeuristic(&rounding);\n");
1808}
1809//#define NEW_ROUNDING
1810// Copy constructor
1811CbcRounding::CbcRounding(const CbcRounding & rhs)
1812        :
1813        CbcHeuristic(rhs),
1814        matrix_(rhs.matrix_),
1815        matrixByRow_(rhs.matrixByRow_),
1816        seed_(rhs.seed_)
1817{
1818#ifdef NEW_ROUNDING
1819    int numberColumns = matrix_.getNumCols();
1820    down_ = CoinCopyOfArray(rhs.down_, numberColumns);
1821    up_ = CoinCopyOfArray(rhs.up_, numberColumns);
1822    equal_ = CoinCopyOfArray(rhs.equal_, numberColumns);
1823#else
1824    down_ = NULL;
1825    up_ = NULL;
1826    equal_ = NULL;
1827#endif
1828}
1829
1830// Assignment operator
1831CbcRounding &
1832CbcRounding::operator=( const CbcRounding & rhs)
1833{
1834    if (this != &rhs) {
1835        CbcHeuristic::operator=(rhs);
1836        matrix_ = rhs.matrix_;
1837        matrixByRow_ = rhs.matrixByRow_;
1838#ifdef NEW_ROUNDING
1839        delete [] down_;
1840        delete [] up_;
1841        delete [] equal_;
1842        int numberColumns = matrix_.getNumCols();
1843        down_ = CoinCopyOfArray(rhs.down_, numberColumns);
1844        up_ = CoinCopyOfArray(rhs.up_, numberColumns);
1845        equal_ = CoinCopyOfArray(rhs.equal_, numberColumns);
1846#else
1847        down_ = NULL;
1848        up_ = NULL;
1849        equal_ = NULL;
1850#endif
1851        seed_ = rhs.seed_;
1852    }
1853    return *this;
1854}
1855
1856// Resets stuff if model changes
1857void
1858CbcRounding::resetModel(CbcModel * model)
1859{
1860    model_ = model;
1861    // Get a copy of original matrix (and by row for rounding);
1862    assert(model_->solver());
1863    matrix_ = *model_->solver()->getMatrixByCol();
1864    matrixByRow_ = *model_->solver()->getMatrixByRow();
1865    validate();
1866}
1867/* Check whether the heuristic should run at all
1868   0 - before cuts at root node (or from doHeuristics)
1869   1 - during cuts at root
1870   2 - after root node cuts
1871   3 - after cuts at other nodes
1872   4 - during cuts at other nodes
1873   8 added if previous heuristic in loop found solution
1874*/
1875bool 
1876CbcRounding::shouldHeurRun(int whereFrom)
1877{
1878  if (whereFrom!=4) {
1879    return CbcHeuristic::shouldHeurRun(whereFrom);
1880  } else {
1881    numCouldRun_++;
1882    return shouldHeurRun_randomChoice();
1883  }
1884}
1885// See if rounding will give solution
1886// Sets value of solution
1887// Assumes rhs for original matrix still okay
1888// At present only works with integers
1889// Fix values if asked for
1890// Returns 1 if solution, 0 if not
1891int
1892CbcRounding::solution(double & solutionValue,
1893                      double * betterSolution)
1894{
1895
1896    numCouldRun_++;
1897    // See if to do
1898    if (!when() || (when() % 10 == 1 && model_->phase() != 1) ||
1899            (when() % 10 == 2 && (model_->phase() != 2 && model_->phase() != 3)))
1900        return 0; // switched off
1901    numRuns_++;
1902#ifdef HEURISTIC_INFORM
1903    printf("Entering heuristic %s - nRuns %d numCould %d when %d\n",
1904           heuristicName(),numRuns_,numCouldRun_,when_);
1905#endif
1906    OsiSolverInterface * solver = model_->solver();
1907    double direction = solver->getObjSense();
1908    double newSolutionValue = direction * solver->getObjValue();
1909    return solution(solutionValue, betterSolution, newSolutionValue);
1910}
1911// See if rounding will give solution
1912// Sets value of solution
1913// Assumes rhs for original matrix still okay
1914// At present only works with integers
1915// Fix values if asked for
1916// Returns 1 if solution, 0 if not
1917int
1918CbcRounding::solution(double & solutionValue,
1919                      double * betterSolution,
1920                      double newSolutionValue)
1921{
1922
1923    // See if to do
1924    if (!when() || (when() % 10 == 1 && model_->phase() != 1) ||
1925            (when() % 10 == 2 && (model_->phase() != 2 && model_->phase() != 3)))
1926        return 0; // switched off
1927    OsiSolverInterface * solver = model_->solver();
1928    const double * lower = solver->getColLower();
1929    const double * upper = solver->getColUpper();
1930    const double * rowLower = solver->getRowLower();
1931    const double * rowUpper = solver->getRowUpper();
1932    const double * solution = solver->getColSolution();
1933    const double * objective = solver->getObjCoefficients();
1934    double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
1935    double primalTolerance;
1936    solver->getDblParam(OsiPrimalTolerance, primalTolerance);
1937    double useTolerance = primalTolerance;
1938
1939    int numberRows = matrix_.getNumRows();
1940    assert (numberRows <= solver->getNumRows());
1941    if (numberRows == 0){
1942       return 0;
1943    }
1944    int numberIntegers = model_->numberIntegers();
1945    const int * integerVariable = model_->integerVariable();
1946    int i;
1947    double direction = solver->getObjSense();
1948    //double newSolutionValue = direction*solver->getObjValue();
1949    int returnCode = 0;
1950    // Column copy
1951    const double * element = matrix_.getElements();
1952    const int * row = matrix_.getIndices();
1953    const CoinBigIndex * columnStart = matrix_.getVectorStarts();
1954    const int * columnLength = matrix_.getVectorLengths();
1955    // Row copy
1956    const double * elementByRow = matrixByRow_.getElements();
1957    const int * column = matrixByRow_.getIndices();
1958    const CoinBigIndex * rowStart = matrixByRow_.getVectorStarts();
1959    const int * rowLength = matrixByRow_.getVectorLengths();
1960
1961    // Get solution array for heuristic solution
1962    int numberColumns = solver->getNumCols();
1963    double * newSolution = new double [numberColumns];
1964    memcpy(newSolution, solution, numberColumns*sizeof(double));
1965
1966    double * rowActivity = new double[numberRows];
1967    memset(rowActivity, 0, numberRows*sizeof(double));
1968    for (i = 0; i < numberColumns; i++) {
1969        int j;
1970        double value = newSolution[i];
1971        if (value < lower[i]) {
1972            value = lower[i];
1973            newSolution[i] = value;
1974        } else if (value > upper[i]) {
1975            value = upper[i];
1976            newSolution[i] = value;
1977        }
1978        if (value) {
1979            for (j = columnStart[i];
1980                    j < columnStart[i] + columnLength[i]; j++) {
1981                int iRow = row[j];
1982                rowActivity[iRow] += value * element[j];
1983            }
1984        }
1985    }
1986    // check was feasible - if not adjust (cleaning may move)
1987    for (i = 0; i < numberRows; i++) {
1988        if (rowActivity[i] < rowLower[i]) {
1989            //assert (rowActivity[i]>rowLower[i]-1000.0*primalTolerance);
1990            rowActivity[i] = rowLower[i];
1991        } else if (rowActivity[i] > rowUpper[i]) {
1992            //assert (rowActivity[i]<rowUpper[i]+1000.0*primalTolerance);
1993            rowActivity[i] = rowUpper[i];
1994        }
1995    }
1996    for (i = 0; i < numberIntegers; i++) {
1997        int iColumn = integerVariable[i];
1998        double value = newSolution[iColumn];
1999        //double thisTolerance = integerTolerance;
2000        if (fabs(floor(value + 0.5) - value) > integerTolerance) {
2001            double below = floor(value);
2002            double newValue = newSolution[iColumn];
2003            double cost = direction * objective[iColumn];
2004            double move;
2005            if (cost > 0.0) {
2006                // try up
2007                move = 1.0 - (value - below);
2008            } else if (cost < 0.0) {
2009                // try down
2010                move = below - value;
2011            } else {
2012                // won't be able to move unless we can grab another variable
2013                double randomNumber = randomNumberGenerator_.randomDouble();
2014                // which way?
2015                if (randomNumber < 0.5)
2016                    move = below - value;
2017                else
2018                    move = 1.0 - (value - below);
2019            }
2020            newValue += move;
2021            newSolution[iColumn] = newValue;
2022            newSolutionValue += move * cost;
2023            int j;
2024            for (j = columnStart[iColumn];
2025                    j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2026                int iRow = row[j];
2027                rowActivity[iRow] += move * element[j];
2028            }
2029        }
2030    }
2031
2032    double penalty = 0.0;
2033    const char * integerType = model_->integerType();
2034    // see if feasible - just using singletons
2035    for (i = 0; i < numberRows; i++) {
2036        double value = rowActivity[i];
2037        double thisInfeasibility = 0.0;
2038        if (value < rowLower[i] - primalTolerance)
2039            thisInfeasibility = value - rowLower[i];
2040        else if (value > rowUpper[i] + primalTolerance)
2041            thisInfeasibility = value - rowUpper[i];
2042        if (thisInfeasibility) {
2043            // See if there are any slacks I can use to fix up
2044            // maybe put in coding for multiple slacks?
2045            double bestCost = 1.0e50;
2046            int k;
2047            int iBest = -1;
2048            double addCost = 0.0;
2049            double newValue = 0.0;
2050            double changeRowActivity = 0.0;
2051            double absInfeasibility = fabs(thisInfeasibility);
2052            for (k = rowStart[i]; k < rowStart[i] + rowLength[i]; k++) {
2053                int iColumn = column[k];
2054                // See if all elements help
2055                if (columnLength[iColumn] == 1) {
2056                    double currentValue = newSolution[iColumn];
2057                    double elementValue = elementByRow[k];
2058                    double lowerValue = lower[iColumn];
2059                    double upperValue = upper[iColumn];
2060                    double gap = rowUpper[i] - rowLower[i];
2061                    double absElement = fabs(elementValue);
2062                    if (thisInfeasibility*elementValue > 0.0) {
2063                        // we want to reduce
2064                        if ((currentValue - lowerValue)*absElement >= absInfeasibility) {
2065                            // possible - check if integer
2066                            double distance = absInfeasibility / absElement;
2067                            double thisCost = -direction * objective[iColumn] * distance;
2068                            if (integerType[iColumn]) {
2069                                distance = ceil(distance - useTolerance);
2070                                if (currentValue - distance >= lowerValue - useTolerance) {
2071                                    if (absInfeasibility - distance*absElement < -gap - useTolerance)
2072                                        thisCost = 1.0e100; // no good
2073                                    else
2074                                        thisCost = -direction * objective[iColumn] * distance;
2075                                } else {
2076                                    thisCost = 1.0e100; // no good
2077                                }
2078                            }
2079                            if (thisCost < bestCost) {
2080                                bestCost = thisCost;
2081                                iBest = iColumn;
2082                                addCost = thisCost;
2083                                newValue = currentValue - distance;
2084                                changeRowActivity = -distance * elementValue;
2085                            }
2086                        }
2087                    } else {
2088                        // we want to increase
2089                        if ((upperValue - currentValue)*absElement >= absInfeasibility) {
2090                            // possible - check if integer
2091                            double distance = absInfeasibility / absElement;
2092                            double thisCost = direction * objective[iColumn] * distance;
2093                            if (integerType[iColumn]) {
2094                                distance = ceil(distance - 1.0e-7);
2095                                assert (currentValue - distance <= upperValue + useTolerance);
2096                                if (absInfeasibility - distance*absElement < -gap - useTolerance)
2097                                    thisCost = 1.0e100; // no good
2098                                else
2099                                    thisCost = direction * objective[iColumn] * distance;
2100                            }
2101                            if (thisCost < bestCost) {
2102                                bestCost = thisCost;
2103                                iBest = iColumn;
2104                                addCost = thisCost;
2105                                newValue = currentValue + distance;
2106                                changeRowActivity = distance * elementValue;
2107                            }
2108                        }
2109                    }
2110                }
2111            }
2112            if (iBest >= 0) {
2113                /*printf("Infeasibility of %g on row %d cost %g\n",
2114                  thisInfeasibility,i,addCost);*/
2115                newSolution[iBest] = newValue;
2116                thisInfeasibility = 0.0;
2117                newSolutionValue += addCost;
2118                rowActivity[i] += changeRowActivity;
2119            }
2120            penalty += fabs(thisInfeasibility);
2121        }
2122    }
2123    if (penalty) {
2124        // see if feasible using any
2125        // first continuous
2126        double penaltyChange = 0.0;
2127        int iColumn;
2128        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
2129            if (integerType[iColumn])
2130                continue;
2131            double currentValue = newSolution[iColumn];
2132            double lowerValue = lower[iColumn];
2133            double upperValue = upper[iColumn];
2134            int j;
2135            int anyBadDown = 0;
2136            int anyBadUp = 0;
2137            double upImprovement = 0.0;
2138            double downImprovement = 0.0;
2139            for (j = columnStart[iColumn];
2140                    j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2141                int iRow = row[j];
2142                if (rowUpper[iRow] > rowLower[iRow]) {
2143                    double value = element[j];
2144                    if (rowActivity[iRow] > rowUpper[iRow] + primalTolerance) {
2145                        // infeasible above
2146                        downImprovement += value;
2147                        upImprovement -= value;
2148                        if (value > 0.0)
2149                            anyBadUp++;
2150                        else
2151                            anyBadDown++;
2152                    } else if (rowActivity[iRow] > rowUpper[iRow] - primalTolerance) {
2153                        // feasible at ub
2154                        if (value > 0.0) {
2155                            upImprovement -= value;
2156                            anyBadUp++;
2157                        } else {
2158                            downImprovement += value;
2159                            anyBadDown++;
2160                        }
2161                    } else if (rowActivity[iRow] > rowLower[iRow] + primalTolerance) {
2162                        // feasible in interior
2163                    } else if (rowActivity[iRow] > rowLower[iRow] - primalTolerance) {
2164                        // feasible at lb
2165                        if (value < 0.0) {
2166                            upImprovement += value;
2167                            anyBadUp++;
2168                        } else {
2169                            downImprovement -= value;
2170                            anyBadDown++;
2171                        }
2172                    } else {
2173                        // infeasible below
2174                        downImprovement -= value;
2175                        upImprovement += value;
2176                        if (value < 0.0)
2177                            anyBadUp++;
2178                        else
2179                            anyBadDown++;
2180                    }
2181                } else {
2182                    // equality row
2183                    double value = element[j];
2184                    if (rowActivity[iRow] > rowUpper[iRow] + primalTolerance) {
2185                        // infeasible above
2186                        downImprovement += value;
2187                        upImprovement -= value;
2188                        if (value > 0.0)
2189                            anyBadUp++;
2190                        else
2191                            anyBadDown++;
2192                    } else if (rowActivity[iRow] < rowLower[iRow] - primalTolerance) {
2193                        // infeasible below
2194                        downImprovement -= value;
2195                        upImprovement += value;
2196                        if (value < 0.0)
2197                            anyBadUp++;
2198                        else
2199                            anyBadDown++;
2200                    } else {
2201                        // feasible - no good
2202                        anyBadUp = -1;
2203                        anyBadDown = -1;
2204                        break;
2205                    }
2206                }
2207            }
2208            // could change tests for anyBad
2209            if (anyBadUp)
2210                upImprovement = 0.0;
2211            if (anyBadDown)
2212                downImprovement = 0.0;
2213            double way = 0.0;
2214            double improvement = 0.0;
2215            if (downImprovement > 0.0 && currentValue > lowerValue) {
2216                way = -1.0;
2217                improvement = downImprovement;
2218            } else if (upImprovement > 0.0 && currentValue < upperValue) {
2219                way = 1.0;
2220                improvement = upImprovement;
2221            }
2222            if (way) {
2223                // can improve
2224                double distance;
2225                if (way > 0.0)
2226                    distance = upperValue - currentValue;
2227                else
2228                    distance = currentValue - lowerValue;
2229                for (j = columnStart[iColumn];
2230                        j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2231                    int iRow = row[j];
2232                    double value = element[j] * way;
2233                    if (rowActivity[iRow] > rowUpper[iRow] + primalTolerance) {
2234                        // infeasible above
2235                        assert (value < 0.0);
2236                        double gap = rowActivity[iRow] - rowUpper[iRow];
2237                        if (gap + value*distance < 0.0)
2238                            distance = -gap / value;
2239                    } else if (rowActivity[iRow] < rowLower[iRow] - primalTolerance) {
2240                        // infeasible below
2241                        assert (value > 0.0);
2242                        double gap = rowActivity[iRow] - rowLower[iRow];
2243                        if (gap + value*distance > 0.0)
2244                            distance = -gap / value;
2245                    } else {
2246                        // feasible
2247                        if (value > 0) {
2248                            double gap = rowActivity[iRow] - rowUpper[iRow];
2249                            if (gap + value*distance > 0.0)
2250                                distance = -gap / value;
2251                        } else {
2252                            double gap = rowActivity[iRow] - rowLower[iRow];
2253                            if (gap + value*distance < 0.0)
2254                                distance = -gap / value;
2255                        }
2256                    }
2257                }
2258                //move
2259                penaltyChange += improvement * distance;
2260                distance *= way;
2261                newSolution[iColumn] += distance;
2262                newSolutionValue += direction * objective[iColumn] * distance;
2263                for (j = columnStart[iColumn];
2264                        j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2265                    int iRow = row[j];
2266                    double value = element[j];
2267                    rowActivity[iRow] += distance * value;
2268                }
2269            }
2270        }
2271        // and now all if improving
2272        double lastChange = penaltyChange ? 1.0 : 0.0;
2273        int numberPasses=0;
2274        while (lastChange > 1.0e-2 && numberPasses < 1000) {
2275            lastChange = 0;
2276            numberPasses++;
2277            for (iColumn = 0; iColumn < numberColumns; iColumn++) {
2278                bool isInteger = (integerType[iColumn] != 0);
2279                double currentValue = newSolution[iColumn];
2280                double lowerValue = lower[iColumn];
2281                double upperValue = upper[iColumn];
2282                int j;
2283                int anyBadDown = 0;
2284                int anyBadUp = 0;
2285                double upImprovement = 0.0;
2286                double downImprovement = 0.0;
2287                for (j = columnStart[iColumn];
2288                        j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2289                    int iRow = row[j];
2290                    double value = element[j];
2291                    if (isInteger) {
2292                        if (value > 0.0) {
2293                            if (rowActivity[iRow] + value > rowUpper[iRow] + primalTolerance)
2294                                anyBadUp++;
2295                            if (rowActivity[iRow] - value < rowLower[iRow] - primalTolerance)
2296                                anyBadDown++;
2297                        } else {
2298                            if (rowActivity[iRow] - value > rowUpper[iRow] + primalTolerance)
2299                                anyBadDown++;
2300                            if (rowActivity[iRow] + value < rowLower[iRow] - primalTolerance)
2301                                anyBadUp++;
2302                        }
2303                    }
2304                    if (rowUpper[iRow] > rowLower[iRow]) {
2305                        if (rowActivity[iRow] > rowUpper[iRow] + primalTolerance) {
2306                            // infeasible above
2307                            downImprovement += value;
2308                            upImprovement -= value;
2309                            if (value > 0.0)
2310                                anyBadUp++;
2311                            else
2312                                anyBadDown++;
2313                        } else if (rowActivity[iRow] > rowUpper[iRow] - primalTolerance) {
2314                            // feasible at ub
2315                            if (value > 0.0) {
2316                                upImprovement -= value;
2317                                anyBadUp++;
2318                            } else {
2319                                downImprovement += value;
2320                                anyBadDown++;
2321                            }
2322                        } else if (rowActivity[iRow] > rowLower[iRow] + primalTolerance) {
2323                            // feasible in interior
2324                        } else if (rowActivity[iRow] > rowLower[iRow] - primalTolerance) {
2325                            // feasible at lb
2326                            if (value < 0.0) {
2327                                upImprovement += value;
2328                                anyBadUp++;
2329                            } else {
2330                                downImprovement -= value;
2331                                anyBadDown++;
2332                            }
2333                        } else {
2334                            // infeasible below
2335                            downImprovement -= value;
2336                            upImprovement += value;
2337                            if (value < 0.0)
2338                                anyBadUp++;
2339                            else
2340                                anyBadDown++;
2341                        }
2342                    } else {
2343                        // equality row
2344                        if (rowActivity[iRow] > rowUpper[iRow] + primalTolerance) {
2345                            // infeasible above
2346                            downImprovement += value;
2347                            upImprovement -= value;
2348                            if (value > 0.0)
2349                                anyBadUp++;
2350                            else
2351                                anyBadDown++;
2352                        } else if (rowActivity[iRow] < rowLower[iRow] - primalTolerance) {
2353                            // infeasible below
2354                            downImprovement -= value;
2355                            upImprovement += value;
2356                            if (value < 0.0)
2357                                anyBadUp++;
2358                            else
2359                                anyBadDown++;
2360                        } else {
2361                            // feasible - no good
2362                            anyBadUp = -1;
2363                            anyBadDown = -1;
2364                            break;
2365                        }
2366                    }
2367                }
2368                // could change tests for anyBad
2369                if (anyBadUp)
2370                    upImprovement = 0.0;
2371                if (anyBadDown)
2372                    downImprovement = 0.0;
2373                double way = 0.0;
2374                double improvement = 0.0;
2375                if (downImprovement > 0.0 && currentValue > lowerValue) {
2376                    way = -1.0;
2377                    improvement = downImprovement;
2378                } else if (upImprovement > 0.0 && currentValue < upperValue) {
2379                    way = 1.0;
2380                    improvement = upImprovement;
2381                }
2382                if (way) {
2383                    // can improve
2384                    double distance = COIN_DBL_MAX;
2385                    for (j = columnStart[iColumn];
2386                            j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2387                        int iRow = row[j];
2388                        double value = element[j] * way;
2389                        if (rowActivity[iRow] > rowUpper[iRow] + primalTolerance) {
2390                            // infeasible above
2391                            assert (value < 0.0);
2392                            double gap = rowActivity[iRow] - rowUpper[iRow];
2393                            if (gap + value*distance < 0.0) {
2394                                // If integer then has to move by 1
2395                                if (!isInteger)
2396                                    distance = -gap / value;
2397                                else
2398                                    distance = CoinMax(-gap / value, 1.0);
2399                            }
2400                        } else if (rowActivity[iRow] < rowLower[iRow] - primalTolerance) {
2401                            // infeasible below
2402                            assert (value > 0.0);
2403                            double gap = rowActivity[iRow] - rowLower[iRow];
2404                            if (gap + value*distance > 0.0) {
2405                                // If integer then has to move by 1
2406                                if (!isInteger)
2407                                    distance = -gap / value;
2408                                else
2409                                    distance = CoinMax(-gap / value, 1.0);
2410                            }
2411                        } else {
2412                            // feasible
2413                            if (value > 0) {
2414                                double gap = rowActivity[iRow] - rowUpper[iRow];
2415                                if (gap + value*distance > 0.0)
2416                                    distance = -gap / value;
2417                            } else {
2418                                double gap = rowActivity[iRow] - rowLower[iRow];
2419                                if (gap + value*distance < 0.0)
2420                                    distance = -gap / value;
2421                            }
2422                        }
2423                    }
2424                    if (isInteger)
2425                        distance = floor(distance + 1.05e-8);
2426                    if (!distance) {
2427                        // should never happen
2428                        //printf("zero distance in CbcRounding - debug\n");
2429                    }
2430                    //move
2431                    lastChange += improvement * distance;
2432                    distance *= way;
2433                    newSolution[iColumn] += distance;
2434                    newSolutionValue += direction * objective[iColumn] * distance;
2435                    for (j = columnStart[iColumn];
2436                            j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2437                        int iRow = row[j];
2438                        double value = element[j];
2439                        rowActivity[iRow] += distance * value;
2440                    }
2441                }
2442            }
2443            penaltyChange += lastChange;
2444        }
2445        penalty -= penaltyChange;
2446        if (penalty < 1.0e-5*fabs(penaltyChange)) {
2447            // recompute
2448            penalty = 0.0;
2449            for (i = 0; i < numberRows; i++) {
2450                double value = rowActivity[i];
2451                if (value < rowLower[i] - primalTolerance)
2452                    penalty += rowLower[i] - value;
2453                else if (value > rowUpper[i] + primalTolerance)
2454                    penalty += value - rowUpper[i];
2455            }
2456        }
2457    }
2458
2459    // Could also set SOS (using random) and repeat
2460    if (!penalty) {
2461        // See if we can do better
2462        //seed_++;
2463        //CoinSeedRandom(seed_);
2464        // Random number between 0 and 1.
2465        double randomNumber = randomNumberGenerator_.randomDouble();
2466        int iPass;
2467        int start[2];
2468        int end[2];
2469        int iRandom = static_cast<int> (randomNumber * (static_cast<double> (numberIntegers)));
2470        start[0] = iRandom;
2471        end[0] = numberIntegers;
2472        start[1] = 0;
2473        end[1] = iRandom;
2474        for (iPass = 0; iPass < 2; iPass++) {
2475            int i;
2476            for (i = start[iPass]; i < end[iPass]; i++) {
2477                int iColumn = integerVariable[i];
2478#ifndef NDEBUG
2479                double value = newSolution[iColumn];
2480                assert (fabs(floor(value + 0.5) - value) < integerTolerance);
2481#endif
2482                double cost = direction * objective[iColumn];
2483                double move = 0.0;
2484                if (cost > 0.0)
2485                    move = -1.0;
2486                else if (cost < 0.0)
2487                    move = 1.0;
2488                while (move) {
2489                    bool good = true;
2490                    double newValue = newSolution[iColumn] + move;
2491                    if (newValue < lower[iColumn] - useTolerance ||
2492                            newValue > upper[iColumn] + useTolerance) {
2493                        move = 0.0;
2494                    } else {
2495                        // see if we can move
2496                        int j;
2497                        for (j = columnStart[iColumn];
2498                                j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2499                            int iRow = row[j];
2500                            double newActivity = rowActivity[iRow] + move * element[j];
2501                            if (newActivity < rowLower[iRow] - primalTolerance ||
2502                                    newActivity > rowUpper[iRow] + primalTolerance) {
2503                                good = false;
2504                                break;
2505                            }
2506                        }
2507                        if (good) {
2508                            newSolution[iColumn] = newValue;
2509                            newSolutionValue += move * cost;
2510                            int j;
2511                            for (j = columnStart[iColumn];
2512                                    j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2513                                int iRow = row[j];
2514                                rowActivity[iRow] += move * element[j];
2515                            }
2516                        } else {
2517                            move = 0.0;
2518                        }
2519                    }
2520                }
2521            }
2522        }
2523        // Just in case of some stupidity
2524        double objOffset = 0.0;
2525        solver->getDblParam(OsiObjOffset, objOffset);
2526        newSolutionValue = -objOffset;
2527        for ( i = 0 ; i < numberColumns ; i++ )
2528            newSolutionValue += objective[i] * newSolution[i];
2529        newSolutionValue *= direction;
2530        //printf("new solution value %g %g\n",newSolutionValue,solutionValue);
2531        if (newSolutionValue < solutionValue) {
2532            // paranoid check
2533            memset(rowActivity, 0, numberRows*sizeof(double));
2534            for (i = 0; i < numberColumns; i++) {
2535                int j;
2536                double value = newSolution[i];
2537                if (value) {
2538                    for (j = columnStart[i];
2539                            j < columnStart[i] + columnLength[i]; j++) {
2540                        int iRow = row[j];
2541                        rowActivity[iRow] += value * element[j];
2542                    }
2543                }
2544            }
2545            // check was approximately feasible
2546            bool feasible = true;
2547            for (i = 0; i < numberRows; i++) {
2548                if (rowActivity[i] < rowLower[i]) {
2549                    if (rowActivity[i] < rowLower[i] - 1000.0*primalTolerance)
2550                        feasible = false;
2551                } else if (rowActivity[i] > rowUpper[i]) {
2552                    if (rowActivity[i] > rowUpper[i] + 1000.0*primalTolerance)
2553                        feasible = false;
2554                }
2555            }
2556            if (feasible) {
2557                // new solution
2558                memcpy(betterSolution, newSolution, numberColumns*sizeof(double));
2559                solutionValue = newSolutionValue;
2560                //printf("** Solution of %g found by rounding\n",newSolutionValue);
2561                returnCode = 1;
2562            } else {
2563                // Can easily happen
2564                //printf("Debug CbcRounding giving bad solution\n");
2565            }
2566        }
2567    }
2568#ifdef NEW_ROUNDING
2569    if (!returnCode) {
2570#ifdef JJF_ZERO
2571        // back to starting point
2572        memcpy(newSolution, solution, numberColumns*sizeof(double));
2573        memset(rowActivity, 0, numberRows*sizeof(double));
2574        for (i = 0; i < numberColumns; i++) {
2575            int j;
2576            double value = newSolution[i];
2577            if (value < lower[i]) {
2578                value = lower[i];
2579                newSolution[i] = value;
2580            } else if (value > upper[i]) {
2581                value = upper[i];
2582                newSolution[i] = value;
2583            }
2584            if (value) {
2585                for (j = columnStart[i];
2586                        j < columnStart[i] + columnLength[i]; j++) {
2587                    int iRow = row[j];
2588                    rowActivity[iRow] += value * element[j];
2589                }
2590            }
2591        }
2592        // check was feasible - if not adjust (cleaning may move)
2593        for (i = 0; i < numberRows; i++) {
2594            if (rowActivity[i] < rowLower[i]) {
2595                //assert (rowActivity[i]>rowLower[i]-1000.0*primalTolerance);
2596                rowActivity[i] = rowLower[i];
2597            } else if (rowActivity[i] > rowUpper[i]) {
2598                //assert (rowActivity[i]<rowUpper[i]+1000.0*primalTolerance);
2599                rowActivity[i] = rowUpper[i];
2600            }
2601        }
2602#endif
2603        int * candidate = new int [numberColumns];
2604        int nCandidate = 0;
2605        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
2606            bool isInteger = (integerType[iColumn] != 0);
2607            if (isInteger) {
2608                double currentValue = newSolution[iColumn];
2609                if (fabs(currentValue - floor(currentValue + 0.5)) > 1.0e-8)
2610                    candidate[nCandidate++] = iColumn;
2611            }
2612        }
2613        if (true) {
2614            // Rounding as in Berthold
2615            while (nCandidate) {
2616                double infeasibility = 1.0e-7;
2617                int iRow = -1;
2618                for (i = 0; i < numberRows; i++) {
2619                    double value = 0.0;
2620                    if (rowActivity[i] < rowLower[i]) {
2621                        value = rowLower[i] - rowActivity[i];
2622                    } else if (rowActivity[i] > rowUpper[i]) {
2623                        value = rowActivity[i] - rowUpper[i];
2624                    }
2625                    if (value > infeasibility) {
2626                        infeasibility = value;
2627                        iRow = i;
2628                    }
2629                }
2630                if (iRow >= 0) {
2631                    // infeasible
2632                } else {
2633                    // feasible
2634                }
2635            }
2636        } else {
2637            // Shifting as in Berthold
2638        }
2639        delete [] candidate;
2640    }
2641#endif
2642    delete [] newSolution;
2643    delete [] rowActivity;
2644    return returnCode;
2645}
2646// update model
2647void CbcRounding::setModel(CbcModel * model)
2648{
2649    model_ = model;
2650    // Get a copy of original matrix (and by row for rounding);
2651    assert(model_->solver());
2652    if (model_->solver()->getNumRows()) {
2653        matrix_ = *model_->solver()->getMatrixByCol();
2654        matrixByRow_ = *model_->solver()->getMatrixByRow();
2655        // make sure model okay for heuristic
2656        validate();
2657    }
2658}
2659// Validate model i.e. sets when_ to 0 if necessary (may be NULL)
2660void
2661CbcRounding::validate()
2662{
2663    if (model_ && (when() % 100) < 10) {
2664        if (model_->numberIntegers() !=
2665                model_->numberObjects() && (model_->numberObjects() ||
2666                                            (model_->specialOptions()&1024) == 0)) {
2667            int numberOdd = 0;
2668            for (int i = 0; i < model_->numberObjects(); i++) {
2669                if (!model_->object(i)->canDoHeuristics())
2670                    numberOdd++;
2671            }
2672            if (numberOdd)
2673                setWhen(0);
2674        }
2675    }
2676#ifdef NEW_ROUNDING
2677    int numberColumns = matrix_.getNumCols();
2678    down_ = new unsigned short [numberColumns];
2679    up_ = new unsigned short [numberColumns];
2680    equal_ = new unsigned short [numberColumns];
2681    // Column copy
2682    const double * element = matrix_.getElements();
2683    const int * row = matrix_.getIndices();
2684    const CoinBigIndex * columnStart = matrix_.getVectorStarts();
2685    const int * columnLength = matrix_.getVectorLengths();
2686    const double * rowLower = model.solver()->getRowLower();
2687    const double * rowUpper = model.solver()->getRowUpper();
2688    for (int i = 0; i < numberColumns; i++) {
2689        int down = 0;
2690        int up = 0;
2691        int equal = 0;
2692        if (columnLength[i] > 65535) {
2693            equal[0] = 65535;
2694            break; // unlikely to work
2695        }
2696        for (CoinBigIndex j = columnStart[i];
2697                j < columnStart[i] + columnLength[i]; j++) {
2698            int iRow = row[j];
2699            if (rowLower[iRow] > -1.0e20 && rowUpper[iRow] < 1.0e20) {
2700                equal++;
2701            } else if (element[j] > 0.0) {
2702                if (rowUpper[iRow] < 1.0e20)
2703                    up++;
2704                else
2705                    down--;
2706            } else {
2707                if (rowLower[iRow] > -1.0e20)
2708                    up++;
2709                else
2710                    down--;
2711            }
2712        }
2713        down_[i] = (unsigned short) down;
2714        up_[i] = (unsigned short) up;
2715        equal_[i] = (unsigned short) equal;
2716    }
2717#else
2718    down_ = NULL;
2719    up_ = NULL;
2720    equal_ = NULL;
2721#endif
2722}
2723
2724// Default Constructor
2725CbcHeuristicPartial::CbcHeuristicPartial()
2726        : CbcHeuristic()
2727{
2728    fixPriority_ = 10000;
2729}
2730
2731// Constructor from model
2732CbcHeuristicPartial::CbcHeuristicPartial(CbcModel & model, int fixPriority, int numberNodes)
2733        : CbcHeuristic(model)
2734{
2735    fixPriority_ = fixPriority;
2736    setNumberNodes(numberNodes);
2737    validate();
2738}
2739
2740// Destructor
2741CbcHeuristicPartial::~CbcHeuristicPartial ()
2742{
2743}
2744
2745// Clone
2746CbcHeuristic *
2747CbcHeuristicPartial::clone() const
2748{
2749    return new CbcHeuristicPartial(*this);
2750}
2751// Create C++ lines to get to current state
2752void
2753CbcHeuristicPartial::generateCpp( FILE * fp)
2754{
2755    CbcHeuristicPartial other;
2756    fprintf(fp, "0#include \"CbcHeuristic.hpp\"\n");
2757    fprintf(fp, "3  CbcHeuristicPartial partial(*cbcModel);\n");
2758    CbcHeuristic::generateCpp(fp, "partial");
2759    if (fixPriority_ != other.fixPriority_)
2760        fprintf(fp, "3  partial.setFixPriority(%d);\n", fixPriority_);
2761    else
2762        fprintf(fp, "4  partial.setFixPriority(%d);\n", fixPriority_);
2763    fprintf(fp, "3  cbcModel->addHeuristic(&partial);\n");
2764}
2765//#define NEW_PARTIAL
2766// Copy constructor
2767CbcHeuristicPartial::CbcHeuristicPartial(const CbcHeuristicPartial & rhs)
2768        :
2769        CbcHeuristic(rhs),
2770        fixPriority_(rhs.fixPriority_)
2771{
2772}
2773
2774// Assignment operator
2775CbcHeuristicPartial &
2776CbcHeuristicPartial::operator=( const CbcHeuristicPartial & rhs)
2777{
2778    if (this != &rhs) {
2779        CbcHeuristic::operator=(rhs);
2780        fixPriority_ = rhs.fixPriority_;
2781    }
2782    return *this;
2783}
2784
2785// Resets stuff if model changes
2786void
2787CbcHeuristicPartial::resetModel(CbcModel * model)
2788{
2789    model_ = model;
2790    // Get a copy of original matrix (and by row for partial);
2791    assert(model_->solver());
2792    validate();
2793}
2794// See if partial will give solution
2795// Sets value of solution
2796// Assumes rhs for original matrix still okay
2797// At present only works with integers
2798// Fix values if asked for
2799// Returns 1 if solution, 0 if not
2800int
2801CbcHeuristicPartial::solution(double & solutionValue,
2802                              double * betterSolution)
2803{
2804    // Return if already done
2805    if (fixPriority_ < 0)
2806        return 0; // switched off
2807#ifdef HEURISTIC_INFORM
2808    printf("Entering heuristic %s - nRuns %d numCould %d when %d\n",
2809           heuristicName(),numRuns_,numCouldRun_,when_);
2810#endif
2811    const double * hotstartSolution = model_->hotstartSolution();
2812    const int * hotstartPriorities = model_->hotstartPriorities();
2813    if (!hotstartSolution)
2814        return 0;
2815    OsiSolverInterface * solver = model_->solver();
2816
2817    int numberIntegers = model_->numberIntegers();
2818    const int * integerVariable = model_->integerVariable();
2819
2820    OsiSolverInterface * newSolver = model_->continuousSolver()->clone();
2821    const double * colLower = newSolver->getColLower();
2822    const double * colUpper = newSolver->getColUpper();
2823
2824    double primalTolerance;
2825    solver->getDblParam(OsiPrimalTolerance, primalTolerance);
2826
2827    int i;
2828    int numberFixed = 0;
2829    int returnCode = 0;
2830
2831    for (i = 0; i < numberIntegers; i++) {
2832        int iColumn = integerVariable[i];
2833        if (abs(hotstartPriorities[iColumn]) <= fixPriority_) {
2834            double value = hotstartSolution[iColumn];
2835            double lower = colLower[iColumn];
2836            double upper = colUpper[iColumn];
2837            value = CoinMax(value, lower);
2838            value = CoinMin(value, upper);
2839            if (fabs(value - floor(value + 0.5)) < 1.0e-8) {
2840                value = floor(value + 0.5);
2841                newSolver->setColLower(iColumn, value);
2842                newSolver->setColUpper(iColumn, value);
2843                numberFixed++;
2844            }
2845        }
2846    }
2847    if (numberFixed > numberIntegers / 5 - 100000000) {
2848#ifdef COIN_DEVELOP
2849        printf("%d integers fixed\n", numberFixed);
2850#endif
2851        returnCode = smallBranchAndBound(newSolver, numberNodes_, betterSolution, solutionValue,
2852                                         model_->getCutoff(), "CbcHeuristicPartial");
2853        if (returnCode < 0)
2854            returnCode = 0; // returned on size
2855        //printf("return code %d",returnCode);
2856        if ((returnCode&2) != 0) {
2857            // could add cut
2858            returnCode &= ~2;
2859            //printf("could add cut with %d elements (if all 0-1)\n",nFix);
2860        } else {
2861            //printf("\n");
2862        }
2863    }
2864    fixPriority_ = -1; // switch off
2865
2866    delete newSolver;
2867    return returnCode;
2868}
2869// update model
2870void CbcHeuristicPartial::setModel(CbcModel * model)
2871{
2872    model_ = model;
2873    assert(model_->solver());
2874    // make sure model okay for heuristic
2875    validate();
2876}
2877// Validate model i.e. sets when_ to 0 if necessary (may be NULL)
2878void
2879CbcHeuristicPartial::validate()
2880{
2881    if (model_ && (when() % 100) < 10) {
2882        if (model_->numberIntegers() !=
2883                model_->numberObjects())
2884            setWhen(0);
2885    }
2886}
2887bool
2888CbcHeuristicPartial::shouldHeurRun(int /*whereFrom*/)
2889{
2890    return true;
2891}
2892
2893// Default Constructor
2894CbcSerendipity::CbcSerendipity()
2895        : CbcHeuristic()
2896{
2897}
2898
2899// Constructor from model
2900CbcSerendipity::CbcSerendipity(CbcModel & model)
2901        : CbcHeuristic(model)
2902{
2903}
2904
2905// Destructor
2906CbcSerendipity::~CbcSerendipity ()
2907{
2908}
2909
2910// Clone
2911CbcHeuristic *
2912CbcSerendipity::clone() const
2913{
2914    return new CbcSerendipity(*this);
2915}
2916// Create C++ lines to get to current state
2917void
2918CbcSerendipity::generateCpp( FILE * fp)
2919{
2920    fprintf(fp, "0#include \"CbcHeuristic.hpp\"\n");
2921    fprintf(fp, "3  CbcSerendipity serendipity(*cbcModel);\n");
2922    CbcHeuristic::generateCpp(fp, "serendipity");
2923    fprintf(fp, "3  cbcModel->addHeuristic(&serendipity);\n");
2924}
2925
2926// Copy constructor
2927CbcSerendipity::CbcSerendipity(const CbcSerendipity & rhs)
2928        :
2929        CbcHeuristic(rhs)
2930{
2931}
2932
2933// Assignment operator
2934CbcSerendipity &
2935CbcSerendipity::operator=( const CbcSerendipity & rhs)
2936{
2937    if (this != &rhs) {
2938        CbcHeuristic::operator=(rhs);
2939    }
2940    return *this;
2941}
2942
2943// Returns 1 if solution, 0 if not
2944int
2945CbcSerendipity::solution(double & solutionValue,
2946                         double * betterSolution)
2947{
2948    if (!model_)
2949        return 0;
2950#ifdef HEURISTIC_INFORM
2951    printf("Entering heuristic %s - nRuns %d numCould %d when %d\n",
2952           heuristicName(),numRuns_,numCouldRun_,when_);
2953#endif
2954    if (!inputSolution_) {
2955        // get information on solver type
2956        OsiAuxInfo * auxInfo = model_->solver()->getAuxiliaryInfo();
2957        OsiBabSolver * auxiliaryInfo = dynamic_cast< OsiBabSolver *> (auxInfo);
2958        if (auxiliaryInfo) {
2959            return auxiliaryInfo->solution(solutionValue, betterSolution, model_->solver()->getNumCols());
2960        } else {
2961            return 0;
2962        }
2963    } else {
2964        int numberColumns = model_->getNumCols();
2965        double value = inputSolution_[numberColumns];
2966        int returnCode = 0;
2967        if (value < solutionValue) {
2968            solutionValue = value;
2969            memcpy(betterSolution, inputSolution_, numberColumns*sizeof(double));
2970            returnCode = 1;
2971        }
2972        delete [] inputSolution_;
2973        inputSolution_ = NULL;
2974        model_ = NULL; // switch off
2975        return returnCode;
2976    }
2977}
2978// update model
2979void CbcSerendipity::setModel(CbcModel * model)
2980{
2981    model_ = model;
2982}
2983// Resets stuff if model changes
2984void
2985CbcSerendipity::resetModel(CbcModel * model)
2986{
2987    model_ = model;
2988}
2989
2990
2991// Default Constructor
2992CbcHeuristicJustOne::CbcHeuristicJustOne()
2993        : CbcHeuristic(),
2994        probabilities_(NULL),
2995        heuristic_(NULL),
2996        numberHeuristics_(0)
2997{
2998}
2999
3000// Constructor from model
3001CbcHeuristicJustOne::CbcHeuristicJustOne(CbcModel & model)
3002        : CbcHeuristic(model),
3003        probabilities_(NULL),
3004        heuristic_(NULL),
3005        numberHeuristics_(0)
3006{
3007}
3008
3009// Destructor
3010CbcHeuristicJustOne::~CbcHeuristicJustOne ()
3011{
3012    for (int i = 0; i < numberHeuristics_; i++)
3013        delete heuristic_[i];
3014    delete [] heuristic_;
3015    delete [] probabilities_;
3016}
3017
3018// Clone
3019CbcHeuristicJustOne *
3020CbcHeuristicJustOne::clone() const
3021{
3022    return new CbcHeuristicJustOne(*this);
3023}
3024
3025// Create C++ lines to get to current state
3026void
3027CbcHeuristicJustOne::generateCpp( FILE * fp)
3028{
3029    CbcHeuristicJustOne other;
3030    fprintf(fp, "0#include \"CbcHeuristicJustOne.hpp\"\n");
3031    fprintf(fp, "3  CbcHeuristicJustOne heuristicJustOne(*cbcModel);\n");
3032    CbcHeuristic::generateCpp(fp, "heuristicJustOne");
3033    fprintf(fp, "3  cbcModel->addHeuristic(&heuristicJustOne);\n");
3034}
3035
3036// Copy constructor
3037CbcHeuristicJustOne::CbcHeuristicJustOne(const CbcHeuristicJustOne & rhs)
3038        :
3039        CbcHeuristic(rhs),
3040        probabilities_(NULL),
3041        heuristic_(NULL),
3042        numberHeuristics_(rhs.numberHeuristics_)
3043{
3044    if (numberHeuristics_) {
3045        probabilities_ = CoinCopyOfArray(rhs.probabilities_, numberHeuristics_);
3046        heuristic_ = new CbcHeuristic * [numberHeuristics_];
3047        for (int i = 0; i < numberHeuristics_; i++)
3048            heuristic_[i] = rhs.heuristic_[i]->clone();
3049    }
3050}
3051
3052// Assignment operator
3053CbcHeuristicJustOne &
3054CbcHeuristicJustOne::operator=( const CbcHeuristicJustOne & rhs)
3055{
3056    if (this != &rhs) {
3057        CbcHeuristic::operator=(rhs);
3058        for (int i = 0; i < numberHeuristics_; i++)
3059            delete heuristic_[i];
3060        delete [] heuristic_;
3061        delete [] probabilities_;
3062        probabilities_ = NULL;
3063        heuristic_ = NULL;
3064        numberHeuristics_ = rhs.numberHeuristics_;
3065        if (numberHeuristics_) {
3066            probabilities_ = CoinCopyOfArray(rhs.probabilities_, numberHeuristics_);
3067            heuristic_ = new CbcHeuristic * [numberHeuristics_];
3068            for (int i = 0; i < numberHeuristics_; i++)
3069                heuristic_[i] = rhs.heuristic_[i]->clone();
3070        }
3071    }
3072    return *this;
3073}
3074// Sets value of solution
3075// Returns 1 if solution, 0 if not
3076int
3077CbcHeuristicJustOne::solution(double & solutionValue,
3078                              double * betterSolution)
3079{
3080#ifdef DIVE_DEBUG
3081    std::cout << "solutionValue = " << solutionValue << std::endl;
3082#endif
3083    ++numCouldRun_;
3084
3085    // test if the heuristic can run
3086    if (!shouldHeurRun_randomChoice() || !numberHeuristics_)
3087        return 0;
3088    double randomNumber = randomNumberGenerator_.randomDouble();
3089    int i;
3090    for (i = 0; i < numberHeuristics_; i++) {
3091        if (randomNumber < probabilities_[i])
3092            break;
3093    }
3094    assert (i < numberHeuristics_);
3095    int returnCode;
3096    //model_->unsetDivingHasRun();
3097#ifdef COIN_DEVELOP
3098    printf("JustOne running %s\n",
3099           heuristic_[i]->heuristicName());
3100#endif
3101    returnCode = heuristic_[i]->solution(solutionValue, betterSolution);
3102#ifdef COIN_DEVELOP
3103    if (returnCode)
3104        printf("JustOne running %s found solution\n",
3105               heuristic_[i]->heuristicName());
3106#endif
3107    return returnCode;
3108}
3109// Resets stuff if model changes
3110void
3111CbcHeuristicJustOne::resetModel(CbcModel * model)
3112{
3113    CbcHeuristic::resetModel(model);
3114    for (int i = 0; i < numberHeuristics_; i++)
3115        heuristic_[i]->resetModel(model);
3116}
3117// update model (This is needed if cliques update matrix etc)
3118void
3119CbcHeuristicJustOne::setModel(CbcModel * model)
3120{
3121    CbcHeuristic::setModel(model);
3122    for (int i = 0; i < numberHeuristics_; i++)
3123        heuristic_[i]->setModel(model);
3124}
3125// Validate model i.e. sets when_ to 0 if necessary (may be NULL)
3126void
3127CbcHeuristicJustOne::validate()
3128{
3129    CbcHeuristic::validate();
3130    for (int i = 0; i < numberHeuristics_; i++)
3131        heuristic_[i]->validate();
3132}
3133// Adds an heuristic with probability
3134void
3135CbcHeuristicJustOne::addHeuristic(const CbcHeuristic * heuristic, double probability)
3136{
3137    CbcHeuristic * thisOne = heuristic->clone();
3138    thisOne->setWhen(-999);
3139    CbcHeuristic ** tempH = CoinCopyOfArrayPartial(heuristic_, numberHeuristics_ + 1,
3140                            numberHeuristics_);
3141    delete [] heuristic_;
3142    heuristic_ = tempH;
3143    heuristic_[numberHeuristics_] = thisOne;
3144    double * tempP = CoinCopyOfArrayPartial(probabilities_, numberHeuristics_ + 1,
3145                                            numberHeuristics_);
3146    delete [] probabilities_;
3147    probabilities_ = tempP;
3148    probabilities_[numberHeuristics_] = probability;
3149    numberHeuristics_++;
3150}
3151// Normalize probabilities
3152void
3153CbcHeuristicJustOne::normalizeProbabilities()
3154{
3155    double sum = 0.0;
3156    for (int i = 0; i < numberHeuristics_; i++)
3157        sum += probabilities_[i];
3158    double multiplier = 1.0 / sum;
3159    sum = 0.0;
3160    for (int i = 0; i < numberHeuristics_; i++) {
3161        sum += probabilities_[i];
3162        probabilities_[i] = sum * multiplier;
3163    }
3164    assert (fabs(probabilities_[numberHeuristics_-1] - 1.0) < 1.0e-5);
3165    probabilities_[numberHeuristics_-1] = 1.000001;
3166}
3167
Note: See TracBrowser for help on using the repository browser.