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

Last change on this file since 2076 was 2072, checked in by unxusr, 5 years ago

fix several memory leaks in fpump, some reported by coverity

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