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

Last change on this file since 1883 was 1883, checked in by stefan, 7 years ago

sync with trunk rev 1882

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