// Copyright (C) 2006, International Business Machines
// Corporation and others.  All Rights Reserved.
// This code is licensed under the terms of the Eclipse Public License (EPL).

#if defined(_MSC_VER)
// Turn off compiler warning about long names
#  pragma warning(disable:4786)
#endif
#include <cassert>
#include <cstdlib>
#include <cmath>
#include <cfloat>
//#define OSI_DEBUG
#include "OsiSolverInterface.hpp"
#include "OsiBranchingObject.hpp"
#include "CoinHelperFunctions.hpp"
#include "CoinPackedMatrix.hpp"
#include "CoinSort.hpp"
#include "CoinError.hpp"
#include "CoinFinite.hpp"

// Default Constructor
OsiObject::OsiObject() 
  :infeasibility_(0.0),
   whichWay_(0),
   numberWays_(2),
   priority_(1000)
{
}


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

// Copy constructor 
OsiObject::OsiObject ( const OsiObject & rhs)
{
  infeasibility_ = rhs.infeasibility_;
  whichWay_ = rhs.whichWay_;
  priority_ = rhs.priority_;
  numberWays_ = rhs.numberWays_;
}

// Assignment operator 
OsiObject & 
OsiObject::operator=( const OsiObject& rhs)
{
  if (this!=&rhs) {
    infeasibility_ = rhs.infeasibility_;
    whichWay_ = rhs.whichWay_;
    priority_ = rhs.priority_;
    numberWays_ = rhs.numberWays_;
  }
  return *this;
}
// Return "up" estimate (default 1.0e-5)
double 
OsiObject::upEstimate() const
{
  return 1.0e-5;
}
// Return "down" estimate (default 1.0e-5)
double 
OsiObject::downEstimate() const
{
  return 1.0e-5;
}
// Column number if single column object -1 otherwise
int 
OsiObject::columnNumber() const
{
  return -1;
}
// Infeasibility - large is 0.5
double 
OsiObject::infeasibility(const OsiSolverInterface * solver, int & preferredWay) const
{
  // Can't guarantee has matrix
  OsiBranchingInformation info(solver,false,false);
  return infeasibility(&info,preferredWay);
}
// This does NOT set mutable stuff
double 
OsiObject::checkInfeasibility(const OsiBranchingInformation * info) const
{
  int way;
  double saveInfeasibility = infeasibility_;
  short int saveWhichWay = whichWay_ ;
  double value = infeasibility(info,way);
  infeasibility_ = saveInfeasibility;
  whichWay_ = saveWhichWay;
  return value;
}

/* For the variable(s) referenced by the object,
   look at the current solution and set bounds to match the solution.
   Returns measure of how much it had to move solution to make feasible
*/
double 
OsiObject::feasibleRegion(OsiSolverInterface * solver) const 
{
  // Can't guarantee has matrix
  OsiBranchingInformation info(solver,false,false);
  return feasibleRegion(solver,&info);
}

// Default Constructor
OsiObject2::OsiObject2() 
  : OsiObject(),
    preferredWay_(-1),
    otherInfeasibility_(0.0)
{
}


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

// Copy constructor 
OsiObject2::OsiObject2 ( const OsiObject2 & rhs)
  : OsiObject(rhs),
    preferredWay_(rhs.preferredWay_),
    otherInfeasibility_ (rhs.otherInfeasibility_)
{
}

// Assignment operator 
OsiObject2 & 
OsiObject2::operator=( const OsiObject2& rhs)
{
  if (this!=&rhs) {
    OsiObject::operator=(rhs);
    preferredWay_ = rhs.preferredWay_;
    otherInfeasibility_ = rhs.otherInfeasibility_;
  }
  return *this;
}
// Default Constructor 
OsiBranchingObject::OsiBranchingObject()
{
  originalObject_=NULL;
  branchIndex_=0;
  value_=0.0;
  numberBranches_=2;
}

// Useful constructor
OsiBranchingObject::OsiBranchingObject (OsiSolverInterface * ,
					 double value)
{
  originalObject_=NULL;
  branchIndex_=0;
  value_=value;
  numberBranches_=2;
}

// Copy constructor 
OsiBranchingObject::OsiBranchingObject ( const OsiBranchingObject & rhs)
{
  originalObject_=rhs.originalObject_;
  branchIndex_=rhs.branchIndex_;
  value_=rhs.value_;
  numberBranches_=rhs.numberBranches_;
}

// Assignment operator 
OsiBranchingObject & 
OsiBranchingObject::operator=( const OsiBranchingObject& rhs)
{
  if (this != &rhs) {
    originalObject_=rhs.originalObject_;
    branchIndex_=rhs.branchIndex_;
    value_=rhs.value_;
    numberBranches_=rhs.numberBranches_;
  }
  return *this;
}

// Destructor 
OsiBranchingObject::~OsiBranchingObject ()
{
}
// For debug
int 
OsiBranchingObject::columnNumber() const
{
  if (originalObject_)
    return originalObject_->columnNumber();
  else
    return -1;
}
/** Default Constructor

*/
OsiBranchingInformation::OsiBranchingInformation ()
  : objectiveValue_(COIN_DBL_MAX),
    cutoff_(COIN_DBL_MAX),
    direction_(COIN_DBL_MAX),
    integerTolerance_(1.0e-7),
    primalTolerance_(1.0e-7),
    timeRemaining_(COIN_DBL_MAX),
    defaultDual_(-1.0),
    solver_(NULL),
    numberColumns_(0),
    lower_(NULL),
    solution_(NULL),
    upper_(NULL),
    hotstartSolution_(NULL),
    pi_(NULL),
    rowActivity_(NULL),
    objective_(NULL),
    rowLower_(NULL),
    rowUpper_(NULL),
    elementByColumn_(NULL),
    columnStart_(NULL),
    columnLength_(NULL),
    row_(NULL),
    usefulRegion_(NULL),
    indexRegion_(NULL),
    numberSolutions_(0),
    numberBranchingSolutions_(0),
    depth_(0),
    owningSolution_(false)
{
}

