// Copyright (C) 2002, International Business Machines
// Corporation and others.  All Rights Reserved.

#include "CoinPragma.hpp"
#include <iostream>
#include <cassert>

#include "CoinIndexedVector.hpp"

#include "ClpSimplex.hpp"
#include "CoinHelperFunctions.hpp"
#include "ClpNonLinearCost.hpp"
//#############################################################################
// Constructors / Destructor / Assignment
//#############################################################################

//-------------------------------------------------------------------
// Default Constructor 
//-------------------------------------------------------------------
ClpNonLinearCost::ClpNonLinearCost () :
  changeCost_(0.0),
  feasibleCost_(0.0),
  infeasibilityWeight_(-1.0),
  largestInfeasibility_(0.0),
  sumInfeasibilities_(0.0),
  averageTheta_(0.0),
  numberRows_(0),
  numberColumns_(0),
  start_(NULL),
  whichRange_(NULL),
  offset_(NULL),
  lower_(NULL),
  cost_(NULL),
  model_(NULL),
  infeasible_(NULL),
  numberInfeasibilities_(-1),
  status_(NULL),
  bound_(NULL),
  cost2_(NULL),
  method_(1),
  convex_(true),
  bothWays_(false)
{

}
//#define VALIDATE
#ifdef VALIDATE
static double * saveLowerV=NULL;
static double * saveUpperV=NULL;
#ifdef NDEBUG
Validate sgould not be set if no debug
#endif
#endif
/* Constructor from simplex.
   This will just set up wasteful arrays for linear, but
   later may do dual analysis and even finding duplicate columns 
*/
ClpNonLinearCost::ClpNonLinearCost ( ClpSimplex * model,int method)
{
  method=2;
  model_ = model;
  numberRows_ = model_->numberRows();
  //if (numberRows_==402) {
  //model_->setLogLevel(63);
  //model_->setMaximumIterations(30000);
  //}
  numberColumns_ = model_->numberColumns();
  // If gub then we need this extra
  int numberExtra = model_->numberExtraRows();
  if (numberExtra)
    method=1;
  int numberTotal1 = numberRows_+numberColumns_;
  int numberTotal = numberTotal1+numberExtra;
  convex_ = true;
  bothWays_ = false;
  method_=method;
  numberInfeasibilities_=0;
  changeCost_=0.0;
  feasibleCost_=0.0;
  infeasibilityWeight_ = -1.0;
  double * cost = model_->costRegion();
  // check if all 0
  int iSequence;
  bool allZero=true;
  for (iSequence=0;iSequence<numberTotal1;iSequence++) {
    if (cost[iSequence]) {
      allZero=false;
      break;
    }
  }
  if (allZero)
    model_->setInfeasibilityCost(1.0);
  double infeasibilityCost = model_->infeasibilityCost();
  sumInfeasibilities_=0.0;
  averageTheta_=0.0;
  largestInfeasibility_=0.0;
  // All arrays NULL to start
  status_ = NULL;
  bound_ = NULL;
  cost2_ = NULL;
  start_ = NULL;
  whichRange_ = NULL;
  offset_ = NULL;
  lower_ = NULL;
  cost_= NULL;
  infeasible_=NULL;

  double * upper = model_->upperRegion();
  double * lower = model_->lowerRegion();

  // See how we are storing things
  bool always4 = (model_->clpMatrix()->
                  generalExpanded(model_,10,iSequence)!=0);
  if (always4)
    method_=1;
  if (CLP_METHOD1) {
    start_ = new int [numberTotal+1];
    whichRange_ = new int [numberTotal];
    offset_ = new int [numberTotal];
    memset(offset_,0,numberTotal*sizeof(int));
    
    
    // First see how much space we need
    int put=0;
    
    // For quadratic we need -inf,0,0,+inf
    for (iSequence=0;iSequence<numberTotal1;iSequence++) {
      if (!always4) {
        if (lower[iSequence]>-COIN_DBL_MAX)
          put++;
        if (upper[iSequence]<COIN_DBL_MAX)
          put++;
        put += 2;
      } else {
        put += 4;
      }
    }
    
    // and for extra
    put += 4*numberExtra;
#ifndef NDEBUG
    int kPut=put;
#endif
    lower_ = new double [put];
    cost_ = new double [put];
    infeasible_ = new unsigned int[(put+31)>>5];
    memset(infeasible_,0,((put+31)>>5)*sizeof(unsigned int));
    
    put=0;
    
    start_[0]=0;
    
    for (iSequence=0;iSequence<numberTotal1;iSequence++) {
      if (!always4) {
        if (lower[iSequence]>-COIN_DBL_MAX) {
          lower_[put] = -COIN_DBL_MAX;
          setInfeasible(put,true);
          cost_[put++] = cost[iSequence]-infeasibilityCost;
        }
        whichRange_[iSequence]=put;
        lower_[put] = lower[iSequence];
        cost_[put++] = cost[iSequence];
        lower_[put] = upper[iSequence];
        cost_[put++] = cost[iSequence]+infeasibilityCost;
        if (upper[iSequence]<COIN_DBL_MAX) {
          lower_[put] = COIN_DBL_MAX;
          setInfeasible(put-1,true);
          cost_[put++] = 1.0e50;
        }
      } else {
        lower_[put] = -COIN_DBL_MAX;
        setInfeasible(put,true);
        cost_[put++] = cost[iSequence]-infeasibilityCost;
        whichRange_[iSequence]=put;
        lower_[put] = lower[iSequence];
        cost_[put++] = cost[iSequence];
        lower_[put] = upper[iSequence];
        cost_[put++] = cost[iSequence]+infeasibilityCost;
        lower_[put] = COIN_DBL_MAX;
        setInfeasible(put-1,true);
        cost_[put++] = 1.0e50;
      }
      start_[iSequence+1]=put;
    }
    for (;iSequence<numberTotal;iSequence++) {
      
      lower_[put] = -COIN_DBL_MAX;
      setInfeasible(put,true);
      put++;
      whichRange_[iSequence]=put;
      lower_[put] = 0.0;
      cost_[put++] = 0.0;
      lower_[put] = 0.0;
      cost_[put++] = 0.0;
      lower_[put] = COIN_DBL_MAX;
      setInfeasible(put-1,true);
      cost_[put++] = 1.0e50;
      start_[iSequence+1]=put;
    }
    assert (put<=kPut);
  }
#ifdef FAST_CLPNON
  // See how we are storing things
  CoinAssert (model_->clpMatrix()->
                  generalExpanded(model_,10,iSequence)==0);
#endif
  if (CLP_METHOD2) {
    assert (!numberExtra);
    bound_ = new double[numberTotal];
    cost2_ = new double[numberTotal];
    status_ = new unsigned char[numberTotal];
#ifdef VALIDATE
    delete [] saveLowerV;
    saveLowerV = CoinCopyOfArray(model_->lowerRegion(),numberTotal);
    delete [] saveUpperV;
    saveUpperV = CoinCopyOfArray(model_->upperRegion(),numberTotal);
#endif
    for (iSequence=0;iSequence<numberTotal;iSequence++) {
      bound_[iSequence]=0.0;
      cost2_[iSequence]=cost[iSequence];
      setInitialStatus(status_[iSequence]);
    }
  }
}
// Refreshes costs always makes row costs zero
void 
ClpNonLinearCost::refreshCosts(const double * columnCosts)
{
  double * cost = model_->costRegion();
  // zero row costs
  memset(cost+numberColumns_,0,numberRows_*sizeof(double));
  // copy column costs
  memcpy(cost,columnCosts,numberColumns_*sizeof(double));
  if ((method_&1)!=0) {
    for (int iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
      int start = start_[iSequence];
      int end = start_[iSequence+1]-1;
      double thisFeasibleCost=cost[iSequence];
      if (infeasible(start)) {
        cost_[start] = thisFeasibleCost-infeasibilityWeight_;
        cost_[start+1] = thisFeasibleCost;
      } else {
        cost_[start] = thisFeasibleCost;
      }
      if (infeasible(end-1)) {
        cost_[end-1] = thisFeasibleCost+infeasibilityWeight_;
      }
    }
  }
  if (CLP_METHOD2) {
    for (int iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
      cost2_[iSequence]=cost[iSequence];
    }
  }
}
ClpNonLinearCost::ClpNonLinearCost(ClpSimplex * model,const int * starts,
		   const double * lowerNon, const double * costNon)
{
#ifndef FAST_CLPNON
  // what about scaling? - only try without it initially
  assert(!model->scalingFlag());
  model_ = model;
  numberRows_ = model_->numberRows();
  numberColumns_ = model_->numberColumns();
  int numberTotal = numberRows_+numberColumns_;
  convex_ = true;
  bothWays_ = true;
  start_ = new int [numberTotal+1];
  whichRange_ = new int [numberTotal];
  offset_ = new int [numberTotal];
  memset(offset_,0,numberTotal*sizeof(int));
  
  double whichWay = model_->optimizationDirection();
  printf("Direction %g\n",whichWay);

  numberInfeasibilities_=0;
  changeCost_=0.0;
  feasibleCost_=0.0;
  double infeasibilityCost = model_->infeasibilityCost();
  infeasibilityWeight_ = infeasibilityCost;;
  largestInfeasibility_=0.0;
  sumInfeasibilities_=0.0;

  int iSequence;
  assert (!model_->rowObjective());
  double * cost = model_->objective();

  // First see how much space we need 
  // Also set up feasible bounds
  int put=starts[numberColumns_];

  double * columnUpper = model_->columnUpper();
  double * columnLower = model_->columnLower();
  for (iSequence=0;iSequence<numberColumns_;iSequence++) {
    if (columnLower[iSequence]>-1.0e20)
      put++;
    if (columnUpper[iSequence]<1.0e20)
      put++;
  }

  double * rowUpper = model_->rowUpper();
  double * rowLower = model_->rowLower();
  for (iSequence=0;iSequence<numberRows_;iSequence++) {
    if (rowLower[iSequence]>-1.0e20)
      put++;
    if (rowUpper[iSequence]<1.0e20)
      put++;
    put +=2;
  }
  lower_ = new double [put];
  cost_ = new double [put];
  infeasible_ = new unsigned int[(put+31)>>5];
  memset(infeasible_,0,((put+31)>>5)*sizeof(unsigned int));

  // now fill in 
  put=0;

  start_[0]=0;
  for (iSequence=0;iSequence<numberTotal;iSequence++) {
    lower_[put] = -COIN_DBL_MAX;
    whichRange_[iSequence]=put+1;
    double thisCost;
    double lowerValue;
    double upperValue;
    if (iSequence>=numberColumns_) {
      // rows
      lowerValue = rowLower[iSequence-numberColumns_];
      upperValue = rowUpper[iSequence-numberColumns_];
      if (lowerValue>-1.0e30) {
	setInfeasible(put,true);
	cost_[put++] = -infeasibilityCost;
	lower_[put] = lowerValue;
      }
      cost_[put++] = 0.0;
      thisCost = 0.0;
    } else {
      // columns - move costs and see if convex
      lowerValue = columnLower[iSequence];
      upperValue = columnUpper[iSequence];
      if (lowerValue>-1.0e30) {
	setInfeasible(put,true);
	cost_[put++] = whichWay*cost[iSequence]-infeasibilityCost;
	lower_[put] = lowerValue;
      }
      int iIndex = starts[iSequence];
      int end = starts[iSequence+1];
      assert (fabs(columnLower[iSequence]-lowerNon[iIndex])<1.0e-8);
      thisCost = -COIN_DBL_MAX;
      for (;iIndex<end;iIndex++) {
	if (lowerNon[iIndex]<columnUpper[iSequence]-1.0e-8) {
	  lower_[put] = lowerNon[iIndex];
	  cost_[put++] = whichWay*costNon[iIndex];
	  // check convexity
	  if (whichWay*costNon[iIndex]<thisCost-1.0e-12)
	    convex_ = false;
	  thisCost = whichWay*costNon[iIndex];
	} else {
	  break;
	}
      }
    }
    lower_[put] = upperValue;
    setInfeasible(put,true);
    cost_[put++] = thisCost+infeasibilityCost;
    if (upperValue<1.0e20) {
      lower_[put] = COIN_DBL_MAX;
      cost_[put++] = 1.0e50;
    }
    int iFirst = start_[iSequence];
    if (lower_[iFirst] != -COIN_DBL_MAX) {
      setInfeasible(iFirst,true);
      whichRange_[iSequence]=iFirst+1;
    } else {
      whichRange_[iSequence]=iFirst;
    }
    start_[iSequence+1]=put;
  }
  // can't handle non-convex at present
  assert(convex_);
  status_ = NULL;
  bound_ = NULL;
  cost2_ = NULL;
  method_ = 1;
#else
  printf("recompile ClpNonLinearCost without FAST_CLPNON\n");
  abort();
#endif
}
//-------------------------------------------------------------------
// Copy constructor 
//-------------------------------------------------------------------
ClpNonLinearCost::ClpNonLinearCost (const ClpNonLinearCost & rhs) :
  changeCost_(0.0),
  feasibleCost_(0.0),
  infeasibilityWeight_(-1.0),
  largestInfeasibility_(0.0),
  sumInfeasibilities_(0.0),
  averageTheta_(0.0),
  numberRows_(rhs.numberRows_),
  numberColumns_(rhs.numberColumns_),
  start_(NULL),
  whichRange_(NULL),
  offset_(NULL),
  lower_(NULL),
  cost_(NULL),
  model_(NULL),
  infeasible_(NULL),
  numberInfeasibilities_(-1),
  status_(NULL),
  bound_(NULL),
  cost2_(NULL),
  method_(rhs.method_),
  convex_(true),
  bothWays_(rhs.bothWays_)
{  
  if (numberRows_) {
    int numberTotal = numberRows_+numberColumns_;
    model_ = rhs.model_;
    numberInfeasibilities_=rhs.numberInfeasibilities_;
    changeCost_ = rhs.changeCost_;
    feasibleCost_ = rhs.feasibleCost_;
    infeasibilityWeight_ = rhs.infeasibilityWeight_;
    largestInfeasibility_ = rhs.largestInfeasibility_;
    sumInfeasibilities_ = rhs.sumInfeasibilities_;
    averageTheta_ = rhs.averageTheta_;
    convex_ = rhs.convex_;
    if (CLP_METHOD1) {
      start_ = new int [numberTotal+1];
      memcpy(start_,rhs.start_,(numberTotal+1)*sizeof(int));
      whichRange_ = new int [numberTotal];
      memcpy(whichRange_,rhs.whichRange_,numberTotal*sizeof(int));
      offset_ = new int [numberTotal];
      memcpy(offset_,rhs.offset_,numberTotal*sizeof(int));
      int numberEntries = start_[numberTotal];
      lower_ = new double [numberEntries];
      memcpy(lower_,rhs.lower_,numberEntries*sizeof(double));
      cost_ = new double [numberEntries];
      memcpy(cost_,rhs.cost_,numberEntries*sizeof(double));
      infeasible_ = new unsigned int[(numberEntries+31)>>5];
      memcpy(infeasible_,rhs.infeasible_,
             ((numberEntries+31)>>5)*sizeof(unsigned int));
    }
    if (CLP_METHOD2) {
      bound_ = CoinCopyOfArray(rhs.bound_,numberTotal);
      cost2_ = CoinCopyOfArray(rhs.cost2_,numberTotal);
      status_ = CoinCopyOfArray(rhs.status_,numberTotal);
    }
  }
}

