source: trunk/Cbc/src/CbcModel.cpp @ 2464

/* $Id: CbcModel.cpp 2464 2019-01-03 19:03:23Z unxusr $ */
// Copyright (C) 2002, International Business Machines
// Corporation and others.  All Rights Reserved.
// This code is licensed under the terms of the Eclipse Public License (EPL).

#if defined(_MSC_VER)
// Turn off compiler warning about long names
#pragma warning(disable : 4786)
#endif

#include "CbcConfig.h"

#include <string>
//#define CBC_DEBUG 1
//#define CHECK_CUT_COUNTS
//#define CHECK_NODE
//#define CHECK_NODE_FULL
//#define NODE_LOG
//#define GLOBAL_CUTS_JUST_POINTERS
#ifdef CGL_DEBUG_GOMORY
extern int gomory_try;
#endif
#include <cassert>
#include <cmath>
#include <cfloat>
#ifdef COIN_HAS_CLP
// include Presolve from Clp
#include "ClpPresolve.hpp"
#include "OsiClpSolverInterface.hpp"
#include "ClpNode.hpp"
#include "ClpDualRowDantzig.hpp"
#include "ClpSimplexPrimal.hpp"
#endif

#include "CbcEventHandler.hpp"

#include "OsiSolverInterface.hpp"
#include "OsiAuxInfo.hpp"
#include "OsiSolverBranch.hpp"
#include "OsiChooseVariable.hpp"
#include "CoinWarmStartBasis.hpp"
#include "CoinPackedMatrix.hpp"
#include "CoinHelperFunctions.hpp"
#include "CbcBranchActual.hpp"
#include "CbcBranchDynamic.hpp"
#include "CbcHeuristic.hpp"
#include "CbcHeuristicFPump.hpp"
#include "CbcHeuristicRINS.hpp"
#include "CbcHeuristicDive.hpp"
#include "CbcModel.hpp"
#include "CbcTreeLocal.hpp"
#include "CbcStatistics.hpp"
#include "CbcStrategy.hpp"
#include "CbcMessage.hpp"
#include "OsiRowCut.hpp"
#include "OsiColCut.hpp"
#include "OsiRowCutDebugger.hpp"
#include "OsiCuts.hpp"
#include "CbcCountRowCut.hpp"
#include "CbcCutGenerator.hpp"
#include "CbcFeasibilityBase.hpp"
#include "CbcFathom.hpp"
#include "CbcFullNodeInfo.hpp"
#ifdef COIN_HAS_NTY
#include "CbcSymmetry.hpp"
#endif
// include Probing
#include "CglProbing.hpp"
#include "CglGomory.hpp"
#include "CglTwomir.hpp"
// include preprocessing
#include "CglPreProcess.hpp"
#include "CglDuplicateRow.hpp"
#include "CglStored.hpp"
#include "CglClique.hpp"
#include "CglKnapsackCover.hpp"

#include "CoinTime.hpp"
#include "CoinMpsIO.hpp"

#include "CbcCompareActual.hpp"
#include "CbcTree.hpp"
// This may be dummy
#include "CbcThread.hpp"
/* Various functions local to CbcModel.cpp */

typedef struct {
  double useCutoff;
  CbcModel *model;
  int switches;
} rootBundle;
static void *doRootCbcThread(void *voidInfo);

namespace {

//-------------------------------------------------------------------
// Returns the greatest common denominator of two
// positive integers, a and b, found using Euclid's algorithm
//-------------------------------------------------------------------
static int gcd(int a, int b)
{
  int remainder = -1;
  // make sure a<=b (will always remain so)
  if (a > b) {
    // Swap a and b
    int temp = a;
    a = b;
    b = temp;
  }
  // if zero then gcd is nonzero (zero may occur in rhs of packed)
  if (!a) {
    if (b) {
      return b;
    } else {
      printf("**** gcd given two zeros!!\n");
      abort();
    }
  }
  while (remainder) {
    remainder = b % a;
    b = a;
    a = remainder;
  }
  return b;
}

#ifdef CHECK_NODE_FULL

/*
  Routine to verify that tree linkage is correct. The invariant that is tested
  is

  reference count = (number of actual references) + (number of branches left)

  The routine builds a set of paired arrays, info and count, by traversing the
  tree. Each CbcNodeInfo is recorded in info, and the number of times it is
  referenced (via the parent field) is recorded in count. Then a final check is
  made to see if the numberPointingToThis_ field agrees.
*/

void verifyTreeNodes(const CbcTree *branchingTree, const CbcModel &model)

{
  if (model.getNodeCount() == 661)
    return;
  printf("*** CHECKING tree after %d nodes\n", model.getNodeCount());

  int j;
  int nNodes = branchingTree->size();
#define MAXINFO 1000
  int *count = new int[MAXINFO];
  CbcNodeInfo **info = new CbcNodeInfo *[MAXINFO];
  int nInfo = 0;
  /*
      Collect all CbcNodeInfo objects in info, by starting from each live node and
      traversing back to the root. Nodes in the live set should have unexplored
      branches remaining.

      TODO: The `while (nodeInfo)' loop could be made to break on reaching a
        common ancester (nodeInfo is found in info[k]). Alternatively, the
        check could change to signal an error if nodeInfo is not found above a
        common ancestor.
    */
  for (j = 0; j < nNodes; j++) {
    CbcNode *node = branchingTree->nodePointer(j);
    if (!node)
      continue;
    CbcNodeInfo *nodeInfo = node->nodeInfo();
    int change = node->nodeInfo()->numberBranchesLeft();
    assert(change);
    while (nodeInfo) {
      int k;
      for (k = 0; k < nInfo; k++) {
        if (nodeInfo == info[k])
          break;
      }
      if (k == nInfo) {
        assert(nInfo < MAXINFO);
        nInfo++;
        info[k] = nodeInfo;
        count[k] = 0;
      }
      nodeInfo = nodeInfo->parent();
    }
  }
  /*
      Walk the info array. For each nodeInfo, look up its parent in info and
      increment the corresponding count.
    */
  for (j = 0; j < nInfo; j++) {
    CbcNodeInfo *nodeInfo = info[j];
    nodeInfo = nodeInfo->parent();
    if (nodeInfo) {
      int k;
      for (k = 0; k < nInfo; k++) {
        if (nodeInfo == info[k])
          break;
      }
      assert(k < nInfo);
      count[k]++;
    }
  }
  /*
      Walk the info array one more time and check that the invariant holds. The
      number of references (numberPointingToThis()) should equal the sum of the
      number of actual references (held in count[]) plus the number of potential
      references (unexplored branches, numberBranchesLeft()).
    */
  for (j = 0; j < nInfo; j++) {
    CbcNodeInfo *nodeInfo = info[j];
    if (nodeInfo) {
      int k;
      for (k = 0; k < nInfo; k++)
        if (nodeInfo == info[k])
          break;
      printf("Nodeinfo %x - %d left, %d count\n",
        nodeInfo,
        nodeInfo->numberBranchesLeft(),
        nodeInfo->numberPointingToThis());
      assert(nodeInfo->numberPointingToThis() == count[k] + nodeInfo->numberBranchesLeft());
    }
  }

  delete[] count;
  delete[] info;

  return;
}

#endif /* CHECK_NODE_FULL */

#ifdef CHECK_CUT_COUNTS

/*
  Routine to verify that cut reference counts are correct.
*/
void verifyCutCounts(const CbcTree *branchingTree, CbcModel &model)

{
  printf("*** CHECKING cuts after %d nodes\n", model.getNodeCount());

  int j;
  int nNodes = branchingTree->size();

  /*
      cut.tempNumber_ exists for the purpose of doing this verification. Clear it
      in all cuts. We traverse the tree by starting from each live node and working
      back to the root. At each CbcNodeInfo, check for cuts.
    */
  for (j = 0; j < nNodes; j++) {
    CbcNode *node = branchingTree->nodePointer(j);
    CbcNodeInfo *nodeInfo = node->nodeInfo();
    assert(node->nodeInfo()->numberBranchesLeft());
    while (nodeInfo) {
      int k;
      for (k = 0; k < nodeInfo->numberCuts(); k++) {
        CbcCountRowCut *cut = nodeInfo->cuts()[k];
        if (cut)
          cut->tempNumber_ = 0;
      }
      nodeInfo = nodeInfo->parent();
    }
  }
  /*
      Walk the live set again, this time collecting the list of cuts in use at each
      node. addCuts1 will collect the cuts in model.addedCuts_. Take into account
      that when we recreate the basis for a node, we compress out the slack cuts.
    */
  for (j = 0; j < nNodes; j++) {
    CoinWarmStartBasis *debugws = model.getEmptyBasis();
    CbcNode *node = branchingTree->nodePointer(j);
    CbcNodeInfo *nodeInfo = node->nodeInfo();
    int change = node->nodeInfo()->numberBranchesLeft();
    printf("Node %d %x (info %x) var %d way %d obj %g", j, node,
      node->nodeInfo(), node->columnNumber(), node->way(),
      node->objectiveValue());

    model.addCuts1(node, debugws);

    int i;
    int numberRowsAtContinuous = model.numberRowsAtContinuous();
    CbcCountRowCut **addedCuts = model.addedCuts();
    for (i = 0; i < model.currentNumberCuts(); i++) {
      CoinWarmStartBasis::Status status = debugws->getArtifStatus(i + numberRowsAtContinuous);
      if (status != CoinWarmStartBasis::basic && addedCuts[i]) {
        addedCuts[i]->tempNumber_ += change;
      }
    }

    while (nodeInfo) {
      nodeInfo = nodeInfo->parent();
      if (nodeInfo)
        printf(" -> %x", nodeInfo);
    }
    printf("\n");
    delete debugws;
  }
  /*
      The moment of truth: We've tallied up the references by direct scan of the  search tree. Check for agreement with the count in the cut.

      TODO: Rewrite to check and print mismatch only when tempNumber_ == 0?
    */
  for (j = 0; j < nNodes; j++) {
    CbcNode *node = branchingTree->nodePointer(j);
    CbcNodeInfo *nodeInfo = node->nodeInfo();
    while (nodeInfo) {
      int k;
      for (k = 0; k < nodeInfo->numberCuts(); k++) {
        CbcCountRowCut *cut = nodeInfo->cuts()[k];
        if (cut && cut->tempNumber_ >= 0) {
          if (cut->tempNumber_ != cut->numberPointingToThis())
            printf("mismatch %x %d %x %d %d\n", nodeInfo, k,
              cut, cut->tempNumber_, cut->numberPointingToThis());
          else
            printf("   match %x %d %x %d %d\n", nodeInfo, k,
              cut, cut->tempNumber_, cut->numberPointingToThis());
          cut->tempNumber_ = -1;
        }
      }
      nodeInfo = nodeInfo->parent();
    }
  }

  return;
}

#endif /* CHECK_CUT_COUNTS */

#ifdef CHECK_CUT_SIZE

/*
  Routine to verify that cut reference counts are correct.
*/
void verifyCutSize(const CbcTree *branchingTree, CbcModel &model)
{

  int j;
  int nNodes = branchingTree->size();
  int totalCuts = 0;

  /*
      cut.tempNumber_ exists for the purpose of doing this verification. Clear it
      in all cuts. We traverse the tree by starting from each live node and working
      back to the root. At each CbcNodeInfo, check for cuts.
    */
  for (j = 0; j < nNodes; j++) {
    CbcNode *node = branchingTree->nodePointer(j);
    CbcNodeInfo *nodeInfo = node->nodeInfo();
    assert(node->nodeInfo()->numberBranchesLeft());
    while (nodeInfo) {
      totalCuts += nodeInfo->numberCuts();
      nodeInfo = nodeInfo->parent();
    }
  }
  printf("*** CHECKING cuts (size) after %d nodes - %d cuts\n", model.getNodeCount(), totalCuts);
  return;
}

#endif /* CHECK_CUT_SIZE */

}

/* End unnamed namespace for CbcModel.cpp */

void CbcModel::analyzeObjective()
/*
  Try to find a minimum change in the objective function. The first scan
  checks that there are no continuous variables with non-zero coefficients,
  and grabs the largest objective coefficient associated with an unfixed
  integer variable. The second scan attempts to scale up the objective
  coefficients to a point where they are sufficiently close to integer that
  we can pretend they are integer, and calculate a gcd over the coefficients
  of interest. This will be the minimum increment for the scaled coefficients.
  The final action is to scale the increment back for the original coefficients
  and install it, if it's better than the existing value.

  John's note: We could do better than this.

  John's second note - apologies for changing s to z
*/
{
  const double *objective = getObjCoefficients();
  const double *lower = getColLower();
  const double *upper = getColUpper();
  /*
      Scan continuous and integer variables to see if continuous
      are cover or network with integral rhs.
    */
  double continuousMultiplier = 1.0;
  double *coeffMultiplier = NULL;
  double largestObj = 0.0;
  double smallestObj = COIN_DBL_MAX;
  {
    const double *rowLower = getRowLower();
    const double *rowUpper = getRowUpper();
    int numberRows = solver_->getNumRows();
    double *rhs = new double[numberRows];
    memset(rhs, 0, numberRows * sizeof(double));
    int iColumn;
    int numberColumns = solver_->getNumCols();
    // Column copy of matrix
    int problemType = -1;
    const double *element = solver_->getMatrixByCol()->getElements();
    const int *row = solver_->getMatrixByCol()->getIndices();
    const CoinBigIndex *columnStart = solver_->getMatrixByCol()->getVectorStarts();
    const int *columnLength = solver_->getMatrixByCol()->getVectorLengths();
    int numberInteger = 0;
    int numberIntegerObj = 0;
    int numberGeneralIntegerObj = 0;
    int numberIntegerWeight = 0;
    int numberContinuousObj = 0;
    double cost = COIN_DBL_MAX;
    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
      if (upper[iColumn] == lower[iColumn]) {
        CoinBigIndex start = columnStart[iColumn];
        CoinBigIndex end = start + columnLength[iColumn];
        for (CoinBigIndex j = start; j < end; j++) {
          int iRow = row[j];
          rhs[iRow] += lower[iColumn] * element[j];
        }
      } else {
        double objValue = objective[iColumn];
        if (solver_->isInteger(iColumn))
          numberInteger++;
        if (objValue) {
          if (!solver_->isInteger(iColumn)) {
            numberContinuousObj++;
          } else {
            largestObj = CoinMax(largestObj, fabs(objValue));
            smallestObj = CoinMin(smallestObj, fabs(objValue));
            numberIntegerObj++;
            if (cost == COIN_DBL_MAX)
              cost = objValue;
            else if (cost != objValue)
              cost = -COIN_DBL_MAX;
            int gap = static_cast<int>(upper[iColumn] - lower[iColumn]);
            if (gap > 1) {
              numberGeneralIntegerObj++;
              numberIntegerWeight += gap;
            }
          }
        }
      }
    }
    int iType = 0;
    if (!numberContinuousObj && numberIntegerObj <= 5 && numberIntegerWeight <= 100 && numberIntegerObj * 3 < numberObjects_ && !parentModel_ && solver_->getNumRows() > 100)
      iType = 1 + 4 + (((moreSpecialOptions_ & 536870912) == 0) ? 2 : 0);
    else if (!numberContinuousObj && numberIntegerObj <= 100 && numberIntegerObj * 5 < numberObjects_ && numberIntegerWeight <= 100 && !parentModel_ && solver_->getNumRows() > 100 && cost != -COIN_DBL_MAX)
      iType = 4 + (((moreSpecialOptions_ & 536870912) == 0) ? 2 : 0);
    else if (!numberContinuousObj && numberIntegerObj <= 100 && numberIntegerObj * 5 < numberObjects_ && !parentModel_ && solver_->getNumRows() > 100 && cost != -COIN_DBL_MAX)
      iType = 8;
    int iTest = getMaximumNodes();
    if (iTest >= 987654320 && iTest < 987654330 && numberObjects_ && !parentModel_) {
      iType = iTest - 987654320;
      printf("Testing %d integer variables out of %d objects (%d integer) have cost of %g - %d continuous\n",
        numberIntegerObj, numberObjects_, numberInteger, cost, numberContinuousObj);
      if (iType == 9)
        exit(77);
      if (numberContinuousObj)
        iType = 0;
    }

    //if (!numberContinuousObj&&(numberIntegerObj<=5||cost!=-COIN_DBL_MAX)&&
    //numberIntegerObj*3<numberObjects_&&!parentModel_&&solver_->getNumRows()>100) {
    if (iType) {
      /*
            A) put high priority on (if none)
            B) create artificial objective (if clp)
            */
      int iPriority = -1;
      for (int i = 0; i < numberObjects_; i++) {
        int k = object_[i]->priority();
        if (iPriority == -1)
          iPriority = k;
        else if (iPriority != k)
          iPriority = -2;
      }
      bool branchOnSatisfied = ((iType & 1) != 0);
      bool createFake = ((iType & 2) != 0);
      bool randomCost = ((iType & 4) != 0);
      if (iPriority >= 0) {
        char general[200];
        if (cost == -COIN_DBL_MAX) {
          sprintf(general, "%d integer variables out of %d objects (%d integer) have costs - high priority",
            numberIntegerObj, numberObjects_, numberInteger);
        } else if (cost == COIN_DBL_MAX) {
          sprintf(general, "No integer variables out of %d objects (%d integer) have costs",
            numberObjects_, numberInteger);
          branchOnSatisfied = false;
        } else {
          sprintf(general, "%d integer variables out of %d objects (%d integer) have cost of %g - high priority",
            numberIntegerObj, numberObjects_, numberInteger, cost);
        }
        messageHandler()->message(CBC_GENERAL,
          messages())
          << general << CoinMessageEol;
        sprintf(general, "branch on satisfied %c create fake objective %c random cost %c",
          branchOnSatisfied ? 'Y' : 'N',
          createFake ? 'Y' : 'N',
          randomCost ? 'Y' : 'N');
        messageHandler()->message(CBC_GENERAL,
          messages())
          << general << CoinMessageEol;
        // switch off clp type branching
        // no ? fastNodeDepth_ = -1;
        int highPriority = (branchOnSatisfied) ? -999 : 100;
        for (int i = 0; i < numberObjects_; i++) {
          CbcSimpleInteger *thisOne = dynamic_cast<CbcSimpleInteger *>(object_[i]);
          object_[i]->setPriority(1000);
          if (thisOne) {
            int iColumn = thisOne->columnNumber();
            if (objective[iColumn])
              thisOne->setPriority(highPriority);
          }
        }
      }
#ifdef COIN_HAS_CLP
      OsiClpSolverInterface *clpSolver
        = dynamic_cast<OsiClpSolverInterface *>(solver_);
      if (clpSolver && createFake) {
        // Create artificial objective to be used when all else fixed
        int numberColumns = clpSolver->getNumCols();
        double *fakeObj = new double[numberColumns];
        // Column copy
        const CoinPackedMatrix *matrixByCol = clpSolver->getMatrixByCol();
        //const double * element = matrixByCol.getElements();
        //const int * row = matrixByCol.getIndices();
        //const CoinBigIndex * columnStart = matrixByCol.getVectorStarts();
        const int *columnLength = matrixByCol->getVectorLengths();
        const double *solution = clpSolver->getColSolution();
#ifdef JJF_ZERO
        int nAtBound = 0;
        for (int i = 0; i < numberColumns; i++) {
          double lowerValue = lower[i];
          double upperValue = upper[i];
          if (clpSolver->isInteger(i)) {
            double lowerValue = lower[i];
            double upperValue = upper[i];
            double value = solution[i];
            if (value < lowerValue + 1.0e-6 || value > upperValue - 1.0e-6)
              nAtBound++;
          }
        }
#endif
        /*
                  Generate a random objective function for problems where the given objective
                  function is not terribly useful. (Nearly feasible, single integer variable,
                  that sort of thing.
                */
        CoinDrand48(true, 1234567);
        for (int i = 0; i < numberColumns; i++) {
          double lowerValue = lower[i];
          double upperValue = upper[i];
          double value = (randomCost) ? ceil((CoinDrand48() + 0.5) * 1000)
                                      : i + 1 + columnLength[i] * 1000;
          value *= 0.001;
          //value += columnLength[i];
          if (lowerValue > -1.0e5 || upperValue < 1.0e5) {
            if (fabs(lowerValue) > fabs(upperValue))
              value = -value;
            if (clpSolver->isInteger(i)) {
              double solValue = solution[i];
              // Better to add in 0.5 or 1.0??
              if (solValue < lowerValue + 1.0e-6)
                value = fabs(value) + 0.5; //fabs(value*1.5);
              else if (solValue > upperValue - 1.0e-6)
                value = -fabs(value) - 0.5; //-fabs(value*1.5);
            }
          } else {
            value = 0.0;
          }
          fakeObj[i] = value;
        }
        // pass to solver
        clpSolver->setFakeObjective(fakeObj);
        delete[] fakeObj;
      }
#endif
    } else if (largestObj < smallestObj * 5.0 && !parentModel_ && !numberContinuousObj && !numberGeneralIntegerObj && numberIntegerObj * 2 < CoinMin(numberColumns, 20)) {
      // up priorities on costed
      int iPriority = -1;
      for (int i = 0; i < numberObjects_; i++) {
        int k = object_[i]->priority();
        if (iPriority == -1)
          iPriority = k;
        else if (iPriority != k)
          iPriority = -2;
      }
      if (iPriority >= 100) {
#if CBC_USEFUL_PRINTING > 1
        printf("Setting variables with obj to high priority\n");
#endif
        for (int i = 0; i < numberObjects_; i++) {
          CbcSimpleInteger *obj = dynamic_cast<CbcSimpleInteger *>(object_[i]);
          if (obj) {
            int iColumn = obj->columnNumber();
            if (objective[iColumn])
              object_[i]->setPriority(iPriority - 1);
          }
        }
      }
    }
    int iRow;
    for (iRow = 0; iRow < numberRows; iRow++) {
      if (rowLower[iRow] > -1.0e20 && fabs(rowLower[iRow] - rhs[iRow] - floor(rowLower[iRow] - rhs[iRow] + 0.5)) > 1.0e-10) {
        continuousMultiplier = 0.0;
        break;
      }
      if (rowUpper[iRow] < 1.0e20 && fabs(rowUpper[iRow] - rhs[iRow] - floor(rowUpper[iRow] - rhs[iRow] + 0.5)) > 1.0e-10) {
        continuousMultiplier = 0.0;
        break;
      }
      // set rhs to limiting value
      if (rowLower[iRow] != rowUpper[iRow]) {
        if (rowLower[iRow] > -1.0e20) {
          if (rowUpper[iRow] < 1.0e20) {
            // no good
            continuousMultiplier = 0.0;
            break;
          } else {
            rhs[iRow] = rowLower[iRow] - rhs[iRow];
            if (problemType < 0)
              problemType = 3; // set cover
            else if (problemType != 3)
              problemType = 4;
          }
        } else {
          rhs[iRow] = rowUpper[iRow] - rhs[iRow];
          if (problemType < 0)
            problemType = 1; // set partitioning <=
          else if (problemType != 1)
            problemType = 4;
        }
      } else {
        rhs[iRow] = rowUpper[iRow] - rhs[iRow];
        if (problemType < 0)
          problemType = 3; // set partitioning ==
        else if (problemType != 2)
          problemType = 2;
      }
      if (fabs(rhs[iRow] - 1.0) > 1.0e-12)
        problemType = 4;
    }
    if (continuousMultiplier) {
      // 1 network, 2 cover, 4 negative cover
      int possible = 7;
      bool unitRhs = true;
      // See which rows could be set cover
      for (iColumn = 0; iColumn < numberColumns; iColumn++) {
        if (upper[iColumn] > lower[iColumn] + 1.0e-8) {
          CoinBigIndex start = columnStart[iColumn];
          CoinBigIndex end = start + columnLength[iColumn];
          for (CoinBigIndex j = start; j < end; j++) {
            double value = element[j];
            if (value == 1.0) {
            } else if (value == -1.0) {
              rhs[row[j]] = -0.5;
            } else {
              rhs[row[j]] = -COIN_DBL_MAX;
            }
          }
        }
      }
      for (iColumn = 0; iColumn < numberColumns; iColumn++) {
        if (upper[iColumn] > lower[iColumn] + 1.0e-8) {
          if (!isInteger(iColumn)) {
            CoinBigIndex start = columnStart[iColumn];
            CoinBigIndex end = start + columnLength[iColumn];
            double rhsValue = 0.0;
            // 1 all ones, -1 all -1s, 2 all +- 1, 3 no good
            int type = 0;
            for (CoinBigIndex j = start; j < end; j++) {
              double value = element[j];
              if (fabs(value) != 1.0) {
                type = 3;
                break;
              } else if (value == 1.0) {
                if (!type)
                  type = 1;
                else if (type != 1)
                  type = 2;
              } else {
                if (!type)
                  type = -1;
                else if (type != -1)
                  type = 2;
              }
              int iRow = row[j];
              if (rhs[iRow] == -COIN_DBL_MAX) {
                type = 3;
                break;
              } else if (rhs[iRow] == -0.5) {
                // different values
                unitRhs = false;
              } else if (rhsValue) {
                if (rhsValue != rhs[iRow])
                  unitRhs = false;
              } else {
                rhsValue = rhs[iRow];
              }
            }
            // if no elements OK
            if (type == 3) {
              // no good
              possible = 0;
              break;
            } else if (type == 2) {
              if (end - start > 2) {
                // no good
                possible = 0;
                break;
              } else {
                // only network
                possible &= 1;
                if (!possible)
                  break;
              }
            } else if (type == 1) {
              // only cover
              possible &= 2;
              if (!possible)
                break;
            } else if (type == -1) {
              // only negative cover
              possible &= 4;
              if (!possible)
                break;
            }
          }
        }
      }
      if ((possible == 2 || possible == 4) && !unitRhs) {
#if COIN_DEVELOP > 1
        printf("XXXXXX Continuous all +1 but different rhs\n");
#endif
        possible = 0;
      }
      // may be all integer
      if (possible != 7) {
        if (!possible)
          continuousMultiplier = 0.0;
        else if (possible == 1)
          continuousMultiplier = 1.0;
        else
          continuousMultiplier = 0.0; // 0.5 was incorrect;
#if COIN_DEVELOP > 1
        if (continuousMultiplier)
          printf("XXXXXX multiplier of %g\n", continuousMultiplier);
#endif
        if (continuousMultiplier == 0.5) {
          coeffMultiplier = new double[numberColumns];
          bool allOne = true;
          for (iColumn = 0; iColumn < numberColumns; iColumn++) {
            coeffMultiplier[iColumn] = 1.0;
            if (upper[iColumn] > lower[iColumn] + 1.0e-8) {
              if (!isInteger(iColumn)) {
                CoinBigIndex start = columnStart[iColumn];
                int iRow = row[start];
                double value = rhs[iRow];
                assert(value >= 0.0);
                if (value != 0.0 && value != 1.0)
                  allOne = false;
                coeffMultiplier[iColumn] = 0.5 * value;
              }
            }
          }
          if (allOne) {
            // back to old way
            delete[] coeffMultiplier;
            coeffMultiplier = NULL;
          }
        }
      } else {
        // all integer
        problemType_ = problemType;
#if COIN_DEVELOP > 1
        printf("Problem type is %d\n", problemType_);
#endif
      }
    }

    // But try again
    if (continuousMultiplier < 1.0) {
      memset(rhs, 0, numberRows * sizeof(double));
      int *count = new int[numberRows];
      memset(count, 0, numberRows * sizeof(int));
      for (iColumn = 0; iColumn < numberColumns; iColumn++) {
        CoinBigIndex start = columnStart[iColumn];
        CoinBigIndex end = start + columnLength[iColumn];
        if (upper[iColumn] == lower[iColumn]) {
          for (CoinBigIndex j = start; j < end; j++) {
            int iRow = row[j];
            rhs[iRow] += lower[iColumn] * element[j];
          }
        } else if (solver_->isInteger(iColumn)) {
          for (CoinBigIndex j = start; j < end; j++) {
            int iRow = row[j];
            if (fabs(element[j] - floor(element[j] + 0.5)) > 1.0e-10)
              rhs[iRow] = COIN_DBL_MAX;
          }
        } else {
          for (CoinBigIndex j = start; j < end; j++) {
            int iRow = row[j];
            count[iRow]++;
            if (fabs(element[j]) != 1.0)
              rhs[iRow] = COIN_DBL_MAX;
          }
        }
      }
      // now look at continuous
      bool allGood = true;
      double direction = solver_->getObjSense();
      int numberObj = 0;
      for (iColumn = 0; iColumn < numberColumns; iColumn++) {
        if (upper[iColumn] > lower[iColumn]) {
          double objValue = objective[iColumn] * direction;
          if (objValue && !solver_->isInteger(iColumn)) {
            numberObj++;
            CoinBigIndex start = columnStart[iColumn];
            CoinBigIndex end = start + columnLength[iColumn];
            if (objValue > 0.0) {
              // wants to be as low as possible
              if (lower[iColumn] < -1.0e10 || fabs(lower[iColumn] - floor(lower[iColumn] + 0.5)) > 1.0e-10) {
                allGood = false;
                break;
              } else if (upper[iColumn] < 1.0e10 && fabs(upper[iColumn] - floor(upper[iColumn] + 0.5)) > 1.0e-10) {
                allGood = false;
                break;
              }
              bool singletonRow = true;
              bool equality = false;
              for (CoinBigIndex j = start; j < end; j++) {
                int iRow = row[j];
                if (count[iRow] > 1)
                  singletonRow = false;
                else if (rowLower[iRow] == rowUpper[iRow])
                  equality = true;
                double rhsValue = rhs[iRow];
                double lowerValue = rowLower[iRow];
                double upperValue = rowUpper[iRow];
                if (rhsValue < 1.0e20) {
                  if (lowerValue > -1.0e20)
                    lowerValue -= rhsValue;
                  if (upperValue < 1.0e20)
                    upperValue -= rhsValue;
                }
                if (fabs(rhsValue) > 1.0e20 || fabs(rhsValue - floor(rhsValue + 0.5)) > 1.0e-10
                  || fabs(element[j]) != 1.0) {
                  // no good
                  allGood = false;
                  break;
                }
                if (element[j] > 0.0) {
                  if (lowerValue > -1.0e20 && fabs(lowerValue - floor(lowerValue + 0.5)) > 1.0e-10) {
                    // no good
                    allGood = false;
                    break;
                  }
                } else {
                  if (upperValue < 1.0e20 && fabs(upperValue - floor(upperValue + 0.5)) > 1.0e-10) {
                    // no good
                    allGood = false;
                    break;
                  }
                }
              }
              if (!singletonRow && end > start + 1 && !equality)
                allGood = false;
              if (!allGood)
                break;
            } else {
              // wants to be as high as possible
              if (upper[iColumn] > 1.0e10 || fabs(upper[iColumn] - floor(upper[iColumn] + 0.5)) > 1.0e-10) {
                allGood = false;
                break;
              } else if (lower[iColumn] > -1.0e10 && fabs(lower[iColumn] - floor(lower[iColumn] + 0.5)) > 1.0e-10) {
                allGood = false;
                break;
              }
              bool singletonRow = true;
              bool equality = false;
              for (CoinBigIndex j = start; j < end; j++) {
                int iRow = row[j];
                if (count[iRow] > 1)
                  singletonRow = false;
                else if (rowLower[iRow] == rowUpper[iRow])
                  equality = true;
                double rhsValue = rhs[iRow];
                double lowerValue = rowLower[iRow];
                double upperValue = rowUpper[iRow];
                if (rhsValue < 1.0e20) {
                  if (lowerValue > -1.0e20)
                    lowerValue -= rhsValue;
                  if (upperValue < 1.0e20)
                    upperValue -= rhsValue;
                }
                if (fabs(rhsValue) > 1.0e20 || fabs(rhsValue - floor(rhsValue + 0.5)) > 1.0e-10
                  || fabs(element[j]) != 1.0) {
                  // no good
                  allGood = false;
                  break;
                }
                if (element[j] < 0.0) {
                  if (lowerValue > -1.0e20 && fabs(lowerValue - floor(lowerValue + 0.5)) > 1.0e-10) {
                    // no good
                    allGood = false;
                    break;
                  }
                } else {
                  if (upperValue < 1.0e20 && fabs(upperValue - floor(upperValue + 0.5)) > 1.0e-10) {
                    // no good
                    allGood = false;
                    break;
                  }
                }
              }
              if (!singletonRow && end > start + 1 && !equality)
                allGood = false;
              if (!allGood)
                break;
            }
          }
        }
      }
      delete[] count;
      if (allGood) {
#if COIN_DEVELOP > 1
        if (numberObj)
          printf("YYYY analysis says all continuous with costs will be integer\n");
#endif
        continuousMultiplier = 1.0;
      }
    }
    delete[] rhs;
  }
  /*
      Take a first scan to see if there are unfixed continuous variables in the
      objective.  If so, the minimum objective change could be arbitrarily small.
      Also pick off the maximum coefficient of an unfixed integer variable.

      If the objective is found to contain only integer variables, set the
      fathoming discipline to strict.
    */
  double maximumCost = 0.0;
  //double trueIncrement=0.0;
  int iColumn;
  int numberColumns = getNumCols();
  double scaleFactor = 1.0; // due to rhs etc
  /*
      Original model did not have integer bounds.
    */
  if ((specialOptions_ & 65536) == 0) {
    /* be on safe side (later look carefully as may be able to
           to get 0.5 say if bounds are multiples of 0.5 */
    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
      if (upper[iColumn] > lower[iColumn] + 1.0e-8) {
        double value;
        value = fabs(lower[iColumn]);
        if (floor(value + 0.5) != value) {
          scaleFactor = CoinMin(scaleFactor, 0.5);
          if (floor(2.0 * value + 0.5) != 2.0 * value) {
            scaleFactor = CoinMin(scaleFactor, 0.25);
            if (floor(4.0 * value + 0.5) != 4.0 * value) {
              scaleFactor = 0.0;
            }
          }
        }
        value = fabs(upper[iColumn]);
        if (floor(value + 0.5) != value) {
          scaleFactor = CoinMin(scaleFactor, 0.5);
          if (floor(2.0 * value + 0.5) != 2.0 * value) {
            scaleFactor = CoinMin(scaleFactor, 0.25);
            if (floor(4.0 * value + 0.5) != 4.0 * value) {
              scaleFactor = 0.0;
            }
          }
        }
      }
    }
  }
  bool possibleMultiple = continuousMultiplier != 0.0 && scaleFactor != 0.0;
  if (possibleMultiple) {
    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
      if (upper[iColumn] > lower[iColumn] + 1.0e-8) {
        maximumCost = CoinMax(maximumCost, fabs(objective[iColumn]));
      }
    }
  }
  setIntParam(CbcModel::CbcFathomDiscipline, possibleMultiple);
  /*
      If a nontrivial increment is possible, try and figure it out. We're looking
      for gcd(c<j>) for all c<j> that are coefficients of unfixed integer
      variables. Since the c<j> might not be integers, try and inflate them
      sufficiently that they look like integers (and we'll deflate the gcd
      later).

      2520.0 is used as it is a nice multiple of 2,3,5,7
    */
  if (possibleMultiple && maximumCost) {
    int increment = 0;
    double multiplier = 2520.0;
    while (10.0 * multiplier * maximumCost < 1.0e8)
      multiplier *= 10.0;
    int bigIntegers = 0; // Count of large costs which are integer
    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
      if (upper[iColumn] > lower[iColumn] + 1.0e-8) {
        double objValue = fabs(objective[iColumn]);
        if (!isInteger(iColumn)) {
          if (!coeffMultiplier)
            objValue *= continuousMultiplier;
          else
            objValue *= coeffMultiplier[iColumn];
        }
        if (objValue) {
          double value = objValue * multiplier;
          if (value < 2.1e9) {
            int nearest = static_cast<int>(floor(value + 0.5));
            if (fabs(value - floor(value + 0.5)) > 1.0e-8) {
              increment = 0;
              break;
            } else if (!increment) {
              increment = nearest;
            } else {
              increment = gcd(increment, nearest);
            }
          } else {
            // large value - may still be multiple of 1.0
            if (fabs(objValue - floor(objValue + 0.5)) > 1.0e-8) {
              increment = 0;
              break;
            } else {
              bigIntegers++;
            }
          }
        }
      }
    }
    delete[] coeffMultiplier;
    /*
          If the increment beats the current value for objective change, install it.
        */
    if (increment) {
      double value = increment;
      double cutoff = getDblParam(CbcModel::CbcCutoffIncrement);
      if (bigIntegers) {
        // allow for 1.0
        increment = gcd(increment, static_cast<int>(multiplier));
        value = increment;
      }
      value /= multiplier;
      value *= scaleFactor;
      //trueIncrement=CoinMax(cutoff,value);;
      if (value * 0.999 > cutoff) {
        messageHandler()->message(CBC_INTEGERINCREMENT,
          messages())
          << value << CoinMessageEol;
        setDblParam(CbcModel::CbcCutoffIncrement, CoinMax(value * 0.999, value - 1.0e-4));
      }
    }
  }

  return;
}

