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

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

#include "CoinIndexedVector.hpp"

#include "ClpSimplex.hpp"
#include "CoinHelperFunctions.hpp"
#include "ClpNonLinearCost.hpp"
#include "ClpEventHandler.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_->clpMatrix()->type() < 15)
    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]);
    }
  }
}
#if 0
// Refresh - assuming regions OK
void 
ClpNonLinearCost::refresh()
{
     int numberTotal = numberRows_ + numberColumns_;
     numberInfeasibilities_ = 0;
     sumInfeasibilities_ = 0.0;
     largestInfeasibility_ = 0.0;
     double infeasibilityCost = model_->infeasibilityCost();
     double primalTolerance = model_->currentPrimalTolerance();
     double * cost = model_->costRegion();
     double * upper = model_->upperRegion();
     double * lower = model_->lowerRegion();
     double * solution = model_->solutionRegion();
     for (int iSequence = 0; iSequence < numberTotal; iSequence++) {
       cost2_[iSequence] = cost[iSequence];
       double value = solution[iSequence];
       double lowerValue = lower[iSequence];
       double upperValue = upper[iSequence];
       if (value - upperValue <= primalTolerance) {
	 if (value - lowerValue >= -primalTolerance) {
	   // feasible
	   status_[iSequence] = static_cast<unsigned char>(CLP_FEASIBLE | (CLP_SAME << 4));
	   bound_[iSequence] = 0.0;
	 } else {
	   // below
	   double infeasibility = lowerValue - value - primalTolerance;
	   sumInfeasibilities_ += infeasibility;
	   largestInfeasibility_ = CoinMax(largestInfeasibility_, infeasibility);
	   cost[iSequence] -= infeasibilityCost;
	   numberInfeasibilities_++;
	   status_[iSequence] = static_cast<unsigned char>(CLP_BELOW_LOWER | (CLP_SAME << 4));
	   bound_[iSequence] = upperValue;
	   upper[iSequence] = lowerValue;
	   lower[iSequence] = -COIN_DBL_MAX;
	 }
       } else {
	 // above
	 double infeasibility = value - upperValue - primalTolerance;
	 sumInfeasibilities_ += infeasibility;
	 largestInfeasibility_ = CoinMax(largestInfeasibility_, infeasibility);
	 cost[iSequence] += infeasibilityCost;
	 numberInfeasibilities_++;
	 status_[iSequence] = static_cast<unsigned char>(CLP_ABOVE_UPPER | (CLP_SAME << 4));
	 bound_[iSequence] = lowerValue;
	 lower[iSequence] = upperValue;
	 upper[iSequence] = COIN_DBL_MAX;
       }
     }
     //     checkInfeasibilities(model_->primalTolerance());
     
}
#endif
// Refresh one- assuming regions OK
void ClpNonLinearCost::refresh(int iSequence)
{
  double infeasibilityCost = model_->infeasibilityCost();
  double primalTolerance = model_->currentPrimalTolerance();
  double *cost = model_->costRegion();
  double *upper = model_->upperRegion();
  double *lower = model_->lowerRegion();
  double *solution = model_->solutionRegion();
  cost2_[iSequence] = cost[iSequence];
  double value = solution[iSequence];
  double lowerValue = lower[iSequence];
  double upperValue = upper[iSequence];
  if (value - upperValue <= primalTolerance) {
    if (value - lowerValue >= -primalTolerance) {
      // feasible
      status_[iSequence] = static_cast< unsigned char >(CLP_FEASIBLE | (CLP_SAME << 4));
      bound_[iSequence] = 0.0;
    } else {
      // below
      cost[iSequence] -= infeasibilityCost;
      status_[iSequence] = static_cast< unsigned char >(CLP_BELOW_LOWER | (CLP_SAME << 4));
      bound_[iSequence] = upperValue;
      upper[iSequence] = lowerValue;
      lower[iSequence] = -COIN_DBL_MAX;
    }
  } else {
    // above
    cost[iSequence] += infeasibilityCost;
    status_[iSequence] = static_cast< unsigned char >(CLP_ABOVE_UPPER | (CLP_SAME << 4));
    bound_[iSequence] = lowerValue;
    lower[iSequence] = upperValue;
    upper[iSequence] = COIN_DBL_MAX;
  }
}
// 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
  CoinMemcpyN(columnCosts, numberColumns_, cost);
  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();
  COIN_DETAIL_PRINT(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];
      CoinMemcpyN(rhs.start_, (numberTotal + 1), start_);
      whichRange_ = new int[numberTotal];
      CoinMemcpyN(rhs.whichRange_, numberTotal, whichRange_);
      offset_ = new int[numberTotal];
      CoinMemcpyN(rhs.offset_, numberTotal, offset_);
      int numberEntries = start_[numberTotal];
      lower_ = new double[numberEntries];
      CoinMemcpyN(rhs.lower_, numberEntries, lower_);
      cost_ = new double[numberEntries];
      CoinMemcpyN(rhs.cost_, numberEntries, cost_);
      infeasible_ = new unsigned int[(numberEntries + 31) >> 5];
      CoinMemcpyN(rhs.infeasible_, ((numberEntries + 31) >> 5), infeasible_);
    }
    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];
        CoinMemcpyN(rhs.start_, (numberTotal + 1), start_);
        whichRange_ = new int[numberTotal];
        CoinMemcpyN(rhs.whichRange_, numberTotal, whichRange_);
        offset_ = new int[numberTotal];
        CoinMemcpyN(rhs.offset_, numberTotal, offset_);
        int numberEntries = start_[numberTotal];
        lower_ = new double[numberEntries];
        CoinMemcpyN(rhs.lower_, numberEntries, lower_);
        cost_ = new double[numberEntries];
        CoinMemcpyN(rhs.cost_, numberEntries, cost_);
        infeasible_ = new unsigned int[(numberEntries + 31) >> 5];
        CoinMemcpyN(rhs.infeasible_, ((numberEntries + 31) >> 5), infeasible_);
      }
      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;
        }
      }
      //#define PRINT_DETAIL7 2
