/* $Id: CbcHeuristic.cpp 1888 2013-04-06 20:52:59Z stefan $ */
// Copyright (C) 2002, 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 "CbcConfig.h"

#include <cassert>
#include <cstdlib>
#include <cmath>
#include <cfloat>

//#define PRINT_DEBUG
#ifdef COIN_HAS_CLP
#include "OsiClpSolverInterface.hpp"
#endif
#include "CbcModel.hpp"
#include "CbcMessage.hpp"
#include "CbcHeuristic.hpp"
#include "CbcHeuristicFPump.hpp"
#include "CbcStrategy.hpp"
#include "CglPreProcess.hpp"
#include "CglGomory.hpp"
#include "CglProbing.hpp"
#include "OsiAuxInfo.hpp"
#include "OsiPresolve.hpp"
#include "CbcBranchActual.hpp"
#include "CbcCutGenerator.hpp"
//==============================================================================

CbcHeuristicNode::CbcHeuristicNode(const CbcHeuristicNode& rhs)
{
    numObjects_ = rhs.numObjects_;
    brObj_ = new CbcBranchingObject*[numObjects_];
    for (int i = 0; i < numObjects_; ++i) {
        brObj_[i] = rhs.brObj_[i]->clone();
    }
}

void
CbcHeuristicNodeList::gutsOfDelete()
{
    for (int i = (static_cast<int>(nodes_.size())) - 1; i >= 0; --i) {
        delete nodes_[i];
    }
}

void
CbcHeuristicNodeList::gutsOfCopy(const CbcHeuristicNodeList& rhs)
{
    append(rhs);
}

CbcHeuristicNodeList::CbcHeuristicNodeList(const CbcHeuristicNodeList& rhs)
{
    gutsOfCopy(rhs);
}

CbcHeuristicNodeList& CbcHeuristicNodeList::operator=
(const CbcHeuristicNodeList & rhs)
{
    if (this != &rhs) {
        gutsOfDelete();
        gutsOfCopy(rhs);
    }
    return *this;
}

CbcHeuristicNodeList::~CbcHeuristicNodeList()
{
    gutsOfDelete();
}

void
CbcHeuristicNodeList::append(CbcHeuristicNode*& node)
{
    nodes_.push_back(node);
    node = NULL;
}

void
CbcHeuristicNodeList::append(const CbcHeuristicNodeList& nodes)
{
    nodes_.reserve(nodes_.size() + nodes.size());
    for (int i = 0; i < nodes.size(); ++i) {
        CbcHeuristicNode* node = new CbcHeuristicNode(*nodes.node(i));
        append(node);
    }
}

//==============================================================================
#define DEFAULT_WHERE ((255-2-16)*(1+256))
// Default Constructor
CbcHeuristic::CbcHeuristic() :
        model_(NULL),
        when_(2),
        numberNodes_(200),
        feasibilityPumpOptions_(-1),
        fractionSmall_(1.0),
        heuristicName_("Unknown"),
        howOften_(1),
        decayFactor_(0.0),
        switches_(0),
        whereFrom_(DEFAULT_WHERE),
        shallowDepth_(1),
        howOftenShallow_(1),
        numInvocationsInShallow_(0),
        numInvocationsInDeep_(0),
        lastRunDeep_(0),
        numRuns_(0),
        minDistanceToRun_(1),
        runNodes_(),
        numCouldRun_(0),
        numberSolutionsFound_(0),
        numberNodesDone_(0),
        inputSolution_(NULL)
{
    // As CbcHeuristic virtual need to modify .cpp if above change
}

// Constructor from model
CbcHeuristic::CbcHeuristic(CbcModel & model) :
        model_(&model),
        when_(2),
        numberNodes_(200),
        feasibilityPumpOptions_(-1),
        fractionSmall_(1.0),
        heuristicName_("Unknown"),
        howOften_(1),
        decayFactor_(0.0),
        switches_(0),
        whereFrom_(DEFAULT_WHERE),
        shallowDepth_(1),
        howOftenShallow_(1),
        numInvocationsInShallow_(0),
        numInvocationsInDeep_(0),
        lastRunDeep_(0),
        numRuns_(0),
        minDistanceToRun_(1),
        runNodes_(),
        numCouldRun_(0),
        numberSolutionsFound_(0),
        numberNodesDone_(0),
        inputSolution_(NULL)
{}

void
CbcHeuristic::gutsOfCopy(const CbcHeuristic & rhs)
{
    model_ = rhs.model_;
    when_ = rhs.when_;
    numberNodes_ = rhs.numberNodes_;
    feasibilityPumpOptions_ = rhs.feasibilityPumpOptions_;
    fractionSmall_ = rhs.fractionSmall_;
    randomNumberGenerator_ = rhs.randomNumberGenerator_;
    heuristicName_ = rhs.heuristicName_;
    howOften_ = rhs.howOften_;
    decayFactor_ = rhs.decayFactor_;
    switches_ = rhs.switches_;
    whereFrom_ = rhs.whereFrom_;
    shallowDepth_ = rhs.shallowDepth_;
    howOftenShallow_ = rhs.howOftenShallow_;
    numInvocationsInShallow_ = rhs.numInvocationsInShallow_;
    numInvocationsInDeep_ = rhs.numInvocationsInDeep_;
    lastRunDeep_ = rhs.lastRunDeep_;
    numRuns_ = rhs.numRuns_;
    numCouldRun_ = rhs.numCouldRun_;
    minDistanceToRun_ = rhs.minDistanceToRun_;
    runNodes_ = rhs.runNodes_;
    numberSolutionsFound_ = rhs.numberSolutionsFound_;
    numberNodesDone_ = rhs.numberNodesDone_;
    if (rhs.inputSolution_) {
        int numberColumns = model_->getNumCols();
        setInputSolution(rhs.inputSolution_, rhs.inputSolution_[numberColumns]);
    }
}
// Copy constructor
CbcHeuristic::CbcHeuristic(const CbcHeuristic & rhs)
{
    inputSolution_ = NULL;
    gutsOfCopy(rhs);
}

// Assignment operator
CbcHeuristic &
CbcHeuristic::operator=( const CbcHeuristic & rhs)
{
    if (this != &rhs) {
        gutsOfDelete();
        gutsOfCopy(rhs);
    }
    return *this;
}

void CbcHeurDebugNodes(CbcModel* model_)
{
    CbcNode* node = model_->currentNode();
    CbcNodeInfo* nodeInfo = node->nodeInfo();
    std::cout << "===============================================================\n";
    while (nodeInfo) {
        const CbcNode* node = nodeInfo->owner();
        printf("nodeinfo: node %i\n", nodeInfo->nodeNumber());
        {
            const CbcIntegerBranchingObject* brPrint =
                dynamic_cast<const CbcIntegerBranchingObject*>(nodeInfo->parentBranch());
            if (!brPrint) {
                printf("    parentBranch: NULL\n");
            } else {
                const double* downBounds = brPrint->downBounds();
                const double* upBounds = brPrint->upBounds();
                int variable = brPrint->variable();
                int way = brPrint->way();
                printf("   parentBranch: var %i downBd [%i,%i] upBd [%i,%i] way %i\n",
                       variable, static_cast<int>(downBounds[0]), static_cast<int>(downBounds[1]),
                       static_cast<int>(upBounds[0]), static_cast<int>(upBounds[1]), way);
            }
        }
        if (! node) {
            printf("    owner: NULL\n");
        } else {
            printf("    owner: node %i depth %i onTree %i active %i",
                   node->nodeNumber(), node->depth(), node->onTree(), node->active());
            const OsiBranchingObject* osibr =
                nodeInfo->owner()->branchingObject();
            const CbcBranchingObject* cbcbr =
                dynamic_cast<const CbcBranchingObject*>(osibr);
            const CbcIntegerBranchingObject* brPrint =
                dynamic_cast<const CbcIntegerBranchingObject*>(cbcbr);
            if (!brPrint) {
                printf("        ownerBranch: NULL\n");
            } else {
                const double* downBounds = brPrint->downBounds();
                const double* upBounds = brPrint->upBounds();
                int variable = brPrint->variable();
                int way = brPrint->way();
                printf("        ownerbranch: var %i downBd [%i,%i] upBd [%i,%i] way %i\n",
                       variable, static_cast<int>(downBounds[0]), static_cast<int>(downBounds[1]),
                       static_cast<int>(upBounds[0]), static_cast<int>(upBounds[1]), way);
            }
        }
        nodeInfo = nodeInfo->parent();
    }
}

void
CbcHeuristic::debugNodes()
{
    CbcHeurDebugNodes(model_);
}

void
CbcHeuristic::printDistanceToNodes()
{
    const CbcNode* currentNode = model_->currentNode();
    if (currentNode != NULL) {
        CbcHeuristicNode* nodeDesc = new CbcHeuristicNode(*model_);
        for (int i = runNodes_.size() - 1; i >= 0; --i) {
            nodeDesc->distance(runNodes_.node(i));
        }
        runNodes_.append(nodeDesc);
    }
}

bool
CbcHeuristic::shouldHeurRun(int whereFrom)
{
    assert (whereFrom >= 0 && whereFrom < 16);
    // take off 8 (code - likes new solution)
    whereFrom &= 7;
    if ((whereFrom_&(1 << whereFrom)) == 0)
        return false;
    // No longer used for original purpose - so use for ever run at all JJF
#ifndef JJF_ONE
    // Don't run if hot start
    if (model_ && model_->hotstartSolution())
        return false;
    else
        return true;
#else
#ifdef JJF_ZERO
    const CbcNode* currentNode = model_->currentNode();
    if (currentNode == NULL) {
        return false;
    }

    debugNodes();
//   return false;

    const int depth = currentNode->depth();
#else
    int depth = model_->currentDepth();
#endif

    const int nodeCount = model_->getNodeCount();  // FIXME: check that this is
    // correct in parallel

    if (nodeCount == 0 || depth <= shallowDepth_) {
        // what to do when we are in the shallow part of the tree
        if (model_->getCurrentPassNumber() == 1) {
            // first time in the node...
            numInvocationsInShallow_ = 0;
        }
        ++numInvocationsInShallow_;
        // Very large howOftenShallow_ will give the original test:
        // (model_->getCurrentPassNumber() != 1)
        //    if ((numInvocationsInShallow_ % howOftenShallow_) != 1) {
        if ((numInvocationsInShallow_ % howOftenShallow_) != 0) {
            return false;
        }
        // LL: should we save these nodes in the list of nodes where the heur was
        // LL: run?
#ifndef JJF_ONE
        if (currentNode != NULL) {
            // Get where we are and create the appropriate CbcHeuristicNode object
            CbcHeuristicNode* nodeDesc = new CbcHeuristicNode(*model_);
            runNodes_.append(nodeDesc);
        }
#endif
    } else {
        // deeper in the tree
        if (model_->getCurrentPassNumber() == 1) {
            // first time in the node...
            ++numInvocationsInDeep_;
        }
        if (numInvocationsInDeep_ - lastRunDeep_ < howOften_) {
            return false;
        }
        if (model_->getCurrentPassNumber() != 1) {
            // Run the heuristic only when first entering the node.
            // LL: I don't think this is right. It should run just before strong
            // LL: branching, I believe.
            return false;
        }
        // Get where we are and create the appropriate CbcHeuristicNode object
        CbcHeuristicNode* nodeDesc = new CbcHeuristicNode(*model_);
        //#ifdef PRINT_DEBUG
#ifndef JJF_ONE
        const double minDistanceToRun = 1.5 * log((double)depth) / log((double)2);
#else
    const double minDistanceToRun = minDistanceToRun_;
#endif
#ifdef PRINT_DEBUG
        double minDistance = nodeDesc->minDistance(runNodes_);
        std::cout << "minDistance = " << minDistance
                  << ", minDistanceToRun = " << minDistanceToRun << std::endl;
#endif
        if (nodeDesc->minDistanceIsSmall(runNodes_, minDistanceToRun)) {
            delete nodeDesc;
            return false;
        }
        runNodes_.append(nodeDesc);
        lastRunDeep_ = numInvocationsInDeep_;
        //    ++lastRunDeep_;
    }
    ++numRuns_;
    return true;
#endif
}