//-------------------------------------------------------------------
// Destructor 
//-------------------------------------------------------------------
ClpNonLinearCost::~ClpNonLinearCost ()
{
  delete [] start_;
  delete [] whichRange_;
  delete [] offset_;
  delete [] lower_;
  delete [] cost_;
  delete [] infeasible_;
  delete [] status_;
  delete [] bound_;
  delete [] cost2_;
}

//----------------------------------------------------------------
// Assignment operator 
//-------------------------------------------------------------------
ClpNonLinearCost &
ClpNonLinearCost::operator=(const ClpNonLinearCost& rhs)
{
  if (this != &rhs) {
    numberRows_ = rhs.numberRows_;
    numberColumns_ = rhs.numberColumns_;
    delete [] start_;
    delete [] whichRange_;
    delete [] offset_;
    delete [] lower_;
    delete [] cost_;
    delete [] infeasible_;
    delete [] status_;
    delete [] bound_;
    delete [] cost2_;
    start_ = NULL;
    whichRange_ = NULL;
    lower_ = NULL;
    cost_= NULL;
    infeasible_=NULL;
    status_ = NULL;
    bound_ = NULL;
    cost2_ = NULL;
    method_=rhs.method_;
    if (numberRows_) {
      int numberTotal = numberRows_+numberColumns_;
      if (CLP_METHOD1) {
        start_ = new int [numberTotal+1];
        memcpy(start_,rhs.start_,(numberTotal+1)*sizeof(int));
        whichRange_ = new int [numberTotal];
        memcpy(whichRange_,rhs.whichRange_,numberTotal*sizeof(int));
        offset_ = new int [numberTotal];
        memcpy(offset_,rhs.offset_,numberTotal*sizeof(int));
        int numberEntries = start_[numberTotal];
        lower_ = new double [numberEntries];
        memcpy(lower_,rhs.lower_,numberEntries*sizeof(double));
        cost_ = new double [numberEntries];
        memcpy(cost_,rhs.cost_,numberEntries*sizeof(double));
        infeasible_ = new unsigned int[(numberEntries+31)>>5];
        memcpy(infeasible_,rhs.infeasible_,
               ((numberEntries+31)>>5)*sizeof(unsigned int));
      }
      if (CLP_METHOD2) {
        bound_ = CoinCopyOfArray(rhs.bound_,numberTotal);
        cost2_ = CoinCopyOfArray(rhs.cost2_,numberTotal);
        status_ = CoinCopyOfArray(rhs.status_,numberTotal);
      }
    }      
    model_ = rhs.model_;
    numberInfeasibilities_=rhs.numberInfeasibilities_;
    changeCost_ = rhs.changeCost_;
    feasibleCost_ = rhs.feasibleCost_;
    infeasibilityWeight_ = rhs.infeasibilityWeight_;
    largestInfeasibility_ = rhs.largestInfeasibility_;
    sumInfeasibilities_ = rhs.sumInfeasibilities_;
    averageTheta_ = rhs.averageTheta_;
    convex_ = rhs.convex_;
    bothWays_ = rhs.bothWays_;
  }
  return *this;
}
// Changes infeasible costs and computes number and cost of infeas
// We will need to re-think objective offsets later
// We will also need a 2 bit per variable array for some
// purpose which will come to me later
void 
ClpNonLinearCost::checkInfeasibilities(double oldTolerance)
{
  numberInfeasibilities_=0;
  double infeasibilityCost = model_->infeasibilityCost();
  changeCost_=0.0;
  largestInfeasibility_ = 0.0;
  sumInfeasibilities_ = 0.0;
  double primalTolerance = model_->currentPrimalTolerance();
  int iSequence;
  double * solution = model_->solutionRegion();
  double * upper = model_->upperRegion();
  double * lower = model_->lowerRegion();
  double * cost = model_->costRegion();
  bool toNearest = oldTolerance<=0.0;
  feasibleCost_=0.0;
  //bool checkCosts = (infeasibilityWeight_ != infeasibilityCost);
  infeasibilityWeight_ = infeasibilityCost;
  int numberTotal = numberColumns_+numberRows_;
  //#define NONLIN_DEBUG
#ifdef NONLIN_DEBUG
  double * saveSolution=NULL;
  double * saveLower=NULL;
  double * saveUpper=NULL;
  unsigned char * saveStatus=NULL;
  if (method_==3) {
    // Save solution as we will be checking
    saveSolution = CoinCopyOfArray(solution,numberTotal);
    saveLower = CoinCopyOfArray(lower,numberTotal);
    saveUpper = CoinCopyOfArray(upper,numberTotal);
    saveStatus = CoinCopyOfArray(model_->statusArray(),numberTotal);
  }
#else
  assert (method_!=3);
#endif
  if (CLP_METHOD1) {
    // nonbasic should be at a valid bound
    for (iSequence=0;iSequence<numberTotal;iSequence++) {
      double lowerValue;
      double upperValue;
      double value=solution[iSequence];
      int iRange;
      // get correct place
      int start = start_[iSequence];
      int end = start_[iSequence+1]-1;
      // correct costs for this infeasibility weight
      // If free then true cost will be first
      double thisFeasibleCost=cost_[start];
      if (infeasible(start)) {
        thisFeasibleCost = cost_[start+1];
        cost_[start] = thisFeasibleCost-infeasibilityCost;
      }
      if (infeasible(end-1)) {
        thisFeasibleCost = cost_[end-2];
        cost_[end-1] = thisFeasibleCost+infeasibilityCost;
      }
      for (iRange=start; iRange<end;iRange++) {
        if (value<lower_[iRange+1]+primalTolerance) {
          // put in better range if infeasible
          if (value>=lower_[iRange+1]-primalTolerance&&infeasible(iRange)&&iRange==start) 
            iRange++;
          whichRange_[iSequence]=iRange;
          break;
        }
      }
      assert(iRange<end);
      lowerValue = lower_[iRange];
      upperValue = lower_[iRange+1];
      ClpSimplex::Status status = model_->getStatus(iSequence);
      if (upperValue==lowerValue && status!=ClpSimplex::isFixed) {
        if (status != ClpSimplex::basic) {
          model_->setStatus(iSequence,ClpSimplex::isFixed);
          status = ClpSimplex::isFixed;
        }
      }
      switch(status) {
        
      case ClpSimplex::basic:
      case ClpSimplex::superBasic:
        // iRange is in correct place
        // slot in here
        if (infeasible(iRange)) {
          if (lower_[iRange]<-1.0e50) {
            //cost_[iRange] = cost_[iRange+1]-infeasibilityCost;
            // possibly below
            lowerValue = lower_[iRange+1];
            if (value-lowerValue<-primalTolerance) {
              value = lowerValue-value-primalTolerance;
#ifndef NDEBUG
              if(value>1.0e15)
                printf("nonlincostb %d %g %g %g\n",
                       iSequence,lowerValue,solution[iSequence],lower_[iRange+2]);
#endif
              sumInfeasibilities_ += value;
              largestInfeasibility_ = CoinMax(largestInfeasibility_,value);
              changeCost_ -= lowerValue*
                (cost_[iRange]-cost[iSequence]);
              numberInfeasibilities_++;
            }
          } else {
            //cost_[iRange] = cost_[iRange-1]+infeasibilityCost;
            // possibly above
            upperValue = lower_[iRange];
            if (value-upperValue>primalTolerance) {
              value = value-upperValue-primalTolerance;
#ifndef NDEBUG
              if(value>1.0e15)
                printf("nonlincostu %d %g %g %g\n",
                       iSequence,lower_[iRange-1],solution[iSequence],upperValue);
#endif
              sumInfeasibilities_ += value;
              largestInfeasibility_ = CoinMax(largestInfeasibility_,value);
              changeCost_ -= upperValue*
                (cost_[iRange]-cost[iSequence]);
              numberInfeasibilities_++;
            }
          }
        }
        //lower[iSequence] = lower_[iRange];
        //upper[iSequence] = lower_[iRange+1];
        //cost[iSequence] = cost_[iRange];
        break;
      case ClpSimplex::isFree:
        //if (toNearest)
        //solution[iSequence] = 0.0;
        break;
      case ClpSimplex::atUpperBound:
        if (!toNearest) {
          // With increasing tolerances - we may be at wrong place
          if (fabs(value-upperValue)>oldTolerance*1.0001) {
            if (fabs(value-lowerValue)<=oldTolerance*1.0001) {
              if  (fabs(value-lowerValue)>primalTolerance)
                solution[iSequence]=lowerValue;
              model_->setStatus(iSequence,ClpSimplex::atLowerBound);
            } else {
              model_->setStatus(iSequence,ClpSimplex::superBasic);
            }
          } else if  (fabs(value-upperValue)>primalTolerance) {
            solution[iSequence]=upperValue;
          }
        } else {
          // Set to nearest and make at upper bound
          int kRange;
          iRange=-1;
          double nearest = COIN_DBL_MAX;
          for (kRange=start; kRange<end;kRange++) {
            if (fabs(lower_[kRange]-value)<nearest) {
              nearest = fabs(lower_[kRange]-value);
              iRange=kRange;
            }
          }
          assert (iRange>=0);
          iRange--;
          whichRange_[iSequence]=iRange;
          solution[iSequence]=lower_[iRange+1];
        }
        break;
      case ClpSimplex::atLowerBound:
        if (!toNearest) {
          // With increasing tolerances - we may be at wrong place
          // below stops compiler error with gcc 3.2!!!
          if (iSequence==-119)
            printf("ZZ %g %g %g %g\n",lowerValue,value,upperValue,oldTolerance);
          if (fabs(value-lowerValue)>oldTolerance*1.0001) {
            if (fabs(value-upperValue)<=oldTolerance*1.0001) {
              if  (fabs(value-upperValue)>primalTolerance)
                solution[iSequence]=upperValue;
              model_->setStatus(iSequence,ClpSimplex::atUpperBound);
            } else {
              model_->setStatus(iSequence,ClpSimplex::superBasic);
            }
          } else if  (fabs(value-lowerValue)>primalTolerance) {
            solution[iSequence]=lowerValue;
          }
        } else {
          // Set to nearest and make at lower bound
          int kRange;
          iRange=-1;
          double nearest = COIN_DBL_MAX;
          for (kRange=start; kRange<end;kRange++) {
            if (fabs(lower_[kRange]-value)<nearest) {
              nearest = fabs(lower_[kRange]-value);
              iRange=kRange;
            }
          }
          assert (iRange>=0);
          whichRange_[iSequence]=iRange;
          solution[iSequence]=lower_[iRange];
        }
        break;
      case ClpSimplex::isFixed:
        if (toNearest) {
          // Set to true fixed
          for (iRange=start; iRange<end;iRange++) {
            if (lower_[iRange]==lower_[iRange+1])
              break;
          }
          if (iRange==end) {
            // Odd - but make sensible
            // Set to nearest and make at lower bound
            int kRange;
            iRange=-1;
            double nearest = COIN_DBL_MAX;
            for (kRange=start; kRange<end;kRange++) {
              if (fabs(lower_[kRange]-value)<nearest) {
                nearest = fabs(lower_[kRange]-value);
                iRange=kRange;
              }
            }
            assert (iRange>=0);
            whichRange_[iSequence]=iRange;
            if (lower_[iRange]!=lower_[iRange+1])
              model_->setStatus(iSequence,ClpSimplex::atLowerBound);
            else
              model_->setStatus(iSequence,ClpSimplex::atUpperBound);
          }
          solution[iSequence]=lower_[iRange];
        }
        break;
      }
      lower[iSequence] = lower_[iRange];
      upper[iSequence] = lower_[iRange+1];
      cost[iSequence] = cost_[iRange];
      feasibleCost_ += thisFeasibleCost*solution[iSequence];
      //assert (iRange==whichRange_[iSequence]);
    }
  }
#ifdef NONLIN_DEBUG
  double saveCost=feasibleCost_;
  if (method_==3) {
    feasibleCost_=0.0;
    // Put back solution as we will be checking
    unsigned char * statusA = model_->statusArray();
    for (iSequence=0;iSequence<numberTotal;iSequence++) {
      double value = solution[iSequence];
      solution[iSequence]=saveSolution[iSequence];
      saveSolution[iSequence]=value;
      value = lower[iSequence];
      lower[iSequence]=saveLower[iSequence];
      saveLower[iSequence]=value;
      value = upper[iSequence];
      upper[iSequence]=saveUpper[iSequence];
      saveUpper[iSequence]=value;
      unsigned char value2 = statusA[iSequence];
      statusA[iSequence]=saveStatus[iSequence];
      saveStatus[iSequence]=value2;
    }
  }
#endif
  if (CLP_METHOD2) {
    // nonbasic should be at a valid bound
    for (iSequence=0;iSequence<numberTotal;iSequence++) {
      double value=solution[iSequence];
      int iStatus = status_[iSequence];
      assert (currentStatus(iStatus)==CLP_SAME);
      double lowerValue=lower[iSequence];
      double upperValue=upper[iSequence];
      double costValue = cost2_[iSequence];
      double trueCost = costValue;
      int iWhere = originalStatus(iStatus);
      if (iWhere==CLP_BELOW_LOWER) {
        lowerValue=upperValue;
        upperValue=bound_[iSequence];
        costValue -= infeasibilityCost;
      } else if (iWhere==CLP_ABOVE_UPPER) {
        upperValue=lowerValue;
        lowerValue=bound_[iSequence];
        costValue += infeasibilityCost;
      }
      // get correct place
      int newWhere=CLP_FEASIBLE;
      ClpSimplex::Status status = model_->getStatus(iSequence);
      if (upperValue==lowerValue && status!=ClpSimplex::isFixed) {
        if (status != ClpSimplex::basic) {
          model_->setStatus(iSequence,ClpSimplex::isFixed);
          status = ClpSimplex::isFixed;
        }
      }
      switch(status) {
        
      case ClpSimplex::basic:
      case ClpSimplex::superBasic:
        if (value-upperValue<=primalTolerance) {
          if (value-lowerValue>=-primalTolerance) {
            // feasible
            //newWhere=CLP_FEASIBLE;
          } else {
            // below
            newWhere=CLP_BELOW_LOWER;
            double infeasibility = lowerValue-value-primalTolerance;
            sumInfeasibilities_ += infeasibility;
            largestInfeasibility_ = CoinMax(largestInfeasibility_,infeasibility);
            costValue = trueCost - infeasibilityCost;
            changeCost_ -= lowerValue*(costValue-cost[iSequence]);
            numberInfeasibilities_++;
          }
        } else {
          // above
          newWhere = CLP_ABOVE_UPPER;
          double infeasibility = value-upperValue-primalTolerance;
          sumInfeasibilities_ += infeasibility;
          largestInfeasibility_ = CoinMax(largestInfeasibility_,infeasibility);
          costValue = trueCost + infeasibilityCost;
          changeCost_ -= upperValue*(costValue-cost[iSequence]);
          numberInfeasibilities_++;
        }
        break;
      case ClpSimplex::isFree:
        break;
      case ClpSimplex::atUpperBound:
        if (!toNearest) {
          // With increasing tolerances - we may be at wrong place
          if (fabs(value-upperValue)>oldTolerance*1.0001) {
            if (fabs(value-lowerValue)<=oldTolerance*1.0001) {
              if  (fabs(value-lowerValue)>primalTolerance) {
                solution[iSequence]=lowerValue;
                value = lowerValue;
              }
              model_->setStatus(iSequence,ClpSimplex::atLowerBound);
            } else {
              if (value<upperValue) {
                if (value>lowerValue) {
                  model_->setStatus(iSequence,ClpSimplex::superBasic);
                } else {
                  // set to lower bound as infeasible
                  solution[iSequence]=lowerValue;
                  value = lowerValue;
                  model_->setStatus(iSequence,ClpSimplex::atLowerBound);
                }
              } else {
                // set to upper bound as infeasible
                solution[iSequence]=upperValue;
                value = upperValue;
              }
            }
          } else if  (fabs(value-upperValue)>primalTolerance) {
            solution[iSequence]=upperValue;
            value = upperValue;
          }
        } else {
          // Set to nearest and make at bound
          if (fabs(value-lowerValue)<fabs(value-upperValue)) {
            solution[iSequence]=lowerValue;
            value = lowerValue;
            model_->setStatus(iSequence,ClpSimplex::atLowerBound);
          } else {
            solution[iSequence]=upperValue;
            value = upperValue;
          }
        }
        break;
      case ClpSimplex::atLowerBound:
        if (!toNearest) {
          // With increasing tolerances - we may be at wrong place
          if (fabs(value-lowerValue)>oldTolerance*1.0001) {
            if (fabs(value-upperValue)<=oldTolerance*1.0001) {
              if  (fabs(value-upperValue)>primalTolerance) {
                solution[iSequence]=upperValue;
                value = upperValue;
              }
              model_->setStatus(iSequence,ClpSimplex::atUpperBound);
            } else {
              if (value<upperValue) {
                if (value>lowerValue) {
                  model_->setStatus(iSequence,ClpSimplex::superBasic);
                } else {
                  // set to lower bound as infeasible
                  solution[iSequence]=lowerValue;
                  value = lowerValue;
                }
              } else {
                // set to upper bound as infeasible
                solution[iSequence]=upperValue;
                value = upperValue;
                model_->setStatus(iSequence,ClpSimplex::atUpperBound);
              }
            }
          } else if  (fabs(value-lowerValue)>primalTolerance) {
            solution[iSequence]=lowerValue;
            value = lowerValue;
          }
        } else {
          // Set to nearest and make at bound
          if (fabs(value-lowerValue)<fabs(value-upperValue)) {
            solution[iSequence]=lowerValue;
            value = lowerValue;
          } else {
            solution[iSequence]=upperValue;
            value = upperValue;
            model_->setStatus(iSequence,ClpSimplex::atUpperBound);
          }
        }
        break;
      case ClpSimplex::isFixed:
        solution[iSequence]=lowerValue;
        value = lowerValue;
        break;
      }
#ifdef NONLIN_DEBUG
      double lo=saveLower[iSequence];
      double up=saveUpper[iSequence];
      double cc=cost[iSequence];
      unsigned char ss=saveStatus[iSequence];
      unsigned char snow=model_->statusArray()[iSequence];
#endif
      if (iWhere!=newWhere) {
        setOriginalStatus(status_[iSequence],newWhere);
        if (newWhere==CLP_BELOW_LOWER) {
          bound_[iSequence]=upperValue;
          upperValue=lowerValue;
          lowerValue=-COIN_DBL_MAX;
          costValue = trueCost - infeasibilityCost;
        } else if (newWhere==CLP_ABOVE_UPPER) {
          bound_[iSequence]=lowerValue;
          lowerValue=upperValue;
          upperValue=COIN_DBL_MAX;
          costValue = trueCost + infeasibilityCost;
        } else {
          costValue = trueCost;
        }
        lower[iSequence] = lowerValue;
        upper[iSequence] = upperValue;
      }
      // always do as other things may change
      cost[iSequence] = costValue;
#ifdef NONLIN_DEBUG
      if (method_==3) {
        assert (ss==snow);
        assert (cc==cost[iSequence]);
        assert (lo==lower[iSequence]);
        assert (up==upper[iSequence]);
        assert (value==saveSolution[iSequence]);
      }
#endif
      feasibleCost_ += trueCost*value;
    }
  }
#ifdef NONLIN_DEBUG
  if (method_==3)
    assert (fabs(saveCost-feasibleCost_)<0.001*(1.0+CoinMax(fabs(saveCost),fabs(feasibleCost_))));
  delete [] saveSolution;
  delete [] saveLower;
  delete [] saveUpper;
  delete [] saveStatus;
#endif
  //feasibleCost_ /= (model_->rhsScale()*model_->objScale());
}
// Puts feasible bounds into lower and upper
void 
ClpNonLinearCost::feasibleBounds()
{
  if (CLP_METHOD2) {
    int iSequence;
    double * upper = model_->upperRegion();
    double * lower = model_->lowerRegion();
    double * cost = model_->costRegion();
    int numberTotal = numberColumns_+numberRows_;
    for (iSequence=0;iSequence<numberTotal;iSequence++) {
      int iStatus = status_[iSequence];
      assert (currentStatus(iStatus)==CLP_SAME);
      double lowerValue=lower[iSequence];
      double upperValue=upper[iSequence];
      double costValue = cost2_[iSequence];
      int iWhere = originalStatus(iStatus);
      if (iWhere==CLP_BELOW_LOWER) {
        lowerValue=upperValue;
        upperValue=bound_[iSequence];
      } else if (iWhere==CLP_ABOVE_UPPER) {
        upperValue=lowerValue;
        lowerValue=bound_[iSequence];
      }
      setOriginalStatus(status_[iSequence],CLP_FEASIBLE);
      lower[iSequence] = lowerValue;
      upper[iSequence] = upperValue;
      cost[iSequence] = costValue;
    }
  }
}
/* Goes through one bound for each variable.
   If array[i]*multiplier>0 goes down, otherwise up.
   The indices are row indices and need converting to sequences
*/
void 
ClpNonLinearCost::goThru(int numberInArray, double multiplier,
	      const int * index, const double * array,
			 double * rhs)
{
  assert (model_!=NULL);
  abort();
  const int * pivotVariable = model_->pivotVariable();
  if (CLP_METHOD1) {
    for (int i=0;i<numberInArray;i++) {
      int iRow = index[i];
      int iSequence = pivotVariable[iRow];
      double alpha = multiplier*array[iRow];
      // get where in bound sequence
      int iRange = whichRange_[iSequence];
      iRange += offset_[iSequence]; //add temporary bias
      double value = model_->solution(iSequence);
      if (alpha>0.0) {
        // down one
        iRange--;
        assert(iRange>=start_[iSequence]);
        rhs[iRow] = value - lower_[iRange];
      } else {
        // up one
        iRange++;
        assert(iRange<start_[iSequence+1]-1);
        rhs[iRow] = lower_[iRange+1] - value;
      }
      offset_[iSequence] = iRange - whichRange_[iSequence];
    }
  }
#ifdef NONLIN_DEBUG
  double * saveRhs = NULL;
  if (method_==3) {
    int numberRows = model_->numberRows();
    saveRhs = CoinCopyOfArray(rhs,numberRows);
  }
#endif
  if (CLP_METHOD2) {
    const double * solution = model_->solutionRegion();
    const double * upper = model_->upperRegion();
    const double * lower = model_->lowerRegion();
    for (int i=0;i<numberInArray;i++) {
      int iRow = index[i];
      int iSequence = pivotVariable[iRow];
      double alpha = multiplier*array[iRow];
      double value = solution[iSequence];
      int iStatus = status_[iSequence];
      int iWhere = currentStatus(iStatus);
      if (iWhere==CLP_SAME)
        iWhere = originalStatus(iStatus);
      if (iWhere==CLP_FEASIBLE) {
        if (alpha>0.0) {
          // going below
          iWhere=CLP_BELOW_LOWER;
          rhs[iRow] = value - lower[iSequence];
        } else {
          // going above
          iWhere=CLP_ABOVE_UPPER;
          rhs[iRow] = upper[iSequence] - value;
        }
      } else if(iWhere==CLP_BELOW_LOWER) {
        assert (alpha<0);
        // going feasible
        iWhere=CLP_FEASIBLE;
        rhs[iRow] = upper[iSequence] - value;
      } else {
        assert (iWhere==CLP_ABOVE_UPPER);
        // going feasible
        iWhere=CLP_FEASIBLE;
        rhs[iRow] = value - lower[iSequence];
      }
#ifdef NONLIN_DEBUG
      if (method_==3)
        assert (rhs[iRow]==saveRhs[iRow]);
#endif
      setCurrentStatus(status_[iSequence],iWhere);
    }
  }
#ifdef NONLIN_DEBUG
  delete [] saveRhs;
#endif
}
/* Takes off last iteration (i.e. offsets closer to 0)
 */