/** Useful constructor
*/
OsiBranchingInformation::OsiBranchingInformation (const OsiSolverInterface * solver,
						  bool /*normalSolver*/,
						  bool owningSolution)
  : timeRemaining_(COIN_DBL_MAX),
    defaultDual_(-1.0),
    solver_(solver),
    hotstartSolution_(NULL),
    usefulRegion_(NULL),
    indexRegion_(NULL),
    numberSolutions_(0),
    numberBranchingSolutions_(0),
    depth_(0),
    owningSolution_(owningSolution)
{
  direction_ = solver_->getObjSense();
  objectiveValue_ = solver_->getObjValue();
  objectiveValue_ *= direction_;
  solver_->getDblParam(OsiDualObjectiveLimit,cutoff_) ;
  cutoff_ *= direction_;
  integerTolerance_ = solver_->getIntegerTolerance();
  solver_->getDblParam(OsiPrimalTolerance,primalTolerance_) ;
  numberColumns_ = solver_->getNumCols();
  lower_ = solver_->getColLower();
  if (owningSolution_)
    solution_ = CoinCopyOfArray(solver_->getColSolution(),numberColumns_);
  else
    solution_ = solver_->getColSolution();
  upper_ = solver_->getColUpper();
  pi_ = solver_->getRowPrice();
  rowActivity_ = solver_->getRowActivity();
  objective_ = solver_->getObjCoefficients();
  rowLower_ = solver_->getRowLower();
  rowUpper_ = solver_->getRowUpper();
  const CoinPackedMatrix* matrix = solver_->getMatrixByCol();
  if (matrix) {
    // Column copy of matrix if matrix exists
    elementByColumn_ = matrix->getElements();
    row_ = matrix->getIndices();
    columnStart_ = matrix->getVectorStarts();
    columnLength_ = matrix->getVectorLengths();
  } else {
    // Matrix does not exist
    elementByColumn_ = NULL;
    row_ = NULL;
    columnStart_ = NULL;
    columnLength_ = NULL;
  }
}
// Copy constructor 
OsiBranchingInformation::OsiBranchingInformation ( const OsiBranchingInformation & rhs)
{
  objectiveValue_ = rhs.objectiveValue_;
  cutoff_ = rhs.cutoff_;
  direction_ = rhs.direction_;
  integerTolerance_ = rhs.integerTolerance_;
  primalTolerance_ = rhs.primalTolerance_;
  timeRemaining_ = rhs.timeRemaining_;
  defaultDual_ = rhs.defaultDual_;
  solver_ = rhs.solver_;
  numberColumns_ = rhs.numberColumns_;
  lower_ = rhs.lower_;
  owningSolution_ = rhs.owningSolution_;
  if (owningSolution_)
    solution_ = CoinCopyOfArray(rhs.solution_,numberColumns_);
  else
    solution_ = rhs.solution_;
  upper_ = rhs.upper_;
  hotstartSolution_ = rhs.hotstartSolution_;
  pi_ = rhs.pi_;
  rowActivity_ = rhs.rowActivity_;
  objective_ = rhs.objective_;
  rowLower_ = rhs.rowLower_;
  rowUpper_ = rhs.rowUpper_;
  elementByColumn_ = rhs.elementByColumn_;
  row_ = rhs.row_;
  columnStart_ = rhs.columnStart_;
  columnLength_ = rhs.columnLength_;
  usefulRegion_ = rhs.usefulRegion_;
  assert (!usefulRegion_);
  indexRegion_ = rhs.indexRegion_;
  numberSolutions_ = rhs.numberSolutions_;
  numberBranchingSolutions_ = rhs.numberBranchingSolutions_;
  depth_ = rhs.depth_;
}

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

// Assignment operator 
OsiBranchingInformation & 
OsiBranchingInformation::operator=( const OsiBranchingInformation& rhs)
{
  if (this!=&rhs) {
    objectiveValue_ = rhs.objectiveValue_;
    cutoff_ = rhs.cutoff_;
    direction_ = rhs.direction_;
    integerTolerance_ = rhs.integerTolerance_;
    primalTolerance_ = rhs.primalTolerance_;
    timeRemaining_ = rhs.timeRemaining_;
    defaultDual_ = rhs.defaultDual_;
    numberColumns_ = rhs.numberColumns_;
    lower_ = rhs.lower_;
    owningSolution_ = rhs.owningSolution_;
    if (owningSolution_) {
      solution_ = CoinCopyOfArray(rhs.solution_,numberColumns_);
      delete [] solution_;
    } else {
      solution_ = rhs.solution_;
    }
    upper_ = rhs.upper_;
    hotstartSolution_ = rhs.hotstartSolution_;
    pi_ = rhs.pi_;
    rowActivity_ = rhs.rowActivity_;
    objective_ = rhs.objective_;
    rowLower_ = rhs.rowLower_;
    rowUpper_ = rhs.rowUpper_;
    elementByColumn_ = rhs.elementByColumn_;
    row_ = rhs.row_;
    columnStart_ = rhs.columnStart_;
    columnLength_ = rhs.columnLength_;
    usefulRegion_ = rhs.usefulRegion_;
    assert (!usefulRegion_);
    indexRegion_ = rhs.indexRegion_;
    numberSolutions_ = rhs.numberSolutions_;
    numberBranchingSolutions_ = rhs.numberBranchingSolutions_;
    depth_ = rhs.depth_;
  }
  return *this;
}

// Destructor 
OsiBranchingInformation::~OsiBranchingInformation ()
{
  if (owningSolution_) 
    delete[] solution_;
}
// Default Constructor 
OsiTwoWayBranchingObject::OsiTwoWayBranchingObject()
  :OsiBranchingObject()
{
  firstBranch_=0;
}

// Useful constructor
OsiTwoWayBranchingObject::OsiTwoWayBranchingObject (OsiSolverInterface * solver, 
						      const OsiObject * object,
						      int way , double value)
  :OsiBranchingObject(solver,value)
{
  originalObject_ = object;
  firstBranch_=way;
}
  

// Copy constructor 
OsiTwoWayBranchingObject::OsiTwoWayBranchingObject ( const OsiTwoWayBranchingObject & rhs) :OsiBranchingObject(rhs)
{
  firstBranch_=rhs.firstBranch_;
}

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

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

/********* Simple Integers *******************************/
/** Default Constructor

  Equivalent to an unspecified binary variable.
*/
OsiSimpleInteger::OsiSimpleInteger ()
  : OsiObject2(),
    originalLower_(0.0),
    originalUpper_(1.0),
    columnNumber_(-1)
{
}

/** Useful constructor

  Loads actual upper & lower bounds for the specified variable.
*/
OsiSimpleInteger::OsiSimpleInteger (const OsiSolverInterface * solver, int iColumn)
  : OsiObject2()
{
  columnNumber_ = iColumn ;
  originalLower_ = solver->getColLower()[columnNumber_] ;
  originalUpper_ = solver->getColUpper()[columnNumber_] ;
}

  
// Useful constructor - passed solver index and original bounds
OsiSimpleInteger::OsiSimpleInteger ( int iColumn, double lower, double upper)
  : OsiObject2()
{
  columnNumber_ = iColumn ;
  originalLower_ = lower;
  originalUpper_ = upper;
}

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

{
  columnNumber_ = rhs.columnNumber_;
  originalLower_ = rhs.originalLower_;
  originalUpper_ = rhs.originalUpper_;
}

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

// Assignment operator 
OsiSimpleInteger & 
OsiSimpleInteger::operator=( const OsiSimpleInteger& rhs)
{
  if (this!=&rhs) {
    OsiObject2::operator=(rhs);
    columnNumber_ = rhs.columnNumber_;
    originalLower_ = rhs.originalLower_;
    originalUpper_ = rhs.originalUpper_;
  }
  return *this;
}

// Destructor 
OsiSimpleInteger::~OsiSimpleInteger ()
{
}
/* Reset variable bounds to their original values.
   
Bounds may be tightened, so it may be good to be able to reset them to
their original values.
*/
void 
OsiSimpleInteger::resetBounds(const OsiSolverInterface * solver) 
{
  originalLower_ = solver->getColLower()[columnNumber_] ;
  originalUpper_ = solver->getColUpper()[columnNumber_] ;
}
// Redoes data when sequence numbers change
void 
OsiSimpleInteger::resetSequenceEtc(int numberColumns, const int * originalColumns)
{
  int i;
  for (i=0;i<numberColumns;i++) {
    if (originalColumns[i]==columnNumber_)
      break;
  }
  if (i<numberColumns)
    columnNumber_=i;
  else
    abort(); // should never happen
}

