/* $Id: CbcLinked.cpp 1902 2013-04-10 16:58:16Z stefan $ */
// Copyright (C) 2006, International Business Machines
// Corporation and others.  All Rights Reserved.
// This code is licensed under the terms of the Eclipse Public License (EPL).

#include "CbcConfig.h"

#include "CoinTime.hpp"

#include "CoinHelperFunctions.hpp"
#include "CoinModel.hpp"
#include "ClpSimplex.hpp"
// returns jColumn (-2 if linear term, -1 if unknown) and coefficient
static
int decodeBit(char * phrase, char * & nextPhrase, double & coefficient, bool ifFirst, const CoinModel & model)
{
    char * pos = phrase;
    // may be leading - (or +)
    char * pos2 = pos;
    double value = 1.0;
    if (*pos2 == '-' || *pos2 == '+')
        pos2++;
    // next terminator * or + or -
    while (*pos2) {
        if (*pos2 == '*') {
            break;
        } else if (*pos2 == '-' || *pos2 == '+') {
            if (pos2 == pos || *(pos2 - 1) != 'e')
                break;
        }
        pos2++;
    }
    // if * must be number otherwise must be name
    if (*pos2 == '*') {
        char * pos3 = pos;
        while (pos3 != pos2) {
            pos3++;
#ifndef NDEBUG
            char x = *(pos3 - 1);
            assert ((x >= '0' && x <= '9') || x == '.' || x == '+' || x == '-' || x == 'e');
#endif
        }
        char saved = *pos2;
        *pos2 = '\0';
        value = atof(pos);
        *pos2 = saved;
        // and down to next
        pos2++;
        pos = pos2;
        while (*pos2) {
            if (*pos2 == '-' || *pos2 == '+')
                break;
            pos2++;
        }
    }
    char saved = *pos2;
    *pos2 = '\0';
    // now name
    // might have + or -
    if (*pos == '+') {
        pos++;
    } else if (*pos == '-') {
        pos++;
        assert (value == 1.0);
        value = - value;
    }
    int jColumn = model.column(pos);
    // must be column unless first when may be linear term
    if (jColumn < 0) {
        if (ifFirst) {
            char * pos3 = pos;
            while (pos3 != pos2) {
                pos3++;
#ifndef NDEBUG
                char x = *(pos3 - 1);
                assert ((x >= '0' && x <= '9') || x == '.' || x == '+' || x == '-' || x == 'e');
#endif
            }
            assert(*pos2 == '\0');
            // keep possible -
            value = value * atof(pos);
            jColumn = -2;
        } else {
            // bad
            *pos2 = saved;
            printf("bad nonlinear term %s\n", phrase);
            abort();
        }
    }
    *pos2 = saved;
    pos = pos2;
    coefficient = value;
    nextPhrase = pos;
    return jColumn;
}
#include "ClpQuadraticObjective.hpp"
#include <cassert>
#if defined(_MSC_VER)
// Turn off compiler warning about long names
#  pragma warning(disable:4786)
#endif
#include "CbcLinked.hpp"
#include "CoinIndexedVector.hpp"
#include "CoinMpsIO.hpp"
//#include "OsiSolverLink.hpp"
//#include "OsiBranchLink.hpp"
#include "ClpPackedMatrix.hpp"
#include "CoinTime.hpp"
#include "CbcModel.hpp"
#include "CbcCutGenerator.hpp"
#include "CglStored.hpp"
#include "CglPreProcess.hpp"
#include "CglGomory.hpp"
#include "CglProbing.hpp"
#include "CglKnapsackCover.hpp"
#include "CglRedSplit.hpp"
#include "CglClique.hpp"
#include "CglFlowCover.hpp"
#include "CglMixedIntegerRounding2.hpp"
#include "CglTwomir.hpp"
#include "CglDuplicateRow.hpp"
#include "CbcHeuristicFPump.hpp"
#include "CbcHeuristic.hpp"
#include "CbcHeuristicLocal.hpp"
#include "CbcHeuristicGreedy.hpp"
#include "ClpLinearObjective.hpp"
#include "CbcBranchActual.hpp"
#include "CbcCompareActual.hpp"
//#############################################################################
// Solve methods
//#############################################################################
void OsiSolverLink::initialSolve()
{
    //writeMps("yy");
    //exit(77);
    specialOptions_ = 0;
    modelPtr_->setWhatsChanged(0);
    if (numberVariables_) {
        CoinPackedMatrix * temp = new CoinPackedMatrix(*matrix_);
        // update all bounds before coefficients
        for (int i = 0; i < numberVariables_; i++ ) {
            info_[i].updateBounds(modelPtr_);
        }
        int updated = updateCoefficients(modelPtr_, temp);
        if (updated || 1) {
            temp->removeGaps(1.0e-14);
            ClpMatrixBase * save = modelPtr_->clpMatrix();
            ClpPackedMatrix * clpMatrix = dynamic_cast<ClpPackedMatrix *> (save);
            assert (clpMatrix);
            if (save->getNumRows() > temp->getNumRows()) {
                // add in cuts
                int numberRows = temp->getNumRows();
                int * which = new int[numberRows];
                for (int i = 0; i < numberRows; i++)
                    which[i] = i;
                save->deleteRows(numberRows, which);
                delete [] which;
                temp->bottomAppendPackedMatrix(*clpMatrix->matrix());
            }
            modelPtr_->replaceMatrix(temp, true);
        } else {
            delete temp;
        }
    }
    if (0) {
        const double * lower = getColLower();
        const double * upper = getColUpper();
        int n = 0;
        for (int i = 84; i < 84 + 16; i++) {
            if (lower[i] + 0.01 < upper[i]) {
                n++;
            }
        }
        if (!n)
            writeMps("sol_query");

    }
    //static int iPass=0;
    //char temp[50];
    //iPass++;
    //sprintf(temp,"cc%d",iPass);
    //writeMps(temp);
    //writeMps("tight");
    //exit(33);
    //printf("wrote cc%d\n",iPass);
    OsiClpSolverInterface::initialSolve();
    int secondaryStatus = modelPtr_->secondaryStatus();
    if (modelPtr_->status() == 0 && (secondaryStatus == 2 || secondaryStatus == 4))
        modelPtr_->cleanup(1);
    //if (!isProvenOptimal())
    //writeMps("yy");
    if (isProvenOptimal() && quadraticModel_ && modelPtr_->numberColumns() == quadraticModel_->numberColumns()) {
        // see if qp can get better solution
        const double * solution = modelPtr_->primalColumnSolution();
        int numberColumns = modelPtr_->numberColumns();
        bool satisfied = true;
        for (int i = 0; i < numberColumns; i++) {
            if (isInteger(i)) {
                double value = solution[i];
                if (fabs(value - floor(value + 0.5)) > 1.0e-6) {
                    satisfied = false;
                    break;
                }
            }
        }
        if (satisfied) {
            ClpSimplex qpTemp(*quadraticModel_);
            double * lower = qpTemp.columnLower();
            double * upper = qpTemp.columnUpper();
            double * lower2 = modelPtr_->columnLower();
            double * upper2 = modelPtr_->columnUpper();
            for (int i = 0; i < numberColumns; i++) {
                if (isInteger(i)) {
                    double value = floor(solution[i] + 0.5);
                    lower[i] = value;
                    upper[i] = value;
                } else {
                    lower[i] = lower2[i];
                    upper[i] = upper2[i];
                }
            }
            //qpTemp.writeMps("bad.mps");
            //modelPtr_->writeMps("bad2.mps");
            //qpTemp.objectiveAsObject()->setActivated(0);
            //qpTemp.primal();
            //qpTemp.objectiveAsObject()->setActivated(1);
            qpTemp.primal();
            //assert (!qpTemp.problemStatus());
            if (qpTemp.objectiveValue() < bestObjectiveValue_ - 1.0e-3 && !qpTemp.problemStatus()) {
                delete [] bestSolution_;
                bestSolution_ = CoinCopyOfArray(qpTemp.primalColumnSolution(), numberColumns);
                bestObjectiveValue_ = qpTemp.objectiveValue();
                printf("better qp objective of %g\n", bestObjectiveValue_);
                // If model has stored then add cut (if convex)
                if (cbcModel_ && (specialOptions2_&4) != 0) {
                    int numberGenerators = cbcModel_->numberCutGenerators();
                    int iGenerator;
                    cbcModel_->lockThread();
                    for (iGenerator = 0; iGenerator < numberGenerators; iGenerator++) {
                        CbcCutGenerator * generator = cbcModel_->cutGenerator(iGenerator);
                        CglCutGenerator * gen = generator->generator();
                        CglStored * gen2 = dynamic_cast<CglStored *> (gen);
                        if (gen2) {
                            // add OA cut
                            double offset;
                            double * gradient = new double [numberColumns+1];
                            memcpy(gradient, qpTemp.objectiveAsObject()->gradient(&qpTemp, bestSolution_, offset, true, 2),
                                   numberColumns*sizeof(double));
                            // assume convex
                            double rhs = 0.0;
                            int * column = new int[numberColumns+1];
                            int n = 0;
                            for (int i = 0; i < numberColumns; i++) {
                                double value = gradient[i];
                                if (fabs(value) > 1.0e-12) {
                                    gradient[n] = value;
                                    rhs += value * solution[i];
                                    column[n++] = i;
                                }
                            }
                            gradient[n] = -1.0;
                            column[n++] = objectiveVariable_;
                            gen2->addCut(-COIN_DBL_MAX, offset + 1.0e-7, n, column, gradient);
                            delete [] gradient;
                            delete [] column;
                            break;
                        }
                    }
                    cbcModel_->unlockThread();
                }
            }
        }
    }
}
//#define WRITE_MATRIX
#ifdef WRITE_MATRIX
static int xxxxxx = 0;
#endif
//-----------------------------------------------------------------------------
void OsiSolverLink::resolve()
{
    if (false) {
        bool takeHint;
        OsiHintStrength strength;
        // Switch off printing if asked to
        getHintParam(OsiDoReducePrint, takeHint, strength);
        if (strength != OsiHintIgnore && takeHint) {
            printf("no printing\n");
        } else {
            printf("printing\n");
        }
    }
    specialOptions_ = 0;
    modelPtr_->setWhatsChanged(0);
    bool allFixed = numberFix_ > 0;
    bool feasible = true;
    if (numberVariables_) {
        CoinPackedMatrix * temp = new CoinPackedMatrix(*matrix_);
        //bool best=true;
        const double * lower = modelPtr_->columnLower();
        const double * upper = modelPtr_->columnUpper();
        // update all bounds before coefficients
        for (int i = 0; i < numberVariables_; i++ ) {
            info_[i].updateBounds(modelPtr_);
            int iColumn = info_[i].variable();
            double lo = lower[iColumn];
            double up = upper[iColumn];
            if (up < lo)
                feasible = false;
        }
        int updated = updateCoefficients(modelPtr_, temp);
        if (updated) {
            temp->removeGaps(1.0e-14);
            ClpMatrixBase * save = modelPtr_->clpMatrix();
            ClpPackedMatrix * clpMatrix = dynamic_cast<ClpPackedMatrix *> (save);
            assert (clpMatrix);
            if (save->getNumRows() > temp->getNumRows()) {
                // add in cuts
                int numberRows = temp->getNumRows();
                int * which = new int[numberRows];
                for (int i = 0; i < numberRows; i++)
                    which[i] = i;
                CoinPackedMatrix * mat = clpMatrix->matrix();
                // for debug
                //mat = new CoinPackedMatrix(*mat);
                mat->deleteRows(numberRows, which);
                delete [] which;
                temp->bottomAppendPackedMatrix(*mat);
                temp->removeGaps(1.0e-14);
            }
            modelPtr_->replaceMatrix(temp, true);
	    modelPtr_->setNewRowCopy(NULL);
	    modelPtr_->setClpScaledMatrix(NULL);
        } else {
            delete temp;
        }
    }
#ifdef WRITE_MATRIX
    {
        xxxxxx++;
        char temp[50];
        sprintf(temp, "bb%d", xxxxxx);
        writeMps(temp);
        printf("wrote bb%d\n", xxxxxx);
    }
#endif
    if (0) {
        const double * lower = getColLower();
        const double * upper = getColUpper();
        int n = 0;
        for (int i = 60; i < 64; i++) {
            if (lower[i]) {
                printf("%d bounds %g %g\n", i, lower[i], upper[i]);
                n++;
            }
        }
        if (n == 1)
            printf("just one?\n");
    }
    // check feasible
    {
        const double * lower = getColLower();
        const double * upper = getColUpper();
        int numberColumns = getNumCols();
        for (int i = 0; i < numberColumns; i++) {
            if (lower[i] > upper[i] + 1.0e-12) {
                feasible = false;
                break;
            }
        }
    }
    if (!feasible)
        allFixed = false;
    if ((specialOptions2_&1) == 0)
        allFixed = false;
    int returnCode = -1;
    // See if in strong branching
    int maxIts = modelPtr_->maximumIterations();
    if (feasible) {
        if (maxIts > 10000) {
            // may do lots of work
            if ((specialOptions2_&1) != 0) {
                // see if fixed
                const double * lower = modelPtr_->columnLower();
                const double * upper = modelPtr_->columnUpper();
                for (int i = 0; i < numberFix_; i++ ) {
                    int iColumn = fixVariables_[i];
                    double lo = lower[iColumn];
                    double up = upper[iColumn];
                    if (up > lo) {
                        allFixed = false;
                        break;
                    }
                }
                returnCode = allFixed ? fathom(allFixed) : 0;
            } else {
                returnCode = 0;
            }
        } else {
            returnCode = 0;
        }
    }
    if (returnCode >= 0) {
        if (returnCode == 0)
            OsiClpSolverInterface::resolve();
        int satisfied = 2;
        const double * solution = getColSolution();
        const double * lower = getColLower();
        const double * upper = getColUpper();
        int numberColumns2 = coinModel_.numberColumns();
        for (int i = 0; i < numberColumns2; i++) {
            if (isInteger(i)) {
                double value = solution[i];
                if (fabs(value - floor(value + 0.5)) > 1.0e-6) {
                    satisfied = 0;
                    break;
                } else if (upper[i] > lower[i]) {
                    satisfied = 1;
                }
            }
        }
        if (isProvenOptimal()) {
            //if (satisfied==2)
            //printf("satisfied %d\n",satisfied);
            if (satisfied && (specialOptions2_&2) != 0) {
                assert (quadraticModel_);
                // look at true objective
#ifndef NDEBUG
                double direction = modelPtr_->optimizationDirection();
                assert (direction == 1.0);
#endif
                double value = - quadraticModel_->objectiveOffset();
                const double * objective = quadraticModel_->objective();
                int i;
                for ( i = 0; i < numberColumns2; i++)
                    value += solution[i] * objective[i];
                // and now rest
                for ( i = 0; i < numberObjects_; i++) {
                    OsiBiLinear * obj = dynamic_cast<OsiBiLinear *> (object_[i]);
                    if (obj) {
                        value += obj->xyCoefficient(solution);
                    }
                }
                if (value < bestObjectiveValue_ - 1.0e-3) {
                    printf("obj of %g\n", value);
                    //modelPtr_->setDualObjectiveLimit(value);
                    delete [] bestSolution_;
                    bestSolution_ = CoinCopyOfArray(modelPtr_->getColSolution(), modelPtr_->getNumCols());
                    bestObjectiveValue_ = value;
                    if (maxIts <= 10000 && cbcModel_) {
                        OsiSolverLink * solver2 = dynamic_cast<OsiSolverLink *> (cbcModel_->solver());
                        assert (solver2);
                        if (solver2 != this) {
                            // in strong branching - need to store in original solver
                            if (value < solver2->bestObjectiveValue_ - 1.0e-3) {
                                delete [] solver2->bestSolution_;
                                solver2->bestSolution_ = CoinCopyOfArray(bestSolution_, modelPtr_->getNumCols());
                                solver2->bestObjectiveValue_ = value;
                            }
                        }
                    }
                    // If model has stored then add cut (if convex)
                    if (cbcModel_ && (specialOptions2_&4) != 0 && quadraticModel_) {
                        int numberGenerators = cbcModel_->numberCutGenerators();
                        int iGenerator;
                        for (iGenerator = 0; iGenerator < numberGenerators; iGenerator++) {
                            CbcCutGenerator * generator = cbcModel_->cutGenerator(iGenerator);
                            CglCutGenerator * gen = generator->generator();
                            CglStored * gen2 = dynamic_cast<CglStored *> (gen);
                            if (gen2) {
                                cbcModel_->lockThread();
                                // add OA cut
                                double offset = 0.0;
                                int numberColumns = quadraticModel_->numberColumns();
                                double * gradient = new double [numberColumns+1];
                                // gradient from bilinear
                                int i;
                                CoinZeroN(gradient, numberColumns + 1);
                                //const double * objective = modelPtr_->objective();
                                assert (objectiveRow_ >= 0);
                                const double * element = originalRowCopy_->getElements();
                                const int * column2 = originalRowCopy_->getIndices();
                                const CoinBigIndex * rowStart = originalRowCopy_->getVectorStarts();
                                //const int * rowLength = originalRowCopy_->getVectorLengths();
                                //int numberColumns2 = coinModel_.numberColumns();
                                for ( i = rowStart[objectiveRow_]; i < rowStart[objectiveRow_+1]; i++)
                                    gradient[column2[i]] = element[i];
                                for ( i = 0; i < numberObjects_; i++) {
                                    OsiBiLinear * obj = dynamic_cast<OsiBiLinear *> (object_[i]);
                                    if (obj) {
                                        int xColumn = obj->xColumn();
                                        int yColumn = obj->yColumn();
                                        if (xColumn != yColumn) {
                                            double coefficient = /* 2.0* */obj->coefficient();
                                            gradient[xColumn] += coefficient * solution[yColumn];
                                            gradient[yColumn] += coefficient * solution[xColumn];
                                            offset += coefficient * solution[xColumn] * solution[yColumn];
                                        } else {
                                            double coefficient = obj->coefficient();
                                            gradient[xColumn] += 2.0 * coefficient * solution[yColumn];
                                            offset += coefficient * solution[xColumn] * solution[yColumn];
                                        }
                                    }
                                }
                                // assume convex
                                double rhs = 0.0;
                                int * column = new int[numberColumns+1];
                                int n = 0;
                                for (int i = 0; i < numberColumns; i++) {
                                    double value = gradient[i];
                                    if (fabs(value) > 1.0e-12) {
                                        gradient[n] = value;
                                        rhs += value * solution[i];
                                        column[n++] = i;
                                    }
                                }
                                gradient[n] = -1.0;
                                column[n++] = objectiveVariable_;
                                gen2->addCut(-COIN_DBL_MAX, offset + 1.0e-4, n, column, gradient);
                                delete [] gradient;
                                delete [] column;
                                cbcModel_->unlockThread();
                                break;
                            }
                        }
                    }
                }
            } else if (satisfied == 2) {
                // is there anything left to do?
                int i;
                int numberContinuous = 0;
                double gap = 0.0;
                for ( i = 0; i < numberObjects_; i++) {
                    OsiBiLinear * obj = dynamic_cast<OsiBiLinear *> (object_[i]);
                    if (obj) {
                        if (obj->xMeshSize() < 1.0 && obj->yMeshSize() < 1.0) {
                            numberContinuous++;
                            int xColumn = obj->xColumn();
                            double gapX = upper[xColumn] - lower[xColumn];
                            int yColumn = obj->yColumn();
                            double gapY = upper[yColumn] - lower[yColumn];
                            gap = CoinMax(gap, CoinMax(gapX, gapY));
                        }
                    }
                }
                if (numberContinuous && 0) {
                    // iterate to get solution and fathom node
                    int numberColumns2 = coinModel_.numberColumns();
                    double * lower2 = CoinCopyOfArray(getColLower(), numberColumns2);
                    double * upper2 = CoinCopyOfArray(getColUpper(), numberColumns2);
                    while (gap > defaultMeshSize_) {
                        gap *= 0.9;
                        const double * solution = getColSolution();
                        double newGap = 0.0;
                        for ( i = 0; i < numberObjects_; i++) {
                            OsiBiLinear * obj = dynamic_cast<OsiBiLinear *> (object_[i]);
                            if (obj && (obj->branchingStrategy()&8) == 0) {
                                if (obj->xMeshSize() < 1.0 && obj->yMeshSize() < 1.0) {
                                    numberContinuous++;
                                    // need to make sure new xy value in range
                                    double xB[3];
                                    double yB[3];
                                    double xybar[4];
                                    obj->getCoefficients(this, xB, yB, xybar);
                                    //double xyTrue = obj->xyCoefficient(solution);
                                    double xyLambda = 0.0;
                                    int firstLambda = obj->firstLambda();
                                    for (int j = 0; j < 4; j++) {
                                        xyLambda += solution[firstLambda+j] * xybar[j];
                                    }
                                    //printf ("x %d, y %d - true %g lambda %g\n",obj->xColumn(),obj->yColumn(),
                                    //  xyTrue,xyLambda);
                                    int xColumn = obj->xColumn();
                                    double gapX = upper[xColumn] - lower[xColumn];
                                    int yColumn = obj->yColumn();
                                    if (gapX > gap) {
                                        double value = solution[xColumn];
                                        double newLower = CoinMax(lower2[xColumn], value - 0.5 * gap);
                                        double newUpper = CoinMin(upper2[xColumn], value + 0.5 * gap);
                                        if (newUpper - newLower < 0.99*gap) {
                                            if (newLower == lower2[xColumn])
                                                newUpper = CoinMin(upper2[xColumn], newLower + gap);
                                            else if (newUpper == upper2[xColumn])
                                                newLower = CoinMax(lower2[xColumn], newUpper - gap);
                                        }
                                        // see if problem
#ifndef NDEBUG
                                        double lambda[4];
#endif
                                        xB[0] = newLower;
                                        xB[1] = newUpper;
                                        xB[2] = value;
                                        yB[2] = solution[yColumn];
                                        xybar[0] = xB[0] * yB[0];
                                        xybar[1] = xB[0] * yB[1];
                                        xybar[2] = xB[1] * yB[0];
                                        xybar[3] = xB[1] * yB[1];
#ifndef NDEBUG
                                        double infeasibility = obj->computeLambdas(xB, yB, xybar, lambda);
                                        assert (infeasibility < 1.0e-9);
#endif
                                        setColLower(xColumn, newLower);
                                        setColUpper(xColumn, newUpper);
                                    }
                                    double gapY = upper[yColumn] - lower[yColumn];
                                    if (gapY > gap) {
                                        double value = solution[yColumn];
                                        double newLower = CoinMax(lower2[yColumn], value - 0.5 * gap);
                                        double newUpper = CoinMin(upper2[yColumn], value + 0.5 * gap);
                                        if (newUpper - newLower < 0.99*gap) {
                                            if (newLower == lower2[yColumn])
                                                newUpper = CoinMin(upper2[yColumn], newLower + gap);
                                            else if (newUpper == upper2[yColumn])
                                                newLower = CoinMax(lower2[yColumn], newUpper - gap);
                                        }
                                        // see if problem
#ifndef NDEBUG
                                        double lambda[4];
#endif
                                        yB[0] = newLower;
                                        yB[1] = newUpper;
                                        xybar[0] = xB[0] * yB[0];
                                        xybar[1] = xB[0] * yB[1];
                                        xybar[2] = xB[1] * yB[0];
                                        xybar[3] = xB[1] * yB[1];
#ifndef NDEBUG
                                        double infeasibility = obj->computeLambdas(xB, yB, xybar, lambda);
                                        assert (infeasibility < 1.0e-9);
#endif
                                        setColLower(yColumn, newLower);
                                        setColUpper(yColumn, newUpper);
                                    }
                                    newGap = CoinMax(newGap, CoinMax(gapX, gapY));
                                }
                            }
                        }
                        printf("solving with gap of %g\n", gap);
                        //OsiClpSolverInterface::resolve();
                        initialSolve();
                        if (!isProvenOptimal())
                            break;
                    }
                    delete [] lower2;
                    delete [] upper2;
                    //if (isProvenOptimal())
                    //writeMps("zz");
                }
            }
            // ???  - try
            // But skip if strong branching
            CbcModel * cbcModel = (modelPtr_->maximumIterations() > 10000) ? cbcModel_ : NULL;
            if ((specialOptions2_&2) != 0) {
                // If model has stored then add cut (if convex)
                // off until I work out problem with ibell3a
                if (cbcModel && (specialOptions2_&4) != 0 && quadraticModel_) {
                    int numberGenerators = cbcModel_->numberCutGenerators();
                    int iGenerator;
                    for (iGenerator = 0; iGenerator < numberGenerators; iGenerator++) {
                        CbcCutGenerator * generator = cbcModel_->cutGenerator(iGenerator);
                        CglCutGenerator * gen = generator->generator();
                        CglTemporary * gen2 = dynamic_cast<CglTemporary *> (gen);
                        if (gen2) {
                            double * solution2 = NULL;
                            int numberColumns = quadraticModel_->numberColumns();
                            int depth = cbcModel_->currentNode() ? cbcModel_->currentNode()->depth() : 0;
                            if (depth < 5) {
                                ClpSimplex qpTemp(*quadraticModel_);
                                double * lower = qpTemp.columnLower();
                                double * upper = qpTemp.columnUpper();
                                double * lower2 = modelPtr_->columnLower();
                                double * upper2 = modelPtr_->columnUpper();
                                for (int i = 0; i < numberColumns; i++) {
                                    lower[i] = lower2[i];
                                    upper[i] = upper2[i];
                                }
                                qpTemp.setLogLevel(modelPtr_->logLevel());
                                qpTemp.primal();
                                assert (!qpTemp.problemStatus());
                                if (qpTemp.objectiveValue() < bestObjectiveValue_ - 1.0e-3 && !qpTemp.problemStatus()) {
                                    solution2 = CoinCopyOfArray(qpTemp.primalColumnSolution(), numberColumns);
                                } else {
                                    printf("QP says expensive - kill\n");
                                    modelPtr_->setProblemStatus(1);
                                    modelPtr_->setObjectiveValue(COIN_DBL_MAX);
                                    break;
                                }
                            }
                            const double * solution = getColSolution();
                            // add OA cut
                            doAOCuts(gen2, solution, solution);
                            if (solution2) {
                                doAOCuts(gen2, solution, solution2);
                                delete [] solution2;
                            }
                            break;
                        }
                    }
                }
            } else if (cbcModel && (specialOptions2_&8) == 8) {
                // convex and nonlinear in constraints
                int numberGenerators = cbcModel_->numberCutGenerators();
                int iGenerator;
                for (iGenerator = 0; iGenerator < numberGenerators; iGenerator++) {
                    CbcCutGenerator * generator = cbcModel_->cutGenerator(iGenerator);
                    CglCutGenerator * gen = generator->generator();
                    CglTemporary * gen2 = dynamic_cast<CglTemporary *> (gen);
                    if (gen2) {
                        const double * solution = getColSolution();
                        const double * rowUpper = getRowUpper();
                        const double * rowLower = getRowLower();
                        const double * element = originalRowCopy_->getElements();
                        const int * column2 = originalRowCopy_->getIndices();
                        const CoinBigIndex * rowStart = originalRowCopy_->getVectorStarts();
                        //const int * rowLength = originalRowCopy_->getVectorLengths();
                        int numberColumns2 = CoinMax(coinModel_.numberColumns(), objectiveVariable_ + 1);
                        double * gradient = new double [numberColumns2];
                        int * column = new int[numberColumns2];
                        //const double * columnLower = modelPtr_->columnLower();
                        //const double * columnUpper = modelPtr_->columnUpper();
                        cbcModel_->lockThread();
                        for (int iNon = 0; iNon < numberNonLinearRows_; iNon++) {
                            int iRow = rowNonLinear_[iNon];
                            bool convex = convex_[iNon] > 0;
                            if (!convex_[iNon])
                                continue; // can't use this row
                            // add OA cuts
                            double offset = 0.0;
                            // gradient from bilinear
                            int i;
                            CoinZeroN(gradient, numberColumns2);
                            //const double * objective = modelPtr_->objective();
                            for ( i = rowStart[iRow]; i < rowStart[iRow+1]; i++)
                                gradient[column2[i]] = element[i];
                            for ( i = startNonLinear_[iNon]; i < startNonLinear_[iNon+1]; i++) {
                                OsiBiLinear * obj = dynamic_cast<OsiBiLinear *> (object_[whichNonLinear_[i]]);
                                assert (obj);
                                int xColumn = obj->xColumn();
                                int yColumn = obj->yColumn();
                                if (xColumn != yColumn) {
                                    double coefficient = /* 2.0* */obj->coefficient();
                                    gradient[xColumn] += coefficient * solution[yColumn];
                                    gradient[yColumn] += coefficient * solution[xColumn];
                                    offset += coefficient * solution[xColumn] * solution[yColumn];
                                } else {
                                    double coefficient = obj->coefficient();
                                    gradient[xColumn] += 2.0 * coefficient * solution[yColumn];
                                    offset += coefficient * solution[xColumn] * solution[yColumn];
                                }
                            }
                            // assume convex
                            double rhs = 0.0;
                            int n = 0;
                            for (int i = 0; i < numberColumns2; i++) {
                                double value = gradient[i];
                                if (fabs(value) > 1.0e-12) {
                                    gradient[n] = value;
                                    rhs += value * solution[i];
                                    column[n++] = i;
                                }
                            }
                            if (iRow == objectiveRow_) {
                                gradient[n] = -1.0;
                                assert (objectiveVariable_ >= 0);
                                rhs -= solution[objectiveVariable_];
                                column[n++] = objectiveVariable_;
                                assert (convex);
                            } else if (convex) {
                                offset += rowUpper[iRow];
                            } else if (!convex) {
                                offset += rowLower[iRow];
                            }
                            if (convex && rhs > offset + 1.0e-5)
                                gen2->addCut(-COIN_DBL_MAX, offset + 1.0e-7, n, column, gradient);
                            else if (!convex && rhs < offset - 1.0e-5)
                                gen2->addCut(offset - 1.0e-7, COIN_DBL_MAX, n, column, gradient);
                        }
                        cbcModel_->unlockThread();
                        delete [] gradient;
                        delete [] column;
                        break;
                    }
                }
            }
        }
    } else {
        modelPtr_->setProblemStatus(1);
        modelPtr_->setObjectiveValue(COIN_DBL_MAX);
    }
}
// Do OA cuts
int
OsiSolverLink::doAOCuts(CglTemporary * cutGen, const double * solution, const double * solution2)
{
    cbcModel_->lockThread();
    // add OA cut
    double offset = 0.0;
    int numberColumns = quadraticModel_->numberColumns();
    double * gradient = new double [numberColumns+1];
    // gradient from bilinear
    int i;
    CoinZeroN(gradient, numberColumns + 1);
    //const double * objective = modelPtr_->objective();
    assert (objectiveRow_ >= 0);
    const double * element = originalRowCopy_->getElements();
    const int * column2 = originalRowCopy_->getIndices();
    const CoinBigIndex * rowStart = originalRowCopy_->getVectorStarts();
    //const int * rowLength = originalRowCopy_->getVectorLengths();
    //int numberColumns2 = coinModel_.numberColumns();
    for ( i = rowStart[objectiveRow_]; i < rowStart[objectiveRow_+1]; i++)
        gradient[column2[i]] = element[i];
    //const double * columnLower = modelPtr_->columnLower();
    //const double * columnUpper = modelPtr_->columnUpper();
    for ( i = 0; i < numberObjects_; i++) {
        OsiBiLinear * obj = dynamic_cast<OsiBiLinear *> (object_[i]);
        if (obj) {
            int xColumn = obj->xColumn();
            int yColumn = obj->yColumn();
            if (xColumn != yColumn) {
                double coefficient = /* 2.0* */obj->coefficient();
                gradient[xColumn] += coefficient * solution2[yColumn];
                gradient[yColumn] += coefficient * solution2[xColumn];
                offset += coefficient * solution2[xColumn] * solution2[yColumn];
            } else {
                double coefficient = obj->coefficient();
                gradient[xColumn] += 2.0 * coefficient * solution2[yColumn];
                offset += coefficient * solution2[xColumn] * solution2[yColumn];
            }
        }
    }
    // assume convex
    double rhs = 0.0;
    int * column = new int[numberColumns+1];
    int n = 0;
    for (int i = 0; i < numberColumns; i++) {
        double value = gradient[i];
        if (fabs(value) > 1.0e-12) {
            gradient[n] = value;
            rhs += value * solution[i];
            column[n++] = i;
        }
    }
    gradient[n] = -1.0;
    assert (objectiveVariable_ >= 0);
    rhs -= solution[objectiveVariable_];
    column[n++] = objectiveVariable_;
    int returnCode = 0;
    if (rhs > offset + 1.0e-5) {
        cutGen->addCut(-COIN_DBL_MAX, offset + 1.0e-7, n, column, gradient);
        //printf("added cut with %d elements\n",n);
        returnCode = 1;
    }
    delete [] gradient;
    delete [] column;
    cbcModel_->unlockThread();
    return returnCode;
}

//#############################################################################
// Constructors, destructors clone and assignment
//#############################################################################

