// Copyright (C) 2002, International Business Machines
// Corporation and others. All Rights Reserved.
#include "CoinPragma.hpp"
#include "ClpSimplex.hpp"
#include "ClpDualRowSteepest.hpp"
#include "CoinIndexedVector.hpp"
#include "ClpFactorization.hpp"
#include "CoinHelperFunctions.hpp"
#include <cstdio>
//#############################################################################
// Constructors / Destructor / Assignment
//#############################################################################
//-------------------------------------------------------------------
// Default Constructor
//-------------------------------------------------------------------
ClpDualRowSteepest::ClpDualRowSteepest (int mode)
: ClpDualRowPivot(),
state_(-1),
mode_(mode),
persistence_(normal),
weights_(NULL),
infeasible_(NULL),
alternateWeights_(NULL),
savedWeights_(NULL),
dubiousWeights_(NULL)
{
type_=2+64*mode;
}
//-------------------------------------------------------------------
// Copy constructor
//-------------------------------------------------------------------
ClpDualRowSteepest::ClpDualRowSteepest (const ClpDualRowSteepest & rhs)
: ClpDualRowPivot(rhs)
{
state_=rhs.state_;
mode_ = rhs.mode_;
persistence_ = rhs.persistence_;
model_ = rhs.model_;
if ((model_&&model_->whatsChanged()&1)!=0) {
int number = model_->numberRows();
if (rhs.savedWeights_)
number = CoinMin(number,rhs.savedWeights_->capacity());
if (rhs.infeasible_) {
infeasible_= new CoinIndexedVector(rhs.infeasible_);
} else {
infeasible_=NULL;
}
if (rhs.weights_) {
weights_= new double[number];
ClpDisjointCopyN(rhs.weights_,number,weights_);
} else {
weights_=NULL;
}
if (rhs.alternateWeights_) {
alternateWeights_= new CoinIndexedVector(rhs.alternateWeights_);
} else {
alternateWeights_=NULL;
}
if (rhs.savedWeights_) {
savedWeights_= new CoinIndexedVector(rhs.savedWeights_);
} else {
savedWeights_=NULL;
}
if (rhs.dubiousWeights_) {
assert(model_);
int number = model_->numberRows();
dubiousWeights_= new int[number];
ClpDisjointCopyN(rhs.dubiousWeights_,number,dubiousWeights_);
} else {
dubiousWeights_=NULL;
}
} else {
infeasible_=NULL;
weights_=NULL;
alternateWeights_=NULL;
savedWeights_=NULL;
dubiousWeights_=NULL;
}
}
//-------------------------------------------------------------------
// Destructor
//-------------------------------------------------------------------
ClpDualRowSteepest::~ClpDualRowSteepest ()
{
delete [] weights_;
delete [] dubiousWeights_;
delete infeasible_;
delete alternateWeights_;
delete savedWeights_;
}
//----------------------------------------------------------------
// Assignment operator
//-------------------------------------------------------------------
ClpDualRowSteepest &
ClpDualRowSteepest::operator=(const ClpDualRowSteepest& rhs)
{
if (this != &rhs) {
ClpDualRowPivot::operator=(rhs);
state_=rhs.state_;
mode_ = rhs.mode_;
persistence_ = rhs.persistence_;
model_ = rhs.model_;
delete [] weights_;
delete [] dubiousWeights_;
delete infeasible_;
delete alternateWeights_;
delete savedWeights_;
assert(model_);
int number = model_->numberRows();
if (rhs.savedWeights_)
number = CoinMin(number,rhs.savedWeights_->capacity());
if (rhs.infeasible_!=NULL) {
infeasible_= new CoinIndexedVector(rhs.infeasible_);
} else {
infeasible_=NULL;
}
if (rhs.weights_!=NULL) {
weights_= new double[number];
ClpDisjointCopyN(rhs.weights_,number,weights_);
} else {
weights_=NULL;
}
if (rhs.alternateWeights_!=NULL) {
alternateWeights_= new CoinIndexedVector(rhs.alternateWeights_);
} else {
alternateWeights_=NULL;
}
if (rhs.savedWeights_!=NULL) {
savedWeights_= new CoinIndexedVector(rhs.savedWeights_);
} else {
savedWeights_=NULL;
}
if (rhs.dubiousWeights_) {
assert(model_);
int number = model_->numberRows();
dubiousWeights_= new int[number];
ClpDisjointCopyN(rhs.dubiousWeights_,number,dubiousWeights_);
} else {
dubiousWeights_=NULL;
}
}
return *this;
}
// Returns pivot row, -1 if none
int
ClpDualRowSteepest::pivotRow()
{
assert(model_);
int i,iRow;
double * infeas = infeasible_->denseVector();
double largest=0.0;
int * index = infeasible_->getIndices();
int number = infeasible_->getNumElements();
const int * pivotVariable =model_->pivotVariable();
int chosenRow=-1;
int lastPivotRow = model_->pivotRow();
double tolerance=model_->currentPrimalTolerance();
// we can't really trust infeasibilities if there is primal error
// this coding has to mimic coding in checkPrimalSolution
double error = CoinMin(1.0e-2,model_->largestPrimalError());
// allow tolerance at least slightly bigger than standard
tolerance = tolerance + error;
// But cap
tolerance = CoinMin(1000.0,tolerance);
tolerance *= tolerance; // as we are using squares
double saveTolerance = tolerance;
double * solution = model_->solutionRegion();
double * lower = model_->lowerRegion();
double * upper = model_->upperRegion();
// do last pivot row one here
//#define COLUMN_BIAS 4.0
//#define FIXED_BIAS 10.0
if (lastPivotRow>=0) {
#ifdef COLUMN_BIAS
int numberColumns = model_->numberColumns();
#endif
int iPivot=pivotVariable[lastPivotRow];
double value = solution[iPivot];
double lower = model_->lower(iPivot);
double upper = model_->upper(iPivot);
if (value>upper+tolerance) {
value -= upper;
value *= value;
#ifdef COLUMN_BIAS
if (iPivot<numberColumns)
value *= COLUMN_BIAS; // bias towards columns
k
#endif
// store square in list
if (infeas[lastPivotRow])
infeas[lastPivotRow] = value; // already there
else
infeasible_->quickAdd(lastPivotRow,value);
} else if (value<lower-tolerance) {
value -= lower;
value *= value;
#ifdef COLUMN_BIAS
if (iPivot<numberColumns)
value *= COLUMN_BIAS; // bias towards columns
#endif
// store square in list
if (infeas[lastPivotRow])
infeas[lastPivotRow] = value; // already there
else
infeasible_->add(lastPivotRow,value);
} else {
// feasible - was it infeasible - if so set tiny
if (infeas[lastPivotRow])
infeas[lastPivotRow] = COIN_INDEXED_REALLY_TINY_ELEMENT;
}
number = infeasible_->getNumElements();
}
if(model_->numberIterations()<model_->lastBadIteration()+200) {
// we can't really trust infeasibilities if there is dual error
if (model_->largestDualError()>model_->largestPrimalError())
tolerance *= CoinMin(model_->largestDualError()/model_->largestPrimalError(),1000.0);
}
int numberWanted;
if (mode_<2 ) {
numberWanted = number+1;
} else if (mode_==2) {
numberWanted = CoinMax(2000,number/8);
} else {
int numberElements = model_->factorization()->numberElements();
double ratio = (double) numberElements/(double) model_->numberRows();
numberWanted = CoinMax(2000,number/8);
if (ratio<1.0) {
numberWanted = CoinMax(2000,number/20);
} else if (ratio>10.0) {
ratio = number * (ratio/80.0);
if (ratio>number)
numberWanted=number+1;
else
numberWanted = CoinMax(2000,(int) ratio);
}
}
if (model_->largestPrimalError()>1.0e-3)
numberWanted = number+1; // be safe
int iPass;
// Setup two passes
int start[4];
start[1]=number;
start[2]=0;
double dstart = ((double) number) * CoinDrand48();
start[0]=(int) dstart;
start[3]=start[0];
//double largestWeight=0.0;
//double smallestWeight=1.0e100;
for (iPass=0;iPass<2;iPass++) {
int end = start[2*iPass+1];
for (i=start[2*iPass];i<end;i++) {
iRow = index[i];
double value = infeas[iRow];
if (value>tolerance) {
//#define OUT_EQ
#ifdef OUT_EQ
{
int iSequence = pivotVariable[iRow];
if (upper[iSequence]==lower[iSequence])
value *= 2.0;
}
#endif
double weight = CoinMin(weights_[iRow],1.0e50);
//largestWeight = CoinMax(largestWeight,weight);
//smallestWeight = CoinMin(smallestWeight,weight);
//double dubious = dubiousWeights_[iRow];
//weight *= dubious;
//if (value>2.0*largest*weight||(value>0.5*largest*weight&&value*largestWeight>dubious*largest*weight)) {
if (value>largest*weight) {
// make last pivot row last resort choice
if (iRow==lastPivotRow) {
if (value*1.0e-10<largest*weight)
continue;
else
value *= 1.0e-10;
}
int iSequence = pivotVariable[iRow];
if (!model_->flagged(iSequence)) {
//#define CLP_DEBUG 3
#ifdef CLP_DEBUG
double value2=0.0;
if (solution[iSequence]>upper[iSequence]+tolerance)
value2=solution[iSequence]-upper[iSequence];
else if (solution[iSequence]<lower[iSequence]-tolerance)
value2=solution[iSequence]-lower[iSequence];
assert(fabs(value2*value2-infeas[iRow])<1.0e-8*CoinMin(value2*value2,infeas[iRow]));
#endif
if (solution[iSequence]>upper[iSequence]+tolerance||
solution[iSequence]<lower[iSequence]-tolerance) {
chosenRow=iRow;
largest=value/weight;
//largestWeight = dubious;
}
} else {
// just to make sure we don't exit before got something
numberWanted++;
}
}
numberWanted--;
if (!numberWanted)
break;
}
}
if (!numberWanted)
break;
}
//printf("smallest %g largest %g\n",smallestWeight,largestWeight);
if (chosenRow<0&& tolerance>saveTolerance) {
// won't line up with checkPrimalSolution - do again
double saveError = model_->largestDualError();
model_->setLargestDualError(0.0);
// can't loop
chosenRow=pivotRow();
model_->setLargestDualError(saveError);
}
return chosenRow;
}
// Updates weights and returns pivot alpha
double
ClpDualRowSteepest::updateWeights(CoinIndexedVector * input,
CoinIndexedVector * spare,
CoinIndexedVector * updatedColumn)
{
#if CLP_DEBUG>2
// Very expensive debug
{
int numberRows = model_->numberRows();
CoinIndexedVector * temp = new CoinIndexedVector();
temp->reserve(numberRows+
model_->factorization()->maximumPivots());
double * array = alternateWeights_->denseVector();
int * which = alternateWeights_->getIndices();
int i;
for (i=0;i<numberRows;i++) {
double value=0.0;
array[i] = 1.0;
which[0] = i;
alternateWeights_->setNumElements(1);
model_->factorization()->updateColumnTranspose(temp,
alternateWeights_);
int number = alternateWeights_->getNumElements();
int j;
for (j=0;j<number;j++) {
int iRow=which[j];
value+=array[iRow]*array[iRow];
array[iRow]=0.0;
}
alternateWeights_->setNumElements(0);
double w = CoinMax(weights_[i],value)*.1;
if (fabs(weights_[i]-value)>w) {
printf("%d old %g, true %g\n",i,weights_[i],value);
weights_[i]=value; // to reduce printout
}
//else
//printf("%d matches %g\n",i,value);
}
delete temp;
}
#endif
assert (input->packedMode());
if (!updatedColumn->packedMode()) {
// I think this means empty
#ifdef COIN_DEVELOP
printf("updatedColumn not packed mode ClpDualRowSteepest::updateWeights\n");
#endif
return 0.0;
}
double alpha=0.0;
if (!model_->factorization()->networkBasis()) {
// clear other region
alternateWeights_->clear();
double norm = 0.0;
int i;
double * work = input->denseVector();
int numberNonZero = input->getNumElements();
int * which = input->getIndices();
double * work2 = spare->denseVector();
int * which2 = spare->getIndices();
// ftran
//permute and move indices into index array
//also compute norm
//int *regionIndex = alternateWeights_->getIndices ( );
const int *permute = model_->factorization()->permute();
//double * region = alternateWeights_->denseVector();
for ( i = 0; i < numberNonZero; i ++ ) {
int iRow = which[i];
double value = work[i];
norm += value*value;
iRow = permute[iRow];
work2[iRow] = value;
which2[i] = iRow;
}
spare->setNumElements ( numberNonZero );
// Only one array active as already permuted
model_->factorization()->updateColumn(spare,spare,true);
// permute back
numberNonZero = spare->getNumElements();
// alternateWeights_ should still be empty
int pivotRow = model_->pivotRow();
#ifdef CLP_DEBUG
if ( model_->logLevel ( ) >4 &&
fabs(norm-weights_[pivotRow])>1.0e-3*(1.0+norm))
printf("on row %d, true weight %g, old %g\n",
pivotRow,sqrt(norm),sqrt(weights_[pivotRow]));
#endif
// could re-initialize here (could be expensive)
norm /= model_->alpha() * model_->alpha();
assert(norm);
// pivot element
alpha=0.0;
double multiplier = 2.0 / model_->alpha();
// look at updated column
work = updatedColumn->denseVector();
numberNonZero = updatedColumn->getNumElements();
which = updatedColumn->getIndices();
int nSave=0;
double * work3 = alternateWeights_->denseVector();
int * which3 = alternateWeights_->getIndices();
const int * pivotColumn = model_->factorization()->pivotColumn();
for (i =0; i < numberNonZero; i++) {
int iRow = which[i];
double theta = work[i];
if (iRow==pivotRow)
alpha = theta;
double devex = weights_[iRow];
work3[nSave]=devex; // save old
which3[nSave++]=iRow;
// transform to match spare
int jRow = pivotColumn[iRow];
double value = work2[jRow];
devex += theta * (theta*norm+value * multiplier);
if (devex < DEVEX_TRY_NORM)
devex = DEVEX_TRY_NORM;
weights_[iRow]=devex;
}
alternateWeights_->setPackedMode(true);
alternateWeights_->setNumElements(nSave);
if (norm < DEVEX_TRY_NORM)
norm = DEVEX_TRY_NORM;
// Try this to make less likely will happen again and stop cycling
//norm *= 1.02;
weights_[pivotRow] = norm;
spare->clear();
#ifdef CLP_DEBUG
spare->checkClear();
#endif
} else {
// clear other region
alternateWeights_->clear();
double norm = 0.0;
int i;
double * work = input->denseVector();
int number = input->getNumElements();
int * which = input->getIndices();
double * work2 = spare->denseVector();
int * which2 = spare->getIndices();
for (i=0;i<number;i++) {
int iRow = which[i];
double value = work[i];
norm += value*value;
work2[iRow]=value;
which2[i]=iRow;
}
spare->setNumElements(number);
// ftran
model_->factorization()->updateColumn(alternateWeights_,spare);
// alternateWeights_ should still be empty
int pivotRow = model_->pivotRow();
#ifdef CLP_DEBUG
if ( model_->logLevel ( ) >4 &&
fabs(norm-weights_[pivotRow])>1.0e-3*(1.0+norm))
printf("on row %d, true weight %g, old %g\n",
pivotRow,sqrt(norm),sqrt(weights_[pivotRow]));
#endif
// could re-initialize here (could be expensive)
norm /= model_->alpha() * model_->alpha();
assert(norm);
//if (norm < DEVEX_TRY_NORM)
//norm = DEVEX_TRY_NORM;
// pivot element
alpha=0.0;
double multiplier = 2.0 / model_->alpha();
// look at updated column
work = updatedColumn->denseVector();
number = updatedColumn->getNumElements();
which = updatedColumn->getIndices();
int nSave=0;
double * work3 = alternateWeights_->denseVector();
int * which3 = alternateWeights_->getIndices();
for (i =0; i < number; i++) {
int iRow = which[i];
double theta = work[i];
if (iRow==pivotRow)
alpha = theta;
double devex = weights_[iRow];
work3[nSave]=devex; // save old
which3[nSave++]=iRow;
double value = work2[iRow];
devex += theta * (theta*norm+value * multiplier);
if (devex < DEVEX_TRY_NORM)
devex = DEVEX_TRY_NORM;
weights_[iRow]=devex;
}
if (!alpha) {
// error - but carry on
alpha =1.0e-50;
}
alternateWeights_->setPackedMode(true);
alternateWeights_->setNumElements(nSave);
if (norm < DEVEX_TRY_NORM)
norm = DEVEX_TRY_NORM;
weights_[pivotRow] = norm;
spare->clear();
}
#ifdef CLP_DEBUG
spare->checkClear();
#endif
return alpha;
}
/* Updates primal solution (and maybe list of candidates)
Uses input vector which it deletes
Computes change in objective function
*/
void
ClpDualRowSteepest::updatePrimalSolution(
CoinIndexedVector * primalUpdate,
double primalRatio,
double & objectiveChange)
{
double * work = primalUpdate->denseVector();
int number = primalUpdate->getNumElements();
int * which = primalUpdate->getIndices();
int i;
double changeObj=0.0;
double tolerance=model_->currentPrimalTolerance();
const int * pivotVariable = model_->pivotVariable();
double * infeas = infeasible_->denseVector();
int pivotRow = model_->pivotRow();
double * solution = model_->solutionRegion();
#ifdef COLUMN_BIAS
int numberColumns = model_->numberColumns();
#endif
if (primalUpdate->packedMode()) {
for (i=0;i<number;i++) {
int iRow=which[i];
int iPivot=pivotVariable[iRow];
double value = solution[iPivot];
double cost = model_->cost(iPivot);
double change = primalRatio*work[i];
work[i]=0.0;
value -= change;
changeObj -= change*cost;
solution[iPivot] = value;
double lower = model_->lower(iPivot);
double upper = model_->upper(iPivot);
// But if pivot row then use value of incoming
// Although it is safer to recompute before next selection
// in case something odd happens
if (iRow==pivotRow) {
iPivot = model_->sequenceIn();
lower = model_->lower(iPivot);
upper = model_->upper(iPivot);
value = model_->valueIncomingDual();
}
if (value<lower-tolerance) {
value -= lower;
value *= value;
#ifdef COLUMN_BIAS
if (iPivot<numberColumns)
value *= COLUMN_BIAS; // bias towards columns
#endif
#ifdef FIXED_BIAS
if (lower==upper)
value *= FIXED_BIAS; // bias towards taking out fixed variables
#endif
// store square in list
if (infeas[iRow])
infeas[iRow] = value; // already there
else
infeasible_->quickAdd(iRow,value);
} else if (value>upper+tolerance) {
value -= upper;
value *= value;
#ifdef COLUMN_BIAS
if (iPivot<numberColumns)
value *= COLUMN_BIAS; // bias towards columns
#endif
#ifdef FIXED_BIAS
if (lower==upper)
value *= FIXED_BIAS; // bias towards taking out fixed variables
#endif
// store square in list
if (infeas[iRow])
infeas[iRow] = value; // already there
else
infeasible_->quickAdd(iRow,value);
} else {
// feasible - was it infeasible - if so set tiny
if (infeas[iRow])
infeas[iRow] = COIN_INDEXED_REALLY_TINY_ELEMENT;
}
}
} else {
for (i=0;i<number;i++) {
int iRow=which[i];
int iPivot=pivotVariable[iRow];
double value = solution[iPivot];
double cost = model_->cost(iPivot);
double change = primalRatio*work[iRow];
value -= change;
changeObj -= change*cost;
solution[iPivot] = value;
double lower = model_->lower(iPivot);
double upper = model_->upper(iPivot);
// But if pivot row then use value of incoming
// Although it is safer to recompute before next selection
// in case something odd happens
if (iRow==pivotRow) {
iPivot = model_->sequenceIn();
lower = model_->lower(iPivot);
upper = model_->upper(iPivot);
value = model_->valueIncomingDual();
}
if (value<lower-tolerance) {
value -= lower;
value *= value;
#ifdef COLUMN_BIAS
if (iPivot<numberColumns)
value *= COLUMN_BIAS; // bias towards columns
#endif
#ifdef FIXED_BIAS
if (lower==upper)
value *= FIXED_BIAS; // bias towards taking out fixed variables
#endif
// store square in list
if (infeas[iRow])
infeas[iRow] = value; // already there
else
infeasible_->quickAdd(iRow,value);
} else if (value>upper+tolerance) {
value -= upper;
value *= value;
#ifdef COLUMN_BIAS
if (iPivot<numberColumns)
value *= COLUMN_BIAS; // bias towards columns
#endif
#ifdef FIXED_BIAS
if (lower==upper)
value *= FIXED_BIAS; // bias towards taking out fixed variables
#endif
// store square in list
if (infeas[iRow])
infeas[iRow] = value; // already there
else
infeasible_->quickAdd(iRow,value);
} else {
// feasible - was it infeasible - if so set tiny
if (infeas[iRow])
infeas[iRow] = COIN_INDEXED_REALLY_TINY_ELEMENT;
}
work[iRow]=0.0;
}
}
primalUpdate->setNumElements(0);
objectiveChange += changeObj;
}
/* Saves any weights round factorization as pivot rows may change
1) before factorization
2) after factorization
3) just redo infeasibilities
4) restore weights
*/
void
ClpDualRowSteepest::saveWeights(ClpSimplex * model,int mode)
{
// alternateWeights_ is defined as indexed but is treated oddly
model_ = model;
int numberRows = model_->numberRows();
int numberColumns = model_->numberColumns();
const int * pivotVariable = model_->pivotVariable();
int i;
if (mode==1) {
if(weights_) {
// Check if size has changed
if (infeasible_->capacity()==numberRows) {
alternateWeights_->clear();
// change from row numbers to sequence numbers
int * which = alternateWeights_->getIndices();
for (i=0;i<numberRows;i++) {
int iPivot=pivotVariable[i];
which[i]=iPivot;
}
state_=1;
} else {
// size has changed - clear everything
delete [] weights_;
weights_=NULL;
delete [] dubiousWeights_;
dubiousWeights_=NULL;
delete infeasible_;
infeasible_=NULL;
delete alternateWeights_;
alternateWeights_=NULL;
delete savedWeights_;
savedWeights_=NULL;
state_=-1;
}
}
} else if (mode==2||mode==4||mode==5) {
// restore
if (!weights_||state_==-1||mode==5) {
// initialize weights
delete [] weights_;
delete alternateWeights_;
weights_ = new double[numberRows];
alternateWeights_ = new CoinIndexedVector();
// enough space so can use it for factorization
alternateWeights_->reserve(numberRows+
model_->factorization()->maximumPivots());
if (mode_!=1||mode==5) {
// initialize to 1.0 (can we do better?)
for (i=0;i<numberRows;i++) {
weights_[i]=1.0;
}
} else {
CoinIndexedVector * temp = new CoinIndexedVector();
temp->reserve(numberRows+
model_->factorization()->maximumPivots());
double * array = alternateWeights_->denseVector();
int * which = alternateWeights_->getIndices();
for (i=0;i<numberRows;i++) {
double value=0.0;
array[0] = 1.0;
which[0] = i;
alternateWeights_->setNumElements(1);
alternateWeights_->setPackedMode(true);
model_->factorization()->updateColumnTranspose(temp,
alternateWeights_);
int number = alternateWeights_->getNumElements();
int j;
for (j=0;j<number;j++) {
value+=array[j]*array[j];
array[j]=0.0;
}
alternateWeights_->setNumElements(0);
weights_[i] = value;
}
delete temp;
}
// create saved weights (not really indexedvector)
savedWeights_ = new CoinIndexedVector();
savedWeights_->reserve(numberRows);
double * array = savedWeights_->denseVector();
int * which = savedWeights_->getIndices();
for (i=0;i<numberRows;i++) {
array[i]=weights_[i];
which[i]=pivotVariable[i];
}
} else {
int * which = alternateWeights_->getIndices();
CoinIndexedVector * rowArray3 = model_->rowArray(3);
rowArray3->clear();
int * back = rowArray3->getIndices();
// In case something went wrong
for (i=0;i<numberRows+numberColumns;i++)
back[i]=-1;
if (mode!=4) {
// save
memcpy(savedWeights_->getIndices(),which,
numberRows*sizeof(int));
memcpy(savedWeights_->denseVector(),weights_,
numberRows*sizeof(double));
} else {
// restore
//memcpy(which,savedWeights_->getIndices(),
// numberRows*sizeof(int));
//memcpy(weights_,savedWeights_->denseVector(),
// numberRows*sizeof(double));
which = savedWeights_->getIndices();
}
// restore (a bit slow - but only every re-factorization)
double * array = savedWeights_->denseVector();
for (i=0;i<numberRows;i++) {
int iSeq=which[i];
back[iSeq]=i;
}
for (i=0;i<numberRows;i++) {
int iPivot=pivotVariable[i];
iPivot = back[iPivot];
if (iPivot>=0) {
weights_[i]=array[iPivot];
if (weights_[i]<DEVEX_TRY_NORM)
weights_[i] = DEVEX_TRY_NORM; // may need to check more
} else {
// odd
weights_[i]=1.0;
}
}
}
state_=0;
// set up infeasibilities
if (!infeasible_) {
infeasible_ = new CoinIndexedVector();
infeasible_->reserve(numberRows);
}
}
if (mode>=2) {
// Get dubious weights
//if (!dubiousWeights_)
//dubiousWeights_=new int[numberRows];
//model_->factorization()->getWeights(dubiousWeights_);
infeasible_->clear();
int iRow;
const int * pivotVariable = model_->pivotVariable();
double tolerance=model_->currentPrimalTolerance();
for (iRow=0;iRow<numberRows;iRow++) {
int iPivot=pivotVariable[iRow];
double value = model_->solution(iPivot);
double lower = model_->lower(iPivot);
double upper = model_->upper(iPivot);
if (value<lower-tolerance) {
value -= lower;
value *= value;
#ifdef COLUMN_BIAS
if (iPivot<numberColumns)
value *= COLUMN_BIAS; // bias towards columns
#endif
#ifdef FIXED_BIAS
if (lower==upper)
value *= FIXED_BIAS; // bias towards taking out fixed variables
#endif
// store square in list
infeasible_->quickAdd(iRow,value);
} else if (value>upper+tolerance) {
value -= upper;
value *= value;
#ifdef COLUMN_BIAS
if (iPivot<numberColumns)
value *= COLUMN_BIAS; // bias towards columns
#endif
#ifdef FIXED_BIAS
if (lower==upper)
value *= FIXED_BIAS; // bias towards taking out fixed variables
#endif
// store square in list
infeasible_->quickAdd(iRow,value);
}
}
}
}
// Gets rid of last update
void
ClpDualRowSteepest::unrollWeights()
{
double * saved = alternateWeights_->denseVector();
int number = alternateWeights_->getNumElements();
int * which = alternateWeights_->getIndices();
int i;
if (alternateWeights_->packedMode()) {
for (i=0;i<number;i++) {
int iRow = which[i];
weights_[iRow]=saved[i];
saved[i]=0.0;
}
} else {
for (i=0;i<number;i++) {
int iRow = which[i];
weights_[iRow]=saved[iRow];
saved[iRow]=0.0;
}
}
alternateWeights_->setNumElements(0);
}
//-------------------------------------------------------------------
// Clone
//-------------------------------------------------------------------
ClpDualRowPivot * ClpDualRowSteepest::clone(bool CopyData) const
{
if (CopyData) {
return new ClpDualRowSteepest(*this);
} else {
return new ClpDualRowSteepest();
}
}
// Gets rid of all arrays
void
ClpDualRowSteepest::clearArrays()
{
if (persistence_==normal) {
delete [] weights_;
weights_=NULL;
delete [] dubiousWeights_;
dubiousWeights_=NULL;
delete infeasible_;
infeasible_ = NULL;
delete alternateWeights_;
alternateWeights_ = NULL;
delete savedWeights_;
savedWeights_ = NULL;
}
state_ =-1;
}
// Returns true if would not find any row
bool
ClpDualRowSteepest::looksOptimal() const
{
int iRow;
const int * pivotVariable = model_->pivotVariable();
double tolerance=model_->currentPrimalTolerance();
// we can't really trust infeasibilities if there is primal error
// this coding has to mimic coding in checkPrimalSolution
double error = CoinMin(1.0e-2,model_->largestPrimalError());
// allow tolerance at least slightly bigger than standard
tolerance = tolerance + error;
// But cap
tolerance = CoinMin(1000.0,tolerance);
int numberRows = model_->numberRows();
int numberInfeasible=0;
for (iRow=0;iRow<numberRows;iRow++) {
int iPivot=pivotVariable[iRow];
double value = model_->solution(iPivot);
double lower = model_->lower(iPivot);
double upper = model_->upper(iPivot);
if (value<lower-tolerance) {
numberInfeasible++;
} else if (value>upper+tolerance) {
numberInfeasible++;
}
}
return (numberInfeasible==0);
}
// Called when maximum pivots changes
void
ClpDualRowSteepest::maximumPivotsChanged()
{
if (alternateWeights_&&
alternateWeights_->capacity()!=model_->numberRows()+
model_->factorization()->maximumPivots()) {
delete alternateWeights_;
alternateWeights_ = new CoinIndexedVector();
// enough space so can use it for factorization
alternateWeights_->reserve(model_->numberRows()+
model_->factorization()->maximumPivots());
}
}
syntax highlighted by Code2HTML, v. 0.9.1