bool
CbcHeuristic::shouldHeurRun_randomChoice()
{
    if (!when_)
        return false;
    int depth = model_->currentDepth();
    // when_ -999 is special marker to force to run
    if (depth != 0 && when_ != -999) {
        const double numerator = depth * depth;
        const double denominator = exp(depth * log(2.0));
        double probability = numerator / denominator;
        double randomNumber = randomNumberGenerator_.randomDouble();
        int when = when_ % 100;
        if (when > 2 && when < 8) {
            /* JJF adjustments
            3 only at root and if no solution
            4 only at root and if this heuristic has not got solution
            5 as 3 but decay more
            6 decay
            7 run up to 2 times if solution found 4 otherwise
            */
            switch (when) {
            case 3:
            default:
                if (model_->bestSolution())
                    probability = -1.0;
                break;
            case 4:
                if (numberSolutionsFound_)
                    probability = -1.0;
                break;
            case 5:
                assert (decayFactor_);
                if (model_->bestSolution()) {
                    probability = -1.0;
                } else if (numCouldRun_ > 1000) {
                    decayFactor_ *= 0.99;
                    probability *= decayFactor_;
                }
                break;
            case 6:
                if (depth >= 3) {
                    if ((numCouldRun_ % howOften_) == 0 &&
                            numberSolutionsFound_*howOften_ < numCouldRun_) {
#ifdef COIN_DEVELOP
                        int old = howOften_;
#endif
                        howOften_ = CoinMin(CoinMax(static_cast<int> (howOften_ * 1.1), howOften_ + 1), 1000000);
#ifdef COIN_DEVELOP
                        printf("Howoften changed from %d to %d for %s\n",
                               old, howOften_, heuristicName_.c_str());
#endif
                    }
                    probability = 1.0 / howOften_;
                    if (model_->bestSolution())
                        probability *= 0.5;
                }
                break;
            case 7:
                if ((model_->bestSolution() && numRuns_ >= 2) || numRuns_ >= 4)
                    probability = -1.0;
                break;
            }
        }
        if (randomNumber > probability)
            return false;

        if (model_->getCurrentPassNumber() > 1)
            return false;
#ifdef COIN_DEVELOP
        printf("Running %s, random %g probability %g\n",
               heuristicName_.c_str(), randomNumber, probability);
#endif
    } else {
#ifdef COIN_DEVELOP
        printf("Running %s, depth %d when %d\n",
               heuristicName_.c_str(), depth, when_);
#endif
    }
    ++numRuns_;
    return true;
}

// Resets stuff if model changes
void
CbcHeuristic::resetModel(CbcModel * model)
{
    model_ = model;
}
// Set seed
void
CbcHeuristic::setSeed(int value)
{
    if (value==0) {
      double time = fabs(CoinGetTimeOfDay());
      while (time>=COIN_INT_MAX)
	time *= 0.5;
      value = static_cast<int>(time);
      char printArray[100];
      sprintf(printArray, "using time of day seed was changed from %d to %d",
	      randomNumberGenerator_.getSeed(), value);
      if (model_)
	model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
	  << printArray
	  << CoinMessageEol;
    }
    randomNumberGenerator_.setSeed(value);
}
// Get seed
int
CbcHeuristic::getSeed() const
{
  return randomNumberGenerator_.getSeed();
}

// Create C++ lines to get to current state
void
CbcHeuristic::generateCpp( FILE * fp, const char * heuristic)
{
    // hard coded as CbcHeuristic virtual
    if (when_ != 2)
        fprintf(fp, "3  %s.setWhen(%d);\n", heuristic, when_);
    else
        fprintf(fp, "4  %s.setWhen(%d);\n", heuristic, when_);
    if (numberNodes_ != 200)
        fprintf(fp, "3  %s.setNumberNodes(%d);\n", heuristic, numberNodes_);
    else
        fprintf(fp, "4  %s.setNumberNodes(%d);\n", heuristic, numberNodes_);
    if (feasibilityPumpOptions_ != -1)
        fprintf(fp, "3  %s.setFeasibilityPumpOptions(%d);\n", heuristic, feasibilityPumpOptions_);
    else
        fprintf(fp, "4  %s.setFeasibilityPumpOptions(%d);\n", heuristic, feasibilityPumpOptions_);
    if (fractionSmall_ != 1.0)
        fprintf(fp, "3  %s.setFractionSmall(%g);\n", heuristic, fractionSmall_);
    else
        fprintf(fp, "4  %s.setFractionSmall(%g);\n", heuristic, fractionSmall_);
    if (heuristicName_ != "Unknown")
        fprintf(fp, "3  %s.setHeuristicName(\"%s\");\n",
                heuristic, heuristicName_.c_str()) ;
    else
        fprintf(fp, "4  %s.setHeuristicName(\"%s\");\n",
                heuristic, heuristicName_.c_str()) ;
    if (decayFactor_ != 0.0)
        fprintf(fp, "3  %s.setDecayFactor(%g);\n", heuristic, decayFactor_);
    else
        fprintf(fp, "4  %s.setDecayFactor(%g);\n", heuristic, decayFactor_);
    if (switches_ != 0)
        fprintf(fp, "3  %s.setSwitches(%d);\n", heuristic, switches_);
    else
        fprintf(fp, "4  %s.setSwitches(%d);\n", heuristic, switches_);
    if (whereFrom_ != DEFAULT_WHERE)
        fprintf(fp, "3  %s.setWhereFrom(%d);\n", heuristic, whereFrom_);
    else
        fprintf(fp, "4  %s.setWhereFrom(%d);\n", heuristic, whereFrom_);
    if (shallowDepth_ != 1)
        fprintf(fp, "3  %s.setShallowDepth(%d);\n", heuristic, shallowDepth_);
    else
        fprintf(fp, "4  %s.setShallowDepth(%d);\n", heuristic, shallowDepth_);
    if (howOftenShallow_ != 1)
        fprintf(fp, "3  %s.setHowOftenShallow(%d);\n", heuristic, howOftenShallow_);
    else
        fprintf(fp, "4  %s.setHowOftenShallow(%d);\n", heuristic, howOftenShallow_);
    if (minDistanceToRun_ != 1)
        fprintf(fp, "3  %s.setMinDistanceToRun(%d);\n", heuristic, minDistanceToRun_);
    else
        fprintf(fp, "4  %s.setMinDistanceToRun(%d);\n", heuristic, minDistanceToRun_);
}
// Destructor
CbcHeuristic::~CbcHeuristic ()
{
    delete [] inputSolution_;
}

// update model
void CbcHeuristic::setModel(CbcModel * model)
{
    model_ = model;
}
/* Clone but ..
   type 0 clone solver, 1 clone continuous solver
   Add 2 to say without integer variables which are at low priority
   Add 4 to say quite likely infeasible so give up easily.*/
OsiSolverInterface *
CbcHeuristic::cloneBut(int type)
{
    OsiSolverInterface * solver;
    if ((type&1) == 0 || !model_->continuousSolver())
        solver = model_->solver()->clone();
    else
        solver = model_->continuousSolver()->clone();
#ifdef COIN_HAS_CLP
    OsiClpSolverInterface * clpSolver
    = dynamic_cast<OsiClpSolverInterface *> (solver);
#endif
    if ((type&2) != 0) {
        int n = model_->numberObjects();
        int priority = model_->continuousPriority();
        if (priority < COIN_INT_MAX) {
            for (int i = 0; i < n; i++) {
                const OsiObject * obj = model_->object(i);
                const CbcSimpleInteger * thisOne =
                    dynamic_cast <const CbcSimpleInteger *> (obj);
                if (thisOne) {
                    int iColumn = thisOne->columnNumber();
                    if (thisOne->priority() >= priority)
                        solver->setContinuous(iColumn);
                }
            }
        }
#ifdef COIN_HAS_CLP
        if (clpSolver) {
            for (int i = 0; i < n; i++) {
                const OsiObject * obj = model_->object(i);
                const CbcSimpleInteger * thisOne =
                    dynamic_cast <const CbcSimpleInteger *> (obj);
                if (thisOne) {
                    int iColumn = thisOne->columnNumber();
                    if (clpSolver->isOptionalInteger(iColumn))
                        clpSolver->setContinuous(iColumn);
                }
            }
        }
#endif
    }
#ifdef COIN_HAS_CLP
    if ((type&4) != 0 && clpSolver) {
        int options = clpSolver->getModelPtr()->moreSpecialOptions();
        clpSolver->getModelPtr()->setMoreSpecialOptions(options | 64);
    }
#endif
    return solver;
}
// Whether to exit at once on gap
bool
CbcHeuristic::exitNow(double bestObjective) const
{
    if ((switches_&2048) != 0) {
        // exit may be forced - but unset for next time
        switches_ &= ~2048;
        if ((switches_&1024) != 0)
            return true;
    } else if ((switches_&1) == 0) {
        return false;
    }
    // See if can stop on gap
    OsiSolverInterface * solver = model_->solver();
    double bestPossibleObjective = solver->getObjValue() * solver->getObjSense();
    double absGap = CoinMax(model_->getAllowableGap(),
                            model_->getHeuristicGap());
    double fracGap = CoinMax(model_->getAllowableFractionGap(),
                             model_->getHeuristicFractionGap());
    double testGap = CoinMax(absGap, fracGap *
                             CoinMax(fabs(bestObjective),
                                     fabs(bestPossibleObjective)));

    if (bestObjective - bestPossibleObjective < testGap
            && model_->getCutoffIncrement() >= 0.0) {
        return true;
    } else {
        return false;
    }
}
#ifdef HISTORY_STATISTICS
extern bool getHistoryStatistics_;
#endif
static double sizeRatio(int numberRowsNow, int numberColumnsNow,
                        int numberRowsStart, int numberColumnsStart)
{
    double valueNow;
    if (numberRowsNow*10 > numberColumnsNow || numberColumnsNow < 200) {
        valueNow = 2 * numberRowsNow + numberColumnsNow;
    } else {
        // long and thin - rows are more important
        if (numberRowsNow*40 > numberColumnsNow)
            valueNow = 10 * numberRowsNow + numberColumnsNow;
        else
            valueNow = 200 * numberRowsNow + numberColumnsNow;
    }
    double valueStart;
    if (numberRowsStart*10 > numberColumnsStart || numberColumnsStart < 200) {
        valueStart = 2 * numberRowsStart + numberColumnsStart;
    } else {
        // long and thin - rows are more important
        if (numberRowsStart*40 > numberColumnsStart)
            valueStart = 10 * numberRowsStart + numberColumnsStart;
        else
            valueStart = 200 * numberRowsStart + numberColumnsStart;
    }
    //printf("sizeProblem Now %g, %d rows, %d columns\nsizeProblem Start %g, %d rows, %d columns\n",
    // valueNow,numberRowsNow,numberColumnsNow,
    // valueStart,numberRowsStart,numberColumnsStart);
    if (10*numberRowsNow < 8*numberRowsStart || 10*numberColumnsNow < 7*numberColumnsStart)
        return valueNow / valueStart;
    else if (10*numberRowsNow < 9*numberRowsStart)
        return 1.1*(valueNow / valueStart);
    else if (numberRowsNow < numberRowsStart)
        return 1.5*(valueNow / valueStart);
    else
        return 2.0*(valueNow / valueStart);
}


