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

Last change on this file since 1943 was 1943, checked in by forrest, 6 years ago

more options, copy statistics structure analysis
start coding of "switch" variables i.e. badly scaled ints or hi/lo
changes to allow more influence on small branch and bound
changes to get correct printout with threads

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