void 
ClpNonLinearCost::goBack(int numberInArray, const int * index, 
	      double * rhs)
{
  assert (model_!=NULL);
  abort();
  const int * pivotVariable = model_->pivotVariable();
  if (CLP_METHOD1) {
    for (int i=0;i<numberInArray;i++) {
      int iRow = index[i];
      int iSequence = pivotVariable[iRow];
      // get where in bound sequence
      int iRange = whichRange_[iSequence];
      // get closer to original
      if (offset_[iSequence]>0) {
        offset_[iSequence]--;
        assert (offset_[iSequence]>=0);
        iRange += offset_[iSequence]; //add temporary bias
        double value = model_->solution(iSequence);
        // up one
        assert(iRange<start_[iSequence+1]-1);
        rhs[iRow] = lower_[iRange+1] - value; // was earlier lower_[iRange]
      } else {
        offset_[iSequence]++;
        assert (offset_[iSequence]<=0);
        iRange += offset_[iSequence]; //add temporary bias
        double value = model_->solution(iSequence);
        // down one
        assert(iRange>=start_[iSequence]);
        rhs[iRow] = value - lower_[iRange]; // was earlier lower_[iRange+1]
      }
    }
  }
#ifdef NONLIN_DEBUG
  double * saveRhs = NULL;
  if (method_==3) {
    int numberRows = model_->numberRows();
    saveRhs = CoinCopyOfArray(rhs,numberRows);
  }
#endif
  if (CLP_METHOD2) {
    const double * solution = model_->solutionRegion();
    const double * upper = model_->upperRegion();
    const double * lower = model_->lowerRegion();
    for (int i=0;i<numberInArray;i++) {
      int iRow = index[i];
      int iSequence = pivotVariable[iRow];
      double value = solution[iSequence];
      int iStatus = status_[iSequence];
      int iWhere = currentStatus(iStatus);
      int original = originalStatus(iStatus);
      assert (iWhere!=CLP_SAME);
      if (iWhere==CLP_FEASIBLE) {
        if (original==CLP_BELOW_LOWER) {
          // going below
          iWhere=CLP_BELOW_LOWER;
          rhs[iRow] = lower[iSequence] - value;
        } else {
          // going above
          iWhere=CLP_ABOVE_UPPER;
          rhs[iRow] = value - upper[iSequence];
        }
      } else if(iWhere==CLP_BELOW_LOWER) {
        // going feasible
        iWhere=CLP_FEASIBLE;
        rhs[iRow] = value - upper[iSequence];
      } else {
        // going feasible
        iWhere=CLP_FEASIBLE;
        rhs[iRow] = lower[iSequence] - value;
      }
#ifdef NONLIN_DEBUG
      if (method_==3)
        assert (rhs[iRow]==saveRhs[iRow]);
#endif
      setCurrentStatus(status_[iSequence],iWhere);
    }
  }
#ifdef NONLIN_DEBUG
  delete [] saveRhs;
#endif
}
void 
ClpNonLinearCost::goBackAll(const CoinIndexedVector * update)
{
  assert (model_!=NULL);
  const int * pivotVariable = model_->pivotVariable();
  int number = update->getNumElements();
  const int * index = update->getIndices();
  if (CLP_METHOD1) {
    for (int i=0;i<number;i++) {
      int iRow = index[i];
      int iSequence = pivotVariable[iRow];
      offset_[iSequence]=0;
    }
#ifdef CLP_DEBUG
    for (i=0;i<numberRows_+numberColumns_;i++) 
      assert(!offset_[i]);
#endif
  }
  if (CLP_METHOD2) {
    for (int i=0;i<number;i++) {
      int iRow = index[i];
      int iSequence = pivotVariable[iRow];
      setSameStatus(status_[iSequence]);
    }
  }
}
void 
ClpNonLinearCost::checkInfeasibilities(int numberInArray, const int * index)
{
  assert (model_!=NULL);
  double primalTolerance = model_->currentPrimalTolerance();
  const int * pivotVariable = model_->pivotVariable();
  if (CLP_METHOD1) {
    for (int i=0;i<numberInArray;i++) {
      int iRow = index[i];
      int iSequence = pivotVariable[iRow];
      // get where in bound sequence
      int iRange;
      int currentRange = whichRange_[iSequence];
      double value = model_->solution(iSequence);
      int start = start_[iSequence];
      int end = start_[iSequence+1]-1;
      for (iRange=start; iRange<end;iRange++) {
        if (value<lower_[iRange+1]+primalTolerance) {
          // put in better range
          if (value>=lower_[iRange+1]-primalTolerance&&infeasible(iRange)&&iRange==start) 
            iRange++;
          break;
        }
      }
      assert(iRange<end);
      assert(model_->getStatus(iSequence)==ClpSimplex::basic);
      double & lower = model_->lowerAddress(iSequence);
      double & upper = model_->upperAddress(iSequence);
      double & cost = model_->costAddress(iSequence);
      whichRange_[iSequence]=iRange;
      if (iRange!=currentRange) {
        if (infeasible(iRange))
          numberInfeasibilities_++;
        if (infeasible(currentRange))
          numberInfeasibilities_--;
      }
      lower = lower_[iRange];
      upper = lower_[iRange+1];
      cost = cost_[iRange];
    }
  }
  if (CLP_METHOD2) {
    double * solution = model_->solutionRegion();
    double * upper = model_->upperRegion();
    double * lower = model_->lowerRegion();
    double * cost = model_->costRegion();
    for (int i=0;i<numberInArray;i++) {
      int iRow = index[i];
      int iSequence = pivotVariable[iRow];
      double value=solution[iSequence];
      int iStatus = status_[iSequence];
      assert (currentStatus(iStatus)==CLP_SAME);
      double lowerValue=lower[iSequence];
      double upperValue=upper[iSequence];
      double costValue = cost2_[iSequence];
      int iWhere = originalStatus(iStatus);
      if (iWhere==CLP_BELOW_LOWER) {
        lowerValue=upperValue;
        upperValue=bound_[iSequence];
        numberInfeasibilities_--;
      } else if (iWhere==CLP_ABOVE_UPPER) {
        upperValue=lowerValue;
        lowerValue=bound_[iSequence];
        numberInfeasibilities_--;
      }
      // get correct place
      int newWhere=CLP_FEASIBLE;
      if (value-upperValue<=primalTolerance) {
        if (value-lowerValue>=-primalTolerance) {
          // feasible
          //newWhere=CLP_FEASIBLE;
        } else {
          // below
          newWhere=CLP_BELOW_LOWER;
          costValue -= infeasibilityWeight_;
          numberInfeasibilities_++;
        }
      } else {
        // above
        newWhere = CLP_ABOVE_UPPER;
        costValue += infeasibilityWeight_;
        numberInfeasibilities_++;
      }
      if (iWhere!=newWhere) {
        setOriginalStatus(status_[iSequence],newWhere);
        if (newWhere==CLP_BELOW_LOWER) {
          bound_[iSequence]=upperValue;
          upperValue=lowerValue;
          lowerValue=-COIN_DBL_MAX;
        } else if (newWhere==CLP_ABOVE_UPPER) {
          bound_[iSequence]=lowerValue;
          lowerValue=upperValue;
          upperValue=COIN_DBL_MAX;
        }
        lower[iSequence] = lowerValue;
        upper[iSequence] = upperValue;
        cost[iSequence] = costValue;
      }
    }
  }
}
/* Puts back correct infeasible costs for each variable
   The input indices are row indices and need converting to sequences
   for costs.
   On input array is empty (but indices exist).  On exit just
   changed costs will be stored as normal CoinIndexedVector
*/
void 
ClpNonLinearCost::checkChanged(int numberInArray, CoinIndexedVector * update)
{
  assert (model_!=NULL);
  double primalTolerance = model_->currentPrimalTolerance();
  const int * pivotVariable = model_->pivotVariable();
  int number=0;
  int * index = update->getIndices();
  double * work = update->denseVector();
  if (CLP_METHOD1) {
    for (int i=0;i<numberInArray;i++) {
      int iRow = index[i];
      int iSequence = pivotVariable[iRow];
      // get where in bound sequence
      int iRange;
      double value = model_->solution(iSequence);
      int start = start_[iSequence];
      int end = start_[iSequence+1]-1;
      for (iRange=start; iRange<end;iRange++) {
        if (value<lower_[iRange+1]+primalTolerance) {
          // put in better range
          if (value>=lower_[iRange+1]-primalTolerance&&infeasible(iRange)&&iRange==start) 
            iRange++;
          break;
        }
      }
      assert(iRange<end);
      assert(model_->getStatus(iSequence)==ClpSimplex::basic);
      int jRange = whichRange_[iSequence];
      if (iRange!=jRange) {
        // changed
        work[iRow] = cost_[jRange]-cost_[iRange];
        index[number++]=iRow;
        double & lower = model_->lowerAddress(iSequence);
        double & upper = model_->upperAddress(iSequence);
        double & cost = model_->costAddress(iSequence);
        whichRange_[iSequence]=iRange;
        if (infeasible(iRange))
          numberInfeasibilities_++;
        if (infeasible(jRange))
          numberInfeasibilities_--;
        lower = lower_[iRange];
        upper = lower_[iRange+1];
        cost = cost_[iRange];
      }
    }
  }
  if (CLP_METHOD2) {
    double * solution = model_->solutionRegion();
    double * upper = model_->upperRegion();
    double * lower = model_->lowerRegion();
    double * cost = model_->costRegion();
    for (int i=0;i<numberInArray;i++) {
      int iRow = index[i];
      int iSequence = pivotVariable[iRow];
      double value=solution[iSequence];
      int iStatus = status_[iSequence];
      assert (currentStatus(iStatus)==CLP_SAME);
      double lowerValue=lower[iSequence];
      double upperValue=upper[iSequence];
      double costValue = cost2_[iSequence];
      int iWhere = originalStatus(iStatus);
      if (iWhere==CLP_BELOW_LOWER) {
        lowerValue=upperValue;
        upperValue=bound_[iSequence];
        numberInfeasibilities_--;
      } else if (iWhere==CLP_ABOVE_UPPER) {
        upperValue=lowerValue;
        lowerValue=bound_[iSequence];
        numberInfeasibilities_--;
      }
      // get correct place
      int newWhere=CLP_FEASIBLE;
      if (value-upperValue<=primalTolerance) {
        if (value-lowerValue>=-primalTolerance) {
          // feasible
          //newWhere=CLP_FEASIBLE;
        } else {
          // below
          newWhere=CLP_BELOW_LOWER;
          costValue -= infeasibilityWeight_;
          numberInfeasibilities_++;
        }
      } else {
        // above
        newWhere = CLP_ABOVE_UPPER;
        costValue += infeasibilityWeight_;
        numberInfeasibilities_++;
      }
      if (iWhere!=newWhere) {
        work[iRow] = cost[iSequence]-costValue;
        index[number++]=iRow;
        setOriginalStatus(status_[iSequence],newWhere);
        if (newWhere==CLP_BELOW_LOWER) {
          bound_[iSequence]=upperValue;
          upperValue=lowerValue;
          lowerValue=-COIN_DBL_MAX;
        } else if (newWhere==CLP_ABOVE_UPPER) {
          bound_[iSequence]=lowerValue;
          lowerValue=upperValue;
          upperValue=COIN_DBL_MAX;
        }
        lower[iSequence] = lowerValue;
        upper[iSequence] = upperValue;
        cost[iSequence] = costValue;
      }
    }
  }
  update->setNumElements(number);
}
/* Sets bounds and cost for one variable - returns change in cost*/
double 
ClpNonLinearCost::setOne(int iSequence, double value)
{
  assert (model_!=NULL);
  double primalTolerance = model_->currentPrimalTolerance();
  // difference in cost
  double difference=0.0;
  if (CLP_METHOD1) {
    // get where in bound sequence
    int iRange;
    int currentRange = whichRange_[iSequence];
    int start = start_[iSequence];
    int end = start_[iSequence+1]-1;
    if (!bothWays_) {
      // If fixed try and get feasible
      if (lower_[start+1]==lower_[start+2]&&fabs(value-lower_[start+1])<1.001*primalTolerance) {
        iRange =start+1;
      } else {
        for (iRange=start; iRange<end;iRange++) {
          if (value<=lower_[iRange+1]+primalTolerance) {
            // put in better range
            if (value>=lower_[iRange+1]-primalTolerance&&infeasible(iRange)&&iRange==start) 
              iRange++;
            break;
          }
        }
      }
    } else {
      // leave in current if possible
      iRange = whichRange_[iSequence];
      if (value<lower_[iRange]-primalTolerance||value>lower_[iRange+1]+primalTolerance) {
        for (iRange=start; iRange<end;iRange++) {
          if (value<lower_[iRange+1]+primalTolerance) {
            // put in better range
            if (value>=lower_[iRange+1]-primalTolerance&&infeasible(iRange)&&iRange==start) 
              iRange++;
            break;
          }
        }
      }
    }
    assert(iRange<end);
    whichRange_[iSequence]=iRange;
    if (iRange!=currentRange) {
      if (infeasible(iRange))
        numberInfeasibilities_++;
      if (infeasible(currentRange))
        numberInfeasibilities_--;
    }
    double & lower = model_->lowerAddress(iSequence);
    double & upper = model_->upperAddress(iSequence);
    double & cost = model_->costAddress(iSequence);
    lower = lower_[iRange];
    upper = lower_[iRange+1];
    ClpSimplex::Status status = model_->getStatus(iSequence);
    if (upper==lower) {
      if (status != ClpSimplex::basic) {
        model_->setStatus(iSequence,ClpSimplex::isFixed);
        status = ClpSimplex::basic; // so will skip
      }
    }
    switch(status) {
      
    case ClpSimplex::basic:
    case ClpSimplex::superBasic:
    case ClpSimplex::isFree:
      break;
    case ClpSimplex::atUpperBound:
    case ClpSimplex::atLowerBound:
    case ClpSimplex::isFixed:
      // set correctly
      if (fabs(value-lower)<=primalTolerance*1.001){
        model_->setStatus(iSequence,ClpSimplex::atLowerBound);
      } else if (fabs(value-upper)<=primalTolerance*1.001){
        model_->setStatus(iSequence,ClpSimplex::atUpperBound);
      } else {
        // set superBasic
        model_->setStatus(iSequence,ClpSimplex::superBasic);
      }
      break;
    }
    difference = cost-cost_[iRange]; 
    cost = cost_[iRange];
  }
  if (CLP_METHOD2) {
    double * upper = model_->upperRegion();
    double * lower = model_->lowerRegion();
    double * cost = model_->costRegion();
    int iStatus = status_[iSequence];
    assert (currentStatus(iStatus)==CLP_SAME);
    double lowerValue=lower[iSequence];
    double upperValue=upper[iSequence];
    double costValue = cost2_[iSequence];
    int iWhere = originalStatus(iStatus);
    if (iWhere==CLP_BELOW_LOWER) {
      lowerValue=upperValue;
      upperValue=bound_[iSequence];
      numberInfeasibilities_--;
    } else if (iWhere==CLP_ABOVE_UPPER) {
      upperValue=lowerValue;
      lowerValue=bound_[iSequence];
      numberInfeasibilities_--;
    }
    // get correct place
    int newWhere=CLP_FEASIBLE;
    if (value-upperValue<=primalTolerance) {
      if (value-lowerValue>=-primalTolerance) {
        // feasible
        //newWhere=CLP_FEASIBLE;
      } else {
        // below
        newWhere=CLP_BELOW_LOWER;
        costValue -= infeasibilityWeight_;
        numberInfeasibilities_++;
      }
    } else {
      // above
      newWhere = CLP_ABOVE_UPPER;
      costValue += infeasibilityWeight_;
      numberInfeasibilities_++;
    }
    if (iWhere!=newWhere) {
      difference = cost[iSequence]-costValue; 
      setOriginalStatus(status_[iSequence],newWhere);
      if (newWhere==CLP_BELOW_LOWER) {
        bound_[iSequence]=upperValue;
        upperValue=lowerValue;
        lowerValue=-COIN_DBL_MAX;
      } else if (newWhere==CLP_ABOVE_UPPER) {
        bound_[iSequence]=lowerValue;
        lowerValue=upperValue;
        upperValue=COIN_DBL_MAX;
      }
      lower[iSequence] = lowerValue;
      upper[iSequence] = upperValue;
      cost[iSequence] = costValue;
    }
    ClpSimplex::Status status = model_->getStatus(iSequence);
    if (upperValue==lowerValue) {
      if (status != ClpSimplex::basic) {
        model_->setStatus(iSequence,ClpSimplex::isFixed);
        status = ClpSimplex::basic; // so will skip
      }
    }
    switch(status) {
      
    case ClpSimplex::basic:
    case ClpSimplex::superBasic:
    case ClpSimplex::isFree:
      break;
    case ClpSimplex::atUpperBound:
    case ClpSimplex::atLowerBound:
    case ClpSimplex::isFixed:
      // set correctly
      if (fabs(value-lowerValue)<=primalTolerance*1.001){
        model_->setStatus(iSequence,ClpSimplex::atLowerBound);
      } else if (fabs(value-upperValue)<=primalTolerance*1.001){
        model_->setStatus(iSequence,ClpSimplex::atUpperBound);
      } else {
        // set superBasic
        model_->setStatus(iSequence,ClpSimplex::superBasic);
      }
      break;
    }
  }
  changeCost_ += value*difference;
  return difference;
}
/* Sets bounds and infeasible cost and true cost for one variable 
   This is for gub and column generation etc */
