limp-cbc-0.3.2.0: cbits/coin/CglZeroHalf.cpp
// $Id: CglZeroHalf.cpp 1123 2013-04-06 20:47:24Z stefan $
// Copyright (C) 2010, 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 <cmath>
#include <cstdio>
#include <cfloat>
#include <cassert>
#include "CoinPragma.hpp"
#include "CglZeroHalf.hpp"
#include "CoinPackedVector.hpp"
#include "CoinSort.hpp"
#include "CoinPackedMatrix.hpp"
//-------------------------------------------------------------
void
CglZeroHalf::generateCuts(const OsiSolverInterface & si, OsiCuts & cs,
const CglTreeInfo info)
{
if (mnz_) {
int cnum=0,cnzcnt=0;
int *cbeg=NULL, *ccnt=NULL,*cind=NULL,*cval=NULL,*crhs=NULL;
char *csense=NULL;
const double * solution = si.getColSolution();
if ((flags_&1)==0) {
// redo bounds
const double * columnLower = si.getColLower();
const double * columnUpper = si.getColUpper();
int numberColumns = si.getNumCols();
for (int iColumn=0;iColumn<numberColumns;iColumn++) {
if (vlb_[iColumn]!=COIN_INT_MAX) {
int ilo,iup;
double lo = columnLower[iColumn];
if (lo<-COIN_INT_MAX)
lo=-COIN_INT_MAX;
ilo= static_cast<int> (ceil(lo));
double up = columnUpper[iColumn];
if (up>COIN_INT_MAX)
up=COIN_INT_MAX;
iup= static_cast<int> (floor(up));
vlb_[iColumn]=ilo;
vub_[iColumn]=iup;
}
}
}
if (true) {
cutInfo_.sep_012_cut(mr_,mc_,mnz_,
mtbeg_,mtcnt_, mtind_, mtval_,
vlb_, vub_,
mrhs_, msense_,
solution,
info.inTree ? false : true,
&cnum,&cnzcnt,
&cbeg,&ccnt,&cind,&cval,&crhs,&csense);
} else {
int k = 4*mr_+2*mnz_;
int * temp = new int[k];
int * mtbeg = temp;
int * mtcnt = mtbeg + mr_;
int * mtind = mtcnt+mr_;
int * mtval = mtind+mnz_;
int * mrhs = mtval+mnz_;
char * msense = reinterpret_cast<char*> (mrhs+mr_);
int i;
k=0;
int kel=0;
for (i=0;i<mr_;i++) {
int kel2=kel;
int rhs = mrhs_[i];
for (int j=mtbeg_[i];j<mtbeg_[i]+mtcnt_[i];j++) {
int iColumn=mtind_[j];
int value=mtval_[j];
if (vlb_[iColumn]<vub_[iColumn]) {
mtind[kel]=mtind_[j];
mtval[kel++]=mtval_[j];
} else {
rhs -= vlb_[iColumn]*value;
}
}
if (kel>kel2) {
mtcnt[k]=kel-kel2;
mtbeg[k]=kel2;
mrhs[k]=rhs;
msense[k++]=msense_[i];
}
}
if (kel) {
cutInfo_.sep_012_cut(k,mc_,kel,
mtbeg,mtcnt, mtind, mtval,
vlb_, vub_,
mrhs, msense,
solution,
info.inTree ? false : true,
&cnum,&cnzcnt,
&cbeg,&ccnt,&cind,&cval,&crhs,&csense);
}
delete [] temp;
}
if (cnum) {
// add cuts
double * element = new double[mc_];
for (int i=0;i<cnum;i++) {
int n = ccnt[i];
int start = cbeg[i];
for (int j=0;j<n;j++)
element[j]=cval[start+j];
OsiRowCut rc;
if (csense[i]=='L') {
rc.setLb(-COIN_DBL_MAX);
rc.setUb(crhs[i]);
} else if (csense[i]=='G') {
rc.setLb(crhs[i]);
rc.setUb(COIN_DBL_MAX);
} else {
abort();
}
rc.setRow(n,cind+start,element,false);
if ((flags_&1)!=0)
rc.setGloballyValid();
//double violation = rc.violated(solution);
//if (violation>1.0e-6)
cs.insert(rc);
//else
//printf("violation of %g\n",violation);
}
delete [] element;
free(cbeg);
free(ccnt);
free(cind);
free(cval);
free(crhs);
free(csense);
}
}
}
//-------------------------------------------------------------------
// Default Constructor
//-------------------------------------------------------------------
CglZeroHalf::CglZeroHalf ()
:
CglCutGenerator(),
mr_(0),
mc_(0),
mnz_(0),
mtbeg_(NULL),
mtcnt_(NULL),
mtind_(NULL),
mtval_(NULL),
vlb_(NULL),
vub_(NULL),
mrhs_(NULL),
msense_(NULL),
flags_(0)
{
// nothing to do here
}
//-------------------------------------------------------------------
// Copy constructor
//-------------------------------------------------------------------
CglZeroHalf::CglZeroHalf (
const CglZeroHalf & source)
:
CglCutGenerator(source),
mtbeg_(NULL),
mtcnt_(NULL),
mtind_(NULL),
mtval_(NULL),
vlb_(NULL),
vub_(NULL),
mrhs_(NULL),
msense_(NULL),
flags_(source.flags_)
{
mr_ = source.mr_;
mc_ = source.mc_;
mnz_ = source.mnz_;
if (mr_) {
mtbeg_ = CoinCopyOfArray(source.mtbeg_,mr_);
mtcnt_ = CoinCopyOfArray(source.mtcnt_,mr_);
mtind_ = CoinCopyOfArray(source.mtind_,mnz_);
mtval_ = CoinCopyOfArray(source.mtval_,mnz_);
vlb_ = CoinCopyOfArray(source.vlb_,mc_);
vub_ = CoinCopyOfArray(source.vub_,mc_);
mrhs_ = CoinCopyOfArray(source.mrhs_,mr_);
msense_ = CoinCopyOfArray(source.msense_,mr_);
}
}
//-------------------------------------------------------------------
// Clone
//-------------------------------------------------------------------
CglCutGenerator *
CglZeroHalf::clone() const
{
return new CglZeroHalf(*this);
}
//-------------------------------------------------------------------
// Destructor
//-------------------------------------------------------------------
CglZeroHalf::~CglZeroHalf ()
{
delete [] mtbeg_;
delete [] mtcnt_;
delete [] mtind_;
delete [] mtval_;
delete [] vlb_;
delete [] vub_;
delete [] mrhs_;
delete [] msense_;
}
//----------------------------------------------------------------
// Assignment operator
//-------------------------------------------------------------------
CglZeroHalf &
CglZeroHalf::operator=(
const CglZeroHalf& rhs)
{
if (this != &rhs) {
CglCutGenerator::operator=(rhs);
delete [] mtbeg_;
delete [] mtcnt_;
delete [] mtind_;
delete [] mtval_;
delete [] vlb_;
delete [] vub_;
delete [] mrhs_;
delete [] msense_;
mr_ = rhs.mr_;
mc_ = rhs.mc_;
mnz_ = rhs.mnz_;
flags_ = rhs.flags_;
cutInfo_=Cgl012Cut();
if (mr_) {
mtbeg_ = CoinCopyOfArray(rhs.mtbeg_,mr_);
mtcnt_ = CoinCopyOfArray(rhs.mtcnt_,mr_);
mtind_ = CoinCopyOfArray(rhs.mtind_,mnz_);
mtval_ = CoinCopyOfArray(rhs.mtval_,mnz_);
vlb_ = CoinCopyOfArray(rhs.vlb_,mc_);
vub_ = CoinCopyOfArray(rhs.vub_,mc_);
mrhs_ = CoinCopyOfArray(rhs.mrhs_,mr_);
msense_ = CoinCopyOfArray(rhs.msense_,mr_);
} else {
mtbeg_ = NULL;
mtcnt_ = NULL;
mtind_ = NULL;
mtval_ = NULL;
vlb_ = NULL;
vub_ = NULL;
mrhs_ = NULL;
msense_ = NULL;
}
}
return *this;
}
void
CglZeroHalf::refreshSolver(OsiSolverInterface * solver)
{
if (!solver||!solver->getNumRows())
return; // no solver
delete [] mtbeg_;
delete [] mtcnt_;
delete [] mtind_;
delete [] mtval_;
delete [] vlb_;
delete [] vub_;
delete [] mrhs_;
delete [] msense_;
mr_ = 0;
mc_ = 0;
mnz_ = 0;
mtbeg_ = NULL;
mtcnt_ = NULL;
mtind_ = NULL;
mtval_ = NULL;
vlb_ = NULL;
vub_ = NULL;
mrhs_ = NULL;
msense_ = NULL;
cutInfo_.free_log_var();
cutInfo_.free_parity_ilp();
cutInfo_.free_ilp();
CoinPackedMatrix rowCopy(*solver->getMatrixByRow());
const int * column = rowCopy.getIndices();
const CoinBigIndex * rowStart = rowCopy.getVectorStarts();
const int * rowLength = rowCopy.getVectorLengths();
const double * rowElements = rowCopy.getElements();
const double * columnLower = solver->getColLower();
const double * columnUpper = solver->getColUpper();
const double * rowLower = solver->getRowLower();
const double * rowUpper = solver->getRowUpper();
int iColumn,iRow;
// count number of possible
int numberColumns = solver->getNumCols();
int numberRows = solver->getNumRows();
vlb_ = new int [numberColumns];
vub_ = new int [numberColumns];
for (iColumn=0;iColumn<numberColumns;iColumn++) {
int ilo,iup;
if (solver->isInteger(iColumn)) {
double lo = columnLower[iColumn];
if (lo<-COIN_INT_MAX)
lo=-COIN_INT_MAX;
ilo= static_cast<int> (ceil(lo));
double up = columnUpper[iColumn];
if (up>COIN_INT_MAX)
up=COIN_INT_MAX;
iup= static_cast<int> (floor(up));
} else {
ilo=COIN_INT_MAX;
iup=-COIN_INT_MAX;
}
vlb_[iColumn]=ilo;
vub_[iColumn]=iup;
}
for (iRow=0;iRow<numberRows;iRow++) {
int n = rowLength[iRow];
bool good=(n>0);
for (CoinBigIndex j=rowStart[iRow];
j<rowStart[iRow]+n;j++) {
int jColumn = column[j];
if (vlb_[jColumn]==COIN_INT_MAX) {
// continuous
good=false;
break;
} else {
double value = rowElements[j];
if (fabs(value-floor(value+0.5))>1.0e-30) {
// not integer coefficient
good=false;
break;
}
}
}
double lo = rowLower[iRow];
double up = rowUpper[iRow];
int iType=1;
double rhs=1.0e20;
if (lo>-1.0e20) {
if (fabs(lo-floor(lo+0.5))>1.0e-12) {
// not integer coefficient
good=false;
}
rhs=fabs(lo);
if (up<1.0e20) {
rhs=CoinMax(fabs(lo),fabs(up));
if (lo!=up)
iType=2; // ranged so make copy
if (fabs(up-floor(up+0.5))>1.0e-12) {
// not integer coefficient
good=false;
}
}
} else if (up<1.0e20) {
rhs=fabs(up);
if (up<1.0e20) {
if (fabs(up-floor(up+0.5))>1.0e-12) {
// not integer coefficient
good=false;
}
}
}
if (good&&rhs<COIN_INT_MAX) {
mr_+=iType;
mnz_ += iType*n;
}
}
if (mnz_) {
mc_ = numberColumns;
mtbeg_ = new int [mr_];
mtcnt_ = new int [mr_];
mtind_ = new int [mnz_];
mtval_ = new int [mnz_];
mrhs_ = new int [mr_];
msense_ = new char [mr_];
mr_=0;
mnz_=0;
for (iRow=0;iRow<numberRows;iRow++) {
int n = rowLength[iRow];
bool good=(n>0);
for (CoinBigIndex j=rowStart[iRow];
j<rowStart[iRow]+n;j++) {
int jColumn = column[j];
if (vlb_[jColumn]==COIN_INT_MAX) {
// continuous
good=false;
break;
} else {
double value = rowElements[j];
if (fabs(value-floor(value+0.5))>1.0e-12) {
// not integer coefficient
good=false;
break;
}
}
}
double lo = rowLower[iRow];
double up = rowUpper[iRow];
int iType=1;
double rhs=1.0e20;
if (lo>-1.0e20) {
if (fabs(lo-floor(lo+0.5))>1.0e-30) {
// not integer coefficient
good=false;
}
rhs=fabs(lo);
if (up<1.0e20) {
rhs=CoinMax(fabs(lo),fabs(up));
if (lo!=up)
iType=2; // ranged so make copy
if (fabs(up-floor(up+0.5))>1.0e-12) {
// not integer coefficient
good=false;
}
}
} else if (up<1.0e20) {
rhs=fabs(up);
if (up<1.0e20) {
if (fabs(up-floor(up+0.5))>1.0e-12) {
// not integer coefficient
good=false;
}
}
}
if (good&&rhs<COIN_INT_MAX) {
mtbeg_[mr_]=mnz_;
for (CoinBigIndex j=rowStart[iRow];
j<rowStart[iRow]+n;j++) {
int jColumn = column[j];
double value = rowElements[j];
assert (fabs(value)<COIN_INT_MAX);
int iValue = static_cast<int> (floor(value+0.5));
if (iValue) {
mtind_[mnz_]=jColumn;
mtval_[mnz_++]=iValue;
}
}
mtcnt_[mr_]=mnz_-mtbeg_[mr_];
if (iType==1) {
if (lo>-1.0e20) {
mrhs_[mr_]=static_cast<int> (floor(lo+0.5));
msense_[mr_]='G';
} else {
mrhs_[mr_]=static_cast<int> (floor(up+0.5));
msense_[mr_]='L';
}
mr_++;
} else {
// ranged!
mrhs_[mr_]=static_cast<int> (floor(lo+0.5));
msense_[mr_]='G';
int k = mnz_-mtbeg_[mr_];
mr_++;
mtbeg_[mr_]=mnz_;
memcpy(mtind_+mnz_,mtind_+mnz_-k,k*sizeof(int));
memcpy(mtval_+mnz_,mtval_+mnz_-k,k*sizeof(int));
mnz_+= mtcnt_[mr_-1];
mtcnt_[mr_]=mnz_-mtbeg_[mr_];
mrhs_[mr_]=static_cast<int> (floor(up+0.5));
msense_[mr_]='L';
mr_++;
}
}
}
cutInfo_.ilp_load(mr_,mc_,mnz_,mtbeg_,mtcnt_,mtind_,mtval_,
vlb_,vub_,mrhs_,msense_);
cutInfo_.alloc_parity_ilp(mr_,mc_,mnz_);
cutInfo_.initialize_log_var();
} else {
// no good
delete [] vlb_;
delete [] vub_;
vlb_ = NULL;
vub_ = NULL;
mr_=0;
mnz_=0;
}
}
// Create C++ lines to get to current state
std::string
CglZeroHalf::generateCpp( FILE * fp)
{
CglZeroHalf other;
fprintf(fp,"0#include \"CglZeroHalf.hpp\"\n");
fprintf(fp,"3 CglZeroHalf zeroHalf;\n");
if (getAggressiveness()!=other.getAggressiveness())
fprintf(fp,"3 zeroHalf.setAggressiveness(%d);\n",getAggressiveness());
else
fprintf(fp,"4 zeroHalf.setAggressiveness(%d);\n",getAggressiveness());
return "zeroHalf";
}
#include <vector>
#include <algorithm>
//bool operator() (cgl_node * x, cgl_node * y) {
bool best(cgl_node * x, cgl_node * y) {
return (x->distanceBack>y->distanceBack);
}
#ifndef CGL_NEW_SHORT
void cglShortestPath(cgl_graph * graph, int source, int maximumLength)
#else
void cglShortestPath(auxiliary_graph * graph, int source, int maximumLength)
#endif
{
int numberNodes=graph->nnodes;
#define HEAP
#ifndef HEAP
int * candidate = new int [numberNodes];
#endif
cgl_node * nodes = graph->nodes;
int i;
for ( i=0;i<numberNodes;i++) {
nodes[i].parentNode=-1;
nodes[i].distanceBack=COIN_INT_MAX;
#ifndef HEAP
candidate[i]=i;
#endif
}
nodes[source].distanceBack=0;
#ifdef HEAP
std::vector <cgl_node *> nodes_;
for ( i=0;i<numberNodes;i++)
nodes_.push_back(nodes+i);
// create heap
std::make_heap(nodes_.begin(), nodes_.end(), best);
#endif
int numberCandidates = numberNodes;
while (numberCandidates>0) {
#ifdef HEAP
cgl_node * bestNode = nodes_.front();
int iNode = bestNode->index;
std::pop_heap(nodes_.begin(), nodes_.end(), best);
nodes_.pop_back();
if (nodes[iNode].distanceBack==COIN_INT_MAX)
break;
numberCandidates--;
#else
int best=-1;
int bestDistance=COIN_INT_MAX;
for (i=0;i<numberCandidates;i++) {
int iNode = candidate[i];
if (nodes[iNode].distanceBack<bestDistance) {
best=i;
bestDistance=nodes[iNode].distanceBack;
}
}
if (best<0)
break; // disconnected graph?
int iNode = candidate[best];
numberCandidates--;
candidate[best]=candidate[numberCandidates];
#endif
cgl_arc * nextNode = nodes[iNode+1].firstArc;
int thisDistance = nodes[iNode].distanceBack;
for (cgl_arc * arc=nodes[iNode].firstArc;arc!=nextNode;arc++) {
int toNode = arc->to;
int dist = arc->length;
if (thisDistance+dist<nodes[toNode].distanceBack) {
nodes[toNode].distanceBack=thisDistance+dist;
nodes[toNode].parentNode=iNode;
#ifdef HEAP
nodes_.push_back(nodes+toNode);
//printf("size %d\n",nodes_.size());
#endif
}
}
}
}