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

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

adding a dubious heuristic

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 138.0 KB
Line 
1/* $Id: CbcHeuristicFPump.cpp 1945 2013-07-29 08:56:04Z 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                if (closestSolution && solver->getObjValue() < closestObjectiveValue) {
1655                    int i;
1656                    const double * objective = solver->getObjCoefficients();
1657                    for (i = 0; i < numberIntegersOrig; i++) {
1658                        int iColumn = integerVariableOrig[i];
1659                        if (objective[iColumn] > 0.0)
1660                            closestSolution[i] = 0;
1661                        else
1662                            closestSolution[i] = 1;
1663                    }
1664                    closestObjectiveValue = solver->getObjValue();
1665                }
1666                // See if we need to think about changing rhs
1667                if ((switches_&12) != 0 && useRhs < 1.0e50) {
1668                    double oldRhs = useRhs;
1669                    bool trying = false;
1670                    if ((switches_&4) != 0 && numberPasses && (numberPasses % 50) == 0) {
1671                        if (solutionValue > 1.0e20) {
1672                            // only if no genuine solution
1673                            double gap = useRhs - continuousObjectiveValue;
1674                            useRhs += 0.1 * gap;
1675                            if (exactMultiple) {
1676                                useRhs = exactMultiple * ceil(useRhs / exactMultiple);
1677                                useRhs = CoinMax(useRhs, oldRhs + exactMultiple);
1678                            }
1679                            trying = true;
1680                        }
1681                    }
1682                    if ((switches_&8) != 0) {
1683                        // Put in new suminf and check
1684                        double largest = newSumInfeas;
1685                        double smallest = newSumInfeas;
1686                        for (int i = 0; i < SIZE_BOBBLE - 1; i++) {
1687                            double value = saveSumInf[i+1];
1688                            saveSumInf[i] = value;
1689                            largest = CoinMax(largest, value);
1690                            smallest = CoinMin(smallest, value);
1691                        }
1692                        saveSumInf[SIZE_BOBBLE-1] = newSumInfeas;
1693                        if (smallest*1.5 > largest && smallest > 2.0) {
1694                            if (bobbleMode == 0) {
1695                                // go closer
1696                                double gap = oldRhs - continuousObjectiveValue;
1697                                useRhs -= 0.4 * gap;
1698                                if (exactMultiple) {
1699                                    double value = floor(useRhs / exactMultiple);
1700                                    useRhs = CoinMin(value * exactMultiple, oldRhs - exactMultiple);
1701                                }
1702                                if (useRhs < continuousObjectiveValue) {
1703                                    // skip decrease
1704                                    bobbleMode = 1;
1705                                    useRhs = oldRhs;
1706                                }
1707                            }
1708                            if (bobbleMode) {
1709                                trying = true;
1710                                // weaken
1711                                if (solutionValue < 1.0e20) {
1712                                    double gap = solutionValue - oldRhs;
1713                                    useRhs += 0.3 * gap;
1714                                } else {
1715                                    double gap = oldRhs - continuousObjectiveValue;
1716                                    useRhs += 0.05 * gap;
1717                                }
1718                                if (exactMultiple) {
1719                                    double value = ceil(useRhs / exactMultiple);
1720                                    useRhs = CoinMin(value * exactMultiple,
1721                                                     solutionValue - exactMultiple);
1722                                }
1723                            }
1724                            bobbleMode++;
1725                            // reset
1726                            CoinFillN(saveSumInf, SIZE_BOBBLE, COIN_DBL_MAX);
1727                        }
1728                    }
1729                    if (useRhs != oldRhs) {
1730                        // tidy up
1731                        if (exactMultiple) {
1732                            double value = floor(useRhs / exactMultiple);
1733                            double bestPossible = ceil(continuousObjectiveValue / exactMultiple);
1734                            useRhs = CoinMax(value, bestPossible) * exactMultiple;
1735                        } else {
1736                            useRhs = CoinMax(useRhs, continuousObjectiveValue);
1737                        }
1738                        int k = solver->getNumRows() - 1;
1739                        solver->setRowUpper(k, useRhs + useOffset);
1740                        bool takeHint;
1741                        OsiHintStrength strength;
1742                        solver->getHintParam(OsiDoDualInResolve, takeHint, strength);
1743                        if (useRhs < oldRhs) {
1744                            solver->setHintParam(OsiDoDualInResolve, true);
1745                            solver->resolve();
1746                        } else if (useRhs > oldRhs) {
1747                            solver->setHintParam(OsiDoDualInResolve, false);
1748                            solver->resolve();
1749                        }
1750                        solver->setHintParam(OsiDoDualInResolve, takeHint);
1751                        if (!solver->isProvenOptimal()) {
1752                            // presumably max time or some such
1753                            exitAll = true;
1754                            break;
1755                        }
1756                    } else if (trying) {
1757                        // doesn't look good
1758                        break;
1759                    }
1760                }
1761            }
1762            // reduce scale factor
1763            scaleFactor *= weightFactor_;
1764        } // END WHILE
1765        // see if rounding worked!
1766        if (roundingObjective < solutionValue) {
1767            if (roundingObjective < solutionValue - 1.0e-6*fabs(roundingObjective)) {
1768                sprintf(pumpPrint, "Rounding solution of %g is better than previous of %g\n",
1769                        roundingObjective, solutionValue);
1770                model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
1771                << pumpPrint
1772                << CoinMessageEol;
1773            }
1774            solutionValue = roundingObjective;
1775            newSolutionValue = solutionValue;
1776            memcpy(betterSolution, roundingSolution, numberColumns*sizeof(double));
1777            solutionFound = true;
1778            numberSolutions++;
1779            if (numberSolutions>=maxSolutions)
1780              exitAll = true;
1781            if (exitNow(roundingObjective))
1782                exitAll = true;
1783        }
1784        if (!solutionFound) {
1785            sprintf(pumpPrint, "No solution found this major pass");
1786            model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
1787            << pumpPrint
1788            << CoinMessageEol;
1789        }
1790        //}
1791        delete solver;
1792        solver = NULL;
1793        for ( j = 0; j < NUMBER_OLD; j++)
1794            delete [] oldSolution[j];
1795        delete [] oldSolution;
1796        delete [] saveObjective;
1797        if (usedColumn && !exitAll) {
1798            OsiSolverInterface * newSolver = cloneBut(3); // was model_->continuousSolver()->clone();
1799#if 0 //def COIN_HAS_CLP
1800            OsiClpSolverInterface * clpSolver
1801              = dynamic_cast<OsiClpSolverInterface *> (newSolver);
1802            if (clpSolver) {
1803              ClpSimplex * simplex = clpSolver->getModelPtr();
1804              simplex->writeMps("start.mps",2,1);
1805            }
1806#endif
1807            const double * colLower = newSolver->getColLower();
1808            const double * colUpper = newSolver->getColUpper();
1809            bool stopBAB = false;
1810            int allowedPass = -1;
1811            if (maximumAllowed > 0)
1812                allowedPass = CoinMax(numberPasses - maximumAllowed, -1);
1813            while (!stopBAB) {
1814                stopBAB = true;
1815                int i;
1816                int nFix = 0;
1817                int nFixI = 0;
1818                int nFixC = 0;
1819                int nFixC2 = 0;
1820                for (i = 0; i < numberIntegersOrig; i++) {
1821                    int iColumn = integerVariableOrig[i];
1822                    //const OsiObject * object = model_->object(i);
1823                    //double originalLower;
1824                    //double originalUpper;
1825                    //getIntegerInformation( object,originalLower, originalUpper);
1826                    //assert(colLower[iColumn]==originalLower);
1827                    //newSolver->setColLower(iColumn,CoinMax(colLower[iColumn],originalLower));
1828                    newSolver->setColLower(iColumn, colLower[iColumn]);
1829                    //assert(colUpper[iColumn]==originalUpper);
1830                    //newSolver->setColUpper(iColumn,CoinMin(colUpper[iColumn],originalUpper));
1831                    newSolver->setColUpper(iColumn, colUpper[iColumn]);
1832                    if (usedColumn[iColumn] <= allowedPass) {
1833                        double value = lastSolution[iColumn];
1834                        double nearest = floor(value + 0.5);
1835                        if (fabs(value - nearest) < 1.0e-7) {
1836                            if (nearest == colLower[iColumn]) {
1837                                newSolver->setColUpper(iColumn, colLower[iColumn]);
1838                                nFix++;
1839                            } else if (nearest == colUpper[iColumn]) {
1840                                newSolver->setColLower(iColumn, colUpper[iColumn]);
1841                                nFix++;
1842                            } else if (fixInternal) {
1843                                newSolver->setColLower(iColumn, nearest);
1844                                newSolver->setColUpper(iColumn, nearest);
1845                                nFix++;
1846                                nFixI++;
1847                            }
1848                        }
1849                    }
1850                }
1851                if (fixContinuous) {
1852                    for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
1853                        if (!newSolver->isInteger(iColumn) && usedColumn[iColumn] <= allowedPass) {
1854                            double value = lastSolution[iColumn];
1855                            if (value < colLower[iColumn] + 1.0e-8) {
1856                                newSolver->setColUpper(iColumn, colLower[iColumn]);
1857                                nFixC++;
1858                            } else if (value > colUpper[iColumn] - 1.0e-8) {
1859                                newSolver->setColLower(iColumn, colUpper[iColumn]);
1860                                nFixC++;
1861                            } else if (fixContinuous == 2) {
1862                                newSolver->setColLower(iColumn, value);
1863                                newSolver->setColUpper(iColumn, value);
1864                                nFixC++;
1865                                nFixC2++;
1866                            }
1867                        }
1868                    }
1869                }
1870                newSolver->initialSolve();
1871                if (!newSolver->isProvenOptimal()) {
1872                    //newSolver->writeMps("bad.mps");
1873                    //assert (newSolver->isProvenOptimal());
1874                    exitAll = true;
1875                    break;
1876                }
1877                sprintf(pumpPrint, "Before mini branch and bound, %d integers at bound fixed and %d continuous",
1878                        nFix, nFixC);
1879                if (nFixC2 + nFixI != 0)
1880                    sprintf(pumpPrint + strlen(pumpPrint), " of which %d were internal integer and %d internal continuous",
1881                            nFixI, nFixC2);
1882                model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
1883                << pumpPrint
1884                << CoinMessageEol;
1885                double saveValue = newSolutionValue;
1886                if (newSolutionValue - model_->getCutoffIncrement()
1887                        > continuousObjectiveValue - 1.0e-7) {
1888                    double saveFraction = fractionSmall_;
1889                    if (numberTries > 1 && !numberBandBsolutions)
1890                        fractionSmall_ *= 0.5;
1891                    // Give branch and bound a bit more freedom
1892                    double cutoff2 = newSolutionValue +
1893                                     CoinMax(model_->getCutoffIncrement(), 1.0e-3);
1894#if 0
1895                      {
1896                        OsiClpSolverInterface * clpSolver
1897                        = dynamic_cast<OsiClpSolverInterface *> (newSolver);
1898                        if (clpSolver) {
1899                            ClpSimplex * simplex = clpSolver->getModelPtr();
1900                            simplex->writeMps("testA.mps",2,1);
1901                        }
1902                      }
1903#endif
1904                    int returnCode2 = smallBranchAndBound(newSolver, numberNodes_, newSolution, newSolutionValue,
1905                                                          cutoff2, "CbcHeuristicLocalAfterFPump");
1906                    fractionSmall_ = saveFraction;
1907                    if (returnCode2 < 0) {
1908                        if (returnCode2 == -2) {
1909                            exitAll = true;
1910                            returnCode = 0;
1911                        } else {
1912                            returnCode2 = 0; // returned on size - try changing
1913                            //#define ROUND_AGAIN
1914#ifdef ROUND_AGAIN
1915                            if (numberTries == 1 && numberPasses > 20 && allowedPass < numberPasses - 1) {
1916                                allowedPass = (numberPasses + allowedPass) >> 1;
1917                                sprintf(pumpPrint,
1918                                        "Fixing all variables which were last changed on pass %d and trying again",
1919                                        allowedPass);
1920                                model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
1921                                << pumpPrint
1922                                << CoinMessageEol;
1923                                stopBAB = false;
1924                                continue;
1925                            }
1926#endif
1927                        }
1928                    }
1929                    if ((returnCode2&2) != 0) {
1930                        // could add cut
1931                        returnCode2 &= ~2;
1932                    }
1933                    if (returnCode2) {
1934                        numberBandBsolutions++;
1935                        // may not have got solution earlier
1936                        returnCode |= 1;
1937                    }
1938                } else {
1939                    // no need
1940                    exitAll = true;
1941                    //returnCode=0;
1942                }
1943                // recompute solution value
1944                if (returnCode && true) {
1945#if 0
1946                      {
1947                        OsiClpSolverInterface * clpSolver
1948                        = dynamic_cast<OsiClpSolverInterface *> (newSolver);
1949                        if (clpSolver) {
1950                            ClpSimplex * simplex = clpSolver->getModelPtr();
1951                            simplex->writeMps("testB.mps",2,1);
1952                        }
1953                      }
1954#endif
1955                    delete newSolver;
1956                    newSolver = cloneBut(3); // was model_->continuousSolver()->clone();
1957                    newSolutionValue = -saveOffset;
1958                    double newSumInfeas = 0.0;
1959                    const double * obj = newSolver->getObjCoefficients();
1960                    for (int i = 0 ; i < numberColumns ; i++ ) {
1961                        if (newSolver->isInteger(i)) {
1962                            double value = newSolution[i];
1963                            double nearest = floor(value + 0.5);
1964                            newSumInfeas += fabs(value - nearest);
1965                        }
1966                        newSolutionValue += obj[i] * newSolution[i];
1967                    }
1968                    newSolutionValue *= direction;
1969                }
1970                bool gotSolution = false;
1971                if (returnCode && newSolutionValue < saveValue) {
1972                    sprintf(pumpPrint, "Mini branch and bound improved solution from %g to %g (%.2f seconds)",
1973                            saveValue, newSolutionValue, model_->getCurrentSeconds());
1974                    model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
1975                    << pumpPrint
1976                    << CoinMessageEol;
1977                    memcpy(betterSolution, newSolution, numberColumns*sizeof(double));
1978                    gotSolution = true;
1979                    if (fixContinuous && nFixC + nFixC2 > 0) {
1980                        // may be able to do even better
1981                        int nFixed = 0;
1982                        const double * lower = model_->solver()->getColLower();
1983                        const double * upper = model_->solver()->getColUpper();
1984                        for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
1985                            double value = newSolution[iColumn];
1986                            if (newSolver->isInteger(iColumn)) {
1987                                value = floor(newSolution[iColumn] + 0.5);
1988                                newSolver->setColLower(iColumn, value);
1989                                newSolver->setColUpper(iColumn, value);
1990                                nFixed++;
1991                            } else {
1992                                newSolver->setColLower(iColumn, lower[iColumn]);
1993                                newSolver->setColUpper(iColumn, upper[iColumn]);
1994                                if (value < lower[iColumn])
1995                                    value = lower[iColumn];
1996                                else if (value > upper[iColumn])
1997                                    value = upper[iColumn];
1998
1999                            }
2000                            newSolution[iColumn] = value;
2001                        }
2002                        newSolver->setColSolution(newSolution);
2003                        //#define CLP_INVESTIGATE2
2004#ifdef CLP_INVESTIGATE2
2005                        {
2006                            // check
2007                            // get row activities
2008                            int numberRows = newSolver->getNumRows();
2009                            double * rowActivity = new double[numberRows];
2010                            memset(rowActivity, 0, numberRows*sizeof(double));
2011                            newSolver->getMatrixByCol()->times(newSolution, rowActivity) ;
2012                            double largestInfeasibility = primalTolerance;
2013                            double sumInfeasibility = 0.0;
2014                            int numberBadRows = 0;
2015                            const double * rowLower = newSolver->getRowLower();
2016                            const double * rowUpper = newSolver->getRowUpper();
2017                            for (i = 0 ; i < numberRows ; i++) {
2018                                double value;
2019                                value = rowLower[i] - rowActivity[i];
2020                                if (value > primalTolerance) {
2021                                    numberBadRows++;
2022                                    largestInfeasibility = CoinMax(largestInfeasibility, value);
2023                                    sumInfeasibility += value;
2024                                }
2025                                value = rowActivity[i] - rowUpper[i];
2026                                if (value > primalTolerance) {
2027                                    numberBadRows++;
2028                                    largestInfeasibility = CoinMax(largestInfeasibility, value);
2029                                    sumInfeasibility += value;
2030                                }
2031                            }
2032                            printf("%d bad rows, largest inf %g sum %g\n",
2033                                   numberBadRows, largestInfeasibility, sumInfeasibility);
2034                            delete [] rowActivity;
2035                        }
2036#endif
2037#ifdef COIN_HAS_CLP
2038                        OsiClpSolverInterface * clpSolver
2039                        = dynamic_cast<OsiClpSolverInterface *> (newSolver);
2040                        if (clpSolver) {
2041                            ClpSimplex * simplex = clpSolver->getModelPtr();
2042                            //simplex->writeBasis("test.bas",true,2);
2043                            //simplex->writeMps("test.mps",2,1);
2044                            if (nFixed*3 > numberColumns*2)
2045                                simplex->allSlackBasis(); // may as well go from all slack
2046                            simplex->primal(1);
2047                            clpSolver->setWarmStart(NULL);
2048                        }
2049#endif
2050                        newSolver->initialSolve();
2051                        if (newSolver->isProvenOptimal()) {
2052                            double value = newSolver->getObjValue() * newSolver->getObjSense();
2053                            if (value < newSolutionValue) {
2054                                //newSolver->writeMps("query","mps");
2055#ifdef JJF_ZERO
2056                                {
2057                                    double saveOffset;
2058                                    newSolver->getDblParam(OsiObjOffset, saveOffset);
2059                                    const double * obj = newSolver->getObjCoefficients();
2060                                    double newTrueSolutionValue = -saveOffset;
2061                                    double newSumInfeas = 0.0;
2062                                    int numberColumns = newSolver->getNumCols();
2063                                    const double * solution = newSolver->getColSolution();
2064                                    for (int  i = 0 ; i < numberColumns ; i++ ) {
2065                                        if (newSolver->isInteger(i)) {
2066                                            double value = solution[i];
2067                                            double nearest = floor(value + 0.5);
2068                                            newSumInfeas += fabs(value - nearest);
2069                                        }
2070                                        if (solution[i])
2071                                            printf("%d obj %g val %g - total %g\n", i, obj[i], solution[i],
2072                                                   newTrueSolutionValue);
2073                                        newTrueSolutionValue += obj[i] * solution[i];
2074                                    }
2075                                    printf("obj %g - inf %g\n", newTrueSolutionValue,
2076                                           newSumInfeas);
2077                                }
2078#endif
2079                                sprintf(pumpPrint, "Freeing continuous variables gives a solution of %g", value);
2080                                model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
2081                                << pumpPrint
2082                                << CoinMessageEol;
2083                                newSolutionValue = value;
2084                                memcpy(betterSolution, newSolver->getColSolution(), numberColumns*sizeof(double));
2085                            }
2086                        } else {
2087                            //newSolver->writeMps("bad3.mps");
2088                          sprintf(pumpPrint, "On closer inspection solution is not valid");
2089                          model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
2090                            << pumpPrint
2091                            << CoinMessageEol;
2092                            exitAll = true;
2093                            break;
2094                        }
2095                    }
2096                } else {
2097                    sprintf(pumpPrint, "Mini branch and bound did not improve solution (%.2f seconds)",
2098                            model_->getCurrentSeconds());
2099                    model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
2100                    << pumpPrint
2101                    << CoinMessageEol;
2102                    if (returnCode && newSolutionValue < saveValue + 1.0e-3 && nFixC + nFixC2) {
2103                        // may be able to do better
2104                        const double * lower = model_->solver()->getColLower();
2105                        const double * upper = model_->solver()->getColUpper();
2106                        for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
2107                            if (newSolver->isInteger(iColumn)) {
2108                                double value = floor(newSolution[iColumn] + 0.5);
2109                                newSolver->setColLower(iColumn, value);
2110                                newSolver->setColUpper(iColumn, value);
2111                            } else {
2112                                newSolver->setColLower(iColumn, lower[iColumn]);
2113                                newSolver->setColUpper(iColumn, upper[iColumn]);
2114                            }
2115                        }
2116                        newSolver->initialSolve();
2117                        if (newSolver->isProvenOptimal()) {
2118                            double value = newSolver->getObjValue() * newSolver->getObjSense();
2119                            if (value < saveValue) {
2120                                sprintf(pumpPrint, "Freeing continuous variables gives a solution of %g", value);
2121                                model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
2122                                << pumpPrint
2123                                << CoinMessageEol;
2124                                newSolutionValue = value;
2125                                memcpy(betterSolution, newSolver->getColSolution(), numberColumns*sizeof(double));
2126                                gotSolution = true;
2127                            }
2128                        }
2129                    }
2130                }
2131                if (gotSolution) {
2132                    if ((accumulate_&1) != 0) {
2133                        model_->incrementUsed(betterSolution); // for local search
2134                    }
2135                    solutionValue = newSolutionValue;
2136                    solutionFound = true;
2137                    numberSolutions++;
2138                    if (numberSolutions>=maxSolutions)
2139                      exitAll = true;
2140                    if (exitNow(newSolutionValue))
2141                        exitAll = true;
2142                    CoinWarmStartBasis * basis =
2143                        dynamic_cast<CoinWarmStartBasis *>(newSolver->getWarmStart()) ;
2144                    if (basis) {
2145                        bestBasis = * basis;
2146                        delete basis;
2147                        int action = model_->dealWithEventHandler(CbcEventHandler::heuristicSolution, newSolutionValue, betterSolution);
2148                        if (action == 0) {
2149                            double * saveOldSolution = CoinCopyOfArray(model_->bestSolution(), numberColumns);
2150                            double saveObjectiveValue = model_->getMinimizationObjValue();
2151                            model_->setBestSolution(betterSolution, numberColumns, newSolutionValue);
2152                            if (saveOldSolution && saveObjectiveValue < model_->getMinimizationObjValue())
2153                                model_->setBestSolution(saveOldSolution, numberColumns, saveObjectiveValue);
2154                            delete [] saveOldSolution;
2155                        }
2156                        if (!action || model_->getCurrentSeconds() > model_->getMaximumSeconds()) {
2157                            exitAll = true; // exit
2158                            break;
2159                        }
2160                    }
2161                }
2162            } // end stopBAB while
2163            delete newSolver;
2164        }
2165        if (solutionFound) finalReturnCode = 1;
2166        cutoff = CoinMin(cutoff, solutionValue - model_->getCutoffIncrement());
2167        if (numberTries >= maximumRetries_ || !solutionFound || exitAll || cutoff < continuousObjectiveValue + 1.0e-7) {
2168            break;
2169        } else {
2170            solutionFound = false;
2171            if (absoluteIncrement_ > 0.0 || relativeIncrement_ > 0.0) {
2172                double gap = relativeIncrement_ * fabs(solutionValue);
2173                double change = CoinMax(gap, absoluteIncrement_);
2174                cutoff = CoinMin(cutoff, solutionValue - change);
2175            } else {
2176                //double weights[10]={0.1,0.1,0.2,0.2,0.2,0.3,0.3,0.3,0.4,0.5};
2177                double weights[10] = {0.1, 0.2, 0.3, 0.3, 0.4, 0.4, 0.4, 0.5, 0.5, 0.6};
2178                cutoff -= weights[CoinMin(numberTries-1, 9)] * (cutoff - continuousObjectiveValue);
2179            }
2180            // But round down
2181            if (exactMultiple)
2182                cutoff = exactMultiple * floor(cutoff / exactMultiple);
2183            if (cutoff < continuousObjectiveValue)
2184                break;
2185            sprintf(pumpPrint, "Round again with cutoff of %g", cutoff);
2186            model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
2187            << pumpPrint
2188            << CoinMessageEol;
2189            if ((accumulate_&3) < 2 && usedColumn) {
2190                for (int i = 0; i < numberColumns; i++)
2191                    usedColumn[i] = -1;
2192            }
2193            averageIterationsPerTry = totalNumberIterations / numberTries;
2194            numberIterationsLastPass =  totalNumberIterations;
2195            totalNumberPasses += numberPasses - 1;
2196        }
2197    }
2198/*
2199  End of the `exitAll' loop.
2200*/
2201#ifdef RAND_RAND
2202    delete [] randomFactor;
2203#endif
2204    delete solver; // probably NULL but do anyway
2205    if (!finalReturnCode && closestSolution && closestObjectiveValue <= 10.0 &&
2206            usedColumn && !model_->maximumSecondsReached()) {
2207        // try a bit of branch and bound
2208        OsiSolverInterface * newSolver = cloneBut(1); // was model_->continuousSolver()->clone();
2209        const double * colLower = newSolver->getColLower();
2210        const double * colUpper = newSolver->getColUpper();
2211        int i;
2212        double rhs = 0.0;
2213        for (i = 0; i < numberIntegersOrig; i++) {
2214            int iColumn = integerVariableOrig[i];
2215            int direction = closestSolution[i];
2216            closestSolution[i] = iColumn;
2217            if (direction == 0) {
2218                // keep close to LB
2219                rhs += colLower[iColumn];
2220                lastSolution[i] = 1.0;
2221            } else {
2222                // keep close to UB
2223                rhs -= colUpper[iColumn];
2224                lastSolution[i] = -1.0;
2225            }
2226        }
2227        newSolver->addRow(numberIntegersOrig, closestSolution,
2228                          lastSolution, -COIN_DBL_MAX, rhs + 10.0);
2229        //double saveValue = newSolutionValue;
2230        //newSolver->writeMps("sub");
2231        int returnCode = smallBranchAndBound(newSolver, numberNodes_, newSolution, newSolutionValue,
2232                                             newSolutionValue, "CbcHeuristicLocalAfterFPump");
2233        if (returnCode < 0)
2234            returnCode = 0; // returned on size
2235        if ((returnCode&2) != 0) {
2236            // could add cut
2237            returnCode &= ~2;
2238        }
2239        if (returnCode) {
2240            //printf("old sol of %g new of %g\n",saveValue,newSolutionValue);
2241            memcpy(betterSolution, newSolution, numberColumns*sizeof(double));
2242            //abort();
2243            solutionValue = newSolutionValue;
2244            solutionFound = true;
2245            numberSolutions++;
2246            if (numberSolutions>=maxSolutions)
2247              exitAll = true;
2248            if (exitNow(newSolutionValue))
2249                exitAll = true;
2250        }
2251        delete newSolver;
2252    }
2253    delete clonedSolver;
2254    delete [] roundingSolution;
2255    delete [] usedColumn;
2256    delete [] lastSolution;
2257    delete [] newSolution;
2258    delete [] closestSolution;
2259    delete [] integerVariable;
2260    delete [] firstPerturbedObjective;
2261    delete [] firstPerturbedSolution;
2262    if (solutionValue == incomingObjective)
2263        sprintf(pumpPrint, "After %.2f seconds - Feasibility pump exiting - took %.2f seconds",
2264                model_->getCurrentSeconds(), CoinCpuTime() - time1);
2265    else if (numberSolutions < maxSolutions)
2266      sprintf(pumpPrint, "After %.2f seconds - Feasibility pump exiting with objective of %g - took %.2f seconds",
2267              model_->getCurrentSeconds(), solutionValue, CoinCpuTime() - time1);
2268    else 
2269        sprintf(pumpPrint, "After %.2f seconds - Feasibility pump exiting with objective of %g (stopping after %d solutions) - took %.2f seconds",
2270                model_->getCurrentSeconds(), solutionValue, 
2271                numberSolutions,CoinCpuTime() - time1);
2272    model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
2273      << pumpPrint
2274      << CoinMessageEol;
2275    if (bestBasis.getNumStructural())
2276        model_->setBestSolutionBasis(bestBasis);
2277    //model_->setMinimizationObjValue(saveBestObjective);
2278    if ((accumulate_&1) != 0 && numberSolutions > 1 && !model_->getSolutionCount()) {
2279        model_->setSolutionCount(1); // for local search
2280        model_->setNumberHeuristicSolutions(1);
2281    }
2282#ifdef COIN_DEVELOP
2283    {
2284        double ncol = model_->solver()->getNumCols();
2285        double nrow = model_->solver()->getNumRows();
2286        printf("XXX total iterations %d ratios - %g %g %g\n",
2287               totalNumberIterations,
2288               static_cast<double> (totalNumberIterations) / nrow,
2289               static_cast<double> (totalNumberIterations) / ncol,
2290               static_cast<double> (totalNumberIterations) / (2*nrow + 2*ncol));
2291    }
2292#endif
2293    return finalReturnCode;
2294}
2295
2296/**************************END MAIN PROCEDURE ***********************************/
2297
2298// update model
2299void CbcHeuristicFPump::setModel(CbcModel * model)
2300{
2301    model_ = model;
2302}
2303
2304/* Rounds solution - down if < downValue
2305   returns 1 if current is a feasible solution
2306*/
2307int
2308CbcHeuristicFPump::rounds(OsiSolverInterface * solver, double * solution,
2309                          //const double * objective,
2310                          int numberIntegers, const int * integerVariable,
2311                          /*char * pumpPrint,*/ int iter,
2312                          /*bool roundExpensive,*/ double downValue, int *flip)
2313{
2314    double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
2315    double primalTolerance ;
2316    solver->getDblParam(OsiPrimalTolerance, primalTolerance) ;
2317
2318    int i;
2319
2320    const double * cost = solver->getObjCoefficients();
2321    int flip_up = 0;
2322    int flip_down  = 0;
2323    double  v = randomNumberGenerator_.randomDouble() * 20.0;
2324    int nn = 10 + static_cast<int> (v);
2325    int nnv = 0;
2326    int * list = new int [nn];
2327    double * val = new double [nn];
2328    for (i = 0; i < nn; i++) val[i] = .001;
2329
2330    const double * rowLower = solver->getRowLower();
2331    const double * rowUpper = solver->getRowUpper();
2332    int numberRows = solver->getNumRows();
2333    if (false && (iter&1) != 0) {
2334        // Do set covering variables
2335        const CoinPackedMatrix * matrixByRow = solver->getMatrixByRow();
2336        const double * elementByRow = matrixByRow->getElements();
2337        const int * column = matrixByRow->getIndices();
2338        const CoinBigIndex * rowStart = matrixByRow->getVectorStarts();
2339        const int * rowLength = matrixByRow->getVectorLengths();
2340        for (i = 0; i < numberRows; i++) {
2341            if (rowLower[i] == 1.0 && rowUpper[i] == 1.0) {
2342                bool cover = true;
2343                double largest = 0.0;
2344                int jColumn = -1;
2345                for (CoinBigIndex k = rowStart[i]; k < rowStart[i] + rowLength[i]; k++) {
2346                    int iColumn = column[k];
2347                    if (elementByRow[k] != 1.0 || !solver->isInteger(iColumn)) {
2348                        cover = false;
2349                        break;
2350                    } else {
2351                        if (solution[iColumn]) {
2352                            double value = solution[iColumn] *
2353                                           (randomNumberGenerator_.randomDouble() + 5.0);
2354                            if (value > largest) {
2355                                largest = value;
2356                                jColumn = iColumn;
2357                            }
2358                        }
2359                    }
2360                }
2361                if (cover) {
2362                    for (CoinBigIndex k = rowStart[i]; k < rowStart[i] + rowLength[i]; k++) {
2363                        int iColumn = column[k];
2364                        if (iColumn == jColumn)
2365                            solution[iColumn] = 1.0;
2366                        else
2367                            solution[iColumn] = 0.0;
2368                    }
2369                }
2370            }
2371        }
2372    }
2373    int numberColumns = solver->getNumCols();
2374#ifdef JJF_ZERO
2375    // Do set covering variables
2376    const CoinPackedMatrix * matrixByRow = solver->getMatrixByRow();
2377    const double * elementByRow = matrixByRow->getElements();
2378    const int * column = matrixByRow->getIndices();
2379    const CoinBigIndex * rowStart = matrixByRow->getVectorStarts();
2380    const int * rowLength = matrixByRow->getVectorLengths();
2381    double * sortTemp = new double[numberColumns];
2382    int * whichTemp = new int [numberColumns];
2383    char * rowTemp = new char [numberRows];
2384    memset(rowTemp, 0, numberRows);
2385    for (i = 0; i < numberColumns; i++)
2386        whichTemp[i] = -1;
2387    int nSOS = 0;
2388    for (i = 0; i < numberRows; i++) {
2389        if (rowLower[i] == 1.0 && rowUpper[i] == 1.0) {
2390            bool cover = true;
2391            for (CoinBigIndex k = rowStart[i]; k < rowStart[i] + rowLength[i]; k++) {
2392                int iColumn = column[k];
2393                if (elementByRow[k] != 1.0 || !solver->isInteger(iColumn)) {
2394                    cover = false;
2395                    break;
2396                }
2397            }
2398            if (cover) {
2399                rowTemp[i] = 1;
2400                nSOS++;
2401                for (CoinBigIndex k = rowStart[i]; k < rowStart[i] + rowLength[i]; k++) {
2402                    int iColumn = column[k];
2403                    double value = solution[iColumn];
2404                    whichTemp[iColumn] = iColumn;
2405                }
2406            }
2407        }
2408    }
2409    if (nSOS) {
2410        // Column copy
2411        const CoinPackedMatrix * matrixByColumn = solver->getMatrixByCol();
2412        //const double * element = matrixByColumn->getElements();
2413        const int * row = matrixByColumn->getIndices();
2414        const CoinBigIndex * columnStart = matrixByColumn->getVectorStarts();
2415        const int * columnLength = matrixByColumn->getVectorLengths();
2416        int nLook = 0;
2417        for (i = 0; i < numberColumns; i++) {
2418            if (whichTemp[i] >= 0) {
2419                whichTemp[nLook] = i;
2420                double value = solution[i];
2421                if (value < 0.5)
2422                    value *= (0.1 * randomNumberGenerator_.randomDouble() + 0.3);
2423                sortTemp[nLook++] = -value;
2424            }
2425        }
2426        CoinSort_2(sortTemp, sortTemp + nLook, whichTemp);
2427        double smallest = 1.0;
2428        int nFix = 0;
2429        int nOne = 0;
2430        for (int j = 0; j < nLook; j++) {
2431            int jColumn = whichTemp[j];
2432            double thisValue = solution[jColumn];
2433            if (!thisValue)
2434                continue;
2435            if (thisValue == 1.0)
2436                nOne++;
2437            smallest = CoinMin(smallest, thisValue);
2438            solution[jColumn] = 1.0;
2439            double largest = 0.0;
2440            for (CoinBigIndex jEl = columnStart[jColumn];
2441                    jEl < columnStart[jColumn] + columnLength[jColumn]; jEl++) {
2442                int jRow = row[jEl];
2443                if (rowTemp[jRow]) {
2444                    for (CoinBigIndex k = rowStart[jRow]; k < rowStart[jRow] + rowLength[jRow]; k++) {
2445                        int iColumn = column[k];
2446                        if (solution[iColumn]) {
2447                            if (iColumn != jColumn) {
2448                                double value = solution[iColumn];
2449                                if (value > largest)
2450                                    largest = value;
2451                                solution[iColumn] = 0.0;
2452                            }
2453                        }
2454                    }
2455                }
2456            }
2457            if (largest > thisValue)
2458                printf("%d was at %g - chosen over a value of %g\n",
2459                       jColumn, thisValue, largest);
2460            nFix++;
2461        }
2462        printf("%d fixed out of %d (%d at one already)\n",
2463               nFix, nLook, nOne);
2464    }
2465    delete [] sortTemp;
2466    delete [] whichTemp;
2467    delete [] rowTemp;
2468#endif
2469    const double * columnLower = solver->getColLower();
2470    const double * columnUpper = solver->getColUpper();
2471    // Check if valid with current solution (allow for 0.99999999s)
2472    double newSumInfeas = 0.0;
2473    int newNumberInfeas = 0;
2474    for (i = 0; i < numberIntegers; i++) {
2475        int iColumn = integerVariable[i];
2476        double value = solution[iColumn];
2477        double round = floor(value + 0.5);
2478        if (fabs(value - round) > primalTolerance) {
2479          newSumInfeas += fabs(value-round);
2480          newNumberInfeas++;
2481        }
2482    }
2483    if (!newNumberInfeas) {
2484        // may be able to use solution even if 0.99999's
2485        double * saveLower = CoinCopyOfArray(columnLower, numberColumns);
2486        double * saveUpper = CoinCopyOfArray(columnUpper, numberColumns);
2487        double * saveSolution = CoinCopyOfArray(solution, numberColumns);
2488        double * tempSolution = CoinCopyOfArray(solution, numberColumns);
2489        CoinWarmStartBasis  * saveBasis =
2490            dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ;
2491        for (i = 0; i < numberIntegers; i++) {
2492            int iColumn = integerVariable[i];
2493            double value = solution[iColumn];
2494            double round = floor(value + 0.5);
2495            solver->setColLower(iColumn, round);
2496            solver->setColUpper(iColumn, round);
2497            tempSolution[iColumn] = round;
2498        }
2499        solver->setColSolution(tempSolution);
2500        delete [] tempSolution;
2501        solver->resolve();
2502        solver->setColLower(saveLower);
2503        solver->setColUpper(saveUpper);
2504        solver->setWarmStart(saveBasis);
2505        delete [] saveLower;
2506        delete [] saveUpper;
2507        delete saveBasis;
2508        if (!solver->isProvenOptimal()) {
2509            solver->setColSolution(saveSolution);
2510        }
2511        delete [] saveSolution;
2512        if (solver->isProvenOptimal()) {
2513            // feasible
2514            delete [] list;
2515            delete [] val;
2516            return 1;
2517        }
2518    }
2519    //double * saveSolution = CoinCopyOfArray(solution,numberColumns);
2520    // return rounded solution
2521    for (i = 0; i < numberIntegers; i++) {
2522        int iColumn = integerVariable[i];
2523        double value = solution[iColumn];
2524        double round = floor(value + primalTolerance);
2525        if (value - round > downValue) round += 1.;
2526#ifndef JJF_ONE
2527        if (round < integerTolerance && cost[iColumn] < -1. + integerTolerance) flip_down++;
2528        if (round > 1. - integerTolerance && cost[iColumn] > 1. - integerTolerance) flip_up++;
2529#else
2530        if (round < columnLower[iColumn] + integerTolerance && cost[iColumn] < -1. + integerTolerance) flip_down++;
2531        if (round > columnUpper[iColumn] - integerTolerance && cost[iColumn] > 1. - integerTolerance) flip_up++;
2532#endif
2533        if (flip_up + flip_down == 0) {
2534            for (int k = 0; k < nn; k++) {
2535                if (fabs(value - round) > val[k]) {
2536                    nnv++;
2537                    for (int j = nn - 2; j >= k; j--) {
2538                        val[j+1] = val[j];
2539                        list[j+1] = list[j];
2540                    }
2541                    val[k] = fabs(value - round);
2542                    list[k] = iColumn;
2543                    break;
2544                }
2545            }
2546        }
2547        solution[iColumn] = round;
2548    }
2549
2550    if (nnv > nn) nnv = nn;
2551    //if (iter != 0)
2552    //sprintf(pumpPrint+strlen(pumpPrint),"up = %5d , down = %5d", flip_up, flip_down);
2553    *flip = flip_up + flip_down;
2554
2555    if (*flip == 0 && iter != 0) {
2556        //sprintf(pumpPrint+strlen(pumpPrint)," -- rand = %4d (%4d) ", nnv, nn);
2557        for (i = 0; i < nnv; i++) {
2558            // was solution[list[i]] = 1. - solution[list[i]]; but does that work for 7>=x>=6
2559            int index = list[i];
2560            double value = solution[index];
2561            if (value <= 1.0) {
2562                solution[index] = 1.0 - value;
2563            } else if (value < columnLower[index] + integerTolerance) {
2564                solution[index] = value + 1.0;
2565            } else if (value > columnUpper[index] - integerTolerance) {
2566                solution[index] = value - 1.0;
2567            } else {
2568                solution[index] = value - 1.0;
2569            }
2570        }
2571        *flip = nnv;
2572    } else {
2573        //sprintf(pumpPrint+strlen(pumpPrint)," ");
2574    }
2575    delete [] list;
2576    delete [] val;
2577    //iter++;
2578
2579    // get row activities
2580    double * rowActivity = new double[numberRows];
2581    memset(rowActivity, 0, numberRows*sizeof(double));
2582    solver->getMatrixByCol()->times(solution, rowActivity) ;
2583    double largestInfeasibility = primalTolerance;
2584    double sumInfeasibility = 0.0;
2585    int numberBadRows = 0;
2586    for (i = 0 ; i < numberRows ; i++) {
2587        double value;
2588        value = rowLower[i] - rowActivity[i];
2589        if (value > primalTolerance) {
2590            numberBadRows++;
2591            largestInfeasibility = CoinMax(largestInfeasibility, value);
2592            sumInfeasibility += value;
2593        }
2594        value = rowActivity[i] - rowUpper[i];
2595        if (value > primalTolerance) {
2596            numberBadRows++;
2597            largestInfeasibility = CoinMax(largestInfeasibility, value);
2598            sumInfeasibility += value;
2599        }
2600    }
2601#ifdef JJF_ZERO
2602    if (largestInfeasibility > primalTolerance && numberBadRows*10 < numberRows) {
2603        // Can we improve by flipping
2604        for (int iPass = 0; iPass < 10; iPass++) {
2605            int numberColumns = solver->getNumCols();
2606            const CoinPackedMatrix * matrixByCol = solver->getMatrixByCol();
2607            const double * element = matrixByCol->getElements();
2608            const int * row = matrixByCol->getIndices();
2609            const CoinBigIndex * columnStart = matrixByCol->getVectorStarts();
2610            const int * columnLength = matrixByCol->getVectorLengths();
2611            double oldSum = sumInfeasibility;
2612            // First improve by moving continuous ones
2613            for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
2614                if (!solver->isInteger(iColumn)) {
2615                    double solValue = solution[iColumn];
2616                    double thetaUp = columnUpper[iColumn] - solValue;
2617                    double improvementUp = 0.0;
2618                    if (thetaUp > primalTolerance) {
2619                        // can go up
2620                        for (CoinBigIndex j = columnStart[iColumn];
2621                                j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2622                            int iRow = row[j];
2623                            double distanceUp = rowUpper[iRow] - rowActivity[iRow];
2624                            double distanceDown = rowLower[iRow] - rowActivity[iRow];
2625                            double el = element[j];
2626                            if (el > 0.0) {
2627                                // positive element
2628                                if (distanceUp > 0.0) {
2629                                    if (thetaUp*el > distanceUp)
2630                                        thetaUp = distanceUp / el;
2631                                } else {
2632                                    improvementUp -= el;
2633                                }
2634                                if (distanceDown > 0.0) {
2635                                    if (thetaUp*el > distanceDown)
2636                                        thetaUp = distanceDown / el;
2637                                    improvementUp += el;
2638                                }
2639                            } else {
2640                                // negative element
2641                                if (distanceDown < 0.0) {
2642                                    if (thetaUp*el < distanceDown)
2643                                        thetaUp = distanceDown / el;
2644                                } else {
2645                                    improvementUp += el;
2646                                }
2647                                if (distanceUp < 0.0) {
2648                                    if (thetaUp*el < distanceUp)
2649                                        thetaUp = distanceUp / el;
2650                                    improvementUp -= el;
2651                                }
2652                            }
2653                        }
2654                    }
2655                    double thetaDown = solValue - columnLower[iColumn];
2656                    double improvementDown = 0.0;
2657                    if (thetaDown > primalTolerance) {
2658                        // can go down
2659                        for (CoinBigIndex j = columnStart[iColumn];
2660                                j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2661                            int iRow = row[j];
2662                            double distanceUp = rowUpper[iRow] - rowActivity[iRow];
2663                            double distanceDown = rowLower[iRow] - rowActivity[iRow];
2664                            double el = -element[j]; // not change in sign form up
2665                            if (el > 0.0) {
2666                                // positive element
2667                                if (distanceUp > 0.0) {
2668                                    if (thetaDown*el > distanceUp)
2669                                        thetaDown = distanceUp / el;
2670                                } else {
2671                                    improvementDown -= el;
2672                                }
2673                                if (distanceDown > 0.0) {
2674                                    if (thetaDown*el > distanceDown)
2675                                        thetaDown = distanceDown / el;
2676                                    improvementDown += el;
2677                                }
2678                            } else {
2679                                // negative element
2680                                if (distanceDown < 0.0) {
2681                                    if (thetaDown*el < distanceDown)
2682                                        thetaDown = distanceDown / el;
2683                                } else {
2684                                    improvementDown += el;
2685                                }
2686                                if (distanceUp < 0.0) {
2687                                    if (thetaDown*el < distanceUp)
2688                                        thetaDown = distanceUp / el;
2689                                    improvementDown -= el;
2690                                }
2691                            }
2692                        }
2693                        if (thetaUp < 1.0e-8)
2694                            improvementUp = 0.0;
2695                        if (thetaDown < 1.0e-8)
2696                            improvementDown = 0.0;
2697                        double theta;
2698                        if (improvementUp >= improvementDown) {
2699                            theta = thetaUp;
2700                        } else {
2701                            improvementUp = improvementDown;
2702                            theta = -thetaDown;
2703                        }
2704                        if (improvementUp > 1.0e-8 && fabs(theta) > 1.0e-8) {
2705                            // Could move
2706                            double oldSum = 0.0;
2707                            double newSum = 0.0;
2708                            solution[iColumn] += theta;
2709                            for (CoinBigIndex j = columnStart[iColumn];
2710                                    j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2711                                int iRow = row[j];
2712                                double lower = rowLower[iRow];
2713                                double upper = rowUpper[iRow];
2714                                double value = rowActivity[iRow];
2715                                if (value > upper)
2716                                    oldSum += value - upper;
2717                                else if (value < lower)
2718                                    oldSum += lower - value;
2719                                value += theta * element[j];
2720                                rowActivity[iRow] = value;
2721                                if (value > upper)
2722                                    newSum += value - upper;
2723                                else if (value < lower)
2724                                    newSum += lower - value;
2725                            }
2726                            assert (newSum <= oldSum);
2727                            sumInfeasibility += newSum - oldSum;
2728                        }
2729                    }
2730                }
2731            }
2732            // Now flip some integers?
2733#ifdef JJF_ZERO
2734            for (i = 0; i < numberIntegers; i++) {
2735                int iColumn = integerVariable[i];
2736                double solValue = solution[iColumn];
2737                assert (fabs(solValue - floor(solValue + 0.5)) < 1.0e-8);
2738                double improvementUp = 0.0;
2739                if (columnUpper[iColumn] >= solValue + 1.0) {
2740                    // can go up
2741                    double oldSum = 0.0;
2742                    double newSum = 0.0;
2743                    for (CoinBigIndex j = columnStart[iColumn];
2744                            j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2745                        int iRow = row[j];
2746                        double lower = rowLower[iRow];
2747                        double upper = rowUpper[iRow];
2748                        double value = rowActivity[iRow];
2749                        if (value > upper)
2750                            oldSum += value - upper;
2751                        else if (value < lower)
2752                            oldSum += lower - value;
2753                        value += element[j];
2754                        if (value > upper)
2755                            newSum += value - upper;
2756                        else if (value < lower)
2757                            newSum += lower - value;
2758                    }
2759                    improvementUp = oldSum - newSum;
2760                }
2761                double improvementDown = 0.0;
2762                if (columnLower[iColumn] <= solValue - 1.0) {
2763                    // can go down
2764                    double oldSum = 0.0;
2765                    double newSum = 0.0;
2766                    for (CoinBigIndex j = columnStart[iColumn];
2767                            j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2768                        int iRow = row[j];
2769                        double lower = rowLower[iRow];
2770                        double upper = rowUpper[iRow];
2771                        double value = rowActivity[iRow];
2772                        if (value > upper)
2773                            oldSum += value - upper;
2774                        else if (value < lower)
2775                            oldSum += lower - value;
2776                        value -= element[j];
2777                        if (value > upper)
2778                            newSum += value - upper;
2779                        else if (value < lower)
2780                            newSum += lower - value;
2781                    }
2782                    improvementDown = oldSum - newSum;
2783                }
2784                double theta;
2785                if (improvementUp >= improvementDown) {
2786                    theta = 1.0;
2787                } else {
2788                    improvementUp = improvementDown;
2789                    theta = -1.0;
2790                }
2791                if (improvementUp > 1.0e-8 && fabs(theta) > 1.0e-8) {
2792                    // Could move
2793                    double oldSum = 0.0;
2794                    double newSum = 0.0;
2795                    solution[iColumn] += theta;
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 += theta * element[j];
2807                        rowActivity[iRow] = value;
2808                        if (value > upper)
2809                            newSum += value - upper;
2810                        else if (value < lower)
2811                            newSum += lower - value;
2812                    }
2813                    assert (newSum <= oldSum);
2814                    sumInfeasibility += newSum - oldSum;
2815                }
2816            }
2817#else
2818            int bestColumn = -1;
2819            double bestImprovement = primalTolerance;
2820            double theta = 0.0;
2821            for (i = 0; i < numberIntegers; i++) {
2822                int iColumn = integerVariable[i];
2823                double solValue = solution[iColumn];
2824                assert (fabs(solValue - floor(solValue + 0.5)) < 1.0e-8);
2825                double improvementUp = 0.0;
2826                if (columnUpper[iColumn] >= solValue + 1.0) {
2827                    // can go up
2828                    double oldSum = 0.0;
2829                    double newSum = 0.0;
2830                    for (CoinBigIndex j = columnStart[iColumn];
2831                            j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2832                        int iRow = row[j];
2833                        double lower = rowLower[iRow];
2834                        double upper = rowUpper[iRow];
2835                        double value = rowActivity[iRow];
2836                        if (value > upper)
2837                            oldSum += value - upper;
2838                        else if (value < lower)
2839                            oldSum += lower - value;
2840                        value += element[j];
2841                        if (value > upper)
2842                            newSum += value - upper;
2843                        else if (value < lower)
2844                            newSum += lower - value;
2845                    }
2846                    improvementUp = oldSum - newSum;
2847                }
2848                double improvementDown = 0.0;
2849                if (columnLower[iColumn] <= solValue - 1.0) {
2850                    // can go down
2851                    double oldSum = 0.0;
2852                    double newSum = 0.0;
2853                    for (CoinBigIndex j = columnStart[iColumn];
2854                            j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2855                        int iRow = row[j];
2856                        double lower = rowLower[iRow];
2857                        double upper = rowUpper[iRow];
2858                        double value = rowActivity[iRow];
2859                        if (value > upper)
2860                            oldSum += value - upper;
2861                        else if (value < lower)
2862                            oldSum += lower - value;
2863                        value -= element[j];
2864                        if (value > upper)
2865                            newSum += value - upper;
2866                        else if (value < lower)
2867                            newSum += lower - value;
2868                    }
2869                    improvementDown = oldSum - newSum;
2870                }
2871                double improvement = CoinMax(improvementUp, improvementDown);
2872                if (improvement > bestImprovement) {
2873                    bestImprovement = improvement;
2874                    bestColumn = iColumn;
2875                    if (improvementUp > improvementDown)
2876                        theta = 1.0;
2877                    else
2878                        theta = -1.0;
2879                }
2880            }
2881            if (bestColumn >= 0) {
2882                // Could move
2883                int iColumn = bestColumn;
2884                double oldSum = 0.0;
2885                double newSum = 0.0;
2886                solution[iColumn] += theta;
2887                for (CoinBigIndex j = columnStart[iColumn];
2888                        j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2889                    int iRow = row[j];
2890                    double lower = rowLower[iRow];
2891                    double upper = rowUpper[iRow];
2892                    double value = rowActivity[iRow];
2893                    if (value > upper)
2894                        oldSum += value - upper;
2895                    else if (value < lower)
2896                        oldSum += lower - value;
2897                    value += theta * element[j];
2898                    rowActivity[iRow] = value;
2899                    if (value > upper)
2900                        newSum += value - upper;
2901                    else if (value < lower)
2902                        newSum += lower - value;
2903                }
2904                assert (newSum <= oldSum);
2905                sumInfeasibility += newSum - oldSum;
2906            }
2907#endif
2908            if (oldSum <= sumInfeasibility + primalTolerance)
2909                break; // no good
2910        }
2911    }
2912    //delete [] saveSolution;
2913#endif
2914    delete [] rowActivity;
2915    return (largestInfeasibility > primalTolerance) ? 0 : 1;
2916}
2917// Set maximum Time (default off) - also sets starttime to current
2918void
2919CbcHeuristicFPump::setMaximumTime(double value)
2920{
2921    startTime_ = CoinCpuTime();
2922    maximumTime_ = value;
2923}
2924
2925# ifdef COIN_HAS_CLP
2926
2927//#############################################################################
2928// Constructors / Destructor / Assignment
2929//#############################################################################
2930
2931//-------------------------------------------------------------------
2932// Default Constructor
2933//-------------------------------------------------------------------
2934CbcDisasterHandler::CbcDisasterHandler (CbcModel * model)
2935        : OsiClpDisasterHandler(),
2936        cbcModel_(model)
2937{
2938    if (model) {
2939        osiModel_
2940        = dynamic_cast<OsiClpSolverInterface *> (model->solver());
2941        if (osiModel_)
2942            setSimplex(osiModel_->getModelPtr());
2943    }
2944}
2945
2946//-------------------------------------------------------------------
2947// Copy constructor
2948//-------------------------------------------------------------------
2949CbcDisasterHandler::CbcDisasterHandler (const CbcDisasterHandler & rhs)
2950        : OsiClpDisasterHandler(rhs),
2951        cbcModel_(rhs.cbcModel_)
2952{
2953}
2954
2955
2956//-------------------------------------------------------------------
2957// Destructor
2958//-------------------------------------------------------------------
2959CbcDisasterHandler::~CbcDisasterHandler ()
2960{
2961}
2962
2963//----------------------------------------------------------------
2964// Assignment operator
2965//-------------------------------------------------------------------
2966CbcDisasterHandler &
2967CbcDisasterHandler::operator=(const CbcDisasterHandler & rhs)
2968{
2969    if (this != &rhs) {
2970        OsiClpDisasterHandler::operator=(rhs);
2971        cbcModel_ = rhs.cbcModel_;
2972    }
2973    return *this;
2974}
2975//-------------------------------------------------------------------
2976// Clone
2977//-------------------------------------------------------------------
2978ClpDisasterHandler * CbcDisasterHandler::clone() const
2979{
2980    return new CbcDisasterHandler(*this);
2981}
2982// Type of disaster 0 can fix, 1 abort
2983int
2984CbcDisasterHandler::typeOfDisaster()
2985{
2986    if (!cbcModel_->parentModel() &&
2987            (cbcModel_->specialOptions()&2048) == 0) {
2988        return 0;
2989    } else {
2990        if (cbcModel_->parentModel())
2991            cbcModel_->setMaximumNodes(0);
2992        return 1;
2993    }
2994}
2995/* set model. */
2996void
2997CbcDisasterHandler::setCbcModel(CbcModel * model)
2998{
2999    cbcModel_ = model;
3000    if (model) {
3001        osiModel_
3002        = dynamic_cast<OsiClpSolverInterface *> (model->solver());
3003        if (osiModel_)
3004            setSimplex(osiModel_->getModelPtr());
3005        else
3006            setSimplex(NULL);
3007    }
3008}
3009#endif
3010
Note: See TracBrowser for help on using the repository browser.