limp-cbc-0.3.2.0: cbits/coin/ClpPresolve.cpp
/* $Id: ClpPresolve.cpp 1931 2013-04-06 20:44:29Z stefan $ */
// Copyright (C) 2002, International Business Machines
// Corporation and others. All Rights Reserved.
// This code is licensed under the terms of the Eclipse Public License (EPL).
//#define PRESOLVE_CONSISTENCY 1
//#define PRESOLVE_DEBUG 1
#include <stdio.h>
#include <cassert>
#include <iostream>
#include "CoinHelperFunctions.hpp"
#include "ClpConfig.h"
#ifdef CLP_HAS_ABC
#include "CoinAbcCommon.hpp"
#endif
#include "CoinPackedMatrix.hpp"
#include "ClpPackedMatrix.hpp"
#include "ClpSimplex.hpp"
#include "ClpSimplexOther.hpp"
#ifndef SLIM_CLP
#include "ClpQuadraticObjective.hpp"
#endif
#include "ClpPresolve.hpp"
#include "CoinPresolveMatrix.hpp"
#include "CoinPresolveEmpty.hpp"
#include "CoinPresolveFixed.hpp"
#include "CoinPresolvePsdebug.hpp"
#include "CoinPresolveSingleton.hpp"
#include "CoinPresolveDoubleton.hpp"
#include "CoinPresolveTripleton.hpp"
#include "CoinPresolveZeros.hpp"
#include "CoinPresolveSubst.hpp"
#include "CoinPresolveForcing.hpp"
#include "CoinPresolveDual.hpp"
#include "CoinPresolveTighten.hpp"
#include "CoinPresolveUseless.hpp"
#include "CoinPresolveDupcol.hpp"
#include "CoinPresolveImpliedFree.hpp"
#include "CoinPresolveIsolated.hpp"
#include "CoinMessage.hpp"
ClpPresolve::ClpPresolve() :
originalModel_(NULL),
presolvedModel_(NULL),
nonLinearValue_(0.0),
originalColumn_(NULL),
originalRow_(NULL),
rowObjective_(NULL),
paction_(0),
ncols_(0),
nrows_(0),
nelems_(0),
#ifdef ABC_INHERIT
numberPasses_(20),
#else
numberPasses_(5),
#endif
substitution_(3),
#ifndef CLP_NO_STD
saveFile_(""),
#endif
presolveActions_(0)
{
}
ClpPresolve::~ClpPresolve()
{
destroyPresolve();
}
// Gets rid of presolve actions (e.g.when infeasible)
void
ClpPresolve::destroyPresolve()
{
const CoinPresolveAction *paction = paction_;
while (paction) {
const CoinPresolveAction *next = paction->next;
delete paction;
paction = next;
}
delete [] originalColumn_;
delete [] originalRow_;
paction_ = NULL;
originalColumn_ = NULL;
originalRow_ = NULL;
delete [] rowObjective_;
rowObjective_ = NULL;
}
/* This version of presolve returns a pointer to a new presolved
model. NULL if infeasible
*/
ClpSimplex *
ClpPresolve::presolvedModel(ClpSimplex & si,
double feasibilityTolerance,
bool keepIntegers,
int numberPasses,
bool dropNames,
bool doRowObjective,
const char * prohibitedRows,
const char * prohibitedColumns)
{
// Check matrix
int checkType = ((si.specialOptions() & 128) != 0) ? 14 : 15;
if (!si.clpMatrix()->allElementsInRange(&si, si.getSmallElementValue(),
1.0e20,checkType))
return NULL;
else
return gutsOfPresolvedModel(&si, feasibilityTolerance, keepIntegers, numberPasses, dropNames,
doRowObjective,
prohibitedRows,
prohibitedColumns);
}
#ifndef CLP_NO_STD
/* This version of presolve updates
model and saves original data to file. Returns non-zero if infeasible
*/
int
ClpPresolve::presolvedModelToFile(ClpSimplex &si, std::string fileName,
double feasibilityTolerance,
bool keepIntegers,
int numberPasses,
bool dropNames,
bool doRowObjective)
{
// Check matrix
if (!si.clpMatrix()->allElementsInRange(&si, si.getSmallElementValue(),
1.0e20))
return 2;
saveFile_ = fileName;
si.saveModel(saveFile_.c_str());
ClpSimplex * model = gutsOfPresolvedModel(&si, feasibilityTolerance, keepIntegers, numberPasses, dropNames,
doRowObjective);
if (model == &si) {
return 0;
} else {
si.restoreModel(saveFile_.c_str());
remove(saveFile_.c_str());
return 1;
}
}
#endif
// Return pointer to presolved model
ClpSimplex *
ClpPresolve::model() const
{
return presolvedModel_;
}
// Return pointer to original model
ClpSimplex *
ClpPresolve::originalModel() const
{
return originalModel_;
}
// Return presolve status (0,1,2)
int
ClpPresolve::presolveStatus() const
{
if (nelems_>=0) {
// feasible (or not done yet)
return 0;
} else {
int presolveStatus = - nelems_;
// If both infeasible and unbounded - say infeasible
if (presolveStatus>2)
presolveStatus = 1;
return presolveStatus;
}
}
void
ClpPresolve::postsolve(bool updateStatus)
{
// Return at once if no presolved model
if (!presolvedModel_)
return;
// Messages
CoinMessages messages = originalModel_->coinMessages();
if (!presolvedModel_->isProvenOptimal()) {
presolvedModel_->messageHandler()->message(COIN_PRESOLVE_NONOPTIMAL,
messages)
<< CoinMessageEol;
}
// this is the size of the original problem
const int ncols0 = ncols_;
const int nrows0 = nrows_;
const CoinBigIndex nelems0 = nelems_;
// this is the reduced problem
int ncols = presolvedModel_->getNumCols();
int nrows = presolvedModel_->getNumRows();
double * acts = NULL;
double * sol = NULL;
unsigned char * rowstat = NULL;
unsigned char * colstat = NULL;
#ifndef CLP_NO_STD
if (saveFile_ == "") {
#endif
// reality check
assert(ncols0 == originalModel_->getNumCols());
assert(nrows0 == originalModel_->getNumRows());
acts = originalModel_->primalRowSolution();
sol = originalModel_->primalColumnSolution();
if (updateStatus) {
// postsolve does not know about fixed
int i;
for (i = 0; i < nrows + ncols; i++) {
if (presolvedModel_->getColumnStatus(i) == ClpSimplex::isFixed)
presolvedModel_->setColumnStatus(i, ClpSimplex::atLowerBound);
}
unsigned char *status = originalModel_->statusArray();
if (!status) {
originalModel_->createStatus();
status = originalModel_->statusArray();
}
rowstat = status + ncols0;
colstat = status;
CoinMemcpyN( presolvedModel_->statusArray(), ncols, colstat);
CoinMemcpyN( presolvedModel_->statusArray() + ncols, nrows, rowstat);
}
#ifndef CLP_NO_STD
} else {
// from file
acts = new double[nrows0];
sol = new double[ncols0];
CoinZeroN(acts, nrows0);
CoinZeroN(sol, ncols0);
if (updateStatus) {
unsigned char *status = new unsigned char [nrows0+ncols0];
rowstat = status + ncols0;
colstat = status;
CoinMemcpyN( presolvedModel_->statusArray(), ncols, colstat);
CoinMemcpyN( presolvedModel_->statusArray() + ncols, nrows, rowstat);
}
}
#endif
// CoinPostsolveMatrix object assumes ownership of sol, acts, colstat;
// will be deleted by ~CoinPostsolveMatrix. delete[] operations below
// cause duplicate free. In case where saveFile == "", as best I can see
// arrays are owned by originalModel_. fix is to
// clear fields in prob to avoid delete[] in ~CoinPostsolveMatrix.
CoinPostsolveMatrix prob(presolvedModel_,
ncols0,
nrows0,
nelems0,
presolvedModel_->getObjSense(),
// end prepost
sol, acts,
colstat, rowstat);
postsolve(prob);
#ifndef CLP_NO_STD
if (saveFile_ != "") {
// From file
assert (originalModel_ == presolvedModel_);
originalModel_->restoreModel(saveFile_.c_str());
remove(saveFile_.c_str());
CoinMemcpyN(acts, nrows0, originalModel_->primalRowSolution());
// delete [] acts;
CoinMemcpyN(sol, ncols0, originalModel_->primalColumnSolution());
// delete [] sol;
if (updateStatus) {
CoinMemcpyN(colstat, nrows0 + ncols0, originalModel_->statusArray());
// delete [] colstat;
}
} else {
#endif
prob.sol_ = 0 ;
prob.acts_ = 0 ;
prob.colstat_ = 0 ;
#ifndef CLP_NO_STD
}
#endif
// put back duals
CoinMemcpyN(prob.rowduals_, nrows_, originalModel_->dualRowSolution());
double maxmin = originalModel_->getObjSense();
if (maxmin < 0.0) {
// swap signs
int i;
double * pi = originalModel_->dualRowSolution();
for (i = 0; i < nrows_; i++)
pi[i] = -pi[i];
}
// Now check solution
double offset;
CoinMemcpyN(originalModel_->objectiveAsObject()->gradient(originalModel_,
originalModel_->primalColumnSolution(), offset, true),
ncols_, originalModel_->dualColumnSolution());
originalModel_->clpMatrix()->transposeTimes(-1.0,
originalModel_->dualRowSolution(),
originalModel_->dualColumnSolution());
memset(originalModel_->primalRowSolution(), 0, nrows_ * sizeof(double));
originalModel_->clpMatrix()->times(1.0,
originalModel_->primalColumnSolution(),
originalModel_->primalRowSolution());
originalModel_->checkSolutionInternal();
if (originalModel_->sumDualInfeasibilities() > 1.0e-1) {
// See if we can fix easily
static_cast<ClpSimplexOther *> (originalModel_)->cleanupAfterPostsolve();
}
// Messages
presolvedModel_->messageHandler()->message(COIN_PRESOLVE_POSTSOLVE,
messages)
<< originalModel_->objectiveValue()
<< originalModel_->sumDualInfeasibilities()
<< originalModel_->numberDualInfeasibilities()
<< originalModel_->sumPrimalInfeasibilities()
<< originalModel_->numberPrimalInfeasibilities()
<< CoinMessageEol;
//originalModel_->objectiveValue_=objectiveValue_;
originalModel_->setNumberIterations(presolvedModel_->numberIterations());
if (!presolvedModel_->status()) {
if (!originalModel_->numberDualInfeasibilities() &&
!originalModel_->numberPrimalInfeasibilities()) {
originalModel_->setProblemStatus( 0);
} else {
originalModel_->setProblemStatus( -1);
// Say not optimal after presolve
originalModel_->setSecondaryStatus(7);
presolvedModel_->messageHandler()->message(COIN_PRESOLVE_NEEDS_CLEANING,
messages)
<< CoinMessageEol;
}
} else {
originalModel_->setProblemStatus( presolvedModel_->status());
// but not if close to feasible
if( originalModel_->sumPrimalInfeasibilities()<1.0e-1) {
originalModel_->setProblemStatus( -1);
// Say not optimal after presolve
originalModel_->setSecondaryStatus(7);
}
}
#ifndef CLP_NO_STD
if (saveFile_ != "")
presolvedModel_ = NULL;
#endif
}
// return pointer to original columns
const int *
ClpPresolve::originalColumns() const
{
return originalColumn_;
}
// return pointer to original rows
const int *
ClpPresolve::originalRows() const
{
return originalRow_;
}
// Set pointer to original model
void
ClpPresolve::setOriginalModel(ClpSimplex * model)
{
originalModel_ = model;
}
#if 0
// A lazy way to restrict which transformations are applied
// during debugging.
static int ATOI(const char *name)
{
return true;
#if PRESOLVE_DEBUG || PRESOLVE_SUMMARY
if (getenv(name)) {
int val = atoi(getenv(name));
printf("%s = %d\n", name, val);
return (val);
} else {
if (strcmp(name, "off"))
return (true);
else
return (false);
}
#else
return (true);
#endif
}
#endif
//#define PRESOLVE_DEBUG 1
#if PRESOLVE_DEBUG
void check_sol(CoinPresolveMatrix *prob, double tol)
{
double *colels = prob->colels_;
int *hrow = prob->hrow_;
int *mcstrt = prob->mcstrt_;
int *hincol = prob->hincol_;
int *hinrow = prob->hinrow_;
int ncols = prob->ncols_;
double * csol = prob->sol_;
double * acts = prob->acts_;
double * clo = prob->clo_;
double * cup = prob->cup_;
int nrows = prob->nrows_;
double * rlo = prob->rlo_;
double * rup = prob->rup_;
int colx;
double * rsol = new double[nrows];
memset(rsol, 0, nrows * sizeof(double));
for (colx = 0; colx < ncols; ++colx) {
if (1) {
CoinBigIndex k = mcstrt[colx];
int nx = hincol[colx];
double solutionValue = csol[colx];
for (int i = 0; i < nx; ++i) {
int row = hrow[k];
double coeff = colels[k];
k++;
rsol[row] += solutionValue * coeff;
}
if (csol[colx] < clo[colx] - tol) {
printf("low CSOL: %d - %g %g %g\n",
colx, clo[colx], csol[colx], cup[colx]);
} else if (csol[colx] > cup[colx] + tol) {
printf("high CSOL: %d - %g %g %g\n",
colx, clo[colx], csol[colx], cup[colx]);
}
}
}
int rowx;
for (rowx = 0; rowx < nrows; ++rowx) {
if (hinrow[rowx]) {
if (fabs(rsol[rowx] - acts[rowx]) > tol)
printf("inacc RSOL: %d - %g %g (acts_ %g) %g\n",
rowx, rlo[rowx], rsol[rowx], acts[rowx], rup[rowx]);
if (rsol[rowx] < rlo[rowx] - tol) {
printf("low RSOL: %d - %g %g %g\n",
rowx, rlo[rowx], rsol[rowx], rup[rowx]);
} else if (rsol[rowx] > rup[rowx] + tol ) {
printf("high RSOL: %d - %g %g %g\n",
rowx, rlo[rowx], rsol[rowx], rup[rowx]);
}
}
}
delete [] rsol;
}
#endif
static int tightenDoubletons2(CoinPresolveMatrix * prob)
{
// column-major representation
const int ncols = prob->ncols_ ;
const CoinBigIndex *const mcstrt = prob->mcstrt_ ;
const int *const hincol = prob->hincol_ ;
const int *const hrow = prob->hrow_ ;
double * colels = prob->colels_ ;
double * cost = prob->cost_ ;
// column type, bounds, solution, and status
const unsigned char *const integerType = prob->integerType_ ;
double * clo = prob->clo_ ;
double * cup = prob->cup_ ;
// row-major representation
//const int nrows = prob->nrows_ ;
const CoinBigIndex *const mrstrt = prob->mrstrt_ ;
const int *const hinrow = prob->hinrow_ ;
const int *const hcol = prob->hcol_ ;
double * rowels = prob->rowels_ ;
// row bounds
double *const rlo = prob->rlo_ ;
double *const rup = prob->rup_ ;
// tolerances
//const double ekkinf2 = PRESOLVE_SMALL_INF ;
//const double ekkinf = ekkinf2*1.0e8 ;
//const double ztolcbarj = prob->ztoldj_ ;
//const CoinRelFltEq relEq(prob->ztolzb_) ;
int numberChanged=0;
double bound[2];
double alpha[2]={0.0,0.0};
double offset=0.0;
for (int icol=0;icol<ncols;icol++) {
if (hincol[icol]==2) {
CoinBigIndex start=mcstrt[icol];
int row0 = hrow[start];
if (hinrow[row0]!=2)
continue;
int row1 = hrow[start+1];
if (hinrow[row1]!=2)
continue;
double element0 = colels[start];
double rowUpper0=rup[row0];
bool swapSigns0=false;
if (rlo[row0]>-1.0e30) {
if (rup[row0]>1.0e30) {
swapSigns0=true;
rowUpper0=-rlo[row0];
element0=-element0;
} else {
// range or equality
continue;
}
} else if (rup[row0]>1.0e30) {
// free
continue;
}
#if 0
// skip here for speed
// skip if no cost (should be able to get rid of)
if (!cost[icol]) {
printf("should be able to get rid of %d with no cost\n",icol);
continue;
}
// skip if negative cost for now
if (cost[icol]<0.0) {
printf("code for negative cost\n");
continue;
}
#endif
double element1 = colels[start+1];
double rowUpper1=rup[row1];
bool swapSigns1=false;
if (rlo[row1]>-1.0e30) {
if (rup[row1]>1.0e30) {
swapSigns1=true;
rowUpper1=-rlo[row1];
element1=-element1;
} else {
// range or equality
continue;
}
} else if (rup[row1]>1.0e30) {
// free
continue;
}
double lowerX=clo[icol];
double upperX=cup[icol];
int otherCol=-1;
CoinBigIndex startRow=mrstrt[row0];
for (CoinBigIndex j=startRow;j<startRow+2;j++) {
int jcol=hcol[j];
if (jcol!=icol) {
alpha[0]=swapSigns0 ? -rowels[j] :rowels[j];
otherCol=jcol;
}
}
startRow=mrstrt[row1];
bool possible=true;
for (CoinBigIndex j=startRow;j<startRow+2;j++) {
int jcol=hcol[j];
if (jcol!=icol) {
if (jcol==otherCol) {
alpha[1]=swapSigns1 ? -rowels[j] :rowels[j];
} else {
possible=false;
}
}
}
if (possible) {
// skip if no cost (should be able to get rid of)
if (!cost[icol]) {
PRESOLVE_DETAIL_PRINT(printf("should be able to get rid of %d with no cost\n",icol));
continue;
}
// skip if negative cost for now
if (cost[icol]<0.0) {
PRESOLVE_DETAIL_PRINT(printf("code for negative cost\n"));
continue;
}
bound[0]=clo[otherCol];
bound[1]=cup[otherCol];
double lowestLowest=COIN_DBL_MAX;
double highestLowest=-COIN_DBL_MAX;
double lowestHighest=COIN_DBL_MAX;
double highestHighest=-COIN_DBL_MAX;
int binding0=0;
int binding1=0;
for (int k=0;k<2;k++) {
bool infLow0=false;
bool infLow1=false;
double sum0=0.0;
double sum1=0.0;
double value=bound[k];
if (fabs(value)<1.0e30) {
sum0+=alpha[0]*value;
sum1+=alpha[1]*value;
} else {
if (alpha[0]>0.0) {
if (value<0.0)
infLow0 =true;
} else if (alpha[0]<0.0) {
if (value>0.0)
infLow0 =true;
}
if (alpha[1]>0.0) {
if (value<0.0)
infLow1 =true;
} else if (alpha[1]<0.0) {
if (value>0.0)
infLow1 =true;
}
}
/* Got sums
*/
double thisLowest0=-COIN_DBL_MAX;
double thisHighest0=COIN_DBL_MAX;
if (element0>0.0) {
// upper bound unless inf&2 !=0
if (!infLow0)
thisHighest0 = (rowUpper0-sum0)/element0;
} else {
// lower bound unless inf&2 !=0
if (!infLow0)
thisLowest0 = (rowUpper0-sum0)/element0;
}
double thisLowest1=-COIN_DBL_MAX;
double thisHighest1=COIN_DBL_MAX;
if (element1>0.0) {
// upper bound unless inf&2 !=0
if (!infLow1)
thisHighest1 = (rowUpper1-sum1)/element1;
} else {
// lower bound unless inf&2 !=0
if (!infLow1)
thisLowest1 = (rowUpper1-sum1)/element1;
}
if (thisLowest0>thisLowest1+1.0e-12) {
if (thisLowest0>lowerX+1.0e-12)
binding0|= 1<<k;
} else if (thisLowest1>thisLowest0+1.0e-12) {
if (thisLowest1>lowerX+1.0e-12)
binding1|= 1<<k;
thisLowest0=thisLowest1;
}
if (thisHighest0<thisHighest1-1.0e-12) {
if (thisHighest0<upperX-1.0e-12)
binding0|= 1<<k;
} else if (thisHighest1<thisHighest0-1.0e-12) {
if (thisHighest1<upperX-1.0e-12)
binding1|= 1<<k;
thisHighest0=thisHighest1;
}
lowestLowest=CoinMin(lowestLowest,thisLowest0);
highestHighest=CoinMax(highestHighest,thisHighest0);
lowestHighest=CoinMin(lowestHighest,thisHighest0);
highestLowest=CoinMax(highestLowest,thisLowest0);
}
// see if any good
//#define PRINT_VALUES
if (!binding0||!binding1) {
PRESOLVE_DETAIL_PRINT(printf("Row redundant for column %d\n",icol));
} else {
#ifdef PRINT_VALUES
printf("Column %d bounds %g,%g lowest %g,%g highest %g,%g\n",
icol,lowerX,upperX,lowestLowest,lowestHighest,
highestLowest,highestHighest);
#endif
// if integer adjust
if (integerType[icol]) {
lowestLowest=ceil(lowestLowest-1.0e-5);
highestLowest=ceil(highestLowest-1.0e-5);
lowestHighest=floor(lowestHighest+1.0e-5);
highestHighest=floor(highestHighest+1.0e-5);
}
// if costed may be able to adjust
if (cost[icol]>=0.0) {
if (highestLowest<upperX&&highestLowest>=lowerX&&highestHighest<1.0e30) {
highestHighest=CoinMin(highestHighest,highestLowest);
}
}
if (cost[icol]<=0.0) {
if (lowestHighest>lowerX&&lowestHighest<=upperX&&lowestHighest>-1.0e30) {
lowestLowest=CoinMax(lowestLowest,lowestHighest);
}
}
#if 1
if (lowestLowest>lowerX+1.0e-8) {
#ifdef PRINT_VALUES
printf("Can increase lower bound on %d from %g to %g\n",
icol,lowerX,lowestLowest);
#endif
lowerX=lowestLowest;
}
if (highestHighest<upperX-1.0e-8) {
#ifdef PRINT_VALUES
printf("Can decrease upper bound on %d from %g to %g\n",
icol,upperX,highestHighest);
#endif
upperX=highestHighest;
}
#endif
// see if we can move costs
double xValue;
double yValue0;
double yValue1;
double newLower=COIN_DBL_MAX;
double newUpper=-COIN_DBL_MAX;
#ifdef PRINT_VALUES
double ranges0[2];
double ranges1[2];
#endif
double costEqual;
double slope[2];
assert (binding0+binding1==3);
// get where equal
xValue=(rowUpper0*element1-rowUpper1*element0)/(alpha[0]*element1-alpha[1]*element0);
yValue0=(rowUpper0-xValue*alpha[0])/element0;
yValue1=(rowUpper1-xValue*alpha[1])/element1;
newLower=CoinMin(newLower,CoinMax(yValue0,yValue1));
newUpper=CoinMax(newUpper,CoinMax(yValue0,yValue1));
double xValueEqual=xValue;
double yValueEqual=yValue0;
costEqual = xValue*cost[otherCol]+yValueEqual*cost[icol];
if (binding0==1) {
#ifdef PRINT_VALUES
ranges0[0]=bound[0];
ranges0[1]=yValue0;
ranges1[0]=yValue0;
ranges1[1]=bound[1];
#endif
// take x 1.0 down
double x=xValue-1.0;
double y=(rowUpper0-x*alpha[0])/element0;
double costTotal = x*cost[otherCol]+y*cost[icol];
slope[0] = costEqual-costTotal;
// take x 1.0 up
x=xValue+1.0;
y=(rowUpper1-x*alpha[1])/element0;
costTotal = x*cost[otherCol]+y*cost[icol];
slope[1] = costTotal-costEqual;
} else {
#ifdef PRINT_VALUES
ranges1[0]=bound[0];
ranges1[1]=yValue0;
ranges0[0]=yValue0;
ranges0[1]=bound[1];
#endif
// take x 1.0 down
double x=xValue-1.0;
double y=(rowUpper1-x*alpha[1])/element0;
double costTotal = x*cost[otherCol]+y*cost[icol];
slope[1] = costEqual-costTotal;
// take x 1.0 up
x=xValue+1.0;
y=(rowUpper0-x*alpha[0])/element0;
costTotal = x*cost[otherCol]+y*cost[icol];
slope[0] = costTotal-costEqual;
}
#ifdef PRINT_VALUES
printf("equal value of %d is %g, value of %d is max(%g,%g) - %g\n",
otherCol,xValue,icol,yValue0,yValue1,CoinMax(yValue0,yValue1));
printf("Cost at equality %g for constraint 0 ranges %g -> %g slope %g for constraint 1 ranges %g -> %g slope %g\n",
costEqual,ranges0[0],ranges0[1],slope[0],ranges1[0],ranges1[1],slope[1]);
#endif
xValue=bound[0];
yValue0=(rowUpper0-xValue*alpha[0])/element0;
yValue1=(rowUpper1-xValue*alpha[1])/element1;
#ifdef PRINT_VALUES
printf("value of %d is %g, value of %d is max(%g,%g) - %g\n",
otherCol,xValue,icol,yValue0,yValue1,CoinMax(yValue0,yValue1));
#endif
newLower=CoinMin(newLower,CoinMax(yValue0,yValue1));
// cost>0 so will be at lower
//double yValueAtBound0=newLower;
newUpper=CoinMax(newUpper,CoinMax(yValue0,yValue1));
xValue=bound[1];
yValue0=(rowUpper0-xValue*alpha[0])/element0;
yValue1=(rowUpper1-xValue*alpha[1])/element1;
#ifdef PRINT_VALUES
printf("value of %d is %g, value of %d is max(%g,%g) - %g\n",
otherCol,xValue,icol,yValue0,yValue1,CoinMax(yValue0,yValue1));
#endif
newLower=CoinMin(newLower,CoinMax(yValue0,yValue1));
// cost>0 so will be at lower
//double yValueAtBound1=newLower;
newUpper=CoinMax(newUpper,CoinMax(yValue0,yValue1));
lowerX=CoinMax(lowerX,newLower-1.0e-12*fabs(newLower));
upperX=CoinMin(upperX,newUpper+1.0e-12*fabs(newUpper));
// Now make duplicate row
// keep row 0 so need to adjust costs so same
#ifdef PRINT_VALUES
printf("Costs for x %g,%g,%g are %g,%g,%g\n",
xValueEqual-1.0,xValueEqual,xValueEqual+1.0,
costEqual-slope[0],costEqual,costEqual+slope[1]);
#endif
double costOther=cost[otherCol]+slope[1];
double costThis=cost[icol]+slope[1]*(element0/alpha[0]);
xValue=xValueEqual;
yValue0=CoinMax((rowUpper0-xValue*alpha[0])/element0,lowerX);
double thisOffset=costEqual-(costOther*xValue+costThis*yValue0);
offset += thisOffset;
#ifdef PRINT_VALUES
printf("new cost at equal %g\n",costOther*xValue+costThis*yValue0+thisOffset);
#endif
xValue=xValueEqual-1.0;
yValue0=CoinMax((rowUpper0-xValue*alpha[0])/element0,lowerX);
#ifdef PRINT_VALUES
printf("new cost at -1 %g\n",costOther*xValue+costThis*yValue0+thisOffset);
#endif
assert(fabs((costOther*xValue+costThis*yValue0+thisOffset)-(costEqual-slope[0]))<1.0e-5);
xValue=xValueEqual+1.0;
yValue0=CoinMax((rowUpper0-xValue*alpha[0])/element0,lowerX);
#ifdef PRINT_VALUES
printf("new cost at +1 %g\n",costOther*xValue+costThis*yValue0+thisOffset);
#endif
assert(fabs((costOther*xValue+costThis*yValue0+thisOffset)-(costEqual+slope[1]))<1.0e-5);
numberChanged++;
// continue;
cost[otherCol] = costOther;
cost[icol] = costThis;
clo[icol]=lowerX;
cup[icol]=upperX;
int startCol[2];
int endCol[2];
startCol[0]=mcstrt[icol];
endCol[0]=startCol[0]+2;
startCol[1]=mcstrt[otherCol];
endCol[1]=startCol[1]+hincol[otherCol];
double values[2]={0.0,0.0};
for (int k=0;k<2;k++) {
for (CoinBigIndex i=startCol[k];i<endCol[k];i++) {
if (hrow[i]==row0)
values[k]=colels[i];
}
for (CoinBigIndex i=startCol[k];i<endCol[k];i++) {
if (hrow[i]==row1)
colels[i]=values[k];
}
}
for (CoinBigIndex i=mrstrt[row1];i<mrstrt[row1]+2;i++) {
if (hcol[i]==icol)
rowels[i]=values[0];
else
rowels[i]=values[1];
}
}
}
}
}
#if ABC_NORMAL_DEBUG>0
if (offset)
printf("Cost offset %g\n",offset);
#endif
return numberChanged;
}
//#define COIN_PRESOLVE_BUG
#ifdef COIN_PRESOLVE_BUG
static int counter=1000000;
static int startEmptyRows=0;
static int startEmptyColumns=0;
static bool break2(CoinPresolveMatrix *prob)
{
int droppedRows = prob->countEmptyRows() - startEmptyRows ;
int droppedColumns = prob->countEmptyCols() - startEmptyColumns;
startEmptyRows=prob->countEmptyRows();
startEmptyColumns=prob->countEmptyCols();
printf("Dropped %d rows and %d columns - current empty %d, %d\n",droppedRows,
droppedColumns,startEmptyRows,startEmptyColumns);
counter--;
if (!counter) {
printf("skipping next and all\n");
}
return (counter<=0);
}
#define possibleBreak if (break2(prob)) break
#define possibleSkip if (!break2(prob))
#else
#define possibleBreak
#define possibleSkip
#endif
#define SOME_PRESOLVE_DETAIL
#ifndef SOME_PRESOLVE_DETAIL
#define printProgress(x,y) {}
#else
#define printProgress(x,y) {if ((presolveActions_ & 0x80000000) != 0) \
printf("%c loop %d %d empty rows, %d empty columns\n",x,y,prob->countEmptyRows(), \
prob->countEmptyCols());}
#endif
// This is the presolve loop.
// It is a separate virtual function so that it can be easily
// customized by subclassing CoinPresolve.
const CoinPresolveAction *ClpPresolve::presolve(CoinPresolveMatrix *prob)
{
// Messages
CoinMessages messages = CoinMessage(prob->messages().language());
paction_ = 0;
prob->maxSubstLevel_ = 3 ;
#ifndef PRESOLVE_DETAIL
if (prob->tuning_) {
#endif
int numberEmptyRows=0;
for ( int i=0;i<prob->nrows_;i++) {
if (!prob->hinrow_[i]) {
PRESOLVE_DETAIL_PRINT(printf("pre_empty row %d\n",i));
//printf("pre_empty row %d\n",i);
numberEmptyRows++;
}
}
int numberEmptyCols=0;
for ( int i=0;i<prob->ncols_;i++) {
if (!prob->hincol_[i]) {
PRESOLVE_DETAIL_PRINT(printf("pre_empty col %d\n",i));
//printf("pre_empty col %d\n",i);
numberEmptyCols++;
}
}
printf("CoinPresolve initial state %d empty rows and %d empty columns\n",
numberEmptyRows,numberEmptyCols);
#ifndef PRESOLVE_DETAIL
}
#endif
prob->status_ = 0; // say feasible
printProgress('A',0);
paction_ = make_fixed(prob, paction_);
paction_ = testRedundant(prob,paction_) ;
printProgress('B',0);
// if integers then switch off dual stuff
// later just do individually
bool doDualStuff = (presolvedModel_->integerInformation() == NULL);
// but allow in some cases
if ((presolveActions_ & 512) != 0)
doDualStuff = true;
if (prob->anyProhibited())
doDualStuff = false;
if (!doDual())
doDualStuff = false;
#if PRESOLVE_CONSISTENCY
// presolve_links_ok(prob->rlink_, prob->mrstrt_, prob->hinrow_, prob->nrows_);
presolve_links_ok(prob, false, true) ;
#endif
if (!prob->status_) {
bool slackSingleton = doSingletonColumn();
slackSingleton = true;
const bool slackd = doSingleton();
const bool doubleton = doDoubleton();
const bool tripleton = doTripleton();
//#define NO_FORCING
#ifndef NO_FORCING
const bool forcing = doForcing();
#endif
const bool ifree = doImpliedFree();
const bool zerocost = doTighten();
const bool dupcol = doDupcol();
const bool duprow = doDuprow();
const bool dual = doDualStuff;
// some things are expensive so just do once (normally)
int i;
// say look at all
if (!prob->anyProhibited()) {
for (i = 0; i < nrows_; i++)
prob->rowsToDo_[i] = i;
prob->numberRowsToDo_ = nrows_;
for (i = 0; i < ncols_; i++)
prob->colsToDo_[i] = i;
prob->numberColsToDo_ = ncols_;
} else {
// some stuff must be left alone
prob->numberRowsToDo_ = 0;
for (i = 0; i < nrows_; i++)
if (!prob->rowProhibited(i))
prob->rowsToDo_[prob->numberRowsToDo_++] = i;
prob->numberColsToDo_ = 0;
for (i = 0; i < ncols_; i++)
if (!prob->colProhibited(i))
prob->colsToDo_[prob->numberColsToDo_++] = i;
}
// transfer costs (may want to do it in OsiPresolve)
// need a transfer back at end of postsolve transferCosts(prob);
int iLoop;
#if PRESOLVE_DEBUG
check_sol(prob, 1.0e0);
#endif
if (dupcol) {
// maybe allow integer columns to be checked
if ((presolveActions_ & 512) != 0)
prob->setPresolveOptions(prob->presolveOptions() | 1);
possibleSkip;
paction_ = dupcol_action::presolve(prob, paction_);
printProgress('C',0);
}
#ifdef ABC_INHERIT
if (doTwoxTwo()) {
possibleSkip;
paction_ = twoxtwo_action::presolve(prob, paction_);
}
#endif
if (duprow) {
possibleSkip;
if (doTwoxTwo()) {
int nTightened=tightenDoubletons2(prob);
if (nTightened)
PRESOLVE_DETAIL_PRINT(printf("%d doubletons tightened\n",
nTightened));
}
paction_ = duprow_action::presolve(prob, paction_);
printProgress('D',0);
}
if (doGubrow()) {
possibleSkip;
paction_ = gubrow_action::presolve(prob, paction_);
printProgress('E',0);
}
if ((presolveActions_ & 16384) != 0)
prob->setPresolveOptions(prob->presolveOptions() | 16384);
// For inaccurate data in implied free
if ((presolveActions_ & 1024) != 0)
prob->setPresolveOptions(prob->presolveOptions() | 0x20000);
// Check number rows dropped
int lastDropped = 0;
prob->pass_ = 0;
#ifdef ABC_INHERIT
int numberRowsStart=nrows_-prob->countEmptyRows();
int numberColumnsStart=ncols_-prob->countEmptyCols();
int numberRowsLeft=numberRowsStart;
int numberColumnsLeft=numberColumnsStart;
bool lastPassWasGood=true;
#if ABC_NORMAL_DEBUG
printf("Original rows,columns %d,%d starting first pass with %d,%d\n",
nrows_,ncols_,numberRowsLeft,numberColumnsLeft);
#endif
#endif
if (numberPasses_<=5)
prob->presolveOptions_ |= 0x10000; // say more lightweight
for (iLoop = 0; iLoop < numberPasses_; iLoop++) {
// See if we want statistics
if ((presolveActions_ & 0x80000000) != 0)
printf("Starting major pass %d after %g seconds with %d rows, %d columns\n", iLoop + 1, CoinCpuTime() - prob->startTime_,
nrows_-prob->countEmptyRows(),
ncols_-prob->countEmptyCols());
#ifdef PRESOLVE_SUMMARY
printf("Starting major pass %d\n", iLoop + 1);
#endif
const CoinPresolveAction * const paction0 = paction_;
// look for substitutions with no fill
//#define IMPLIED 3
#ifdef IMPLIED
int fill_level = 3;
#define IMPLIED2 99
#if IMPLIED!=3
#if IMPLIED>2&&IMPLIED<11
fill_level = IMPLIED;
COIN_DETAIL_PRINT(printf("** fill_level == %d !\n", fill_level));
#endif
#if IMPLIED>11&&IMPLIED<21
fill_level = -(IMPLIED - 10);
COIN_DETAIL_PRINT(printf("** fill_level == %d !\n", fill_level));
#endif
#endif
#else
int fill_level = prob->maxSubstLevel_;
#endif
int whichPass = 0;
while (1) {
whichPass++;
prob->pass_++;
const CoinPresolveAction * const paction1 = paction_;
if (slackd) {
bool notFinished = true;
while (notFinished) {
possibleBreak;
paction_ = slack_doubleton_action::presolve(prob, paction_,
notFinished);
}
printProgress('F',iLoop+1);
if (prob->status_)
break;
}
if (dual && whichPass == 1) {
// this can also make E rows so do one bit here
possibleBreak;
paction_ = remove_dual_action::presolve(prob, paction_);
if (prob->status_)
break;
printProgress('G',iLoop+1);
}
if (doubleton) {
possibleBreak;
paction_ = doubleton_action::presolve(prob, paction_);
if (prob->status_)
break;
printProgress('H',iLoop+1);
}
if (tripleton) {
possibleBreak;
paction_ = tripleton_action::presolve(prob, paction_);
if (prob->status_)
break;
printProgress('I',iLoop+1);
}
if (zerocost) {
possibleBreak;
paction_ = do_tighten_action::presolve(prob, paction_);
if (prob->status_)
break;
printProgress('J',iLoop+1);
}
#ifndef NO_FORCING
if (forcing) {
possibleBreak;
paction_ = forcing_constraint_action::presolve(prob, paction_);
if (prob->status_)
break;
printProgress('K',iLoop+1);
}
#endif
if (ifree && (whichPass % 5) == 1) {
possibleBreak;
paction_ = implied_free_action::presolve(prob, paction_, fill_level);
if (prob->status_)
break;
printProgress('L',iLoop+1);
}
#if PRESOLVE_DEBUG
check_sol(prob, 1.0e0);
#endif
#if PRESOLVE_CONSISTENCY
// presolve_links_ok(prob->rlink_, prob->mrstrt_, prob->hinrow_,
// prob->nrows_);
presolve_links_ok(prob, false, true) ;
#endif
//#if PRESOLVE_DEBUG
// presolve_no_zeros(prob->mcstrt_, prob->colels_, prob->hincol_,
// prob->ncols_);
//#endif
//#if PRESOLVE_CONSISTENCY
// prob->consistent();
//#endif
#if PRESOLVE_CONSISTENCY
presolve_no_zeros(prob, true, false) ;
presolve_consistent(prob, true) ;
#endif
{
// set up for next pass
// later do faster if many changes i.e. memset and memcpy
const int * count = prob->hinrow_;
const int * nextToDo = prob->nextRowsToDo_;
int * toDo = prob->rowsToDo_;
int nNext = prob->numberNextRowsToDo_;
int n = 0;
for (int i = 0; i < nNext; i++) {
int index = nextToDo[i];
prob->unsetRowChanged(index);
if (count[index])
toDo[n++] = index;
}
prob->numberRowsToDo_ = n;
prob->numberNextRowsToDo_ = 0;
count = prob->hincol_;
nextToDo = prob->nextColsToDo_;
toDo = prob->colsToDo_;
nNext = prob->numberNextColsToDo_;
n = 0;
for (int i = 0; i < nNext; i++) {
int index = nextToDo[i];
prob->unsetColChanged(index);
if (count[index])
toDo[n++] = index;
}
prob->numberColsToDo_ = n;
prob->numberNextColsToDo_ = 0;
}
if (paction_ == paction1 && fill_level > 0)
break;
}
// say look at all
int i;
if (!prob->anyProhibited()) {
const int * count = prob->hinrow_;
int * toDo = prob->rowsToDo_;
int n = 0;
for (int i = 0; i < nrows_; i++) {
prob->unsetRowChanged(i);
if (count[i])
toDo[n++] = i;
}
prob->numberRowsToDo_ = n;
prob->numberNextRowsToDo_ = 0;
count = prob->hincol_;
toDo = prob->colsToDo_;
n = 0;
for (int i = 0; i < ncols_; i++) {
prob->unsetColChanged(i);
if (count[i])
toDo[n++] = i;
}
prob->numberColsToDo_ = n;
prob->numberNextColsToDo_ = 0;
} else {
// some stuff must be left alone
prob->numberRowsToDo_ = 0;
for (i = 0; i < nrows_; i++)
if (!prob->rowProhibited(i))
prob->rowsToDo_[prob->numberRowsToDo_++] = i;
prob->numberColsToDo_ = 0;
for (i = 0; i < ncols_; i++)
if (!prob->colProhibited(i))
prob->colsToDo_[prob->numberColsToDo_++] = i;
}
// now expensive things
// this caused world.mps to run into numerical difficulties
#ifdef PRESOLVE_SUMMARY
printf("Starting expensive\n");
#endif
if (dual) {
int itry;
for (itry = 0; itry < 5; itry++) {
possibleBreak;
paction_ = remove_dual_action::presolve(prob, paction_);
if (prob->status_)
break;
printProgress('M',iLoop+1);
const CoinPresolveAction * const paction2 = paction_;
if (ifree) {
#ifdef IMPLIED
#if IMPLIED2 ==0
int fill_level = 0; // switches off substitution
#elif IMPLIED2!=99
int fill_level = IMPLIED2;
#endif
#endif
if ((itry & 1) == 0) {
possibleBreak;
paction_ = implied_free_action::presolve(prob, paction_, fill_level);
}
if (prob->status_)
break;
printProgress('N',iLoop+1);
}
if (paction_ == paction2)
break;
}
} else if (ifree) {
// just free
#ifdef IMPLIED
#if IMPLIED2 ==0
int fill_level = 0; // switches off substitution
#elif IMPLIED2!=99
int fill_level = IMPLIED2;
#endif
#endif
possibleBreak;
paction_ = implied_free_action::presolve(prob, paction_, fill_level);
if (prob->status_)
break;
printProgress('O',iLoop+1);
}
#if PRESOLVE_DEBUG
check_sol(prob, 1.0e0);
#endif
if (dupcol) {
// maybe allow integer columns to be checked
if ((presolveActions_ & 512) != 0)
prob->setPresolveOptions(prob->presolveOptions() | 1);
possibleBreak;
paction_ = dupcol_action::presolve(prob, paction_);
if (prob->status_)
break;
printProgress('P',iLoop+1);
}
#if PRESOLVE_DEBUG
check_sol(prob, 1.0e0);
#endif
if (duprow) {
possibleBreak;
paction_ = duprow_action::presolve(prob, paction_);
if (prob->status_)
break;
printProgress('Q',iLoop+1);
}
// Marginally slower on netlib if this call is enabled.
// paction_ = testRedundant(prob,paction_) ;
#if PRESOLVE_DEBUG
check_sol(prob, 1.0e0);
#endif
bool stopLoop = false;
{
int * hinrow = prob->hinrow_;
int numberDropped = 0;
for (i = 0; i < nrows_; i++)
if (!hinrow[i])
numberDropped++;
prob->messageHandler()->message(COIN_PRESOLVE_PASS,
messages)
<< numberDropped << iLoop + 1
<< CoinMessageEol;
//printf("%d rows dropped after pass %d\n",numberDropped,
// iLoop+1);
if (numberDropped == lastDropped)
stopLoop = true;
else
lastDropped = numberDropped;
}
// Do this here as not very loopy
if (slackSingleton) {
// On most passes do not touch costed slacks
if (paction_ != paction0 && !stopLoop) {
possibleBreak;
paction_ = slack_singleton_action::presolve(prob, paction_, NULL);
} else {
// do costed if Clp (at end as ruins rest of presolve)
possibleBreak;
paction_ = slack_singleton_action::presolve(prob, paction_, rowObjective_);
stopLoop = true;
}
printProgress('R',iLoop+1);
}
#if PRESOLVE_DEBUG
check_sol(prob, 1.0e0);
#endif
if (paction_ == paction0 || stopLoop)
break;
#ifdef ABC_INHERIT
// see whether to stop anyway
int numberRowsNow=nrows_-prob->countEmptyRows();
int numberColumnsNow=ncols_-prob->countEmptyCols();
#if ABC_NORMAL_DEBUG
printf("Original rows,columns %d,%d - last %d,%d end of pass %d has %d,%d\n",
nrows_,ncols_,numberRowsLeft,numberColumnsLeft,iLoop+1,numberRowsNow,
numberColumnsNow);
#endif
int rowsDeleted=numberRowsLeft-numberRowsNow;
int columnsDeleted=numberColumnsLeft-numberColumnsNow;
if (iLoop>15) {
if (rowsDeleted*100<numberRowsStart&&
columnsDeleted*100<numberColumnsStart)
break;
lastPassWasGood=true;
} else if (rowsDeleted*100<numberRowsStart&&rowsDeleted<500&&
columnsDeleted*100<numberColumnsStart&&columnsDeleted<500) {
if (!lastPassWasGood)
break;
else
lastPassWasGood=false;
} else {
lastPassWasGood=true;
}
numberRowsLeft=numberRowsNow;
numberColumnsLeft=numberColumnsNow;
#endif
}
}
prob->presolveOptions_ &= ~0x10000;
if (!prob->status_) {
paction_ = drop_zero_coefficients(prob, paction_);
#if PRESOLVE_DEBUG
check_sol(prob, 1.0e0);
#endif
paction_ = drop_empty_cols_action::presolve(prob, paction_);
paction_ = drop_empty_rows_action::presolve(prob, paction_);
#if PRESOLVE_DEBUG
check_sol(prob, 1.0e0);
#endif
}
if (prob->status_) {
if (prob->status_ == 1)
prob->messageHandler()->message(COIN_PRESOLVE_INFEAS,
messages)
<< prob->feasibilityTolerance_
<< CoinMessageEol;
else if (prob->status_ == 2)
prob->messageHandler()->message(COIN_PRESOLVE_UNBOUND,
messages)
<< CoinMessageEol;
else
prob->messageHandler()->message(COIN_PRESOLVE_INFEASUNBOUND,
messages)
<< CoinMessageEol;
// get rid of data
destroyPresolve();
}
return (paction_);
}
void check_djs(CoinPostsolveMatrix *prob);
// We could have implemented this by having each postsolve routine
// directly call the next one, but this may make it easier to add debugging checks.
void ClpPresolve::postsolve(CoinPostsolveMatrix &prob)
{
{
// Check activities
double *colels = prob.colels_;
int *hrow = prob.hrow_;
CoinBigIndex *mcstrt = prob.mcstrt_;
int *hincol = prob.hincol_;
int *link = prob.link_;
int ncols = prob.ncols_;
char *cdone = prob.cdone_;
double * csol = prob.sol_;
int nrows = prob.nrows_;
int colx;
double * rsol = prob.acts_;
memset(rsol, 0, nrows * sizeof(double));
for (colx = 0; colx < ncols; ++colx) {
if (cdone[colx]) {
CoinBigIndex k = mcstrt[colx];
int nx = hincol[colx];
double solutionValue = csol[colx];
for (int i = 0; i < nx; ++i) {
int row = hrow[k];
double coeff = colels[k];
k = link[k];
rsol[row] += solutionValue * coeff;
}
}
}
}
if (prob.maxmin_<0) {
//for (int i=0;i<presolvedModel_->numberRows();i++)
//prob.rowduals_[i]=-prob.rowduals_[i];
for (int i=0;i<ncols_;i++) {
prob.cost_[i]=-prob.cost_[i];
//prob.rcosts_[i]=-prob.rcosts_[i];
}
prob.maxmin_=1.0;
}
const CoinPresolveAction *paction = paction_;
//#define PRESOLVE_DEBUG 1
#if PRESOLVE_DEBUG
// Check only works after first one
int checkit = -1;
#endif
while (paction) {
#if PRESOLVE_DEBUG
printf("POSTSOLVING %s\n", paction->name());
#endif
paction->postsolve(&prob);
#if PRESOLVE_DEBUG
# if 0
/*
This check fails (on exmip1 (!) in osiUnitTest) because clp
enters postsolve with a solution that seems to have incorrect
status for a logical. You can see similar behaviour with
column status --- incorrect entering postsolve.
-- lh, 111207 --
*/
{
int nr = 0;
int i;
for (i = 0; i < prob.nrows_; i++) {
if ((prob.rowstat_[i] & 7) == 1) {
nr++;
} else if ((prob.rowstat_[i] & 7) == 2) {
// at ub
assert (prob.acts_[i] > prob.rup_[i] - 1.0e-6);
} else if ((prob.rowstat_[i] & 7) == 3) {
// at lb
assert (prob.acts_[i] < prob.rlo_[i] + 1.0e-6);
}
}
int nc = 0;
for (i = 0; i < prob.ncols_; i++) {
if ((prob.colstat_[i] & 7) == 1)
nc++;
}
printf("%d rows (%d basic), %d cols (%d basic)\n", prob.nrows_, nr, prob.ncols_, nc);
}
# endif // if 0
checkit++;
if (prob.colstat_ && checkit > 0) {
presolve_check_nbasic(&prob) ;
presolve_check_sol(&prob, 2, 2, 1) ;
}
#endif
paction = paction->next;
#if PRESOLVE_DEBUG
// check_djs(&prob);
if (checkit > 0)
presolve_check_reduced_costs(&prob) ;
#endif
}
#if PRESOLVE_DEBUG
if (prob.colstat_) {
presolve_check_nbasic(&prob) ;
presolve_check_sol(&prob, 2, 2, 1) ;
}
#endif
#undef PRESOLVE_DEBUG
#if 0 && PRESOLVE_DEBUG
for (i = 0; i < ncols0; i++) {
if (!cdone[i]) {
printf("!cdone[%d]\n", i);
abort();
}
}
for (i = 0; i < nrows0; i++) {
if (!rdone[i]) {
printf("!rdone[%d]\n", i);
abort();
}
}
for (i = 0; i < ncols0; i++) {
if (sol[i] < -1e10 || sol[i] > 1e10)
printf("!!!%d %g\n", i, sol[i]);
}
#endif
#if 0 && PRESOLVE_DEBUG
// debug check: make sure we ended up with same original matrix
{
int identical = 1;
for (int i = 0; i < ncols0; i++) {
PRESOLVEASSERT(hincol[i] == &prob->mcstrt0[i+1] - &prob->mcstrt0[i]);
CoinBigIndex kcs0 = &prob->mcstrt0[i];
CoinBigIndex kcs = mcstrt[i];
int n = hincol[i];
for (int k = 0; k < n; k++) {
CoinBigIndex k1 = presolve_find_row1(&prob->hrow0[kcs0+k], kcs, kcs + n, hrow);
if (k1 == kcs + n) {
printf("ROW %d NOT IN COL %d\n", &prob->hrow0[kcs0+k], i);
abort();
}
if (colels[k1] != &prob->dels0[kcs0+k])
printf("BAD COLEL[%d %d %d]: %g\n",
k1, i, &prob->hrow0[kcs0+k], colels[k1] - &prob->dels0[kcs0+k]);
if (kcs0 + k != k1)
identical = 0;
}
}
printf("identical? %d\n", identical);
}
#endif
}
static inline double getTolerance(const ClpSimplex *si, ClpDblParam key)
{
double tol;
if (! si->getDblParam(key, tol)) {
CoinPresolveAction::throwCoinError("getDblParam failed",
"CoinPrePostsolveMatrix::CoinPrePostsolveMatrix");
}
return (tol);
}
// Assumptions:
// 1. nrows>=m.getNumRows()
// 2. ncols>=m.getNumCols()
//
// In presolve, these values are equal.
// In postsolve, they may be inequal, since the reduced problem
// may be smaller, but we need room for the large problem.
// ncols may be larger than si.getNumCols() in postsolve,
// this at that point si will be the reduced problem,
// but we need to reserve enough space for the original problem.
CoinPrePostsolveMatrix::CoinPrePostsolveMatrix(const ClpSimplex * si,
int ncols_in,
int nrows_in,
CoinBigIndex nelems_in,
double bulkRatio)
: ncols_(si->getNumCols()),
nrows_(si->getNumRows()),
nelems_(si->getNumElements()),
ncols0_(ncols_in),
nrows0_(nrows_in),
bulkRatio_(bulkRatio),
mcstrt_(new CoinBigIndex[ncols_in+1]),
hincol_(new int[ncols_in+1]),
cost_(new double[ncols_in]),
clo_(new double[ncols_in]),
cup_(new double[ncols_in]),
rlo_(new double[nrows_in]),
rup_(new double[nrows_in]),
originalColumn_(new int[ncols_in]),
originalRow_(new int[nrows_in]),
ztolzb_(getTolerance(si, ClpPrimalTolerance)),
ztoldj_(getTolerance(si, ClpDualTolerance)),
maxmin_(si->getObjSense()),
sol_(NULL),
rowduals_(NULL),
acts_(NULL),
rcosts_(NULL),
colstat_(NULL),
rowstat_(NULL),
handler_(NULL),
defaultHandler_(false),
messages_()
{
bulk0_ = static_cast<CoinBigIndex> (bulkRatio_ * nelems_in);
hrow_ = new int [bulk0_];
colels_ = new double[bulk0_];
si->getDblParam(ClpObjOffset, originalOffset_);
int ncols = si->getNumCols();
int nrows = si->getNumRows();
setMessageHandler(si->messageHandler()) ;
ClpDisjointCopyN(si->getColLower(), ncols, clo_);
ClpDisjointCopyN(si->getColUpper(), ncols, cup_);
//ClpDisjointCopyN(si->getObjCoefficients(), ncols, cost_);
double offset;
ClpDisjointCopyN(si->objectiveAsObject()->gradient(si, si->getColSolution(), offset, true), ncols, cost_);
ClpDisjointCopyN(si->getRowLower(), nrows, rlo_);
ClpDisjointCopyN(si->getRowUpper(), nrows, rup_);
int i;
for (i = 0; i < ncols_in; i++)
originalColumn_[i] = i;
for (i = 0; i < nrows_in; i++)
originalRow_[i] = i;
sol_ = NULL;
rowduals_ = NULL;
acts_ = NULL;
rcosts_ = NULL;
colstat_ = NULL;
rowstat_ = NULL;
}
// I am not familiar enough with CoinPackedMatrix to be confident
// that I will implement a row-ordered version of toColumnOrderedGapFree
// properly.
static bool isGapFree(const CoinPackedMatrix& matrix)
{
const CoinBigIndex * start = matrix.getVectorStarts();
const int * length = matrix.getVectorLengths();
int i = matrix.getSizeVectorLengths() - 1;
// Quick check
if (matrix.getNumElements() == start[i]) {
return true;
} else {
for (i = matrix.getSizeVectorLengths() - 1; i >= 0; --i) {
if (start[i+1] - start[i] != length[i])
break;
}
return (! (i >= 0));
}
}
#if PRESOLVE_DEBUG
static void matrix_bounds_ok(const double *lo, const double *up, int n)
{
int i;
for (i = 0; i < n; i++) {
PRESOLVEASSERT(lo[i] <= up[i]);
PRESOLVEASSERT(lo[i] < PRESOLVE_INF);
PRESOLVEASSERT(-PRESOLVE_INF < up[i]);
}
}
#endif
CoinPresolveMatrix::CoinPresolveMatrix(int ncols0_in,
double /*maxmin*/,
// end prepost members
ClpSimplex * si,
// rowrep
int nrows_in,
CoinBigIndex nelems_in,
bool doStatus,
double nonLinearValue,
double bulkRatio) :
CoinPrePostsolveMatrix(si,
ncols0_in, nrows_in, nelems_in, bulkRatio),
clink_(new presolvehlink[ncols0_in+1]),
rlink_(new presolvehlink[nrows_in+1]),
dobias_(0.0),
// temporary init
integerType_(new unsigned char[ncols0_in]),
tuning_(false),
startTime_(0.0),
feasibilityTolerance_(0.0),
status_(-1),
colsToDo_(new int [ncols0_in]),
numberColsToDo_(0),
nextColsToDo_(new int[ncols0_in]),
numberNextColsToDo_(0),
rowsToDo_(new int [nrows_in]),
numberRowsToDo_(0),
nextRowsToDo_(new int[nrows_in]),
numberNextRowsToDo_(0),
presolveOptions_(0)
{
const int bufsize = bulk0_;
nrows_ = si->getNumRows() ;
// Set up change bits etc
rowChanged_ = new unsigned char[nrows_];
memset(rowChanged_, 0, nrows_);
colChanged_ = new unsigned char[ncols_];
memset(colChanged_, 0, ncols_);
CoinPackedMatrix * m = si->matrix();
// The coefficient matrix is a big hunk of stuff.
// Do the copy here to try to avoid running out of memory.
const CoinBigIndex * start = m->getVectorStarts();
const int * row = m->getIndices();
const double * element = m->getElements();
int icol, nel = 0;
mcstrt_[0] = 0;
ClpDisjointCopyN(m->getVectorLengths(), ncols_, hincol_);
if (si->getObjSense() < 0.0) {
for (int i=0;i<ncols_;i++)
cost_[i]=-cost_[i];
maxmin_=1.0;
}
for (icol = 0; icol < ncols_; icol++) {
CoinBigIndex j;
for (j = start[icol]; j < start[icol] + hincol_[icol]; j++) {
hrow_[nel] = row[j];
if (fabs(element[j]) > ZTOLDP)
colels_[nel++] = element[j];
}
mcstrt_[icol+1] = nel;
hincol_[icol] = nel - mcstrt_[icol];
}
// same thing for row rep
CoinPackedMatrix * mRow = new CoinPackedMatrix();
mRow->setExtraGap(0.0);
mRow->setExtraMajor(0.0);
mRow->reverseOrderedCopyOf(*m);
//mRow->removeGaps();
//mRow->setExtraGap(0.0);
// Now get rid of matrix
si->createEmptyMatrix();
double * el = mRow->getMutableElements();
int * ind = mRow->getMutableIndices();
CoinBigIndex * strt = mRow->getMutableVectorStarts();
int * len = mRow->getMutableVectorLengths();
// Do carefully to save memory
rowels_ = new double[bulk0_];
ClpDisjointCopyN(el, nelems_, rowels_);
mRow->nullElementArray();
delete [] el;
hcol_ = new int[bulk0_];
ClpDisjointCopyN(ind, nelems_, hcol_);
mRow->nullIndexArray();
delete [] ind;
mrstrt_ = new CoinBigIndex[nrows_in+1];
ClpDisjointCopyN(strt, nrows_, mrstrt_);
mRow->nullStartArray();
mrstrt_[nrows_] = nelems_;
delete [] strt;
hinrow_ = new int[nrows_in+1];
ClpDisjointCopyN(len, nrows_, hinrow_);
if (nelems_ > nel) {
nelems_ = nel;
// Clean any small elements
int irow;
nel = 0;
CoinBigIndex start = 0;
for (irow = 0; irow < nrows_; irow++) {
CoinBigIndex j;
for (j = start; j < start + hinrow_[irow]; j++) {
hcol_[nel] = hcol_[j];
if (fabs(rowels_[j]) > ZTOLDP)
rowels_[nel++] = rowels_[j];
}
start = mrstrt_[irow+1];
mrstrt_[irow+1] = nel;
hinrow_[irow] = nel - mrstrt_[irow];
}
}
delete mRow;
if (si->integerInformation()) {
CoinMemcpyN(reinterpret_cast<unsigned char *> (si->integerInformation()), ncols_, integerType_);
} else {
ClpFillN<unsigned char>(integerType_, ncols_, static_cast<unsigned char> (0));
}
#ifndef SLIM_CLP
#ifndef NO_RTTI
ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(si->objectiveAsObject()));
#else
ClpQuadraticObjective * quadraticObj = NULL;
if (si->objectiveAsObject()->type() == 2)
quadraticObj = (static_cast< ClpQuadraticObjective*>(si->objectiveAsObject()));
#endif
#endif
// Set up prohibited bits if needed
if (nonLinearValue) {
anyProhibited_ = true;
for (icol = 0; icol < ncols_; icol++) {
int j;
bool nonLinearColumn = false;
if (cost_[icol] == nonLinearValue)
nonLinearColumn = true;
for (j = mcstrt_[icol]; j < mcstrt_[icol+1]; j++) {
if (colels_[j] == nonLinearValue) {
nonLinearColumn = true;
setRowProhibited(hrow_[j]);
}
}
if (nonLinearColumn)
setColProhibited(icol);
}
#ifndef SLIM_CLP
} else if (quadraticObj) {
CoinPackedMatrix * quadratic = quadraticObj->quadraticObjective();
//const int * columnQuadratic = quadratic->getIndices();
//const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
const int * columnQuadraticLength = quadratic->getVectorLengths();
//double * quadraticElement = quadratic->getMutableElements();
int numberColumns = quadratic->getNumCols();
anyProhibited_ = true;
for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
if (columnQuadraticLength[iColumn]) {
setColProhibited(iColumn);
//printf("%d prohib\n",iColumn);
}
}
#endif
} else {
anyProhibited_ = false;
}
if (doStatus) {
// allow for status and solution
sol_ = new double[ncols_];
CoinMemcpyN(si->primalColumnSolution(), ncols_, sol_);;
acts_ = new double [nrows_];
CoinMemcpyN(si->primalRowSolution(), nrows_, acts_);
if (!si->statusArray())
si->createStatus();
colstat_ = new unsigned char [nrows_+ncols_];
CoinMemcpyN(si->statusArray(), (nrows_ + ncols_), colstat_);
rowstat_ = colstat_ + ncols_;
}
// the original model's fields are now unneeded - free them
si->resize(0, 0);
#if PRESOLVE_DEBUG
matrix_bounds_ok(rlo_, rup_, nrows_);
matrix_bounds_ok(clo_, cup_, ncols_);
#endif
#if 0
for (i = 0; i < nrows; ++i)
printf("NR: %6d\n", hinrow[i]);
for (int i = 0; i < ncols; ++i)
printf("NC: %6d\n", hincol[i]);
#endif
presolve_make_memlists(/*mcstrt_,*/ hincol_, clink_, ncols_);
presolve_make_memlists(/*mrstrt_,*/ hinrow_, rlink_, nrows_);
// this allows last col/row to expand up to bufsize-1 (22);
// this must come after the calls to presolve_prefix
mcstrt_[ncols_] = bufsize - 1;
mrstrt_[nrows_] = bufsize - 1;
// Allocate useful arrays
initializeStuff();
#if PRESOLVE_CONSISTENCY
//consistent(false);
presolve_consistent(this, false) ;
#endif
}
// avoid compiler warnings
#if PRESOLVE_SUMMARY > 0
void CoinPresolveMatrix::update_model(ClpSimplex * si,
int nrows0, int ncols0,
CoinBigIndex nelems0)
#else
void CoinPresolveMatrix::update_model(ClpSimplex * si,
int /*nrows0*/,
int /*ncols0*/,
CoinBigIndex /*nelems0*/)
#endif
{
if (si->getObjSense() < 0.0) {
for (int i=0;i<ncols_;i++)
cost_[i]=-cost_[i];
dobias_=-dobias_;
}
si->loadProblem(ncols_, nrows_, mcstrt_, hrow_, colels_, hincol_,
clo_, cup_, cost_, rlo_, rup_);
//delete [] si->integerInformation();
int numberIntegers = 0;
for (int i = 0; i < ncols_; i++) {
if (integerType_[i])
numberIntegers++;
}
if (numberIntegers)
si->copyInIntegerInformation(reinterpret_cast<const char *> (integerType_));
else
si->copyInIntegerInformation(NULL);
#if PRESOLVE_SUMMARY
printf("NEW NCOL/NROW/NELS: %d(-%d) %d(-%d) %d(-%d)\n",
ncols_, ncols0 - ncols_,
nrows_, nrows0 - nrows_,
si->getNumElements(), nelems0 - si->getNumElements());
#endif
si->setDblParam(ClpObjOffset, originalOffset_ - dobias_);
if (si->getObjSense() < 0.0) {
// put back
for (int i=0;i<ncols_;i++)
cost_[i]=-cost_[i];
dobias_=-dobias_;
maxmin_=-1.0;
}
}
//////////////// POSTSOLVE
CoinPostsolveMatrix::CoinPostsolveMatrix(ClpSimplex* si,
int ncols0_in,
int nrows0_in,
CoinBigIndex nelems0,
double maxmin,
// end prepost members
double *sol_in,
double *acts_in,
unsigned char *colstat_in,
unsigned char *rowstat_in) :
CoinPrePostsolveMatrix(si,
ncols0_in, nrows0_in, nelems0, 2.0),
free_list_(0),
// link, free_list, maxlink
maxlink_(bulk0_),
link_(new int[/*maxlink*/ bulk0_]),
cdone_(new char[ncols0_]),
rdone_(new char[nrows0_in])
{
bulk0_ = maxlink_ ;
nrows_ = si->getNumRows() ;
ncols_ = si->getNumCols() ;
sol_ = sol_in;
rowduals_ = NULL;
acts_ = acts_in;
rcosts_ = NULL;
colstat_ = colstat_in;
rowstat_ = rowstat_in;
// this is the *reduced* model, which is probably smaller
int ncols1 = ncols_ ;
int nrows1 = nrows_ ;
const CoinPackedMatrix * m = si->matrix();
const CoinBigIndex nelemsr = m->getNumElements();
if (m->getNumElements() && !isGapFree(*m)) {
// Odd - gaps
CoinPackedMatrix mm(*m);
mm.removeGaps();
mm.setExtraGap(0.0);
ClpDisjointCopyN(mm.getVectorStarts(), ncols1, mcstrt_);
CoinZeroN(mcstrt_ + ncols1, ncols0_ - ncols1);
mcstrt_[ncols1] = nelems0; // ?? (should point to end of bulk store -- lh --)
ClpDisjointCopyN(mm.getVectorLengths(), ncols1, hincol_);
ClpDisjointCopyN(mm.getIndices(), nelemsr, hrow_);
ClpDisjointCopyN(mm.getElements(), nelemsr, colels_);
} else {
// No gaps
ClpDisjointCopyN(m->getVectorStarts(), ncols1, mcstrt_);
CoinZeroN(mcstrt_ + ncols1, ncols0_ - ncols1);
mcstrt_[ncols1] = nelems0; // ?? (should point to end of bulk store -- lh --)
ClpDisjointCopyN(m->getVectorLengths(), ncols1, hincol_);
ClpDisjointCopyN(m->getIndices(), nelemsr, hrow_);
ClpDisjointCopyN(m->getElements(), nelemsr, colels_);
}
#if 0 && PRESOLVE_DEBUG
presolve_check_costs(model, &colcopy);
#endif
// This determines the size of the data structure that contains
// the matrix being postsolved. Links are taken from the free_list
// to recreate matrix entries that were presolved away,
// and links are added to the free_list when entries created during
// presolve are discarded. There is never a need to gc this list.
// Naturally, it should contain
// exactly nelems0 entries "in use" when postsolving is done,
// but I don't know whether the matrix could temporarily get
// larger during postsolving. Substitution into more than two
// rows could do that, in principle. I am being very conservative
// here by reserving much more than the amount of space I probably need.
// If this guess is wrong, check_free_list may be called.
// int bufsize = 2*nelems0;
memset(cdone_, -1, ncols0_);
memset(rdone_, -1, nrows0_);
rowduals_ = new double[nrows0_];
ClpDisjointCopyN(si->getRowPrice(), nrows1, rowduals_);
rcosts_ = new double[ncols0_];
ClpDisjointCopyN(si->getReducedCost(), ncols1, rcosts_);
if (maxmin < 0.0) {
// change so will look as if minimize
int i;
for (i = 0; i < nrows1; i++)
rowduals_[i] = - rowduals_[i];
for (i = 0; i < ncols1; i++) {
rcosts_[i] = - rcosts_[i];
}
}
//ClpDisjointCopyN(si->getRowUpper(), nrows1, rup_);
//ClpDisjointCopyN(si->getRowLower(), nrows1, rlo_);
ClpDisjointCopyN(si->getColSolution(), ncols1, sol_);
si->setDblParam(ClpObjOffset, originalOffset_);
for (int j = 0; j < ncols1; j++) {
#ifdef COIN_SLOW_PRESOLVE
if (hincol_[j]) {
#endif
CoinBigIndex kcs = mcstrt_[j];
CoinBigIndex kce = kcs + hincol_[j];
for (CoinBigIndex k = kcs; k < kce; ++k) {
link_[k] = k + 1;
}
link_[kce-1] = NO_LINK ;
#ifdef COIN_SLOW_PRESOLVE
}
#endif
}
{
int ml = maxlink_;
for (CoinBigIndex k = nelemsr; k < ml; ++k)
link_[k] = k + 1;
if (ml)
link_[ml-1] = NO_LINK;
}
free_list_ = nelemsr;
# if PRESOLVE_DEBUG || PRESOLVE_CONSISTENCY
/*
These are used to track the action of postsolve transforms during debugging.
*/
CoinFillN(cdone_, ncols1, PRESENT_IN_REDUCED) ;
CoinZeroN(cdone_ + ncols1, ncols0_in - ncols1) ;
CoinFillN(rdone_, nrows1, PRESENT_IN_REDUCED) ;
CoinZeroN(rdone_ + nrows1, nrows0_in - nrows1) ;
# endif
}
/* This is main part of Presolve */
ClpSimplex *
ClpPresolve::gutsOfPresolvedModel(ClpSimplex * originalModel,
double feasibilityTolerance,
bool keepIntegers,
int numberPasses,
bool dropNames,
bool doRowObjective,
const char * prohibitedRows,
const char * prohibitedColumns)
{
ncols_ = originalModel->getNumCols();
nrows_ = originalModel->getNumRows();
nelems_ = originalModel->getNumElements();
numberPasses_ = numberPasses;
double maxmin = originalModel->getObjSense();
originalModel_ = originalModel;
delete [] originalColumn_;
originalColumn_ = new int[ncols_];
delete [] originalRow_;
originalRow_ = new int[nrows_];
// and fill in case returns early
int i;
for (i = 0; i < ncols_; i++)
originalColumn_[i] = i;
for (i = 0; i < nrows_; i++)
originalRow_[i] = i;
delete [] rowObjective_;
if (doRowObjective) {
rowObjective_ = new double [nrows_];
memset(rowObjective_, 0, nrows_ * sizeof(double));
} else {
rowObjective_ = NULL;
}
// result is 0 - okay, 1 infeasible, -1 go round again, 2 - original model
int result = -1;
// User may have deleted - its their responsibility
presolvedModel_ = NULL;
// Messages
CoinMessages messages = originalModel->coinMessages();
// Only go round 100 times even if integer preprocessing
int totalPasses = 100;
while (result == -1) {
#ifndef CLP_NO_STD
// make new copy
if (saveFile_ == "") {
#endif
delete presolvedModel_;
#ifndef CLP_NO_STD
// So won't get names
int lengthNames = originalModel->lengthNames();
originalModel->setLengthNames(0);
#endif
presolvedModel_ = new ClpSimplex(*originalModel);
#ifndef CLP_NO_STD
originalModel->setLengthNames(lengthNames);
presolvedModel_->dropNames();
} else {
presolvedModel_ = originalModel;
if (dropNames)
presolvedModel_->dropNames();
}
#endif
// drop integer information if wanted
if (!keepIntegers)
presolvedModel_->deleteIntegerInformation();
totalPasses--;
double ratio = 2.0;
if (substitution_ > 3)
ratio = substitution_;
else if (substitution_ == 2)
ratio = 1.5;
CoinPresolveMatrix prob(ncols_,
maxmin,
presolvedModel_,
nrows_, nelems_, true, nonLinearValue_, ratio);
if (prohibitedRows) {
prob.setAnyProhibited();
for (int i=0;i<nrows_;i++) {
if (prohibitedRows[i])
prob.setRowProhibited(i);
}
}
if (prohibitedColumns) {
prob.setAnyProhibited();
for (int i=0;i<ncols_;i++) {
if (prohibitedColumns[i])
prob.setColProhibited(i);
}
}
prob.setMaximumSubstitutionLevel(substitution_);
if (doRowObjective)
memset(rowObjective_, 0, nrows_ * sizeof(double));
// See if we want statistics
if ((presolveActions_ & 0x80000000) != 0)
prob.statistics();
// make sure row solution correct
{
double *colels = prob.colels_;
int *hrow = prob.hrow_;
CoinBigIndex *mcstrt = prob.mcstrt_;
int *hincol = prob.hincol_;
int ncols = prob.ncols_;
double * csol = prob.sol_;
double * acts = prob.acts_;
int nrows = prob.nrows_;
int colx;
memset(acts, 0, nrows * sizeof(double));
for (colx = 0; colx < ncols; ++colx) {
double solutionValue = csol[colx];
for (int i = mcstrt[colx]; i < mcstrt[colx] + hincol[colx]; ++i) {
int row = hrow[i];
double coeff = colels[i];
acts[row] += solutionValue * coeff;
}
}
}
// move across feasibility tolerance
prob.feasibilityTolerance_ = feasibilityTolerance;
// Do presolve
paction_ = presolve(&prob);
// Get rid of useful arrays
prob.deleteStuff();
result = 0;
bool fixInfeasibility = (prob.presolveOptions_&16384)!=0;
bool hasSolution = (prob.presolveOptions_&32768)!=0;
if (prob.status_ == 0 && paction_ && (!hasSolution || !fixInfeasibility)) {
// Looks feasible but double check to see if anything slipped through
int n = prob.ncols_;
double * lo = prob.clo_;
double * up = prob.cup_;
int i;
for (i = 0; i < n; i++) {
if (up[i] < lo[i]) {
if (up[i] < lo[i] - feasibilityTolerance && !fixInfeasibility) {
// infeasible
prob.status_ = 1;
} else {
up[i] = lo[i];
}
}
}
n = prob.nrows_;
lo = prob.rlo_;
up = prob.rup_;
for (i = 0; i < n; i++) {
if (up[i] < lo[i]) {
if (up[i] < lo[i] - feasibilityTolerance && !fixInfeasibility) {
// infeasible
prob.status_ = 1;
} else {
up[i] = lo[i];
}
}
}
}
if (prob.status_ == 0 && paction_) {
// feasible
prob.update_model(presolvedModel_, nrows_, ncols_, nelems_);
// copy status and solution
CoinMemcpyN( prob.sol_, prob.ncols_, presolvedModel_->primalColumnSolution());
CoinMemcpyN( prob.acts_, prob.nrows_, presolvedModel_->primalRowSolution());
CoinMemcpyN( prob.colstat_, prob.ncols_, presolvedModel_->statusArray());
CoinMemcpyN( prob.rowstat_, prob.nrows_, presolvedModel_->statusArray() + prob.ncols_);
if (fixInfeasibility && hasSolution) {
// Looks feasible but double check to see if anything slipped through
int n = prob.ncols_;
double * lo = prob.clo_;
double * up = prob.cup_;
double * rsol = prob.acts_;
//memset(prob.acts_,0,prob.nrows_*sizeof(double));
presolvedModel_->matrix()->times(prob.sol_,rsol);
int i;
for (i = 0; i < n; i++) {
double gap=up[i]-lo[i];
if (rsol[i]<lo[i]-feasibilityTolerance&&fabs(rsol[i]-lo[i])<1.0e-3) {
lo[i]=rsol[i];
if (gap<1.0e5)
up[i]=lo[i]+gap;
} else if (rsol[i]>up[i]+feasibilityTolerance&&fabs(rsol[i]-up[i])<1.0e-3) {
up[i]=rsol[i];
if (gap<1.0e5)
lo[i]=up[i]-gap;
}
if (up[i] < lo[i]) {
up[i] = lo[i];
}
}
}
int n = prob.nrows_;
double * lo = prob.rlo_;
double * up = prob.rup_;
for (i = 0; i < n; i++) {
if (up[i] < lo[i]) {
if (up[i] < lo[i] - feasibilityTolerance && !fixInfeasibility) {
// infeasible
prob.status_ = 1;
} else {
up[i] = lo[i];
}
}
}
delete [] prob.sol_;
delete [] prob.acts_;
delete [] prob.colstat_;
prob.sol_ = NULL;
prob.acts_ = NULL;
prob.colstat_ = NULL;
int ncolsNow = presolvedModel_->getNumCols();
CoinMemcpyN(prob.originalColumn_, ncolsNow, originalColumn_);
#ifndef SLIM_CLP
#ifndef NO_RTTI
ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(originalModel->objectiveAsObject()));
#else
ClpQuadraticObjective * quadraticObj = NULL;
if (originalModel->objectiveAsObject()->type() == 2)
quadraticObj = (static_cast< ClpQuadraticObjective*>(originalModel->objectiveAsObject()));
#endif
if (quadraticObj) {
// set up for subset
char * mark = new char [ncols_];
memset(mark, 0, ncols_);
CoinPackedMatrix * quadratic = quadraticObj->quadraticObjective();
//const int * columnQuadratic = quadratic->getIndices();
//const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
const int * columnQuadraticLength = quadratic->getVectorLengths();
//double * quadraticElement = quadratic->getMutableElements();
int numberColumns = quadratic->getNumCols();
ClpQuadraticObjective * newObj = new ClpQuadraticObjective(*quadraticObj,
ncolsNow,
originalColumn_);
// and modify linear and check
double * linear = newObj->linearObjective();
CoinMemcpyN(presolvedModel_->objective(), ncolsNow, linear);
int iColumn;
for ( iColumn = 0; iColumn < numberColumns; iColumn++) {
if (columnQuadraticLength[iColumn])
mark[iColumn] = 1;
}
// and new
quadratic = newObj->quadraticObjective();
columnQuadraticLength = quadratic->getVectorLengths();
int numberColumns2 = quadratic->getNumCols();
for ( iColumn = 0; iColumn < numberColumns2; iColumn++) {
if (columnQuadraticLength[iColumn])
mark[originalColumn_[iColumn]] = 0;
}
presolvedModel_->setObjective(newObj);
delete newObj;
// final check
for ( iColumn = 0; iColumn < numberColumns; iColumn++)
if (mark[iColumn])
printf("Quadratic column %d modified - may be okay\n", iColumn);
delete [] mark;
}
#endif
delete [] prob.originalColumn_;
prob.originalColumn_ = NULL;
int nrowsNow = presolvedModel_->getNumRows();
CoinMemcpyN(prob.originalRow_, nrowsNow, originalRow_);
delete [] prob.originalRow_;
prob.originalRow_ = NULL;
#ifndef CLP_NO_STD
if (!dropNames && originalModel->lengthNames()) {
// Redo names
int iRow;
std::vector<std::string> rowNames;
rowNames.reserve(nrowsNow);
for (iRow = 0; iRow < nrowsNow; iRow++) {
int kRow = originalRow_[iRow];
rowNames.push_back(originalModel->rowName(kRow));
}
int iColumn;
std::vector<std::string> columnNames;
columnNames.reserve(ncolsNow);
for (iColumn = 0; iColumn < ncolsNow; iColumn++) {
int kColumn = originalColumn_[iColumn];
columnNames.push_back(originalModel->columnName(kColumn));
}
presolvedModel_->copyNames(rowNames, columnNames);
} else {
presolvedModel_->setLengthNames(0);
}
#endif
if (rowObjective_) {
int iRow;
#ifndef NDEBUG
int k = -1;
#endif
int nObj = 0;
for (iRow = 0; iRow < nrowsNow; iRow++) {
int kRow = originalRow_[iRow];
#ifndef NDEBUG
assert (kRow > k);
k = kRow;
#endif
rowObjective_[iRow] = rowObjective_[kRow];
if (rowObjective_[iRow])
nObj++;
}
if (nObj) {
printf("%d costed slacks\n", nObj);
presolvedModel_->setRowObjective(rowObjective_);
}
}
/* now clean up integer variables. This can modify original
Don't do if dupcol added columns together */
int i;
const char * information = presolvedModel_->integerInformation();
if ((prob.presolveOptions_ & 0x80000000) == 0 && information) {
int numberChanges = 0;
double * lower0 = originalModel_->columnLower();
double * upper0 = originalModel_->columnUpper();
double * lower = presolvedModel_->columnLower();
double * upper = presolvedModel_->columnUpper();
for (i = 0; i < ncolsNow; i++) {
if (!information[i])
continue;
int iOriginal = originalColumn_[i];
double lowerValue0 = lower0[iOriginal];
double upperValue0 = upper0[iOriginal];
double lowerValue = ceil(lower[i] - 1.0e-5);
double upperValue = floor(upper[i] + 1.0e-5);
lower[i] = lowerValue;
upper[i] = upperValue;
if (lowerValue > upperValue) {
numberChanges++;
presolvedModel_->messageHandler()->message(COIN_PRESOLVE_COLINFEAS,
messages)
<< iOriginal
<< lowerValue
<< upperValue
<< CoinMessageEol;
result = 1;
} else {
if (lowerValue > lowerValue0 + 1.0e-8) {
lower0[iOriginal] = lowerValue;
numberChanges++;
}
if (upperValue < upperValue0 - 1.0e-8) {
upper0[iOriginal] = upperValue;
numberChanges++;
}
}
}
if (numberChanges) {
presolvedModel_->messageHandler()->message(COIN_PRESOLVE_INTEGERMODS,
messages)
<< numberChanges
<< CoinMessageEol;
if (!result && totalPasses > 0) {
result = -1; // round again
const CoinPresolveAction *paction = paction_;
while (paction) {
const CoinPresolveAction *next = paction->next;
delete paction;
paction = next;
}
paction_ = NULL;
}
}
}
} else if (prob.status_) {
// infeasible or unbounded
result = 1;
// Put status in nelems_!
nelems_ = - prob.status_;
originalModel->setProblemStatus(prob.status_);
} else {
// no changes - model needs restoring after Lou's changes
#ifndef CLP_NO_STD
if (saveFile_ == "") {
#endif
delete presolvedModel_;
presolvedModel_ = new ClpSimplex(*originalModel);
// but we need to remove gaps
ClpPackedMatrix* clpMatrix =
dynamic_cast< ClpPackedMatrix*>(presolvedModel_->clpMatrix());
if (clpMatrix) {
clpMatrix->getPackedMatrix()->removeGaps();
}
#ifndef CLP_NO_STD
} else {
presolvedModel_ = originalModel;
}
presolvedModel_->dropNames();
#endif
// drop integer information if wanted
if (!keepIntegers)
presolvedModel_->deleteIntegerInformation();
result = 2;
}
}
if (result == 0 || result == 2) {
int nrowsAfter = presolvedModel_->getNumRows();
int ncolsAfter = presolvedModel_->getNumCols();
CoinBigIndex nelsAfter = presolvedModel_->getNumElements();
presolvedModel_->messageHandler()->message(COIN_PRESOLVE_STATS,
messages)
<< nrowsAfter << -(nrows_ - nrowsAfter)
<< ncolsAfter << -(ncols_ - ncolsAfter)
<< nelsAfter << -(nelems_ - nelsAfter)
<< CoinMessageEol;
} else {
destroyPresolve();
if (presolvedModel_ != originalModel_)
delete presolvedModel_;
presolvedModel_ = NULL;
}
return presolvedModel_;
}