// Infeasibility - large is 0.5
double 
OsiSimpleInteger::infeasibility(const OsiBranchingInformation * info, int & whichWay) const
{
  double value = info->solution_[columnNumber_];
  value = CoinMax(value, info->lower_[columnNumber_]);
  value = CoinMin(value, info->upper_[columnNumber_]);
  double nearest = floor(value+(1.0-0.5));
  if (nearest>value) { 
    whichWay=1;
  } else {
    whichWay=0;
  }
  infeasibility_ = fabs(value-nearest);
  double returnValue = infeasibility_;
  if (infeasibility_<=info->integerTolerance_) {
    otherInfeasibility_ = 1.0;
    returnValue = 0.0;
  } else if (info->defaultDual_<0.0) {
    otherInfeasibility_ = 1.0-infeasibility_;
  } else {
    const double * pi = info->pi_;
    const double * activity = info->rowActivity_;
    const double * lower = info->rowLower_;
    const double * upper = info->rowUpper_;
    const double * element = info->elementByColumn_;
    const int * row = info->row_;
    const CoinBigIndex * columnStart = info->columnStart_;
    const int * columnLength = info->columnLength_;
    double direction = info->direction_;
    double downMovement = value - floor(value);
    double upMovement = 1.0-downMovement;
    double valueP = info->objective_[columnNumber_]*direction;
    CoinBigIndex start = columnStart[columnNumber_];
    CoinBigIndex end = start + columnLength[columnNumber_];
    double upEstimate = 0.0;
    double downEstimate = 0.0;
    if (valueP>0.0)
      upEstimate = valueP*upMovement;
    else
      downEstimate -= valueP*downMovement;
    double tolerance = info->primalTolerance_;
    for (CoinBigIndex j=start;j<end;j++) {
      int iRow = row[j];
      if (lower[iRow]<-1.0e20) 
	assert (pi[iRow]<=1.0e-4);
      if (upper[iRow]>1.0e20) 
	assert (pi[iRow]>=-1.0e-4);
      valueP = pi[iRow]*direction;
      double el2 = element[j];
      double value2 = valueP*el2;
      double u=0.0;
      double d=0.0;
      if (value2>0.0)
	u = value2;
      else
	d = -value2;
      // if up makes infeasible then make at least default
      double newUp = activity[iRow] + upMovement*el2;
      if (newUp>upper[iRow]+tolerance||newUp<lower[iRow]-tolerance)
	u = CoinMax(u,info->defaultDual_);
      upEstimate += u*upMovement;
      // if down makes infeasible then make at least default
      double newDown = activity[iRow] - downMovement*el2;
      if (newDown>upper[iRow]+tolerance||newDown<lower[iRow]-tolerance)
	d = CoinMax(d,info->defaultDual_);
      downEstimate += d*downMovement;
    }
    if (downEstimate>=upEstimate) {
      infeasibility_ = CoinMax(1.0e-12,upEstimate);
      otherInfeasibility_ = CoinMax(1.0e-12,downEstimate);
      whichWay = 1;
    } else {
      infeasibility_ = CoinMax(1.0e-12,downEstimate);
      otherInfeasibility_ = CoinMax(1.0e-12,upEstimate);
      whichWay = 0;
    }
    returnValue = infeasibility_;
  }
  if (preferredWay_>=0&&returnValue)
    whichWay = preferredWay_;
  whichWay_ = static_cast<short int>(whichWay) ;
  return returnValue;
}

// This looks at solution and sets bounds to contain solution
/** More precisely: it first forces the variable within the existing
    bounds, and then tightens the bounds to fix the variable at the
    nearest integer value.
*/
double
OsiSimpleInteger::feasibleRegion(OsiSolverInterface * solver,
				 const OsiBranchingInformation * info) const
{
  double value = info->solution_[columnNumber_];
  double newValue = CoinMax(value, info->lower_[columnNumber_]);
  newValue = CoinMin(newValue, info->upper_[columnNumber_]);
  newValue = floor(newValue+0.5);
  solver->setColLower(columnNumber_,newValue);
  solver->setColUpper(columnNumber_,newValue);
  return fabs(value-newValue);
}
/* Column number if single column object -1 otherwise,
   so returns >= 0
   Used by heuristics
*/
int 
OsiSimpleInteger::columnNumber() const
{
  return columnNumber_;
}
// Creates a branching object
OsiBranchingObject * 
OsiSimpleInteger::createBranch(OsiSolverInterface * solver, const OsiBranchingInformation * info, int way) const 
{
  double value = info->solution_[columnNumber_];
  value = CoinMax(value, info->lower_[columnNumber_]);
  value = CoinMin(value, info->upper_[columnNumber_]);
  assert (info->upper_[columnNumber_]>info->lower_[columnNumber_]);
#ifndef NDEBUG
  double nearest = floor(value+0.5);
  assert (fabs(value-nearest)>info->integerTolerance_);
#endif
  OsiBranchingObject * branch = new OsiIntegerBranchingObject(solver,this,way,
					     value);
  return branch;
}
// Return "down" estimate
double 
OsiSimpleInteger::downEstimate() const
{
  if (whichWay_)
    return 1.0-infeasibility_;
  else
    return infeasibility_;
}
// Return "up" estimate
double 
OsiSimpleInteger::upEstimate() const
{
  if (!whichWay_)
    return 1.0-infeasibility_;
  else
    return infeasibility_;
}

// Default Constructor 
OsiIntegerBranchingObject::OsiIntegerBranchingObject()
  :OsiTwoWayBranchingObject()
{
  down_[0] = 0.0;
  down_[1] = 0.0;
  up_[0] = 0.0;
  up_[1] = 0.0;
}

// Useful constructor
OsiIntegerBranchingObject::OsiIntegerBranchingObject (OsiSolverInterface * solver, 
						      const OsiSimpleInteger * object,
						      int way , double value)
  :OsiTwoWayBranchingObject(solver,object, way, value)
{
  int iColumn = object->columnNumber();
  down_[0] = solver->getColLower()[iColumn];
  down_[1] = floor(value_);
  up_[0] = ceil(value_);
  up_[1] = solver->getColUpper()[iColumn];
}
/* Create a standard floor/ceiling branch object
   Specifies a simple two-way branch in a more flexible way. One arm of the
   branch will be lb <= x <= downUpperBound, the other upLowerBound <= x <= ub.
   Specify way = -1 to set the object state to perform the down arm first,
   way = 1 for the up arm.
*/
OsiIntegerBranchingObject::OsiIntegerBranchingObject (OsiSolverInterface * solver, 
						      const OsiSimpleInteger * object,
						      int way , double value, double downUpperBound, 
						      double upLowerBound) 
  :OsiTwoWayBranchingObject(solver,object, way, value)
{
  int iColumn = object->columnNumber();
  down_[0] = solver->getColLower()[iColumn];
  down_[1] = downUpperBound;
  up_[0] = upLowerBound;
  up_[1] = solver->getColUpper()[iColumn];
}
  

// Copy constructor 
OsiIntegerBranchingObject::OsiIntegerBranchingObject ( const OsiIntegerBranchingObject & rhs) :OsiTwoWayBranchingObject(rhs)
{
  down_[0] = rhs.down_[0];
  down_[1] = rhs.down_[1];
  up_[0] = rhs.up_[0];
  up_[1] = rhs.up_[1];
}

// Assignment operator 
OsiIntegerBranchingObject & 
OsiIntegerBranchingObject::operator=( const OsiIntegerBranchingObject& rhs)
{
  if (this != &rhs) {
    OsiTwoWayBranchingObject::operator=(rhs);
    down_[0] = rhs.down_[0];
    down_[1] = rhs.down_[1];
    up_[0] = rhs.up_[0];
    up_[1] = rhs.up_[1];
  }
  return *this;
}
OsiBranchingObject * 
OsiIntegerBranchingObject::clone() const
{ 
  return (new OsiIntegerBranchingObject(*this));
}


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

