source: stable/2.9/Cbc/src/CbcHeuristicFPump.cpp @ 2186

Last change on this file since 2186 was 2186, checked in by stefan, 3 years ago

sync with trunk

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