/*
saveModel called (carved out of) BranchandBound
*/
void CbcModel::saveModel(OsiSolverInterface *saveSolver, double *checkCutoffForRestart, bool *feasible)
{
  if (saveSolver && (specialOptions_ & 32768) != 0) {
    // See if worth trying reduction
    *checkCutoffForRestart = getCutoff();
    bool tryNewSearch = solverCharacteristics_->reducedCostsAccurate() && (*checkCutoffForRestart < 1.0e20);
    int numberColumns = getNumCols();
    if (tryNewSearch) {
#if CBC_USEFUL_PRINTING > 1
      printf("after %d nodes, cutoff %g - looking\n",
        numberNodes_, getCutoff());
#endif
      saveSolver->resolve();
      double direction = saveSolver->getObjSense();
      double gap = *checkCutoffForRestart - saveSolver->getObjValue() * direction;
      double tolerance;
      saveSolver->getDblParam(OsiDualTolerance, tolerance);
      if (gap <= 0.0)
        gap = tolerance;
      gap += 100.0 * tolerance;
      double integerTolerance = getDblParam(CbcIntegerTolerance);

      const double *lower = saveSolver->getColLower();
      const double *upper = saveSolver->getColUpper();
      const double *solution = saveSolver->getColSolution();
      const double *reducedCost = saveSolver->getReducedCost();

      int numberFixed = 0;
      int numberFixed2 = 0;
      for (int i = 0; i < numberIntegers_; i++) {
        int iColumn = integerVariable_[i];
        double djValue = direction * reducedCost[iColumn];
        if (upper[iColumn] - lower[iColumn] > integerTolerance) {
          if (solution[iColumn] < lower[iColumn] + integerTolerance && djValue > gap) {
            saveSolver->setColUpper(iColumn, lower[iColumn]);
            numberFixed++;
          } else if (solution[iColumn] > upper[iColumn] - integerTolerance && -djValue > gap) {
            saveSolver->setColLower(iColumn, upper[iColumn]);
            numberFixed++;
          }
        } else {
          numberFixed2++;
        }
      }
#ifdef COIN_DEVELOP
      /*
              We're debugging. (specialOptions 1)
            */
      if ((specialOptions_ & 1) != 0) {
        const OsiRowCutDebugger *debugger = saveSolver->getRowCutDebugger();
        if (debugger) {
          printf("Contains optimal\n");
          OsiSolverInterface *temp = saveSolver->clone();
          const double *solution = debugger->optimalSolution();
          const double *lower = temp->getColLower();
          const double *upper = temp->getColUpper();
          int n = temp->getNumCols();
          for (int i = 0; i < n; i++) {
            if (temp->isInteger(i)) {
              double value = floor(solution[i] + 0.5);
              assert(value >= lower[i] && value <= upper[i]);
              temp->setColLower(i, value);
              temp->setColUpper(i, value);
            }
          }
          temp->writeMps("reduced_fix");
          delete temp;
          saveSolver->writeMps("reduced");
        } else {
          abort();
        }
      }
      printf("Restart could fix %d integers (%d already fixed)\n",
        numberFixed + numberFixed2, numberFixed2);
#endif
      numberFixed += numberFixed2;
      if (numberFixed * 20 < numberColumns)
        tryNewSearch = false;
    }
    if (tryNewSearch) {
      // back to solver without cuts?
      OsiSolverInterface *solver2 = continuousSolver_->clone();
      const double *lower = saveSolver->getColLower();
      const double *upper = saveSolver->getColUpper();
      for (int i = 0; i < numberIntegers_; i++) {
        int iColumn = integerVariable_[i];
        solver2->setColLower(iColumn, lower[iColumn]);
        solver2->setColUpper(iColumn, upper[iColumn]);
      }
      // swap
      delete saveSolver;
      saveSolver = solver2;
      double *newSolution = new double[numberColumns];
      double objectiveValue = *checkCutoffForRestart;
      CbcSerendipity heuristic(*this);
      if (bestSolution_)
        heuristic.setInputSolution(bestSolution_, bestObjective_);
      heuristic.setFractionSmall(0.9);
      heuristic.setFeasibilityPumpOptions(1008013);
      // Use numberNodes to say how many are original rows
      heuristic.setNumberNodes(continuousSolver_->getNumRows());
#ifdef COIN_DEVELOP
      if (continuousSolver_->getNumRows() < saveSolver->getNumRows())
        printf("%d rows added ZZZZZ\n",
          solver_->getNumRows() - continuousSolver_->getNumRows());
#endif
      int returnCode = heuristic.smallBranchAndBound(saveSolver,
        -1, newSolution,
        objectiveValue,
        *checkCutoffForRestart, "Reduce");
      if (returnCode < 0) {
#ifdef COIN_DEVELOP
        printf("Restart - not small enough to do search after fixing\n");
#endif
        delete[] newSolution;
      } else {
        if ((returnCode & 1) != 0) {
          // increment number of solutions so other heuristics can test
          numberSolutions_++;
          numberHeuristicSolutions_++;
          lastHeuristic_ = NULL;
          setBestSolution(CBC_ROUNDING, objectiveValue, newSolution);
        }
        delete[] newSolution;
        *feasible = false; // stop search
      }
#if 0 // probably not needed def CBC_THREAD
            if (master_) {
                lockThread();
                if (parallelMode() > 0) {
                    while (master_->waitForThreadsInTree(0)) {
                        lockThread();
                        double dummyBest;
                        tree_->cleanTree(this, -COIN_DBL_MAX, dummyBest) ;
                        //unlockThread();
                    }
                }
                master_->waitForThreadsInTree(2);
                delete master_;
                master_ = NULL;
                masterThread_ = NULL;
            }
#endif
    }
  }
}
/*
Adds integers, called from BranchandBound()
*/
void CbcModel::AddIntegers()
{
  int numberColumns = continuousSolver_->getNumCols();
  int numberRows = continuousSolver_->getNumRows();
  int numberOriginalIntegers = numberIntegers_;
  int *del = new int[CoinMax(numberColumns, numberRows)];
  int *original = new int[numberColumns];
  char *possibleRow = new char[numberRows];
  {
    const CoinPackedMatrix *rowCopy = continuousSolver_->getMatrixByRow();
    const int *column = rowCopy->getIndices();
    const int *rowLength = rowCopy->getVectorLengths();
    const CoinBigIndex *rowStart = rowCopy->getVectorStarts();
    const double *rowLower = continuousSolver_->getRowLower();
    const double *rowUpper = continuousSolver_->getRowUpper();
    const double *element = rowCopy->getElements();
    for (int i = 0; i < numberRows; i++) {
      int nLeft = 0;
      bool possible = false;
      if (rowLower[i] < -1.0e20) {
        double value = rowUpper[i];
        if (fabs(value - floor(value + 0.5)) < 1.0e-8)
          possible = true;
      } else if (rowUpper[i] > 1.0e20) {
        double value = rowLower[i];
        if (fabs(value - floor(value + 0.5)) < 1.0e-8)
          possible = true;
      } else {
        double value = rowUpper[i];
        if (rowLower[i] == rowUpper[i] && fabs(value - floor(value + 0.5)) < 1.0e-8)
          possible = true;
      }
      double allSame = (possible) ? 0.0 : -1.0;
      for (CoinBigIndex j = rowStart[i];
           j < rowStart[i] + rowLength[i]; j++) {
        int iColumn = column[j];
        if (continuousSolver_->isInteger(iColumn)) {
          if (fabs(element[j]) != 1.0)
            possible = false;
        } else {
          nLeft++;
          if (!allSame) {
            allSame = fabs(element[j]);
          } else if (allSame > 0.0) {
            if (allSame != fabs(element[j]))
              allSame = -1.0;
          }
        }
      }
      if (nLeft == rowLength[i] && allSame > 0.0)
        possibleRow[i] = 2;
      else if (possible || !nLeft)
        possibleRow[i] = 1;
      else
        possibleRow[i] = 0;
    }
  }
  int nDel = 0;
  for (int i = 0; i < numberColumns; i++) {
    original[i] = i;
    if (continuousSolver_->isInteger(i))
      del[nDel++] = i;
  }
  {
    // we must not exclude current best solution (rounding errors)
    // also not if large values
    const int *row = continuousSolver_->getMatrixByCol()->getIndices();
    const CoinBigIndex *columnStart = continuousSolver_->getMatrixByCol()->getVectorStarts();
    const int *columnLength = continuousSolver_->getMatrixByCol()->getVectorLengths();
    const double *solution = continuousSolver_->getColSolution();
    for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
      if (!continuousSolver_->isInteger(iColumn)) {
        double value = bestSolution_ ? bestSolution_[iColumn] : 0.0;
        double value2 = solution[iColumn];
        if (fabs(value - floor(value + 0.5)) > 1.0e-8 || fabs(value2) > 1.0e3) {
          CoinBigIndex start = columnStart[iColumn];
          CoinBigIndex end = start + columnLength[iColumn];
          for (CoinBigIndex j = start; j < end; j++) {
            int iRow = row[j];
            possibleRow[iRow] = 0;
          }
        }
      }
    }
  }
  int nExtra = 0;
  OsiSolverInterface *copy1 = continuousSolver_->clone();
  int nPass = 0;
  while (nDel && nPass < 10) {
    nPass++;
    OsiSolverInterface *copy2 = copy1->clone();
    int nLeft = 0;
    for (int i = 0; i < nDel; i++)
      original[del[i]] = -1;
    for (int i = 0; i < numberColumns; i++) {
      int kOrig = original[i];
      if (kOrig >= 0)
        original[nLeft++] = kOrig;
    }
    assert(nLeft == numberColumns - nDel);
    copy2->deleteCols(nDel, del);
    numberColumns = copy2->getNumCols();
    const CoinPackedMatrix *rowCopy = copy2->getMatrixByRow();
    numberRows = rowCopy->getNumRows();
    const int *column = rowCopy->getIndices();
    const int *rowLength = rowCopy->getVectorLengths();
    const CoinBigIndex *rowStart = rowCopy->getVectorStarts();
    const double *rowLower = copy2->getRowLower();
    const double *rowUpper = copy2->getRowUpper();
    const double *element = rowCopy->getElements();
    const CoinPackedMatrix *columnCopy = copy2->getMatrixByCol();
    const int *columnLength = columnCopy->getVectorLengths();
    nDel = 0;
    // Could do gcd stuff on ones with costs
    for (int i = 0; i < numberRows; i++) {
      if (!rowLength[i]) {
        del[nDel++] = i;
      } else if (possibleRow[i]) {
        if (rowLength[i] == 1) {
          CoinBigIndex k = rowStart[i];
          int iColumn = column[k];
          if (!copy2->isInteger(iColumn)) {
            double mult = 1.0 / fabs(element[k]);
            if (rowLower[i] < -1.0e20) {
              // treat rhs as multiple of 1 unless elements all same
              double value = ((possibleRow[i] == 2) ? rowUpper[i] : 1.0) * mult;
              if (fabs(value - floor(value + 0.5)) < 1.0e-8) {
                del[nDel++] = i;
                if (columnLength[iColumn] == 1) {
                  copy2->setInteger(iColumn);
                  int kOrig = original[iColumn];
                  setOptionalInteger(kOrig);
                }
              }
            } else if (rowUpper[i] > 1.0e20) {
              // treat rhs as multiple of 1 unless elements all same
              double value = ((possibleRow[i] == 2) ? rowLower[i] : 1.0) * mult;
              if (fabs(value - floor(value + 0.5)) < 1.0e-8) {
                del[nDel++] = i;
                if (columnLength[iColumn] == 1) {
                  copy2->setInteger(iColumn);
                  int kOrig = original[iColumn];
                  setOptionalInteger(kOrig);
                }
              }
            } else {
              // treat rhs as multiple of 1 unless elements all same
              double value = ((possibleRow[i] == 2) ? rowUpper[i] : 1.0) * mult;
              if (rowLower[i] == rowUpper[i] && fabs(value - floor(value + 0.5)) < 1.0e-8) {
                del[nDel++] = i;
                copy2->setInteger(iColumn);
                int kOrig = original[iColumn];
                setOptionalInteger(kOrig);
              }
            }
          }
        } else {
          // only if all singletons
          bool possible = false;
          if (rowLower[i] < -1.0e20) {
            double value = rowUpper[i];
            if (fabs(value - floor(value + 0.5)) < 1.0e-8)
              possible = true;
          } else if (rowUpper[i] > 1.0e20) {
            double value = rowLower[i];
            if (fabs(value - floor(value + 0.5)) < 1.0e-8)
              possible = true;
          } else {
            double value = rowUpper[i];
            if (rowLower[i] == rowUpper[i] && fabs(value - floor(value + 0.5)) < 1.0e-8)
              possible = true;
          }
          if (possible) {
            for (CoinBigIndex j = rowStart[i];
                 j < rowStart[i] + rowLength[i]; j++) {
              int iColumn = column[j];
              if (columnLength[iColumn] != 1 || fabs(element[j]) != 1.0) {
                possible = false;
                break;
              }
            }
            if (possible) {
              for (CoinBigIndex j = rowStart[i];
                   j < rowStart[i] + rowLength[i]; j++) {
                int iColumn = column[j];
                if (!copy2->isInteger(iColumn)) {
                  copy2->setInteger(iColumn);
                  int kOrig = original[iColumn];
                  setOptionalInteger(kOrig);
                }
              }
              del[nDel++] = i;
            }
          }
        }
      }
    }
    if (nDel) {
      copy2->deleteRows(nDel, del);
      // pack down possible
      int n = 0;
      for (int i = 0; i < nDel; i++)
        possibleRow[del[i]] = -1;
      for (int i = 0; i < numberRows; i++) {
        if (possibleRow[i] >= 0)
          possibleRow[n++] = possibleRow[i];
      }
    }
    if (nDel != numberRows) {
      nDel = 0;
      for (int i = 0; i < numberColumns; i++) {
        if (copy2->isInteger(i)) {
          del[nDel++] = i;
          nExtra++;
        }
      }
    } else {
      nDel = 0;
    }
    delete copy1;
    copy1 = copy2->clone();
    delete copy2;
  }
  // See if what's left is a network
  bool couldBeNetwork = false;
  if (copy1->getNumRows() && copy1->getNumCols()) {
#ifdef COIN_HAS_CLP
    OsiClpSolverInterface *clpSolver
      = dynamic_cast<OsiClpSolverInterface *>(copy1);
    if (false && clpSolver) {
      numberRows = clpSolver->getNumRows();
      char *rotate = new char[numberRows];
      int n = clpSolver->getModelPtr()->findNetwork(rotate, 1.0);
      delete[] rotate;
#if CBC_USEFUL_PRINTING > 1
      printf("INTA network %d rows out of %d\n", n, numberRows);
#endif
      if (CoinAbs(n) == numberRows) {
        couldBeNetwork = true;
        for (int i = 0; i < numberRows; i++) {
          if (!possibleRow[i]) {
            couldBeNetwork = false;
#if CBC_USEFUL_PRINTING > 1
            printf("but row %d is bad\n", i);
#endif
            break;
          }
        }
      }
    } else
#endif
    {
      numberColumns = copy1->getNumCols();
      numberRows = copy1->getNumRows();
      const double *rowLower = copy1->getRowLower();
      const double *rowUpper = copy1->getRowUpper();
      couldBeNetwork = true;
      for (int i = 0; i < numberRows; i++) {
        if (rowLower[i] > -1.0e20 && fabs(rowLower[i] - floor(rowLower[i] + 0.5)) > 1.0e-12) {
          couldBeNetwork = false;
          break;
        }
        if (rowUpper[i] < 1.0e20 && fabs(rowUpper[i] - floor(rowUpper[i] + 0.5)) > 1.0e-12) {
          couldBeNetwork = false;
          break;
        }
        if (possibleRow[i] == 0) {
          couldBeNetwork = false;
          break;
        }
      }
      if (couldBeNetwork) {
        const CoinPackedMatrix *matrixByCol = copy1->getMatrixByCol();
        const double *element = matrixByCol->getElements();
        //const int * row = matrixByCol->getIndices();
        const CoinBigIndex *columnStart = matrixByCol->getVectorStarts();
        const int *columnLength = matrixByCol->getVectorLengths();
        for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
          CoinBigIndex start = columnStart[iColumn];
          CoinBigIndex end = start + columnLength[iColumn];
          if (end > start + 2) {
            couldBeNetwork = false;
            break;
          }
          int type = 0;
          for (CoinBigIndex j = start; j < end; j++) {
            double value = element[j];
            if (fabs(value) != 1.0) {
              couldBeNetwork = false;
              break;
            } else if (value == 1.0) {
              if ((type & 1) == 0)
                type |= 1;
              else
                type = 7;
            } else if (value == -1.0) {
              if ((type & 2) == 0)
                type |= 2;
              else
                type = 7;
            }
          }
          if (type > 3) {
            couldBeNetwork = false;
            break;
          }
        }
      }
    }
  }
  if (couldBeNetwork) {
    for (int i = 0; i < numberColumns; i++)
      setOptionalInteger(original[i]);
  }
  if (nExtra || couldBeNetwork) {
    numberColumns = copy1->getNumCols();
    numberRows = copy1->getNumRows();
    if (!numberColumns || !numberRows) {
      int numberColumns = solver_->getNumCols();
      for (int i = 0; i < numberColumns; i++)
        assert(solver_->isInteger(i));
    }
#if CBC_USEFUL_PRINTING > 1
    if (couldBeNetwork || nExtra)
      printf("INTA %d extra integers, %d left%s\n", nExtra,
        numberColumns,
        couldBeNetwork ? ", all network" : "");
#endif
    findIntegers(true, 2);
    convertToDynamic();
  }
