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

Last change on this file since 1432 was 1432, checked in by bjarni, 11 years ago

Added extra return at end of each source file where needed, to remove possible linefeed conflicts (NightlyBuild? errors)

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