source: stable/2.8/Cbc/src/CbcHeuristic.cpp @ 1865

Last change on this file since 1865 was 1865, checked in by forrest, 6 years ago

cuts off when user asks for them off

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