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

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

new event plus flexible output format

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 138.8 KB
Line 
1/* $Id: CbcHeuristicFPump.cpp 1954 2013-08-06 19:19:32Z 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                          bool takeHint;
836                          OsiHintStrength strength;
837                          solver->getHintParam(OsiDoDualInResolve, takeHint, strength);
838                          //solver->setHintParam(OsiDoReducePrint, false, OsiHintTry);
839                          solver->setHintParam(OsiDoDualInResolve, false, OsiHintDo);
840                          //solver->setHintParam(OsiDoScale, false, OsiHintDo);
841                          solver->resolve();
842                          solver->setHintParam(OsiDoDualInResolve, takeHint, strength);
843                          if (solver->isProvenOptimal()) {
844                            memcpy(newSolution,solver->getColSolution(),
845                                   numberColumns*sizeof(double));
846                            newSolutionValue = -saveOffset;
847                            for (  i = 0 ; i < numberColumns ; i++ ) {
848                              newSolutionValue += saveObjective[i] * newSolution[i];
849                            }
850                            newSolutionValue *= direction;
851                            sprintf(pumpPrint, "Relaxing continuous gives %g", newSolutionValue);
852                            //#define DEBUG_BEST
853#ifdef DEBUG_BEST
854                            {
855                              int numberColumns=solver->getNumCols();
856                              FILE * fp = fopen("solution.data2","wb");
857                              printf("Solution data on file solution.data2\n");
858                              size_t numberWritten;
859                              numberWritten=fwrite(&numberColumns,sizeof(int),1,fp);
860                              assert (numberWritten==1);
861                              numberWritten=fwrite(&newSolutionValue,sizeof(double),1,fp);
862                              assert (numberWritten==1);
863                              numberWritten=fwrite(newSolution,sizeof(double),numberColumns,fp);
864                              assert (numberWritten==numberColumns);
865                              fclose(fp);
866                              const double * rowLower = solver->getRowLower();
867                              const double * rowUpper = solver->getRowUpper();
868                              const double * columnLower = solver->getColLower();
869                              const double * columnUpper = solver->getColUpper();
870                              int numberRows = solver->getNumRows() ;
871                              double *rowActivity = new double[numberRows] ;
872                              memset(rowActivity, 0, numberRows*sizeof(double)) ;
873                              const double * element = solver->getMatrixByCol()->getElements();
874                              const int * row = solver->getMatrixByCol()->getIndices();
875                              const CoinBigIndex * columnStart = solver->getMatrixByCol()->getVectorStarts();
876                              const int * columnLength = solver->getMatrixByCol()->getVectorLengths();
877                              double largestAway=0.0;
878                              int away=-1;
879                              double saveOffset;
880                              solver->getDblParam(OsiObjOffset, saveOffset);
881                              double newSolutionValue = -saveOffset;
882                              const double * objective = solver->getObjCoefficients();
883                              for ( int iColumn=0 ; iColumn<numberColumns ; ++iColumn ) {
884                                double value=newSolution[iColumn];
885                                CoinBigIndex start = columnStart[iColumn];
886                                CoinBigIndex end = start + columnLength[iColumn];
887                                for (CoinBigIndex j = start; j < end; j++) {
888                                  int iRow = row[j];
889                                  if (iRow==1996)
890                                    printf("fp col %d val %g el %g old y %g\n",
891                                  iColumn,value,element[j],rowActivity[iRow]);
892                                  rowActivity[iRow] += value * element[j];
893                                }
894                                newSolutionValue += objective[iColumn] * newSolution[iColumn];
895                                if (solver->isInteger(iColumn)) {
896                                  double intValue = floor(value+0.5);
897                                  if (fabs(value-intValue)>largestAway) {
898                                    largestAway=fabs(value-intValue);
899                                    away=iColumn;
900                                  }
901                                }
902                              }
903                              printf("Largest away from int at column %d was %g - obj %g\n",away,
904                                     largestAway,newSolutionValue);
905                              double largestInfeasibility=0.0;
906                              for (int i = 0 ; i < numberRows ; i++) {
907#if 0 //def CLP_INVESTIGATE
908                                double inf;
909                                inf = rowLower[i] - rowActivity[i];
910                                if (inf > primalTolerance)
911                                  printf("Row %d inf %g sum %g %g <= %g <= %g\n",
912                                         i, inf, rowSum[i], rowLower[i], rowActivity[i], rowUpper[i]);
913                                inf = rowActivity[i] - rowUpper[i];
914                                if (inf > primalTolerance)
915                                  printf("Row %d inf %g %g <= %g <= %g\n",
916                                         i, inf, rowLower[i], rowActivity[i], rowUpper[i]);
917#endif
918                                double infeasibility = CoinMax(rowActivity[i]-rowUpper[i],
919                                                               rowLower[i]-rowActivity[i]);
920                                if (infeasibility>largestInfeasibility) {
921                                  largestInfeasibility = infeasibility;
922                                  printf("Binf of %g on row %d\n",
923                                         infeasibility,i);
924                                }
925                              }
926                              delete [] rowActivity ;
927                              printf("Blargest infeasibility is %g - obj %g\n", largestInfeasibility,newSolutionValue);
928                            }
929#endif
930                          } else {
931                            sprintf(pumpPrint,"Infeasible when relaxing continuous!\n");
932                          }
933                          model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
934                            << pumpPrint
935                            << CoinMessageEol;
936                        }
937                    }
938                    if (returnCode && newSolutionValue < saveValue) {
939                        memcpy(betterSolution, newSolution, numberColumns*sizeof(double));
940                        solutionFound = true;
941                        if (exitNow(newSolutionValue))
942                            exitAll = true;
943                        CoinWarmStartBasis * basis =
944                            dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ;
945                        if (basis) {
946                            bestBasis = * basis;
947                            delete basis;
948                            int action = model_->dealWithEventHandler(CbcEventHandler::heuristicSolution, newSolutionValue, betterSolution);
949                            if (action == 0) {
950                                double * saveOldSolution = CoinCopyOfArray(model_->bestSolution(), numberColumns);
951                                double saveObjectiveValue = model_->getMinimizationObjValue();
952                                model_->setBestSolution(betterSolution, numberColumns, newSolutionValue);
953                                if (saveOldSolution && saveObjectiveValue < model_->getMinimizationObjValue())
954                                    model_->setBestSolution(saveOldSolution, numberColumns, saveObjectiveValue);
955                                delete [] saveOldSolution;
956                            }
957                            if (action == 0 || model_->maximumSecondsReached()) {
958                                exitAll = true; // exit
959                                break;
960                            }
961                        }
962                        if ((accumulate_&1) != 0) {
963                            model_->incrementUsed(betterSolution); // for local search
964                        }
965                        solutionValue = newSolutionValue;
966                        solutionFound = true;
967                        numberSolutions++;
968                        if (numberSolutions>=maxSolutions)
969                          exitAll = true;
970                        if (general && saveValue != newSolutionValue) {
971                            sprintf(pumpPrint, "Cleaned solution of %g", solutionValue);
972                            model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
973                            << pumpPrint
974                            << CoinMessageEol;
975                        }
976                        if (exitNow(newSolutionValue))
977                            exitAll = true;
978                    } else {
979                        sprintf(pumpPrint, "Mini branch and bound could not fix general integers");
980                        model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
981                        << pumpPrint
982                        << CoinMessageEol;
983                    }
984                } else {
985                    sprintf(pumpPrint, "After further testing solution no better than previous of %g", solutionValue);
986                    model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
987                    << pumpPrint
988                    << CoinMessageEol;
989                    //newLineNeeded=false;
990                    returnCode = 0;
991                }
992                break;
993            } else {
994                // SOLUTION IS not INTEGER
995                // 1. check for loop
996                bool matched;
997                for (int k = NUMBER_OLD - 1; k > 0; k--) {
998                    double * b = oldSolution[k];
999                    matched = true;
1000                    for (i = 0; i < numberIntegers; i++) {
1001                        int iColumn = integerVariable[i];
1002                        if (newSolution[iColumn] != b[iColumn]) {
1003                            matched = false;
1004                            break;
1005                        }
1006                    }
1007                    if (matched) break;
1008                }
1009                int numberPerturbed = 0;
1010                if (matched || numberPasses % 100 == 0) {
1011                    // perturbation
1012                    //sprintf(pumpPrint+strlen(pumpPrint)," perturbation applied");
1013                    //newLineNeeded=true;
1014                    double factorX[10] = {0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0};
1015                    double factor = 1.0;
1016                    double target = -1.0;
1017                    double * randomX = new double [numberIntegers];
1018                    for (i = 0; i < numberIntegers; i++)
1019                        randomX[i] = CoinMax(0.0, randomNumberGenerator_.randomDouble() - 0.3);
1020                    for (int k = 0; k < 10; k++) {
1021#ifdef COIN_DEVELOP_x
1022                        printf("kpass %d\n", k);
1023#endif
1024                        int numberX[10] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
1025                        for (i = 0; i < numberIntegers; i++) {
1026                            int iColumn = integerVariable[i];
1027                            double value = randomX[i];
1028                            double difference = fabs(solution[iColumn] - newSolution[iColumn]);
1029                            for (int j = 0; j < 10; j++) {
1030                                if (difference + value*factorX[j] > 0.5)
1031                                    numberX[j]++;
1032                            }
1033                        }
1034                        if (target < 0.0) {
1035                            if (numberX[9] <= 200)
1036                                break; // not very many changes
1037                            target = CoinMax(200.0, CoinMin(0.05 * numberX[9], 1000.0));
1038                        }
1039                        int iX = -1;
1040                        int iBand = -1;
1041                        for (i = 0; i < 10; i++) {
1042#ifdef COIN_DEVELOP_x
1043                            printf("** %d changed at %g\n", numberX[i], factorX[i]);
1044#endif
1045                            if (numberX[i] >= target && numberX[i] < 2.0*target && iX < 0)
1046                                iX = i;
1047                            if (iBand<0 && numberX[i]>target) {
1048                                iBand = i;
1049                                factor = factorX[i];
1050                            }
1051                        }
1052                        if (iX >= 0) {
1053                            factor = factorX[iX];
1054                            break;
1055                        } else {
1056                            assert (iBand >= 0);
1057                            double hi = factor;
1058                            double lo = (iBand > 0) ? factorX[iBand-1] : 0.0;
1059                            double diff = (hi - lo) / 9.0;
1060                            for (i = 0; i < 10; i++) {
1061                                factorX[i] = lo;
1062                                lo += diff;
1063                            }
1064                        }
1065                    }
1066                    for (i = 0; i < numberIntegers; i++) {
1067                        int iColumn = integerVariable[i];
1068                        double value = randomX[i];
1069                        double difference = fabs(solution[iColumn] - newSolution[iColumn]);
1070                        if (difference + value*factor > 0.5) {
1071                            numberPerturbed++;
1072                            if (newSolution[iColumn] < lower[iColumn] + primalTolerance) {
1073                                newSolution[iColumn] += 1.0;
1074                            } else if (newSolution[iColumn] > upper[iColumn] - primalTolerance) {
1075                                newSolution[iColumn] -= 1.0;
1076                            } else {
1077                                // general integer
1078                                if (difference + value > 0.75)
1079                                    newSolution[iColumn] += 1.0;
1080                                else
1081                                    newSolution[iColumn] -= 1.0;
1082                            }
1083                        }
1084                    }
1085                    delete [] randomX;
1086                } else {
1087                    for (j = NUMBER_OLD - 1; j > 0; j--) {
1088                        for (i = 0; i < numberColumns; i++) oldSolution[j][i] = oldSolution[j-1][i];
1089                    }
1090                    for (j = 0; j < numberColumns; j++) oldSolution[0][j] = newSolution[j];
1091                }
1092
1093                // 2. update the objective function based on the new rounded solution
1094                double offset = 0.0;
1095                double costValue = (1.0 - scaleFactor) * solver->getObjSense();
1096                int numberChanged = 0;
1097                const double * oldObjective = solver->getObjCoefficients();
1098                for (i = 0; i < numberColumns; i++) {
1099                    // below so we can keep original code and allow for objective
1100                    int iColumn = i;
1101                    // Special code for "artificials"
1102                    if (direction*saveObjective[iColumn] >= artificialCost_) {
1103                        //solver->setObjCoeff(iColumn,scaleFactor*saveObjective[iColumn]);
1104                        solver->setObjCoeff(iColumn, (artificialFactor*saveObjective[iColumn]) / artificialCost_);
1105                    }
1106                    if (!solver->isBinary(iColumn) && !doGeneral)
1107                        continue;
1108                    // deal with fixed variables (i.e., upper=lower)
1109                    if (fabs(lower[iColumn] - upper[iColumn]) < primalTolerance || !solver->isInteger(iColumn)) {
1110                        //if (lower[iColumn] > 1. - primalTolerance) solver->setObjCoeff(iColumn,-costValue);
1111                        //else                                       solver->setObjCoeff(iColumn,costValue);
1112                        continue;
1113                    }
1114                    double newValue = 0.0;
1115                    if (newSolution[iColumn] < lower[iColumn] + primalTolerance) {
1116                      newValue = costValue + scaleFactor * saveObjective[iColumn];
1117                    } else {
1118                      if (newSolution[iColumn] > upper[iColumn] - primalTolerance) {
1119                        newValue = -costValue + scaleFactor * saveObjective[iColumn];
1120                      }
1121                    }
1122#ifdef RAND_RAND
1123                    if (!offRandom)
1124                        newValue *= randomFactor[iColumn];
1125#endif
1126                    if (newValue != oldObjective[iColumn]) {
1127                        numberChanged++;
1128                    }
1129                    solver->setObjCoeff(iColumn, newValue);
1130                    offset += costValue * newSolution[iColumn];
1131                }
1132                if (numberPasses==1 && !totalNumberPasses && (model_->specialOptions()&8388608)!=0) {
1133                  // doing multiple solvers - make a real difference - flip 5%
1134                  for (i = 0; i < numberIntegers; i++) {
1135                    int iColumn = integerVariable[i];
1136                    double value = floor(newSolution[iColumn]+0.5);
1137                    if (fabs(value-solution[iColumn])>primalTolerance) {
1138                      value = randomNumberGenerator_.randomDouble();
1139                      if(value<0.05) {
1140                        //printf("Flipping %d - random %g\n",iColumn,value);
1141                        solver->setObjCoeff(iColumn,-solver->getObjCoefficients()[iColumn]);
1142                      }
1143                    }
1144                  }
1145                }
1146                solver->setDblParam(OsiObjOffset, -offset);
1147                if (!general && false) {
1148                    // Solve in two goes - first keep satisfied ones fixed
1149                    double * saveLower = new double [numberIntegers];
1150                    double * saveUpper = new double [numberIntegers];
1151                    for (i = 0; i < numberIntegers; i++) {
1152                        int iColumn = integerVariable[i];
1153                        saveLower[i] = COIN_DBL_MAX;
1154                        saveUpper[i] = -COIN_DBL_MAX;
1155                        if (solution[iColumn] < lower[iColumn] + primalTolerance) {
1156                            saveUpper[i] = upper[iColumn];
1157                            solver->setColUpper(iColumn, lower[iColumn]);
1158                        } else if (solution[iColumn] > upper[iColumn] - primalTolerance) {
1159                            saveLower[i] = lower[iColumn];
1160                            solver->setColLower(iColumn, upper[iColumn]);
1161                        }
1162                    }
1163                    solver->resolve();
1164                    if (!solver->isProvenOptimal()) {
1165                        // presumably max time or some such
1166                        exitAll = true;
1167                        break;
1168                    }
1169                    for (i = 0; i < numberIntegers; i++) {
1170                        int iColumn = integerVariable[i];
1171                        if (saveLower[i] != COIN_DBL_MAX)
1172                            solver->setColLower(iColumn, saveLower[i]);
1173                        if (saveUpper[i] != -COIN_DBL_MAX)
1174                            solver->setColUpper(iColumn, saveUpper[i]);
1175                        saveUpper[i] = -COIN_DBL_MAX;
1176                    }
1177                    memcpy(newSolution, solution, numberColumns*sizeof(double));
1178                    int flip;
1179                    returnCode = rounds(solver, newSolution,/*saveObjective,*/
1180                                        numberIntegers, integerVariable,
1181                                        /*pumpPrint,*/numberPasses,
1182                                        /*roundExpensive_,*/defaultRounding_, &flip);
1183                    numberPasses++;
1184                    if (returnCode) {
1185                        // solution - but may not be better
1186                        // Compute using dot product
1187                        double newSolutionValue = -saveOffset;
1188                        for (  i = 0 ; i < numberColumns ; i++ )
1189                            newSolutionValue += saveObjective[i] * newSolution[i];
1190                        newSolutionValue *= direction;
1191                        sprintf(pumpPrint, "Intermediate solution found of %g", newSolutionValue);
1192                        model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
1193                        << pumpPrint
1194                        << CoinMessageEol;
1195                        if (newSolutionValue < solutionValue) {
1196                            memcpy(betterSolution, newSolution, numberColumns*sizeof(double));
1197                            CoinWarmStartBasis * basis =
1198                                dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ;
1199                            solutionFound = true;
1200                            numberSolutions++;
1201                            if (numberSolutions>=maxSolutions)
1202                              exitAll = true;
1203                            if (exitNow(newSolutionValue))
1204                                exitAll = true;
1205                            if (basis) {
1206                                bestBasis = * basis;
1207                                delete basis;
1208                                int action = model_->dealWithEventHandler(CbcEventHandler::heuristicSolution, newSolutionValue, betterSolution);
1209                                if (!action) {
1210                                    double * saveOldSolution = CoinCopyOfArray(model_->bestSolution(), numberColumns);
1211                                    double saveObjectiveValue = model_->getMinimizationObjValue();
1212                                    model_->setBestSolution(betterSolution, numberColumns, newSolutionValue);
1213                                    if (saveOldSolution && saveObjectiveValue < model_->getMinimizationObjValue())
1214                                        model_->setBestSolution(saveOldSolution, numberColumns, saveObjectiveValue);
1215                                    delete [] saveOldSolution;
1216                                }
1217                                if (!action || model_->maximumSecondsReached()) {
1218                                    exitAll = true; // exit
1219                                    break;
1220                                }
1221                            }
1222                            if ((accumulate_&1) != 0) {
1223                                model_->incrementUsed(betterSolution); // for local search
1224                            }
1225                            solutionValue = newSolutionValue;
1226                            solutionFound = true;
1227                            numberSolutions++;
1228                            if (numberSolutions>=maxSolutions)
1229                              exitAll = true;
1230                            if (exitNow(newSolutionValue))
1231                                exitAll = true;
1232                        } else {
1233                            returnCode = 0;
1234                        }
1235                    }
1236                }
1237                int numberIterations = 0;
1238                if (!doGeneral) {
1239                    // faster to do from all slack!!!!
1240                    if (allSlack) {
1241                        CoinWarmStartBasis dummy;
1242                        solver->setWarmStart(&dummy);
1243                    }
1244#ifdef COIN_DEVELOP
1245                    printf("%d perturbed out of %d columns (%d changed)\n", numberPerturbed, numberColumns, numberChanged);
1246#endif
1247                    bool takeHint;
1248                    OsiHintStrength strength;
1249                    solver->getHintParam(OsiDoDualInResolve, takeHint, strength);
1250                    if (dualPass && numberChanged > 2) {
1251                        solver->setHintParam(OsiDoDualInResolve, true); // dual may be better
1252                        if (dualPass == 1 && 2*numberChanged < numberColumns &&
1253                                (numberChanged < 5000 || 6*numberChanged < numberColumns)) {
1254                            // but we need to make infeasible
1255                            CoinWarmStartBasis * basis =
1256                                dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ;
1257                            if (basis) {
1258                                // modify
1259                                const double * lower = solver->getColLower();
1260                                const double * upper = solver->getColUpper();
1261                                double * solution = CoinCopyOfArray(solver->getColSolution(),
1262                                                                    numberColumns);
1263                                const double * objective = solver->getObjCoefficients();
1264                                int nChanged = 0;
1265                                for (i = 0; i < numberIntegersOrig; i++) {
1266                                    int iColumn = integerVariableOrig[i];
1267#ifdef RAND_RAND
1268                                    if (nChanged > numberChanged)
1269                                        break;
1270#endif
1271                                    if (objective[iColumn] > 0.0) {
1272                                        if (basis->getStructStatus(iColumn) ==
1273                                                CoinWarmStartBasis::atUpperBound) {
1274                                            solution[iColumn] = lower[iColumn];
1275                                            basis->setStructStatus(iColumn, CoinWarmStartBasis::atLowerBound);
1276                                            nChanged++;
1277                                        }
1278                                    } else if (objective[iColumn] < 0.0) {
1279                                        if (basis->getStructStatus(iColumn) ==
1280                                                CoinWarmStartBasis::atLowerBound) {
1281                                            solution[iColumn] = upper[iColumn];
1282                                            basis->setStructStatus(iColumn, CoinWarmStartBasis::atUpperBound);
1283                                            nChanged++;
1284                                        }
1285                                    }
1286                                }
1287                                if (!nChanged) {
1288                                    for (i = 0; i < numberIntegersOrig; i++) {
1289                                        int iColumn = integerVariableOrig[i];
1290                                        if (objective[iColumn] > 0.0) {
1291                                            if (basis->getStructStatus(iColumn) ==
1292                                                    CoinWarmStartBasis::basic) {
1293                                                solution[iColumn] = lower[iColumn];
1294                                                basis->setStructStatus(iColumn, CoinWarmStartBasis::atLowerBound);
1295                                                break;
1296                                            }
1297                                        } else if (objective[iColumn] < 0.0) {
1298                                            if (basis->getStructStatus(iColumn) ==
1299                                                    CoinWarmStartBasis::basic) {
1300                                                solution[iColumn] = upper[iColumn];
1301                                                basis->setStructStatus(iColumn, CoinWarmStartBasis::atUpperBound);
1302                                                break;
1303                                            }
1304                                        }
1305                                    }
1306                                }
1307                                solver->setColSolution(solution);
1308                                delete [] solution;
1309                                solver->setWarmStart(basis);
1310                                delete basis;
1311                            }
1312                        } else {
1313                            // faster to do from all slack!!!! ???
1314                            CoinWarmStartBasis dummy;
1315                            solver->setWarmStart(&dummy);
1316                        }
1317                    }
1318                    if (numberTries > 1 && numberPasses == 1 && firstPerturbedObjective) {
1319                        // Modify to use convex combination
1320                        // use basis from first time
1321                        solver->setWarmStart(&saveBasis);
1322                        // and objective
1323                        if (secondPassOpt < 3 || (secondPassOpt >= 4 && secondPassOpt < 6))
1324                            solver->setObjective(firstPerturbedObjective);
1325                        // and solution
1326                        solver->setColSolution(firstPerturbedSolution);
1327                        //if (secondPassOpt==2||secondPassOpt==5||
1328                        if (firstPerturbedValue > cutoff)
1329                            solver->setHintParam(OsiDoDualInResolve, true); // dual may be better
1330                    }
1331                    solver->resolve();
1332                    if (!solver->isProvenOptimal()) {
1333                        // presumably max time or some such
1334                        exitAll = true;
1335                        break;
1336                    }
1337                    solver->setHintParam(OsiDoDualInResolve, takeHint);
1338                    newTrueSolutionValue = -saveOffset;
1339                    newSumInfeas = 0.0;
1340                    newNumberInfeas = 0;
1341                    {
1342                        const double * newSolution = solver->getColSolution();
1343                        for (int iColumn = 0 ; iColumn < numberColumns ; iColumn++ ) {
1344                            if (solver->isInteger(iColumn)) {
1345                                double value = newSolution[iColumn];
1346                                double nearest = floor(value + 0.5);
1347                                newSumInfeas += fabs(value - nearest);
1348                                if (fabs(value - nearest) > 1.0e-6) {
1349                                    newNumberInfeas++;
1350                                }
1351                            }
1352                            newTrueSolutionValue += saveObjective[iColumn] * newSolution[iColumn];
1353                        }
1354                        newTrueSolutionValue *= direction;
1355                        if (numberPasses == 1 && secondPassOpt) {
1356                            if (numberTries == 1 || secondPassOpt > 3) {
1357                                // save basis
1358                                CoinWarmStartBasis * basis =
1359                                    dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ;
1360                                if (basis) {
1361                                    saveBasis = * basis;
1362                                    delete basis;
1363                                }
1364                                delete [] firstPerturbedObjective;
1365                                delete [] firstPerturbedSolution;
1366                                firstPerturbedObjective = CoinCopyOfArray(solver->getObjCoefficients(), numberColumns);
1367                                firstPerturbedSolution = CoinCopyOfArray(solver->getColSolution(), numberColumns);
1368                                firstPerturbedValue = newTrueSolutionValue;
1369                            }
1370                        }
1371                        if (newNumberInfeas && newNumberInfeas < 15) {
1372#ifdef JJF_ZERO
1373                            roundingObjective = solutionValue;
1374                            OsiSolverInterface * saveSolver = model_->swapSolver(solver);
1375                            double * currentObjective =
1376                                CoinCopyOfArray(solver->getObjCoefficients(), numberColumns);
1377                            solver->setObjective(saveObjective);
1378                            double saveOffset2;
1379                            solver->getDblParam(OsiObjOffset, saveOffset2);
1380                            solver->setDblParam(OsiObjOffset, saveOffset);
1381                            int ifSol = roundingHeuristic.solution(roundingObjective, roundingSolution);
1382                            solver->setObjective(currentObjective);
1383                            solver->setDblParam(OsiObjOffset, saveOffset2);
1384                            delete [] currentObjective;
1385                            model_->swapSolver(saveSolver);
1386                            if (ifSol > 0)
1387                                abort();
1388#endif
1389                            int numberRows = solver->getNumRows();
1390                            double * rowActivity = new double[numberRows];
1391                            memset(rowActivity, 0, numberRows*sizeof(double));
1392                            int * which = new int[newNumberInfeas];
1393                            int * stack = new int[newNumberInfeas+1];
1394                            double * baseValue = new double[newNumberInfeas];
1395                            int * whichRow = new int[numberRows];
1396                            double * rowValue = new double[numberRows];
1397                            memset(rowValue, 0, numberRows*sizeof(double));
1398                            int nRow = 0;
1399                            // Column copy
1400                            const double * element = solver->getMatrixByCol()->getElements();
1401                            const int * row = solver->getMatrixByCol()->getIndices();
1402                            const CoinBigIndex * columnStart = solver->getMatrixByCol()->getVectorStarts();
1403                            const int * columnLength = solver->getMatrixByCol()->getVectorLengths();
1404                            int n = 0;
1405                            double contrib = 0.0;
1406                            for (  i = 0 ; i < numberColumns ; i++ ) {
1407                                double value = newSolution[i];
1408                                if (solver->isInteger(i)) {
1409                                    double nearest = floor(value + 0.5);
1410                                    if (fabs(value - nearest) > 1.0e-6) {
1411                                        //printf("Column %d value %g\n",i,value);
1412                                        for (CoinBigIndex j = columnStart[i];
1413                                                j < columnStart[i] + columnLength[i]; j++) {
1414                                            int iRow = row[j];
1415                                            //printf("row %d element %g\n",iRow,element[j]);
1416                                            if (!rowValue[iRow]) {
1417                                                rowValue[iRow] = 1.0;
1418                                                whichRow[nRow++] = iRow;
1419                                            }
1420                                        }
1421                                        baseValue[n] = floor(value);
1422                                        contrib += saveObjective[i] * value;
1423                                        value = 0.0;
1424                                        stack[n] = 0;
1425                                        which[n++] = i;
1426                                    }
1427                                }
1428                                for (CoinBigIndex j = columnStart[i];
1429                                        j < columnStart[i] + columnLength[i]; j++) {
1430                                    int iRow = row[j];
1431                                    rowActivity[iRow] += value * element[j];
1432                                }
1433                            }
1434                            if (newNumberInfeas < 15) {
1435                                stack[n] = newNumberInfeas + 100;
1436                                int iStack = n;
1437                                memset(rowValue, 0, numberRows*sizeof(double));
1438                                const double * rowLower = solver->getRowLower();
1439                                const double * rowUpper = solver->getRowUpper();
1440                                while (iStack >= 0) {
1441                                    double contrib2 = 0.0;
1442                                    // Could do faster
1443                                    for (int k = 0 ; k < n ; k++ ) {
1444                                        i = which[k];
1445                                        double value = baseValue[k] + stack[k];
1446                                        contrib2 += saveObjective[i] * value;
1447                                        for (CoinBigIndex j = columnStart[i];
1448                                                j < columnStart[i] + columnLength[i]; j++) {
1449                                            int iRow = row[j];
1450                                            rowValue[iRow] += value * element[j];
1451                                        }
1452                                    }
1453                                    // check if feasible
1454                                    bool feasible = true;
1455                                    for (int k = 0; k < nRow; k++) {
1456                                        i = whichRow[k];
1457                                        double value = rowValue[i] + rowActivity[i];
1458                                        rowValue[i] = 0.0;
1459                                        if (value < rowLower[i] - 1.0e-7 ||
1460                                                value > rowUpper[i] + 1.0e-7)
1461                                            feasible = false;
1462                                    }
1463                                    if (feasible) {
1464                                        double newObj = newTrueSolutionValue * direction;
1465                                        newObj += contrib2 - contrib;
1466                                        newObj *= direction;
1467#ifdef COIN_DEVELOP
1468                                        printf("FFFeasible! - obj %g\n", newObj);
1469#endif
1470                                        if (newObj < roundingObjective - 1.0e-6) {
1471#ifdef COIN_DEVELOP
1472                                            printf("FBetter\n");
1473#endif
1474                                            roundingObjective = newObj;
1475                                            memcpy(roundingSolution, newSolution, numberColumns*sizeof(double));
1476                                            for (int k = 0 ; k < n ; k++ ) {
1477                                                i = which[k];
1478                                                double value = baseValue[k] + stack[k];
1479                                                roundingSolution[i] = value;
1480                                            }
1481                                        }
1482                                    }
1483                                    while (iStack >= 0 && stack[iStack]) {
1484                                        stack[iStack]--;
1485                                        iStack--;
1486                                    }
1487                                    if (iStack >= 0) {
1488                                        stack[iStack] = 1;
1489                                        iStack = n;
1490                                        stack[n] = 1;
1491                                    }
1492                                }
1493                            }
1494                            delete [] rowActivity;
1495                            delete [] which;
1496                            delete [] stack;
1497                            delete [] baseValue;
1498                            delete [] whichRow;
1499                            delete [] rowValue;
1500                        }
1501                    }
1502                    if (true) {
1503                        OsiSolverInterface * saveSolver = model_->swapSolver(clonedSolver);
1504                        clonedSolver->setColSolution(solver->getColSolution());
1505                        CbcRounding heuristic1(*model_);
1506                        heuristic1.setHeuristicName("rounding in feaspump!");
1507                        heuristic1.setWhen(1);
1508                        roundingObjective = CoinMin(roundingObjective, solutionValue);
1509                        double testSolutionValue = newTrueSolutionValue;
1510                        int returnCode = heuristic1.solution(roundingObjective,
1511                                                             roundingSolution,
1512                                                             testSolutionValue) ;
1513                        if (returnCode == 1) {
1514#ifdef COIN_DEVELOP
1515                            printf("rounding obj of %g?\n", roundingObjective);
1516#endif
1517                            //roundingObjective = newSolutionValue;
1518                            //} else {
1519                            //roundingObjective = COIN_DBL_MAX;
1520                        }
1521                        model_->swapSolver(saveSolver);
1522                    }
1523                    if (!solver->isProvenOptimal()) {
1524                        // presumably max time or some such
1525                        exitAll = true;
1526                        break;
1527                    }
1528                    // in case very dubious solver
1529                    lower = solver->getColLower();
1530                    upper = solver->getColUpper();
1531                    solution = solver->getColSolution();
1532                    numberIterations = solver->getIterationCount();
1533                } else {
1534                    int * addStart = new int[2*general+1];
1535                    int * addIndex = new int[4*general];
1536                    double * addElement = new double[4*general];
1537                    double * addLower = new double[2*general];
1538                    double * addUpper = new double[2*general];
1539                    double * obj = new double[general];
1540                    int nAdd = 0;
1541                    for (i = 0; i < numberIntegers; i++) {
1542                        int iColumn = integerVariable[i];
1543                        if (newSolution[iColumn] > lower[iColumn] + primalTolerance &&
1544                                newSolution[iColumn] < upper[iColumn] - primalTolerance) {
1545                            assert (upper[iColumn] - lower[iColumn] > 1.00001);
1546                            obj[nAdd] = 1.0;
1547                            addLower[nAdd] = 0.0;
1548                            addUpper[nAdd] = COIN_DBL_MAX;
1549                            nAdd++;
1550                        }
1551                    }
1552                    OsiSolverInterface * solver2 = solver;
1553                    if (nAdd) {
1554                        CoinZeroN(addStart, nAdd + 1);
1555                        solver2 = solver->clone();
1556                        solver2->addCols(nAdd, addStart, NULL, NULL, addLower, addUpper, obj);
1557                        // feasible solution
1558                        double * sol = new double[nAdd+numberColumns];
1559                        memcpy(sol, solution, numberColumns*sizeof(double));
1560                        // now rows
1561                        int nAdd = 0;
1562                        int nEl = 0;
1563                        int nAddRow = 0;
1564                        for (i = 0; i < numberIntegers; i++) {
1565                            int iColumn = integerVariable[i];
1566                            if (newSolution[iColumn] > lower[iColumn] + primalTolerance &&
1567                                    newSolution[iColumn] < upper[iColumn] - primalTolerance) {
1568                                addLower[nAddRow] = -newSolution[iColumn];;
1569                                addUpper[nAddRow] = COIN_DBL_MAX;
1570                                addIndex[nEl] = iColumn;
1571                                addElement[nEl++] = -1.0;
1572                                addIndex[nEl] = numberColumns + nAdd;
1573                                addElement[nEl++] = 1.0;
1574                                nAddRow++;
1575                                addStart[nAddRow] = nEl;
1576                                addLower[nAddRow] = newSolution[iColumn];;
1577                                addUpper[nAddRow] = COIN_DBL_MAX;
1578                                addIndex[nEl] = iColumn;
1579                                addElement[nEl++] = 1.0;
1580                                addIndex[nEl] = numberColumns + nAdd;
1581                                addElement[nEl++] = 1.0;
1582                                nAddRow++;
1583                                addStart[nAddRow] = nEl;
1584                                sol[nAdd+numberColumns] = fabs(sol[iColumn] - newSolution[iColumn]);
1585                                nAdd++;
1586                            }
1587                        }
1588                        solver2->setColSolution(sol);
1589                        delete [] sol;
1590                        solver2->addRows(nAddRow, addStart, addIndex, addElement, addLower, addUpper);
1591                    }
1592                    delete [] addStart;
1593                    delete [] addIndex;
1594                    delete [] addElement;
1595                    delete [] addLower;
1596                    delete [] addUpper;
1597                    delete [] obj;
1598                    solver2->resolve();
1599                    if (!solver2->isProvenOptimal()) {
1600                        // presumably max time or some such
1601                        exitAll = true;
1602                        break;
1603                    }
1604                    //assert (solver2->isProvenOptimal());
1605                    if (nAdd) {
1606                        solver->setColSolution(solver2->getColSolution());
1607                        numberIterations = solver2->getIterationCount();
1608                        delete solver2;
1609                    } else {
1610                        numberIterations = solver->getIterationCount();
1611                    }
1612                    lower = solver->getColLower();
1613                    upper = solver->getColUpper();
1614                    solution = solver->getColSolution();
1615                    newTrueSolutionValue = -saveOffset;
1616                    newSumInfeas = 0.0;
1617                    newNumberInfeas = 0;
1618                    {
1619                        const double * newSolution = solver->getColSolution();
1620                        for (  i = 0 ; i < numberColumns ; i++ ) {
1621                            if (solver->isInteger(i)) {
1622                                double value = newSolution[i];
1623                                double nearest = floor(value + 0.5);
1624                                newSumInfeas += fabs(value - nearest);
1625                                if (fabs(value - nearest) > 1.0e-6)
1626                                    newNumberInfeas++;
1627                            }
1628                            newTrueSolutionValue += saveObjective[i] * newSolution[i];
1629                        }
1630                        newTrueSolutionValue *= direction;
1631                    }
1632                }
1633                if (lastMove != 1000000) {
1634                    if (newSumInfeas < lastSumInfeas) {
1635                        lastMove = numberPasses;
1636                        lastSumInfeas = newSumInfeas;
1637                    } else if (newSumInfeas > lastSumInfeas + 1.0e-5) {
1638                        lastMove = 1000000; // going up
1639                    }
1640                }
1641                totalNumberIterations += numberIterations;
1642                if (solver->getNumRows() < 3000)
1643                    sprintf(pumpPrint, "Pass %3d: suminf. %10.5f (%d) obj. %g iterations %d",
1644                            numberPasses + totalNumberPasses,
1645                            newSumInfeas, newNumberInfeas,
1646                            newTrueSolutionValue, numberIterations);
1647                else
1648                    sprintf(pumpPrint, "Pass %3d: (%.2f seconds) suminf. %10.5f (%d) obj. %g iterations %d", numberPasses + totalNumberPasses,
1649                            model_->getCurrentSeconds(), newSumInfeas, newNumberInfeas,
1650                            newTrueSolutionValue, numberIterations);
1651                model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
1652                << pumpPrint
1653                << CoinMessageEol;
1654                CbcEventHandler *eventHandler = model_->getEventHandler() ;
1655                if (eventHandler) {
1656                  typedef struct {
1657                    double newSumInfeas;
1658                    double trueSolutionValue;
1659                    double spareDouble[2];
1660                    OsiSolverInterface * solver;
1661                    void * sparePointer[2];
1662                    int numberPasses;
1663                    int totalNumberPasses;
1664                    int numberInfeas;
1665                    int numberIterations;
1666                    int spareInt[3];
1667                  } HeurPass;
1668                  HeurPass temp;
1669                  temp.solver=solver;
1670                  temp.newSumInfeas = newSumInfeas;
1671                  temp.trueSolutionValue = newTrueSolutionValue;
1672                  temp.numberPasses=numberPasses;
1673                  temp.totalNumberPasses=totalNumberPasses;
1674                  temp.numberInfeas=newNumberInfeas;
1675                  temp.numberIterations=numberIterations;
1676                  CbcEventHandler::CbcAction status = 
1677                    eventHandler->event(CbcEventHandler::heuristicPass,
1678                                        &temp);
1679                  if (status==CbcEventHandler::killSolution) {
1680                    exitAll = true;
1681                    break;
1682                  }
1683                }
1684                if (closestSolution && solver->getObjValue() < closestObjectiveValue) {
1685                    int i;
1686                    const double * objective = solver->getObjCoefficients();
1687                    for (i = 0; i < numberIntegersOrig; i++) {
1688                        int iColumn = integerVariableOrig[i];
1689                        if (objective[iColumn] > 0.0)
1690                            closestSolution[i] = 0;
1691                        else
1692                            closestSolution[i] = 1;
1693                    }
1694                    closestObjectiveValue = solver->getObjValue();
1695                }
1696                // See if we need to think about changing rhs
1697                if ((switches_&12) != 0 && useRhs < 1.0e50) {
1698                    double oldRhs = useRhs;
1699                    bool trying = false;
1700                    if ((switches_&4) != 0 && numberPasses && (numberPasses % 50) == 0) {
1701                        if (solutionValue > 1.0e20) {
1702                            // only if no genuine solution
1703                            double gap = useRhs - continuousObjectiveValue;
1704                            useRhs += 0.1 * gap;
1705                            if (exactMultiple) {
1706                                useRhs = exactMultiple * ceil(useRhs / exactMultiple);
1707                                useRhs = CoinMax(useRhs, oldRhs + exactMultiple);
1708                            }
1709                            trying = true;
1710                        }
1711                    }
1712                    if ((switches_&8) != 0) {
1713                        // Put in new suminf and check
1714                        double largest = newSumInfeas;
1715                        double smallest = newSumInfeas;
1716                        for (int i = 0; i < SIZE_BOBBLE - 1; i++) {
1717                            double value = saveSumInf[i+1];
1718                            saveSumInf[i] = value;
1719                            largest = CoinMax(largest, value);
1720                            smallest = CoinMin(smallest, value);
1721                        }
1722                        saveSumInf[SIZE_BOBBLE-1] = newSumInfeas;
1723                        if (smallest*1.5 > largest && smallest > 2.0) {
1724                            if (bobbleMode == 0) {
1725                                // go closer
1726                                double gap = oldRhs - continuousObjectiveValue;
1727                                useRhs -= 0.4 * gap;
1728                                if (exactMultiple) {
1729                                    double value = floor(useRhs / exactMultiple);
1730                                    useRhs = CoinMin(value * exactMultiple, oldRhs - exactMultiple);
1731                                }
1732                                if (useRhs < continuousObjectiveValue) {
1733                                    // skip decrease
1734                                    bobbleMode = 1;
1735                                    useRhs = oldRhs;
1736                                }
1737                            }
1738                            if (bobbleMode) {
1739                                trying = true;
1740                                // weaken
1741                                if (solutionValue < 1.0e20) {
1742                                    double gap = solutionValue - oldRhs;
1743                                    useRhs += 0.3 * gap;
1744                                } else {
1745                                    double gap = oldRhs - continuousObjectiveValue;
1746                                    useRhs += 0.05 * gap;
1747                                }
1748                                if (exactMultiple) {
1749                                    double value = ceil(useRhs / exactMultiple);
1750                                    useRhs = CoinMin(value * exactMultiple,
1751                                                     solutionValue - exactMultiple);
1752                                }
1753                            }
1754                            bobbleMode++;
1755                            // reset
1756                            CoinFillN(saveSumInf, SIZE_BOBBLE, COIN_DBL_MAX);
1757                        }
1758                    }
1759                    if (useRhs != oldRhs) {
1760                        // tidy up
1761                        if (exactMultiple) {
1762                            double value = floor(useRhs / exactMultiple);
1763                            double bestPossible = ceil(continuousObjectiveValue / exactMultiple);
1764                            useRhs = CoinMax(value, bestPossible) * exactMultiple;
1765                        } else {
1766                            useRhs = CoinMax(useRhs, continuousObjectiveValue);
1767                        }
1768                        int k = solver->getNumRows() - 1;
1769                        solver->setRowUpper(k, useRhs + useOffset);
1770                        bool takeHint;
1771                        OsiHintStrength strength;
1772                        solver->getHintParam(OsiDoDualInResolve, takeHint, strength);
1773                        if (useRhs < oldRhs) {
1774                            solver->setHintParam(OsiDoDualInResolve, true);
1775                            solver->resolve();
1776                        } else if (useRhs > oldRhs) {
1777                            solver->setHintParam(OsiDoDualInResolve, false);
1778                            solver->resolve();
1779                        }
1780                        solver->setHintParam(OsiDoDualInResolve, takeHint);
1781                        if (!solver->isProvenOptimal()) {
1782                            // presumably max time or some such
1783                            exitAll = true;
1784                            break;
1785                        }
1786                    } else if (trying) {
1787                        // doesn't look good
1788                        break;
1789                    }
1790                }
1791            }
1792            // reduce scale factor
1793            scaleFactor *= weightFactor_;
1794        } // END WHILE
1795        // see if rounding worked!
1796        if (roundingObjective < solutionValue) {
1797            if (roundingObjective < solutionValue - 1.0e-6*fabs(roundingObjective)) {
1798                sprintf(pumpPrint, "Rounding solution of %g is better than previous of %g\n",
1799                        roundingObjective, solutionValue);
1800                model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
1801                << pumpPrint
1802                << CoinMessageEol;
1803            }
1804            solutionValue = roundingObjective;
1805            newSolutionValue = solutionValue;
1806            memcpy(betterSolution, roundingSolution, numberColumns*sizeof(double));
1807            solutionFound = true;
1808            numberSolutions++;
1809            if (numberSolutions>=maxSolutions)
1810              exitAll = true;
1811            if (exitNow(roundingObjective))
1812                exitAll = true;
1813        }
1814        if (!solutionFound) {
1815            sprintf(pumpPrint, "No solution found this major pass");
1816            model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
1817            << pumpPrint
1818            << CoinMessageEol;
1819        }
1820        //}
1821        delete solver;
1822        solver = NULL;
1823        for ( j = 0; j < NUMBER_OLD; j++)
1824            delete [] oldSolution[j];
1825        delete [] oldSolution;
1826        delete [] saveObjective;
1827        if (usedColumn && !exitAll) {
1828            OsiSolverInterface * newSolver = cloneBut(3); // was model_->continuousSolver()->clone();
1829#if 0 //def COIN_HAS_CLP
1830            OsiClpSolverInterface * clpSolver
1831              = dynamic_cast<OsiClpSolverInterface *> (newSolver);
1832            if (clpSolver) {
1833              ClpSimplex * simplex = clpSolver->getModelPtr();
1834              simplex->writeMps("start.mps",2,1);
1835            }
1836#endif
1837            const double * colLower = newSolver->getColLower();
1838            const double * colUpper = newSolver->getColUpper();
1839            bool stopBAB = false;
1840            int allowedPass = -1;
1841            if (maximumAllowed > 0)
1842                allowedPass = CoinMax(numberPasses - maximumAllowed, -1);
1843            while (!stopBAB) {
1844                stopBAB = true;
1845                int i;
1846                int nFix = 0;
1847                int nFixI = 0;
1848                int nFixC = 0;
1849                int nFixC2 = 0;
1850                for (i = 0; i < numberIntegersOrig; i++) {
1851                    int iColumn = integerVariableOrig[i];
1852                    //const OsiObject * object = model_->object(i);
1853                    //double originalLower;
1854                    //double originalUpper;
1855                    //getIntegerInformation( object,originalLower, originalUpper);
1856                    //assert(colLower[iColumn]==originalLower);
1857                    //newSolver->setColLower(iColumn,CoinMax(colLower[iColumn],originalLower));
1858                    newSolver->setColLower(iColumn, colLower[iColumn]);
1859                    //assert(colUpper[iColumn]==originalUpper);
1860                    //newSolver->setColUpper(iColumn,CoinMin(colUpper[iColumn],originalUpper));
1861                    newSolver->setColUpper(iColumn, colUpper[iColumn]);
1862                    if (usedColumn[iColumn] <= allowedPass) {
1863                        double value = lastSolution[iColumn];
1864                        double nearest = floor(value + 0.5);
1865                        if (fabs(value - nearest) < 1.0e-7) {
1866                            if (nearest == colLower[iColumn]) {
1867                                newSolver->setColUpper(iColumn, colLower[iColumn]);
1868                                nFix++;
1869                            } else if (nearest == colUpper[iColumn]) {
1870                                newSolver->setColLower(iColumn, colUpper[iColumn]);
1871                                nFix++;
1872                            } else if (fixInternal) {
1873                                newSolver->setColLower(iColumn, nearest);
1874                                newSolver->setColUpper(iColumn, nearest);
1875                                nFix++;
1876                                nFixI++;
1877                            }
1878                        }
1879                    }
1880                }
1881                if (fixContinuous) {
1882                    for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
1883                        if (!newSolver->isInteger(iColumn) && usedColumn[iColumn] <= allowedPass) {
1884                            double value = lastSolution[iColumn];
1885                            if (value < colLower[iColumn] + 1.0e-8) {
1886                                newSolver->setColUpper(iColumn, colLower[iColumn]);
1887                                nFixC++;
1888                            } else if (value > colUpper[iColumn] - 1.0e-8) {
1889                                newSolver->setColLower(iColumn, colUpper[iColumn]);
1890                                nFixC++;
1891                            } else if (fixContinuous == 2) {
1892                                newSolver->setColLower(iColumn, value);
1893                                newSolver->setColUpper(iColumn, value);
1894                                nFixC++;
1895                                nFixC2++;
1896                            }
1897                        }
1898                    }
1899                }
1900                newSolver->initialSolve();
1901                if (!newSolver->isProvenOptimal()) {
1902                    //newSolver->writeMps("bad.mps");
1903                    //assert (newSolver->isProvenOptimal());
1904                    exitAll = true;
1905                    break;
1906                }
1907                sprintf(pumpPrint, "Before mini branch and bound, %d integers at bound fixed and %d continuous",
1908                        nFix, nFixC);
1909                if (nFixC2 + nFixI != 0)
1910                    sprintf(pumpPrint + strlen(pumpPrint), " of which %d were internal integer and %d internal continuous",
1911                            nFixI, nFixC2);
1912                model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
1913                << pumpPrint
1914                << CoinMessageEol;
1915                double saveValue = newSolutionValue;
1916                if (newSolutionValue - model_->getCutoffIncrement()
1917                        > continuousObjectiveValue - 1.0e-7) {
1918                    double saveFraction = fractionSmall_;
1919                    if (numberTries > 1 && !numberBandBsolutions)
1920                        fractionSmall_ *= 0.5;
1921                    // Give branch and bound a bit more freedom
1922                    double cutoff2 = newSolutionValue +
1923                                     CoinMax(model_->getCutoffIncrement(), 1.0e-3);
1924#if 0
1925                      {
1926                        OsiClpSolverInterface * clpSolver
1927                        = dynamic_cast<OsiClpSolverInterface *> (newSolver);
1928                        if (clpSolver) {
1929                            ClpSimplex * simplex = clpSolver->getModelPtr();
1930                            simplex->writeMps("testA.mps",2,1);
1931                        }
1932                      }
1933#endif
1934                    int returnCode2 = smallBranchAndBound(newSolver, numberNodes_, newSolution, newSolutionValue,
1935                                                          cutoff2, "CbcHeuristicLocalAfterFPump");
1936                    fractionSmall_ = saveFraction;
1937                    if (returnCode2 < 0) {
1938                        if (returnCode2 == -2) {
1939                            exitAll = true;
1940                            returnCode = 0;
1941                        } else {
1942                            returnCode2 = 0; // returned on size - try changing
1943                            //#define ROUND_AGAIN
1944#ifdef ROUND_AGAIN
1945                            if (numberTries == 1 && numberPasses > 20 && allowedPass < numberPasses - 1) {
1946                                allowedPass = (numberPasses + allowedPass) >> 1;
1947                                sprintf(pumpPrint,
1948                                        "Fixing all variables which were last changed on pass %d and trying again",
1949                                        allowedPass);
1950                                model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
1951                                << pumpPrint
1952                                << CoinMessageEol;
1953                                stopBAB = false;
1954                                continue;
1955                            }
1956#endif
1957                        }
1958                    }
1959                    if ((returnCode2&2) != 0) {
1960                        // could add cut
1961                        returnCode2 &= ~2;
1962                    }
1963                    if (returnCode2) {
1964                        numberBandBsolutions++;
1965                        // may not have got solution earlier
1966                        returnCode |= 1;
1967                    }
1968                } else {
1969                    // no need
1970                    exitAll = true;
1971                    //returnCode=0;
1972                }
1973                // recompute solution value
1974                if (returnCode && true) {
1975#if 0
1976                      {
1977                        OsiClpSolverInterface * clpSolver
1978                        = dynamic_cast<OsiClpSolverInterface *> (newSolver);
1979                        if (clpSolver) {
1980                            ClpSimplex * simplex = clpSolver->getModelPtr();
1981                            simplex->writeMps("testB.mps",2,1);
1982                        }
1983                      }
1984#endif
1985                    delete newSolver;
1986                    newSolver = cloneBut(3); // was model_->continuousSolver()->clone();
1987                    newSolutionValue = -saveOffset;
1988                    double newSumInfeas = 0.0;
1989                    const double * obj = newSolver->getObjCoefficients();
1990                    for (int i = 0 ; i < numberColumns ; i++ ) {
1991                        if (newSolver->isInteger(i)) {
1992                            double value = newSolution[i];
1993                            double nearest = floor(value + 0.5);
1994                            newSumInfeas += fabs(value - nearest);
1995                        }
1996                        newSolutionValue += obj[i] * newSolution[i];
1997                    }
1998                    newSolutionValue *= direction;
1999                }
2000                bool gotSolution = false;
2001                if (returnCode && newSolutionValue < saveValue) {
2002                    sprintf(pumpPrint, "Mini branch and bound improved solution from %g to %g (%.2f seconds)",
2003                            saveValue, newSolutionValue, model_->getCurrentSeconds());
2004                    model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
2005                    << pumpPrint
2006                    << CoinMessageEol;
2007                    memcpy(betterSolution, newSolution, numberColumns*sizeof(double));
2008                    gotSolution = true;
2009                    if (fixContinuous && nFixC + nFixC2 > 0) {
2010                        // may be able to do even better
2011                        int nFixed = 0;
2012                        const double * lower = model_->solver()->getColLower();
2013                        const double * upper = model_->solver()->getColUpper();
2014                        for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
2015                            double value = newSolution[iColumn];
2016                            if (newSolver->isInteger(iColumn)) {
2017                                value = floor(newSolution[iColumn] + 0.5);
2018                                newSolver->setColLower(iColumn, value);
2019                                newSolver->setColUpper(iColumn, value);
2020                                nFixed++;
2021                            } else {
2022                                newSolver->setColLower(iColumn, lower[iColumn]);
2023                                newSolver->setColUpper(iColumn, upper[iColumn]);
2024                                if (value < lower[iColumn])
2025                                    value = lower[iColumn];
2026                                else if (value > upper[iColumn])
2027                                    value = upper[iColumn];
2028
2029                            }
2030                            newSolution[iColumn] = value;
2031                        }
2032                        newSolver->setColSolution(newSolution);
2033                        //#define CLP_INVESTIGATE2
2034#ifdef CLP_INVESTIGATE2
2035                        {
2036                            // check
2037                            // get row activities
2038                            int numberRows = newSolver->getNumRows();
2039                            double * rowActivity = new double[numberRows];
2040                            memset(rowActivity, 0, numberRows*sizeof(double));
2041                            newSolver->getMatrixByCol()->times(newSolution, rowActivity) ;
2042                            double largestInfeasibility = primalTolerance;
2043                            double sumInfeasibility = 0.0;
2044                            int numberBadRows = 0;
2045                            const double * rowLower = newSolver->getRowLower();
2046                            const double * rowUpper = newSolver->getRowUpper();
2047                            for (i = 0 ; i < numberRows ; i++) {
2048                                double value;
2049                                value = rowLower[i] - rowActivity[i];
2050                                if (value > primalTolerance) {
2051                                    numberBadRows++;
2052                                    largestInfeasibility = CoinMax(largestInfeasibility, value);
2053                                    sumInfeasibility += value;
2054                                }
2055                                value = rowActivity[i] - rowUpper[i];
2056                                if (value > primalTolerance) {
2057                                    numberBadRows++;
2058                                    largestInfeasibility = CoinMax(largestInfeasibility, value);
2059                                    sumInfeasibility += value;
2060                                }
2061                            }
2062                            printf("%d bad rows, largest inf %g sum %g\n",
2063                                   numberBadRows, largestInfeasibility, sumInfeasibility);
2064                            delete [] rowActivity;
2065                        }
2066#endif
2067#ifdef COIN_HAS_CLP
2068                        OsiClpSolverInterface * clpSolver
2069                        = dynamic_cast<OsiClpSolverInterface *> (newSolver);
2070                        if (clpSolver) {
2071                            ClpSimplex * simplex = clpSolver->getModelPtr();
2072                            //simplex->writeBasis("test.bas",true,2);
2073                            //simplex->writeMps("test.mps",2,1);
2074                            if (nFixed*3 > numberColumns*2)
2075                                simplex->allSlackBasis(); // may as well go from all slack
2076                            simplex->primal(1);
2077                            clpSolver->setWarmStart(NULL);
2078                        }
2079#endif
2080                        newSolver->initialSolve();
2081                        if (newSolver->isProvenOptimal()) {
2082                            double value = newSolver->getObjValue() * newSolver->getObjSense();
2083                            if (value < newSolutionValue) {
2084                                //newSolver->writeMps("query","mps");
2085#ifdef JJF_ZERO
2086                                {
2087                                    double saveOffset;
2088                                    newSolver->getDblParam(OsiObjOffset, saveOffset);
2089                                    const double * obj = newSolver->getObjCoefficients();
2090                                    double newTrueSolutionValue = -saveOffset;
2091                                    double newSumInfeas = 0.0;
2092                                    int numberColumns = newSolver->getNumCols();
2093                                    const double * solution = newSolver->getColSolution();
2094                                    for (int  i = 0 ; i < numberColumns ; i++ ) {
2095                                        if (newSolver->isInteger(i)) {
2096                                            double value = solution[i];
2097                                            double nearest = floor(value + 0.5);
2098                                            newSumInfeas += fabs(value - nearest);
2099                                        }
2100                                        if (solution[i])
2101                                            printf("%d obj %g val %g - total %g\n", i, obj[i], solution[i],
2102                                                   newTrueSolutionValue);
2103                                        newTrueSolutionValue += obj[i] * solution[i];
2104                                    }
2105                                    printf("obj %g - inf %g\n", newTrueSolutionValue,
2106                                           newSumInfeas);
2107                                }
2108#endif
2109                                sprintf(pumpPrint, "Freeing continuous variables gives a solution of %g", value);
2110                                model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
2111                                << pumpPrint
2112                                << CoinMessageEol;
2113                                newSolutionValue = value;
2114                                memcpy(betterSolution, newSolver->getColSolution(), numberColumns*sizeof(double));
2115                            }
2116                        } else {
2117                            //newSolver->writeMps("bad3.mps");
2118                          sprintf(pumpPrint, "On closer inspection solution is not valid");
2119                          model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
2120                            << pumpPrint
2121                            << CoinMessageEol;
2122                            exitAll = true;
2123                            break;
2124                        }
2125                    }
2126                } else {
2127                    sprintf(pumpPrint, "Mini branch and bound did not improve solution (%.2f seconds)",
2128                            model_->getCurrentSeconds());
2129                    model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
2130                    << pumpPrint
2131                    << CoinMessageEol;
2132                    if (returnCode && newSolutionValue < saveValue + 1.0e-3 && nFixC + nFixC2) {
2133                        // may be able to do better
2134                        const double * lower = model_->solver()->getColLower();
2135                        const double * upper = model_->solver()->getColUpper();
2136                        for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
2137                            if (newSolver->isInteger(iColumn)) {
2138                                double value = floor(newSolution[iColumn] + 0.5);
2139                                newSolver->setColLower(iColumn, value);
2140                                newSolver->setColUpper(iColumn, value);
2141                            } else {
2142                                newSolver->setColLower(iColumn, lower[iColumn]);
2143                                newSolver->setColUpper(iColumn, upper[iColumn]);
2144                            }
2145                        }
2146                        newSolver->initialSolve();
2147                        if (newSolver->isProvenOptimal()) {
2148                            double value = newSolver->getObjValue() * newSolver->getObjSense();
2149                            if (value < saveValue) {
2150                                sprintf(pumpPrint, "Freeing continuous variables gives a solution of %g", value);
2151                                model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
2152                                << pumpPrint
2153                                << CoinMessageEol;
2154                                newSolutionValue = value;
2155                                memcpy(betterSolution, newSolver->getColSolution(), numberColumns*sizeof(double));
2156                                gotSolution = true;
2157                            }
2158                        }
2159                    }
2160                }
2161                if (gotSolution) {
2162                    if ((accumulate_&1) != 0) {
2163                        model_->incrementUsed(betterSolution); // for local search
2164                    }
2165                    solutionValue = newSolutionValue;
2166                    solutionFound = true;
2167                    numberSolutions++;
2168                    if (numberSolutions>=maxSolutions)
2169                      exitAll = true;
2170                    if (exitNow(newSolutionValue))
2171                        exitAll = true;
2172                    CoinWarmStartBasis * basis =
2173                        dynamic_cast<CoinWarmStartBasis *>(newSolver->getWarmStart()) ;
2174                    if (basis) {
2175                        bestBasis = * basis;
2176                        delete basis;
2177                        int action = model_->dealWithEventHandler(CbcEventHandler::heuristicSolution, newSolutionValue, betterSolution);
2178                        if (action == 0) {
2179                            double * saveOldSolution = CoinCopyOfArray(model_->bestSolution(), numberColumns);
2180                            double saveObjectiveValue = model_->getMinimizationObjValue();
2181                            model_->setBestSolution(betterSolution, numberColumns, newSolutionValue);
2182                            if (saveOldSolution && saveObjectiveValue < model_->getMinimizationObjValue())
2183                                model_->setBestSolution(saveOldSolution, numberColumns, saveObjectiveValue);
2184                            delete [] saveOldSolution;
2185                        }
2186                        if (!action || model_->getCurrentSeconds() > model_->getMaximumSeconds()) {
2187                            exitAll = true; // exit
2188                            break;
2189                        }
2190                    }
2191                }
2192            } // end stopBAB while
2193            delete newSolver;
2194        }
2195        if (solutionFound) finalReturnCode = 1;
2196        cutoff = CoinMin(cutoff, solutionValue - model_->getCutoffIncrement());
2197        if (numberTries >= maximumRetries_ || !solutionFound || exitAll || cutoff < continuousObjectiveValue + 1.0e-7) {
2198            break;
2199        } else {
2200            solutionFound = false;
2201            if (absoluteIncrement_ > 0.0 || relativeIncrement_ > 0.0) {
2202                double gap = relativeIncrement_ * fabs(solutionValue);
2203                double change = CoinMax(gap, absoluteIncrement_);
2204                cutoff = CoinMin(cutoff, solutionValue - change);
2205            } else {
2206                //double weights[10]={0.1,0.1,0.2,0.2,0.2,0.3,0.3,0.3,0.4,0.5};
2207                double weights[10] = {0.1, 0.2, 0.3, 0.3, 0.4, 0.4, 0.4, 0.5, 0.5, 0.6};
2208                cutoff -= weights[CoinMin(numberTries-1, 9)] * (cutoff - continuousObjectiveValue);
2209            }
2210            // But round down
2211            if (exactMultiple)
2212                cutoff = exactMultiple * floor(cutoff / exactMultiple);
2213            if (cutoff < continuousObjectiveValue)
2214                break;
2215            sprintf(pumpPrint, "Round again with cutoff of %g", cutoff);
2216            model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
2217            << pumpPrint
2218            << CoinMessageEol;
2219            if ((accumulate_&3) < 2 && usedColumn) {
2220                for (int i = 0; i < numberColumns; i++)
2221                    usedColumn[i] = -1;
2222            }
2223            averageIterationsPerTry = totalNumberIterations / numberTries;
2224            numberIterationsLastPass =  totalNumberIterations;
2225            totalNumberPasses += numberPasses - 1;
2226        }
2227    }
2228/*
2229  End of the `exitAll' loop.
2230*/
2231#ifdef RAND_RAND
2232    delete [] randomFactor;
2233#endif
2234    delete solver; // probably NULL but do anyway
2235    if (!finalReturnCode && closestSolution && closestObjectiveValue <= 10.0 &&
2236            usedColumn && !model_->maximumSecondsReached()) {
2237        // try a bit of branch and bound
2238        OsiSolverInterface * newSolver = cloneBut(1); // was model_->continuousSolver()->clone();
2239        const double * colLower = newSolver->getColLower();
2240        const double * colUpper = newSolver->getColUpper();
2241        int i;
2242        double rhs = 0.0;
2243        for (i = 0; i < numberIntegersOrig; i++) {
2244            int iColumn = integerVariableOrig[i];
2245            int direction = closestSolution[i];
2246            closestSolution[i] = iColumn;
2247            if (direction == 0) {
2248                // keep close to LB
2249                rhs += colLower[iColumn];
2250                lastSolution[i] = 1.0;
2251            } else {
2252                // keep close to UB
2253                rhs -= colUpper[iColumn];
2254                lastSolution[i] = -1.0;
2255            }
2256        }
2257        newSolver->addRow(numberIntegersOrig, closestSolution,
2258                          lastSolution, -COIN_DBL_MAX, rhs + 10.0);
2259        //double saveValue = newSolutionValue;
2260        //newSolver->writeMps("sub");
2261        int returnCode = smallBranchAndBound(newSolver, numberNodes_, newSolution, newSolutionValue,
2262                                             newSolutionValue, "CbcHeuristicLocalAfterFPump");
2263        if (returnCode < 0)
2264            returnCode = 0; // returned on size
2265        if ((returnCode&2) != 0) {
2266            // could add cut
2267            returnCode &= ~2;
2268        }
2269        if (returnCode) {
2270            //printf("old sol of %g new of %g\n",saveValue,newSolutionValue);
2271            memcpy(betterSolution, newSolution, numberColumns*sizeof(double));
2272            //abort();
2273            solutionValue = newSolutionValue;
2274            solutionFound = true;
2275            numberSolutions++;
2276            if (numberSolutions>=maxSolutions)
2277              exitAll = true;
2278            if (exitNow(newSolutionValue))
2279                exitAll = true;
2280        }
2281        delete newSolver;
2282    }
2283    delete clonedSolver;
2284    delete [] roundingSolution;
2285    delete [] usedColumn;
2286    delete [] lastSolution;
2287    delete [] newSolution;
2288    delete [] closestSolution;
2289    delete [] integerVariable;
2290    delete [] firstPerturbedObjective;
2291    delete [] firstPerturbedSolution;
2292    if (solutionValue == incomingObjective)
2293        sprintf(pumpPrint, "After %.2f seconds - Feasibility pump exiting - took %.2f seconds",
2294                model_->getCurrentSeconds(), CoinCpuTime() - time1);
2295    else if (numberSolutions < maxSolutions)
2296      sprintf(pumpPrint, "After %.2f seconds - Feasibility pump exiting with objective of %g - took %.2f seconds",
2297              model_->getCurrentSeconds(), solutionValue, CoinCpuTime() - time1);
2298    else 
2299        sprintf(pumpPrint, "After %.2f seconds - Feasibility pump exiting with objective of %g (stopping after %d solutions) - took %.2f seconds",
2300                model_->getCurrentSeconds(), solutionValue, 
2301                numberSolutions,CoinCpuTime() - time1);
2302    model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
2303      << pumpPrint
2304      << CoinMessageEol;
2305    if (bestBasis.getNumStructural())
2306        model_->setBestSolutionBasis(bestBasis);
2307    //model_->setMinimizationObjValue(saveBestObjective);
2308    if ((accumulate_&1) != 0 && numberSolutions > 1 && !model_->getSolutionCount()) {
2309        model_->setSolutionCount(1); // for local search
2310        model_->setNumberHeuristicSolutions(1);
2311    }
2312#ifdef COIN_DEVELOP
2313    {
2314        double ncol = model_->solver()->getNumCols();
2315        double nrow = model_->solver()->getNumRows();
2316        printf("XXX total iterations %d ratios - %g %g %g\n",
2317               totalNumberIterations,
2318               static_cast<double> (totalNumberIterations) / nrow,
2319               static_cast<double> (totalNumberIterations) / ncol,
2320               static_cast<double> (totalNumberIterations) / (2*nrow + 2*ncol));
2321    }
2322#endif
2323    return finalReturnCode;
2324}
2325
2326/**************************END MAIN PROCEDURE ***********************************/
2327
2328// update model
2329void CbcHeuristicFPump::setModel(CbcModel * model)
2330{
2331    model_ = model;
2332}
2333
2334/* Rounds solution - down if < downValue
2335   returns 1 if current is a feasible solution
2336*/
2337int
2338CbcHeuristicFPump::rounds(OsiSolverInterface * solver, double * solution,
2339                          //const double * objective,
2340                          int numberIntegers, const int * integerVariable,
2341                          /*char * pumpPrint,*/ int iter,
2342                          /*bool roundExpensive,*/ double downValue, int *flip)
2343{
2344    double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
2345    double primalTolerance ;
2346    solver->getDblParam(OsiPrimalTolerance, primalTolerance) ;
2347
2348    int i;
2349
2350    const double * cost = solver->getObjCoefficients();
2351    int flip_up = 0;
2352    int flip_down  = 0;
2353    double  v = randomNumberGenerator_.randomDouble() * 20.0;
2354    int nn = 10 + static_cast<int> (v);
2355    int nnv = 0;
2356    int * list = new int [nn];
2357    double * val = new double [nn];
2358    for (i = 0; i < nn; i++) val[i] = .001;
2359
2360    const double * rowLower = solver->getRowLower();
2361    const double * rowUpper = solver->getRowUpper();
2362    int numberRows = solver->getNumRows();
2363    if (false && (iter&1) != 0) {
2364        // Do set covering variables
2365        const CoinPackedMatrix * matrixByRow = solver->getMatrixByRow();
2366        const double * elementByRow = matrixByRow->getElements();
2367        const int * column = matrixByRow->getIndices();
2368        const CoinBigIndex * rowStart = matrixByRow->getVectorStarts();
2369        const int * rowLength = matrixByRow->getVectorLengths();
2370        for (i = 0; i < numberRows; i++) {
2371            if (rowLower[i] == 1.0 && rowUpper[i] == 1.0) {
2372                bool cover = true;
2373                double largest = 0.0;
2374                int jColumn = -1;
2375                for (CoinBigIndex k = rowStart[i]; k < rowStart[i] + rowLength[i]; k++) {
2376                    int iColumn = column[k];
2377                    if (elementByRow[k] != 1.0 || !solver->isInteger(iColumn)) {
2378                        cover = false;
2379                        break;
2380                    } else {
2381                        if (solution[iColumn]) {
2382                            double value = solution[iColumn] *
2383                                           (randomNumberGenerator_.randomDouble() + 5.0);
2384                            if (value > largest) {
2385                                largest = value;
2386                                jColumn = iColumn;
2387                            }
2388                        }
2389                    }
2390                }
2391                if (cover) {
2392                    for (CoinBigIndex k = rowStart[i]; k < rowStart[i] + rowLength[i]; k++) {
2393                        int iColumn = column[k];
2394                        if (iColumn == jColumn)
2395                            solution[iColumn] = 1.0;
2396                        else
2397                            solution[iColumn] = 0.0;
2398                    }
2399                }
2400            }
2401        }
2402    }
2403    int numberColumns = solver->getNumCols();
2404#ifdef JJF_ZERO
2405    // Do set covering variables
2406    const CoinPackedMatrix * matrixByRow = solver->getMatrixByRow();
2407    const double * elementByRow = matrixByRow->getElements();
2408    const int * column = matrixByRow->getIndices();
2409    const CoinBigIndex * rowStart = matrixByRow->getVectorStarts();
2410    const int * rowLength = matrixByRow->getVectorLengths();
2411    double * sortTemp = new double[numberColumns];
2412    int * whichTemp = new int [numberColumns];
2413    char * rowTemp = new char [numberRows];
2414    memset(rowTemp, 0, numberRows);
2415    for (i = 0; i < numberColumns; i++)
2416        whichTemp[i] = -1;
2417    int nSOS = 0;
2418    for (i = 0; i < numberRows; i++) {
2419        if (rowLower[i] == 1.0 && rowUpper[i] == 1.0) {
2420            bool cover = true;
2421            for (CoinBigIndex k = rowStart[i]; k < rowStart[i] + rowLength[i]; k++) {
2422                int iColumn = column[k];
2423                if (elementByRow[k] != 1.0 || !solver->isInteger(iColumn)) {
2424                    cover = false;
2425                    break;
2426                }
2427            }
2428            if (cover) {
2429                rowTemp[i] = 1;
2430                nSOS++;
2431                for (CoinBigIndex k = rowStart[i]; k < rowStart[i] + rowLength[i]; k++) {
2432                    int iColumn = column[k];
2433                    double value = solution[iColumn];
2434                    whichTemp[iColumn] = iColumn;
2435                }
2436            }
2437        }
2438    }
2439    if (nSOS) {
2440        // Column copy
2441        const CoinPackedMatrix * matrixByColumn = solver->getMatrixByCol();
2442        //const double * element = matrixByColumn->getElements();
2443        const int * row = matrixByColumn->getIndices();
2444        const CoinBigIndex * columnStart = matrixByColumn->getVectorStarts();
2445        const int * columnLength = matrixByColumn->getVectorLengths();
2446        int nLook = 0;
2447        for (i = 0; i < numberColumns; i++) {
2448            if (whichTemp[i] >= 0) {
2449                whichTemp[nLook] = i;
2450                double value = solution[i];
2451                if (value < 0.5)
2452                    value *= (0.1 * randomNumberGenerator_.randomDouble() + 0.3);
2453                sortTemp[nLook++] = -value;
2454            }
2455        }
2456        CoinSort_2(sortTemp, sortTemp + nLook, whichTemp);
2457        double smallest = 1.0;
2458        int nFix = 0;
2459        int nOne = 0;
2460        for (int j = 0; j < nLook; j++) {
2461            int jColumn = whichTemp[j];
2462            double thisValue = solution[jColumn];
2463            if (!thisValue)
2464                continue;
2465            if (thisValue == 1.0)
2466                nOne++;
2467            smallest = CoinMin(smallest, thisValue);
2468            solution[jColumn] = 1.0;
2469            double largest = 0.0;
2470            for (CoinBigIndex jEl = columnStart[jColumn];
2471                    jEl < columnStart[jColumn] + columnLength[jColumn]; jEl++) {
2472                int jRow = row[jEl];
2473                if (rowTemp[jRow]) {
2474                    for (CoinBigIndex k = rowStart[jRow]; k < rowStart[jRow] + rowLength[jRow]; k++) {
2475                        int iColumn = column[k];
2476                        if (solution[iColumn]) {
2477                            if (iColumn != jColumn) {
2478                                double value = solution[iColumn];
2479                                if (value > largest)
2480                                    largest = value;
2481                                solution[iColumn] = 0.0;
2482                            }
2483                        }
2484                    }
2485                }
2486            }
2487            if (largest > thisValue)
2488                printf("%d was at %g - chosen over a value of %g\n",
2489                       jColumn, thisValue, largest);
2490            nFix++;
2491        }
2492        printf("%d fixed out of %d (%d at one already)\n",
2493               nFix, nLook, nOne);
2494    }
2495    delete [] sortTemp;
2496    delete [] whichTemp;
2497    delete [] rowTemp;
2498#endif
2499    const double * columnLower = solver->getColLower();
2500    const double * columnUpper = solver->getColUpper();
2501    // Check if valid with current solution (allow for 0.99999999s)
2502    double newSumInfeas = 0.0;
2503    int newNumberInfeas = 0;
2504    for (i = 0; i < numberIntegers; i++) {
2505        int iColumn = integerVariable[i];
2506        double value = solution[iColumn];
2507        double round = floor(value + 0.5);
2508        if (fabs(value - round) > primalTolerance) {
2509          newSumInfeas += fabs(value-round);
2510          newNumberInfeas++;
2511        }
2512    }
2513    if (!newNumberInfeas) {
2514        // may be able to use solution even if 0.99999's
2515        double * saveLower = CoinCopyOfArray(columnLower, numberColumns);
2516        double * saveUpper = CoinCopyOfArray(columnUpper, numberColumns);
2517        double * saveSolution = CoinCopyOfArray(solution, numberColumns);
2518        double * tempSolution = CoinCopyOfArray(solution, numberColumns);
2519        CoinWarmStartBasis  * saveBasis =
2520            dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ;
2521        for (i = 0; i < numberIntegers; i++) {
2522            int iColumn = integerVariable[i];
2523            double value = solution[iColumn];
2524            double round = floor(value + 0.5);
2525            solver->setColLower(iColumn, round);
2526            solver->setColUpper(iColumn, round);
2527            tempSolution[iColumn] = round;
2528        }
2529        solver->setColSolution(tempSolution);
2530        delete [] tempSolution;
2531        solver->resolve();
2532        solver->setColLower(saveLower);
2533        solver->setColUpper(saveUpper);
2534        solver->setWarmStart(saveBasis);
2535        delete [] saveLower;
2536        delete [] saveUpper;
2537        delete saveBasis;
2538        if (!solver->isProvenOptimal()) {
2539            solver->setColSolution(saveSolution);
2540        }
2541        delete [] saveSolution;
2542        if (solver->isProvenOptimal()) {
2543            // feasible
2544            delete [] list;
2545            delete [] val;
2546            return 1;
2547        }
2548    }
2549    //double * saveSolution = CoinCopyOfArray(solution,numberColumns);
2550    // return rounded solution
2551    for (i = 0; i < numberIntegers; i++) {
2552        int iColumn = integerVariable[i];
2553        double value = solution[iColumn];
2554        double round = floor(value + primalTolerance);
2555        if (value - round > downValue) round += 1.;
2556#ifndef JJF_ONE
2557        if (round < integerTolerance && cost[iColumn] < -1. + integerTolerance) flip_down++;
2558        if (round > 1. - integerTolerance && cost[iColumn] > 1. - integerTolerance) flip_up++;
2559#else
2560        if (round < columnLower[iColumn] + integerTolerance && cost[iColumn] < -1. + integerTolerance) flip_down++;
2561        if (round > columnUpper[iColumn] - integerTolerance && cost[iColumn] > 1. - integerTolerance) flip_up++;
2562#endif
2563        if (flip_up + flip_down == 0) {
2564            for (int k = 0; k < nn; k++) {
2565                if (fabs(value - round) > val[k]) {
2566                    nnv++;
2567                    for (int j = nn - 2; j >= k; j--) {
2568                        val[j+1] = val[j];
2569                        list[j+1] = list[j];
2570                    }
2571                    val[k] = fabs(value - round);
2572                    list[k] = iColumn;
2573                    break;
2574                }
2575            }
2576        }
2577        solution[iColumn] = round;
2578    }
2579
2580    if (nnv > nn) nnv = nn;
2581    //if (iter != 0)
2582    //sprintf(pumpPrint+strlen(pumpPrint),"up = %5d , down = %5d", flip_up, flip_down);
2583    *flip = flip_up + flip_down;
2584
2585    if (*flip == 0 && iter != 0) {
2586        //sprintf(pumpPrint+strlen(pumpPrint)," -- rand = %4d (%4d) ", nnv, nn);
2587        for (i = 0; i < nnv; i++) {
2588            // was solution[list[i]] = 1. - solution[list[i]]; but does that work for 7>=x>=6
2589            int index = list[i];
2590            double value = solution[index];
2591            if (value <= 1.0) {
2592                solution[index] = 1.0 - value;
2593            } else if (value < columnLower[index] + integerTolerance) {
2594                solution[index] = value + 1.0;
2595            } else if (value > columnUpper[index] - integerTolerance) {
2596                solution[index] = value - 1.0;
2597            } else {
2598                solution[index] = value - 1.0;
2599            }
2600        }
2601        *flip = nnv;
2602    } else {
2603        //sprintf(pumpPrint+strlen(pumpPrint)," ");
2604    }
2605    delete [] list;
2606    delete [] val;
2607    //iter++;
2608
2609    // get row activities
2610    double * rowActivity = new double[numberRows];
2611    memset(rowActivity, 0, numberRows*sizeof(double));
2612    solver->getMatrixByCol()->times(solution, rowActivity) ;
2613    double largestInfeasibility = primalTolerance;
2614    double sumInfeasibility = 0.0;
2615    int numberBadRows = 0;
2616    for (i = 0 ; i < numberRows ; i++) {
2617        double value;
2618        value = rowLower[i] - rowActivity[i];
2619        if (value > primalTolerance) {
2620            numberBadRows++;
2621            largestInfeasibility = CoinMax(largestInfeasibility, value);
2622            sumInfeasibility += value;
2623        }
2624        value = rowActivity[i] - rowUpper[i];
2625        if (value > primalTolerance) {
2626            numberBadRows++;
2627            largestInfeasibility = CoinMax(largestInfeasibility, value);
2628            sumInfeasibility += value;
2629        }
2630    }
2631#ifdef JJF_ZERO
2632    if (largestInfeasibility > primalTolerance && numberBadRows*10 < numberRows) {
2633        // Can we improve by flipping
2634        for (int iPass = 0; iPass < 10; iPass++) {
2635            int numberColumns = solver->getNumCols();
2636            const CoinPackedMatrix * matrixByCol = solver->getMatrixByCol();
2637            const double * element = matrixByCol->getElements();
2638            const int * row = matrixByCol->getIndices();
2639            const CoinBigIndex * columnStart = matrixByCol->getVectorStarts();
2640            const int * columnLength = matrixByCol->getVectorLengths();
2641            double oldSum = sumInfeasibility;
2642            // First improve by moving continuous ones
2643            for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
2644                if (!solver->isInteger(iColumn)) {
2645                    double solValue = solution[iColumn];
2646                    double thetaUp = columnUpper[iColumn] - solValue;
2647                    double improvementUp = 0.0;
2648                    if (thetaUp > primalTolerance) {
2649                        // can go up
2650                        for (CoinBigIndex j = columnStart[iColumn];
2651                                j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2652                            int iRow = row[j];
2653                            double distanceUp = rowUpper[iRow] - rowActivity[iRow];
2654                            double distanceDown = rowLower[iRow] - rowActivity[iRow];
2655                            double el = element[j];
2656                            if (el > 0.0) {
2657                                // positive element
2658                                if (distanceUp > 0.0) {
2659                                    if (thetaUp*el > distanceUp)
2660                                        thetaUp = distanceUp / el;
2661                                } else {
2662                                    improvementUp -= el;
2663                                }
2664                                if (distanceDown > 0.0) {
2665                                    if (thetaUp*el > distanceDown)
2666                                        thetaUp = distanceDown / el;
2667                                    improvementUp += el;
2668                                }
2669                            } else {
2670                                // negative element
2671                                if (distanceDown < 0.0) {
2672                                    if (thetaUp*el < distanceDown)
2673                                        thetaUp = distanceDown / el;
2674                                } else {
2675                                    improvementUp += el;
2676                                }
2677                                if (distanceUp < 0.0) {
2678                                    if (thetaUp*el < distanceUp)
2679                                        thetaUp = distanceUp / el;
2680                                    improvementUp -= el;
2681                                }
2682                            }
2683                        }
2684                    }
2685                    double thetaDown = solValue - columnLower[iColumn];
2686                    double improvementDown = 0.0;
2687                    if (thetaDown > primalTolerance) {
2688                        // can go down
2689                        for (CoinBigIndex j = columnStart[iColumn];
2690                                j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2691                            int iRow = row[j];
2692                            double distanceUp = rowUpper[iRow] - rowActivity[iRow];
2693                            double distanceDown = rowLower[iRow] - rowActivity[iRow];
2694                            double el = -element[j]; // not change in sign form up
2695                            if (el > 0.0) {
2696                                // positive element
2697                                if (distanceUp > 0.0) {
2698                                    if (thetaDown*el > distanceUp)
2699                                        thetaDown = distanceUp / el;
2700                                } else {
2701                                    improvementDown -= el;
2702                                }
2703                                if (distanceDown > 0.0) {
2704                                    if (thetaDown*el > distanceDown)
2705                                        thetaDown = distanceDown / el;
2706                                    improvementDown += el;
2707                                }
2708                            } else {
2709                                // negative element
2710                                if (distanceDown < 0.0) {
2711                                    if (thetaDown*el < distanceDown)
2712                                        thetaDown = distanceDown / el;
2713                                } else {
2714                                    improvementDown += el;
2715                                }
2716                                if (distanceUp < 0.0) {
2717                                    if (thetaDown*el < distanceUp)
2718                                        thetaDown = distanceUp / el;
2719                                    improvementDown -= el;
2720                                }
2721                            }
2722                        }
2723                        if (thetaUp < 1.0e-8)
2724                            improvementUp = 0.0;
2725                        if (thetaDown < 1.0e-8)
2726                            improvementDown = 0.0;
2727                        double theta;
2728                        if (improvementUp >= improvementDown) {
2729                            theta = thetaUp;
2730                        } else {
2731                            improvementUp = improvementDown;
2732                            theta = -thetaDown;
2733                        }
2734                        if (improvementUp > 1.0e-8 && fabs(theta) > 1.0e-8) {
2735                            // Could move
2736                            double oldSum = 0.0;
2737                            double newSum = 0.0;
2738                            solution[iColumn] += theta;
2739                            for (CoinBigIndex j = columnStart[iColumn];
2740                                    j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2741                                int iRow = row[j];
2742                                double lower = rowLower[iRow];
2743                                double upper = rowUpper[iRow];
2744                                double value = rowActivity[iRow];
2745                                if (value > upper)
2746                                    oldSum += value - upper;
2747                                else if (value < lower)
2748                                    oldSum += lower - value;
2749                                value += theta * element[j];
2750                                rowActivity[iRow] = value;
2751                                if (value > upper)
2752                                    newSum += value - upper;
2753                                else if (value < lower)
2754                                    newSum += lower - value;
2755                            }
2756                            assert (newSum <= oldSum);
2757                            sumInfeasibility += newSum - oldSum;
2758                        }
2759                    }
2760                }
2761            }
2762            // Now flip some integers?
2763#ifdef JJF_ZERO
2764            for (i = 0; i < numberIntegers; i++) {
2765                int iColumn = integerVariable[i];
2766                double solValue = solution[iColumn];
2767                assert (fabs(solValue - floor(solValue + 0.5)) < 1.0e-8);
2768                double improvementUp = 0.0;
2769                if (columnUpper[iColumn] >= solValue + 1.0) {
2770                    // can go up
2771                    double oldSum = 0.0;
2772                    double newSum = 0.0;
2773                    for (CoinBigIndex j = columnStart[iColumn];
2774                            j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2775                        int iRow = row[j];
2776                        double lower = rowLower[iRow];
2777                        double upper = rowUpper[iRow];
2778                        double value = rowActivity[iRow];
2779                        if (value > upper)
2780                            oldSum += value - upper;
2781                        else if (value < lower)
2782                            oldSum += lower - value;
2783                        value += element[j];
2784                        if (value > upper)
2785                            newSum += value - upper;
2786                        else if (value < lower)
2787                            newSum += lower - value;
2788                    }
2789                    improvementUp = oldSum - newSum;
2790                }
2791                double improvementDown = 0.0;
2792                if (columnLower[iColumn] <= solValue - 1.0) {
2793                    // can go down
2794                    double oldSum = 0.0;
2795                    double newSum = 0.0;
2796                    for (CoinBigIndex j = columnStart[iColumn];
2797                            j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2798                        int iRow = row[j];
2799                        double lower = rowLower[iRow];
2800                        double upper = rowUpper[iRow];
2801                        double value = rowActivity[iRow];
2802                        if (value > upper)
2803                            oldSum += value - upper;
2804                        else if (value < lower)
2805                            oldSum += lower - value;
2806                        value -= element[j];
2807                        if (value > upper)
2808                            newSum += value - upper;
2809                        else if (value < lower)
2810                            newSum += lower - value;
2811                    }
2812                    improvementDown = oldSum - newSum;
2813                }
2814                double theta;
2815                if (improvementUp >= improvementDown) {
2816                    theta = 1.0;
2817                } else {
2818                    improvementUp = improvementDown;
2819                    theta = -1.0;
2820                }
2821                if (improvementUp > 1.0e-8 && fabs(theta) > 1.0e-8) {
2822                    // Could move
2823                    double oldSum = 0.0;
2824                    double newSum = 0.0;
2825                    solution[iColumn] += theta;
2826                    for (CoinBigIndex j = columnStart[iColumn];
2827                            j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2828                        int iRow = row[j];
2829                        double lower = rowLower[iRow];
2830                        double upper = rowUpper[iRow];
2831                        double value = rowActivity[iRow];
2832                        if (value > upper)
2833                            oldSum += value - upper;
2834                        else if (value < lower)
2835                            oldSum += lower - value;
2836                        value += theta * element[j];
2837                        rowActivity[iRow] = value;
2838                        if (value > upper)
2839                            newSum += value - upper;
2840                        else if (value < lower)
2841                            newSum += lower - value;
2842                    }
2843                    assert (newSum <= oldSum);
2844                    sumInfeasibility += newSum - oldSum;
2845                }
2846            }
2847#else
2848            int bestColumn = -1;
2849            double bestImprovement = primalTolerance;
2850            double theta = 0.0;
2851            for (i = 0; i < numberIntegers; i++) {
2852                int iColumn = integerVariable[i];
2853                double solValue = solution[iColumn];
2854                assert (fabs(solValue - floor(solValue + 0.5)) < 1.0e-8);
2855                double improvementUp = 0.0;
2856                if (columnUpper[iColumn] >= solValue + 1.0) {
2857                    // can go up
2858                    double oldSum = 0.0;
2859                    double newSum = 0.0;
2860                    for (CoinBigIndex j = columnStart[iColumn];
2861                            j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2862                        int iRow = row[j];
2863                        double lower = rowLower[iRow];
2864                        double upper = rowUpper[iRow];
2865                        double value = rowActivity[iRow];
2866                        if (value > upper)
2867                            oldSum += value - upper;
2868                        else if (value < lower)
2869                            oldSum += lower - value;
2870                        value += element[j];
2871                        if (value > upper)
2872                            newSum += value - upper;
2873                        else if (value < lower)
2874                            newSum += lower - value;
2875                    }
2876                    improvementUp = oldSum - newSum;
2877                }
2878                double improvementDown = 0.0;
2879                if (columnLower[iColumn] <= solValue - 1.0) {
2880                    // can go down
2881                    double oldSum = 0.0;
2882                    double newSum = 0.0;
2883                    for (CoinBigIndex j = columnStart[iColumn];
2884                            j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2885                        int iRow = row[j];
2886                        double lower = rowLower[iRow];
2887                        double upper = rowUpper[iRow];
2888                        double value = rowActivity[iRow];
2889                        if (value > upper)
2890                            oldSum += value - upper;
2891                        else if (value < lower)
2892                            oldSum += lower - value;
2893                        value -= element[j];
2894                        if (value > upper)
2895                            newSum += value - upper;
2896                        else if (value < lower)
2897                            newSum += lower - value;
2898                    }
2899                    improvementDown = oldSum - newSum;
2900                }
2901                double improvement = CoinMax(improvementUp, improvementDown);
2902                if (improvement > bestImprovement) {
2903                    bestImprovement = improvement;
2904                    bestColumn = iColumn;
2905                    if (improvementUp > improvementDown)
2906                        theta = 1.0;
2907                    else
2908                        theta = -1.0;
2909                }
2910            }
2911            if (bestColumn >= 0) {
2912                // Could move
2913                int iColumn = bestColumn;
2914                double oldSum = 0.0;
2915                double newSum = 0.0;
2916                solution[iColumn] += theta;
2917                for (CoinBigIndex j = columnStart[iColumn];
2918                        j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2919                    int iRow = row[j];
2920                    double lower = rowLower[iRow];
2921                    double upper = rowUpper[iRow];
2922                    double value = rowActivity[iRow];
2923                    if (value > upper)
2924                        oldSum += value - upper;
2925                    else if (value < lower)
2926                        oldSum += lower - value;
2927                    value += theta * element[j];
2928                    rowActivity[iRow] = value;
2929                    if (value > upper)
2930                        newSum += value - upper;
2931                    else if (value < lower)
2932                        newSum += lower - value;
2933                }
2934                assert (newSum <= oldSum);
2935                sumInfeasibility += newSum - oldSum;
2936            }
2937#endif
2938            if (oldSum <= sumInfeasibility + primalTolerance)
2939                break; // no good
2940        }
2941    }
2942    //delete [] saveSolution;
2943#endif
2944    delete [] rowActivity;
2945    return (largestInfeasibility > primalTolerance) ? 0 : 1;
2946}
2947// Set maximum Time (default off) - also sets starttime to current
2948void
2949CbcHeuristicFPump::setMaximumTime(double value)
2950{
2951    startTime_ = CoinCpuTime();
2952    maximumTime_ = value;
2953}
2954
2955# ifdef COIN_HAS_CLP
2956
2957//#############################################################################
2958// Constructors / Destructor / Assignment
2959//#############################################################################
2960
2961//-------------------------------------------------------------------
2962// Default Constructor
2963//-------------------------------------------------------------------
2964CbcDisasterHandler::CbcDisasterHandler (CbcModel * model)
2965        : OsiClpDisasterHandler(),
2966        cbcModel_(model)
2967{
2968    if (model) {
2969        osiModel_
2970        = dynamic_cast<OsiClpSolverInterface *> (model->solver());
2971        if (osiModel_)
2972            setSimplex(osiModel_->getModelPtr());
2973    }
2974}
2975
2976//-------------------------------------------------------------------
2977// Copy constructor
2978//-------------------------------------------------------------------
2979CbcDisasterHandler::CbcDisasterHandler (const CbcDisasterHandler & rhs)
2980        : OsiClpDisasterHandler(rhs),
2981        cbcModel_(rhs.cbcModel_)
2982{
2983}
2984
2985
2986//-------------------------------------------------------------------
2987// Destructor
2988//-------------------------------------------------------------------
2989CbcDisasterHandler::~CbcDisasterHandler ()
2990{
2991}
2992
2993//----------------------------------------------------------------
2994// Assignment operator
2995//-------------------------------------------------------------------
2996CbcDisasterHandler &
2997CbcDisasterHandler::operator=(const CbcDisasterHandler & rhs)
2998{
2999    if (this != &rhs) {
3000        OsiClpDisasterHandler::operator=(rhs);
3001        cbcModel_ = rhs.cbcModel_;
3002    }
3003    return *this;
3004}
3005//-------------------------------------------------------------------
3006// Clone
3007//-------------------------------------------------------------------
3008ClpDisasterHandler * CbcDisasterHandler::clone() const
3009{
3010    return new CbcDisasterHandler(*this);
3011}
3012// Type of disaster 0 can fix, 1 abort
3013int
3014CbcDisasterHandler::typeOfDisaster()
3015{
3016    if (!cbcModel_->parentModel() &&
3017            (cbcModel_->specialOptions()&2048) == 0) {
3018        return 0;
3019    } else {
3020        if (cbcModel_->parentModel())
3021            cbcModel_->setMaximumNodes(0);
3022        return 1;
3023    }
3024}
3025/* set model. */
3026void
3027CbcDisasterHandler::setCbcModel(CbcModel * model)
3028{
3029    cbcModel_ = model;
3030    if (model) {
3031        osiModel_
3032        = dynamic_cast<OsiClpSolverInterface *> (model->solver());
3033        if (osiModel_)
3034            setSimplex(osiModel_->getModelPtr());
3035        else
3036            setSimplex(NULL);
3037    }
3038}
3039#endif
3040
Note: See TracBrowser for help on using the repository browser.