//-------------------------------------------------------------------
// Default Constructor
//-------------------------------------------------------------------
OsiSolverLink::OsiSolverLink ()
        : CbcOsiSolver()
{
    gutsOfDestructor(true);
}
#ifdef JJF_ZERO
/* returns
   sequence of nonlinear or
   -1 numeric
   -2 not found
   -3 too many terms
*/
static int getVariable(const CoinModel & model, char * expression,
                       int & linear)
{
    int non = -1;
    linear = -1;
    if (strcmp(expression, "Numeric")) {
        // function
        char * first = strchr(expression, '*');
        int numberColumns = model.numberColumns();
        int j;
        if (first) {
            *first = '\0';
            for (j = 0; j < numberColumns; j++) {
                if (!strcmp(expression, model.columnName(j))) {
                    linear = j;
                    memmove(expression, first + 1, strlen(first + 1) + 1);
                    break;
                }
            }
        }
        // find nonlinear
        for (j = 0; j < numberColumns; j++) {
            const char * name = model.columnName(j);
            first = strstr(expression, name);
            if (first) {
                if (first != expression && isalnum(*(first - 1)))
                    continue; // not real match
                first += strlen(name);
                if (!isalnum(*first)) {
                    // match
                    non = j;
                    // but check no others
                    j++;
                    for (; j < numberColumns; j++) {
                        const char * name = model.columnName(j);
                        first = strstr(expression, name);
                        if (first) {
                            if (isalnum(*(first - 1)))
                                continue; // not real match
                            first += strlen(name);
                            if (!isalnum(*first)) {
                                // match - ouch
                                non = -3;
                                break;
                            }
                        }
                    }
                    break;
                }
            }
        }
        if (non == -1)
            non = -2;
    }
    return non;
}
#endif
/* This creates from a coinModel object

   if errors.then number of sets is -1

   This creates linked ordered sets information.  It assumes -

   for product terms syntax is yy*f(zz)
   also just f(zz) is allowed
   and even a constant

   modelObject not const as may be changed as part of process.
*/
OsiSolverLink::OsiSolverLink ( CoinModel & coinModel)
        : CbcOsiSolver()
{
    gutsOfDestructor(true);
    load(coinModel);
}
// need bounds
static void fakeBounds(OsiSolverInterface * solver, int column, double maximumValue,
                       CoinModel * model1, CoinModel * model2)
{
    double lo = solver->getColLower()[column];
    if (lo < -maximumValue) {
        solver->setColLower(column, -maximumValue);
        if (model1)
            model1->setColLower(column, -maximumValue);
        if (model2)
            model2->setColLower(column, -maximumValue);
    }
    double up = solver->getColUpper()[column];
    if (up > maximumValue) {
        solver->setColUpper(column, maximumValue);
        if (model1)
            model1->setColUpper(column, maximumValue);
        if (model2)
            model2->setColUpper(column, maximumValue);
    }
}
void OsiSolverLink::load ( CoinModel & coinModelOriginal, bool tightenBounds, int logLevel)
{
    // first check and set up arrays
    int numberColumns = coinModelOriginal.numberColumns();
    int numberRows = coinModelOriginal.numberRows();
    // List of nonlinear entries
    int * which = new int[numberColumns];
    numberVariables_ = 0;
    //specialOptions2_=0;
    int iColumn;
    int numberErrors = 0;
    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
        CoinModelLink triple = coinModelOriginal.firstInColumn(iColumn);
        bool linear = true;
        int n = 0;
        // See if quadratic objective
        const char * expr = coinModelOriginal.getColumnObjectiveAsString(iColumn);
        if (strcmp(expr, "Numeric")) {
            linear = false;
        }
        while (triple.row() >= 0) {
            int iRow = triple.row();
            const char * expr = coinModelOriginal.getElementAsString(iRow, iColumn);
            if (strcmp(expr, "Numeric")) {
                linear = false;
            }
            triple = coinModelOriginal.next(triple);
            n++;
        }
        if (!linear) {
            which[numberVariables_++] = iColumn;
        }
    }
    // return if nothing
    if (!numberVariables_) {
        delete [] which;
        coinModel_ = coinModelOriginal;
        int nInt = 0;
        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
            if (coinModel_.isInteger(iColumn))
                nInt++;
        }
        printf("There are %d integers\n", nInt);
        loadFromCoinModel(coinModelOriginal, true);
        OsiObject ** objects = new OsiObject * [nInt];
        nInt = 0;
        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
            if (coinModel_.isInteger(iColumn)) {
                objects[nInt] = new OsiSimpleInteger(this, iColumn);
                objects[nInt]->setPriority(integerPriority_);
                nInt++;
            }
        }
        addObjects(nInt, objects);
        int i;
        for (i = 0; i < nInt; i++)
            delete objects[i];
        delete [] objects;
        return;
    } else {
        coinModel_ = coinModelOriginal;
        // arrays for tightening bounds
        int * freeRow = new int [numberRows];
        CoinZeroN(freeRow, numberRows);
        int * tryColumn = new int [numberColumns];
        CoinZeroN(tryColumn, numberColumns);
        int nBi = 0;
        int numberQuadratic = 0;
        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
            // See if quadratic objective
            const char * expr = coinModel_.getColumnObjectiveAsString(iColumn);
            if (strcmp(expr, "Numeric")) {
                // check if value*x+-value*y....
                assert (strlen(expr) < 20000);
                tryColumn[iColumn] = 1;
                char temp[20000];
                strcpy(temp, expr);
                char * pos = temp;
                bool ifFirst = true;
                double linearTerm = 0.0;
                while (*pos) {
                    double value;
                    int jColumn = decodeBit(pos, pos, value, ifFirst, coinModel_);
                    // must be column unless first when may be linear term
                    if (jColumn >= 0) {
                        tryColumn[jColumn] = 1;
                        numberQuadratic++;
                        nBi++;
                    } else if (jColumn == -2) {
                        linearTerm = value;
                    } else {
                        printf("bad nonlinear term %s\n", temp);
                        abort();
                    }
                    ifFirst = false;
                }
                coinModelOriginal.setObjective(iColumn, linearTerm);
            }
        }
        int iRow;
        int saveNBi = nBi;
        for (iRow = 0; iRow < numberRows; iRow++) {
            CoinModelLink triple = coinModel_.firstInRow(iRow);
            while (triple.column() >= 0) {
                int iColumn = triple.column();
                const char *  el = coinModel_.getElementAsString(iRow, iColumn);
                if (strcmp("Numeric", el)) {
                    // check if value*x+-value*y....
                    assert (strlen(el) < 20000);
                    char temp[20000];
                    strcpy(temp, el);
                    char * pos = temp;
                    bool ifFirst = true;
                    double linearTerm = 0.0;
                    tryColumn[iColumn] = 1;
                    freeRow[iRow] = 1;
                    while (*pos) {
                        double value;
                        int jColumn = decodeBit(pos, pos, value, ifFirst, coinModel_);
                        // must be column unless first when may be linear term
                        if (jColumn >= 0) {
                            tryColumn[jColumn] = 1;
                            nBi++;
                        } else if (jColumn == -2) {
                            linearTerm = value;
                        } else {
                            printf("bad nonlinear term %s\n", temp);
                            abort();
                        }
                        ifFirst = false;
                    }
                    coinModelOriginal.setElement(iRow, iColumn, linearTerm);
                }
                triple = coinModel_.next(triple);
            }
        }
        if (!nBi)
            exit(1);
        bool quadraticObjective = false;
        int nInt = 0;
        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
            if (coinModel_.isInteger(iColumn))
                nInt++;
        }
        printf("There are %d bilinear and %d integers\n", nBi, nInt);
        loadFromCoinModel(coinModelOriginal, true);
        CoinModel coinModel = coinModelOriginal;
        if (tightenBounds && numberColumns < 100) {
            // first fake bounds
            for (iColumn = 0; iColumn < numberColumns; iColumn++) {
                if (tryColumn[iColumn]) {
                    fakeBounds(this, iColumn, defaultBound_, &coinModel, &coinModel_);
                }
            }
            ClpSimplex tempModel(*modelPtr_);
            int nDelete = 0;
            for (iRow = 0; iRow < numberRows; iRow++) {
                if (freeRow[iRow])
                    freeRow[nDelete++] = iRow;
            }
            tempModel.deleteRows(nDelete, freeRow);
            tempModel.setOptimizationDirection(1.0);
            if (logLevel < 3) {
                tempModel.setLogLevel(0);
                tempModel.setSpecialOptions(32768);
            }
            double * objective = tempModel.objective();
            CoinZeroN(objective, numberColumns);
            // now up and down
            double * columnLower = modelPtr_->columnLower();
            double * columnUpper = modelPtr_->columnUpper();
            const double * solution = tempModel.primalColumnSolution();
            for (iColumn = 0; iColumn < numberColumns; iColumn++) {
                if (tryColumn[iColumn]) {
                    objective[iColumn] = 1.0;
                    tempModel.primal(1);
                    if (solution[iColumn] > columnLower[iColumn] + 1.0e-3) {
                        double value = solution[iColumn];
                        if (coinModel_.isInteger(iColumn))
                            value = ceil(value - 0.9e-3);
                        if (logLevel > 1)
                            printf("lower bound on %d changed from %g to %g\n", iColumn, columnLower[iColumn], value);
                        columnLower[iColumn] = value;
                        coinModel_.setColumnLower(iColumn, value);
                        coinModel.setColumnLower(iColumn, value);
                    }
                    objective[iColumn] = -1.0;
                    tempModel.primal(1);
                    if (solution[iColumn] < columnUpper[iColumn] - 1.0e-3) {
                        double value = solution[iColumn];
                        if (coinModel_.isInteger(iColumn))
                            value = floor(value + 0.9e-3);
                        if (logLevel > 1)
                            printf("upper bound on %d changed from %g to %g\n", iColumn, columnUpper[iColumn], value);
                        columnUpper[iColumn] = value;
                        coinModel_.setColumnUpper(iColumn, value);
                        coinModel.setColumnUpper(iColumn, value);
                    }
                    objective[iColumn] = 0.0;
                }
            }
        }
        delete [] freeRow;
        delete [] tryColumn;
        CoinBigIndex * startQuadratic = NULL;
        int * columnQuadratic = NULL;
        double * elementQuadratic = NULL;
        if ( saveNBi == nBi) {
            printf("all bilinearity in objective\n");
            specialOptions2_ |= 2;
            quadraticObjective = true;
            // save copy as quadratic model
            quadraticModel_ = new ClpSimplex(*modelPtr_);
            startQuadratic = new CoinBigIndex [numberColumns+1];
            columnQuadratic = new int [numberQuadratic];
            elementQuadratic = new double [numberQuadratic];
            numberQuadratic = 0;
        }
        //if (quadraticObjective||((specialOptions2_&8)!=0&&saveNBi)) {
        if (saveNBi) {
            // add in objective as constraint
            objectiveVariable_ = numberColumns;
            objectiveRow_ = coinModel.numberRows();
            coinModel.addColumn(0, NULL, NULL, -COIN_DBL_MAX, COIN_DBL_MAX, 1.0);
            int * column = new int[numberColumns+1];
            double * element = new double[numberColumns+1];
            double * objective = coinModel.objectiveArray();
            int n = 0;
            for (int i = 0; i < numberColumns; i++) {
                if (objective[i]) {
                    column[n] = i;
                    element[n++] = objective[i];
                    objective[i] = 0.0;
                }
            }
            column[n] = objectiveVariable_;
            element[n++] = -1.0;
            double offset = - coinModel.objectiveOffset();
            //assert (!offset); // get sign right if happens
            printf("***** offset %g\n", offset);
            coinModel.setObjectiveOffset(0.0);
            double lowerBound = -COIN_DBL_MAX;
            coinModel.addRow(n, column, element, lowerBound, offset);
            delete [] column;
            delete [] element;
        }
        OsiObject ** objects = new OsiObject * [nBi+nInt];
        char * marked = new char [numberColumns];
        memset(marked, 0, numberColumns);
        // statistics I-I I-x x-x
        int stats[3] = {0, 0, 0};
        double * sort = new double [nBi];
        nBi = nInt;
        const OsiObject ** justBi = const_cast<const OsiObject **> (objects + nInt);
        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
            if (quadraticObjective)
                startQuadratic[iColumn] = numberQuadratic;
            // See if quadratic objective
            const char * expr = coinModel_.getColumnObjectiveAsString(iColumn);
            if (strcmp(expr, "Numeric")) {
                // need bounds
                fakeBounds(this, iColumn, defaultBound_, &coinModel, &coinModel_);
                // value*x*y
                char temp[20000];
                strcpy(temp, expr);
                char * pos = temp;
                bool ifFirst = true;
                while (*pos) {
                    double value;
                    int jColumn = decodeBit(pos, pos, value, ifFirst, coinModel_);
                    // must be column unless first when may be linear term
                    if (jColumn >= 0) {
                        if (quadraticObjective) {
                            columnQuadratic[numberQuadratic] = jColumn;
                            if (jColumn == iColumn)
                                elementQuadratic[numberQuadratic++] = 2.0 * value; // convention
                            else
                                elementQuadratic[numberQuadratic++] = 1.0 * value; // convention
                        }
                        // need bounds
                        fakeBounds(this, jColumn, defaultBound_, &coinModel, &coinModel_);
                        double meshI = coinModel_.isInteger(iColumn) ? 1.0 : 0.0;
                        if (meshI)
                            marked[iColumn] = 1;
                        double meshJ = coinModel_.isInteger(jColumn) ? 1.0 : 0.0;
                        if (meshJ)
                            marked[jColumn] = 1;
                        // stats etc
                        if (meshI) {
                            if (meshJ)
                                stats[0]++;
                            else
                                stats[1]++;
                        } else {
                            if (meshJ)
                                stats[1]++;
                            else
                                stats[2]++;
                        }
                        if (iColumn <= jColumn)
                            sort[nBi-nInt] = iColumn + numberColumns * jColumn;
                        else
                            sort[nBi-nInt] = jColumn + numberColumns * iColumn;
                        if (!meshJ && !meshI) {
                            meshI = defaultMeshSize_;
                            meshJ = 0.0;
                        }
                        OsiBiLinear * newObj = new OsiBiLinear(&coinModel, iColumn, jColumn, objectiveRow_, value, meshI, meshJ,
                                                               nBi - nInt, justBi);
                        newObj->setPriority(biLinearPriority_);
                        objects[nBi++] = newObj;
                    } else if (jColumn == -2) {
                    } else {
                        printf("bad nonlinear term %s\n", temp);
                        abort();
                    }
                    ifFirst = false;
                }
            }
        }
        // stats
        printf("There were %d I-I, %d I-x and %d x-x bilinear in objective\n",
               stats[0], stats[1], stats[2]);
        if (quadraticObjective) {
            startQuadratic[numberColumns] = numberQuadratic;
            quadraticModel_->loadQuadraticObjective(numberColumns, startQuadratic,
                                                    columnQuadratic, elementQuadratic);
            delete [] startQuadratic;
            delete [] columnQuadratic;
            delete [] elementQuadratic;
        }
        for (iRow = 0; iRow < numberRows; iRow++) {
            CoinModelLink triple = coinModel_.firstInRow(iRow);
            while (triple.column() >= 0) {
                int iColumn = triple.column();
                const char *  el = coinModel_.getElementAsString(iRow, iColumn);
                if (strcmp("Numeric", el)) {
                    // need bounds
                    fakeBounds(this, iColumn, defaultBound_, &coinModel, &coinModel_);
                    // value*x*y
                    char temp[20000];
                    strcpy(temp, el);
                    char * pos = temp;
                    bool ifFirst = true;
                    while (*pos) {
                        double value;
                        int jColumn = decodeBit(pos, pos, value, ifFirst, coinModel_);
                        // must be column unless first when may be linear term
                        if (jColumn >= 0) {
                            // need bounds
                            fakeBounds(this, jColumn, defaultBound_, &coinModel, &coinModel_);
                            double meshI = coinModel_.isInteger(iColumn) ? 1.0 : 0.0;
                            if (meshI)
                                marked[iColumn] = 1;
                            double meshJ = coinModel_.isInteger(jColumn) ? 1.0 : 0.0;
                            if (meshJ)
                                marked[jColumn] = 1;
                            // stats etc
                            if (meshI) {
                                if (meshJ)
                                    stats[0]++;
                                else
                                    stats[1]++;
                            } else {
                                if (meshJ)
                                    stats[1]++;
                                else
                                    stats[2]++;
                            }
                            if (iColumn <= jColumn)
                                sort[nBi-nInt] = iColumn + numberColumns * jColumn;
                            else
                                sort[nBi-nInt] = jColumn + numberColumns * iColumn;
                            if (!meshJ && !meshI) {
                                meshI = defaultMeshSize_;
                                meshJ = 0.0;
                            }
                            OsiBiLinear * newObj = new OsiBiLinear(&coinModel, iColumn, jColumn, iRow, value, meshI, meshJ,
                                                                   nBi - nInt, justBi);
                            newObj->setPriority(biLinearPriority_);
                            objects[nBi++] = newObj;
                        } else if (jColumn == -2) {
                        } else {
                            printf("bad nonlinear term %s\n", temp);
                            abort();
                        }
                        ifFirst = false;
                    }
                }
                triple = coinModel_.next(triple);
            }
        }
        {
            // stats
            std::sort(sort, sort + nBi - nInt);
            int nDiff = 0;
            double last = -1.0;
            for (int i = 0; i < nBi - nInt; i++) {
                if (sort[i] != last)
                    nDiff++;
                last = sort[i];
            }
            delete [] sort;
            printf("There were %d I-I, %d I-x and %d x-x bilinear in total of which %d were duplicates\n",
                   stats[0], stats[1], stats[2], nBi - nInt - nDiff);
        }
        // reload with all bilinear stuff
        loadFromCoinModel(coinModel, true);
        //exit(77);
        nInt = 0;
        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
            if (coinModel_.isInteger(iColumn)) {
                objects[nInt] = new OsiSimpleInteger(this, iColumn);
                if (marked[iColumn])
                    objects[nInt]->setPriority(integerPriority_);
                else
                    objects[nInt]->setPriority(integerPriority_);
                nInt++;
            }
        }
        nInt = nBi;
        delete [] marked;
        if (numberErrors) {
            // errors
            gutsOfDestructor();
            numberVariables_ = -1;
        } else {
            addObjects(nInt, objects);
            int i;
            for (i = 0; i < nInt; i++)
                delete objects[i];
            delete [] objects;
            // Now do dummy bound stuff
            matrix_ = new CoinPackedMatrix(*getMatrixByCol());
            info_ = new OsiLinkedBound [numberVariables_];
            for ( i = 0; i < numberVariables_; i++) {
                info_[i] = OsiLinkedBound(this, which[i], 0, NULL, NULL, NULL);
            }
            // Do row copy but just part
            int numberRows2 = objectiveRow_ >= 0 ? numberRows + 1 : numberRows;
            int * whichRows = new int [numberRows2];
            int * whichColumns = new int [numberColumns];
            CoinIotaN(whichRows, numberRows2, 0);
            CoinIotaN(whichColumns, numberColumns, 0);
            originalRowCopy_ = new CoinPackedMatrix(*getMatrixByRow(),
                                                    numberRows2, whichRows,
                                                    numberColumns, whichColumns);
            delete [] whichColumns;
            numberNonLinearRows_ = 0;
            CoinZeroN(whichRows, numberRows2);
            for ( i = 0; i < numberObjects_; i++) {
                OsiBiLinear * obj = dynamic_cast<OsiBiLinear *> (object_[i]);
                if (obj) {
                    int xyRow = obj->xyRow();
                    assert (xyRow >= 0 && xyRow < numberRows2); // even if obj we should move
                    whichRows[xyRow]++;
                }
            }
            int * pos = new int [numberRows2];
            int n = 0;
            for (i = 0; i < numberRows2; i++) {
                if (whichRows[i]) {
                    pos[numberNonLinearRows_] = n;
                    n += whichRows[i];
                    whichRows[i] = numberNonLinearRows_;
                    numberNonLinearRows_++;
                } else {
                    whichRows[i] = -1;
                }
            }
            startNonLinear_ = new int [numberNonLinearRows_+1];
            memcpy(startNonLinear_, pos, numberNonLinearRows_*sizeof(int));
            startNonLinear_[numberNonLinearRows_] = n;
            rowNonLinear_ = new int [numberNonLinearRows_];
            convex_ = new int [numberNonLinearRows_];
            // do row numbers now
            numberNonLinearRows_ = 0;
            for (i = 0; i < numberRows2; i++) {
                if (whichRows[i] >= 0) {
                    rowNonLinear_[numberNonLinearRows_++] = i;
                }
            }
            whichNonLinear_ = new int [n];
            for ( i = 0; i < numberObjects_; i++) {
                OsiBiLinear * obj = dynamic_cast<OsiBiLinear *> (object_[i]);
                if (obj) {
                    int xyRow = obj->xyRow();
                    int k = whichRows[xyRow];
                    int put = pos[k];
                    pos[k]++;
                    whichNonLinear_[put] = i;
                }
            }
            delete [] pos;
            delete [] whichRows;
            analyzeObjects();
        }
    }
    // See if there are any quadratic bounds
    int nQ = 0;
    const CoinPackedMatrix * rowCopy = getMatrixByRow();
    //const double * element = rowCopy->getElements();
    //const int * column = rowCopy->getIndices();
    //const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
    const int * rowLength = rowCopy->getVectorLengths();
    const double * rowLower = getRowLower();
    const double * rowUpper = getRowUpper();
    for (int iObject = 0; iObject < numberObjects_; iObject++) {
        OsiBiLinear * obj = dynamic_cast<OsiBiLinear *> (object_[iObject]);
        if (obj) {
            int xyRow = obj->xyRow();
            if (rowLength[xyRow] == 4 && false) {
                // we have simple bound
                nQ++;
                double coefficient = obj->coefficient();
                double lo = rowLower[xyRow];
                double up = rowUpper[xyRow];
                if (coefficient != 1.0) {
                    printf("*** double check code here\n");
                    if (coefficient < 0.0) {
                        double temp = lo;
                        lo = - up;
                        up = - temp;
                        coefficient = - coefficient;
                    }
                    if (lo > -1.0e20)
                        lo /= coefficient;
                    if (up < 1.0e20)
                        up /= coefficient;
                    setRowLower(xyRow, lo);
                    setRowUpper(xyRow, up);
                    // we also need to change elements in matrix_
                }
                int type = 0;
                if (lo == up) {
                    // good news
                    type = 3;
                    coefficient = lo;
                } else if (lo < -1.0e20) {
                    assert (up < 1.0e20);
                    coefficient = up;
                    type = 1;
                    // can we make equality?
                } else if (up > 1.0e20) {
                    coefficient = lo;
                    type = 2;
                    // can we make equality?
                } else {
                    // we would need extra code
                    abort();
                }
                obj->setBoundType(type);
                obj->setCoefficient(coefficient);
                // can do better if integer?
                assert (!isInteger(obj->xColumn()));
                assert (!isInteger(obj->yColumn()));
            }
        }
    }
    delete [] which;
    if ((specialOptions2_&16) != 0)
        addTighterConstraints();
}
// Add reformulated bilinear constraints
void
OsiSolverLink::addTighterConstraints()
{
    // This is first attempt - for now get working on trimloss
    int numberW = 0;
    int * xW = new int[numberObjects_];
    int * yW = new int[numberObjects_];
    // Points to firstlambda
    int * wW = new int[numberObjects_];
    // Coefficient
    double * alphaW = new double[numberObjects_];
    // Objects
    OsiBiLinear ** objW = new OsiBiLinear * [numberObjects_];
    int numberColumns = getNumCols();
    int firstLambda = numberColumns;
    // set up list (better to rethink and do properly as column ordered)
    int * list = new int[numberColumns];
    memset(list, 0, numberColumns*sizeof(int));
    int i;
    for ( i = 0; i < numberObjects_; i++) {
        OsiBiLinear * obj = dynamic_cast<OsiBiLinear *> (object_[i]);
        if (obj) {
            //obj->setBranchingStrategy(4); // ***** temp
            objW[numberW] = obj;
            xW[numberW] = obj->xColumn();
            yW[numberW] = obj->yColumn();
            list[xW[numberW]] = 1;
            list[yW[numberW]] = 1;
            wW[numberW] = obj->firstLambda();
            firstLambda = CoinMin(firstLambda, obj->firstLambda());
            alphaW[numberW] = obj->coefficient();
            //assert (alphaW[numberW]==1.0); // fix when occurs
            numberW++;
        }
    }
    int nList = 0;
    for (i = 0; i < numberColumns; i++) {
        if (list[i])
            list[nList++] = i;
    }
    // set up mark array
    char * mark = new char [firstLambda*firstLambda];
    memset(mark, 0, firstLambda*firstLambda);
    for (i = 0; i < numberW; i++) {
        int x = xW[i];
        int y = yW[i];
        mark[x*firstLambda+y] = 1;
        mark[y*firstLambda+x] = 1;
    }
    int numberRows2 = originalRowCopy_->getNumRows();
    int * addColumn = new int [numberColumns];
    double * addElement = new double [numberColumns];
    int * addW = new int [numberColumns];
    assert (objectiveRow_ < 0); // fix when occurs
    for (int iRow = 0; iRow < numberRows2; iRow++) {
        for (int iList = 0; iList < nList; iList++) {
            int kColumn = list[iList];
#ifndef NDEBUG
            const double * columnLower = getColLower();
#endif
            //const double * columnUpper = getColUpper();
            const double * rowLower = getRowLower();
            const double * rowUpper = getRowUpper();
            const CoinPackedMatrix * rowCopy = getMatrixByRow();
            const double * element = rowCopy->getElements();
            const int * column = rowCopy->getIndices();
            const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
            const int * rowLength = rowCopy->getVectorLengths();
            CoinBigIndex j;
            int numberElements = rowLength[iRow];
            int n = 0;
            for (j = rowStart[iRow]; j < rowStart[iRow] + numberElements; j++) {
                int iColumn = column[j];
                if (iColumn >= firstLambda) {
                    // no good
                    n = -1;
                    break;
                }
                if (mark[iColumn*firstLambda+kColumn])
                    n++;
            }
            if (n == numberElements) {
                printf("can add row %d\n", iRow);
                assert (columnLower[kColumn] >= 0); // might be able to fix
                n = 0;
                for (j = rowStart[iRow]; j < rowStart[iRow] + numberElements; j++) {
                    int xColumn = kColumn;
                    int yColumn = column[j];
                    int k;
                    for (k = 0; k < numberW; k++) {
                        if ((xW[k] == yColumn && yW[k] == xColumn) ||
                                (yW[k] == yColumn && xW[k] == xColumn))
                            break;
                    }
                    assert (k < numberW);
                    if (xW[k] != xColumn) {
                        int temp = xColumn;
                        xColumn = yColumn;
                        yColumn = temp;
                    }
                    addW[n/4] = k;
                    int start = wW[k];
                    double value = element[j];
                    for (int kk = 0; kk < 4; kk++) {
                        // Dummy value
                        addElement[n] = value;
                        addColumn[n++] = start + kk;
                    }
                }
                addColumn[n++] = kColumn;
                double lo = rowLower[iRow];
                double up = rowUpper[iRow];
                if (lo > -1.0e20) {
                    // and tell object
                    for (j = 0; j < n - 1; j += 4) {
                        int iObject = addW[j/4];
                        objW[iObject]->addExtraRow(matrix_->getNumRows(), addElement[j]);
                    }
                    addElement[n-1] = -lo;
                    if (lo == up)
                        addRow(n, addColumn, addElement, 0.0, 0.0);
                    else
                        addRow(n, addColumn, addElement, 0.0, COIN_DBL_MAX);
                    matrix_->appendRow(n, addColumn, addElement);
                }
                if (up<1.0e20 && up>lo) {
                    // and tell object
                    for (j = 0; j < n - 1; j += 4) {
                        int iObject = addW[j/4];
                        objW[iObject]->addExtraRow(matrix_->getNumRows(), addElement[j]);
                    }
                    addElement[n-1] = -up;
                    addRow(n, addColumn, addElement, -COIN_DBL_MAX, 0.0);
                    matrix_->appendRow(n, addColumn, addElement);
                }
            }
        }
    }
#ifdef JJF_ZERO
    // possibly do bounds
    for (int iColumn = 0; iColumn < firstLambda; iColumn++) {
        for (int iList = 0; iList < nList; iList++) {
            int kColumn = list[iList];
            const double * columnLower = getColLower();
            const double * columnUpper = getColUpper();
            if (mark[iColumn*firstLambda+kColumn]) {
                printf("can add column %d\n", iColumn);
                assert (columnLower[kColumn] >= 0); // might be able to fix
                int xColumn = kColumn;
                int yColumn = iColumn;
                int k;
                for (k = 0; k < numberW; k++) {
                    if ((xW[k] == yColumn && yW[k] == xColumn) ||
                            (yW[k] == yColumn && xW[k] == xColumn))
                        break;
                }
                assert (k < numberW);
                if (xW[k] != xColumn) {
                    int temp = xColumn;
                    xColumn = yColumn;
                    yColumn = temp;
                }
                int start = wW[k];
                int n = 0;
                for (int kk = 0; kk < 4; kk++) {
                    // Dummy value
                    addElement[n] = 1.0e-19;
                    addColumn[n++] = start + kk;
                }
                // Tell object about this
                objW[k]->addExtraRow(matrix_->getNumRows(), 1.0);
                addColumn[n++] = kColumn;
                double lo = columnLower[iColumn];
                double up = columnUpper[iColumn];
                if (lo > -1.0e20) {
                    addElement[n-1] = -lo;
                    if (lo == up)
                        addRow(n, addColumn, addElement, 0.0, 0.0);
                    else
                        addRow(n, addColumn, addElement, 0.0, COIN_DBL_MAX);
                    matrix_->appendRow(n, addColumn, addElement);
                }
                if (up<1.0e20 && up>lo) {
                    addElement[n-1] = -up;
                    addRow(n, addColumn, addElement, -COIN_DBL_MAX, 0.0);
                    matrix_->appendRow(n, addColumn, addElement);
                }
            }
        }
    }
#endif
    delete [] xW;
    delete [] yW;
    delete [] wW;
    delete [] alphaW;
    delete [] addColumn;
    delete [] addElement;
    delete [] addW;
    delete [] mark;
    delete [] list;
    delete [] objW;
}
// Set all biLinear priorities on x-x variables
void
OsiSolverLink::setBiLinearPriorities(int value, double meshSize)
{
    OsiObject ** newObject = new OsiObject * [numberObjects_];
    int numberOdd = 0;
    int i;
    for ( i = 0; i < numberObjects_; i++) {
        OsiBiLinear * obj = dynamic_cast<OsiBiLinear *> (object_[i]);
        if (obj) {
            if (obj->xMeshSize() < 1.0 && obj->yMeshSize() < 1.0) {
                double oldSatisfied = CoinMax(obj->xSatisfied(),
                                              obj->ySatisfied());
                OsiBiLinear * objNew = new OsiBiLinear(*obj);
                newObject[numberOdd++] = objNew;
                objNew->setXSatisfied(0.5*meshSize);
                obj->setXOtherSatisfied(0.5*meshSize);
                objNew->setXOtherSatisfied(oldSatisfied);
                objNew->setXMeshSize(meshSize);
                objNew->setYSatisfied(0.5*meshSize);
                obj->setYOtherSatisfied(0.5*meshSize);
                objNew->setYOtherSatisfied(oldSatisfied);
                objNew->setYMeshSize(meshSize);
                objNew->setXYSatisfied(0.25*meshSize);
                objNew->setPriority(value);
                objNew->setBranchingStrategy(8);
            }
        }
    }
    addObjects(numberOdd, newObject);
    for (i = 0; i < numberOdd; i++)
        delete newObject[i];
    delete [] newObject;
}
/* Set options and priority on all or some biLinear variables
   1 - on I-I
   2 - on I-x
   4 - on x-x
      or combinations.
      -1 means leave (for priority value and strategy value)
*/
void
OsiSolverLink::setBranchingStrategyOnVariables(int strategyValue, int priorityValue,
        int mode)
{
    int i;
    for ( i = 0; i < numberObjects_; i++) {
        OsiBiLinear * obj = dynamic_cast<OsiBiLinear *> (object_[i]);
        if (obj) {
            bool change = false;
            if (obj->xMeshSize() < 1.0 && obj->yMeshSize() < 1.0 && (mode&4) != 0)
                change = true;
            else if (((obj->xMeshSize() == 1.0 && obj->yMeshSize() < 1.0) ||
                      (obj->xMeshSize() < 1.0 && obj->yMeshSize() == 1.0)) && (mode&2) != 0)
                change = true;
            else if (obj->xMeshSize() == 1.0 && obj->yMeshSize() == 1.0 && (mode&1) != 0)
                change = true;
            else if (obj->xMeshSize() > 1.0 || obj->yMeshSize() > 1.0)
                abort();
            if (change) {
                if (strategyValue >= 0)
                    obj->setBranchingStrategy(strategyValue);
                if (priorityValue >= 0)
                    obj->setPriority(priorityValue);
            }
        }
    }
}

