/* $Id: CbcHeuristicFPump.cpp 1883 2013-04-06 13:33:15Z stefan $ */
// Copyright (C) 2004, International Business Machines
// Corporation and others.  All Rights Reserved.
// This code is licensed under the terms of the Eclipse Public License (EPL).

#if defined(_MSC_VER)
// Turn off compiler warning about long names
#  pragma warning(disable:4786)
#endif
#include <cassert>
#include <cstdlib>
#include <cmath>
#include <cfloat>

#include "OsiSolverInterface.hpp"
#include "CbcModel.hpp"
#ifdef COIN_HAS_CLP
#include "OsiClpSolverInterface.hpp"
#endif
#include "CbcMessage.hpp"
#include "CbcHeuristicFPump.hpp"
#include "CbcBranchActual.hpp"
#include "CbcBranchDynamic.hpp"
#include "CoinHelperFunctions.hpp"
#include "CoinWarmStartBasis.hpp"
#include "CoinTime.hpp"
#include "CbcEventHandler.hpp"


// Default Constructor
CbcHeuristicFPump::CbcHeuristicFPump()
        : CbcHeuristic(),
        startTime_(0.0),
        maximumTime_(0.0),
        fakeCutoff_(COIN_DBL_MAX),
        absoluteIncrement_(0.0),
        relativeIncrement_(0.0),
        defaultRounding_(0.49999),
        initialWeight_(0.0),
        weightFactor_(0.1),
        artificialCost_(COIN_DBL_MAX),
        iterationRatio_(0.0),
        reducedCostMultiplier_(1.0),
        maximumPasses_(100),
        maximumRetries_(1),
        accumulate_(0),
        fixOnReducedCosts_(1),
        roundExpensive_(false)
{
    setWhen(1);
}

// Constructor from model
CbcHeuristicFPump::CbcHeuristicFPump(CbcModel & model,
                                     double downValue, bool roundExpensive)
        : CbcHeuristic(model),
        startTime_(0.0),
        maximumTime_(0.0),
        fakeCutoff_(COIN_DBL_MAX),
        absoluteIncrement_(0.0),
        relativeIncrement_(0.0),
        defaultRounding_(downValue),
        initialWeight_(0.0),
        weightFactor_(0.1),
        artificialCost_(COIN_DBL_MAX),
        iterationRatio_(0.0),
        reducedCostMultiplier_(1.0),
        maximumPasses_(100),
        maximumRetries_(1),
        accumulate_(0),
        fixOnReducedCosts_(1),
        roundExpensive_(roundExpensive)
{
    setWhen(1);
}

// Destructor
CbcHeuristicFPump::~CbcHeuristicFPump ()
{
}

// Clone
CbcHeuristic *
CbcHeuristicFPump::clone() const
{
    return new CbcHeuristicFPump(*this);
}
// Create C++ lines to get to current state
void
CbcHeuristicFPump::generateCpp( FILE * fp)
{
    CbcHeuristicFPump other;
    fprintf(fp, "0#include \"CbcHeuristicFPump.hpp\"\n");
    fprintf(fp, "3  CbcHeuristicFPump heuristicFPump(*cbcModel);\n");
    CbcHeuristic::generateCpp(fp, "heuristicFPump");
    if (maximumPasses_ != other.maximumPasses_)
        fprintf(fp, "3  heuristicFPump.setMaximumPasses(%d);\n", maximumPasses_);
    else
        fprintf(fp, "4  heuristicFPump.setMaximumPasses(%d);\n", maximumPasses_);
    if (maximumRetries_ != other.maximumRetries_)
        fprintf(fp, "3  heuristicFPump.setMaximumRetries(%d);\n", maximumRetries_);
    else
        fprintf(fp, "4  heuristicFPump.setMaximumRetries(%d);\n", maximumRetries_);
    if (accumulate_ != other.accumulate_)
        fprintf(fp, "3  heuristicFPump.setAccumulate(%d);\n", accumulate_);
    else
        fprintf(fp, "4  heuristicFPump.setAccumulate(%d);\n", accumulate_);
    if (fixOnReducedCosts_ != other.fixOnReducedCosts_)
        fprintf(fp, "3  heuristicFPump.setFixOnReducedCosts(%d);\n", fixOnReducedCosts_);
    else
        fprintf(fp, "4  heuristicFPump.setFixOnReducedCosts(%d);\n", fixOnReducedCosts_);
    if (maximumTime_ != other.maximumTime_)
        fprintf(fp, "3  heuristicFPump.setMaximumTime(%g);\n", maximumTime_);
    else
        fprintf(fp, "4  heuristicFPump.setMaximumTime(%g);\n", maximumTime_);
    if (fakeCutoff_ != other.fakeCutoff_)
        fprintf(fp, "3  heuristicFPump.setFakeCutoff(%g);\n", fakeCutoff_);
    else
        fprintf(fp, "4  heuristicFPump.setFakeCutoff(%g);\n", fakeCutoff_);
    if (absoluteIncrement_ != other.absoluteIncrement_)
        fprintf(fp, "3  heuristicFPump.setAbsoluteIncrement(%g);\n", absoluteIncrement_);
    else
        fprintf(fp, "4  heuristicFPump.setAbsoluteIncrement(%g);\n", absoluteIncrement_);
    if (relativeIncrement_ != other.relativeIncrement_)
        fprintf(fp, "3  heuristicFPump.setRelativeIncrement(%g);\n", relativeIncrement_);
    else
        fprintf(fp, "4  heuristicFPump.setRelativeIncrement(%g);\n", relativeIncrement_);
    if (defaultRounding_ != other.defaultRounding_)
        fprintf(fp, "3  heuristicFPump.setDefaultRounding(%g);\n", defaultRounding_);
    else
        fprintf(fp, "4  heuristicFPump.setDefaultRounding(%g);\n", defaultRounding_);
    if (initialWeight_ != other.initialWeight_)
        fprintf(fp, "3  heuristicFPump.setInitialWeight(%g);\n", initialWeight_);
    else
        fprintf(fp, "4  heuristicFPump.setInitialWeight(%g);\n", initialWeight_);
    if (weightFactor_ != other.weightFactor_)
        fprintf(fp, "3  heuristicFPump.setWeightFactor(%g);\n", weightFactor_);
    else
        fprintf(fp, "4  heuristicFPump.setWeightFactor(%g);\n", weightFactor_);
    if (artificialCost_ != other.artificialCost_)
        fprintf(fp, "3  heuristicFPump.setArtificialCost(%g);\n", artificialCost_);
    else
        fprintf(fp, "4  heuristicFPump.setArtificialCost(%g);\n", artificialCost_);
    if (iterationRatio_ != other.iterationRatio_)
        fprintf(fp, "3  heuristicFPump.setIterationRatio(%g);\n", iterationRatio_);
    else
        fprintf(fp, "4  heuristicFPump.setIterationRatio(%g);\n", iterationRatio_);
    if (reducedCostMultiplier_ != other.reducedCostMultiplier_)
        fprintf(fp, "3  heuristicFPump.setReducedCostMultiplier(%g);\n", reducedCostMultiplier_);
    else
        fprintf(fp, "4  heuristicFPump.setReducedCostMultiplier(%g);\n", reducedCostMultiplier_);
    fprintf(fp, "3  cbcModel->addHeuristic(&heuristicFPump);\n");
}

// Copy constructor
CbcHeuristicFPump::CbcHeuristicFPump(const CbcHeuristicFPump & rhs)
        :
        CbcHeuristic(rhs),
        startTime_(rhs.startTime_),
        maximumTime_(rhs.maximumTime_),
        fakeCutoff_(rhs.fakeCutoff_),
        absoluteIncrement_(rhs.absoluteIncrement_),
        relativeIncrement_(rhs.relativeIncrement_),
        defaultRounding_(rhs.defaultRounding_),
        initialWeight_(rhs.initialWeight_),
        weightFactor_(rhs.weightFactor_),
        artificialCost_(rhs.artificialCost_),
        iterationRatio_(rhs.iterationRatio_),
        reducedCostMultiplier_(rhs.reducedCostMultiplier_),
        maximumPasses_(rhs.maximumPasses_),
        maximumRetries_(rhs.maximumRetries_),
        accumulate_(rhs.accumulate_),
        fixOnReducedCosts_(rhs.fixOnReducedCosts_),
        roundExpensive_(rhs.roundExpensive_)
{
}

// Assignment operator
CbcHeuristicFPump &
CbcHeuristicFPump::operator=( const CbcHeuristicFPump & rhs)
{
    if (this != &rhs) {
        CbcHeuristic::operator=(rhs);
        startTime_ = rhs.startTime_;
        maximumTime_ = rhs.maximumTime_;
        fakeCutoff_ = rhs.fakeCutoff_;
        absoluteIncrement_ = rhs.absoluteIncrement_;
        relativeIncrement_ = rhs.relativeIncrement_;
        defaultRounding_ = rhs.defaultRounding_;
        initialWeight_ = rhs.initialWeight_;
        weightFactor_ = rhs.weightFactor_;
        artificialCost_ = rhs.artificialCost_;
        iterationRatio_ = rhs.iterationRatio_;
        reducedCostMultiplier_ = rhs.reducedCostMultiplier_;
        maximumPasses_ = rhs.maximumPasses_;
        maximumRetries_ = rhs.maximumRetries_;
        accumulate_ = rhs.accumulate_;
        fixOnReducedCosts_ = rhs.fixOnReducedCosts_;
        roundExpensive_ = rhs.roundExpensive_;
    }
    return *this;
}

// Resets stuff if model changes
void
CbcHeuristicFPump::resetModel(CbcModel * )
{
}

/**************************BEGIN MAIN PROCEDURE ***********************************/