// Do mini branch and bound (return 1 if solution)
int
CbcHeuristic::smallBranchAndBound(OsiSolverInterface * solver, int numberNodes,
                                  double * newSolution, double & newSolutionValue,
                                  double cutoff, std::string name) const
{
    // size before
    int shiftRows = 0;
    if (numberNodes < 0)
        shiftRows = solver->getNumRows() - numberNodes_;
    int numberRowsStart = solver->getNumRows() - shiftRows;
    int numberColumnsStart = solver->getNumCols();
#ifdef CLP_INVESTIGATE
    printf("%s has %d rows, %d columns\n",
           name.c_str(), solver->getNumRows(), solver->getNumCols());
#endif
    // Use this fraction
    double fractionSmall = fractionSmall_;
    double before = 2 * numberRowsStart + numberColumnsStart;
    if (before > 40000.0) {
        // fairly large - be more conservative
        double multiplier = 1.0 - 0.3 * CoinMin(100000.0, before - 40000.0) / 100000.0;
        if (multiplier < 1.0) {
            fractionSmall *= multiplier;
#ifdef CLP_INVESTIGATE
            printf("changing fractionSmall from %g to %g for %s\n",
                   fractionSmall_, fractionSmall, name.c_str());
#endif
        }
    }
#ifdef COIN_HAS_CLP
    OsiClpSolverInterface * osiclp = dynamic_cast< OsiClpSolverInterface*> (solver);
    if (osiclp && (osiclp->specialOptions()&65536) == 0) {
        // go faster stripes
        if (osiclp->getNumRows() < 300 && osiclp->getNumCols() < 500) {
            osiclp->setupForRepeatedUse(2, 0);
        } else {
            osiclp->setupForRepeatedUse(0, 0);
        }
        // Turn this off if you get problems
        // Used to be automatically set
        osiclp->setSpecialOptions(osiclp->specialOptions() | (128 + 64 - 128));
        ClpSimplex * lpSolver = osiclp->getModelPtr();
        lpSolver->setSpecialOptions(lpSolver->specialOptions() | 0x01000000); // say is Cbc (and in branch and bound)
        lpSolver->setSpecialOptions(lpSolver->specialOptions() |
                                    (/*16384+*/4096 + 512 + 128));
    }
#endif
#ifdef HISTORY_STATISTICS
    getHistoryStatistics_ = false;
#endif
#ifdef COIN_DEVELOP
    int status = 0;
#endif
    int logLevel = model_->logLevel();
#define LEN_PRINT 250
    char generalPrint[LEN_PRINT];
    // Do presolve to see if possible
    int numberColumns = solver->getNumCols();
    char * reset = NULL;
    int returnCode = 1;
    int saveModelOptions = model_->specialOptions();
    //assert ((saveModelOptions&2048) == 0);
    model_->setSpecialOptions(saveModelOptions | 2048);
    {
        int saveLogLevel = solver->messageHandler()->logLevel();
        if (saveLogLevel == 1) 
            solver->messageHandler()->setLogLevel(0);
        OsiPresolve * pinfo = new OsiPresolve();
        int presolveActions = 0;
        // Allow dual stuff on integers
        presolveActions = 1;
        // Do not allow all +1 to be tampered with
        //if (allPlusOnes)
        //presolveActions |= 2;
        // allow transfer of costs
        // presolveActions |= 4;
        pinfo->setPresolveActions(presolveActions);
        OsiSolverInterface * presolvedModel = pinfo->presolvedModel(*solver, 1.0e-8, true, 2);
        delete pinfo;
        // see if too big

        if (presolvedModel) {
            int afterRows = presolvedModel->getNumRows();
            int afterCols = presolvedModel->getNumCols();
            //#define COIN_DEVELOP
#ifdef COIN_DEVELOP_z
            if (numberNodes < 0) {
                solver->writeMpsNative("before.mps", NULL, NULL, 2, 1);
                presolvedModel->writeMpsNative("after1.mps", NULL, NULL, 2, 1);
            }
#endif
            delete presolvedModel;
            double ratio = sizeRatio(afterRows - shiftRows, afterCols,
                                     numberRowsStart, numberColumnsStart);
            double after = 2 * afterRows + afterCols;
            if (ratio > fractionSmall && after > 300 && numberNodes >= 0) {
                // Need code to try again to compress further using used
                const int * used =  model_->usedInSolution();
                int maxUsed = 0;
                int iColumn;
                const double * lower = solver->getColLower();
                const double * upper = solver->getColUpper();
                for (iColumn = 0; iColumn < numberColumns; iColumn++) {
                    if (upper[iColumn] > lower[iColumn]) {
                        if (solver->isBinary(iColumn))
                            maxUsed = CoinMax(maxUsed, used[iColumn]);
                    }
                }
                if (maxUsed) {
                    reset = new char [numberColumns];
                    int nFix = 0;
                    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
                        reset[iColumn] = 0;
                        if (upper[iColumn] > lower[iColumn]) {
                            if (solver->isBinary(iColumn) && used[iColumn] == maxUsed) {
                                bool setValue = true;
                                if (maxUsed == 1) {
                                    double randomNumber = randomNumberGenerator_.randomDouble();
                                    if (randomNumber > 0.3)
                                        setValue = false;
                                }
                                if (setValue) {
                                    reset[iColumn] = 1;
                                    solver->setColLower(iColumn, 1.0);
                                    nFix++;
                                }
                            }
                        }
                    }
                    pinfo = new OsiPresolve();
                    presolveActions = 0;
                    // Allow dual stuff on integers
                    presolveActions = 1;
                    // Do not allow all +1 to be tampered with
                    //if (allPlusOnes)
                    //presolveActions |= 2;
                    // allow transfer of costs
                    // presolveActions |= 4;
                    pinfo->setPresolveActions(presolveActions);
                    presolvedModel = pinfo->presolvedModel(*solver, 1.0e-8, true, 2);
                    delete pinfo;
                    if (presolvedModel) {
                        // see if too big
                        int afterRows2 = presolvedModel->getNumRows();
                        int afterCols2 = presolvedModel->getNumCols();
                        delete presolvedModel;
                        double ratio = sizeRatio(afterRows2 - shiftRows, afterCols2,
                                                 numberRowsStart, numberColumnsStart);
                        double after = 2 * afterRows2 + afterCols2;
                        if (ratio > fractionSmall && (after > 300 || numberNodes < 0)) {
                            sprintf(generalPrint, "Full problem %d rows %d columns, reduced to %d rows %d columns - %d fixed gives %d, %d - still too large",
                                    solver->getNumRows(), solver->getNumCols(),
                                    afterRows, afterCols, nFix, afterRows2, afterCols2);
                            // If much too big - give up
                            if (ratio > 0.75)
                                returnCode = -1;
                        } else {
                            sprintf(generalPrint, "Full problem %d rows %d columns, reduced to %d rows %d columns - %d fixed gives %d, %d - ok now",
                                    solver->getNumRows(), solver->getNumCols(),
                                    afterRows, afterCols, nFix, afterRows2, afterCols2);
                        }
                        model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
                        << generalPrint
                        << CoinMessageEol;
                    } else {
                        returnCode = 2; // infeasible
                    }
                }
            } else if (ratio > fractionSmall && after > 300 && numberNodes >=0) {
                returnCode = -1;
            }
        } else {
            returnCode = 2; // infeasible
        }
        solver->messageHandler()->setLogLevel(saveLogLevel);
    }
    if (returnCode == 2 || returnCode == -1) {
        model_->setSpecialOptions(saveModelOptions);
        delete [] reset;
#ifdef HISTORY_STATISTICS
        getHistoryStatistics_ = true;
#endif
        //printf("small no good\n");
        return returnCode;
    }
    // Reduce printout
    bool takeHint;
    OsiHintStrength strength;
    solver->getHintParam(OsiDoReducePrint, takeHint, strength);
    solver->setHintParam(OsiDoReducePrint, true, OsiHintTry);
    solver->setHintParam(OsiDoPresolveInInitial, false, OsiHintTry);
    double signedCutoff = cutoff*solver->getObjSense();
    solver->setDblParam(OsiDualObjectiveLimit, signedCutoff);
    solver->initialSolve();
    if (solver->isProvenOptimal()) {
        CglPreProcess process;
	OsiSolverInterface * solver2 = NULL;
	if ((model_->moreSpecialOptions()&65536)!=0)
	  process.setOptions(2+4+8); // no cuts
	/* Do not try and produce equality cliques and
	   do up to 2 passes (normally) 5 if restart */
	int numberPasses = 2;
	if (numberNodes < 0) {
	  numberPasses = 5;
	  // Say some rows cuts
	  int numberRows = solver->getNumRows();
	  if (numberNodes_ < numberRows && true /* think */) {
	    char * type = new char[numberRows];
	    memset(type, 0, numberNodes_);
	    memset(type + numberNodes_, 1, numberRows - numberNodes_);
	    process.passInRowTypes(type, numberRows);
	    delete [] type;
	  }
	}
	if (logLevel <= 1)
	  process.messageHandler()->setLogLevel(0);
	if (!solver->defaultHandler()&&
	    solver->messageHandler()->logLevel(0)!=-1000)
	  process.passInMessageHandler(solver->messageHandler());
	solver2 = process.preProcessNonDefault(*solver, false,
					       numberPasses);
	  if (!solver2) {
            if (logLevel > 1)
	      printf("Pre-processing says infeasible\n");
            returnCode = 2; // so will be infeasible
	  } else {
#ifdef COIN_DEVELOP_z
            if (numberNodes < 0) {
	      solver2->writeMpsNative("after2.mps", NULL, NULL, 2, 1);
            }
#endif
            // see if too big
            double ratio = sizeRatio(solver2->getNumRows() - shiftRows, solver2->getNumCols(),
                                     numberRowsStart, numberColumnsStart);
            double after = 2 * solver2->getNumRows() + solver2->getNumCols();
            if (ratio > fractionSmall && (after > 300 || numberNodes < 0)) {
                sprintf(generalPrint, "Full problem %d rows %d columns, reduced to %d rows %d columns - too large",
                        solver->getNumRows(), solver->getNumCols(),
                        solver2->getNumRows(), solver2->getNumCols());
                model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
                << generalPrint
                << CoinMessageEol;
                returnCode = -1;
                //printf("small no good2\n");
            } else {
                sprintf(generalPrint, "Full problem %d rows %d columns, reduced to %d rows %d columns",
                        solver->getNumRows(), solver->getNumCols(),
                        solver2->getNumRows(), solver2->getNumCols());
                model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
                << generalPrint
                << CoinMessageEol;
            }
            if (returnCode == 1) {
                solver2->resolve();
                CbcModel model(*solver2);
		// move seed across
		model.randomNumberGenerator()->setSeed(model_->randomNumberGenerator()->getSeed());
                if (numberNodes >= 0) {
                    // normal
                    model.setSpecialOptions(saveModelOptions | 2048);
                    if (logLevel <= 1 && feasibilityPumpOptions_ != -3)
                        model.setLogLevel(0);
                    else
                        model.setLogLevel(logLevel);
                    // No small fathoming
                    model.setFastNodeDepth(-1);
                    model.setCutoff(signedCutoff);
		    model.setStrongStrategy(0);
		    // Don't do if original fraction > 1.0 and too large
		    if (fractionSmall_>1.0 && fractionSmall_ < 1000000.0) {
		      /* 1.4 means -1 nodes if >.4
			 2.4 means -1 nodes if >.5 and 0 otherwise
			 3.4 means -1 nodes if >.6 and 0 or 5
			 4.4 means -1 nodes if >.7 and 0, 5 or 10
		      */
		      double fraction = fractionSmall_-floor(fractionSmall_);
		      if (ratio>fraction) {
			int type = static_cast<int>(floor(fractionSmall_*0.1));
			int over = static_cast<int>(ceil(ratio-fraction));
			int maxNodes[]={-1,0,5,10};
			if (type>over)
			  numberNodes=maxNodes[type-over];
			else
			  numberNodes=-1;
		      }
		    }
                    model.setMaximumNodes(numberNodes);
                    model.solver()->setHintParam(OsiDoReducePrint, true, OsiHintTry);
		    if ((saveModelOptions&2048) == 0)
		      model.setMoreSpecialOptions(model_->moreSpecialOptions());
		    // off conflict analysis
		    model.setMoreSpecialOptions(model.moreSpecialOptions()&~4194304);
		    
                    // Lightweight
                    CbcStrategyDefaultSubTree strategy(model_, 1, 5, 1, 0);
                    model.setStrategy(strategy);
                    model.solver()->setIntParam(OsiMaxNumIterationHotStart, 10);
                    model.setMaximumCutPassesAtRoot(CoinMin(20, CoinAbs(model_->getMaximumCutPassesAtRoot())));
                    model.setMaximumCutPasses(CoinMin(10, model_->getMaximumCutPasses()));
                } else {
                    model.setSpecialOptions(saveModelOptions);
                    model_->messageHandler()->message(CBC_RESTART, model_->messages())
                    << solver2->getNumRows() << solver2->getNumCols()
                    << CoinMessageEol;
                    // going for full search and copy across more stuff
                    model.gutsOfCopy(*model_, 2);
		    assert (!model_->heuristicModel());
		    model_->setHeuristicModel(&model);
                    for (int i = 0; i < model.numberCutGenerators(); i++) {
		        CbcCutGenerator * generator = model.cutGenerator(i);
			CglGomory * gomory = dynamic_cast<CglGomory *>
			  (generator->generator());
			if (gomory&&gomory->originalSolver()) 
			  gomory->passInOriginalSolver(model.solver());
                        generator->setTiming(true);
                        // Turn on if was turned on
                        int iOften = model_->cutGenerator(i)->howOften();
#ifdef CLP_INVESTIGATE
                        printf("Gen %d often %d %d\n",
                               i, generator->howOften(),
                               iOften);
#endif
                        if (iOften > 0)
                            generator->setHowOften(iOften % 1000000);
                        if (model_->cutGenerator(i)->howOftenInSub() == -200)
                            generator->setHowOften(-100);
                    }
                    model.setCutoff(signedCutoff);
                    // make sure can't do nested search! but allow heuristics
                    model.setSpecialOptions((model.specialOptions()&(~(512 + 2048))) | 1024);
                    bool takeHint;
                    OsiHintStrength strength;
                    // Switch off printing if asked to
                    model_->solver()->getHintParam(OsiDoReducePrint, takeHint, strength);
                    model.solver()->setHintParam(OsiDoReducePrint, takeHint, strength);
		    // no cut generators if none in parent
                    CbcStrategyDefault 
		      strategy(model_->numberCutGenerators() ? 1 : -1, 
			       model_->numberStrong(),
			       model_->numberBeforeTrust());
                    // Set up pre-processing - no
                    strategy.setupPreProcessing(0); // was (4);
                    model.setStrategy(strategy);
                    //model.solver()->writeMps("crunched");
                    int numberCuts = process.cuts().sizeRowCuts();
                    if (numberCuts) {
                        // add in cuts
                        CglStored cuts = process.cuts();
                        model.addCutGenerator(&cuts, 1, "Stored from first");
			model.cutGenerator(model.numberCutGenerators()-1)->setGlobalCuts(true);
                    }
                }
                // Do search
                if (logLevel > 1)
                    model_->messageHandler()->message(CBC_START_SUB, model_->messages())
                    << name
                    << model.getMaximumNodes()
                    << CoinMessageEol;
                // probably faster to use a basis to get integer solutions
                model.setSpecialOptions(model.specialOptions() | 2);
#ifdef CBC_THREAD
                if (model_->getNumberThreads() > 0 && (model_->getThreadMode()&4) != 0) {
                    // See if at root node
                    bool atRoot = model_->getNodeCount() == 0;
                    int passNumber = model_->getCurrentPassNumber();
                    if (atRoot && passNumber == 1)
                        model.setNumberThreads(model_->getNumberThreads());
                }
#endif
                model.setParentModel(*model_);
		model.setMaximumSolutions(model_->getMaximumSolutions()); 
		model.setOriginalColumns(process.originalColumns());
                model.setSearchStrategy(-1);
                // If no feasibility pump then insert a lightweight one
                if (feasibilityPumpOptions_ >= 0 || feasibilityPumpOptions_ == -2) {
		    CbcHeuristicFPump * fpump = NULL;
                    for (int i = 0; i < model.numberHeuristics(); i++) {
                        CbcHeuristicFPump* pump =
                            dynamic_cast<CbcHeuristicFPump*>(model.heuristic(i));
                        if (pump) {
                            fpump = pump;
                            break;
			}
                    }
                    if (!fpump) {
                        CbcHeuristicFPump heuristic4;
			// use any cutoff
			heuristic4.setFakeCutoff(0.5*COIN_DBL_MAX);
			if (fractionSmall_<=1.0) 
			  heuristic4.setMaximumPasses(10);
                        int pumpTune = feasibilityPumpOptions_;
			if (pumpTune==-2)
			  pumpTune = 4; // proximity
                        if (pumpTune > 0) {
                            /*
                            >=10000000 for using obj
                            >=1000000 use as accumulate switch
                            >=1000 use index+1 as number of large loops
                            >=100 use 0.05 objvalue as increment
                            %100 == 10,20 etc for experimentation
                            1 == fix ints at bounds, 2 fix all integral ints, 3 and continuous at bounds
                            4 and static continuous, 5 as 3 but no internal integers
                            6 as 3 but all slack basis!
                            */
                            double value = solver2->getObjSense() * solver2->getObjValue();
                            int w = pumpTune / 10;
                            int ix = w % 10;
                            w /= 10;
                            int c = w % 10;
                            w /= 10;
                            int r = w;
                            int accumulate = r / 1000;
                            r -= 1000 * accumulate;
                            if (accumulate >= 10) {
                                int which = accumulate / 10;
                                accumulate -= 10 * which;
                                which--;
                                // weights and factors
                                double weight[] = {0.1, 0.1, 0.5, 0.5, 1.0, 1.0, 5.0, 5.0};
                                double factor[] = {0.1, 0.5, 0.1, 0.5, 0.1, 0.5, 0.1, 0.5};
                                heuristic4.setInitialWeight(weight[which]);
                                heuristic4.setWeightFactor(factor[which]);
                            }
                            // fake cutoff
                            if (c) {
                                double cutoff;
                                solver2->getDblParam(OsiDualObjectiveLimit, cutoff);
                                cutoff = CoinMin(cutoff, value + 0.1 * fabs(value) * c);
                                heuristic4.setFakeCutoff(cutoff);
                            }
                            if (r) {
                                // also set increment
                                //double increment = (0.01*i+0.005)*(fabs(value)+1.0e-12);
                                double increment = 0.0;
                                heuristic4.setAbsoluteIncrement(increment);
                                heuristic4.setAccumulate(accumulate);
                                heuristic4.setMaximumRetries(r + 1);
                            }
                            pumpTune = pumpTune % 100;
                            if (pumpTune == 6)
                                pumpTune = 13;
                            if (pumpTune != 13)
                                pumpTune = pumpTune % 10;
                            heuristic4.setWhen(pumpTune);
                            if (ix) {
                                heuristic4.setFeasibilityPumpOptions(ix*10);
                            }
                        }
                        model.addHeuristic(&heuristic4, "feasibility pump", 0);
	
                    }
		} else if (feasibilityPumpOptions_==-3) {
		  // add all (except this)
		  for (int i = 0; i < model_->numberHeuristics(); i++) {
		    if (strcmp(heuristicName(),model_->heuristic(i)->heuristicName()))
		      model.addHeuristic(model_->heuristic(i)); 
		  }
		}
                //printf("sol %x\n",inputSolution_);
                if (inputSolution_) {
                    // translate and add a serendipity heuristic
                    int numberColumns = solver2->getNumCols();
                    const int * which = process.originalColumns();
                    OsiSolverInterface * solver3 = solver2->clone();
                    for (int i = 0; i < numberColumns; i++) {
                        if (solver3->isInteger(i)) {
                            int k = which[i];
                            double value = inputSolution_[k];
                            //if (value)
                            //printf("orig col %d now %d val %g\n",
                            //       k,i,value);
                            solver3->setColLower(i, value);
                            solver3->setColUpper(i, value);
                        }
                    }
                    solver3->setDblParam(OsiDualObjectiveLimit, COIN_DBL_MAX);
                    solver3->resolve();
                    if (!solver3->isProvenOptimal()) {
                        // Try just setting nonzeros
                        OsiSolverInterface * solver4 = solver2->clone();
                        for (int i = 0; i < numberColumns; i++) {
                            if (solver4->isInteger(i)) {
                                int k = which[i];
                                double value = floor(inputSolution_[k] + 0.5);
                                if (value) {
                                    solver3->setColLower(i, value);
                                    solver3->setColUpper(i, value);
                                }
                            }
                        }
                        solver4->setDblParam(OsiDualObjectiveLimit, COIN_DBL_MAX);
                        solver4->resolve();
                        int nBad = -1;
                        if (solver4->isProvenOptimal()) {
                            nBad = 0;
                            const double * solution = solver4->getColSolution();
                            for (int i = 0; i < numberColumns; i++) {
                                if (solver4->isInteger(i)) {
                                    double value = floor(solution[i] + 0.5);
                                    if (fabs(value - solution[i]) > 1.0e-6)
                                        nBad++;
                                }
                            }
                        }
                        if (nBad) {
                            delete solver4;
                        } else {
                            delete solver3;
                            solver3 = solver4;
                        }
                    }
                    if (solver3->isProvenOptimal()) {
                        // good
                        CbcSerendipity heuristic(model);
                        double value = solver3->getObjSense() * solver3->getObjValue();
                        heuristic.setInputSolution(solver3->getColSolution(), value);
                        value = value + 1.0e-7*(1.0 + fabs(value));
			value *= solver3->getObjSense();
                        model.setCutoff(value);
                        model.addHeuristic(&heuristic, "Previous solution", 0);
                        //printf("added seren\n");
                    } else {
                        double value = model_->getMinimizationObjValue();
                        value = value + 1.0e-7*(1.0 + fabs(value));
			value *= solver3->getObjSense();
                        model.setCutoff(value);
#ifdef CLP_INVESTIGATE
                        printf("NOT added seren\n");
                        solver3->writeMps("bad_seren");
                        solver->writeMps("orig_seren");
#endif
                    }
                    delete solver3;
                }
                if (model_->searchStrategy() == 2) {
                    model.setNumberStrong(5);
                    model.setNumberBeforeTrust(5);
                }
                if (model.getNumCols()) {
                    if (numberNodes >= 0) {
                        setCutAndHeuristicOptions(model);
                        // not too many iterations
                        model.setMaximumNumberIterations(100*(numberNodes + 10));
                        // Not fast stuff
                        model.setFastNodeDepth(-1);
                    } else if (model.fastNodeDepth() >= 1000000) {
                        // already set
                        model.setFastNodeDepth(model.fastNodeDepth() - 1000000);
                    }
                    model.setWhenCuts(999998);
#define ALWAYS_DUAL
#ifdef ALWAYS_DUAL
		    OsiSolverInterface * solverD = model.solver();
		    bool takeHint;
		    OsiHintStrength strength;
		    solverD->getHintParam(OsiDoDualInResolve, takeHint, strength);
		    solverD->setHintParam(OsiDoDualInResolve, true, OsiHintDo);
#endif
		    model.passInEventHandler(model_->getEventHandler());
		    // say model_ is sitting there
		    int saveOptions = model_->specialOptions();
		    model_->setSpecialOptions(saveOptions|1048576);
                    model.branchAndBound();
		    model_->setHeuristicModel(NULL);
		    model_->setSpecialOptions(saveOptions);
#ifdef ALWAYS_DUAL
		    solverD = model.solver();
		    solverD->setHintParam(OsiDoDualInResolve, takeHint, strength);
#endif
		    numberNodesDone_ = model.getNodeCount();
#ifdef COIN_DEVELOP
                    printf("sub branch %d nodes, %d iterations - max %d\n",
                           model.getNodeCount(), model.getIterationCount(),
                           100*(numberNodes + 10));
#endif
                    if (numberNodes < 0) {
                        model_->incrementIterationCount(model.getIterationCount());
                        model_->incrementNodeCount(model.getNodeCount());
			// update best solution (in case ctrl-c)
			// !!! not a good idea - think a bit harder
			//model_->setMinimizationObjValue(model.getMinimizationObjValue());
                        for (int iGenerator = 0; iGenerator < model.numberCutGenerators(); iGenerator++) {
                            CbcCutGenerator * generator = model.cutGenerator(iGenerator);
                            sprintf(generalPrint,
                                    "%s was tried %d times and created %d cuts of which %d were active after adding rounds of cuts (%.3f seconds)",
                                    generator->cutGeneratorName(),
                                    generator->numberTimesEntered(),
                                    generator->numberCutsInTotal() +
                                    generator->numberColumnCuts(),
                                    generator->numberCutsActive(),
                                    generator->timeInCutGenerator());
                            CglStored * stored = dynamic_cast<CglStored*>(generator->generator());
                            if (stored && !generator->numberCutsInTotal())
                                continue;
#ifndef CLP_INVESTIGATE
                            CglImplication * implication = dynamic_cast<CglImplication*>(generator->generator());
                            if (implication)
                                continue;
#endif
                            model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
                            << generalPrint
                            << CoinMessageEol;
                        }
                    }
                } else {
                    // empty model
                    model.setMinimizationObjValue(model.solver()->getObjSense()*model.solver()->getObjValue());
                }
                if (logLevel > 1)
                    model_->messageHandler()->message(CBC_END_SUB, model_->messages())
                    << name
                    << CoinMessageEol;
                if (model.getMinimizationObjValue() < CoinMin(cutoff, 1.0e30)) {
                    // solution
                    if (model.getNumCols())
                        returnCode = model.isProvenOptimal() ? 3 : 1;
                    else
                        returnCode = 3;
                    // post process
#ifdef COIN_HAS_CLP
                    OsiClpSolverInterface * clpSolver = dynamic_cast< OsiClpSolverInterface*> (model.solver());
                    if (clpSolver) {
                        ClpSimplex * lpSolver = clpSolver->getModelPtr();
                        lpSolver->setSpecialOptions(lpSolver->specialOptions() | 0x01000000); // say is Cbc (and in branch and bound)
                    }
#endif
		    //if (fractionSmall_ < 1000000.0) 
		      process.postProcess(*model.solver());
                    if (solver->isProvenOptimal() && solver->getObjValue()*solver->getObjSense() < cutoff) {
                        // Solution now back in solver
                        int numberColumns = solver->getNumCols();
                        memcpy(newSolution, solver->getColSolution(),
                               numberColumns*sizeof(double));
                        newSolutionValue = model.getMinimizationObjValue();
                    } else {
                        // odd - but no good
                        returnCode = 0; // so will be infeasible
                    }
                } else {
                    // no good
                    returnCode = model.isProvenInfeasible() ? 2 : 0; // so will be infeasible
                }
                int totalNumberIterations = model.getIterationCount() +
                                            process.numberIterationsPre() +
                                            process.numberIterationsPost();
                if (totalNumberIterations > 100*(numberNodes + 10)
		    && fractionSmall_ < 1000000.0) {
                    // only allow smaller problems
                    fractionSmall = fractionSmall_;
                    fractionSmall_ *= 0.9;
#ifdef CLP_INVESTIGATE
                    printf("changing fractionSmall from %g to %g for %s as %d iterations\n",
                           fractionSmall, fractionSmall_, name.c_str(), totalNumberIterations);
#endif
                }
                if (model.status() == 5)
                   model_->sayEventHappened();
#ifdef COIN_DEVELOP
                if (model.isProvenInfeasible())
                    status = 1;
                else if (model.isProvenOptimal())
                    status = 2;
#endif
            }
        }
    } else {
        returnCode = 2; // infeasible finished
    }
    model_->setSpecialOptions(saveModelOptions);
    model_->setLogLevel(logLevel);
    if (returnCode == 1 || returnCode == 2) {
        OsiSolverInterface * solverC = model_->continuousSolver();
        if (false && solverC) {
            const double * lower = solver->getColLower();
            const double * upper = solver->getColUpper();
            const double * lowerC = solverC->getColLower();
            const double * upperC = solverC->getColUpper();
            bool good = true;
            for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
                if (solverC->isInteger(iColumn)) {
                    if (lower[iColumn] > lowerC[iColumn] &&
                            upper[iColumn] < upperC[iColumn]) {
                        good = false;
                        printf("CUT - can't add\n");
                        break;
                    }
                }
            }
            if (good) {
                double * cut = new double [numberColumns];
                int * which = new int [numberColumns];
                double rhs = -1.0;
                int n = 0;
                for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
                    if (solverC->isInteger(iColumn)) {
                        if (lower[iColumn] == upperC[iColumn]) {
                            rhs += lower[iColumn];
                            cut[n] = 1.0;
                            which[n++] = iColumn;
                        } else if (upper[iColumn] == lowerC[iColumn]) {
                            rhs -= upper[iColumn];
                            cut[n] = -1.0;
                            which[n++] = iColumn;
                        }
                    }
                }
                printf("CUT has %d entries\n", n);
                OsiRowCut newCut;
                newCut.setLb(-COIN_DBL_MAX);
                newCut.setUb(rhs);
                newCut.setRow(n, which, cut, false);
                model_->makeGlobalCut(newCut);
                delete [] cut;
                delete [] which;
            }
        }
