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

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

coding for getting postprocessed model

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