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

Last change on this file since 2080 was 2040, checked in by forrest, 5 years ago

fixes for odd SOS

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