#ifdef COIN_DEVELOP
        if (status == 1)
            printf("heuristic could add cut because infeasible (%s)\n", heuristicName_.c_str());
        else if (status == 2)
            printf("heuristic could add cut because optimal (%s)\n", heuristicName_.c_str());
#endif
    }
    if (reset) {
        for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
            if (reset[iColumn])
                solver->setColLower(iColumn, 0.0);
        }
        delete [] reset;
    }
#ifdef HISTORY_STATISTICS
    getHistoryStatistics_ = true;
#endif
    solver->setHintParam(OsiDoReducePrint, takeHint, strength);
    return returnCode;
}
// Set input solution
void
CbcHeuristic::setInputSolution(const double * solution, double objValue)
{
    delete [] inputSolution_;
    inputSolution_ = NULL;
    if (model_ && solution) {
        int numberColumns = model_->getNumCols();
        inputSolution_ = new double [numberColumns+1];
        memcpy(inputSolution_, solution, numberColumns*sizeof(double));
        inputSolution_[numberColumns] = objValue;
    }
}

//##############################################################################

inline int compare3BranchingObjects(const CbcBranchingObject* br0,
                                    const CbcBranchingObject* br1)
{
    const int t0 = br0->type();
    const int t1 = br1->type();
    if (t0 < t1) {
        return -1;
    }
    if (t0 > t1) {
        return 1;
    }
    return br0->compareOriginalObject(br1);
}