// Say convex (should work it out)
void
OsiSolverLink::sayConvex(bool convex)
{
    specialOptions2_ |= 4;
    if (convex_) {
        for (int iNon = 0; iNon < numberNonLinearRows_; iNon++) {
            convex_[iNon] = convex ? 1 : -1;
        }
    }
}
// Set all mesh sizes on x-x variables
void
OsiSolverLink::setMeshSizes(double value)
{
    int i;
    for ( i = 0; i < numberObjects_; i++) {
        OsiBiLinear * obj = dynamic_cast<OsiBiLinear *> (object_[i]);
        if (obj) {
            if (obj->xMeshSize() < 1.0 && obj->yMeshSize() < 1.0) {
#ifdef JJF_ZERO
                numberContinuous++;
                int xColumn = obj->xColumn();
                double gapX = upper[xColumn] - lower[xColumn];
                int yColumn = obj->yColumn();
                double gapY = upper[yColumn] - lower[yColumn];
                gap = CoinMax(gap, CoinMax(gapX, gapY));
#endif
                obj->setMeshSizes(this, value, value);
            }
        }
    }
}
/* Solves nonlinear problem from CoinModel using SLP - may be used as crash
   for other algorithms when number of iterations small.
   Also exits if all problematical variables are changing
   less than deltaTolerance
   Returns solution array
*/
double *
OsiSolverLink::nonlinearSLP(int numberPasses, double deltaTolerance)
{
    if (!coinModel_.numberRows()) {
        printf("Model not set up or nonlinear arrays not created!\n");
        return NULL;
    }
    // first check and set up arrays
    int numberColumns = coinModel_.numberColumns();
    int numberRows = coinModel_.numberRows();
    char * markNonlinear = new char [numberColumns+numberRows];
    CoinZeroN(markNonlinear, numberColumns + numberRows);
    // List of nonlinear entries
    int * listNonLinearColumn = new int[numberColumns];
    // List of nonlinear constraints
    int * whichRow = new int [numberRows];
    CoinZeroN(whichRow, numberRows);
    int numberNonLinearColumns = 0;
    int iColumn;
    CoinModel coinModel = coinModel_;
    //const CoinModelHash * stringArray = coinModel.stringArray();
    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
        CoinModelLink triple = coinModel.firstInColumn(iColumn);
        bool linear = true;
        int n = 0;
        // See if nonlinear objective
        const char * expr = coinModel.getColumnObjectiveAsString(iColumn);
        if (strcmp(expr, "Numeric")) {
            linear = false;
            // try and see which columns
            assert (strlen(expr) < 20000);
            char temp[20000];
            strcpy(temp, expr);
            char * pos = temp;
            bool ifFirst = true;
            while (*pos) {
                double value;
                int jColumn = decodeBit(pos, pos, value, ifFirst, coinModel_);
                // must be column unless first when may be linear term
                if (jColumn >= 0) {
                    markNonlinear[jColumn] = 1;
                } else if (jColumn != -2) {
                    printf("bad nonlinear term %s\n", temp);
                    abort();
                }
                ifFirst = false;
            }
        }
        while (triple.row() >= 0) {
            int iRow = triple.row();
            const char * expr = coinModel.getElementAsString(iRow, iColumn);
            if (strcmp(expr, "Numeric")) {
                linear = false;
                whichRow[iRow]++;
                // try and see which columns
                assert (strlen(expr) < 20000);
                char temp[20000];
                strcpy(temp, expr);
                char * pos = temp;
                bool ifFirst = true;
                while (*pos) {
                    double value;
                    int jColumn = decodeBit(pos, pos, value, ifFirst, coinModel_);
                    // must be column unless first when may be linear term
                    if (jColumn >= 0) {
                        markNonlinear[jColumn] = 1;
                    } else if (jColumn != -2) {
                        printf("bad nonlinear term %s\n", temp);
                        abort();
                    }
                    ifFirst = false;
                }
            }
            triple = coinModel.next(triple);
            n++;
        }
        if (!linear) {
            markNonlinear[iColumn] = 1;
        }
    }
    //int xxxx[]={3,2,0,4,3,0};
    //double initialSolution[6];
    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
        if (markNonlinear[iColumn]) {
            // put in something
            double lower = coinModel.columnLower(iColumn);
            double upper = CoinMin(coinModel.columnUpper(iColumn), lower + 1000.0);
            coinModel.associateElement(coinModel.columnName(iColumn), 0.5*(lower + upper));
            //coinModel.associateElement(coinModel.columnName(iColumn),xxxx[iColumn]);
            listNonLinearColumn[numberNonLinearColumns++] = iColumn;
            //initialSolution[iColumn]=xxxx[iColumn];
        }
    }
    // if nothing just solve
    if (!numberNonLinearColumns) {
        delete [] listNonLinearColumn;
        delete [] whichRow;
        delete [] markNonlinear;
        ClpSimplex tempModel;
        tempModel.loadProblem(coinModel, true);
        tempModel.initialSolve();
        double * solution = CoinCopyOfArray(tempModel.getColSolution(), numberColumns);
        return solution;
    }
    // Create artificials
    ClpSimplex tempModel;
    tempModel.loadProblem(coinModel, true);
    const double * rowLower = tempModel.rowLower();
    const double * rowUpper = tempModel.rowUpper();
    bool takeAll = false;
    int iRow;
    int numberArtificials = 0;
    for (iRow = 0; iRow < numberRows; iRow++) {
        if (whichRow[iRow] || takeAll) {
            if (rowLower[iRow] > -1.0e30)
                numberArtificials++;
            if (rowUpper[iRow] < 1.0e30)
                numberArtificials++;
        }
    }
    CoinBigIndex * startArtificial = new CoinBigIndex [numberArtificials+1];
    int * rowArtificial = new int [numberArtificials];
    double * elementArtificial = new double [numberArtificials];
    double * objectiveArtificial = new double [numberArtificials];
    numberArtificials = 0;
    startArtificial[0] = 0;
    double artificialCost = 1.0e9;
    for (iRow = 0; iRow < numberRows; iRow++) {
        if (whichRow[iRow] || takeAll) {
            if (rowLower[iRow] > -1.0e30) {
                rowArtificial[numberArtificials] = iRow;
                elementArtificial[numberArtificials] = 1.0;
                objectiveArtificial[numberArtificials] = artificialCost;
                numberArtificials++;
                startArtificial[numberArtificials] = numberArtificials;
            }
            if (rowUpper[iRow] < 1.0e30) {
                rowArtificial[numberArtificials] = iRow;
                elementArtificial[numberArtificials] = -1.0;
                objectiveArtificial[numberArtificials] = artificialCost;
                numberArtificials++;
                startArtificial[numberArtificials] = numberArtificials;
            }
        }
    }
    // Get first solution
    int numberColumnsSmall = numberColumns;
    ClpSimplex model;
    model.loadProblem(coinModel, true);
    model.addColumns(numberArtificials, NULL, NULL, objectiveArtificial,
                     startArtificial, rowArtificial, elementArtificial);
    double * columnLower = model.columnLower();
    double * columnUpper = model.columnUpper();
    double * trueLower = new double[numberNonLinearColumns];
    double * trueUpper = new double[numberNonLinearColumns];
    int jNon;
    for (jNon = 0; jNon < numberNonLinearColumns; jNon++) {
        iColumn = listNonLinearColumn[jNon];
        trueLower[jNon] = columnLower[iColumn];
        trueUpper[jNon] = columnUpper[iColumn];
        //columnLower[iColumn]=initialSolution[iColumn];
        //columnUpper[iColumn]=initialSolution[iColumn];
    }
    model.initialSolve();
    //model.writeMps("bad.mps");
    // redo number of columns
    numberColumns = model.numberColumns();
    int * last[3];
    double * solution = model.primalColumnSolution();

    double * trust = new double[numberNonLinearColumns];
    for (jNon = 0; jNon < numberNonLinearColumns; jNon++) {
        iColumn = listNonLinearColumn[jNon];
        trust[jNon] = 0.5;
        if (solution[iColumn] < trueLower[jNon])
            solution[iColumn] = trueLower[jNon];
        else if (solution[iColumn] > trueUpper[jNon])
            solution[iColumn] = trueUpper[jNon];
    }
    int iPass;
    double lastObjective = 1.0e31;
    double * saveSolution = new double [numberColumns];
    double * saveRowSolution = new double [numberRows];
    memset(saveRowSolution, 0, numberRows*sizeof(double));
    double * savePi = new double [numberRows];
    double * safeSolution = new double [numberColumns];
    unsigned char * saveStatus = new unsigned char[numberRows+numberColumns];
    double targetDrop = 1.0e31;
    //double objectiveOffset;
    //model.getDblParam(ClpObjOffset,objectiveOffset);
    // 1 bound up, 2 up, -1 bound down, -2 down, 0 no change
    for (iPass = 0; iPass < 3; iPass++) {
        last[iPass] = new int[numberNonLinearColumns];
        for (jNon = 0; jNon < numberNonLinearColumns; jNon++)
            last[iPass][jNon] = 0;
    }
    // goodMove +1 yes, 0 no, -1 last was bad - just halve gaps, -2 do nothing
    int goodMove = -2;
    char * statusCheck = new char[numberColumns];
    double * changeRegion = new double [numberColumns];
    int logLevel = 63;
    double dualTolerance = model.dualTolerance();
    double primalTolerance = model.primalTolerance();
    int lastGoodMove = 1;
    for (iPass = 0; iPass < numberPasses; iPass++) {
        lastGoodMove = goodMove;
        columnLower = model.columnLower();
        columnUpper = model.columnUpper();
        solution = model.primalColumnSolution();
        double * rowActivity = model.primalRowSolution();
        // redo objective
        ClpSimplex tempModel;
        // load new values
        for (jNon = 0; jNon < numberNonLinearColumns; jNon++) {
            iColumn = listNonLinearColumn[jNon];
            coinModel.associateElement(coinModel.columnName(iColumn), solution[iColumn]);
        }
        tempModel.loadProblem(coinModel);
        double objectiveOffset;
        tempModel.getDblParam(ClpObjOffset, objectiveOffset);
        double objValue = -objectiveOffset;
        const double * objective = tempModel.objective();
        for (iColumn = 0; iColumn < numberColumnsSmall; iColumn++)
            objValue += solution[iColumn] * objective[iColumn];
        double * rowActivity2 = tempModel.primalRowSolution();
        const double * rowLower2 = tempModel.rowLower();
        const double * rowUpper2 = tempModel.rowUpper();
        memset(rowActivity2, 0, numberRows*sizeof(double));
        tempModel.times(1.0, solution, rowActivity2);
        for (iRow = 0; iRow < numberRows; iRow++) {
            if (rowActivity2[iRow] < rowLower2[iRow] - primalTolerance)
                objValue += (rowLower2[iRow] - rowActivity2[iRow] - primalTolerance) * artificialCost;
            else if (rowActivity2[iRow] > rowUpper2[iRow] + primalTolerance)
                objValue -= (rowUpper2[iRow] - rowActivity2[iRow] + primalTolerance) * artificialCost;
        }
        double theta = -1.0;
        double maxTheta = COIN_DBL_MAX;
        if (objValue <= lastObjective + 1.0e-15*fabs(lastObjective) || !iPass)
            goodMove = 1;
        else
            goodMove = -1;
        //maxTheta=1.0;
        if (iPass) {
            int jNon = 0;
            for (iColumn = 0; iColumn < numberColumns; iColumn++) {
                changeRegion[iColumn] = solution[iColumn] - saveSolution[iColumn];
                double alpha = changeRegion[iColumn];
                double oldValue = saveSolution[iColumn];
                if (markNonlinear[iColumn] == 0) {
                    // linear
                    if (alpha < -1.0e-15) {
                        // variable going towards lower bound
                        double bound = columnLower[iColumn];
                        oldValue -= bound;
                        if (oldValue + maxTheta*alpha < 0.0) {
                            maxTheta = CoinMax(0.0, oldValue / (-alpha));
                        }
                    } else if (alpha > 1.0e-15) {
                        // variable going towards upper bound
                        double bound = columnUpper[iColumn];
                        oldValue = bound - oldValue;
                        if (oldValue - maxTheta*alpha < 0.0) {
                            maxTheta = CoinMax(0.0, oldValue / alpha);
                        }
                    }
                } else {
                    // nonlinear
                    if (alpha < -1.0e-15) {
                        // variable going towards lower bound
                        double bound = trueLower[jNon];
                        oldValue -= bound;
                        if (oldValue + maxTheta*alpha < 0.0) {
                            maxTheta = CoinMax(0.0, oldValue / (-alpha));
                        }
                    } else if (alpha > 1.0e-15) {
                        // variable going towards upper bound
                        double bound = trueUpper[jNon];
                        oldValue = bound - oldValue;
                        if (oldValue - maxTheta*alpha < 0.0) {
                            maxTheta = CoinMax(0.0, oldValue / alpha);
                        }
                    }
                    jNon++;
                }
            }
            // make sure both accurate
            memset(rowActivity, 0, numberRows*sizeof(double));
            model.times(1.0, solution, rowActivity);
            memset(saveRowSolution, 0, numberRows*sizeof(double));
            model.times(1.0, saveSolution, saveRowSolution);
            for (int iRow = 0; iRow < numberRows; iRow++) {
                double alpha = rowActivity[iRow] - saveRowSolution[iRow];
                double oldValue = saveRowSolution[iRow];
                if (alpha < -1.0e-15) {
                    // variable going towards lower bound
                    double bound = rowLower[iRow];
                    oldValue -= bound;
                    if (oldValue + maxTheta*alpha < 0.0) {
                        maxTheta = CoinMax(0.0, oldValue / (-alpha));
                    }
                } else if (alpha > 1.0e-15) {
                    // variable going towards upper bound
                    double bound = rowUpper[iRow];
                    oldValue = bound - oldValue;
                    if (oldValue - maxTheta*alpha < 0.0) {
                        maxTheta = CoinMax(0.0, oldValue / alpha);
                    }
                }
            }
        } else {
            for (iColumn = 0; iColumn < numberColumns; iColumn++) {
                changeRegion[iColumn] = 0.0;
                saveSolution[iColumn] = solution[iColumn];
            }
            memcpy(saveRowSolution, rowActivity, numberRows*sizeof(double));
        }
        if (goodMove >= 0) {
            //theta = CoinMin(theta2,maxTheta);
            theta = maxTheta;
            if (theta > 0.0 && theta <= 1.0) {
                // update solution
                double lambda = 1.0 - theta;
                for (iColumn = 0; iColumn < numberColumns; iColumn++)
                    solution[iColumn] = lambda * saveSolution[iColumn]
                                        + theta * solution[iColumn];
                memset(rowActivity, 0, numberRows*sizeof(double));
                model.times(1.0, solution, rowActivity);
                if (lambda > 0.999) {
                    memcpy(model.dualRowSolution(), savePi, numberRows*sizeof(double));
                    memcpy(model.statusArray(), saveStatus, numberRows + numberColumns);
                }
                // redo rowActivity
                memset(rowActivity, 0, numberRows*sizeof(double));
                model.times(1.0, solution, rowActivity);
            }
        }
        // load new values
        for (jNon = 0; jNon < numberNonLinearColumns; jNon++) {
            iColumn = listNonLinearColumn[jNon];
            coinModel.associateElement(coinModel.columnName(iColumn), solution[iColumn]);
        }
        double * sol2 = CoinCopyOfArray(model.primalColumnSolution(), numberColumns);
        unsigned char * status2 = CoinCopyOfArray(model.statusArray(), numberColumns);
        model.loadProblem(coinModel);
        model.addColumns(numberArtificials, NULL, NULL, objectiveArtificial,
                         startArtificial, rowArtificial, elementArtificial);
        memcpy(model.primalColumnSolution(), sol2, numberColumns*sizeof(double));
        memcpy(model.statusArray(), status2, numberColumns);
        delete [] sol2;
        delete [] status2;
        columnLower = model.columnLower();
        columnUpper = model.columnUpper();
        solution = model.primalColumnSolution();
        rowActivity = model.primalRowSolution();
        int * temp = last[2];
        last[2] = last[1];
        last[1] = last[0];
        last[0] = temp;
        for (jNon = 0; jNon < numberNonLinearColumns; jNon++) {
            iColumn = listNonLinearColumn[jNon];
            double change = solution[iColumn] - saveSolution[iColumn];
            if (change < -1.0e-5) {
                if (fabs(change + trust[jNon]) < 1.0e-5)
                    temp[jNon] = -1;
                else
                    temp[jNon] = -2;
            } else if (change > 1.0e-5) {
                if (fabs(change - trust[jNon]) < 1.0e-5)
                    temp[jNon] = 1;
                else
                    temp[jNon] = 2;
            } else {
                temp[jNon] = 0;
            }
        }
        // goodMove +1 yes, 0 no, -1 last was bad - just halve gaps, -2 do nothing
        double maxDelta = 0.0;
        if (goodMove >= 0) {
            if (objValue <= lastObjective + 1.0e-15*fabs(lastObjective))
                goodMove = 1;
            else
                goodMove = 0;
        } else {
            maxDelta = 1.0e10;
        }
        double maxGap = 0.0;
        int numberSmaller = 0;
        int numberSmaller2 = 0;
        int numberLarger = 0;
        for (jNon = 0; jNon < numberNonLinearColumns; jNon++) {
            iColumn = listNonLinearColumn[jNon];
            maxDelta = CoinMax(maxDelta,
                               fabs(solution[iColumn] - saveSolution[iColumn]));
            if (goodMove > 0) {
                if (last[0][jNon]*last[1][jNon] < 0) {
                    // halve
                    trust[jNon] *= 0.5;
                    numberSmaller2++;
                } else {
                    if (last[0][jNon] == last[1][jNon] &&
                            last[0][jNon] == last[2][jNon])
                        trust[jNon] = CoinMin(1.5 * trust[jNon], 1.0e6);
                    numberLarger++;
                }
            } else if (goodMove != -2 && trust[jNon] > 10.0*deltaTolerance) {
                trust[jNon] *= 0.2;
                numberSmaller++;
            }
            maxGap = CoinMax(maxGap, trust[jNon]);
        }
#ifdef CLP_DEBUG
        if (logLevel&32)
            std::cout << "largest gap is " << maxGap << " "
                      << numberSmaller + numberSmaller2 << " reduced ("
                      << numberSmaller << " badMove ), "
                      << numberLarger << " increased" << std::endl;
#endif
        if (iPass > 10000) {
            for (jNon = 0; jNon < numberNonLinearColumns; jNon++)
                trust[jNon] *= 0.0001;
        }
        printf("last good %d goodMove %d\n", lastGoodMove, goodMove);
        if (goodMove > 0) {
            double drop = lastObjective - objValue;
            printf("Pass %d, objective %g - drop %g maxDelta %g\n", iPass, objValue, drop, maxDelta);
            if (iPass > 20 && drop < 1.0e-12*fabs(objValue) && lastGoodMove > 0)
                drop = 0.999e-4; // so will exit
            if (maxDelta < deltaTolerance && drop < 1.0e-4 && goodMove && theta<0.99999 && lastGoodMove>0) {
                if (logLevel > 1)
                    std::cout << "Exiting as maxDelta < tolerance and small drop" << std::endl;
                break;
            }
        } else if (!numberSmaller && iPass > 1) {
            if (logLevel > 1)
                std::cout << "Exiting as all gaps small" << std::endl;
            break;
        }
        if (!iPass)
            goodMove = 1;
        targetDrop = 0.0;
        double * r = model.dualColumnSolution();
        for (jNon = 0; jNon < numberNonLinearColumns; jNon++) {
            iColumn = listNonLinearColumn[jNon];
            columnLower[iColumn] = CoinMax(solution[iColumn]
                                           - trust[jNon],
                                           trueLower[jNon]);
            columnUpper[iColumn] = CoinMin(solution[iColumn]
                                           + trust[jNon],
                                           trueUpper[jNon]);
        }
        if (iPass) {
            // get reduced costs
            model.matrix()->transposeTimes(savePi,
                                           model.dualColumnSolution());
            const double * objective = model.objective();
            for (jNon = 0; jNon < numberNonLinearColumns; jNon++) {
                iColumn = listNonLinearColumn[jNon];
                double dj = objective[iColumn] - r[iColumn];
                r[iColumn] = dj;
                if (dj < -dualTolerance)
                    targetDrop -= dj * (columnUpper[iColumn] - solution[iColumn]);
                else if (dj > dualTolerance)
                    targetDrop -= dj * (columnLower[iColumn] - solution[iColumn]);
            }
        } else {
            memset(r, 0, numberColumns*sizeof(double));
        }
#ifdef JJF_ZERO
        for (jNon = 0; jNon < numberNonLinearColumns; jNon++) {
            iColumn = listNonLinearColumn[jNon];
            if (statusCheck[iColumn] == 'L' && r[iColumn] < -1.0e-4) {
                columnLower[iColumn] = CoinMax(solution[iColumn],
                                               trueLower[jNon]);
                columnUpper[iColumn] = CoinMin(solution[iColumn]
                                               + trust[jNon],
                                               trueUpper[jNon]);
            } else if (statusCheck[iColumn] == 'U' && r[iColumn] > 1.0e-4) {
                columnLower[iColumn] = CoinMax(solution[iColumn]
                                               - trust[jNon],
                                               trueLower[jNon]);
                columnUpper[iColumn] = CoinMin(solution[iColumn],
                                               trueUpper[jNon]);
            } else {
                columnLower[iColumn] = CoinMax(solution[iColumn]
                                               - trust[jNon],
                                               trueLower[jNon]);
                columnUpper[iColumn] = CoinMin(solution[iColumn]
                                               + trust[jNon],
                                               trueUpper[jNon]);
            }
        }
#endif
        if (goodMove > 0) {
            memcpy(saveSolution, solution, numberColumns*sizeof(double));
            memcpy(saveRowSolution, rowActivity, numberRows*sizeof(double));
            memcpy(savePi, model.dualRowSolution(), numberRows*sizeof(double));
            memcpy(saveStatus, model.statusArray(), numberRows + numberColumns);

#ifdef CLP_DEBUG
            if (logLevel&32)
                std::cout << "Pass - " << iPass
                          << ", target drop is " << targetDrop
                          << std::endl;
#endif
            lastObjective = objValue;
            if (targetDrop < CoinMax(1.0e-8, CoinMin(1.0e-6, 1.0e-6*fabs(objValue))) && lastGoodMove && iPass > 3) {
                if (logLevel > 1)
                    printf("Exiting on target drop %g\n", targetDrop);
                break;
            }
#ifdef CLP_DEBUG
            {
                double * r = model.dualColumnSolution();
                for (jNon = 0; jNon < numberNonLinearColumns; jNon++) {
                    iColumn = listNonLinearColumn[jNon];
                    if (logLevel&32)
                        printf("Trust %d %g - solution %d %g obj %g dj %g state %c - bounds %g %g\n",
                               jNon, trust[jNon], iColumn, solution[iColumn], objective[iColumn],
                               r[iColumn], statusCheck[iColumn], columnLower[iColumn],
                               columnUpper[iColumn]);
                }
            }
#endif
            model.scaling(false);
            model.primal(1);
            for (jNon = 0; jNon < numberNonLinearColumns; jNon++) {
                iColumn = listNonLinearColumn[jNon];
                printf("%d bounds etc %g %g %g\n", iColumn, columnLower[iColumn], solution[iColumn], columnUpper[iColumn]);
            }
            char temp[20];
            sprintf(temp, "pass%d.mps", iPass);
            //model.writeMps(temp);
#ifdef CLP_DEBUG
            if (model.status()) {
                model.writeMps("xx.mps");
            }
#endif
            if (model.status() == 1) {
                // not feasible ! - backtrack and exit
                // use safe solution
                memcpy(solution, safeSolution, numberColumns*sizeof(double));
                memcpy(saveSolution, solution, numberColumns*sizeof(double));
                memset(rowActivity, 0, numberRows*sizeof(double));
                model.times(1.0, solution, rowActivity);
                memcpy(saveRowSolution, rowActivity, numberRows*sizeof(double));
                memcpy(model.dualRowSolution(), savePi, numberRows*sizeof(double));
                memcpy(model.statusArray(), saveStatus, numberRows + numberColumns);
                for (jNon = 0; jNon < numberNonLinearColumns; jNon++) {
                    iColumn = listNonLinearColumn[jNon];
                    columnLower[iColumn] = CoinMax(solution[iColumn]
                                                   - trust[jNon],
                                                   trueLower[jNon]);
                    columnUpper[iColumn] = CoinMin(solution[iColumn]
                                                   + trust[jNon],
                                                   trueUpper[jNon]);
                }
                break;
            } else {
                // save in case problems
                memcpy(safeSolution, solution, numberColumns*sizeof(double));
            }
            goodMove = 1;
        } else {
            // bad pass - restore solution
#ifdef CLP_DEBUG
            if (logLevel&32)
                printf("Backtracking\n");
#endif
            // load old values
            for (jNon = 0; jNon < numberNonLinearColumns; jNon++) {
                iColumn = listNonLinearColumn[jNon];
                coinModel.associateElement(coinModel.columnName(iColumn), saveSolution[iColumn]);
            }
            model.loadProblem(coinModel);
            model.addColumns(numberArtificials, NULL, NULL, objectiveArtificial,
                             startArtificial, rowArtificial, elementArtificial);
            solution = model.primalColumnSolution();
            rowActivity = model.primalRowSolution();
            memcpy(solution, saveSolution, numberColumns*sizeof(double));
            memcpy(rowActivity, saveRowSolution, numberRows*sizeof(double));
            memcpy(model.dualRowSolution(), savePi, numberRows*sizeof(double));
            memcpy(model.statusArray(), saveStatus, numberRows + numberColumns);
            columnLower = model.columnLower();
            columnUpper = model.columnUpper();
            for (jNon = 0; jNon < numberNonLinearColumns; jNon++) {
                iColumn = listNonLinearColumn[jNon];
                columnLower[iColumn] = solution[iColumn];
                columnUpper[iColumn] = solution[iColumn];
            }
            model.primal(1);
            //model.writeMps("xx.mps");
            iPass--;
            goodMove = -1;
        }
    }
    // restore solution
    memcpy(solution, saveSolution, numberColumns*sizeof(double));
    delete [] statusCheck;
    delete [] savePi;
    delete [] saveStatus;
    // load new values
    for (jNon = 0; jNon < numberNonLinearColumns; jNon++) {
        iColumn = listNonLinearColumn[jNon];
        coinModel.associateElement(coinModel.columnName(iColumn), solution[iColumn]);
    }
    double * sol2 = CoinCopyOfArray(model.primalColumnSolution(), numberColumns);
    unsigned char * status2 = CoinCopyOfArray(model.statusArray(), numberColumns);
    model.loadProblem(coinModel);
    model.addColumns(numberArtificials, NULL, NULL, objectiveArtificial,
                     startArtificial, rowArtificial, elementArtificial);
    memcpy(model.primalColumnSolution(), sol2, numberColumns*sizeof(double));
    memcpy(model.statusArray(), status2, numberColumns);
    delete [] sol2;
    delete [] status2;
    columnLower = model.columnLower();
    columnUpper = model.columnUpper();
    solution = model.primalColumnSolution();
    for (jNon = 0; jNon < numberNonLinearColumns; jNon++) {
        iColumn = listNonLinearColumn[jNon];
        columnLower[iColumn] = CoinMax(solution[iColumn],
                                       trueLower[jNon]);
        columnUpper[iColumn] = CoinMin(solution[iColumn],
                                       trueUpper[jNon]);
    }
    model.primal(1);
    for (jNon = 0; jNon < numberNonLinearColumns; jNon++) {
        iColumn = listNonLinearColumn[jNon];
        columnLower[iColumn] = trueLower[jNon];
        columnUpper[iColumn] = trueUpper[jNon];
    }
    delete [] saveSolution;
    delete [] safeSolution;
    delete [] saveRowSolution;
    for (iPass = 0; iPass < 3; iPass++)
        delete [] last[iPass];
    delete [] trust;
    delete [] trueUpper;
    delete [] trueLower;
    delete [] changeRegion;
    delete [] startArtificial;
    delete [] rowArtificial;
    delete [] elementArtificial;
    delete [] objectiveArtificial;
    delete [] listNonLinearColumn;
    delete [] whichRow;
    delete [] markNonlinear;
    return CoinCopyOfArray(solution, coinModel.numberColumns());
}
/* Solve linearized quadratic objective branch and bound.
   Return cutoff and OA cut
*/
double
OsiSolverLink::linearizedBAB(CglStored * cut)
{
    double bestObjectiveValue = COIN_DBL_MAX;
    if (quadraticModel_) {
        ClpSimplex * qp = new ClpSimplex(*quadraticModel_);
        // bounds
        int numberColumns = qp->numberColumns();
        double * lower = qp->columnLower();
        double * upper = qp->columnUpper();
        const double * lower2 = getColLower();
        const double * upper2 = getColUpper();
        for (int i = 0; i < numberColumns; i++) {
            lower[i] = CoinMax(lower[i], lower2[i]);
            upper[i] = CoinMin(upper[i], upper2[i]);
        }
        qp->nonlinearSLP(20, 1.0e-5);
        qp->primal();
        OsiSolverLinearizedQuadratic solver2(qp);
        const double * solution = NULL;
        // Reduce printout
        solver2.setHintParam(OsiDoReducePrint, true, OsiHintTry);
        CbcModel model2(solver2);
        // Now do requested saves and modifications
        CbcModel * cbcModel = & model2;
        OsiSolverInterface * osiModel = model2.solver();
        OsiClpSolverInterface * osiclpModel = dynamic_cast< OsiClpSolverInterface*> (osiModel);
        ClpSimplex * clpModel = osiclpModel->getModelPtr();

        // Set changed values

        CglProbing probing;
        probing.setMaxProbe(10);
        probing.setMaxLook(10);
        probing.setMaxElements(200);
        probing.setMaxProbeRoot(50);
        probing.setMaxLookRoot(10);
        probing.setRowCuts(3);
        probing.setUsingObjective(true);
        cbcModel->addCutGenerator(&probing, -1, "Probing", true, false, false, -100, -1, -1);
        cbcModel->cutGenerator(0)->setTiming(true);

        CglGomory gomory;
        gomory.setLimitAtRoot(512);
        cbcModel->addCutGenerator(&gomory, -98, "Gomory", true, false, false, -100, -1, -1);
        cbcModel->cutGenerator(1)->setTiming(true);

        CglKnapsackCover knapsackCover;
        cbcModel->addCutGenerator(&knapsackCover, -98, "KnapsackCover", true, false, false, -100, -1, -1);
        cbcModel->cutGenerator(2)->setTiming(true);

        CglClique clique;
        clique.setStarCliqueReport(false);
        clique.setRowCliqueReport(false);
        clique.setMinViolation(0.1);
        cbcModel->addCutGenerator(&clique, -98, "Clique", true, false, false, -100, -1, -1);
        cbcModel->cutGenerator(3)->setTiming(true);

        CglMixedIntegerRounding2 mixedIntegerRounding2;
        cbcModel->addCutGenerator(&mixedIntegerRounding2, -98, "MixedIntegerRounding2", true, false, false, -100, -1, -1);
        cbcModel->cutGenerator(4)->setTiming(true);

        CglFlowCover flowCover;
        cbcModel->addCutGenerator(&flowCover, -98, "FlowCover", true, false, false, -100, -1, -1);
        cbcModel->cutGenerator(5)->setTiming(true);

        CglTwomir twomir;
        twomir.setMaxElements(250);
        cbcModel->addCutGenerator(&twomir, -99, "Twomir", true, false, false, -100, -1, -1);
        cbcModel->cutGenerator(6)->setTiming(true);
        // For now - switch off most heuristics (because CglPreProcess is bad with QP)
#ifndef JJF_ONE
        CbcHeuristicFPump heuristicFPump(*cbcModel);
        heuristicFPump.setWhen(13);
        heuristicFPump.setMaximumPasses(20);
        heuristicFPump.setMaximumRetries(7);
        heuristicFPump.setAbsoluteIncrement(4332.64);
        cbcModel->addHeuristic(&heuristicFPump);
        heuristicFPump.setInitialWeight(1);

        CbcHeuristicLocal heuristicLocal(*cbcModel);
        heuristicLocal.setSearchType(1);
        cbcModel->addHeuristic(&heuristicLocal);

        CbcHeuristicGreedyCover heuristicGreedyCover(*cbcModel);
        cbcModel->addHeuristic(&heuristicGreedyCover);

        CbcHeuristicGreedyEquality heuristicGreedyEquality(*cbcModel);
        cbcModel->addHeuristic(&heuristicGreedyEquality);
#endif

        CbcRounding rounding(*cbcModel);
        rounding.setHeuristicName("rounding");
        cbcModel->addHeuristic(&rounding);

        cbcModel->setNumberBeforeTrust(5);
        cbcModel->setSpecialOptions(2);
        cbcModel->messageHandler()->setLogLevel(1);
        cbcModel->setMaximumCutPassesAtRoot(-100);
        cbcModel->setMaximumCutPasses(1);
        cbcModel->setMinimumDrop(0.05);
        // For branchAndBound this may help
        clpModel->defaultFactorizationFrequency();
        clpModel->setDualBound(1.0001e+08);
        clpModel->setPerturbation(50);
        osiclpModel->setSpecialOptions(193);
        osiclpModel->messageHandler()->setLogLevel(0);
        osiclpModel->setIntParam(OsiMaxNumIterationHotStart, 100);
        osiclpModel->setHintParam(OsiDoReducePrint, true, OsiHintTry);
        // You can save some time by switching off message building
        // clpModel->messagesPointer()->setDetailMessages(100,10000,(int *) NULL);

        // Solve

        cbcModel->initialSolve();
        if (clpModel->tightenPrimalBounds() != 0) {
            std::cout << "Problem is infeasible - tightenPrimalBounds!" << std::endl;
            delete qp;
            return COIN_DBL_MAX;
        }
        clpModel->dual();  // clean up
        cbcModel->initialSolve();
        cbcModel->branchAndBound();
        OsiSolverLinearizedQuadratic * solver3 = dynamic_cast<OsiSolverLinearizedQuadratic *> (model2.solver());
        assert (solver3);
        solution = solver3->bestSolution();
        bestObjectiveValue = solver3->bestObjectiveValue();
        setBestObjectiveValue(bestObjectiveValue);
        setBestSolution(solution, solver3->getNumCols());
        // if convex
        if ((specialOptions2()&4) != 0) {
            if (cbcModel_)
                cbcModel_->lockThread();
            // add OA cut
            double offset;
            double * gradient = new double [numberColumns+1];
            memcpy(gradient, qp->objectiveAsObject()->gradient(qp, solution, offset, true, 2),
                   numberColumns*sizeof(double));
            double rhs = 0.0;
            int * column = new int[numberColumns+1];
            int n = 0;
            for (int i = 0; i < numberColumns; i++) {
                double value = gradient[i];
                if (fabs(value) > 1.0e-12) {
                    gradient[n] = value;
                    rhs += value * solution[i];
                    column[n++] = i;
                }
            }
            gradient[n] = -1.0;
            column[n++] = numberColumns;
            cut->addCut(-COIN_DBL_MAX, offset + 1.0e-7, n, column, gradient);
            delete [] gradient;
            delete [] column;
            if (cbcModel_)
                cbcModel_->unlockThread();
        }
        delete qp;
        printf("obj %g\n", bestObjectiveValue);
    }
    return bestObjectiveValue;
}
/* Solves nonlinear problem from CoinModel using SLP - and then tries to get
   heuristic solution
   Returns solution array
*/
double *
OsiSolverLink::heuristicSolution(int numberPasses, double deltaTolerance, int mode)
{
    // get a solution
    CoinModel tempModel = coinModel_;
    ClpSimplex * temp = approximateSolution(tempModel, numberPasses, deltaTolerance);
    int numberColumns = coinModel_.numberColumns();
    double * solution = CoinCopyOfArray(temp->primalColumnSolution(), numberColumns);
    delete temp;
    if (mode == 0) {
        return solution;
    } else if (mode == 2) {
        const double * lower = getColLower();
        const double * upper = getColUpper();
        for (int iObject = 0; iObject < numberObjects_; iObject++) {
            OsiSimpleInteger * obj = dynamic_cast<OsiSimpleInteger *> (object_[iObject]);
            if (obj && (obj->priority() < biLinearPriority_ || biLinearPriority_ <= 0)) {
                int iColumn = obj->columnNumber();
                double value = solution[iColumn];
                value = floor(value + 0.5);
                if (fabs(value - solution[iColumn]) > 0.01) {
                    setColLower(iColumn, CoinMax(lower[iColumn], value - CoinMax(defaultBound_, 0.0)));
                    setColUpper(iColumn, CoinMin(upper[iColumn], value + CoinMax(defaultBound_, 1.0)));
                } else {
                    // could fix to integer
                    setColLower(iColumn, CoinMax(lower[iColumn], value - CoinMax(defaultBound_, 0.0)));
                    setColUpper(iColumn, CoinMin(upper[iColumn], value + CoinMax(defaultBound_, 0.0)));
                }
            }
        }
        return solution;
    }
    OsiClpSolverInterface newSolver;
    if (mode == 1) {
        // round all with priority < biLinearPriority_
        setFixedPriority(biLinearPriority_);
        // ? should we save and restore coin model
        tempModel = coinModel_;
        // solve modified problem
        char * mark = new char[numberColumns];
        memset(mark, 0, numberColumns);
        for (int iObject = 0; iObject < numberObjects_; iObject++) {
            OsiSimpleInteger * obj = dynamic_cast<OsiSimpleInteger *> (object_[iObject]);
            if (obj && obj->priority() < biLinearPriority_) {
                int iColumn = obj->columnNumber();
                double value = solution[iColumn];
                value = ceil(value - 1.0e-7);
                tempModel.associateElement(coinModel_.columnName(iColumn), value);
                mark[iColumn] = 1;
            }
            OsiBiLinear * objB = dynamic_cast<OsiBiLinear *> (object_[iObject]);
            if (objB) {
                // if one or both continuous then fix one
                if (objB->xMeshSize() < 1.0) {
                    int xColumn = objB->xColumn();
                    double value = solution[xColumn];
                    tempModel.associateElement(coinModel_.columnName(xColumn), value);
                    mark[xColumn] = 1;
                } else if (objB->yMeshSize() < 1.0) {
                    int yColumn = objB->yColumn();
                    double value = solution[yColumn];
                    tempModel.associateElement(coinModel_.columnName(yColumn), value);
                    mark[yColumn] = 1;
                }
            }
        }
        CoinModel * reOrdered = tempModel.reorder(mark);
        assert (reOrdered);
        tempModel = *reOrdered;
        delete reOrdered;
        delete [] mark;
        newSolver.loadFromCoinModel(tempModel, true);
        for (int iObject = 0; iObject < numberObjects_; iObject++) {
            OsiSimpleInteger * obj = dynamic_cast<OsiSimpleInteger *> (object_[iObject]);
            if (obj && obj->priority() < biLinearPriority_) {
                int iColumn = obj->columnNumber();
                double value = solution[iColumn];
                value = ceil(value - 1.0e-7);
                newSolver.setColLower(iColumn, value);
                newSolver.setColUpper(iColumn, value);
            }
            OsiBiLinear * objB = dynamic_cast<OsiBiLinear *> (object_[iObject]);
            if (objB) {
                // if one or both continuous then fix one
                if (objB->xMeshSize() < 1.0) {
                    int xColumn = objB->xColumn();
                    double value = solution[xColumn];
                    newSolver.setColLower(xColumn, value);
                    newSolver.setColUpper(xColumn, value);
                } else if (objB->yMeshSize() < 1.0) {
                    int yColumn = objB->yColumn();
                    double value = solution[yColumn];
                    newSolver.setColLower(yColumn, value);
                    newSolver.setColUpper(yColumn, value);
                }
            }
        }
    }
    CbcModel model(newSolver);
    // Now do requested saves and modifications
    CbcModel * cbcModel = & model;
    OsiSolverInterface * osiModel = model.solver();
    OsiClpSolverInterface * osiclpModel = dynamic_cast< OsiClpSolverInterface*> (osiModel);
    ClpSimplex * clpModel = osiclpModel->getModelPtr();
    CglProbing probing;
    probing.setMaxProbe(10);
    probing.setMaxLook(10);
    probing.setMaxElements(200);
    probing.setMaxProbeRoot(50);
    probing.setMaxLookRoot(10);
    probing.setRowCuts(3);
    probing.setRowCuts(0);
    probing.setUsingObjective(true);
    cbcModel->addCutGenerator(&probing, -1, "Probing", true, false, false, -100, -1, -1);

    CglGomory gomory;
    gomory.setLimitAtRoot(512);
    cbcModel->addCutGenerator(&gomory, -98, "Gomory", true, false, false, -100, -1, -1);

    CglKnapsackCover knapsackCover;
    cbcModel->addCutGenerator(&knapsackCover, -98, "KnapsackCover", true, false, false, -100, -1, -1);

    CglClique clique;
    clique.setStarCliqueReport(false);
    clique.setRowCliqueReport(false);
    clique.setMinViolation(0.1);
    cbcModel->addCutGenerator(&clique, -98, "Clique", true, false, false, -100, -1, -1);
    CglMixedIntegerRounding2 mixedIntegerRounding2;
    cbcModel->addCutGenerator(&mixedIntegerRounding2, -98, "MixedIntegerRounding2", true, false, false, -100, -1, -1);

    CglFlowCover flowCover;
    cbcModel->addCutGenerator(&flowCover, -98, "FlowCover", true, false, false, -100, -1, -1);

    CglTwomir twomir;
    twomir.setMaxElements(250);
    cbcModel->addCutGenerator(&twomir, -99, "Twomir", true, false, false, -100, -1, -1);
    cbcModel->cutGenerator(6)->setTiming(true);

    CbcHeuristicFPump heuristicFPump(*cbcModel);
    heuristicFPump.setWhen(1);
    heuristicFPump.setMaximumPasses(20);
    heuristicFPump.setDefaultRounding(0.5);
    cbcModel->addHeuristic(&heuristicFPump);

    CbcRounding rounding(*cbcModel);
    cbcModel->addHeuristic(&rounding);

    CbcHeuristicLocal heuristicLocal(*cbcModel);
    heuristicLocal.setSearchType(1);
    cbcModel->addHeuristic(&heuristicLocal);

    CbcHeuristicGreedyCover heuristicGreedyCover(*cbcModel);
    cbcModel->addHeuristic(&heuristicGreedyCover);

    CbcHeuristicGreedyEquality heuristicGreedyEquality(*cbcModel);
    cbcModel->addHeuristic(&heuristicGreedyEquality);

    CbcCompareDefault compare;
    cbcModel->setNodeComparison(compare);
    cbcModel->setNumberBeforeTrust(5);
    cbcModel->setSpecialOptions(2);
    cbcModel->messageHandler()->setLogLevel(1);
    cbcModel->setMaximumCutPassesAtRoot(-100);
    cbcModel->setMaximumCutPasses(1);
    cbcModel->setMinimumDrop(0.05);
    clpModel->setNumberIterations(1);
    // For branchAndBound this may help
    clpModel->defaultFactorizationFrequency();
    clpModel->setDualBound(6.71523e+07);
    clpModel->setPerturbation(50);
    osiclpModel->setSpecialOptions(193);
    osiclpModel->messageHandler()->setLogLevel(0);
    osiclpModel->setIntParam(OsiMaxNumIterationHotStart, 100);
    osiclpModel->setHintParam(OsiDoReducePrint, true, OsiHintTry);
    // You can save some time by switching off message building
    // clpModel->messagesPointer()->setDetailMessages(100,10000,(int *) NULL);
    // Solve

    cbcModel->initialSolve();
    //double cutoff = model_->getCutoff();
    if (!cbcModel_)
        cbcModel->setCutoff(1.0e50);
    else
        cbcModel->setCutoff(cbcModel_->getCutoff());
    int saveLogLevel = clpModel->logLevel();
    clpModel->setLogLevel(0);
#ifndef NDEBUG
    int returnCode = 0;
#endif
    if (clpModel->tightenPrimalBounds() != 0) {
        clpModel->setLogLevel(saveLogLevel);
#ifndef NDEBUG
        returnCode = -1; // infeasible//std::cout<<"Problem is infeasible - tightenPrimalBounds!"<<std::endl;
#endif
        //clpModel->writeMps("infeas2.mps");
    } else {
        clpModel->setLogLevel(saveLogLevel);
        clpModel->dual();  // clean up
        // compute some things using problem size
        cbcModel->setMinimumDrop(CoinMin(5.0e-2,
                                         fabs(cbcModel->getMinimizationObjValue())*1.0e-3 + 1.0e-4));
        if (cbcModel->getNumCols() < 500)
            cbcModel->setMaximumCutPassesAtRoot(-100); // always do 100 if possible
        else if (cbcModel->getNumCols() < 5000)
            cbcModel->setMaximumCutPassesAtRoot(100); // use minimum drop
        else
            cbcModel->setMaximumCutPassesAtRoot(20);
        cbcModel->setMaximumCutPasses(1);
        // Hand coded preprocessing
        CglPreProcess process;
        OsiSolverInterface * saveSolver = cbcModel->solver()->clone();
        // Tell solver we are in Branch and Cut
        saveSolver->setHintParam(OsiDoInBranchAndCut, true, OsiHintDo) ;
        // Default set of cut generators
        CglProbing generator1;
        generator1.setUsingObjective(true);
        generator1.setMaxPass(3);
        generator1.setMaxProbeRoot(saveSolver->getNumCols());
        generator1.setMaxElements(100);
        generator1.setMaxLookRoot(50);
        generator1.setRowCuts(3);
        // Add in generators
        process.addCutGenerator(&generator1);
        process.messageHandler()->setLogLevel(cbcModel->logLevel());
        OsiSolverInterface * solver2 =
            process.preProcessNonDefault(*saveSolver, 0, 10);
        // Tell solver we are not in Branch and Cut
        saveSolver->setHintParam(OsiDoInBranchAndCut, false, OsiHintDo) ;
        if (solver2)
            solver2->setHintParam(OsiDoInBranchAndCut, false, OsiHintDo) ;
        if (!solver2) {
            std::cout << "Pre-processing says infeasible!" << std::endl;
            delete saveSolver;
#ifndef NDEBUG
            returnCode = -1;
#endif
        } else {
            std::cout << "processed model has " << solver2->getNumRows()
                      << " rows, " << solver2->getNumCols()
                      << " and " << solver2->getNumElements() << std::endl;
            // we have to keep solver2 so pass clone
            solver2 = solver2->clone();
            //solver2->writeMps("intmodel");
            cbcModel->assignSolver(solver2);
            cbcModel->initialSolve();
            cbcModel->branchAndBound();
            // For best solution
            int numberColumns = newSolver.getNumCols();
            if (cbcModel->getMinimizationObjValue() < 1.0e50) {
                // post process
                process.postProcess(*cbcModel->solver());
                // Solution now back in saveSolver
                cbcModel->assignSolver(saveSolver);
                memcpy(cbcModel->bestSolution(), cbcModel->solver()->getColSolution(),
                       numberColumns*sizeof(double));
                // put back in original solver
                newSolver.setColSolution(cbcModel->bestSolution());
            } else {
                delete saveSolver;
            }
        }
    }
    assert (!returnCode);
    abort();
    return solution;
}
// Analyze constraints to see which are convex (quadratic)
void
OsiSolverLink::analyzeObjects()
{
    // space for starts
    int numberColumns = coinModel_.numberColumns();
    int * start = new int [numberColumns+1];
    const double * rowLower = getRowLower();
    const double * rowUpper = getRowUpper();
    for (int iNon = 0; iNon < numberNonLinearRows_; iNon++) {
        int iRow = rowNonLinear_[iNon];
        int numberElements = startNonLinear_[iNon+1] - startNonLinear_[iNon];
        // triplet arrays
        int * iColumn = new int [2*numberElements+1];
        int * jColumn = new int [2*numberElements];
        double * element = new double [2*numberElements];
        int i;
        int n = 0;
        for ( i = startNonLinear_[iNon]; i < startNonLinear_[iNon+1]; i++) {
            OsiBiLinear * obj = dynamic_cast<OsiBiLinear *> (object_[whichNonLinear_[i]]);
            assert (obj);
            int xColumn = obj->xColumn();
            int yColumn = obj->yColumn();
            double coefficient = obj->coefficient();
            if (xColumn != yColumn) {
                iColumn[n] = xColumn;
                jColumn[n] = yColumn;
                element[n++] = coefficient;
                iColumn[n] = yColumn;
                jColumn[n] = xColumn;
                element[n++] = coefficient;
            } else {
                iColumn[n] = xColumn;
                jColumn[n] = xColumn;
                element[n++] = coefficient;
            }
        }
        // First sort in column order
        CoinSort_3(iColumn, iColumn + n, jColumn, element);
        // marker at end
        iColumn[n] = numberColumns;
        int lastI = iColumn[0];
        // compute starts
        start[0] = 0;
        for (i = 1; i < n + 1; i++) {
            if (iColumn[i] != lastI) {
                while (lastI < iColumn[i]) {
                    start[lastI+1] = i;
                    lastI++;
                }
                lastI = iColumn[i];
            }
        }
        // -1 unknown, 0 convex, 1 nonconvex
        int status = -1;
        int statusNegative = -1;
        int numberLong = 0; // number with >2 elements
        for (int k = 0; k < numberColumns; k++) {
            int first = start[k];
            int last = start[k+1];
            if (last > first) {
                int j;
                double diagonal = 0.0;
                int whichK = -1;
                for (j = first; j < last; j++) {
                    if (jColumn[j] == k) {
                        diagonal = element[j];
                        status = diagonal > 0 ? 0 : 1;
                        statusNegative = diagonal < 0 ? 0 : 1;
                        whichK = (j == first) ? j + 1 : j - 1;
                        break;
                    }
                }
                if (last == first + 1) {
                    // just one entry
                    if (!diagonal) {
                        // one off diagonal - not positive semi definite
                        status = 1;
                        statusNegative = 1;
                    }
                } else if (diagonal) {
                    if (last == first + 2) {
                        // other column and element
                        double otherElement = element[whichK];;
                        int otherColumn = jColumn[whichK];
                        double otherDiagonal = 0.0;
                        // check 2x2 determinant - unless past and 2 long
                        if (otherColumn > i || start[otherColumn+1] > start[otherColumn] + 2) {
                            for (j = start[otherColumn]; j < start[otherColumn+1]; j++) {
                                if (jColumn[j] == otherColumn) {
                                    otherDiagonal = element[j];
                                    break;
                                }
                            }
                            // determinant
                            double determinant = diagonal * otherDiagonal - otherElement * otherElement;
                            if (determinant < -1.0e-12) {
                                // not positive semi definite
                                status = 1;
                                statusNegative = 1;
                            } else if (start[otherColumn+1] > start[otherColumn] + 2 && determinant < 1.0e-12) {
                                // not positive semi definite
                                status = 1;
                                statusNegative = 1;
                            }
                        }
                    } else {
                        numberLong++;
                    }
                }
            }
        }
        if ((status == 0 || statusNegative == 0) && numberLong) {
            // need to do more work
            //printf("Needs more work\n");
        }
        assert (status > 0 || statusNegative > 0);
        if (!status) {
            convex_[iNon] = 1;
            // equality may be ok
            if (rowUpper[iRow] < 1.0e20)
                specialOptions2_ |= 8;
            else
                convex_[iNon] = 0;
        } else if (!statusNegative) {
            convex_[iNon] = -1;
            // equality may be ok
            if (rowLower[iRow] > -1.0e20)
                specialOptions2_ |= 8;
            else
                convex_[iNon] = 0;
        } else {
            convex_[iNon] = 0;
        }
        //printf("Convexity of row %d is %d\n",iRow,convex_[iNon]);
        delete [] iColumn;
        delete [] jColumn;
        delete [] element;
    }
    delete [] start;
}
//-------------------------------------------------------------------
// Clone
//-------------------------------------------------------------------
OsiSolverInterface *
OsiSolverLink::clone(bool /*copyData*/) const
{
    //assert (copyData);
    OsiSolverLink * newModel = new OsiSolverLink(*this);
    return newModel;
}


