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

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

add random seed setting

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