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

Last change on this file since 1942 was 1888, checked in by stefan, 6 years ago

sync with trunk rev 1887

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