//-------------------------------------------------------------------
// Copy constructor
//-------------------------------------------------------------------
OsiSolverLink::OsiSolverLink (
    const OsiSolverLink & rhs)
        : OsiSolverInterface(rhs),
        CbcOsiSolver(rhs)
{
    gutsOfDestructor(true);
    gutsOfCopy(rhs);
    // something odd happens - try this
    OsiSolverInterface::operator=(rhs);
}

//-------------------------------------------------------------------
// Destructor
//-------------------------------------------------------------------
OsiSolverLink::~OsiSolverLink ()
{
    gutsOfDestructor();
}

//-------------------------------------------------------------------
// Assignment operator
//-------------------------------------------------------------------
OsiSolverLink &
OsiSolverLink::operator=(const OsiSolverLink & rhs)
{
    if (this != &rhs) {
        gutsOfDestructor();
        CbcOsiSolver::operator=(rhs);
        gutsOfCopy(rhs);
    }
    return *this;
}
void
OsiSolverLink::gutsOfDestructor(bool justNullify)
{
    if (!justNullify) {
        delete matrix_;
        delete originalRowCopy_;
        delete [] info_;
        delete [] bestSolution_;
        delete quadraticModel_;
        delete [] startNonLinear_;
        delete [] rowNonLinear_;
        delete [] convex_;
        delete [] whichNonLinear_;
        delete [] fixVariables_;
    }
    matrix_ = NULL;
    originalRowCopy_ = NULL;
    quadraticModel_ = NULL;
    numberNonLinearRows_ = 0;
    startNonLinear_ = NULL;
    rowNonLinear_ = NULL;
    convex_ = NULL;
    whichNonLinear_ = NULL;
    info_ = NULL;
    fixVariables_ = NULL;
    numberVariables_ = 0;
    specialOptions2_ = 0;
    objectiveRow_ = -1;
    objectiveVariable_ = -1;
    bestSolution_ = NULL;
    bestObjectiveValue_ = 1.0e100;
    defaultMeshSize_ = 1.0e-4;
    defaultBound_ = 1.0e5;
    integerPriority_ = 1000;
    biLinearPriority_ = 10000;
    numberFix_ = 0;
}
void
OsiSolverLink::gutsOfCopy(const OsiSolverLink & rhs)
{
    coinModel_ = rhs.coinModel_;
    numberVariables_ = rhs.numberVariables_;
    numberNonLinearRows_ = rhs.numberNonLinearRows_;
    specialOptions2_ = rhs.specialOptions2_;
    objectiveRow_ = rhs.objectiveRow_;
    objectiveVariable_ = rhs.objectiveVariable_;
    bestObjectiveValue_ = rhs.bestObjectiveValue_;
    defaultMeshSize_ = rhs.defaultMeshSize_;
    defaultBound_ = rhs.defaultBound_;
    integerPriority_ = rhs.integerPriority_;
    biLinearPriority_ = rhs.biLinearPriority_;
    numberFix_ = rhs.numberFix_;
    if (numberVariables_) {
        if (rhs.matrix_)
            matrix_ = new CoinPackedMatrix(*rhs.matrix_);
        else
            matrix_ = NULL;
        if (rhs.originalRowCopy_)
            originalRowCopy_ = new CoinPackedMatrix(*rhs.originalRowCopy_);
        else
            originalRowCopy_ = NULL;
        info_ = new OsiLinkedBound [numberVariables_];
        for (int i = 0; i < numberVariables_; i++) {
            info_[i] = OsiLinkedBound(rhs.info_[i]);
        }
        if (rhs.bestSolution_) {
            bestSolution_ = CoinCopyOfArray(rhs.bestSolution_, modelPtr_->getNumCols());
        } else {
            bestSolution_ = NULL;
        }
    }
    if (numberNonLinearRows_) {
        startNonLinear_ = CoinCopyOfArray(rhs.startNonLinear_, numberNonLinearRows_ + 1);
        rowNonLinear_ = CoinCopyOfArray(rhs.rowNonLinear_, numberNonLinearRows_);
        convex_ = CoinCopyOfArray(rhs.convex_, numberNonLinearRows_);
        int numberEntries = startNonLinear_[numberNonLinearRows_];
        whichNonLinear_ = CoinCopyOfArray(rhs.whichNonLinear_, numberEntries);
    }
    if (rhs.quadraticModel_) {
        quadraticModel_ = new ClpSimplex(*rhs.quadraticModel_);
    } else {
        quadraticModel_ = NULL;
    }
    fixVariables_ = CoinCopyOfArray(rhs.fixVariables_, numberFix_);
}
// Add a bound modifier
void
OsiSolverLink::addBoundModifier(bool upperBoundAffected, bool useUpperBound, int whichVariable, int whichVariableAffected,
                                double multiplier)
{
    bool found = false;
    int i;
    for ( i = 0; i < numberVariables_; i++) {
        if (info_[i].variable() == whichVariable) {
            found = true;
            break;
        }
    }
    if (!found) {
        // add in
        OsiLinkedBound * temp = new OsiLinkedBound [numberVariables_+1];
        for (int i = 0; i < numberVariables_; i++)
            temp[i] = info_[i];
        delete [] info_;
        info_ = temp;
        info_[numberVariables_++] = OsiLinkedBound(this, whichVariable, 0, NULL, NULL, NULL);
    }
    info_[i].addBoundModifier(upperBoundAffected, useUpperBound, whichVariableAffected, multiplier);
}
// Update coefficients
int
OsiSolverLink::updateCoefficients(ClpSimplex * solver, CoinPackedMatrix * matrix)
{
    double * lower = solver->columnLower();
    double * upper = solver->columnUpper();
    double * objective = solver->objective();
    int numberChanged = 0;
    for (int iObject = 0; iObject < numberObjects_; iObject++) {
        OsiBiLinear * obj = dynamic_cast<OsiBiLinear *> (object_[iObject]);
        if (obj) {
            numberChanged += obj->updateCoefficients(lower, upper, objective, matrix, &basis_);
        }
    }
    return numberChanged;
}
// Set best solution found internally
void
OsiSolverLink::setBestSolution(const double * solution, int numberColumns)
{
    delete [] bestSolution_;
    int numberColumnsThis = modelPtr_->numberColumns();
    bestSolution_ = new double [numberColumnsThis];
    CoinZeroN(bestSolution_, numberColumnsThis);
    memcpy(bestSolution_, solution, CoinMin(numberColumns, numberColumnsThis)*sizeof(double));
}
/* Two tier integer problem where when set of variables with priority
   less than this are fixed the problem becomes an easier integer problem
*/
void
OsiSolverLink::setFixedPriority(int priorityValue)
{
    delete [] fixVariables_;
    fixVariables_ = NULL;
    numberFix_ = 0;
    int i;
    for ( i = 0; i < numberObjects_; i++) {
        OsiSimpleInteger * obj = dynamic_cast<OsiSimpleInteger *> (object_[i]);
        if (obj) {
#ifndef NDEBUG
            int iColumn = obj->columnNumber();
            assert (iColumn >= 0);
#endif
            if (obj->priority() < priorityValue)
                numberFix_++;
        }
    }
    if (numberFix_) {
        specialOptions2_ |= 1;
        fixVariables_ = new int [numberFix_];
        numberFix_ = 0;
        // need to make sure coinModel_ is correct
        int numberColumns = coinModel_.numberColumns();
        char * highPriority = new char [numberColumns];
        CoinZeroN(highPriority, numberColumns);
        for ( i = 0; i < numberObjects_; i++) {
            OsiSimpleInteger * obj = dynamic_cast<OsiSimpleInteger *> (object_[i]);
            if (obj) {
                int iColumn = obj->columnNumber();
                assert (iColumn >= 0);
                if (iColumn < numberColumns) {
                    if (obj->priority() < priorityValue) {
                        object_[i] = new OsiSimpleFixedInteger(*obj);
                        delete obj;
                        fixVariables_[numberFix_++] = iColumn;
                        highPriority[iColumn] = 1;
                    }
                }
            }
        }
        CoinModel * newModel = coinModel_.reorder(highPriority);
        if (newModel) {
            coinModel_ = * newModel;
        } else {
            printf("Unable to use priorities\n");
            delete [] fixVariables_;
            fixVariables_ = NULL;
            numberFix_ = 0;
        }
        delete newModel;
        delete [] highPriority;
    }
}
// Gets correct form for a quadratic row - user to delete
CoinPackedMatrix *
OsiSolverLink::quadraticRow(int rowNumber, double * linearRow) const
{
    int numberColumns = coinModel_.numberColumns();
    CoinZeroN(linearRow, numberColumns);
    int numberElements = 0;
#ifndef NDEBUG
    int numberRows = coinModel_.numberRows();
    assert (rowNumber >= 0 && rowNumber < numberRows);
#endif
    CoinModelLink triple = coinModel_.firstInRow(rowNumber);
    while (triple.column() >= 0) {
        int iColumn = triple.column();
        const char * expr = coinModel_.getElementAsString(rowNumber, iColumn);
        if (strcmp(expr, "Numeric")) {
            // try and see which columns
            assert (strlen(expr) < 20000);
            char temp[20000];
            strcpy(temp, expr);
            char * pos = temp;
            bool ifFirst = true;
            while (*pos) {
                double value;
                int jColumn = decodeBit(pos, pos, value, ifFirst, coinModel_);
                // must be column unless first when may be linear term
                if (jColumn >= 0) {
                    numberElements++;
                } else if (jColumn == -2) {
                    linearRow[iColumn] = value;
                } else {
                    printf("bad nonlinear term %s\n", temp);
                    abort();
                }
                ifFirst = false;
            }
        } else {
            linearRow[iColumn] = coinModel_.getElement(rowNumber, iColumn);
        }
        triple = coinModel_.next(triple);
    }
    if (!numberElements) {
        return NULL;
    } else {
        int * column = new int[numberElements];
        int * column2 = new int[numberElements];
        double * element = new double[numberElements];
        numberElements = 0;
        CoinModelLink triple = coinModel_.firstInRow(rowNumber);
        while (triple.column() >= 0) {
            int iColumn = triple.column();
            const char * expr = coinModel_.getElementAsString(rowNumber, iColumn);
            if (strcmp(expr, "Numeric")) {
                // try and see which columns
                assert (strlen(expr) < 20000);
                char temp[20000];
                strcpy(temp, expr);
                char * pos = temp;
                bool ifFirst = true;
                while (*pos) {
                    double value;
                    int jColumn = decodeBit(pos, pos, value, ifFirst, coinModel_);
                    // must be column unless first when may be linear term
                    if (jColumn >= 0) {
                        column[numberElements] = iColumn;
                        column2[numberElements] = jColumn;
                        element[numberElements++] = value;
                    } else if (jColumn != -2) {
                        printf("bad nonlinear term %s\n", temp);
                        abort();
                    }
                    ifFirst = false;
                }
            }
            triple = coinModel_.next(triple);
        }
        return new CoinPackedMatrix(true, column2, column, element, numberElements);
    }
}
/*
  Problem specific
  Returns -1 if node fathomed and no solution
  0 if did nothing
  1 if node fathomed and solution
  allFixed is true if all LinkedBound variables are fixed
*/
int
OsiSolverLink::fathom(bool allFixed)
{
    int returnCode = 0;
    if (allFixed) {
        // solve anyway
        OsiClpSolverInterface::resolve();
        if (!isProvenOptimal()) {
            printf("cutoff before fathoming\n");
            return -1;
        }
        // all fixed so we can reformulate
        OsiClpSolverInterface newSolver;
        // set values
        const double * lower = modelPtr_->columnLower();
        const double * upper = modelPtr_->columnUpper();
        int i;
        for (i = 0; i < numberFix_; i++ ) {
            int iColumn = fixVariables_[i];
            double lo = lower[iColumn];
#ifndef NDEBUG
            double up = upper[iColumn];
            assert (lo == up);
#endif
            //printf("column %d fixed to %g\n",iColumn,lo);
            coinModel_.associateElement(coinModel_.columnName(iColumn), lo);
        }
        newSolver.loadFromCoinModel(coinModel_, true);
        for (i = 0; i < numberFix_; i++ ) {
            int iColumn = fixVariables_[i];
            newSolver.setColLower(iColumn, lower[iColumn]);
            newSolver.setColUpper(iColumn, lower[iColumn]);
        }
        // see if everything with objective fixed
        const double * objective = modelPtr_->objective();
        int numberColumns = newSolver.getNumCols();
        bool zeroObjective = true;
        double sum = 0.0;
        for (i = 0; i < numberColumns; i++) {
            if (upper[i] > lower[i] && objective[i]) {
                zeroObjective = false;
                break;
            } else {
                sum += lower[i] * objective[i];
            }
        }
        int fake[] = {5, 4, 3, 2, 0, 0, 0};
        bool onOptimalPath = true;
        for (i = 0; i < 7; i++) {
            if (static_cast<int> (upper[i]) != fake[i])
                onOptimalPath = false;
        }
        if (onOptimalPath)
            printf("possible\n");
        if (zeroObjective) {
            // randomize objective
            ClpSimplex * clpModel = newSolver.getModelPtr();
            const double * element = clpModel->matrix()->getMutableElements();
            //const int * row = clpModel->matrix()->getIndices();
            const CoinBigIndex * columnStart = clpModel->matrix()->getVectorStarts();
            const int * columnLength = clpModel->matrix()->getVectorLengths();
            double * objective = clpModel->objective();
            for (i = 0; i < numberColumns; i++) {
                if (clpModel->isInteger(i)) {
                    double value = 0.0;
                    for (int j = columnStart[i]; j < columnStart[i] + columnLength[i]; j++) {
                        value += fabs(element[j]);
                    }
                    objective[i] = value;
                }
            }
        }
        //newSolver.writeMps("xx");
        CbcModel model(newSolver);
        // Now do requested saves and modifications
        CbcModel * cbcModel = & model;
        OsiSolverInterface * osiModel = model.solver();
        OsiClpSolverInterface * osiclpModel = dynamic_cast< OsiClpSolverInterface*> (osiModel);
        ClpSimplex * clpModel = osiclpModel->getModelPtr();
        CglProbing probing;
        probing.setMaxProbe(10);
        probing.setMaxLook(10);
        probing.setMaxElements(200);
        probing.setMaxProbeRoot(50);
        probing.setMaxLookRoot(10);
        probing.setRowCuts(3);
        probing.setRowCuts(0);
        probing.setUsingObjective(true);
        cbcModel->addCutGenerator(&probing, -1, "Probing", true, false, false, -100, -1, -1);

        CglGomory gomory;
        gomory.setLimitAtRoot(512);
        cbcModel->addCutGenerator(&gomory, -98, "Gomory", true, false, false, -100, -1, -1);

        CglKnapsackCover knapsackCover;
        cbcModel->addCutGenerator(&knapsackCover, -98, "KnapsackCover", true, false, false, -100, -1, -1);

        CglClique clique;
        clique.setStarCliqueReport(false);
        clique.setRowCliqueReport(false);
        clique.setMinViolation(0.1);
        cbcModel->addCutGenerator(&clique, -98, "Clique", true, false, false, -100, -1, -1);
        CglMixedIntegerRounding2 mixedIntegerRounding2;
        cbcModel->addCutGenerator(&mixedIntegerRounding2, -98, "MixedIntegerRounding2", true, false, false, -100, -1, -1);

        CglFlowCover flowCover;
        cbcModel->addCutGenerator(&flowCover, -98, "FlowCover", true, false, false, -100, -1, -1);

        CglTwomir twomir;
        twomir.setMaxElements(250);
        cbcModel->addCutGenerator(&twomir, -99, "Twomir", true, false, false, -100, -1, -1);
        cbcModel->cutGenerator(6)->setTiming(true);

        CbcHeuristicFPump heuristicFPump(*cbcModel);
        heuristicFPump.setWhen(1);
        heuristicFPump.setMaximumPasses(20);
        heuristicFPump.setDefaultRounding(0.5);
        cbcModel->addHeuristic(&heuristicFPump);

        CbcRounding rounding(*cbcModel);
        cbcModel->addHeuristic(&rounding);

        CbcHeuristicLocal heuristicLocal(*cbcModel);
        heuristicLocal.setSearchType(1);
        cbcModel->addHeuristic(&heuristicLocal);

        CbcHeuristicGreedyCover heuristicGreedyCover(*cbcModel);
        cbcModel->addHeuristic(&heuristicGreedyCover);

        CbcHeuristicGreedyEquality heuristicGreedyEquality(*cbcModel);
        cbcModel->addHeuristic(&heuristicGreedyEquality);

        CbcCompareDefault compare;
        cbcModel->setNodeComparison(compare);
        cbcModel->setNumberBeforeTrust(5);
        cbcModel->setSpecialOptions(2);
        cbcModel->messageHandler()->setLogLevel(1);
        cbcModel->setMaximumCutPassesAtRoot(-100);
        cbcModel->setMaximumCutPasses(1);
        cbcModel->setMinimumDrop(0.05);
        clpModel->setNumberIterations(1);
        // For branchAndBound this may help
        clpModel->defaultFactorizationFrequency();
        clpModel->setDualBound(6.71523e+07);
        clpModel->setPerturbation(50);
        osiclpModel->setSpecialOptions(193);
        osiclpModel->messageHandler()->setLogLevel(0);
        osiclpModel->setIntParam(OsiMaxNumIterationHotStart, 100);
        osiclpModel->setHintParam(OsiDoReducePrint, true, OsiHintTry);
        // You can save some time by switching off message building
        // clpModel->messagesPointer()->setDetailMessages(100,10000,(int *) NULL);
        // Solve

        cbcModel->initialSolve();
        //double cutoff = model_->getCutoff();
        if (zeroObjective || !cbcModel_)
            cbcModel->setCutoff(1.0e50);
        else
            cbcModel->setCutoff(cbcModel_->getCutoff());
        // to change exits
        bool isFeasible = false;
        int saveLogLevel = clpModel->logLevel();
        clpModel->setLogLevel(0);
        if (clpModel->tightenPrimalBounds() != 0) {
            clpModel->setLogLevel(saveLogLevel);
            returnCode = -1; // infeasible//std::cout<<"Problem is infeasible - tightenPrimalBounds!"<<std::endl;
        } else {
            clpModel->setLogLevel(saveLogLevel);
            clpModel->dual();  // clean up
            // compute some things using problem size
            cbcModel->setMinimumDrop(CoinMin(5.0e-2,
                                             fabs(cbcModel->getMinimizationObjValue())*1.0e-3 + 1.0e-4));
            if (cbcModel->getNumCols() < 500)
                cbcModel->setMaximumCutPassesAtRoot(-100); // always do 100 if possible
            else if (cbcModel->getNumCols() < 5000)
                cbcModel->setMaximumCutPassesAtRoot(100); // use minimum drop
            else
                cbcModel->setMaximumCutPassesAtRoot(20);
            cbcModel->setMaximumCutPasses(1);
            // Hand coded preprocessing
            CglPreProcess process;
            OsiSolverInterface * saveSolver = cbcModel->solver()->clone();
            // Tell solver we are in Branch and Cut
            saveSolver->setHintParam(OsiDoInBranchAndCut, true, OsiHintDo) ;
            // Default set of cut generators
            CglProbing generator1;
            generator1.setUsingObjective(true);
            generator1.setMaxPass(3);
            generator1.setMaxProbeRoot(saveSolver->getNumCols());
            generator1.setMaxElements(100);
            generator1.setMaxLookRoot(50);
            generator1.setRowCuts(3);
            // Add in generators
            process.addCutGenerator(&generator1);
            process.messageHandler()->setLogLevel(cbcModel->logLevel());
            OsiSolverInterface * solver2 =
                process.preProcessNonDefault(*saveSolver, 0, 10);
            // Tell solver we are not in Branch and Cut
            saveSolver->setHintParam(OsiDoInBranchAndCut, false, OsiHintDo) ;
            if (solver2)
                solver2->setHintParam(OsiDoInBranchAndCut, false, OsiHintDo) ;
            if (!solver2) {
                std::cout << "Pre-processing says infeasible!" << std::endl;
                delete saveSolver;
                returnCode = -1;
            } else {
                std::cout << "processed model has " << solver2->getNumRows()
                          << " rows, " << solver2->getNumCols()
                          << " and " << solver2->getNumElements() << std::endl;
                // we have to keep solver2 so pass clone
                solver2 = solver2->clone();
                //solver2->writeMps("intmodel");
                cbcModel->assignSolver(solver2);
                cbcModel->initialSolve();
                if (zeroObjective) {
                    cbcModel->setMaximumSolutions(1); // just getting a solution
#ifdef JJF_ZERO
                    OsiClpSolverInterface * osiclpModel = dynamic_cast< OsiClpSolverInterface*> (cbcModel->solver());
                    ClpSimplex * clpModel = osiclpModel->getModelPtr();
                    const double * element = clpModel->matrix()->getMutableElements();
                    //const int * row = clpModel->matrix()->getIndices();
                    const CoinBigIndex * columnStart = clpModel->matrix()->getVectorStarts();
                    const int * columnLength = clpModel->matrix()->getVectorLengths();
                    int n = clpModel->numberColumns();
                    int * sort2 = new int[n];
                    int * pri = new int[n];
                    double * sort = new double[n];
                    int i;
                    int nint = 0;
                    for (i = 0; i < n; i++) {
                        if (clpModel->isInteger(i)) {
                            double largest = 0.0;
                            for (int j = columnStart[i]; j < columnStart[i] + columnLength[i]; j++) {
                                largest = CoinMax(largest, fabs(element[j]));
                            }
                            sort2[nint] = nint;
                            sort[nint++] = -largest;
                        }
                    }
                    CoinSort_2(sort, sort + nint, sort2);
                    int kpri = 1;
                    double last = sort[0];
                    for (i = 0; i < nint; i++) {
                        if (sort[i] != last) {
                            kpri++;
                            last = sort[i];
                        }
                        pri[sort2[i]] = kpri;
                    }
                    cbcModel->passInPriorities(pri, false);
                    delete [] sort;
                    delete [] sort2;
                    delete [] pri;
#endif
                }
                cbcModel->branchAndBound();
                // For best solution
                int numberColumns = newSolver.getNumCols();
                if (cbcModel->getMinimizationObjValue() < 1.0e50) {
                    // post process
                    process.postProcess(*cbcModel->solver());
                    // Solution now back in saveSolver
                    cbcModel->assignSolver(saveSolver);
                    memcpy(cbcModel->bestSolution(), cbcModel->solver()->getColSolution(),
                           numberColumns*sizeof(double));
                    // put back in original solver
                    newSolver.setColSolution(cbcModel->bestSolution());
                    isFeasible = true;
                } else {
                    delete saveSolver;
                }
            }
            //const double * solution = newSolver.getColSolution();
            if (isFeasible && cbcModel->getMinimizationObjValue() < 1.0e50) {
                int numberColumns = this->getNumCols();
                int i;
                const double * solution = cbcModel->bestSolution();
                int numberColumns2 = newSolver.getNumCols();
                for (i = 0; i < numberColumns2; i++) {
                    double value = solution[i];
                    assert (fabs(value - floor(value + 0.5)) < 0.0001);
                    value = floor(value + 0.5);
                    this->setColLower(i, value);
                    this->setColUpper(i, value);
                }
                for (; i < numberColumns; i++) {
                    this->setColLower(i, 0.0);
                    this->setColUpper(i, 1.1);
                }
                // but take off cuts
                int numberRows = getNumRows();
                int numberRows2 = cbcModel_->continuousSolver()->getNumRows();

                for (i = numberRows2; i < numberRows; i++)
                    setRowBounds(i, -COIN_DBL_MAX, COIN_DBL_MAX);
                initialSolve();
                //if (!isProvenOptimal())
                //getModelPtr()->writeMps("bad.mps");
                if (isProvenOptimal()) {
                    delete [] bestSolution_;
                    bestSolution_ = CoinCopyOfArray(modelPtr_->getColSolution(), modelPtr_->getNumCols());
                    bestObjectiveValue_ = modelPtr_->objectiveValue();
                    printf("BB best value %g\n", bestObjectiveValue_);
                    returnCode = 1;
                } else {
                    printf("*** WHY BAD SOL\n");
                    returnCode = -1;
                }
            } else {
                modelPtr_->setProblemStatus(1);
                modelPtr_->setObjectiveValue(COIN_DBL_MAX);
                returnCode = -1;
            }
        }
    }
    return returnCode;
}
//#############################################################################
// Constructors, destructors  and assignment
//#############################################################################

//-------------------------------------------------------------------
// Default Constructor
//-------------------------------------------------------------------
OsiLinkedBound::OsiLinkedBound ()
{
    model_ = NULL;
    variable_ = -1;
    numberAffected_ = 0;
    maximumAffected_ = numberAffected_;
    affected_ = NULL;
}
// Useful Constructor
OsiLinkedBound::OsiLinkedBound(OsiSolverInterface * model, int variable,
                               int numberAffected, const int * positionL,
                               const int * positionU, const double * multiplier)
{
    model_ = model;
    variable_ = variable;
    numberAffected_ = 2 * numberAffected;
    maximumAffected_ = numberAffected_;
    if (numberAffected_) {
        affected_ = new boundElementAction[numberAffected_];
        int n = 0;
        for (int i = 0; i < numberAffected; i++) {
            // LB
            boundElementAction action;
            action.affect = 2;
            action.ubUsed = 0;
            action.type = 0;
            action.affected = positionL[i];
            action.multiplier = multiplier[i];
            affected_[n++] = action;
            // UB
            action.affect = 2;
            action.ubUsed = 1;
            action.type = 0;
            action.affected = positionU[i];
            action.multiplier = multiplier[i];
            affected_[n++] = action;
        }
    } else {
        affected_ = NULL;
    }
}

//-------------------------------------------------------------------
// Copy constructor
//-------------------------------------------------------------------
OsiLinkedBound::OsiLinkedBound (
    const OsiLinkedBound & rhs)
{
    model_ = rhs.model_;
    variable_ = rhs.variable_;
    numberAffected_ = rhs.numberAffected_;
    maximumAffected_ = rhs.maximumAffected_;
    if (numberAffected_) {
        affected_ = new boundElementAction[maximumAffected_];
        memcpy(affected_, rhs.affected_, numberAffected_*sizeof(boundElementAction));
    } else {
        affected_ = NULL;
    }
}

//-------------------------------------------------------------------
// Destructor
//-------------------------------------------------------------------
OsiLinkedBound::~OsiLinkedBound ()
{
    delete [] affected_;
}

//-------------------------------------------------------------------
// Assignment operator
//-------------------------------------------------------------------
OsiLinkedBound &
OsiLinkedBound::operator=(const OsiLinkedBound & rhs)
{
    if (this != &rhs) {
        delete [] affected_;
        model_ = rhs.model_;
        variable_ = rhs.variable_;
        numberAffected_ = rhs.numberAffected_;
        maximumAffected_ = rhs.maximumAffected_;
        if (numberAffected_) {
            affected_ = new boundElementAction[maximumAffected_];
            memcpy(affected_, rhs.affected_, numberAffected_*sizeof(boundElementAction));
        } else {
            affected_ = NULL;
        }
    }
    return *this;
}
// Add a bound modifier
void
OsiLinkedBound::addBoundModifier(bool upperBoundAffected, bool useUpperBound, int whichVariable,
                                 double multiplier)
{
    if (numberAffected_ == maximumAffected_) {
        maximumAffected_ = maximumAffected_ + 10 + maximumAffected_ / 4;
        boundElementAction * temp = new boundElementAction[maximumAffected_];
        memcpy(temp, affected_, numberAffected_*sizeof(boundElementAction));
        delete [] affected_;
        affected_ = temp;
    }
    boundElementAction action;
    action.affect = static_cast<unsigned char>(upperBoundAffected ? 1 : 0);
    action.ubUsed = static_cast<unsigned char>(useUpperBound ? 1 : 0);
    action.type = 2;
    action.affected = static_cast<short int>(whichVariable);
    action.multiplier = multiplier;
    affected_[numberAffected_++] = action;

}
// Update other bounds
void
OsiLinkedBound::updateBounds(ClpSimplex * solver)
{
    double * lower = solver->columnLower();
    double * upper = solver->columnUpper();
    double lo = lower[variable_];
    double up = upper[variable_];
    // printf("bounds for %d are %g and %g\n",variable_,lo,up);
    for (int j = 0; j < numberAffected_; j++) {
        if (affected_[j].affect < 2) {
            double multiplier = affected_[j].multiplier;
            assert (affected_[j].type == 2);
            int iColumn = affected_[j].affected;
            double useValue = (affected_[j].ubUsed) ? up : lo;
            if (affected_[j].affect == 0)
                lower[iColumn] = CoinMin(upper[iColumn], CoinMax(lower[iColumn], multiplier * useValue));
            else
                upper[iColumn] = CoinMax(lower[iColumn], CoinMin(upper[iColumn], multiplier * useValue));
        }
    }
}
#ifdef JJF_ZERO
// Add an element modifier
void
OsiLinkedBound::addCoefficientModifier(bool useUpperBound, int position,
                                       double multiplier)
{
    if (numberAffected_ == maximumAffected_) {
        maximumAffected_ = maximumAffected_ + 10 + maximumAffected_ / 4;
        boundElementAction * temp = new boundElementAction[maximumAffected_];
        memcpy(temp, affected_, numberAffected_*sizeof(boundElementAction));
        delete [] affected_;
        affected_ = temp;
    }
    boundElementAction action;
    action.affect = 2;
    action.ubUsed = useUpperBound ? 1 : 0;
    action.type = 0;
    action.affected = position;
    action.multiplier = multiplier;
    affected_[numberAffected_++] = action;

}
// Update coefficients
void
OsiLinkedBound::updateCoefficients(ClpSimplex * solver, CoinPackedMatrix * matrix)
{
    double * lower = solver->columnLower();
    double * upper = solver->columnUpper();
    double * element = matrix->getMutableElements();
    double lo = lower[variable_];
    double up = upper[variable_];
    // printf("bounds for %d are %g and %g\n",variable_,lo,up);
    for (int j = 0; j < numberAffected_; j++) {
        if (affected_[j].affect == 2) {
            double multiplier = affected_[j].multiplier;
            assert (affected_[j].type == 0);
            int position = affected_[j].affected;
            //double old = element[position];
            if (affected_[j].ubUsed)
                element[position] = multiplier * up;
            else
                element[position] = multiplier * lo;
            //if ( old != element[position])
            //printf("change at %d from %g to %g\n",position,old,element[position]);
        }
    }
}
#endif
// Default Constructor
CbcHeuristicDynamic3::CbcHeuristicDynamic3()
        : CbcHeuristic()
{
}

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

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

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

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

// Returns 1 if solution, 0 if not
int
CbcHeuristicDynamic3::solution(double & solutionValue,
                               double * betterSolution)
{
    if (!model_)
        return 0;
    OsiSolverLink * clpSolver
    = dynamic_cast<OsiSolverLink *> (model_->solver());
    assert (clpSolver);
    double newSolutionValue = clpSolver->bestObjectiveValue();
    const double * solution = clpSolver->bestSolution();
    if (newSolutionValue < solutionValue && solution) {
        int numberColumns = clpSolver->getNumCols();
        // new solution
        memcpy(betterSolution, solution, numberColumns*sizeof(double));
        solutionValue = newSolutionValue;
        return 1;
    } else {
        return 0;
    }
}
// update model
void CbcHeuristicDynamic3::setModel(CbcModel * model)
{
    model_ = model;
}
// Resets stuff if model changes
void
CbcHeuristicDynamic3::resetModel(CbcModel * model)
{
    model_ = model;
}
#include <cassert>
#include <cmath>
#include <cfloat>
//#define CBC_DEBUG

#include "OsiSolverInterface.hpp"
//#include "OsiBranchLink.hpp"
#include "CoinError.hpp"
#include "CoinHelperFunctions.hpp"
#include "CoinPackedMatrix.hpp"
#include "CoinWarmStartBasis.hpp"

// Default Constructor
OsiOldLink::OsiOldLink ()
        : OsiSOS(),
        numberLinks_(0)
{
}

// Useful constructor (which are indices)
OsiOldLink::OsiOldLink (const OsiSolverInterface * /*solver*/,  int numberMembers,
                        int numberLinks, int first , const double * weights, int /*identifier*/)
        : OsiSOS(),
        numberLinks_(numberLinks)
{
    numberMembers_ = numberMembers;
    members_ = NULL;
    sosType_ = 1;
    if (numberMembers_) {
        weights_ = new double[numberMembers_];
        members_ = new int[numberMembers_*numberLinks_];
        if (weights) {
            memcpy(weights_, weights, numberMembers_*sizeof(double));
        } else {
            for (int i = 0; i < numberMembers_; i++)
                weights_[i] = i;
        }
        // weights must be increasing
        int i;
#ifndef NDEBUG
        for (i = 1; i < numberMembers_; i++)
            assert (weights_[i] > weights_[i-1] + 1.0e-12);
#endif
        for (i = 0; i < numberMembers_*numberLinks_; i++) {
            members_[i] = first + i;
        }
    } else {
        weights_ = NULL;
    }
}