// See if feasibility pump will give better solution
// Sets value of solution
// Returns 1 if solution, 0 if not
int
CbcHeuristicFPump::solution(double & solutionValue,
                            double * betterSolution)
{
    startTime_ = CoinCpuTime();
    numCouldRun_++;
    double incomingObjective = solutionValue;
#define LEN_PRINT 200
    char pumpPrint[LEN_PRINT];
    pumpPrint[0] = '\0';
/*
  Decide if we want to run. Standard values for when are described in
  CbcHeuristic.hpp. If we're off, or running only at root and this isn't the
  root, bail out.

  The double test (against phase, then atRoot and passNumber) has a fair bit
  of redundancy, but the results will differ depending on whether we're
  actually at the root of the main search tree or at the root of a small tree
  (recursive call to branchAndBound).

  FPump also supports some exotic values (11 -- 15) for when, described in
  CbcHeuristicFPump.hpp.
*/
    if (!when() || (when() == 1 && model_->phase() != 1))
        return 0; // switched off
    // See if at root node
    bool atRoot = model_->getNodeCount() == 0;
    int passNumber = model_->getCurrentPassNumber();
    // just do once
    if (!atRoot)
        return 0;
    int options = feasibilityPumpOptions_;
    if ((options % 1000000) > 0) {
        int kOption = options / 1000000;
        options = options % 1000000;
        /*
          Add 10 to do even if solution
          1 - do after cuts
          2 - do after cuts (not before)
          3 - not used do after every cut round (and after cuts)
          k not used do after every (k-2)th round
        */
        if (kOption < 10 && model_->getSolutionCount())
            return 0;
        if (model_->getSolutionCount())
            kOption = kOption % 10;
        bool good;
        if (kOption == 1) {
            good = (passNumber == 999999);
        } else if (kOption == 2) {
            good = (passNumber == 999999);
            passNumber = 2; // so won't run before
            //} else if (kOption==3) {
            //good = true;
        } else {
            //good = (((passNumber-1)%(kOption-2))==0);
            good = false;
        }
        if (passNumber != 1 && !good)
            return 0;
    } else {
        if (passNumber != 1)
            return 0;
    }
    // loop round doing repeated pumps
    double cutoff;
    model_->solver()->getDblParam(OsiDualObjectiveLimit, cutoff);
    double direction = model_->solver()->getObjSense();
    cutoff *= direction;
    int numberBandBsolutions = 0;
    double firstCutoff = fabs(cutoff);
    cutoff = CoinMin(cutoff, solutionValue);
    // check plausible and space for rounded solution
    int numberColumns = model_->getNumCols();
    int numberIntegers = model_->numberIntegers();
    const int * integerVariableOrig = model_->integerVariable();
    double iterationLimit = -1.0;
    //iterationRatio_=1.0;
    if (iterationRatio_ > 0.0)
        iterationLimit = (2 * model_->solver()->getNumRows() + 2 * numberColumns) *
                         iterationRatio_;
    int totalNumberIterations = 0;
    int averageIterationsPerTry = -1;
    int numberIterationsLastPass = 0;
    // 1. initially check 0-1
/*
  I'm skeptical of the above comment, but it's likely accurate as the default.
  Bit 4 or bit 8 needs to be set in order to consider working with general
  integers.
*/
    int i, j;
    int general = 0;
    int * integerVariable = new int[numberIntegers];
    const double * lower = model_->solver()->getColLower();
    const double * upper = model_->solver()->getColUpper();
    bool doGeneral = (accumulate_ & 4) != 0;
    j = 0;
/*
  Scan the objects, recording the columns and counting general integers.

  Seems like the NDEBUG tests could be made into an applicability test. If
  a scan of the objects reveals complex objects, just clean up and return
  failure.
*/
    for (i = 0; i < numberIntegers; i++) {
        int iColumn = integerVariableOrig[i];
#ifndef NDEBUG
        const OsiObject * object = model_->object(i);
        const CbcSimpleInteger * integerObject =
            dynamic_cast<const  CbcSimpleInteger *> (object);
        const OsiSimpleInteger * integerObject2 =
            dynamic_cast<const  OsiSimpleInteger *> (object);
        assert(integerObject || integerObject2);
#endif
        if (upper[iColumn] - lower[iColumn] > 1.000001) {
            general++;
            if (doGeneral)
                integerVariable[j++] = iColumn;
        } else {
            integerVariable[j++] = iColumn;
        }
    }
/*
  If 2/3 of integers are general integers, and we're not going to work with
  them, might as well go home.

  The else case is unclear to me. We reach it if general integers are less than
  2/3 of the total, or if either of bit 4 or 8 is set. But only bit 8 is used
  in the decision. (Let manyGen = 1 if more than 2/3 of integers are general
  integers. Then a k-map on manyGen, bit4, and bit8 shows it clearly.)

  So there's something odd here. In the case where bit4 = 1 and bit8 = 0,
  we've included general integers in integerVariable, but we're not going to
  process them.
*/
    if (general*3 > 2*numberIntegers && !doGeneral) {
        delete [] integerVariable;
        return 0;
    } else if ((accumulate_&4) == 0) {
        doGeneral = false;
        j = 0;
        for (i = 0; i < numberIntegers; i++) {
            int iColumn = integerVariableOrig[i];
            if (upper[iColumn] - lower[iColumn] < 1.000001)
                integerVariable[j++] = iColumn;
        }
    }
    if (!general)
        doGeneral = false;
#ifdef CLP_INVESTIGATE
    if (doGeneral)
        printf("DOing general with %d out of %d\n", general, numberIntegers);
#endif
/*
  This `closest solution' will satisfy integrality, but violate some other
  constraints?
*/
    // For solution closest to feasible if none found
    int * closestSolution = general ? NULL : new int[numberIntegers];
    double closestObjectiveValue = COIN_DBL_MAX;

    int numberIntegersOrig = numberIntegers;
    numberIntegers = j;
    double * newSolution = new double [numberColumns];
    double newSolutionValue = COIN_DBL_MAX;
    int maxSolutions = model_->getMaximumSolutions();
    int numberSolutions=0;
    bool solutionFound = false;
    int * usedColumn = NULL;
    double * lastSolution = NULL;
    int fixContinuous = 0;
    bool fixInternal = false;
    bool allSlack = false;
    if (when_ >= 21 && when_ <= 25) {
        when_ -= 10;
        allSlack = true;
    }
    double time1 = CoinCpuTime();
/*
  Obtain a relaxed lp solution.
*/
    model_->solver()->resolve();
    if (!model_->solver()->isProvenOptimal()) {
        // presumably max time or some such
        return 0;
    }
    numRuns_++;
    if (cutoff < 1.0e50 && false) {
        // Fix on djs
        double direction = model_->solver()->getObjSense() ;
        double gap = cutoff - model_->solver()->getObjValue() * direction ;
        double tolerance;
        model_->solver()->getDblParam(OsiDualTolerance, tolerance) ;
        if (gap > 0.0) {
            gap += 100.0 * tolerance;
            int nFix = model_->solver()->reducedCostFix(gap);
            printf("dj fixing fixed %d variables\n", nFix);
        }
    }
/*
  I have no idea why we're doing this, except perhaps that saveBasis will be
  automagically deleted on exit from the routine.
*/
    CoinWarmStartBasis saveBasis;
    CoinWarmStartBasis * basis =
        dynamic_cast<CoinWarmStartBasis *>(model_->solver()->getWarmStart()) ;
    if (basis) {
        saveBasis = * basis;
        delete basis;
    }
    double continuousObjectiveValue = model_->solver()->getObjValue() * model_->solver()->getObjSense();
    double * firstPerturbedObjective = NULL;
    double * firstPerturbedSolution = NULL;
    double firstPerturbedValue = COIN_DBL_MAX;
    if (when_ >= 11 && when_ <= 15) {
        fixInternal = when_ > 11 && when_ < 15;
        if (when_ < 13)
            fixContinuous = 0;
        else if (when_ != 14)
            fixContinuous = 1;
        else
            fixContinuous = 2;
        when_ = 1;
        if ((accumulate_&1) != 0) {
            usedColumn = new int [numberColumns];
            for (int i = 0; i < numberColumns; i++)
                usedColumn[i] = -1;
        }
        lastSolution = CoinCopyOfArray(model_->solver()->getColSolution(), numberColumns);
    }
    int finalReturnCode = 0;
    int totalNumberPasses = 0;
    int numberTries = 0;
    CoinWarmStartBasis bestBasis;
    bool exitAll = false;
    //double saveBestObjective = model_->getMinimizationObjValue();
    OsiSolverInterface * solver = NULL;
    double artificialFactor = 0.00001;
    // also try rounding!
    double * roundingSolution = new double[numberColumns];
    double roundingObjective = solutionValue;
    CbcRounding roundingHeuristic(*model_);
    int dualPass = 0;
    int secondPassOpt = 0;
#define RAND_RAND
#ifdef RAND_RAND
    int offRandom = 0;
#endif
    int maximumAllowed = -1;
    bool moreIterations = false;
    if (options > 0) {
        if (options >= 1000)
            maximumAllowed = options / 1000;
        int options2 = (options % 1000) / 100;
#ifdef RAND_RAND
        offRandom = options2 & 1;
#endif
        moreIterations = (options2 & 2) != 0;
        secondPassOpt = (options / 10) % 10;
        /* 1 to 7 - re-use solution
           8 use dual and current solution(ish)
           9 use dual and allslack
           1 - primal and mod obj
           2 - dual and mod obj
           3 - primal and no mod obj
           add 3 to redo current solution
        */
        if (secondPassOpt >= 8) {
            dualPass = secondPassOpt - 7;
            secondPassOpt = 0;
        }
    }
    // Number of passes to do
    int maximumPasses = maximumPasses_;
#ifdef COIN_HAS_CLP
    {
      OsiClpSolverInterface * clpSolver
        = dynamic_cast<OsiClpSolverInterface *> (model_->solver());
      if (clpSolver ) {
	if (maximumPasses == 30) {
	  if (clpSolver->fakeObjective())
            maximumPasses = 100; // feasibility problem?
	}
	randomNumberGenerator_.randomize();
	if (model_->getRandomSeed()!=-1)
	  clpSolver->getModelPtr()->setRandomSeed(randomNumberGenerator_.getSeed());
	clpSolver->getModelPtr()->randomNumberGenerator()->randomize();
      }
    }
#endif
#ifdef RAND_RAND
    double * randomFactor = new double [numberColumns];
    for (int i = 0; i < numberColumns; i++) {
        double value = floor(1.0e3 * randomNumberGenerator_.randomDouble());
        randomFactor[i] = 1.0 + value * 1.0e-4;
    }
#endif
    // guess exact multiple of objective
    double exactMultiple = model_->getCutoffIncrement();
    exactMultiple *= 2520;
    if (fabs(exactMultiple / 0.999 - floor(exactMultiple / 0.999 + 0.5)) < 1.0e-9)
        exactMultiple /= 2520.0 * 0.999;
    else if (fabs(exactMultiple - floor(exactMultiple + 0.5)) < 1.0e-9)
        exactMultiple /= 2520.0;
    else
        exactMultiple = 0.0;
    // check for rounding errors (only for integral case)
    if (fabs(exactMultiple - floor(exactMultiple + 0.5)) < 1.0e-8)
        exactMultiple = floor(exactMultiple + 0.5);
    //printf("exact multiple %g\n",exactMultiple);
    // Clone solver for rounding
    OsiSolverInterface * clonedSolver = cloneBut(2); // wasmodel_->solver()->clone();
    while (!exitAll) {
        // Cutoff rhs
        double useRhs = COIN_DBL_MAX;
        double useOffset = 0.0;
        int numberPasses = 0;
        artificialFactor *= 10.0;
        int lastMove = (!numberTries) ? -10 : 1000000;
        double lastSumInfeas = COIN_DBL_MAX;
        numberTries++;
        // Clone solver - otherwise annoys root node computations
        solver = cloneBut(2); // was model_->solver()->clone();
#ifdef COIN_HAS_CLP
        {
            OsiClpSolverInterface * clpSolver
            = dynamic_cast<OsiClpSolverInterface *> (solver);
            if (clpSolver) {
                // better to clean up using primal?
                ClpSimplex * lp = clpSolver->getModelPtr();
                int options = lp->specialOptions();
                lp->setSpecialOptions(options | 8192);
                //lp->setSpecialOptions(options|0x01000000);
#ifdef CLP_INVESTIGATE
                clpSolver->setHintParam(OsiDoReducePrint, false, OsiHintTry);
                lp->setLogLevel(CoinMax(1, lp->logLevel()));
#endif
            }
        }
#endif
        if (CoinMin(fakeCutoff_, cutoff) < 1.0e50) {
            // Fix on djs
            double direction = solver->getObjSense() ;
            double gap = CoinMin(fakeCutoff_, cutoff) - solver->getObjValue() * direction ;
            double tolerance;
            solver->getDblParam(OsiDualTolerance, tolerance) ;
            if (gap > 0.0 && (fixOnReducedCosts_ == 1 || (numberTries == 1 && fixOnReducedCosts_ == 2))) {
                gap += 100.0 * tolerance;
                gap *= reducedCostMultiplier_;
                int nFix = solver->reducedCostFix(gap);
                if (nFix) {
                    sprintf(pumpPrint, "Reduced cost fixing fixed %d variables on major pass %d", nFix, numberTries);
                    model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
                    << pumpPrint
                    << CoinMessageEol;
                    //pumpPrint[0]='\0';
                }
            }
        }
        // if cutoff exists then add constraint
        bool useCutoff = (fabs(cutoff) < 1.0e20 && (fakeCutoff_ != COIN_DBL_MAX || numberTries > 1));
        // but there may be a close one
        if (firstCutoff < 2.0*solutionValue && numberTries == 1 && CoinMin(cutoff, fakeCutoff_) < 1.0e20)
            useCutoff = true;
        if (useCutoff) {
            double rhs = CoinMin(cutoff, fakeCutoff_);
            const double * objective = solver->getObjCoefficients();
            int numberColumns = solver->getNumCols();
            int * which = new int[numberColumns];
            double * els = new double[numberColumns];
            int nel = 0;
            for (int i = 0; i < numberColumns; i++) {
                double value = objective[i];
                if (value) {
                    which[nel] = i;
                    els[nel++] = direction * value;
                }
            }
            solver->getDblParam(OsiObjOffset, useOffset);
#ifdef COIN_DEVELOP
            if (useOffset)
                printf("CbcHeuristicFPump obj offset %g\n", useOffset);
#endif
            useOffset *= direction;
            // Tweak rhs and save
            useRhs = rhs;
#ifdef JJF_ZERO
            double tempValue = 60.0 * useRhs;
            if (fabs(tempValue - floor(tempValue + 0.5)) < 1.0e-7 && rhs != fakeCutoff_) {
                // add a little
                useRhs += 1.0e-5;
            }
#endif
            solver->addRow(nel, which, els, -COIN_DBL_MAX, useRhs + useOffset);
            delete [] which;
            delete [] els;
            bool takeHint;
            OsiHintStrength strength;
            solver->getHintParam(OsiDoDualInResolve, takeHint, strength);
            solver->setHintParam(OsiDoDualInResolve, true, OsiHintDo);
            solver->resolve();
            solver->setHintParam(OsiDoDualInResolve, takeHint, strength);
            if (!solver->isProvenOptimal()) {
                // presumably max time or some such
                exitAll = true;
                break;
            }
        }
        solver->setDblParam(OsiDualObjectiveLimit, 1.0e50);
        solver->resolve();
        // Solver may not be feasible
        if (!solver->isProvenOptimal()) {
            exitAll = true;
            break;
        }
        const double * lower = solver->getColLower();
        const double * upper = solver->getColUpper();
        const double * solution = solver->getColSolution();
        if (lastSolution)
            memcpy(lastSolution, solution, numberColumns*sizeof(double));
        double primalTolerance;
        solver->getDblParam(OsiPrimalTolerance, primalTolerance);

        // 2 space for last rounded solutions
#define NUMBER_OLD 4
        double ** oldSolution = new double * [NUMBER_OLD];
        for (j = 0; j < NUMBER_OLD; j++) {
            oldSolution[j] = new double[numberColumns];
            for (i = 0; i < numberColumns; i++) oldSolution[j][i] = -COIN_DBL_MAX;
        }

        // 3. Replace objective with an initial 0-valued objective
        double * saveObjective = new double [numberColumns];
        memcpy(saveObjective, solver->getObjCoefficients(), numberColumns*sizeof(double));
        for (i = 0; i < numberColumns; i++) {
            solver->setObjCoeff(i, 0.0);
        }
        bool finished = false;
        double direction = solver->getObjSense();
        int returnCode = 0;
        bool takeHint;
        OsiHintStrength strength;
        solver->getHintParam(OsiDoDualInResolve, takeHint, strength);
        solver->setHintParam(OsiDoDualInResolve, false);
        //solver->messageHandler()->setLogLevel(0);

        // 4. Save objective offset so we can see progress
        double saveOffset;
        solver->getDblParam(OsiObjOffset, saveOffset);
        // Get amount for original objective
        double scaleFactor = 0.0;
#ifdef COIN_DEVELOP
        double largestCost = 0.0;
        int nArtificial = 0;
#endif
        for (i = 0; i < numberColumns; i++) {
            double value = saveObjective[i];
            scaleFactor += value * value;
#ifdef COIN_DEVELOP
            largestCost = CoinMax(largestCost, fabs(value));
            if (value*direction >= artificialCost_)
                nArtificial++;
#endif
        }
        if (scaleFactor)
            scaleFactor = (initialWeight_ * sqrt(static_cast<double> (numberIntegers))) / sqrt(scaleFactor);
#ifdef CLP_INVESTIGATE
#ifdef COIN_DEVELOP
        if (scaleFactor || nArtificial)
            printf("Using %g fraction of original objective (decay %g) - largest %g - %d artificials\n", scaleFactor, weightFactor_,
                   largestCost, nArtificial);
#else
        if (scaleFactor)
            printf("Using %g fraction of original objective (decay %g)\n",
                   scaleFactor, weightFactor_);
#endif
#endif
        // This is an array of sums of infeasibilities so can see if "bobbling"
#define SIZE_BOBBLE 20
        double saveSumInf[SIZE_BOBBLE];
        CoinFillN(saveSumInf, SIZE_BOBBLE, COIN_DBL_MAX);
        // 0 before doing anything
        int bobbleMode = 0;
        // 5. MAIN WHILE LOOP
        //bool newLineNeeded=false;
/*
  finished occurs exactly twice in this routine: immediately above, where it's
  set to false, and here in the loop condition.
*/
        while (!finished) {
            double newTrueSolutionValue = 0.0;
            double newSumInfeas = 0.0;
            int newNumberInfeas = 0;
            returnCode = 0;
            if (model_->maximumSecondsReached()) {
                exitAll = true;
                break;
            }
            // see what changed
            if (usedColumn) {
                for (i = 0; i < numberColumns; i++) {
                    if (fabs(solution[i] - lastSolution[i]) > 1.0e-8)
                        usedColumn[i] = numberPasses;
                    lastSolution[i] = solution[i];
                }
            }
            if (averageIterationsPerTry >= 0) {
                int n = totalNumberIterations - numberIterationsLastPass;
                double perPass = totalNumberIterations /
                                 (totalNumberPasses + numberPasses + 1.0e-5);
                perPass /= (solver->getNumRows() + numberColumns);
                double test = moreIterations ? 0.3 : 0.05;
                if (n > CoinMax(20000, 3*averageIterationsPerTry)
                        && (switches_&2) == 0 && maximumPasses<200 && perPass>test) {
                    exitAll = true;
                }
            }
            // Exit on exact total number if maximumPasses large
            if ((maximumPasses >= 200 || (switches_&2) != 0)
                    && numberPasses + totalNumberPasses >=
                    maximumPasses)
                exitAll = true;
            bool exitThis = false;
            if (iterationLimit < 0.0) {
                if (numberPasses >= maximumPasses) {
                    // If going well then keep going if maximumPasses small
                    if (lastMove < numberPasses - 4 || lastMove == 1000000)
                        exitThis = true;
                    if (maximumPasses > 20 || numberPasses >= 40)
                        exitThis = true;
                }
            }
            if (iterationLimit > 0.0 && totalNumberIterations > iterationLimit
                    && numberPasses > 15) {
                // exiting on iteration count
                exitAll = true;
            } else if (maximumPasses<30 && numberPasses>100) {
                // too many passes anyway
                exitAll = true;
            }
            if (maximumTime_ > 0.0 && CoinCpuTime() >= startTime_ + maximumTime_) {
                exitAll = true;
                // force exit
                switches_ |= 2048;
            }
            if (exitAll || exitThis)
                break;
            memcpy(newSolution, solution, numberColumns*sizeof(double));
            int flip;
            if (numberPasses == 0 && false) {
                // always use same seed
                randomNumberGenerator_.setSeed(987654321);
            }
            returnCode = rounds(solver, newSolution,/*saveObjective,*/
                                numberIntegers, integerVariable,
                                /*pumpPrint,*/numberPasses,
                                /*roundExpensive_,*/defaultRounding_, &flip);
            if (numberPasses == 0 && false) {
                // Make sure random will be different
                for (i = 1; i < numberTries; i++)
                    randomNumberGenerator_.randomDouble();
            }
            numberPasses++;
            if (returnCode) {
                // SOLUTION IS INTEGER
                // Put back correct objective
                for (i = 0; i < numberColumns; i++)
                    solver->setObjCoeff(i, saveObjective[i]);

                // solution - but may not be better
                // Compute using dot product
                solver->setDblParam(OsiObjOffset, saveOffset);
                newSolutionValue = -saveOffset;
                for (  i = 0 ; i < numberColumns ; i++ )
                    newSolutionValue += saveObjective[i] * newSolution[i];
                newSolutionValue *= direction;
                sprintf(pumpPrint, "Solution found of %g", newSolutionValue);
                model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
                << pumpPrint
                << CoinMessageEol;
                //newLineNeeded=false;
                if (newSolutionValue < solutionValue) {
                    double saveValue = solutionValue;
                    if (!doGeneral) {
                        int numberLeft = 0;
                        for (i = 0; i < numberIntegersOrig; i++) {
                            int iColumn = integerVariableOrig[i];
                            double value = floor(newSolution[iColumn] + 0.5);
                            if (solver->isBinary(iColumn)) {
                                solver->setColLower(iColumn, value);
                                solver->setColUpper(iColumn, value);
                            } else {
                                if (fabs(value - newSolution[iColumn]) > 1.0e-7)
                                    numberLeft++;
                            }
                        }
                        if (numberLeft) {
                            sprintf(pumpPrint, "Branch and bound needed to clear up %d general integers", numberLeft);
                            model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
                            << pumpPrint
                            << CoinMessageEol;
                            returnCode = smallBranchAndBound(solver, numberNodes_, newSolution, newSolutionValue,
                                                             solutionValue, "CbcHeuristicFpump");
                            if (returnCode < 0) {
                                if (returnCode == -2)
                                    exitAll = true;
                                returnCode = 0; // returned on size or event
                            }
                            if ((returnCode&2) != 0) {
                                // could add cut
                                returnCode &= ~2;
                            }
                            if (returnCode != 1)
                                newSolutionValue = saveValue;
                            if (returnCode && newSolutionValue < saveValue)
                                numberBandBsolutions++;
                        }
                    }
                    if (returnCode && newSolutionValue < saveValue) {
                        memcpy(betterSolution, newSolution, numberColumns*sizeof(double));
                        solutionFound = true;
                        if (exitNow(newSolutionValue))
                            exitAll = true;
                        CoinWarmStartBasis * basis =
                            dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ;
                        if (basis) {
                            bestBasis = * basis;
                            delete basis;
                            int action = model_->dealWithEventHandler(CbcEventHandler::heuristicSolution, newSolutionValue, betterSolution);
                            if (action == 0) {
                                double * saveOldSolution = CoinCopyOfArray(model_->bestSolution(), numberColumns);
                                double saveObjectiveValue = model_->getMinimizationObjValue();
                                model_->setBestSolution(betterSolution, numberColumns, newSolutionValue);
                                if (saveOldSolution && saveObjectiveValue < model_->getMinimizationObjValue())
                                    model_->setBestSolution(saveOldSolution, numberColumns, saveObjectiveValue);
                                delete [] saveOldSolution;
                            }
                            if (action == 0 || model_->maximumSecondsReached()) {
                                exitAll = true; // exit
                                break;
                            }
                        }
                        if ((accumulate_&1) != 0) {
                            model_->incrementUsed(betterSolution); // for local search
                        }
                        solutionValue = newSolutionValue;
                        solutionFound = true;
			numberSolutions++;
			if (numberSolutions>=maxSolutions)
			  exitAll = true;
                        if (general && saveValue != newSolutionValue) {
                            sprintf(pumpPrint, "Cleaned solution of %g", solutionValue);
                            model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
                            << pumpPrint
                            << CoinMessageEol;
                        }
                        if (exitNow(newSolutionValue))
                            exitAll = true;
                    } else {
                        sprintf(pumpPrint, "Mini branch and bound could not fix general integers");
                        model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
                        << pumpPrint
                        << CoinMessageEol;
                    }
                } else {
                    sprintf(pumpPrint, "After further testing solution no better than previous of %g", solutionValue);
                    model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
                    << pumpPrint
                    << CoinMessageEol;
                    //newLineNeeded=false;
                    returnCode = 0;
                }
                break;
            } else {
                // SOLUTION IS not INTEGER
                // 1. check for loop
                bool matched;
                for (int k = NUMBER_OLD - 1; k > 0; k--) {
                    double * b = oldSolution[k];
                    matched = true;
                    for (i = 0; i < numberIntegers; i++) {
                        int iColumn = integerVariable[i];
                        if (newSolution[iColumn] != b[iColumn]) {
                            matched = false;
                            break;
                        }
                    }
                    if (matched) break;
                }
                int numberPerturbed = 0;
                if (matched || numberPasses % 100 == 0) {
                    // perturbation
                    //sprintf(pumpPrint+strlen(pumpPrint)," perturbation applied");
                    //newLineNeeded=true;
                    double factorX[10] = {0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0};
                    double factor = 1.0;
                    double target = -1.0;
                    double * randomX = new double [numberIntegers];
                    for (i = 0; i < numberIntegers; i++)
                        randomX[i] = CoinMax(0.0, randomNumberGenerator_.randomDouble() - 0.3);
                    for (int k = 0; k < 10; k++) {
#ifdef COIN_DEVELOP_x
                        printf("kpass %d\n", k);
#endif
                        int numberX[10] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
                        for (i = 0; i < numberIntegers; i++) {
                            int iColumn = integerVariable[i];
                            double value = randomX[i];
                            double difference = fabs(solution[iColumn] - newSolution[iColumn]);
                            for (int j = 0; j < 10; j++) {
                                if (difference + value*factorX[j] > 0.5)
                                    numberX[j]++;
                            }
                        }
                        if (target < 0.0) {
                            if (numberX[9] <= 200)
                                break; // not very many changes
                            target = CoinMax(200.0, CoinMin(0.05 * numberX[9], 1000.0));
                        }
                        int iX = -1;
                        int iBand = -1;
                        for (i = 0; i < 10; i++) {
#ifdef COIN_DEVELOP_x
                            printf("** %d changed at %g\n", numberX[i], factorX[i]);
#endif
                            if (numberX[i] >= target && numberX[i] < 2.0*target && iX < 0)
                                iX = i;
                            if (iBand<0 && numberX[i]>target) {
                                iBand = i;
                                factor = factorX[i];
                            }
                        }
                        if (iX >= 0) {
                            factor = factorX[iX];
                            break;
                        } else {
                            assert (iBand >= 0);
                            double hi = factor;
                            double lo = (iBand > 0) ? factorX[iBand-1] : 0.0;
                            double diff = (hi - lo) / 9.0;
                            for (i = 0; i < 10; i++) {
                                factorX[i] = lo;
                                lo += diff;
                            }
                        }
                    }
                    for (i = 0; i < numberIntegers; i++) {
                        int iColumn = integerVariable[i];
                        double value = randomX[i];
                        double difference = fabs(solution[iColumn] - newSolution[iColumn]);
                        if (difference + value*factor > 0.5) {
                            numberPerturbed++;
                            if (newSolution[iColumn] < lower[iColumn] + primalTolerance) {
                                newSolution[iColumn] += 1.0;
                            } else if (newSolution[iColumn] > upper[iColumn] - primalTolerance) {
                                newSolution[iColumn] -= 1.0;
                            } else {
                                // general integer
                                if (difference + value > 0.75)
                                    newSolution[iColumn] += 1.0;
                                else
                                    newSolution[iColumn] -= 1.0;
                            }
                        }
                    }
                    delete [] randomX;
                } else {
                    for (j = NUMBER_OLD - 1; j > 0; j--) {
                        for (i = 0; i < numberColumns; i++) oldSolution[j][i] = oldSolution[j-1][i];
                    }
                    for (j = 0; j < numberColumns; j++) oldSolution[0][j] = newSolution[j];
                }

                // 2. update the objective function based on the new rounded solution
                double offset = 0.0;
                double costValue = (1.0 - scaleFactor) * solver->getObjSense();
                int numberChanged = 0;
                const double * oldObjective = solver->getObjCoefficients();
                for (i = 0; i < numberColumns; i++) {
                    // below so we can keep original code and allow for objective
                    int iColumn = i;
                    // Special code for "artificials"
                    if (direction*saveObjective[iColumn] >= artificialCost_) {
                        //solver->setObjCoeff(iColumn,scaleFactor*saveObjective[iColumn]);
                        solver->setObjCoeff(iColumn, (artificialFactor*saveObjective[iColumn]) / artificialCost_);
                    }
                    if (!solver->isBinary(iColumn) && !doGeneral)
                        continue;
                    // deal with fixed variables (i.e., upper=lower)
                    if (fabs(lower[iColumn] - upper[iColumn]) < primalTolerance || !solver->isInteger(iColumn)) {
                        //if (lower[iColumn] > 1. - primalTolerance) solver->setObjCoeff(iColumn,-costValue);
                        //else                                       solver->setObjCoeff(iColumn,costValue);
                        continue;
                    }
                    double newValue = 0.0;
                    if (newSolution[iColumn] < lower[iColumn] + primalTolerance) {
                        newValue = costValue + scaleFactor * saveObjective[iColumn];
                    } else {
                        if (newSolution[iColumn] > upper[iColumn] - primalTolerance) {
                            newValue = -costValue + scaleFactor * saveObjective[iColumn];
                        }
                    }
#ifdef RAND_RAND
                    if (!offRandom)
                        newValue *= randomFactor[iColumn];
#endif
                    if (newValue != oldObjective[iColumn]) {
                        numberChanged++;
                    }
                    solver->setObjCoeff(iColumn, newValue);
                    offset += costValue * newSolution[iColumn];
                }
		if (numberPasses==1 && !totalNumberPasses && (model_->specialOptions()&8388608)!=0) {
		  // doing multiple solvers - make a real difference - flip 5%
		  for (i = 0; i < numberIntegers; i++) {
		    int iColumn = integerVariable[i];
		    double value = floor(newSolution[iColumn]+0.5);
		    if (fabs(value-solution[iColumn])>primalTolerance) {
		      value = randomNumberGenerator_.randomDouble();
		      if(value<0.05) {
			//printf("Flipping %d - random %g\n",iColumn,value);
			solver->setObjCoeff(iColumn,-solver->getObjCoefficients()[iColumn]);
		      }
		    }
		  }
		}
                solver->setDblParam(OsiObjOffset, -offset);
                if (!general && false) {
                    // Solve in two goes - first keep satisfied ones fixed
                    double * saveLower = new double [numberIntegers];
                    double * saveUpper = new double [numberIntegers];
                    for (i = 0; i < numberIntegers; i++) {
                        int iColumn = integerVariable[i];
                        saveLower[i] = COIN_DBL_MAX;
                        saveUpper[i] = -COIN_DBL_MAX;
                        if (solution[iColumn] < lower[iColumn] + primalTolerance) {
                            saveUpper[i] = upper[iColumn];
                            solver->setColUpper(iColumn, lower[iColumn]);
                        } else if (solution[iColumn] > upper[iColumn] - primalTolerance) {
                            saveLower[i] = lower[iColumn];
                            solver->setColLower(iColumn, upper[iColumn]);
                        }
                    }
                    solver->resolve();
                    if (!solver->isProvenOptimal()) {
                        // presumably max time or some such
                        exitAll = true;
                        break;
                    }
                    for (i = 0; i < numberIntegers; i++) {
                        int iColumn = integerVariable[i];
                        if (saveLower[i] != COIN_DBL_MAX)
                            solver->setColLower(iColumn, saveLower[i]);
                        if (saveUpper[i] != -COIN_DBL_MAX)
                            solver->setColUpper(iColumn, saveUpper[i]);
                        saveUpper[i] = -COIN_DBL_MAX;
                    }
                    memcpy(newSolution, solution, numberColumns*sizeof(double));
                    int flip;
                    returnCode = rounds(solver, newSolution,/*saveObjective,*/
                                        numberIntegers, integerVariable,
                                        /*pumpPrint,*/numberPasses,
                                        /*roundExpensive_,*/defaultRounding_, &flip);
                    numberPasses++;
                    if (returnCode) {
                        // solution - but may not be better
                        // Compute using dot product
                        double newSolutionValue = -saveOffset;
                        for (  i = 0 ; i < numberColumns ; i++ )
                            newSolutionValue += saveObjective[i] * newSolution[i];
                        newSolutionValue *= direction;
                        sprintf(pumpPrint, "Intermediate solution found of %g", newSolutionValue);
                        model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
                        << pumpPrint
                        << CoinMessageEol;
                        if (newSolutionValue < solutionValue) {
                            memcpy(betterSolution, newSolution, numberColumns*sizeof(double));
                            CoinWarmStartBasis * basis =
                                dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ;
                            solutionFound = true;
			    numberSolutions++;
			    if (numberSolutions>=maxSolutions)
			      exitAll = true;
                            if (exitNow(newSolutionValue))
                                exitAll = true;
                            if (basis) {
                                bestBasis = * basis;
                                delete basis;
                                int action = model_->dealWithEventHandler(CbcEventHandler::heuristicSolution, newSolutionValue, betterSolution);
                                if (!action) {
                                    double * saveOldSolution = CoinCopyOfArray(model_->bestSolution(), numberColumns);
                                    double saveObjectiveValue = model_->getMinimizationObjValue();
                                    model_->setBestSolution(betterSolution, numberColumns, newSolutionValue);
                                    if (saveOldSolution && saveObjectiveValue < model_->getMinimizationObjValue())
                                        model_->setBestSolution(saveOldSolution, numberColumns, saveObjectiveValue);
                                    delete [] saveOldSolution;
                                }
                                if (!action || model_->maximumSecondsReached()) {
                                    exitAll = true; // exit
                                    break;
                                }
                            }
                            if ((accumulate_&1) != 0) {
                                model_->incrementUsed(betterSolution); // for local search
                            }
                            solutionValue = newSolutionValue;
                            solutionFound = true;
			    numberSolutions++;
			    if (numberSolutions>=maxSolutions)
			      exitAll = true;
                            if (exitNow(newSolutionValue))
                                exitAll = true;
                        } else {
                            returnCode = 0;
                        }
                    }
                }
                int numberIterations = 0;
                if (!doGeneral) {
                    // faster to do from all slack!!!!
                    if (allSlack) {
                        CoinWarmStartBasis dummy;
                        solver->setWarmStart(&dummy);
                    }
#ifdef COIN_DEVELOP
                    printf("%d perturbed out of %d columns (%d changed)\n", numberPerturbed, numberColumns, numberChanged);
#endif
                    bool takeHint;
                    OsiHintStrength strength;
                    solver->getHintParam(OsiDoDualInResolve, takeHint, strength);
                    if (dualPass && numberChanged > 2) {
                        solver->setHintParam(OsiDoDualInResolve, true); // dual may be better
                        if (dualPass == 1 && 2*numberChanged < numberColumns &&
                                (numberChanged < 5000 || 6*numberChanged < numberColumns)) {
                            // but we need to make infeasible
                            CoinWarmStartBasis * basis =
                                dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ;
                            if (basis) {
                                // modify
                                const double * lower = solver->getColLower();
                                const double * upper = solver->getColUpper();
                                double * solution = CoinCopyOfArray(solver->getColSolution(),
                                                                    numberColumns);
                                const double * objective = solver->getObjCoefficients();
                                int nChanged = 0;
                                for (i = 0; i < numberIntegersOrig; i++) {
                                    int iColumn = integerVariableOrig[i];
#ifdef RAND_RAND
                                    if (nChanged > numberChanged)
                                        break;
#endif
                                    if (objective[iColumn] > 0.0) {
                                        if (basis->getStructStatus(iColumn) ==
                                                CoinWarmStartBasis::atUpperBound) {
                                            solution[iColumn] = lower[iColumn];
                                            basis->setStructStatus(iColumn, CoinWarmStartBasis::atLowerBound);
                                            nChanged++;
                                        }
                                    } else if (objective[iColumn] < 0.0) {
                                        if (basis->getStructStatus(iColumn) ==
                                                CoinWarmStartBasis::atLowerBound) {
                                            solution[iColumn] = upper[iColumn];
                                            basis->setStructStatus(iColumn, CoinWarmStartBasis::atUpperBound);
                                            nChanged++;
                                        }
                                    }
                                }
                                if (!nChanged) {
                                    for (i = 0; i < numberIntegersOrig; i++) {
                                        int iColumn = integerVariableOrig[i];
                                        if (objective[iColumn] > 0.0) {
                                            if (basis->getStructStatus(iColumn) ==
                                                    CoinWarmStartBasis::basic) {
                                                solution[iColumn] = lower[iColumn];
                                                basis->setStructStatus(iColumn, CoinWarmStartBasis::atLowerBound);
                                                break;
                                            }
                                        } else if (objective[iColumn] < 0.0) {
                                            if (basis->getStructStatus(iColumn) ==
                                                    CoinWarmStartBasis::basic) {
                                                solution[iColumn] = upper[iColumn];
                                                basis->setStructStatus(iColumn, CoinWarmStartBasis::atUpperBound);
                                                break;
                                            }
                                        }
                                    }
                                }
                                solver->setColSolution(solution);
                                delete [] solution;
                                solver->setWarmStart(basis);
                                delete basis;
                            }
                        } else {
                            // faster to do from all slack!!!! ???
                            CoinWarmStartBasis dummy;
                            solver->setWarmStart(&dummy);
                        }
                    }
                    if (numberTries > 1 && numberPasses == 1 && firstPerturbedObjective) {
                        // Modify to use convex combination
                        // use basis from first time
                        solver->setWarmStart(&saveBasis);
                        // and objective
                        if (secondPassOpt < 3 || (secondPassOpt >= 4 && secondPassOpt < 6))
                            solver->setObjective(firstPerturbedObjective);
                        // and solution
                        solver->setColSolution(firstPerturbedSolution);
                        //if (secondPassOpt==2||secondPassOpt==5||
                        if (firstPerturbedValue > cutoff)
                            solver->setHintParam(OsiDoDualInResolve, true); // dual may be better
                    }
                    solver->resolve();
                    if (!solver->isProvenOptimal()) {
                        // presumably max time or some such
                        exitAll = true;
                        break;
                    }
                    solver->setHintParam(OsiDoDualInResolve, takeHint);
                    newTrueSolutionValue = -saveOffset;
                    newSumInfeas = 0.0;
                    newNumberInfeas = 0;
                    {
                        const double * newSolution = solver->getColSolution();
                        for (  i = 0 ; i < numberColumns ; i++ ) {
                            if (solver->isInteger(i)) {
                                double value = newSolution[i];
                                double nearest = floor(value + 0.5);
                                newSumInfeas += fabs(value - nearest);
                                if (fabs(value - nearest) > 1.0e-6)
                                    newNumberInfeas++;
                            }
                            newTrueSolutionValue += saveObjective[i] * newSolution[i];
                        }
                        newTrueSolutionValue *= direction;
                        if (numberPasses == 1 && secondPassOpt) {
                            if (numberTries == 1 || secondPassOpt > 3) {
                                // save basis
                                CoinWarmStartBasis * basis =
                                    dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ;
                                if (basis) {
                                    saveBasis = * basis;
                                    delete basis;
                                }
                                delete [] firstPerturbedObjective;
                                delete [] firstPerturbedSolution;
                                firstPerturbedObjective = CoinCopyOfArray(solver->getObjCoefficients(), numberColumns);
                                firstPerturbedSolution = CoinCopyOfArray(solver->getColSolution(), numberColumns);
                                firstPerturbedValue = newTrueSolutionValue;
                            }
                        }
                        if (newNumberInfeas && newNumberInfeas < 15) {
#ifdef JJF_ZERO
                            roundingObjective = solutionValue;
                            OsiSolverInterface * saveSolver = model_->swapSolver(solver);
                            double * currentObjective =
                                CoinCopyOfArray(solver->getObjCoefficients(), numberColumns);
                            solver->setObjective(saveObjective);
                            double saveOffset2;
                            solver->getDblParam(OsiObjOffset, saveOffset2);
                            solver->setDblParam(OsiObjOffset, saveOffset);
                            int ifSol = roundingHeuristic.solution(roundingObjective, roundingSolution);
                            solver->setObjective(currentObjective);
                            solver->setDblParam(OsiObjOffset, saveOffset2);
                            delete [] currentObjective;
                            model_->swapSolver(saveSolver);
                            if (ifSol > 0)
                                abort();
#endif
                            int numberRows = solver->getNumRows();
                            double * rowActivity = new double[numberRows];
                            memset(rowActivity, 0, numberRows*sizeof(double));
                            int * which = new int[newNumberInfeas];
                            int * stack = new int[newNumberInfeas+1];
                            double * baseValue = new double[newNumberInfeas];
                            int * whichRow = new int[numberRows];
                            double * rowValue = new double[numberRows];
                            memset(rowValue, 0, numberRows*sizeof(double));
                            int nRow = 0;
                            // Column copy
                            const double * element = solver->getMatrixByCol()->getElements();
                            const int * row = solver->getMatrixByCol()->getIndices();
                            const CoinBigIndex * columnStart = solver->getMatrixByCol()->getVectorStarts();
                            const int * columnLength = solver->getMatrixByCol()->getVectorLengths();
                            int n = 0;
                            double contrib = 0.0;
                            for (  i = 0 ; i < numberColumns ; i++ ) {
                                double value = newSolution[i];
                                if (solver->isInteger(i)) {
                                    double nearest = floor(value + 0.5);
                                    if (fabs(value - nearest) > 1.0e-6) {
                                        //printf("Column %d value %g\n",i,value);
                                        for (CoinBigIndex j = columnStart[i];
                                                j < columnStart[i] + columnLength[i]; j++) {
                                            int iRow = row[j];
                                            //printf("row %d element %g\n",iRow,element[j]);
                                            if (!rowValue[iRow]) {
                                                rowValue[iRow] = 1.0;
                                                whichRow[nRow++] = iRow;
                                            }
                                        }
                                        baseValue[n] = floor(value);
                                        contrib += saveObjective[i] * value;
                                        value = 0.0;
                                        stack[n] = 0;
                                        which[n++] = i;
                                    }
                                }
                                for (CoinBigIndex j = columnStart[i];
                                        j < columnStart[i] + columnLength[i]; j++) {
                                    int iRow = row[j];
                                    rowActivity[iRow] += value * element[j];
                                }
                            }
                            if (newNumberInfeas < 15) {
                                stack[n] = newNumberInfeas + 100;
                                int iStack = n;
                                memset(rowValue, 0, numberRows*sizeof(double));
                                const double * rowLower = solver->getRowLower();
                                const double * rowUpper = solver->getRowUpper();
                                while (iStack >= 0) {
                                    double contrib2 = 0.0;
                                    // Could do faster
                                    for (int k = 0 ; k < n ; k++ ) {
                                        i = which[k];
                                        double value = baseValue[k] + stack[k];
                                        contrib2 += saveObjective[i] * value;
                                        for (CoinBigIndex j = columnStart[i];
                                                j < columnStart[i] + columnLength[i]; j++) {
                                            int iRow = row[j];
                                            rowValue[iRow] += value * element[j];
                                        }
                                    }
                                    // check if feasible
                                    bool feasible = true;
                                    for (int k = 0; k < nRow; k++) {
                                        i = whichRow[k];
                                        double value = rowValue[i] + rowActivity[i];
                                        rowValue[i] = 0.0;
                                        if (value < rowLower[i] - 1.0e-7 ||
                                                value > rowUpper[i] + 1.0e-7)
                                            feasible = false;
                                    }
                                    if (feasible) {
                                        double newObj = newTrueSolutionValue * direction;
                                        newObj += contrib2 - contrib;
                                        newObj *= direction;
#ifdef COIN_DEVELOP
                                        printf("FFFeasible! - obj %g\n", newObj);
#endif
                                        if (newObj < roundingObjective - 1.0e-6) {
#ifdef COIN_DEVELOP
                                            printf("FBetter\n");
#endif
                                            roundingObjective = newObj;
                                            memcpy(roundingSolution, newSolution, numberColumns*sizeof(double));
                                            for (int k = 0 ; k < n ; k++ ) {
                                                i = which[k];
                                                double value = baseValue[k] + stack[k];
                                                roundingSolution[i] = value;
                                            }
                                        }
                                    }
                                    while (iStack >= 0 && stack[iStack]) {
                                        stack[iStack]--;
                                        iStack--;
                                    }
                                    if (iStack >= 0) {
                                        stack[iStack] = 1;
                                        iStack = n;
                                        stack[n] = 1;
                                    }
                                }
                            }
                            delete [] rowActivity;
                            delete [] which;
                            delete [] stack;
                            delete [] baseValue;
                            delete [] whichRow;
                            delete [] rowValue;
                        }
                    }
                    if (true) {
                        OsiSolverInterface * saveSolver = model_->swapSolver(clonedSolver);
                        clonedSolver->setColSolution(solver->getColSolution());
                        CbcRounding heuristic1(*model_);
                        heuristic1.setHeuristicName("rounding in feaspump!");
                        heuristic1.setWhen(1);
                        roundingObjective = CoinMin(roundingObjective, solutionValue);
                        double testSolutionValue = newTrueSolutionValue;
                        int returnCode = heuristic1.solution(roundingObjective,
                                                             roundingSolution,
                                                             testSolutionValue) ;
                        if (returnCode == 1) {
#ifdef COIN_DEVELOP
                            printf("rounding obj of %g?\n", roundingObjective);
#endif
                            //roundingObjective = newSolutionValue;
                            //} else {
                            //roundingObjective = COIN_DBL_MAX;
                        }
                        model_->swapSolver(saveSolver);
                    }
                    if (!solver->isProvenOptimal()) {
                        // presumably max time or some such
                        exitAll = true;
                        break;
                    }
                    // in case very dubious solver
                    lower = solver->getColLower();
                    upper = solver->getColUpper();
                    solution = solver->getColSolution();
                    numberIterations = solver->getIterationCount();
                } else {
                    int * addStart = new int[2*general+1];
                    int * addIndex = new int[4*general];
                    double * addElement = new double[4*general];
                    double * addLower = new double[2*general];
                    double * addUpper = new double[2*general];
                    double * obj = new double[general];
                    int nAdd = 0;
                    for (i = 0; i < numberIntegers; i++) {
                        int iColumn = integerVariable[i];
                        if (newSolution[iColumn] > lower[iColumn] + primalTolerance &&
                                newSolution[iColumn] < upper[iColumn] - primalTolerance) {
                            assert (upper[iColumn] - lower[iColumn] > 1.00001);
                            obj[nAdd] = 1.0;
                            addLower[nAdd] = 0.0;
                            addUpper[nAdd] = COIN_DBL_MAX;
                            nAdd++;
                        }
                    }
                    OsiSolverInterface * solver2 = solver;
                    if (nAdd) {
                        CoinZeroN(addStart, nAdd + 1);
                        solver2 = solver->clone();
                        solver2->addCols(nAdd, addStart, NULL, NULL, addLower, addUpper, obj);
                        // feasible solution
                        double * sol = new double[nAdd+numberColumns];
                        memcpy(sol, solution, numberColumns*sizeof(double));
                        // now rows
                        int nAdd = 0;
                        int nEl = 0;
                        int nAddRow = 0;
                        for (i = 0; i < numberIntegers; i++) {
                            int iColumn = integerVariable[i];
                            if (newSolution[iColumn] > lower[iColumn] + primalTolerance &&
                                    newSolution[iColumn] < upper[iColumn] - primalTolerance) {
                                addLower[nAddRow] = -newSolution[iColumn];;
                                addUpper[nAddRow] = COIN_DBL_MAX;
                                addIndex[nEl] = iColumn;
                                addElement[nEl++] = -1.0;
                                addIndex[nEl] = numberColumns + nAdd;
                                addElement[nEl++] = 1.0;
                                nAddRow++;
                                addStart[nAddRow] = nEl;
                                addLower[nAddRow] = newSolution[iColumn];;
                                addUpper[nAddRow] = COIN_DBL_MAX;
                                addIndex[nEl] = iColumn;
                                addElement[nEl++] = 1.0;
                                addIndex[nEl] = numberColumns + nAdd;
                                addElement[nEl++] = 1.0;
                                nAddRow++;
                                addStart[nAddRow] = nEl;
                                sol[nAdd+numberColumns] = fabs(sol[iColumn] - newSolution[iColumn]);
                                nAdd++;
                            }
                        }
                        solver2->setColSolution(sol);
                        delete [] sol;
                        solver2->addRows(nAddRow, addStart, addIndex, addElement, addLower, addUpper);
                    }
                    delete [] addStart;
                    delete [] addIndex;
                    delete [] addElement;
                    delete [] addLower;
                    delete [] addUpper;
                    delete [] obj;
                    solver2->resolve();
                    if (!solver2->isProvenOptimal()) {
                        // presumably max time or some such
                        exitAll = true;
                        break;
                    }
                    //assert (solver2->isProvenOptimal());
                    if (nAdd) {
                        solver->setColSolution(solver2->getColSolution());
                        numberIterations = solver2->getIterationCount();
                        delete solver2;
                    } else {
                        numberIterations = solver->getIterationCount();
                    }
                    lower = solver->getColLower();
                    upper = solver->getColUpper();
                    solution = solver->getColSolution();
                    newTrueSolutionValue = -saveOffset;
                    newSumInfeas = 0.0;
                    newNumberInfeas = 0;
                    {
                        const double * newSolution = solver->getColSolution();
                        for (  i = 0 ; i < numberColumns ; i++ ) {
                            if (solver->isInteger(i)) {
                                double value = newSolution[i];
                                double nearest = floor(value + 0.5);
                                newSumInfeas += fabs(value - nearest);
                                if (fabs(value - nearest) > 1.0e-6)
                                    newNumberInfeas++;
                            }
                            newTrueSolutionValue += saveObjective[i] * newSolution[i];
                        }
                        newTrueSolutionValue *= direction;
                    }
                }
                if (lastMove != 1000000) {
                    if (newSumInfeas < lastSumInfeas) {
                        lastMove = numberPasses;
                        lastSumInfeas = newSumInfeas;
                    } else if (newSumInfeas > lastSumInfeas + 1.0e-5) {
                        lastMove = 1000000; // going up
                    }
                }
                totalNumberIterations += numberIterations;
                if (solver->getNumRows() < 3000)
                    sprintf(pumpPrint, "Pass %3d: suminf. %10.5f (%d) obj. %g iterations %d",
                            numberPasses + totalNumberPasses,
                            newSumInfeas, newNumberInfeas,
                            newTrueSolutionValue, numberIterations);
                else
                    sprintf(pumpPrint, "Pass %3d: (%.2f seconds) suminf. %10.5f (%d) obj. %g iterations %d", numberPasses + totalNumberPasses,
                            model_->getCurrentSeconds(), newSumInfeas, newNumberInfeas,
                            newTrueSolutionValue, numberIterations);
                model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
                << pumpPrint
                << CoinMessageEol;
                if (closestSolution && solver->getObjValue() < closestObjectiveValue) {
                    int i;
                    const double * objective = solver->getObjCoefficients();
                    for (i = 0; i < numberIntegersOrig; i++) {
                        int iColumn = integerVariableOrig[i];
                        if (objective[iColumn] > 0.0)
                            closestSolution[i] = 0;
                        else
                            closestSolution[i] = 1;
                    }
                    closestObjectiveValue = solver->getObjValue();
                }
                // See if we need to think about changing rhs
                if ((switches_&12) != 0 && useRhs < 1.0e50) {
                    double oldRhs = useRhs;
                    bool trying = false;
                    if ((switches_&4) != 0 && numberPasses && (numberPasses % 50) == 0) {
                        if (solutionValue > 1.0e20) {
                            // only if no genuine solution
                            double gap = useRhs - continuousObjectiveValue;
                            useRhs += 0.1 * gap;
                            if (exactMultiple) {
                                useRhs = exactMultiple * ceil(useRhs / exactMultiple);
                                useRhs = CoinMax(useRhs, oldRhs + exactMultiple);
                            }
                            trying = true;
                        }
                    }
                    if ((switches_&8) != 0) {
                        // Put in new suminf and check
                        double largest = newSumInfeas;
                        double smallest = newSumInfeas;
                        for (int i = 0; i < SIZE_BOBBLE - 1; i++) {
                            double value = saveSumInf[i+1];
                            saveSumInf[i] = value;
                            largest = CoinMax(largest, value);
                            smallest = CoinMin(smallest, value);
                        }
                        saveSumInf[SIZE_BOBBLE-1] = newSumInfeas;
                        if (smallest*1.5 > largest && smallest > 2.0) {
                            if (bobbleMode == 0) {
                                // go closer
                                double gap = oldRhs - continuousObjectiveValue;
                                useRhs -= 0.4 * gap;
                                if (exactMultiple) {
                                    double value = floor(useRhs / exactMultiple);
                                    useRhs = CoinMin(value * exactMultiple, oldRhs - exactMultiple);
                                }
                                if (useRhs < continuousObjectiveValue) {
                                    // skip decrease
                                    bobbleMode = 1;
                                    useRhs = oldRhs;
                                }
                            }
                            if (bobbleMode) {
                                trying = true;
                                // weaken
                                if (solutionValue < 1.0e20) {
                                    double gap = solutionValue - oldRhs;
                                    useRhs += 0.3 * gap;
                                } else {
                                    double gap = oldRhs - continuousObjectiveValue;
                                    useRhs += 0.05 * gap;
                                }
                                if (exactMultiple) {
                                    double value = ceil(useRhs / exactMultiple);
                                    useRhs = CoinMin(value * exactMultiple,
                                                     solutionValue - exactMultiple);
                                }
                            }
                            bobbleMode++;
                            // reset
                            CoinFillN(saveSumInf, SIZE_BOBBLE, COIN_DBL_MAX);
                        }
                    }
                    if (useRhs != oldRhs) {
                        // tidy up
                        if (exactMultiple) {
                            double value = floor(useRhs / exactMultiple);
                            double bestPossible = ceil(continuousObjectiveValue / exactMultiple);
                            useRhs = CoinMax(value, bestPossible) * exactMultiple;
                        } else {
                            useRhs = CoinMax(useRhs, continuousObjectiveValue);
                        }
                        int k = solver->getNumRows() - 1;
                        solver->setRowUpper(k, useRhs + useOffset);
                        bool takeHint;
                        OsiHintStrength strength;
                        solver->getHintParam(OsiDoDualInResolve, takeHint, strength);
                        if (useRhs < oldRhs) {
                            solver->setHintParam(OsiDoDualInResolve, true);
                            solver->resolve();
                        } else if (useRhs > oldRhs) {
                            solver->setHintParam(OsiDoDualInResolve, false);
                            solver->resolve();
                        }
                        solver->setHintParam(OsiDoDualInResolve, takeHint);
                        if (!solver->isProvenOptimal()) {
                            // presumably max time or some such
                            exitAll = true;
                            break;
                        }
                    } else if (trying) {
                        // doesn't look good
                        break;
                    }
                }
            }
            // reduce scale factor
            scaleFactor *= weightFactor_;
        } // END WHILE
        // see if rounding worked!
        if (roundingObjective < solutionValue) {
            if (roundingObjective < solutionValue - 1.0e-6*fabs(roundingObjective)) {
                sprintf(pumpPrint, "Rounding solution of %g is better than previous of %g\n",
                        roundingObjective, solutionValue);
                model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
                << pumpPrint
                << CoinMessageEol;
            }
            solutionValue = roundingObjective;
            newSolutionValue = solutionValue;
            memcpy(betterSolution, roundingSolution, numberColumns*sizeof(double));
            solutionFound = true;
	    numberSolutions++;
	    if (numberSolutions>=maxSolutions)
	      exitAll = true;
            if (exitNow(roundingObjective))
                exitAll = true;
        }
        if (!solutionFound) {
            sprintf(pumpPrint, "No solution found this major pass");
            model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
            << pumpPrint
            << CoinMessageEol;
        }
        //}
        delete solver;
        solver = NULL;
        for ( j = 0; j < NUMBER_OLD; j++)
            delete [] oldSolution[j];
        delete [] oldSolution;
        delete [] saveObjective;
        if (usedColumn && !exitAll) {
            OsiSolverInterface * newSolver = cloneBut(3); // was model_->continuousSolver()->clone();
#if 0 //def COIN_HAS_CLP
	    OsiClpSolverInterface * clpSolver
	      = dynamic_cast<OsiClpSolverInterface *> (newSolver);
	    if (clpSolver) {
	      ClpSimplex * simplex = clpSolver->getModelPtr();
	      simplex->writeMps("start.mps",2,1);
	    }
#endif
            const double * colLower = newSolver->getColLower();
            const double * colUpper = newSolver->getColUpper();
            bool stopBAB = false;
            int allowedPass = -1;
            if (maximumAllowed > 0)
                allowedPass = CoinMax(numberPasses - maximumAllowed, -1);
            while (!stopBAB) {
                stopBAB = true;
                int i;
                int nFix = 0;
                int nFixI = 0;
                int nFixC = 0;
                int nFixC2 = 0;
                for (i = 0; i < numberIntegersOrig; i++) {
                    int iColumn = integerVariableOrig[i];
                    //const OsiObject * object = model_->object(i);
                    //double originalLower;
                    //double originalUpper;
                    //getIntegerInformation( object,originalLower, originalUpper);
                    //assert(colLower[iColumn]==originalLower);
                    //newSolver->setColLower(iColumn,CoinMax(colLower[iColumn],originalLower));
                    newSolver->setColLower(iColumn, colLower[iColumn]);
                    //assert(colUpper[iColumn]==originalUpper);
                    //newSolver->setColUpper(iColumn,CoinMin(colUpper[iColumn],originalUpper));
                    newSolver->setColUpper(iColumn, colUpper[iColumn]);
                    if (usedColumn[iColumn] <= allowedPass) {
                        double value = lastSolution[iColumn];
                        double nearest = floor(value + 0.5);
                        if (fabs(value - nearest) < 1.0e-7) {
                            if (nearest == colLower[iColumn]) {
                                newSolver->setColUpper(iColumn, colLower[iColumn]);
                                nFix++;
                            } else if (nearest == colUpper[iColumn]) {
                                newSolver->setColLower(iColumn, colUpper[iColumn]);
                                nFix++;
                            } else if (fixInternal) {
                                newSolver->setColLower(iColumn, nearest);
                                newSolver->setColUpper(iColumn, nearest);
                                nFix++;
                                nFixI++;
                            }
                        }
                    }
                }
                if (fixContinuous) {
                    for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
                        if (!newSolver->isInteger(iColumn) && usedColumn[iColumn] <= allowedPass) {
                            double value = lastSolution[iColumn];
                            if (value < colLower[iColumn] + 1.0e-8) {
                                newSolver->setColUpper(iColumn, colLower[iColumn]);
                                nFixC++;
                            } else if (value > colUpper[iColumn] - 1.0e-8) {
                                newSolver->setColLower(iColumn, colUpper[iColumn]);
                                nFixC++;
                            } else if (fixContinuous == 2) {
                                newSolver->setColLower(iColumn, value);
                                newSolver->setColUpper(iColumn, value);
                                nFixC++;
                                nFixC2++;
                            }
                        }
                    }
                }
                newSolver->initialSolve();
                if (!newSolver->isProvenOptimal()) {
                    //newSolver->writeMps("bad.mps");
                    //assert (newSolver->isProvenOptimal());
                    exitAll = true;
                    break;
                }
                sprintf(pumpPrint, "Before mini branch and bound, %d integers at bound fixed and %d continuous",
                        nFix, nFixC);
                if (nFixC2 + nFixI != 0)
                    sprintf(pumpPrint + strlen(pumpPrint), " of which %d were internal integer and %d internal continuous",
                            nFixI, nFixC2);
                model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
                << pumpPrint
                << CoinMessageEol;
                double saveValue = newSolutionValue;
                if (newSolutionValue - model_->getCutoffIncrement()
                        > continuousObjectiveValue - 1.0e-7) {
                    double saveFraction = fractionSmall_;
                    if (numberTries > 1 && !numberBandBsolutions)
                        fractionSmall_ *= 0.5;
                    // Give branch and bound a bit more freedom
                    double cutoff2 = newSolutionValue +
                                     CoinMax(model_->getCutoffIncrement(), 1.0e-3);
#if 0
		      {
                        OsiClpSolverInterface * clpSolver
                        = dynamic_cast<OsiClpSolverInterface *> (newSolver);
                        if (clpSolver) {
                            ClpSimplex * simplex = clpSolver->getModelPtr();
                            simplex->writeMps("testA.mps",2,1);
			}
		      }
#endif
                    int returnCode2 = smallBranchAndBound(newSolver, numberNodes_, newSolution, newSolutionValue,
                                                          cutoff2, "CbcHeuristicLocalAfterFPump");
                    fractionSmall_ = saveFraction;
                    if (returnCode2 < 0) {
                        if (returnCode2 == -2) {
                            exitAll = true;
                            returnCode = 0;
                        } else {
                            returnCode2 = 0; // returned on size - try changing
                            //#define ROUND_AGAIN
#ifdef ROUND_AGAIN
                            if (numberTries == 1 && numberPasses > 20 && allowedPass < numberPasses - 1) {
                                allowedPass = (numberPasses + allowedPass) >> 1;
                                sprintf(pumpPrint,
                                        "Fixing all variables which were last changed on pass %d and trying again",
                                        allowedPass);
                                model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
                                << pumpPrint
                                << CoinMessageEol;
                                stopBAB = false;
                                continue;
                            }
#endif
                        }
                    }
                    if ((returnCode2&2) != 0) {
                        // could add cut
                        returnCode2 &= ~2;
                    }
                    if (returnCode2) {
                        numberBandBsolutions++;
			// may not have got solution earlier
			returnCode |= 1;
		    }
                } else {
                    // no need
                    exitAll = true;
                    //returnCode=0;
                }
                // recompute solution value
                if (returnCode && true) {
#if 0
		      {
                        OsiClpSolverInterface * clpSolver
                        = dynamic_cast<OsiClpSolverInterface *> (newSolver);
                        if (clpSolver) {
                            ClpSimplex * simplex = clpSolver->getModelPtr();
                            simplex->writeMps("testB.mps",2,1);
			}
		      }
#endif
                    delete newSolver;
                    newSolver = cloneBut(3); // was model_->continuousSolver()->clone();
                    newSolutionValue = -saveOffset;
                    double newSumInfeas = 0.0;
                    const double * obj = newSolver->getObjCoefficients();
                    for (int i = 0 ; i < numberColumns ; i++ ) {
                        if (newSolver->isInteger(i)) {
                            double value = newSolution[i];
                            double nearest = floor(value + 0.5);
                            newSumInfeas += fabs(value - nearest);
                        }
                        newSolutionValue += obj[i] * newSolution[i];
                    }
                    newSolutionValue *= direction;
                }
                bool gotSolution = false;
                if (returnCode && newSolutionValue < saveValue) {
                    sprintf(pumpPrint, "Mini branch and bound improved solution from %g to %g (%.2f seconds)",
                            saveValue, newSolutionValue, model_->getCurrentSeconds());
                    model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
                    << pumpPrint
                    << CoinMessageEol;
                    memcpy(betterSolution, newSolution, numberColumns*sizeof(double));
                    gotSolution = true;
                    if (fixContinuous && nFixC + nFixC2 > 0) {
                        // may be able to do even better
                        int nFixed = 0;
                        const double * lower = model_->solver()->getColLower();
                        const double * upper = model_->solver()->getColUpper();
                        for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
                            double value = newSolution[iColumn];
                            if (newSolver->isInteger(iColumn)) {
                                value = floor(newSolution[iColumn] + 0.5);
                                newSolver->setColLower(iColumn, value);
                                newSolver->setColUpper(iColumn, value);
                                nFixed++;
                            } else {
                                newSolver->setColLower(iColumn, lower[iColumn]);
                                newSolver->setColUpper(iColumn, upper[iColumn]);
                                if (value < lower[iColumn])
                                    value = lower[iColumn];
                                else if (value > upper[iColumn])
                                    value = upper[iColumn];

                            }
                            newSolution[iColumn] = value;
                        }
                        newSolver->setColSolution(newSolution);
                        //#define CLP_INVESTIGATE2
#ifdef CLP_INVESTIGATE2
                        {
                            // check
                            // get row activities
                            int numberRows = newSolver->getNumRows();
                            double * rowActivity = new double[numberRows];
                            memset(rowActivity, 0, numberRows*sizeof(double));
                            newSolver->getMatrixByCol()->times(newSolution, rowActivity) ;
                            double largestInfeasibility = primalTolerance;
                            double sumInfeasibility = 0.0;
                            int numberBadRows = 0;
                            const double * rowLower = newSolver->getRowLower();
                            const double * rowUpper = newSolver->getRowUpper();
                            for (i = 0 ; i < numberRows ; i++) {
                                double value;
                                value = rowLower[i] - rowActivity[i];
                                if (value > primalTolerance) {
                                    numberBadRows++;
                                    largestInfeasibility = CoinMax(largestInfeasibility, value);
                                    sumInfeasibility += value;
                                }
                                value = rowActivity[i] - rowUpper[i];
                                if (value > primalTolerance) {
                                    numberBadRows++;
                                    largestInfeasibility = CoinMax(largestInfeasibility, value);
                                    sumInfeasibility += value;
                                }
                            }
                            printf("%d bad rows, largest inf %g sum %g\n",
                                   numberBadRows, largestInfeasibility, sumInfeasibility);
                            delete [] rowActivity;
                        }
#endif
#ifdef COIN_HAS_CLP
                        OsiClpSolverInterface * clpSolver
                        = dynamic_cast<OsiClpSolverInterface *> (newSolver);
                        if (clpSolver) {
                            ClpSimplex * simplex = clpSolver->getModelPtr();
                            //simplex->writeBasis("test.bas",true,2);
                            //simplex->writeMps("test.mps",2,1);
                            if (nFixed*3 > numberColumns*2)
                                simplex->allSlackBasis(); // may as well go from all slack
                            simplex->primal(1);
                            clpSolver->setWarmStart(NULL);
                        }
#endif
                        newSolver->initialSolve();
                        if (newSolver->isProvenOptimal()) {
                            double value = newSolver->getObjValue() * newSolver->getObjSense();
                            if (value < newSolutionValue) {
                                //newSolver->writeMps("query","mps");
#ifdef JJF_ZERO
                                {
                                    double saveOffset;
                                    newSolver->getDblParam(OsiObjOffset, saveOffset);
                                    const double * obj = newSolver->getObjCoefficients();
                                    double newTrueSolutionValue = -saveOffset;
                                    double newSumInfeas = 0.0;
                                    int numberColumns = newSolver->getNumCols();
                                    const double * solution = newSolver->getColSolution();
                                    for (int  i = 0 ; i < numberColumns ; i++ ) {
                                        if (newSolver->isInteger(i)) {
                                            double value = solution[i];
                                            double nearest = floor(value + 0.5);
                                            newSumInfeas += fabs(value - nearest);
                                        }
                                        if (solution[i])
                                            printf("%d obj %g val %g - total %g\n", i, obj[i], solution[i],
                                                   newTrueSolutionValue);
                                        newTrueSolutionValue += obj[i] * solution[i];
                                    }
                                    printf("obj %g - inf %g\n", newTrueSolutionValue,
                                           newSumInfeas);
                                }
#endif
                                sprintf(pumpPrint, "Freeing continuous variables gives a solution of %g", value);
                                model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
                                << pumpPrint
                                << CoinMessageEol;
                                newSolutionValue = value;
                                memcpy(betterSolution, newSolver->getColSolution(), numberColumns*sizeof(double));
                            }
                        } else {
                            //newSolver->writeMps("bad3.mps");
			  sprintf(pumpPrint, "On closer inspection solution is not valid");
			  model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
			    << pumpPrint
			    << CoinMessageEol;
                            exitAll = true;
                            break;
                        }
                    }
                } else {
                    sprintf(pumpPrint, "Mini branch and bound did not improve solution (%.2f seconds)",
                            model_->getCurrentSeconds());
                    model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
                    << pumpPrint
                    << CoinMessageEol;
                    if (returnCode && newSolutionValue < saveValue + 1.0e-3 && nFixC + nFixC2) {
                        // may be able to do better
                        const double * lower = model_->solver()->getColLower();
                        const double * upper = model_->solver()->getColUpper();
                        for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
                            if (newSolver->isInteger(iColumn)) {
                                double value = floor(newSolution[iColumn] + 0.5);
                                newSolver->setColLower(iColumn, value);
                                newSolver->setColUpper(iColumn, value);
                            } else {
                                newSolver->setColLower(iColumn, lower[iColumn]);
                                newSolver->setColUpper(iColumn, upper[iColumn]);
                            }
                        }
                        newSolver->initialSolve();
                        if (newSolver->isProvenOptimal()) {
                            double value = newSolver->getObjValue() * newSolver->getObjSense();
                            if (value < saveValue) {
                                sprintf(pumpPrint, "Freeing continuous variables gives a solution of %g", value);
                                model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
                                << pumpPrint
                                << CoinMessageEol;
                                newSolutionValue = value;
                                memcpy(betterSolution, newSolver->getColSolution(), numberColumns*sizeof(double));
                                gotSolution = true;
                            }
                        }
                    }
                }
                if (gotSolution) {
                    if ((accumulate_&1) != 0) {
                        model_->incrementUsed(betterSolution); // for local search
                    }
                    solutionValue = newSolutionValue;
                    solutionFound = true;
		    numberSolutions++;
		    if (numberSolutions>=maxSolutions)
		      exitAll = true;
                    if (exitNow(newSolutionValue))
                        exitAll = true;
                    CoinWarmStartBasis * basis =
                        dynamic_cast<CoinWarmStartBasis *>(newSolver->getWarmStart()) ;
                    if (basis) {
                        bestBasis = * basis;
                        delete basis;
                        int action = model_->dealWithEventHandler(CbcEventHandler::heuristicSolution, newSolutionValue, betterSolution);
                        if (action == 0) {
                            double * saveOldSolution = CoinCopyOfArray(model_->bestSolution(), numberColumns);
                            double saveObjectiveValue = model_->getMinimizationObjValue();
                            model_->setBestSolution(betterSolution, numberColumns, newSolutionValue);
                            if (saveOldSolution && saveObjectiveValue < model_->getMinimizationObjValue())
                                model_->setBestSolution(saveOldSolution, numberColumns, saveObjectiveValue);
                            delete [] saveOldSolution;
                        }
                        if (!action || model_->getCurrentSeconds() > model_->getMaximumSeconds()) {
                            exitAll = true; // exit
                            break;
                        }
                    }
                }
            } // end stopBAB while
            delete newSolver;
        }
        if (solutionFound) finalReturnCode = 1;
        cutoff = CoinMin(cutoff, solutionValue - model_->getCutoffIncrement());
        if (numberTries >= maximumRetries_ || !solutionFound || exitAll || cutoff < continuousObjectiveValue + 1.0e-7) {
            break;
        } else {
            solutionFound = false;
            if (absoluteIncrement_ > 0.0 || relativeIncrement_ > 0.0) {
                double gap = relativeIncrement_ * fabs(solutionValue);
                double change = CoinMax(gap, absoluteIncrement_);
                cutoff = CoinMin(cutoff, solutionValue - change);
            } else {
                //double weights[10]={0.1,0.1,0.2,0.2,0.2,0.3,0.3,0.3,0.4,0.5};
                double weights[10] = {0.1, 0.2, 0.3, 0.3, 0.4, 0.4, 0.4, 0.5, 0.5, 0.6};
                cutoff -= weights[CoinMin(numberTries-1, 9)] * (cutoff - continuousObjectiveValue);
            }
            // But round down
            if (exactMultiple)
                cutoff = exactMultiple * floor(cutoff / exactMultiple);
            if (cutoff < continuousObjectiveValue)
                break;
            sprintf(pumpPrint, "Round again with cutoff of %g", cutoff);
            model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
            << pumpPrint
            << CoinMessageEol;
            if ((accumulate_&3) < 2 && usedColumn) {
                for (int i = 0; i < numberColumns; i++)
                    usedColumn[i] = -1;
            }
            averageIterationsPerTry = totalNumberIterations / numberTries;
            numberIterationsLastPass =  totalNumberIterations;
            totalNumberPasses += numberPasses - 1;
        }
    }