//==============================================================================

inline bool compareBranchingObjects(const CbcBranchingObject* br0,
                                    const CbcBranchingObject* br1)
{
    return compare3BranchingObjects(br0, br1) < 0;
}

//==============================================================================

void
CbcHeuristicNode::gutsOfConstructor(CbcModel& model)
{
    //  CbcHeurDebugNodes(&model);
    CbcNode* node = model.currentNode();
    brObj_ = new CbcBranchingObject*[node->depth()];
    CbcNodeInfo* nodeInfo = node->nodeInfo();
    int cnt = 0;
    while (nodeInfo->parentBranch() != NULL) {
        const OsiBranchingObject* br = nodeInfo->parentBranch();
        const CbcBranchingObject* cbcbr = dynamic_cast<const CbcBranchingObject*>(br);
        if (! cbcbr) {
            throw CoinError("CbcHeuristicNode can be used only with CbcBranchingObjects.\n",
                            "gutsOfConstructor",
                            "CbcHeuristicNode",
                            __FILE__, __LINE__);
        }
        brObj_[cnt] = cbcbr->clone();
        brObj_[cnt]->previousBranch();
        ++cnt;
        nodeInfo = nodeInfo->parent();
    }
    std::sort(brObj_, brObj_ + cnt, compareBranchingObjects);
    if (cnt <= 1) {
        numObjects_ = cnt;
    } else {
        numObjects_ = 0;
        CbcBranchingObject* br = NULL; // What should this be?
        for (int i = 1; i < cnt; ++i) {
            if (compare3BranchingObjects(brObj_[numObjects_], brObj_[i]) == 0) {
                int comp = brObj_[numObjects_]->compareBranchingObject(brObj_[i], br != 0);
                switch (comp) {
                case CbcRangeSame: // the same range
                case CbcRangeDisjoint: // disjoint decisions
                    // should not happen! we are on a chain!
                    abort();
                case CbcRangeSubset: // brObj_[numObjects_] is a subset of brObj_[i]
                    delete brObj_[i];
                    break;
                case CbcRangeSuperset: // brObj_[i] is a subset of brObj_[numObjects_]
                    delete brObj_[numObjects_];
                    brObj_[numObjects_] = brObj_[i];
                    break;
                case CbcRangeOverlap: // overlap
                    delete brObj_[i];
                    delete brObj_[numObjects_];
                    brObj_[numObjects_] = br;
                    break;
                }
                continue;
            } else {
                brObj_[++numObjects_] = brObj_[i];
            }
        }
        ++numObjects_;
    }
}

//==============================================================================

CbcHeuristicNode::CbcHeuristicNode(CbcModel& model)
{
    gutsOfConstructor(model);
}

//==============================================================================

double
CbcHeuristicNode::distance(const CbcHeuristicNode* node) const
{

    const double disjointWeight = 1;
    const double overlapWeight = 0.4;
    const double subsetWeight = 0.2;
    int countDisjointWeight = 0;
    int countOverlapWeight = 0;
    int countSubsetWeight = 0;
    int i = 0;
    int j = 0;
    double dist = 0.0;
#ifdef PRINT_DEBUG
    printf(" numObjects_ = %i, node->numObjects_ = %i\n",
           numObjects_, node->numObjects_);
#endif
    while ( i < numObjects_ && j < node->numObjects_) {
        CbcBranchingObject* br0 = brObj_[i];
        const CbcBranchingObject* br1 = node->brObj_[j];
#ifdef PRINT_DEBUG
        const CbcIntegerBranchingObject* brPrint0 =
            dynamic_cast<const CbcIntegerBranchingObject*>(br0);
        const double* downBounds = brPrint0->downBounds();
        const double* upBounds = brPrint0->upBounds();
        int variable = brPrint0->variable();
        int way = brPrint0->way();
        printf("   br0: var %i downBd [%i,%i] upBd [%i,%i] way %i\n",
               variable, static_cast<int>(downBounds[0]), static_cast<int>(downBounds[1]),
               static_cast<int>(upBounds[0]), static_cast<int>(upBounds[1]), way);
        const CbcIntegerBranchingObject* brPrint1 =
            dynamic_cast<const CbcIntegerBranchingObject*>(br1);
        downBounds = brPrint1->downBounds();
        upBounds = brPrint1->upBounds();
        variable = brPrint1->variable();
        way = brPrint1->way();
        printf("   br1: var %i downBd [%i,%i] upBd [%i,%i] way %i\n",
               variable, static_cast<int>(downBounds[0]), static_cast<int>(downBounds[1]),
               static_cast<int>(upBounds[0]), static_cast<int>(upBounds[1]), way);
#endif
        const int brComp = compare3BranchingObjects(br0, br1);
        if (brComp < 0) {
            dist += subsetWeight;
            countSubsetWeight++;
            ++i;
        } else if (brComp > 0) {
            dist += subsetWeight;
            countSubsetWeight++;
            ++j;
        } else {
            const int comp = br0->compareBranchingObject(br1, false);
            switch (comp) {
            case CbcRangeSame:
                // do nothing
                break;
            case CbcRangeDisjoint: // disjoint decisions
                dist += disjointWeight;
                countDisjointWeight++;
                break;
            case CbcRangeSubset: // subset one way or another
            case CbcRangeSuperset:
                dist += subsetWeight;
                countSubsetWeight++;
                break;
            case CbcRangeOverlap: // overlap
                dist += overlapWeight;
                countOverlapWeight++;
                break;
            }
            ++i;
            ++j;
        }
    }
    dist += subsetWeight * (numObjects_ - i + node->numObjects_ - j);
    countSubsetWeight += (numObjects_ - i + node->numObjects_ - j);
    COIN_DETAIL_PRINT(printf("subset = %i, overlap = %i, disjoint = %i\n", countSubsetWeight,
			     countOverlapWeight, countDisjointWeight));
    return dist;
}

//==============================================================================

CbcHeuristicNode::~CbcHeuristicNode()
{
    for (int i = 0; i < numObjects_; ++i) {
        delete brObj_[i];
    }
    delete [] brObj_;
}

//==============================================================================

double
CbcHeuristicNode::minDistance(const CbcHeuristicNodeList& nodeList) const
{
    double minDist = COIN_DBL_MAX;
    for (int i = nodeList.size() - 1; i >= 0; --i) {
        minDist = CoinMin(minDist, distance(nodeList.node(i)));
    }
    return minDist;
}

//==============================================================================

bool
CbcHeuristicNode::minDistanceIsSmall(const CbcHeuristicNodeList& nodeList,
                                     const double threshold) const
{
    for (int i = nodeList.size() - 1; i >= 0; --i) {
        if (distance(nodeList.node(i)) >= threshold) {
            continue;
        } else {
            return true;
        }
    }
    return false;
}

//==============================================================================

double
CbcHeuristicNode::avgDistance(const CbcHeuristicNodeList& nodeList) const
{
    if (nodeList.size() == 0) {
        return COIN_DBL_MAX;
    }
    double sumDist = 0;
    for (int i = nodeList.size() - 1; i >= 0; --i) {
        sumDist += distance(nodeList.node(i));
    }
    return sumDist / nodeList.size();
}

//##############################################################################

// Default Constructor
CbcRounding::CbcRounding()
        : CbcHeuristic()
{
    // matrix and row copy will automatically be empty
    seed_ = 7654321;
    down_ = NULL;
    up_ = NULL;
    equal_ = NULL;
    //whereFrom_ |= 16; // allow more often
}

// Constructor from model
CbcRounding::CbcRounding(CbcModel & model)
        : CbcHeuristic(model)
{
    // Get a copy of original matrix (and by row for rounding);
    assert(model.solver());
    if (model.solver()->getNumRows()) {
        matrix_ = *model.solver()->getMatrixByCol();
        matrixByRow_ = *model.solver()->getMatrixByRow();
        validate();
    }
    down_ = NULL;
    up_ = NULL;
    equal_ = NULL;
    seed_ = 7654321;
    //whereFrom_ |= 16; // allow more often
}

// Destructor
CbcRounding::~CbcRounding ()
{
    delete [] down_;
    delete [] up_;
    delete [] equal_;
}

// Clone
CbcHeuristic *
CbcRounding::clone() const
{
    return new CbcRounding(*this);
}
// Create C++ lines to get to current state
void
CbcRounding::generateCpp( FILE * fp)
{
    CbcRounding other;
    fprintf(fp, "0#include \"CbcHeuristic.hpp\"\n");
    fprintf(fp, "3  CbcRounding rounding(*cbcModel);\n");
    CbcHeuristic::generateCpp(fp, "rounding");
    if (seed_ != other.seed_)
        fprintf(fp, "3  rounding.setSeed(%d);\n", seed_);
    else
        fprintf(fp, "4  rounding.setSeed(%d);\n", seed_);
    fprintf(fp, "3  cbcModel->addHeuristic(&rounding);\n");
}
//#define NEW_ROUNDING
// Copy constructor
CbcRounding::CbcRounding(const CbcRounding & rhs)
        :
        CbcHeuristic(rhs),
        matrix_(rhs.matrix_),
        matrixByRow_(rhs.matrixByRow_),
        seed_(rhs.seed_)
{
#ifdef NEW_ROUNDING
    int numberColumns = matrix_.getNumCols();
    down_ = CoinCopyOfArray(rhs.down_, numberColumns);
    up_ = CoinCopyOfArray(rhs.up_, numberColumns);
    equal_ = CoinCopyOfArray(rhs.equal_, numberColumns);
#else
    down_ = NULL;
    up_ = NULL;
    equal_ = NULL;
#endif
}

// Assignment operator
CbcRounding &
CbcRounding::operator=( const CbcRounding & rhs)
{
    if (this != &rhs) {
        CbcHeuristic::operator=(rhs);
        matrix_ = rhs.matrix_;
        matrixByRow_ = rhs.matrixByRow_;
#ifdef NEW_ROUNDING
        delete [] down_;
        delete [] up_;
        delete [] equal_;
        int numberColumns = matrix_.getNumCols();
        down_ = CoinCopyOfArray(rhs.down_, numberColumns);
        up_ = CoinCopyOfArray(rhs.up_, numberColumns);
        equal_ = CoinCopyOfArray(rhs.equal_, numberColumns);
#else
        down_ = NULL;
        up_ = NULL;
        equal_ = NULL;
#endif
        seed_ = rhs.seed_;
    }
    return *this;
}