#if CBC_USEFUL_PRINTING > 1
  if (!couldBeNetwork && copy1->getNumCols() && copy1->getNumRows()) {
    printf("INTA %d rows and %d columns remain\n",
      copy1->getNumRows(), copy1->getNumCols());
    if (copy1->getNumCols() < 200) {
      copy1->writeMps("moreint");
      printf("INTA Written remainder to moreint.mps.gz %d rows %d cols\n",
        copy1->getNumRows(), copy1->getNumCols());
    }
  }
#endif
  delete copy1;
  delete[] del;
  delete[] original;
  delete[] possibleRow;
  // double check increment
  analyzeObjective();
  // If any changes - tell code
  if (numberOriginalIntegers < numberIntegers_)
    synchronizeModel();
}
/**
  \todo
  Normally, it looks like we enter here from command dispatch in the main
  routine, after calling the solver for an initial solution
  (CbcModel::initialSolve, which simply calls the solver's initialSolve
  routine.) The first thing we do is call resolve. Presumably there are
  circumstances where this is nontrivial? There's also a call from
  CbcModel::originalModel (tied up with integer presolve), which should be
  checked.

*/

#ifdef CONFLICT_CUTS
#if PRINT_CONFLICT == 1
static int numberConflictCuts = 0;
static int lastNumberConflictCuts = 0;
static double lengthConflictCuts = 0.0;
#endif
#endif
/*
  The overall flow can be divided into three stages:
    * Prep: Check that the lp relaxation remains feasible at the root. If so,
      do all the setup for B&C.
    * Process the root node: Generate cuts, apply heuristics, and in general do
      the best we can to resolve the problem without B&C.
    * Do B&C search until we hit a limit or exhaust the search tree.

  Keep in mind that in general there is no node in the search tree that
  corresponds to the active subproblem. The active subproblem is represented
  by the current state of the model,  of the solver, and of the constraint
  system held by the solver.
*/
#ifdef COIN_HAS_CPX
#include "OsiCpxSolverInterface.hpp"
#include "cplex.h"
#endif
void CbcModel::branchAndBound(int doStatistics)

{
  if (!parentModel_) {
    /*
        Capture a time stamp before we start (unless set).
      */
    if (!dblParam_[CbcStartSeconds]) {
      if (!useElapsedTime())
        dblParam_[CbcStartSeconds] = CoinCpuTime();
      else
        dblParam_[CbcStartSeconds] = CoinGetTimeOfDay();
    }
  }
  dblParam_[CbcSmallestChange] = COIN_DBL_MAX;
  dblParam_[CbcSumChange] = 0.0;
  dblParam_[CbcLargestChange] = 0.0;
  intParam_[CbcNumberBranches] = 0;
  double lastBestPossibleObjective = -COIN_DBL_MAX;
  // when to check for restart
  int nextCheckRestart = 50;
  // Force minimization !!!!
  bool flipObjective = (solver_->getObjSense() < 0.0);
  if (flipObjective)
    flipModel();
  dblParam_[CbcOptimizationDirection] = 1.0; // was solver_->getObjSense();
  strongInfo_[0] = 0;
  strongInfo_[1] = 0;
  strongInfo_[2] = 0;
  strongInfo_[3] = 0;
  strongInfo_[4] = 0;
  strongInfo_[5] = 0;
  strongInfo_[6] = 0;
  numberStrongIterations_ = 0;
  currentNode_ = NULL;
  // See if should do cuts old way
  if (parallelMode() < 0) {
    specialOptions_ |= 4096 + 8192;
  } else if (parallelMode() > 0) {
    specialOptions_ |= 4096;
  }
  int saveMoreSpecialOptions = moreSpecialOptions_;
  if (dynamic_cast<CbcTreeLocal *>(tree_))
    specialOptions_ |= 4096 + 8192;
#ifdef COIN_HAS_CLP
  {
    OsiClpSolverInterface *clpSolver
      = dynamic_cast<OsiClpSolverInterface *>(solver_);
    if (clpSolver) {
      // pass in disaster handler
      CbcDisasterHandler handler(this);
      clpSolver->passInDisasterHandler(&handler);
      // Initialise solvers seed (unless users says not)
      if ((specialOptions_ & 4194304) == 0)
        clpSolver->getModelPtr()->setRandomSeed(1234567);
#ifdef JJF_ZERO
      // reduce factorization frequency
      int frequency = clpSolver->getModelPtr()->factorizationFrequency();
      clpSolver->getModelPtr()->setFactorizationFrequency(CoinMin(frequency, 120));
#endif
    }
  }
#endif
  // original solver (only set if pre-processing)
  OsiSolverInterface *originalSolver = NULL;
  int numberOriginalObjects = numberObjects_;
  OsiObject **originalObject = NULL;
  // Save whether there were any objects
  bool noObjects = (numberObjects_ == 0);
  // Set up strategies
  /*
      See if the user has supplied a strategy object and deal with it if present.
      The call to setupOther will set numberStrong_ and numberBeforeTrust_, and
      perform integer preprocessing, if requested.

      We need to hang on to a pointer to solver_. setupOther will assign a
      preprocessed solver to model, but will instruct assignSolver not to trash the
      existing one.
    */
  if (strategy_) {
    // May do preprocessing
    originalSolver = solver_;
    strategy_->setupOther(*this);
    if (strategy_->preProcessState()) {
      // pre-processing done
      if (strategy_->preProcessState() < 0) {
        // infeasible (or unbounded)
        status_ = 0;
        if (!solver_->isProvenDualInfeasible()) {
          handler_->message(CBC_INFEAS, messages_) << CoinMessageEol;
          secondaryStatus_ = 1;
        } else {
          handler_->message(CBC_UNBOUNDED,
            messages_)
            << CoinMessageEol;
          secondaryStatus_ = 7;
        }
        originalContinuousObjective_ = COIN_DBL_MAX;
        if (flipObjective)
          flipModel();
        return;
      } else if (numberObjects_ && object_) {
        numberOriginalObjects = numberObjects_;
        // redo sequence
        numberIntegers_ = 0;
        int numberColumns = getNumCols();
        int nOrig = originalSolver->getNumCols();
        CglPreProcess *process = strategy_->process();
        assert(process);
        const int *originalColumns = process->originalColumns();
        // allow for cliques etc
        nOrig = CoinMax(nOrig, originalColumns[numberColumns - 1] + 1);
        // try and redo debugger
        OsiRowCutDebugger *debugger = const_cast<OsiRowCutDebugger *>(solver_->getRowCutDebuggerAlways());
        if (debugger) {
          if (numberColumns <= debugger->numberColumns())
            debugger->redoSolution(numberColumns, originalColumns);
          else
            debugger = NULL; // no idea how to handle (SOS?)
        }
        // User-provided solution might have been best. Synchronise.
        if (bestSolution_) {
          // need to redo - in case no better found in BAB
          // just get integer part right
          for (int i = 0; i < numberColumns; i++) {
            int jColumn = originalColumns[i];
            bestSolution_[i] = bestSolution_[jColumn];
          }
        }
        originalObject = object_;
        // object number or -1
        int *temp = new int[nOrig];
        int iColumn;
        for (iColumn = 0; iColumn < nOrig; iColumn++)
          temp[iColumn] = -1;
        int iObject;
        int nNonInt = 0;
        for (iObject = 0; iObject < numberOriginalObjects; iObject++) {
          iColumn = originalObject[iObject]->columnNumber();
          if (iColumn < 0) {
            nNonInt++;
          } else {
            temp[iColumn] = iObject;
          }
        }
        int numberNewIntegers = 0;
        int numberOldIntegers = 0;
        int numberOldOther = 0;
        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
          int jColumn = originalColumns[iColumn];
          if (temp[jColumn] >= 0) {
            int iObject = temp[jColumn];
            CbcSimpleInteger *obj = dynamic_cast<CbcSimpleInteger *>(originalObject[iObject]);
            if (obj)
              numberOldIntegers++;
            else
              numberOldOther++;
          } else if (isInteger(iColumn)) {
            numberNewIntegers++;
          }
        }
        /*
                  Allocate an array to hold the indices of the integer variables.
                  Make a large enough array for all objects
                */
        numberObjects_ = numberNewIntegers + numberOldIntegers + numberOldOther + nNonInt;
        object_ = new OsiObject *[numberObjects_];
        delete[] integerVariable_;
        integerVariable_ = new int[numberNewIntegers + numberOldIntegers];
        /*
                  Walk the variables again, filling in the indices and creating objects for
                  the integer variables. Initially, the objects hold the index and upper &
                  lower bounds.
                */
        numberIntegers_ = 0;
        int n = originalColumns[numberColumns - 1] + 1;
        int *backward = new int[n];
        int i;
        for (i = 0; i < n; i++)
          backward[i] = -1;
        for (i = 0; i < numberColumns; i++)
          backward[originalColumns[i]] = i;
        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
          int jColumn = originalColumns[iColumn];
          if (temp[jColumn] >= 0) {
            int iObject = temp[jColumn];
            CbcSimpleInteger *obj = dynamic_cast<CbcSimpleInteger *>(originalObject[iObject]);
            if (obj) {
              object_[numberIntegers_] = originalObject[iObject]->clone();
              // redo ids etc
              //object_[numberIntegers_]->resetSequenceEtc(numberColumns,originalColumns);
              object_[numberIntegers_]->resetSequenceEtc(numberColumns, backward);
              integerVariable_[numberIntegers_++] = iColumn;
            }
          } else if (isInteger(iColumn)) {
            object_[numberIntegers_] = new CbcSimpleInteger(this, iColumn);
            integerVariable_[numberIntegers_++] = iColumn;
          }
        }
        delete[] backward;
        numberObjects_ = numberIntegers_;
        // Now append other column stuff
        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
          int jColumn = originalColumns[iColumn];
          if (temp[jColumn] >= 0) {
            int iObject = temp[jColumn];
            CbcSimpleInteger *obj = dynamic_cast<CbcSimpleInteger *>(originalObject[iObject]);
            if (!obj) {
              object_[numberObjects_] = originalObject[iObject]->clone();
              // redo ids etc
              CbcObject *obj = dynamic_cast<CbcObject *>(object_[numberObjects_]);
              assert(obj);
              obj->redoSequenceEtc(this, numberColumns, originalColumns);
              numberObjects_++;
            }
          }
        }
        // now append non column stuff
        for (iObject = 0; iObject < numberOriginalObjects; iObject++) {
          iColumn = originalObject[iObject]->columnNumber();
          if (iColumn < 0) {
            // already has column numbers changed
            object_[numberObjects_] = originalObject[iObject]->clone();
#ifdef JJF_ZERO
            // redo ids etc
            CbcObject *obj = dynamic_cast<CbcObject *>(object_[numberObjects_]);
            assert(obj);
            obj->redoSequenceEtc(this, numberColumns, originalColumns);
#endif
            numberObjects_++;
          }
        }
        delete[] temp;
        if (!numberObjects_)
          handler_->message(CBC_NOINT, messages_) << CoinMessageEol;
      } else {
        int numberColumns = getNumCols();
        CglPreProcess *process = strategy_->process();
        assert(process);
        const int *originalColumns = process->originalColumns();
        // try and redo debugger
        OsiRowCutDebugger *debugger = const_cast<OsiRowCutDebugger *>(solver_->getRowCutDebuggerAlways());
        if (debugger)
          debugger->redoSolution(numberColumns, originalColumns);
      }
    } else {
      //no preprocessing
      originalSolver = NULL;
    }
    strategy_->setupCutGenerators(*this);
    strategy_->setupHeuristics(*this);
    // Set strategy print level to models
    strategy_->setupPrinting(*this, handler_->logLevel());
  }
  eventHappened_ = false;
  CbcEventHandler *eventHandler = getEventHandler();
  if (eventHandler)
    eventHandler->setModel(this);
#define CLIQUE_ANALYSIS
#ifdef CLIQUE_ANALYSIS
  // set up for probing
  // If we're doing clever stuff with cliques, additional info here.
  if (!parentModel_)
    probingInfo_ = new CglTreeProbingInfo(solver_);
  else
    probingInfo_ = NULL;
#else
  probingInfo_ = NULL;
#endif

  // Try for dominated columns
  if ((specialOptions_ & 64) != 0) {
    CglDuplicateRow dupcuts(solver_);
    dupcuts.setMode(2);
    CglStored *storedCuts = dupcuts.outDuplicates(solver_);
    if (storedCuts) {
      COIN_DETAIL_PRINT(printf("adding dup cuts\n"));
      addCutGenerator(storedCuts, 1, "StoredCuts from dominated",
        true, false, false, -200);
    }
  }
  if (!nodeCompare_)
    nodeCompare_ = new CbcCompareDefault();
  ;
  // See if hot start wanted
  CbcCompareBase *saveCompare = NULL;
  // User supplied hotstart. Adapt for preprocessing.
  if (hotstartSolution_) {
    if (strategy_ && strategy_->preProcessState() > 0) {
      CglPreProcess *process = strategy_->process();
      assert(process);
      int n = solver_->getNumCols();
      const int *originalColumns = process->originalColumns();
      // columns should be in order ... but
      double *tempS = new double[n];
      for (int i = 0; i < n; i++) {
        int iColumn = originalColumns[i];
        tempS[i] = hotstartSolution_[iColumn];
      }
      delete[] hotstartSolution_;
      hotstartSolution_ = tempS;
      if (hotstartPriorities_) {
        int *tempP = new int[n];
        for (int i = 0; i < n; i++) {
          int iColumn = originalColumns[i];
          tempP[i] = hotstartPriorities_[iColumn];
        }
        delete[] hotstartPriorities_;
        hotstartPriorities_ = tempP;
      }
    }
    saveCompare = nodeCompare_;
    // depth first
    nodeCompare_ = new CbcCompareDepth();
  }
  if (!problemFeasibility_)
    problemFeasibility_ = new CbcFeasibilityBase();
#ifdef CBC_DEBUG
  std::string problemName;
  solver_->getStrParam(OsiProbName, problemName);
  printf("Problem name - %s\n", problemName.c_str());
  solver_->setHintParam(OsiDoReducePrint, false, OsiHintDo, 0);
#endif
  /*
      Assume we're done, and see if we're proven wrong.
    */
  status_ = 0;
  secondaryStatus_ = 0;
  phase_ = 0;
  /*
      Scan the variables, noting the integer variables. Create an
      CbcSimpleInteger object for each integer variable.
    */
  findIntegers(false);
  // Say not dynamic pseudo costs
  ownership_ &= ~0x40000000;
  // If dynamic pseudo costs then do
  if (numberBeforeTrust_)
    convertToDynamic();
  // Set up char array to say if integer (speed)
  delete[] integerInfo_;
  {
    int n = solver_->getNumCols();
    integerInfo_ = new char[n];
    for (int i = 0; i < n; i++) {
      if (solver_->isInteger(i))
        integerInfo_[i] = 1;
      else
        integerInfo_[i] = 0;
    }
  }
  if (preferredWay_) {
    // set all unset ones
    for (int iObject = 0; iObject < numberObjects_; iObject++) {
      CbcObject *obj = dynamic_cast<CbcObject *>(object_[iObject]);
      if (obj && !obj->preferredWay())
        obj->setPreferredWay(preferredWay_);
    }
  }
  /*
      Ensure that objects on the lists of OsiObjects, heuristics, and cut
      generators attached to this model all refer to this model.
    */
  synchronizeModel();
  if (!solverCharacteristics_) {
    OsiBabSolver *solverCharacteristics = dynamic_cast<OsiBabSolver *>(solver_->getAuxiliaryInfo());
    if (solverCharacteristics) {
      solverCharacteristics_ = solverCharacteristics;
    } else {
      // replace in solver
      OsiBabSolver defaultC;
      solver_->setAuxiliaryInfo(&defaultC);
      solverCharacteristics_ = dynamic_cast<OsiBabSolver *>(solver_->getAuxiliaryInfo());
    }
  }

  solverCharacteristics_->setSolver(solver_);
  // Set so we can tell we are in initial phase in resolve
  continuousObjective_ = -COIN_DBL_MAX;
  /*
      Solve the relaxation.

      Apparently there are circumstances where this will be non-trivial --- i.e.,
      we've done something since initialSolve that's trashed the solution to the
      continuous relaxation.
    */
  /* Tell solver we are in Branch and Cut
       Could use last parameter for subtle differences */
  solver_->setHintParam(OsiDoInBranchAndCut, true, OsiHintDo, NULL);
#ifdef COIN_HAS_CLP
  {
    OsiClpSolverInterface *clpSolver
      = dynamic_cast<OsiClpSolverInterface *>(solver_);
    if (clpSolver) {
      ClpSimplex *clpSimplex = clpSolver->getModelPtr();
      if ((specialOptions_ & 32) == 0) {
        // take off names (unless going to be saving)
        if (numberAnalyzeIterations_ >= 0 || (-numberAnalyzeIterations_ & 64) == 0)
          clpSimplex->dropNames();
      }
      // no crunch if mostly continuous
      if ((clpSolver->specialOptions() & (1 + 8)) != (1 + 8)) {
        int numberColumns = solver_->getNumCols();
        if (numberColumns > 1000 && numberIntegers_ * 4 < numberColumns)
          clpSolver->setSpecialOptions(clpSolver->specialOptions() & (~1));
      }
      //#define NO_CRUNCH
#ifdef NO_CRUNCH
      printf("TEMP switching off crunch\n");
      int iOpt = clpSolver->specialOptions();
      iOpt &= ~1;
      iOpt |= 65536;
      clpSolver->setSpecialOptions(iOpt);
#endif
    }
  }