/*
  End of the `exitAll' loop.
*/
#ifdef RAND_RAND
    delete [] randomFactor;
#endif
    delete solver; // probably NULL but do anyway
    if (!finalReturnCode && closestSolution && closestObjectiveValue <= 10.0 &&
            usedColumn && !model_->maximumSecondsReached()) {
        // try a bit of branch and bound
        OsiSolverInterface * newSolver = cloneBut(1); // was model_->continuousSolver()->clone();
        const double * colLower = newSolver->getColLower();
        const double * colUpper = newSolver->getColUpper();
        int i;
        double rhs = 0.0;
        for (i = 0; i < numberIntegersOrig; i++) {
            int iColumn = integerVariableOrig[i];
            int direction = closestSolution[i];
            closestSolution[i] = iColumn;
            if (direction == 0) {
                // keep close to LB
                rhs += colLower[iColumn];
                lastSolution[i] = 1.0;
            } else {
                // keep close to UB
                rhs -= colUpper[iColumn];
                lastSolution[i] = -1.0;
            }
        }
        newSolver->addRow(numberIntegersOrig, closestSolution,
                          lastSolution, -COIN_DBL_MAX, rhs + 10.0);
        //double saveValue = newSolutionValue;
        //newSolver->writeMps("sub");
        int returnCode = smallBranchAndBound(newSolver, numberNodes_, newSolution, newSolutionValue,
                                             newSolutionValue, "CbcHeuristicLocalAfterFPump");
        if (returnCode < 0)
            returnCode = 0; // returned on size
        if ((returnCode&2) != 0) {
            // could add cut
            returnCode &= ~2;
        }
        if (returnCode) {
            //printf("old sol of %g new of %g\n",saveValue,newSolutionValue);
            memcpy(betterSolution, newSolution, numberColumns*sizeof(double));
            //abort();
            solutionValue = newSolutionValue;
            solutionFound = true;
	    numberSolutions++;
	    if (numberSolutions>=maxSolutions)
	      exitAll = true;
            if (exitNow(newSolutionValue))
                exitAll = true;
        }
        delete newSolver;
    }
    delete clonedSolver;
    delete [] roundingSolution;
    delete [] usedColumn;
    delete [] lastSolution;
    delete [] newSolution;
    delete [] closestSolution;
    delete [] integerVariable;
    delete [] firstPerturbedObjective;
    delete [] firstPerturbedSolution;
    if (solutionValue == incomingObjective)
        sprintf(pumpPrint, "After %.2f seconds - Feasibility pump exiting - took %.2f seconds",
                model_->getCurrentSeconds(), CoinCpuTime() - time1);
    else if (numberSolutions < maxSolutions)
      sprintf(pumpPrint, "After %.2f seconds - Feasibility pump exiting with objective of %g - took %.2f seconds",
	      model_->getCurrentSeconds(), solutionValue, CoinCpuTime() - time1);
    else 
        sprintf(pumpPrint, "After %.2f seconds - Feasibility pump exiting with objective of %g (stopping after %d solutions) - took %.2f seconds",
                model_->getCurrentSeconds(), solutionValue, 
		numberSolutions,CoinCpuTime() - time1);
    model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
      << pumpPrint
      << CoinMessageEol;
    if (bestBasis.getNumStructural())
        model_->setBestSolutionBasis(bestBasis);
    //model_->setMinimizationObjValue(saveBestObjective);
    if ((accumulate_&1) != 0 && numberSolutions > 1 && !model_->getSolutionCount()) {
        model_->setSolutionCount(1); // for local search
        model_->setNumberHeuristicSolutions(1);
    }