// Resets stuff if model changes
void
CbcRounding::resetModel(CbcModel * model)
{
    model_ = model;
    // Get a copy of original matrix (and by row for rounding);
    assert(model_->solver());
    matrix_ = *model_->solver()->getMatrixByCol();
    matrixByRow_ = *model_->solver()->getMatrixByRow();
    validate();
}
// See if rounding will give solution
// Sets value of solution
// Assumes rhs for original matrix still okay
// At present only works with integers
// Fix values if asked for
// Returns 1 if solution, 0 if not
int
CbcRounding::solution(double & solutionValue,
                      double * betterSolution)
{

    numCouldRun_++;
    // See if to do
    if (!when() || (when() % 10 == 1 && model_->phase() != 1) ||
            (when() % 10 == 2 && (model_->phase() != 2 && model_->phase() != 3)))
        return 0; // switched off
    numRuns_++;
    OsiSolverInterface * solver = model_->solver();
    double direction = solver->getObjSense();
    double newSolutionValue = direction * solver->getObjValue();
    return solution(solutionValue, betterSolution, newSolutionValue);
}
// See if rounding will give solution
// Sets value of solution
// Assumes rhs for original matrix still okay
// At present only works with integers
// Fix values if asked for
// Returns 1 if solution, 0 if not
int
CbcRounding::solution(double & solutionValue,
                      double * betterSolution,
                      double newSolutionValue)
{

    // See if to do
    if (!when() || (when() % 10 == 1 && model_->phase() != 1) ||
            (when() % 10 == 2 && (model_->phase() != 2 && model_->phase() != 3)))
        return 0; // switched off
    OsiSolverInterface * solver = model_->solver();
    const double * lower = solver->getColLower();
    const double * upper = solver->getColUpper();
    const double * rowLower = solver->getRowLower();
    const double * rowUpper = solver->getRowUpper();
    const double * solution = solver->getColSolution();
    const double * objective = solver->getObjCoefficients();
    double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
    double primalTolerance;
    solver->getDblParam(OsiPrimalTolerance, primalTolerance);

    int numberRows = matrix_.getNumRows();
    assert (numberRows <= solver->getNumRows());
    int numberIntegers = model_->numberIntegers();
    const int * integerVariable = model_->integerVariable();
    int i;
    double direction = solver->getObjSense();
    //double newSolutionValue = direction*solver->getObjValue();
    int returnCode = 0;
    // Column copy
    const double * element = matrix_.getElements();
    const int * row = matrix_.getIndices();
    const CoinBigIndex * columnStart = matrix_.getVectorStarts();
    const int * columnLength = matrix_.getVectorLengths();
    // Row copy
    const double * elementByRow = matrixByRow_.getElements();
    const int * column = matrixByRow_.getIndices();
    const CoinBigIndex * rowStart = matrixByRow_.getVectorStarts();
    const int * rowLength = matrixByRow_.getVectorLengths();

    // Get solution array for heuristic solution
    int numberColumns = solver->getNumCols();
    double * newSolution = new double [numberColumns];
    memcpy(newSolution, solution, numberColumns*sizeof(double));

    double * rowActivity = new double[numberRows];
    memset(rowActivity, 0, numberRows*sizeof(double));
    for (i = 0; i < numberColumns; i++) {
        int j;
        double value = newSolution[i];
        if (value < lower[i]) {
            value = lower[i];
            newSolution[i] = value;
        } else if (value > upper[i]) {
            value = upper[i];
            newSolution[i] = value;
        }
        if (value) {
            for (j = columnStart[i];
                    j < columnStart[i] + columnLength[i]; j++) {
                int iRow = row[j];
                rowActivity[iRow] += value * element[j];
            }
        }
    }
    // check was feasible - if not adjust (cleaning may move)
    for (i = 0; i < numberRows; i++) {
        if (rowActivity[i] < rowLower[i]) {
            //assert (rowActivity[i]>rowLower[i]-1000.0*primalTolerance);
            rowActivity[i] = rowLower[i];
        } else if (rowActivity[i] > rowUpper[i]) {
            //assert (rowActivity[i]<rowUpper[i]+1000.0*primalTolerance);
            rowActivity[i] = rowUpper[i];
        }
    }
    for (i = 0; i < numberIntegers; i++) {
        int iColumn = integerVariable[i];
        double value = newSolution[iColumn];
        if (fabs(floor(value + 0.5) - value) > integerTolerance) {
            double below = floor(value);
            double newValue = newSolution[iColumn];
            double cost = direction * objective[iColumn];
            double move;
            if (cost > 0.0) {
                // try up
                move = 1.0 - (value - below);
            } else if (cost < 0.0) {
                // try down
                move = below - value;
            } else {
                // won't be able to move unless we can grab another variable
                double randomNumber = randomNumberGenerator_.randomDouble();
                // which way?
                if (randomNumber < 0.5)
                    move = below - value;
                else
                    move = 1.0 - (value - below);
            }
            newValue += move;
            newSolution[iColumn] = newValue;
            newSolutionValue += move * cost;
            int j;
            for (j = columnStart[iColumn];
                    j < columnStart[iColumn] + columnLength[iColumn]; j++) {
                int iRow = row[j];
                rowActivity[iRow] += move * element[j];
            }
        }
    }

    double penalty = 0.0;
    const char * integerType = model_->integerType();
    // see if feasible - just using singletons
    for (i = 0; i < numberRows; i++) {
        double value = rowActivity[i];
        double thisInfeasibility = 0.0;
        if (value < rowLower[i] - primalTolerance)
            thisInfeasibility = value - rowLower[i];
        else if (value > rowUpper[i] + primalTolerance)
            thisInfeasibility = value - rowUpper[i];
        if (thisInfeasibility) {
            // See if there are any slacks I can use to fix up
            // maybe put in coding for multiple slacks?
            double bestCost = 1.0e50;
            int k;
            int iBest = -1;
            double addCost = 0.0;
            double newValue = 0.0;
            double changeRowActivity = 0.0;
            double absInfeasibility = fabs(thisInfeasibility);
            for (k = rowStart[i]; k < rowStart[i] + rowLength[i]; k++) {
                int iColumn = column[k];
                // See if all elements help
                if (columnLength[iColumn] == 1) {
                    double currentValue = newSolution[iColumn];
                    double elementValue = elementByRow[k];
                    double lowerValue = lower[iColumn];
                    double upperValue = upper[iColumn];
                    double gap = rowUpper[i] - rowLower[i];
                    double absElement = fabs(elementValue);
                    if (thisInfeasibility*elementValue > 0.0) {
                        // we want to reduce
                        if ((currentValue - lowerValue)*absElement >= absInfeasibility) {
                            // possible - check if integer
                            double distance = absInfeasibility / absElement;
                            double thisCost = -direction * objective[iColumn] * distance;
                            if (integerType[iColumn]) {
                                distance = ceil(distance - primalTolerance);
                                if (currentValue - distance >= lowerValue - primalTolerance) {
                                    if (absInfeasibility - distance*absElement < -gap - primalTolerance)
                                        thisCost = 1.0e100; // no good
                                    else
                                        thisCost = -direction * objective[iColumn] * distance;
                                } else {
                                    thisCost = 1.0e100; // no good
                                }
                            }
                            if (thisCost < bestCost) {
                                bestCost = thisCost;
                                iBest = iColumn;
                                addCost = thisCost;
                                newValue = currentValue - distance;
                                changeRowActivity = -distance * elementValue;
                            }
                        }
                    } else {
                        // we want to increase
                        if ((upperValue - currentValue)*absElement >= absInfeasibility) {
                            // possible - check if integer
                            double distance = absInfeasibility / absElement;
                            double thisCost = direction * objective[iColumn] * distance;
                            if (integerType[iColumn]) {
                                distance = ceil(distance - 1.0e-7);
                                assert (currentValue - distance <= upperValue + primalTolerance);
                                if (absInfeasibility - distance*absElement < -gap - primalTolerance)
                                    thisCost = 1.0e100; // no good
                                else
                                    thisCost = direction * objective[iColumn] * distance;
                            }
                            if (thisCost < bestCost) {
                                bestCost = thisCost;
                                iBest = iColumn;
                                addCost = thisCost;
                                newValue = currentValue + distance;
                                changeRowActivity = distance * elementValue;
                            }
                        }
                    }
                }
            }
            if (iBest >= 0) {
                /*printf("Infeasibility of %g on row %d cost %g\n",
                  thisInfeasibility,i,addCost);*/
                newSolution[iBest] = newValue;
                thisInfeasibility = 0.0;
                newSolutionValue += addCost;
                rowActivity[i] += changeRowActivity;
            }
            penalty += fabs(thisInfeasibility);
        }
    }
    if (penalty) {
        // see if feasible using any
        // first continuous
        double penaltyChange = 0.0;
        int iColumn;
        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
            if (integerType[iColumn])
                continue;
            double currentValue = newSolution[iColumn];
            double lowerValue = lower[iColumn];
            double upperValue = upper[iColumn];
            int j;
            int anyBadDown = 0;
            int anyBadUp = 0;
            double upImprovement = 0.0;
            double downImprovement = 0.0;
            for (j = columnStart[iColumn];
                    j < columnStart[iColumn] + columnLength[iColumn]; j++) {
                int iRow = row[j];
                if (rowUpper[iRow] > rowLower[iRow]) {
                    double value = element[j];
                    if (rowActivity[iRow] > rowUpper[iRow] + primalTolerance) {
                        // infeasible above
                        downImprovement += value;
                        upImprovement -= value;
                        if (value > 0.0)
                            anyBadUp++;
                        else
                            anyBadDown++;
                    } else if (rowActivity[iRow] > rowUpper[iRow] - primalTolerance) {
                        // feasible at ub
                        if (value > 0.0) {
                            upImprovement -= value;
                            anyBadUp++;
                        } else {
                            downImprovement += value;
                            anyBadDown++;
                        }
                    } else if (rowActivity[iRow] > rowLower[iRow] + primalTolerance) {
                        // feasible in interior
                    } else if (rowActivity[iRow] > rowLower[iRow] - primalTolerance) {
                        // feasible at lb
                        if (value < 0.0) {
                            upImprovement += value;
                            anyBadUp++;
                        } else {
                            downImprovement -= value;
                            anyBadDown++;
                        }
                    } else {
                        // infeasible below
                        downImprovement -= value;
                        upImprovement += value;
                        if (value < 0.0)
                            anyBadUp++;
                        else
                            anyBadDown++;
                    }
                } else {
                    // equality row
                    double value = element[j];
                    if (rowActivity[iRow] > rowUpper[iRow] + primalTolerance) {
                        // infeasible above
                        downImprovement += value;
                        upImprovement -= value;
                        if (value > 0.0)
                            anyBadUp++;
                        else
                            anyBadDown++;
                    } else if (rowActivity[iRow] < rowLower[iRow] - primalTolerance) {
                        // infeasible below
                        downImprovement -= value;
                        upImprovement += value;
                        if (value < 0.0)
                            anyBadUp++;
                        else
                            anyBadDown++;
                    } else {
                        // feasible - no good
                        anyBadUp = -1;
                        anyBadDown = -1;
                        break;
                    }
                }
            }
            // could change tests for anyBad
            if (anyBadUp)
                upImprovement = 0.0;
            if (anyBadDown)
                downImprovement = 0.0;
            double way = 0.0;
            double improvement = 0.0;
            if (downImprovement > 0.0 && currentValue > lowerValue) {
                way = -1.0;
                improvement = downImprovement;
            } else if (upImprovement > 0.0 && currentValue < upperValue) {
                way = 1.0;
                improvement = upImprovement;
            }
            if (way) {
                // can improve
                double distance;
                if (way > 0.0)
                    distance = upperValue - currentValue;
                else
                    distance = currentValue - lowerValue;
                for (j = columnStart[iColumn];
                        j < columnStart[iColumn] + columnLength[iColumn]; j++) {
                    int iRow = row[j];
                    double value = element[j] * way;
                    if (rowActivity[iRow] > rowUpper[iRow] + primalTolerance) {
                        // infeasible above
                        assert (value < 0.0);
                        double gap = rowActivity[iRow] - rowUpper[iRow];
                        if (gap + value*distance < 0.0)
                            distance = -gap / value;
                    } else if (rowActivity[iRow] < rowLower[iRow] - primalTolerance) {
                        // infeasible below
                        assert (value > 0.0);
                        double gap = rowActivity[iRow] - rowLower[iRow];
                        if (gap + value*distance > 0.0)
                            distance = -gap / value;
                    } else {
                        // feasible
                        if (value > 0) {
                            double gap = rowActivity[iRow] - rowUpper[iRow];
                            if (gap + value*distance > 0.0)
                                distance = -gap / value;
                        } else {
                            double gap = rowActivity[iRow] - rowLower[iRow];
                            if (gap + value*distance < 0.0)
                                distance = -gap / value;
                        }
                    }
                }
                //move
                penaltyChange += improvement * distance;
                distance *= way;
                newSolution[iColumn] += distance;
                newSolutionValue += direction * objective[iColumn] * distance;
                for (j = columnStart[iColumn];
                        j < columnStart[iColumn] + columnLength[iColumn]; j++) {
                    int iRow = row[j];
                    double value = element[j];
                    rowActivity[iRow] += distance * value;
                }
            }
        }
        // and now all if improving
        double lastChange = penaltyChange ? 1.0 : 0.0;
        while (lastChange > 1.0e-2) {
            lastChange = 0;
            for (iColumn = 0; iColumn < numberColumns; iColumn++) {
                bool isInteger = (integerType[iColumn] != 0);
                double currentValue = newSolution[iColumn];
                double lowerValue = lower[iColumn];
                double upperValue = upper[iColumn];
                int j;
                int anyBadDown = 0;
                int anyBadUp = 0;
                double upImprovement = 0.0;
                double downImprovement = 0.0;
                for (j = columnStart[iColumn];
                        j < columnStart[iColumn] + columnLength[iColumn]; j++) {
                    int iRow = row[j];
                    double value = element[j];
                    if (isInteger) {
                        if (value > 0.0) {
                            if (rowActivity[iRow] + value > rowUpper[iRow] + primalTolerance)
                                anyBadUp++;
                            if (rowActivity[iRow] - value < rowLower[iRow] - primalTolerance)
                                anyBadDown++;
                        } else {
                            if (rowActivity[iRow] - value > rowUpper[iRow] + primalTolerance)
                                anyBadDown++;
                            if (rowActivity[iRow] + value < rowLower[iRow] - primalTolerance)
                                anyBadUp++;
                        }
                    }
                    if (rowUpper[iRow] > rowLower[iRow]) {
                        if (rowActivity[iRow] > rowUpper[iRow] + primalTolerance) {
                            // infeasible above
                            downImprovement += value;
                            upImprovement -= value;
                            if (value > 0.0)
                                anyBadUp++;
                            else
                                anyBadDown++;
                        } else if (rowActivity[iRow] > rowUpper[iRow] - primalTolerance) {
                            // feasible at ub
                            if (value > 0.0) {
                                upImprovement -= value;
                                anyBadUp++;
                            } else {
                                downImprovement += value;
                                anyBadDown++;
                            }
                        } else if (rowActivity[iRow] > rowLower[iRow] + primalTolerance) {
                            // feasible in interior
                        } else if (rowActivity[iRow] > rowLower[iRow] - primalTolerance) {
                            // feasible at lb
                            if (value < 0.0) {
                                upImprovement += value;
                                anyBadUp++;
                            } else {
                                downImprovement -= value;
                                anyBadDown++;
                            }
                        } else {
                            // infeasible below
                            downImprovement -= value;
                            upImprovement += value;
                            if (value < 0.0)
                                anyBadUp++;
                            else
                                anyBadDown++;
                        }
                    } else {
                        // equality row
                        if (rowActivity[iRow] > rowUpper[iRow] + primalTolerance) {
                            // infeasible above
                            downImprovement += value;
                            upImprovement -= value;
                            if (value > 0.0)
                                anyBadUp++;
                            else
                                anyBadDown++;
                        } else if (rowActivity[iRow] < rowLower[iRow] - primalTolerance) {
                            // infeasible below
                            downImprovement -= value;
                            upImprovement += value;
                            if (value < 0.0)
                                anyBadUp++;
                            else
                                anyBadDown++;
                        } else {
                            // feasible - no good
                            anyBadUp = -1;
                            anyBadDown = -1;
                            break;
                        }
                    }
                }
                // could change tests for anyBad
                if (anyBadUp)
                    upImprovement = 0.0;
                if (anyBadDown)
                    downImprovement = 0.0;
                double way = 0.0;
                double improvement = 0.0;
                if (downImprovement > 0.0 && currentValue > lowerValue) {
                    way = -1.0;
                    improvement = downImprovement;
                } else if (upImprovement > 0.0 && currentValue < upperValue) {
                    way = 1.0;
                    improvement = upImprovement;
                }
                if (way) {
                    // can improve
                    double distance = COIN_DBL_MAX;
                    for (j = columnStart[iColumn];
                            j < columnStart[iColumn] + columnLength[iColumn]; j++) {
                        int iRow = row[j];
                        double value = element[j] * way;
                        if (rowActivity[iRow] > rowUpper[iRow] + primalTolerance) {
                            // infeasible above
                            assert (value < 0.0);
                            double gap = rowActivity[iRow] - rowUpper[iRow];
                            if (gap + value*distance < 0.0) {
                                // If integer then has to move by 1
                                if (!isInteger)
                                    distance = -gap / value;
                                else
                                    distance = CoinMax(-gap / value, 1.0);
                            }
                        } else if (rowActivity[iRow] < rowLower[iRow] - primalTolerance) {
                            // infeasible below
                            assert (value > 0.0);
                            double gap = rowActivity[iRow] - rowLower[iRow];
                            if (gap + value*distance > 0.0) {
                                // If integer then has to move by 1
                                if (!isInteger)
                                    distance = -gap / value;
                                else
                                    distance = CoinMax(-gap / value, 1.0);
                            }
                        } else {
                            // feasible
                            if (value > 0) {
                                double gap = rowActivity[iRow] - rowUpper[iRow];
                                if (gap + value*distance > 0.0)
                                    distance = -gap / value;
                            } else {
                                double gap = rowActivity[iRow] - rowLower[iRow];
                                if (gap + value*distance < 0.0)
                                    distance = -gap / value;
                            }
                        }
                    }
                    if (isInteger)
                        distance = floor(distance + 1.05e-8);
                    if (!distance) {
                        // should never happen
                        //printf("zero distance in CbcRounding - debug\n");
                    }
                    //move
                    lastChange += improvement * distance;
                    distance *= way;
                    newSolution[iColumn] += distance;
                    newSolutionValue += direction * objective[iColumn] * distance;
                    for (j = columnStart[iColumn];
                            j < columnStart[iColumn] + columnLength[iColumn]; j++) {
                        int iRow = row[j];
                        double value = element[j];
                        rowActivity[iRow] += distance * value;
                    }
                }
            }
            penaltyChange += lastChange;
        }
        penalty -= penaltyChange;
        if (penalty < 1.0e-5*fabs(penaltyChange)) {
            // recompute
            penalty = 0.0;
            for (i = 0; i < numberRows; i++) {
                double value = rowActivity[i];
                if (value < rowLower[i] - primalTolerance)
                    penalty += rowLower[i] - value;
                else if (value > rowUpper[i] + primalTolerance)
                    penalty += value - rowUpper[i];
            }
        }
    }

    // Could also set SOS (using random) and repeat
    if (!penalty) {
        // See if we can do better
        //seed_++;
        //CoinSeedRandom(seed_);
        // Random number between 0 and 1.
        double randomNumber = randomNumberGenerator_.randomDouble();
        int iPass;
        int start[2];
        int end[2];
        int iRandom = static_cast<int> (randomNumber * (static_cast<double> (numberIntegers)));
        start[0] = iRandom;
        end[0] = numberIntegers;
        start[1] = 0;
        end[1] = iRandom;
        for (iPass = 0; iPass < 2; iPass++) {
            int i;
            for (i = start[iPass]; i < end[iPass]; i++) {
                int iColumn = integerVariable[i];
#ifndef NDEBUG
                double value = newSolution[iColumn];
                assert (fabs(floor(value + 0.5) - value) < integerTolerance);
#endif
                double cost = direction * objective[iColumn];
                double move = 0.0;
                if (cost > 0.0)
                    move = -1.0;
                else if (cost < 0.0)
                    move = 1.0;
                while (move) {
                    bool good = true;
                    double newValue = newSolution[iColumn] + move;
                    if (newValue < lower[iColumn] - primalTolerance ||
                            newValue > upper[iColumn] + primalTolerance) {
                        move = 0.0;
                    } else {
                        // see if we can move
                        int j;
                        for (j = columnStart[iColumn];
                                j < columnStart[iColumn] + columnLength[iColumn]; j++) {
                            int iRow = row[j];
                            double newActivity = rowActivity[iRow] + move * element[j];
                            if (newActivity < rowLower[iRow] - primalTolerance ||
                                    newActivity > rowUpper[iRow] + primalTolerance) {
                                good = false;
                                break;
                            }
                        }
                        if (good) {
                            newSolution[iColumn] = newValue;
                            newSolutionValue += move * cost;
                            int j;
                            for (j = columnStart[iColumn];
                                    j < columnStart[iColumn] + columnLength[iColumn]; j++) {
                                int iRow = row[j];
                                rowActivity[iRow] += move * element[j];
                            }
                        } else {
                            move = 0.0;
                        }
                    }
                }
            }
        }
        // Just in case of some stupidity
        double objOffset = 0.0;
        solver->getDblParam(OsiObjOffset, objOffset);
        newSolutionValue = -objOffset;
        for ( i = 0 ; i < numberColumns ; i++ )
            newSolutionValue += objective[i] * newSolution[i];
        newSolutionValue *= direction;
        //printf("new solution value %g %g\n",newSolutionValue,solutionValue);
        if (newSolutionValue < solutionValue) {
            // paranoid check
            memset(rowActivity, 0, numberRows*sizeof(double));
            for (i = 0; i < numberColumns; i++) {
                int j;
                double value = newSolution[i];
                if (value) {
                    for (j = columnStart[i];
                            j < columnStart[i] + columnLength[i]; j++) {
                        int iRow = row[j];
                        rowActivity[iRow] += value * element[j];
                    }
                }
            }
            // check was approximately feasible
            bool feasible = true;
            for (i = 0; i < numberRows; i++) {
                if (rowActivity[i] < rowLower[i]) {
                    if (rowActivity[i] < rowLower[i] - 1000.0*primalTolerance)
                        feasible = false;
                } else if (rowActivity[i] > rowUpper[i]) {
                    if (rowActivity[i] > rowUpper[i] + 1000.0*primalTolerance)
                        feasible = false;
                }
            }
            if (feasible) {
                // new solution
                memcpy(betterSolution, newSolution, numberColumns*sizeof(double));
                solutionValue = newSolutionValue;
                //printf("** Solution of %g found by rounding\n",newSolutionValue);
                returnCode = 1;
            } else {
                // Can easily happen
                //printf("Debug CbcRounding giving bad solution\n");
            }
        }
    }