// Useful constructor (which are indices)
OsiOldLink::OsiOldLink (const OsiSolverInterface * /*solver*/,  int numberMembers,
                        int numberLinks, int /*sosType*/, const int * which ,
                        const double * weights, int /*identifier*/)
        : OsiSOS(),
        numberLinks_(numberLinks)
{
    numberMembers_ = numberMembers;
    members_ = NULL;
    sosType_ = 1;
    if (numberMembers_) {
        weights_ = new double[numberMembers_];
        members_ = new int[numberMembers_*numberLinks_];
        if (weights) {
            memcpy(weights_, weights, numberMembers_*sizeof(double));
        } else {
            for (int i = 0; i < numberMembers_; i++)
                weights_[i] = i;
        }
        // weights must be increasing
        int i;
#ifndef NDEBUG
        for (i = 1; i < numberMembers_; i++)
            assert (weights_[i] > weights_[i-1] + 1.0e-12);
#endif
        for (i = 0; i < numberMembers_*numberLinks_; i++) {
            members_[i] = which[i];
        }
    } else {
        weights_ = NULL;
    }
}

// Copy constructor
OsiOldLink::OsiOldLink ( const OsiOldLink & rhs)
        : OsiSOS(rhs)
{
    numberLinks_ = rhs.numberLinks_;
    if (numberMembers_) {
        delete [] members_;
        members_ = CoinCopyOfArray(rhs.members_, numberMembers_ * numberLinks_);
    }
}

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

// Assignment operator
OsiOldLink &
OsiOldLink::operator=( const OsiOldLink & rhs)
{
    if (this != &rhs) {
        OsiSOS::operator=(rhs);
        delete [] members_;
        numberLinks_ = rhs.numberLinks_;
        if (numberMembers_) {
            members_ = CoinCopyOfArray(rhs.members_, numberMembers_ * numberLinks_);
        } else {
            members_ = NULL;
        }
    }
    return *this;
}

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

// Infeasibility - large is 0.5
double
OsiOldLink::infeasibility(const OsiBranchingInformation * info, int & whichWay) const
{
    int j;
    int firstNonZero = -1;
    int lastNonZero = -1;
    const double * solution = info->solution_;
    //const double * lower = info->lower_;
    const double * upper = info->upper_;
    double integerTolerance = info->integerTolerance_;
    double weight = 0.0;
    double sum = 0.0;

    // check bounds etc
    double lastWeight = -1.0e100;
    int base = 0;
    for (j = 0; j < numberMembers_; j++) {
        for (int k = 0; k < numberLinks_; k++) {
            int iColumn = members_[base+k];
            if (lastWeight >= weights_[j] - 1.0e-7)
                throw CoinError("Weights too close together in OsiLink", "infeasibility", "OsiLink");
            lastWeight = weights_[j];
            double value = CoinMax(0.0, solution[iColumn]);
            sum += value;
            if (value > integerTolerance && upper[iColumn]) {
                // Possibly due to scaling a fixed variable might slip through
                if (value > upper[iColumn] + 1.0e-8) {
#ifdef OSI_DEBUG
                    printf("** Variable %d (%d) has value %g and upper bound of %g\n",
                           iColumn, j, value, upper[iColumn]);
#endif
                }
                value = CoinMin(value, upper[iColumn]);
                weight += weights_[j] * value;
                if (firstNonZero < 0)
                    firstNonZero = j;
                lastNonZero = j;
            }
        }
        base += numberLinks_;
    }
    double valueInfeasibility;
    whichWay = 1;
    whichWay_ = 1;
    if (lastNonZero - firstNonZero >= sosType_) {
        // find where to branch
        assert (sum > 0.0);
        weight /= sum;
        valueInfeasibility = lastNonZero - firstNonZero + 1;
        valueInfeasibility *= 0.5 / static_cast<double> (numberMembers_);
        //#define DISTANCE
#ifdef DISTANCE
        assert (sosType_ == 1); // code up
        /* may still be satisfied.
           For LOS type 2 we might wish to move coding around
           and keep initial info in model_ for speed
        */
        int iWhere;
        bool possible = false;
        for (iWhere = firstNonZero; iWhere <= lastNonZero; iWhere++) {
            if (fabs(weight - weights_[iWhere]) < 1.0e-8) {
                possible = true;
                break;
            }
        }
        if (possible) {
            // One could move some of this (+ arrays) into model_
            const CoinPackedMatrix * matrix = solver->getMatrixByCol();
            const double * element = matrix->getMutableElements();
            const int * row = matrix->getIndices();
            const CoinBigIndex * columnStart = matrix->getVectorStarts();
            const int * columnLength = matrix->getVectorLengths();
            const double * rowSolution = solver->getRowActivity();
            const double * rowLower = solver->getRowLower();
            const double * rowUpper = solver->getRowUpper();
            int numberRows = matrix->getNumRows();
            double * array = new double [numberRows];
            CoinZeroN(array, numberRows);
            int * which = new int [numberRows];
            int n = 0;
            int base = numberLinks_ * firstNonZero;
            for (j = firstNonZero; j <= lastNonZero; j++) {
                for (int k = 0; k < numberLinks_; k++) {
                    int iColumn = members_[base+k];
                    double value = CoinMax(0.0, solution[iColumn]);
                    if (value > integerTolerance && upper[iColumn]) {
                        value = CoinMin(value, upper[iColumn]);
                        for (int j = columnStart[iColumn]; j < columnStart[iColumn] + columnLength[iColumn]; j++) {
                            int iRow = row[j];
                            double a = array[iRow];
                            if (a) {
                                a += value * element[j];
                                if (!a)
                                    a = 1.0e-100;
                            } else {
                                which[n++] = iRow;
                                a = value * element[j];
                                assert (a);
                            }
                            array[iRow] = a;
                        }
                    }
                }
                base += numberLinks_;
            }
            base = numberLinks_ * iWhere;
            for (int k = 0; k < numberLinks_; k++) {
                int iColumn = members_[base+k];
                const double value = 1.0;
                for (int j = columnStart[iColumn]; j < columnStart[iColumn] + columnLength[iColumn]; j++) {
                    int iRow = row[j];
                    double a = array[iRow];
                    if (a) {
                        a -= value * element[j];
                        if (!a)
                            a = 1.0e-100;
                    } else {
                        which[n++] = iRow;
                        a = -value * element[j];
                        assert (a);
                    }
                    array[iRow] = a;
                }
            }
            for (j = 0; j < n; j++) {
                int iRow = which[j];
                // moving to point will increase row solution by this
                double distance = array[iRow];
                if (distance > 1.0e-8) {
                    if (distance + rowSolution[iRow] > rowUpper[iRow] + 1.0e-8) {
                        possible = false;
                        break;
                    }
                } else if (distance < -1.0e-8) {
                    if (distance + rowSolution[iRow] < rowLower[iRow] - 1.0e-8) {
                        possible = false;
                        break;
                    }
                }
            }
            for (j = 0; j < n; j++)
                array[which[j]] = 0.0;
            delete [] array;
            delete [] which;
            if (possible) {
                valueInfeasibility = 0.0;
                printf("possible %d %d %d\n", firstNonZero, lastNonZero, iWhere);
            }
        }
#endif
    } else {
        valueInfeasibility = 0.0; // satisfied
    }
    infeasibility_ = valueInfeasibility;
    otherInfeasibility_ = 1.0 - valueInfeasibility;
    return valueInfeasibility;
}

// This looks at solution and sets bounds to contain solution
double
OsiOldLink::feasibleRegion(OsiSolverInterface * solver, const OsiBranchingInformation * info) const
{
    int j;
    int firstNonZero = -1;
    int lastNonZero = -1;
    const double * solution = info->solution_;
    const double * upper = info->upper_;
    double integerTolerance = info->integerTolerance_;
    double weight = 0.0;
    double sum = 0.0;

    int base = 0;
    for (j = 0; j < numberMembers_; j++) {
        for (int k = 0; k < numberLinks_; k++) {
            int iColumn = members_[base+k];
            double value = CoinMax(0.0, solution[iColumn]);
            sum += value;
            if (value > integerTolerance && upper[iColumn]) {
                weight += weights_[j] * value;
                if (firstNonZero < 0)
                    firstNonZero = j;
                lastNonZero = j;
            }
        }
        base += numberLinks_;
    }
#ifdef DISTANCE
    if (lastNonZero - firstNonZero > sosType_ - 1) {
        /* may still be satisfied.
           For LOS type 2 we might wish to move coding around
           and keep initial info in model_ for speed
        */
        int iWhere;
        bool possible = false;
        for (iWhere = firstNonZero; iWhere <= lastNonZero; iWhere++) {
            if (fabs(weight - weights_[iWhere]) < 1.0e-8) {
                possible = true;
                break;
            }
        }
        if (possible) {
            // One could move some of this (+ arrays) into model_
            const CoinPackedMatrix * matrix = solver->getMatrixByCol();
            const double * element = matrix->getMutableElements();
            const int * row = matrix->getIndices();
            const CoinBigIndex * columnStart = matrix->getVectorStarts();
            const int * columnLength = matrix->getVectorLengths();
            const double * rowSolution = solver->getRowActivity();
            const double * rowLower = solver->getRowLower();
            const double * rowUpper = solver->getRowUpper();
            int numberRows = matrix->getNumRows();
            double * array = new double [numberRows];
            CoinZeroN(array, numberRows);
            int * which = new int [numberRows];
            int n = 0;
            int base = numberLinks_ * firstNonZero;
            for (j = firstNonZero; j <= lastNonZero; j++) {
                for (int k = 0; k < numberLinks_; k++) {
                    int iColumn = members_[base+k];
                    double value = CoinMax(0.0, solution[iColumn]);
                    if (value > integerTolerance && upper[iColumn]) {
                        value = CoinMin(value, upper[iColumn]);
                        for (int j = columnStart[iColumn]; j < columnStart[iColumn] + columnLength[iColumn]; j++) {
                            int iRow = row[j];
                            double a = array[iRow];
                            if (a) {
                                a += value * element[j];
                                if (!a)
                                    a = 1.0e-100;
                            } else {
                                which[n++] = iRow;
                                a = value * element[j];
                                assert (a);
                            }
                            array[iRow] = a;
                        }
                    }
                }
                base += numberLinks_;
            }
            base = numberLinks_ * iWhere;
            for (int k = 0; k < numberLinks_; k++) {
                int iColumn = members_[base+k];
                const double value = 1.0;
                for (int j = columnStart[iColumn]; j < columnStart[iColumn] + columnLength[iColumn]; j++) {
                    int iRow = row[j];
                    double a = array[iRow];
                    if (a) {
                        a -= value * element[j];
                        if (!a)
                            a = 1.0e-100;
                    } else {
                        which[n++] = iRow;
                        a = -value * element[j];
                        assert (a);
                    }
                    array[iRow] = a;
                }
            }
            for (j = 0; j < n; j++) {
                int iRow = which[j];
                // moving to point will increase row solution by this
                double distance = array[iRow];
                if (distance > 1.0e-8) {
                    if (distance + rowSolution[iRow] > rowUpper[iRow] + 1.0e-8) {
                        possible = false;
                        break;
                    }
                } else if (distance < -1.0e-8) {
                    if (distance + rowSolution[iRow] < rowLower[iRow] - 1.0e-8) {
                        possible = false;
                        break;
                    }
                }
            }
            for (j = 0; j < n; j++)
                array[which[j]] = 0.0;
            delete [] array;
            delete [] which;
            if (possible) {
                printf("possible feas region %d %d %d\n", firstNonZero, lastNonZero, iWhere);
                firstNonZero = iWhere;
                lastNonZero = iWhere;
            }
        }
    }
#else
    assert (lastNonZero - firstNonZero < sosType_) ;
#endif
    base = 0;
    for (j = 0; j < firstNonZero; j++) {
        for (int k = 0; k < numberLinks_; k++) {
            int iColumn = members_[base+k];
            solver->setColUpper(iColumn, 0.0);
        }
        base += numberLinks_;
    }
    // skip
    base += numberLinks_;
    for (j = lastNonZero + 1; j < numberMembers_; j++) {
        for (int k = 0; k < numberLinks_; k++) {
            int iColumn = members_[base+k];
            solver->setColUpper(iColumn, 0.0);
        }
        base += numberLinks_;
    }
    // go to coding as in OsiSOS
    abort();
    return -1.0;
}

// Redoes data when sequence numbers change
void
OsiOldLink::resetSequenceEtc(int numberColumns, const int * originalColumns)
{
    int n2 = 0;
    for (int j = 0; j < numberMembers_*numberLinks_; j++) {
        int iColumn = members_[j];
        int i;
#ifdef JJF_ZERO
        for (i = 0; i < numberColumns; i++) {
            if (originalColumns[i] == iColumn)
                break;
        }
#else
        i = originalColumns[iColumn];
#endif
        if (i >= 0 && i < numberColumns) {
            members_[n2] = i;
            weights_[n2++] = weights_[j];
        }
    }
    if (n2 < numberMembers_) {
        printf("** SOS number of members reduced from %d to %d!\n", numberMembers_, n2 / numberLinks_);
        numberMembers_ = n2 / numberLinks_;
    }
}

// Creates a branching object
OsiBranchingObject *
OsiOldLink::createBranch(OsiSolverInterface * solver, const OsiBranchingInformation * info, int way) const
{
    int j;
    const double * solution = info->solution_;
    double tolerance = info->primalTolerance_;
    const double * upper = info->upper_;
    int firstNonFixed = -1;
    int lastNonFixed = -1;
    int firstNonZero = -1;
    int lastNonZero = -1;
    double weight = 0.0;
    double sum = 0.0;
    int base = 0;
    for (j = 0; j < numberMembers_; j++) {
        for (int k = 0; k < numberLinks_; k++) {
            int iColumn = members_[base+k];
            if (upper[iColumn]) {
                double value = CoinMax(0.0, solution[iColumn]);
                sum += value;
                if (firstNonFixed < 0)
                    firstNonFixed = j;
                lastNonFixed = j;
                if (value > tolerance) {
                    weight += weights_[j] * value;
                    if (firstNonZero < 0)
                        firstNonZero = j;
                    lastNonZero = j;
                }
            }
        }
        base += numberLinks_;
    }
    assert (lastNonZero - firstNonZero >= sosType_) ;
    // find where to branch
    assert (sum > 0.0);
    weight /= sum;
    int iWhere;
    double separator = 0.0;
    for (iWhere = firstNonZero; iWhere < lastNonZero; iWhere++)
        if (weight < weights_[iWhere+1])
            break;
    if (sosType_ == 1) {
        // SOS 1
        separator = 0.5 * (weights_[iWhere] + weights_[iWhere+1]);
    } else {
        // SOS 2
        if (iWhere == firstNonFixed)
            iWhere++;;
        if (iWhere == lastNonFixed - 1)
            iWhere = lastNonFixed - 2;
        separator = weights_[iWhere+1];
    }
    // create object
    OsiBranchingObject * branch;
    branch = new OsiOldLinkBranchingObject(solver, this, way, separator);
    return branch;
}
OsiOldLinkBranchingObject::OsiOldLinkBranchingObject()
        : OsiSOSBranchingObject()
{
}

// Useful constructor
OsiOldLinkBranchingObject::OsiOldLinkBranchingObject (OsiSolverInterface * solver,
        const OsiOldLink * set,
        int way ,
        double separator)
        : OsiSOSBranchingObject(solver, set, way, separator)
{
}

// Copy constructor
OsiOldLinkBranchingObject::OsiOldLinkBranchingObject ( const OsiOldLinkBranchingObject & rhs) : OsiSOSBranchingObject(rhs)
{
}

// Assignment operator
OsiOldLinkBranchingObject &
OsiOldLinkBranchingObject::operator=( const OsiOldLinkBranchingObject & rhs)
{
    if (this != &rhs) {
        OsiSOSBranchingObject::operator=(rhs);
    }
    return *this;
}
OsiBranchingObject *
OsiOldLinkBranchingObject::clone() const
{
    return (new OsiOldLinkBranchingObject(*this));
}


// Destructor
OsiOldLinkBranchingObject::~OsiOldLinkBranchingObject ()
{
}
double
OsiOldLinkBranchingObject::branch(OsiSolverInterface * solver)
{
    const OsiOldLink * set =
        dynamic_cast <const OsiOldLink *>(originalObject_) ;
    assert (set);
    int way = (!branchIndex_) ? (2 * firstBranch_ - 1) : -(2 * firstBranch_ - 1);
    branchIndex_++;
    int numberMembers = set->numberMembers();
    const int * which = set->members();
    const double * weights = set->weights();
    int numberLinks = set->numberLinks();
    //const double * lower = info->lower_;
    //const double * upper = solver->getColUpper();
    // *** for way - up means fix all those in down section
    if (way < 0) {
        int i;
        for ( i = 0; i < numberMembers; i++) {
            if (weights[i] > value_)
                break;
        }
        assert (i < numberMembers);
        int base = i * numberLinks;;
        for (; i < numberMembers; i++) {
            for (int k = 0; k < numberLinks; k++) {
                int iColumn = which[base+k];
                solver->setColUpper(iColumn, 0.0);
            }
            base += numberLinks;
        }
    } else {
        int i;
        int base = 0;
        for ( i = 0; i < numberMembers; i++) {
            if (weights[i] >= value_) {
                break;
            } else {
                for (int k = 0; k < numberLinks; k++) {
                    int iColumn = which[base+k];
                    solver->setColUpper(iColumn, 0.0);
                }
                base += numberLinks;
            }
        }
        assert (i < numberMembers);
    }
    return 0.0;
}
// Print what would happen
void
OsiOldLinkBranchingObject::print(const OsiSolverInterface * solver)
{
    const OsiOldLink * set =
        dynamic_cast <const OsiOldLink *>(originalObject_) ;
    assert (set);
    int way = (!branchIndex_) ? (2 * firstBranch_ - 1) : -(2 * firstBranch_ - 1);
    int numberMembers = set->numberMembers();
    int numberLinks = set->numberLinks();
    const double * weights = set->weights();
    const int * which = set->members();
    const double * upper = solver->getColUpper();
    int first = numberMembers;
    int last = -1;
    int numberFixed = 0;
    int numberOther = 0;
    int i;
    int base = 0;
    for ( i = 0; i < numberMembers; i++) {
        for (int k = 0; k < numberLinks; k++) {
            int iColumn = which[base+k];
            double bound = upper[iColumn];
            if (bound) {
                first = CoinMin(first, i);
                last = CoinMax(last, i);
            }
        }
        base += numberLinks;
    }
    // *** for way - up means fix all those in down section
    base = 0;
    if (way < 0) {
        printf("SOS Down");
        for ( i = 0; i < numberMembers; i++) {
            if (weights[i] > value_)
                break;
            for (int k = 0; k < numberLinks; k++) {
                int iColumn = which[base+k];
                double bound = upper[iColumn];
                if (bound)
                    numberOther++;
            }
            base += numberLinks;
        }
        assert (i < numberMembers);
        for (; i < numberMembers; i++) {
            for (int k = 0; k < numberLinks; k++) {
                int iColumn = which[base+k];
                double bound = upper[iColumn];
                if (bound)
                    numberFixed++;
            }
            base += numberLinks;
        }
    } else {
        printf("SOS Up");
        for ( i = 0; i < numberMembers; i++) {
            if (weights[i] >= value_)
                break;
            for (int k = 0; k < numberLinks; k++) {
                int iColumn = which[base+k];
                double bound = upper[iColumn];
                if (bound)
                    numberFixed++;
            }
            base += numberLinks;
        }
        assert (i < numberMembers);
        for (; i < numberMembers; i++) {
            for (int k = 0; k < numberLinks; k++) {
                int iColumn = which[base+k];
                double bound = upper[iColumn];
                if (bound)
                    numberOther++;
            }
            base += numberLinks;
        }
    }
    assert ((numberFixed % numberLinks) == 0);
    assert ((numberOther % numberLinks) == 0);
    printf(" - at %g, free range %d (%g) => %d (%g), %d would be fixed, %d other way\n",
           value_, first, weights[first], last, weights[last], numberFixed / numberLinks,
           numberOther / numberLinks);
}
// Default Constructor
OsiBiLinear::OsiBiLinear ()
        : OsiObject2(),
        coefficient_(0.0),
        xMeshSize_(0.0),
        yMeshSize_(0.0),
        xSatisfied_(1.0e-6),
        ySatisfied_(1.0e-6),
        xOtherSatisfied_(0.0),
        yOtherSatisfied_(0.0),
        xySatisfied_(1.0e-6),
        xyBranchValue_(0.0),
        xColumn_(-1),
        yColumn_(-1),
        firstLambda_(-1),
        branchingStrategy_(0),
        boundType_(0),
        xRow_(-1),
        yRow_(-1),
        xyRow_(-1),
        convexity_(-1),
        numberExtraRows_(0),
        multiplier_(NULL),
        extraRow_(NULL),
        chosen_(-1)
{
}

// Useful constructor
OsiBiLinear::OsiBiLinear (OsiSolverInterface * solver, int xColumn,
                          int yColumn, int xyRow, double coefficient,
                          double xMesh, double yMesh,
                          int numberExistingObjects, const OsiObject ** objects )
        : OsiObject2(),
        coefficient_(coefficient),
        xMeshSize_(xMesh),
        yMeshSize_(yMesh),
        xSatisfied_(1.0e-6),
        ySatisfied_(1.0e-6),
        xOtherSatisfied_(0.0),
        yOtherSatisfied_(0.0),
        xySatisfied_(1.0e-6),
        xyBranchValue_(0.0),
        xColumn_(xColumn),
        yColumn_(yColumn),
        firstLambda_(-1),
        branchingStrategy_(0),
        boundType_(0),
        xRow_(-1),
        yRow_(-1),
        xyRow_(xyRow),
        convexity_(-1),
        numberExtraRows_(0),
        multiplier_(NULL),
        extraRow_(NULL),
        chosen_(-1)
{
    double columnLower[4];
    double columnUpper[4];
    double objective[4];
    double rowLower[3];
    double rowUpper[3];
    CoinBigIndex starts[5];
    int index[16];
    double element[16];
    int i;
    starts[0] = 0;
    // rows
    int numberRows = solver->getNumRows();
    // convexity
    rowLower[0] = 1.0;
    rowUpper[0] = 1.0;
    convexity_ = numberRows;
    starts[1] = 0;
    // x
    rowLower[1] = 0.0;
    rowUpper[1] = 0.0;
    index[0] = xColumn_;
    element[0] = -1.0;
    xRow_ = numberRows + 1;
    starts[2] = 1;
    int nAdd = 2;
    if (xColumn_ != yColumn_) {
        rowLower[2] = 0.0;
        rowUpper[2] = 0.0;
        index[1] = yColumn;
        element[1] = -1.0;
        nAdd = 3;
        yRow_ = numberRows + 2;
        starts[3] = 2;
    } else {
        yRow_ = -1;
        branchingStrategy_ = 1;
    }
    // may be objective
    assert (xyRow_ >= -1);
    solver->addRows(nAdd, starts, index, element, rowLower, rowUpper);
    int n = 0;
    // order is LxLy, LxUy, UxLy and UxUy
    firstLambda_ = solver->getNumCols();
    // bit sloppy as theoretically could be infeasible but otherwise need to do more work
    double xB[2];
    double yB[2];
    const double * lower = solver->getColLower();
    const double * upper = solver->getColUpper();
    xB[0] = lower[xColumn_];
    xB[1] = upper[xColumn_];
    yB[0] = lower[yColumn_];
    yB[1] = upper[yColumn_];
    if (xMeshSize_ != floor(xMeshSize_)) {
        // not integral
        xSatisfied_ = CoinMax(xSatisfied_, 0.51 * xMeshSize_);
        if (!yMeshSize_) {
            xySatisfied_ = CoinMax(xySatisfied_, xSatisfied_ * CoinMax(fabs(yB[0]), fabs(yB[1])));
        }
    }
    if (yMeshSize_ != floor(yMeshSize_)) {
        // not integral
        ySatisfied_ = CoinMax(ySatisfied_, 0.51 * yMeshSize_);
        if (!xMeshSize_) {
            xySatisfied_ = CoinMax(xySatisfied_, ySatisfied_ * CoinMax(fabs(xB[0]), fabs(xB[1])));
        }
    }
    // adjust
    double distance;
    double steps;
    if (xMeshSize_) {
        distance = xB[1] - xB[0];
        steps = floor ((distance + 0.5 * xMeshSize_) / xMeshSize_);
        distance = xB[0] + xMeshSize_ * steps;
        if (fabs(xB[1] - distance) > xSatisfied_) {
            printf("bad x mesh %g %g %g -> %g\n", xB[0], xMeshSize_, xB[1], distance);
            //double newValue = CoinMax(fabs(xB[1]-distance),xMeshSize_);
            //printf("xSatisfied increased to %g\n",newValue);
            //xSatisfied_ = newValue;
            //xB[1]=distance;
            //solver->setColUpper(xColumn_,distance);
        }
    }
    if (yMeshSize_) {
        distance = yB[1] - yB[0];
        steps = floor ((distance + 0.5 * yMeshSize_) / yMeshSize_);
        distance = yB[0] + yMeshSize_ * steps;
        if (fabs(yB[1] - distance) > ySatisfied_) {
            printf("bad y mesh %g %g %g -> %g\n", yB[0], yMeshSize_, yB[1], distance);
            //double newValue = CoinMax(fabs(yB[1]-distance),yMeshSize_);
            //printf("ySatisfied increased to %g\n",newValue);
            //ySatisfied_ = newValue;
            //yB[1]=distance;
            //solver->setColUpper(yColumn_,distance);
        }
    }
    for (i = 0; i < 4; i++) {
        double x = (i < 2) ? xB[0] : xB[1];
        double y = ((i & 1) == 0) ? yB[0] : yB[1];
        columnLower[i] = 0.0;
        columnUpper[i] = 2.0;
        objective[i] = 0.0;
        double value;
        // xy
        value = coefficient_ * x * y;
        if (xyRow_ >= 0) {
            if (fabs(value) < 1.0e-19)
                value = 1.0e-19;
            element[n] = value;
            index[n++] = xyRow_;
        } else {
            objective[i] = value;
        }
        // convexity
        value = 1.0;
        element[n] = value;
        index[n++] = 0 + numberRows;
        // x
        value = x;
        if (fabs(value) < 1.0e-19)
            value = 1.0e-19;
        element[n] = value;
        index[n++] = 1 + numberRows;
        if (xColumn_ != yColumn_) {
            // y
            value = y;
            if (fabs(value) < 1.0e-19)
                value = 1.0e-19;
            element[n] = value;
            index[n++] = 2 + numberRows;
        }
        starts[i+1] = n;
    }
    solver->addCols(4, starts, index, element, columnLower, columnUpper, objective);
    // At least one has to have a mesh
    if (!xMeshSize_ && (!yMeshSize_ || yRow_ < 0)) {
        printf("one of x and y must have a mesh size\n");
        abort();
    } else if (yRow_ >= 0) {
        if (!xMeshSize_)
            branchingStrategy_ = 2;
        else if (!yMeshSize_)
            branchingStrategy_ = 1;
    }
    // Now add constraints to link in x and or y to existing ones.
    bool xDone = false;
    bool yDone = false;
    // order is LxLy, LxUy, UxLy and UxUy
    for (i = numberExistingObjects - 1; i >= 0; i--) {
        const OsiObject * obj = objects[i];
        const OsiBiLinear * obj2 =
            dynamic_cast <const OsiBiLinear *>(obj) ;
        if (obj2) {
            if (xColumn_ == obj2->xColumn_ && !xDone) {
                // make sure y equal
                double rhs = 0.0;
                CoinBigIndex starts[2];
                int index[4];
                double element[4] = {1.0, 1.0, -1.0, -1.0};
                starts[0] = 0;
                starts[1] = 4;
                index[0] = firstLambda_ + 0;
                index[1] = firstLambda_ + 1;
                index[2] = obj2->firstLambda_ + 0;
                index[3] = obj2->firstLambda_ + 1;
                solver->addRows(1, starts, index, element, &rhs, &rhs);
                xDone = true;
            }
            if (yColumn_ == obj2->yColumn_ && yRow_ >= 0 && !yDone) {
                // make sure x equal
                double rhs = 0.0;
                CoinBigIndex starts[2];
                int index[4];
                double element[4] = {1.0, 1.0, -1.0, -1.0};
                starts[0] = 0;
                starts[1] = 4;
                index[0] = firstLambda_ + 0;
                index[1] = firstLambda_ + 2;
                index[2] = obj2->firstLambda_ + 0;
                index[3] = obj2->firstLambda_ + 2;
                solver->addRows(1, starts, index, element, &rhs, &rhs);
                yDone = true;
            }
        }
    }
}
// Set sizes and other stuff
void
OsiBiLinear::setMeshSizes(const OsiSolverInterface * solver, double x, double y)
{
    xMeshSize_ = x;
    yMeshSize_ = y;
    double xB[2];
    double yB[2];
    const double * lower = solver->getColLower();
    const double * upper = solver->getColUpper();
    xB[0] = lower[xColumn_];
    xB[1] = upper[xColumn_];
    yB[0] = lower[yColumn_];
    yB[1] = upper[yColumn_];
    if (xMeshSize_ != floor(xMeshSize_)) {
        // not integral
        xSatisfied_ = CoinMax(xSatisfied_, 0.51 * xMeshSize_);
        if (!yMeshSize_) {
            xySatisfied_ = CoinMax(xySatisfied_, xSatisfied_ * CoinMax(fabs(yB[0]), fabs(yB[1])));
        }
    }
    if (yMeshSize_ != floor(yMeshSize_)) {
        // not integral
        ySatisfied_ = CoinMax(ySatisfied_, 0.51 * yMeshSize_);
        if (!xMeshSize_) {
            xySatisfied_ = CoinMax(xySatisfied_, ySatisfied_ * CoinMax(fabs(xB[0]), fabs(xB[1])));
        }
    }
}
// Useful constructor
OsiBiLinear::OsiBiLinear (CoinModel * coinModel, int xColumn,
                          int yColumn, int xyRow, double coefficient,
                          double xMesh, double yMesh,
                          int numberExistingObjects, const OsiObject ** objects )
        : OsiObject2(),
        coefficient_(coefficient),
        xMeshSize_(xMesh),
        yMeshSize_(yMesh),
        xSatisfied_(1.0e-6),
        ySatisfied_(1.0e-6),
        xOtherSatisfied_(0.0),
        yOtherSatisfied_(0.0),
        xySatisfied_(1.0e-6),
        xyBranchValue_(0.0),
        xColumn_(xColumn),
        yColumn_(yColumn),
        firstLambda_(-1),
        branchingStrategy_(0),
        boundType_(0),
        xRow_(-1),
        yRow_(-1),
        xyRow_(xyRow),
        convexity_(-1),
        numberExtraRows_(0),
        multiplier_(NULL),
        extraRow_(NULL),
        chosen_(-1)
{
    double columnLower[4];
    double columnUpper[4];
    double objective[4];
    double rowLower[3];
    double rowUpper[3];
    CoinBigIndex starts[5];
    int index[16];
    double element[16];
    int i;
    starts[0] = 0;
    // rows
    int numberRows = coinModel->numberRows();
    // convexity
    rowLower[0] = 1.0;
    rowUpper[0] = 1.0;
    convexity_ = numberRows;
    starts[1] = 0;
    // x
    rowLower[1] = 0.0;
    rowUpper[1] = 0.0;
    index[0] = xColumn_;
    element[0] = -1.0;
    xRow_ = numberRows + 1;
    starts[2] = 1;
    int nAdd = 2;
    if (xColumn_ != yColumn_) {
        rowLower[2] = 0.0;
        rowUpper[2] = 0.0;
        index[1] = yColumn;
        element[1] = -1.0;
        nAdd = 3;
        yRow_ = numberRows + 2;
        starts[3] = 2;
    } else {
        yRow_ = -1;
        branchingStrategy_ = 1;
    }
    // may be objective
    assert (xyRow_ >= -1);
    for (i = 0; i < nAdd; i++) {
        CoinBigIndex iStart = starts[i];
        coinModel->addRow(starts[i+1] - iStart, index + iStart, element + iStart, rowLower[i], rowUpper[i]);
    }
    int n = 0;
    // order is LxLy, LxUy, UxLy and UxUy
    firstLambda_ = coinModel->numberColumns();
    // bit sloppy as theoretically could be infeasible but otherwise need to do more work
    double xB[2];
    double yB[2];
    const double * lower = coinModel->columnLowerArray();
    const double * upper = coinModel->columnUpperArray();
    xB[0] = lower[xColumn_];
    xB[1] = upper[xColumn_];
    yB[0] = lower[yColumn_];
    yB[1] = upper[yColumn_];
    if (xMeshSize_ != floor(xMeshSize_)) {
        // not integral
        xSatisfied_ = CoinMax(xSatisfied_, 0.51 * xMeshSize_);
        if (!yMeshSize_) {
            xySatisfied_ = CoinMax(xySatisfied_, xSatisfied_ * CoinMax(fabs(yB[0]), fabs(yB[1])));
        }
    }
    if (yMeshSize_ != floor(yMeshSize_)) {
        // not integral
        ySatisfied_ = CoinMax(ySatisfied_, 0.51 * yMeshSize_);
        if (!xMeshSize_) {
            xySatisfied_ = CoinMax(xySatisfied_, ySatisfied_ * CoinMax(fabs(xB[0]), fabs(xB[1])));
        }
    }
    // adjust
    double distance;
    double steps;
    if (xMeshSize_) {
        distance = xB[1] - xB[0];
        steps = floor ((distance + 0.5 * xMeshSize_) / xMeshSize_);
        distance = xB[0] + xMeshSize_ * steps;
        if (fabs(xB[1] - distance) > xSatisfied_) {
            printf("bad x mesh %g %g %g -> %g\n", xB[0], xMeshSize_, xB[1], distance);
            //double newValue = CoinMax(fabs(xB[1]-distance),xMeshSize_);
            //printf("xSatisfied increased to %g\n",newValue);
            //xSatisfied_ = newValue;
            //xB[1]=distance;
            //coinModel->setColUpper(xColumn_,distance);
        }
    }
    if (yMeshSize_) {
        distance = yB[1] - yB[0];
        steps = floor ((distance + 0.5 * yMeshSize_) / yMeshSize_);
        distance = yB[0] + yMeshSize_ * steps;
        if (fabs(yB[1] - distance) > ySatisfied_) {
            printf("bad y mesh %g %g %g -> %g\n", yB[0], yMeshSize_, yB[1], distance);
            //double newValue = CoinMax(fabs(yB[1]-distance),yMeshSize_);
            //printf("ySatisfied increased to %g\n",newValue);
            //ySatisfied_ = newValue;
            //yB[1]=distance;
            //coinModel->setColUpper(yColumn_,distance);
        }
    }
    for (i = 0; i < 4; i++) {
        double x = (i < 2) ? xB[0] : xB[1];
        double y = ((i & 1) == 0) ? yB[0] : yB[1];
        columnLower[i] = 0.0;
        columnUpper[i] = 2.0;
        objective[i] = 0.0;
        double value;
        // xy
        value = coefficient_ * x * y;
        if (xyRow_ >= 0) {
            if (fabs(value) < 1.0e-19)
                value = 1.0e-19;
            element[n] = value;
            index[n++] = xyRow_;
        } else {
            objective[i] = value;
        }
        // convexity
        value = 1.0;
        element[n] = value;
        index[n++] = 0 + numberRows;
        // x
        value = x;
        if (fabs(value) < 1.0e-19)
            value = 1.0e-19;
        element[n] = value;
        index[n++] = 1 + numberRows;
        if (xColumn_ != yColumn_) {
            // y
            value = y;
            if (fabs(value) < 1.0e-19)
                value = 1.0e-19;
            element[n] = value;
            index[n++] = 2 + numberRows;
        }
        starts[i+1] = n;
    }
    for (i = 0; i < 4; i++) {
        CoinBigIndex iStart = starts[i];
        coinModel->addColumn(starts[i+1] - iStart, index + iStart, element + iStart, columnLower[i],
                             columnUpper[i], objective[i]);
    }
    // At least one has to have a mesh
    if (!xMeshSize_ && (!yMeshSize_ || yRow_ < 0)) {
        printf("one of x and y must have a mesh size\n");
        abort();
    } else if (yRow_ >= 0) {
        if (!xMeshSize_)
            branchingStrategy_ = 2;
        else if (!yMeshSize_)
            branchingStrategy_ = 1;
    }
    // Now add constraints to link in x and or y to existing ones.
    bool xDone = false;
    bool yDone = false;
    // order is LxLy, LxUy, UxLy and UxUy
    for (i = numberExistingObjects - 1; i >= 0; i--) {
        const OsiObject * obj = objects[i];
        const OsiBiLinear * obj2 =
            dynamic_cast <const OsiBiLinear *>(obj) ;
        if (obj2) {
            if (xColumn_ == obj2->xColumn_ && !xDone) {
                // make sure y equal
                double rhs = 0.0;
                int index[4];
                double element[4] = {1.0, 1.0, -1.0, -1.0};
                index[0] = firstLambda_ + 0;
                index[1] = firstLambda_ + 1;
                index[2] = obj2->firstLambda_ + 0;
                index[3] = obj2->firstLambda_ + 1;
                coinModel->addRow(4, index, element, rhs, rhs);
                xDone = true;
            }
            if (yColumn_ == obj2->yColumn_ && yRow_ >= 0 && !yDone) {
                // make sure x equal
                double rhs = 0.0;
                int index[4];
                double element[4] = {1.0, 1.0, -1.0, -1.0};
                index[0] = firstLambda_ + 0;
                index[1] = firstLambda_ + 2;
                index[2] = obj2->firstLambda_ + 0;
                index[3] = obj2->firstLambda_ + 2;
                coinModel->addRow(4, index, element, rhs, rhs);
                yDone = true;
            }
        }
    }
}