#ifdef COIN_DEVELOP
    {
        double ncol = model_->solver()->getNumCols();
        double nrow = model_->solver()->getNumRows();
        printf("XXX total iterations %d ratios - %g %g %g\n",
               totalNumberIterations,
               static_cast<double> (totalNumberIterations) / nrow,
               static_cast<double> (totalNumberIterations) / ncol,
               static_cast<double> (totalNumberIterations) / (2*nrow + 2*ncol));
    }
#endif
    return finalReturnCode;
}

/**************************END MAIN PROCEDURE ***********************************/

// update model
void CbcHeuristicFPump::setModel(CbcModel * model)
{
    model_ = model;
}

/* Rounds solution - down if < downValue
   returns 1 if current is a feasible solution
*/
int
CbcHeuristicFPump::rounds(OsiSolverInterface * solver, double * solution,
                          //const double * objective,
                          int numberIntegers, const int * integerVariable,
                          /*char * pumpPrint,*/ int iter,
                          /*bool roundExpensive,*/ double downValue, int *flip)
{
    double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
    double primalTolerance ;
    solver->getDblParam(OsiPrimalTolerance, primalTolerance) ;

    int i;

    const double * cost = solver->getObjCoefficients();
    int flip_up = 0;
    int flip_down  = 0;
    double  v = randomNumberGenerator_.randomDouble() * 20.0;
    int nn = 10 + static_cast<int> (v);
    int nnv = 0;
    int * list = new int [nn];
    double * val = new double [nn];
    for (i = 0; i < nn; i++) val[i] = .001;

    const double * rowLower = solver->getRowLower();
    const double * rowUpper = solver->getRowUpper();
    int numberRows = solver->getNumRows();
    if (false && (iter&1) != 0) {
        // Do set covering variables
        const CoinPackedMatrix * matrixByRow = solver->getMatrixByRow();
        const double * elementByRow = matrixByRow->getElements();
        const int * column = matrixByRow->getIndices();
        const CoinBigIndex * rowStart = matrixByRow->getVectorStarts();
        const int * rowLength = matrixByRow->getVectorLengths();
        for (i = 0; i < numberRows; i++) {
            if (rowLower[i] == 1.0 && rowUpper[i] == 1.0) {
                bool cover = true;
                double largest = 0.0;
                int jColumn = -1;
                for (CoinBigIndex k = rowStart[i]; k < rowStart[i] + rowLength[i]; k++) {
                    int iColumn = column[k];
                    if (elementByRow[k] != 1.0 || !solver->isInteger(iColumn)) {
                        cover = false;
                        break;
                    } else {
                        if (solution[iColumn]) {
                            double value = solution[iColumn] *
                                           (randomNumberGenerator_.randomDouble() + 5.0);
                            if (value > largest) {
                                largest = value;
                                jColumn = iColumn;
                            }
                        }
                    }
                }
                if (cover) {
                    for (CoinBigIndex k = rowStart[i]; k < rowStart[i] + rowLength[i]; k++) {
                        int iColumn = column[k];
                        if (iColumn == jColumn)
                            solution[iColumn] = 1.0;
                        else
                            solution[iColumn] = 0.0;
                    }
                }
            }
        }
    }
    int numberColumns = solver->getNumCols();
#ifdef JJF_ZERO
    // Do set covering variables
    const CoinPackedMatrix * matrixByRow = solver->getMatrixByRow();
    const double * elementByRow = matrixByRow->getElements();
    const int * column = matrixByRow->getIndices();
    const CoinBigIndex * rowStart = matrixByRow->getVectorStarts();
    const int * rowLength = matrixByRow->getVectorLengths();
    double * sortTemp = new double[numberColumns];
    int * whichTemp = new int [numberColumns];
    char * rowTemp = new char [numberRows];
    memset(rowTemp, 0, numberRows);
    for (i = 0; i < numberColumns; i++)
        whichTemp[i] = -1;
    int nSOS = 0;
    for (i = 0; i < numberRows; i++) {
        if (rowLower[i] == 1.0 && rowUpper[i] == 1.0) {
            bool cover = true;
            for (CoinBigIndex k = rowStart[i]; k < rowStart[i] + rowLength[i]; k++) {
                int iColumn = column[k];
                if (elementByRow[k] != 1.0 || !solver->isInteger(iColumn)) {
                    cover = false;
                    break;
                }
            }
            if (cover) {
                rowTemp[i] = 1;
                nSOS++;
                for (CoinBigIndex k = rowStart[i]; k < rowStart[i] + rowLength[i]; k++) {
                    int iColumn = column[k];
                    double value = solution[iColumn];
                    whichTemp[iColumn] = iColumn;
                }
            }
        }
    }
    if (nSOS) {
        // Column copy
        const CoinPackedMatrix * matrixByColumn = solver->getMatrixByCol();
        //const double * element = matrixByColumn->getElements();
        const int * row = matrixByColumn->getIndices();
        const CoinBigIndex * columnStart = matrixByColumn->getVectorStarts();
        const int * columnLength = matrixByColumn->getVectorLengths();
        int nLook = 0;
        for (i = 0; i < numberColumns; i++) {
            if (whichTemp[i] >= 0) {
                whichTemp[nLook] = i;
                double value = solution[i];
                if (value < 0.5)
                    value *= (0.1 * randomNumberGenerator_.randomDouble() + 0.3);
                sortTemp[nLook++] = -value;
            }
        }
        CoinSort_2(sortTemp, sortTemp + nLook, whichTemp);
        double smallest = 1.0;
        int nFix = 0;
        int nOne = 0;
        for (int j = 0; j < nLook; j++) {
            int jColumn = whichTemp[j];
            double thisValue = solution[jColumn];
            if (!thisValue)
                continue;
            if (thisValue == 1.0)
                nOne++;
            smallest = CoinMin(smallest, thisValue);
            solution[jColumn] = 1.0;
            double largest = 0.0;
            for (CoinBigIndex jEl = columnStart[jColumn];
                    jEl < columnStart[jColumn] + columnLength[jColumn]; jEl++) {
                int jRow = row[jEl];
                if (rowTemp[jRow]) {
                    for (CoinBigIndex k = rowStart[jRow]; k < rowStart[jRow] + rowLength[jRow]; k++) {
                        int iColumn = column[k];
                        if (solution[iColumn]) {
                            if (iColumn != jColumn) {
                                double value = solution[iColumn];
                                if (value > largest)
                                    largest = value;
                                solution[iColumn] = 0.0;
                            }
                        }
                    }
                }
            }
            if (largest > thisValue)
                printf("%d was at %g - chosen over a value of %g\n",
                       jColumn, thisValue, largest);
            nFix++;
        }
        printf("%d fixed out of %d (%d at one already)\n",
               nFix, nLook, nOne);
    }
    delete [] sortTemp;
    delete [] whichTemp;
    delete [] rowTemp;
#endif
    const double * columnLower = solver->getColLower();
    const double * columnUpper = solver->getColUpper();
    // Check if valid with current solution (allow for 0.99999999s)
    for (i = 0; i < numberIntegers; i++) {
        int iColumn = integerVariable[i];
        double value = solution[iColumn];
        double round = floor(value + 0.5);
        if (fabs(value - round) > primalTolerance)
            break;
    }
    if (i == numberIntegers) {
        // may be able to use solution even if 0.99999's
        double * saveLower = CoinCopyOfArray(columnLower, numberColumns);
        double * saveUpper = CoinCopyOfArray(columnUpper, numberColumns);
        double * saveSolution = CoinCopyOfArray(solution, numberColumns);
        double * tempSolution = CoinCopyOfArray(solution, numberColumns);
        CoinWarmStartBasis  * saveBasis =
            dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ;
        for (i = 0; i < numberIntegers; i++) {
            int iColumn = integerVariable[i];
            double value = solution[iColumn];
            double round = floor(value + 0.5);
            solver->setColLower(iColumn, round);
            solver->setColUpper(iColumn, round);
            tempSolution[iColumn] = round;
        }
        solver->setColSolution(tempSolution);
        delete [] tempSolution;
        solver->resolve();
        solver->setColLower(saveLower);
        solver->setColUpper(saveUpper);
        solver->setWarmStart(saveBasis);
        delete [] saveLower;
        delete [] saveUpper;
        delete saveBasis;
        if (!solver->isProvenOptimal()) {
            solver->setColSolution(saveSolution);
        }
        delete [] saveSolution;
        if (solver->isProvenOptimal()) {
            // feasible
            delete [] list;
            delete [] val;
            return 1;
        }
    }
    //double * saveSolution = CoinCopyOfArray(solution,numberColumns);
    // return rounded solution
    for (i = 0; i < numberIntegers; i++) {
        int iColumn = integerVariable[i];
        double value = solution[iColumn];
        double round = floor(value + primalTolerance);
        if (value - round > downValue) round += 1.;
#ifndef JJF_ONE
        if (round < integerTolerance && cost[iColumn] < -1. + integerTolerance) flip_down++;
        if (round > 1. - integerTolerance && cost[iColumn] > 1. - integerTolerance) flip_up++;
#else
        if (round < columnLower[iColumn] + integerTolerance && cost[iColumn] < -1. + integerTolerance) flip_down++;
        if (round > columnUpper[iColumn] - integerTolerance && cost[iColumn] > 1. - integerTolerance) flip_up++;
#endif
        if (flip_up + flip_down == 0) {
            for (int k = 0; k < nn; k++) {
                if (fabs(value - round) > val[k]) {
                    nnv++;
                    for (int j = nn - 2; j >= k; j--) {
                        val[j+1] = val[j];
                        list[j+1] = list[j];
                    }
                    val[k] = fabs(value - round);
                    list[k] = iColumn;
                    break;
                }
            }
        }
        solution[iColumn] = round;
    }

    if (nnv > nn) nnv = nn;
    //if (iter != 0)
    //sprintf(pumpPrint+strlen(pumpPrint),"up = %5d , down = %5d", flip_up, flip_down);
    *flip = flip_up + flip_down;

    if (*flip == 0 && iter != 0) {
        //sprintf(pumpPrint+strlen(pumpPrint)," -- rand = %4d (%4d) ", nnv, nn);
        for (i = 0; i < nnv; i++) {
            // was solution[list[i]] = 1. - solution[list[i]]; but does that work for 7>=x>=6
            int index = list[i];
            double value = solution[index];
            if (value <= 1.0) {
                solution[index] = 1.0 - value;
            } else if (value < columnLower[index] + integerTolerance) {
                solution[index] = value + 1.0;
            } else if (value > columnUpper[index] - integerTolerance) {
                solution[index] = value - 1.0;
            } else {
                solution[index] = value - 1.0;
            }
        }
        *flip = nnv;
    } else {
        //sprintf(pumpPrint+strlen(pumpPrint)," ");
    }
    delete [] list;
    delete [] val;
    //iter++;

    // get row activities
    double * rowActivity = new double[numberRows];
    memset(rowActivity, 0, numberRows*sizeof(double));
    solver->getMatrixByCol()->times(solution, rowActivity) ;
    double largestInfeasibility = primalTolerance;
    double sumInfeasibility = 0.0;
    int numberBadRows = 0;
    for (i = 0 ; i < numberRows ; i++) {
        double value;
        value = rowLower[i] - rowActivity[i];
        if (value > primalTolerance) {
            numberBadRows++;
            largestInfeasibility = CoinMax(largestInfeasibility, value);
            sumInfeasibility += value;
        }
        value = rowActivity[i] - rowUpper[i];
        if (value > primalTolerance) {
            numberBadRows++;
            largestInfeasibility = CoinMax(largestInfeasibility, value);
            sumInfeasibility += value;
        }
    }
#ifdef JJF_ZERO
    if (largestInfeasibility > primalTolerance && numberBadRows*10 < numberRows) {
        // Can we improve by flipping
        for (int iPass = 0; iPass < 10; iPass++) {
            int numberColumns = solver->getNumCols();
            const CoinPackedMatrix * matrixByCol = solver->getMatrixByCol();
            const double * element = matrixByCol->getElements();
            const int * row = matrixByCol->getIndices();
            const CoinBigIndex * columnStart = matrixByCol->getVectorStarts();
            const int * columnLength = matrixByCol->getVectorLengths();
            double oldSum = sumInfeasibility;
            // First improve by moving continuous ones
            for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
                if (!solver->isInteger(iColumn)) {
                    double solValue = solution[iColumn];
                    double thetaUp = columnUpper[iColumn] - solValue;
                    double improvementUp = 0.0;
                    if (thetaUp > primalTolerance) {
                        // can go up
                        for (CoinBigIndex j = columnStart[iColumn];
                                j < columnStart[iColumn] + columnLength[iColumn]; j++) {
                            int iRow = row[j];
                            double distanceUp = rowUpper[iRow] - rowActivity[iRow];
                            double distanceDown = rowLower[iRow] - rowActivity[iRow];
                            double el = element[j];
                            if (el > 0.0) {
                                // positive element
                                if (distanceUp > 0.0) {
                                    if (thetaUp*el > distanceUp)
                                        thetaUp = distanceUp / el;
                                } else {
                                    improvementUp -= el;
                                }
                                if (distanceDown > 0.0) {
                                    if (thetaUp*el > distanceDown)
                                        thetaUp = distanceDown / el;
                                    improvementUp += el;
                                }
                            } else {
                                // negative element
                                if (distanceDown < 0.0) {
                                    if (thetaUp*el < distanceDown)
                                        thetaUp = distanceDown / el;
                                } else {
                                    improvementUp += el;
                                }
                                if (distanceUp < 0.0) {
                                    if (thetaUp*el < distanceUp)
                                        thetaUp = distanceUp / el;
                                    improvementUp -= el;
                                }
                            }
                        }
                    }
                    double thetaDown = solValue - columnLower[iColumn];
                    double improvementDown = 0.0;
                    if (thetaDown > primalTolerance) {
                        // can go down
                        for (CoinBigIndex j = columnStart[iColumn];
                                j < columnStart[iColumn] + columnLength[iColumn]; j++) {
                            int iRow = row[j];
                            double distanceUp = rowUpper[iRow] - rowActivity[iRow];
                            double distanceDown = rowLower[iRow] - rowActivity[iRow];
                            double el = -element[j]; // not change in sign form up
                            if (el > 0.0) {
                                // positive element
                                if (distanceUp > 0.0) {
                                    if (thetaDown*el > distanceUp)
                                        thetaDown = distanceUp / el;
                                } else {
                                    improvementDown -= el;
                                }
                                if (distanceDown > 0.0) {
                                    if (thetaDown*el > distanceDown)
                                        thetaDown = distanceDown / el;
                                    improvementDown += el;
                                }
                            } else {
                                // negative element
                                if (distanceDown < 0.0) {
                                    if (thetaDown*el < distanceDown)
                                        thetaDown = distanceDown / el;
                                } else {
                                    improvementDown += el;
                                }
                                if (distanceUp < 0.0) {
                                    if (thetaDown*el < distanceUp)
                                        thetaDown = distanceUp / el;
                                    improvementDown -= el;
                                }
                            }
                        }
                        if (thetaUp < 1.0e-8)
                            improvementUp = 0.0;
                        if (thetaDown < 1.0e-8)
                            improvementDown = 0.0;
                        double theta;
                        if (improvementUp >= improvementDown) {
                            theta = thetaUp;
                        } else {
                            improvementUp = improvementDown;
                            theta = -thetaDown;
                        }
                        if (improvementUp > 1.0e-8 && fabs(theta) > 1.0e-8) {
                            // Could move
                            double oldSum = 0.0;
                            double newSum = 0.0;
                            solution[iColumn] += theta;
                            for (CoinBigIndex j = columnStart[iColumn];
                                    j < columnStart[iColumn] + columnLength[iColumn]; j++) {
                                int iRow = row[j];
                                double lower = rowLower[iRow];
                                double upper = rowUpper[iRow];
                                double value = rowActivity[iRow];
                                if (value > upper)
                                    oldSum += value - upper;
                                else if (value < lower)
                                    oldSum += lower - value;
                                value += theta * element[j];
                                rowActivity[iRow] = value;
                                if (value > upper)
                                    newSum += value - upper;
                                else if (value < lower)
                                    newSum += lower - value;
                            }
                            assert (newSum <= oldSum);
                            sumInfeasibility += newSum - oldSum;
                        }
                    }
                }
            }
            // Now flip some integers?
#ifdef JJF_ZERO
            for (i = 0; i < numberIntegers; i++) {
                int iColumn = integerVariable[i];
                double solValue = solution[iColumn];
                assert (fabs(solValue - floor(solValue + 0.5)) < 1.0e-8);
                double improvementUp = 0.0;
                if (columnUpper[iColumn] >= solValue + 1.0) {
                    // can go up
                    double oldSum = 0.0;
                    double newSum = 0.0;
                    for (CoinBigIndex j = columnStart[iColumn];
                            j < columnStart[iColumn] + columnLength[iColumn]; j++) {
                        int iRow = row[j];
                        double lower = rowLower[iRow];
                        double upper = rowUpper[iRow];
                        double value = rowActivity[iRow];
                        if (value > upper)
                            oldSum += value - upper;
                        else if (value < lower)
                            oldSum += lower - value;
                        value += element[j];
                        if (value > upper)
                            newSum += value - upper;
                        else if (value < lower)
                            newSum += lower - value;
                    }
                    improvementUp = oldSum - newSum;
                }
                double improvementDown = 0.0;
                if (columnLower[iColumn] <= solValue - 1.0) {
                    // can go down
                    double oldSum = 0.0;
                    double newSum = 0.0;
                    for (CoinBigIndex j = columnStart[iColumn];
                            j < columnStart[iColumn] + columnLength[iColumn]; j++) {
                        int iRow = row[j];
                        double lower = rowLower[iRow];
                        double upper = rowUpper[iRow];
                        double value = rowActivity[iRow];
                        if (value > upper)
                            oldSum += value - upper;
                        else if (value < lower)
                            oldSum += lower - value;
                        value -= element[j];
                        if (value > upper)
                            newSum += value - upper;
                        else if (value < lower)
                            newSum += lower - value;
                    }
                    improvementDown = oldSum - newSum;
                }
                double theta;
                if (improvementUp >= improvementDown) {
                    theta = 1.0;
                } else {
                    improvementUp = improvementDown;
                    theta = -1.0;
                }
                if (improvementUp > 1.0e-8 && fabs(theta) > 1.0e-8) {
                    // Could move
                    double oldSum = 0.0;
                    double newSum = 0.0;
                    solution[iColumn] += theta;
                    for (CoinBigIndex j = columnStart[iColumn];
                            j < columnStart[iColumn] + columnLength[iColumn]; j++) {
                        int iRow = row[j];
                        double lower = rowLower[iRow];
                        double upper = rowUpper[iRow];
                        double value = rowActivity[iRow];
                        if (value > upper)
                            oldSum += value - upper;
                        else if (value < lower)
                            oldSum += lower - value;
                        value += theta * element[j];
                        rowActivity[iRow] = value;
                        if (value > upper)
                            newSum += value - upper;
                        else if (value < lower)
                            newSum += lower - value;
                    }
                    assert (newSum <= oldSum);
                    sumInfeasibility += newSum - oldSum;
                }
            }
#else
            int bestColumn = -1;
            double bestImprovement = primalTolerance;
            double theta = 0.0;
            for (i = 0; i < numberIntegers; i++) {
                int iColumn = integerVariable[i];
                double solValue = solution[iColumn];
                assert (fabs(solValue - floor(solValue + 0.5)) < 1.0e-8);
                double improvementUp = 0.0;
                if (columnUpper[iColumn] >= solValue + 1.0) {
                    // can go up
                    double oldSum = 0.0;
                    double newSum = 0.0;
                    for (CoinBigIndex j = columnStart[iColumn];
                            j < columnStart[iColumn] + columnLength[iColumn]; j++) {
                        int iRow = row[j];
                        double lower = rowLower[iRow];
                        double upper = rowUpper[iRow];
                        double value = rowActivity[iRow];
                        if (value > upper)
                            oldSum += value - upper;
                        else if (value < lower)
                            oldSum += lower - value;
                        value += element[j];
                        if (value > upper)
                            newSum += value - upper;
                        else if (value < lower)
                            newSum += lower - value;
                    }
                    improvementUp = oldSum - newSum;
                }
                double improvementDown = 0.0;
                if (columnLower[iColumn] <= solValue - 1.0) {
                    // can go down
                    double oldSum = 0.0;
                    double newSum = 0.0;
                    for (CoinBigIndex j = columnStart[iColumn];
                            j < columnStart[iColumn] + columnLength[iColumn]; j++) {
                        int iRow = row[j];
                        double lower = rowLower[iRow];
                        double upper = rowUpper[iRow];
                        double value = rowActivity[iRow];
                        if (value > upper)
                            oldSum += value - upper;
                        else if (value < lower)
                            oldSum += lower - value;
                        value -= element[j];
                        if (value > upper)
                            newSum += value - upper;
                        else if (value < lower)
                            newSum += lower - value;
                    }
                    improvementDown = oldSum - newSum;
                }
                double improvement = CoinMax(improvementUp, improvementDown);
                if (improvement > bestImprovement) {
                    bestImprovement = improvement;
                    bestColumn = iColumn;
                    if (improvementUp > improvementDown)
                        theta = 1.0;
                    else
                        theta = -1.0;
                }
            }
            if (bestColumn >= 0) {
                // Could move
                int iColumn = bestColumn;
                double oldSum = 0.0;
                double newSum = 0.0;
                solution[iColumn] += theta;
                for (CoinBigIndex j = columnStart[iColumn];
                        j < columnStart[iColumn] + columnLength[iColumn]; j++) {
                    int iRow = row[j];
                    double lower = rowLower[iRow];
                    double upper = rowUpper[iRow];
                    double value = rowActivity[iRow];
                    if (value > upper)
                        oldSum += value - upper;
                    else if (value < lower)
                        oldSum += lower - value;
                    value += theta * element[j];
                    rowActivity[iRow] = value;
                    if (value > upper)
                        newSum += value - upper;
                    else if (value < lower)
                        newSum += lower - value;
                }
                assert (newSum <= oldSum);
                sumInfeasibility += newSum - oldSum;
            }
#endif
            if (oldSum <= sumInfeasibility + primalTolerance)
                break; // no good
        }
    }
    //delete [] saveSolution;
