limp-cbc-0.3.2.0: cbits/coin/CglTwomir.cpp
// $Id: CglTwomir.cpp 1132 2013-04-25 18:57:12Z 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).
#include <cstdlib>
#include <cstdio>
#include <cmath>
#include <cfloat>
#include <cassert>
#include <iostream>
#include "CoinPragma.hpp"
#include "CoinHelperFunctions.hpp"
#include "CoinPackedVector.hpp"
#include "CoinPackedMatrix.hpp"
#include "CoinIndexedVector.hpp"
#include "OsiRowCutDebugger.hpp"
#include "CoinFactorization.hpp"
#include "CoinWarmStartBasis.hpp"
#include "CoinFinite.hpp"
#include "CoinPackedMatrix.hpp"
#include "OsiRowCutDebugger.hpp"
#include "CoinWarmStartBasis.hpp"
#include "CglTwomir.hpp"
class CoinWarmStartBasis;
#define COIN_HAS_CLP_TWOMIR
#ifdef COIN_HAS_CLP_TWOMIR
#include "OsiClpSolverInterface.hpp"
#endif
#undef DGG_DEBUG_DGG
//#define DGG_DEBUG_DGG 1
//#define CGL_DEBUG
//#define CGL_DEBUG_ZERO
#define q_max data->cparams.q_max
#define q_min data->cparams.q_min
#define t_max data->cparams.t_max
#define t_min data->cparams.t_min
#define a_max data->cparams.a_max
#define max_elements data->cparams.max_elements
#ifdef CGL_DEBUG
// Declarations and defines for debug build.
#define talk true
namespace {
const OsiSolverInterface *six ;
}
void write_cut( DGG_constraint_t *cut){ //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
printf("2mir_test: cut: !!!!!!!!!!!!!!!!!!!!!!!***********************************\n");
for (int i=0; i<cut->nz; i++)
{ printf(" %12.10f x[%d] ", cut->coeff[i],cut->index[i]);}
printf(" >= %12.10f \n", cut->rhs);
}
void testus( DGG_constraint_t *cut){ //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
// const OsiRowCutDebugger debugg(*six, "flugpl.mps");
const OsiRowCutDebugger debugg(*six, "egout.mps");
const OsiRowCutDebugger *debugger = &debugg;
if (debugger&&!debugger->onOptimalPath(*six))
return;
OsiRowCut rowcut;
rowcut.setRow(cut->nz, cut->index, cut->coeff);
rowcut.setUb(DBL_MAX);
rowcut.setLb(cut->rhs);
if(debugger->invalidCut(rowcut)){
write_cut(cut);
//assert(0);
}
}
#else // CGL_DEBUG
#define talk false
#endif // CGL_DEBUG
//-------------------------------------------------------------------
// Generate cuts
//-------------------------------------------------------------------
void CglTwomir::generateCuts(const OsiSolverInterface & si, OsiCuts & cs,
const CglTreeInfo info )
{
# ifdef CGL_DEBUG
//!!!!!!!!!!!!!!!!!!
six = &si;
# endif
const double * colUpper = si.getColUpper();
const double * colLower = si.getColLower();
const OsiSolverInterface * useSolver;
#ifdef COIN_HAS_CLP_TWOMIR
double * objective = NULL;
OsiClpSolverInterface * clpSolver = dynamic_cast<OsiClpSolverInterface *>(originalSolver_);
int numberOriginalRows;
int numberColumns=si.getNumCols();
int twomirType=0;
if (!clpSolver) {
#endif
useSolver=&si;
// Temp - check if free variables
int ncol = useSolver->getNumCols();
int numberFree=0;
for (int i=0;i<ncol;i++) {
if (colLower[i]<-1.0e20&&colUpper[i]>1.0e20)
numberFree++;
}
if (numberFree) {
#ifdef COIN_DEVELOP
if (!info.inTree&&!info.pass)
printf("CglTwoMir - %d free variables - returning\n",numberFree);
#endif
return;
}
#ifdef COIN_HAS_CLP_TWOMIR
} else {
useSolver = originalSolver_;
assert (twomirType_);
// check simplex is plausible
if (!clpSolver->getNumRows()||numberColumns!=clpSolver->getNumCols()) {
delete originalSolver_;
originalSolver_=si.clone();
clpSolver = dynamic_cast<OsiClpSolverInterface *>(originalSolver_);
assert (clpSolver);
useSolver = originalSolver_;
}
ClpSimplex * simplex = clpSolver->getModelPtr();
numberOriginalRows = simplex->numberRows();
int numberRows = si.getNumRows();
assert (numberOriginalRows<=numberRows);
// only do if different (unless type 2x)
twomirType = twomirType_%10;
int whenToDo = twomirType_/10;
if (whenToDo==2 ||(numberRows>numberOriginalRows && whenToDo==1
&& (info.options&512)==0) ||
((info.options&1024)!=0 && (info.options&512)==0)) {
// bounds
const double * solution = si.getColSolution();
memcpy(simplex->columnLower(),colLower,numberColumns*sizeof(double));
memcpy(simplex->columnUpper(),colUpper,numberColumns*sizeof(double));
for (int i=0;i<numberColumns;i++) {
if (colLower[i]<-1.0e20&&colUpper[i]>1.0e20) {
double lower=-COIN_DBL_MAX;
double upper=COIN_DBL_MAX;
if (solution[i]>0.0)
lower=-1.0e10;
else
upper=1.0e10;
originalSolver_->setColLower(i,lower);
originalSolver_->setColUpper(i,upper);
}
}
double * obj = simplex->objective();
objective = CoinCopyOfArray(obj,numberColumns);
const double * pi = si.getRowPrice();
const CoinPackedMatrix * rowCopy = si.getMatrixByRow();
const int * column = rowCopy->getIndices();
const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
const int * rowLength = rowCopy->getVectorLengths();
const double * rowElements = rowCopy->getElements();
const double * rowLower = si.getRowLower();
const double * rowUpper = si.getRowUpper();
int numberCopy;
int numberAdd;
double * rowLower2 = NULL;
double * rowUpper2 = NULL;
int * column2 = NULL;
CoinBigIndex * rowStart2 = NULL;
double * rowElements2 = NULL;
char * copy = new char [numberRows-numberOriginalRows];
memset(copy,0,numberRows-numberOriginalRows);
if (twomirType==2) {
numberCopy=0;
numberAdd=0;
for (int iRow=numberOriginalRows;iRow<numberRows;iRow++) {
bool simple = true;
for (int k=rowStart[iRow];
k<rowStart[iRow]+rowLength[iRow];k++) {
double value = rowElements[k];
if (value!=floor(value+0.5)) {
simple=false;
break;
}
}
if (simple) {
numberCopy++;
numberAdd+=rowLength[iRow];
copy[iRow-numberOriginalRows]=1;
}
}
if (numberCopy) {
//printf("Using %d rows out of %d\n",numberCopy,numberRows-numberOriginalRows);
rowLower2 = new double [numberCopy];
rowUpper2 = new double [numberCopy];
rowStart2 = new CoinBigIndex [numberCopy+1];
rowStart2[0]=0;
column2 = new int [numberAdd];
rowElements2 = new double [numberAdd];
}
}
numberCopy=0;
numberAdd=0;
const double * rowSolution = si.getRowActivity();
double offset=0.0;
for (int iRow=numberOriginalRows;iRow<numberRows;iRow++) {
if (!copy[iRow-numberOriginalRows]) {
double value = pi[iRow];
offset += rowSolution[iRow]*value;
for (int k=rowStart[iRow];
k<rowStart[iRow]+rowLength[iRow];k++) {
int iColumn=column[k];
obj[iColumn] -= value*rowElements[k];
}
} else {
rowLower2[numberCopy]=rowLower[iRow];
rowUpper2[numberCopy]=rowUpper[iRow];
for (int k=rowStart[iRow];
k<rowStart[iRow]+rowLength[iRow];k++) {
column2[numberAdd]=column[k];
rowElements2[numberAdd++]=rowElements[k];
}
numberCopy++;
rowStart2[numberCopy]=numberAdd;
}
}
#if 0
CoinThreadRandom randomNumberGenerator;
const double * solution = si.getColSolution();
for (int i=0;i<numberColumns;i++) {
if (intVar[i]==1) {
double randomNumber = randomNumberGenerator.randomDouble();
//randomNumber = 0.001*floor(randomNumber*1000.0);
if (solution[i]>0.5)
obj[i] -= randomNumber*0.001*fabs(obj[i]);
else
obj[i] += randomNumber*0.001*fabs(obj[i]);
}
}
#endif
if (numberCopy) {
clpSolver->addRows(numberCopy,
rowStart2,column2,rowElements2,
rowLower2,rowUpper2);
delete [] rowLower2 ;
delete [] rowUpper2 ;
delete [] column2 ;
delete [] rowStart2 ;
delete [] rowElements2 ;
}
delete [] copy;
memcpy(simplex->primalColumnSolution(),si.getColSolution(),
numberColumns*sizeof(double));
CoinWarmStart * warmstart = si.getWarmStart();
CoinWarmStartBasis* warm =
dynamic_cast<CoinWarmStartBasis*>(warmstart);
warm->resize(simplex->numberRows(),numberColumns);
clpSolver->setBasis(*warm);
delete warm;
simplex->setDualObjectiveLimit(COIN_DBL_MAX);
simplex->setLogLevel(0);
clpSolver->resolve();
//printf("Trying - %d its status %d objs %g %g - with offset %g\n",
// simplex->numberIterations(),simplex->status(),
// simplex->objectiveValue(),si.getObjValue(),simplex->objectiveValue()+offset);
//simplex->setLogLevel(0);
if (simplex->status()) {
//printf("BAD status %d\n",simplex->status());
//clpSolver->writeMps("clp");
//si.writeMps("si");
delete [] objective;
objective=NULL;
useSolver=&si;
}
}
}
#endif
useSolver->getStrParam(OsiProbName,probname_) ;
int numberRowCutsBefore = cs.sizeRowCuts();
DGG_list_t cut_list;
DGG_list_init (&cut_list);
DGG_data_t* data = DGG_getData(reinterpret_cast<const void *> (useSolver));
// Note that the lhs variables are hash defines to data->cparams.*
q_max = q_max_;
q_min = q_min_;
t_max = t_max_;
t_min = t_min_;
a_max = a_max_;
max_elements = info.inTree ? max_elements_ : max_elements_root_;
data->gomory_threshold = info.inTree ? away_ : awayAtRoot_;
if (!info.inTree) {
//const CoinPackedMatrix * columnCopy = useSolver->getMatrixByCol();
//int numberColumns=columnCopy->getNumCols();
if (!info.pass||(info.options&32)!=0) {
max_elements=useSolver->getNumCols();
//} else {
//int numberRows=columnCopy.getNumRows();
//int numberElements=columnCopy->getNumElements();
//if (max_elements>500&&numberElements>10*numberColumns)
//max_elements=numberColumns;
}
}
if (!do_mir_) t_max = t_min - 1;
if (!do_2mir_) q_max = q_min - 1;
if (do_tab_ && info.level < 1 && info.pass < 6)
DGG_generateTabRowCuts( &cut_list, data, reinterpret_cast<const void *> (useSolver) );
if (do_form_)
DGG_generateFormulationCuts( &cut_list, data, reinterpret_cast<const void *> (useSolver),
info.formulation_rows,
randomNumberGenerator_);
#ifdef CGL_DEBUG
const OsiRowCutDebugger debugg(si,probname_.c_str()) ;
const OsiRowCutDebugger *debugger = &debugg;
if (debugger&&!debugger->onOptimalPath(si))
debugger = NULL;
else
{if(talk) printf ("2mir_test: debug success\n");}
#endif
int i;
for ( i=0; i<cut_list.n; i++){
DGG_constraint_t *cut = cut_list.c[i];
OsiRowCut rowcut;
if (cut->nz<max_elements) {
// See if any zero coefficients!!!!!!!
int nZero=0;
for( int i=0; i < cut->nz; i++) {
if (!cut->coeff[i])
nZero++;
}
#ifdef CGL_DEBUG_ZERO
if (nZero) {
printf("Cut ");
for( int i=0; i < cut->nz; i++) {
printf("%d %g ",cut->index[i],cut->coeff[i]);
}
printf("\n");
}
#endif
if (nZero) {
#ifdef CGL_DEBUG_ZERO
printf("TwoMir cut had %d zero coefficients!\n",nZero);
#endif
} else {
#ifdef CBC_CHECK_CUT
double rhs = cut->rhs;
int * cutIndex = cut->index;
double * packed = cut->coeff;
int i,number2=cut->nz;
int number=0;
double largest=0.0;
double smallest=1.0e30;
const double *colUpper = useSolver->getColUpper();
const double *colLower = useSolver->getColLower();
bool goodCut=true;
for (i=0;i<number2;i++) {
double value=fabs(packed[i]);
if (value<1.0e-9) {
int iColumn = cutIndex[i];
if (colUpper[iColumn]-colLower[iColumn]<100.0) {
// weaken cut
if (packed[i]>0.0)
rhs -= value*colUpper[iColumn];
else
rhs += value*colLower[iColumn];
} else {
// throw away
goodCut=false;
break;
}
} else {
int iColumn = cutIndex[i];
if (colUpper[iColumn]!=colLower[iColumn]) {
largest=CoinMax(largest,value);
smallest=CoinMin(smallest,value);
cutIndex[number]=cutIndex[i];
packed[number++]=packed[i];
} else {
// fixed so subtract out
rhs -= packed[i]*colLower[iColumn];
}
}
}
if (largest<5.0e9*smallest&&goodCut) {
rowcut.setRow(number, cutIndex, packed);
rowcut.setUb(COIN_DBL_MAX);
rowcut.setLb(rhs);
cs.insert(rowcut);
}
#else
rowcut.setRow(cut->nz, cut->index, cut->coeff);
rowcut.setUb(DBL_MAX);
rowcut.setLb(cut->rhs);
cs.insert(rowcut);
#endif
}
#ifdef CGL_DEBUG
if (debugger) {
if (debugger->invalidCut(rowcut)) {
write_cut(cut);
printf ("2mir_test: debug failed, mayday, mayday **********************************\n");}
//assert(0);
}
//assert(!debugger->invalidCut(rowcut));
#endif
}
}
for ( i=0; i<cut_list.n; i++)
DGG_freeConstraint (cut_list.c[i]);
DGG_list_free (&cut_list);
DGG_freeData (data);
if (!info.inTree&&((info.options&4)==4||((info.options&8)&&!info.pass))) {
int numberRowCutsAfter = cs.sizeRowCuts();
for (int i=numberRowCutsBefore;i<numberRowCutsAfter;i++) {
int length = cs.rowCutPtr(i)->row().getNumElements();
if (length<=max_elements_)
cs.rowCutPtr(i)->setGloballyValid();
}
}
#ifdef COIN_HAS_CLP_TWOMIR
if (objective) {
int numberRowCutsAfter = cs.sizeRowCuts();
ClpSimplex * simplex = clpSolver->getModelPtr();
memcpy(simplex->objective(),objective,numberColumns*sizeof(double));
delete [] objective;
// take out locally useless cuts
const double * solution = si.getColSolution();
double primalTolerance = 1.0e-7;
for (int k = numberRowCutsAfter - 1; k >= numberRowCutsBefore; k--) {
const OsiRowCut * thisCut = cs.rowCutPtr(k) ;
double sum = 0.0;
int n = thisCut->row().getNumElements();
const int * column = thisCut->row().getIndices();
const double * element = thisCut->row().getElements();
assert (n);
for (int i = 0; i < n; i++) {
double value = element[i];
sum += value * solution[column[i]];
}
if (sum > thisCut->ub() + primalTolerance) {
sum = sum - thisCut->ub();
} else if (sum < thisCut->lb() - primalTolerance) {
sum = thisCut->lb() - sum;
} else {
sum = 0.0;
}
if (!sum) {
// take out
cs.eraseRowCut(k);
}
}
#ifdef CLP_INVESTIGATE2
printf("OR %p pass %d inTree %c - %d cuts (but %d deleted)\n",
originalSolver_,info.pass,info.inTree?'Y':'N',
numberRowCutsAfter-numberRowCutsBefore,
numberRowCutsAfter-cs.sizeRowCuts());
#endif
}
int numberRowCutsAfter = cs.sizeRowCuts();
if (!info.inTree) {
for (int i=numberRowCutsBefore;i<numberRowCutsAfter;i++) {
cs.rowCutPtr(i)->setGloballyValid();
}
}
if (twomirType==2) {
// back to original
int numberRows = clpSolver->getNumRows();
if (numberRows>numberOriginalRows) {
int numberDelete = numberRows-numberOriginalRows;
int * delRow = new int [numberDelete];
for (int i=0;i<numberDelete;i++)
delRow[i]=i+numberOriginalRows;
clpSolver->deleteRows(numberDelete,delRow);
delete [] delRow;
}
}
#endif
}
//-------------------------------------------------------------------
// Default Constructor
//-------------------------------------------------------------------
CglTwomir::CglTwomir () :
CglCutGenerator(),
probname_(),
randomNumberGenerator_(987654321),originalSolver_(NULL),
away_(0.0005),awayAtRoot_(0.0005),twomirType_(0),
do_mir_(true), do_2mir_(true), do_tab_(true), do_form_(true),
t_min_(1), t_max_(1), q_min_(1), q_max_(1), a_max_(2),max_elements_(50000),
max_elements_root_(50000),form_nrows_(0) {}
//-------------------------------------------------------------------
// Copy constructor
//-------------------------------------------------------------------
CglTwomir::CglTwomir (const CglTwomir & source) :
CglCutGenerator(source),
randomNumberGenerator_(source.randomNumberGenerator_),
originalSolver_(NULL),
away_(source.away_),
awayAtRoot_(source.awayAtRoot_),
twomirType_(source.twomirType_),
do_mir_(source.do_mir_),
do_2mir_(source.do_2mir_),
do_tab_(source.do_tab_),
do_form_(source.do_form_),
t_min_(source.t_min_),
t_max_(source.t_max_),
q_min_(source.q_min_),
q_max_(source.q_max_),
a_max_(source.a_max_),
max_elements_(source.max_elements_),
max_elements_root_(source.max_elements_root_),
form_nrows_(source.form_nrows_)
{
probname_ = source.probname_ ;
if (source.originalSolver_)
originalSolver_ = source.originalSolver_->clone();
}
//-------------------------------------------------------------------
// Clone
//-------------------------------------------------------------------
CglCutGenerator *
CglTwomir::clone() const
{
return new CglTwomir(*this);
}
//-------------------------------------------------------------------
// Destructor
//-------------------------------------------------------------------
CglTwomir::~CglTwomir ()
{
delete originalSolver_;
}
//----------------------------------------------------------------
// Assignment operator
//-------------------------------------------------------------------
CglTwomir &
CglTwomir::operator=(const CglTwomir& rhs)
{
if (this != &rhs) {
CglCutGenerator::operator=(rhs);
randomNumberGenerator_ = rhs.randomNumberGenerator_;
away_=rhs.away_;
awayAtRoot_=rhs.awayAtRoot_;
twomirType_ = rhs.twomirType_;
delete originalSolver_;
if (rhs.originalSolver_)
originalSolver_ = rhs.originalSolver_->clone();
else
originalSolver_=NULL;
do_mir_=rhs.do_mir_;
do_2mir_=rhs.do_2mir_;
do_tab_=rhs.do_tab_;
do_form_=rhs.do_form_;
t_min_=rhs.t_min_;
t_max_=rhs.t_max_;
q_min_=rhs.q_min_;
q_max_=rhs.q_max_;
a_max_=rhs.a_max_;
max_elements_=rhs.max_elements_;
max_elements_root_ = rhs.max_elements_root_;
form_nrows_=rhs.form_nrows_;
}
return *this;
}
// Pass in a copy of original solver (clone it)
void
CglTwomir::passInOriginalSolver(OsiSolverInterface * solver)
{
delete originalSolver_;
if (solver) {
if (!twomirType_)
twomirType_=1;
originalSolver_ = solver->clone();
originalSolver_->setHintParam(OsiDoDualInResolve, false, OsiHintDo);
// Temp - check if free variables
const double *colUpper = originalSolver_->getColUpper();
const double *colLower = originalSolver_->getColLower();
int ncol = originalSolver_->getNumCols();
int numberFree=0;
for (int i=0;i<ncol;i++) {
if (colLower[i]<-1.0e20&&colUpper[i]>1.0e20)
numberFree++;
}
if (numberFree)
printf("CglTwoMir - %d free variables - take care\n",numberFree);
} else {
twomirType_=0;
originalSolver_=NULL;
}
}
int DGG_freeData( DGG_data_t *data )
{
free(data->info);
free(data->lb);
free(data->ub);
free(data->x);
free(data->rc);
free(data);
return 0;
}
DGG_data_t* DGG_getData(const void *osi_ptr )
{
DGG_data_t *data = NULL;
const OsiSolverInterface *si = reinterpret_cast<const OsiSolverInterface *> (osi_ptr);
data = reinterpret_cast<DGG_data_t*> (malloc( sizeof(DGG_data_t)) );
/* retrieve basis information */
CoinWarmStart *startbasis = si->getWarmStart();
const CoinWarmStartBasis *basis = dynamic_cast<const CoinWarmStartBasis*>(startbasis);
/* retrieve bounds information */
const double *colUpper = si->getColUpper();
const double *colLower = si->getColLower();
const double *rowUpper = si->getRowUpper();
const double *rowLower = si->getRowLower();
const double *redCost = si->getReducedCost();
const double *dualVal = si->getRowPrice();
/* retrieve current optimal solution */
const double *colSolut = si->getColSolution();
/* retrieve the matrix in row format */
const CoinPackedMatrix *rowMatrixPtr = si->getMatrixByRow();
const int *rowBeg = 0, *rowCnt = 0, *rowInd = 0;
const double *rowMat;
rowBeg = rowMatrixPtr->getVectorStarts();
rowCnt = rowMatrixPtr->getVectorLengths();
rowMat = rowMatrixPtr->getElements();
rowInd = rowMatrixPtr->getIndices();
/* set number of columns and number of rows */
data->ncol = si->getNumCols();
data->nrow = si->getNumRows();
/* set ninteger */
data->ninteger = 0;
/* allocate memory for the arrays in 'data' */
data->info = reinterpret_cast<int*> (malloc( sizeof(int)*(data->ncol+data->nrow)) );
data->lb = reinterpret_cast<double*> (malloc( sizeof(double)*(data->ncol+data->nrow)) );
data->ub = reinterpret_cast<double*> (malloc( sizeof(double)*(data->ncol+data->nrow)) );
data->x = reinterpret_cast<double*> (malloc( sizeof(double)*(data->ncol+data->nrow)) );
data->rc = reinterpret_cast<double*> (malloc( sizeof(double)*(data->ncol+data->nrow)) );
memset(data->info, 0, sizeof(int)*(data->ncol+data->nrow));
/* set parameters for column variables */
data->nbasic_col = 0;
for(int i=0; i < data->ncol; i++){
/* is variable basic */
if ( basis->getStructStatus(i) == CoinWarmStartBasis::basic ){
data->nbasic_col++;
DGG_setIsBasic(data,i);
}
#if DGG_DEBUG_DGG
{
int error = 0;
if ( basis->getStructStatus(i) != CoinWarmStartBasis::basic )
if ( fabs(colSolut[i] - colUpper[i]) > DGG_BOUND_THRESH )
if ( fabs(colSolut[i] - colLower[i]) > DGG_BOUND_THRESH ){
fprintf(stdout, "WARNING!!!! : ");
fprintf(stdout, "variable %d non-basic, lb = %f, ub = %f, x = %f\n",
i, colLower[i], colUpper[i], colSolut[i]);
error+=1;
}
if (error)
fprintf(stdout, "\nFOUND %d errors. BYE.\n", error);
}
#endif
/* set variable bounds*/
data->lb[i] = colLower[i];
data->ub[i] = colUpper[i];
/* is variable integer */
if ( si->isInteger(i) ){
data->ninteger++;
DGG_setIsInteger(data,i);
/* tighten variable bounds*/
data->lb[i] = ceil(colLower[i]);
data->ub[i] = floor(colUpper[i]);
}
/* set x value */
data->x[i] = colSolut[i];
/* WARNING: remember to set rc!! Its not set!! */
data->rc[i] = redCost[i];
}
/* set parameters for row variables */
/* slack variables (row variables) work as follows:
for a ranged constraint, b_dw < ax < b_up, define a variable s so that
1) if b_up is not infinity:
ax + s = b_up, 0 < s < b_up - b_dw
2) if b_up is infinity:
ax - s = b_dw, 0 < s < b_up - b_dw
*/
{
int i,j;
double activity;
data->nbasic_row = 0;
for(i=0, j=data->ncol; i < data->nrow; i++, j++){
/* check if the row is an equality constraint */
if ( fabs( rowUpper[i] - rowLower[i] ) <= DGG_BOUND_THRESH )
DGG_setEqualityConstraint(data,j);
/* check if the row is bounded above/below and define variable bounds */
if ( rowUpper[i] < COIN_DBL_MAX )
DGG_setIsConstraintBoundedAbove(data,j);
if ( rowLower[i] > -1*COIN_DBL_MAX )
DGG_setIsConstraintBoundedBelow(data,j);
data->lb[j] = 0.0;
if (DGG_isConstraintBoundedAbove(data,j) && DGG_isConstraintBoundedBelow(data,j))
data->ub[j] = rowUpper[i] - rowLower[i];
else
data->ub[j] = COIN_DBL_MAX;
/* compute row activity. for this we need to go to the row in question,
and multiply all the coefficients times their respective variables.
For the moment, we will store the inverse of this value in
the 'x' field (since it is in fact a partial computation of it) */
{
int k;
activity = 0.0;
for(k=rowBeg[i]; k < rowBeg[i]+rowCnt[i]; k++)
activity += rowMat[k]*colSolut[rowInd[k]];
}
/* compute x value */
if ( DGG_isConstraintBoundedAbove(data,j) )
data->x[j] = rowUpper[i] - activity;
else
data->x[j] = activity - rowLower[i];
if ( data->x[j] < -DGG_NULL_SLACK ){
#if DGG_DEBUG_DGG
int k;
double norm = 0.0, min = DBL_MAX, amin = DBL_MAX, max = DBL_MIN;
printf("** warning: row %d has negative slack!\n", i);
for(k=rowBeg[i]; k < rowBeg[i]+rowCnt[i]; k++){
norm += rowMat[k]*rowMat[k];
if ( fabs(rowMat[k]) < amin ) amin = fabs(rowMat[k]);
if ( rowMat[k] < min ) min = rowMat[k];
if ( rowMat[k] > max ) max = rowMat[k];
}
norm = sqrt(norm);
printf("min = %f amin = %f max = %f\n", min, amin, max);
printf("rlower = %f activity = %f\n", rowLower[i], activity);
printf("norm = %f (b-ax) = %f\n", norm, (rowLower[i] - activity));
printf("steepn = %f\n", (rowLower[i] - activity)/norm);
#endif
}
data->rc[j] = dualVal[i];
#if DGG_DEBUG_SOLVER
DGG_IF_EXIT( !DGG_isConstraintBoundedAbove(data,j) && !DGG_isConstraintBoundedBelow(data,j),
1, "some row is not bounded above or below");
#endif
/* is variable basic */
if ( basis->getArtifStatus(i) == CoinWarmStartBasis::basic ){
data->nbasic_row++;
DGG_setIsBasic(data,j);
}
/* is variable integer. For this we need to go to the row in question,
and check that the rhs is integer, and that all of the coefficients
and variables participating in the constraint are also integer. */
{
int k;
if( DGG_isConstraintBoundedAbove(data,j)) {
if ( frac_part(rowUpper[i]) > DGG_INTEGRALITY_THRESH )
goto DONE_ROW;
}
else
if ( frac_part(rowLower[i]) > DGG_INTEGRALITY_THRESH )
goto DONE_ROW;
for(k=rowBeg[i]; k < rowBeg[i]+rowCnt[i]; k++)
if ( frac_part(rowMat[k]) > DGG_INTEGRALITY_THRESH || !DGG_isInteger(data, rowInd[k]))
goto DONE_ROW;
DGG_setIsInteger(data, j);
data->ninteger++;
}
DONE_ROW:;
/* set variable bounds: careful!! Later, remember to adjust
the INFINITY to a DGG standard (to deal with neq solvers). */
/* WARNING: remember to set rc!! Its not set!! */
}
}
/* CLEANUP */
delete basis;
return data;
}
DGG_constraint_t*
DGG_getSlackExpression(const void *osi_ptr, DGG_data_t* data, int row_index)
{
DGG_constraint_t *row = 0;
int i,j;
/* retrieve the matrix in row format */
const OsiSolverInterface *si = reinterpret_cast<const OsiSolverInterface *> (osi_ptr);
const CoinPackedMatrix *rowMatrixPtr = si->getMatrixByRow();
const int *rowBeg = 0, *rowCnt = 0, *rowInd = 0;
const double *rowMat;
const double *rowUpper;
const double *rowLower;
row = DGG_newConstraint(data->ncol);
rowBeg = rowMatrixPtr->getVectorStarts();
rowCnt = rowMatrixPtr->getVectorLengths();
rowMat = rowMatrixPtr->getElements();
rowInd = rowMatrixPtr->getIndices();
rowUpper = si->getRowUpper();
rowLower = si->getRowLower();
#if DGG_DEBUG_DGG
if ( row_index < 0 || row_index > data->nrow )
DGG_THROW(0, "bad row index");
#endif
/* copy the information into the row ADT */
row->nz = rowCnt[row_index];
for(j=0, i=rowBeg[row_index]; i < rowBeg[row_index]+rowCnt[row_index]; i++, j++){
row->coeff[j] = rowMat[i];
row->index[j] = rowInd[i];
if (DGG_isConstraintBoundedAbove (data, data->ncol + row_index))
row->coeff[j] = -row->coeff[j];
}
row->sense = '?';
if ( DGG_isConstraintBoundedAbove(data, data->ncol + row_index) )
row->rhs = rowUpper[row_index];
else
row->rhs = -rowLower[row_index];
return row;
}
int
DGG_getTableauConstraint( int index, const void *osi_ptr, DGG_data_t *data,
DGG_constraint_t* tabrow,
const int * colIsBasic,
const int * /*rowIsBasic*/,
CoinFactorization & factorization,
int mode )
{
#if DGG_DEBUG_DGG
/* ensure that the index corresponds to a basic variable */
if ( !DGG_isBasic(data, index) )
DGG_THROW(1, "index is non-basic");
/* ensure that the index corresponds to a column variable */
if ( index < 0 || index > (data->ncol - 1) )
DGG_THROW(1, "index not a column variable");
#endif
/* obtain pointer to solver interface */
const OsiSolverInterface *si = reinterpret_cast<const OsiSolverInterface *> (osi_ptr);
DGG_TEST(!si, 1, "null OsiSolverInterfave");
/* obtain address of the LP matrix */
const CoinPackedMatrix *colMatrixPtr = si->getMatrixByCol();
const int* colBeg = colMatrixPtr->getVectorStarts();
const int* colCnt = colMatrixPtr->getVectorLengths();
const int* colInd = colMatrixPtr->getIndices();
const double* colMat = colMatrixPtr->getElements();
/* obtain row right-hand-sides */
const double *rowUpper = si->getRowUpper();
const double *rowLower = si->getRowLower();
/* allocate memory for constraint in non-sparse form */
int nz = 0;
double *value = NULL, rhs = 0.0;
value = reinterpret_cast<double*>(malloc(sizeof(double)*(data->nrow+data->ncol)));
memset(value, 0, sizeof(double)*(data->nrow+data->ncol));
/* obtain the tableau row coefficients for all variables */
/* note: we could speed this up by only computing non-basic variables */
{
int i, j, cnt = 0;
double one = 1.0;
CoinIndexedVector work;
CoinIndexedVector array;
work.reserve(data->nrow);
array.reserve(data->nrow);
array.setVector(1,&colIsBasic[index],&one);
factorization.updateColumnTranspose ( &work, &array );
int * arrayRows = array.getIndices();
double *arrayElements = array.denseVector();
cnt = array.getNumElements();
/* compute column (structural) variable coefficients */
for(j = 0; j < data->ncol; j++) {
value[j] = 0.0;
for(i=colBeg[j]; i < colBeg[j]+colCnt[j]; i++)
value[j] += colMat[i]*arrayElements[ colInd[i] ];
}
#if DGG_DEBUG_SOLVER
/* check pivot */
if ( fabs(value[index] - 1.0) > DGG_INTEGRALITY_THRESH )
DGG_THROW(1, "pivot is not one");
#endif
/* compute row variable (slack/logical) variable coefficients */
for(j = 0; j < cnt; j++){
if ( DGG_isEqualityConstraint(data,data->ncol + arrayRows[j]) && !mode )
value[ data->ncol + arrayRows[j] ] = 0.0;
else if ( DGG_isConstraintBoundedAbove(data, data->ncol + arrayRows[j]) )
value[ data->ncol + arrayRows[j] ] = arrayElements[ arrayRows[j] ];
else
value[ data->ncol + arrayRows[j] ] = -1*arrayElements[ arrayRows[j] ];
}
/* compute rhs */
rhs = 0.0;
for(i=0; i < cnt; i++){
if ( DGG_isConstraintBoundedAbove(data,data->ncol + arrayRows[i]) )
rhs += arrayElements[arrayRows[i]]*rowUpper[arrayRows[i]];
else
rhs += arrayElements[arrayRows[i]]*rowLower[arrayRows[i]];
}
/* free 'work' and 'array' ?? do the empty functions do it?? they are not
cleared in CglGomory. Is that due to a mistake? Is it done on purpose? */
/*
work.empty();
array.empty();
*/
}
/* count non-zeroes */
nz = 0;
int j;
for( j=0; j<data->ncol+data->nrow; j++){
if ( fabs(value[j]) > DGG_MIN_TABLEAU_COEFFICIENT )
nz += 1;
}
/* put in sparse constraint format */
/* technical issue: should we let max_nz == nz or should we instead
set max_nz == (nrow+ncol). The advantage of the latter approach
is that later, when we substitute the slacks, the denser column
will not require us to re-allocate memory */
tabrow->max_nz = nz;
if (tabrow->coeff)
free(tabrow->coeff);
if (tabrow->index)
free(tabrow->index);
tabrow->coeff = reinterpret_cast<double*> (malloc(sizeof(double)*nz));
tabrow->index = reinterpret_cast<int*> (malloc(sizeof(int)*nz));
tabrow->nz = 0;
for( j = 0; j < data->ncol + data->nrow; j++)
if ( fabs(value[j]) > DGG_MIN_TABLEAU_COEFFICIENT ){
tabrow->coeff[tabrow->nz] = value[j];
tabrow->index[tabrow->nz] = j;
tabrow->nz += 1;
}
tabrow->sense = 'E';
tabrow->rhs = rhs;
/* CLEANUP */
free(value);
return 0;
}
int
DGG_getFormulaConstraint( int da_row,
const void *osi_ptr,
DGG_data_t *data,
DGG_constraint_t* form_row )
{
/* ensure that the da_row corresponds to a row */
if ( data->nrow <= da_row || 0> da_row) DGG_THROW(1, "row out of range...");
/* obtain pointer to solver interface */
const OsiSolverInterface *si = reinterpret_cast<const OsiSolverInterface *> (osi_ptr);
//DGG_TEST(!si, 1, "null OsiSolverInterfave");
/* obtain address of the LP matrix */
const CoinPackedMatrix *rowMatrixPtr = si->getMatrixByRow();
const int* rowBeg = rowMatrixPtr->getVectorStarts();
const int* rowCnt = rowMatrixPtr->getVectorLengths();
const int* rowInd = rowMatrixPtr->getIndices();
const double* rowMat = rowMatrixPtr->getElements();
/* obtain row right-hand-sides */
const double *rowUpper = si->getRowUpper();
const double *rowLower = si->getRowLower();
int nz = rowCnt[da_row];
form_row->nz = nz;
form_row->max_nz = nz+1;
int i;
for( i=0; i < nz; i++) form_row->coeff[i] = rowMat[rowBeg[da_row]+i];
for( i=0; i < nz; i++) form_row->index[i] = rowInd[rowBeg[da_row]+i];
if ( DGG_isConstraintBoundedAbove(data,data->ncol + da_row) ){
form_row->rhs = rowUpper[da_row];
form_row->sense = 'L';
}
else{
form_row->rhs = rowLower[da_row];
form_row->sense = 'G';
}
if ( DGG_isEqualityConstraint(data,data->ncol + da_row) )
form_row->sense = 'E';
/* now add slack/surplus if there is one */
if ( DGG_isEqualityConstraint(data,data->ncol + da_row) == 0 ){
form_row->index[nz] = data->ncol + da_row;
if ( DGG_isConstraintBoundedAbove(data, data->ncol + da_row) )
form_row->coeff[nz] = 1;
else
form_row->coeff[nz] = -1;
form_row->nz +=1;
}
return 0;
}
//---------------------------------------------------------------
//---------------------------------------------------------------
//---------------------------------------------------------------
//---------------------------------------------------------------
//---------------------------------------------------------------
/******************** CONSTRAINT ADTs *****************************************/
DGG_constraint_t* DGG_newConstraint(int max_arrays)
{
DGG_constraint_t *c = NULL;
if (max_arrays <= 0) return NULL;
c = reinterpret_cast<DGG_constraint_t*> (malloc(sizeof(DGG_constraint_t)));
c->nz = 0;
c->max_nz = max_arrays;
c->rhs = 0.0;
c->sense = '?';
c->coeff = NULL;
c->index = NULL;
c->coeff = reinterpret_cast<double*>(malloc(sizeof(double)*max_arrays));
c->index = reinterpret_cast<int*>(malloc(sizeof(int)*max_arrays));
return c;
}
void DGG_freeConstraint(DGG_constraint_t *c)
{
if (c == NULL) return;
if (c->coeff) free(c->coeff);
if (c->index) free(c->index);
free(c);
}
DGG_constraint_t *DGG_copyConstraint(DGG_constraint_t* c)
{
DGG_constraint_t *nc = NULL;
if (!c || c->max_nz <= 0) return nc;
nc = DGG_newConstraint(c->max_nz);
if (nc == NULL) return nc;
nc->nz = c->nz;
nc->rhs = c->rhs;
nc->sense = c->sense;
memcpy(nc->coeff, c->coeff, sizeof(double)*nc->nz);
memcpy(nc->index, c->index, sizeof(int)*nc->nz);
return nc;
}
void DGG_scaleConstraint(DGG_constraint_t *c, int t)
{
int i;
c->rhs *= t;
if (t < 0){
if (c->sense == 'G') c->sense = 'L';
else if (c->sense == 'L') c->sense = 'G';
}
for(i=0; i<c->nz; i++) c->coeff[i] *= t;
}
void DGG_list_init (DGG_list_t *l)
{
l->n = 0;
l->c = NULL;
l->ctype = NULL;
l->alpha = NULL;
}
void DGG_list_free(DGG_list_t *l)
{
if (l->c != NULL) free (l->c);
if (l->ctype != NULL) free (l->ctype);
if (l->alpha != NULL) free (l->alpha);
}
int DGG_list_addcut (DGG_list_t *l, DGG_constraint_t *cut, int ctype, double alpha)
{
l->n ++;
l->c = reinterpret_cast<DGG_constraint_t **>(realloc (l->c, l->n * sizeof(DGG_constraint_t *)));
l->ctype = reinterpret_cast<int *>(realloc (l->ctype, l->n * sizeof (int)));
l->alpha = reinterpret_cast<double *>(realloc (l->alpha, l->n * sizeof (double)));
if (l->c == NULL || l->ctype == NULL || l->alpha == NULL){
printf ("No memory, bailing out\n");
return -1;
}
l->c[l->n - 1] = cut;
l->ctype[l->n - 1] = ctype;
l->alpha[l->n - 1] = alpha;
return 0;
}
void DGG_list_delcut (DGG_list_t *l, int i)
{
if (i >= l->n && i < 0) return;
DGG_freeConstraint (l->c[i]);
l->c[i] = l->c[l->n - 1];
l->ctype[i] = l->ctype[l->n - 1];
l->alpha[i] = l->alpha[l->n - 1];
l->n --;
}
/******************* CONSTRAINT MANIPULATION **********************************/
/* VARIABLES CLOSE TO UPPER BOUNDS:
we will substitute: x' = (u - x); hence the constraint will change from ax ~ b
to -ax' ~ b - au note: the new bounds of x' will be, 0 <= x' <= u - l
VARIABLES CLOSE TO LOWER BOUNDS:
we will substitute: x' = (x - l); hence, the constraint will change from
ax ~ b to ax' ~ b - al. note: some variable lower bounds may have changed
when doing the complement in the previous stage - this must be taken into
account. note: the new bounds of x' will be, 0 <= x' <= u - l */
int DGG_transformConstraint( DGG_data_t *data,
double **x_out,
double **rc_out,
char **isint_out,
DGG_constraint_t *constraint )
{
double half;
double *px = reinterpret_cast<double*> (malloc( sizeof(double)*constraint->max_nz ));
double *rc = reinterpret_cast<double*> (malloc( sizeof(double)*constraint->max_nz ));
char *pi = reinterpret_cast<char*> (malloc( sizeof(char) *constraint->max_nz ));
{
int i, idx;
for(i=0; i < constraint->nz; i++){
idx = constraint->index[i];
px[i] = data->x[idx];
rc[i] = data->rc[idx];
pi[i] = static_cast<char>(DGG_isInteger(data, idx));
half = (data->ub[idx] - data->lb[idx]) / 2;
if ( data->ub[idx] - data->x[idx] < half ){
px[i] = data->ub[idx] - data->x[idx];
if (fabs(px[i]) <= DGG_BOUND_THRESH)
px[i] = 0.0;
constraint->rhs -= constraint->coeff[i]*data->ub[idx];
constraint->coeff[i] *= -1;
}
else {
px[i] = data->x[idx] - data->lb[idx];
if (fabs(px[i]) <= DGG_BOUND_THRESH)
px[i] = 0.0;
constraint->rhs -= constraint->coeff[i]*data->lb[idx];
}
}
}
*x_out = px;
*rc_out = rc;
*isint_out = pi;
#if DGG_DEBUG_DGG
DGG_TEST(DGG_isConstraintViolated(data, constraint), 1, "bad transformation");
#endif
return 0;
}
int DGG_unTransformConstraint( DGG_data_t *data,
DGG_constraint_t *constraint )
{
int i, idx;
double half;
for(i=0; i < constraint->nz; i++){
idx = constraint->index[i];
half = (data->ub[idx] - data->lb[idx]) / 2;
if ( data->ub[idx] - data->x[idx] < half ){
constraint->rhs -= constraint->coeff[i]*data->ub[idx];
constraint->coeff[i] *= -1;
}
else
constraint->rhs += constraint->coeff[i]*data->lb[idx];
}
return 0;
}
int
DGG_substituteSlacks( const void *solver_ptr,
DGG_data_t *data,
DGG_constraint_t *cut )
{
int i,j, lnz;
double *lcut, lrhs;
DGG_constraint_t *row=NULL;
/* lcut will store all the column coefficients. allocate space and init. */
lcut = reinterpret_cast<double*>(malloc(sizeof(double)*data->ncol));
memset(lcut, 0, sizeof(double)*data->ncol);
/* initialize lrhs */
lrhs = cut->rhs;
/* set coefficients in lcut */
/* technical: we could speed this up by re-using allocated memory
for row->coeff and row->index */
for(i=0; i < cut->nz; i++){
if ( cut->index[i] < data->ncol )
lcut[ cut->index[i] ] += cut->coeff[i];
else{
row = DGG_getSlackExpression(solver_ptr, data, (cut->index[i] - data->ncol));
for(j=0; j < row->nz; j++)
lcut[ row->index[j] ] += row->coeff[j]*cut->coeff[i];
lrhs -= row->rhs*cut->coeff[i];
DGG_freeConstraint(row);
}
}
/* count nz in new constraint */
lnz = 0;
for(i=0; i < data->ncol; i++)
if ( fabs(lcut[i]) > DGG_MIN_TABLEAU_COEFFICIENT )
lnz += 1;
/* free row->coeff and row->index, and re-allocate */
free(cut->coeff); cut->coeff = 0;
free(cut->index); cut->index = 0;
cut->nz = lnz;
cut->max_nz = lnz;
if (lnz)
{
cut->coeff = reinterpret_cast<double*> (malloc( sizeof(double)*lnz ));
cut->index = reinterpret_cast<int*> (malloc( sizeof(int)*lnz ));
}
/* set new constraint */
lnz = 0;
for(i=0; i < data->ncol; i++){
if ( fabs(lcut[i]) > DGG_MIN_TABLEAU_COEFFICIENT ){
cut->coeff[lnz] = lcut[i];
cut->index[lnz] = i;
lnz += 1;
}
}
cut->rhs = lrhs;
free(lcut);
return 0;
}
int DGG_nicefyConstraint( const void * /*solver_ptr*/,
DGG_data_t *data,
DGG_constraint_t *cut)
{
double min_coef = COIN_DBL_MAX, max_coef = COIN_DBL_MIN;
DGG_TEST(cut->sense == 'L', 1, "can't nicefy an L constraint");
int i;
for( i=0; i<cut->nz; i++) // first clean out noise
if( fabs(cut->coeff[i]) < DGG_NICEFY_MIN_ABSVALUE)
cut->coeff[i] = 0;
for( i=0; i<cut->nz; i++){
if( DGG_isInteger(data, cut->index[i])){// look at integral vars.
double aht = ABOV(cut->coeff[i]);
double ub = data->ub[ cut->index[i]];
if(aht < DGG_NICEFY_MIN_FIX){// coefficient = integer + epsylon
cut->coeff[i] = floor( cut->coeff[i]);
double ahtu = aht * ub;
if(ahtu<DGG_NICEFY_MAX_PADDING)
cut->rhs -= ahtu;// safely remove the fractional part
else
cut->coeff[i] += DGG_NICEFY_MIN_FIX; // inflate the fractional part
}
else
if (1-aht < DGG_NICEFY_MIN_FIX) // coefficient = integer - epsylon
cut->coeff[i] = ceil( cut->coeff[i]);
}// done with integers
else // now look at continuous variables
if ( cut->coeff[i] < DGG_NICEFY_MIN_ABSVALUE) // delete all negative and noise
cut->coeff[i] = 0.0;
else
if(cut->coeff[i] < DGG_NICEFY_MIN_FIX) {// coefficient = epsylon
double au = cut->coeff[i] * data->ub[ cut->index[i]];
if(au<DGG_NICEFY_MAX_PADDING){ // safely remove the variable
cut->coeff[i] = 0.0;
cut->rhs -= au;
}
else
cut->coeff[i] = DGG_NICEFY_MIN_FIX; // inflate the coefficient
}// done with continuous variables too
double abs_coef = fabs(cut->coeff[i]);
min_coef = DGG_MIN(min_coef, abs_coef);
max_coef = DGG_MAX(max_coef, abs_coef);
}
cut->sense = 'G';
/*
if ( max_coef > DGG_NICEFY_MAX_RATIO*min_coef ) // kill the cut if numbers are all over the place
cut->nz = 0;
*/
return 0;
}
/******************* CUT GENERATION *******************************************/
int
DGG_generateTabRowCuts( DGG_list_t *cut_list,
DGG_data_t *data,
const void *solver_ptr )
{
int k, rval = 0;
DGG_constraint_t *base = NULL;
int nc = cut_list->n;
base = DGG_newConstraint(data->ncol + data->nrow);
if(talk) printf ("2mir_test: generating tab row cuts\n");
/* allocate memory for basic column/row indicators */
int *rowIsBasic = 0, *colIsBasic = 0;
rowIsBasic = reinterpret_cast<int*>(malloc(sizeof(int)*data->nrow));
colIsBasic = reinterpret_cast<int*>(malloc(sizeof(int)*data->ncol));
/* initialize the IsBasic arrays with -1 / 1 values indicating
where the basic rows and columns are. NOTE: WE could do this
only once and keep it in osi_data at the expense of space!! */
int i;
for( i=0; i<data->ncol; i++){
if ( DGG_isBasic(data,i) ) colIsBasic[i] = 1;
else colIsBasic[i] = -1;
}
for( i=0; i<data->nrow; i++){
if ( DGG_isBasic(data,i+data->ncol) ) rowIsBasic[i] = 1;
else rowIsBasic[i] = -1;
}
/* obtain factorization */
CoinFactorization factorization;
/* obtain address of the LP matrix */
const OsiSolverInterface *si = reinterpret_cast<const OsiSolverInterface *> (solver_ptr);
const CoinPackedMatrix *colMatrixPtr = si->getMatrixByCol();
rval = factorization.factorize(*colMatrixPtr, rowIsBasic, colIsBasic);
/* 0 = okay. -1 = singular. -2 = too many in basis. -99 = memory. */
DGG_TEST2(rval, 1, "factorization error = %d", rval);
for(k=0; k<data->ncol; k++){
if (!(DGG_isBasic(data, k) && DGG_isInteger(data,k))) continue;
double frac = frac_part (data->x[k]);
if (frac < data->gomory_threshold || frac > 1-data->gomory_threshold) continue;
base->nz = 0;
rval = DGG_getTableauConstraint(k, solver_ptr, data, base,
colIsBasic,rowIsBasic,factorization,0);
DGG_CHECKRVAL(rval, rval);
if (base->nz == 0){
printf ("2mir_test: why does constraint not exist ?\n");
continue;
}
if (base->nz > 500) continue;
rval = DGG_generateCutsFromBase(base, cut_list, data, solver_ptr);
DGG_CHECKRVAL(rval, rval);
}
free(rowIsBasic);
free(colIsBasic);
if(talk) printf ("2mir_test: generated %d tab cuts\n", cut_list->n - nc); fflush (stdout);
DGG_freeConstraint(base);
return rval;
}
int DGG_generateFormulationCuts( DGG_list_t *cut_list,
DGG_data_t *data,
const void *solver_ptr,
int nrows,
CoinThreadRandom & generator)
{
int k, rval = 0;
DGG_constraint_t *base = NULL;
int num_rows = (data->nrow < nrows) ? data->nrow : nrows;
int nc = cut_list->n;
base = DGG_newConstraint(data->ncol + data->nrow);
if(talk) printf ("2mir_test: generating form row cuts %d\n", num_rows);
for(k=0; k<num_rows; k++) {
base->nz = 0;
rval = DGG_getFormulaConstraint(k, solver_ptr, data, base);
DGG_CHECKRVAL1(rval, rval);
//printf ("generating formulation for row %d\n", k);
rval = DGG_generateFormulationCutsFromBase(base, data->x[data->ncol+k],
cut_list, data, solver_ptr,
generator);
DGG_CHECKRVAL1(rval, rval);
if (base->nz == 0){
#ifdef COIN_DEVELOP
printf ("why does constraint not exist ?\n");
#endif
continue;
}
}
CLEANUP:
if(talk) printf ("2mir_test: generated %d form cuts\n", cut_list->n - nc); fflush (stdout);
DGG_freeConstraint(base);
return rval;
}
int DGG_generateFormulationCutsFromBase( DGG_constraint_t *base,
double slack,
DGG_list_t *cut_list,
DGG_data_t *data,
const void *solver_ptr,
CoinThreadRandom & generator)
{
int i, p, rval;
int int_skala;
double skala;
int num_inlist = 0;
int* skala_list = reinterpret_cast<int*> (malloc( sizeof(int)*base->nz ));
char *isint = NULL;
double *xout = NULL, *rcout = NULL;
DGG_constraint_t *scaled_base = NULL;
int tot_int = 0;
double prob_choose = 0.0;
rval = DGG_transformConstraint(data, &xout, &rcout, &isint, base);
DGG_CHECKRVAL1(rval, rval);
for(p = 0; p < base->nz; p++) if(isint[p]) tot_int ++;
if (tot_int == 0) goto CLEANUP;
prob_choose = 5.0/tot_int;
for(p = 0; p < base->nz; p++) {
if(isint[p]) if(generator.randomDouble() < prob_choose){
if(xout[p]<0.01) continue;
skala =fabs(base->coeff[p]);
if(skala<0.01) continue;
// check if slack is too large
if (fabs(slack/skala) > 0.5) continue;
scaled_base = DGG_copyConstraint(base);
DGG_CHECKRVAL1((scaled_base == NULL),-1);
if(base->sense == 'L') {
skala = -skala;
scaled_base->sense = 'G';
}
int_skala = int(100*skala);
for(i = 0; i< num_inlist; i++)
if(int_skala == skala_list[i])
goto END_LOOP;
skala_list[num_inlist++] = int_skala;
scaled_base->rhs = base->rhs/skala;
for(i = 0; i<base->nz; i++)
scaled_base->coeff[i] = base->coeff[i] / skala;
rval = DGG_unTransformConstraint(data, scaled_base);
DGG_CHECKRVAL1(rval, rval);
rval = DGG_generateCutsFromBase(scaled_base, cut_list,
data, solver_ptr);
DGG_CHECKRVAL1(rval, rval);
END_LOOP:
DGG_freeConstraint(scaled_base);
scaled_base = NULL;
}
}
CLEANUP:
if (isint) free(isint);
if (xout) free(xout);
if (rcout) free(rcout);
if (skala_list) free(skala_list);
if (scaled_base != NULL) DGG_freeConstraint (scaled_base);
return rval;
}
int
DGG_generateCutsFromBase( DGG_constraint_t *orig_base,
DGG_list_t *cut_list,
DGG_data_t *data,
const void *solver_ptr )
{
int rval = 0;
int t;
double *x = NULL, *rc = NULL;
char *isint = NULL;
DGG_constraint_t *base = NULL;
bool not_nicefied = true;
int new_pos = cut_list->n;
// DGG_constraint_t *keep_origbase = DGG_copyConstraint(orig_base); //for debug only ------
if (orig_base->sense == 'L') return 0;
if (orig_base->nz == 0) return 0;
rval = DGG_transformConstraint(data, &x, &rc, &isint, orig_base);
double frac = frac_part(orig_base->rhs);
//printf ("frac = %.7f, r %.7f, fr %.7f\n", frac, orig_base->rhs, floor(orig_base->rhs));
if (rval || frac < data->gomory_threshold || frac > 1-data->gomory_threshold){
free (x); free (rc); free (isint);
return 0;
}
int min_t = t_min;
int min_q = q_min;
if (orig_base->sense == 'G' && min_t < 1) min_t = 1;
if (orig_base->sense == 'G' && min_q < 1) min_q = 1;
if (min_q > 0 && min_t > 0 ) {
not_nicefied = false;
rval = DGG_nicefyConstraint(solver_ptr, data, orig_base);
DGG_CHECKRVAL(rval, rval);
if (orig_base->nz == 0){
if(talk) printf ("2mir_test: Nicefy returns empty constraint\n"); rval = 0; goto CLEANUP;
}
}
for(t = min_t; t <= t_max ; t++){
if (t == 0) continue;
base = DGG_copyConstraint(orig_base);
DGG_TEST(!base, 1, "error making copy of base");
DGG_scaleConstraint (base, t);
if(not_nicefied){
rval = DGG_nicefyConstraint(solver_ptr, data, base);
DGG_CHECKRVAL(rval, rval);
if (base->nz == 0){
if(talk) printf ("2mir_test: Nicefy returns empty constraint\n"); goto MIR_DONE;
}
}
if ( DGG_isBaseTrivial(data, base) ) goto MIR_DONE;
rval = DGG_addMirToList(base, isint, x, cut_list, data, orig_base);
DGG_CHECKRVAL(rval, rval);
MIR_DONE:
DGG_freeConstraint(base);
}
for( t = min_q; t <= q_max; t++ ){
if (t == 0) continue;
base = DGG_copyConstraint(orig_base);
DGG_TEST(!base, 1, "error making copy of base");
DGG_scaleConstraint (base, t);
if(not_nicefied){
rval = DGG_nicefyConstraint(solver_ptr, data, base);
DGG_CHECKRVAL(rval, rval);
if (base->nz == 0){
if(talk) printf ("2mir_test: Nicefy returns empty constraint\n"); goto TWOMIR_DONE;
}
}
if ( DGG_isBaseTrivial(data, base) ) goto TWOMIR_DONE;
rval = DGG_add2stepToList(base, isint, x, rc, cut_list, data, orig_base);
DGG_CHECKRVAL(rval, rval);
TWOMIR_DONE:
DGG_freeConstraint(base);
}
int i;
for ( i = cut_list->n-1; i>=new_pos; i--){
DGG_constraint_t *lcut = cut_list->c[i];
rval = DGG_unTransformConstraint(data, lcut);
DGG_CHECKRVAL(rval, rval);
rval = DGG_substituteSlacks(solver_ptr, data, lcut);
DGG_CHECKRVAL(rval, rval);
if ( !DGG_isCutDesirable(lcut, data) ){
DGG_list_delcut (cut_list, i);
continue;
}
//else testus(lcut);//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
/*
if ( data->opt_x && DGG_cutsOffPoint(data->opt_x, lcut) ){
//DGG_cutDisplay_sparse(data, data->opt_x, lcut, stdout);
DGG_TEST(1,1, "new cut is infeasible for optimal solution\n");
}
*/
}
CLEANUP:
if (x) free(x);
if (rc) free (rc);
if (isint) free(isint);
return 0;
}
int
DGG_addMirToList ( DGG_constraint_t *base, char *isint, double * /*x*/,
DGG_list_t *list, DGG_data_t * /*data*/,
DGG_constraint_t * /*orig_base*/ )
{
int rval = 0;
DGG_constraint_t *cut = NULL;
rval = DGG_buildMir(isint, base, &cut);
DGG_CHECKRVAL(rval, rval);
DGG_list_addcut(list, cut, DGG_TMIR_CUT, 0.0);
return 0;
}
int
DGG_add2stepToList ( DGG_constraint_t *base, char *isint, double * /*x*/,
double *rc, DGG_list_t *list, DGG_data_t *data,
DGG_constraint_t * /*orig_base*/ )
{
int rval;
DGG_constraint_t *cut = NULL;
int i;
double norm_val, best_norm_val, best_norm_alpha=-1.0;
double rc_val, best_rc_val, best_rc_alpha=-1.0;
double vht, bht, alpha;
best_rc_val = best_norm_val = COIN_DBL_MAX;
bht = ABOV(base->rhs);
double best_rc = 0;
for(i=0; i<base->nz; i++) if (isint[i]) best_rc = CoinMax(best_rc, fabs(rc[i]));
double rc_cutoff = best_rc / 10;
for(i=0; i<base->nz; i++){
if (!isint[i]) continue;
if (fabs(rc[i]) <= rc_cutoff) continue; //too unimportant
vht = ABOV(base->coeff[i]);
if(vht >= bht) continue; // too big
if(vht < bht/a_max) continue; // too small
alpha = vht;
int kk = 1;
while ( !DGG_is2stepValid(alpha, bht) && bht/alpha <= a_max) {
alpha = vht/kk;
kk++;
if (kk>1000)
break;
}
if ( !DGG_is2stepValid(alpha, bht) ) continue;
rval = DGG_build2step(alpha, isint, base, &cut);
DGG_CHECKRVAL(rval, rval);
rc_val = COIN_DBL_MAX; // this gives a lower bound on obj. fn. improvement
for(i=0; i<cut->nz; i++) if(cut->coeff[i]> 1E-6){
rc_val = CoinMin(rc_val, fabs(rc[i])/cut->coeff[i]);
}
rc_val *= cut->rhs;
norm_val = 0; // this is the square of the L2 norm
for(i=0; i<cut->nz; i++) if(cut->coeff[i]> 1E-6){
norm_val += (cut->coeff[i]*cut->coeff[i]);
}
norm_val /= cut->rhs * cut->rhs;
if (rc_val < best_rc_val ) {
best_rc_val = rc_val; best_rc_alpha = alpha; }
if (norm_val < best_norm_val ) {
best_norm_val = norm_val; best_norm_alpha = alpha; }
DGG_freeConstraint(cut);
}
if( best_rc_val> 1E-6 && best_rc_alpha != -1.0){
rval = DGG_build2step(best_rc_alpha, isint, base, &cut);
DGG_CHECKRVAL(rval, rval);
DGG_list_addcut(list, cut, DGG_2STEP_CUT, best_rc_alpha);
}
else if (best_norm_alpha != -1.0){
rval = DGG_build2step(best_norm_alpha, isint, base, &cut);
DGG_CHECKRVAL(rval, rval);
DGG_list_addcut(list, cut, DGG_2STEP_CUT, best_norm_alpha);
}
return 0;
}
int DGG_buildMir( char *isint,
DGG_constraint_t *base,
DGG_constraint_t **cut_out )
{
int i, lnz = 0;
double b = (base->rhs);
double bht = ABOV(b);
double bup = ceil(b);
DGG_constraint_t *tmir = NULL;
DGG_TEST( base->sense == 'L', 1, "this form not valid for L");
DGG_TEST( base->nz == 0, 1, "base must have some coefficients\n");
tmir = DGG_newConstraint( base->nz );
tmir->sense = 'G';
tmir->rhs = bht * bup;
for(i=0; i<base->nz; i++){
double v = base->coeff[i];
if (!isint[i]) {
if (v > 0.0) tmir->coeff[lnz] = v;
else tmir->coeff[lnz] = 0.0;
}
else {
double vht = ABOV(v);
DGG_IF_EXIT( vht<0, 1, "negative vht");
tmir->coeff[lnz] = bht * floor(v) + DGG_MIN(bht,vht);
}
tmir->index[lnz] = base->index[i];
lnz += 1;
}
tmir->nz = lnz;
*cut_out = tmir;
return 0;
}
int DGG_build2step( double alpha,
char *isint,
DGG_constraint_t *base,
DGG_constraint_t **cut_out )
{
DGG_constraint_t *tmir = 0;
int i, lnz = 0;
double vht, bht, bup, rho, tau, k;
double b = (base->rhs);
DGG_TEST( base->sense == 'L', 1, "this form not valid for L");
DGG_TEST( base->nz == 0, 1, "base must have some coefficients\n");
bht = ABOV(b);
bup = ceil(b);
tau = ceil(bht/alpha);
rho = bht - alpha*floor(bht/alpha);
/* ensure bht > alpha > 0 */
DGG_TEST3( (bht <= alpha) || (alpha <= 0.0), 1, "bad alpha (%f) / bht (%f) pair", alpha, bht);
/* ensure that we are not in a limiting case */
DGG_TEST( DGG_is_a_multiple_of_b(alpha, bht), 1, "can't generate simple 2mir in limiting case");
/* ensure that rho is not zero */
DGG_TEST2( rho < DGG_MIN_RHO, 1, "rho (%f) too small", rho);
/* initialize constraint */
tmir = DGG_newConstraint( base->nz );
tmir->rhs = bup*tau*rho;
tmir->sense = 'G';
/* compute cut coefficients */
for(i=0; i<base->nz; i++){
double v = base->coeff[i];
if (!isint[i]) {
if (v > 0.0) tmir->coeff[lnz] = v;
else tmir->coeff[lnz] = 0.0;
}
else {
vht = v - floor(v);
DGG_IF_EXIT( vht < 0.0, 1, "negative vht");
k = DGG_MIN(tau-1,floor(vht/alpha));
tmir->coeff[lnz] = floor(v)*tau*rho + k*rho + DGG_MIN(rho,vht-k*alpha);
}
tmir->index[lnz] = base->index[i];
lnz += 1;
}
tmir->nz = lnz;
*cut_out = tmir;
return 0;
}
/******************* TEST / DEBUGGING ROUTINES ********************************/
/* DGG_is2stepValid:
checks that:
bht > alpha > 0
(1/alpha) >= tau > (bht/alpha)
*/
int DGG_is2stepValid(double alpha, double bht)
{
/* d */
double tau;
/* ensure that alpha is not null or negative */
if ( alpha < DGG_MIN_ALPHA )
return 0;
/* compute tau and tau_lim */
tau = ceil( bht / alpha );
/* make sure alpha is not a divisor of bht */
if ( DGG_is_a_multiple_of_b(alpha, bht) )
return 0;
/* page 15, definition 12 */
/* check if alpha is admissible for simple-2-step-tmir */
if ( (bht > alpha) && (alpha > 0.0) )
if ( (1/alpha) >= tau )
return 1;
/* not admissible */
return 0;
}
/* checks that its worth doing a 1MIR on the constraint. More precisely,
- Is the RHS null?
- Are there any integer variables set at fractional values? */
int DGG_isBaseTrivial(DGG_data_t *d, DGG_constraint_t* c)
{
/* is rhs sufficiently fractional */
if ( frac_part(ABOV(c->rhs)) < d->gomory_threshold )
return 1;
if ( (1.0 - frac_part(ABOV(c->rhs))) < d->gomory_threshold )
return 1;
return 0;
}
/* tests lhs vs rhs of a constraint */
int DGG_isConstraintViolated(DGG_data_t *d, DGG_constraint_t *c)
{
double lhs = DGG_cutLHS(c, d->x);
double rhs = c->rhs;
/* compare LHS and RHS */
if (c->sense == 'G')
if ( lhs > (rhs - DGG_NULL_SLACK) )
return 0;
if (c->sense == 'L')
if ( lhs < (rhs + DGG_NULL_SLACK) )
return 0;
if (c->sense == 'E')
if ( fabs(lhs - rhs) < DGG_NULL_SLACK )
return 0;
return 0;
}
double DGG_cutLHS(DGG_constraint_t *c, double *x)
{
int i;
double lhs = 0.0;
for(i=0; i < c->nz; i++)
lhs += c->coeff[i]*x[c->index[i]];
return lhs;
}
int DGG_isCutDesirable(DGG_constraint_t *c, DGG_data_t *d)
{
double lhs, rhs;
lhs = DGG_cutLHS(c, d->x);
rhs = c->rhs;
if (c->nz > 500) return 0;
/* if the cut is not violated, return 0 */
if (c->sense == 'G')
if ( lhs > (rhs - DGG_NULL_SLACK) )
return 0;
if (c->sense == 'L')
if ( lhs < (rhs + DGG_NULL_SLACK) )
return 0;
if (c->sense == 'E')
if ( fabs(lhs - rhs) < DGG_NULL_SLACK )
return 0;
return 1;
}
/******************** SIMPLE MACROS AND FUNCTIONS *****************************/
int DGG_is_even(double vht, double bht, int tau, int q)
{
double v2 = V2I(bht, tau, q);
if ( vht > v2 )
return 1;
return 0;
}
double frac_part(double value)
{
return value-floor(value);
}
int DGG_is_a_multiple_of_b(double a, double b)
{
double c = b/a;
if ( (b - a*floor(c)) < DGG_MIN_RHO )
return 1;
return 0;
}
int DGG_cutsOffPoint(double *x, DGG_constraint_t *cut)
{
int i;
double LHS = 0.0;
for(i=0; i < cut->nz; i++)
LHS += cut->coeff[i]*(x[ cut->index[i] ]);
//fprintf(stdout, "LHS = %f, SENSE = %c, RHS = %f\n", LHS, cut->sense, cut->rhs);
if ( cut->sense == 'E' )
if ( fabs(LHS - cut->rhs) > DGG_NULL_SLACK )
goto BAD;
if (cut->sense == 'G' )
if ( (cut->rhs - LHS) > DGG_NULL_SLACK )
goto BAD;
if (cut->sense == 'L' )
if ( (LHS - cut->rhs) > DGG_NULL_SLACK )
goto BAD;
return 0;
BAD:
fprintf(stdout, "LHS = %f, SENSE = %c, RHS = %f\n", LHS, cut->sense, cut->rhs);
DGG_TEST(1, 1, "found a bad cut!");
return 0;
}
// Returns true if needs optimal basis to do cuts
bool
CglTwomir::needsOptimalBasis() const
{
return true;
}
// Away stuff
void CglTwomir::setAway(double value)
{
if (value>0.0&&value<=0.5)
away_=value;
}
double CglTwomir::getAway() const
{
return away_;
}
// Away stuff at root
void CglTwomir::setAwayAtRoot(double value)
{
if (value>0.0&&value<=0.5)
awayAtRoot_=value;
}
double CglTwomir::getAwayAtRoot() const
{
return awayAtRoot_;
}
// This can be used to refresh any information
void
CglTwomir::refreshSolver(OsiSolverInterface * solver)
{
if (originalSolver_) {
delete originalSolver_;
originalSolver_ = solver->clone();
}
}
// Create C++ lines to get to current state
std::string
CglTwomir::generateCpp( FILE * fp)
{
CglTwomir other;
fprintf(fp,"0#include \"CglTwomir.hpp\"\n");
fprintf(fp,"3 CglTwomir twomir;\n");
if (t_min_!=other.t_min_||t_max_!=other.t_max_)
fprintf(fp,"3 twomir.setMirScale(%d,%d);\n",t_min_,t_max_);
else
fprintf(fp,"4 twomir.setMirScale(%d,%d);\n",t_min_,t_max_);
if (q_min_!=other.q_min_||q_max_!=other.q_max_)
fprintf(fp,"3 twomir.setTwomirScale(%d,%d);\n",q_min_,q_max_);
else
fprintf(fp,"4 twomir.setTwomirScale(%d,%d);\n",q_min_,q_max_);
if (do_mir_!=other.do_mir_||do_2mir_!=other.do_2mir_||
do_tab_!=other.do_tab_||do_form_!=other.do_form_)
fprintf(fp,"3 twomir.setCutTypes(%s,%s,%s,%s);\n",
do_mir_ ? "true" : "false",
do_2mir_ ? "true" : "false",
do_tab_ ? "true" : "false",
do_form_ ? "true" : "false");
else
fprintf(fp,"4 twomir.setCutTypes(%s,%s,%s,%s);\n",
do_mir_ ? "true" : "false",
do_2mir_ ? "true" : "false",
do_tab_ ? "true" : "false",
do_form_ ? "true" : "false");
if (a_max_!=other.a_max_)
fprintf(fp,"3 twomir.setAMax(%d);\n",a_max_);
else
fprintf(fp,"4 twomir.setAMax(%d);\n",a_max_);
if (max_elements_!=other.max_elements_)
fprintf(fp,"3 twomir.setMaxElements(%d);\n",max_elements_);
else
fprintf(fp,"4 twomir.setMaxElements(%d);\n",max_elements_);
if (max_elements_root_!=other.max_elements_root_)
fprintf(fp,"3 twomir.setMaxElementsRoot(%d);\n",max_elements_root_);
else
fprintf(fp,"4 twomir.setMaxElementsRoot(%d);\n",max_elements_root_);
if (getAggressiveness()!=other.getAggressiveness())
fprintf(fp,"3 twomir.setAggressiveness(%d);\n",getAggressiveness());
else
fprintf(fp,"4 twomir.setAggressiveness(%d);\n",getAggressiveness());
return "twomir";
}