/*
  Perform a branch by adjusting the bounds of the specified variable. Note
  that each arm of the branch advances the object to the next arm by
  advancing the value of branchIndex_.

  Providing new values for the variable's lower and upper bounds for each
  branching direction gives a little bit of additional flexibility and will
  be easily extensible to multi-way branching.
  Returns change in guessed objective on next branch
*/
double
OsiIntegerBranchingObject::branch(OsiSolverInterface * solver)
{
  const OsiSimpleInteger * obj =
    dynamic_cast <const OsiSimpleInteger *>(originalObject_) ;
  assert (obj);
  int iColumn = obj->columnNumber();
  double olb,oub ;
  olb = solver->getColLower()[iColumn] ;
  oub = solver->getColUpper()[iColumn] ;
  int way = (!branchIndex_) ? (2*firstBranch_-1) : -(2*firstBranch_-1);
  if (0) {
    printf("branching %s on %d bounds %g %g / %g %g\n",
	   (way==-1) ? "down" :"up",iColumn,
	   down_[0],down_[1],up_[0],up_[1]);
    const double * lower = solver->getColLower();
    const double * upper = solver->getColUpper();
    for (int i=0;i<8;i++) 
      printf(" [%d (%g,%g)]",i,lower[i],upper[i]);
    printf("\n");
  }
  if (way<0) {
#ifdef OSI_DEBUG
  { double olb,oub ;
    olb = solver->getColLower()[iColumn] ;
    oub = solver->getColUpper()[iColumn] ;
    printf("branching down on var %d: [%g,%g] => [%g,%g]\n",
	   iColumn,olb,oub,down_[0],down_[1]) ; }
#endif
    solver->setColLower(iColumn,down_[0]);
    solver->setColUpper(iColumn,down_[1]);
  } else {
#ifdef OSI_DEBUG
  { double olb,oub ;
    olb = solver->getColLower()[iColumn] ;
    oub = solver->getColUpper()[iColumn] ;
    printf("branching up on var %d: [%g,%g] => [%g,%g]\n",
	   iColumn,olb,oub,up_[0],up_[1]) ; }
#endif
    solver->setColLower(iColumn,up_[0]);
    solver->setColUpper(iColumn,up_[1]);
  }
  double nlb = solver->getColLower()[iColumn];
  if (nlb<olb) {
#ifndef NDEBUG
    printf("bad lb change for column %d from %g to %g\n",iColumn,olb,nlb);
#endif
    solver->setColLower(iColumn,olb);
  }
  double nub = solver->getColUpper()[iColumn];
  if (nub>oub) {
#ifndef NDEBUG
    printf("bad ub change for column %d from %g to %g\n",iColumn,oub,nub);
#endif
    solver->setColUpper(iColumn,oub);
  }
#ifndef NDEBUG
  if (nlb<olb+1.0e-8&&nub>oub-1.0e-8)
    printf("bad null change for column %d - bounds %g,%g\n",iColumn,olb,oub);
#endif
  branchIndex_++;
  return 0.0;
}
// Print what would happen  
void
OsiIntegerBranchingObject::print(const OsiSolverInterface * solver)
{
  const OsiSimpleInteger * obj =
    dynamic_cast <const OsiSimpleInteger *>(originalObject_) ;
  assert (obj);
  int iColumn = obj->columnNumber();
  int way = (!branchIndex_) ? (2*firstBranch_-1) : -(2*firstBranch_-1);
  if (way<0) {
  { double olb,oub ;
    olb = solver->getColLower()[iColumn] ;
    oub = solver->getColUpper()[iColumn] ;
    printf("OsiInteger would branch down on var %d : [%g,%g] => [%g,%g]\n",
	   iColumn,olb,oub,down_[0],down_[1]) ; }
  } else {
  { double olb,oub ;
    olb = solver->getColLower()[iColumn] ;
    oub = solver->getColUpper()[iColumn] ;
    printf("OsiInteger would branch up on var %d : [%g,%g] => [%g,%g]\n",
	   iColumn,olb,oub,up_[0],up_[1]) ; }
  }
}
// Default Constructor 
OsiSOS::OsiSOS ()
  : OsiObject2(),
    members_(NULL),
    weights_(NULL),
    numberMembers_(0),
    sosType_(-1),
    integerValued_(false)
{
}

// Useful constructor (which are indices)
OsiSOS::OsiSOS (const OsiSolverInterface * ,  int numberMembers,
	   const int * which, const double * weights, int type)
  : OsiObject2(),
    numberMembers_(numberMembers),
    sosType_(type)
{
  integerValued_ = type==1; // not strictly true - should check problem
  if (numberMembers_) {
    members_ = new int[numberMembers_];
    weights_ = new double[numberMembers_];
    memcpy(members_,which,numberMembers_*sizeof(int));
    if (weights) {
      memcpy(weights_,weights,numberMembers_*sizeof(double));
    } else {
      for (int i=0;i<numberMembers_;i++)
        weights_[i]=i;
    }
    // sort so weights increasing
    CoinSort_2(weights_,weights_+numberMembers_,members_);
    double last = -COIN_DBL_MAX;
    int i;
    for (i=0;i<numberMembers_;i++) {
      double possible = CoinMax(last+1.0e-10,weights_[i]);
      weights_[i] = possible;
      last=possible;
    }
  } else {
    members_ = NULL;
    weights_ = NULL;
  }
  assert (sosType_>0&&sosType_<3);
}

// Copy constructor 
OsiSOS::OsiSOS ( const OsiSOS & rhs)
  :OsiObject2(rhs)
{
  numberMembers_ = rhs.numberMembers_;
  sosType_ = rhs.sosType_;
  integerValued_ = rhs.integerValued_;
  if (numberMembers_) {
    members_ = new int[numberMembers_];
    weights_ = new double[numberMembers_];
    memcpy(members_,rhs.members_,numberMembers_*sizeof(int));
    memcpy(weights_,rhs.weights_,numberMembers_*sizeof(double));
  } else {
    members_ = NULL;
    weights_ = NULL;
  }
}

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

// Assignment operator 
OsiSOS & 
OsiSOS::operator=( const OsiSOS& rhs)
{
  if (this!=&rhs) {
    OsiObject2::operator=(rhs);
    delete [] members_;
    delete [] weights_;
    numberMembers_ = rhs.numberMembers_;
    sosType_ = rhs.sosType_;
    integerValued_ = rhs.integerValued_;
    if (numberMembers_) {
      members_ = new int[numberMembers_];
      weights_ = new double[numberMembers_];
      memcpy(members_,rhs.members_,numberMembers_*sizeof(int));
      memcpy(weights_,rhs.weights_,numberMembers_*sizeof(double));
    } else {
      members_ = NULL;
      weights_ = NULL;
    }
  }
  return *this;
}