// Copy constructor
OsiBiLinear::OsiBiLinear ( const OsiBiLinear & rhs)
        : OsiObject2(rhs),
        coefficient_(rhs.coefficient_),
        xMeshSize_(rhs.xMeshSize_),
        yMeshSize_(rhs.yMeshSize_),
        xSatisfied_(rhs.xSatisfied_),
        ySatisfied_(rhs.ySatisfied_),
        xOtherSatisfied_(rhs.xOtherSatisfied_),
        yOtherSatisfied_(rhs.yOtherSatisfied_),
        xySatisfied_(rhs.xySatisfied_),
        xyBranchValue_(rhs.xyBranchValue_),
        xColumn_(rhs.xColumn_),
        yColumn_(rhs.yColumn_),
        firstLambda_(rhs.firstLambda_),
        branchingStrategy_(rhs.branchingStrategy_),
        boundType_(rhs.boundType_),
        xRow_(rhs.xRow_),
        yRow_(rhs.yRow_),
        xyRow_(rhs.xyRow_),
        convexity_(rhs.convexity_),
        numberExtraRows_(rhs.numberExtraRows_),
        multiplier_(NULL),
        extraRow_(NULL),
        chosen_(rhs.chosen_)
{
    if (numberExtraRows_) {
        multiplier_ = CoinCopyOfArray(rhs.multiplier_, numberExtraRows_);
        extraRow_ = CoinCopyOfArray(rhs.extraRow_, numberExtraRows_);
    }
}

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

// Assignment operator
OsiBiLinear &
OsiBiLinear::operator=( const OsiBiLinear & rhs)
{
    if (this != &rhs) {
        OsiObject2::operator=(rhs);
        coefficient_ = rhs.coefficient_;
        xMeshSize_ = rhs.xMeshSize_;
        yMeshSize_ = rhs.yMeshSize_;
        xSatisfied_ = rhs.xSatisfied_;
        ySatisfied_ = rhs.ySatisfied_;
        xOtherSatisfied_ = rhs.xOtherSatisfied_;
        yOtherSatisfied_ = rhs.yOtherSatisfied_;
        xySatisfied_ = rhs.xySatisfied_;
        xyBranchValue_ = rhs.xyBranchValue_;
        xColumn_ = rhs.xColumn_;
        yColumn_ = rhs.yColumn_;
        firstLambda_ = rhs.firstLambda_;
        branchingStrategy_ = rhs.branchingStrategy_;
        boundType_ = rhs.boundType_;
        xRow_ = rhs.xRow_;
        yRow_ = rhs.yRow_;
        xyRow_ = rhs.xyRow_;
        convexity_ = rhs.convexity_;
        numberExtraRows_ = rhs.numberExtraRows_;
        delete [] multiplier_;
        delete [] extraRow_;
        if (numberExtraRows_) {
            multiplier_ = CoinCopyOfArray(rhs.multiplier_, numberExtraRows_);
            extraRow_ = CoinCopyOfArray(rhs.extraRow_, numberExtraRows_);
        } else {
            multiplier_ = NULL;
            extraRow_ = NULL;
        }
        chosen_ = rhs.chosen_;
    }
    return *this;
}

// Destructor
OsiBiLinear::~OsiBiLinear ()
{
    delete [] multiplier_;
    delete [] extraRow_;
}
// Adds in data for extra row with variable coefficients
void
OsiBiLinear::addExtraRow(int row, double multiplier)
{
    int * tempI = new int [numberExtraRows_+1];
    double * tempD = new double [numberExtraRows_+1];
    memcpy(tempI, extraRow_, numberExtraRows_*sizeof(int));
    memcpy(tempD, multiplier_, numberExtraRows_*sizeof(double));
    tempI[numberExtraRows_] = row;
    tempD[numberExtraRows_] = multiplier;
    if (numberExtraRows_)
        assert (row > tempI[numberExtraRows_-1]);
    numberExtraRows_++;
    delete [] extraRow_;
    extraRow_ = tempI;
    delete [] multiplier_;
    multiplier_ = tempD;
}
static bool testCoarse = true;
// Infeasibility - large is 0.5
double
OsiBiLinear::infeasibility(const OsiBranchingInformation * info, int & whichWay) const
{
    // order is LxLy, LxUy, UxLy and UxUy
    double xB[2];
    double yB[2];
    xB[0] = info->lower_[xColumn_];
    xB[1] = info->upper_[xColumn_];
    yB[0] = info->lower_[yColumn_];
    yB[1] = info->upper_[yColumn_];
#ifdef JJF_ZERO
    if (info->lower_[1] <= 43.0 && info->upper_[1] >= 43.0) {
        if (info->lower_[4] <= 49.0 && info->upper_[4] >= 49.0) {
            if (info->lower_[2] <= 16.0 && info->upper_[2] >= 16.0) {
                if (info->lower_[3] <= 19.0 && info->upper_[3] >= 19.0) {
                    printf("feas %g %g %g %g  p %g t %g\n",
                           info->solution_[1],
                           info->solution_[2],
                           info->solution_[3],
                           info->solution_[4],
                           info->solution_[0],
                           info->solution_[5]);
                }
            }
        }
    }
#endif
    double x = info->solution_[xColumn_];
    x = CoinMax(x, xB[0]);
    x = CoinMin(x, xB[1]);
    double y = info->solution_[yColumn_];
    y = CoinMax(y, yB[0]);
    y = CoinMin(y, yB[1]);
    int j;
#ifndef NDEBUG
    double xLambda = 0.0;
    double yLambda = 0.0;
    if ((branchingStrategy_&4) == 0) {
        for (j = 0; j < 4; j++) {
            int iX = j >> 1;
            int iY = j & 1;
            xLambda += xB[iX] * info->solution_[firstLambda_+j];
            if (yRow_ >= 0)
                yLambda += yB[iY] * info->solution_[firstLambda_+j];
        }
    } else {
        const double * element = info->elementByColumn_;
        const int * row = info->row_;
        const CoinBigIndex * columnStart = info->columnStart_;
        const int * columnLength = info->columnLength_;
        for (j = 0; j < 4; j++) {
            int iColumn = firstLambda_ + j;
            int iStart = columnStart[iColumn];
            int iEnd = iStart + columnLength[iColumn];
            int k = iStart;
            double sol = info->solution_[iColumn];
            for (; k < iEnd; k++) {
                if (xRow_ == row[k])
                    xLambda += element[k] * sol;
                if (yRow_ == row[k])
                    yLambda += element[k] * sol;
            }
        }
    }
    assert (fabs(x - xLambda) < 1.0e-1);
    if (yRow_ >= 0)
        assert (fabs(y - yLambda) < 1.0e-1);
#endif
    // If x or y not satisfied then branch on that
    double distance;
    double steps;
    bool xSatisfied;
    double xNew = xB[0];
    if (xMeshSize_) {
        if (x < 0.5*(xB[0] + xB[1])) {
            distance = x - xB[0];
            steps = floor ((distance + 0.5 * xMeshSize_) / xMeshSize_);
            xNew = xB[0] + steps * xMeshSize_;
            assert (xNew <= xB[1] + xSatisfied_);
            xSatisfied =  (fabs(xNew - x) < xSatisfied_);
        } else {
            distance = xB[1] - x;
            steps = floor ((distance + 0.5 * xMeshSize_) / xMeshSize_);
            xNew = xB[1] - steps * xMeshSize_;
            assert (xNew >= xB[0] - xSatisfied_);
            xSatisfied =  (fabs(xNew - x) < xSatisfied_);
        }
        // but if first coarse grid then only if gap small
        if (testCoarse && (branchingStrategy_&8) != 0 && xSatisfied &&
                xB[1] - xB[0] >= xMeshSize_) {
            // but allow if fine grid would allow
            if (fabs(xNew - x) >= xOtherSatisfied_ && fabs(yB[0] - y) > yOtherSatisfied_
                    && fabs(yB[1] - y) > yOtherSatisfied_) {
                xNew = 0.5 * (xB[0] + xB[1]);
                x = xNew;
                xSatisfied = false;
            }
        }
    } else {
        xSatisfied = true;
    }
    bool ySatisfied;
    double yNew = yB[0];
    if (yMeshSize_) {
        if (y < 0.5*(yB[0] + yB[1])) {
            distance = y - yB[0];
            steps = floor ((distance + 0.5 * yMeshSize_) / yMeshSize_);
            yNew = yB[0] + steps * yMeshSize_;
            assert (yNew <= yB[1] + ySatisfied_);
            ySatisfied =  (fabs(yNew - y) < ySatisfied_);
        } else {
            distance = yB[1] - y;
            steps = floor ((distance + 0.5 * yMeshSize_) / yMeshSize_);
            yNew = yB[1] - steps * yMeshSize_;
            assert (yNew >= yB[0] - ySatisfied_);
            ySatisfied =  (fabs(yNew - y) < ySatisfied_);
        }
        // but if first coarse grid then only if gap small
        if (testCoarse && (branchingStrategy_&8) != 0 && ySatisfied &&
                yB[1] - yB[0] >= yMeshSize_) {
            // but allow if fine grid would allow
            if (fabs(yNew - y) >= yOtherSatisfied_ && fabs(xB[0] - x) > xOtherSatisfied_
                    && fabs(xB[1] - x) > xOtherSatisfied_) {
                yNew = 0.5 * (yB[0] + yB[1]);
                y = yNew;
                ySatisfied = false;
            }
        }
    } else {
        ySatisfied = true;
    }
    /* There are several possibilities
       1 - one or both are unsatisfied and branching strategy tells us what to do
       2 - both are unsatisfied and branching strategy is 0
       3 - both are satisfied but xy is not
           3a one has bounds within satisfied_ - other does not
       (or neither have but branching strategy tells us what to do)
       3b neither do - and branching strategy does not tell us
       3c both do - treat as feasible knowing another copy of object will fix
       4 - both are satisfied and xy is satisfied - as 3c
    */
    chosen_ = -1;
    xyBranchValue_ = COIN_DBL_MAX;
    whichWay_ = 0;
    double xyTrue = x * y;
    double xyLambda = 0.0;
    if ((branchingStrategy_&4) == 0) {
        for (j = 0; j < 4; j++) {
            int iX = j >> 1;
            int iY = j & 1;
            xyLambda += xB[iX] * yB[iY] * info->solution_[firstLambda_+j];
        }
    } else {
        if (xyRow_ >= 0) {
            const double * element = info->elementByColumn_;
            const int * row = info->row_;
            const CoinBigIndex * columnStart = info->columnStart_;
            const int * columnLength = info->columnLength_;
            for (j = 0; j < 4; j++) {
                int iColumn = firstLambda_ + j;
                int iStart = columnStart[iColumn];
                int iEnd = iStart + columnLength[iColumn];
                int k = iStart;
                double sol = info->solution_[iColumn];
                for (; k < iEnd; k++) {
                    if (xyRow_ == row[k])
                        xyLambda += element[k] * sol;
                }
            }
        } else {
            // objective
            const double * objective = info->objective_;
            for (j = 0; j < 4; j++) {
                int iColumn = firstLambda_ + j;
                double sol = info->solution_[iColumn];
                xyLambda += objective[iColumn] * sol;
            }
        }
        xyLambda /= coefficient_;
    }
    if (0) {
        // only true with positive values
        // see if all convexification constraints OK with true
        assert (xyTrue + 1.0e-5 > xB[0]*y + yB[0]*x - xB[0]*yB[0]);
        assert (xyTrue + 1.0e-5 > xB[1]*y + yB[1]*x - xB[1]*yB[1]);
        assert (xyTrue - 1.0e-5 < xB[1]*y + yB[0]*x - xB[1]*yB[0]);
        assert (xyTrue - 1.0e-5 < xB[0]*y + yB[1]*x - xB[0]*yB[1]);
        // see if all convexification constraints OK with lambda version
#ifndef JJF_ONE
        assert (xyLambda + 1.0e-5 > xB[0]*y + yB[0]*x - xB[0]*yB[0]);
        assert (xyLambda + 1.0e-5 > xB[1]*y + yB[1]*x - xB[1]*yB[1]);
        assert (xyLambda - 1.0e-5 < xB[1]*y + yB[0]*x - xB[1]*yB[0]);
        assert (xyLambda - 1.0e-5 < xB[0]*y + yB[1]*x - xB[0]*yB[1]);
#endif
        // see if other bound stuff true
        assert (xyLambda + 1.0e-5 > xB[0]*y);
        assert (xyLambda + 1.0e-5 > yB[0]*x);
        assert (xyLambda - 1.0e-5 < xB[1]*y);
        assert (xyLambda - 1.0e-5 < yB[1]*x);
#define SIZE 2
        if (yColumn_ == xColumn_ + SIZE) {
#if SIZE==6
            double bMax = 2200.0;
            double bMin = bMax - 100.0;
            double b[] = {330.0, 360.0, 380.0, 430.0, 490.0, 530.0};
#elif SIZE==2
            double bMax = 1900.0;
            double bMin = bMax - 200.0;
            double b[] = {460.0, 570.0};
#else
            abort();
#endif
            double sum = 0.0;
            double sum2 = 0.0;
            int m = xColumn_;
            double x = info->solution_[m];
            double xB[2];
            double yB[2];
            xB[0] = info->lower_[m];
            xB[1] = info->upper_[m];
            for (int i = 0; i < SIZE*SIZE; i += SIZE) {
                int n = i + SIZE + m;
                double y = info->solution_[n];
                yB[0] = info->lower_[n];
                yB[1] = info->upper_[n];
                int firstLambda = SIZE * SIZE + 2 * SIZE + 4 * i + 4 * m;
                double xyLambda = 0.0;
                for (int j = 0; j < 4; j++) {
                    int iX = j >> 1;
                    int iY = j & 1;
                    xyLambda += xB[iX] * yB[iY] * info->solution_[firstLambda+j];
                }
                sum += xyLambda * b[i/SIZE];
                double xyTrue = x * y;
                sum2 += xyTrue * b[i/SIZE];
            }
            if (sum > bMax*x + 1.0e-5 || sum < bMin*x - 1.0e-5) {
                //if (sum<bMax*x+1.0e-5&&sum>bMin*x-1.0e-5) {
                printf("bmin*x %g b*w %g bmax*x %g (true) %g\n", bMin*x, sum, bMax*x, sum2);
                printf("m %d lb %g value %g up %g\n",
                       m, xB[0], x, xB[1]);
                sum = 0.0;
                for (int i = 0; i < SIZE*SIZE; i += SIZE) {
                    int n = i + SIZE + m;
                    double y = info->solution_[n];
                    yB[0] = info->lower_[n];
                    yB[1] = info->upper_[n];
                    printf("n %d lb %g value %g up %g\n",
                           n, yB[0], y, yB[1]);
                    int firstLambda = SIZE * SIZE + 2 * SIZE + 4 * i + m * 4;
                    double xyLambda = 0.0;
                    for (int j = 0; j < 4; j++) {
                        int iX = j >> 1;
                        int iY = j & 1;
                        xyLambda += xB[iX] * yB[iY] * info->solution_[firstLambda+j];
                        printf("j %d l %d new xylambda %g ", j, firstLambda + j, xyLambda);
                    }
                    sum += xyLambda * b[i/SIZE];
                    printf(" - sum now %g\n", sum);
                }
            }
            if (sum2 > bMax*x + 1.0e-5 || sum2 < bMin*x - 1.0e-5) {
                printf("bmin*x %g b*x*y %g bmax*x %g (estimate) %g\n", bMin*x, sum2, bMax*x, sum);
                printf("m %d lb %g value %g up %g\n",
                       m, xB[0], x, xB[1]);
                sum2 = 0.0;
                for (int i = 0; i < SIZE*SIZE; i += SIZE) {
                    int n = i + SIZE + m;
                    double y = info->solution_[n];
                    yB[0] = info->lower_[n];
                    yB[1] = info->upper_[n];
                    printf("n %d lb %g value %g up %g\n",
                           n, yB[0], y, yB[1]);
                    double xyTrue = x * y;
                    sum2 += xyTrue * b[i/SIZE];
                    printf("xyTrue %g - sum now %g\n", xyTrue, sum2);
                }
            }
        }
    }
    // If pseudo shadow prices then see what would happen
    //double pseudoEstimate = 0.0;
    if (info->defaultDual_ >= 0.0) {
        // If we move to xy then we move by coefficient * (xyTrue-xyLambda) on row xyRow_
        double move = xyTrue - xyLambda;
        assert (xyRow_ >= 0);
        if (boundType_ == 0) {
            move *= coefficient_;
            move *= info->pi_[xyRow_];
            move = CoinMax(move, 0.0);
        } else if (boundType_ == 1) {
            // if OK then say satisfied
        } else if (boundType_ == 2) {
        } else {
            // == row so move x and y not xy
        }
    }
    if ((branchingStrategy_&16) != 0) {
        // always treat as satisfied!!
        xSatisfied = true;
        ySatisfied = true;
        xyTrue = xyLambda;
    }
    if ( !xSatisfied) {
        if (!ySatisfied) {
            if ((branchingStrategy_&3) == 0) {
                // If pseudo shadow prices then see what would happen
                if (info->defaultDual_ >= 0.0) {
                    // need coding here
                    if (fabs(x - xNew) > fabs(y - yNew)) {
                        chosen_ = 0;
                        xyBranchValue_ = x;
                    } else {
                        chosen_ = 1;
                        xyBranchValue_ = y;
                    }
                } else {
                    if (fabs(x - xNew) > fabs(y - yNew)) {
                        chosen_ = 0;
                        xyBranchValue_ = x;
                    } else {
                        chosen_ = 1;
                        xyBranchValue_ = y;
                    }
                }
            } else if ((branchingStrategy_&3) == 1) {
                chosen_ = 0;
                xyBranchValue_ = x;
            } else {
                chosen_ = 1;
                xyBranchValue_ = y;
            }
        } else {
            // y satisfied
            chosen_ = 0;
            xyBranchValue_ = x;
        }
    } else {
        // x satisfied
        if (!ySatisfied) {
            chosen_ = 1;
            xyBranchValue_ = y;
        } else {
            /*
            3 - both are satisfied but xy is not
            3a one has bounds within satisfied_ - other does not
            (or neither have but branching strategy tells us what to do)
            3b neither do - and branching strategy does not tell us
            3c both do - treat as feasible knowing another copy of object will fix
              4 - both are satisfied and xy is satisfied - as 3c
            */
            if (fabs(xyLambda - xyTrue) < xySatisfied_ || (xB[0] == xB[1] && yB[0] == yB[1])) {
                // satisfied
#ifdef JJF_ZERO
                printf("all satisfied true %g lambda %g\n",
                       xyTrue, xyLambda);
                printf("x %d (%g,%g,%g) y %d (%g,%g,%g)\n",
                       xColumn_, xB[0], x, xB[1],
                       yColumn_, yB[0], y, yB[1]);
#endif
            } else {
                // May be infeasible - check
                bool feasible = true;
                if (xB[0] == xB[1] && yB[0] == yB[1]) {
                    double lambda[4];
                    computeLambdas(info->solver_, lambda);
                    for (int j = 0; j < 4; j++) {
                        int iColumn = firstLambda_ + j;
                        if (info->lower_[iColumn] > lambda[j] + 1.0e-5 ||
                                info->upper_[iColumn] < lambda[j] - 1.0e-5)
                            feasible = false;
                    }
                }
                if (testCoarse && (branchingStrategy_&8) != 0 && xB[1] - xB[0] < 1.0001*xSatisfied_ &&
                        yB[1] - yB[0] < 1.0001*ySatisfied_)
                    feasible = true;
                if (feasible) {
                    if (xB[1] - xB[0] >= xSatisfied_ && xMeshSize_) {
                        if (yB[1] - yB[0] >= ySatisfied_ && yMeshSize_) {
                            if ((branchingStrategy_&3) == 0) {
                                // If pseudo shadow prices then see what would happen
                                if (info->defaultDual_ >= 0.0) {
                                    // need coding here
                                    if (xB[1] - xB[0] > yB[1] - yB[0]) {
                                        chosen_ = 0;
                                        xyBranchValue_ = 0.5 * (xB[0] + xB[1]);
                                    } else {
                                        chosen_ = 1;
                                        xyBranchValue_ = 0.5 * (yB[0] + yB[1]);
                                    }
                                } else {
                                    if (xB[1] - xB[0] > yB[1] - yB[0]) {
                                        chosen_ = 0;
                                        xyBranchValue_ = 0.5 * (xB[0] + xB[1]);
                                    } else {
                                        chosen_ = 1;
                                        xyBranchValue_ = 0.5 * (yB[0] + yB[1]);
                                    }
                                }
                            } else if ((branchingStrategy_&3) == 1) {
                                chosen_ = 0;
                                xyBranchValue_ = 0.5 * (xB[0] + xB[1]);
                            } else {
                                chosen_ = 1;
                                xyBranchValue_ = 0.5 * (yB[0] + yB[1]);
                            }
                        } else {
                            // y satisfied
                            chosen_ = 0;
                            xyBranchValue_ = 0.5 * (xB[0] + xB[1]);
                        }
                    } else if (yB[1] - yB[0] >= ySatisfied_ && yMeshSize_) {
                        chosen_ = 1;
                        xyBranchValue_ = 0.5 * (yB[0] + yB[1]);
                    } else {
                        // treat as satisfied unless no coefficient tightening
                        if ((branchingStrategy_&4) != 0) {
                            chosen_ = 0; // fix up in branch
                            xyBranchValue_ = x;
                        }
                    }
                } else {
                    // node not feasible!!!
                    chosen_ = 0;
                    infeasibility_ = COIN_DBL_MAX;
                    otherInfeasibility_ = COIN_DBL_MAX;
                    whichWay = whichWay_;
                    return infeasibility_;
                }
            }
        }
    }
    if (chosen_ == -1) {
        infeasibility_ = 0.0;
    } else if (chosen_ == 0) {
        infeasibility_ = CoinMax(fabs(xyBranchValue_ - x), 1.0e-12);
        //assert (xyBranchValue_>=info->lower_[xColumn_]&&xyBranchValue_<=info->upper_[xColumn_]);
    } else {
        infeasibility_ = CoinMax(fabs(xyBranchValue_ - y), 1.0e-12);
        //assert (xyBranchValue_>=info->lower_[yColumn_]&&xyBranchValue_<=info->upper_[yColumn_]);
    }
    if (info->defaultDual_ < 0.0) {
        // not using pseudo shadow prices
        otherInfeasibility_ = 1.0 - infeasibility_;
    } else {
        abort();
    }
    if (infeasibility_) {
        bool fixed = true;
        for (int j = 0; j < 4; j++) {
            int iColumn = firstLambda_ + j;
            //if (info->lower_[iColumn]) printf("lower of %g on %d\n",info->lower_[iColumn],iColumn);
            if (info->lower_[iColumn] < info->upper_[iColumn])
                fixed = false;
        }
        if (fixed) {
            //printf("must be tolerance problem - xy true %g lambda %g\n",xyTrue,xyLambda);
            chosen_ = -1;
            infeasibility_ = 0.0;
        }
    }
    whichWay = whichWay_;
    //if (infeasibility_&&priority_==10)
    //printf("x %d %g %g %g, y %d %g %g %g\n",xColumn_,xB[0],x,xB[1],yColumn_,yB[0],y,yB[1]);
    return infeasibility_;
}
// Sets infeasibility and other when pseudo shadow prices
void
OsiBiLinear::getPseudoShadow(const OsiBranchingInformation * info)
{
    // order is LxLy, LxUy, UxLy and UxUy
    double xB[2];
    double yB[2];
    xB[0] = info->lower_[xColumn_];
    xB[1] = info->upper_[xColumn_];
    yB[0] = info->lower_[yColumn_];
    yB[1] = info->upper_[yColumn_];
    double x = info->solution_[xColumn_];
    x = CoinMax(x, xB[0]);
    x = CoinMin(x, xB[1]);
    double y = info->solution_[yColumn_];
    y = CoinMax(y, yB[0]);
    y = CoinMin(y, yB[1]);
    int j;
    double xyTrue = x * y;
    double xyLambda = 0.0;
    if ((branchingStrategy_&4) == 0) {
        for (j = 0; j < 4; j++) {
            int iX = j >> 1;
            int iY = j & 1;
            xyLambda += xB[iX] * yB[iY] * info->solution_[firstLambda_+j];
        }
    } else {
        if (xyRow_ >= 0) {
            const double * element = info->elementByColumn_;
            const int * row = info->row_;
            const CoinBigIndex * columnStart = info->columnStart_;
            const int * columnLength = info->columnLength_;
            for (j = 0; j < 4; j++) {
                int iColumn = firstLambda_ + j;
                int iStart = columnStart[iColumn];
                int iEnd = iStart + columnLength[iColumn];
                int k = iStart;
                double sol = info->solution_[iColumn];
                for (; k < iEnd; k++) {
                    if (xyRow_ == row[k])
                        xyLambda += element[k] * sol;
                }
            }
        } else {
            // objective
            const double * objective = info->objective_;
            for (j = 0; j < 4; j++) {
                int iColumn = firstLambda_ + j;
                double sol = info->solution_[iColumn];
                xyLambda += objective[iColumn] * sol;
            }
        }
        xyLambda /= coefficient_;
    }
    assert (info->defaultDual_ >= 0.0);
    // If we move to xy then we move by coefficient * (xyTrue-xyLambda) on row xyRow_
    double movement = xyTrue - xyLambda;
    infeasibility_ = 0.0;
    const double * pi = info->pi_;
    const double * activity = info->rowActivity_;
    const double * lower = info->rowLower_;
    const double * upper = info->rowUpper_;
    double tolerance = info->primalTolerance_;
    double direction = info->direction_;
    bool infeasible = false;
    if (xyRow_ >= 0) {
        assert (!boundType_);
        if (lower[xyRow_] < -1.0e20)
            assert (pi[xyRow_] <= 1.0e-3);
        if (upper[xyRow_] > 1.0e20)
            assert (pi[xyRow_] >= -1.0e-3);
        double valueP = pi[xyRow_] * direction;
        // if move makes infeasible then make at least default
        double newValue = activity[xyRow_] + movement * coefficient_;
        if (newValue > upper[xyRow_] + tolerance || newValue < lower[xyRow_] - tolerance) {
            infeasibility_ += fabs(movement * coefficient_) * CoinMax(fabs(valueP), info->defaultDual_);
            infeasible = true;
        }
    } else {
        // objective
        assert (movement > -1.0e-7);
        infeasibility_ += movement;
    }
    for (int i = 0; i < numberExtraRows_; i++) {
        int iRow = extraRow_[i];
        if (lower[iRow] < -1.0e20)
            assert (pi[iRow] <= 1.0e-3);
        if (upper[iRow] > 1.0e20)
            assert (pi[iRow] >= -1.0e-3);
        double valueP = pi[iRow] * direction;
        // if move makes infeasible then make at least default
        double newValue = activity[iRow] + movement * multiplier_[i];
        if (newValue > upper[iRow] + tolerance || newValue < lower[iRow] - tolerance) {
            infeasibility_ += fabs(movement * multiplier_[i]) * CoinMax(fabs(valueP), info->defaultDual_);
            infeasible = true;
        }
    }
    if (infeasibility_ < info->integerTolerance_) {
        if (!infeasible)
            infeasibility_ = 0.0;
        else
            infeasibility_ = info->integerTolerance_;
    }
    otherInfeasibility_ = CoinMax(1.0e-12, infeasibility_ * 10.0);
}
// Gets sum of movements to correct value
double
OsiBiLinear::getMovement(const OsiBranchingInformation * info)
{
    // order is LxLy, LxUy, UxLy and UxUy
    double xB[2];
    double yB[2];
    xB[0] = info->lower_[xColumn_];
    xB[1] = info->upper_[xColumn_];
    yB[0] = info->lower_[yColumn_];
    yB[1] = info->upper_[yColumn_];
    double x = info->solution_[xColumn_];
    x = CoinMax(x, xB[0]);
    x = CoinMin(x, xB[1]);
    double y = info->solution_[yColumn_];
    y = CoinMax(y, yB[0]);
    y = CoinMin(y, yB[1]);
    int j;
    double xyTrue = x * y;
    double xyLambda = 0.0;
    if ((branchingStrategy_&4) == 0) {
        for (j = 0; j < 4; j++) {
            int iX = j >> 1;
            int iY = j & 1;
            xyLambda += xB[iX] * yB[iY] * info->solution_[firstLambda_+j];
        }
    } else {
        if (xyRow_ >= 0) {
            const double * element = info->elementByColumn_;
            const int * row = info->row_;
            const CoinBigIndex * columnStart = info->columnStart_;
            const int * columnLength = info->columnLength_;
            for (j = 0; j < 4; j++) {
                int iColumn = firstLambda_ + j;
                int iStart = columnStart[iColumn];
                int iEnd = iStart + columnLength[iColumn];
                int k = iStart;
                double sol = info->solution_[iColumn];
                for (; k < iEnd; k++) {
                    if (xyRow_ == row[k])
                        xyLambda += element[k] * sol;
                }
            }
        } else {
            // objective
            const double * objective = info->objective_;
            for (j = 0; j < 4; j++) {
                int iColumn = firstLambda_ + j;
                double sol = info->solution_[iColumn];
                xyLambda += objective[iColumn] * sol;
            }
        }
        xyLambda /= coefficient_;
    }
    // If we move to xy then we move by coefficient * (xyTrue-xyLambda) on row xyRow_
    double movement = xyTrue - xyLambda;
    double mesh = CoinMax(xMeshSize_, yMeshSize_);
    if (fabs(movement) < xySatisfied_ && (xB[1] - xB[0] < mesh || yB[1] - yB[0] < mesh))
        return 0.0; // say feasible
    const double * activity = info->rowActivity_;
    const double * lower = info->rowLower_;
    const double * upper = info->rowUpper_;
    double tolerance = info->primalTolerance_;
    double  infeasibility = 0.0;
    if (xyRow_ >= 0) {
        assert (!boundType_);
        // if move makes infeasible
        double newValue = activity[xyRow_] + movement * coefficient_;
        if (newValue > upper[xyRow_] + tolerance)
            infeasibility += newValue - upper[xyRow_];
        else if (newValue < lower[xyRow_] - tolerance)
            infeasibility += lower[xyRow_] - newValue;
    } else {
        // objective
        assert (movement > -1.0e-7);
        infeasibility += movement;
    }
    for (int i = 0; i < numberExtraRows_; i++) {
        int iRow = extraRow_[i];
        // if move makes infeasible
        double newValue = activity[iRow] + movement * multiplier_[i];
        if (newValue > upper[iRow] + tolerance)
            infeasibility += newValue - upper[iRow];
        else if (newValue < lower[iRow] - tolerance)
            infeasibility += lower[iRow] - newValue;
    }
    return infeasibility;
}
// This looks at solution and sets bounds to contain solution
double
OsiBiLinear::feasibleRegion(OsiSolverInterface * solver, const OsiBranchingInformation * info) const
{
    // If another object has finer mesh ignore this
    if ((branchingStrategy_&8) != 0)
        return 0.0;
    // order is LxLy, LxUy, UxLy and UxUy
    double xB[2];
    double yB[2];
    xB[0] = info->lower_[xColumn_];
    xB[1] = info->upper_[xColumn_];
    yB[0] = info->lower_[yColumn_];
    yB[1] = info->upper_[yColumn_];
    double x = info->solution_[xColumn_];
    double y = info->solution_[yColumn_];
    int j;
#ifndef NDEBUG
    double xLambda = 0.0;
    double yLambda = 0.0;
    if ((branchingStrategy_&4) == 0) {
        for (j = 0; j < 4; j++) {
            int iX = j >> 1;
            int iY = j & 1;
            xLambda += xB[iX] * info->solution_[firstLambda_+j];
            if (yRow_ >= 0)
                yLambda += yB[iY] * info->solution_[firstLambda_+j];
        }
    } else {
        const double * element = info->elementByColumn_;
        const int * row = info->row_;
        const CoinBigIndex * columnStart = info->columnStart_;
        const int * columnLength = info->columnLength_;
        for (j = 0; j < 4; j++) {
            int iColumn = firstLambda_ + j;
            int iStart = columnStart[iColumn];
            int iEnd = iStart + columnLength[iColumn];
            int k = iStart;
            double sol = info->solution_[iColumn];
            for (; k < iEnd; k++) {
                if (xRow_ == row[k])
                    xLambda += element[k] * sol;
                if (yRow_ == row[k])
                    yLambda += element[k] * sol;
            }
        }
    }
    if (yRow_ < 0)
        yLambda = xLambda;
#ifdef JJF_ZERO
    if (fabs(x - xLambda) > 1.0e-4 ||
            fabs(y - yLambda) > 1.0e-4)
        printf("feasibleregion x %d given %g lambda %g y %d given %g lambda %g\n",
               xColumn_, x, xLambda,
               yColumn_, y, yLambda);
#endif
#endif
    double infeasibility = 0.0;
    double distance;
    double steps;
    double xNew = x;
    if (xMeshSize_) {
        distance = x - xB[0];
        if (x < 0.5*(xB[0] + xB[1])) {
            distance = x - xB[0];
            steps = floor ((distance + 0.5 * xMeshSize_) / xMeshSize_);
            xNew = xB[0] + steps * xMeshSize_;
            assert (xNew <= xB[1] + xSatisfied_);
        } else {
            distance = xB[1] - x;
            steps = floor ((distance + 0.5 * xMeshSize_) / xMeshSize_);
            xNew = xB[1] - steps * xMeshSize_;
            assert (xNew >= xB[0] - xSatisfied_);
        }
        if (xMeshSize_ < 1.0 && fabs(xNew - x) <= xSatisfied_) {
            double lo = CoinMax(xB[0], x - 0.5 * xSatisfied_);
            double up = CoinMin(xB[1], x + 0.5 * xSatisfied_);
            solver->setColLower(xColumn_, lo);
            solver->setColUpper(xColumn_, up);
        } else {
            infeasibility +=  fabs(xNew - x);
            solver->setColLower(xColumn_, xNew);
            solver->setColUpper(xColumn_, xNew);
        }
    }
    double yNew = y;
    if (yMeshSize_) {
        distance = y - yB[0];
        if (y < 0.5*(yB[0] + yB[1])) {
            distance = y - yB[0];
            steps = floor ((distance + 0.5 * yMeshSize_) / yMeshSize_);
            yNew = yB[0] + steps * yMeshSize_;
            assert (yNew <= yB[1] + ySatisfied_);
        } else {
            distance = yB[1] - y;
            steps = floor ((distance + 0.5 * yMeshSize_) / yMeshSize_);
            yNew = yB[1] - steps * yMeshSize_;
            assert (yNew >= yB[0] - ySatisfied_);
        }
        if (yMeshSize_ < 1.0 && fabs(yNew - y) <= ySatisfied_) {
            double lo = CoinMax(yB[0], y - 0.5 * ySatisfied_);
            double up = CoinMin(yB[1], y + 0.5 * ySatisfied_);
            solver->setColLower(yColumn_, lo);
            solver->setColUpper(yColumn_, up);
        } else {
            infeasibility +=  fabs(yNew - y);
            solver->setColLower(yColumn_, yNew);
            solver->setColUpper(yColumn_, yNew);
        }
    }
    if (0) {
        // temp
        solver->setColLower(xColumn_, x);
        solver->setColUpper(xColumn_, x);
        solver->setColLower(yColumn_, y);
        solver->setColUpper(yColumn_, y);
    }
    if ((branchingStrategy_&4)) {
        // fake to make correct
        double lambda[4];
        computeLambdas(solver, lambda);
        for (int j = 0; j < 4; j++) {
            int iColumn = firstLambda_ + j;
            double value = lambda[j];
            solver->setColLower(iColumn, value);
            solver->setColUpper(iColumn, value);
        }
    }
    double xyTrue = xNew * yNew;
    double xyLambda = 0.0;
    for (j = 0; j < 4; j++) {
        int iX = j >> 1;
        int iY = j & 1;
        xyLambda += xB[iX] * yB[iY] * info->solution_[firstLambda_+j];
    }
    infeasibility += fabs(xyTrue - xyLambda);
    return infeasibility;
}
// Returns true value of single xyRow coefficient
double
OsiBiLinear::xyCoefficient(const double * solution) const
{
    // If another object has finer mesh ignore this
    if ((branchingStrategy_&8) != 0)
        return 0.0;
    double x = solution[xColumn_];
    double y = solution[yColumn_];
    //printf("x (%d,%g) y (%d,%g) x*y*coefficient %g\n",
    // xColumn_,x,yColumn_,y,x*y*coefficient_);
    return x*y*coefficient_;
}

// Redoes data when sequence numbers change
void
OsiBiLinear::resetSequenceEtc(int numberColumns, const int * originalColumns)
{
    int i;
#ifdef JJF_ZERO
    for (i = 0; i < numberColumns; i++) {
        if (originalColumns[i] == firstLambda_)
            break;
    }
#else
    i = originalColumns[firstLambda_];
#endif
    if (i >= 0 && i < numberColumns) {
        firstLambda_ = i;
        for (int j = 0; j < 4; j++) {
            assert (originalColumns[j+i] - firstLambda_ == j);
        }
    } else {
        printf("lost set\n");
        abort();
    }
    // rows will be out anyway
    abort();
}