#ifdef NEW_ROUNDING
    if (!returnCode) {
#ifdef JJF_ZERO
        // back to starting point
        memcpy(newSolution, solution, numberColumns*sizeof(double));
        memset(rowActivity, 0, numberRows*sizeof(double));
        for (i = 0; i < numberColumns; i++) {
            int j;
            double value = newSolution[i];
            if (value < lower[i]) {
                value = lower[i];
                newSolution[i] = value;
            } else if (value > upper[i]) {
                value = upper[i];
                newSolution[i] = value;
            }
            if (value) {
                for (j = columnStart[i];
                        j < columnStart[i] + columnLength[i]; j++) {
                    int iRow = row[j];
                    rowActivity[iRow] += value * element[j];
                }
            }
        }
        // check was feasible - if not adjust (cleaning may move)
        for (i = 0; i < numberRows; i++) {
            if (rowActivity[i] < rowLower[i]) {
                //assert (rowActivity[i]>rowLower[i]-1000.0*primalTolerance);
                rowActivity[i] = rowLower[i];
            } else if (rowActivity[i] > rowUpper[i]) {
                //assert (rowActivity[i]<rowUpper[i]+1000.0*primalTolerance);
                rowActivity[i] = rowUpper[i];
            }
        }
#endif
        int * candidate = new int [numberColumns];
        int nCandidate = 0;
        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
            bool isInteger = (integerType[iColumn] != 0);
            if (isInteger) {
                double currentValue = newSolution[iColumn];
                if (fabs(currentValue - floor(currentValue + 0.5)) > 1.0e-8)
                    candidate[nCandidate++] = iColumn;
            }
        }
        if (true) {
            // Rounding as in Berthold
            while (nCandidate) {
                double infeasibility = 1.0e-7;
                int iRow = -1;
                for (i = 0; i < numberRows; i++) {
                    double value = 0.0;
                    if (rowActivity[i] < rowLower[i]) {
                        value = rowLower[i] - rowActivity[i];
                    } else if (rowActivity[i] > rowUpper[i]) {
                        value = rowActivity[i] - rowUpper[i];
                    }
                    if (value > infeasibility) {
                        infeasibility = value;
                        iRow = i;
                    }
                }
                if (iRow >= 0) {
                    // infeasible
                } else {
                    // feasible
                }
            }
        } else {
            // Shifting as in Berthold
        }
        delete [] candidate;
    }
#endif
    delete [] newSolution;
    delete [] rowActivity;
    return returnCode;
}
// update model
void CbcRounding::setModel(CbcModel * model)
{
    model_ = model;
    // Get a copy of original matrix (and by row for rounding);
    assert(model_->solver());
    if (model_->solver()->getNumRows()) {
        matrix_ = *model_->solver()->getMatrixByCol();
        matrixByRow_ = *model_->solver()->getMatrixByRow();
        // make sure model okay for heuristic
        validate();
    }
}
// Validate model i.e. sets when_ to 0 if necessary (may be NULL)
void
CbcRounding::validate()
{
    if (model_ && (when() % 100) < 10) {
        if (model_->numberIntegers() !=
                model_->numberObjects() && (model_->numberObjects() ||
                                            (model_->specialOptions()&1024) == 0)) {
            int numberOdd = 0;
            for (int i = 0; i < model_->numberObjects(); i++) {
                if (!model_->object(i)->canDoHeuristics())
                    numberOdd++;
            }
            if (numberOdd)
                setWhen(0);
        }
    }
#ifdef NEW_ROUNDING
    int numberColumns = matrix_.getNumCols();
    down_ = new unsigned short [numberColumns];
    up_ = new unsigned short [numberColumns];
    equal_ = new unsigned short [numberColumns];
    // Column copy
    const double * element = matrix_.getElements();
    const int * row = matrix_.getIndices();
    const CoinBigIndex * columnStart = matrix_.getVectorStarts();
    const int * columnLength = matrix_.getVectorLengths();
    const double * rowLower = model.solver()->getRowLower();
    const double * rowUpper = model.solver()->getRowUpper();
    for (int i = 0; i < numberColumns; i++) {
        int down = 0;
        int up = 0;
        int equal = 0;
        if (columnLength[i] > 65535) {
            equal[0] = 65535;
            break; // unlikely to work
        }
        for (CoinBigIndex j = columnStart[i];
                j < columnStart[i] + columnLength[i]; j++) {
            int iRow = row[j];
            if (rowLower[iRow] > -1.0e20 && rowUpper[iRow] < 1.0e20) {
                equal++;
            } else if (element[j] > 0.0) {
                if (rowUpper[iRow] < 1.0e20)
                    up++;
                else
                    down--;
            } else {
                if (rowLower[iRow] > -1.0e20)
                    up++;
                else
                    down--;
            }
        }
        down_[i] = (unsigned short) down;
        up_[i] = (unsigned short) up;
        equal_[i] = (unsigned short) equal;
    }
#else
    down_ = NULL;
    up_ = NULL;
    equal_ = NULL;
#endif
}

// Default Constructor
CbcHeuristicPartial::CbcHeuristicPartial()
        : CbcHeuristic()
{
    fixPriority_ = 10000;
}