#endif
  bool feasible;
  numberSolves_ = 0;
  {
    // check
    int numberOdd = 0;
    for (int i = 0; i < numberObjects_; i++) {
      CbcSimpleInteger *obj = dynamic_cast<CbcSimpleInteger *>(object_[i]);
      if (!obj)
        numberOdd++;
    }
    if (numberOdd)
      moreSpecialOptions_ |= 1073741824;
  }
  // If NLP then we assume already solved outside branchAndbound
  if (!solverCharacteristics_->solverType() || solverCharacteristics_->solverType() == 4) {
    feasible = resolve(NULL, 0) != 0;
  } else {
    // pick up given status
    feasible = (solver_->isProvenOptimal() && !solver_->isDualObjectiveLimitReached());
  }
  if (problemFeasibility_->feasible(this, 0) < 0) {
    feasible = false; // pretend infeasible
  }
  numberSavedSolutions_ = 0;
  int saveNumberStrong = numberStrong_;
  int saveNumberBeforeTrust = numberBeforeTrust_;
  /*
      If the linear relaxation of the root is infeasible, bail out now. Otherwise,
      continue with processing the root node.
    */
  if (!feasible) {
    status_ = 0;
    if (!solver_->isProvenDualInfeasible()) {
      handler_->message(CBC_INFEAS, messages_) << CoinMessageEol;
      secondaryStatus_ = 1;
    } else {
      handler_->message(CBC_UNBOUNDED, messages_) << CoinMessageEol;
      secondaryStatus_ = 7;
    }
    originalContinuousObjective_ = COIN_DBL_MAX;
    if (bestSolution_ && ((specialOptions_ & 8388608) == 0 || (specialOptions_ & 2048) != 0)) {
      // best solution found by various heuristics - set solution
      char general[200];
      sprintf(general, "Solution of %g already found by heuristic",
        bestObjective_);
      messageHandler()->message(CBC_GENERAL,
        messages())
        << general << CoinMessageEol;
      setCutoff(1.0e50); // As best solution should be worse than cutoff
      // change cutoff as constraint if wanted
      if (cutoffRowNumber_ >= 0) {
        if (solver_->getNumRows() > cutoffRowNumber_)
          solver_->setRowUpper(cutoffRowNumber_, 1.0e50);
      }
      // also in continuousSolver_
      if (continuousSolver_) {
        // Solvers know about direction
        double direction = solver_->getObjSense();
        continuousSolver_->setDblParam(OsiDualObjectiveLimit, 1.0e50 * direction);
      } else {
        continuousSolver_ = solver_->clone();
      }
      phase_ = 5;
      double increment = getDblParam(CbcModel::CbcCutoffIncrement);
      if ((specialOptions_ & 4) == 0)
        bestObjective_ += 100.0 * increment + 1.0e-3; // only set if we are going to solve
      setBestSolution(CBC_END_SOLUTION, bestObjective_, bestSolution_, 1);
      continuousSolver_->resolve();
      if (!continuousSolver_->isProvenOptimal()) {
        continuousSolver_->messageHandler()->setLogLevel(2);
        continuousSolver_->initialSolve();
      }
      delete solver_;
      solverCharacteristics_ = NULL;
      solver_ = continuousSolver_;
      setPointers(solver_);
      continuousSolver_ = NULL;
    }
    solverCharacteristics_ = NULL;
    if (flipObjective)
      flipModel();
    return;
  } else if (!numberObjects_) {
    // nothing to do
    // Undo preprocessing performed during BaB.
    if (strategy_ && strategy_->preProcessState() > 0) {
      // undo preprocessing
      CglPreProcess *process = strategy_->process();
      assert(process);
      int n = originalSolver->getNumCols();
      if (bestSolution_) {
        delete[] bestSolution_;
        bestSolution_ = new double[n];
        process->postProcess(*solver_);
      }
      strategy_->deletePreProcess();
      // Solution now back in originalSolver
      delete solver_;
      solver_ = originalSolver;
      if (bestSolution_) {
        bestObjective_ = solver_->getObjValue() * solver_->getObjSense();
        memcpy(bestSolution_, solver_->getColSolution(), n * sizeof(double));
      }
      // put back original objects if there were any
      if (originalObject) {
        int iColumn;
        assert(ownObjects_);
        for (iColumn = 0; iColumn < numberObjects_; iColumn++)
          delete object_[iColumn];
        delete[] object_;
        numberObjects_ = numberOriginalObjects;
        object_ = originalObject;
        delete[] integerVariable_;
        numberIntegers_ = 0;
        for (iColumn = 0; iColumn < n; iColumn++) {
          if (solver_->isInteger(iColumn))
            numberIntegers_++;
        }
        integerVariable_ = new int[numberIntegers_];
        numberIntegers_ = 0;
        for (iColumn = 0; iColumn < n; iColumn++) {
          if (solver_->isInteger(iColumn))
            integerVariable_[numberIntegers_++] = iColumn;
        }
      }
    }
    if (flipObjective)
      flipModel();
    solverCharacteristics_ = NULL;
    bestObjective_ = solver_->getObjValue() * solver_->getObjSense();
    int numberColumns = solver_->getNumCols();
    delete[] bestSolution_;
    bestSolution_ = new double[numberColumns];
    CoinCopyN(solver_->getColSolution(), numberColumns, bestSolution_);
    return;
  }
  /*
      See if we're using the Osi side of the branching hierarchy. If so, either
      convert existing CbcObjects to OsiObjects, or generate them fresh. In the
      first case, CbcModel owns the objects on the object_ list. In the second
      case, the solver holds the objects and object_ simply points to the
      solver's list.

      080417 The conversion code here (the block protected by `if (obj)') cannot
      possibly be correct. On the Osi side, descent is OsiObject -> OsiObject2 ->
      all other Osi object classes. On the Cbc side, it's OsiObject -> CbcObject
      -> all other Cbc object classes. It's structurally impossible for any Osi
      object to descend from CbcObject. The only thing I can see is that this is
      really dead code, and object detection is now handled from the Osi side.
    */
  // Convert to Osi if wanted
  //OsiBranchingInformation * persistentInfo = NULL;
  if (branchingMethod_ && branchingMethod_->chooseMethod()) {
    //persistentInfo = new OsiBranchingInformation(solver_);
    if (numberOriginalObjects) {
      for (int iObject = 0; iObject < numberObjects_; iObject++) {
        CbcObject *obj = dynamic_cast<CbcObject *>(object_[iObject]);
        if (obj) {
          CbcSimpleInteger *obj2 = dynamic_cast<CbcSimpleInteger *>(obj);
          if (obj2) {
            // back to Osi land
            object_[iObject] = obj2->osiObject();
            delete obj;
          } else {
            OsiSimpleInteger *obj3 = dynamic_cast<OsiSimpleInteger *>(obj);
            if (!obj3) {
              OsiSOS *obj4 = dynamic_cast<OsiSOS *>(obj);
              if (!obj4) {
                CbcSOS *obj5 = dynamic_cast<CbcSOS *>(obj);
                if (obj5) {
                  // back to Osi land
                  object_[iObject] = obj5->osiObject(solver_);
                } else {
                  printf("Code up CbcObject type in Osi land\n");
                  abort();
                }
              }
            }
          }
        }
      }
      // and add to solver
      //if (!solver_->numberObjects()) {
      solver_->addObjects(numberObjects_, object_);
      //} else {
      //if (solver_->numberObjects()!=numberOriginalObjects) {
      //printf("should have trapped that solver has objects before\n");
      //abort();
      //}
      //}
    } else {
      /*
              As of 080104, findIntegersAndSOS is misleading --- the default OSI
              implementation finds only integers.
            */
      // do from solver
      deleteObjects(false);
      solver_->findIntegersAndSOS(false);
      numberObjects_ = solver_->numberObjects();
      object_ = solver_->objects();
      ownObjects_ = false;
    }
    branchingMethod_->chooseMethod()->setSolver(solver_);
  }
  // take off heuristics if have to (some do not work with SOS, for example)
  // object should know what's safe.
  {
    int numberOdd = 0;
    int numberSOS = 0;
    for (int i = 0; i < numberObjects_; i++) {
      if (!object_[i]->canDoHeuristics())
        numberOdd++;
      CbcSOS *obj = dynamic_cast<CbcSOS *>(object_[i]);
      if (obj)
        numberSOS++;
    }
    if (numberOdd) {
      if (numberHeuristics_ && (specialOptions_ & 1024) == 0) {
        int k = 0;
        for (int i = 0; i < numberHeuristics_; i++) {
          if (!heuristic_[i]->canDealWithOdd())
            delete heuristic_[i];
          else
            heuristic_[k++] = heuristic_[i];
        }
        if (!k) {
          delete[] heuristic_;
          heuristic_ = NULL;
        }
        numberHeuristics_ = k;
        handler_->message(CBC_HEURISTICS_OFF, messages_) << numberOdd << CoinMessageEol;
      }
      // If odd switch off AddIntegers
      specialOptions_ &= ~65536;
      // switch off fast nodes for now
      fastNodeDepth_ = -1;
      moreSpecialOptions_ &= ~33554432; // no diving
    } else if (numberSOS) {
      specialOptions_ |= 128; // say can do SOS in dynamic mode
      // switch off fast nodes for now
      fastNodeDepth_ = -1;
      moreSpecialOptions_ &= ~33554432; // no diving
    }
    if (numberThreads_ > 0) {
      /* switch off fast nodes for now
               Trouble is that by time mini bab finishes code is
               looking at a different node
             */
      fastNodeDepth_ = -1;
    }
  }
  // Save objective (just so user can access it)
  originalContinuousObjective_ = solver_->getObjValue() * solver_->getObjSense();
  bestPossibleObjective_ = originalContinuousObjective_;
  sumChangeObjective1_ = 0.0;
  sumChangeObjective2_ = 0.0;
  /*
      OsiRowCutDebugger knows an optimal answer for a subset of MIP problems.
      Assuming it recognises the problem, when called upon it will check a cut to
      see if it cuts off the optimal answer.
    */
  // If debugger exists set specialOptions_ bit
  if (solver_->getRowCutDebuggerAlways()) {
    specialOptions_ |= 1;
  }

#ifdef CBC_DEBUG
  if ((specialOptions_ & 1) == 0)
    solver_->activateRowCutDebugger(problemName.c_str());
  if (solver_->getRowCutDebuggerAlways())
    specialOptions_ |= 1;
#endif

  /*
      Begin setup to process a feasible root node.
    */
  bestObjective_ = CoinMin(bestObjective_, 1.0e50);
  if (!bestSolution_) {
    numberSolutions_ = 0;
    numberHeuristicSolutions_ = 0;
  }
  stateOfSearch_ = 0;
  // Everything is minimization
  {
    // needed to sync cutoffs
    double value;
    solver_->getDblParam(OsiDualObjectiveLimit, value);
    dblParam_[CbcCurrentCutoff] = value * solver_->getObjSense();
  }
  double cutoff = getCutoff();
  double direction = solver_->getObjSense();
  dblParam_[CbcOptimizationDirection] = direction;
  if (cutoff < 1.0e20 && direction < 0.0)
    messageHandler()->message(CBC_CUTOFF_WARNING1,
      messages())
      << cutoff << -cutoff << CoinMessageEol;
  if (cutoff > bestObjective_)
    cutoff = bestObjective_;
  setCutoff(cutoff);
  /*
      We probably already have a current solution, but just in case ...
    */
  int numberColumns = getNumCols();
  if (!currentSolution_)
    currentSolution_ = new double[numberColumns];
  testSolution_ = currentSolution_;
  /*
      Create a copy of the solver, thus capturing the original (root node)
      constraint system (aka the continuous system).
    */
  delete continuousSolver_;
  continuousSolver_ = solver_->clone();
#ifdef CONFLICT_CUTS
  if ((moreSpecialOptions_ & 4194304) != 0) {
#ifdef COIN_HAS_CLP
    OsiClpSolverInterface *clpSolver
      = dynamic_cast<OsiClpSolverInterface *>(solver_);
    if (clpSolver) {
      int specialOptions = clpSolver->getModelPtr()->specialOptions();
      // 2097152 switches on rays in crunch
      if (!parentModel_)
        clpSolver->getModelPtr()->setSpecialOptions(specialOptions | 32 | 2097152);
      else
        clpSolver->getModelPtr()->setSpecialOptions(specialOptions & ~(32 | 2097152));
    }
  }
#endif
#endif
#ifdef COIN_HAS_NTY
  // maybe allow on fix and restart later
  if ((moreSpecialOptions2_ & (128 | 256)) != 0 && !parentModel_) {
    symmetryInfo_ = new CbcSymmetry();
    symmetryInfo_->setupSymmetry(*continuousSolver_);
    int numberGenerators = symmetryInfo_->statsOrbits(this, 0);
    if (!symmetryInfo_->numberUsefulOrbits() && (moreSpecialOptions2_ & (128 | 256)) != (128 | 256)) {
      delete symmetryInfo_;
      symmetryInfo_ = NULL;
      moreSpecialOptions2_ &= ~(128 | 256);
    }
    if ((moreSpecialOptions2_ & (128 | 256)) == (128 | 256)) {
      //moreSpecialOptions2_ &= ~256;
    }
  }
#endif

  // add cutoff as constraint if wanted
  if (cutoffRowNumber_ == -2) {
    if (!parentModel_) {
      int numberColumns = solver_->getNumCols();
      double *obj = CoinCopyOfArray(solver_->getObjCoefficients(), numberColumns);
      int *indices = new int[numberColumns];
      int n = 0;
      for (int i = 0; i < numberColumns; i++) {
        if (obj[i]) {
          indices[n] = i;
          obj[n++] = obj[i];
        }
      }
      if (n) {
        double cutoff = getCutoff();
        // relax a little bit
        cutoff += 1.0e-4;
        double offset;
        solver_->getDblParam(OsiObjOffset, offset);
        cutoffRowNumber_ = solver_->getNumRows();
        solver_->addRow(n, indices, obj, -COIN_DBL_MAX, CoinMin(cutoff, 1.0e25) + offset);
      } else {
        // no objective!
        cutoffRowNumber_ = -1;
      }
      delete[] indices;
      delete[] obj;
    } else {
      // switch off
      cutoffRowNumber_ = -1;
    }
  }
  numberRowsAtContinuous_ = getNumRows();
  solver_->saveBaseModel();
  /*
      Check the objective to see if we can deduce a nontrivial increment. If
      it's better than the current value for CbcCutoffIncrement, it'll be
      installed.
    */
  if (solverCharacteristics_->reducedCostsAccurate())
    analyzeObjective();
  {
    // may be able to change cutoff now
    double cutoff = getCutoff();
    double increment = getDblParam(CbcModel::CbcCutoffIncrement);
    if (cutoff > bestObjective_ - increment) {
      cutoff = bestObjective_ - increment;
      setCutoff(cutoff);
    }
  }
#ifdef COIN_HAS_CLP
  // Possible save of pivot method
  ClpDualRowPivot *savePivotMethod = NULL;
  {
    // pass tolerance and increment to solver
    OsiClpSolverInterface *clpSolver
      = dynamic_cast<OsiClpSolverInterface *>(solver_);
    if (clpSolver)
      clpSolver->setStuff(getIntegerTolerance(), getCutoffIncrement());
#ifdef CLP_RESOLVE
    if ((moreSpecialOptions_ & 1048576) != 0 && !parentModel_ && clpSolver) {
      resolveClp(clpSolver, 0);
    }
#endif
  }
#endif
  /*
      Set up for cut generation. addedCuts_ holds the cuts which are relevant for
      the active subproblem. whichGenerator will be used to record the generator
      that produced a given cut.
    */
#define INITIAL_MAXIMUM_WHICH 1000
  maximumWhich_ = INITIAL_MAXIMUM_WHICH;
  delete[] whichGenerator_;
  whichGenerator_ = new int[maximumWhich_];
  memset(whichGenerator_, 0, maximumWhich_ * sizeof(int));
  maximumNumberCuts_ = 0;
  currentNumberCuts_ = 0;
  delete[] addedCuts_;
  addedCuts_ = NULL;
  OsiObject **saveObjects = NULL;
  maximumRows_ = numberRowsAtContinuous_;
  currentDepth_ = 0;
  workingBasis_.resize(maximumRows_, numberColumns);
  /*
      Set up an empty heap and associated data structures to hold the live set
      (problems which require further exploration).
    */
  CbcCompareDefault *compareActual
    = dynamic_cast<CbcCompareDefault *>(nodeCompare_);
  if (compareActual) {
    compareActual->setBestPossible(direction * solver_->getObjValue());
    compareActual->setCutoff(getCutoff());
#ifdef JJF_ZERO
    if (false && !numberThreads_ && !parentModel_) {
      printf("CbcTreeArray ? threads ? parentArray\n");
      // Setup new style tree
      delete tree_;
      tree_ = new CbcTreeArray();
    }
#endif
  }
  tree_->setComparison(*nodeCompare_);
  /*
      Used to record the path from a node to the root of the search tree, so that
      we can then traverse from the root to the node when restoring a subproblem.
    */
  maximumDepth_ = 10;
  delete[] walkback_;
  walkback_ = new CbcNodeInfo *[maximumDepth_];
  lastDepth_ = 0;
  delete[] lastNodeInfo_;
  lastNodeInfo_ = new CbcNodeInfo *[maximumDepth_];
  delete[] lastNumberCuts_;
  lastNumberCuts_ = new int[maximumDepth_];
  maximumCuts_ = 100;
  lastNumberCuts2_ = 0;
  delete[] lastCut_;
  lastCut_ = new const OsiRowCut *[maximumCuts_];
  /*
      Used to generate bound edits for CbcPartialNodeInfo.
    */
  double *lowerBefore = new double[numberColumns];
  double *upperBefore = new double[numberColumns];
  /*
    Set up to run heuristics and generate cuts at the root node. The heavy
    lifting is hidden inside the calls to doHeuristicsAtRoot and solveWithCuts.

    To start, tell cut generators they can be a bit more aggressive at the
    root node.

    QUESTION: phase_ = 0 is documented as `initial solve', phase = 1 as `solve
        with cuts at root'. Is phase_ = 1 the correct indication when
        doHeurisiticsAtRoot is called to run heuristics outside of the main
        cut / heurisitc / reoptimise loop in solveWithCuts?

      Generate cuts at the root node and reoptimise. solveWithCuts does the heavy
      lifting. It will iterate a generate/reoptimise loop (including reduced cost
      fixing) until no cuts are generated, the change in objective falls off,  or
      the limit on the number of rounds of cut generation is exceeded.

      At the end of all this, any cuts will be recorded in cuts and also
      installed in the solver's constraint system. We'll have reoptimised, and
      removed any slack cuts (numberOldActiveCuts_ and numberNewCuts_ have been
      adjusted accordingly).

      Tell cut generators they can be a bit more aggressive at root node

      TODO: Why don't we make a copy of the solution after solveWithCuts?
      TODO: If numberUnsatisfied == 0, don't we have a solution?
    */
  phase_ = 1;
  int iCutGenerator;
  for (iCutGenerator = 0; iCutGenerator < numberCutGenerators_; iCutGenerator++) {
    // If parallel switch off global cuts
    if (numberThreads_) {
      generator_[iCutGenerator]->setGlobalCuts(false);
      generator_[iCutGenerator]->setGlobalCutsAtRoot(false);
    }
    CglCutGenerator *generator = generator_[iCutGenerator]->generator();
    generator->setAggressiveness(generator->getAggressiveness() + 100);
    if (!generator->canDoGlobalCuts())
      generator->setGlobalCuts(false);
  }
  OsiCuts cuts;
  int anyAction = -1;
  numberOldActiveCuts_ = 0;
  numberNewCuts_ = 0;
  // Array to mark solution
  delete[] usedInSolution_;
  usedInSolution_ = new int[numberColumns];
  CoinZeroN(usedInSolution_, numberColumns);
  /*
      For printing totals and for CbcNode (numberNodes_)
    */
  numberIterations_ = 0;
  numberNodes_ = 0;
  numberNodes2_ = 0;
  maximumStatistics_ = 0;
  maximumDepthActual_ = 0;
  numberDJFixed_ = 0.0;
  if (!parentModel_) {
    if ((specialOptions_ & 262144) != 0) {
      // create empty stored cuts
      //storedRowCuts_ = new CglStored(solver_->getNumCols());
    } else if ((specialOptions_ & 524288) != 0 && storedRowCuts_) {
      // tighten and set best solution
      // A) tight bounds on integer variables
      /*
                storedRowCuts_ are coming in from outside, probably for nonlinear.
              John was unsure about origin.
            */
      const double *lower = solver_->getColLower();
      const double *upper = solver_->getColUpper();
      const double *tightLower = storedRowCuts_->tightLower();
      const double *tightUpper = storedRowCuts_->tightUpper();
      int nTightened = 0;
      for (int i = 0; i < numberIntegers_; i++) {
        int iColumn = integerVariable_[i];
        if (tightLower[iColumn] > lower[iColumn]) {
          nTightened++;
          solver_->setColLower(iColumn, tightLower[iColumn]);
        }
        if (tightUpper[iColumn] < upper[iColumn]) {
          nTightened++;
          solver_->setColUpper(iColumn, tightUpper[iColumn]);
        }
      }
      if (nTightened)
        COIN_DETAIL_PRINT(printf("%d tightened by alternate cuts\n", nTightened));
      if (storedRowCuts_->bestObjective() < bestObjective_) {
        // B) best solution
        double objValue = storedRowCuts_->bestObjective();
        setBestSolution(CBC_SOLUTION, objValue,
          storedRowCuts_->bestSolution());
        // Do heuristics
        // Allow RINS
        for (int i = 0; i < numberHeuristics_; i++) {
          CbcHeuristicRINS *rins
            = dynamic_cast<CbcHeuristicRINS *>(heuristic_[i]);
          if (rins) {
            rins->setLastNode(-100);
          }
        }
      }
    }
  }
#ifdef SWITCH_VARIABLES
  // see if any switching variables
  if (numberIntegers_ < solver_->getNumCols())
    findSwitching();