// Destructor 
OsiSOS::~OsiSOS ()
{
  delete [] members_;
  delete [] weights_;
}

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

  // check bounds etc
  double lastWeight=-1.0e100;
  for (j=0;j<numberMembers_;j++) {
    int iColumn = members_[j];
    if (lastWeight>=weights_[j]-1.0e-12)
      throw CoinError("Weights too close together in SOS","infeasibility","OsiSOS");
    lastWeight = weights_[j];
    if (upper[iColumn]) {
      double value = CoinMax(0.0,solution[iColumn]);
      if (value>integerTolerance) {
	// Possibly due to scaling a fixed variable might slip through
#ifdef COIN_DEVELOP
	if (value>upper[iColumn]+10.0*primalTolerance)
	  printf("** Variable %d (%d) has value %g and upper bound of %g\n",
		 iColumn,j,value,upper[iColumn]);
#endif
	if (value>upper[iColumn]) {
	  value=upper[iColumn];
	} 
	sum += value;
	weight += weights_[j]*value;
	if (firstNonZero<0)
	  firstNonZero=j;
	lastNonZero=j;
      }
      if (firstNonFixed<0)
	firstNonFixed=j;
      lastNonFixed=j;
    }
  }
  whichWay=1;
  whichWay_=1;
  if (lastNonZero-firstNonZero>=sosType_) {
    // find where to branch
    assert (sum>0.0);
    // probably best to use pseudo duals
    double value = lastNonZero-firstNonZero+1;
    value *= 0.5/static_cast<double> (numberMembers_);
    infeasibility_=value;
    otherInfeasibility_=1.0-value;
    if (info->defaultDual_>=0.0) {
      // Using pseudo shadow prices
      weight /= sum;
      int iWhere;
      for (iWhere=firstNonZero;iWhere<lastNonZero;iWhere++) 
	if (weight<weights_[iWhere+1])
	  break;
      assert (iWhere!=lastNonZero);
      /* Complicated - infeasibility is being used for branching so we
	 don't want estimate of satisfying set but of each way on branch.
	 So let us suppose that all on side being fixed to 0 goes to closest
      */
      int lastDown=iWhere;
      int firstUp=iWhere+1;
      if (sosType_==2) {
	// SOS 2 - choose nearest
	if (weight-weights_[iWhere]>=weights_[iWhere+1]-weight)
	  lastDown++;
	// But make sure OK
	if (lastDown==firstNonFixed) {
	  lastDown ++;
	} else if (lastDown==lastNonFixed) {
	  lastDown --;
	} 
	firstUp=lastDown;
      }
      // Now get current contribution and compute weight for end points
      double weightDown = 0.0;
      double weightUp = 0.0;
      const double * element = info->elementByColumn_;
      const int * row = info->row_;
      const CoinBigIndex * columnStart = info->columnStart_;
      const int * columnLength = info->columnLength_;
      double direction = info->direction_;
      const double * objective = info->objective_;
      // Compute where we would move to
      double objValue=0.0;
      double * useful = info->usefulRegion_;
      int * index = info->indexRegion_;
      int n=0;
      for (j=firstNonZero;j<=lastNonZero;j++) {
	int iColumn = members_[j];
	double multiplier = solution[iColumn];
	if (j>=lastDown)
	  weightDown += multiplier;
	if (j<=firstUp)
	  weightUp += multiplier;
	if (multiplier>0.0) {
	  objValue += objective[iColumn]*multiplier;
	  CoinBigIndex start = columnStart[iColumn];
	  CoinBigIndex end = start + columnLength[iColumn];
	  for (CoinBigIndex j=start;j<end;j++) {
	    int iRow = row[j];
	    double value = element[j]*multiplier;
	    if (useful[iRow]) {
	      value += useful[iRow];
	      if (!value)
		value = 1.0e-100;
	    } else {
	      assert (value);
	      index[n++]=iRow;
	    }
	    useful[iRow] = value;
	  }
	}
      }
      if (sosType_==2) 
	assert (fabs(weightUp+weightDown-sum-solution[members_[lastDown]])<1.0e-4);
      int startX[2];
      int endX[2];
      startX[0]=firstNonZero;
      startX[1]=firstUp;
      endX[0]=lastDown;
      endX[1]=lastNonZero;
      double fakeSolution[2];
      int check[2];
      fakeSolution[0]=weightDown;
      check[0]=members_[lastDown];
      fakeSolution[1]=weightUp;
      check[1]=members_[firstUp];
      const double * pi = info->pi_;
      const double * activity = info->rowActivity_;
      const double * lower = info->rowLower_;
      const double * upper = info->rowUpper_;
      int numberRows = info->solver_->getNumRows();
      double * useful2 = useful+numberRows;
      int * index2 = index+numberRows;
      for (int i=0;i<2;i++) {
	double obj=0.0;
	int n2=0;
	for (j=startX[i];j<=endX[i];j++) {
	  int iColumn = members_[j];
	  double multiplier = solution[iColumn];
	  if (iColumn==check[i])
	    multiplier=fakeSolution[i];
	  if (multiplier>0.0) {
	    obj += objective[iColumn]*multiplier;
	    CoinBigIndex start = columnStart[iColumn];
	    CoinBigIndex end = start + columnLength[iColumn];
	    for (CoinBigIndex j=start;j<end;j++) {
	      int iRow = row[j];
	      double value = element[j]*multiplier;
	      if (useful2[iRow]) {
		value += useful2[iRow];
		if (!value)
		  value = 1.0e-100;
	      } else {
		assert (value);
		index2[n2++]=iRow;
	      }
	      useful2[iRow] = value;
	    }
	  }
	}
	// movement in objective
	obj = (obj-objValue) * direction;
	double estimate = (obj>0.0) ? obj : 0.0;
	for (j=0;j<n;j++) {
	  int iRow = index[j];
	  // movement
	  double movement = useful2[iRow]-useful[iRow];
	  useful[iRow]=0.0;
	  useful2[iRow]=0.0;
	  double valueP = pi[iRow]*direction;
	  if (lower[iRow]<-1.0e20) 
	    assert (valueP<=1.0e-4);
	  if (upper[iRow]>1.0e20) 
	    assert (valueP>=-1.0e-4);
	  double value2 = valueP*movement;
	  double thisEstimate = (value2>0.0) ? value2 : 0;
	  // if makes infeasible then make at least default
	  double newValue = activity[iRow] + movement;
	  if (newValue>upper[iRow]+primalTolerance||newValue<lower[iRow]-primalTolerance)
	    thisEstimate = CoinMax(thisEstimate,info->defaultDual_);
	  estimate += thisEstimate;
	}
	for (j=0;j<n2;j++) {
	  int iRow = index2[j];
	  // movement
	  double movement = useful2[iRow]-useful[iRow];
	  useful[iRow]=0.0;
	  useful2[iRow]=0.0;
	  if (movement) {
	    double valueP = pi[iRow]*direction;
	    if (lower[iRow]<-1.0e20) 
	      assert (valueP<=1.0e-4);
	    if (upper[iRow]>1.0e20) 
	      assert (valueP>=-1.0e-4);
	    double value2 = valueP*movement;
	    double thisEstimate = (value2>0.0) ? value2 : 0;
	    // if makes infeasible then make at least default
	    double newValue = activity[iRow] + movement;
	    if (newValue>upper[iRow]+primalTolerance||newValue<lower[iRow]-primalTolerance)
	      thisEstimate = CoinMax(thisEstimate,info->defaultDual_);
	    estimate += thisEstimate;
	  }
	}
	// store in fakeSolution
	fakeSolution[i]=estimate;
      }
      double downEstimate = fakeSolution[0];
      double upEstimate = fakeSolution[1];
      if (downEstimate>=upEstimate) {
	infeasibility_ = CoinMax(1.0e-12,upEstimate);
	otherInfeasibility_ = CoinMax(1.0e-12,downEstimate);
	whichWay = 1;
      } else {
	infeasibility_ = CoinMax(1.0e-12,downEstimate);
	otherInfeasibility_ = CoinMax(1.0e-12,upEstimate);
	whichWay = 0;
      }
      whichWay_=static_cast<short>(whichWay);
      value=infeasibility_;
    }
    return value;
  } else {
    infeasibility_=0.0;
    otherInfeasibility_=1.0;
    return 0.0; // satisfied
  }
}