#endif
    delete [] rowActivity;
    return (largestInfeasibility > primalTolerance) ? 0 : 1;
}
// Set maximum Time (default off) - also sets starttime to current
void
CbcHeuristicFPump::setMaximumTime(double value)
{
    startTime_ = CoinCpuTime();
    maximumTime_ = value;
}

# ifdef COIN_HAS_CLP

//#############################################################################
// Constructors / Destructor / Assignment
//#############################################################################

//-------------------------------------------------------------------
// Default Constructor
//-------------------------------------------------------------------
CbcDisasterHandler::CbcDisasterHandler (CbcModel * model)
        : OsiClpDisasterHandler(),
        cbcModel_(model)
{
    if (model) {
        osiModel_
        = dynamic_cast<OsiClpSolverInterface *> (model->solver());
        if (osiModel_)
            setSimplex(osiModel_->getModelPtr());
    }
}

//-------------------------------------------------------------------
// Copy constructor
//-------------------------------------------------------------------
CbcDisasterHandler::CbcDisasterHandler (const CbcDisasterHandler & rhs)
        : OsiClpDisasterHandler(rhs),
        cbcModel_(rhs.cbcModel_)
{
}


//-------------------------------------------------------------------
// Destructor
//-------------------------------------------------------------------
CbcDisasterHandler::~CbcDisasterHandler ()
{
}

