// Copyright (C) 2004, International Business Machines
// Corporation and others. All Rights Reserved.
#include <cstdio>
#include "CoinPragma.hpp"
#include "CoinIndexedVector.hpp"
#include "CoinHelperFunctions.hpp"
#include "ClpSimplex.hpp"
#include "ClpFactorization.hpp"
#include "ClpQuadraticObjective.hpp"
#include "ClpNonLinearCost.hpp"
// at end to get min/max!
#include "ClpDynamicExampleMatrix.hpp"
#include "ClpMessage.hpp"
//#define CLP_DEBUG
//#define CLP_DEBUG_PRINT
//#############################################################################
// Constructors / Destructor / Assignment
//#############################################################################
//-------------------------------------------------------------------
// Default Constructor
//-------------------------------------------------------------------
ClpDynamicExampleMatrix::ClpDynamicExampleMatrix ()
: ClpDynamicMatrix(),
numberColumns_(0),
startColumnGen_(NULL),
rowGen_(NULL),
elementGen_(NULL),
costGen_(NULL),
fullStartGen_(NULL),
dynamicStatusGen_(NULL),
idGen_(NULL),
columnLowerGen_(NULL),
columnUpperGen_(NULL)
{
setType(25);
}
//-------------------------------------------------------------------
// Copy constructor
//-------------------------------------------------------------------
ClpDynamicExampleMatrix::ClpDynamicExampleMatrix (const ClpDynamicExampleMatrix & rhs)
: ClpDynamicMatrix(rhs)
{
numberColumns_ = rhs.numberColumns_;
startColumnGen_ = ClpCopyOfArray(rhs.startColumnGen_,numberColumns_+1);
CoinBigIndex numberElements = startColumnGen_[numberColumns_];
rowGen_ = ClpCopyOfArray(rhs.rowGen_,numberElements);;
elementGen_ = ClpCopyOfArray(rhs.elementGen_,numberElements);;
costGen_ = ClpCopyOfArray(rhs.costGen_,numberColumns_);
fullStartGen_ = ClpCopyOfArray(rhs.fullStartGen_,numberSets_+1);
dynamicStatusGen_ = ClpCopyOfArray(rhs.dynamicStatusGen_,numberColumns_);
idGen_ = ClpCopyOfArray(rhs.idGen_,maximumGubColumns_);
columnLowerGen_ = ClpCopyOfArray(rhs.columnLowerGen_,numberColumns_);
columnUpperGen_ = ClpCopyOfArray(rhs.columnUpperGen_,numberColumns_);
}
/* This is the real constructor*/
ClpDynamicExampleMatrix::ClpDynamicExampleMatrix(ClpSimplex * model, int numberSets,
int numberGubColumns, const int * starts,
const double * lower, const double * upper,
const CoinBigIndex * startColumn, const int * row,
const double * element, const double * cost,
const double * columnLower, const double * columnUpper,
const unsigned char * status,
const unsigned char * dynamicStatus,
int numberIds,const int *ids)
: ClpDynamicMatrix(model,numberSets,0,NULL,lower,upper,NULL,NULL,NULL,NULL,NULL,NULL,
NULL,NULL)
{
setType(25);
numberColumns_ = numberGubColumns;
// start with safe values - then experiment
maximumGubColumns_=numberColumns_;
maximumElements_ = startColumn[numberColumns_];
// delete odd stuff created by ClpDynamicMatrix constructor
delete [] startSet_;
startSet_ = new int [numberSets_];
delete [] next_;
next_ = new int [maximumGubColumns_];
delete [] row_;
delete [] element_;
delete [] startColumn_;
delete [] cost_;
delete [] columnLower_;
delete [] columnUpper_;
delete [] dynamicStatus_;
delete [] status_;
delete [] id_;
// and size correctly
row_ = new int [maximumElements_];
element_ = new double [maximumElements_];
startColumn_ = new CoinBigIndex [maximumGubColumns_+1];
// say no columns yet
numberGubColumns_=0;
startColumn_[0]=0;
cost_ = new double[maximumGubColumns_];
dynamicStatus_ = new unsigned char [maximumGubColumns_];
memset(dynamicStatus_,0,maximumGubColumns_);
id_ = new int[maximumGubColumns_];
if (columnLower)
columnLower_ = new double[maximumGubColumns_];
else
columnLower_ = NULL;
if (columnUpper)
columnUpper_ = new double[maximumGubColumns_];
else
columnUpper_=NULL;
// space for ids
idGen_ = new int [maximumGubColumns_];
int iSet;
for (iSet=0;iSet<numberSets_;iSet++)
startSet_[iSet]=-1;
// This starts code specific to this storage method
CoinBigIndex i;
fullStartGen_ = ClpCopyOfArray(starts,numberSets_+1);
startColumnGen_ = ClpCopyOfArray(startColumn,numberColumns_+1);
CoinBigIndex numberElements = startColumnGen_[numberColumns_];
rowGen_ = ClpCopyOfArray(row,numberElements);
elementGen_ = new double[numberElements];
for (i=0;i<numberElements;i++)
elementGen_[i]=element[i];
costGen_ = new double[numberColumns_];
for (i=0;i<numberColumns_;i++) {
costGen_[i]=cost[i];
// I don't think I need sorted but ...
CoinSort_2(rowGen_+startColumnGen_[i],rowGen_+startColumnGen_[i+1],elementGen_+startColumnGen_[i]);
}
if (columnLower) {
columnLowerGen_ = new double[numberColumns_];
for (i=0;i<numberColumns_;i++) {
columnLowerGen_[i]=columnLower[i];
if (columnLowerGen_[i]) {
printf("Non-zero lower bounds not allowed - subtract from model\n");
abort();
}
}
} else {
columnLowerGen_=NULL;
}
if (columnUpper) {
columnUpperGen_ = new double[numberColumns_];
for (i=0;i<numberColumns_;i++)
columnUpperGen_[i]=columnUpper[i];
} else {
columnUpperGen_=NULL;
}
// end specific coding
if (columnUpper_) {
// set all upper bounds so we have enough space
double * columnUpper = model->columnUpper();
for(i=firstDynamic_;i<lastDynamic_;i++)
columnUpper[i]=1.0e10;
}
if (status) {
status_ = ClpCopyOfArray(status,numberSets_);
assert (dynamicStatus);
memcpy(dynamicStatus_,dynamicStatus,numberIds);
assert (numberIds);
} else {
assert (!numberIds);
status_= new unsigned char [numberSets_];
memset(status_,0,numberSets_);
for (i=0;i<numberSets_;i++) {
// make slack key
setStatus(i,ClpSimplex::basic);
}
}
dynamicStatusGen_ = new unsigned char [numberColumns_];
memset(dynamicStatusGen_,0,numberColumns_); // for clarity
for (i=0;i<numberColumns_;i++)
setDynamicStatusGen(i,atLowerBound);
// Populate with enough columns
if (!numberIds) {
// This could be made more sophisticated
for (iSet=0;iSet<numberSets_;iSet++) {
int sequence = fullStartGen_[iSet];
CoinBigIndex start = startColumnGen_[sequence];
addColumn(startColumnGen_[sequence+1]-start,
rowGen_+start,
elementGen_+start,
costGen_[sequence],
columnLowerGen_ ? columnLowerGen_[sequence] : 0,
columnUpperGen_ ? columnUpperGen_[sequence] : 1.0e30,
iSet,getDynamicStatusGen(sequence));
idGen_[iSet]=sequence; // say which one in
setDynamicStatusGen(sequence,inSmall);
}
} else {
// put back old ones
int * set = new int[numberColumns_];
for (iSet=0;iSet<numberSets_;iSet++) {
for (CoinBigIndex j=fullStartGen_[iSet];j<fullStartGen_[iSet+1];j++)
set[j]=iSet;
}
for (int i=0;i<numberIds;i++) {
int sequence = ids[i];
CoinBigIndex start = startColumnGen_[sequence];
addColumn(startColumnGen_[sequence+1]-start,
rowGen_+start,
elementGen_+start,
costGen_[sequence],
columnLowerGen_ ? columnLowerGen_[sequence] : 0,
columnUpperGen_ ? columnUpperGen_[sequence] : 1.0e30,
set[sequence],getDynamicStatus(i));
idGen_[iSet]=sequence; // say which one in
setDynamicStatusGen(sequence,inSmall);
}
delete [] set;
}
if (!status) {
gubCrash();
} else {
initialProblem();
}
}
// This constructor just takes over ownership
ClpDynamicExampleMatrix::ClpDynamicExampleMatrix(ClpSimplex * model, int numberSets,
int numberGubColumns, int * starts,
const double * lower, const double * upper,
int * startColumn, int * row,
double * element, double * cost,
double * columnLower, double * columnUpper,
const unsigned char * status,
const unsigned char * dynamicStatus,
int numberIds,const int *ids)
: ClpDynamicMatrix(model,numberSets,0,NULL,lower,upper,NULL,NULL,NULL,NULL,NULL,NULL,
NULL,NULL)
{
setType(25);
numberColumns_ = numberGubColumns;
// start with safe values - then experiment
maximumGubColumns_=numberColumns_;
maximumElements_ = startColumn[numberColumns_];
// delete odd stuff created by ClpDynamicMatrix constructor
delete [] startSet_;
startSet_ = new int [numberSets_];
delete [] next_;
next_ = new int [maximumGubColumns_];
delete [] row_;
delete [] element_;
delete [] startColumn_;
delete [] cost_;
delete [] columnLower_;
delete [] columnUpper_;
delete [] dynamicStatus_;
delete [] status_;
delete [] id_;
// and size correctly
row_ = new int [maximumElements_];
element_ = new double [maximumElements_];
startColumn_ = new CoinBigIndex [maximumGubColumns_+1];
// say no columns yet
numberGubColumns_=0;
startColumn_[0]=0;
cost_ = new double[maximumGubColumns_];
dynamicStatus_ = new unsigned char [maximumGubColumns_];
memset(dynamicStatus_,0,maximumGubColumns_);
id_ = new int[maximumGubColumns_];
if (columnLower)
columnLower_ = new double[maximumGubColumns_];
else
columnLower_ = NULL;
if (columnUpper)
columnUpper_ = new double[maximumGubColumns_];
else
columnUpper_=NULL;
// space for ids
idGen_ = new int [maximumGubColumns_];
int iSet;
for (iSet=0;iSet<numberSets_;iSet++)
startSet_[iSet]=-1;
// This starts code specific to this storage method
CoinBigIndex i;
fullStartGen_ = starts;
startColumnGen_ = startColumn;
rowGen_ = row;
elementGen_ = element;
costGen_ = cost;
for (i=0;i<numberColumns_;i++) {
// I don't think I need sorted but ...
CoinSort_2(rowGen_+startColumnGen_[i],rowGen_+startColumnGen_[i+1],elementGen_+startColumnGen_[i]);
}
if (columnLower) {
columnLowerGen_ = columnLower;
for (i=0;i<numberColumns_;i++) {
if (columnLowerGen_[i]) {
printf("Non-zero lower bounds not allowed - subtract from model\n");
abort();
}
}
} else {
columnLowerGen_=NULL;
}
if (columnUpper) {
columnUpperGen_ = columnUpper;
} else {
columnUpperGen_=NULL;
}
// end specific coding
if (columnUpper_) {
// set all upper bounds so we have enough space
double * columnUpper = model->columnUpper();
for(i=firstDynamic_;i<lastDynamic_;i++)
columnUpper[i]=1.0e10;
}
if (status) {
status_ = ClpCopyOfArray(status,numberSets_);
assert (dynamicStatus);
memcpy(dynamicStatus_,dynamicStatus,numberIds);
assert (numberIds);
} else {
assert (!numberIds);
status_= new unsigned char [numberSets_];
memset(status_,0,numberSets_);
for (i=0;i<numberSets_;i++) {
// make slack key
setStatus(i,ClpSimplex::basic);
}
}
dynamicStatusGen_ = new unsigned char [numberColumns_];
memset(dynamicStatusGen_,0,numberColumns_); // for clarity
for (i=0;i<numberColumns_;i++)
setDynamicStatusGen(i,atLowerBound);
// Populate with enough columns
if (!numberIds) {
// This could be made more sophisticated
for (iSet=0;iSet<numberSets_;iSet++) {
int sequence = fullStartGen_[iSet];
CoinBigIndex start = startColumnGen_[sequence];
addColumn(startColumnGen_[sequence+1]-start,
rowGen_+start,
elementGen_+start,
costGen_[sequence],
columnLowerGen_ ? columnLowerGen_[sequence] : 0,
columnUpperGen_ ? columnUpperGen_[sequence] : 1.0e30,
iSet,getDynamicStatusGen(sequence));
idGen_[iSet]=sequence; // say which one in
setDynamicStatusGen(sequence,inSmall);
}
} else {
// put back old ones
int * set = new int[numberColumns_];
for (iSet=0;iSet<numberSets_;iSet++) {
for (CoinBigIndex j=fullStartGen_[iSet];j<fullStartGen_[iSet+1];j++)
set[j]=iSet;
}
for (int i=0;i<numberIds;i++) {
int sequence = ids[i];
int iSet = set[sequence];
CoinBigIndex start = startColumnGen_[sequence];
addColumn(startColumnGen_[sequence+1]-start,
rowGen_+start,
elementGen_+start,
costGen_[sequence],
columnLowerGen_ ? columnLowerGen_[sequence] : 0,
columnUpperGen_ ? columnUpperGen_[sequence] : 1.0e30,
iSet,getDynamicStatus(i));
idGen_[i]=sequence; // say which one in
setDynamicStatusGen(sequence,inSmall);
}
delete [] set;
}
if (!status) {
gubCrash();
} else {
initialProblem();
}
}
//-------------------------------------------------------------------
// Destructor
//-------------------------------------------------------------------
ClpDynamicExampleMatrix::~ClpDynamicExampleMatrix ()
{
delete [] startColumnGen_;
delete [] rowGen_;
delete [] elementGen_;
delete [] costGen_;
delete [] fullStartGen_;
delete [] dynamicStatusGen_;
delete [] idGen_;
delete [] columnLowerGen_;
delete [] columnUpperGen_;
}
//----------------------------------------------------------------
// Assignment operator
//-------------------------------------------------------------------
ClpDynamicExampleMatrix &
ClpDynamicExampleMatrix::operator=(const ClpDynamicExampleMatrix& rhs)
{
if (this != &rhs) {
ClpDynamicMatrix::operator=(rhs);
numberColumns_ = rhs.numberColumns_;
delete [] startColumnGen_;
delete [] rowGen_;
delete [] elementGen_;
delete [] costGen_;
delete [] fullStartGen_;
delete [] dynamicStatusGen_;
delete [] idGen_;
delete [] columnLowerGen_;
delete [] columnUpperGen_;
startColumnGen_ = ClpCopyOfArray(rhs.startColumnGen_,numberColumns_+1);
CoinBigIndex numberElements = startColumnGen_[numberColumns_];
rowGen_ = ClpCopyOfArray(rhs.rowGen_,numberElements);;
elementGen_ = ClpCopyOfArray(rhs.elementGen_,numberElements);;
costGen_ = ClpCopyOfArray(rhs.costGen_,numberColumns_);
fullStartGen_ = ClpCopyOfArray(rhs.fullStartGen_,numberSets_+1);
dynamicStatusGen_ = ClpCopyOfArray(rhs.dynamicStatusGen_,numberColumns_);
idGen_ = ClpCopyOfArray(rhs.idGen_,maximumGubColumns_);
columnLowerGen_ = ClpCopyOfArray(rhs.columnLowerGen_,numberColumns_);
columnUpperGen_ = ClpCopyOfArray(rhs.columnUpperGen_,numberColumns_);
}
return *this;
}
//-------------------------------------------------------------------
// Clone
//-------------------------------------------------------------------
ClpMatrixBase * ClpDynamicExampleMatrix::clone() const
{
return new ClpDynamicExampleMatrix(*this);
}
// Partial pricing
void
ClpDynamicExampleMatrix::partialPricing(ClpSimplex * model, double startFraction, double endFraction,
int & bestSequence, int & numberWanted)
{
numberWanted=currentWanted_;
assert(!model->rowScale());
if (!numberSets_) {
// no gub
ClpPackedMatrix::partialPricing(model,startFraction,endFraction,bestSequence,numberWanted);
} else {
// and do some proportion of full set
int startG2 = (int) (startFraction*numberSets_);
int endG2 = (int) (endFraction*numberSets_+0.1);
endG2 = CoinMin(endG2,numberSets_);
//printf("gub price - set start %d end %d\n",
// startG2,endG2);
double tolerance=model->currentDualTolerance();
double * reducedCost = model->djRegion();
const double * duals = model->dualRowSolution();
double bestDj;
int numberRows = model->numberRows();
int slackOffset = lastDynamic_+numberRows;
int structuralOffset = slackOffset+numberSets_;
int structuralOffset2 = structuralOffset+maximumGubColumns_;
// If nothing found yet can go all the way to end
int endAll = endG2;
if (bestSequence<0&&!startG2)
endAll = numberSets_;
if (bestSequence>=0) {
if (bestSequence!=savedBestSequence_)
bestDj = fabs(reducedCost[bestSequence]); // dj from slacks or permanent
else
bestDj = savedBestDj_;
} else {
bestDj=tolerance;
}
int saveSequence = bestSequence;
double djMod=0.0;
double bestDjMod=0.0;
//printf("iteration %d start %d end %d - wanted %d\n",model->numberIterations(),
// startG2,endG2,numberWanted);
int bestSet=-1;
int minSet = minimumObjectsScan_<0 ? 5 : minimumObjectsScan_;
int minNeg = minimumGoodReducedCosts_<0 ? 5 : minimumGoodReducedCosts_;
for (int iSet = startG2;iSet<endAll;iSet++) {
if (numberWanted+minNeg<originalWanted_&&iSet>startG2+minSet) {
// give up
numberWanted=0;
break;
} else if (iSet==endG2&&bestSequence>=0) {
break;
}
int gubRow = toIndex_[iSet];
if (gubRow>=0) {
djMod = duals[gubRow+numberStaticRows_]; // have I got sign right?
} else {
int iBasic = keyVariable_[iSet];
if (iBasic>=numberColumns_) {
djMod = 0.0; // set not in
} else {
// get dj without
djMod=0.0;
for (CoinBigIndex j=startColumn_[iBasic];
j<startColumn_[iBasic+1];j++) {
int jRow=row_[j];
djMod -= duals[jRow]*element_[j];
}
djMod += cost_[iBasic];
// See if gub slack possible - dj is djMod
if (getStatus(iSet)==ClpSimplex::atLowerBound) {
double value = -djMod;
if (value>tolerance) {
numberWanted--;
if (value>bestDj) {
// check flagged variable and correct dj
if (!flagged(iSet)) {
bestDj=value;
bestSequence = slackOffset+iSet;
bestDjMod = djMod;
bestSet = iSet;
} else {
// just to make sure we don't exit before got something
numberWanted++;
abort();
}
}
}
} else if (getStatus(iSet)==ClpSimplex::atUpperBound) {
double value = djMod;
if (value>tolerance) {
numberWanted--;
if (value>bestDj) {
// check flagged variable and correct dj
if (!flagged(iSet)) {
bestDj=value;
bestSequence = slackOffset+iSet;
bestDjMod = djMod;
bestSet = iSet;
} else {
// just to make sure we don't exit before got something
numberWanted++;
abort();
}
}
}
}
}
}
// do ones in small
int iSequence= startSet_[iSet];
while (iSequence>=0) {
DynamicStatus status = getDynamicStatus(iSequence);
if (status==atLowerBound||status==atUpperBound) {
double value=cost_[iSequence]-djMod;
for (CoinBigIndex j=startColumn_[iSequence];
j<startColumn_[iSequence+1];j++) {
int jRow=row_[j];
value -= duals[jRow]*element_[j];
}
// change sign if at lower bound
if (status==atLowerBound)
value = -value;
if (value>tolerance) {
numberWanted--;
if (value>bestDj) {
// check flagged variable and correct dj
if (!flagged(iSequence)) {
bestDj=value;
bestSequence = structuralOffset+iSequence;
bestDjMod = djMod;
bestSet = iSet;
} else {
// just to make sure we don't exit before got something
numberWanted++;
}
}
}
}
iSequence = next_[iSequence]; //onto next in set
}
// and now get best by column generation
// If no upper bounds we may not need status test
for (iSequence=fullStartGen_[iSet];iSequence<fullStartGen_[iSet+1];iSequence++) {
DynamicStatus status = getDynamicStatusGen(iSequence);
assert (status!=atUpperBound&&status!=soloKey);
if (status==atLowerBound) {
double value=costGen_[iSequence]-djMod;
for (CoinBigIndex j=startColumnGen_[iSequence];
j<startColumnGen_[iSequence+1];j++) {
int jRow=rowGen_[j];
value -= duals[jRow]*elementGen_[j];
}
// change sign as at lower bound
value = -value;
if (value>tolerance) {
numberWanted--;
if (value>bestDj) {
// check flagged variable and correct dj
if (!flaggedGen(iSequence)) {
bestDj=value;
bestSequence = structuralOffset2+iSequence;
bestDjMod = djMod;
bestSet = iSet;
} else {
// just to make sure we don't exit before got something
numberWanted++;
}
}
}
}
}
if (numberWanted<=0) {
numberWanted=0;
break;
}
}
if (bestSequence!=saveSequence) {
savedBestGubDual_ = bestDjMod;
savedBestDj_ = bestDj;
savedBestSequence_ = bestSequence;
savedBestSet_ = bestSet;
}
// Do packed part before gub
// always???
// Resize so just do to gub
numberActiveColumns_=firstDynamic_;
int saveMinNeg=minimumGoodReducedCosts_;
if (bestSequence>=0)
minimumGoodReducedCosts_=-2;
currentWanted_=numberWanted;
ClpPackedMatrix::partialPricing(model,startFraction,endFraction,bestSequence,numberWanted);
numberActiveColumns_=matrix_->getNumCols();
minimumGoodReducedCosts_=saveMinNeg;
// See if may be finished
if (!startG2&&bestSequence<0)
infeasibilityWeight_=model_->infeasibilityCost();
else if (bestSequence>=0)
infeasibilityWeight_=-1.0;
currentWanted_=numberWanted;
}
}
/* Creates a variable. This is called after partial pricing and may modify matrix.
May update bestSequence.
*/
void
ClpDynamicExampleMatrix::createVariable(ClpSimplex * model, int & bestSequence)
{
int numberRows = model->numberRows();
int slackOffset = lastDynamic_+numberRows;
int structuralOffset = slackOffset+numberSets_;
int bestSequence2=savedBestSequence_-structuralOffset;
if (bestSequence2>=0) {
// See if needs new
if (bestSequence2>=maximumGubColumns_) {
bestSequence2 -= maximumGubColumns_;
int sequence = addColumn(startColumnGen_[bestSequence2+1]-startColumnGen_[bestSequence2],
rowGen_+startColumnGen_[bestSequence2],
elementGen_+startColumnGen_[bestSequence2],
costGen_[bestSequence2],
columnLowerGen_ ? columnLowerGen_[bestSequence2] : 0,
columnUpperGen_ ? columnUpperGen_[bestSequence2] : 1.0e30,
savedBestSet_,getDynamicStatusGen(bestSequence2));
savedBestSequence_ = structuralOffset + sequence;
idGen_[sequence]=bestSequence2;
setDynamicStatusGen(bestSequence2,inSmall);
}
}
ClpDynamicMatrix::createVariable(model,bestSequence);
// clear for next iteration
savedBestSequence_=-1;
}
/* If addColumn forces compression then this allows descendant to know what to do.
If >=0 then entry stayed in, if -1 then entry went out to lower bound.of zero.
Entries at upper bound (really nonzero) never go out (at present).
*/
void
ClpDynamicExampleMatrix::packDown(const int * in, int numberToPack)
{
int put=0;
for (int i=0;i<numberToPack;i++) {
int id = idGen_[i];
if (in[i]>=0) {
// stays in
assert (put==in[i]); // true for now
idGen_[put++]=id;
} else {
// out to lower bound
setDynamicStatusGen(id,atLowerBound);
}
}
assert (put==numberGubColumns_);
}
syntax highlighted by Code2HTML, v. 0.9.1