// This looks at solution and sets bounds to contain solution
double
OsiSOS::feasibleRegion(OsiSolverInterface * solver, const OsiBranchingInformation * info) const
{
  int j;
  int firstNonZero=-1;
  int lastNonZero = -1;
  const double * solution = info->solution_;
  //const double * lower = info->lower_;
  const double * upper = info->upper_;
  double sum =0.0;
  // Find largest one or pair
  double movement=0.0;
  if (sosType_==1) {
    for (j=0;j<numberMembers_;j++) {
      int iColumn = members_[j];
      double value = CoinMax(0.0,solution[iColumn]);
      if (value>sum&&upper[iColumn]) {
	firstNonZero=j;
	sum=value;
      }
    }
    lastNonZero=firstNonZero;
  } else {
    // type 2
    for (j=1;j<numberMembers_;j++) {
      int iColumn = members_[j];
      int jColumn = members_[j-1];
      double value1 = CoinMax(0.0,solution[iColumn]);
      double value0 = CoinMax(0.0,solution[jColumn]);
      double value = value0+value1;
      if (value>sum) {
	if (upper[iColumn]||upper[jColumn]) {
	  firstNonZero=upper[jColumn] ? j-1 : j;
	  lastNonZero=upper[iColumn] ? j : j-1;
	  sum=value;
	}
      }
    }
  }
  for (j=0;j<numberMembers_;j++) {
    if (j<firstNonZero||j>lastNonZero) {
      int iColumn = members_[j];
      double value = CoinMax(0.0,solution[iColumn]);
      movement += value;
      solver->setColUpper(iColumn,0.0);
    }
  }
  return movement;
}
// Redoes data when sequence numbers change
void 
OsiSOS::resetSequenceEtc(int numberColumns, const int * originalColumns)
{
  int n2=0;
  for (int j=0;j<numberMembers_;j++) {
    int iColumn = members_[j];
    int i;
    for (i=0;i<numberColumns;i++) {
      if (originalColumns[i]==iColumn)
        break;
    }
    if (i<numberColumns) {
      members_[n2]=i;
      weights_[n2++]=weights_[j];
    }
  }
  if (n2<numberMembers_) {
    printf("** SOS number of members reduced from %d to %d!\n",numberMembers_,n2);
    numberMembers_=n2;
  }
}
// Return "down" estimate
double 
OsiSOS::downEstimate() const
{
  if (whichWay_)
    return otherInfeasibility_;
  else
    return infeasibility_;
}
// Return "up" estimate
double 
OsiSOS::upEstimate() const
{
  if (!whichWay_)
    return otherInfeasibility_;
  else
    return infeasibility_;
}

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

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

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

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


// Destructor 
OsiSOSBranchingObject::~OsiSOSBranchingObject ()
{
}
double
OsiSOSBranchingObject::branch(OsiSolverInterface * solver)
{
  const OsiSOS * set =
    dynamic_cast <const OsiSOS *>(originalObject_) ;
  assert (set);
  int way = (!branchIndex_) ? (2*firstBranch_-1) : -(2*firstBranch_-1);
  branchIndex_++;
  int numberMembers = set->numberMembers();
  const int * which = set->members();
  const double * weights = set->weights();
  //const double * lower = solver->getColLower();
  //const double * upper = solver->getColUpper();
  // *** for way - up means fix all those in down section
  if (way<0) {
    int i;
    for ( i=0;i<numberMembers;i++) {
      if (weights[i] > value_)
	break;
    }
    assert (i<numberMembers);
    for (;i<numberMembers;i++) 
      solver->setColUpper(which[i],0.0);
  } else {
    int i;
    for ( i=0;i<numberMembers;i++) {
      if (weights[i] >= value_)
	break;
      else
	solver->setColUpper(which[i],0.0);
    }
    assert (i<numberMembers);
  }
  return 0.0;
}
// Print what would happen  
void
OsiSOSBranchingObject::print(const OsiSolverInterface * solver)
{
  const OsiSOS * set =
    dynamic_cast <const OsiSOS *>(originalObject_) ;
  assert (set);
  int way = (!branchIndex_) ? (2*firstBranch_-1) : -(2*firstBranch_-1);
  int numberMembers = set->numberMembers();
  const int * which = set->members();
  const double * weights = set->weights();
  //const double * lower = solver->getColLower();
  const double * upper = solver->getColUpper();
  int first=numberMembers;
  int last=-1;
  int numberFixed=0;
  int numberOther=0;
  int i;
  for ( i=0;i<numberMembers;i++) {
    double bound = upper[which[i]];
    if (bound) {
      first = CoinMin(first,i);
      last = CoinMax(last,i);
    }
  }
  // *** for way - up means fix all those in down section
  if (way<0) {
    printf("SOS Down");
    for ( i=0;i<numberMembers;i++) {
      double bound = upper[which[i]];
      if (weights[i] > value_)
	break;
      else if (bound)
	numberOther++;
    }
    assert (i<numberMembers);
    for (;i<numberMembers;i++) {
      double bound = upper[which[i]];
      if (bound)
	numberFixed++;
    }
  } else {
    printf("SOS Up");
    for ( i=0;i<numberMembers;i++) {
      double bound = upper[which[i]];
      if (weights[i] >= value_)
	break;
      else if (bound)
	numberFixed++;
    }
    assert (i<numberMembers);
    for (;i<numberMembers;i++) {
      double bound = upper[which[i]];
      if (bound)
	numberOther++;
    }
  }
  printf(" - at %g, free range %d (%g) => %d (%g), %d would be fixed, %d other way\n",
	 value_,which[first],weights[first],which[last],weights[last],numberFixed,numberOther);
}
/** Default Constructor

*/
OsiLotsize::OsiLotsize ()
  : OsiObject2(),
    columnNumber_(-1),
    rangeType_(0),
    numberRanges_(0),
    largestGap_(0),
    bound_(NULL),
    range_(0)
{
}

/** Useful constructor

  Loads actual upper & lower bounds for the specified variable.
*/
OsiLotsize::OsiLotsize (const OsiSolverInterface * , 
				    int iColumn, int numberPoints,
			const double * points, bool range)
  : OsiObject2()
{
  assert (numberPoints>0);
  columnNumber_ = iColumn ;
  // sort ranges
  int * sort = new int[numberPoints];
  double * weight = new double [numberPoints];
  int i;
  if (range) {
    rangeType_=2;
  } else {
    rangeType_=1;
  }
  for (i=0;i<numberPoints;i++) {
    sort[i]=i;
    weight[i]=points[i*rangeType_];
  }
  CoinSort_2(weight,weight+numberPoints,sort);
  numberRanges_=1;
  largestGap_=0;
  if (rangeType_==1) {
    bound_ = new double[numberPoints+1];
    bound_[0]=weight[0];
    for (i=1;i<numberPoints;i++) {
      if (weight[i]!=weight[i-1]) 
	bound_[numberRanges_++]=weight[i];
    }
    // and for safety
    bound_[numberRanges_]=bound_[numberRanges_-1];
    for (i=1;i<numberRanges_;i++) {
      largestGap_ = CoinMax(largestGap_,bound_[i]-bound_[i-1]);
    }
  } else {
    bound_ = new double[2*numberPoints+2];
    bound_[0]=points[sort[0]*2];
    bound_[1]=points[sort[0]*2+1];
    double hi=bound_[1];
    assert (hi>=bound_[0]);
    for (i=1;i<numberPoints;i++) {
      double thisLo =points[sort[i]*2];
      double thisHi =points[sort[i]*2+1];
      assert (thisHi>=thisLo);
      if (thisLo>hi) {
	bound_[2*numberRanges_]=thisLo;
	bound_[2*numberRanges_+1]=thisHi;
	numberRanges_++;
	hi=thisHi;
      } else {
	//overlap
	hi=CoinMax(hi,thisHi);
	bound_[2*numberRanges_-1]=hi;
      }
    }
    // and for safety
    bound_[2*numberRanges_]=bound_[2*numberRanges_-2];
    bound_[2*numberRanges_+1]=bound_[2*numberRanges_-1];
    for (i=1;i<numberRanges_;i++) {
      largestGap_ = CoinMax(largestGap_,bound_[2*i]-bound_[2*i-1]);
    }
  }
  delete [] sort;
  delete [] weight;
  range_=0;
}