//----------------------------------------------------------------
// Assignment operator
//-------------------------------------------------------------------
CbcDisasterHandler &
CbcDisasterHandler::operator=(const CbcDisasterHandler & rhs)
{
    if (this != &rhs) {
        OsiClpDisasterHandler::operator=(rhs);
        cbcModel_ = rhs.cbcModel_;
    }
    return *this;
}
//-------------------------------------------------------------------
// Clone
//-------------------------------------------------------------------
ClpDisasterHandler * CbcDisasterHandler::clone() const
{
    return new CbcDisasterHandler(*this);
}
// Type of disaster 0 can fix, 1 abort
int
CbcDisasterHandler::typeOfDisaster()
{
    if (!cbcModel_->parentModel() &&
            (cbcModel_->specialOptions()&2048) == 0) {
        return 0;
    } else {
        if (cbcModel_->parentModel())
            cbcModel_->setMaximumNodes(0);
        return 1;
    }
}
/* set model. */
void
CbcDisasterHandler::setCbcModel(CbcModel * model)
{
    cbcModel_ = model;
    if (model) {
        osiModel_
        = dynamic_cast<OsiClpSolverInterface *> (model->solver());
        if (osiModel_)
            setSimplex(osiModel_->getModelPtr());
        else
            setSimplex(NULL);
    }
}
#endif