source: trunk/Cbc/src/CbcHeuristicFPump.cpp @ 1871

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

fix for mini b&b success

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 132.8 KB
Line 
1/* $Id: CbcHeuristicFPump.cpp 1871 2013-02-09 11:55:54Z forrest $ */
2// Copyright (C) 2004, 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#include <cassert>
11#include <cstdlib>
12#include <cmath>
13#include <cfloat>
14
15#include "OsiSolverInterface.hpp"
16#include "CbcModel.hpp"
17#ifdef COIN_HAS_CLP
18#include "OsiClpSolverInterface.hpp"
19#endif
20#include "CbcMessage.hpp"
21#include "CbcHeuristicFPump.hpp"
22#include "CbcBranchActual.hpp"
23#include "CbcBranchDynamic.hpp"
24#include "CoinHelperFunctions.hpp"
25#include "CoinWarmStartBasis.hpp"
26#include "CoinTime.hpp"
27#include "CbcEventHandler.hpp"
28
29
30// Default Constructor
31CbcHeuristicFPump::CbcHeuristicFPump()
32        : CbcHeuristic(),
33        startTime_(0.0),
34        maximumTime_(0.0),
35        fakeCutoff_(COIN_DBL_MAX),
36        absoluteIncrement_(0.0),
37        relativeIncrement_(0.0),
38        defaultRounding_(0.49999),
39        initialWeight_(0.0),
40        weightFactor_(0.1),
41        artificialCost_(COIN_DBL_MAX),
42        iterationRatio_(0.0),
43        reducedCostMultiplier_(1.0),
44        maximumPasses_(100),
45        maximumRetries_(1),
46        accumulate_(0),
47        fixOnReducedCosts_(1),
48        roundExpensive_(false)
49{
50    setWhen(1);
51}
52
53// Constructor from model
54CbcHeuristicFPump::CbcHeuristicFPump(CbcModel & model,
55                                     double downValue, bool roundExpensive)
56        : CbcHeuristic(model),
57        startTime_(0.0),
58        maximumTime_(0.0),
59        fakeCutoff_(COIN_DBL_MAX),
60        absoluteIncrement_(0.0),
61        relativeIncrement_(0.0),
62        defaultRounding_(downValue),
63        initialWeight_(0.0),
64        weightFactor_(0.1),
65        artificialCost_(COIN_DBL_MAX),
66        iterationRatio_(0.0),
67        reducedCostMultiplier_(1.0),
68        maximumPasses_(100),
69        maximumRetries_(1),
70        accumulate_(0),
71        fixOnReducedCosts_(1),
72        roundExpensive_(roundExpensive)
73{
74    setWhen(1);
75}
76
77// Destructor
78CbcHeuristicFPump::~CbcHeuristicFPump ()
79{
80}
81
82// Clone
83CbcHeuristic *
84CbcHeuristicFPump::clone() const
85{
86    return new CbcHeuristicFPump(*this);
87}
88// Create C++ lines to get to current state
89void
90CbcHeuristicFPump::generateCpp( FILE * fp)
91{
92    CbcHeuristicFPump other;
93    fprintf(fp, "0#include \"CbcHeuristicFPump.hpp\"\n");
94    fprintf(fp, "3  CbcHeuristicFPump heuristicFPump(*cbcModel);\n");
95    CbcHeuristic::generateCpp(fp, "heuristicFPump");
96    if (maximumPasses_ != other.maximumPasses_)
97        fprintf(fp, "3  heuristicFPump.setMaximumPasses(%d);\n", maximumPasses_);
98    else
99        fprintf(fp, "4  heuristicFPump.setMaximumPasses(%d);\n", maximumPasses_);
100    if (maximumRetries_ != other.maximumRetries_)
101        fprintf(fp, "3  heuristicFPump.setMaximumRetries(%d);\n", maximumRetries_);
102    else
103        fprintf(fp, "4  heuristicFPump.setMaximumRetries(%d);\n", maximumRetries_);
104    if (accumulate_ != other.accumulate_)
105        fprintf(fp, "3  heuristicFPump.setAccumulate(%d);\n", accumulate_);
106    else
107        fprintf(fp, "4  heuristicFPump.setAccumulate(%d);\n", accumulate_);
108    if (fixOnReducedCosts_ != other.fixOnReducedCosts_)
109        fprintf(fp, "3  heuristicFPump.setFixOnReducedCosts(%d);\n", fixOnReducedCosts_);
110    else
111        fprintf(fp, "4  heuristicFPump.setFixOnReducedCosts(%d);\n", fixOnReducedCosts_);
112    if (maximumTime_ != other.maximumTime_)
113        fprintf(fp, "3  heuristicFPump.setMaximumTime(%g);\n", maximumTime_);
114    else
115        fprintf(fp, "4  heuristicFPump.setMaximumTime(%g);\n", maximumTime_);
116    if (fakeCutoff_ != other.fakeCutoff_)
117        fprintf(fp, "3  heuristicFPump.setFakeCutoff(%g);\n", fakeCutoff_);
118    else
119        fprintf(fp, "4  heuristicFPump.setFakeCutoff(%g);\n", fakeCutoff_);
120    if (absoluteIncrement_ != other.absoluteIncrement_)
121        fprintf(fp, "3  heuristicFPump.setAbsoluteIncrement(%g);\n", absoluteIncrement_);
122    else
123        fprintf(fp, "4  heuristicFPump.setAbsoluteIncrement(%g);\n", absoluteIncrement_);
124    if (relativeIncrement_ != other.relativeIncrement_)
125        fprintf(fp, "3  heuristicFPump.setRelativeIncrement(%g);\n", relativeIncrement_);
126    else
127        fprintf(fp, "4  heuristicFPump.setRelativeIncrement(%g);\n", relativeIncrement_);
128    if (defaultRounding_ != other.defaultRounding_)
129        fprintf(fp, "3  heuristicFPump.setDefaultRounding(%g);\n", defaultRounding_);
130    else
131        fprintf(fp, "4  heuristicFPump.setDefaultRounding(%g);\n", defaultRounding_);
132    if (initialWeight_ != other.initialWeight_)
133        fprintf(fp, "3  heuristicFPump.setInitialWeight(%g);\n", initialWeight_);
134    else
135        fprintf(fp, "4  heuristicFPump.setInitialWeight(%g);\n", initialWeight_);
136    if (weightFactor_ != other.weightFactor_)
137        fprintf(fp, "3  heuristicFPump.setWeightFactor(%g);\n", weightFactor_);
138    else
139        fprintf(fp, "4  heuristicFPump.setWeightFactor(%g);\n", weightFactor_);
140    if (artificialCost_ != other.artificialCost_)
141        fprintf(fp, "3  heuristicFPump.setArtificialCost(%g);\n", artificialCost_);
142    else
143        fprintf(fp, "4  heuristicFPump.setArtificialCost(%g);\n", artificialCost_);
144    if (iterationRatio_ != other.iterationRatio_)
145        fprintf(fp, "3  heuristicFPump.setIterationRatio(%g);\n", iterationRatio_);
146    else
147        fprintf(fp, "4  heuristicFPump.setIterationRatio(%g);\n", iterationRatio_);
148    if (reducedCostMultiplier_ != other.reducedCostMultiplier_)
149        fprintf(fp, "3  heuristicFPump.setReducedCostMultiplier(%g);\n", reducedCostMultiplier_);
150    else
151        fprintf(fp, "4  heuristicFPump.setReducedCostMultiplier(%g);\n", reducedCostMultiplier_);
152    fprintf(fp, "3  cbcModel->addHeuristic(&heuristicFPump);\n");
153}
154
155// Copy constructor
156CbcHeuristicFPump::CbcHeuristicFPump(const CbcHeuristicFPump & rhs)
157        :
158        CbcHeuristic(rhs),
159        startTime_(rhs.startTime_),
160        maximumTime_(rhs.maximumTime_),
161        fakeCutoff_(rhs.fakeCutoff_),
162        absoluteIncrement_(rhs.absoluteIncrement_),
163        relativeIncrement_(rhs.relativeIncrement_),
164        defaultRounding_(rhs.defaultRounding_),
165        initialWeight_(rhs.initialWeight_),
166        weightFactor_(rhs.weightFactor_),
167        artificialCost_(rhs.artificialCost_),
168        iterationRatio_(rhs.iterationRatio_),
169        reducedCostMultiplier_(rhs.reducedCostMultiplier_),
170        maximumPasses_(rhs.maximumPasses_),
171        maximumRetries_(rhs.maximumRetries_),
172        accumulate_(rhs.accumulate_),
173        fixOnReducedCosts_(rhs.fixOnReducedCosts_),
174        roundExpensive_(rhs.roundExpensive_)
175{
176}
177
178// Assignment operator
179CbcHeuristicFPump &
180CbcHeuristicFPump::operator=( const CbcHeuristicFPump & rhs)
181{
182    if (this != &rhs) {
183        CbcHeuristic::operator=(rhs);
184        startTime_ = rhs.startTime_;
185        maximumTime_ = rhs.maximumTime_;
186        fakeCutoff_ = rhs.fakeCutoff_;
187        absoluteIncrement_ = rhs.absoluteIncrement_;
188        relativeIncrement_ = rhs.relativeIncrement_;
189        defaultRounding_ = rhs.defaultRounding_;
190        initialWeight_ = rhs.initialWeight_;
191        weightFactor_ = rhs.weightFactor_;
192        artificialCost_ = rhs.artificialCost_;
193        iterationRatio_ = rhs.iterationRatio_;
194        reducedCostMultiplier_ = rhs.reducedCostMultiplier_;
195        maximumPasses_ = rhs.maximumPasses_;
196        maximumRetries_ = rhs.maximumRetries_;
197        accumulate_ = rhs.accumulate_;
198        fixOnReducedCosts_ = rhs.fixOnReducedCosts_;
199        roundExpensive_ = rhs.roundExpensive_;
200    }
201    return *this;
202}
203
204// Resets stuff if model changes
205void
206CbcHeuristicFPump::resetModel(CbcModel * )
207{
208}
209
210/**************************BEGIN MAIN PROCEDURE ***********************************/
211
212// See if feasibility pump will give better solution
213// Sets value of solution
214// Returns 1 if solution, 0 if not
215int
216CbcHeuristicFPump::solution(double & solutionValue,
217                            double * betterSolution)
218{
219    startTime_ = CoinCpuTime();
220    numCouldRun_++;
221    double incomingObjective = solutionValue;
222#define LEN_PRINT 200
223    char pumpPrint[LEN_PRINT];
224    pumpPrint[0] = '\0';
225/*
226  Decide if we want to run. Standard values for when are described in
227  CbcHeuristic.hpp. If we're off, or running only at root and this isn't the
228  root, bail out.
229
230  The double test (against phase, then atRoot and passNumber) has a fair bit
231  of redundancy, but the results will differ depending on whether we're
232  actually at the root of the main search tree or at the root of a small tree
233  (recursive call to branchAndBound).
234
235  FPump also supports some exotic values (11 -- 15) for when, described in
236  CbcHeuristicFPump.hpp.
237*/
238    if (!when() || (when() == 1 && model_->phase() != 1))
239        return 0; // switched off
240    // See if at root node
241    bool atRoot = model_->getNodeCount() == 0;
242    int passNumber = model_->getCurrentPassNumber();
243    // just do once
244    if (!atRoot)
245        return 0;
246    int options = feasibilityPumpOptions_;
247    if ((options % 1000000) > 0) {
248        int kOption = options / 1000000;
249        options = options % 1000000;
250        /*
251          Add 10 to do even if solution
252          1 - do after cuts
253          2 - do after cuts (not before)
254          3 - not used do after every cut round (and after cuts)
255          k not used do after every (k-2)th round
256        */
257        if (kOption < 10 && model_->getSolutionCount())
258            return 0;
259        if (model_->getSolutionCount())
260            kOption = kOption % 10;
261        bool good;
262        if (kOption == 1) {
263            good = (passNumber == 999999);
264        } else if (kOption == 2) {
265            good = (passNumber == 999999);
266            passNumber = 2; // so won't run before
267            //} else if (kOption==3) {
268            //good = true;
269        } else {
270            //good = (((passNumber-1)%(kOption-2))==0);
271            good = false;
272        }
273        if (passNumber != 1 && !good)
274            return 0;
275    } else {
276        if (passNumber != 1)
277            return 0;
278    }
279    // loop round doing repeated pumps
280    double cutoff;
281    model_->solver()->getDblParam(OsiDualObjectiveLimit, cutoff);
282    double direction = model_->solver()->getObjSense();
283    cutoff *= direction;
284    int numberBandBsolutions = 0;
285    double firstCutoff = fabs(cutoff);
286    cutoff = CoinMin(cutoff, solutionValue);
287    // check plausible and space for rounded solution
288    int numberColumns = model_->getNumCols();
289    int numberIntegers = model_->numberIntegers();
290    const int * integerVariableOrig = model_->integerVariable();
291    double iterationLimit = -1.0;
292    //iterationRatio_=1.0;
293    if (iterationRatio_ > 0.0)
294        iterationLimit = (2 * model_->solver()->getNumRows() + 2 * numberColumns) *
295                         iterationRatio_;
296    int totalNumberIterations = 0;
297    int averageIterationsPerTry = -1;
298    int numberIterationsLastPass = 0;
299    // 1. initially check 0-1
300/*
301  I'm skeptical of the above comment, but it's likely accurate as the default.
302  Bit 4 or bit 8 needs to be set in order to consider working with general
303  integers.
304*/
305    int i, j;
306    int general = 0;
307    int * integerVariable = new int[numberIntegers];
308    const double * lower = model_->solver()->getColLower();
309    const double * upper = model_->solver()->getColUpper();
310    bool doGeneral = (accumulate_ & 4) != 0;
311    j = 0;
312/*
313  Scan the objects, recording the columns and counting general integers.
314
315  Seems like the NDEBUG tests could be made into an applicability test. If
316  a scan of the objects reveals complex objects, just clean up and return
317  failure.
318*/
319    for (i = 0; i < numberIntegers; i++) {
320        int iColumn = integerVariableOrig[i];
321#ifndef NDEBUG
322        const OsiObject * object = model_->object(i);
323        const CbcSimpleInteger * integerObject =
324            dynamic_cast<const  CbcSimpleInteger *> (object);
325        const OsiSimpleInteger * integerObject2 =
326            dynamic_cast<const  OsiSimpleInteger *> (object);
327        assert(integerObject || integerObject2);
328#endif
329        if (upper[iColumn] - lower[iColumn] > 1.000001) {
330            general++;
331            if (doGeneral)
332                integerVariable[j++] = iColumn;
333        } else {
334            integerVariable[j++] = iColumn;
335        }
336    }
337/*
338  If 2/3 of integers are general integers, and we're not going to work with
339  them, might as well go home.
340
341  The else case is unclear to me. We reach it if general integers are less than
342  2/3 of the total, or if either of bit 4 or 8 is set. But only bit 8 is used
343  in the decision. (Let manyGen = 1 if more than 2/3 of integers are general
344  integers. Then a k-map on manyGen, bit4, and bit8 shows it clearly.)
345
346  So there's something odd here. In the case where bit4 = 1 and bit8 = 0,
347  we've included general integers in integerVariable, but we're not going to
348  process them.
349*/
350    if (general*3 > 2*numberIntegers && !doGeneral) {
351        delete [] integerVariable;
352        return 0;
353    } else if ((accumulate_&4) == 0) {
354        doGeneral = false;
355        j = 0;
356        for (i = 0; i < numberIntegers; i++) {
357            int iColumn = integerVariableOrig[i];
358            if (upper[iColumn] - lower[iColumn] < 1.000001)
359                integerVariable[j++] = iColumn;
360        }
361    }
362    if (!general)
363        doGeneral = false;
364#ifdef CLP_INVESTIGATE
365    if (doGeneral)
366        printf("DOing general with %d out of %d\n", general, numberIntegers);
367#endif
368/*
369  This `closest solution' will satisfy integrality, but violate some other
370  constraints?
371*/
372    // For solution closest to feasible if none found
373    int * closestSolution = general ? NULL : new int[numberIntegers];
374    double closestObjectiveValue = COIN_DBL_MAX;
375
376    int numberIntegersOrig = numberIntegers;
377    numberIntegers = j;
378    double * newSolution = new double [numberColumns];
379    double newSolutionValue = COIN_DBL_MAX;
380    int maxSolutions = model_->getMaximumSolutions();
381    int numberSolutions=0;
382    bool solutionFound = false;
383    int * usedColumn = NULL;
384    double * lastSolution = NULL;
385    int fixContinuous = 0;
386    bool fixInternal = false;
387    bool allSlack = false;
388    if (when_ >= 21 && when_ <= 25) {
389        when_ -= 10;
390        allSlack = true;
391    }
392    double time1 = CoinCpuTime();
393/*
394  Obtain a relaxed lp solution.
395*/
396    model_->solver()->resolve();
397    if (!model_->solver()->isProvenOptimal()) {
398        // presumably max time or some such
399        return 0;
400    }
401    numRuns_++;
402    if (cutoff < 1.0e50 && false) {
403        // Fix on djs
404        double direction = model_->solver()->getObjSense() ;
405        double gap = cutoff - model_->solver()->getObjValue() * direction ;
406        double tolerance;
407        model_->solver()->getDblParam(OsiDualTolerance, tolerance) ;
408        if (gap > 0.0) {
409            gap += 100.0 * tolerance;
410            int nFix = model_->solver()->reducedCostFix(gap);
411            printf("dj fixing fixed %d variables\n", nFix);
412        }
413    }
414/*
415  I have no idea why we're doing this, except perhaps that saveBasis will be
416  automagically deleted on exit from the routine.
417*/
418    CoinWarmStartBasis saveBasis;
419    CoinWarmStartBasis * basis =
420        dynamic_cast<CoinWarmStartBasis *>(model_->solver()->getWarmStart()) ;
421    if (basis) {
422        saveBasis = * basis;
423        delete basis;
424    }
425    double continuousObjectiveValue = model_->solver()->getObjValue() * model_->solver()->getObjSense();
426    double * firstPerturbedObjective = NULL;
427    double * firstPerturbedSolution = NULL;
428    double firstPerturbedValue = COIN_DBL_MAX;
429    if (when_ >= 11 && when_ <= 15) {
430        fixInternal = when_ > 11 && when_ < 15;
431        if (when_ < 13)
432            fixContinuous = 0;
433        else if (when_ != 14)
434            fixContinuous = 1;
435        else
436            fixContinuous = 2;
437        when_ = 1;
438        if ((accumulate_&1) != 0) {
439            usedColumn = new int [numberColumns];
440            for (int i = 0; i < numberColumns; i++)
441                usedColumn[i] = -1;
442        }
443        lastSolution = CoinCopyOfArray(model_->solver()->getColSolution(), numberColumns);
444    }
445    int finalReturnCode = 0;
446    int totalNumberPasses = 0;
447    int numberTries = 0;
448    CoinWarmStartBasis bestBasis;
449    bool exitAll = false;
450    //double saveBestObjective = model_->getMinimizationObjValue();
451    OsiSolverInterface * solver = NULL;
452    double artificialFactor = 0.00001;
453    // also try rounding!
454    double * roundingSolution = new double[numberColumns];
455    double roundingObjective = solutionValue;
456    CbcRounding roundingHeuristic(*model_);
457    int dualPass = 0;
458    int secondPassOpt = 0;
459#define RAND_RAND
460#ifdef RAND_RAND
461    int offRandom = 0;
462#endif
463    int maximumAllowed = -1;
464    bool moreIterations = false;
465    if (options > 0) {
466        if (options >= 1000)
467            maximumAllowed = options / 1000;
468        int options2 = (options % 1000) / 100;
469#ifdef RAND_RAND
470        offRandom = options2 & 1;
471#endif
472        moreIterations = (options2 & 2) != 0;
473        secondPassOpt = (options / 10) % 10;
474        /* 1 to 7 - re-use solution
475           8 use dual and current solution(ish)
476           9 use dual and allslack
477           1 - primal and mod obj
478           2 - dual and mod obj
479           3 - primal and no mod obj
480           add 3 to redo current solution
481        */
482        if (secondPassOpt >= 8) {
483            dualPass = secondPassOpt - 7;
484            secondPassOpt = 0;
485        }
486    }
487    // Number of passes to do
488    int maximumPasses = maximumPasses_;
489#ifdef COIN_HAS_CLP
490    {
491      OsiClpSolverInterface * clpSolver
492        = dynamic_cast<OsiClpSolverInterface *> (model_->solver());
493      if (clpSolver ) {
494        if (maximumPasses == 30) {
495          if (clpSolver->fakeObjective())
496            maximumPasses = 100; // feasibility problem?
497        }
498        if (model_->getRandomSeed()!=-1)
499          clpSolver->getModelPtr()->setRandomSeed(randomNumberGenerator_.getSeed());
500      }
501    }
502#endif
503#ifdef RAND_RAND
504    double * randomFactor = new double [numberColumns];
505    for (int i = 0; i < numberColumns; i++) {
506        double value = floor(1.0e3 * randomNumberGenerator_.randomDouble());
507        randomFactor[i] = 1.0 + value * 1.0e-4;
508    }
509#endif
510    // guess exact multiple of objective
511    double exactMultiple = model_->getCutoffIncrement();
512    exactMultiple *= 2520;
513    if (fabs(exactMultiple / 0.999 - floor(exactMultiple / 0.999 + 0.5)) < 1.0e-9)
514        exactMultiple /= 2520.0 * 0.999;
515    else if (fabs(exactMultiple - floor(exactMultiple + 0.5)) < 1.0e-9)
516        exactMultiple /= 2520.0;
517    else
518        exactMultiple = 0.0;
519    // check for rounding errors (only for integral case)
520    if (fabs(exactMultiple - floor(exactMultiple + 0.5)) < 1.0e-8)
521        exactMultiple = floor(exactMultiple + 0.5);
522    //printf("exact multiple %g\n",exactMultiple);
523    // Clone solver for rounding
524    OsiSolverInterface * clonedSolver = cloneBut(2); // wasmodel_->solver()->clone();
525    while (!exitAll) {
526        // Cutoff rhs
527        double useRhs = COIN_DBL_MAX;
528        double useOffset = 0.0;
529        int numberPasses = 0;
530        artificialFactor *= 10.0;
531        int lastMove = (!numberTries) ? -10 : 1000000;
532        double lastSumInfeas = COIN_DBL_MAX;
533        numberTries++;
534        // Clone solver - otherwise annoys root node computations
535        solver = cloneBut(2); // was model_->solver()->clone();
536#ifdef COIN_HAS_CLP
537        {
538            OsiClpSolverInterface * clpSolver
539            = dynamic_cast<OsiClpSolverInterface *> (solver);
540            if (clpSolver) {
541                // better to clean up using primal?
542                ClpSimplex * lp = clpSolver->getModelPtr();
543                int options = lp->specialOptions();
544                lp->setSpecialOptions(options | 8192);
545                //lp->setSpecialOptions(options|0x01000000);
546#ifdef CLP_INVESTIGATE
547                clpSolver->setHintParam(OsiDoReducePrint, false, OsiHintTry);
548                lp->setLogLevel(CoinMax(1, lp->logLevel()));
549#endif
550            }
551        }
552#endif
553        if (CoinMin(fakeCutoff_, cutoff) < 1.0e50) {
554            // Fix on djs
555            double direction = solver->getObjSense() ;
556            double gap = CoinMin(fakeCutoff_, cutoff) - solver->getObjValue() * direction ;
557            double tolerance;
558            solver->getDblParam(OsiDualTolerance, tolerance) ;
559            if (gap > 0.0 && (fixOnReducedCosts_ == 1 || (numberTries == 1 && fixOnReducedCosts_ == 2))) {
560                gap += 100.0 * tolerance;
561                gap *= reducedCostMultiplier_;
562                int nFix = solver->reducedCostFix(gap);
563                if (nFix) {
564                    sprintf(pumpPrint, "Reduced cost fixing fixed %d variables on major pass %d", nFix, numberTries);
565                    model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
566                    << pumpPrint
567                    << CoinMessageEol;
568                    //pumpPrint[0]='\0';
569                }
570            }
571        }
572        // if cutoff exists then add constraint
573        bool useCutoff = (fabs(cutoff) < 1.0e20 && (fakeCutoff_ != COIN_DBL_MAX || numberTries > 1));
574        // but there may be a close one
575        if (firstCutoff < 2.0*solutionValue && numberTries == 1 && CoinMin(cutoff, fakeCutoff_) < 1.0e20)
576            useCutoff = true;
577        if (useCutoff) {
578            double rhs = CoinMin(cutoff, fakeCutoff_);
579            const double * objective = solver->getObjCoefficients();
580            int numberColumns = solver->getNumCols();
581            int * which = new int[numberColumns];
582            double * els = new double[numberColumns];
583            int nel = 0;
584            for (int i = 0; i < numberColumns; i++) {
585                double value = objective[i];
586                if (value) {
587                    which[nel] = i;
588                    els[nel++] = direction * value;
589                }
590            }
591            solver->getDblParam(OsiObjOffset, useOffset);
592#ifdef COIN_DEVELOP
593            if (useOffset)
594                printf("CbcHeuristicFPump obj offset %g\n", useOffset);
595#endif
596            useOffset *= direction;
597            // Tweak rhs and save
598            useRhs = rhs;
599#ifdef JJF_ZERO
600            double tempValue = 60.0 * useRhs;
601            if (fabs(tempValue - floor(tempValue + 0.5)) < 1.0e-7 && rhs != fakeCutoff_) {
602                // add a little
603                useRhs += 1.0e-5;
604            }
605#endif
606            solver->addRow(nel, which, els, -COIN_DBL_MAX, useRhs + useOffset);
607            delete [] which;
608            delete [] els;
609            bool takeHint;
610            OsiHintStrength strength;
611            solver->getHintParam(OsiDoDualInResolve, takeHint, strength);
612            solver->setHintParam(OsiDoDualInResolve, true, OsiHintDo);
613            solver->resolve();
614            solver->setHintParam(OsiDoDualInResolve, takeHint, strength);
615            if (!solver->isProvenOptimal()) {
616                // presumably max time or some such
617                exitAll = true;
618                break;
619            }
620        }
621        solver->setDblParam(OsiDualObjectiveLimit, 1.0e50);
622        solver->resolve();
623        // Solver may not be feasible
624        if (!solver->isProvenOptimal()) {
625            exitAll = true;
626            break;
627        }
628        const double * lower = solver->getColLower();
629        const double * upper = solver->getColUpper();
630        const double * solution = solver->getColSolution();
631        if (lastSolution)
632            memcpy(lastSolution, solution, numberColumns*sizeof(double));
633        double primalTolerance;
634        solver->getDblParam(OsiPrimalTolerance, primalTolerance);
635
636        // 2 space for last rounded solutions
637#define NUMBER_OLD 4
638        double ** oldSolution = new double * [NUMBER_OLD];
639        for (j = 0; j < NUMBER_OLD; j++) {
640            oldSolution[j] = new double[numberColumns];
641            for (i = 0; i < numberColumns; i++) oldSolution[j][i] = -COIN_DBL_MAX;
642        }
643
644        // 3. Replace objective with an initial 0-valued objective
645        double * saveObjective = new double [numberColumns];
646        memcpy(saveObjective, solver->getObjCoefficients(), numberColumns*sizeof(double));
647        for (i = 0; i < numberColumns; i++) {
648            solver->setObjCoeff(i, 0.0);
649        }
650        bool finished = false;
651        double direction = solver->getObjSense();
652        int returnCode = 0;
653        bool takeHint;
654        OsiHintStrength strength;
655        solver->getHintParam(OsiDoDualInResolve, takeHint, strength);
656        solver->setHintParam(OsiDoDualInResolve, false);
657        //solver->messageHandler()->setLogLevel(0);
658
659        // 4. Save objective offset so we can see progress
660        double saveOffset;
661        solver->getDblParam(OsiObjOffset, saveOffset);
662        // Get amount for original objective
663        double scaleFactor = 0.0;
664#ifdef COIN_DEVELOP
665        double largestCost = 0.0;
666        int nArtificial = 0;
667#endif
668        for (i = 0; i < numberColumns; i++) {
669            double value = saveObjective[i];
670            scaleFactor += value * value;
671#ifdef COIN_DEVELOP
672            largestCost = CoinMax(largestCost, fabs(value));
673            if (value*direction >= artificialCost_)
674                nArtificial++;
675#endif
676        }
677        if (scaleFactor)
678            scaleFactor = (initialWeight_ * sqrt(static_cast<double> (numberIntegers))) / sqrt(scaleFactor);
679#ifdef CLP_INVESTIGATE
680#ifdef COIN_DEVELOP
681        if (scaleFactor || nArtificial)
682            printf("Using %g fraction of original objective (decay %g) - largest %g - %d artificials\n", scaleFactor, weightFactor_,
683                   largestCost, nArtificial);
684#else
685        if (scaleFactor)
686            printf("Using %g fraction of original objective (decay %g)\n",
687                   scaleFactor, weightFactor_);
688#endif
689#endif
690        // This is an array of sums of infeasibilities so can see if "bobbling"
691#define SIZE_BOBBLE 20
692        double saveSumInf[SIZE_BOBBLE];
693        CoinFillN(saveSumInf, SIZE_BOBBLE, COIN_DBL_MAX);
694        // 0 before doing anything
695        int bobbleMode = 0;
696        // 5. MAIN WHILE LOOP
697        //bool newLineNeeded=false;
698/*
699  finished occurs exactly twice in this routine: immediately above, where it's
700  set to false, and here in the loop condition.
701*/
702        while (!finished) {
703            double newTrueSolutionValue = 0.0;
704            double newSumInfeas = 0.0;
705            int newNumberInfeas = 0;
706            returnCode = 0;
707            if (model_->maximumSecondsReached()) {
708                exitAll = true;
709                break;
710            }
711            // see what changed
712            if (usedColumn) {
713                for (i = 0; i < numberColumns; i++) {
714                    if (fabs(solution[i] - lastSolution[i]) > 1.0e-8)
715                        usedColumn[i] = numberPasses;
716                    lastSolution[i] = solution[i];
717                }
718            }
719            if (averageIterationsPerTry >= 0) {
720                int n = totalNumberIterations - numberIterationsLastPass;
721                double perPass = totalNumberIterations /
722                                 (totalNumberPasses + numberPasses + 1.0e-5);
723                perPass /= (solver->getNumRows() + numberColumns);
724                double test = moreIterations ? 0.3 : 0.05;
725                if (n > CoinMax(20000, 3*averageIterationsPerTry)
726                        && (switches_&2) == 0 && maximumPasses<200 && perPass>test) {
727                    exitAll = true;
728                }
729            }
730            // Exit on exact total number if maximumPasses large
731            if ((maximumPasses >= 200 || (switches_&2) != 0)
732                    && numberPasses + totalNumberPasses >=
733                    maximumPasses)
734                exitAll = true;
735            bool exitThis = false;
736            if (iterationLimit < 0.0) {
737                if (numberPasses >= maximumPasses) {
738                    // If going well then keep going if maximumPasses small
739                    if (lastMove < numberPasses - 4 || lastMove == 1000000)
740                        exitThis = true;
741                    if (maximumPasses > 20 || numberPasses >= 40)
742                        exitThis = true;
743                }
744            }
745            if (iterationLimit > 0.0 && totalNumberIterations > iterationLimit
746                    && numberPasses > 15) {
747                // exiting on iteration count
748                exitAll = true;
749            } else if (maximumPasses<30 && numberPasses>100) {
750                // too many passes anyway
751                exitAll = true;
752            }
753            if (maximumTime_ > 0.0 && CoinCpuTime() >= startTime_ + maximumTime_) {
754                exitAll = true;
755                // force exit
756                switches_ |= 2048;
757            }
758            if (exitAll || exitThis)
759                break;
760            memcpy(newSolution, solution, numberColumns*sizeof(double));
761            int flip;
762            if (numberPasses == 0 && false) {
763                // always use same seed
764                randomNumberGenerator_.setSeed(987654321);
765            }
766            returnCode = rounds(solver, newSolution,/*saveObjective,*/
767                                numberIntegers, integerVariable,
768                                /*pumpPrint,*/numberPasses,
769                                /*roundExpensive_,*/defaultRounding_, &flip);
770            if (numberPasses == 0 && false) {
771                // Make sure random will be different
772                for (i = 1; i < numberTries; i++)
773                    randomNumberGenerator_.randomDouble();
774            }
775            numberPasses++;
776            if (returnCode) {
777                // SOLUTION IS INTEGER
778                // Put back correct objective
779                for (i = 0; i < numberColumns; i++)
780                    solver->setObjCoeff(i, saveObjective[i]);
781
782                // solution - but may not be better
783                // Compute using dot product
784                solver->setDblParam(OsiObjOffset, saveOffset);
785                newSolutionValue = -saveOffset;
786                for (  i = 0 ; i < numberColumns ; i++ )
787                    newSolutionValue += saveObjective[i] * newSolution[i];
788                newSolutionValue *= direction;
789                sprintf(pumpPrint, "Solution found of %g", newSolutionValue);
790                model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
791                << pumpPrint
792                << CoinMessageEol;
793                //newLineNeeded=false;
794                if (newSolutionValue < solutionValue) {
795                    double saveValue = solutionValue;
796                    if (!doGeneral) {
797                        int numberLeft = 0;
798                        for (i = 0; i < numberIntegersOrig; i++) {
799                            int iColumn = integerVariableOrig[i];
800                            double value = floor(newSolution[iColumn] + 0.5);
801                            if (solver->isBinary(iColumn)) {
802                                solver->setColLower(iColumn, value);
803                                solver->setColUpper(iColumn, value);
804                            } else {
805                                if (fabs(value - newSolution[iColumn]) > 1.0e-7)
806                                    numberLeft++;
807                            }
808                        }
809                        if (numberLeft) {
810                            sprintf(pumpPrint, "Branch and bound needed to clear up %d general integers", numberLeft);
811                            model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
812                            << pumpPrint
813                            << CoinMessageEol;
814                            returnCode = smallBranchAndBound(solver, numberNodes_, newSolution, newSolutionValue,
815                                                             solutionValue, "CbcHeuristicFpump");
816                            if (returnCode < 0) {
817                                if (returnCode == -2)
818                                    exitAll = true;
819                                returnCode = 0; // returned on size or event
820                            }
821                            if ((returnCode&2) != 0) {
822                                // could add cut
823                                returnCode &= ~2;
824                            }
825                            if (returnCode != 1)
826                                newSolutionValue = saveValue;
827                            if (returnCode && newSolutionValue < saveValue)
828                                numberBandBsolutions++;
829                        }
830                    }
831                    if (returnCode && newSolutionValue < saveValue) {
832                        memcpy(betterSolution, newSolution, numberColumns*sizeof(double));
833                        solutionFound = true;
834                        if (exitNow(newSolutionValue))
835                            exitAll = true;
836                        CoinWarmStartBasis * basis =
837                            dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ;
838                        if (basis) {
839                            bestBasis = * basis;
840                            delete basis;
841                            int action = model_->dealWithEventHandler(CbcEventHandler::heuristicSolution, newSolutionValue, betterSolution);
842                            if (action == 0) {
843                                double * saveOldSolution = CoinCopyOfArray(model_->bestSolution(), numberColumns);
844                                double saveObjectiveValue = model_->getMinimizationObjValue();
845                                model_->setBestSolution(betterSolution, numberColumns, newSolutionValue);
846                                if (saveOldSolution && saveObjectiveValue < model_->getMinimizationObjValue())
847                                    model_->setBestSolution(saveOldSolution, numberColumns, saveObjectiveValue);
848                                delete [] saveOldSolution;
849                            }
850                            if (action == 0 || model_->maximumSecondsReached()) {
851                                exitAll = true; // exit
852                                break;
853                            }
854                        }
855                        if ((accumulate_&1) != 0) {
856                            model_->incrementUsed(betterSolution); // for local search
857                        }
858                        solutionValue = newSolutionValue;
859                        solutionFound = true;
860                        numberSolutions++;
861                        if (numberSolutions>=maxSolutions)
862                          exitAll = true;
863                        if (general && saveValue != newSolutionValue) {
864                            sprintf(pumpPrint, "Cleaned solution of %g", solutionValue);
865                            model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
866                            << pumpPrint
867                            << CoinMessageEol;
868                        }
869                        if (exitNow(newSolutionValue))
870                            exitAll = true;
871                    } else {
872                        sprintf(pumpPrint, "Mini branch and bound could not fix general integers");
873                        model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
874                        << pumpPrint
875                        << CoinMessageEol;
876                    }
877                } else {
878                    sprintf(pumpPrint, "After further testing solution no better than previous of %g", solutionValue);
879                    model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
880                    << pumpPrint
881                    << CoinMessageEol;
882                    //newLineNeeded=false;
883                    returnCode = 0;
884                }
885                break;
886            } else {
887                // SOLUTION IS not INTEGER
888                // 1. check for loop
889                bool matched;
890                for (int k = NUMBER_OLD - 1; k > 0; k--) {
891                    double * b = oldSolution[k];
892                    matched = true;
893                    for (i = 0; i < numberIntegers; i++) {
894                        int iColumn = integerVariable[i];
895                        if (newSolution[iColumn] != b[iColumn]) {
896                            matched = false;
897                            break;
898                        }
899                    }
900                    if (matched) break;
901                }
902                int numberPerturbed = 0;
903                if (matched || numberPasses % 100 == 0) {
904                    // perturbation
905                    //sprintf(pumpPrint+strlen(pumpPrint)," perturbation applied");
906                    //newLineNeeded=true;
907                    double factorX[10] = {0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0};
908                    double factor = 1.0;
909                    double target = -1.0;
910                    double * randomX = new double [numberIntegers];
911                    for (i = 0; i < numberIntegers; i++)
912                        randomX[i] = CoinMax(0.0, randomNumberGenerator_.randomDouble() - 0.3);
913                    for (int k = 0; k < 10; k++) {
914#ifdef COIN_DEVELOP_x
915                        printf("kpass %d\n", k);
916#endif
917                        int numberX[10] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
918                        for (i = 0; i < numberIntegers; i++) {
919                            int iColumn = integerVariable[i];
920                            double value = randomX[i];
921                            double difference = fabs(solution[iColumn] - newSolution[iColumn]);
922                            for (int j = 0; j < 10; j++) {
923                                if (difference + value*factorX[j] > 0.5)
924                                    numberX[j]++;
925                            }
926                        }
927                        if (target < 0.0) {
928                            if (numberX[9] <= 200)
929                                break; // not very many changes
930                            target = CoinMax(200.0, CoinMin(0.05 * numberX[9], 1000.0));
931                        }
932                        int iX = -1;
933                        int iBand = -1;
934                        for (i = 0; i < 10; i++) {
935#ifdef COIN_DEVELOP_x
936                            printf("** %d changed at %g\n", numberX[i], factorX[i]);
937#endif
938                            if (numberX[i] >= target && numberX[i] < 2.0*target && iX < 0)
939                                iX = i;
940                            if (iBand<0 && numberX[i]>target) {
941                                iBand = i;
942                                factor = factorX[i];
943                            }
944                        }
945                        if (iX >= 0) {
946                            factor = factorX[iX];
947                            break;
948                        } else {
949                            assert (iBand >= 0);
950                            double hi = factor;
951                            double lo = (iBand > 0) ? factorX[iBand-1] : 0.0;
952                            double diff = (hi - lo) / 9.0;
953                            for (i = 0; i < 10; i++) {
954                                factorX[i] = lo;
955                                lo += diff;
956                            }
957                        }
958                    }
959                    for (i = 0; i < numberIntegers; i++) {
960                        int iColumn = integerVariable[i];
961                        double value = randomX[i];
962                        double difference = fabs(solution[iColumn] - newSolution[iColumn]);
963                        if (difference + value*factor > 0.5) {
964                            numberPerturbed++;
965                            if (newSolution[iColumn] < lower[iColumn] + primalTolerance) {
966                                newSolution[iColumn] += 1.0;
967                            } else if (newSolution[iColumn] > upper[iColumn] - primalTolerance) {
968                                newSolution[iColumn] -= 1.0;
969                            } else {
970                                // general integer
971                                if (difference + value > 0.75)
972                                    newSolution[iColumn] += 1.0;
973                                else
974                                    newSolution[iColumn] -= 1.0;
975                            }
976                        }
977                    }
978                    delete [] randomX;
979                } else {
980                    for (j = NUMBER_OLD - 1; j > 0; j--) {
981                        for (i = 0; i < numberColumns; i++) oldSolution[j][i] = oldSolution[j-1][i];
982                    }
983                    for (j = 0; j < numberColumns; j++) oldSolution[0][j] = newSolution[j];
984                }
985
986                // 2. update the objective function based on the new rounded solution
987                double offset = 0.0;
988                double costValue = (1.0 - scaleFactor) * solver->getObjSense();
989                int numberChanged = 0;
990                const double * oldObjective = solver->getObjCoefficients();
991                for (i = 0; i < numberColumns; i++) {
992                    // below so we can keep original code and allow for objective
993                    int iColumn = i;
994                    // Special code for "artificials"
995                    if (direction*saveObjective[iColumn] >= artificialCost_) {
996                        //solver->setObjCoeff(iColumn,scaleFactor*saveObjective[iColumn]);
997                        solver->setObjCoeff(iColumn, (artificialFactor*saveObjective[iColumn]) / artificialCost_);
998                    }
999                    if (!solver->isBinary(iColumn) && !doGeneral)
1000                        continue;
1001                    // deal with fixed variables (i.e., upper=lower)
1002                    if (fabs(lower[iColumn] - upper[iColumn]) < primalTolerance || !solver->isInteger(iColumn)) {
1003                        //if (lower[iColumn] > 1. - primalTolerance) solver->setObjCoeff(iColumn,-costValue);
1004                        //else                                       solver->setObjCoeff(iColumn,costValue);
1005                        continue;
1006                    }
1007                    double newValue = 0.0;
1008                    if (newSolution[iColumn] < lower[iColumn] + primalTolerance) {
1009                        newValue = costValue + scaleFactor * saveObjective[iColumn];
1010                    } else {
1011                        if (newSolution[iColumn] > upper[iColumn] - primalTolerance) {
1012                            newValue = -costValue + scaleFactor * saveObjective[iColumn];
1013                        }
1014                    }
1015#ifdef RAND_RAND
1016                    if (!offRandom)
1017                        newValue *= randomFactor[iColumn];
1018#endif
1019                    if (newValue != oldObjective[iColumn]) {
1020                        numberChanged++;
1021                    }
1022                    solver->setObjCoeff(iColumn, newValue);
1023                    offset += costValue * newSolution[iColumn];
1024                }
1025                solver->setDblParam(OsiObjOffset, -offset);
1026                if (!general && false) {
1027                    // Solve in two goes - first keep satisfied ones fixed
1028                    double * saveLower = new double [numberIntegers];
1029                    double * saveUpper = new double [numberIntegers];
1030                    for (i = 0; i < numberIntegers; i++) {
1031                        int iColumn = integerVariable[i];
1032                        saveLower[i] = COIN_DBL_MAX;
1033                        saveUpper[i] = -COIN_DBL_MAX;
1034                        if (solution[iColumn] < lower[iColumn] + primalTolerance) {
1035                            saveUpper[i] = upper[iColumn];
1036                            solver->setColUpper(iColumn, lower[iColumn]);
1037                        } else if (solution[iColumn] > upper[iColumn] - primalTolerance) {
1038                            saveLower[i] = lower[iColumn];
1039                            solver->setColLower(iColumn, upper[iColumn]);
1040                        }
1041                    }
1042                    solver->resolve();
1043                    if (!solver->isProvenOptimal()) {
1044                        // presumably max time or some such
1045                        exitAll = true;
1046                        break;
1047                    }
1048                    for (i = 0; i < numberIntegers; i++) {
1049                        int iColumn = integerVariable[i];
1050                        if (saveLower[i] != COIN_DBL_MAX)
1051                            solver->setColLower(iColumn, saveLower[i]);
1052                        if (saveUpper[i] != -COIN_DBL_MAX)
1053                            solver->setColUpper(iColumn, saveUpper[i]);
1054                        saveUpper[i] = -COIN_DBL_MAX;
1055                    }
1056                    memcpy(newSolution, solution, numberColumns*sizeof(double));
1057                    int flip;
1058                    returnCode = rounds(solver, newSolution,/*saveObjective,*/
1059                                        numberIntegers, integerVariable,
1060                                        /*pumpPrint,*/numberPasses,
1061                                        /*roundExpensive_,*/defaultRounding_, &flip);
1062                    numberPasses++;
1063                    if (returnCode) {
1064                        // solution - but may not be better
1065                        // Compute using dot product
1066                        double newSolutionValue = -saveOffset;
1067                        for (  i = 0 ; i < numberColumns ; i++ )
1068                            newSolutionValue += saveObjective[i] * newSolution[i];
1069                        newSolutionValue *= direction;
1070                        sprintf(pumpPrint, "Intermediate solution found of %g", newSolutionValue);
1071                        model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
1072                        << pumpPrint
1073                        << CoinMessageEol;
1074                        if (newSolutionValue < solutionValue) {
1075                            memcpy(betterSolution, newSolution, numberColumns*sizeof(double));
1076                            CoinWarmStartBasis * basis =
1077                                dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ;
1078                            solutionFound = true;
1079                            numberSolutions++;
1080                            if (numberSolutions>=maxSolutions)
1081                              exitAll = true;
1082                            if (exitNow(newSolutionValue))
1083                                exitAll = true;
1084                            if (basis) {
1085                                bestBasis = * basis;
1086                                delete basis;
1087                                int action = model_->dealWithEventHandler(CbcEventHandler::heuristicSolution, newSolutionValue, betterSolution);
1088                                if (!action) {
1089                                    double * saveOldSolution = CoinCopyOfArray(model_->bestSolution(), numberColumns);
1090                                    double saveObjectiveValue = model_->getMinimizationObjValue();
1091                                    model_->setBestSolution(betterSolution, numberColumns, newSolutionValue);
1092                                    if (saveOldSolution && saveObjectiveValue < model_->getMinimizationObjValue())
1093                                        model_->setBestSolution(saveOldSolution, numberColumns, saveObjectiveValue);
1094                                    delete [] saveOldSolution;
1095                                }
1096                                if (!action || model_->maximumSecondsReached()) {
1097                                    exitAll = true; // exit
1098                                    break;
1099                                }
1100                            }
1101                            if ((accumulate_&1) != 0) {
1102                                model_->incrementUsed(betterSolution); // for local search
1103                            }
1104                            solutionValue = newSolutionValue;
1105                            solutionFound = true;
1106                            numberSolutions++;
1107                            if (numberSolutions>=maxSolutions)
1108                              exitAll = true;
1109                            if (exitNow(newSolutionValue))
1110                                exitAll = true;
1111                        } else {
1112                            returnCode = 0;
1113                        }
1114                    }
1115                }
1116                int numberIterations = 0;
1117                if (!doGeneral) {
1118                    // faster to do from all slack!!!!
1119                    if (allSlack) {
1120                        CoinWarmStartBasis dummy;
1121                        solver->setWarmStart(&dummy);
1122                    }
1123#ifdef COIN_DEVELOP
1124                    printf("%d perturbed out of %d columns (%d changed)\n", numberPerturbed, numberColumns, numberChanged);
1125#endif
1126                    bool takeHint;
1127                    OsiHintStrength strength;
1128                    solver->getHintParam(OsiDoDualInResolve, takeHint, strength);
1129                    if (dualPass && numberChanged > 2) {
1130                        solver->setHintParam(OsiDoDualInResolve, true); // dual may be better
1131                        if (dualPass == 1 && 2*numberChanged < numberColumns &&
1132                                (numberChanged < 5000 || 6*numberChanged < numberColumns)) {
1133                            // but we need to make infeasible
1134                            CoinWarmStartBasis * basis =
1135                                dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ;
1136                            if (basis) {
1137                                // modify
1138                                const double * lower = solver->getColLower();
1139                                const double * upper = solver->getColUpper();
1140                                double * solution = CoinCopyOfArray(solver->getColSolution(),
1141                                                                    numberColumns);
1142                                const double * objective = solver->getObjCoefficients();
1143                                int nChanged = 0;
1144                                for (i = 0; i < numberIntegersOrig; i++) {
1145                                    int iColumn = integerVariableOrig[i];
1146#ifdef RAND_RAND
1147                                    if (nChanged > numberChanged)
1148                                        break;
1149#endif
1150                                    if (objective[iColumn] > 0.0) {
1151                                        if (basis->getStructStatus(iColumn) ==
1152                                                CoinWarmStartBasis::atUpperBound) {
1153                                            solution[iColumn] = lower[iColumn];
1154                                            basis->setStructStatus(iColumn, CoinWarmStartBasis::atLowerBound);
1155                                            nChanged++;
1156                                        }
1157                                    } else if (objective[iColumn] < 0.0) {
1158                                        if (basis->getStructStatus(iColumn) ==
1159                                                CoinWarmStartBasis::atLowerBound) {
1160                                            solution[iColumn] = upper[iColumn];
1161                                            basis->setStructStatus(iColumn, CoinWarmStartBasis::atUpperBound);
1162                                            nChanged++;
1163                                        }
1164                                    }
1165                                }
1166                                if (!nChanged) {
1167                                    for (i = 0; i < numberIntegersOrig; i++) {
1168                                        int iColumn = integerVariableOrig[i];
1169                                        if (objective[iColumn] > 0.0) {
1170                                            if (basis->getStructStatus(iColumn) ==
1171                                                    CoinWarmStartBasis::basic) {
1172                                                solution[iColumn] = lower[iColumn];
1173                                                basis->setStructStatus(iColumn, CoinWarmStartBasis::atLowerBound);
1174                                                break;
1175                                            }
1176                                        } else if (objective[iColumn] < 0.0) {
1177                                            if (basis->getStructStatus(iColumn) ==
1178                                                    CoinWarmStartBasis::basic) {
1179                                                solution[iColumn] = upper[iColumn];
1180                                                basis->setStructStatus(iColumn, CoinWarmStartBasis::atUpperBound);
1181                                                break;
1182                                            }
1183                                        }
1184                                    }
1185                                }
1186                                solver->setColSolution(solution);
1187                                delete [] solution;
1188                                solver->setWarmStart(basis);
1189                                delete basis;
1190                            }
1191                        } else {
1192                            // faster to do from all slack!!!! ???
1193                            CoinWarmStartBasis dummy;
1194                            solver->setWarmStart(&dummy);
1195                        }
1196                    }
1197                    if (numberTries > 1 && numberPasses == 1 && firstPerturbedObjective) {
1198                        // Modify to use convex combination
1199                        // use basis from first time
1200                        solver->setWarmStart(&saveBasis);
1201                        // and objective
1202                        if (secondPassOpt < 3 || (secondPassOpt >= 4 && secondPassOpt < 6))
1203                            solver->setObjective(firstPerturbedObjective);
1204                        // and solution
1205                        solver->setColSolution(firstPerturbedSolution);
1206                        //if (secondPassOpt==2||secondPassOpt==5||
1207                        if (firstPerturbedValue > cutoff)
1208                            solver->setHintParam(OsiDoDualInResolve, true); // dual may be better
1209                    }
1210                    solver->resolve();
1211                    if (!solver->isProvenOptimal()) {
1212                        // presumably max time or some such
1213                        exitAll = true;
1214                        break;
1215                    }
1216                    solver->setHintParam(OsiDoDualInResolve, takeHint);
1217                    newTrueSolutionValue = -saveOffset;
1218                    newSumInfeas = 0.0;
1219                    newNumberInfeas = 0;
1220                    {
1221                        const double * newSolution = solver->getColSolution();
1222                        for (  i = 0 ; i < numberColumns ; i++ ) {
1223                            if (solver->isInteger(i)) {
1224                                double value = newSolution[i];
1225                                double nearest = floor(value + 0.5);
1226                                newSumInfeas += fabs(value - nearest);
1227                                if (fabs(value - nearest) > 1.0e-6)
1228                                    newNumberInfeas++;
1229                            }
1230                            newTrueSolutionValue += saveObjective[i] * newSolution[i];
1231                        }
1232                        newTrueSolutionValue *= direction;
1233                        if (numberPasses == 1 && secondPassOpt) {
1234                            if (numberTries == 1 || secondPassOpt > 3) {
1235                                // save basis
1236                                CoinWarmStartBasis * basis =
1237                                    dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ;
1238                                if (basis) {
1239                                    saveBasis = * basis;
1240                                    delete basis;
1241                                }
1242                                delete [] firstPerturbedObjective;
1243                                delete [] firstPerturbedSolution;
1244                                firstPerturbedObjective = CoinCopyOfArray(solver->getObjCoefficients(), numberColumns);
1245                                firstPerturbedSolution = CoinCopyOfArray(solver->getColSolution(), numberColumns);
1246                                firstPerturbedValue = newTrueSolutionValue;
1247                            }
1248                        }
1249                        if (newNumberInfeas && newNumberInfeas < 15) {
1250#ifdef JJF_ZERO
1251                            roundingObjective = solutionValue;
1252                            OsiSolverInterface * saveSolver = model_->swapSolver(solver);
1253                            double * currentObjective =
1254                                CoinCopyOfArray(solver->getObjCoefficients(), numberColumns);
1255                            solver->setObjective(saveObjective);
1256                            double saveOffset2;
1257                            solver->getDblParam(OsiObjOffset, saveOffset2);
1258                            solver->setDblParam(OsiObjOffset, saveOffset);
1259                            int ifSol = roundingHeuristic.solution(roundingObjective, roundingSolution);
1260                            solver->setObjective(currentObjective);
1261                            solver->setDblParam(OsiObjOffset, saveOffset2);
1262                            delete [] currentObjective;
1263                            model_->swapSolver(saveSolver);
1264                            if (ifSol > 0)
1265                                abort();
1266#endif
1267                            int numberRows = solver->getNumRows();
1268                            double * rowActivity = new double[numberRows];
1269                            memset(rowActivity, 0, numberRows*sizeof(double));
1270                            int * which = new int[newNumberInfeas];
1271                            int * stack = new int[newNumberInfeas+1];
1272                            double * baseValue = new double[newNumberInfeas];
1273                            int * whichRow = new int[numberRows];
1274                            double * rowValue = new double[numberRows];
1275                            memset(rowValue, 0, numberRows*sizeof(double));
1276                            int nRow = 0;
1277                            // Column copy
1278                            const double * element = solver->getMatrixByCol()->getElements();
1279                            const int * row = solver->getMatrixByCol()->getIndices();
1280                            const CoinBigIndex * columnStart = solver->getMatrixByCol()->getVectorStarts();
1281                            const int * columnLength = solver->getMatrixByCol()->getVectorLengths();
1282                            int n = 0;
1283                            double contrib = 0.0;
1284                            for (  i = 0 ; i < numberColumns ; i++ ) {
1285                                double value = newSolution[i];
1286                                if (solver->isInteger(i)) {
1287                                    double nearest = floor(value + 0.5);
1288                                    if (fabs(value - nearest) > 1.0e-6) {
1289                                        //printf("Column %d value %g\n",i,value);
1290                                        for (CoinBigIndex j = columnStart[i];
1291                                                j < columnStart[i] + columnLength[i]; j++) {
1292                                            int iRow = row[j];
1293                                            //printf("row %d element %g\n",iRow,element[j]);
1294                                            if (!rowValue[iRow]) {
1295                                                rowValue[iRow] = 1.0;
1296                                                whichRow[nRow++] = iRow;
1297                                            }
1298                                        }
1299                                        baseValue[n] = floor(value);
1300                                        contrib += saveObjective[i] * value;
1301                                        value = 0.0;
1302                                        stack[n] = 0;
1303                                        which[n++] = i;
1304                                    }
1305                                }
1306                                for (CoinBigIndex j = columnStart[i];
1307                                        j < columnStart[i] + columnLength[i]; j++) {
1308                                    int iRow = row[j];
1309                                    rowActivity[iRow] += value * element[j];
1310                                }
1311                            }
1312                            if (newNumberInfeas < 15) {
1313                                stack[n] = newNumberInfeas + 100;
1314                                int iStack = n;
1315                                memset(rowValue, 0, numberRows*sizeof(double));
1316                                const double * rowLower = solver->getRowLower();
1317                                const double * rowUpper = solver->getRowUpper();
1318                                while (iStack >= 0) {
1319                                    double contrib2 = 0.0;
1320                                    // Could do faster
1321                                    for (int k = 0 ; k < n ; k++ ) {
1322                                        i = which[k];
1323                                        double value = baseValue[k] + stack[k];
1324                                        contrib2 += saveObjective[i] * value;
1325                                        for (CoinBigIndex j = columnStart[i];
1326                                                j < columnStart[i] + columnLength[i]; j++) {
1327                                            int iRow = row[j];
1328                                            rowValue[iRow] += value * element[j];
1329                                        }
1330                                    }
1331                                    // check if feasible
1332                                    bool feasible = true;
1333                                    for (int k = 0; k < nRow; k++) {
1334                                        i = whichRow[k];
1335                                        double value = rowValue[i] + rowActivity[i];
1336                                        rowValue[i] = 0.0;
1337                                        if (value < rowLower[i] - 1.0e-7 ||
1338                                                value > rowUpper[i] + 1.0e-7)
1339                                            feasible = false;
1340                                    }
1341                                    if (feasible) {
1342                                        double newObj = newTrueSolutionValue * direction;
1343                                        newObj += contrib2 - contrib;
1344                                        newObj *= direction;
1345#ifdef COIN_DEVELOP
1346                                        printf("FFFeasible! - obj %g\n", newObj);
1347#endif
1348                                        if (newObj < roundingObjective - 1.0e-6) {
1349#ifdef COIN_DEVELOP
1350                                            printf("FBetter\n");
1351#endif
1352                                            roundingObjective = newObj;
1353                                            memcpy(roundingSolution, newSolution, numberColumns*sizeof(double));
1354                                            for (int k = 0 ; k < n ; k++ ) {
1355                                                i = which[k];
1356                                                double value = baseValue[k] + stack[k];
1357                                                roundingSolution[i] = value;
1358                                            }
1359                                        }
1360                                    }
1361                                    while (iStack >= 0 && stack[iStack]) {
1362                                        stack[iStack]--;
1363                                        iStack--;
1364                                    }
1365                                    if (iStack >= 0) {
1366                                        stack[iStack] = 1;
1367                                        iStack = n;
1368                                        stack[n] = 1;
1369                                    }
1370                                }
1371                            }
1372                            delete [] rowActivity;
1373                            delete [] which;
1374                            delete [] stack;
1375                            delete [] baseValue;
1376                            delete [] whichRow;
1377                            delete [] rowValue;
1378                        }
1379                    }
1380                    if (true) {
1381                        OsiSolverInterface * saveSolver = model_->swapSolver(clonedSolver);
1382                        clonedSolver->setColSolution(solver->getColSolution());
1383                        CbcRounding heuristic1(*model_);
1384                        heuristic1.setHeuristicName("rounding in feaspump!");
1385                        heuristic1.setWhen(1);
1386                        roundingObjective = CoinMin(roundingObjective, solutionValue);
1387                        double testSolutionValue = newTrueSolutionValue;
1388                        int returnCode = heuristic1.solution(roundingObjective,
1389                                                             roundingSolution,
1390                                                             testSolutionValue) ;
1391                        if (returnCode == 1) {
1392#ifdef COIN_DEVELOP
1393                            printf("rounding obj of %g?\n", roundingObjective);
1394#endif
1395                            //roundingObjective = newSolutionValue;
1396                            //} else {
1397                            //roundingObjective = COIN_DBL_MAX;
1398                        }
1399                        model_->swapSolver(saveSolver);
1400                    }
1401                    if (!solver->isProvenOptimal()) {
1402                        // presumably max time or some such
1403                        exitAll = true;
1404                        break;
1405                    }
1406                    // in case very dubious solver
1407                    lower = solver->getColLower();
1408                    upper = solver->getColUpper();
1409                    solution = solver->getColSolution();
1410                    numberIterations = solver->getIterationCount();
1411                } else {
1412                    int * addStart = new int[2*general+1];
1413                    int * addIndex = new int[4*general];
1414                    double * addElement = new double[4*general];
1415                    double * addLower = new double[2*general];
1416                    double * addUpper = new double[2*general];
1417                    double * obj = new double[general];
1418                    int nAdd = 0;
1419                    for (i = 0; i < numberIntegers; i++) {
1420                        int iColumn = integerVariable[i];
1421                        if (newSolution[iColumn] > lower[iColumn] + primalTolerance &&
1422                                newSolution[iColumn] < upper[iColumn] - primalTolerance) {
1423                            assert (upper[iColumn] - lower[iColumn] > 1.00001);
1424                            obj[nAdd] = 1.0;
1425                            addLower[nAdd] = 0.0;
1426                            addUpper[nAdd] = COIN_DBL_MAX;
1427                            nAdd++;
1428                        }
1429                    }
1430                    OsiSolverInterface * solver2 = solver;
1431                    if (nAdd) {
1432                        CoinZeroN(addStart, nAdd + 1);
1433                        solver2 = solver->clone();
1434                        solver2->addCols(nAdd, addStart, NULL, NULL, addLower, addUpper, obj);
1435                        // feasible solution
1436                        double * sol = new double[nAdd+numberColumns];
1437                        memcpy(sol, solution, numberColumns*sizeof(double));
1438                        // now rows
1439                        int nAdd = 0;
1440                        int nEl = 0;
1441                        int nAddRow = 0;
1442                        for (i = 0; i < numberIntegers; i++) {
1443                            int iColumn = integerVariable[i];
1444                            if (newSolution[iColumn] > lower[iColumn] + primalTolerance &&
1445                                    newSolution[iColumn] < upper[iColumn] - primalTolerance) {
1446                                addLower[nAddRow] = -newSolution[iColumn];;
1447                                addUpper[nAddRow] = COIN_DBL_MAX;
1448                                addIndex[nEl] = iColumn;
1449                                addElement[nEl++] = -1.0;
1450                                addIndex[nEl] = numberColumns + nAdd;
1451                                addElement[nEl++] = 1.0;
1452                                nAddRow++;
1453                                addStart[nAddRow] = nEl;
1454                                addLower[nAddRow] = newSolution[iColumn];;
1455                                addUpper[nAddRow] = COIN_DBL_MAX;
1456                                addIndex[nEl] = iColumn;
1457                                addElement[nEl++] = 1.0;
1458                                addIndex[nEl] = numberColumns + nAdd;
1459                                addElement[nEl++] = 1.0;
1460                                nAddRow++;
1461                                addStart[nAddRow] = nEl;
1462                                sol[nAdd+numberColumns] = fabs(sol[iColumn] - newSolution[iColumn]);
1463                                nAdd++;
1464                            }
1465                        }
1466                        solver2->setColSolution(sol);
1467                        delete [] sol;
1468                        solver2->addRows(nAddRow, addStart, addIndex, addElement, addLower, addUpper);
1469                    }
1470                    delete [] addStart;
1471                    delete [] addIndex;
1472                    delete [] addElement;
1473                    delete [] addLower;
1474                    delete [] addUpper;
1475                    delete [] obj;
1476                    solver2->resolve();
1477                    if (!solver2->isProvenOptimal()) {
1478                        // presumably max time or some such
1479                        exitAll = true;
1480                        break;
1481                    }
1482                    //assert (solver2->isProvenOptimal());
1483                    if (nAdd) {
1484                        solver->setColSolution(solver2->getColSolution());
1485                        numberIterations = solver2->getIterationCount();
1486                        delete solver2;
1487                    } else {
1488                        numberIterations = solver->getIterationCount();
1489                    }
1490                    lower = solver->getColLower();
1491                    upper = solver->getColUpper();
1492                    solution = solver->getColSolution();
1493                    newTrueSolutionValue = -saveOffset;
1494                    newSumInfeas = 0.0;
1495                    newNumberInfeas = 0;
1496                    {
1497                        const double * newSolution = solver->getColSolution();
1498                        for (  i = 0 ; i < numberColumns ; i++ ) {
1499                            if (solver->isInteger(i)) {
1500                                double value = newSolution[i];
1501                                double nearest = floor(value + 0.5);
1502                                newSumInfeas += fabs(value - nearest);
1503                                if (fabs(value - nearest) > 1.0e-6)
1504                                    newNumberInfeas++;
1505                            }
1506                            newTrueSolutionValue += saveObjective[i] * newSolution[i];
1507                        }
1508                        newTrueSolutionValue *= direction;
1509                    }
1510                }
1511                if (lastMove != 1000000) {
1512                    if (newSumInfeas < lastSumInfeas) {
1513                        lastMove = numberPasses;
1514                        lastSumInfeas = newSumInfeas;
1515                    } else if (newSumInfeas > lastSumInfeas + 1.0e-5) {
1516                        lastMove = 1000000; // going up
1517                    }
1518                }
1519                totalNumberIterations += numberIterations;
1520                if (solver->getNumRows() < 3000)
1521                    sprintf(pumpPrint, "Pass %3d: suminf. %10.5f (%d) obj. %g iterations %d",
1522                            numberPasses + totalNumberPasses,
1523                            newSumInfeas, newNumberInfeas,
1524                            newTrueSolutionValue, numberIterations);
1525                else
1526                    sprintf(pumpPrint, "Pass %3d: (%.2f seconds) suminf. %10.5f (%d) obj. %g iterations %d", numberPasses + totalNumberPasses,
1527                            model_->getCurrentSeconds(), newSumInfeas, newNumberInfeas,
1528                            newTrueSolutionValue, numberIterations);
1529                model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
1530                << pumpPrint
1531                << CoinMessageEol;
1532                if (closestSolution && solver->getObjValue() < closestObjectiveValue) {
1533                    int i;
1534                    const double * objective = solver->getObjCoefficients();
1535                    for (i = 0; i < numberIntegersOrig; i++) {
1536                        int iColumn = integerVariableOrig[i];
1537                        if (objective[iColumn] > 0.0)
1538                            closestSolution[i] = 0;
1539                        else
1540                            closestSolution[i] = 1;
1541                    }
1542                    closestObjectiveValue = solver->getObjValue();
1543                }
1544                // See if we need to think about changing rhs
1545                if ((switches_&12) != 0 && useRhs < 1.0e50) {
1546                    double oldRhs = useRhs;
1547                    bool trying = false;
1548                    if ((switches_&4) != 0 && numberPasses && (numberPasses % 50) == 0) {
1549                        if (solutionValue > 1.0e20) {
1550                            // only if no genuine solution
1551                            double gap = useRhs - continuousObjectiveValue;
1552                            useRhs += 0.1 * gap;
1553                            if (exactMultiple) {
1554                                useRhs = exactMultiple * ceil(useRhs / exactMultiple);
1555                                useRhs = CoinMax(useRhs, oldRhs + exactMultiple);
1556                            }
1557                            trying = true;
1558                        }
1559                    }
1560                    if ((switches_&8) != 0) {
1561                        // Put in new suminf and check
1562                        double largest = newSumInfeas;
1563                        double smallest = newSumInfeas;
1564                        for (int i = 0; i < SIZE_BOBBLE - 1; i++) {
1565                            double value = saveSumInf[i+1];
1566                            saveSumInf[i] = value;
1567                            largest = CoinMax(largest, value);
1568                            smallest = CoinMin(smallest, value);
1569                        }
1570                        saveSumInf[SIZE_BOBBLE-1] = newSumInfeas;
1571                        if (smallest*1.5 > largest && smallest > 2.0) {
1572                            if (bobbleMode == 0) {
1573                                // go closer
1574                                double gap = oldRhs - continuousObjectiveValue;
1575                                useRhs -= 0.4 * gap;
1576                                if (exactMultiple) {
1577                                    double value = floor(useRhs / exactMultiple);
1578                                    useRhs = CoinMin(value * exactMultiple, oldRhs - exactMultiple);
1579                                }
1580                                if (useRhs < continuousObjectiveValue) {
1581                                    // skip decrease
1582                                    bobbleMode = 1;
1583                                    useRhs = oldRhs;
1584                                }
1585                            }
1586                            if (bobbleMode) {
1587                                trying = true;
1588                                // weaken
1589                                if (solutionValue < 1.0e20) {
1590                                    double gap = solutionValue - oldRhs;
1591                                    useRhs += 0.3 * gap;
1592                                } else {
1593                                    double gap = oldRhs - continuousObjectiveValue;
1594                                    useRhs += 0.05 * gap;
1595                                }
1596                                if (exactMultiple) {
1597                                    double value = ceil(useRhs / exactMultiple);
1598                                    useRhs = CoinMin(value * exactMultiple,
1599                                                     solutionValue - exactMultiple);
1600                                }
1601                            }
1602                            bobbleMode++;
1603                            // reset
1604                            CoinFillN(saveSumInf, SIZE_BOBBLE, COIN_DBL_MAX);
1605                        }
1606                    }
1607                    if (useRhs != oldRhs) {
1608                        // tidy up
1609                        if (exactMultiple) {
1610                            double value = floor(useRhs / exactMultiple);
1611                            double bestPossible = ceil(continuousObjectiveValue / exactMultiple);
1612                            useRhs = CoinMax(value, bestPossible) * exactMultiple;
1613                        } else {
1614                            useRhs = CoinMax(useRhs, continuousObjectiveValue);
1615                        }
1616                        int k = solver->getNumRows() - 1;
1617                        solver->setRowUpper(k, useRhs + useOffset);
1618                        bool takeHint;
1619                        OsiHintStrength strength;
1620                        solver->getHintParam(OsiDoDualInResolve, takeHint, strength);
1621                        if (useRhs < oldRhs) {
1622                            solver->setHintParam(OsiDoDualInResolve, true);
1623                            solver->resolve();
1624                        } else if (useRhs > oldRhs) {
1625                            solver->setHintParam(OsiDoDualInResolve, false);
1626                            solver->resolve();
1627                        }
1628                        solver->setHintParam(OsiDoDualInResolve, takeHint);
1629                        if (!solver->isProvenOptimal()) {
1630                            // presumably max time or some such
1631                            exitAll = true;
1632                            break;
1633                        }
1634                    } else if (trying) {
1635                        // doesn't look good
1636                        break;
1637                    }
1638                }
1639            }
1640            // reduce scale factor
1641            scaleFactor *= weightFactor_;
1642        } // END WHILE
1643        // see if rounding worked!
1644        if (roundingObjective < solutionValue) {
1645            if (roundingObjective < solutionValue - 1.0e-6*fabs(roundingObjective)) {
1646                sprintf(pumpPrint, "Rounding solution of %g is better than previous of %g\n",
1647                        roundingObjective, solutionValue);
1648                model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
1649                << pumpPrint
1650                << CoinMessageEol;
1651            }
1652            solutionValue = roundingObjective;
1653            newSolutionValue = solutionValue;
1654            memcpy(betterSolution, roundingSolution, numberColumns*sizeof(double));
1655            solutionFound = true;
1656            numberSolutions++;
1657            if (numberSolutions>=maxSolutions)
1658              exitAll = true;
1659            if (exitNow(roundingObjective))
1660                exitAll = true;
1661        }
1662        if (!solutionFound) {
1663            sprintf(pumpPrint, "No solution found this major pass");
1664            model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
1665            << pumpPrint
1666            << CoinMessageEol;
1667        }
1668        //}
1669        delete solver;
1670        solver = NULL;
1671        for ( j = 0; j < NUMBER_OLD; j++)
1672            delete [] oldSolution[j];
1673        delete [] oldSolution;
1674        delete [] saveObjective;
1675        if (usedColumn && !exitAll) {
1676            OsiSolverInterface * newSolver = cloneBut(3); // was model_->continuousSolver()->clone();
1677#if 0 //def COIN_HAS_CLP
1678            OsiClpSolverInterface * clpSolver
1679              = dynamic_cast<OsiClpSolverInterface *> (newSolver);
1680            if (clpSolver) {
1681              ClpSimplex * simplex = clpSolver->getModelPtr();
1682              simplex->writeMps("start.mps",2,1);
1683            }
1684#endif
1685            const double * colLower = newSolver->getColLower();
1686            const double * colUpper = newSolver->getColUpper();
1687            bool stopBAB = false;
1688            int allowedPass = -1;
1689            if (maximumAllowed > 0)
1690                allowedPass = CoinMax(numberPasses - maximumAllowed, -1);
1691            while (!stopBAB) {
1692                stopBAB = true;
1693                int i;
1694                int nFix = 0;
1695                int nFixI = 0;
1696                int nFixC = 0;
1697                int nFixC2 = 0;
1698                for (i = 0; i < numberIntegersOrig; i++) {
1699                    int iColumn = integerVariableOrig[i];
1700                    //const OsiObject * object = model_->object(i);
1701                    //double originalLower;
1702                    //double originalUpper;
1703                    //getIntegerInformation( object,originalLower, originalUpper);
1704                    //assert(colLower[iColumn]==originalLower);
1705                    //newSolver->setColLower(iColumn,CoinMax(colLower[iColumn],originalLower));
1706                    newSolver->setColLower(iColumn, colLower[iColumn]);
1707                    //assert(colUpper[iColumn]==originalUpper);
1708                    //newSolver->setColUpper(iColumn,CoinMin(colUpper[iColumn],originalUpper));
1709                    newSolver->setColUpper(iColumn, colUpper[iColumn]);
1710                    if (usedColumn[iColumn] <= allowedPass) {
1711                        double value = lastSolution[iColumn];
1712                        double nearest = floor(value + 0.5);
1713                        if (fabs(value - nearest) < 1.0e-7) {
1714                            if (nearest == colLower[iColumn]) {
1715                                newSolver->setColUpper(iColumn, colLower[iColumn]);
1716                                nFix++;
1717                            } else if (nearest == colUpper[iColumn]) {
1718                                newSolver->setColLower(iColumn, colUpper[iColumn]);
1719                                nFix++;
1720                            } else if (fixInternal) {
1721                                newSolver->setColLower(iColumn, nearest);
1722                                newSolver->setColUpper(iColumn, nearest);
1723                                nFix++;
1724                                nFixI++;
1725                            }
1726                        }
1727                    }
1728                }
1729                if (fixContinuous) {
1730                    for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
1731                        if (!newSolver->isInteger(iColumn) && usedColumn[iColumn] <= allowedPass) {
1732                            double value = lastSolution[iColumn];
1733                            if (value < colLower[iColumn] + 1.0e-8) {
1734                                newSolver->setColUpper(iColumn, colLower[iColumn]);
1735                                nFixC++;
1736                            } else if (value > colUpper[iColumn] - 1.0e-8) {
1737                                newSolver->setColLower(iColumn, colUpper[iColumn]);
1738                                nFixC++;
1739                            } else if (fixContinuous == 2) {
1740                                newSolver->setColLower(iColumn, value);
1741                                newSolver->setColUpper(iColumn, value);
1742                                nFixC++;
1743                                nFixC2++;
1744                            }
1745                        }
1746                    }
1747                }
1748                newSolver->initialSolve();
1749                if (!newSolver->isProvenOptimal()) {
1750                    //newSolver->writeMps("bad.mps");
1751                    //assert (newSolver->isProvenOptimal());
1752                    exitAll = true;
1753                    break;
1754                }
1755                sprintf(pumpPrint, "Before mini branch and bound, %d integers at bound fixed and %d continuous",
1756                        nFix, nFixC);
1757                if (nFixC2 + nFixI != 0)
1758                    sprintf(pumpPrint + strlen(pumpPrint), " of which %d were internal integer and %d internal continuous",
1759                            nFixI, nFixC2);
1760                model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
1761                << pumpPrint
1762                << CoinMessageEol;
1763                double saveValue = newSolutionValue;
1764                if (newSolutionValue - model_->getCutoffIncrement()
1765                        > continuousObjectiveValue - 1.0e-7) {
1766                    double saveFraction = fractionSmall_;
1767                    if (numberTries > 1 && !numberBandBsolutions)
1768                        fractionSmall_ *= 0.5;
1769                    // Give branch and bound a bit more freedom
1770                    double cutoff2 = newSolutionValue +
1771                                     CoinMax(model_->getCutoffIncrement(), 1.0e-3);
1772#if 0
1773                      {
1774                        OsiClpSolverInterface * clpSolver
1775                        = dynamic_cast<OsiClpSolverInterface *> (newSolver);
1776                        if (clpSolver) {
1777                            ClpSimplex * simplex = clpSolver->getModelPtr();
1778                            simplex->writeMps("testA.mps",2,1);
1779                        }
1780                      }
1781#endif
1782                    int returnCode2 = smallBranchAndBound(newSolver, numberNodes_, newSolution, newSolutionValue,
1783                                                          cutoff2, "CbcHeuristicLocalAfterFPump");
1784                    fractionSmall_ = saveFraction;
1785                    if (returnCode2 < 0) {
1786                        if (returnCode2 == -2) {
1787                            exitAll = true;
1788                            returnCode = 0;
1789                        } else {
1790                            returnCode2 = 0; // returned on size - try changing
1791                            //#define ROUND_AGAIN
1792#ifdef ROUND_AGAIN
1793                            if (numberTries == 1 && numberPasses > 20 && allowedPass < numberPasses - 1) {
1794                                allowedPass = (numberPasses + allowedPass) >> 1;
1795                                sprintf(pumpPrint,
1796                                        "Fixing all variables which were last changed on pass %d and trying again",
1797                                        allowedPass);
1798                                model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
1799                                << pumpPrint
1800                                << CoinMessageEol;
1801                                stopBAB = false;
1802                                continue;
1803                            }
1804#endif
1805                        }
1806                    }
1807                    if ((returnCode2&2) != 0) {
1808                        // could add cut
1809                        returnCode2 &= ~2;
1810                    }
1811                    if (returnCode2) {
1812                        numberBandBsolutions++;
1813                        // may not have got solution earlier
1814                        returnCode |= 1;
1815                    }
1816                } else {
1817                    // no need
1818                    exitAll = true;
1819                    //returnCode=0;
1820                }
1821                // recompute solution value
1822                if (returnCode && true) {
1823#if 0
1824                      {
1825                        OsiClpSolverInterface * clpSolver
1826                        = dynamic_cast<OsiClpSolverInterface *> (newSolver);
1827                        if (clpSolver) {
1828                            ClpSimplex * simplex = clpSolver->getModelPtr();
1829                            simplex->writeMps("testB.mps",2,1);
1830                        }
1831                      }
1832#endif
1833                    delete newSolver;
1834                    newSolver = cloneBut(3); // was model_->continuousSolver()->clone();
1835                    newSolutionValue = -saveOffset;
1836                    double newSumInfeas = 0.0;
1837                    const double * obj = newSolver->getObjCoefficients();
1838                    for (int i = 0 ; i < numberColumns ; i++ ) {
1839                        if (newSolver->isInteger(i)) {
1840                            double value = newSolution[i];
1841                            double nearest = floor(value + 0.5);
1842                            newSumInfeas += fabs(value - nearest);
1843                        }
1844                        newSolutionValue += obj[i] * newSolution[i];
1845                    }
1846                    newSolutionValue *= direction;
1847                }
1848                bool gotSolution = false;
1849                if (returnCode && newSolutionValue < saveValue) {
1850                    sprintf(pumpPrint, "Mini branch and bound improved solution from %g to %g (%.2f seconds)",
1851                            saveValue, newSolutionValue, model_->getCurrentSeconds());
1852                    model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
1853                    << pumpPrint
1854                    << CoinMessageEol;
1855                    memcpy(betterSolution, newSolution, numberColumns*sizeof(double));
1856                    gotSolution = true;
1857                    if (fixContinuous && nFixC + nFixC2 > 0) {
1858                        // may be able to do even better
1859                        int nFixed = 0;
1860                        const double * lower = model_->solver()->getColLower();
1861                        const double * upper = model_->solver()->getColUpper();
1862                        for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
1863                            double value = newSolution[iColumn];
1864                            if (newSolver->isInteger(iColumn)) {
1865                                value = floor(newSolution[iColumn] + 0.5);
1866                                newSolver->setColLower(iColumn, value);
1867                                newSolver->setColUpper(iColumn, value);
1868                                nFixed++;
1869                            } else {
1870                                newSolver->setColLower(iColumn, lower[iColumn]);
1871                                newSolver->setColUpper(iColumn, upper[iColumn]);
1872                                if (value < lower[iColumn])
1873                                    value = lower[iColumn];
1874                                else if (value > upper[iColumn])
1875                                    value = upper[iColumn];
1876
1877                            }
1878                            newSolution[iColumn] = value;
1879                        }
1880                        newSolver->setColSolution(newSolution);
1881                        //#define CLP_INVESTIGATE2
1882#ifdef CLP_INVESTIGATE2
1883                        {
1884                            // check
1885                            // get row activities
1886                            int numberRows = newSolver->getNumRows();
1887                            double * rowActivity = new double[numberRows];
1888                            memset(rowActivity, 0, numberRows*sizeof(double));
1889                            newSolver->getMatrixByCol()->times(newSolution, rowActivity) ;
1890                            double largestInfeasibility = primalTolerance;
1891                            double sumInfeasibility = 0.0;
1892                            int numberBadRows = 0;
1893                            const double * rowLower = newSolver->getRowLower();
1894                            const double * rowUpper = newSolver->getRowUpper();
1895                            for (i = 0 ; i < numberRows ; i++) {
1896                                double value;
1897                                value = rowLower[i] - rowActivity[i];
1898                                if (value > primalTolerance) {
1899                                    numberBadRows++;
1900                                    largestInfeasibility = CoinMax(largestInfeasibility, value);
1901                                    sumInfeasibility += value;
1902                                }
1903                                value = rowActivity[i] - rowUpper[i];
1904                                if (value > primalTolerance) {
1905                                    numberBadRows++;
1906                                    largestInfeasibility = CoinMax(largestInfeasibility, value);
1907                                    sumInfeasibility += value;
1908                                }
1909                            }
1910                            printf("%d bad rows, largest inf %g sum %g\n",
1911                                   numberBadRows, largestInfeasibility, sumInfeasibility);
1912                            delete [] rowActivity;
1913                        }
1914#endif
1915#ifdef COIN_HAS_CLP
1916                        OsiClpSolverInterface * clpSolver
1917                        = dynamic_cast<OsiClpSolverInterface *> (newSolver);
1918                        if (clpSolver) {
1919                            ClpSimplex * simplex = clpSolver->getModelPtr();
1920                            //simplex->writeBasis("test.bas",true,2);
1921                            //simplex->writeMps("test.mps",2,1);
1922                            if (nFixed*3 > numberColumns*2)
1923                                simplex->allSlackBasis(); // may as well go from all slack
1924                            simplex->primal(1);
1925                            clpSolver->setWarmStart(NULL);
1926                        }
1927#endif
1928                        newSolver->initialSolve();
1929                        if (newSolver->isProvenOptimal()) {
1930                            double value = newSolver->getObjValue() * newSolver->getObjSense();
1931                            if (value < newSolutionValue) {
1932                                //newSolver->writeMps("query","mps");
1933#ifdef JJF_ZERO
1934                                {
1935                                    double saveOffset;
1936                                    newSolver->getDblParam(OsiObjOffset, saveOffset);
1937                                    const double * obj = newSolver->getObjCoefficients();
1938                                    double newTrueSolutionValue = -saveOffset;
1939                                    double newSumInfeas = 0.0;
1940                                    int numberColumns = newSolver->getNumCols();
1941                                    const double * solution = newSolver->getColSolution();
1942                                    for (int  i = 0 ; i < numberColumns ; i++ ) {
1943                                        if (newSolver->isInteger(i)) {
1944                                            double value = solution[i];
1945                                            double nearest = floor(value + 0.5);
1946                                            newSumInfeas += fabs(value - nearest);
1947                                        }
1948                                        if (solution[i])
1949                                            printf("%d obj %g val %g - total %g\n", i, obj[i], solution[i],
1950                                                   newTrueSolutionValue);
1951                                        newTrueSolutionValue += obj[i] * solution[i];
1952                                    }
1953                                    printf("obj %g - inf %g\n", newTrueSolutionValue,
1954                                           newSumInfeas);
1955                                }
1956#endif
1957                                sprintf(pumpPrint, "Freeing continuous variables gives a solution of %g", value);
1958                                model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
1959                                << pumpPrint
1960                                << CoinMessageEol;
1961                                newSolutionValue = value;
1962                                memcpy(betterSolution, newSolver->getColSolution(), numberColumns*sizeof(double));
1963                            }
1964                        } else {
1965                            //newSolver->writeMps("bad3.mps");
1966                          sprintf(pumpPrint, "On closer inspection solution is not valid");
1967                          model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
1968                            << pumpPrint
1969                            << CoinMessageEol;
1970                            exitAll = true;
1971                            break;
1972                        }
1973                    }
1974                } else {
1975                    sprintf(pumpPrint, "Mini branch and bound did not improve solution (%.2f seconds)",
1976                            model_->getCurrentSeconds());
1977                    model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
1978                    << pumpPrint
1979                    << CoinMessageEol;
1980                    if (returnCode && newSolutionValue < saveValue + 1.0e-3 && nFixC + nFixC2) {
1981                        // may be able to do better
1982                        const double * lower = model_->solver()->getColLower();
1983                        const double * upper = model_->solver()->getColUpper();
1984                        for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
1985                            if (newSolver->isInteger(iColumn)) {
1986                                double value = floor(newSolution[iColumn] + 0.5);
1987                                newSolver->setColLower(iColumn, value);
1988                                newSolver->setColUpper(iColumn, value);
1989                            } else {
1990                                newSolver->setColLower(iColumn, lower[iColumn]);
1991                                newSolver->setColUpper(iColumn, upper[iColumn]);
1992                            }
1993                        }
1994                        newSolver->initialSolve();
1995                        if (newSolver->isProvenOptimal()) {
1996                            double value = newSolver->getObjValue() * newSolver->getObjSense();
1997                            if (value < saveValue) {
1998                                sprintf(pumpPrint, "Freeing continuous variables gives a solution of %g", value);
1999                                model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
2000                                << pumpPrint
2001                                << CoinMessageEol;
2002                                newSolutionValue = value;
2003                                memcpy(betterSolution, newSolver->getColSolution(), numberColumns*sizeof(double));
2004                                gotSolution = true;
2005                            }
2006                        }
2007                    }
2008                }
2009                if (gotSolution) {
2010                    if ((accumulate_&1) != 0) {
2011                        model_->incrementUsed(betterSolution); // for local search
2012                    }
2013                    solutionValue = newSolutionValue;
2014                    solutionFound = true;
2015                    numberSolutions++;
2016                    if (numberSolutions>=maxSolutions)
2017                      exitAll = true;
2018                    if (exitNow(newSolutionValue))
2019                        exitAll = true;
2020                    CoinWarmStartBasis * basis =
2021                        dynamic_cast<CoinWarmStartBasis *>(newSolver->getWarmStart()) ;
2022                    if (basis) {
2023                        bestBasis = * basis;
2024                        delete basis;
2025                        int action = model_->dealWithEventHandler(CbcEventHandler::heuristicSolution, newSolutionValue, betterSolution);
2026                        if (action == 0) {
2027                            double * saveOldSolution = CoinCopyOfArray(model_->bestSolution(), numberColumns);
2028                            double saveObjectiveValue = model_->getMinimizationObjValue();
2029                            model_->setBestSolution(betterSolution, numberColumns, newSolutionValue);
2030                            if (saveOldSolution && saveObjectiveValue < model_->getMinimizationObjValue())
2031                                model_->setBestSolution(saveOldSolution, numberColumns, saveObjectiveValue);
2032                            delete [] saveOldSolution;
2033                        }
2034                        if (!action || model_->getCurrentSeconds() > model_->getMaximumSeconds()) {
2035                            exitAll = true; // exit
2036                            break;
2037                        }
2038                    }
2039                }
2040            } // end stopBAB while
2041            delete newSolver;
2042        }
2043        if (solutionFound) finalReturnCode = 1;
2044        cutoff = CoinMin(cutoff, solutionValue - model_->getCutoffIncrement());
2045        if (numberTries >= maximumRetries_ || !solutionFound || exitAll || cutoff < continuousObjectiveValue + 1.0e-7) {
2046            break;
2047        } else {
2048            solutionFound = false;
2049            if (absoluteIncrement_ > 0.0 || relativeIncrement_ > 0.0) {
2050                double gap = relativeIncrement_ * fabs(solutionValue);
2051                double change = CoinMax(gap, absoluteIncrement_);
2052                cutoff = CoinMin(cutoff, solutionValue - change);
2053            } else {
2054                //double weights[10]={0.1,0.1,0.2,0.2,0.2,0.3,0.3,0.3,0.4,0.5};
2055                double weights[10] = {0.1, 0.2, 0.3, 0.3, 0.4, 0.4, 0.4, 0.5, 0.5, 0.6};
2056                cutoff -= weights[CoinMin(numberTries-1, 9)] * (cutoff - continuousObjectiveValue);
2057            }
2058            // But round down
2059            if (exactMultiple)
2060                cutoff = exactMultiple * floor(cutoff / exactMultiple);
2061            if (cutoff < continuousObjectiveValue)
2062                break;
2063            sprintf(pumpPrint, "Round again with cutoff of %g", cutoff);
2064            model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
2065            << pumpPrint
2066            << CoinMessageEol;
2067            if ((accumulate_&3) < 2 && usedColumn) {
2068                for (int i = 0; i < numberColumns; i++)
2069                    usedColumn[i] = -1;
2070            }
2071            averageIterationsPerTry = totalNumberIterations / numberTries;
2072            numberIterationsLastPass =  totalNumberIterations;
2073            totalNumberPasses += numberPasses - 1;
2074        }
2075    }
2076/*
2077  End of the `exitAll' loop.
2078*/
2079#ifdef RAND_RAND
2080    delete [] randomFactor;
2081#endif
2082    delete solver; // probably NULL but do anyway
2083    if (!finalReturnCode && closestSolution && closestObjectiveValue <= 10.0 &&
2084            usedColumn && !model_->maximumSecondsReached()) {
2085        // try a bit of branch and bound
2086        OsiSolverInterface * newSolver = cloneBut(1); // was model_->continuousSolver()->clone();
2087        const double * colLower = newSolver->getColLower();
2088        const double * colUpper = newSolver->getColUpper();
2089        int i;
2090        double rhs = 0.0;
2091        for (i = 0; i < numberIntegersOrig; i++) {
2092            int iColumn = integerVariableOrig[i];
2093            int direction = closestSolution[i];
2094            closestSolution[i] = iColumn;
2095            if (direction == 0) {
2096                // keep close to LB
2097                rhs += colLower[iColumn];
2098                lastSolution[i] = 1.0;
2099            } else {
2100                // keep close to UB
2101                rhs -= colUpper[iColumn];
2102                lastSolution[i] = -1.0;
2103            }
2104        }
2105        newSolver->addRow(numberIntegersOrig, closestSolution,
2106                          lastSolution, -COIN_DBL_MAX, rhs + 10.0);
2107        //double saveValue = newSolutionValue;
2108        //newSolver->writeMps("sub");
2109        int returnCode = smallBranchAndBound(newSolver, numberNodes_, newSolution, newSolutionValue,
2110                                             newSolutionValue, "CbcHeuristicLocalAfterFPump");
2111        if (returnCode < 0)
2112            returnCode = 0; // returned on size
2113        if ((returnCode&2) != 0) {
2114            // could add cut
2115            returnCode &= ~2;
2116        }
2117        if (returnCode) {
2118            //printf("old sol of %g new of %g\n",saveValue,newSolutionValue);
2119            memcpy(betterSolution, newSolution, numberColumns*sizeof(double));
2120            //abort();
2121            solutionValue = newSolutionValue;
2122            solutionFound = true;
2123            numberSolutions++;
2124            if (numberSolutions>=maxSolutions)
2125              exitAll = true;
2126            if (exitNow(newSolutionValue))
2127                exitAll = true;
2128        }
2129        delete newSolver;
2130    }
2131    delete clonedSolver;
2132    delete [] roundingSolution;
2133    delete [] usedColumn;
2134    delete [] lastSolution;
2135    delete [] newSolution;
2136    delete [] closestSolution;
2137    delete [] integerVariable;
2138    delete [] firstPerturbedObjective;
2139    delete [] firstPerturbedSolution;
2140    if (solutionValue == incomingObjective)
2141        sprintf(pumpPrint, "After %.2f seconds - Feasibility pump exiting - took %.2f seconds",
2142                model_->getCurrentSeconds(), CoinCpuTime() - time1);
2143    else if (numberSolutions < maxSolutions)
2144      sprintf(pumpPrint, "After %.2f seconds - Feasibility pump exiting with objective of %g - took %.2f seconds",
2145              model_->getCurrentSeconds(), solutionValue, CoinCpuTime() - time1);
2146    else 
2147        sprintf(pumpPrint, "After %.2f seconds - Feasibility pump exiting with objective of %g (stopping after %d solutions) - took %.2f seconds",
2148                model_->getCurrentSeconds(), solutionValue, 
2149                numberSolutions,CoinCpuTime() - time1);
2150    model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
2151      << pumpPrint
2152      << CoinMessageEol;
2153    if (bestBasis.getNumStructural())
2154        model_->setBestSolutionBasis(bestBasis);
2155    //model_->setMinimizationObjValue(saveBestObjective);
2156    if ((accumulate_&1) != 0 && numberSolutions > 1 && !model_->getSolutionCount()) {
2157        model_->setSolutionCount(1); // for local search
2158        model_->setNumberHeuristicSolutions(1);
2159    }
2160#ifdef COIN_DEVELOP
2161    {
2162        double ncol = model_->solver()->getNumCols();
2163        double nrow = model_->solver()->getNumRows();
2164        printf("XXX total iterations %d ratios - %g %g %g\n",
2165               totalNumberIterations,
2166               static_cast<double> (totalNumberIterations) / nrow,
2167               static_cast<double> (totalNumberIterations) / ncol,
2168               static_cast<double> (totalNumberIterations) / (2*nrow + 2*ncol));
2169    }
2170#endif
2171    return finalReturnCode;
2172}
2173
2174/**************************END MAIN PROCEDURE ***********************************/
2175
2176// update model
2177void CbcHeuristicFPump::setModel(CbcModel * model)
2178{
2179    model_ = model;
2180}
2181
2182/* Rounds solution - down if < downValue
2183   returns 1 if current is a feasible solution
2184*/
2185int
2186CbcHeuristicFPump::rounds(OsiSolverInterface * solver, double * solution,
2187                          //const double * objective,
2188                          int numberIntegers, const int * integerVariable,
2189                          /*char * pumpPrint,*/ int iter,
2190                          /*bool roundExpensive,*/ double downValue, int *flip)
2191{
2192    double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
2193    double primalTolerance ;
2194    solver->getDblParam(OsiPrimalTolerance, primalTolerance) ;
2195
2196    int i;
2197
2198    const double * cost = solver->getObjCoefficients();
2199    int flip_up = 0;
2200    int flip_down  = 0;
2201    double  v = randomNumberGenerator_.randomDouble() * 20.0;
2202    int nn = 10 + static_cast<int> (v);
2203    int nnv = 0;
2204    int * list = new int [nn];
2205    double * val = new double [nn];
2206    for (i = 0; i < nn; i++) val[i] = .001;
2207
2208    const double * rowLower = solver->getRowLower();
2209    const double * rowUpper = solver->getRowUpper();
2210    int numberRows = solver->getNumRows();
2211    if (false && (iter&1) != 0) {
2212        // Do set covering variables
2213        const CoinPackedMatrix * matrixByRow = solver->getMatrixByRow();
2214        const double * elementByRow = matrixByRow->getElements();
2215        const int * column = matrixByRow->getIndices();
2216        const CoinBigIndex * rowStart = matrixByRow->getVectorStarts();
2217        const int * rowLength = matrixByRow->getVectorLengths();
2218        for (i = 0; i < numberRows; i++) {
2219            if (rowLower[i] == 1.0 && rowUpper[i] == 1.0) {
2220                bool cover = true;
2221                double largest = 0.0;
2222                int jColumn = -1;
2223                for (CoinBigIndex k = rowStart[i]; k < rowStart[i] + rowLength[i]; k++) {
2224                    int iColumn = column[k];
2225                    if (elementByRow[k] != 1.0 || !solver->isInteger(iColumn)) {
2226                        cover = false;
2227                        break;
2228                    } else {
2229                        if (solution[iColumn]) {
2230                            double value = solution[iColumn] *
2231                                           (randomNumberGenerator_.randomDouble() + 5.0);
2232                            if (value > largest) {
2233                                largest = value;
2234                                jColumn = iColumn;
2235                            }
2236                        }
2237                    }
2238                }
2239                if (cover) {
2240                    for (CoinBigIndex k = rowStart[i]; k < rowStart[i] + rowLength[i]; k++) {
2241                        int iColumn = column[k];
2242                        if (iColumn == jColumn)
2243                            solution[iColumn] = 1.0;
2244                        else
2245                            solution[iColumn] = 0.0;
2246                    }
2247                }
2248            }
2249        }
2250    }
2251    int numberColumns = solver->getNumCols();
2252#ifdef JJF_ZERO
2253    // Do set covering variables
2254    const CoinPackedMatrix * matrixByRow = solver->getMatrixByRow();
2255    const double * elementByRow = matrixByRow->getElements();
2256    const int * column = matrixByRow->getIndices();
2257    const CoinBigIndex * rowStart = matrixByRow->getVectorStarts();
2258    const int * rowLength = matrixByRow->getVectorLengths();
2259    double * sortTemp = new double[numberColumns];
2260    int * whichTemp = new int [numberColumns];
2261    char * rowTemp = new char [numberRows];
2262    memset(rowTemp, 0, numberRows);
2263    for (i = 0; i < numberColumns; i++)
2264        whichTemp[i] = -1;
2265    int nSOS = 0;
2266    for (i = 0; i < numberRows; i++) {
2267        if (rowLower[i] == 1.0 && rowUpper[i] == 1.0) {
2268            bool cover = true;
2269            for (CoinBigIndex k = rowStart[i]; k < rowStart[i] + rowLength[i]; k++) {
2270                int iColumn = column[k];
2271                if (elementByRow[k] != 1.0 || !solver->isInteger(iColumn)) {
2272                    cover = false;
2273                    break;
2274                }
2275            }
2276            if (cover) {
2277                rowTemp[i] = 1;
2278                nSOS++;
2279                for (CoinBigIndex k = rowStart[i]; k < rowStart[i] + rowLength[i]; k++) {
2280                    int iColumn = column[k];
2281                    double value = solution[iColumn];
2282                    whichTemp[iColumn] = iColumn;
2283                }
2284            }
2285        }
2286    }
2287    if (nSOS) {
2288        // Column copy
2289        const CoinPackedMatrix * matrixByColumn = solver->getMatrixByCol();
2290        //const double * element = matrixByColumn->getElements();
2291        const int * row = matrixByColumn->getIndices();
2292        const CoinBigIndex * columnStart = matrixByColumn->getVectorStarts();
2293        const int * columnLength = matrixByColumn->getVectorLengths();
2294        int nLook = 0;
2295        for (i = 0; i < numberColumns; i++) {
2296            if (whichTemp[i] >= 0) {
2297                whichTemp[nLook] = i;
2298                double value = solution[i];
2299                if (value < 0.5)
2300                    value *= (0.1 * randomNumberGenerator_.randomDouble() + 0.3);
2301                sortTemp[nLook++] = -value;
2302            }
2303        }
2304        CoinSort_2(sortTemp, sortTemp + nLook, whichTemp);
2305        double smallest = 1.0;
2306        int nFix = 0;
2307        int nOne = 0;
2308        for (int j = 0; j < nLook; j++) {
2309            int jColumn = whichTemp[j];
2310            double thisValue = solution[jColumn];
2311            if (!thisValue)
2312                continue;
2313            if (thisValue == 1.0)
2314                nOne++;
2315            smallest = CoinMin(smallest, thisValue);
2316            solution[jColumn] = 1.0;
2317            double largest = 0.0;
2318            for (CoinBigIndex jEl = columnStart[jColumn];
2319                    jEl < columnStart[jColumn] + columnLength[jColumn]; jEl++) {
2320                int jRow = row[jEl];
2321                if (rowTemp[jRow]) {
2322                    for (CoinBigIndex k = rowStart[jRow]; k < rowStart[jRow] + rowLength[jRow]; k++) {
2323                        int iColumn = column[k];
2324                        if (solution[iColumn]) {
2325                            if (iColumn != jColumn) {
2326                                double value = solution[iColumn];
2327                                if (value > largest)
2328                                    largest = value;
2329                                solution[iColumn] = 0.0;
2330                            }
2331                        }
2332                    }
2333                }
2334            }
2335            if (largest > thisValue)
2336                printf("%d was at %g - chosen over a value of %g\n",
2337                       jColumn, thisValue, largest);
2338            nFix++;
2339        }
2340        printf("%d fixed out of %d (%d at one already)\n",
2341               nFix, nLook, nOne);
2342    }
2343    delete [] sortTemp;
2344    delete [] whichTemp;
2345    delete [] rowTemp;
2346#endif
2347    const double * columnLower = solver->getColLower();
2348    const double * columnUpper = solver->getColUpper();
2349    // Check if valid with current solution (allow for 0.99999999s)
2350    for (i = 0; i < numberIntegers; i++) {
2351        int iColumn = integerVariable[i];
2352        double value = solution[iColumn];
2353        double round = floor(value + 0.5);
2354        if (fabs(value - round) > primalTolerance)
2355            break;
2356    }
2357    if (i == numberIntegers) {
2358        // may be able to use solution even if 0.99999's
2359        double * saveLower = CoinCopyOfArray(columnLower, numberColumns);
2360        double * saveUpper = CoinCopyOfArray(columnUpper, numberColumns);
2361        double * saveSolution = CoinCopyOfArray(solution, numberColumns);
2362        double * tempSolution = CoinCopyOfArray(solution, numberColumns);
2363        CoinWarmStartBasis  * saveBasis =
2364            dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ;
2365        for (i = 0; i < numberIntegers; i++) {
2366            int iColumn = integerVariable[i];
2367            double value = solution[iColumn];
2368            double round = floor(value + 0.5);
2369            solver->setColLower(iColumn, round);
2370            solver->setColUpper(iColumn, round);
2371            tempSolution[iColumn] = round;
2372        }
2373        solver->setColSolution(tempSolution);
2374        delete [] tempSolution;
2375        solver->resolve();
2376        solver->setColLower(saveLower);
2377        solver->setColUpper(saveUpper);
2378        solver->setWarmStart(saveBasis);
2379        delete [] saveLower;
2380        delete [] saveUpper;
2381        delete saveBasis;
2382        if (!solver->isProvenOptimal()) {
2383            solver->setColSolution(saveSolution);
2384        }
2385        delete [] saveSolution;
2386        if (solver->isProvenOptimal()) {
2387            // feasible
2388            delete [] list;
2389            delete [] val;
2390            return 1;
2391        }
2392    }
2393    //double * saveSolution = CoinCopyOfArray(solution,numberColumns);
2394    // return rounded solution
2395    for (i = 0; i < numberIntegers; i++) {
2396        int iColumn = integerVariable[i];
2397        double value = solution[iColumn];
2398        double round = floor(value + primalTolerance);
2399        if (value - round > downValue) round += 1.;
2400#ifndef JJF_ONE
2401        if (round < integerTolerance && cost[iColumn] < -1. + integerTolerance) flip_down++;
2402        if (round > 1. - integerTolerance && cost[iColumn] > 1. - integerTolerance) flip_up++;
2403#else
2404        if (round < columnLower[iColumn] + integerTolerance && cost[iColumn] < -1. + integerTolerance) flip_down++;
2405        if (round > columnUpper[iColumn] - integerTolerance && cost[iColumn] > 1. - integerTolerance) flip_up++;
2406#endif
2407        if (flip_up + flip_down == 0) {
2408            for (int k = 0; k < nn; k++) {
2409                if (fabs(value - round) > val[k]) {
2410                    nnv++;
2411                    for (int j = nn - 2; j >= k; j--) {
2412                        val[j+1] = val[j];
2413                        list[j+1] = list[j];
2414                    }
2415                    val[k] = fabs(value - round);
2416                    list[k] = iColumn;
2417                    break;
2418                }
2419            }
2420        }
2421        solution[iColumn] = round;
2422    }
2423
2424    if (nnv > nn) nnv = nn;
2425    //if (iter != 0)
2426    //sprintf(pumpPrint+strlen(pumpPrint),"up = %5d , down = %5d", flip_up, flip_down);
2427    *flip = flip_up + flip_down;
2428
2429    if (*flip == 0 && iter != 0) {
2430        //sprintf(pumpPrint+strlen(pumpPrint)," -- rand = %4d (%4d) ", nnv, nn);
2431        for (i = 0; i < nnv; i++) {
2432            // was solution[list[i]] = 1. - solution[list[i]]; but does that work for 7>=x>=6
2433            int index = list[i];
2434            double value = solution[index];
2435            if (value <= 1.0) {
2436                solution[index] = 1.0 - value;
2437            } else if (value < columnLower[index] + integerTolerance) {
2438                solution[index] = value + 1.0;
2439            } else if (value > columnUpper[index] - integerTolerance) {
2440                solution[index] = value - 1.0;
2441            } else {
2442                solution[index] = value - 1.0;
2443            }
2444        }
2445        *flip = nnv;
2446    } else {
2447        //sprintf(pumpPrint+strlen(pumpPrint)," ");
2448    }
2449    delete [] list;
2450    delete [] val;
2451    //iter++;
2452
2453    // get row activities
2454    double * rowActivity = new double[numberRows];
2455    memset(rowActivity, 0, numberRows*sizeof(double));
2456    solver->getMatrixByCol()->times(solution, rowActivity) ;
2457    double largestInfeasibility = primalTolerance;
2458    double sumInfeasibility = 0.0;
2459    int numberBadRows = 0;
2460    for (i = 0 ; i < numberRows ; i++) {
2461        double value;
2462        value = rowLower[i] - rowActivity[i];
2463        if (value > primalTolerance) {
2464            numberBadRows++;
2465            largestInfeasibility = CoinMax(largestInfeasibility, value);
2466            sumInfeasibility += value;
2467        }
2468        value = rowActivity[i] - rowUpper[i];
2469        if (value > primalTolerance) {
2470            numberBadRows++;
2471            largestInfeasibility = CoinMax(largestInfeasibility, value);
2472            sumInfeasibility += value;
2473        }
2474    }
2475#ifdef JJF_ZERO
2476    if (largestInfeasibility > primalTolerance && numberBadRows*10 < numberRows) {
2477        // Can we improve by flipping
2478        for (int iPass = 0; iPass < 10; iPass++) {
2479            int numberColumns = solver->getNumCols();
2480            const CoinPackedMatrix * matrixByCol = solver->getMatrixByCol();
2481            const double * element = matrixByCol->getElements();
2482            const int * row = matrixByCol->getIndices();
2483            const CoinBigIndex * columnStart = matrixByCol->getVectorStarts();
2484            const int * columnLength = matrixByCol->getVectorLengths();
2485            double oldSum = sumInfeasibility;
2486            // First improve by moving continuous ones
2487            for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
2488                if (!solver->isInteger(iColumn)) {
2489                    double solValue = solution[iColumn];
2490                    double thetaUp = columnUpper[iColumn] - solValue;
2491                    double improvementUp = 0.0;
2492                    if (thetaUp > primalTolerance) {
2493                        // can go up
2494                        for (CoinBigIndex j = columnStart[iColumn];
2495                                j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2496                            int iRow = row[j];
2497                            double distanceUp = rowUpper[iRow] - rowActivity[iRow];
2498                            double distanceDown = rowLower[iRow] - rowActivity[iRow];
2499                            double el = element[j];
2500                            if (el > 0.0) {
2501                                // positive element
2502                                if (distanceUp > 0.0) {
2503                                    if (thetaUp*el > distanceUp)
2504                                        thetaUp = distanceUp / el;
2505                                } else {
2506                                    improvementUp -= el;
2507                                }
2508                                if (distanceDown > 0.0) {
2509                                    if (thetaUp*el > distanceDown)
2510                                        thetaUp = distanceDown / el;
2511                                    improvementUp += el;
2512                                }
2513                            } else {
2514                                // negative element
2515                                if (distanceDown < 0.0) {
2516                                    if (thetaUp*el < distanceDown)
2517                                        thetaUp = distanceDown / el;
2518                                } else {
2519                                    improvementUp += el;
2520                                }
2521                                if (distanceUp < 0.0) {
2522                                    if (thetaUp*el < distanceUp)
2523                                        thetaUp = distanceUp / el;
2524                                    improvementUp -= el;
2525                                }
2526                            }
2527                        }
2528                    }
2529                    double thetaDown = solValue - columnLower[iColumn];
2530                    double improvementDown = 0.0;
2531                    if (thetaDown > primalTolerance) {
2532                        // can go down
2533                        for (CoinBigIndex j = columnStart[iColumn];
2534                                j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2535                            int iRow = row[j];
2536                            double distanceUp = rowUpper[iRow] - rowActivity[iRow];
2537                            double distanceDown = rowLower[iRow] - rowActivity[iRow];
2538                            double el = -element[j]; // not change in sign form up
2539                            if (el > 0.0) {
2540                                // positive element
2541                                if (distanceUp > 0.0) {
2542                                    if (thetaDown*el > distanceUp)
2543                                        thetaDown = distanceUp / el;
2544                                } else {
2545                                    improvementDown -= el;
2546                                }
2547                                if (distanceDown > 0.0) {
2548                                    if (thetaDown*el > distanceDown)
2549                                        thetaDown = distanceDown / el;
2550                                    improvementDown += el;
2551                                }
2552                            } else {
2553                                // negative element
2554                                if (distanceDown < 0.0) {
2555                                    if (thetaDown*el < distanceDown)
2556                                        thetaDown = distanceDown / el;
2557                                } else {
2558                                    improvementDown += el;
2559                                }
2560                                if (distanceUp < 0.0) {
2561                                    if (thetaDown*el < distanceUp)
2562                                        thetaDown = distanceUp / el;
2563                                    improvementDown -= el;
2564                                }
2565                            }
2566                        }
2567                        if (thetaUp < 1.0e-8)
2568                            improvementUp = 0.0;
2569                        if (thetaDown < 1.0e-8)
2570                            improvementDown = 0.0;
2571                        double theta;
2572                        if (improvementUp >= improvementDown) {
2573                            theta = thetaUp;
2574                        } else {
2575                            improvementUp = improvementDown;
2576                            theta = -thetaDown;
2577                        }
2578                        if (improvementUp > 1.0e-8 && fabs(theta) > 1.0e-8) {
2579                            // Could move
2580                            double oldSum = 0.0;
2581                            double newSum = 0.0;
2582                            solution[iColumn] += theta;
2583                            for (CoinBigIndex j = columnStart[iColumn];
2584                                    j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2585                                int iRow = row[j];
2586                                double lower = rowLower[iRow];
2587                                double upper = rowUpper[iRow];
2588                                double value = rowActivity[iRow];
2589                                if (value > upper)
2590                                    oldSum += value - upper;
2591                                else if (value < lower)
2592                                    oldSum += lower - value;
2593                                value += theta * element[j];
2594                                rowActivity[iRow] = value;
2595                                if (value > upper)
2596                                    newSum += value - upper;
2597                                else if (value < lower)
2598                                    newSum += lower - value;
2599                            }
2600                            assert (newSum <= oldSum);
2601                            sumInfeasibility += newSum - oldSum;
2602                        }
2603                    }
2604                }
2605            }
2606            // Now flip some integers?
2607#ifdef JJF_ZERO
2608            for (i = 0; i < numberIntegers; i++) {
2609                int iColumn = integerVariable[i];
2610                double solValue = solution[iColumn];
2611                assert (fabs(solValue - floor(solValue + 0.5)) < 1.0e-8);
2612                double improvementUp = 0.0;
2613                if (columnUpper[iColumn] >= solValue + 1.0) {
2614                    // can go up
2615                    double oldSum = 0.0;
2616                    double newSum = 0.0;
2617                    for (CoinBigIndex j = columnStart[iColumn];
2618                            j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2619                        int iRow = row[j];
2620                        double lower = rowLower[iRow];
2621                        double upper = rowUpper[iRow];
2622                        double value = rowActivity[iRow];
2623                        if (value > upper)
2624                            oldSum += value - upper;
2625                        else if (value < lower)
2626                            oldSum += lower - value;
2627                        value += element[j];
2628                        if (value > upper)
2629                            newSum += value - upper;
2630                        else if (value < lower)
2631                            newSum += lower - value;
2632                    }
2633                    improvementUp = oldSum - newSum;
2634                }
2635                double improvementDown = 0.0;
2636                if (columnLower[iColumn] <= solValue - 1.0) {
2637                    // can go down
2638                    double oldSum = 0.0;
2639                    double newSum = 0.0;
2640                    for (CoinBigIndex j = columnStart[iColumn];
2641                            j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2642                        int iRow = row[j];
2643                        double lower = rowLower[iRow];
2644                        double upper = rowUpper[iRow];
2645                        double value = rowActivity[iRow];
2646                        if (value > upper)
2647                            oldSum += value - upper;
2648                        else if (value < lower)
2649                            oldSum += lower - value;
2650                        value -= element[j];
2651                        if (value > upper)
2652                            newSum += value - upper;
2653                        else if (value < lower)
2654                            newSum += lower - value;
2655                    }
2656                    improvementDown = oldSum - newSum;
2657                }
2658                double theta;
2659                if (improvementUp >= improvementDown) {
2660                    theta = 1.0;
2661                } else {
2662                    improvementUp = improvementDown;
2663                    theta = -1.0;
2664                }
2665                if (improvementUp > 1.0e-8 && fabs(theta) > 1.0e-8) {
2666                    // Could move
2667                    double oldSum = 0.0;
2668                    double newSum = 0.0;
2669                    solution[iColumn] += theta;
2670                    for (CoinBigIndex j = columnStart[iColumn];
2671                            j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2672                        int iRow = row[j];
2673                        double lower = rowLower[iRow];
2674                        double upper = rowUpper[iRow];
2675                        double value = rowActivity[iRow];
2676                        if (value > upper)
2677                            oldSum += value - upper;
2678                        else if (value < lower)
2679                            oldSum += lower - value;
2680                        value += theta * element[j];
2681                        rowActivity[iRow] = value;
2682                        if (value > upper)
2683                            newSum += value - upper;
2684                        else if (value < lower)
2685                            newSum += lower - value;
2686                    }
2687                    assert (newSum <= oldSum);
2688                    sumInfeasibility += newSum - oldSum;
2689                }
2690            }
2691#else
2692            int bestColumn = -1;
2693            double bestImprovement = primalTolerance;
2694            double theta = 0.0;
2695            for (i = 0; i < numberIntegers; i++) {
2696                int iColumn = integerVariable[i];
2697                double solValue = solution[iColumn];
2698                assert (fabs(solValue - floor(solValue + 0.5)) < 1.0e-8);
2699                double improvementUp = 0.0;
2700                if (columnUpper[iColumn] >= solValue + 1.0) {
2701                    // can go up
2702                    double oldSum = 0.0;
2703                    double newSum = 0.0;
2704                    for (CoinBigIndex j = columnStart[iColumn];
2705                            j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2706                        int iRow = row[j];
2707                        double lower = rowLower[iRow];
2708                        double upper = rowUpper[iRow];
2709                        double value = rowActivity[iRow];
2710                        if (value > upper)
2711                            oldSum += value - upper;
2712                        else if (value < lower)
2713                            oldSum += lower - value;
2714                        value += element[j];
2715                        if (value > upper)
2716                            newSum += value - upper;
2717                        else if (value < lower)
2718                            newSum += lower - value;
2719                    }
2720                    improvementUp = oldSum - newSum;
2721                }
2722                double improvementDown = 0.0;
2723                if (columnLower[iColumn] <= solValue - 1.0) {
2724                    // can go down
2725                    double oldSum = 0.0;
2726                    double newSum = 0.0;
2727                    for (CoinBigIndex j = columnStart[iColumn];
2728                            j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2729                        int iRow = row[j];
2730                        double lower = rowLower[iRow];
2731                        double upper = rowUpper[iRow];
2732                        double value = rowActivity[iRow];
2733                        if (value > upper)
2734                            oldSum += value - upper;
2735                        else if (value < lower)
2736                            oldSum += lower - value;
2737                        value -= element[j];
2738                        if (value > upper)
2739                            newSum += value - upper;
2740                        else if (value < lower)
2741                            newSum += lower - value;
2742                    }
2743                    improvementDown = oldSum - newSum;
2744                }
2745                double improvement = CoinMax(improvementUp, improvementDown);
2746                if (improvement > bestImprovement) {
2747                    bestImprovement = improvement;
2748                    bestColumn = iColumn;
2749                    if (improvementUp > improvementDown)
2750                        theta = 1.0;
2751                    else
2752                        theta = -1.0;
2753                }
2754            }
2755            if (bestColumn >= 0) {
2756                // Could move
2757                int iColumn = bestColumn;
2758                double oldSum = 0.0;
2759                double newSum = 0.0;
2760                solution[iColumn] += theta;
2761                for (CoinBigIndex j = columnStart[iColumn];
2762                        j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2763                    int iRow = row[j];
2764                    double lower = rowLower[iRow];
2765                    double upper = rowUpper[iRow];
2766                    double value = rowActivity[iRow];
2767                    if (value > upper)
2768                        oldSum += value - upper;
2769                    else if (value < lower)
2770                        oldSum += lower - value;
2771                    value += theta * element[j];
2772                    rowActivity[iRow] = value;
2773                    if (value > upper)
2774                        newSum += value - upper;
2775                    else if (value < lower)
2776                        newSum += lower - value;
2777                }
2778                assert (newSum <= oldSum);
2779                sumInfeasibility += newSum - oldSum;
2780            }
2781#endif
2782            if (oldSum <= sumInfeasibility + primalTolerance)
2783                break; // no good
2784        }
2785    }
2786    //delete [] saveSolution;
2787#endif
2788    delete [] rowActivity;
2789    return (largestInfeasibility > primalTolerance) ? 0 : 1;
2790}
2791// Set maximum Time (default off) - also sets starttime to current
2792void
2793CbcHeuristicFPump::setMaximumTime(double value)
2794{
2795    startTime_ = CoinCpuTime();
2796    maximumTime_ = value;
2797}
2798
2799# ifdef COIN_HAS_CLP
2800
2801//#############################################################################
2802// Constructors / Destructor / Assignment
2803//#############################################################################
2804
2805//-------------------------------------------------------------------
2806// Default Constructor
2807//-------------------------------------------------------------------
2808CbcDisasterHandler::CbcDisasterHandler (CbcModel * model)
2809        : OsiClpDisasterHandler(),
2810        cbcModel_(model)
2811{
2812    if (model) {
2813        osiModel_
2814        = dynamic_cast<OsiClpSolverInterface *> (model->solver());
2815        if (osiModel_)
2816            setSimplex(osiModel_->getModelPtr());
2817    }
2818}
2819
2820//-------------------------------------------------------------------
2821// Copy constructor
2822//-------------------------------------------------------------------
2823CbcDisasterHandler::CbcDisasterHandler (const CbcDisasterHandler & rhs)
2824        : OsiClpDisasterHandler(rhs),
2825        cbcModel_(rhs.cbcModel_)
2826{
2827}
2828
2829
2830//-------------------------------------------------------------------
2831// Destructor
2832//-------------------------------------------------------------------
2833CbcDisasterHandler::~CbcDisasterHandler ()
2834{
2835}
2836
2837//----------------------------------------------------------------
2838// Assignment operator
2839//-------------------------------------------------------------------
2840CbcDisasterHandler &
2841CbcDisasterHandler::operator=(const CbcDisasterHandler & rhs)
2842{
2843    if (this != &rhs) {
2844        OsiClpDisasterHandler::operator=(rhs);
2845        cbcModel_ = rhs.cbcModel_;
2846    }
2847    return *this;
2848}
2849//-------------------------------------------------------------------
2850// Clone
2851//-------------------------------------------------------------------
2852ClpDisasterHandler * CbcDisasterHandler::clone() const
2853{
2854    return new CbcDisasterHandler(*this);
2855}
2856// Type of disaster 0 can fix, 1 abort
2857int
2858CbcDisasterHandler::typeOfDisaster()
2859{
2860    if (!cbcModel_->parentModel() &&
2861            (cbcModel_->specialOptions()&2048) == 0) {
2862        return 0;
2863    } else {
2864        if (cbcModel_->parentModel())
2865            cbcModel_->setMaximumNodes(0);
2866        return 1;
2867    }
2868}
2869/* set model. */
2870void
2871CbcDisasterHandler::setCbcModel(CbcModel * model)
2872{
2873    cbcModel_ = model;
2874    if (model) {
2875        osiModel_
2876        = dynamic_cast<OsiClpSolverInterface *> (model->solver());
2877        if (osiModel_)
2878            setSimplex(osiModel_->getModelPtr());
2879        else
2880            setSimplex(NULL);
2881    }
2882}
2883#endif
2884
Note: See TracBrowser for help on using the repository browser.