// Creates a branching object
OsiBranchingObject *
OsiBiLinear::createBranch(OsiSolverInterface * solver, const OsiBranchingInformation * /*info*/, int way) const
{
    // create object
    OsiBranchingObject * branch;
    assert (chosen_ == 0 || chosen_ == 1);
    //if (chosen_==0)
    //assert (xyBranchValue_>=info->lower_[xColumn_]&&xyBranchValue_<=info->upper_[xColumn_]);
    //else
    //assert (xyBranchValue_>=info->lower_[yColumn_]&&xyBranchValue_<=info->upper_[yColumn_]);
    branch = new OsiBiLinearBranchingObject(solver, this, way, xyBranchValue_, chosen_);
    return branch;
}
// Does work of branching
void
OsiBiLinear::newBounds(OsiSolverInterface * solver, int way, short xOrY, double separator) const
{
    int iColumn;
    double mesh;
    double satisfied;
    if (xOrY == 0) {
        iColumn = xColumn_;
        mesh = xMeshSize_;
        satisfied = xSatisfied_;
    } else {
        iColumn = yColumn_;
        mesh = yMeshSize_;
        satisfied = ySatisfied_;
    }
    const double * columnLower = solver->getColLower();
    const double * columnUpper = solver->getColUpper();
    double lower = columnLower[iColumn];
    double distance;
    double steps;
    double zNew = separator;
    distance = separator - lower;
    assert (mesh);
    steps = floor ((distance + 0.5 * mesh) / mesh);
    if (mesh < 1.0)
        zNew = lower + steps * mesh;
    if (zNew > columnUpper[iColumn] - satisfied)
        zNew = 0.5 * (columnUpper[iColumn] - lower);
    double oldUpper = columnUpper[iColumn] ;
    double oldLower = columnLower[iColumn] ;
#ifndef NDEBUG
    int nullChange = 0;
#endif
    if (way < 0) {
        if (zNew > separator && mesh < 1.0)
            zNew -= mesh;
        double oldUpper = columnUpper[iColumn] ;
        if (zNew + satisfied >= oldUpper)
            zNew = 0.5 * (oldUpper + oldLower);
        if (mesh == 1.0)
            zNew = floor(separator);
#ifndef NDEBUG
        if (oldUpper < zNew + 1.0e-8)
            nullChange = -1;
#endif
        solver->setColUpper(iColumn, zNew);
    } else {
        if (zNew < separator && mesh < 1.0)
            zNew += mesh;
        if (zNew - satisfied <= oldLower)
            zNew = 0.5 * (oldUpper + oldLower);
        if (mesh == 1.0)
            zNew = ceil(separator);
#ifndef NDEBUG
        if (oldLower > zNew - 1.0e-8)
            nullChange = 1;
#endif
        solver->setColLower(iColumn, zNew);
    }
    if ((branchingStrategy_&4) != 0
            && columnLower[xColumn_] == columnUpper[xColumn_]
            && columnLower[yColumn_] == columnUpper[yColumn_]) {
        // fake to make correct
        double lambda[4];
        computeLambdas(solver, lambda);
        for (int j = 0; j < 4; j++) {
            int iColumn = firstLambda_ + j;
            double value = lambda[j];
#ifndef NDEBUG
            if (fabs(value - columnLower[iColumn]) > 1.0e-5 ||
                    fabs(value - columnUpper[iColumn]) > 1.0e-5)
                nullChange = 0;
#endif
            solver->setColLower(iColumn, value);
            solver->setColUpper(iColumn, value);
        }
    }
#ifndef NDEBUG
    if (nullChange)
        printf("null change on column%s %d - bounds %g,%g\n", nullChange > 0 ? "Lower" : "Upper",
               iColumn, oldLower, oldUpper);
#endif
#ifdef JJF_ZERO
    // always free up lambda
    for (int i = firstLambda_; i < firstLambda_ + 4; i++) {
        solver->setColLower(i, 0.0);
        solver->setColUpper(i, 2.0);
    }
#endif
    double xB[3];
    xB[0] = columnLower[xColumn_];
    xB[1] = columnUpper[xColumn_];
    double yB[3];
    yB[0] = columnLower[yColumn_];
    yB[1] = columnUpper[yColumn_];
    if (false && (branchingStrategy_&4) != 0 && yRow_ >= 0 &&
            xMeshSize_ == 1.0 && yMeshSize_ == 1.0) {
        if ((xB[1] - xB[0])*(yB[1] - yB[0]) < 40) {
            // try looking at all solutions
            double lower[4];
            double upper[4];
            double lambda[4];
            int i;
            double lowerLambda[4];
            double upperLambda[4];
            for (i = 0; i < 4; i++) {
                lower[i] = CoinMax(0.0, columnLower[firstLambda_+i]);
                upper[i] = CoinMin(1.0, columnUpper[firstLambda_+i]);
                lowerLambda[i] = 1.0;
                upperLambda[i] = 0.0;
            }
            // get coefficients
            double xybar[4];
            getCoefficients(solver, xB, yB, xybar);
            double x, y;
            for (x = xB[0]; x <= xB[1]; x++) {
                xB[2] = x;
                for (y = yB[0]; y <= yB[1]; y++) {
                    yB[2] = y;
                    computeLambdas(xB, yB, xybar, lambda);
                    for (i = 0; i < 4; i++) {
                        lowerLambda[i] = CoinMin(lowerLambda[i], lambda[i]);
                        upperLambda[i] = CoinMax(upperLambda[i], lambda[i]);
                    }
                }
            }
            double change = 0.0;;
            for (i = 0; i < 4; i++) {
                if (lowerLambda[i] > lower[i] + 1.0e-12) {
                    solver->setColLower(firstLambda_ + i, lowerLambda[i]);
                    change += lowerLambda[i] - lower[i];
                }
                if (upperLambda[i] < upper[i] - 1.0e-12) {
                    solver->setColUpper(firstLambda_ + i, upperLambda[i]);
                    change -= upperLambda[i] - upper[i];
                }
            }
            if (change > 1.0e-5)
                printf("change of %g\n", change);
        }
    }
    if (boundType_) {
        assert (!xMeshSize_ || !yMeshSize_);
        if (xMeshSize_) {
            // can tighten bounds on y
            if ((boundType_&1) != 0) {
                if (xB[0]*yB[1] > coefficient_) {
                    // tighten upper bound on y
                    solver->setColUpper(yColumn_, coefficient_ / xB[0]);
                }
            }
            if ((boundType_&2) != 0) {
                if (xB[1]*yB[0] < coefficient_) {
                    // tighten lower bound on y
                    solver->setColLower(yColumn_, coefficient_ / xB[1]);
                }
            }
        } else {
            // can tighten bounds on x
            if ((boundType_&1) != 0) {
                if (yB[0]*xB[1] > coefficient_) {
                    // tighten upper bound on x
                    solver->setColUpper(xColumn_, coefficient_ / yB[0]);
                }
            }
            if ((boundType_&2) != 0) {
                if (yB[1]*xB[0] < coefficient_) {
                    // tighten lower bound on x
                    solver->setColLower(xColumn_, coefficient_ / yB[1]);
                }
            }
        }
    }
}
// Compute lambdas if coefficients not changing
void
OsiBiLinear::computeLambdas(const OsiSolverInterface * solver, double lambda[4]) const
{
    // fix so correct
    double xB[3], yB[3];
    double xybar[4];
    getCoefficients(solver, xB, yB, xybar);
    double x, y;
    x = solver->getColLower()[xColumn_];
    assert(x == solver->getColUpper()[xColumn_]);
    xB[2] = x;
    y = solver->getColLower()[yColumn_];
    assert(y == solver->getColUpper()[yColumn_]);
    yB[2] = y;
    computeLambdas(xB, yB, xybar, lambda);
    assert (xyRow_ >= 0);
}
// Get LU coefficients from matrix
void
OsiBiLinear::getCoefficients(const OsiSolverInterface * solver, double xB[2], double yB[2],
                             double xybar[4]) const
{
    const CoinPackedMatrix * matrix = solver->getMatrixByCol();
    const double * element = matrix->getElements();
    const double * objective = solver->getObjCoefficients();
    const int * row = matrix->getIndices();
    const CoinBigIndex * columnStart = matrix->getVectorStarts();
    const int * columnLength = matrix->getVectorLengths();
    // order is LxLy, LxUy, UxLy and UxUy
    int j;
    double multiplier = (boundType_ == 0) ? 1.0 / coefficient_ : 1.0;
    if (yRow_ >= 0) {
        for (j = 0; j < 4; j++) {
            int iColumn = firstLambda_ + j;
            int iStart = columnStart[iColumn];
            int iEnd = iStart + columnLength[iColumn];
            int k = iStart;
            double x = 0.0;
            double y = 0.0;
            xybar[j] = 0.0;
            for (; k < iEnd; k++) {
                if (xRow_ == row[k])
                    x = element[k];
                if (yRow_ == row[k])
                    y = element[k];
                if (xyRow_ == row[k])
                    xybar[j] = element[k] * multiplier;
            }
            if (xyRow_ < 0)
                xybar[j] = objective[iColumn] * multiplier;
            if (j == 0)
                xB[0] = x;
            else if (j == 1)
                yB[1] = y;
            else if (j == 2)
                yB[0] = y;
            else if (j == 3)
                xB[1] = x;
            assert (fabs(xybar[j] - x*y) < 1.0e-4);
        }
    } else {
        // x==y
        for (j = 0; j < 4; j++) {
            int iColumn = firstLambda_ + j;
            int iStart = columnStart[iColumn];
            int iEnd = iStart + columnLength[iColumn];
            int k = iStart;
            double x = 0.0;
            xybar[j] = 0.0;
            for (; k < iEnd; k++) {
                if (xRow_ == row[k])
                    x = element[k];
                if (xyRow_ == row[k])
                    xybar[j] = element[k] * multiplier;
            }
            if (xyRow_ < 0)
                xybar[j] = objective[iColumn] * multiplier;
            if (j == 0) {
                xB[0] = x;
                yB[0] = x;
            } else if (j == 2) {
                xB[1] = x;
                yB[1] = x;
            }
        }
        assert (fabs(xybar[0] - xB[0]*yB[0]) < 1.0e-4);
        assert (fabs(xybar[1] - xB[0]*yB[1]) < 1.0e-4);
        assert (fabs(xybar[2] - xB[1]*yB[0]) < 1.0e-4);
        assert (fabs(xybar[3] - xB[1]*yB[1]) < 1.0e-4);
    }
}
// Compute lambdas (third entry in each .B is current value)
double
OsiBiLinear::computeLambdas(const double xB[3], const double yB[3], const double xybar[4], double lambda[4]) const
{
    // fake to make correct
    double x = xB[2];
    double y = yB[2];
    // order is LxLy, LxUy, UxLy and UxUy
    // l0 + l1 = this
    double rhs1 = (xB[1] - x) / (xB[1] - xB[0]);
    // l0 + l2 = this
    double rhs2 = (yB[1] - y) / (yB[1] - yB[0]);
    // For xy (taking out l3)
    double rhs3 = xB[1] * yB[1] - x * y;
    double a0 = xB[1] * yB[1] - xB[0] * yB[0];
    double a1 = xB[1] * yB[1] - xB[0] * yB[1];
    double a2 = xB[1] * yB[1] - xB[1] * yB[0];
    // divide through to get l0 coefficient
    rhs3 /= a0;
    a1 /= a0;
    a2 /= a0;
    // subtract out l0
    double b[2][2];
    double rhs[2];
    // first for l1 and l2
    b[0][0] = 1.0 - a1;
    b[0][1] = -a2;
    rhs[0] = rhs1 - rhs3;
    // second
    b[1][0] = -a1;
    b[1][1] = 1.0 - a2;
    rhs[1] = rhs2 - rhs3;
    if (fabs(b[0][0]) > fabs(b[0][1])) {
        double sub = b[1][0] / b[0][0];
        b[1][1] -= sub * b[0][1];
        rhs[1] -= sub * rhs[0];
        assert (fabs(b[1][1]) > 1.0e-12);
        lambda[2] = rhs[1] / b[1][1];
        lambda[0] = rhs2 - lambda[2];
        lambda[1] = rhs1 - lambda[0];
    } else {
        double sub = b[1][1] / b[0][1];
        b[1][0] -= sub * b[0][0];
        rhs[1] -= sub * rhs[0];
        assert (fabs(b[1][0]) > 1.0e-12);
        lambda[1] = rhs[1] / b[1][0];
        lambda[0] = rhs1 - lambda[1];
        lambda[2] = rhs2 - lambda[0];
    }
    lambda[3] = 1.0 - (lambda[0] + lambda[1] + lambda[2]);
    double infeasibility = 0.0;
    double xy = 0.0;
    for (int j = 0; j < 4; j++) {
        double value = lambda[j];
        if (value > 1.0) {
            infeasibility += value - 1.0;
            value = 1.0;
        }
        if (value < 0.0) {
            infeasibility -= value;
            value = 0.0;
        }
        lambda[j] = value;
        xy += xybar[j] * value;
    }
    assert (fabs(xy - x*y) < 1.0e-4);
    return infeasibility;
}
// Updates coefficients
int
OsiBiLinear::updateCoefficients(const double * lower, const double * upper, double * objective,
                                CoinPackedMatrix * matrix, CoinWarmStartBasis * basis) const
{
    // Return if no updates
    if ((branchingStrategy_&4) != 0)
        return 0;
    int numberUpdated = 0;
    double * element = matrix->getMutableElements();
    const int * row = matrix->getIndices();
    const CoinBigIndex * columnStart = matrix->getVectorStarts();
    const int * columnLength = matrix->getVectorLengths();
    // order is LxLy, LxUy, UxLy and UxUy
    double xB[2];
    double yB[2];
    xB[0] = lower[xColumn_];
    xB[1] = upper[xColumn_];
    yB[0] = lower[yColumn_];
    yB[1] = upper[yColumn_];
    //printf("x %d (%g,%g) y %d (%g,%g)\n",
    // xColumn_,xB[0],xB[1],
    // yColumn_,yB[0],yB[1]);
    CoinWarmStartBasis::Status status[4];
    int numStruct = basis ? basis->getNumStructural() - firstLambda_ : 0;
    double coefficient = (boundType_ == 0) ? coefficient_ : 1.0;
    for (int j = 0; j < 4; j++) {
        status[j] = (j < numStruct) ? basis->getStructStatus(j + firstLambda_) : CoinWarmStartBasis::atLowerBound;
        int iX = j >> 1;
        double x = xB[iX];
        int iY = j & 1;
        double y = yB[iY];
        CoinBigIndex k = columnStart[j+firstLambda_];
        CoinBigIndex last = k + columnLength[j+firstLambda_];
        double value;
        // xy
        value = coefficient * x * y;
        if (xyRow_ >= 0) {
            assert (row[k] == xyRow_);
#if BI_PRINT > 1
            printf("j %d xy (%d,%d) coeff from %g to %g\n", j, xColumn_, yColumn_, element[k], value);
#endif
            element[k++] = value;
        } else {
            // objective
            objective[j+firstLambda_] = value;
        }
        numberUpdated++;
        // convexity
        assert (row[k] == convexity_);
        k++;
        // x
        value = x;
#if BI_PRINT > 1
        printf("j %d x (%d) coeff from %g to %g\n", j, xColumn_, element[k], value);
#endif
        assert (row[k] == xRow_);
        element[k++] = value;
        numberUpdated++;
        if (yRow_ >= 0) {
            // y
            value = y;
#if BI_PRINT > 1
            printf("j %d y (%d) coeff from %g to %g\n", j, yColumn_, element[k], value);
#endif
            assert (row[k] == yRow_);
            element[k++] = value;
            numberUpdated++;
        }
        // Do extra rows
        for (int i = 0; i < numberExtraRows_; i++) {
            int iRow = extraRow_[i];
            for (; k < last; k++) {
                if (row[k] == iRow)
                    break;
            }
            assert (k < last);
            element[k++] = x * y * multiplier_[i];
        }
    }

    if (xB[0] == xB[1]) {
        if (yB[0] == yB[1]) {
            // only one basic
            bool first = true;
            for (int j = 0; j < 4; j++) {
                if (status[j] == CoinWarmStartBasis::basic) {
                    if (first) {
                        first = false;
                    } else {
                        basis->setStructStatus(j + firstLambda_, CoinWarmStartBasis::atLowerBound);
#if BI_PRINT
                        printf("zapping %d (x=%d,y=%d)\n", j, xColumn_, yColumn_);
#endif
                    }
                }
            }
        } else {
            if (status[0] == CoinWarmStartBasis::basic &&
                    status[2] == CoinWarmStartBasis::basic) {
                basis->setStructStatus(2 + firstLambda_, CoinWarmStartBasis::atLowerBound);
#if BI_PRINT
                printf("zapping %d (x=%d,y=%d)\n", 2, xColumn_, yColumn_);
#endif
            }
            if (status[1] == CoinWarmStartBasis::basic &&
                    status[3] == CoinWarmStartBasis::basic) {
                basis->setStructStatus(3 + firstLambda_, CoinWarmStartBasis::atLowerBound);
#if BI_PRINT
                printf("zapping %d (x=%d,y=%d)\n", 3, xColumn_, yColumn_);
#endif
            }
        }
    } else if (yB[0] == yB[1]) {
        if (status[0] == CoinWarmStartBasis::basic &&
                status[1] == CoinWarmStartBasis::basic) {
            basis->setStructStatus(1 + firstLambda_, CoinWarmStartBasis::atLowerBound);
#if BI_PRINT
            printf("zapping %d (x=%d,y=%d)\n", 1, xColumn_, yColumn_);
#endif
        }
        if (status[2] == CoinWarmStartBasis::basic &&
                status[3] == CoinWarmStartBasis::basic) {
            basis->setStructStatus(3 + firstLambda_, CoinWarmStartBasis::atLowerBound);
#if BI_PRINT
            printf("zapping %d (x=%d,y=%d)\n", 3, xColumn_, yColumn_);
#endif
        }
    }
    return numberUpdated;
}
// This does NOT set mutable stuff
double
OsiBiLinear::checkInfeasibility(const OsiBranchingInformation * info) const
{
    // If another object has finer mesh ignore this
    if ((branchingStrategy_&8) != 0)
        return 0.0;
    int way;
    double saveInfeasibility = infeasibility_;
    short int saveWhichWay = whichWay_;
    double saveXyBranchValue = xyBranchValue_;
    short saveChosen = chosen_;
    double value = infeasibility(info, way);
    infeasibility_ = saveInfeasibility;
    whichWay_ = saveWhichWay;
    xyBranchValue_ = saveXyBranchValue;
    chosen_ = saveChosen;
    return value;
}
OsiBiLinearBranchingObject::OsiBiLinearBranchingObject()
        : OsiTwoWayBranchingObject(),
        chosen_(0)
{
}

// Useful constructor
OsiBiLinearBranchingObject::OsiBiLinearBranchingObject (OsiSolverInterface * solver,
        const OsiBiLinear * set,
        int way ,
        double separator,
        int chosen)
        : OsiTwoWayBranchingObject(solver, set, way, separator),
        chosen_(static_cast<short int>(chosen))
{
    assert (chosen_ >= 0 && chosen_ < 2);
}

// Copy constructor
OsiBiLinearBranchingObject::OsiBiLinearBranchingObject ( const OsiBiLinearBranchingObject & rhs)
        : OsiTwoWayBranchingObject(rhs),
        chosen_(rhs.chosen_)
{
}

// Assignment operator
OsiBiLinearBranchingObject &
OsiBiLinearBranchingObject::operator=( const OsiBiLinearBranchingObject & rhs)
{
    if (this != &rhs) {
        OsiTwoWayBranchingObject::operator=(rhs);
        chosen_ = rhs.chosen_;
    }
    return *this;
}
OsiBranchingObject *
OsiBiLinearBranchingObject::clone() const
{
    return (new OsiBiLinearBranchingObject(*this));
}


// Destructor
OsiBiLinearBranchingObject::~OsiBiLinearBranchingObject ()
{
}
double
OsiBiLinearBranchingObject::branch(OsiSolverInterface * solver)
{
    const OsiBiLinear * set =
        dynamic_cast <const OsiBiLinear *>(originalObject_) ;
    assert (set);
    int way = (!branchIndex_) ? (2 * firstBranch_ - 1) : -(2 * firstBranch_ - 1);
    branchIndex_++;
    set->newBounds(solver, way, chosen_, value_);
    return 0.0;
}
/* Return true if branch should only bound variables
 */
bool
OsiBiLinearBranchingObject::boundBranch() const
{
    const OsiBiLinear * set =
        dynamic_cast <const OsiBiLinear *>(originalObject_) ;
    assert (set);
    return (set->branchingStrategy()&4) != 0;
}
// Print what would happen
void
OsiBiLinearBranchingObject::print(const OsiSolverInterface * /*solver*/)
{
    const OsiBiLinear * set =
        dynamic_cast <const OsiBiLinear *>(originalObject_) ;
    assert (set);
    int way = (!branchIndex_) ? (2 * firstBranch_ - 1) : -(2 * firstBranch_ - 1);
    int iColumn = (chosen_ == 1) ? set->xColumn() : set->yColumn();
    printf("OsiBiLinear would branch %s on %c variable %d from value %g\n",
           (way < 0) ? "down" : "up",
           (chosen_ == 0) ? 'X' : 'Y', iColumn, value_);
}
// Default Constructor
OsiBiLinearEquality::OsiBiLinearEquality ()
        : OsiBiLinear(),
        numberPoints_(0)
{
}

// Useful constructor
OsiBiLinearEquality::OsiBiLinearEquality (OsiSolverInterface * solver, int xColumn,
        int yColumn, int xyRow, double rhs,
        double xMesh)
        : OsiBiLinear(),
        numberPoints_(0)
{
    double xB[2];
    double yB[2];
    const double * lower = solver->getColLower();
    const double * upper = solver->getColUpper();
    xColumn_ = xColumn;
    yColumn_ = yColumn;
    xyRow_ = xyRow;
    coefficient_ = rhs;
    xB[0] = lower[xColumn_];
    xB[1] = upper[xColumn_];
    yB[0] = lower[yColumn_];
    yB[1] = upper[yColumn_];
    if (xB[1]*yB[1] < coefficient_ + 1.0e-12 ||
            xB[0]*yB[0] > coefficient_ - 1.0e-12) {
        printf("infeasible row - reformulate\n");
        abort();
    }
    // reduce range of x if possible
    if (yB[0]*xB[1] > coefficient_ + 1.0e12) {
        xB[1] = coefficient_ / yB[0];
        solver->setColUpper(xColumn_, xB[1]);
    }
    if (yB[1]*xB[0] < coefficient_ - 1.0e12) {
        xB[0] = coefficient_ / yB[1];
        solver->setColLower(xColumn_, xB[0]);
    }
    // See how many points
    numberPoints_ = static_cast<int> ((xB[1] - xB[0] + 0.5 * xMesh) / xMesh);
    // redo exactly
    xMeshSize_ = (xB[1] - xB[0]) / static_cast<double> (numberPoints_);
    numberPoints_++;
    //#define KEEPXY
#ifndef KEEPXY
    // Take out xyRow
    solver->setRowLower(xyRow_, 0.0);
    solver->setRowUpper(xyRow_, 0.0);
#else
    // make >=
    solver->setRowLower(xyRow_, coefficient_ - 0.05);
    solver->setRowUpper(xyRow_, COIN_DBL_MAX);
#endif
    double rowLower[3];
    double rowUpper[3];
#ifndef KEEPXY
    double * columnLower = new double [numberPoints_];
    double * columnUpper = new double [numberPoints_];
    double * objective = new double [numberPoints_];
    CoinBigIndex *starts = new CoinBigIndex[numberPoints_+1];
    int * index = new int[3*numberPoints_];
    double * element = new double [3*numberPoints_];
#else
    double * columnLower = new double [numberPoints_+2];
    double * columnUpper = new double [numberPoints_+2];
    double * objective = new double [numberPoints_+2];
    CoinBigIndex *starts = new CoinBigIndex[numberPoints_+3];
    int * index = new int[4*numberPoints_+2];
    double * element = new double [4*numberPoints_+2];
#endif
    int i;
    starts[0] = 0;
    // rows
    int numberRows = solver->getNumRows();
    // convexity
    rowLower[0] = 1.0;
    rowUpper[0] = 1.0;
    convexity_ = numberRows;
    starts[1] = 0;
    // x
    rowLower[1] = 0.0;
    rowUpper[1] = 0.0;
    index[0] = xColumn_;
    element[0] = -1.0;
    xRow_ = numberRows + 1;
    starts[2] = 1;
    rowLower[2] = 0.0;
    rowUpper[2] = 0.0;
    index[1] = yColumn;
    element[1] = -1.0;
    yRow_ = numberRows + 2;
    starts[3] = 2;
    solver->addRows(3, starts, index, element, rowLower, rowUpper);
    int n = 0;
    firstLambda_ = solver->getNumCols();
    double x = xB[0];
    assert(xColumn_ != yColumn_);
    for (i = 0; i < numberPoints_; i++) {
        double y = coefficient_ / x;
        columnLower[i] = 0.0;
        columnUpper[i] = 2.0;
        objective[i] = 0.0;
        double value;
#ifdef KEEPXY
        // xy
        value = coefficient_;
        element[n] = value;
        index[n++] = xyRow_;
#endif
        // convexity
        value = 1.0;
        element[n] = value;
        index[n++] = 0 + numberRows;
        // x
        value = x;
        if (fabs(value) < 1.0e-19)
            value = 1.0e-19;
        element[n] = value;
        index[n++] = 1 + numberRows;
        // y
        value = y;
        if (fabs(value) < 1.0e-19)
            value = 1.0e-19;
        element[n] = value;
        index[n++] = 2 + numberRows;
        starts[i+1] = n;
        x += xMeshSize_;
    }
#ifdef KEEPXY
    // costed slacks
    columnLower[numberPoints_] = 0.0;
    columnUpper[numberPoints_] = xMeshSize_;
    objective[numberPoints_] = 1.0e3;;
    // convexity
    element[n] = 1.0;
    index[n++] = 0 + numberRows;
    starts[numberPoints_+1] = n;
    columnLower[numberPoints_+1] = 0.0;
    columnUpper[numberPoints_+1] = xMeshSize_;
    objective[numberPoints_+1] = 1.0e3;;
    // convexity
    element[n] = -1.0;
    index[n++] = 0 + numberRows;
    starts[numberPoints_+2] = n;
    solver->addCols(numberPoints_ + 2, starts, index, element, columnLower, columnUpper, objective);
#else
    solver->addCols(numberPoints_, starts, index, element, columnLower, columnUpper, objective);
#endif
    delete [] columnLower;
    delete [] columnUpper;
    delete [] objective;
    delete [] starts;
    delete [] index;
    delete [] element;
}

// Copy constructor
OsiBiLinearEquality::OsiBiLinearEquality ( const OsiBiLinearEquality & rhs)
        : OsiBiLinear(rhs),
        numberPoints_(rhs.numberPoints_)
{
}

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

// Assignment operator
OsiBiLinearEquality &
OsiBiLinearEquality::operator=( const OsiBiLinearEquality & rhs)
{
    if (this != &rhs) {
        OsiBiLinear::operator=(rhs);
        numberPoints_ = rhs.numberPoints_;
    }
    return *this;
}

// Destructor
OsiBiLinearEquality::~OsiBiLinearEquality ()
{
}
// Possible improvement
double
OsiBiLinearEquality::improvement(const OsiSolverInterface * solver) const
{
    const double * pi = solver->getRowPrice();
    int i;
    const double * solution = solver->getColSolution();
    printf(" for x %d y %d - pi %g %g\n", xColumn_, yColumn_, pi[xRow_], pi[yRow_]);
    for (i = 0; i < numberPoints_; i++) {
        if (fabs(solution[i+firstLambda_]) > 1.0e-7)
            printf("(%d %g) ", i, solution[i+firstLambda_]);
    }
    printf("\n");
    return 0.0;
}
/* change grid
   if type 0 then use solution and make finer
   if 1 then back to original
*/
double
OsiBiLinearEquality::newGrid(OsiSolverInterface * solver, int type) const
{
    CoinPackedMatrix * matrix = solver->getMutableMatrixByCol();
    if (!matrix) {
        printf("Unable to modify matrix\n");
        abort();
    }
    double * element = matrix->getMutableElements();
#ifndef NDEBUG
    const int * row = matrix->getIndices();
#endif
    const CoinBigIndex * columnStart = matrix->getVectorStarts();
    //const int * columnLength = matrix->getVectorLengths();
    // get original bounds
    double xB[2];
    const double * lower = solver->getColLower();
    const double * upper = solver->getColUpper();
    xB[0] = lower[xColumn_];
    xB[1] = upper[xColumn_];
    assert (fabs((xB[1] - xB[0]) - xMeshSize_*(numberPoints_ - 1)) < 1.0e-7);
    double mesh = 0.0;
    int i;
    if (type == 0) {
        const double * solution = solver->getColSolution();
        int first = -1;
        int last = -1;
        double xValue = 0.0;
        double step = 0.0;
        for (i = 0; i < numberPoints_; i++) {
            int iColumn = i + firstLambda_;
            if (fabs(solution[iColumn]) > 1.0e-7) {
                int k = columnStart[iColumn] + 1;
                xValue += element[k] * solution[iColumn];
                if (first == -1) {
                    first = i;
                    step = -element[k];
                } else {
                    step += element[k];
                }
                last = i;
            }
        }
        if (last > first + 1) {
            printf("not adjacent - presuming small djs\n");
        }
        // new step size
        assert (numberPoints_ > 2);
        step = CoinMax((1.5 * step) / static_cast<double> (numberPoints_ - 1), 0.5 * step);
        xB[0] = CoinMax(xB[0], xValue - 0.5 * step);
        xB[1] = CoinMin(xB[1], xValue + 0.5 * step);
        // and now divide these
        mesh = (xB[1] - xB[0]) / static_cast<double> (numberPoints_ - 1);
    } else {
        // back to original
        mesh = xMeshSize_;
    }
    double x = xB[0];
    for (i = 0; i < numberPoints_; i++) {
        int iColumn = i + firstLambda_;
        double y = coefficient_ / x;
        //assert (columnLength[iColumn]==3); - could have cuts
        int k = columnStart[iColumn];
#ifdef KEEPXY
        // xy
        assert (row[k] == xyRow_);
        k++;
#endif
        assert (row[k] == convexity_);
        k++;
        double value;
        // x
        value = x;
        assert (row[k] == xRow_);
        assert (fabs(value) > 1.0e-10);
        element[k++] = value;
        // y
        value = y;
        assert (row[k] == yRow_);
        assert (fabs(value) > 1.0e-10);
        element[k++] = value;
        x += mesh;
    }
    return mesh;
}
/** Default Constructor

  Equivalent to an unspecified binary variable.
*/
OsiSimpleFixedInteger::OsiSimpleFixedInteger ()
        : OsiSimpleInteger()
{
}

/** Useful constructor

  Loads actual upper & lower bounds for the specified variable.
*/
OsiSimpleFixedInteger::OsiSimpleFixedInteger (const OsiSolverInterface * solver, int iColumn)
        : OsiSimpleInteger(solver, iColumn)
{
}


// Useful constructor - passed solver index and original bounds
OsiSimpleFixedInteger::OsiSimpleFixedInteger ( int iColumn, double lower, double upper)
        : OsiSimpleInteger(iColumn, lower, upper)
{
}

// Useful constructor - passed simple integer
OsiSimpleFixedInteger::OsiSimpleFixedInteger ( const OsiSimpleInteger &rhs)
        : OsiSimpleInteger(rhs)
{
}

// Copy constructor
OsiSimpleFixedInteger::OsiSimpleFixedInteger ( const OsiSimpleFixedInteger & rhs)
        : OsiSimpleInteger(rhs)

{
}

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

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

// Destructor
OsiSimpleFixedInteger::~OsiSimpleFixedInteger ()
{
}
// Infeasibility - large is 0.5
double
OsiSimpleFixedInteger::infeasibility(const OsiBranchingInformation * info, int & whichWay) const
{
    double value = info->solution_[columnNumber_];
    value = CoinMax(value, info->lower_[columnNumber_]);
    value = CoinMin(value, info->upper_[columnNumber_]);
    double nearest = floor(value + (1.0 - 0.5));
    if (nearest > value) {
        whichWay = 1;
    } else {
        whichWay = 0;
    }
    infeasibility_ = fabs(value - nearest);
    bool satisfied = false;
    if (infeasibility_ <= info->integerTolerance_) {
        otherInfeasibility_ = 1.0;
        satisfied = true;
        if (info->lower_[columnNumber_] != info->upper_[columnNumber_])
            infeasibility_ = 1.0e-5;
        else
            infeasibility_ = 0.0;
    } else if (info->defaultDual_ < 0.0) {
        otherInfeasibility_ = 1.0 - infeasibility_;
    } else {
        const double * pi = info->pi_;
        const double * activity = info->rowActivity_;
        const double * lower = info->rowLower_;
        const double * upper = info->rowUpper_;
        const double * element = info->elementByColumn_;
        const int * row = info->row_;
        const CoinBigIndex * columnStart = info->columnStart_;
        const int * columnLength = info->columnLength_;
        double direction = info->direction_;
        double downMovement = value - floor(value);
        double upMovement = 1.0 - downMovement;
        double valueP = info->objective_[columnNumber_] * direction;
        CoinBigIndex start = columnStart[columnNumber_];
        CoinBigIndex end = start + columnLength[columnNumber_];
        double upEstimate = 0.0;
        double downEstimate = 0.0;
        if (valueP > 0.0)
            upEstimate = valueP * upMovement;
        else
            downEstimate -= valueP * downMovement;
        double tolerance = info->primalTolerance_;
        for (CoinBigIndex j = start; j < end; j++) {
            int iRow = row[j];
            if (lower[iRow] < -1.0e20)
                assert (pi[iRow] <= 1.0e-3);
            if (upper[iRow] > 1.0e20)
                assert (pi[iRow] >= -1.0e-3);
            valueP = pi[iRow] * direction;
            double el2 = element[j];
            double value2 = valueP * el2;
            double u = 0.0;
            double d = 0.0;
            if (value2 > 0.0)
                u = value2;
            else
                d = -value2;
            // if up makes infeasible then make at least default
            double newUp = activity[iRow] + upMovement * el2;
            if (newUp > upper[iRow] + tolerance || newUp < lower[iRow] - tolerance)
                u = CoinMax(u, info->defaultDual_);
            upEstimate += u * upMovement * fabs(el2);
            // if down makes infeasible then make at least default
            double newDown = activity[iRow] - downMovement * el2;
            if (newDown > upper[iRow] + tolerance || newDown < lower[iRow] - tolerance)
                d = CoinMax(d, info->defaultDual_);
            downEstimate += d * downMovement * fabs(el2);
        }
        if (downEstimate >= upEstimate) {
            infeasibility_ = CoinMax(1.0e-12, upEstimate);
            otherInfeasibility_ = CoinMax(1.0e-12, downEstimate);
            whichWay = 1;
        } else {
            infeasibility_ = CoinMax(1.0e-12, downEstimate);
            otherInfeasibility_ = CoinMax(1.0e-12, upEstimate);
            whichWay = 0;
        }
    }
    if (preferredWay_ >= 0 && !satisfied)
        whichWay = preferredWay_;
    whichWay_ = static_cast<short int>(whichWay);
    return infeasibility_;
}
// Creates a branching object
OsiBranchingObject *
OsiSimpleFixedInteger::createBranch(OsiSolverInterface * solver, const OsiBranchingInformation * info, int way) const
{
    double value = info->solution_[columnNumber_];
    value = CoinMax(value, info->lower_[columnNumber_]);
    value = CoinMin(value, info->upper_[columnNumber_]);
    assert (info->upper_[columnNumber_] > info->lower_[columnNumber_]);
    double nearest = floor(value + 0.5);
    double integerTolerance = info->integerTolerance_;
    if (fabs(value - nearest) < integerTolerance) {
        // adjust value
        if (nearest != info->upper_[columnNumber_])
            value = nearest + 2.0 * integerTolerance;
        else
            value = nearest - 2.0 * integerTolerance;
    }
    OsiBranchingObject * branch = new OsiIntegerBranchingObject(solver, this, way,
            value);
    return branch;
}

#include <cstdlib>
#include <cstdio>
#include <cmath>
#include <cfloat>
#include <cassert>
#include <iostream>
//#define CGL_DEBUG 2
#include "CoinHelperFunctions.hpp"
#include "CoinPackedVector.hpp"
#include "OsiRowCutDebugger.hpp"
#include "CoinWarmStartBasis.hpp"
//#include "CglTemporary.hpp"
#include "CoinFinite.hpp"
//-------------------------------------------------------------------
// Generate Stored cuts
//-------------------------------------------------------------------
void
CglTemporary::generateCuts(const OsiSolverInterface & si, OsiCuts & cs,
                           const CglTreeInfo /*info*/)
{
    // Get basic problem information
    const double * solution = si.getColSolution();
    int numberRowCuts = cuts_.sizeRowCuts();
    for (int i = 0; i < numberRowCuts; i++) {
        const OsiRowCut * rowCutPointer = cuts_.rowCutPtr(i);
        double violation = rowCutPointer->violated(solution);
        if (violation >= requiredViolation_)
            cs.insert(*rowCutPointer);
    }
    // delete
    cuts_ = OsiCuts();
}

//-------------------------------------------------------------------
// Default Constructor
//-------------------------------------------------------------------
CglTemporary::CglTemporary ()
        :
        CglStored()
{
}

//-------------------------------------------------------------------
// Copy constructor
//-------------------------------------------------------------------
CglTemporary::CglTemporary (const CglTemporary & source) :
        CglStored(source)
{
}

//-------------------------------------------------------------------
// Clone
//-------------------------------------------------------------------
CglCutGenerator *
CglTemporary::clone() const
{
    return new CglTemporary(*this);
}

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

//----------------------------------------------------------------
// Assignment operator
//-------------------------------------------------------------------
CglTemporary &
CglTemporary::operator=(const CglTemporary & rhs)
{
    if (this != &rhs) {
        CglStored::operator=(rhs);
    }
    return *this;
}
void checkQP(ClpSimplex * /*model*/)
{
#ifdef JJF_ZERO
    printf("Checking quadratic model %x\n", model);
    if (model) {
        ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(model->objectiveAsObject()));
        assert (quadraticObj);
        CoinPackedMatrix * quadraticObjective = quadraticObj->quadraticObjective();
        int numberColumns = quadraticObj->numberColumns();
        const int * columnQuadratic = quadraticObjective->getIndices();
        const CoinBigIndex * columnQuadraticStart = quadraticObjective->getVectorStarts();
        const int * columnQuadraticLength = quadraticObjective->getVectorLengths();
        //const double * quadraticElement = quadraticObjective->getElements();
        for (int i = 0; i < numberColumns; i++) {
            for (int j = columnQuadraticStart[i]; j < columnQuadraticStart[i] + columnQuadraticLength[i]; j++)
                assert (columnQuadratic[j] >= 0 && columnQuadratic[j] < 1000);
        }
    }