// Constructor from model
CbcHeuristicPartial::CbcHeuristicPartial(CbcModel & model, int fixPriority, int numberNodes)
        : CbcHeuristic(model)
{
    fixPriority_ = fixPriority;
    setNumberNodes(numberNodes);
    validate();
}

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

// Clone
CbcHeuristic *
CbcHeuristicPartial::clone() const
{
    return new CbcHeuristicPartial(*this);
}
// Create C++ lines to get to current state
void
CbcHeuristicPartial::generateCpp( FILE * fp)
{
    CbcHeuristicPartial other;
    fprintf(fp, "0#include \"CbcHeuristic.hpp\"\n");
    fprintf(fp, "3  CbcHeuristicPartial partial(*cbcModel);\n");
    CbcHeuristic::generateCpp(fp, "partial");
    if (fixPriority_ != other.fixPriority_)
        fprintf(fp, "3  partial.setFixPriority(%d);\n", fixPriority_);
    else
        fprintf(fp, "4  partial.setFixPriority(%d);\n", fixPriority_);
    fprintf(fp, "3  cbcModel->addHeuristic(&partial);\n");
}
//#define NEW_PARTIAL
// Copy constructor
CbcHeuristicPartial::CbcHeuristicPartial(const CbcHeuristicPartial & rhs)
        :
        CbcHeuristic(rhs),
        fixPriority_(rhs.fixPriority_)
{
}

// Assignment operator
CbcHeuristicPartial &
CbcHeuristicPartial::operator=( const CbcHeuristicPartial & rhs)
{
    if (this != &rhs) {
        CbcHeuristic::operator=(rhs);
        fixPriority_ = rhs.fixPriority_;
    }
    return *this;
}

// Resets stuff if model changes
void
CbcHeuristicPartial::resetModel(CbcModel * model)
{
    model_ = model;
    // Get a copy of original matrix (and by row for partial);
    assert(model_->solver());
    validate();
}
// See if partial will give solution
// Sets value of solution
// Assumes rhs for original matrix still okay
// At present only works with integers
// Fix values if asked for
// Returns 1 if solution, 0 if not
int
CbcHeuristicPartial::solution(double & solutionValue,
                              double * betterSolution)
{
    // Return if already done
    if (fixPriority_ < 0)
        return 0; // switched off
    const double * hotstartSolution = model_->hotstartSolution();
    const int * hotstartPriorities = model_->hotstartPriorities();
    if (!hotstartSolution)
        return 0;
    OsiSolverInterface * solver = model_->solver();

    int numberIntegers = model_->numberIntegers();
    const int * integerVariable = model_->integerVariable();

    OsiSolverInterface * newSolver = model_->continuousSolver()->clone();
    const double * colLower = newSolver->getColLower();
    const double * colUpper = newSolver->getColUpper();

    double primalTolerance;
    solver->getDblParam(OsiPrimalTolerance, primalTolerance);

    int i;
    int numberFixed = 0;
    int returnCode = 0;

    for (i = 0; i < numberIntegers; i++) {
        int iColumn = integerVariable[i];
        if (abs(hotstartPriorities[iColumn]) <= fixPriority_) {
            double value = hotstartSolution[iColumn];
            double lower = colLower[iColumn];
            double upper = colUpper[iColumn];
            value = CoinMax(value, lower);
            value = CoinMin(value, upper);
            if (fabs(value - floor(value + 0.5)) < 1.0e-8) {
                value = floor(value + 0.5);
                newSolver->setColLower(iColumn, value);
                newSolver->setColUpper(iColumn, value);
                numberFixed++;
            }
        }
    }
    if (numberFixed > numberIntegers / 5 - 100000000) {
#ifdef COIN_DEVELOP
        printf("%d integers fixed\n", numberFixed);
#endif
        returnCode = smallBranchAndBound(newSolver, numberNodes_, betterSolution, solutionValue,
                                         model_->getCutoff(), "CbcHeuristicPartial");
        if (returnCode < 0)
            returnCode = 0; // returned on size
        //printf("return code %d",returnCode);
        if ((returnCode&2) != 0) {
            // could add cut
            returnCode &= ~2;
            //printf("could add cut with %d elements (if all 0-1)\n",nFix);
        } else {
            //printf("\n");
        }
    }
    fixPriority_ = -1; // switch off

    delete newSolver;
    return returnCode;
}
// update model
void CbcHeuristicPartial::setModel(CbcModel * model)
{
    model_ = model;
    assert(model_->solver());
    // make sure model okay for heuristic
    validate();
}
// Validate model i.e. sets when_ to 0 if necessary (may be NULL)
void
CbcHeuristicPartial::validate()
{
    if (model_ && (when() % 100) < 10) {
        if (model_->numberIntegers() !=
                model_->numberObjects())
            setWhen(0);
    }
}
bool
CbcHeuristicPartial::shouldHeurRun(int /*whereFrom*/)
{
    return true;
}

// Default Constructor
CbcSerendipity::CbcSerendipity()
        : CbcHeuristic()
{
}

// Constructor from model
CbcSerendipity::CbcSerendipity(CbcModel & model)
        : CbcHeuristic(model)
{
}

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

// Clone
CbcHeuristic *
CbcSerendipity::clone() const
{
    return new CbcSerendipity(*this);
}
// Create C++ lines to get to current state
void
CbcSerendipity::generateCpp( FILE * fp)
{
    fprintf(fp, "0#include \"CbcHeuristic.hpp\"\n");
    fprintf(fp, "3  CbcSerendipity serendipity(*cbcModel);\n");
    CbcHeuristic::generateCpp(fp, "serendipity");
    fprintf(fp, "3  cbcModel->addHeuristic(&serendipity);\n");
}

// Copy constructor
CbcSerendipity::CbcSerendipity(const CbcSerendipity & rhs)
        :
        CbcHeuristic(rhs)
{
}

// Assignment operator
CbcSerendipity &
CbcSerendipity::operator=( const CbcSerendipity & rhs)
{
    if (this != &rhs) {
        CbcHeuristic::operator=(rhs);
    }
    return *this;
}

// Returns 1 if solution, 0 if not
int
CbcSerendipity::solution(double & solutionValue,
                         double * betterSolution)
{
    if (!model_)
        return 0;
    if (!inputSolution_) {
        // get information on solver type
        OsiAuxInfo * auxInfo = model_->solver()->getAuxiliaryInfo();
        OsiBabSolver * auxiliaryInfo = dynamic_cast< OsiBabSolver *> (auxInfo);
        if (auxiliaryInfo) {
            return auxiliaryInfo->solution(solutionValue, betterSolution, model_->solver()->getNumCols());
        } else {
            return 0;
        }
    } else {
        int numberColumns = model_->getNumCols();
        double value = inputSolution_[numberColumns];
        int returnCode = 0;
        if (value < solutionValue) {
            solutionValue = value;
            memcpy(betterSolution, inputSolution_, numberColumns*sizeof(double));
            returnCode = 1;
        }
        delete [] inputSolution_;
        inputSolution_ = NULL;
        model_ = NULL; // switch off
        return returnCode;
    }
}
// update model
void CbcSerendipity::setModel(CbcModel * model)
{
    model_ = model;
}
// Resets stuff if model changes
void
CbcSerendipity::resetModel(CbcModel * model)
{
    model_ = model;
}


// Default Constructor
CbcHeuristicJustOne::CbcHeuristicJustOne()
        : CbcHeuristic(),
        probabilities_(NULL),
        heuristic_(NULL),
        numberHeuristics_(0)
{
}

// Constructor from model
CbcHeuristicJustOne::CbcHeuristicJustOne(CbcModel & model)
        : CbcHeuristic(model),
        probabilities_(NULL),
        heuristic_(NULL),
        numberHeuristics_(0)
{
}

// Destructor
CbcHeuristicJustOne::~CbcHeuristicJustOne ()
{
    for (int i = 0; i < numberHeuristics_; i++)
        delete heuristic_[i];
    delete [] heuristic_;
    delete [] probabilities_;
}

// Clone
CbcHeuristicJustOne *
CbcHeuristicJustOne::clone() const
{
    return new CbcHeuristicJustOne(*this);
}

// Create C++ lines to get to current state
void
CbcHeuristicJustOne::generateCpp( FILE * fp)
{
    CbcHeuristicJustOne other;
    fprintf(fp, "0#include \"CbcHeuristicJustOne.hpp\"\n");
    fprintf(fp, "3  CbcHeuristicJustOne heuristicJustOne(*cbcModel);\n");
    CbcHeuristic::generateCpp(fp, "heuristicJustOne");
    fprintf(fp, "3  cbcModel->addHeuristic(&heuristicJustOne);\n");
}

// Copy constructor
CbcHeuristicJustOne::CbcHeuristicJustOne(const CbcHeuristicJustOne & rhs)
        :
        CbcHeuristic(rhs),
        probabilities_(NULL),
        heuristic_(NULL),
        numberHeuristics_(rhs.numberHeuristics_)
{
    if (numberHeuristics_) {
        probabilities_ = CoinCopyOfArray(rhs.probabilities_, numberHeuristics_);
        heuristic_ = new CbcHeuristic * [numberHeuristics_];
        for (int i = 0; i < numberHeuristics_; i++)
            heuristic_[i] = rhs.heuristic_[i]->clone();
    }
}

// Assignment operator
CbcHeuristicJustOne &
CbcHeuristicJustOne::operator=( const CbcHeuristicJustOne & rhs)
{
    if (this != &rhs) {
        CbcHeuristic::operator=(rhs);
        for (int i = 0; i < numberHeuristics_; i++)
            delete heuristic_[i];
        delete [] heuristic_;
        delete [] probabilities_;
        probabilities_ = NULL;
        heuristic_ = NULL;
        numberHeuristics_ = rhs.numberHeuristics_;
        if (numberHeuristics_) {
            probabilities_ = CoinCopyOfArray(rhs.probabilities_, numberHeuristics_);
            heuristic_ = new CbcHeuristic * [numberHeuristics_];
            for (int i = 0; i < numberHeuristics_; i++)
                heuristic_[i] = rhs.heuristic_[i]->clone();
        }
    }
    return *this;
}
// Sets value of solution
// Returns 1 if solution, 0 if not
int
CbcHeuristicJustOne::solution(double & solutionValue,
                              double * betterSolution)
{
#ifdef DIVE_DEBUG
    std::cout << "solutionValue = " << solutionValue << std::endl;
#endif
    ++numCouldRun_;

    // test if the heuristic can run
    if (!shouldHeurRun_randomChoice() || !numberHeuristics_)
        return 0;
    double randomNumber = randomNumberGenerator_.randomDouble();
    int i;
    for (i = 0; i < numberHeuristics_; i++) {
        if (randomNumber < probabilities_[i])
            break;
    }
    assert (i < numberHeuristics_);
    int returnCode;
    //model_->unsetDivingHasRun();
#ifdef COIN_DEVELOP
    printf("JustOne running %s\n",
           heuristic_[i]->heuristicName());
#endif
    returnCode = heuristic_[i]->solution(solutionValue, betterSolution);
#ifdef COIN_DEVELOP
    if (returnCode)
        printf("JustOne running %s found solution\n",
               heuristic_[i]->heuristicName());
#endif
    return returnCode;
}
// Resets stuff if model changes
void
CbcHeuristicJustOne::resetModel(CbcModel * model)
{
    CbcHeuristic::resetModel(model);
    for (int i = 0; i < numberHeuristics_; i++)
        heuristic_[i]->resetModel(model);
}
// update model (This is needed if cliques update matrix etc)
void
CbcHeuristicJustOne::setModel(CbcModel * model)
{
    CbcHeuristic::setModel(model);
    for (int i = 0; i < numberHeuristics_; i++)
        heuristic_[i]->setModel(model);
}
// Validate model i.e. sets when_ to 0 if necessary (may be NULL)
void
CbcHeuristicJustOne::validate()
{
    CbcHeuristic::validate();
    for (int i = 0; i < numberHeuristics_; i++)
        heuristic_[i]->validate();
}
// Adds an heuristic with probability
void
CbcHeuristicJustOne::addHeuristic(const CbcHeuristic * heuristic, double probability)
{
    CbcHeuristic * thisOne = heuristic->clone();
    thisOne->setWhen(-999);
    CbcHeuristic ** tempH = CoinCopyOfArrayPartial(heuristic_, numberHeuristics_ + 1,
                            numberHeuristics_);
    delete [] heuristic_;
    heuristic_ = tempH;
    heuristic_[numberHeuristics_] = thisOne;
    double * tempP = CoinCopyOfArrayPartial(probabilities_, numberHeuristics_ + 1,
                                            numberHeuristics_);
    delete [] probabilities_;
    probabilities_ = tempP;
    probabilities_[numberHeuristics_] = probability;
    numberHeuristics_++;
}
// Normalize probabilities
void
CbcHeuristicJustOne::normalizeProbabilities()
{
    double sum = 0.0;
    for (int i = 0; i < numberHeuristics_; i++)
        sum += probabilities_[i];
    double multiplier = 1.0 / sum;
    sum = 0.0;
    for (int i = 0; i < numberHeuristics_; i++) {
        sum += probabilities_[i];
        probabilities_[i] = sum * multiplier;
    }
    assert (fabs(probabilities_[numberHeuristics_-1] - 1.0) < 1.0e-5);
    probabilities_[numberHeuristics_-1] = 1.000001;
}