// 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
syntax highlighted by Code2HTML, v. 0.9.1