#endif
}
//#############################################################################
// Solve methods
//#############################################################################
void OsiSolverLinearizedQuadratic::initialSolve()
{
    OsiClpSolverInterface::initialSolve();
    int secondaryStatus = modelPtr_->secondaryStatus();
    if (modelPtr_->status() == 0 && (secondaryStatus == 2 || secondaryStatus == 4))
        modelPtr_->cleanup(1);
    if (isProvenOptimal() && modelPtr_->numberColumns() == quadraticModel_->numberColumns()) {
        // see if qp can get better solution
        const double * solution = modelPtr_->primalColumnSolution();
        int numberColumns = modelPtr_->numberColumns();
        bool satisfied = true;
        for (int i = 0; i < numberColumns; i++) {
            if (isInteger(i)) {
                double value = solution[i];
                if (fabs(value - floor(value + 0.5)) > 1.0e-6) {
                    satisfied = false;
                    break;
                }
            }
        }
        if (satisfied) {
            checkQP(quadraticModel_);
            ClpSimplex qpTemp(*quadraticModel_);
            checkQP(&qpTemp);
            double * lower = qpTemp.columnLower();
            double * upper = qpTemp.columnUpper();
            double * lower2 = modelPtr_->columnLower();
            double * upper2 = modelPtr_->columnUpper();
            for (int i = 0; i < numberColumns; i++) {
                if (isInteger(i)) {
                    double value = floor(solution[i] + 0.5);
                    lower[i] = value;
                    upper[i] = value;
                } else {
                    lower[i] = lower2[i];
                    upper[i] = upper2[i];
                }
            }
            //qpTemp.writeMps("bad.mps");
            //modelPtr_->writeMps("bad2.mps");
            //qpTemp.objectiveAsObject()->setActivated(0);
            //qpTemp.primal();
            //qpTemp.objectiveAsObject()->setActivated(1);
            qpTemp.primal();
            //assert (!qpTemp.problemStatus());
            if (qpTemp.objectiveValue() < bestObjectiveValue_ && !qpTemp.problemStatus()) {
                delete [] bestSolution_;
                bestSolution_ = CoinCopyOfArray(qpTemp.primalColumnSolution(), numberColumns);
                bestObjectiveValue_ = qpTemp.objectiveValue();
                printf("better qp objective of %g\n", bestObjectiveValue_);
            }
        }
    }
}
//#############################################################################
// Constructors, destructors clone and assignment
//#############################################################################
//-------------------------------------------------------------------
// Default Constructor
//-------------------------------------------------------------------
OsiSolverLinearizedQuadratic::OsiSolverLinearizedQuadratic ()
        : OsiClpSolverInterface()
{
    bestObjectiveValue_ = COIN_DBL_MAX;
    bestSolution_ = NULL;
    specialOptions3_ = 0;
    quadraticModel_ = NULL;
}
OsiSolverLinearizedQuadratic::OsiSolverLinearizedQuadratic ( ClpSimplex * quadraticModel)
        : OsiClpSolverInterface(new ClpSimplex(*quadraticModel), true)
{
    bestObjectiveValue_ = COIN_DBL_MAX;
    bestSolution_ = NULL;
    specialOptions3_ = 0;
    quadraticModel_ = new ClpSimplex(*quadraticModel);
    // linearize
    int numberColumns = modelPtr_->numberColumns();
    const double * solution = modelPtr_->primalColumnSolution();
    // Replace objective
    ClpObjective * trueObjective = modelPtr_->objectiveAsObject();
    ClpObjective * objective = new ClpLinearObjective(NULL, numberColumns);
    modelPtr_->setObjectivePointer(objective);
    double offset;
    double saveOffset = modelPtr_->objectiveOffset();
    memcpy(modelPtr_->objective(), trueObjective->gradient(modelPtr_, solution, offset, true, 2),
           numberColumns*sizeof(double));
    modelPtr_->setObjectiveOffset(saveOffset + offset);
    delete trueObjective;
    checkQP(quadraticModel_);
}
//-------------------------------------------------------------------
// Clone
//-------------------------------------------------------------------
OsiSolverInterface *
OsiSolverLinearizedQuadratic::clone(bool /*copyData*/) const
{
    //assert (copyData);
    return new OsiSolverLinearizedQuadratic(*this);
}


//-------------------------------------------------------------------
// Copy constructor
//-------------------------------------------------------------------
OsiSolverLinearizedQuadratic::OsiSolverLinearizedQuadratic (
    const OsiSolverLinearizedQuadratic & rhs)
        : OsiSolverInterface(rhs)
        , OsiClpSolverInterface(rhs)
{
    bestObjectiveValue_ = rhs.bestObjectiveValue_;
    if (rhs.bestSolution_) {
        bestSolution_ = CoinCopyOfArray(rhs.bestSolution_, modelPtr_->numberColumns());
    } else {
        bestSolution_ = NULL;
    }
    specialOptions3_ = rhs.specialOptions3_;
    if (rhs.quadraticModel_) {
        quadraticModel_ = new ClpSimplex(*rhs.quadraticModel_);
    } else {
        quadraticModel_ = NULL;
    }
    checkQP(rhs.quadraticModel_);
    checkQP(quadraticModel_);
}

//-------------------------------------------------------------------
// Destructor
//-------------------------------------------------------------------
OsiSolverLinearizedQuadratic::~OsiSolverLinearizedQuadratic ()
{
    delete [] bestSolution_;
    delete quadraticModel_;
}

//-------------------------------------------------------------------
// Assignment operator
//-------------------------------------------------------------------
OsiSolverLinearizedQuadratic &
OsiSolverLinearizedQuadratic::operator=(const OsiSolverLinearizedQuadratic & rhs)
{
    if (this != &rhs) {
        delete [] bestSolution_;
        delete quadraticModel_;
        OsiClpSolverInterface::operator=(rhs);
        bestObjectiveValue_ = rhs.bestObjectiveValue_;
        if (rhs.bestSolution_) {
            bestSolution_ = CoinCopyOfArray(rhs.bestSolution_, modelPtr_->numberColumns());
        } else {
            bestSolution_ = NULL;
        }
        specialOptions3_ = rhs.specialOptions3_;
        if (rhs.quadraticModel_) {
            quadraticModel_ = new ClpSimplex(*rhs.quadraticModel_);
        } else {
            quadraticModel_ = NULL;
        }
        checkQP(rhs.quadraticModel_);
        checkQP(quadraticModel_);
    }
    return *this;
}
/* Expands out all possible combinations for a knapsack
   If buildObj NULL then just computes space needed - returns number elements
   On entry numberOutput is maximum allowed, on exit it is number needed or
   -1 (as will be number elements) if maximum exceeded.  numberOutput will have at
   least space to return values which reconstruct input.
   Rows returned will be original rows but no entries will be returned for
   any rows all of whose entries are in knapsack.  So up to user to allow for this.
   If reConstruct >=0 then returns number of entrie which make up item "reConstruct"
   in expanded knapsack.  Values in buildRow and buildElement;
*/
int
CoinModel::expandKnapsack(int knapsackRow, int & numberOutput, double * buildObj, CoinBigIndex * buildStart,
                          int * buildRow, double * buildElement, int reConstruct) const
{
    /* mark rows
       -2 in knapsack and other variables
       -1 not involved
       0 only in knapsack
    */
    int * markRow = new int [numberRows_];
    int iRow;
    int iColumn;
    int * whichColumn = new int [numberColumns_];
    for (iColumn = 0; iColumn < numberColumns_; iColumn++)
        whichColumn[iColumn] = -1;
    int numJ = 0;
    for (iRow = 0; iRow < numberRows_; iRow++)
        markRow[iRow] = -1;
    CoinModelLink triple;
    triple = firstInRow(knapsackRow);
    while (triple.column() >= 0) {
        int iColumn = triple.column();
#ifndef NDEBUG
        const char *  el = getElementAsString(knapsackRow, iColumn);
        assert (!strcmp("Numeric", el));
#endif
        whichColumn[iColumn] = numJ;
        numJ++;
        triple = next(triple);
    }
    for (iRow = 0; iRow < numberRows_; iRow++) {
        triple = firstInRow(iRow);
        int type = -3;
        while (triple.column() >= 0) {
            int iColumn = triple.column();
            if (whichColumn[iColumn] >= 0) {
                if (type == -3)
                    type = 0;
                else if (type != 0)
                    type = -2;
            } else {
                if (type == -3)
                    type = -1;
                else if (type == 0)
                    type = -2;
            }
            triple = next(triple);
        }
        if (type == -3)
            type = -1;
        markRow[iRow] = type;
    }
    int * bound = new int [numberColumns_+1];
    int * whichRow = new int [numberRows_];
    ClpSimplex tempModel;
    CoinModel tempModel2(*this);
    tempModel.loadProblem(tempModel2);
    int * stack = new int [numberColumns_+1];
    double * size = new double [numberColumns_+1];
    double * rhsOffset = new double[numberRows_];
    int * build = new int[numberColumns_];
    int maxNumber = numberOutput;
    numJ = 0;
    double minSize = getRowLower(knapsackRow);
    double maxSize = getRowUpper(knapsackRow);
    double offset = 0.0;
    triple = firstInRow(knapsackRow);
    while (triple.column() >= 0) {
        iColumn = triple.column();
        double lowerColumn = columnLower(iColumn);
        double upperColumn = columnUpper(iColumn);
        double gap = upperColumn - lowerColumn;
        if (gap > 1.0e8)
            gap = 1.0e8;
        assert (fabs(floor(gap + 0.5) - gap) < 1.0e-5);
        whichColumn[numJ] = iColumn;
        bound[numJ] = static_cast<int> (gap);
        size[numJ++] = triple.value();
        offset += triple.value() * lowerColumn;
        triple = next(triple);
    }
    int jRow;
    for (iRow = 0; iRow < numberRows_; iRow++)
        whichRow[iRow] = iRow;
    ClpSimplex smallModel(&tempModel, numberRows_, whichRow, numJ, whichColumn, true, true, true);
    // modify rhs to allow for nonzero lower bounds
    double * rowLower = smallModel.rowLower();
    double * rowUpper = smallModel.rowUpper();
    const double * columnLower = smallModel.columnLower();
    //const double * columnUpper = smallModel.columnUpper();
    const CoinPackedMatrix * matrix = smallModel.matrix();
    const double * element = matrix->getElements();
    const int * row = matrix->getIndices();
    const CoinBigIndex * columnStart = matrix->getVectorStarts();
    const int * columnLength = matrix->getVectorLengths();
    const double * objective = smallModel.objective();
    double objectiveOffset = 0.0;
    CoinZeroN(rhsOffset, numberRows_);
    for (iColumn = 0; iColumn < numJ; iColumn++) {
        double lower = columnLower[iColumn];
        if (lower) {
            objectiveOffset += objective[iColumn];
            for (CoinBigIndex j = columnStart[iColumn];
                    j < columnStart[iColumn] + columnLength[iColumn]; j++) {
                double value = element[j] * lower;
                int kRow = row[j];
                rhsOffset[kRow] += value;
                if (rowLower[kRow] > -1.0e20)
                    rowLower[kRow] -= value;
                if (rowUpper[kRow] < 1.0e20)
                    rowUpper[kRow] -= value;
            }
        }
    }
    // relax
    for (jRow = 0; jRow < numberRows_; jRow++) {
        if (markRow[jRow] == 0 && knapsackRow != jRow) {
            if (rowLower[jRow] > -1.0e20)
                rowLower[jRow] -= 1.0e-7;
            if (rowUpper[jRow] < 1.0e20)
                rowUpper[jRow] += 1.0e-7;
        } else {
            rowLower[jRow] = -COIN_DBL_MAX;
            rowUpper[jRow] = COIN_DBL_MAX;
        }
    }
    double * rowActivity = smallModel.primalRowSolution();
    CoinZeroN(rowActivity, numberRows_);
    maxSize -= offset;
    minSize -= offset;
    // now generate
    int i;
    int iStack = numJ;
    for (i = 0; i < numJ; i++) {
        stack[i] = 0;
    }
    double tooMuch = 10.0 * maxSize;
    stack[numJ] = 1;
    size[numJ] = tooMuch;
    bound[numJ] = 0;
    double sum = tooMuch;
    numberOutput = 0;
    int nelCreate = 0;
    /* typeRun is - 0 for initial sizes
                    1 for build
      	  2 for reconstruct
    */
    int typeRun = buildObj ? 1 : 0;
    if (reConstruct >= 0) {
        assert (buildRow && buildElement);
        typeRun = 2;
    }
    if (typeRun == 1)
        buildStart[0] = 0;
    while (iStack >= 0) {
        if (sum >= minSize && sum <= maxSize) {
            double checkSize = 0.0;
            bool good = true;
            int nRow = 0;
            double obj = objectiveOffset;
            // nRow is zero? CoinZeroN(rowActivity,nRow);
            for (iColumn = 0; iColumn < numJ; iColumn++) {
                int iValue = stack[iColumn];
                if (iValue > bound[iColumn]) {
                    good = false;
                    break;
                } else if (iValue) {
                    obj += objective[iColumn] * iValue;
                    for (CoinBigIndex j = columnStart[iColumn];
                            j < columnStart[iColumn] + columnLength[iColumn]; j++) {
                        double value = element[j] * iValue;
                        int kRow = row[j];
                        if (rowActivity[kRow]) {
                            rowActivity[kRow] += value;
                            if (!rowActivity[kRow])
                                rowActivity[kRow] = 1.0e-100;
                        } else {
                            build[nRow++] = kRow;
                            rowActivity[kRow] = value;
                        }
                    }
                }
            }
            if (good) {
                for (jRow = 0; jRow < nRow; jRow++) {
                    int kRow = build[jRow];
                    double value = rowActivity[kRow];
                    if (value > rowUpper[kRow] || value < rowLower[kRow]) {
                        good = false;
                        break;
                    }
                }
            }
            if (good) {
                if (typeRun == 1) {
                    buildObj[numberOutput] = obj;
                    for (jRow = 0; jRow < nRow; jRow++) {
                        int kRow = build[jRow];
                        double value = rowActivity[kRow];
                        if (markRow[kRow] < 0 && fabs(value) > 1.0e-13) {
                            buildElement[nelCreate] = value;
                            buildRow[nelCreate++] = kRow;
                        }
                    }
                    buildStart[numberOutput+1] = nelCreate;
                } else if (!typeRun) {
                    for (jRow = 0; jRow < nRow; jRow++) {
                        int kRow = build[jRow];
                        double value = rowActivity[kRow];
                        if (markRow[kRow] < 0 && fabs(value) > 1.0e-13) {
                            nelCreate++;
                        }
                    }
                }
                if (typeRun == 2 && reConstruct == numberOutput) {
                    // build and exit
                    nelCreate = 0;
                    for (iColumn = 0; iColumn < numJ; iColumn++) {
                        int iValue = stack[iColumn];
                        if (iValue) {
                            buildRow[nelCreate] = whichColumn[iColumn];
                            buildElement[nelCreate++] = iValue;
                        }
                    }
                    numberOutput = 1;
                    for (i = 0; i < numJ; i++) {
                        bound[i] = 0;
                    }
                    break;
                }
                numberOutput++;
                if (numberOutput > maxNumber) {
                    nelCreate = -1;
                    numberOutput = -1;
                    for (i = 0; i < numJ; i++) {
                        bound[i] = 0;
                    }
                    break;
                } else if (typeRun == 1 && numberOutput == maxNumber) {
                    // On second run
                    for (i = 0; i < numJ; i++) {
                        bound[i] = 0;
                    }
                    break;
                }
                for (int j = 0; j < numJ; j++) {
                    checkSize += stack[j] * size[j];
                }
                assert (fabs(sum - checkSize) < 1.0e-3);
            }
            for (jRow = 0; jRow < nRow; jRow++) {
                int kRow = build[jRow];
                rowActivity[kRow] = 0.0;
            }
        }
        if (sum > maxSize || stack[iStack] > bound[iStack]) {
            sum -= size[iStack] * stack[iStack];
            stack[iStack--] = 0;
            if (iStack >= 0) {
                stack[iStack] ++;
                sum += size[iStack];
            }
        } else {
            // must be less
            // add to last possible
            iStack = numJ - 1;
            sum += size[iStack];
            stack[iStack]++;
        }
    }
    //printf("%d will be created\n",numberOutput);
    delete [] whichColumn;
    delete [] whichRow;
    delete [] bound;
    delete [] stack;
    delete [] size;
    delete [] rhsOffset;
    delete [] build;
    delete [] markRow;
    return nelCreate;
}
#include "ClpConstraint.hpp"
#include "ClpConstraintLinear.hpp"
#include "ClpConstraintQuadratic.hpp"
#ifdef COIN_HAS_ASL
//#include "ClpAmplObjective.hpp"
#endif
/* Return an approximate solution to a CoinModel.
    Lots of bounds may be odd to force a solution.
    mode = 0 just tries to get a continuous solution
*/
ClpSimplex *
approximateSolution(CoinModel & coinModel,
                    int numberPasses, double deltaTolerance,
                    int /*mode*/)
{
#ifndef JJF_ONE
    //#ifdef COIN_HAS_ASL
    // matrix etc will be changed
    CoinModel coinModel2 = coinModel;
    if (coinModel2.moreInfo()) {
        // for now just ampl objective
        ClpSimplex * model = new ClpSimplex();
        model->loadProblem(coinModel2);
        int numberConstraints;
        ClpConstraint ** constraints = NULL;
        int type = model->loadNonLinear(coinModel2.moreInfo(),
                                        numberConstraints, constraints);
        if (type == 1 || type == 3) {
            model->nonlinearSLP(numberPasses, deltaTolerance);
        } else if (type == 2 || type == 4) {
            model->nonlinearSLP(numberConstraints, constraints,
                                numberPasses, deltaTolerance);
        } else {
            printf("error or linear - fix %d\n", type);
        }
        //exit(66);
        return model;
    }
    // first check and set up arrays
    int numberColumns = coinModel.numberColumns();
    int numberRows = coinModel.numberRows();
    // List of nonlinear rows
    int * which = new int[numberRows];
    bool testLinear = false;
    int numberConstraints = 0;
    int iColumn;
    bool linearObjective = true;
    int maximumQuadraticElements = 0;
    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
        // See if quadratic objective
        const char * expr = coinModel.getColumnObjectiveAsString(iColumn);
        if (strcmp(expr, "Numeric")) {
            linearObjective = false;
            // check if value*x+-value*y....
            assert (strlen(expr) < 20000);
            char temp[20000];
            strcpy(temp, expr);
            char * pos = temp;
            bool ifFirst = true;
            while (*pos) {
                double value;
                int jColumn = decodeBit(pos, pos, value, ifFirst, coinModel);
                // must be column unless first when may be linear term
                if (jColumn >= 0) {
                    maximumQuadraticElements++;
                } else if (jColumn != -2) {
                    printf("bad nonlinear term %s\n", temp);
                    abort();
                }
                ifFirst = false;
            }
        }
    }
    if (!linearObjective) {
        // zero objective
        for (iColumn = 0; iColumn < numberColumns; iColumn++)
            coinModel2.setObjective(iColumn, 0.0);
    }
    int iRow;
    for (iRow = 0; iRow < numberRows; iRow++) {
        int numberQuadratic = 0;
        bool linear = true;
        CoinModelLink triple = coinModel.firstInRow(iRow);
        while (triple.column() >= 0) {
            int iColumn = triple.column();
            const char *  expr = coinModel.getElementAsString(iRow, iColumn);
            if (strcmp("Numeric", expr)) {
                linear = false;
                // check if value*x+-value*y....
                assert (strlen(expr) < 20000);
                char temp[20000];
                strcpy(temp, expr);
                char * pos = temp;
                bool ifFirst = true;
                while (*pos) {
                    double value;
                    int jColumn = decodeBit(pos, pos, value, ifFirst, coinModel);
                    // must be column unless first when may be linear term
                    if (jColumn >= 0) {
                        numberQuadratic++;
                    } else if (jColumn != -2) {
                        printf("bad nonlinear term %s\n", temp);
                        abort();
                    }
                    ifFirst = false;
                }
            }
            triple = coinModel.next(triple);
        }
        if (!linear || testLinear) {
            CoinModelLink triple = coinModel.firstInRow(iRow);
            while (triple.column() >= 0) {
                int iColumn = triple.column();
                coinModel2.setElement(iRow, iColumn, 0.0);
                triple = coinModel.next(triple);
            }
            which[numberConstraints++] = iRow;
            maximumQuadraticElements = CoinMax(maximumQuadraticElements, numberQuadratic);
        }
    }
    ClpSimplex * model = new ClpSimplex();
    // return if nothing
    if (!numberConstraints && linearObjective) {
        delete [] which;
        model->loadProblem(coinModel);
        model->dual();
        return model;
    }
    // space for quadratic
    // allow for linear term
    maximumQuadraticElements += numberColumns;
    CoinBigIndex * startQuadratic = new CoinBigIndex [numberColumns+1];
    int * columnQuadratic = new int [maximumQuadraticElements];
    double * elementQuadratic = new double [maximumQuadraticElements];
    ClpConstraint ** constraints = new ClpConstraint * [numberConstraints];
    double * linearTerm = new double [numberColumns];
    int saveNumber = numberConstraints;
    numberConstraints = 0;
    ClpQuadraticObjective * quadObj = NULL;
    if (!linearObjective) {
        int numberQuadratic = 0;
        CoinZeroN(linearTerm, numberColumns);
        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
            startQuadratic[iColumn] = numberQuadratic;
            // See if quadratic objective
            const char * expr = coinModel.getColumnObjectiveAsString(iColumn);
            if (strcmp(expr, "Numeric")) {
                // value*x*y
                char temp[20000];
                strcpy(temp, expr);
                char * pos = temp;
                bool ifFirst = true;
                while (*pos) {
                    double value;
                    int jColumn = decodeBit(pos, pos, value, ifFirst, coinModel);
                    // must be column unless first when may be linear term
                    if (jColumn >= 0) {
                        columnQuadratic[numberQuadratic] = jColumn;
                        if (jColumn != iColumn)
                            elementQuadratic[numberQuadratic++] = 1.0 * value; // convention
                        else if (jColumn == iColumn)
                            elementQuadratic[numberQuadratic++] = 2.0 * value; // convention
                    } else if (jColumn == -2) {
                        linearTerm[iColumn] = value;
                    } else {
                        printf("bad nonlinear term %s\n", temp);
                        abort();
                    }
                    ifFirst = false;
                }
            } else {
                // linear part
                linearTerm[iColumn] = coinModel.getColumnObjective(iColumn);
            }
        }
        startQuadratic[numberColumns] = numberQuadratic;
        quadObj = new ClpQuadraticObjective(linearTerm, numberColumns,
                                            startQuadratic, columnQuadratic, elementQuadratic);
    }
    int iConstraint;
    for (iConstraint = 0; iConstraint < saveNumber; iConstraint++) {
        iRow = which[iConstraint];
        if (iRow >= 0) {
            int numberQuadratic = 0;
            int lastColumn = -1;
            int largestColumn = -1;
            CoinZeroN(linearTerm, numberColumns);
            CoinModelLink triple = coinModel.firstInRow(iRow);
            while (triple.column() >= 0) {
                int iColumn = triple.column();
                while (lastColumn < iColumn) {
                    startQuadratic[lastColumn+1] = numberQuadratic;
                    lastColumn++;
                }
                const char *  expr = coinModel.getElementAsString(iRow, iColumn);
                if (strcmp("Numeric", expr)) {
                    largestColumn = CoinMax(largestColumn, iColumn);
                    // value*x*y
                    char temp[20000];
                    strcpy(temp, expr);
                    char * pos = temp;
                    bool ifFirst = true;
                    while (*pos) {
                        double value;
                        int jColumn = decodeBit(pos, pos, value, ifFirst, coinModel);
                        // must be column unless first when may be linear term
                        if (jColumn >= 0) {
                            columnQuadratic[numberQuadratic] = jColumn;
                            if (jColumn == iColumn)
                                elementQuadratic[numberQuadratic++] = 2.0 * value; // convention
                            else
                                elementQuadratic[numberQuadratic++] = 1.0 * value; // convention
                            largestColumn = CoinMax(largestColumn, jColumn);
                        } else if (jColumn == -2) {
                            linearTerm[iColumn] = value;
                            // and put in as row -1
                            columnQuadratic[numberQuadratic] = -1;
                            if (jColumn == iColumn)
                                elementQuadratic[numberQuadratic++] = 2.0 * value; // convention
                            else
                                elementQuadratic[numberQuadratic++] = 1.0 * value; // convention
                            largestColumn = CoinMax(largestColumn, iColumn);
                        } else {
                            printf("bad nonlinear term %s\n", temp);
                            abort();
                        }
                        ifFirst = false;
                    }
                } else {
                    // linear part
                    linearTerm[iColumn] = coinModel.getElement(iRow, iColumn);
                    // and put in as row -1
                    columnQuadratic[numberQuadratic] = -1;
                    elementQuadratic[numberQuadratic++] = linearTerm[iColumn];
                    if (linearTerm[iColumn])
                        largestColumn = CoinMax(largestColumn, iColumn);
                }
                triple = coinModel.next(triple);
            }
            while (lastColumn < numberColumns) {
                startQuadratic[lastColumn+1] = numberQuadratic;
                lastColumn++;
            }
            // here we create ClpConstraint
            if (testLinear) {
                int n = 0;
                int * indices = new int[numberColumns];
                for (int j = 0; j < numberColumns; j++) {
                    if (linearTerm[j]) {
                        linearTerm[n] = linearTerm[j];
                        indices[n++] = j;
                    }
                }
                /// Constructor from constraint
                constraints[numberConstraints++] = new ClpConstraintLinear(iRow, n, numberColumns,
                        indices, linearTerm);
                delete [] indices;
            } else {
                constraints[numberConstraints++] = new ClpConstraintQuadratic(iRow, largestColumn + 1, numberColumns,
                        startQuadratic, columnQuadratic, elementQuadratic);
            }
        }
    }
    delete [] startQuadratic;
    delete [] columnQuadratic;
    delete [] elementQuadratic;
    delete [] linearTerm;
    delete [] which;
    model->loadProblem(coinModel2);
    if (quadObj)
        model->setObjective(quadObj);
    delete quadObj;
#ifndef NDEBUG
    int returnCode;
    if (numberConstraints) {
        returnCode = model->nonlinearSLP(numberConstraints, constraints,
                                         numberPasses, deltaTolerance);
        for (iConstraint = 0; iConstraint < saveNumber; iConstraint++)
            delete constraints[iConstraint];
    } else {
        returnCode = model->nonlinearSLP(numberPasses, deltaTolerance);
    }
    assert (!returnCode);
#else
    if (numberConstraints) {
        model->nonlinearSLP(numberConstraints, constraints, numberPasses, deltaTolerance);
        for (iConstraint = 0; iConstraint < saveNumber; iConstraint++)
            delete constraints[iConstraint];
    } else {
        model->nonlinearSLP(numberPasses, deltaTolerance);
    }
#endif
    delete [] constraints;
    return model;
#else
    printf("loadNonLinear needs ampl\n");
    abort();
    return NULL;
#endif
}
OsiChooseStrongSubset::OsiChooseStrongSubset() :
        OsiChooseStrong(),
        numberObjectsToUse_(0)
{
}

OsiChooseStrongSubset::OsiChooseStrongSubset(const OsiSolverInterface * solver) :
        OsiChooseStrong(solver),
        numberObjectsToUse_(-1)
{
}

OsiChooseStrongSubset::OsiChooseStrongSubset(const OsiChooseStrongSubset & rhs)
        : OsiChooseStrong(rhs)
{
    numberObjectsToUse_ = -1;
}

OsiChooseStrongSubset &
OsiChooseStrongSubset::operator=(const OsiChooseStrongSubset & rhs)
{
    if (this != &rhs) {
        OsiChooseStrong::operator=(rhs);
        numberObjectsToUse_ = -1;
    }
    return *this;
}


OsiChooseStrongSubset::~OsiChooseStrongSubset ()
{
}

// Clone
OsiChooseVariable *
OsiChooseStrongSubset::clone() const
{
    return new OsiChooseStrongSubset(*this);
}
// Initialize
int
OsiChooseStrongSubset::setupList ( OsiBranchingInformation *info, bool initialize)
{
    assert (solver_ == info->solver_);
    // Only has to work with Clp
    OsiSolverInterface * solverA = const_cast<OsiSolverInterface *> (solver_);
    OsiSolverLink * solver = dynamic_cast<OsiSolverLink *> (solverA);
    assert (solver);
    int numberObjects = solver->numberObjects();
    if (numberObjects > pseudoCosts_.numberObjects()) {
        // redo useful arrays
        pseudoCosts_.initialize(numberObjects);
    }
    int numObj = numberObjects;
    if (numberObjectsToUse_ < 0) {
        // Sort objects so bilinear at end
        OsiObject ** sorted = new OsiObject * [numberObjects];
        OsiObject ** objects = solver->objects();
        numObj = 0;
        int numberBiLinear = 0;
        int i;
        for (i = 0; i < numberObjects; i++) {
            OsiObject * obj = objects[i];
            OsiBiLinear * objB = dynamic_cast<OsiBiLinear *> (obj);
            if (!objB)
                objects[numObj++] = obj;
            else
                sorted[numberBiLinear++] = obj;
        }
        numberObjectsToUse_ = numObj;
        for (i = 0; i < numberBiLinear; i++)
            objects[numObj++] = sorted[i];
        delete [] sorted;
        // See if any master objects
        for (i = 0; i < numberObjectsToUse_; i++) {
            OsiUsesBiLinear * obj = dynamic_cast<OsiUsesBiLinear *> (objects[i]);
            if (obj)
                obj->addBiLinearObjects(solver);
        }
    }
    solver->setNumberObjects(numberObjectsToUse_);
    numObj = numberObjectsToUse_;
    // Use shadow prices
    //info->defaultDual_=0.0;
    int numberUnsatisfied = OsiChooseStrong::setupList ( info, initialize);
    solver->setNumberObjects(numberObjects);
    numObj = numberObjects;
    return numberUnsatisfied;
}
/* Choose a variable
   Returns -
   -1 Node is infeasible
   0  Normal termination - we have a candidate
   1  All looks satisfied - no candidate
   2  We can change the bound on a variable - but we also have a strong branching candidate
   3  We can change the bound on a variable - but we have a non-strong branching candidate
   4  We can change the bound on a variable - no other candidates
   We can pick up branch from whichObject() and whichWay()
   We can pick up a forced branch (can change bound) from whichForcedObject() and whichForcedWay()
   If we have a solution then we can pick up from goodObjectiveValue() and goodSolution()
*/
int
OsiChooseStrongSubset::chooseVariable( OsiSolverInterface * solver, OsiBranchingInformation *info, bool fixVariables)
{
    //int numberObjects = solver->numberObjects();
    //solver->setNumberObjects(numberObjectsToUse_);
    //numberObjects_=numberObjectsToUse_;
    // Use shadow prices
    //info->defaultDual_=0.0;
    int returnCode = OsiChooseStrong::chooseVariable(solver, info, fixVariables);
    //solver->setNumberObjects(numberObjects);
    //numberObjects_=numberObjects;
    return returnCode;
}
/** Default Constructor

  Equivalent to an unspecified binary variable.
*/
OsiUsesBiLinear::OsiUsesBiLinear ()
        : OsiSimpleInteger(),
        numberBiLinear_(0),
        type_(0),
        objects_(NULL)
{
}

/** Useful constructor

  Loads actual upper & lower bounds for the specified variable.
*/
OsiUsesBiLinear::OsiUsesBiLinear (const OsiSolverInterface * solver, int iColumn, int type)
        : OsiSimpleInteger(solver, iColumn),
        numberBiLinear_(0),
        type_(type),
        objects_(NULL)
{
    if (type_) {
        assert(originalLower_ == floor(originalLower_ + 0.5));
        assert(originalUpper_ == floor(originalUpper_ + 0.5));
    }
}


// Useful constructor - passed solver index and original bounds
OsiUsesBiLinear::OsiUsesBiLinear ( int iColumn, double lower, double upper, int type)
        : OsiSimpleInteger(iColumn, lower, upper),
        numberBiLinear_(0),
        type_(type),
        objects_(NULL)
{
    if (type_) {
        assert(originalLower_ == floor(originalLower_ + 0.5));
        assert(originalUpper_ == floor(originalUpper_ + 0.5));
    }
}

// Useful constructor - passed simple integer
OsiUsesBiLinear::OsiUsesBiLinear ( const OsiSimpleInteger &rhs, int type)
        : OsiSimpleInteger(rhs),
        numberBiLinear_(0),
        type_(type),
        objects_(NULL)
{
    if (type_) {
        assert(originalLower_ == floor(originalLower_ + 0.5));
        assert(originalUpper_ == floor(originalUpper_ + 0.5));
    }
}

// Copy constructor
OsiUsesBiLinear::OsiUsesBiLinear ( const OsiUsesBiLinear & rhs)
        : OsiSimpleInteger(rhs),
        numberBiLinear_(0),
        type_(rhs.type_),
        objects_(NULL)

{
}

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

// Assignment operator
OsiUsesBiLinear &
OsiUsesBiLinear::operator=( const OsiUsesBiLinear & rhs)
{
    if (this != &rhs) {
        OsiSimpleInteger::operator=(rhs);
        delete [] objects_;
        numberBiLinear_ = 0;
        type_ = rhs.type_;
        objects_ = NULL;
    }
    return *this;
}

// Destructor
OsiUsesBiLinear::~OsiUsesBiLinear ()
{
    delete [] objects_;
}
// Infeasibility - large is 0.5
double
OsiUsesBiLinear::infeasibility(const OsiBranchingInformation * info, int & whichWay) const
{
    assert (type_ == 0); // just continuous for now
    double value = info->solution_[columnNumber_];
    value = CoinMax(value, info->lower_[columnNumber_]);
    value = CoinMin(value, info->upper_[columnNumber_]);
    infeasibility_ = 0.0;
    for (int i = 0; i < numberBiLinear_; i++) {
        OsiBiLinear * obj = dynamic_cast<OsiBiLinear *> (objects_[i]);
        assert (obj);
        //obj->getPseudoShadow(info);
        //infeasibility_ += objects_[i]->infeasibility(info,whichWay);
        infeasibility_ += obj->getMovement(info);
    }
    bool satisfied = false;
    whichWay = -1;
    if (!infeasibility_) {
        otherInfeasibility_ = 1.0;
        satisfied = true;
        infeasibility_ = 0.0;
    } else {
        otherInfeasibility_ = 10.0 * infeasibility_;
        if (value - info->lower_[columnNumber_] >
                info->upper_[columnNumber_] - value)
            whichWay = 1;
        else
            whichWay = -1;
    }
    if (preferredWay_ >= 0 && !satisfied)
        whichWay = preferredWay_;
    whichWay_ = static_cast<short int>(whichWay);
    return infeasibility_;
}
// Creates a branching object
OsiBranchingObject *
OsiUsesBiLinear::createBranch(OsiSolverInterface * solver, const OsiBranchingInformation * info, int way) const
{
    double value = info->solution_[columnNumber_];
    value = CoinMax(value, info->lower_[columnNumber_]);
    value = CoinMin(value, info->upper_[columnNumber_]);
    assert (info->upper_[columnNumber_] > info->lower_[columnNumber_]);
    double nearest = floor(value + 0.5);
    double integerTolerance = info->integerTolerance_;
    if (fabs(value - nearest) < integerTolerance) {
        // adjust value
        if (nearest != info->upper_[columnNumber_])
            value = nearest + 2.0 * integerTolerance;
        else
            value = nearest - 2.0 * integerTolerance;
    }
    OsiBranchingObject * branch = new OsiIntegerBranchingObject(solver, this, way,
            value, value, value);
    return branch;
}
// This looks at solution and sets bounds to contain solution
/** More precisely: it first forces the variable within the existing
    bounds, and then tightens the bounds to fix the variable at the
    nearest integer value.
*/
double
OsiUsesBiLinear::feasibleRegion(OsiSolverInterface * solver,
                                const OsiBranchingInformation * info) const
{
    double value = info->solution_[columnNumber_];
    double newValue = CoinMax(value, info->lower_[columnNumber_]);
    newValue = CoinMin(newValue, info->upper_[columnNumber_]);
    solver->setColLower(columnNumber_, newValue);
    solver->setColUpper(columnNumber_, newValue);
    return fabs(value - newValue);
}
// Add all bi-linear objects
void
OsiUsesBiLinear::addBiLinearObjects(OsiSolverLink * solver)
{
    delete [] objects_;
    numberBiLinear_ = 0;
    OsiObject ** objects = solver->objects();
    int i;
    int numberObjects = solver->numberObjects();
    for (i = 0; i < numberObjects; i++) {
        OsiObject * obj = objects[i];
        OsiBiLinear * objB = dynamic_cast<OsiBiLinear *> (obj);
        if (objB) {
            if (objB->xColumn() == columnNumber_ || objB->yColumn() == columnNumber_)
                numberBiLinear_++;
        }
    }
    if (numberBiLinear_) {
        objects_ = new OsiObject * [numberBiLinear_];
        numberBiLinear_ = 0;
        for (i = 0; i < numberObjects; i++) {
            OsiObject * obj = objects[i];
            OsiBiLinear * objB = dynamic_cast<OsiBiLinear *> (obj);
            if (objB) {
                if (objB->xColumn() == columnNumber_ || objB->yColumn() == columnNumber_)
                    objects_[numberBiLinear_++] = obj;;
            }
        }
    } else {
        objects_ = NULL;
    }
}