#endif
  /*
      Run heuristics at the root. This is the only opportunity to run FPump; it
      will be removed from the heuristics list by doHeuristicsAtRoot.
    */
  // See if multiple runs wanted
  CbcModel **rootModels = NULL;
  if (!parentModel_ && multipleRootTries_ % 100) {
    double rootTimeCpu = CoinCpuTime();
    double startTimeRoot = CoinGetTimeOfDay();
    int numberRootThreads = 1;
    /* undocumented fine tuning
         aabbcc where cc is number of tries
         bb if nonzero is number of threads
         aa if nonzero just do heuristics
      */
    int numberModels = multipleRootTries_ % 100;
#ifdef CBC_THREAD
    numberRootThreads = (multipleRootTries_ / 100) % 100;
    if (!numberRootThreads) {
      if (numberThreads_ < 2)
        numberRootThreads = numberModels;
      else
        numberRootThreads = CoinMin(numberThreads_, numberModels);
    }
#endif
    int otherOptions = (multipleRootTries_ / 10000) % 100;
    rootModels = new CbcModel *[numberModels];
    unsigned int newSeed = randomSeed_;
    if (newSeed == 0) {
      double time = fabs(CoinGetTimeOfDay());
      while (time >= COIN_INT_MAX)
        time *= 0.5;
      newSeed = static_cast<unsigned int>(time);
    } else if (newSeed < 0) {
      newSeed = 123456789;
#ifdef COIN_HAS_CLP
      OsiClpSolverInterface *clpSolver
        = dynamic_cast<OsiClpSolverInterface *>(solver_);
      if (clpSolver) {
        newSeed += clpSolver->getModelPtr()->randomNumberGenerator()->getSeed();
      }
#endif
    }
    CoinWarmStartBasis *basis = dynamic_cast<CoinWarmStartBasis *>(solver_->getEmptyWarmStart());
    for (int i = 0; i < numberModels; i++) {
      rootModels[i] = new CbcModel(*this);
      rootModels[i]->setNumberThreads(0);
      rootModels[i]->setMaximumNodes(otherOptions ? -1 : 0);
      rootModels[i]->setRandomSeed(newSeed + 10000000 * i);
      rootModels[i]->randomNumberGenerator()->setSeed(newSeed + 50000000 * i);
      rootModels[i]->setMultipleRootTries(0);
#ifdef COIN_HAS_NTY
      rootModels[i]->zapSymmetry();
      rootModels[i]->moreSpecialOptions2_ &= ~(128 | 256); // off nauty
#endif
      // use seed
      rootModels[i]->setSpecialOptions(specialOptions_ | (4194304 | 8388608));
      rootModels[i]->setMoreSpecialOptions(moreSpecialOptions_ & (~(134217728 | 4194304)));
      rootModels[i]->setMoreSpecialOptions2(moreSpecialOptions2_ & (~(128 | 256)));
      rootModels[i]->solver_->setWarmStart(basis);
#ifdef COIN_HAS_CLP
      OsiClpSolverInterface *clpSolver
        = dynamic_cast<OsiClpSolverInterface *>(rootModels[i]->solver_);
#define NEW_RANDOM_BASIS
#ifdef NEW_RANDOM_BASIS
      if (i == 0)
        continue;
#endif
      if (clpSolver) {
        ClpSimplex *simplex = clpSolver->getModelPtr();
        if (defaultHandler_)
          simplex->setDefaultMessageHandler();
        simplex->setRandomSeed(newSeed + 20000000 * i);
        simplex->allSlackBasis();
        int logLevel = simplex->logLevel();
        if (logLevel == 1)
          simplex->setLogLevel(0);
        if (i != 0) {
#ifdef NEW_RANDOM_BASIS
          int numberRows = simplex->numberRows();
          int throwOut = 20; //2+numberRows/100;
          for (int iThrow = 0; iThrow < throwOut; iThrow++) {
            double random = simplex->randomNumberGenerator()->randomDouble();
            int iStart = static_cast<int>(random * numberRows);
            for (int j = iStart; j < numberRows; j++) {
              if (simplex->getRowStatus(j) != ClpSimplex::basic) {
                simplex->setRowStatus(j, ClpSimplex::basic);
                break;
              }
            }
          }
          clpSolver->setWarmStart(NULL);
#else
            double random = simplex->randomNumberGenerator()->randomDouble();
            int bias = static_cast<int>(random * (numberIterations / 4));
            simplex->setMaximumIterations(numberIterations / 2 + bias);
            simplex->primal();
            simplex->setMaximumIterations(COIN_INT_MAX);
            simplex->dual();
#endif
        } else {
#ifndef NEW_RANDOM_BASIS
          simplex->primal();
#endif
#endif
        }
#ifdef NEW_RANDOM_BASIS
        simplex->setLogLevel(logLevel);
        clpSolver->setWarmStart(NULL);
#endif
      }
      for (int j = 0; j < numberHeuristics_; j++)
        rootModels[i]->heuristic_[j]->setSeed(rootModels[i]->heuristic_[j]->getSeed() + 100000000 * i);
      for (int j = 0; j < numberCutGenerators_; j++)
        rootModels[i]->generator_[j]->generator()->refreshSolver(rootModels[i]->solver_);
    }
    delete basis;
#ifdef CBC_THREAD
    if (numberRootThreads == 1) {
#endif
      for (int iModel = 0; iModel < numberModels; iModel++) {
        doRootCbcThread(rootModels[iModel]);
        // see if solved at root node
        if (rootModels[iModel]->getMaximumNodes()) {
          feasible = false;
          break;
        }
      }
#ifdef CBC_THREAD
    } else {
      Coin_pthread_t *threadId = new Coin_pthread_t[numberRootThreads];
      for (int kModel = 0; kModel < numberModels; kModel += numberRootThreads) {
        bool finished = false;
        for (int iModel = kModel; iModel < CoinMin(numberModels, kModel + numberRootThreads); iModel++) {
          pthread_create(&(threadId[iModel - kModel].thr), NULL,
            doRootCbcThread,
            rootModels[iModel]);
        }
        // wait
        for (int iModel = kModel; iModel < CoinMin(numberModels, kModel + numberRootThreads); iModel++) {
          pthread_join(threadId[iModel - kModel].thr, NULL);
        }
        // see if solved at root node
        for (int iModel = kModel; iModel < CoinMin(numberModels, kModel + numberRootThreads); iModel++) {
          if (rootModels[iModel]->getMaximumNodes())
            finished = true;
        }
        if (finished) {
          feasible = false;
          break;
        }
      }
      delete[] threadId;
    }
#endif
    // sort solutions
    int *which = new int[numberModels];
    double *value = new double[numberModels];
    int numberSolutions = 0;
    for (int iModel = 0; iModel < numberModels; iModel++) {
      if (rootModels[iModel]->bestSolution()) {
        which[numberSolutions] = iModel;
        value[numberSolutions++] = -rootModels[iModel]->getMinimizationObjValue();
      }
    }
    char general[100];
    rootTimeCpu = CoinCpuTime() - rootTimeCpu;
    if (numberRootThreads == 1)
      sprintf(general, "Multiple root solvers took a total of %.2f seconds\n",
        rootTimeCpu);
    else
      sprintf(general, "Multiple root solvers took a total of %.2f seconds (%.2f elapsed)\n",
        rootTimeCpu, CoinGetTimeOfDay() - startTimeRoot);
    messageHandler()->message(CBC_GENERAL,
      messages())
      << general << CoinMessageEol;
    CoinSort_2(value, value + numberSolutions, which);
    // to get name
    CbcHeuristicRINS dummyHeuristic;
    dummyHeuristic.setHeuristicName("Multiple root solvers");
    lastHeuristic_ = &dummyHeuristic;
    for (int i = 0; i < numberSolutions; i++) {
      double objValue = -value[i];
      if (objValue < getCutoff()) {
        int iModel = which[i];
        setBestSolution(CBC_ROUNDING, objValue,
          rootModels[iModel]->bestSolution());
      }
    }
    lastHeuristic_ = NULL;
    delete[] which;
    delete[] value;
  }
  // Do heuristics
  if (numberObjects_ && !rootModels)
    doHeuristicsAtRoot();
  if (solverCharacteristics_->solutionAddsCuts()) {
    // With some heuristics solver needs a resolve here
    solver_->resolve();
    if (!isProvenOptimal()) {
      solver_->initialSolve();
    }
  }
  /*
      Grepping through the code, it would appear that this is a command line
      debugging hook.  There's no obvious place in the code where this is set to
      a negative value.

      User hook, says John.
    */
  if (intParam_[CbcMaxNumNode] < 0
    || numberSolutions_ >= getMaximumSolutions())
    eventHappened_ = true; // stop as fast as possible
  stoppedOnGap_ = false;
  // See if can stop on gap
  bestPossibleObjective_ = solver_->getObjValue() * solver_->getObjSense();
  if (canStopOnGap()) {
    if (bestPossibleObjective_ < getCutoff())
      stoppedOnGap_ = true;
    feasible = false;
    //eventHappened_=true; // stop as fast as possible
  }
  /*
      Set up for statistics collection, if requested. Standard values are
      documented in CbcModel.hpp. The magic number 100 will trigger a dump of
      CbcSimpleIntegerDynamicPseudoCost objects (no others). Looks like another
      command line debugging hook.
    */
  statistics_ = NULL;
  // Do on switch
  if (doStatistics > 0 && doStatistics <= 100) {
    maximumStatistics_ = 10000;
    statistics_ = new CbcStatistics *[maximumStatistics_];
    memset(statistics_, 0, maximumStatistics_ * sizeof(CbcStatistics *));
  }
  // See if we can add integers
  if (noObjects && numberIntegers_ < solver_->getNumCols() && (specialOptions_ & 65536) != 0 && !parentModel_ && false) {
    int numberIntegers1 = 0;
    int numberColumns = solver_->getNumCols();
    for (int i = 0; i < numberColumns; i++) {
      if (solver_->isInteger(i))
        numberIntegers1++;
    }
    AddIntegers();
    // make sure in sync
    int numberIntegers2 = 0;
    for (int i = 0; i < numberColumns; i++) {
      if (solver_->isInteger(i))
        numberIntegers2++;
    }
    if (numberIntegers1 < numberIntegers2) {
      findIntegers(true, 2);
      convertToDynamic();
    }
  }

  /*
      Do an initial round of cut generation for the root node. Depending on the
      type of underlying solver, we may want to do this even if the initial query
      to the objects indicates they're satisfied.

      solveWithCuts does the heavy lifting. It will iterate a generate/reoptimise
      loop (including reduced cost fixing) until no cuts are generated, the
      change in objective falls off,  or the limit on the number of rounds of cut
      generation is exceeded.

      At the end of all this, any cuts will be recorded in cuts and also
      installed in the solver's constraint system. We'll have reoptimised, and
      removed any slack cuts (numberOldActiveCuts_ and numberNewCuts_ have been
      adjusted accordingly).
    */
  int iObject;
  int numberUnsatisfied = 0;
  delete[] currentSolution_;
  currentSolution_ = new double[numberColumns];
  testSolution_ = currentSolution_;
  memcpy(currentSolution_, solver_->getColSolution(),
    numberColumns * sizeof(double));
  // point to useful information
  OsiBranchingInformation usefulInfo = usefulInformation();

  for (iObject = 0; iObject < numberObjects_; iObject++) {
    double infeasibility = object_[iObject]->checkInfeasibility(&usefulInfo);
    if (infeasibility)
      numberUnsatisfied++;
  }
  // replace solverType
  double *tightBounds = NULL;
  if (solverCharacteristics_->tryCuts()) {

    if (numberUnsatisfied) {
      // User event
      if (!eventHappened_ && feasible) {
        if (rootModels) {
          // for fixings
          int numberColumns = solver_->getNumCols();
          tightBounds = new double[2 * numberColumns];
          {
            const double *lower = solver_->getColLower();
            const double *upper = solver_->getColUpper();
            for (int i = 0; i < numberColumns; i++) {
              tightBounds[2 * i + 0] = lower[i];
              tightBounds[2 * i + 1] = upper[i];
            }
          }
          int numberModels = multipleRootTries_ % 100;
          const OsiSolverInterface **solvers = new const OsiSolverInterface *[numberModels];
          int numberRows = continuousSolver_->getNumRows();
          int maxCuts = 0;
          for (int i = 0; i < numberModels; i++) {
            solvers[i] = rootModels[i]->solver();
            const double *lower = solvers[i]->getColLower();
            const double *upper = solvers[i]->getColUpper();
            for (int j = 0; j < numberColumns; j++) {
              tightBounds[2 * j + 0] = CoinMax(lower[j], tightBounds[2 * j + 0]);
              tightBounds[2 * j + 1] = CoinMin(upper[j], tightBounds[2 * j + 1]);
            }
            int numberRows2 = solvers[i]->getNumRows();
            assert(numberRows2 >= numberRows);
            maxCuts += numberRows2 - numberRows;
            // accumulate statistics
            for (int j = 0; j < numberCutGenerators_; j++) {
              generator_[j]->addStatistics(rootModels[i]->cutGenerator(j));
            }
          }
          for (int j = 0; j < numberCutGenerators_; j++) {
            generator_[j]->scaleBackStatistics(numberModels);
          }
          //CbcRowCuts rowCut(maxCuts);
          const OsiRowCutDebugger *debugger = NULL;
          if ((specialOptions_ & 1) != 0)
            debugger = solver_->getRowCutDebugger();
          for (int iModel = 0; iModel < numberModels; iModel++) {
            int numberRows2 = solvers[iModel]->getNumRows();
            const CoinPackedMatrix *rowCopy = solvers[iModel]->getMatrixByRow();
            const int *rowLength = rowCopy->getVectorLengths();
            const double *elements = rowCopy->getElements();
            const int *column = rowCopy->getIndices();
            const CoinBigIndex *rowStart = rowCopy->getVectorStarts();
            const double *rowLower = solvers[iModel]->getRowLower();
            const double *rowUpper = solvers[iModel]->getRowUpper();
            for (int iRow = numberRows; iRow < numberRows2; iRow++) {
              OsiRowCut rc;
              rc.setLb(rowLower[iRow]);
              rc.setUb(rowUpper[iRow]);
              CoinBigIndex start = rowStart[iRow];
              rc.setRow(rowLength[iRow], column + start, elements + start, false);
              if (debugger)
                CoinAssert(!debugger->invalidCut(rc));
              globalCuts_.addCutIfNotDuplicate(rc);
            }
            //int cutsAdded=globalCuts_.numberCuts()-numberCuts;
            //numberCuts += cutsAdded;
            //printf("Model %d gave %d cuts (out of %d possible)\n",
            //     iModel,cutsAdded,numberRows2-numberRows);
          }
          // normally replace global cuts
          //if (!globalCuts_.())
          //globalCuts_=rowCutrowCut.addCuts(globalCuts_);
          //rowCut.addCuts(globalCuts_);
          int nTightened = 0;
          assert(feasible);
          {
            double tolerance = 1.0e-5;
            const double *lower = solver_->getColLower();
            const double *upper = solver_->getColUpper();
            for (int i = 0; i < numberColumns; i++) {
              if (tightBounds[2 * i + 0] > tightBounds[2 * i + 1] + 1.0e-9) {
                feasible = false;
                char general[200];
                sprintf(general, "Solvers give infeasible bounds on %d %g,%g was %g,%g - search finished\n",
                  i, tightBounds[2 * i + 0], tightBounds[2 * i + 1], lower[i], upper[i]);
                messageHandler()->message(CBC_GENERAL, messages())
                  << general << CoinMessageEol;
                break;
              }
              double oldLower = lower[i];
              double oldUpper = upper[i];
              if (tightBounds[2 * i + 0] > oldLower + tolerance) {
                nTightened++;
                solver_->setColLower(i, tightBounds[2 * i + 0]);
              }
              if (tightBounds[2 * i + 1] < oldUpper - tolerance) {
                nTightened++;
                solver_->setColUpper(i, tightBounds[2 * i + 1]);
              }
            }
          }
          delete[] tightBounds;
          tightBounds = NULL;
          char printBuffer[200];
          sprintf(printBuffer, "%d solvers added %d different cuts out of pool of %d",
            numberModels, globalCuts_.sizeRowCuts(), maxCuts);
          messageHandler()->message(CBC_GENERAL, messages())
            << printBuffer << CoinMessageEol;
          if (nTightened) {
            sprintf(printBuffer, "%d bounds were tightened",
              nTightened);
            messageHandler()->message(CBC_GENERAL, messages())
              << printBuffer << CoinMessageEol;
          }
          delete[] solvers;
        }
        if (!parentModel_ && (moreSpecialOptions_ & 67108864) != 0) {
          // load cuts from file
          FILE *fp = fopen("global.cuts", "rb");
          if (fp) {
            size_t nRead;
            int numberColumns = solver_->getNumCols();
            int numCols;
            nRead = fread(&numCols, sizeof(int), 1, fp);
            if (nRead != 1)
              throw("Error in fread");
            if (numberColumns != numCols) {
              printf("Mismatch on columns %d %d\n", numberColumns, numCols);
              fclose(fp);
            } else {
              // If rootModel just do some
              double threshold = -1.0;
              if (!multipleRootTries_)
                threshold = 0.5;
              int initialCuts = 0;
              int initialGlobal = globalCuts_.sizeRowCuts();
              double *elements = new double[numberColumns + 2];
              int *indices = new int[numberColumns];
              int numberEntries = 1;
              while (numberEntries > 0) {
                nRead = fread(&numberEntries, sizeof(int), 1, fp);
                if (nRead != 1)
                  throw("Error in fread");
                double randomNumber = randomNumberGenerator_.randomDouble();
                if (numberEntries > 0) {
                  initialCuts++;
                  nRead = fread(elements, sizeof(double), numberEntries + 2, fp);
                  if (nRead != static_cast<size_t>(numberEntries + 2))
                    throw("Error in fread");
                  nRead = fread(indices, sizeof(int), numberEntries, fp);
                  if (nRead != static_cast<size_t>(numberEntries))
                    throw("Error in fread");
                  if (randomNumber > threshold) {
                    OsiRowCut rc;
                    rc.setLb(elements[numberEntries]);
                    rc.setUb(elements[numberEntries + 1]);
                    rc.setRow(numberEntries, indices, elements,
                      false);
                    rc.setGloballyValidAsInteger(2);
                    globalCuts_.addCutIfNotDuplicate(rc);
                  }
                }
              }
              fclose(fp);
              // fixes
              int nTightened = 0;
              fp = fopen("global.fix", "rb");
              if (fp) {
                nRead = fread(indices, sizeof(int), 2, fp);
                if (nRead != 2)
                  throw("Error in fread");
                if (numberColumns != indices[0]) {
                  printf("Mismatch on columns %d %d\n", numberColumns,
                    indices[0]);
                } else {
                  indices[0] = 1;
                  while (indices[0] >= 0) {
                    nRead = fread(indices, sizeof(int), 2, fp);
                    if (nRead != 2)
                      throw("Error in fread");
                    int iColumn = indices[0];
                    if (iColumn >= 0) {
                      nTightened++;
                      nRead = fread(elements, sizeof(double), 4, fp);
                      if (nRead != 4)
                        throw("Error in fread");
                      solver_->setColLower(iColumn, elements[0]);
                      solver_->setColUpper(iColumn, elements[1]);
                    }
                  }
                }
              }
              if (fp)
                fclose(fp);
              char printBuffer[200];
              sprintf(printBuffer, "%d cuts read in of which %d were unique, %d bounds tightened",
                initialCuts,
                globalCuts_.sizeRowCuts() - initialGlobal, nTightened);
              messageHandler()->message(CBC_GENERAL, messages())
                << printBuffer << CoinMessageEol;
              delete[] elements;
              delete[] indices;
            }
          }
        }
        if (feasible)
          feasible = solveWithCuts(cuts, maximumCutPassesAtRoot_,
            NULL);
        if (multipleRootTries_ && (moreSpecialOptions_ & 134217728) != 0) {
          FILE *fp = NULL;
          size_t nRead;
          int numberColumns = solver_->getNumCols();
          int initialCuts = 0;
          if ((moreSpecialOptions_ & 134217728) != 0) {
            // append so go down to end
            fp = fopen("global.cuts", "r+b");
            if (fp) {
              int numCols;
              nRead = fread(&numCols, sizeof(int), 1, fp);
              if (nRead != 1)
                throw("Error in fread");
              if (numberColumns != numCols) {
                printf("Mismatch on columns %d %d\n", numberColumns, numCols);
                fclose(fp);
                fp = NULL;
              }
            }
          }
          double *elements = new double[numberColumns + 2];
          int *indices = new int[numberColumns];
          if (fp) {
            int numberEntries = 1;
            while (numberEntries > 0) {
              fpos_t position;
              fgetpos(fp, &position);
              nRead = fread(&numberEntries, sizeof(int), 1, fp);
              if (nRead != 1)
                throw("Error in fread");
              if (numberEntries > 0) {
                initialCuts++;
                nRead = fread(elements, sizeof(double), numberEntries + 2, fp);
                if (nRead != static_cast<size_t>(numberEntries + 2))
                  throw("Error in fread");
                nRead = fread(indices, sizeof(int), numberEntries, fp);
                if (nRead != static_cast<size_t>(numberEntries))
                  throw("Error in fread");
              } else {
                // end
                fsetpos(fp, &position);
              }
            }
          } else {
            fp = fopen("global.cuts", "wb");
            size_t nWrite;
            nWrite = fwrite(&numberColumns, sizeof(int), 1, fp);
            if (nWrite != 1)
              throw("Error in fwrite");
          }
          size_t nWrite;
          // now append binding cuts
          int numberC = continuousSolver_->getNumRows();
          int numberRows = solver_->getNumRows();
          printf("Saving %d cuts (up from %d)\n",
            initialCuts + numberRows - numberC, initialCuts);
          const double *rowLower = solver_->getRowLower();
          const double *rowUpper = solver_->getRowUpper();
          // Row copy
          CoinPackedMatrix matrixByRow(*solver_->getMatrixByRow());
          const double *elementByRow = matrixByRow.getElements();
          const int *column = matrixByRow.getIndices();
          const CoinBigIndex *rowStart = matrixByRow.getVectorStarts();
          const int *rowLength = matrixByRow.getVectorLengths();
          for (int iRow = numberC; iRow < numberRows; iRow++) {
            int n = rowLength[iRow];
            assert(n);
            CoinBigIndex start = rowStart[iRow];
            memcpy(elements, elementByRow + start, n * sizeof(double));
            memcpy(indices, column + start, n * sizeof(int));
            elements[n] = rowLower[iRow];
            elements[n + 1] = rowUpper[iRow];
            nWrite = fwrite(&n, sizeof(int), 1, fp);
            if (nWrite != 1)
              throw("Error in fwrite");
            nWrite = fwrite(elements, sizeof(double), n + 2, fp);
            if (nWrite != static_cast<size_t>(n + 2))
              throw("Error in fwrite");
            nWrite = fwrite(indices, sizeof(int), n, fp);
            if (nWrite != static_cast<size_t>(n))
              throw("Error in fwrite");
          }
          // eof marker
          int eofMarker = -1;
          nWrite = fwrite(&eofMarker, sizeof(int), 1, fp);
          if (nWrite != 1)
            throw("Error in fwrite");
          fclose(fp);
          // do tighter bounds (? later extra to original columns)
          int nTightened = 0;
          const double *lower = solver_->getColLower();
          const double *upper = solver_->getColUpper();
          const double *originalLower = continuousSolver_->getColLower();
          const double *originalUpper = continuousSolver_->getColUpper();
          double tolerance = 1.0e-5;
          for (int i = 0; i < numberColumns; i++) {
            if (lower[i] > originalLower[i] + tolerance) {
              nTightened++;
            }
            if (upper[i] < originalUpper[i] - tolerance) {
              nTightened++;
            }
          }
          if (nTightened) {
            fp = fopen("global.fix", "wb");
            size_t nWrite;
            indices[0] = numberColumns;
            if (originalColumns_)
              indices[1] = COIN_INT_MAX;
            else
              indices[1] = -1;
            nWrite = fwrite(indices, sizeof(int), 2, fp);
            if (nWrite != 2)
              throw("Error in fwrite");
            for (int i = 0; i < numberColumns; i++) {
              int nTightened = 0;
              if (lower[i] > originalLower[i] + tolerance) {
                nTightened++;
              }
              if (upper[i] < originalUpper[i] - tolerance) {
                nTightened++;
              }
              if (nTightened) {
                indices[0] = i;
                if (originalColumns_)
                  indices[1] = originalColumns_[i];
                elements[0] = lower[i];
                elements[1] = upper[i];
                elements[2] = originalLower[i];
                elements[3] = originalUpper[i];
                nWrite = fwrite(indices, sizeof(int), 2, fp);
                if (nWrite != 2)
                  throw("Error in fwrite");
                nWrite = fwrite(elements, sizeof(double), 4, fp);
                if (nWrite != 4)
                  throw("Error in fwrite");
              }
            }
            // eof marker
            indices[0] = -1;
            nWrite = fwrite(indices, sizeof(int), 2, fp);
            if (nWrite != 2)
              throw("Error in fwrite");
            fclose(fp);
          }
          delete[] elements;
          delete[] indices;
        }
        if ((specialOptions_ & 524288) != 0 && !parentModel_
          && storedRowCuts_) {
          if (feasible) {
            /* pick up stuff and try again
                        add cuts, maybe keep around
                        do best solution and if so new heuristics
                        obviously tighten bounds
                        */
            // A and B probably done on entry
            // A) tight bounds on integer variables
            const double *lower = solver_->getColLower();
            const double *upper = solver_->getColUpper();
            const double *tightLower = storedRowCuts_->tightLower();
            const double *tightUpper = storedRowCuts_->tightUpper();
            int nTightened = 0;
            for (int i = 0; i < numberIntegers_; i++) {
              int iColumn = integerVariable_[i];
              if (tightLower[iColumn] > lower[iColumn]) {
                nTightened++;
                solver_->setColLower(iColumn, tightLower[iColumn]);
              }
              if (tightUpper[iColumn] < upper[iColumn]) {
                nTightened++;
                solver_->setColUpper(iColumn, tightUpper[iColumn]);
              }
            }
            if (nTightened)
              COIN_DETAIL_PRINT(printf("%d tightened by alternate cuts\n", nTightened));
            if (storedRowCuts_->bestObjective() < bestObjective_) {
              // B) best solution
              double objValue = storedRowCuts_->bestObjective();
              setBestSolution(CBC_SOLUTION, objValue,
                storedRowCuts_->bestSolution());
              // Do heuristics
              // Allow RINS
              for (int i = 0; i < numberHeuristics_; i++) {
                CbcHeuristicRINS *rins
                  = dynamic_cast<CbcHeuristicRINS *>(heuristic_[i]);
                if (rins) {
                  rins->setLastNode(-100);
                }
              }
              doHeuristicsAtRoot();
            }
#ifdef JJF_ZERO
            int nCuts = storedRowCuts_->sizeRowCuts();
            // add to global list
            for (int i = 0; i < nCuts; i++) {
              OsiRowCut newCut(*storedRowCuts_->rowCutPointer(i));
              newCut.setGloballyValidAsInteger(2);
              newCut.mutableRow().setTestForDuplicateIndex(false);
              globalCuts_.insert(newCut);
            }
#else
          addCutGenerator(storedRowCuts_, -99, "Stored from previous run",
            true, false, false, -200);
#endif
            // Set cuts as active
            delete[] addedCuts_;
            maximumNumberCuts_ = cuts.sizeRowCuts();
            if (maximumNumberCuts_) {
              addedCuts_ = new CbcCountRowCut *[maximumNumberCuts_];
            } else {
              addedCuts_ = NULL;
            }
            for (int i = 0; i < maximumNumberCuts_; i++)
              addedCuts_[i] = new CbcCountRowCut(*cuts.rowCutPtr(i),
                NULL, -1, -1, 2);
            COIN_DETAIL_PRINT(printf("size %d\n", cuts.sizeRowCuts()));
            cuts = OsiCuts();
            currentNumberCuts_ = maximumNumberCuts_;
            feasible = solveWithCuts(cuts, maximumCutPassesAtRoot_,
              NULL);
            for (int i = 0; i < maximumNumberCuts_; i++)
              delete addedCuts_[i];
          }
          delete storedRowCuts_;
          storedRowCuts_ = NULL;
        }
      } else {
        feasible = false;
      }
    } else if (solverCharacteristics_->solutionAddsCuts() || solverCharacteristics_->alwaysTryCutsAtRootNode()) {
      // may generate cuts and turn the solution
      //to an infeasible one
      feasible = solveWithCuts(cuts, 2,
        NULL);
    }
  }
  if (rootModels) {
    int numberModels = multipleRootTries_ % 100;
    for (int i = 0; i < numberModels; i++)
      delete rootModels[i];
    delete[] rootModels;
  }
  // check extra info on feasibility
  if (!solverCharacteristics_->mipFeasible())
    feasible = false;
  // If max nodes==0 - don't do strong branching
  if (!getMaximumNodes()) {
    if (feasible)
      feasible = false;
    else
      setMaximumNodes(1); //allow to stop on success
  }
  topOfTree_ = NULL;
#ifdef CLP_RESOLVE
  if ((moreSpecialOptions_ & 2097152) != 0 && !parentModel_ && feasible) {
    OsiClpSolverInterface *clpSolver
      = dynamic_cast<OsiClpSolverInterface *>(solver_);
    if (clpSolver)
      resolveClp(clpSolver, 0);
  }
#endif
  // make cut generators less aggressive
  for (iCutGenerator = 0; iCutGenerator < numberCutGenerators_; iCutGenerator++) {
    CglCutGenerator *generator = generator_[iCutGenerator]->generator();
    generator->setAggressiveness(generator->getAggressiveness() - 100);
  }
  currentNumberCuts_ = numberNewCuts_;
  if (solverCharacteristics_->solutionAddsCuts()) {
    // With some heuristics solver needs a resolve here (don't know if this is bug in heuristics)
    solver_->resolve();
    if (!isProvenOptimal()) {
      solver_->initialSolve();
    }
  }
  // See if can stop on gap
  bestPossibleObjective_ = solver_->getObjValue() * solver_->getObjSense();
  if (canStopOnGap()) {
    if (bestPossibleObjective_ < getCutoff())
      stoppedOnGap_ = true;
    feasible = false;
  }
  // User event
  if (eventHappened_)
    feasible = false;
#if defined(COIN_HAS_CLP) && defined(COIN_HAS_CPX)
  /*
      This is the notion of using Cbc stuff to get going, then calling cplex to
      finish off.
    */
  if (feasible && (specialOptions_ & 16384) != 0 && fastNodeDepth_ == -2 && !parentModel_) {
    // Use Cplex to do search!
    double time1 = CoinCpuTime();
    OsiClpSolverInterface *clpSolver
      = dynamic_cast<OsiClpSolverInterface *>(solver_);
    OsiCpxSolverInterface cpxSolver;
    double direction = clpSolver->getObjSense();
    cpxSolver.setObjSense(direction);
    // load up cplex
    const CoinPackedMatrix *matrix = continuousSolver_->getMatrixByCol();
    const double *rowLower = continuousSolver_->getRowLower();
    const double *rowUpper = continuousSolver_->getRowUpper();
    const double *columnLower = continuousSolver_->getColLower();
    const double *columnUpper = continuousSolver_->getColUpper();
    const double *objective = continuousSolver_->getObjCoefficients();
    cpxSolver.loadProblem(*matrix, columnLower, columnUpper,
      objective, rowLower, rowUpper);
    double *setSol = new double[numberIntegers_];
    int *setVar = new int[numberIntegers_];
    // cplex doesn't know about objective offset
    double offset = clpSolver->getModelPtr()->objectiveOffset();
    for (int i = 0; i < numberIntegers_; i++) {
      int iColumn = integerVariable_[i];
      cpxSolver.setInteger(iColumn);
      if (bestSolution_) {
        setSol[i] = bestSolution_[iColumn];
        setVar[i] = iColumn;
      }
    }
    CPXENVptr env = cpxSolver.getEnvironmentPtr();
    CPXLPptr lpPtr = cpxSolver.getLpPtr(OsiCpxSolverInterface::KEEPCACHED_ALL);
    cpxSolver.switchToMIP();
    if (bestSolution_) {
#if 0
            CPXcopymipstart(env, lpPtr, numberIntegers_, setVar, setSol);
#else
        int zero = 0;
        CPXaddmipstarts(env, lpPtr, 1, numberIntegers_, &zero, setVar, setSol, NULL, NULL);
#endif
    }
    if (clpSolver->getNumRows() > continuousSolver_->getNumRows() && false) {
      // add cuts
      const CoinPackedMatrix *matrix = clpSolver->getMatrixByRow();
      const double *rhs = clpSolver->getRightHandSide();
      const char *rowSense = clpSolver->getRowSense();
      const double *elementByRow = matrix->getElements();
      const int *column = matrix->getIndices();
      const CoinBigIndex *rowStart = matrix->getVectorStarts();
      const int *rowLength = matrix->getVectorLengths();
      int nStart = continuousSolver_->getNumRows();
      int nRows = clpSolver->getNumRows();
      int size = rowStart[nRows - 1] + rowLength[nRows - 1] - rowStart[nStart];
      int nAdd = 0;
      double *rmatval = new double[size];
      int *rmatind = new int[size];
      int *rmatbeg = new int[nRows - nStart + 1];
      size = 0;
      rmatbeg[0] = 0;
      for (int i = nStart; i < nRows; i++) {
        for (int k = rowStart[i]; k < rowStart[i] + rowLength[i]; k++) {
          rmatind[size] = column[k];
          rmatval[size++] = elementByRow[k];
        }
        nAdd++;
        rmatbeg[nAdd] = size;
      }
      CPXaddlazyconstraints(env, lpPtr, nAdd, size,
        rhs, rowSense, rmatbeg,
        rmatind, rmatval, NULL);
      CPXsetintparam(env, CPX_PARAM_REDUCE,
        // CPX_PREREDUCE_NOPRIMALORDUAL (0)
        CPX_PREREDUCE_PRIMALONLY);
    }
    if (getCutoff() < 1.0e50) {
      double useCutoff = getCutoff() + offset;
      if (bestObjective_ < 1.0e50)
        useCutoff = bestObjective_ + offset + 1.0e-7;
      cpxSolver.setDblParam(OsiDualObjectiveLimit, useCutoff * direction);
      if (direction > 0.0)
        CPXsetdblparam(env, CPX_PARAM_CUTUP, useCutoff); // min
      else
        CPXsetdblparam(env, CPX_PARAM_CUTLO, useCutoff); // max
    }
    CPXsetdblparam(env, CPX_PARAM_EPGAP, dblParam_[CbcAllowableFractionGap]);
    delete[] setSol;
    delete[] setVar;
    char printBuffer[200];
    if (offset) {
      sprintf(printBuffer, "Add %g to all Cplex messages for true objective",
        -offset);
      messageHandler()->message(CBC_GENERAL, messages())
        << printBuffer << CoinMessageEol;
      cpxSolver.setDblParam(OsiObjOffset, offset);
    }
    cpxSolver.branchAndBound();
    double timeTaken = CoinCpuTime() - time1;
    sprintf(printBuffer, "Cplex took %g seconds",
      timeTaken);
    messageHandler()->message(CBC_GENERAL, messages())
      << printBuffer << CoinMessageEol;
    numberExtraNodes_ = CPXgetnodecnt(env, lpPtr);
    numberExtraIterations_ = CPXgetmipitcnt(env, lpPtr);
    double value = cpxSolver.getObjValue() * direction;
    if (cpxSolver.isProvenOptimal() && value <= getCutoff()) {
      feasible = true;
      clpSolver->setWarmStart(NULL);
      // try and do solution
      double *newSolution = CoinCopyOfArray(cpxSolver.getColSolution(),
        getNumCols());
      setBestSolution(CBC_STRONGSOL, value, newSolution);
      delete[] newSolution;
    }
    feasible = false;
  }
