limp-cbc-0.3.2.0: cbits/coin/CglLandP.cpp
// Copyright (C) 2005-2009, Pierre Bonami and others. All Rights Reserved.
// Author: Pierre Bonami
// Tepper School of Business
// Carnegie Mellon University, Pittsburgh, PA 15213
// Date: 07/21/05
//
// $Id: CglLandP.cpp 1123 2013-04-06 20:47:24Z stefan $
//
// This code is licensed under the terms of the Eclipse Public License (EPL).
//---------------------------------------------------------------------------
#include "CglLandP.hpp"
#include "CglLandPSimplex.hpp"
#define INT_INFEAS(value) fabs(value - floor(value+0.5))
#include "CglConfig.h"
#ifdef COIN_HAS_OSICLP
#include "OsiClpSolverInterface.hpp"
#endif
#define CLONE_SI //Solver is cloned between two cuts
#include "CoinTime.hpp"
#include "CglGomory.hpp"
#include "CoinFactorization.hpp"
#include <fstream>
namespace LAP
{
//Setup output messages
LapMessages::LapMessages( )
:CoinMessages(LAP_MESSAGES_DUMMY_END)
{
strcpy(source_,"Lap");
addMessage(BEGIN_ROUND,CoinOneMessage( 1, 2,"Starting %s round %d variable considered for separation."));
addMessage(END_ROUND,CoinOneMessage(2, 2,"End ouf %s round %d cut generated in %g seconds."));
addMessage(DURING_SEP,CoinOneMessage(3,1,"After %g seconds, separated %d cuts."));
addMessage(CUT_REJECTED, CoinOneMessage(4,1,"Cut rejected for %s."));
addMessage(CUT_FAILED,CoinOneMessage(5,1,"Generation failed."));
addMessage(CUT_GAP, CoinOneMessage(7,1,"CUTGAP after %i pass objective is %g"));
addMessage(LAP_CUT_FAILED_DO_MIG, CoinOneMessage(3006,1,"Failed to generate a cut generate a Gomory cut instead"));
}
}
using namespace LAP;
CglLandP::Parameters::Parameters():
CglParam(),
pivotLimit(20),
pivotLimitInTree(10),
maxCutPerRound(5000),
failedPivotLimit(1),
degeneratePivotLimit(0),
extraCutsLimit(5),
pivotTol(1e-4),
away(5e-4),
timeLimit(COIN_DBL_MAX),
singleCutTimeLimit(COIN_DBL_MAX),
rhsWeight(1.),
useTableauRow(true),
modularize(false),
strengthen(true),
countMistakenRc(false),
sepSpace(Fractional),
perturb(true),
normalization(Unweighted),
rhsWeightType(Fixed),
lhs_norm(L1),
generateExtraCuts(none),
pivotSelection(mostNegativeRc)
{
EPS = 1e-08;
}
CglLandP::Parameters::Parameters(const Parameters &other):
CglParam(other),
pivotLimit(other.pivotLimit),
pivotLimitInTree(other.pivotLimitInTree),
maxCutPerRound(other.maxCutPerRound),
failedPivotLimit(other.failedPivotLimit),
degeneratePivotLimit(other.degeneratePivotLimit),
extraCutsLimit(other.extraCutsLimit),
pivotTol(other.pivotTol),
away(other.away),
timeLimit(other.timeLimit),
singleCutTimeLimit(other.singleCutTimeLimit),
rhsWeight(other.rhsWeight),
useTableauRow(other.useTableauRow),
modularize(other.modularize),
strengthen(other.strengthen),
countMistakenRc(other.countMistakenRc),
sepSpace(other.sepSpace),
perturb(other.perturb),
normalization(other.normalization),
rhsWeightType(other.rhsWeightType),
lhs_norm(other.lhs_norm),
generateExtraCuts(other.generateExtraCuts),
pivotSelection(other.pivotSelection)
{}
CglLandP::Parameters & CglLandP::Parameters::operator=(const Parameters &other)
{
if (this != &other)
{
CglParam::operator=(other);
pivotLimit = other.pivotLimit;
pivotLimitInTree = other.pivotLimitInTree;
maxCutPerRound = other.maxCutPerRound;
failedPivotLimit = other.failedPivotLimit;
degeneratePivotLimit = other.failedPivotLimit;
extraCutsLimit = other.extraCutsLimit;
pivotTol = other.pivotTol;
away = other.away;
timeLimit = other.timeLimit;
singleCutTimeLimit = other.singleCutTimeLimit;
rhsWeight = other.rhsWeight;
useTableauRow = other.useTableauRow;
modularize = other.modularize;
strengthen = other.strengthen;
countMistakenRc = other.countMistakenRc;
sepSpace = other.sepSpace;
perturb = other.perturb;
normalization = other.normalization;
rhsWeightType = other.rhsWeightType;
lhs_norm = other.lhs_norm;
generateExtraCuts = other.generateExtraCuts;
pivotSelection = other.pivotSelection;
}
return *this;
}
CglLandP::CachedData::CachedData(int nBasics, int nNonBasics):
basics_(NULL), nonBasics_(NULL), nBasics_(nBasics),
nNonBasics_(nNonBasics), basis_(NULL), colsol_(NULL),
slacks_(NULL), integers_(NULL), solver_(NULL)
{
if (nBasics_>0)
{
basics_ = new int[nBasics_];
integers_ = new bool [nNonBasics_ + nBasics_];
}
if (nNonBasics_>0)
nonBasics_ = new int[nNonBasics_];
if (nBasics_ + nNonBasics_ > 0)
{
colsol_ = new double[nBasics_ + nNonBasics_];
slacks_ = &colsol_[nNonBasics_];
}
}
CglLandP::CachedData::CachedData(const CachedData &source):
basics_(NULL), nonBasics_(NULL), nBasics_(source.nBasics_),
nNonBasics_(source.nNonBasics_), basis_(NULL),
colsol_(NULL), slacks_(NULL), integers_(NULL), solver_(NULL)
{
if (nBasics_>0)
{
basics_ = new int[nBasics_];
CoinCopyN(source.basics_, nBasics_, basics_);
integers_ = new bool [nNonBasics_ + nBasics_];
CoinCopyN(source.integers_, nBasics_ + nNonBasics_, integers_);
}
if (nNonBasics_>0)
{
nonBasics_ = new int[nNonBasics_];
CoinCopyN(source.nonBasics_, nBasics_, nonBasics_);
}
if (nBasics_ + nNonBasics_ > 0)
{
colsol_ = new double[nBasics_ + nNonBasics_];
slacks_ = &colsol_[nNonBasics_];
CoinCopyN(source.colsol_, nBasics_ + nNonBasics_, colsol_);
}
if (source.basis_!=NULL)
basis_ = new CoinWarmStartBasis(*source.basis_);
if (source.solver_!=NULL)
solver_ = source.solver_->clone();
}
CglLandP::CachedData& CglLandP::CachedData::operator=(const CachedData &source)
{
if (this != &source)
{
nBasics_ = source.nBasics_;
nNonBasics_ = source.nNonBasics_;
if (basics_ == NULL) delete [] basics_;
basics_ = NULL;
if (nonBasics_ == NULL) delete [] nonBasics_;
nonBasics_ = NULL;
if (basis_ == NULL) delete [] basis_;
basis_ = NULL;
if (colsol_ == NULL) delete [] colsol_;
colsol_ = NULL;
if (slacks_ == NULL) delete [] slacks_;
slacks_ = NULL;
if (integers_ == NULL) delete [] integers_;
integers_ = NULL;
if (nBasics_>0)
{
basics_ = new int[nBasics_];
CoinCopyN(source.basics_, nBasics_, basics_);
integers_ = new bool [nBasics_ + nNonBasics_];
CoinCopyN(source.integers_, nBasics_ + nNonBasics_, integers_);
}
if (nNonBasics_>0)
{
nonBasics_ = new int[nNonBasics_];
CoinCopyN(source.nonBasics_, nBasics_, nonBasics_);
}
if (nBasics_ + nNonBasics_ > 0)
{
colsol_ = new double[nBasics_ + nNonBasics_];
slacks_ = &colsol_[nNonBasics_];
CoinCopyN(source.colsol_, nBasics_ + nNonBasics_, colsol_);
}
if (source.basis_!=NULL)
basis_ = new CoinWarmStartBasis(*source.basis_);
delete solver_;
if (source.solver_)
solver_ = source.solver_->clone();
}
return *this;
}
void
CglLandP::CachedData::getData(const OsiSolverInterface &si)
{
int nBasics = si.getNumRows();
int nNonBasics = si.getNumCols();
if (basis_ != NULL)
delete basis_;
basis_ = dynamic_cast<CoinWarmStartBasis *> (si.getWarmStart());
if (!basis_)
throw NoBasisError();
if (nBasics_ > 0 || nBasics != nBasics_)
{
delete [] basics_;
basics_ = NULL;
}
if (basics_ == NULL)
{
basics_ = new int[nBasics];
nBasics_ = nBasics;
}
if (nNonBasics_ > 0 || nNonBasics != nNonBasics_)
{
delete [] nonBasics_;
nonBasics_ = NULL;
}
if (nonBasics_ == NULL)
{
nonBasics_ = new int[nNonBasics];
nNonBasics_ = nNonBasics;
}
int n = nBasics + nNonBasics;
if ( nBasics_ + nNonBasics_ > 0 || nBasics_ + nNonBasics_ != n)
{
delete [] colsol_;
delete [] integers_;
integers_ = NULL;
colsol_ = NULL;
slacks_ = NULL;
}
if (colsol_ == NULL)
{
colsol_ = new double[n];
slacks_ = &colsol_[nNonBasics];
}
if (integers_ == NULL)
{
integers_ = new bool[n];
}
const double * rowLower = si.getRowLower();
const double * rowUpper = si.getRowUpper();
//determine which slacks are integer
const CoinPackedMatrix * m = si.getMatrixByCol();
const double * elems = m->getElements();
const int * inds = m->getIndices();
const CoinBigIndex * starts = m->getVectorStarts();
const int * lengths = m->getVectorLengths();
// int numElems = m->getNumElements();
int numCols = m->getNumCols();
assert(numCols == nNonBasics_);
// int numRows = m->getNumRows();
CoinFillN(integers_ ,n, true);
for (int i = 0 ; i < numCols ; i++)
{
if (si.isContinuous(i))
integers_[i] = false;
}
bool * integerSlacks = integers_ + numCols;
for (int i = 0 ; i < nBasics ; i++)
{
if (rowLower[i] > -1e50 && INT_INFEAS(rowLower[i]) > 1e-15)
integerSlacks[i] = false;
if (rowUpper[i] < 1e50 && INT_INFEAS(rowUpper[i]) > 1e-15)
integerSlacks[i] = false;
}
for (int i = 0 ; i < numCols ; i++)
{
CoinBigIndex end = starts[i] + lengths[i];
if (integers_[i])
{
for (CoinBigIndex k=starts[i] ; k < end; k++)
{
if (integerSlacks[inds[k]] && INT_INFEAS(elems[k])>1e-15 )
integerSlacks[inds[k]] = false;
}
}
else
{
for (CoinBigIndex k=starts[i] ; k < end; k++)
{
if (integerSlacks[inds[k]])
integerSlacks[inds[k]] = false;
}
}
}
CoinCopyN(si.getColSolution(), si.getNumCols(), colsol_);
CoinCopyN(si.getRowActivity(), si.getNumRows(), slacks_);
for (int i = 0 ; i < si.getNumRows() ; i++)
{
slacks_[i]*=-1;
if (rowLower[i]>-1e50)
{
slacks_[i] += rowLower[i];
}
else
{
slacks_[i] += rowUpper[i];
}
}
//Now get the fill the arrays;
nNonBasics = 0;
nBasics = 0;
//For having the index variables correctly ordered we need to access to OsiSimplexInterface
{
OsiSolverInterface * ncSi = (const_cast<OsiSolverInterface *>(&si));
ncSi->enableSimplexInterface(0);
ncSi->getBasics(basics_);
// Save enabled solver
solver_ = si.clone();
#ifdef COIN_HAS_OSICLP
OsiClpSolverInterface * clpSi = dynamic_cast<OsiClpSolverInterface *>(solver_);
const OsiClpSolverInterface * clpSiRhs = dynamic_cast<const OsiClpSolverInterface *>(&si);
if (clpSi)
clpSi->getModelPtr()->copyEnabledStuff(clpSiRhs->getModelPtr());;
#endif
ncSi->disableSimplexInterface();
}
int numStructural = basis_->getNumStructural();
for (int i = 0 ; i < numStructural ; i++)
{
if (basis_->getStructStatus(i)== CoinWarmStartBasis::basic)
{
nBasics++;
//Basically do nothing
}
else
{
nonBasics_[nNonBasics++] = i;
}
}
int numArtificial = basis_->getNumArtificial();
for (int i = 0 ; i < numArtificial ; i++)
{
if (basis_->getArtifStatus(i)== CoinWarmStartBasis::basic)
{
//Just check number of basics
nBasics++;
}
else
{
nonBasics_[nNonBasics++] = i + basis_->getNumStructural();
}
}
}
void
CglLandP::CachedData::clean(){
if (basics_!=NULL)
delete [] basics_;
basics_ = NULL;
if (nonBasics_!=NULL)
delete [] nonBasics_;
nonBasics_ = NULL;
if (colsol_ != NULL)
delete [] colsol_;
colsol_ = NULL;
delete basis_;
basis_ = NULL;
if (integers_)
delete [] integers_;
integers_ = NULL;
nBasics_ = 0;
nNonBasics_ = 0;
delete solver_;
solver_ = NULL;
}
CglLandP::CachedData::~CachedData()
{
if (basics_!=NULL)
delete [] basics_;
if (nonBasics_!=NULL)
delete [] nonBasics_;
if (colsol_ != NULL)
delete [] colsol_;
delete basis_;
if (integers_)
delete [] integers_;
delete solver_;
}
CglLandP::CglLandP(const CglLandP::Parameters ¶ms,
const LAP::Validator &validator):
params_(params), cached_(), validator_(validator), numcols_(-1),
originalColLower_(NULL), originalColUpper_(NULL),
canLift_(false),
extraCuts_()
{
handler_ = new CoinMessageHandler();
handler_->setLogLevel(0);
messages_ = LapMessages();
}
CglLandP::~CglLandP()
{
delete handler_;
if (originalColLower_ != NULL)
delete [] originalColLower_;
if (originalColUpper_ != NULL)
delete [] originalColUpper_;
}
CglLandP::CglLandP(const CglLandP & source):
CglCutGenerator(source),
params_(source.params_), cached_(source.cached_),
validator_(source.validator_), numcols_(source.numcols_),
originalColLower_(NULL), originalColUpper_(NULL),
canLift_(source.canLift_),
extraCuts_(source.extraCuts_)
{
handler_ = new CoinMessageHandler();
handler_->setLogLevel(source.handler_->logLevel());
messages_ = LapMessages();
if (numcols_ != -1)
{
assert(numcols_ > 0);
assert(originalColLower_!=NULL);
assert(originalColUpper_!=NULL);
originalColLower_ = new double[numcols_];
originalColUpper_ = new double[numcols_];
CoinCopyN(source.originalColLower_,numcols_,originalColLower_);
CoinCopyN(source.originalColUpper_,numcols_,originalColUpper_);
}
}
/** Assignment operator */
CglLandP& CglLandP::operator=(const CglLandP &rhs)
{
if (this != &rhs)
{
params_ = rhs.params_;
cached_ = rhs.cached_;
validator_ = rhs.validator_;
extraCuts_ = rhs.extraCuts_;
}
return *this;
}
CglCutGenerator *
CglLandP::clone() const
{
return new CglLandP(*this);
}
extern double restaurationTime;
struct cutsCos
{
int i;
int j;
double angle;
cutsCos(int i_, int j_ , double angle_):i(i_), j(j_), angle(angle_)
{
}
bool operator<(const cutsCos&other)const
{
return angle > other.angle;
}
};
void
CglLandP::scanExtraCuts(OsiCuts& cs, const double * colsol) const
{
int numAdded = 0;
for (int i = extraCuts_.sizeRowCuts() - 1; i > -1 ; i--)
{
double violation = extraCuts_.rowCut(i).violated(colsol);
if (violation > 0.)
{
cs.insert(extraCuts_.rowCut(i));
numAdded++;
// std::cout<<"A cut computed in a previous iteration is violated by "<<violation<<"."<<std::endl;
//extraCuts_.eraseRowCut(i);
}
}
// std::cout<<"Added "<<numAdded<<" previously generated cuts."<<std::endl;
}
void
CglLandP::generateCuts(const OsiSolverInterface & si, OsiCuts & cs,
const CglTreeInfo info )
{
if ((info.pass == 0) && !info.inTree)
{
numrows_ = si.getNumRows();
}
// scanExtraCuts(cs, si.getColSolution());
Parameters params = params_;
params.rhsWeight = numrows_ + 2;
handler_->message(CUT_GAP, messages_)<<info.pass<<si.getObjValue() <<CoinMessageEol;
if (info.inTree) //put lower pivot limit
{
params.pivotLimit = std::min(params.pivotLimit, params.pivotLimitInTree);
params.countMistakenRc = true;
}
if (params.timeLimit < 0)
{
params.pivotLimit = 0;
}
assert(si.basisIsAvailable());
#ifdef APPEND_ROW
OsiSolverInterface * t_si = si.clone();
if (params.modularize)
{
int new_idx = si.getNumCols();
int v_idx[1] = {new_idx};
double v_val[1] = {-1};
CoinPackedVector v(1, v_idx, v_val, false);
t_si->addCol(CoinPackedVector(), 0, 1, 0);
t_si->setInteger(new_idx);
t_si->addRow(v,0, 0);
t_si->resolve();
}
#else
const OsiSolverInterface * t_si = &si;
#endif
cached_.getData(*t_si);
CglLandPSimplex landpSi(*t_si, cached_, params, validator_);
if (params.generateExtraCuts == CglLandP::AllViolatedMigs)
{
landpSi.genThisBasisMigs(cached_, params);
}
landpSi.setLogLevel(handler_->logLevel());
int nCut = 0;
std::vector<int> indices;
getSortedFractionalIndices(indices,cached_, params);
#ifndef NDEBUG
int numrows = si.getNumRows();
#endif
#ifdef DO_STAT
//Get informations on current optimum
{
OsiSolverInterface * gapTester = si.clone();
gapTester->resolve();
roundsStats_.analyseOptimalBasis(gapTester,info.pass, numrows_);
delete gapTester;
}
#endif
params_.timeLimit += CoinCpuTime();
CoinRelFltEq eq(1e-04);
for (unsigned int i = 0; i < indices.size() && nCut < params.maxCutPerRound &&
nCut < cached_.nBasics_ ; i++)
{
//Check for time limit
int iRow = indices[i];
assert(iRow < numrows);
OsiRowCut cut;
int code=1;
OsiSolverInterface * ncSi = NULL;
if (params.pivotLimit != 0)
{
ncSi = t_si->clone();
landpSi.setSi(ncSi);
ncSi->setDblParam(OsiDualObjectiveLimit, COIN_DBL_MAX);
ncSi->messageHandler()->setLogLevel(0);
}
int generated = 0;
if (params.pivotLimit == 0)
{
generated = landpSi.generateMig(iRow, cut, params);
}
else
{
generated = landpSi.optimize(iRow, cut, cached_, params);
if (params.generateExtraCuts == CglLandP::AllViolatedMigs)
{
landpSi.genThisBasisMigs(cached_, params);
}
landpSi.resetSolver(cached_.basis_);
}
code = 0;
if (generated)
code = validator_(cut, cached_.colsol_, si, params, originalColLower_, originalColUpper_);
if (!generated || code)
{
if (params.pivotLimit !=0)
{
handler_->message(LAP_CUT_FAILED_DO_MIG, messages_)<<validator_.failureString(code)<<CoinMessageEol;
landpSi.freeSi();
OsiSolverInterface * ncSi = t_si->clone();
landpSi.setSi(ncSi);
params.pivotLimit = 0;
if (landpSi.optimize(iRow, cut, cached_, params))
{
code = validator_(cut, cached_.colsol_, si, params, originalColLower_, originalColUpper_);
}
params.pivotLimit = params_.pivotLimit;
}
}
if (params.pivotLimit != 0)
{
landpSi.freeSi();
}
if (code)
{
handler_->message(CUT_REJECTED, messages_)<<
validator_.failureString(code)<<CoinMessageEol;
}
else
{
if (canLift_)
{
cut.setGloballyValid(true);
}
cs.insertIfNotDuplicate(cut, eq);
//cs.insert(cut);
{
//std::cout<<"Violation "<<cut.violated(cached_.colsol_)<<std::endl;
nCut++;
}
}
}
Cuts& extra = landpSi.extraCuts();
for (int i = 0 ; i < cached_.nNonBasics_; i++)
{
OsiRowCut * cut = extra.rowCut(i);
if (cut == NULL) continue;
int code = validator_(*cut, cached_.colsol_, si, params,
originalColLower_, originalColUpper_);
if (code)
{
handler_->message(LAP_CUT_FAILED_DO_MIG, messages_)
<<validator_.failureString(code)<<CoinMessageEol;
}
else
{
cs.insertIfNotDuplicate(*cut, eq);
{
nCut++;
}
}
delete cut;
}
landpSi.outPivInfo(nCut);
params_.timeLimit -= CoinCpuTime();
cached_.clean();
#ifdef APPEND_ROW
assert(t_si != &si);
delete t_si;
#endif
}
template < class S, class T, class U >
class StableCompare
{
public:
inline bool operator()(const CoinTriple<S,T,U>& t1,
const CoinTriple<S,T,U>& t2) const
{
return (t1.third < t2.third) ||
((t1.third == t2.third) && (t1.second < t2.second));
}
};
template <class T1,class T2>
struct StableExternalComp
{
const std::vector<T1> &vec_1_;
const std::vector<T2> &vec_2_;
StableExternalComp(const std::vector<T1> &vec_1,
const std::vector<T2> &vec_2):
vec_1_(vec_1),
vec_2_(vec_2)
{
}
CoinRelFltEq eq;
bool operator()(int i, int j)
{
bool result = (vec_1_[i] < vec_1_[j]) ||
( ((vec_1_[i]== vec_1_[j]))
&& (vec_2_[i] < vec_2_[j]));
return result;
}
};
void
CglLandP::getSortedFractionalIndices(std::vector<int> &frac_indices,
const CachedData &data,
const CglLandP::Parameters & params) const
{
std::vector<int> colIndices;
std::vector<double> values;
std::vector<int> indices;
for (int i = 0 ; i < data.nBasics_ ; i++)
{
const int& iCol = data.basics_[i];
if (iCol >= data.nNonBasics_ ||
!data.integers_[iCol] ||
INT_INFEAS(data.colsol_[iCol]) <= params.away)
continue;
const double value = INT_INFEAS(data.colsol_[iCol]);
frac_indices.push_back(i);
indices.push_back(static_cast<int>(values.size()));
values.push_back(- value);
colIndices.push_back(iCol);
}
std::sort(indices.begin(), indices.end(),StableExternalComp<double, int>(values,colIndices));
colIndices = frac_indices;
for (unsigned int i = 0; i < indices.size() ; i++)
{
frac_indices[i] = colIndices[indices[i]];
}
}