// Copy constructor 
OsiLotsize::OsiLotsize ( const OsiLotsize & rhs)
  :OsiObject2(rhs)

{
  columnNumber_ = rhs.columnNumber_;
  rangeType_ = rhs.rangeType_;
  numberRanges_ = rhs.numberRanges_;
  range_ = rhs.range_;
  largestGap_ = rhs.largestGap_;
  if (numberRanges_) {
    assert (rangeType_>0&&rangeType_<3);
    bound_= new double [(numberRanges_+1)*rangeType_];
    memcpy(bound_,rhs.bound_,(numberRanges_+1)*rangeType_*sizeof(double));
  } else {
    bound_=NULL;
  }
}

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

// Assignment operator 
OsiLotsize & 
OsiLotsize::operator=( const OsiLotsize& rhs)
{
  if (this!=&rhs) {
    OsiObject2::operator=(rhs);
    columnNumber_ = rhs.columnNumber_;
    rangeType_ = rhs.rangeType_;
    numberRanges_ = rhs.numberRanges_;
    largestGap_ = rhs.largestGap_;
    delete [] bound_;
    range_ = rhs.range_;
    if (numberRanges_) {
      assert (rangeType_>0&&rangeType_<3);
      bound_= new double [(numberRanges_+1)*rangeType_];
      memcpy(bound_,rhs.bound_,(numberRanges_+1)*rangeType_*sizeof(double));
    } else {
      bound_=NULL;
    }
  }
  return *this;
}

// Destructor 
OsiLotsize::~OsiLotsize ()
{
  delete [] bound_;
}
/* Finds range of interest so value is feasible in range range_ or infeasible 
   between hi[range_] and lo[range_+1].  Returns true if feasible.
*/
bool 
OsiLotsize::findRange(double value, double integerTolerance) const
{
  assert (range_>=0&&range_<numberRanges_+1);
  int iLo;
  int iHi;
  double infeasibility=0.0;
  if (rangeType_==1) {
    if (value<bound_[range_]-integerTolerance) {
      iLo=0;
      iHi=range_-1;
    } else if (value<bound_[range_]+integerTolerance) {
      return true;
    } else if (value<bound_[range_+1]-integerTolerance) {
      return false;
    } else {
      iLo=range_+1;
      iHi=numberRanges_-1;
    }
    // check lo and hi
    bool found=false;
    if (value>bound_[iLo]-integerTolerance&&value<bound_[iLo+1]+integerTolerance) {
      range_=iLo;
      found=true;
    } else if (value>bound_[iHi]-integerTolerance&&value<bound_[iHi+1]+integerTolerance) {
      range_=iHi;
      found=true;
    } else {
      range_ = (iLo+iHi)>>1;
    }
    //points
    while (!found) {
      if (value<bound_[range_]) {
	if (value>=bound_[range_-1]) {
	  // found
	  range_--;
	  break;
	} else {
	  iHi = range_;
	}
      } else {
	if (value<bound_[range_+1]) {
	  // found
	  break;
	} else {
	  iLo = range_;
	}
      }
      range_ = (iLo+iHi)>>1;
    }
    if (value-bound_[range_]<=bound_[range_+1]-value) {
      infeasibility = value-bound_[range_];
    } else {
      infeasibility = bound_[range_+1]-value;
      if (infeasibility<integerTolerance)
	range_++;
    }
    return (infeasibility<integerTolerance);
  } else {
    // ranges
    if (value<bound_[2*range_]-integerTolerance) {
      iLo=0;
      iHi=range_-1;
    } else if (value<bound_[2*range_+1]+integerTolerance) {
      return true;
    } else if (value<bound_[2*range_+2]-integerTolerance) {
      return false;
    } else {
      iLo=range_+1;
      iHi=numberRanges_-1;
    }
    // check lo and hi
    bool found=false;
    if (value>bound_[2*iLo]-integerTolerance&&value<bound_[2*iLo+2]-integerTolerance) {
      range_=iLo;
      found=true;
    } else if (value>=bound_[2*iHi]-integerTolerance) {
      range_=iHi;
      found=true;
    } else {
      range_ = (iLo+iHi)>>1;
    }
    //points
    while (!found) {
      if (value<bound_[2*range_]) {
	if (value>=bound_[2*range_-2]) {
	  // found
	  range_--;
	  break;
	} else {
	  iHi = range_;
	}
      } else {
	if (value<bound_[2*range_+2]) {
	  // found
	  break;
	} else {
	  iLo = range_;
	}
      }
      range_ = (iLo+iHi)>>1;
    }
    if (value>=bound_[2*range_]-integerTolerance&&value<=bound_[2*range_+1]+integerTolerance)
      infeasibility=0.0;
    else if (value-bound_[2*range_+1]<bound_[2*range_+2]-value) {
      infeasibility = value-bound_[2*range_+1];
    } else {
      infeasibility = bound_[2*range_+2]-value;
    }
    return (infeasibility<integerTolerance);
  }
}
/* Returns floor and ceiling
 */
void 
OsiLotsize::floorCeiling(double & floorLotsize, double & ceilingLotsize, double value,
			 double tolerance) const
{
  bool feasible=findRange(value,tolerance);
  if (rangeType_==1) {
    floorLotsize=bound_[range_];
    ceilingLotsize=bound_[range_+1];
    // may be able to adjust
    if (feasible&&fabs(value-floorLotsize)>fabs(value-ceilingLotsize)) {
      floorLotsize=bound_[range_+1];
      ceilingLotsize=bound_[range_+2];
    }
  } else {
    // ranges
    assert (value>=bound_[2*range_+1]);
    floorLotsize=bound_[2*range_+1];
    ceilingLotsize=bound_[2*range_+2];
  }
}

// Infeasibility - large is 0.5
double 
OsiLotsize::infeasibility(const OsiBranchingInformation * info, int & preferredWay) const
{
  const double * solution = info->solution_;
  const double * lower = info->lower_;
  const double * upper = info->upper_;
  double value = solution[columnNumber_];
  value = CoinMax(value, lower[columnNumber_]);
  value = CoinMin(value, upper[columnNumber_]);
  double integerTolerance = info->integerTolerance_;
  /*printf("%d %g %g %g %g\n",columnNumber_,value,lower[columnNumber_],
    solution[columnNumber_],upper[columnNumber_]);*/
  assert (value>=bound_[0]-integerTolerance
          &&value<=bound_[rangeType_*numberRanges_-1]+integerTolerance);
  infeasibility_=0.0;
  bool feasible = findRange(value,integerTolerance);
  if (!feasible) {
    if (rangeType_==1) {
      if (value-bound_[range_]<bound_[range_+1]-value) {
	preferredWay=-1;
	infeasibility_ = value-bound_[range_];
	otherInfeasibility_ = bound_[range_+1] - value ;
      } else {
	preferredWay=1;
	infeasibility_ = bound_[range_+1]-value;
	otherInfeasibility_ = value-bound_[range_];
      }
    } else {
      // ranges
      if (value-bound_[2*range_+1]<bound_[2*range_+2]-value) {
	preferredWay=-1;
	infeasibility_ = value-bound_[2*range_+1];
	otherInfeasibility_ = bound_[2*range_+2]-value;
      } else {
	preferredWay=1;
	infeasibility_ = bound_[2*range_+2]-value;
	otherInfeasibility_ = value-bound_[2*range_+1];
      }
    }
  } else {
    // always satisfied
    preferredWay=-1;
    otherInfeasibility_ = 1.0;
  }
  if (infeasibility_<integerTolerance)
    infeasibility_=0.0;
  else
    infeasibility_ /= largestGap_;
  return infeasibility_;
}
/* Column number if single column object -1 otherwise,
   so returns >= 0
   Used by heuristics
*/
int 
OsiLotsize::columnNumber() const
{
  return columnNumber_;
}
/* Set bounds to contain the current solution.
   More precisely, for the variable associated with this object, take the
   value given in the current solution, force it within the current bounds
   if required, then set the bounds to fix the variable at the integer
   nearest the solution value.  Returns amount it had to move variable.
*/
double 
OsiLotsize::feasibleRegion(OsiSolverInterface * solver, const OsiBranchingInformation * info) const
{
  const double * lower = solver->getColLower();
  const double * upper = solver->getColUpper();
  const double * solution = info->solution_;
  double value = solution[columnNumber_];
  value = CoinMax(value, lower[columnNumber_]);
  value = CoinMin(value, upper[columnNumber_]);
  findRange(value,info->integerTolerance_);
  double nearest;
  if (rangeType_==1) {
    nearest = bound_[range_];
    solver->setColLower(columnNumber_,nearest);
    solver->setColUpper(columnNumber_,nearest);
  } else {
    // ranges
    solver->setColLower(columnNumber_,bound_[2*range_]);
    solver->setColUpper(columnNumber_,bound_[2*range_+1]);
    if (value>bound_[2*range_+1]) 
      nearest=bound_[2*range_+1];
    else if (value<bound_[2*range_]) 
      nearest = bound_[2*range_];
    else
      nearest = value;
  }
  // Scaling may have moved it a bit
  // Lotsizing variables could be a lot larger
#ifndef NDEBUG
  assert (fabs(value-nearest)<=(100.0+10.0*fabs(nearest))*info->integerTolerance_);
#endif
  return fabs(value-nearest);
}