#endif
  if (!parentModel_ && (moreSpecialOptions_ & 268435456) != 0) {
    // try idiotic idea
    CbcObject *obj = new CbcIdiotBranch(this);
    obj->setPriority(1); // temp
    addObjects(1, &obj);
    delete obj;
  }

  /*
      A hook to use clp to quickly explore some part of the tree.
    */
  if (fastNodeDepth_ == 1000 && /*!parentModel_*/ (specialOptions_ & 2048) == 0) {
    fastNodeDepth_ = -1;
    CbcObject *obj = new CbcFollowOn(this);
    obj->setPriority(1);
    addObjects(1, &obj);
    delete obj;
  }
  int saveNumberSolves = numberSolves_;
  int saveNumberIterations = numberIterations_;
  if ((fastNodeDepth_ >= 0 || (moreSpecialOptions_ & 33554432) != 0)
    && /*!parentModel_*/ (specialOptions_ & 2048) == 0) {
    // add in a general depth object doClp
    int type = (fastNodeDepth_ <= 100) ? fastNodeDepth_ : -(fastNodeDepth_ - 100);
    if ((moreSpecialOptions_ & 33554432) != 0)
      type = 12;
    else
      fastNodeDepth_ += 1000000; // mark as done
    CbcObject *obj = new CbcGeneralDepth(this, type);
    addObjects(1, &obj);
    delete obj;
    // fake number of objects
    numberObjects_--;
    if (parallelMode() < -1) {
      // But make sure position is correct
      OsiObject *obj2 = object_[numberObjects_];
      obj = dynamic_cast<CbcObject *>(obj2);
      assert(obj);
      obj->setPosition(numberObjects_);
    }
  }
#ifdef COIN_HAS_CLP
#ifdef NO_CRUNCH
  if (true) {
    OsiClpSolverInterface *clpSolver
      = dynamic_cast<OsiClpSolverInterface *>(solver_);
    if (clpSolver && !parentModel_) {
      ClpSimplex *clpSimplex = clpSolver->getModelPtr();
      clpSimplex->setSpecialOptions(clpSimplex->specialOptions() | 131072);
      //clpSimplex->startPermanentArrays();
      clpSimplex->setPersistenceFlag(2);
    }
  }
#endif
#endif
  // Save copy of solver
  OsiSolverInterface *saveSolver = NULL;
  if (!parentModel_ && (specialOptions_ & (512 + 32768)) != 0)
    saveSolver = solver_->clone();
  double checkCutoffForRestart = 1.0e100;
  saveModel(saveSolver, &checkCutoffForRestart, &feasible);
  if ((specialOptions_ & 262144) != 0 && !parentModel_) {
    // Save stuff and return!
    storedRowCuts_->saveStuff(bestObjective_, bestSolution_,
      solver_->getColLower(),
      solver_->getColUpper());
    delete[] lowerBefore;
    delete[] upperBefore;
    delete saveSolver;
    if (flipObjective)
      flipModel();
    return;
  }
  /*
      We've taken the continuous relaxation as far as we can. Time to branch.
      The first order of business is to actually create a node. chooseBranch
      currently uses strong branching to evaluate branch object candidates,
      unless forced back to simple branching. If chooseBranch concludes that a
      branching candidate is monotone (anyAction == -1) or infeasible (anyAction
      == -2) when forced to integer values, it returns here immediately.

      Monotone variables trigger a call to resolve(). If the problem remains
      feasible, try again to choose a branching variable. At the end of the loop,
      resolved == true indicates that some variables were fixed.

      Loss of feasibility will result in the deletion of newNode.
    */

  bool resolved = false;
  CbcNode *newNode = NULL;
  numberFixedAtRoot_ = 0;
  numberFixedNow_ = 0;
  if (!parentModel_ && (moreSpecialOptions2_ & 2) != 0) {
#ifdef COIN_HAS_CLP
    OsiClpSolverInterface *clpSolver
      = dynamic_cast<OsiClpSolverInterface *>(solver_);
    if (clpSolver) {
      if (getCutoff() > 1.0e20) {
        printf("Zapping costs\n");
        int numberColumns = solver_->getNumCols();
        double *zeroCost = new double[numberColumns];
        // could make small random
        memset(zeroCost, 0, numberColumns * sizeof(double));
        solver_->setObjective(zeroCost);
        double objValue = solver_->getObjValue();
        solver_->setDblParam(OsiObjOffset, -objValue);
        clpSolver->getModelPtr()->setObjectiveValue(objValue);
        delete[] zeroCost;
      } else {
        moreSpecialOptions2_ &= ~2;
      }
    } else {
#endif
      moreSpecialOptions2_ &= ~2;
#ifdef COIN_HAS_CLP
    }
#endif
  }
  int numberIterationsAtContinuous = numberIterations_;
  //solverCharacteristics_->setSolver(solver_);
  if (feasible) {
    // mark all cuts as globally valid
    int numberCuts = cuts.sizeRowCuts();
    resizeWhichGenerator(0, numberCuts);
    for (int i = 0; i < numberCuts; i++) {
      cuts.rowCutPtr(i)->setGloballyValid();
      whichGenerator_[i] = 20000 + (whichGenerator_[i] % 10000);
    }
#define HOTSTART -1
#if HOTSTART < 0
    if (bestSolution_ && !parentModel_ && !hotstartSolution_ && (moreSpecialOptions_ & 1024) != 0 && (specialOptions_ & 2048) == 0) {
      // Set priorities so only branch on ones we need to
      // use djs and see if only few branches needed
#ifndef NDEBUG
      double integerTolerance = getIntegerTolerance();
#endif
      bool possible = true;
      const double *saveLower = continuousSolver_->getColLower();
      const double *saveUpper = continuousSolver_->getColUpper();
      for (int i = 0; i < numberObjects_; i++) {
        const CbcSimpleInteger *thisOne = dynamic_cast<const CbcSimpleInteger *>(object_[i]);
        if (thisOne) {
          int iColumn = thisOne->columnNumber();
          if (saveUpper[iColumn] > saveLower[iColumn] + 1.5) {
            possible = false;
            break;
          }
        } else {
          possible = false;
          break;
        }
      }
      if (possible) {
        OsiSolverInterface *solver = continuousSolver_->clone();
        int numberColumns = solver->getNumCols();
        for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
          double value = bestSolution_[iColumn];
          value = CoinMax(value, saveLower[iColumn]);
          value = CoinMin(value, saveUpper[iColumn]);
          value = floor(value + 0.5);
          if (solver->isInteger(iColumn)) {
            solver->setColLower(iColumn, value);
            solver->setColUpper(iColumn, value);
          }
        }
        solver->setHintParam(OsiDoDualInResolve, false, OsiHintTry);
        // objlim and all slack
        double direction = solver->getObjSense();
        solver->setDblParam(OsiDualObjectiveLimit, 1.0e50 * direction);
        CoinWarmStartBasis *basis = dynamic_cast<CoinWarmStartBasis *>(solver->getEmptyWarmStart());
        solver->setWarmStart(basis);
        delete basis;
        bool changed = true;
        hotstartPriorities_ = new int[numberColumns];
        for (int iColumn = 0; iColumn < numberColumns; iColumn++)
          hotstartPriorities_[iColumn] = 1;
        while (changed) {
          changed = false;
          solver->resolve();
          if (!solver->isProvenOptimal()) {
            possible = false;
            break;
          }
          const double *dj = solver->getReducedCost();
          const double *colLower = solver->getColLower();
          const double *colUpper = solver->getColUpper();
          const double *solution = solver->getColSolution();
          int nAtLbNatural = 0;
          int nAtUbNatural = 0;
          int nZeroDj = 0;
          int nForced = 0;
          for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
            double value = solution[iColumn];
            value = CoinMax(value, saveLower[iColumn]);
            value = CoinMin(value, saveUpper[iColumn]);
            if (solver->isInteger(iColumn)) {
              assert(fabs(value - solution[iColumn]) <= integerTolerance);
              if (hotstartPriorities_[iColumn] == 1) {
                if (dj[iColumn] < -1.0e-6) {
                  // negative dj
                  if (saveUpper[iColumn] == colUpper[iColumn]) {
                    nAtUbNatural++;
                    hotstartPriorities_[iColumn] = 2;
                    solver->setColLower(iColumn, saveLower[iColumn]);
                    solver->setColUpper(iColumn, saveUpper[iColumn]);
                  } else {
                    nForced++;
                  }
                } else if (dj[iColumn] > 1.0e-6) {
                  // positive dj
                  if (saveLower[iColumn] == colLower[iColumn]) {
                    nAtLbNatural++;
                    hotstartPriorities_[iColumn] = 2;
                    solver->setColLower(iColumn, saveLower[iColumn]);
                    solver->setColUpper(iColumn, saveUpper[iColumn]);
                  } else {
                    nForced++;
                  }
                } else {
                  // zero dj
                  nZeroDj++;
                }
              }
            }
          }
#if CBC_USEFUL_PRINTING > 1
          printf("%d forced, %d naturally at lower, %d at upper - %d zero dj\n",
            nForced, nAtLbNatural, nAtUbNatural, nZeroDj);
#endif
          if (nAtLbNatural || nAtUbNatural) {
            changed = true;
          } else {
            if (nForced + nZeroDj > 5000 || (nForced + nZeroDj) * 2 > numberIntegers_)
              possible = false;
          }
        }
        delete solver;
      }
      if (possible) {
        setHotstartSolution(bestSolution_);
        if (!saveCompare) {
          // create depth first comparison
          saveCompare = nodeCompare_;
          // depth first
          nodeCompare_ = new CbcCompareDepth();
          tree_->setComparison(*nodeCompare_);
        }
      } else {
        delete[] hotstartPriorities_;
        hotstartPriorities_ = NULL;
      }
    }
#endif
#if HOTSTART > 0
    if (hotstartSolution_ && !hotstartPriorities_) {
      // Set up hot start
      OsiSolverInterface *solver = solver_->clone();
      double direction = solver_->getObjSense();
      int numberColumns = solver->getNumCols();
      double *saveLower = CoinCopyOfArray(solver->getColLower(), numberColumns);
      double *saveUpper = CoinCopyOfArray(solver->getColUpper(), numberColumns);
      // move solution
      solver->setColSolution(hotstartSolution_);
      // point to useful information
      const double *saveSolution = testSolution_;
      testSolution_ = solver->getColSolution();
      OsiBranchingInformation usefulInfo = usefulInformation();
      testSolution_ = saveSolution;
      /*
            Run through the objects and use feasibleRegion() to set variable bounds
            so as to fix the variables specified in the objects at their value in this
            solution. Since the object list contains (at least) one object for every
            integer variable, this has the effect of fixing all integer variables.
            */
      for (int i = 0; i < numberObjects_; i++)
        object_[i]->feasibleRegion(solver, &usefulInfo);
      solver->resolve();
      assert(solver->isProvenOptimal());
      double gap = CoinMax((solver->getObjValue() - solver_->getObjValue()) * direction, 0.0);
      const double *dj = solver->getReducedCost();
      const double *colLower = solver->getColLower();
      const double *colUpper = solver->getColUpper();
      const double *solution = solver->getColSolution();
      int nAtLbNatural = 0;
      int nAtUbNatural = 0;
      int nAtLbNaturalZero = 0;
      int nAtUbNaturalZero = 0;
      int nAtLbFixed = 0;
      int nAtUbFixed = 0;
      int nAtOther = 0;
      int nAtOtherNatural = 0;
      int nNotNeeded = 0;
      delete[] hotstartSolution_;
      hotstartSolution_ = new double[numberColumns];
      delete[] hotstartPriorities_;
      hotstartPriorities_ = new int[numberColumns];
      int *order = (int *)saveUpper;
      int nFix = 0;
      double bestRatio = COIN_DBL_MAX;
      for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
        double value = solution[iColumn];
        value = CoinMax(value, saveLower[iColumn]);
        value = CoinMin(value, saveUpper[iColumn]);
        double sortValue = COIN_DBL_MAX;
        if (solver->isInteger(iColumn)) {
          assert(fabs(value - solution[iColumn]) <= 1.0e-5);
          double value2 = floor(value + 0.5);
          if (dj[iColumn] < -1.0e-6) {
            // negative dj
            //assert (value2==colUpper[iColumn]);
            if (saveUpper[iColumn] == colUpper[iColumn]) {
              nAtUbNatural++;
              sortValue = 0.0;
              double value = -dj[iColumn];
              if (value > gap)
                nFix++;
              else if (gap < value * bestRatio)
                bestRatio = gap / value;
              if (saveLower[iColumn] != colLower[iColumn]) {
                nNotNeeded++;
                sortValue = 1.0e20;
              }
            } else if (saveLower[iColumn] == colUpper[iColumn]) {
              nAtLbFixed++;
              sortValue = dj[iColumn];
            } else {
              nAtOther++;
              sortValue = 0.0;
              if (saveLower[iColumn] != colLower[iColumn] && saveUpper[iColumn] != colUpper[iColumn]) {
                nNotNeeded++;
                sortValue = 1.0e20;
              }
            }
          } else if (dj[iColumn] > 1.0e-6) {
            // positive dj
            //assert (value2==colLower[iColumn]);
            if (saveLower[iColumn] == colLower[iColumn]) {
              nAtLbNatural++;
              sortValue = 0.0;
              double value = dj[iColumn];
              if (value > gap)
                nFix++;
              else if (gap < value * bestRatio)
                bestRatio = gap / value;
              if (saveUpper[iColumn] != colUpper[iColumn]) {
                nNotNeeded++;
                sortValue = 1.0e20;
              }
            } else if (saveUpper[iColumn] == colLower[iColumn]) {
              nAtUbFixed++;
              sortValue = -dj[iColumn];
            } else {
              nAtOther++;
              sortValue = 0.0;
              if (saveLower[iColumn] != colLower[iColumn] && saveUpper[iColumn] != colUpper[iColumn]) {
                nNotNeeded++;
                sortValue = 1.0e20;
              }
            }
          } else {
            // zero dj
            if (value2 == saveUpper[iColumn]) {
              nAtUbNaturalZero++;
              sortValue = 0.0;
              if (saveLower[iColumn] != colLower[iColumn]) {
                nNotNeeded++;
                sortValue = 1.0e20;
              }
            } else if (value2 == saveLower[iColumn]) {
              nAtLbNaturalZero++;
              sortValue = 0.0;
            } else {
              nAtOtherNatural++;
              sortValue = 0.0;
              if (saveLower[iColumn] != colLower[iColumn] && saveUpper[iColumn] != colUpper[iColumn]) {
                nNotNeeded++;
                sortValue = 1.0e20;
              }
            }
          }
#if HOTSTART == 3
          sortValue = -fabs(dj[iColumn]);
#endif
        }
        hotstartSolution_[iColumn] = value;
        saveLower[iColumn] = sortValue;
        order[iColumn] = iColumn;
      }
      COIN_DETAIL_PRINT(printf("** can fix %d columns - best ratio for others is %g on gap of %g\n",
        nFix, bestRatio, gap));
      int nNeg = 0;
      CoinSort_2(saveLower, saveLower + numberColumns, order);
      for (int i = 0; i < numberColumns; i++) {
        if (saveLower[i] < 0.0) {
          nNeg++;
#if HOTSTART == 2 || HOTSTART == 3
          // swap sign ?
          saveLower[i] = -saveLower[i];
#endif
        }
      }
      CoinSort_2(saveLower, saveLower + nNeg, order);
      for (int i = 0; i < numberColumns; i++) {
#if HOTSTART == 1
        hotstartPriorities_[order[i]] = 100;
#else
          hotstartPriorities_[order[i]] = -(i + 1);
#endif
      }
      COIN_DETAIL_PRINT(printf("nAtLbNat %d,nAtUbNat %d,nAtLbNatZero %d,nAtUbNatZero %d,nAtLbFixed %d,nAtUbFixed %d,nAtOther %d,nAtOtherNat %d, useless %d %d\n",
        nAtLbNatural,
        nAtUbNatural,
        nAtLbNaturalZero,
        nAtUbNaturalZero,
        nAtLbFixed,
        nAtUbFixed,
        nAtOther,
        nAtOtherNatural, nNotNeeded, nNeg));
      delete[] saveLower;
      delete[] saveUpper;
      if (!saveCompare) {
        // create depth first comparison
        saveCompare = nodeCompare_;
        // depth first
        nodeCompare_ = new CbcCompareDepth();
        tree_->setComparison(*nodeCompare_);
      }
    }
#endif
    newNode = new CbcNode;
    // Set objective value (not so obvious if NLP etc)
    setObjectiveValue(newNode, NULL);
    anyAction = -1;
    // To make depth available we may need a fake node
    CbcNode fakeNode;
    if (!currentNode_) {
      // Not true if sub trees assert (!numberNodes_);
      currentNode_ = &fakeNode;
    }
    phase_ = 3;
    // only allow 1000 passes
    int numberPassesLeft = 1000;
    // This is first crude step
    if (numberAnalyzeIterations_ && !parentModel_) {
      delete[] analyzeResults_;
      //int numberColumns = solver_->getNumCols();
      analyzeResults_ = new double[5 * numberIntegers_];
      numberFixedAtRoot_ = newNode->analyze(this, analyzeResults_);
      if (numberFixedAtRoot_ > 0) {
        COIN_DETAIL_PRINT(printf("%d fixed by analysis\n", numberFixedAtRoot_));
        setPointers(solver_);
        numberFixedNow_ = numberFixedAtRoot_;
      } else if (numberFixedAtRoot_ < 0) {
        COIN_DETAIL_PRINT(printf("analysis found to be infeasible\n"));
        anyAction = -2;
        delete newNode;
        newNode = NULL;
        feasible = false;
      }
    }
    OsiSolverBranch *branches = NULL;
    if (feasible)
      anyAction = chooseBranch(newNode, numberPassesLeft, NULL, cuts, resolved,
        NULL, NULL, NULL, branches);
    if (anyAction == -2 || newNode->objectiveValue() >= cutoff) {
      if (anyAction != -2) {
        // zap parent nodeInfo
#ifdef COIN_DEVELOP
        printf("zapping CbcNodeInfo %x\n", newNode->nodeInfo()->parent());
#endif
        if (newNode->nodeInfo())
          newNode->nodeInfo()->nullParent();
      }
      delete newNode;
      newNode = NULL;
      feasible = false;
    }
  }
  if (newNode && probingInfo_) {
    int number01 = probingInfo_->numberIntegers();
    //const fixEntry * entry = probingInfo_->fixEntries();
    const int *toZero = probingInfo_->toZero();
    //const int * toOne = probingInfo_->toOne();
    //const int * integerVariable = probingInfo_->integerVariable();
    if (toZero[number01]) {
      CglTreeProbingInfo info(*probingInfo_);
      if ((moreSpecialOptions2_ & 64) != 0 && !parentModel_) {
        /*
                Marginal idea. Further exploration probably good. Build some extra
                cliques from probing info. Not quite worth the effort?
              */
        CglProbing generator1;
        generator1.setUsingObjective(false);
        generator1.setMaxPass(1);
        generator1.setMaxPassRoot(1);
        generator1.setMaxLook(100);
        generator1.setRowCuts(3);
        generator1.setMaxElements(300);
        generator1.setMaxProbeRoot(solver_->getNumCols());
        CoinThreadRandom randomGenerator;
        //CglTreeProbingInfo info(solver_);
        info.level = 0;
        info.formulation_rows = solver_->getNumRows();
        info.inTree = false;
        if (parentModel_) {
          info.parentSolver = parentModel_->continuousSolver();
          // indicate if doing full search
          info.hasParent = ((specialOptions_ & 67108864) == 0) ? 1 : 2;
        } else {
          info.hasParent = 0;
          info.parentSolver = NULL;
        }
        info.originalColumns = originalColumns();
        info.randomNumberGenerator = &randomGenerator;
        info.pass = 4;
        generator1.setMode(8);
        OsiCuts cs;
        generator1.generateCutsAndModify(*solver_, cs, &info);
        // very clunky
        OsiSolverInterface *temp = generator1.cliqueModel(solver_, 2);
        CglPreProcess dummy;
        OsiSolverInterface *newSolver = dummy.cliqueIt(*temp, 0.0001);
        delete temp;
        OsiSolverInterface *fake = NULL;
        if (newSolver) {
#if 0
                int numberCliques = generator1.numberCliques();
                cliqueEntry * entry = generator1.cliqueEntry();
                cliqueType * type = new cliqueType [numberCliques];
                int * start = new int [numberCliques+1];
                start[numberCliques]=2*numberCliques;
                int n=0;
                for (int i=0;i<numberCliques;i++) {
                  start[i]=2*i;
                  setOneFixesInCliqueEntry(entry[2*i],true);
                  setOneFixesInCliqueEntry(entry[2*i+1],true);
                  type[i]=0;
                }
                fake = info.analyze(*solver_, 1,numberCliques,start,
                                    entry,type);
                delete [] type;
                delete [] entry;
#else
        fake = info.analyze(*newSolver, 1, -1);
#endif
          delete newSolver;
        } else {
          fake = info.analyze(*solver_, 1);
        }
        if (fake) {
          //fake->writeMps("fake");
          CglFakeClique cliqueGen(fake);
          cliqueGen.setStarCliqueReport(false);
          cliqueGen.setRowCliqueReport(false);
          cliqueGen.setMinViolation(0.1);
          addCutGenerator(&cliqueGen, 1, "Fake cliques", true, false, false, -200);
          generator_[numberCutGenerators_ - 1]->setTiming(true);
          for (int i = 0; i < numberCutGenerators_; i++) {
            CglKnapsackCover *cutGen = dynamic_cast<CglKnapsackCover *>(generator_[i]->generator());
            if (cutGen) {
              cutGen->createCliques(*fake, 2, 200, false);
            }
          }
        }
      }
      if (probingInfo_->packDown()) {
#if CBC_USEFUL_PRINTING > 1
        printf("%d implications on %d 0-1\n", toZero[number01], number01);
#endif
        // Create a cut generator that remembers implications discovered at root.
        CglImplication implication(probingInfo_);
        addCutGenerator(&implication, 1, "ImplicationCuts", true, false, false, -200);
        generator_[numberCutGenerators_ - 1]->setGlobalCuts(true);
        generator_[numberCutGenerators_ - 1]->setTiming(true);
      } else {
        delete probingInfo_;
        probingInfo_ = NULL;
      }
    } else {
      delete probingInfo_;

      probingInfo_ = NULL;
    }
  }
  /*
      At this point, the root subproblem is infeasible or fathomed by bound
      (newNode == NULL), or we're live with an objective value that satisfies the
      current objective cutoff.
    */
  assert(!newNode || newNode->objectiveValue() <= cutoff);
  // Save address of root node as we don't want to delete it
  /*
      The common case is that the lp relaxation is feasible but doesn't satisfy
      integrality (i.e., newNode->branchingObject(), indicating we've been able to
      select a branching variable). Remove any cuts that have gone slack due to
      forcing monotone variables. Then tack on an CbcFullNodeInfo object and full
      basis (via createInfo()) and stash the new cuts in the nodeInfo (via
      addCuts()). If, by some miracle, we have an integral solution at the root
      (newNode->branchingObject() is NULL), takeOffCuts() will ensure that the solver holds
      a valid solution for use by setBestSolution().
    */
  CoinWarmStartBasis *lastws = NULL;
  if (feasible && newNode->branchingObject()) {
    if (resolved) {
      takeOffCuts(cuts, false, NULL);
#ifdef CHECK_CUT_COUNTS
      {
        printf("Number of rows after chooseBranch fix (root)"
               "(active only) %d\n",
          numberRowsAtContinuous_ + numberNewCuts_ + numberOldActiveCuts_);
        const CoinWarmStartBasis *debugws = dynamic_cast<const CoinWarmStartBasis *>(solver_->getWarmStart());
        debugws->print();
        delete debugws;
      }
#endif
    }
    //newNode->createInfo(this,NULL,NULL,NULL,NULL,0,0) ;
    //newNode->nodeInfo()->addCuts(cuts,
    //                   newNode->numberBranches(),whichGenerator_) ;
    if (lastws)
      delete lastws;
    lastws = dynamic_cast<CoinWarmStartBasis *>(solver_->getWarmStart());
  }
  /*
      Continuous data to be used later
    */
  continuousObjective_ = solver_->getObjValue() * solver_->getObjSense();
  continuousInfeasibilities_ = 0;
  if (newNode) {
    continuousObjective_ = newNode->objectiveValue();
    delete[] continuousSolution_;
    continuousSolution_ = CoinCopyOfArray(solver_->getColSolution(),
      numberColumns);
    continuousInfeasibilities_ = newNode->numberUnsatisfied();
  }
  /*
      Bound may have changed so reset in objects
    */
  {
    int i;
    for (i = 0; i < numberObjects_; i++)
      object_[i]->resetBounds(solver_);
  }
  /*
      Feasible? Then we should have either a live node prepped for future
      expansion (indicated by variable() >= 0), or (miracle of miracles) an
      integral solution at the root node.

      initializeInfo sets the reference counts in the nodeInfo object.  Since
      this node is still live, push it onto the heap that holds the live set.
    */
  if (newNode) {
    if (newNode->branchingObject()) {
      newNode->initializeInfo();
      tree_->push(newNode);
      // save pointer to root node - so can pick up bounds
      if (!topOfTree_)
        topOfTree_ = dynamic_cast<CbcFullNodeInfo *>(newNode->nodeInfo());
      if (statistics_) {
        if (numberNodes2_ == maximumStatistics_) {
          maximumStatistics_ = 2 * maximumStatistics_;
          CbcStatistics **temp = new CbcStatistics *[maximumStatistics_];
          memset(temp, 0, maximumStatistics_ * sizeof(CbcStatistics *));
          memcpy(temp, statistics_, numberNodes2_ * sizeof(CbcStatistics *));
          delete[] statistics_;
          statistics_ = temp;
        }
        assert(!statistics_[numberNodes2_]);
        statistics_[numberNodes2_] = new CbcStatistics(newNode, this);
      }
      numberNodes2_++;
#ifdef CHECK_NODE
      printf("Node %x on tree\n", newNode);
#endif
    } else {
      // continuous is integer
      double objectiveValue = newNode->objectiveValue();
      setBestSolution(CBC_SOLUTION, objectiveValue,
        solver_->getColSolution());
      if (eventHandler) {
        // we are stopping anyway so no need to test return code
        eventHandler->event(CbcEventHandler::solution);
      }
      delete newNode;
      newNode = NULL;
    }
  }

  if (printFrequency_ <= 0) {
    printFrequency_ = 1000;
    if (getNumCols() > 2000)
      printFrequency_ = 100;
  }
  /*
      It is possible that strong branching fixes one variable and then the code goes round
      again and again.  This can take too long.  So we need to warn user - just once.
    */
  numberLongStrong_ = 0;
  CbcNode *createdNode = NULL;
