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

Last change on this file since 1570 was 1570, checked in by forrest, 8 years ago

modifications for printing

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