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

Last change on this file since 2280 was 2280, checked in by forrest, 3 years ago

allow heuristics to see if integers are 'optional'

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