#ifdef CBC_THREAD
  if ((specialOptions_ & 2048) != 0)
    numberThreads_ = 0;
  if (numberThreads_) {
    nodeCompare_->sayThreaded(); // need to use addresses
    master_ = new CbcBaseModel(*this,
      (parallelMode() < -1) ? 1 : 0);
    masterThread_ = master_->masterThread();
  }
#endif
#ifdef COIN_HAS_CLP
  {
    OsiClpSolverInterface *clpSolver
      = dynamic_cast<OsiClpSolverInterface *>(solver_);
    if (clpSolver && !parentModel_) {
      clpSolver->computeLargestAway();
    }
  }
#endif
  /*
      At last, the actual branch-and-cut search loop, which will iterate until
      the live set is empty or we hit some limit (integrality gap, time, node
      count, etc.). The overall flow is to rebuild a subproblem, reoptimise using
      solveWithCuts(), choose a branching pattern with chooseBranch(), and finally
      add the node to the live set.

      The first action is to winnow the live set to remove nodes which are worse
      than the current objective cutoff.
    */
  if (solver_->getRowCutDebuggerAlways()) {
    OsiRowCutDebugger *debuggerX = const_cast<OsiRowCutDebugger *>(solver_->getRowCutDebuggerAlways());
    const OsiRowCutDebugger *debugger = solver_->getRowCutDebugger();
    if (!debugger) {
      // infeasible!!
      printf("before search\n");
      const double *lower = solver_->getColLower();
      const double *upper = solver_->getColUpper();
      const double *solution = debuggerX->optimalSolution();
      int numberColumns = solver_->getNumCols();
      for (int i = 0; i < numberColumns; i++) {
        if (solver_->isInteger(i)) {
          if (solution[i] < lower[i] - 1.0e-6 || solution[i] > upper[i] + 1.0e-6)
            printf("**** ");
          printf("%d %g <= %g <= %g\n",
            i, lower[i], solution[i], upper[i]);
        }
      }
      //abort();
    }
  }
  {
    // may be able to change cutoff now
    double cutoff = getCutoff();
    double increment = getDblParam(CbcModel::CbcCutoffIncrement);
    if (cutoff > bestObjective_ - increment) {
      cutoff = bestObjective_ - increment;
      setCutoff(cutoff);
    }
  }
#ifdef CBC_THREAD
  bool goneParallel = false;
#endif
#define MAX_DEL_NODE 1
  CbcNode *delNode[MAX_DEL_NODE + 1];
  int nDeleteNode = 0;
  // For Printing etc when parallel
  int lastEvery1000 = 0;
  int lastPrintEvery = 0;
  int numberConsecutiveInfeasible = 0;
#define PERTURB_IN_FATHOM
#ifdef PERTURB_IN_FATHOM
  // allow in fathom
  if ((moreSpecialOptions_ & 262144) != 0)
    specialOptions_ |= 131072;