// Creates a branching object
// Creates a branching object
OsiBranchingObject * 
OsiLotsize::createBranch(OsiSolverInterface * solver, const OsiBranchingInformation * info, int way) const 
{
  const double * solution = info->solution_;
  const double * lower = solver->getColLower();
  const double * upper = solver->getColUpper();
  double value = solution[columnNumber_];
  value = CoinMax(value, lower[columnNumber_]);
  value = CoinMin(value, upper[columnNumber_]);
  assert (!findRange(value,info->integerTolerance_));
  return new OsiLotsizeBranchingObject(solver,this,way,
					     value);
}

  
/*
  Bounds may be tightened, so it may be good to be able to refresh the local
  copy of the original bounds.
 */
void 
OsiLotsize::resetBounds(const OsiSolverInterface * )
{
}
// Return "down" estimate
double 
OsiLotsize::downEstimate() const
{
  if (whichWay_)
    return otherInfeasibility_;
  else
    return infeasibility_;
}
// Return "up" estimate
double 
OsiLotsize::upEstimate() const
{
  if (!whichWay_)
    return otherInfeasibility_;
  else
    return infeasibility_;
}
// Redoes data when sequence numbers change
void 
OsiLotsize::resetSequenceEtc(int numberColumns, const int * originalColumns)
{
  int i;
  for (i=0;i<numberColumns;i++) {
    if (originalColumns[i]==columnNumber_)
      break;
  }
  if (i<numberColumns)
    columnNumber_=i;
  else
    abort(); // should never happen
}


// Default Constructor 
OsiLotsizeBranchingObject::OsiLotsizeBranchingObject()
  :OsiTwoWayBranchingObject()
{
  down_[0] = 0.0;
  down_[1] = 0.0;
  up_[0] = 0.0;
  up_[1] = 0.0;
}

// Useful constructor
OsiLotsizeBranchingObject::OsiLotsizeBranchingObject (OsiSolverInterface * solver, 
						      const OsiLotsize * originalObject, 
						      int way , double value)
  :OsiTwoWayBranchingObject(solver,originalObject,way,value)
{
  int iColumn = originalObject->columnNumber();
  down_[0] = solver->getColLower()[iColumn];
  double integerTolerance = solver->getIntegerTolerance();
  originalObject->floorCeiling(down_[1],up_[0],value,integerTolerance);
  up_[1] = solver->getColUpper()[iColumn];
}

// Copy constructor 
OsiLotsizeBranchingObject::OsiLotsizeBranchingObject ( const OsiLotsizeBranchingObject & rhs) :OsiTwoWayBranchingObject(rhs)
{
  down_[0] = rhs.down_[0];
  down_[1] = rhs.down_[1];
  up_[0] = rhs.up_[0];
  up_[1] = rhs.up_[1];
}

// Assignment operator 
OsiLotsizeBranchingObject & 
OsiLotsizeBranchingObject::operator=( const OsiLotsizeBranchingObject& rhs)
{
  if (this != &rhs) {
    OsiTwoWayBranchingObject::operator=(rhs);
    down_[0] = rhs.down_[0];
    down_[1] = rhs.down_[1];
    up_[0] = rhs.up_[0];
    up_[1] = rhs.up_[1];
  }
  return *this;
}
OsiBranchingObject * 
OsiLotsizeBranchingObject::clone() const
{ 
  return (new OsiLotsizeBranchingObject(*this));
}


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

/*
  Perform a branch by adjusting the bounds of the specified variable. Note
  that each arm of the branch advances the object to the next arm by
  advancing the value of way_.

  Providing new values for the variable's lower and upper bounds for each
  branching direction gives a little bit of additional flexibility and will
  be easily extensible to multi-way branching.
*/
double
OsiLotsizeBranchingObject::branch(OsiSolverInterface * solver)
{
  const OsiLotsize * obj =
    dynamic_cast <const OsiLotsize *>(originalObject_) ;
  assert (obj);
  int iColumn = obj->columnNumber();
  int way = (!branchIndex_) ? (2*firstBranch_-1) : -(2*firstBranch_-1);
  if (way<0) {
#ifdef OSI_DEBUG
  { double olb,oub ;
    olb = solver->getColLower()[iColumn] ;
    oub = solver->getColUpper()[iColumn] ;
    printf("branching down on var %d: [%g,%g] => [%g,%g]\n",
	   iColumn,olb,oub,down_[0],down_[1]) ; }
#endif
    solver->setColLower(iColumn,down_[0]);
    solver->setColUpper(iColumn,down_[1]);
  } else {
#ifdef OSI_DEBUG
  { double olb,oub ;
    olb = solver->getColLower()[iColumn] ;
    oub = solver->getColUpper()[iColumn] ;
    printf("branching up on var %d: [%g,%g] => [%g,%g]\n",
	   iColumn,olb,oub,up_[0],up_[1]) ; }
#endif
    solver->setColLower(iColumn,up_[0]);
    solver->setColUpper(iColumn,up_[1]);
  }
  branchIndex_++;
  return 0.0;
}
// Print
void
OsiLotsizeBranchingObject::print(const OsiSolverInterface * solver)
{
  const OsiLotsize * obj =
    dynamic_cast <const OsiLotsize *>(originalObject_) ;
  assert (obj);
  int iColumn = obj->columnNumber();
  int way = (!branchIndex_) ? (2*firstBranch_-1) : -(2*firstBranch_-1);
  if (way<0) {
  { double olb,oub ;
    olb = solver->getColLower()[iColumn] ;
    oub = solver->getColUpper()[iColumn] ;
    printf("branching down on var %d: [%g,%g] => [%g,%g]\n",
	   iColumn,olb,oub,down_[0],down_[1]) ; }
  } else {
  { double olb,oub ;
    olb = solver->getColLower()[iColumn] ;
    oub = solver->getColUpper()[iColumn] ;
    printf("branching up on var %d: [%g,%g] => [%g,%g]\n",
	   iColumn,olb,oub,up_[0],up_[1]) ; }
  }
}