#if PRINT_DETAIL7 > 1
      printf("NL %d sol %g bounds %g %g\n",
        iSequence, solution[iSequence],
        lowerValue, upperValue);
#endif
      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
#if PRINT_DETAIL7
              printf("**NL %d sol %g below %g\n",
                iSequence, solution[iSequence],
                lowerValue);
#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
#if PRINT_DETAIL7
              printf("**NL %d sol %g above %g\n",
                iSequence, 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) {
    //#define CLP_NON_JUST_BASIC
#ifndef CLP_NON_JUST_BASIC
    // nonbasic should be at a valid bound
    for (iSequence = 0; iSequence < numberTotal; iSequence++) {
#else
    const int *pivotVariable = model_->pivotVariable();
    for (int i = 0; i < numberRows_; i++) {
      int iSequence = pivotVariable[i];
#endif
      double value = solution[iSequence];
      unsigned char 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;
            assert(fabs(lowerValue) < 1.0e100);
            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;
        assert(lowerValue <= 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++) {
      unsigned char 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];
        assert(fabs(lowerValue) < 1.0e100);
      } 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];
      unsigned char 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];
      unsigned char 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];
      unsigned char 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_--;
        assert(fabs(lowerValue) < 1.0e100);
      } 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;
          assert(fabs(lowerValue) < 1.0e100);
          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];
      unsigned char 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_--;
        assert(fabs(lowerValue) < 1.0e100);
      } 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_++;
          assert(fabs(lowerValue) < 1.0e100);
        }
      } 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();
    unsigned char iStatus = status_[iSequence];
    assert(currentStatus(iStatus) == CLP_SAME);
    double lowerValue = lower[iSequence];
    double upperValue = upper[iSequence];
    double costValue = cost2_[iSequence];
    int iWhere = originalStatus(iStatus);
#undef CLP_USER_DRIVEN
#ifdef CLP_USER_DRIVEN
    clpUserStruct info;
    info.type = 3;
    info.sequence = iSequence;
    info.value = value;
    info.change = 0;
    info.lower = lowerValue;
    info.upper = upperValue;
    info.cost = costValue;
#endif
    if (iWhere == CLP_BELOW_LOWER) {
      lowerValue = upperValue;
      upperValue = bound_[iSequence];
      numberInfeasibilities_--;
      assert(fabs(lowerValue) < 1.0e100);
    } 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_++;
        assert(fabs(lowerValue) < 1.0e100);
      }
    } 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;
    }
#ifdef CLP_USER_DRIVEN
    model_->eventHandler()->eventWithInfo(ClpEventHandler::pivotRow,
      &info);
    if (info.change) {
    }
#endif
  }
  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) {
    bound_[sequence] = 0.0;
    cost2_[sequence] = costValue;
    setInitialStatus(status_[sequence]);
  }
}
/* 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;
#ifdef CLP_USER_DRIVEN
  double saveLower = model_->lowerRegion()[iSequence];
  double saveUpper = model_->upperRegion()[iSequence];
  double saveCost = model_->costRegion()[iSequence];
#endif
  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();
    unsigned char 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_--;
      assert(fabs(lowerValue) < 1.0e100);
    } 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_++;
        assert(fabs(lowerValue) < 1.0e100);
      }
    } 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;
#ifdef CLP_USER_DRIVEN
  //if (iSequence>=model_->numberColumns)
  if (difference || saveCost != model_->costRegion()[iSequence]) {
    printf("Out %d (%d) %g <= %g <= %g cost %g x=>x %g < %g cost %g\n",
      iSequence, iSequence - model_->numberColumns(),
      saveLower, value, saveUpper, saveCost,
      model_->lowerRegion()[iSequence],
      model_->upperRegion()[iSequence],
      model_->costRegion()[iSequence]);
  }
#endif
  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];
      assert(fabs(lowerValue) < 1.0e100);
    } 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];
      assert(fabs(lowerValue) < 1.0e100);
      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

/* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
*/