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

Last change on this file since 1802 was 1802, checked in by forrest, 7 years ago

add Proximity heuristic (Fischetti and Monaci) - shouldn't break anything

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