limp-cbc-0.3.2.0: cbits/coin/ClpQuadraticObjective.cpp
/* $Id: ClpQuadraticObjective.cpp 1723 2011-04-17 15:07:10Z forrest $ */
// Copyright (C) 2003, International Business Machines
// Corporation and others. All Rights Reserved.
// This code is licensed under the terms of the Eclipse Public License (EPL).
#include "CoinPragma.hpp"
#include "CoinHelperFunctions.hpp"
#include "CoinIndexedVector.hpp"
#include "ClpFactorization.hpp"
#include "ClpSimplex.hpp"
#include "ClpQuadraticObjective.hpp"
//#############################################################################
// Constructors / Destructor / Assignment
//#############################################################################
//-------------------------------------------------------------------
// Default Constructor
//-------------------------------------------------------------------
ClpQuadraticObjective::ClpQuadraticObjective ()
: ClpObjective()
{
type_ = 2;
objective_ = NULL;
quadraticObjective_ = NULL;
gradient_ = NULL;
numberColumns_ = 0;
numberExtendedColumns_ = 0;
activated_ = 0;
fullMatrix_ = false;
}
//-------------------------------------------------------------------
// Useful Constructor
//-------------------------------------------------------------------
ClpQuadraticObjective::ClpQuadraticObjective (const double * objective ,
int numberColumns,
const CoinBigIndex * start,
const int * column, const double * element,
int numberExtendedColumns)
: ClpObjective()
{
type_ = 2;
numberColumns_ = numberColumns;
if (numberExtendedColumns >= 0)
numberExtendedColumns_ = CoinMax(numberColumns_, numberExtendedColumns);
else
numberExtendedColumns_ = numberColumns_;
if (objective) {
objective_ = new double [numberExtendedColumns_];
CoinMemcpyN(objective, numberColumns_, objective_);
memset(objective_ + numberColumns_, 0, (numberExtendedColumns_ - numberColumns_)*sizeof(double));
} else {
objective_ = new double [numberExtendedColumns_];
memset(objective_, 0, numberExtendedColumns_ * sizeof(double));
}
if (start)
quadraticObjective_ = new CoinPackedMatrix(true, numberColumns, numberColumns,
start[numberColumns], element, column, start, NULL);
else
quadraticObjective_ = NULL;
gradient_ = NULL;
activated_ = 1;
fullMatrix_ = false;
}
//-------------------------------------------------------------------
// Copy constructor
//-------------------------------------------------------------------
ClpQuadraticObjective::ClpQuadraticObjective (const ClpQuadraticObjective & rhs,
int type)
: ClpObjective(rhs)
{
numberColumns_ = rhs.numberColumns_;
numberExtendedColumns_ = rhs.numberExtendedColumns_;
fullMatrix_ = rhs.fullMatrix_;
if (rhs.objective_) {
objective_ = new double [numberExtendedColumns_];
CoinMemcpyN(rhs.objective_, numberExtendedColumns_, objective_);
} else {
objective_ = NULL;
}
if (rhs.gradient_) {
gradient_ = new double [numberExtendedColumns_];
CoinMemcpyN(rhs.gradient_, numberExtendedColumns_, gradient_);
} else {
gradient_ = NULL;
}
if (rhs.quadraticObjective_) {
// see what type of matrix wanted
if (type == 0) {
// just copy
quadraticObjective_ = new CoinPackedMatrix(*rhs.quadraticObjective_);
} else if (type == 1) {
// expand to full symmetric
fullMatrix_ = true;
const int * columnQuadratic1 = rhs.quadraticObjective_->getIndices();
const CoinBigIndex * columnQuadraticStart1 = rhs.quadraticObjective_->getVectorStarts();
const int * columnQuadraticLength1 = rhs.quadraticObjective_->getVectorLengths();
const double * quadraticElement1 = rhs.quadraticObjective_->getElements();
CoinBigIndex * columnQuadraticStart2 = new CoinBigIndex [numberExtendedColumns_+1];
int * columnQuadraticLength2 = new int [numberExtendedColumns_];
int iColumn;
int numberColumns = rhs.quadraticObjective_->getNumCols();
int numberBelow = 0;
int numberAbove = 0;
int numberDiagonal = 0;
CoinZeroN(columnQuadraticLength2, numberExtendedColumns_);
for (iColumn = 0; iColumn < numberColumns; iColumn++) {
for (CoinBigIndex j = columnQuadraticStart1[iColumn];
j < columnQuadraticStart1[iColumn] + columnQuadraticLength1[iColumn]; j++) {
int jColumn = columnQuadratic1[j];
if (jColumn > iColumn) {
numberBelow++;
columnQuadraticLength2[jColumn]++;
columnQuadraticLength2[iColumn]++;
} else if (jColumn == iColumn) {
numberDiagonal++;
columnQuadraticLength2[iColumn]++;
} else {
numberAbove++;
}
}
}
if (numberAbove > 0) {
if (numberAbove == numberBelow) {
// already done
quadraticObjective_ = new CoinPackedMatrix(*rhs.quadraticObjective_);
delete [] columnQuadraticStart2;
delete [] columnQuadraticLength2;
} else {
printf("number above = %d, number below = %d, error\n",
numberAbove, numberBelow);
abort();
}
} else {
int numberElements = numberDiagonal + 2 * numberBelow;
int * columnQuadratic2 = new int [numberElements];
double * quadraticElement2 = new double [numberElements];
columnQuadraticStart2[0] = 0;
numberElements = 0;
for (iColumn = 0; iColumn < numberColumns; iColumn++) {
int n = columnQuadraticLength2[iColumn];
columnQuadraticLength2[iColumn] = 0;
numberElements += n;
columnQuadraticStart2[iColumn+1] = numberElements;
}
for (iColumn = 0; iColumn < numberColumns; iColumn++) {
for (CoinBigIndex j = columnQuadraticStart1[iColumn];
j < columnQuadraticStart1[iColumn] + columnQuadraticLength1[iColumn]; j++) {
int jColumn = columnQuadratic1[j];
if (jColumn > iColumn) {
// put in two places
CoinBigIndex put = columnQuadraticLength2[jColumn] + columnQuadraticStart2[jColumn];
columnQuadraticLength2[jColumn]++;
quadraticElement2[put] = quadraticElement1[j];
columnQuadratic2[put] = iColumn;
put = columnQuadraticLength2[iColumn] + columnQuadraticStart2[iColumn];
columnQuadraticLength2[iColumn]++;
quadraticElement2[put] = quadraticElement1[j];
columnQuadratic2[put] = jColumn;
} else if (jColumn == iColumn) {
CoinBigIndex put = columnQuadraticLength2[iColumn] + columnQuadraticStart2[iColumn];
columnQuadraticLength2[iColumn]++;
quadraticElement2[put] = quadraticElement1[j];
columnQuadratic2[put] = iColumn;
} else {
abort();
}
}
}
// Now create
quadraticObjective_ =
new CoinPackedMatrix (true,
rhs.numberExtendedColumns_,
rhs.numberExtendedColumns_,
numberElements,
quadraticElement2,
columnQuadratic2,
columnQuadraticStart2,
columnQuadraticLength2, 0.0, 0.0);
delete [] columnQuadraticStart2;
delete [] columnQuadraticLength2;
delete [] columnQuadratic2;
delete [] quadraticElement2;
}
} else {
fullMatrix_ = false;
abort(); // code when needed
}
} else {
quadraticObjective_ = NULL;
}
}
/* Subset constructor. Duplicates are allowed
and order is as given.
*/
ClpQuadraticObjective::ClpQuadraticObjective (const ClpQuadraticObjective &rhs,
int numberColumns,
const int * whichColumn)
: ClpObjective(rhs)
{
fullMatrix_ = rhs.fullMatrix_;
objective_ = NULL;
int extra = rhs.numberExtendedColumns_ - rhs.numberColumns_;
numberColumns_ = 0;
numberExtendedColumns_ = numberColumns + extra;
if (numberColumns > 0) {
// check valid lists
int numberBad = 0;
int i;
for (i = 0; i < numberColumns; i++)
if (whichColumn[i] < 0 || whichColumn[i] >= rhs.numberColumns_)
numberBad++;
if (numberBad)
throw CoinError("bad column list", "subset constructor",
"ClpQuadraticObjective");
numberColumns_ = numberColumns;
objective_ = new double[numberExtendedColumns_];
for (i = 0; i < numberColumns_; i++)
objective_[i] = rhs.objective_[whichColumn[i]];
CoinMemcpyN(rhs.objective_ + rhs.numberColumns_, (numberExtendedColumns_ - numberColumns_),
objective_ + numberColumns_);
if (rhs.gradient_) {
gradient_ = new double[numberExtendedColumns_];
for (i = 0; i < numberColumns_; i++)
gradient_[i] = rhs.gradient_[whichColumn[i]];
CoinMemcpyN(rhs.gradient_ + rhs.numberColumns_, (numberExtendedColumns_ - numberColumns_),
gradient_ + numberColumns_);
} else {
gradient_ = NULL;
}
} else {
gradient_ = NULL;
objective_ = NULL;
}
if (rhs.quadraticObjective_) {
quadraticObjective_ = new CoinPackedMatrix(*rhs.quadraticObjective_,
numberColumns, whichColumn,
numberColumns, whichColumn);
} else {
quadraticObjective_ = NULL;
}
}
//-------------------------------------------------------------------
// Destructor
//-------------------------------------------------------------------
ClpQuadraticObjective::~ClpQuadraticObjective ()
{
delete [] objective_;
delete [] gradient_;
delete quadraticObjective_;
}
//----------------------------------------------------------------
// Assignment operator
//-------------------------------------------------------------------
ClpQuadraticObjective &
ClpQuadraticObjective::operator=(const ClpQuadraticObjective& rhs)
{
if (this != &rhs) {
fullMatrix_ = rhs.fullMatrix_;
delete quadraticObjective_;
quadraticObjective_ = NULL;
delete [] objective_;
delete [] gradient_;
ClpObjective::operator=(rhs);
numberColumns_ = rhs.numberColumns_;
numberExtendedColumns_ = rhs.numberExtendedColumns_;
if (rhs.objective_) {
objective_ = new double [numberExtendedColumns_];
CoinMemcpyN(rhs.objective_, numberExtendedColumns_, objective_);
} else {
objective_ = NULL;
}
if (rhs.gradient_) {
gradient_ = new double [numberExtendedColumns_];
CoinMemcpyN(rhs.gradient_, numberExtendedColumns_, gradient_);
} else {
gradient_ = NULL;
}
if (rhs.quadraticObjective_) {
quadraticObjective_ = new CoinPackedMatrix(*rhs.quadraticObjective_);
} else {
quadraticObjective_ = NULL;
}
}
return *this;
}
// Returns gradient
double *
ClpQuadraticObjective::gradient(const ClpSimplex * model,
const double * solution, double & offset, bool refresh,
int includeLinear)
{
offset = 0.0;
bool scaling = false;
if (model && (model->rowScale() ||
model->objectiveScale() != 1.0 || model->optimizationDirection() != 1.0))
scaling = true;
const double * cost = NULL;
if (model)
cost = model->costRegion();
if (!cost) {
// not in solve
cost = objective_;
scaling = false;
}
if (!scaling) {
if (!quadraticObjective_ || !solution || !activated_) {
return objective_;
} else {
if (refresh || !gradient_) {
if (!gradient_)
gradient_ = new double[numberExtendedColumns_];
const int * columnQuadratic = quadraticObjective_->getIndices();
const CoinBigIndex * columnQuadraticStart = quadraticObjective_->getVectorStarts();
const int * columnQuadraticLength = quadraticObjective_->getVectorLengths();
const double * quadraticElement = quadraticObjective_->getElements();
offset = 0.0;
// use current linear cost region
if (includeLinear == 1)
CoinMemcpyN(cost, numberExtendedColumns_, gradient_);
else if (includeLinear == 2)
CoinMemcpyN(objective_, numberExtendedColumns_, gradient_);
else
memset(gradient_, 0, numberExtendedColumns_ * sizeof(double));
if (activated_) {
if (!fullMatrix_) {
int iColumn;
for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
double valueI = solution[iColumn];
CoinBigIndex j;
for (j = columnQuadraticStart[iColumn];
j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
int jColumn = columnQuadratic[j];
double valueJ = solution[jColumn];
double elementValue = quadraticElement[j];
if (iColumn != jColumn) {
offset += valueI * valueJ * elementValue;
//if (fabs(valueI*valueJ*elementValue)>1.0e-12)
//printf("%d %d %g %g %g -> %g\n",
// iColumn,jColumn,valueI,valueJ,elementValue,
// valueI*valueJ*elementValue);
double gradientI = valueJ * elementValue;
double gradientJ = valueI * elementValue;
gradient_[iColumn] += gradientI;
gradient_[jColumn] += gradientJ;
} else {
offset += 0.5 * valueI * valueI * elementValue;
//if (fabs(valueI*valueI*elementValue)>1.0e-12)
//printf("XX %d %g %g -> %g\n",
// iColumn,valueI,elementValue,
// 0.5*valueI*valueI*elementValue);
double gradientI = valueI * elementValue;
gradient_[iColumn] += gradientI;
}
}
}
} else {
// full matrix
int iColumn;
offset *= 2.0;
for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
CoinBigIndex j;
double value = 0.0;
double current = gradient_[iColumn];
for (j = columnQuadraticStart[iColumn];
j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
int jColumn = columnQuadratic[j];
double valueJ = solution[jColumn] * quadraticElement[j];
value += valueJ;
}
offset += value * solution[iColumn];
gradient_[iColumn] = current + value;
}
offset *= 0.5;
}
}
}
if (model)
offset *= model->optimizationDirection() * model->objectiveScale();
return gradient_;
}
} else {
// do scaling
assert (solution);
// for now only if half
assert (!fullMatrix_);
if (refresh || !gradient_) {
if (!gradient_)
gradient_ = new double[numberExtendedColumns_];
double direction = model->optimizationDirection() * model->objectiveScale();
// direction is actually scale out not scale in
//if (direction)
//direction = 1.0/direction;
const int * columnQuadratic = quadraticObjective_->getIndices();
const CoinBigIndex * columnQuadraticStart = quadraticObjective_->getVectorStarts();
const int * columnQuadraticLength = quadraticObjective_->getVectorLengths();
const double * quadraticElement = quadraticObjective_->getElements();
int iColumn;
const double * columnScale = model->columnScale();
// use current linear cost region (already scaled)
if (includeLinear == 1) {
CoinMemcpyN(model->costRegion(), numberExtendedColumns_, gradient_);
} else if (includeLinear == 2) {
memset(gradient_ + numberColumns_, 0, (numberExtendedColumns_ - numberColumns_)*sizeof(double));
if (!columnScale) {
for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
gradient_[iColumn] = objective_[iColumn] * direction;
}
} else {
for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
gradient_[iColumn] = objective_[iColumn] * direction * columnScale[iColumn];
}
}
} else {
memset(gradient_, 0, numberExtendedColumns_ * sizeof(double));
}
if (!columnScale) {
if (activated_) {
for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
double valueI = solution[iColumn];
CoinBigIndex j;
for (j = columnQuadraticStart[iColumn];
j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
int jColumn = columnQuadratic[j];
double valueJ = solution[jColumn];
double elementValue = quadraticElement[j];
elementValue *= direction;
if (iColumn != jColumn) {
offset += valueI * valueJ * elementValue;
double gradientI = valueJ * elementValue;
double gradientJ = valueI * elementValue;
gradient_[iColumn] += gradientI;
gradient_[jColumn] += gradientJ;
} else {
offset += 0.5 * valueI * valueI * elementValue;
double gradientI = valueI * elementValue;
gradient_[iColumn] += gradientI;
}
}
}
}
} else {
// scaling
if (activated_) {
for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
double valueI = solution[iColumn];
double scaleI = columnScale[iColumn] * direction;
CoinBigIndex j;
for (j = columnQuadraticStart[iColumn];
j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
int jColumn = columnQuadratic[j];
double valueJ = solution[jColumn];
double elementValue = quadraticElement[j];
double scaleJ = columnScale[jColumn];
elementValue *= scaleI * scaleJ;
if (iColumn != jColumn) {
offset += valueI * valueJ * elementValue;
double gradientI = valueJ * elementValue;
double gradientJ = valueI * elementValue;
gradient_[iColumn] += gradientI;
gradient_[jColumn] += gradientJ;
} else {
offset += 0.5 * valueI * valueI * elementValue;
double gradientI = valueI * elementValue;
gradient_[iColumn] += gradientI;
}
}
}
}
}
}
if (model)
offset *= model->optimizationDirection();
return gradient_;
}
}
//-------------------------------------------------------------------
// Clone
//-------------------------------------------------------------------
ClpObjective * ClpQuadraticObjective::clone() const
{
return new ClpQuadraticObjective(*this);
}
/* Subset clone. Duplicates are allowed
and order is as given.
*/
ClpObjective *
ClpQuadraticObjective::subsetClone (int numberColumns,
const int * whichColumns) const
{
return new ClpQuadraticObjective(*this, numberColumns, whichColumns);
}
// Resize objective
void
ClpQuadraticObjective::resize(int newNumberColumns)
{
if (numberColumns_ != newNumberColumns) {
int newExtended = newNumberColumns + (numberExtendedColumns_ - numberColumns_);
int i;
double * newArray = new double[newExtended];
if (objective_)
CoinMemcpyN(objective_, CoinMin(newExtended, numberExtendedColumns_), newArray);
delete [] objective_;
objective_ = newArray;
for (i = numberColumns_; i < newNumberColumns; i++)
objective_[i] = 0.0;
if (gradient_) {
newArray = new double[newExtended];
if (gradient_)
CoinMemcpyN(gradient_, CoinMin(newExtended, numberExtendedColumns_), newArray);
delete [] gradient_;
gradient_ = newArray;
for (i = numberColumns_; i < newNumberColumns; i++)
gradient_[i] = 0.0;
}
if (quadraticObjective_) {
if (newNumberColumns < numberColumns_) {
int * which = new int[numberColumns_-newNumberColumns];
int i;
for (i = newNumberColumns; i < numberColumns_; i++)
which[i-newNumberColumns] = i;
quadraticObjective_->deleteRows(numberColumns_ - newNumberColumns, which);
quadraticObjective_->deleteCols(numberColumns_ - newNumberColumns, which);
delete [] which;
} else {
quadraticObjective_->setDimensions(newNumberColumns, newNumberColumns);
}
}
numberColumns_ = newNumberColumns;
numberExtendedColumns_ = newExtended;
}
}
// Delete columns in objective
void
ClpQuadraticObjective::deleteSome(int numberToDelete, const int * which)
{
int newNumberColumns = numberColumns_ - numberToDelete;
int newExtended = numberExtendedColumns_ - numberToDelete;
if (objective_) {
int i ;
char * deleted = new char[numberColumns_];
int numberDeleted = 0;
memset(deleted, 0, numberColumns_ * sizeof(char));
for (i = 0; i < numberToDelete; i++) {
int j = which[i];
if (j >= 0 && j < numberColumns_ && !deleted[j]) {
numberDeleted++;
deleted[j] = 1;
}
}
newNumberColumns = numberColumns_ - numberDeleted;
newExtended = numberExtendedColumns_ - numberDeleted;
double * newArray = new double[newExtended];
int put = 0;
for (i = 0; i < numberColumns_; i++) {
if (!deleted[i]) {
newArray[put++] = objective_[i];
}
}
delete [] objective_;
objective_ = newArray;
delete [] deleted;
CoinMemcpyN(objective_ + numberColumns_, (numberExtendedColumns_ - numberColumns_),
objective_ + newNumberColumns);
}
if (gradient_) {
int i ;
char * deleted = new char[numberColumns_];
int numberDeleted = 0;
memset(deleted, 0, numberColumns_ * sizeof(char));
for (i = 0; i < numberToDelete; i++) {
int j = which[i];
if (j >= 0 && j < numberColumns_ && !deleted[j]) {
numberDeleted++;
deleted[j] = 1;
}
}
newNumberColumns = numberColumns_ - numberDeleted;
newExtended = numberExtendedColumns_ - numberDeleted;
double * newArray = new double[newExtended];
int put = 0;
for (i = 0; i < numberColumns_; i++) {
if (!deleted[i]) {
newArray[put++] = gradient_[i];
}
}
delete [] gradient_;
gradient_ = newArray;
delete [] deleted;
CoinMemcpyN(gradient_ + numberColumns_, (numberExtendedColumns_ - numberColumns_),
gradient_ + newNumberColumns);
}
numberColumns_ = newNumberColumns;
numberExtendedColumns_ = newExtended;
if (quadraticObjective_) {
quadraticObjective_->deleteCols(numberToDelete, which);
quadraticObjective_->deleteRows(numberToDelete, which);
}
}
// Load up quadratic objective
void
ClpQuadraticObjective::loadQuadraticObjective(const int numberColumns, const CoinBigIndex * start,
const int * column, const double * element, int numberExtended)
{
fullMatrix_ = false;
delete quadraticObjective_;
quadraticObjective_ = new CoinPackedMatrix(true, numberColumns, numberColumns,
start[numberColumns], element, column, start, NULL);
numberColumns_ = numberColumns;
if (numberExtended > numberExtendedColumns_) {
if (objective_) {
// make correct size
double * newArray = new double[numberExtended];
CoinMemcpyN(objective_, numberColumns_, newArray);
delete [] objective_;
objective_ = newArray;
memset(objective_ + numberColumns_, 0, (numberExtended - numberColumns_)*sizeof(double));
}
if (gradient_) {
// make correct size
double * newArray = new double[numberExtended];
CoinMemcpyN(gradient_, numberColumns_, newArray);
delete [] gradient_;
gradient_ = newArray;
memset(gradient_ + numberColumns_, 0, (numberExtended - numberColumns_)*sizeof(double));
}
numberExtendedColumns_ = numberExtended;
} else {
numberExtendedColumns_ = numberColumns_;
}
}
void
ClpQuadraticObjective::loadQuadraticObjective ( const CoinPackedMatrix& matrix)
{
delete quadraticObjective_;
quadraticObjective_ = new CoinPackedMatrix(matrix);
}
// Get rid of quadratic objective
void
ClpQuadraticObjective::deleteQuadraticObjective()
{
delete quadraticObjective_;
quadraticObjective_ = NULL;
}
/* Returns reduced gradient.Returns an offset (to be added to current one).
*/
double
ClpQuadraticObjective::reducedGradient(ClpSimplex * model, double * region,
bool useFeasibleCosts)
{
int numberRows = model->numberRows();
int numberColumns = model->numberColumns();
//work space
CoinIndexedVector * workSpace = model->rowArray(0);
CoinIndexedVector arrayVector;
arrayVector.reserve(numberRows + 1);
int iRow;
#ifdef CLP_DEBUG
workSpace->checkClear();
#endif
double * array = arrayVector.denseVector();
int * index = arrayVector.getIndices();
int number = 0;
const double * costNow = gradient(model, model->solutionRegion(), offset_,
true, useFeasibleCosts ? 2 : 1);
double * cost = model->costRegion();
const int * pivotVariable = model->pivotVariable();
for (iRow = 0; iRow < numberRows; iRow++) {
int iPivot = pivotVariable[iRow];
double value;
if (iPivot < numberColumns)
value = costNow[iPivot];
else if (!useFeasibleCosts)
value = cost[iPivot];
else
value = 0.0;
if (value) {
array[iRow] = value;
index[number++] = iRow;
}
}
arrayVector.setNumElements(number);
// Btran basic costs
model->factorization()->updateColumnTranspose(workSpace, &arrayVector);
double * work = workSpace->denseVector();
ClpFillN(work, numberRows, 0.0);
// now look at dual solution
double * rowReducedCost = region + numberColumns;
double * dual = rowReducedCost;
const double * rowCost = cost + numberColumns;
for (iRow = 0; iRow < numberRows; iRow++) {
dual[iRow] = array[iRow];
}
double * dj = region;
ClpDisjointCopyN(costNow, numberColumns, dj);
model->transposeTimes(-1.0, dual, dj);
for (iRow = 0; iRow < numberRows; iRow++) {
// slack
double value = dual[iRow];
value += rowCost[iRow];
rowReducedCost[iRow] = value;
}
return offset_;
}
/* Returns step length which gives minimum of objective for
solution + theta * change vector up to maximum theta.
arrays are numberColumns+numberRows
*/
double
ClpQuadraticObjective::stepLength(ClpSimplex * model,
const double * solution,
const double * change,
double maximumTheta,
double & currentObj,
double & predictedObj,
double & thetaObj)
{
const double * cost = model->costRegion();
bool inSolve = true;
if (!cost) {
// not in solve
cost = objective_;
inSolve = false;
}
double delta = 0.0;
double linearCost = 0.0;
int numberRows = model->numberRows();
int numberColumns = model->numberColumns();
int numberTotal = numberColumns;
if (inSolve)
numberTotal += numberRows;
currentObj = 0.0;
thetaObj = 0.0;
for (int iColumn = 0; iColumn < numberTotal; iColumn++) {
delta += cost[iColumn] * change[iColumn];
linearCost += cost[iColumn] * solution[iColumn];
}
if (!activated_ || !quadraticObjective_) {
currentObj = linearCost;
thetaObj = currentObj + delta * maximumTheta;
if (delta < 0.0) {
return maximumTheta;
} else {
COIN_DETAIL_PRINT(printf("odd linear direction %g\n", delta));
return 0.0;
}
}
assert (model);
bool scaling = false;
if ((model->rowScale() ||
model->objectiveScale() != 1.0 || model->optimizationDirection() != 1.0) && inSolve)
scaling = true;
const int * columnQuadratic = quadraticObjective_->getIndices();
const CoinBigIndex * columnQuadraticStart = quadraticObjective_->getVectorStarts();
const int * columnQuadraticLength = quadraticObjective_->getVectorLengths();
const double * quadraticElement = quadraticObjective_->getElements();
double a = 0.0;
double b = delta;
double c = 0.0;
if (!scaling) {
if (!fullMatrix_) {
int iColumn;
for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
double valueI = solution[iColumn];
double changeI = change[iColumn];
CoinBigIndex j;
for (j = columnQuadraticStart[iColumn];
j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
int jColumn = columnQuadratic[j];
double valueJ = solution[jColumn];
double changeJ = change[jColumn];
double elementValue = quadraticElement[j];
if (iColumn != jColumn) {
a += changeI * changeJ * elementValue;
b += (changeI * valueJ + changeJ * valueI) * elementValue;
c += valueI * valueJ * elementValue;
} else {
a += 0.5 * changeI * changeI * elementValue;
b += changeI * valueI * elementValue;
c += 0.5 * valueI * valueI * elementValue;
}
}
}
} else {
// full matrix stored
int iColumn;
for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
double valueI = solution[iColumn];
double changeI = change[iColumn];
CoinBigIndex j;
for (j = columnQuadraticStart[iColumn];
j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
int jColumn = columnQuadratic[j];
double valueJ = solution[jColumn];
double changeJ = change[jColumn];
double elementValue = quadraticElement[j];
valueJ *= elementValue;
a += changeI * changeJ * elementValue;
b += changeI * valueJ;
c += valueI * valueJ;
}
}
a *= 0.5;
c *= 0.5;
}
} else {
// scaling
// for now only if half
assert (!fullMatrix_);
const double * columnScale = model->columnScale();
double direction = model->optimizationDirection() * model->objectiveScale();
// direction is actually scale out not scale in
if (direction)
direction = 1.0 / direction;
if (!columnScale) {
for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
double valueI = solution[iColumn];
double changeI = change[iColumn];
CoinBigIndex j;
for (j = columnQuadraticStart[iColumn];
j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
int jColumn = columnQuadratic[j];
double valueJ = solution[jColumn];
double changeJ = change[jColumn];
double elementValue = quadraticElement[j];
elementValue *= direction;
if (iColumn != jColumn) {
a += changeI * changeJ * elementValue;
b += (changeI * valueJ + changeJ * valueI) * elementValue;
c += valueI * valueJ * elementValue;
} else {
a += 0.5 * changeI * changeI * elementValue;
b += changeI * valueI * elementValue;
c += 0.5 * valueI * valueI * elementValue;
}
}
}
} else {
// scaling
for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
double valueI = solution[iColumn];
double changeI = change[iColumn];
double scaleI = columnScale[iColumn] * direction;
CoinBigIndex j;
for (j = columnQuadraticStart[iColumn];
j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
int jColumn = columnQuadratic[j];
double valueJ = solution[jColumn];
double changeJ = change[jColumn];
double elementValue = quadraticElement[j];
elementValue *= scaleI * columnScale[jColumn];
if (iColumn != jColumn) {
a += changeI * changeJ * elementValue;
b += (changeI * valueJ + changeJ * valueI) * elementValue;
c += valueI * valueJ * elementValue;
} else {
a += 0.5 * changeI * changeI * elementValue;
b += changeI * valueI * elementValue;
c += 0.5 * valueI * valueI * elementValue;
}
}
}
}
}
double theta;
//printf("Current cost %g\n",c+linearCost);
currentObj = c + linearCost;
thetaObj = currentObj + a * maximumTheta * maximumTheta + b * maximumTheta;
// minimize a*x*x + b*x + c
if (a <= 0.0) {
theta = maximumTheta;
} else {
theta = -0.5 * b / a;
}
predictedObj = currentObj + a * theta * theta + b * theta;
if (b > 0.0) {
if (model->messageHandler()->logLevel() & 32)
printf("a %g b %g c %g => %g\n", a, b, c, theta);
b = 0.0;
}
return CoinMin(theta, maximumTheta);
}
// Return objective value (without any ClpModel offset) (model may be NULL)
double
ClpQuadraticObjective::objectiveValue(const ClpSimplex * model, const double * solution) const
{
bool scaling = false;
if (model && (model->rowScale() ||
model->objectiveScale() != 1.0))
scaling = true;
const double * cost = NULL;
if (model)
cost = model->costRegion();
if (!cost) {
// not in solve
cost = objective_;
scaling = false;
}
double linearCost = 0.0;
int numberColumns = model->numberColumns();
int numberTotal = numberColumns;
double currentObj = 0.0;
for (int iColumn = 0; iColumn < numberTotal; iColumn++) {
linearCost += cost[iColumn] * solution[iColumn];
}
if (!activated_ || !quadraticObjective_) {
currentObj = linearCost;
return currentObj;
}
assert (model);
const int * columnQuadratic = quadraticObjective_->getIndices();
const CoinBigIndex * columnQuadraticStart = quadraticObjective_->getVectorStarts();
const int * columnQuadraticLength = quadraticObjective_->getVectorLengths();
const double * quadraticElement = quadraticObjective_->getElements();
double c = 0.0;
if (!scaling) {
if (!fullMatrix_) {
int iColumn;
for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
double valueI = solution[iColumn];
CoinBigIndex j;
for (j = columnQuadraticStart[iColumn];
j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
int jColumn = columnQuadratic[j];
double valueJ = solution[jColumn];
double elementValue = quadraticElement[j];
if (iColumn != jColumn) {
c += valueI * valueJ * elementValue;
} else {
c += 0.5 * valueI * valueI * elementValue;
}
}
}
} else {
// full matrix stored
int iColumn;
for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
double valueI = solution[iColumn];
CoinBigIndex j;
for (j = columnQuadraticStart[iColumn];
j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
int jColumn = columnQuadratic[j];
double valueJ = solution[jColumn];
double elementValue = quadraticElement[j];
valueJ *= elementValue;
c += valueI * valueJ;
}
}
c *= 0.5;
}
} else {
// scaling
// for now only if half
assert (!fullMatrix_);
const double * columnScale = model->columnScale();
double direction = model->objectiveScale();
// direction is actually scale out not scale in
if (direction)
direction = 1.0 / direction;
if (!columnScale) {
for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
double valueI = solution[iColumn];
CoinBigIndex j;
for (j = columnQuadraticStart[iColumn];
j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
int jColumn = columnQuadratic[j];
double valueJ = solution[jColumn];
double elementValue = quadraticElement[j];
elementValue *= direction;
if (iColumn != jColumn) {
c += valueI * valueJ * elementValue;
} else {
c += 0.5 * valueI * valueI * elementValue;
}
}
}
} else {
// scaling
for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
double valueI = solution[iColumn];
double scaleI = columnScale[iColumn] * direction;
CoinBigIndex j;
for (j = columnQuadraticStart[iColumn];
j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
int jColumn = columnQuadratic[j];
double valueJ = solution[jColumn];
double elementValue = quadraticElement[j];
elementValue *= scaleI * columnScale[jColumn];
if (iColumn != jColumn) {
c += valueI * valueJ * elementValue;
} else {
c += 0.5 * valueI * valueI * elementValue;
}
}
}
}
}
currentObj = c + linearCost;
return currentObj;
}
// Scale objective
void
ClpQuadraticObjective::reallyScale(const double * columnScale)
{
const int * columnQuadratic = quadraticObjective_->getIndices();
const CoinBigIndex * columnQuadraticStart = quadraticObjective_->getVectorStarts();
const int * columnQuadraticLength = quadraticObjective_->getVectorLengths();
double * quadraticElement = quadraticObjective_->getMutableElements();
for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
double scaleI = columnScale[iColumn];
objective_[iColumn] *= scaleI;
CoinBigIndex j;
for (j = columnQuadraticStart[iColumn];
j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
int jColumn = columnQuadratic[j];
quadraticElement[j] *= scaleI * columnScale[jColumn];
}
}
}
/* Given a zeroed array sets nonlinear columns to 1.
Returns number of nonlinear columns
*/
int
ClpQuadraticObjective::markNonlinear(char * which)
{
int iColumn;
const int * columnQuadratic = quadraticObjective_->getIndices();
const CoinBigIndex * columnQuadraticStart = quadraticObjective_->getVectorStarts();
const int * columnQuadraticLength = quadraticObjective_->getVectorLengths();
for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
CoinBigIndex j;
for (j = columnQuadraticStart[iColumn];
j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
int jColumn = columnQuadratic[j];
which[jColumn] = 1;
which[iColumn] = 1;
}
}
int numberNonLinearColumns = 0;
for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
if(which[iColumn])
numberNonLinearColumns++;
}
return numberNonLinearColumns;
}