#endif
  while (true) {
    lockThread();
#ifdef COIN_HAS_CLP
    // See if we want dantzig row choice
    goToDantzig(100, savePivotMethod);
#endif
    //#define REPORT_DYNAMIC 2
#if REPORT_DYNAMIC
    if (numberNodes_ && !parentModel_ && (tree_->empty() || (numberNodes_ % 10000) == 0)) {
      // Get average up and down costs
      double averageUp = 0.0;
      double averageDown = 0.0;
      int numberUp = 0;
      int numberDown = 0;
      int minTimesDown = COIN_INT_MAX;
      int maxTimesDown = 0;
      int neverBranchedDown = 0;
      int infeasibleTimesDown = 0;
      int minTimesUp = COIN_INT_MAX;
      int maxTimesUp = 0;
      int infeasibleTimesUp = 0;
      int neverBranchedUp = 0;
      int neverBranched = 0;
      int i;
      int numberInts = 0;
      bool endOfSearch = tree_->empty();
      int numberUp2 = 0;
      int numberDown2 = 0;
      for (i = 0; i < numberObjects_; i++) {
        OsiObject *object = object_[i];
        CbcSimpleIntegerDynamicPseudoCost *dynamicObject = dynamic_cast<CbcSimpleIntegerDynamicPseudoCost *>(object);
        if (dynamicObject) {
          numberInts++;
          if (dynamicObject->numberTimesUp() || dynamicObject->numberTimesDown()) {
            int nUp = 0;
            int nDown = 0;
            double up = 0.0;
            double down = 0.0;
            if (dynamicObject->numberTimesUp()) {
              numberUp++;
              nUp = dynamicObject->numberTimesUp();
              minTimesUp = CoinMin(minTimesUp, nUp);
              maxTimesUp = CoinMax(maxTimesUp, nUp);
              up = dynamicObject->upDynamicPseudoCost();
              averageUp += up;
              numberUp2 += nUp;
              infeasibleTimesUp += dynamicObject->numberTimesUpInfeasible();
            } else {
              neverBranchedUp++;
            }
            if (dynamicObject->numberTimesDown()) {
              numberDown++;
              nDown = dynamicObject->numberTimesDown();
              minTimesDown = CoinMin(minTimesDown, nDown);
              maxTimesDown = CoinMax(maxTimesDown, nDown);
              down = dynamicObject->downDynamicPseudoCost();
              averageDown += down;
              numberDown2 += dynamicObject->numberTimesDown();
              infeasibleTimesDown += dynamicObject->numberTimesDownInfeasible();
            } else {
              neverBranchedDown++;
            }
#if REPORT_DYNAMIC > 1
#if REPORT_DYNAMIC == 2
            if (endOfSearch && numberIntegers_ < 400) {
#elif REPORT_DYNAMIC == 3
            if (endOfSearch) {
#else
            {
#endif
              dynamicObject->print(0, 0.0);
            }
#endif
          } else {
            neverBranched++;
#if REPORT_DYNAMIC > 2
#if REPORT_DYNAMIC == 3
            if (endOfSearch && numberIntegers_ < 400) {
#elif REPORT_DYNAMIC == 4
            if (endOfSearch) {
#else
            {
#endif
              printf("col %d - never branched on\n", dynamicObject->columnNumber());
            }
#endif
          }
        }
      }
      if (numberUp)
        averageUp /= static_cast<double>(numberUp);
      else
        averageUp = 0.0;
      if (numberDown)
        averageDown /= static_cast<double>(numberDown);
      else
        averageDown = 0.0;
      printf("Report for %d variables (%d never branched on) after %d nodes - total solves down %d up %d\n",
        numberInts, neverBranched, numberNodes_, numberDown2, numberUp2);
      if ((neverBranchedDown || neverBranchedUp) && endOfSearch)
        printf("odd %d never branched down and %d never branched up\n",
          neverBranchedDown, neverBranchedUp);
      printf("down average %g times (%d infeasible) average increase %g min/max times (%d,%d)\n",
        static_cast<double>(numberDown2) / numberDown, infeasibleTimesDown, averageDown,
        minTimesDown, maxTimesDown);
      printf("up average %g times (%d infeasible) average increase %g min/max times (%d,%d)\n",
        static_cast<double>(numberUp2) / numberUp, infeasibleTimesUp, averageUp,
        minTimesUp, maxTimesUp);
    }
#endif
    if (tree_->empty()) {
#ifdef CBC_THREAD
      if (parallelMode() > 0 && master_) {
        int anyLeft = master_->waitForThreadsInTree(0);
        if (!anyLeft) {
          master_->stopThreads(-1);
          break;
        }
      } else {
        break;
      }
#else
    break;
#endif
    } else {
      unlockThread();
    }
    // If done 50/100 nodes see if worth trying reduction
    if (numberNodes_ >= nextCheckRestart) {
      if (nextCheckRestart < 100)
        nextCheckRestart = 100;
      else
        nextCheckRestart = COIN_INT_MAX;
#ifdef COIN_HAS_CLP
      OsiClpSolverInterface *clpSolver
        = dynamic_cast<OsiClpSolverInterface *>(solver_);
      if (clpSolver && ((specialOptions_ & 131072) == 0) && true) {
        ClpSimplex *simplex = clpSolver->getModelPtr();
        int perturbation = simplex->perturbation();
#if CBC_USEFUL_PRINTING > 1
        printf("Testing its n,s %d %d solves n,s %d %d - pert %d\n",
          numberIterations_, saveNumberIterations,
          numberSolves_, saveNumberSolves, perturbation);
#endif
        if (perturbation == 50 && (numberIterations_ - saveNumberIterations) < 8 * (numberSolves_ - saveNumberSolves)) {
          // switch off perturbation
          simplex->setPerturbation(100);
#if CBC_USEFUL_PRINTING > 1
          printf("Perturbation switched off\n");
#endif
        }
      }
#endif
      /*
              Decide if we want to do a restart.
            */
      if (saveSolver && (specialOptions_ & (512 + 32768)) != 0) {
        bool tryNewSearch = solverCharacteristics_->reducedCostsAccurate() && (getCutoff() < 1.0e20 && getCutoff() < checkCutoffForRestart);
        int numberColumns = getNumCols();
        if (tryNewSearch) {
          // adding increment back allows current best - tiny bit weaker
          checkCutoffForRestart = getCutoff() + getCutoffIncrement();
#if CBC_USEFUL_PRINTING > 1
          printf("after %d nodes, cutoff %g - looking\n",
            numberNodes_, getCutoff());
#endif
          saveSolver->resolve();
          double direction = saveSolver->getObjSense();
          double gap = checkCutoffForRestart - saveSolver->getObjValue() * direction;
          double tolerance;
          saveSolver->getDblParam(OsiDualTolerance, tolerance);
          if (gap <= 0.0)
            gap = tolerance;
          gap += 100.0 * tolerance;
          double integerTolerance = getDblParam(CbcIntegerTolerance);

          const double *lower = saveSolver->getColLower();
          const double *upper = saveSolver->getColUpper();
          const double *solution = saveSolver->getColSolution();
          const double *reducedCost = saveSolver->getReducedCost();

          int numberFixed = 0;
          int numberFixed2 = 0;
#ifdef COIN_DEVELOP
          printf("gap %g\n", gap);
#endif
          for (int i = 0; i < numberIntegers_; i++) {
            int iColumn = integerVariable_[i];
            double djValue = direction * reducedCost[iColumn];
            if (upper[iColumn] - lower[iColumn] > integerTolerance) {
              if (solution[iColumn] < lower[iColumn] + integerTolerance && djValue > gap) {
                //printf("%d to lb on dj of %g - bounds %g %g\n",
                //     iColumn,djValue,lower[iColumn],upper[iColumn]);
                saveSolver->setColUpper(iColumn, lower[iColumn]);
                numberFixed++;
              } else if (solution[iColumn] > upper[iColumn] - integerTolerance && -djValue > gap) {
                //printf("%d to ub on dj of %g - bounds %g %g\n",
                //     iColumn,djValue,lower[iColumn],upper[iColumn]);
                saveSolver->setColLower(iColumn, upper[iColumn]);
                numberFixed++;
              }
            } else {
              //printf("%d has dj of %g - already fixed to %g\n",
              //     iColumn,djValue,lower[iColumn]);
              numberFixed2++;
            }
          }
#ifdef COIN_DEVELOP
          if ((specialOptions_ & 1) != 0) {
            const OsiRowCutDebugger *debugger = saveSolver->getRowCutDebugger();
            if (debugger) {
              printf("Contains optimal\n");
              OsiSolverInterface *temp = saveSolver->clone();
              const double *solution = debugger->optimalSolution();
              const double *lower = temp->getColLower();
              const double *upper = temp->getColUpper();
              int n = temp->getNumCols();
              for (int i = 0; i < n; i++) {
                if (temp->isInteger(i)) {
                  double value = floor(solution[i] + 0.5);
                  assert(value >= lower[i] && value <= upper[i]);
                  temp->setColLower(i, value);
                  temp->setColUpper(i, value);
                }
              }
              temp->writeMps("reduced_fix");
              delete temp;
              saveSolver->writeMps("reduced");
            } else {
              abort();
            }
          }
          printf("Restart could fix %d integers (%d already fixed)\n",
            numberFixed + numberFixed2, numberFixed2);
#endif
          numberFixed += numberFixed2;
          if (numberFixed * 10 < numberColumns && numberFixed * 4 < numberIntegers_)
            tryNewSearch = false;
        }
#ifdef CONFLICT_CUTS
        // temporary
        //if ((moreSpecialOptions_&4194304)!=0)
        //tryNewSearch=false;
#endif
        if (tryNewSearch) {
          // back to solver without cuts?
          OsiSolverInterface *solver2 = saveSolver->clone();
          const double *lower = saveSolver->getColLower();
          const double *upper = saveSolver->getColUpper();
          for (int i = 0; i < numberIntegers_; i++) {
            int iColumn = integerVariable_[i];
            solver2->setColLower(iColumn, lower[iColumn]);
            solver2->setColUpper(iColumn, upper[iColumn]);
          }
          // swap
          delete saveSolver;
          saveSolver = solver2;
          double *newSolution = new double[numberColumns];
          double objectiveValue = checkCutoffForRestart;
          // Save the best solution so far.
          CbcSerendipity heuristic(*this);
          if (bestSolution_)
            heuristic.setInputSolution(bestSolution_, bestObjective_);
          // Magic number
          heuristic.setFractionSmall(0.8);
          // `pumpTune' to stand-alone solver for explanations.
          heuristic.setFeasibilityPumpOptions(1008013);
          // Use numberNodes to say how many are original rows
          heuristic.setNumberNodes(continuousSolver_->getNumRows());
#ifdef COIN_DEVELOP
          if (continuousSolver_->getNumRows() < solver_->getNumRows())
            printf("%d rows added ZZZZZ\n",
              solver_->getNumRows() - continuousSolver_->getNumRows());
#endif
          int returnCode = heuristic.smallBranchAndBound(saveSolver,
            -1, newSolution,
            objectiveValue,
            checkCutoffForRestart, "Reduce");
          if (returnCode < 0) {
#ifdef COIN_DEVELOP
            printf("Restart - not small enough to do search after fixing\n");
#endif
            delete[] newSolution;
          } else {
            // 1 for sol'n, 2 for finished, 3 for both
            if ((returnCode & 1) != 0) {
              // increment number of solutions so other heuristics can test
              numberSolutions_++;
              numberHeuristicSolutions_++;
              lastHeuristic_ = NULL;
              setBestSolution(CBC_ROUNDING, objectiveValue, newSolution);
            }
            delete[] newSolution;
#ifdef CBC_THREAD
            if (master_) {
              lockThread();
              if (parallelMode() > 0) {
                while (master_->waitForThreadsInTree(0)) {
                  lockThread();
                  double dummyBest;
                  tree_->cleanTree(this, -COIN_DBL_MAX, dummyBest);
                  //unlockThread();
                }
              } else {
                double dummyBest;
                tree_->cleanTree(this, -COIN_DBL_MAX, dummyBest);
              }
              master_->waitForThreadsInTree(2);
              delete master_;
              master_ = NULL;
              masterThread_ = NULL;
            }
#endif
            if (tree_->size()) {
              double dummyBest;
              tree_->cleanTree(this, -COIN_DBL_MAX, dummyBest);
            }
            break;
          }
        }
        delete saveSolver;
        saveSolver = NULL;
      }
    }
    /*
          Check for abort on limits: node count, solution count, time, integrality gap.
        */
    if (!(numberNodes_ < intParam_[CbcMaxNumNode] && numberSolutions_ < intParam_[CbcMaxNumSol] && !maximumSecondsReached() && !stoppedOnGap_ && !eventHappened_ && (maximumNumberIterations_ < 0 || numberIterations_ < maximumNumberIterations_))) {
      // out of loop
      break;
    }
#ifdef BONMIN
    assert(!solverCharacteristics_->solutionAddsCuts() || solverCharacteristics_->mipFeasible());
#endif
// Sets percentage of time when we try diving. Diving requires a bit of heap reorganisation, because
// we need to replace the comparison function to dive, and that requires reordering to retain the
// heap property.
#define DIVE_WHEN 1000
#define DIVE_STOP 2000
    int kNode = numberNodes_ % 4000;
    if (numberNodes_ < 100000 && kNode > DIVE_WHEN && kNode <= DIVE_STOP) {
      if (!parallelMode()) {
        if (kNode == DIVE_WHEN + 1 || numberConsecutiveInfeasible > 1) {
          CbcCompareDefault *compare = dynamic_cast<CbcCompareDefault *>(nodeCompare_);
          // Don't interfere if user has replaced the compare function.
          if (compare) {
            //printf("Redoing tree\n");
            compare->startDive(this);
            numberConsecutiveInfeasible = 0;
          }
        }
      }
    }
    // replace current cutoff?
    if (cutoff > getCutoff()) {
      double newCutoff = getCutoff();
      if (analyzeResults_) {
        // see if we could fix any (more)
        int n = 0;
        double *newLower = analyzeResults_;
        double *objLower = newLower + numberIntegers_;
        double *newUpper = objLower + numberIntegers_;
        double *objUpper = newUpper + numberIntegers_;
        for (int i = 0; i < numberIntegers_; i++) {
          if (objLower[i] > newCutoff) {
            n++;
            if (objUpper[i] > newCutoff) {
              newCutoff = -COIN_DBL_MAX;
              break;
            }
            // add as global cut
            objLower[i] = -COIN_DBL_MAX;
            OsiRowCut rc;
            rc.setLb(newLower[i]);
            rc.setUb(COIN_DBL_MAX);
            double one = 1.0;
            rc.setRow(1, integerVariable_ + i, &one, false);
            rc.setGloballyValidAsInteger(2);
            globalCuts_.addCutIfNotDuplicate(rc);
          } else if (objUpper[i] > newCutoff) {
            n++;
            // add as global cut
            objUpper[i] = -COIN_DBL_MAX;
            OsiRowCut rc;
            rc.setLb(-COIN_DBL_MAX);
            rc.setUb(newUpper[i]);
            double one = 1.0;
            rc.setRow(1, integerVariable_ + i, &one, false);
            rc.setGloballyValidAsInteger(2);
            globalCuts_.addCutIfNotDuplicate(rc);
          }
        }
        if (newCutoff == -COIN_DBL_MAX) {
          COIN_DETAIL_PRINT(printf("Root analysis says finished\n"));
        } else if (n > numberFixedNow_) {
          COIN_DETAIL_PRINT(printf("%d more fixed by analysis - now %d\n", n - numberFixedNow_, n));
          numberFixedNow_ = n;
        }
      }
      if (eventHandler) {
        if (!eventHandler->event(CbcEventHandler::solution)) {
          eventHappened_ = true; // exit
        }
        newCutoff = getCutoff();
      }
      lockThread();
      /*
              Clean the tree to reflect the new solution, then see if the
              node comparison predicate wants to make any changes. If so,
              call setComparison for the side effect of rebuilding the heap.
            */
      tree_->cleanTree(this, newCutoff, bestPossibleObjective_);
      if (nodeCompare_->newSolution(this) || nodeCompare_->newSolution(this, continuousObjective_, continuousInfeasibilities_)) {
        tree_->setComparison(*nodeCompare_);
      }
      if (tree_->empty()) {
        continue;
      }
      unlockThread();
    }
    cutoff = getCutoff();
    /*
            Periodic activities: Opportunities to
            + tweak the nodeCompare criteria,
            + check if we've closed the integrality gap enough to quit,
            + print a summary line to let the user know we're working
        */
    if (numberNodes_ >= lastEvery1000) {
      lockThread();
#ifdef COIN_HAS_CLP
      // See if we want dantzig row choice
      goToDantzig(1000, savePivotMethod);
#endif
      lastEvery1000 = numberNodes_ + 1000;
      bool redoTree = nodeCompare_->every1000Nodes(this, numberNodes_);
#ifdef CHECK_CUT_SIZE
      verifyCutSize(tree_, *this);
#endif
      // redo tree if requested
      if (redoTree)
        tree_->setComparison(*nodeCompare_);
      unlockThread();
    }
    // Had hotstart before, now switched off
    if (saveCompare && !hotstartSolution_) {
      // hotstart switched off
      delete nodeCompare_; // off depth first
      nodeCompare_ = saveCompare;
      saveCompare = NULL;
      // redo tree
      lockThread();
      tree_->setComparison(*nodeCompare_);
      unlockThread();
    }
    if (numberNodes_ >= lastPrintEvery) {
      lastPrintEvery = numberNodes_ + printFrequency_;
      lockThread();
      int nNodes = tree_->size();

      //MODIF PIERRE
      bestPossibleObjective_ = tree_->getBestPossibleObjective();
#ifdef CBC_THREAD
      if (parallelMode() > 0 && master_) {
        // need to adjust for ones not on tree
        int numberThreads = master_->numberThreads();
        for (int i = 0; i < numberThreads; i++) {
          CbcThread *child = master_->child(i);
          if (child->node()) {
            // adjust
            double value = child->node()->objectiveValue();
            bestPossibleObjective_ = CoinMin(bestPossibleObjective_, value);
          }
        }
      }
#endif
      unlockThread();
#if CBC_USEFUL_PRINTING > 1
      if (getCutoff() < 1.0e20) {
        if (fabs(getCutoff() - (bestObjective_ - getCutoffIncrement())) > 1.0e-6 && !parentModel_)
          printf("model cutoff in status %g, best %g, increment %g\n",
            getCutoff(), bestObjective_, getCutoffIncrement());
        assert(getCutoff() < bestObjective_ - getCutoffIncrement() + 1.0e-6 + 1.0e-10 * fabs(bestObjective_));
      }
#endif
      if (!intParam_[CbcPrinting]) {
        // Parallel may not have any nodes
        if (!nNodes)
          bestPossibleObjective_ = lastBestPossibleObjective;
        else
          lastBestPossibleObjective = bestPossibleObjective_;
        messageHandler()->message(CBC_STATUS, messages())
          << numberNodes_ << CoinMax(nNodes, 1) << bestObjective_ << bestPossibleObjective_
          << getCurrentSeconds()
          << CoinMessageEol;
      } else if (intParam_[CbcPrinting] == 1) {
        messageHandler()->message(CBC_STATUS2, messages())
          << numberNodes_ << nNodes << bestObjective_ << bestPossibleObjective_
          << tree_->lastDepth() << tree_->lastUnsatisfied()
          << tree_->lastObjective() << numberIterations_
          << getCurrentSeconds()
          << CoinMessageEol;
      } else if (!numberExtraIterations_) {
        messageHandler()->message(CBC_STATUS2, messages())
          << numberNodes_ << nNodes << bestObjective_ << bestPossibleObjective_
          << tree_->lastDepth() << tree_->lastUnsatisfied() << numberIterations_
          << getCurrentSeconds()
          << CoinMessageEol;
      } else {
        messageHandler()->message(CBC_STATUS3, messages())
          << numberNodes_ << numberFathoms_ << numberExtraNodes_ << nNodes
          << bestObjective_ << bestPossibleObjective_
          << tree_->lastDepth() << tree_->lastUnsatisfied() << numberIterations_ << numberExtraIterations_
          << getCurrentSeconds()
          << CoinMessageEol;
      }
#ifdef COIN_HAS_NTY
      if (symmetryInfo_)
        symmetryInfo_->statsOrbits(this, 1);
#endif
#if PRINT_CONFLICT == 1
      if (numberConflictCuts > lastNumberConflictCuts) {
        double length = lengthConflictCuts / numberConflictCuts;
        printf("%d new conflict cuts - total %d - average length %g\n",
          numberConflictCuts - lastNumberConflictCuts,
          numberConflictCuts, length);
        lastNumberConflictCuts = numberConflictCuts;
      }
#endif
      if (eventHandler && !eventHandler->event(CbcEventHandler::treeStatus)) {
        eventHappened_ = true; // exit
      }
    }
    // See if can stop on gap
    if (canStopOnGap()) {
      stoppedOnGap_ = true;
    }

#ifdef CHECK_NODE_FULL
    verifyTreeNodes(tree_, *this);
#endif
#ifdef CHECK_CUT_COUNTS
    verifyCutCounts(tree_, *this);
#endif
    /*
          Now we come to the meat of the loop. To create the active subproblem, we'll
          pop the most promising node in the live set, rebuild the subproblem it
          represents, and then execute the current arm of the branch to create the
          active subproblem.
        */
    CbcNode *node = NULL;
#ifdef CBC_THREAD
    if (!parallelMode() || parallelMode() == -1) {
#endif
      node = tree_->bestNode(cutoff);
      // Possible one on tree worse than cutoff
      // Weird comparison function can leave ineligible nodes on tree
      if (!node || node->objectiveValue() > cutoff)
        continue;
      // Do main work of solving node here
      doOneNode(this, node, createdNode);
#ifdef JJF_ZERO
      if (node) {
        if (createdNode) {
          printf("Node %d depth %d, created %d depth %d\n",
            node->nodeNumber(), node->depth(),
            createdNode->nodeNumber(), createdNode->depth());
        } else {
          printf("Node %d depth %d,  no created node\n",
            node->nodeNumber(), node->depth());
        }
      } else if (createdNode) {
        printf("Node exhausted, created %d depth %d\n",
          createdNode->nodeNumber(), createdNode->depth());
      } else {
        printf("Node exhausted,  no created node\n");
        numberConsecutiveInfeasible = 2;
      }
#endif
      //if (createdNode)
      //numberConsecutiveInfeasible=0;
      //else
      //numberConsecutiveInfeasible++;
#ifdef CBC_THREAD
    } else if (parallelMode() > 0) {
      //lockThread();
      //node = tree_->bestNode(cutoff) ;
      // Possible one on tree worse than cutoff
      if (true || !node || node->objectiveValue() > cutoff) {
        assert(master_);
        if (master_) {
          int anyLeft = master_->waitForThreadsInTree(1);
          // may need to go round again
          if (anyLeft) {
            continue;
          } else {
            master_->stopThreads(-1);
          }
        }
      }
      //unlockThread();
    } else {
      // Deterministic parallel
      if ((tree_->size() < CoinMax(numberThreads_, 8) || hotstartSolution_) && !goneParallel) {
        node = tree_->bestNode(cutoff);
        // Possible one on tree worse than cutoff
        if (!node || node->objectiveValue() > cutoff)
          continue;
        // Do main work of solving node here
        doOneNode(this, node, createdNode);
        assert(createdNode);
        if (!createdNode->active()) {
          delete createdNode;
          createdNode = NULL;
        } else {
          // Say one more pointing to this
          node->nodeInfo()->increment();
          tree_->push(createdNode);
        }
        if (node->active()) {
          assert(node->nodeInfo());
          if (node->nodeInfo()->numberBranchesLeft()) {
            tree_->push(node);
          } else {
            node->setActive(false);
          }
        } else {
          if (node->nodeInfo()) {
            if (!node->nodeInfo()->numberBranchesLeft())
              node->nodeInfo()->allBranchesGone(); // can clean up
            // So will delete underlying stuff
            node->setActive(true);
          }
          delNode[nDeleteNode++] = node;
          node = NULL;
        }
        if (nDeleteNode >= MAX_DEL_NODE) {
          for (int i = 0; i < nDeleteNode; i++) {
            //printf("trying to del %d %x\n",i,delNode[i]);
            delete delNode[i];
            //printf("done to del %d %x\n",i,delNode[i]);
          }
          nDeleteNode = 0;
        }
      } else {
        // Split and solve
        master_->deterministicParallel();
        goneParallel = true;
      }
    }
#endif
  }
  if (nDeleteNode) {
    for (int i = 0; i < nDeleteNode; i++) {
      delete delNode[i];
    }
    nDeleteNode = 0;
  }
#ifdef CBC_THREAD
  if (master_) {
    master_->stopThreads(-1);
    master_->waitForThreadsInTree(2);
    // adjust time to allow for children on some systems
    //dblParam_[CbcStartSeconds] -= CoinCpuTimeJustChildren();
  }
#endif
  /*
      End of the non-abort actions. The next block of code is executed if we've
      aborted because we hit one of the limits. Clean up by deleting the live set
      and break out of the node processing loop. Note that on an abort, node may
      have been pushed back onto the tree for further processing, in which case
      it'll be deleted in cleanTree. We need to check.
    */
  if (!(numberNodes_ < intParam_[CbcMaxNumNode] && numberSolutions_ < intParam_[CbcMaxNumSol] && !maximumSecondsReached() && !stoppedOnGap_ && !eventHappened_ && (maximumNumberIterations_ < 0 || numberIterations_ < maximumNumberIterations_))) {
    if (tree_->size()) {
      double dummyBest;
      tree_->cleanTree(this, -COIN_DBL_MAX, dummyBest);
#if 0 // Does not seem to be needed def CBC_THREAD
            if (parallelMode() > 0 && master_) {
              // see if any dangling nodes
              int numberThreads = master_->numberThreads();
              for (int i=0;i<numberThreads;i++) {
                CbcThread * child = master_->child(i);
                //if (child->createdNode())
                //printf("CHILD_NODE %p\n",child->createdNode());
                delete child->createdNode();
              }
            }
#endif
    }
    delete nextRowCut_;
    /* order is important here:
         * maximumSecondsReached() should be checked before eventHappened_ and
         * isNodeLimitReached() should be checked after eventHappened_
         * reason is, that at timelimit, eventHappened_ is set to true to make Cbc stop fast
         *   and if Ctrl+C is hit, then the nodelimit is set to -1 to make Cbc stop
         */
    if (stoppedOnGap_) {
      messageHandler()->message(CBC_GAP, messages())
        << bestObjective_ - bestPossibleObjective_
        << dblParam_[CbcAllowableGap]
        << dblParam_[CbcAllowableFractionGap] * 100.0
        << CoinMessageEol;
      secondaryStatus_ = 2;
      status_ = 0;
    } else if (maximumSecondsReached()) {
      handler_->message(CBC_MAXTIME, messages_) << CoinMessageEol;
      secondaryStatus_ = 4;
      status_ = 1;
    } else if (numberSolutions_ >= intParam_[CbcMaxNumSol]) {
      handler_->message(CBC_MAXSOLS, messages_) << CoinMessageEol;
      secondaryStatus_ = 6;
      status_ = 1;
    } else if (isNodeLimitReached()) {
      handler_->message(CBC_MAXNODES, messages_) << CoinMessageEol;
      secondaryStatus_ = 3;
      status_ = 1;
    } else if (maximumNumberIterations_ >= 0 && numberIterations_ >= maximumNumberIterations_) {
      handler_->message(CBC_MAXITERS, messages_) << CoinMessageEol;
      secondaryStatus_ = 8;
      status_ = 1;
    } else {
      handler_->message(CBC_EVENT, messages_) << CoinMessageEol;
      secondaryStatus_ = 5;
      status_ = 5;
    }
  }
#ifdef CBC_THREAD
  if (master_) {
    delete master_;
    master_ = NULL;
    masterThread_ = NULL;
  }
#endif
  /*
      That's it, we've exhausted the search tree, or broken out of the loop because
      we hit some limit on evaluation.

      We may have got an intelligent tree so give it one more chance
    */
  // Tell solver we are not in Branch and Cut
  solver_->setHintParam(OsiDoInBranchAndCut, false, OsiHintDo, NULL);
  tree_->endSearch();
  //  If we did any sub trees - did we give up on any?
  if (numberStoppedSubTrees_)
    status_ = 1;
  numberNodes_ += numberExtraNodes_;
  numberIterations_ += numberExtraIterations_;
  if (eventHandler) {
    eventHandler->event(CbcEventHandler::endSearch);
  }
  if (!status_) {
    // Set best possible unless stopped on gap
    if (secondaryStatus_ != 2)
      bestPossibleObjective_ = bestObjective_;
    handler_->message(CBC_END_GOOD, messages_)
      << bestObjective_ << numberIterations_ << numberNodes_ << getCurrentSeconds()
      << CoinMessageEol;
  } else {
    handler_->message(CBC_END, messages_)
      << bestObjective_ << bestPossibleObjective_
      << numberIterations_ << numberNodes_ << getCurrentSeconds()
      << CoinMessageEol;
  }
  if ((moreSpecialOptions_ & 4194304) != 0) {
    // Conflict cuts
    int numberCuts = globalCuts_.sizeRowCuts();
    int nConflict = 0;
    double sizeConflict = 0.0;
    for (int i = 0; i < numberCuts; i++) {
      OsiRowCut2 *cut = globalCuts_.cut(i);
      if (cut->whichRow() == 1) {
        nConflict++;
        sizeConflict += cut->row().getNumElements();
      }
    }
    if (nConflict) {
      sizeConflict /= nConflict;
      char general[200];
      sprintf(general, "%d conflict cuts generated - average length %g",
        nConflict, sizeConflict);
      messageHandler()->message(CBC_GENERAL,
        messages())
        << general << CoinMessageEol;
    }
  }
  if (numberStrongIterations_)
    handler_->message(CBC_STRONG_STATS, messages_)
      << strongInfo_[0] << numberStrongIterations_ << strongInfo_[2]
      << strongInfo_[1] << CoinMessageEol;
  if (!numberExtraNodes_)
    handler_->message(CBC_OTHER_STATS, messages_)
      << maximumDepthActual_
      << numberDJFixed_ << CoinMessageEol;
  else
    handler_->message(CBC_OTHER_STATS2, messages_)
      << maximumDepthActual_
      << numberDJFixed_ << numberFathoms_ << numberExtraNodes_ << numberExtraIterations_
      << CoinMessageEol;
#ifdef COIN_HAS_NTY
  if (symmetryInfo_)
    symmetryInfo_->statsOrbits(this, 1);
#endif
  if (doStatistics == 100) {
    for (int i = 0; i < numberObjects_; i++) {
      CbcSimpleIntegerDynamicPseudoCost *obj = dynamic_cast<CbcSimpleIntegerDynamicPseudoCost *>(object_[i]);
      if (obj)
        obj->print();
    }
  }
  if (statistics_) {
    // report in some way
    int *lookup = new int[numberObjects_];
    int i;
    for (i = 0; i < numberObjects_; i++)
      lookup[i] = -1;
    bool goodIds = false; //true;
    for (i = 0; i < numberObjects_; i++) {
      int iColumn = object_[i]->columnNumber();
      if (iColumn >= 0 && iColumn < numberColumns) {
        if (lookup[i] == -1) {
          lookup[i] = iColumn;
        } else {
          goodIds = false;
          break;
        }
      } else {
        goodIds = false;
        break;
      }
    }
    if (!goodIds) {
      delete[] lookup;
      lookup = NULL;
    }
    if (doStatistics >= 3) {
      printf("  node parent depth column   value                    obj      inf\n");
      for (i = 0; i < numberNodes2_; i++) {
        statistics_[i]->print(lookup);
      }
    }
    if (doStatistics > 1) {
      // Find last solution
      int k;
      for (k = numberNodes2_ - 1; k >= 0; k--) {
        if (statistics_[k]->endingObjective() != COIN_DBL_MAX && !statistics_[k]->endingInfeasibility())
          break;
      }
      if (k >= 0) {
        int depth = statistics_[k]->depth();
        int *which = new int[depth + 1];
        for (i = depth; i >= 0; i--) {
          which[i] = k;
          k = statistics_[k]->parentNode();
        }
        printf("  node parent depth column   value                    obj      inf\n");
        for (i = 0; i <= depth; i++) {
          statistics_[which[i]]->print(lookup);
        }
        delete[] which;
      }
    }
    // now summary
    int maxDepth = 0;
    double averageSolutionDepth = 0.0;
    int numberSolutions = 0;
    double averageCutoffDepth = 0.0;
    double averageSolvedDepth = 0.0;
    int numberCutoff = 0;
    int numberDown = 0;
    int numberFirstDown = 0;
    double averageInfDown = 0.0;
    double averageObjDown = 0.0;
    int numberCutoffDown = 0;
    int numberUp = 0;
    int numberFirstUp = 0;
    double averageInfUp = 0.0;
    double averageObjUp = 0.0;
    int numberCutoffUp = 0;
    double averageNumberIterations1 = 0.0;
    double averageValue = 0.0;
    for (i = 0; i < numberNodes2_; i++) {
      int depth = statistics_[i]->depth();
      int way = statistics_[i]->way();
      double value = statistics_[i]->value();
      double startingObjective = statistics_[i]->startingObjective();
      int startingInfeasibility = statistics_[i]->startingInfeasibility();
      double endingObjective = statistics_[i]->endingObjective();
      int endingInfeasibility = statistics_[i]->endingInfeasibility();
      maxDepth = CoinMax(depth, maxDepth);
      // Only for completed
      averageNumberIterations1 += statistics_[i]->numberIterations();
      averageValue += value;
      if (endingObjective != COIN_DBL_MAX && !endingInfeasibility) {
        numberSolutions++;
        averageSolutionDepth += depth;
      }
      if (endingObjective == COIN_DBL_MAX) {
        numberCutoff++;
        averageCutoffDepth += depth;
        if (way < 0) {
          numberDown++;
          numberCutoffDown++;
          if (way == -1)
            numberFirstDown++;
        } else {
          numberUp++;
          numberCutoffUp++;
          if (way == 1)
            numberFirstUp++;
        }
      } else {
        averageSolvedDepth += depth;
        if (way < 0) {
          numberDown++;
          averageInfDown += startingInfeasibility - endingInfeasibility;
          averageObjDown += endingObjective - startingObjective;
          if (way == -1)
            numberFirstDown++;
        } else {
          numberUp++;
          averageInfUp += startingInfeasibility - endingInfeasibility;
          averageObjUp += endingObjective - startingObjective;
          if (way == 1)
            numberFirstUp++;
        }
      }
    }
    // Now print
    if (numberSolutions)
      averageSolutionDepth /= static_cast<double>(numberSolutions);
    int numberSolved = numberNodes2_ - numberCutoff;
    double averageNumberIterations2 = numberIterations_ - averageNumberIterations1
      - numberIterationsAtContinuous;
    if (numberCutoff) {
      averageCutoffDepth /= static_cast<double>(numberCutoff);
      averageNumberIterations2 /= static_cast<double>(numberCutoff);
    }
    if (numberNodes2_)
      averageValue /= static_cast<double>(numberNodes2_);
    if (numberSolved) {
      averageNumberIterations1 /= static_cast<double>(numberSolved);
      averageSolvedDepth /= static_cast<double>(numberSolved);
    }
    printf("%d solution(s) were found (by branching) at an average depth of %g\n",
      numberSolutions, averageSolutionDepth);
    printf("average value of variable being branched on was %g\n",
      averageValue);
    printf("%d nodes were cutoff at an average depth of %g with iteration count of %g\n",
      numberCutoff, averageCutoffDepth, averageNumberIterations2);
    printf("%d nodes were solved at an average depth of %g with iteration count of %g\n",
      numberSolved, averageSolvedDepth, averageNumberIterations1);
    if (numberDown) {
      averageInfDown /= static_cast<double>(numberDown);
      averageObjDown /= static_cast<double>(numberDown);
    }
    printf("Down %d nodes (%d first, %d second) - %d cutoff, rest decrease numinf %g increase obj %g\n",
      numberDown, numberFirstDown, numberDown - numberFirstDown, numberCutoffDown,
      averageInfDown, averageObjDown);
    if (numberUp) {
      averageInfUp /= static_cast<double>(numberUp);
      averageObjUp /= static_cast<double>(numberUp);
    }
    printf("Up %d nodes (%d first, %d second) - %d cutoff, rest decrease numinf %g increase obj %g\n",
      numberUp, numberFirstUp, numberUp - numberFirstUp, numberCutoffUp,
      averageInfUp, averageObjUp);
    for (i = 0; i < numberNodes2_; i++)
      delete statistics_[i];
    delete[] statistics_;
    statistics_ = NULL;
    maximumStatistics_ = 0;
    delete[] lookup;
  }
  /*
      If we think we have a solution, restore and confirm it with a call to
      setBestSolution().  We need to reset the cutoff value so as not to fathom
      the solution on bounds.  Note that calling setBestSolution( ..., true)
      leaves the continuousSolver_ bounds vectors fixed at the solution value.

      Running resolve() here is a failsafe --- setBestSolution has already
      reoptimised using the continuousSolver_. If for some reason we fail to
      prove optimality, run the problem again after instructing the solver to
      tell us more.

      If all looks good, replace solver_ with continuousSolver_, so that the
      outside world will be able to obtain information about the solution using
      public methods.

      Don't replace if we are trying to save cuts
    */
  if (bestSolution_ && (solverCharacteristics_->solverType() < 2 || solverCharacteristics_->solverType() == 4) && ((specialOptions_ & 8388608) == 0 || (specialOptions_ & 2048) != 0)) {
    setCutoff(1.0e50); // As best solution should be worse than cutoff
    // change cutoff as constraint if wanted
    if (cutoffRowNumber_ >= 0) {
      if (solver_->getNumRows() > cutoffRowNumber_)
        solver_->setRowUpper(cutoffRowNumber_, 1.0e50);
    }
    // also in continuousSolver_
    if (continuousSolver_) {
      // Solvers know about direction
      double direction = solver_->getObjSense();
      continuousSolver_->setDblParam(OsiDualObjectiveLimit, 1.0e50 * direction);
    }
    phase_ = 5;
    double increment = getDblParam(CbcModel::CbcCutoffIncrement);
    if ((specialOptions_ & 4) == 0)
      bestObjective_ += 100.0 * increment + 1.0e-3; // only set if we are going to solve
    setBestSolution(CBC_END_SOLUTION, bestObjective_, bestSolution_, 1);
    continuousSolver_->resolve();
    // Deal with funny variables
    if ((moreSpecialOptions2_ & 32768) != 0)
      cleanBounds(continuousSolver_, NULL);
    if (!continuousSolver_->isProvenOptimal()) {
      continuousSolver_->messageHandler()->setLogLevel(2);
      continuousSolver_->initialSolve();
    }
    delete solver_;
    // above deletes solverCharacteristics_
    solverCharacteristics_ = NULL;
    solver_ = continuousSolver_;
    setPointers(solver_);
    continuousSolver_ = NULL;
  }
  /*
      Clean up dangling objects. continuousSolver_ may already be toast.
    */
  delete lastws;
  if (saveObjects) {
    for (int i = 0; i < numberObjects_; i++)
      delete saveObjects[i];
    delete[] saveObjects;
  }
  numberStrong_ = saveNumberStrong;
  numberBeforeTrust_ = saveNumberBeforeTrust;
  delete[] whichGenerator_;
  whichGenerator_ = NULL;
  delete[] lowerBefore;
  delete[] upperBefore;
  delete[] walkback_;
  walkback_ = NULL;
  delete[] lastNodeInfo_;
  lastNodeInfo_ = NULL;
  delete[] lastNumberCuts_;
  lastNumberCuts_ = NULL;
  delete[] lastCut_;
  lastCut_ = NULL;
  delete[] addedCuts_;
  addedCuts_ = NULL;
  //delete persistentInfo;
  // Get rid of characteristics
  solverCharacteristics_ = NULL;
  if (continuousSolver_) {
    delete continuousSolver_;
    continuousSolver_ = NULL;
  }
  /*
      Destroy global cuts by replacing with an empty OsiCuts object.
    */
  globalCuts_ = CbcRowCuts();
  delete globalConflictCuts_;
  globalConflictCuts_ = NULL;
  if (!bestSolution_ && (specialOptions_ & 8388608) == 0 && false) {
    // make sure lp solver is infeasible
    int numberColumns = solver_->getNumCols();
    const double *columnLower = solver_->getColLower();
    int iColumn;
    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
      if (solver_->isInteger(iColumn))
        solver_->setColUpper(iColumn, columnLower[iColumn]);
    }
    solver_->initialSolve();
  }
#ifdef COIN_HAS_CLP
  {
    OsiClpSolverInterface *clpSolver
      = dynamic_cast<OsiClpSolverInterface *>(solver_);
    if (clpSolver) {
      // Possible restore of pivot method
      if (savePivotMethod) {
        // model may have changed
        savePivotMethod->setModel(NULL);
        clpSolver->getModelPtr()->setDualRowPivotAlgorithm(*savePivotMethod);
        delete savePivotMethod;
      }
      clpSolver->setLargestAway(-1.0);
    }
  }
#endif
  if ((fastNodeDepth_ >= 1000000 || (moreSpecialOptions_ & 33554432) != 0)
    && !parentModel_) {
    // delete object off end
    delete object_[numberObjects_];
    if ((moreSpecialOptions_ & 33554432) == 0)
      fastNodeDepth_ -= 1000000;
  }
  delete saveSolver;
  // Undo preprocessing performed during BaB.
  if (strategy_ && strategy_->preProcessState() > 0) {
    // undo preprocessing
    CglPreProcess *process = strategy_->process();
    assert(process);
    int n = originalSolver->getNumCols();
    if (bestSolution_) {
      delete[] bestSolution_;
      bestSolution_ = new double[n];
      process->postProcess(*solver_);
    }
    strategy_->deletePreProcess();
    // Solution now back in originalSolver
    delete solver_;
    solver_ = originalSolver;
    if (bestSolution_) {
      bestObjective_ = solver_->getObjValue() * solver_->getObjSense();
      memcpy(bestSolution_, solver_->getColSolution(), n * sizeof(double));
    }
    // put back original objects if there were any
    if (originalObject) {
      int iColumn;
      assert(ownObjects_);
      for (iColumn = 0; iColumn < numberObjects_; iColumn++)
        delete object_[iColumn];
      delete[] object_;
      numberObjects_ = numberOriginalObjects;
      object_ = originalObject;
      delete[] integerVariable_;
      numberIntegers_ = 0;
      for (iColumn = 0; iColumn < n; iColumn++) {
        if (solver_->isInteger(iColumn))
          numberIntegers_++;
      }
      integerVariable_ = new int[numberIntegers_];
      numberIntegers_ = 0;
      for (iColumn = 0; iColumn < n; iColumn++) {
        if (solver_->isInteger(iColumn))
          integerVariable_[numberIntegers_++] = iColumn;
      }
    }
  }
  if (flipObjective)
    flipModel();
#ifdef COIN_HAS_CLP
  {
    OsiClpSolverInterface *clpSolver
      = dynamic_cast<OsiClpSolverInterface *>(solver_);
    if (clpSolver)
      clpSolver->setFakeObjective(reinterpret_cast<double *>(NULL));
  }
#endif
  moreSpecialOptions_ = saveMoreSpecialOptions;
  return;
}

// Solve the initial LP relaxation
void CbcModel::initialSolve()
{
  assert(solver_);
  // Double check optimization directions line up
  dblParam_[CbcOptimizationDirection] = solver_->getObjSense();
  // Check if bounds are all integral (as may get messed up later)
  checkModel();
  if (!solverCharacteristics_) {
    OsiBabSolver *solverCharacteristics = dynamic_cast<OsiBabSolver *>(solver_->getAuxiliaryInfo());
    if (solverCharacteristics) {
      solverCharacteristics_ = solverCharacteristics;
    } else {
      // replace in solver
      OsiBabSolver defaultC;
      solver_->setAuxiliaryInfo(&defaultC);
      solverCharacteristics_ = dynamic_cast<OsiBabSolver *>(solver_->getAuxiliaryInfo());
    }
  }
  solverCharacteristics_->setSolver(solver_);
  solver_->setHintParam(OsiDoInBranchAndCut, true, OsiHintDo, NULL);
  solver_->initialSolve();
  solver_->setHintParam(OsiDoInBranchAndCut, false, OsiHintDo, NULL);
  if (!solver_->isProvenOptimal())
    solver_->resolve();
  // But set up so Jon Lee will be happy
  status_ = -1;
  secondaryStatus_ = -1;
  originalContinuousObjective_ = solver_->getObjValue() * solver_->getObjSense();
  bestPossibleObjective_ = originalContinuousObjective_;
  if (solver_->isProvenDualInfeasible())
    originalContinuousObjective_ = -COIN_DBL_MAX;