void 
ClpNonLinearCost::setOne(int sequence, double solutionValue, double lowerValue, double upperValue,
			 double costValue)
{
  if (CLP_METHOD1) {
    int iRange=-1;
    int start = start_[sequence];
    double infeasibilityCost = model_->infeasibilityCost();
    cost_[start] = costValue-infeasibilityCost;
    lower_[start+1]=lowerValue;
    cost_[start+1] = costValue;
    lower_[start+2]=upperValue;
    cost_[start+2] = costValue+infeasibilityCost;
    double primalTolerance = model_->currentPrimalTolerance();
    if (solutionValue-lowerValue>=-primalTolerance) {
      if (solutionValue-upperValue<=primalTolerance) {
        iRange=start+1;
      } else {
        iRange=start+2;
      }
    } else {
      iRange = start;
    }
    model_->costRegion()[sequence]=cost_[iRange];
    whichRange_[sequence]=iRange;
  }
  if (CLP_METHOD2) {
    abort(); // may never work
  }
  
}
/* Sets bounds and cost for outgoing variable 
   may change value
   Returns direction */
int 
ClpNonLinearCost::setOneOutgoing(int iSequence, double & value)
{
  assert (model_!=NULL);
  double primalTolerance = model_->currentPrimalTolerance();
  // difference in cost
  double difference=0.0;
  int direction=0;
  if (CLP_METHOD1) {
    // get where in bound sequence
    int iRange;
    int currentRange = whichRange_[iSequence];
    int start = start_[iSequence];
    int end = start_[iSequence+1]-1;
    // Set perceived direction out
    if (value<=lower_[currentRange]+1.001*primalTolerance) {
      direction=1;
    } else if (value>=lower_[currentRange+1]-1.001*primalTolerance) {
      direction=-1;
    } else {
      // odd
      direction=0;
    }
    // If fixed try and get feasible
    if (lower_[start+1]==lower_[start+2]&&fabs(value-lower_[start+1])<1.001*primalTolerance) {
      iRange =start+1;
    } else {
      // See if exact
      for (iRange=start; iRange<end;iRange++) {
        if (value==lower_[iRange+1]) {
          // put in better range
          if (infeasible(iRange)&&iRange==start) 
            iRange++;
          break;
        }
      }
      if (iRange==end) {
        // not exact
        for (iRange=start; iRange<end;iRange++) {
          if (value<=lower_[iRange+1]+primalTolerance) {
            // put in better range
            if (value>=lower_[iRange+1]-primalTolerance&&infeasible(iRange)&&iRange==start) 
              iRange++;
            break;
          }
        }
      }
    }
    assert(iRange<end);
    whichRange_[iSequence]=iRange;
    if (iRange!=currentRange) {
      if (infeasible(iRange))
        numberInfeasibilities_++;
      if (infeasible(currentRange))
        numberInfeasibilities_--;
    }
    double & lower = model_->lowerAddress(iSequence);
    double & upper = model_->upperAddress(iSequence);
    double & cost = model_->costAddress(iSequence);
    lower = lower_[iRange];
    upper = lower_[iRange+1];
    if (upper==lower) {
      value=upper;
    } else {
      // set correctly
      if (fabs(value-lower)<=primalTolerance*1.001){
        value = CoinMin(value,lower+primalTolerance);
      } else if (fabs(value-upper)<=primalTolerance*1.001){
        value = CoinMax(value,upper-primalTolerance);
      } else {
        //printf("*** variable wandered off bound %g %g %g!\n",
        //     lower,value,upper);
        if (value-lower<=upper-value) 
          value = lower+primalTolerance;
        else 
          value = upper-primalTolerance;
      }
    }
    difference = cost-cost_[iRange]; 
    cost = cost_[iRange];
  }
  if (CLP_METHOD2) {
    double * upper = model_->upperRegion();
    double * lower = model_->lowerRegion();
    double * cost = model_->costRegion();
    int iStatus = status_[iSequence];
    assert (currentStatus(iStatus)==CLP_SAME);
    double lowerValue=lower[iSequence];
    double upperValue=upper[iSequence];
    double costValue = cost2_[iSequence];
    // Set perceived direction out
    if (value<=lowerValue+1.001*primalTolerance) {
      direction=1;
    } else if (value>=upperValue-1.001*primalTolerance) {
      direction=-1;
    } else {
      // odd
      direction=0;
    }
    int iWhere = originalStatus(iStatus);
    if (iWhere==CLP_BELOW_LOWER) {
      lowerValue=upperValue;
      upperValue=bound_[iSequence];
      numberInfeasibilities_--;
    } else if (iWhere==CLP_ABOVE_UPPER) {
      upperValue=lowerValue;
      lowerValue=bound_[iSequence];
      numberInfeasibilities_--;
    }
    // get correct place
    // If fixed give benefit of doubt
    if (lowerValue==upperValue)
      value=lowerValue;
    int newWhere=CLP_FEASIBLE;
    if (value-upperValue<=primalTolerance) {
      if (value-lowerValue>=-primalTolerance) {
        // feasible
        //newWhere=CLP_FEASIBLE;
      } else {
        // below
        newWhere=CLP_BELOW_LOWER;
        costValue -= infeasibilityWeight_;
        numberInfeasibilities_++;
      }
    } else {
      // above
      newWhere = CLP_ABOVE_UPPER;
      costValue += infeasibilityWeight_;
      numberInfeasibilities_++;
    }
    if (iWhere!=newWhere) {
      difference = cost[iSequence]-costValue; 
      setOriginalStatus(status_[iSequence],newWhere);
      if (newWhere==CLP_BELOW_LOWER) {
        bound_[iSequence]=upperValue;
        upper[iSequence]=lowerValue;
        lower[iSequence]=-COIN_DBL_MAX;
      } else if (newWhere==CLP_ABOVE_UPPER) {
        bound_[iSequence]=lowerValue;
        lower[iSequence]=upperValue;
        upper[iSequence]=COIN_DBL_MAX;
      } else {
        lower[iSequence] = lowerValue;
        upper[iSequence] = upperValue;
      }
      cost[iSequence] = costValue;
    }
    // set correctly
    if (fabs(value-lowerValue)<=primalTolerance*1.001){
      value = CoinMin(value,lowerValue+primalTolerance);
    } else if (fabs(value-upperValue)<=primalTolerance*1.001){
      value = CoinMax(value,upperValue-primalTolerance);
    } else {
      //printf("*** variable wandered off bound %g %g %g!\n",
      //     lowerValue,value,upperValue);
      if (value-lowerValue<=upperValue-value) 
        value = lowerValue+primalTolerance;
      else 
        value = upperValue-primalTolerance;
    }
  }
  changeCost_ += value*difference;
  return direction;
}
// Returns nearest bound
double 
ClpNonLinearCost::nearest(int iSequence, double solutionValue)
{
  assert (model_!=NULL);
  double nearest=0.0;
  if (CLP_METHOD1) {
    // get where in bound sequence
    int iRange;
    int start = start_[iSequence];
    int end = start_[iSequence+1];
    int jRange=-1;
    nearest=COIN_DBL_MAX;
    for (iRange=start; iRange<end;iRange++) {
      if (fabs(solutionValue-lower_[iRange])<nearest) {
        jRange=iRange;
        nearest=fabs(solutionValue-lower_[iRange]);
      }
    }
    assert(jRange<end);
    nearest= lower_[jRange];
  }
  if (CLP_METHOD2) {
    const double * upper = model_->upperRegion();
    const double * lower = model_->lowerRegion();
    double lowerValue=lower[iSequence];
    double upperValue=upper[iSequence];
    int iWhere = originalStatus(status_[iSequence]);
    if (iWhere==CLP_BELOW_LOWER) {
      lowerValue=upperValue;
      upperValue=bound_[iSequence];
    } else if (iWhere==CLP_ABOVE_UPPER) {
      upperValue=lowerValue;
      lowerValue=bound_[iSequence];
    }
    if (fabs(solutionValue-lowerValue)<fabs(solutionValue-upperValue))
      nearest = lowerValue;
    else
      nearest = upperValue;
  }
  return nearest;
}
/// Feasible cost with offset and direction (i.e. for reporting)
double 
ClpNonLinearCost::feasibleReportCost() const
{ 
  double value;
  model_->getDblParam(ClpObjOffset,value);
  return (feasibleCost_+model_->objectiveAsObject()->nonlinearOffset())*model_->optimizationDirection()/
    (model_->objectiveScale()*model_->rhsScale())-value;
}
// Get rid of real costs (just for moment)
void 
ClpNonLinearCost::zapCosts()
{
  int iSequence;
  double infeasibilityCost = model_->infeasibilityCost();
  // zero out all costs
  int numberTotal = numberColumns_+numberRows_;
  if (CLP_METHOD1) {
    int n = start_[numberTotal];
    memset(cost_,0,n*sizeof(double));
    for (iSequence=0;iSequence<numberTotal;iSequence++) {
      int start = start_[iSequence];
      int end = start_[iSequence+1]-1;
      // correct costs for this infeasibility weight
      if (infeasible(start)) {
        cost_[start] = -infeasibilityCost;
      }
      if (infeasible(end-1)) {
        cost_[end-1] = infeasibilityCost;
      }
    }
  }
  if (CLP_METHOD2) {
  }
}
#ifdef VALIDATE
// For debug
void 
ClpNonLinearCost::validate()
{
  double primalTolerance = model_->currentPrimalTolerance();
  int iSequence;
  const double * solution = model_->solutionRegion();
  const double * upper = model_->upperRegion();
  const double * lower = model_->lowerRegion();
  const double * cost = model_->costRegion();
  double infeasibilityCost = model_->infeasibilityCost();
  int numberTotal = numberRows_+numberColumns_;
  int numberInfeasibilities=0;
  double sumInfeasibilities = 0.0;
  
  for (iSequence=0;iSequence<numberTotal;iSequence++) {
    double value=solution[iSequence];
    int iStatus = status_[iSequence];
    assert (currentStatus(iStatus)==CLP_SAME);
    double lowerValue=lower[iSequence];
    double upperValue=upper[iSequence];
    double costValue = cost2_[iSequence];
    int iWhere = originalStatus(iStatus);
    if (iWhere==CLP_BELOW_LOWER) {
      lowerValue=upperValue;
      upperValue=bound_[iSequence];
      costValue -= infeasibilityCost;
      assert (value<=lowerValue-primalTolerance);
      numberInfeasibilities++;
      sumInfeasibilities += lowerValue-value-primalTolerance;
      assert (model_->getStatus(iSequence)==ClpSimplex::basic);
    } else if (iWhere==CLP_ABOVE_UPPER) {
      upperValue=lowerValue;
      lowerValue=bound_[iSequence];
      costValue += infeasibilityCost;
      assert (value>=upperValue+primalTolerance);
      numberInfeasibilities++;
      sumInfeasibilities += value -upperValue-primalTolerance;
      assert (model_->getStatus(iSequence)==ClpSimplex::basic);
    } else {
      assert (value>=lowerValue-primalTolerance&&value<=upperValue+primalTolerance);
    }
    assert (lowerValue==saveLowerV[iSequence]);
    assert (upperValue==saveUpperV[iSequence]);
    assert (costValue==cost[iSequence]);
  }
  if (numberInfeasibilities)
    printf("JJ %d infeasibilities summing to %g\n",
           numberInfeasibilities,sumInfeasibilities);
}
#endif