source: branches/sandbox/Cbc/src/CbcHeuristic.cpp @ 1286

Last change on this file since 1286 was 1286, checked in by EdwinStraver, 10 years ago

Changed formatting using AStyle -A4 -p

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