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

Last change on this file since 1952 was 1952, checked in by tkr, 6 years ago

Adding log level check before printing message

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