limp-cbc-0.3.2.0: cbits/coin/Cgl012cut.cpp
// $Id: Cgl012cut.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).
/** @file 012cut.c Definition file for C coded 0-1/2 separator */
#include "CoinFinite.hpp"
#include "CoinTime.hpp"
#include "Cgl012cut.hpp"
#include "CglZeroHalf.hpp"
int MAX_CUTS = 1000;
//#define PRINT_TABU
//#define PRINT_CUTS
//#define PRINT_TIME
//#define TIME
/* #define TIME */
#undef TIME
#define TRUE 1
#define FALSE 0
#define ODD 1
#define EVEN 0
#define NONE -1
#define BOTH 2
#define IN 1
#define OUT 0
#define ADD 1
#define DEL 0
#define LOWER_BOUND 0
#define UPPER_BOUND 1
#define EPS 0.0001 /* small tolerance */
//#define EPS 0.000001 /* small tolerance */
#define ZERO 0.000001 /* estimated accuracy for doubles */
//#define ZERO 0.0001 /* estimated accuracy for doubles */
#define INF 1000000000.0
#define IINF 1000000000
#define MAX_SLACK 1.0
#define MAX_LOSS 1.0
#define MAX_CYCLE_WEIGHT 1.0
#define MIN_VIOLATION 0.001
#define MIN_SCORE_RANGE 10.0
#define MAX_SCORE_RANGE ZERO /* 1.0 */
#define ISCALE 10000
//#define MAX_CUTS 10
#define MAX_CUT_POOL 10000
#define MAX_CUT_COD 10000
#define MAX_ITER_POOL 100
#define CLEAN_THRESH 0.9
#define MANY_IT_ZERO 10
#define mod2(I) ( I % 2 == 0 ? 0 : 1 )
#ifdef TIME
static float tot_basic_sep_time = 0.0; /* total time spent for basic
separation */
static float avg_basic_sep_time; /* average time per iteration spent for basic
separation */
static float total_time = 0.0; /* total time spent in the separation */
static float prep_time = 0.0; /* time spent for the definition of the
parity ILP data structure */
static float weak_time = 0.0; /* time spent for the construction of the
separation graph by weakening */
static float aux_time = 0.0; /* time spent for the definition of the
auxiliary graph */
static float path_time = 0.0; /* time spent in the computation of the
shortest paths */
static float cycle_time = 0.0; /* time spent in the determination of the
shortest cycles */
static float cut_time = 0.0; /* time spent in the determination of the
violated cuts */
static float bw_time = 0.0; /* time spent in best_weakening */
static float coef_time = 0.0; /* time spent in the initial computation
of coef in get_cut */
static float pool_time = 0.0; /* time spent for the addition and
extraction of cuts from the pool */
static int cut_ncalls = 0; /* number of calls to get_cut */
static float tabu_time = 0.0; /* time spent within tabu search */
float ti, tf, td, tti, ttf, tsi, tsf, tii, tff, tpi, tpf, ttabi, ttabf;
void
second_(float *t) {*t=CoinCpuTime();}
#endif
/* #endif */
/* global data structures */
#define CGGGGG
#ifndef CGGGGG
static ilp *inp_ilp; /* input ILP data structure */
static parity_ilp *p_ilp; /* parity ILP data structure */
#endif
#ifdef PRINT_CUTS
/* utility subroutines */
void print_int_vect(char *s,int *v,int n)
{
int i;
printf("integer vector %s:",s);
for ( i = 0; i < n; i++ ) printf(" %d",v[i]);
printf("\n");
}
void print_short_int_vect(char *s,short int *v,int n)
{
int i;
printf("short integer vector %s:",s);
for ( i = 0; i < n; i++ ) printf(" %d",v[i]);
printf("\n");
}
void print_double_vect(char *s,double *v,int n)
{
int i;
printf("double vector %s:",s);
for ( i = 0; i < n; i++ ) printf(" %f",v[i]);
printf("\n");
}
#endif
void alloc_error(char *s)
{
printf("\n Warning: Not enough memory to allocate %s\n",s);
printf("\n Cannot proceed with 0-1/2 cut separation\n");
exit(FALSE);
}
/* double2int: compute the integer equivalent of a double */
int double2int(double x)
{
if ( x > IINF ) return (IINF);
if ( x < - IINF ) return (- IINF);
if ( x < ZERO && x > - ZERO ) return(0);
if ( x > 0.0 ) return(static_cast<int> (x + ZERO));
return(static_cast<int> (x - ZERO));
}
/* gcd: compute the greatest common divisor of two integers */
int gcd(int a,int b)
{
int c;
if ( a < 0 ) a = - a;
if ( b < 0 ) b = - b;
if ( a < b ) { c = a; a = b; b = c; }
while ( b != 0 ) {
c = a % b; a = b; b = c;
}
return(a);
}
/* ILP data structures subroutines */
/* ilp_load: load the input ILP into an internal data structure */
void Cgl012Cut::ilp_load(
int mr, /* number of rows in the ILP matrix */
int mc, /* number of columns in the ILP matrix */
int mnz, /* number of nonzero's in the ILP matrix */
int *mtbeg, /* starting position of each row in arrays mtind and mtval */
int *mtcnt, /* number of entries of each row in arrays mtind and mtval */
int *mtind, /* column indices of the nonzero entries of the ILP matrix */
int *mtval, /* values of the nonzero entries of the ILP matrix */
int *vlb, /* lower bounds on the variables */
int *vub, /* upper bounds on the variables */
int *mrhs, /* right hand sides of the constraints */
char *msense /* senses of the constraints: 'L', 'G' or 'E' */
)
{
inp_ilp = reinterpret_cast<ilp *> (calloc(1,sizeof(ilp)));
if ( inp_ilp == NULL ) alloc_error(const_cast<char*>("inp_ilp"));
inp_ilp->mr = mr; inp_ilp->mc = mc; inp_ilp->mnz = mnz;
inp_ilp->mtbeg = mtbeg; inp_ilp->mtcnt = mtcnt;
inp_ilp->mtind = mtind; inp_ilp->mtval = mtval;
inp_ilp->vlb = vlb; inp_ilp->vub = vub;
inp_ilp->mrhs = mrhs; inp_ilp->msense = msense;
}
void Cgl012Cut::free_ilp()
{
free(inp_ilp);
inp_ilp=NULL;
}
#ifdef PRINT_CUTS
void Cgl012Cut::print_constr(int i /* constraint to be printed */)
{
printf("\n content of constraint %d: nzcnt = %d, rhs = %d, sense = %c, slack = %f\n",
i, inp_ilp->mtcnt[i], inp_ilp->mrhs[i], inp_ilp->msense[i], p_ilp->slack[i]);
print_int_vect(const_cast<char*>("ind"),inp_ilp->mtind + inp_ilp->mtbeg[i],inp_ilp->mtcnt[i]);
print_int_vect(const_cast<char*>("val"),inp_ilp->mtval + inp_ilp->mtbeg[i],inp_ilp->mtcnt[i]);
}
#endif
/* alloc_parity_ilp: allocate the memory for the parity ILP data structure */
void Cgl012Cut::alloc_parity_ilp(
int mr, /* number of rows in the ILP matrix */
int mc, /* number of columns in the ILP matrix */
int mnz /* number of nonzero's in the ILP matrix */
)
{
p_ilp = reinterpret_cast<parity_ilp *> (calloc(1,sizeof(parity_ilp)));
if ( p_ilp == NULL ) alloc_error(const_cast<char*>("p_ilp"));
p_ilp->mtbeg = reinterpret_cast<int *> (calloc(mr,sizeof(int)));
if ( p_ilp->mtbeg == NULL ) alloc_error(const_cast<char*>("p_ilp->mtbeg"));
p_ilp->mtcnt = reinterpret_cast<int *> (calloc(mr,sizeof(int)));
if ( p_ilp->mtcnt == NULL ) alloc_error(const_cast<char*>("p_ilp->mtcnt"));
p_ilp->mtind = reinterpret_cast<int *> (calloc(mnz,sizeof(int)));
if ( p_ilp->mtind == NULL ) alloc_error(const_cast<char*>("p_ilp->mtind"));
p_ilp->mrhs = reinterpret_cast<short int *> (calloc(mr,sizeof(short int)));
if ( p_ilp->mrhs== NULL ) alloc_error(const_cast<char*>("p_ilp->mrhs"));
p_ilp->xstar = reinterpret_cast<double *> (calloc(mc,sizeof(double)));
if ( p_ilp->xstar== NULL ) alloc_error(const_cast<char*>("p_ilp->xstar"));
p_ilp->slack = reinterpret_cast<double *> (calloc(mr,sizeof(double)));
if ( p_ilp->slack == NULL ) alloc_error(const_cast<char*>("p_ilp->slack"));
p_ilp->row_to_delete = reinterpret_cast<short int *> (calloc(mr,sizeof(short int)));
if ( p_ilp->row_to_delete == NULL ) alloc_error(const_cast<char*>("p_ilp->row_to_delete"));
p_ilp->col_to_delete = reinterpret_cast<short int *> (calloc(mc,sizeof(short int)));
if ( p_ilp->col_to_delete == NULL ) alloc_error(const_cast<char*>("p_ilp->col_to_delete"));
p_ilp->gcd = reinterpret_cast<int *> (calloc(mr,sizeof(int)));
if ( p_ilp->gcd == NULL ) alloc_error(const_cast<char*>("p_ilp->gcd"));
p_ilp->possible_weak = reinterpret_cast<short int *> (calloc(mc,sizeof(short int)));
if ( p_ilp->possible_weak == NULL ) alloc_error(const_cast<char*>("p_ilp->possible_weak"));
p_ilp->type_even_weak = reinterpret_cast<short int *> (calloc(mc,sizeof(short int)));
if ( p_ilp->type_even_weak == NULL ) alloc_error(const_cast<char*>("p_ilp->type_even_weak"));
p_ilp->type_odd_weak = reinterpret_cast<short int *> (calloc(mc,sizeof(short int)));
if ( p_ilp->type_odd_weak == NULL ) alloc_error(const_cast<char*>("p_ilp->type_odd_weak"));
p_ilp->loss_even_weak = reinterpret_cast<double *> (calloc(mc,sizeof(double)));
if ( p_ilp->loss_even_weak == NULL ) alloc_error(const_cast<char*>("p_ilp->loss_even_weak"));
p_ilp->loss_odd_weak = reinterpret_cast<double *> (calloc(mc,sizeof(double)));
if ( p_ilp->loss_odd_weak == NULL ) alloc_error(const_cast<char*>("p_ilp->loss_odd_weak"));
p_ilp->min_loss_by_weak = reinterpret_cast<double *> (calloc(mc,sizeof(double)));
if ( p_ilp->min_loss_by_weak == NULL ) alloc_error(const_cast<char*>("p_ilp->min_loss_by_weak"));
p_ilp->mr=mr;
p_ilp->mc=mc;
p_ilp->mnz=mnz;
}
#ifdef PRINT_CUTS
void Cgl012Cut::print_parity_ilp()
{
printf("\n content of parity_ilp data structure: mc = %d, mr = %d, mnz = %d\n",
p_ilp->mc,p_ilp->mr,p_ilp->mnz);
print_int_vect(const_cast<char*>("mtbeg"),p_ilp->mtbeg,p_ilp->mr);
print_int_vect(const_cast<char*>("mtcnt"),p_ilp->mtcnt,p_ilp->mr);
print_int_vect(const_cast<char*>("mtind"),p_ilp->mtind,p_ilp->mnz);
print_short_int_vect(const_cast<char*>("mrhs"),p_ilp->mrhs,p_ilp->mr);
print_double_vect(const_cast<char*>("xstar"),p_ilp->xstar,p_ilp->mc);
print_double_vect(const_cast<char*>("slack"),p_ilp->slack,p_ilp->mr);
print_short_int_vect(const_cast<char*>("row_to_delete"),p_ilp->row_to_delete,p_ilp->mr);
print_short_int_vect(const_cast<char*>("col_to_delete"),p_ilp->col_to_delete,p_ilp->mc);
print_int_vect(const_cast<char*>("gcd"),p_ilp->gcd,p_ilp->mr);
print_short_int_vect(const_cast<char*>("possible_weak"),p_ilp->possible_weak,p_ilp->mc);
print_short_int_vect(const_cast<char*>("type_even_weak"),p_ilp->type_even_weak,p_ilp->mc);
print_short_int_vect(const_cast<char*>("type_odd_weak"),p_ilp->type_odd_weak,p_ilp->mc);
print_double_vect(const_cast<char*>("loss_even_weak"),p_ilp->loss_even_weak,p_ilp->mc);
print_double_vect(const_cast<char*>("loss_odd_weak"),p_ilp->loss_odd_weak,p_ilp->mc);
print_double_vect(const_cast<char*>("min_loss_by_weak"),p_ilp->min_loss_by_weak,p_ilp->mc);
}
#endif
void Cgl012Cut::free_parity_ilp()
{
if (p_ilp) {
free(p_ilp->mtbeg);
free(p_ilp->mtcnt);
free(p_ilp->mtind);
free(p_ilp->mrhs);
free(p_ilp->xstar);
free(p_ilp->slack);
free(p_ilp->row_to_delete);
free(p_ilp->col_to_delete);
free(p_ilp->gcd);
free(p_ilp->possible_weak);
free(p_ilp->type_even_weak);
free(p_ilp->type_odd_weak);
free(p_ilp->loss_even_weak);
free(p_ilp->loss_odd_weak);
free(p_ilp->min_loss_by_weak);
free(p_ilp);
p_ilp=NULL;
}
}
/* alloc_info_weak: allocate memory for the weakening info data structure */
info_weak *alloc_info_weak(int nweak /* number of variables to be weakened */)
{
info_weak *i_weak;
i_weak = reinterpret_cast<info_weak *> (calloc(1,sizeof(info_weak)));
if ( i_weak == NULL ) alloc_error(const_cast<char*>("i_weak"));
if ( nweak > 0 ) {
i_weak->var = reinterpret_cast<int *> (calloc(nweak,sizeof(int)));
if ( i_weak->var == NULL ) alloc_error(const_cast<char*>("i_weak->var"));
i_weak->type = reinterpret_cast<short int *> (calloc(nweak,sizeof(short int)));
if ( i_weak->type == NULL ) alloc_error(const_cast<char*>("i_weak->type"));
}
return(i_weak);
}
#ifdef PRINT_CUTS
void print_info_weak(info_weak *i_weak)
{
printf("\n content of info_weak: nweak = %d\n",i_weak->nweak);
if ( i_weak->nweak > 0 ) {
print_int_vect(const_cast<char*>("var"),i_weak->var,i_weak->nweak);
print_short_int_vect(const_cast<char*>("type"),i_weak->type,i_weak->nweak);
}
}
#endif
void free_info_weak(info_weak *i_weak)
{
if ( i_weak->nweak > 0 ) {
free(i_weak->var);
free(i_weak->type);
}
free(i_weak);
}
/* get_parity_ilp: construct an internal data structure containing all the
information which can be useful for 0-1/2 cut separation */
void Cgl012Cut::get_parity_ilp()
{
int i, j, h, ij, aij, cnti, cnttot, begi, begh, ofsj, gcdi, ubj, lbj;
double slacki, xstarj, loss_upper, loss_lower;
short int parity_col_removed, equalih;
/* allocate the memory for the parity ILP data structure */
//alloc_parity_ilp(inp_ilp->mr,inp_ilp->mc,inp_ilp->mnz);
p_ilp->mr = inp_ilp->mr;
p_ilp->mc = inp_ilp->mc;
/* mark the variables equal to their lower/upper bound */
parity_col_removed = 0;
for ( j = 0; j < inp_ilp->mc; j++ ) {
xstarj = p_ilp->xstar[j] = inp_ilp->xstar[j];
ubj = inp_ilp->vub[j]; lbj = inp_ilp->vlb[j];
if ( xstarj > static_cast<double> (ubj - ZERO) ) {
/* variable at its upper bound */
p_ilp->col_to_delete[j] = TRUE;
if ( mod2(ubj) == ODD ) {
p_ilp->possible_weak[j] = ODD;
p_ilp->type_odd_weak[j] = UPPER_BOUND;
p_ilp->loss_odd_weak[j] = 0.0;
if ( parity_col_removed == EVEN ) parity_col_removed = ODD;
else parity_col_removed = EVEN;
}
else {
p_ilp->possible_weak[j] = EVEN;
p_ilp->type_even_weak[j] = UPPER_BOUND;
p_ilp->loss_even_weak[j] = 0.0;
}
p_ilp->min_loss_by_weak[j] = 0.0;
}
else if ( xstarj < static_cast<double> (lbj) + ZERO ) {
/* variable at its lower bound */
p_ilp->col_to_delete[j] = TRUE;
if ( mod2(lbj) == ODD ) {
p_ilp->possible_weak[j] = ODD;
p_ilp->type_odd_weak[j] = LOWER_BOUND;
p_ilp->loss_odd_weak[j] = 0.0;
p_ilp->min_loss_by_weak[j] = 0.0;
if ( parity_col_removed == EVEN ) parity_col_removed = ODD;
else parity_col_removed = EVEN;
}
else {
p_ilp->possible_weak[j] = EVEN;
p_ilp->type_even_weak[j] = LOWER_BOUND;
p_ilp->loss_even_weak[j] = 0.0;
p_ilp->min_loss_by_weak[j] = 0.0;
}
p_ilp->min_loss_by_weak[j] = 0.0;
}
else {
/* variable neither at its lower nor at its upper bound */
p_ilp->col_to_delete[j] = FALSE;
loss_upper = static_cast<double> (ubj) - xstarj;
loss_lower = xstarj - static_cast<double> (lbj);
if ( ( loss_upper > MAX_LOSS ) && ( loss_lower > MAX_LOSS ) )
/* no weakening for the variable */
p_ilp->possible_weak[j] = NONE;
else if ( loss_upper > MAX_LOSS ) {
/* lower weakening only */
if ( mod2(lbj) == EVEN ) {
p_ilp->possible_weak[j] = EVEN;
p_ilp->type_even_weak[j] = LOWER_BOUND;
p_ilp->loss_even_weak[j] = loss_lower;
}
else {
p_ilp->possible_weak[j] = ODD;
p_ilp->type_odd_weak[j] = LOWER_BOUND;
p_ilp->loss_odd_weak[j] = loss_lower;
}
}
else if ( loss_lower > MAX_LOSS ) {
/* upper weakening only */
if ( mod2(ubj) == EVEN ) {
p_ilp->possible_weak[j] = EVEN;
p_ilp->type_even_weak[j] = UPPER_BOUND;
p_ilp->loss_even_weak[j] = loss_upper;
}
else {
p_ilp->possible_weak[j] = ODD;
p_ilp->type_odd_weak[j] = UPPER_BOUND;
p_ilp->loss_odd_weak[j] = loss_upper;
}
}
else if ( mod2(ubj) == mod2(lbj) ) {
/* lower and upper bound have the same parity:
choose the best weakening */
if ( mod2(ubj) == EVEN ) {
p_ilp->possible_weak[j] = EVEN;
if ( loss_lower <= loss_upper ) {
p_ilp->type_even_weak[j] = LOWER_BOUND;
p_ilp->loss_even_weak[j] = loss_lower;
}
else {
p_ilp->type_even_weak[j] = UPPER_BOUND;
p_ilp->loss_even_weak[j] = loss_upper;
}
}
else {
p_ilp->possible_weak[j] = ODD;
if ( loss_lower <= loss_upper ) {
p_ilp->type_odd_weak[j] = LOWER_BOUND;
p_ilp->loss_odd_weak[j] = loss_lower;
}
else {
p_ilp->type_odd_weak[j] = UPPER_BOUND;
p_ilp->loss_odd_weak[j] = loss_upper;
}
}
}
else {
/* lower and upper bound have different parities:
consider both weakenings */
p_ilp->possible_weak[j] = BOTH;
if ( mod2(ubj) == EVEN ) {
p_ilp->type_even_weak[j] = UPPER_BOUND;
p_ilp->loss_even_weak[j] = loss_upper;
p_ilp->type_odd_weak[j] = LOWER_BOUND;
p_ilp->loss_odd_weak[j] = loss_lower;
}
else {
p_ilp->type_even_weak[j] = LOWER_BOUND;
p_ilp->loss_even_weak[j] = loss_lower;
p_ilp->type_odd_weak[j] = UPPER_BOUND;
p_ilp->loss_odd_weak[j] = loss_upper;
}
}
if ( loss_upper > loss_lower ) p_ilp->min_loss_by_weak[j] = loss_lower;
else p_ilp->min_loss_by_weak[j] = loss_upper;
}
}
/* scan the constraints and delete those which are trivially useless
in the 0-1/2 cut separation */
cnttot = 0;
for ( i = 0; i < inp_ilp->mr; i++ ) {
begi = inp_ilp->mtbeg[i];
/* compute the row slack and the GCD of the entries of the row */
slacki = static_cast<double> (inp_ilp->mrhs[i]);
gcdi = inp_ilp->mrhs[i];
for ( ofsj = 0; ofsj < inp_ilp->mtcnt[i]; ofsj++ ) {
ij = begi + ofsj;
j = inp_ilp->mtind[ij];
aij = inp_ilp->mtval[ij];
slacki -= static_cast<double> (aij ) * ( inp_ilp->xstar[j] );
gcdi = gcd(gcdi,aij);
}
if ( inp_ilp->msense[i] == 'G' ) slacki = -slacki;
if ( slacki < -ZERO || ( inp_ilp->msense[i] == 'E' && slacki > ZERO ) ) {
#ifdef COIN_DEVELOP
printf("\n Warning: constraint %d in the model is violated:\n",i);
printf("\n 0-1/2 cut separation is not possible\n");
printf("\nnumber of nonzero's %d\n",inp_ilp->mtcnt[i]);
printf("nonzero's (col,coef,xstar) ");
for (ofsj=0;ofsj<inp_ilp->mtcnt[i];ofsj++) printf("(%d,%d,%f) ",
inp_ilp->mtind[begi+ofsj], inp_ilp->mtval[begi+ofsj],
inp_ilp->xstar[inp_ilp->mtind[begi+ofsj]]);
printf("\n");
printf("sense %c and rhs %d and slack %.5e\n",inp_ilp->msense[i],inp_ilp->mrhs[i], slacki);
#endif
//exit(0);
slacki = INF;
}
p_ilp->slack[i] = slacki;
/* mark the rows with slack greater than the maximum allowed */
if ( slacki > MAX_SLACK - EPS ) p_ilp->row_to_delete[i] = TRUE;
else p_ilp->row_to_delete[i] = FALSE;
/* store the odd entries in the (possibly scaled) row i */
//if ( gcdi != 1 )
//printf("Warning: constraint %d with nonprime coefficients\n",i);
p_ilp->gcd[i] = gcdi;
p_ilp->mrhs[i] = mod2(( inp_ilp->mrhs[i] / gcdi ));
p_ilp->mtbeg[i] = cnttot;
cnti = 0;
for ( ofsj = 0; ofsj < inp_ilp->mtcnt[i]; ofsj++ ) {
ij = begi + ofsj;
j = inp_ilp->mtind[ij];
aij = mod2(( inp_ilp->mtval[ij] / gcdi ));
if ( aij == ODD ) {
if ( ! p_ilp->col_to_delete[j] ) {
p_ilp->mtind[cnttot] = j;
cnti++; cnttot++;
}
else if ( p_ilp->possible_weak[j] == ODD ) {
if ( p_ilp->mrhs[i] == EVEN ) p_ilp->mrhs[i] = ODD;
else p_ilp->mrhs[i] = EVEN;
}
}
}
p_ilp->mtcnt[i] = cnti;
if ( cnti == 0 ) /* (scaled) row with even entries only */
p_ilp->row_to_delete[i] = TRUE;
else {
if ( cnti == 1 && slacki < EPS ) { /* the row could be deleted */
#ifdef PRINT
printf("get_parity_ilp: row %d could be deleted since it\n",i);
printf("has only one odd entry and is tight, but it is not ...\n");
#endif
}
}
}
p_ilp->mnz = cnttot;
#ifdef REDUCTION
/* remove identical rows in the parity matrix */
/* very trivial implementation */
for ( i = 0; i < p_ilp->mr; i++ )
for ( h = i+1; h < p_ilp->mr; h++ )
if ( ( p_ilp->mrhs[i] == p_ilp->mrhs[h] ) &&
( p_ilp->mtcnt[i] == p_ilp->mtcnt[h] ) &&
( ! p_ilp->row_to_delete[i] ) &&
( ! p_ilp->row_to_delete[h] ) ) {
begi = p_ilp->mtbeg[i]; begh = p_ilp->mtbeg[h];
equalih = TRUE;
for ( ofsj = 0; ofsj < p_ilp->mtcnt[i]; ofsj++ )
/* the check assumes the indexes of the columns associated
with each row are ordered in p_ilp->mtind[] ... */
if ( p_ilp->mtind[begi+ofsj] != p_ilp->mtind[begh+ofsj] ) {
equalih = FALSE;
break;
}
if ( equalih ) {
if ( p_ilp->slack[h] > p_ilp->slack[i] )
p_ilp->row_to_delete[h] = TRUE;
else
p_ilp->row_to_delete[i] = TRUE;
}
}
/* check for the existence of separate connected components in the
parity matrix row intersection graph */
/* not implemented so far - if ever, the availability of the parity
matrix in column form also would be really convenient */
#endif
}
/* separation graph subroutines */
#define SG_EDGE_INDEX(s_graph,J,K) ( ((J) < (K)) ? ( (s_graph->nnodes * (J)) - (((J)+1)*(J)/2) + (K) - (J) -1 ) : ( (s_graph->nnodes * (K)) - (((K)+1)*(K)/2) + (J) - (K) -1 ) )
/* initialize_sep_graph: allocate and initialize the data structure
to contain the information associated with a separation graph */
separation_graph *Cgl012Cut::initialize_sep_graph()
{
int maxnodes, maxedges, nnodes, j, jk;
int *nodes, *ind;
separation_graph *s_graph;
s_graph = reinterpret_cast<separation_graph *> (calloc(1,sizeof(separation_graph)));
if ( s_graph == NULL ) alloc_error(const_cast<char*>("s_graph"));
maxnodes = p_ilp->mc + 1;
nnodes = 0;
nodes = reinterpret_cast<int *> (calloc(maxnodes,sizeof(int)));
if ( nodes == NULL ) alloc_error(const_cast<char*>("nodes"));
ind = reinterpret_cast<int *> (calloc(maxnodes,sizeof(int)));
if ( ind == NULL ) alloc_error(const_cast<char*>("ind"));
for ( j = 0; j < p_ilp->mc; j++ )
if ( ! p_ilp->col_to_delete[j] ) {
/* variable not removed from the separation problem */
nodes[nnodes] = j;
ind[j] = nnodes;
nnodes++;
}
/* take into account the special node */
nodes[nnodes] = maxnodes - 1;
ind[maxnodes-1] = nnodes;
nnodes++;
s_graph->nnodes = nnodes;
s_graph->nedges = 0;
s_graph->nodes = reinterpret_cast<int *> (malloc(nnodes*sizeof(int)));
if ( s_graph->nodes == NULL ) alloc_error(const_cast<char*>("s_graph->nodes"));
for ( j = 0; j < nnodes; j++ ) s_graph->nodes[j] = nodes[j];
free(nodes);
s_graph->ind = reinterpret_cast<int *> (malloc(maxnodes*sizeof(int)));
if ( s_graph->ind == NULL ) alloc_error(const_cast<char*>("s_graph->ind"));
for ( j = 0; j < maxnodes; j++ ) s_graph->ind[j] = ind[j];
free(ind);
maxedges = (nnodes * (nnodes - 1)) / 2;
s_graph->even_adj_list = reinterpret_cast<edge **> (malloc(maxedges*sizeof(edge *)));
if ( s_graph->even_adj_list == NULL ) alloc_error(const_cast<char*>("s_graph->even_adj_list"));
s_graph->odd_adj_list = reinterpret_cast<edge **> (malloc(maxedges*sizeof(edge *)));
if ( s_graph->odd_adj_list == NULL ) alloc_error(const_cast<char*>("s_graph->odd_adj_list"));
for ( jk = 0; jk < maxedges; jk++ )
s_graph->even_adj_list[jk] = s_graph->odd_adj_list[jk] = NULL;
return(s_graph);
}
/* update_weight_sep_graph: consider a new edge obtained from the
(weakened) parity ILP and (possibly) add it to the separation graph */
separation_graph *update_weight_sep_graph(
int j, int k, /* endpoints of the new edge */
double weight, /* weight of the new edge */
short int parity, /* parity of the new edge */
int i, /* constraint associated with the new edge */
info_weak *i_weak, /* information associated with the weakening */
separation_graph *s_graph /* separation graph to be updated */
)
{
int indj, indk, indjk;
edge *old_edge, *new_edge;
indj = s_graph->ind[j]; indk = s_graph->ind[k];
indjk = SG_EDGE_INDEX(s_graph,indj,indk);
if ( parity == EVEN ) old_edge = s_graph->even_adj_list[indjk];
else old_edge = s_graph->odd_adj_list[indjk];
if ( old_edge == NULL ) {
/* edge is not in the graph */
new_edge = reinterpret_cast<edge *> (calloc(1,sizeof(edge)));
if ( new_edge == NULL ) alloc_error(const_cast<char*>("new_edge"));
new_edge->endpoint1 = indj; new_edge->endpoint2 = indk;
new_edge->weight = weight; new_edge->parity = parity;
new_edge->constr = i; new_edge->weak = i_weak;
(s_graph->nedges)++;
if ( parity == EVEN ) s_graph->even_adj_list[indjk] = new_edge;
else s_graph->odd_adj_list[indjk] = new_edge;
}
else {
/* edge is already in the graph */
if ( old_edge->weight > weight ) {
/* replace the old edge */
old_edge->weight = weight; old_edge->constr = i;
free_info_weak(old_edge->weak);
old_edge->weak = i_weak;
}
else {
/* keep the old edge */
free_info_weak(i_weak);
}
}
return(s_graph);
}
#ifdef PRINT_CUTS
void print_edge(edge *e)
{
printf("\n content of edge: endpoint1 = %d, endpoint2 = %d, weight = %f, parity = %d, constr = %d\n",
e->endpoint1,e->endpoint2,e->weight,e->parity,e->constr);
print_info_weak(e->weak);
}
#endif
void free_edge(edge *e)
{
if ( e->weak != NULL )
free_info_weak(e->weak);
free(e);
}
#ifdef PRINT_CUTS
void print_sep_graph(separation_graph *s_graph)
{
int nnodes, maxedges, jk;
nnodes = s_graph->nnodes;
maxedges = (nnodes * (nnodes - 1)) / 2;
printf("\n content of separation_graph: nnodes = %d, nedges = %d\n",
nnodes, s_graph->nedges);
print_int_vect(const_cast<char*>("nodes"),s_graph->nodes,nnodes);
print_int_vect(const_cast<char*>("ind"),s_graph->ind,nnodes);
for ( jk = 0; jk < maxedges; jk++ ) {
if ( s_graph->even_adj_list[jk] != NULL )
print_edge(s_graph->even_adj_list[jk]);
if ( s_graph->odd_adj_list[jk] != NULL )
print_edge(s_graph->odd_adj_list[jk]);
}
}
#endif
void free_sep_graph(separation_graph *s_graph)
{
int nnodes, maxedges, jk;
nnodes = s_graph->nnodes;
maxedges = (nnodes * (nnodes - 1)) / 2;
for ( jk = 0; jk < maxedges; jk++ ) {
if ( s_graph->even_adj_list[jk] != NULL )
free_edge(s_graph->even_adj_list[jk]);
if ( s_graph->odd_adj_list[jk] != NULL )
free_edge(s_graph->odd_adj_list[jk]);
}
free(s_graph->nodes);
free(s_graph->ind);
free(s_graph->even_adj_list);
free(s_graph->odd_adj_list);
free(s_graph);
}
/* auxiliary graph subroutines - depend on the shortest path code used */
#ifndef CGL_NEW_SHORT
// will error if we get here
#include "Cgldikbd.c"
#endif
#define AG_TWIN1(J) 2 * J
#define AG_TWIN2(J) 2 * J + 1
#define AG_MATE(J) 2 * static_cast<int> ( J / 2 ) + ( J % 2 == EVEN ? 1 : 0 )
#define AG_TYPE(J,K) ( (J % 2) == (K % 2) ? EVEN : ODD )
#define SG_ORIG(J) static_cast<int> (J / 2)
/* define_aux_graph: construct the auxiliary graph for the shortest
path computation - the data structure is based on that used by
Cherkassky, Goldberg and Radzik's shortest path codes */
auxiliary_graph *define_aux_graph(separation_graph *s_graph /* input separation graph */)
{
int j, k, indjk, auxj1, auxj2, auxk1, auxk2, noutj, totoutj, narcs;
edge *s_edge;
auxiliary_graph *a_graph;
a_graph = reinterpret_cast<auxiliary_graph *> (calloc(1,sizeof(auxiliary_graph)));
if ( a_graph == NULL ) alloc_error(const_cast<char*>("a_graph"));
a_graph->nnodes = 2 * s_graph->nnodes;
a_graph->narcs = 4 * s_graph->nedges;
#ifndef CGL_NEW_SHORT
a_graph->nodes = reinterpret_cast<node *> (calloc((a_graph->nnodes + 1),sizeof(node)));
#else
a_graph->nodes = reinterpret_cast<cgl_node *> (calloc((a_graph->nnodes + 1),sizeof(cgl_node)));
#endif
if ( a_graph->nodes == NULL ) alloc_error(const_cast<char*>("a_graph->nodes"));
#ifndef CGL_NEW_SHORT
a_graph->arcs = reinterpret_cast<arc *> (calloc(((a_graph->narcs) + 1),sizeof(arc)));
#else
a_graph->arcs = reinterpret_cast<cgl_arc *> (calloc(((a_graph->narcs) + 1),sizeof(cgl_arc)));
#endif
if ( a_graph->arcs == NULL ) alloc_error(const_cast<char*>("a_graph->arcs"));
narcs = 0;
for ( j = 0; j < s_graph->nnodes; j++ ) {
/* count the number of edges incident with j in the separation graph */
totoutj = 0;
for ( k = 0; k < s_graph->nnodes; k++ )
if ( k != j ) {
indjk = SG_EDGE_INDEX(s_graph,j,k);
if ( s_graph->even_adj_list[indjk] != NULL ) totoutj++;
if ( s_graph->odd_adj_list[indjk] != NULL ) totoutj++;
}
auxj1 = AG_TWIN1(j); auxj2 = AG_TWIN2(j);
a_graph->nodes[auxj1].index = auxj1;
a_graph->nodes[auxj2].index = auxj2;
#ifndef CGL_NEW_SHORT
a_graph->nodes[auxj1].first = &(a_graph->arcs[narcs]);
a_graph->nodes[auxj2].first = &(a_graph->arcs[narcs+totoutj]);
#else
a_graph->nodes[auxj1].firstArc = &(a_graph->arcs[narcs]);
a_graph->nodes[auxj2].firstArc = &(a_graph->arcs[narcs+totoutj]);
#endif
/* add the edges as arcs outgoing from j to the auxiliary graph */
noutj = 0;
for ( k = 0; k < s_graph->nnodes; k++ ) {
if ( k != j ) {
auxk1 = AG_TWIN1(k); auxk2 = AG_TWIN2(k);
indjk = SG_EDGE_INDEX(s_graph,j,k);
s_edge = s_graph->even_adj_list[indjk];
if ( s_edge != NULL ) {
/* there is an even edge between j and k */
#ifndef CGL_NEW_SHORT
a_graph->arcs[narcs].len = a_graph->arcs[narcs+totoutj].len =
(int) (s_edge->weight * ISCALE);
a_graph->arcs[narcs].head = &(a_graph->nodes[auxk1]);
a_graph->arcs[narcs+totoutj].head = &(a_graph->nodes[auxk2]);
#else
a_graph->arcs[narcs].length = a_graph->arcs[narcs+totoutj].length =
static_cast<int> (s_edge->weight * ISCALE);
a_graph->arcs[narcs].to = auxk1;
a_graph->arcs[narcs+totoutj].to = auxk2;
#endif
narcs++; noutj++;
}
s_edge = s_graph->odd_adj_list[indjk];
if ( s_edge != NULL ) {
/* there is an odd edge between j and k */
#ifndef CGL_NEW_SHORT
a_graph->arcs[narcs].len = a_graph->arcs[narcs+totoutj].len =
(int) (s_edge->weight * ISCALE);
a_graph->arcs[narcs].head = &(a_graph->nodes[auxk2]);
a_graph->arcs[narcs+totoutj].head = &(a_graph->nodes[auxk1]);
#else
a_graph->arcs[narcs].length = a_graph->arcs[narcs+totoutj].length =
static_cast<int> (s_edge->weight * ISCALE);
a_graph->arcs[narcs].to = auxk2;
a_graph->arcs[narcs+totoutj].to = auxk1;
#endif
/* this looks really useless - to be removed ...
if ( noutj == 0 ) {
a_graph->nodes[auxj1].first = &(a_graph->arcs[narcs]);
a_graph->nodes[auxj2].first = &(a_graph->arcs[narcs+totoutj]);
}
... */
narcs++; noutj++;
}
}
}
narcs += totoutj;
}
#ifndef CGL_NEW_SHORT
a_graph->nodes[a_graph->nnodes].first = &(a_graph->arcs[narcs]);
#else
a_graph->nodes[a_graph->nnodes].firstArc = &(a_graph->arcs[narcs]);
#endif
return(a_graph);
}
/* cancel_node_aux_graph: remove the node j in the separation graph
from the auxiliary graph - all the outgoing arc lengths are set to
a large value */
auxiliary_graph *cancel_node_aux_graph(
int j, /* index of the node in the separation graph */
auxiliary_graph *a_graph /* auxiliary graph to be updated */
)
{
int auxj1, auxj2;
#ifndef CGL_NEW_SHORT
arc *arc_ptr;
#else
cgl_arc *arc_ptr;
#endif
auxj1 = AG_TWIN1(j); auxj2 = AG_TWIN2(j);
#ifndef CGL_NEW_SHORT
for ( arc_ptr = a_graph->nodes[auxj1].first;
arc_ptr < a_graph->nodes[auxj1+1].first;
arc_ptr++ )
(*arc_ptr).len = ISCALE;
for ( arc_ptr = a_graph->nodes[auxj2].first;
arc_ptr < a_graph->nodes[auxj2+1].first;
arc_ptr++ )
(*arc_ptr).len = ISCALE;
#else
for ( arc_ptr = a_graph->nodes[auxj1].firstArc;
arc_ptr < a_graph->nodes[auxj1+1].firstArc;
arc_ptr++ )
(*arc_ptr).length = ISCALE;
for ( arc_ptr = a_graph->nodes[auxj2].firstArc;
arc_ptr < a_graph->nodes[auxj2+1].firstArc;
arc_ptr++ )
(*arc_ptr).length = ISCALE;
#endif
return(a_graph);
}
#ifdef PRINT_CUTS
#ifndef CGL_NEW_SHORT
void print_node(node *n)
{
printf("\n content of node (addr = %d): first = %d, dist = %d, parent = %d, next = %d, prev = %d, status = %d\n",
(int)n,(int)(*n).first,(*n).dist,(int)(*n).parent,(int)(*n).next,(int)(*n).prev,
(*n).status);
}
void print_arc(arc *a)
{
printf("\n content of arc (addr = %d): len = %d, head = %d\n",
(int)a,(*a).len,(int)(*a).head);
}
void print_node_vect(char *s,node *v,int n)
{
int i;
printf("node vector %s:",s);
for ( i = 0; i < n; i++ ) print_node(&v[i]);
printf("\n");
}
void print_arc_vect(char *s,arc *v,int n)
{
int i;
printf("arc vector %s:",s);
for ( i = 0; i < n; i++ ) print_arc(&v[i]);
printf("\n");
}
#else
void print_node(cgl_node *n)
{
printf("\n content of node (addr = %p): first = %p, dist = %d, parent = %p\n",
reinterpret_cast<void *>(n),
reinterpret_cast<void *>((*n).firstArc),
(*n).distanceBack,
static_cast<int>((*n).parentNode));
}
void print_arc(cgl_arc *a)
{
printf("\n content of arc (addr = %p): len = %d, head = %d\n",
reinterpret_cast<void *>(a),(*a).length,static_cast<int>((*a).to));
}
void print_node_vect(char *s,cgl_node *v,int n)
{
int i;
printf("node vector %s:",s);
for ( i = 0; i < n; i++ ) print_node(&v[i]);
printf("\n");
}
void print_arc_vect(char *s,cgl_arc *v,int n)
{
int i;
printf("arc vector %s:",s);
for ( i = 0; i < n; i++ ) print_arc(&v[i]);
printf("\n");
}
#endif
void print_aux_graph(auxiliary_graph *a_graph)
{
printf("\n content of auxiliary graph: nnodes = %d, narcs = %d\n",
a_graph->nnodes,a_graph->narcs);
print_node_vect(const_cast<char*>("nodes"),a_graph->nodes,a_graph->nnodes);
print_arc_vect(const_cast<char*>("nodes"),a_graph->arcs,a_graph->narcs);
}
#endif
void free_aux_graph(auxiliary_graph *a_graph)
{
free(a_graph->nodes);
free(a_graph->arcs);
free(a_graph);
}
/* odd cycles management subroutines */
/* simple_cycle: check whether a given cycle is simple
(and therefore may correspond to a non-dominated ineq.) */
short int simple_cycle(cycle *s_cyc /* cycle to be checked */)
{
int i, e, maxnodes;
int *cnt;
maxnodes = 0;
for ( e = 0; e < s_cyc->length; e++ ) {
if (!s_cyc->edge_list[e]) {
// bad
maxnodes=-1;
abort();//break;
}
i = s_cyc->edge_list[e]->endpoint1;
if ( i > maxnodes ) maxnodes = i;
i = s_cyc->edge_list[e]->endpoint2;
if ( i > maxnodes ) maxnodes = i;
}
if (maxnodes<0)
return FALSE;
cnt = reinterpret_cast<int *> (calloc(maxnodes+1,sizeof(int)));
if ( cnt == NULL ) alloc_error(const_cast<char*>("cnt"));
//for ( i = 0; i <= maxnodes; i++ ) cnt[i] = 0;
for ( e = 0; e < s_cyc->length; e++ ) {
i = s_cyc->edge_list[e]->endpoint1;
cnt[i]++;
if ( cnt[i] > 2 ) {
free(cnt);
return(FALSE);
}
i = s_cyc->edge_list[e]->endpoint2;
cnt[i]++;
if ( cnt[i] > 2 ) {
free(cnt);
return(FALSE);
}
}
free(cnt);
return(TRUE);
}
/* same_cycle: check whether two cycles are identical
(assumes the first nodes of the cycles coincide) */
short int same_cycle(cycle *s_cyc1, cycle *s_cyc2 /* cycles to be compared */)
{
int e, eb;
short int same;
if ( s_cyc1->length != s_cyc2->length ) return(FALSE);
/* check the cycles in the same direction ... */
same = TRUE;
for ( e = 0; e < s_cyc1->length; e++ ) {
if ( s_cyc1->edge_list[e] != s_cyc2->edge_list[e] ) {
same = FALSE;
break;
}
}
if ( same ) return(TRUE);
/* ... and in reverse direction */
same = TRUE;
for ( e = 0, eb = s_cyc2->length - 1; e < s_cyc1->length; e++, eb-- ) {
if ( s_cyc1->edge_list[e] != s_cyc2->edge_list[eb] ) {
same = FALSE;
break;
}
}
if ( same ) return(TRUE);
return(FALSE);
}
#ifdef PRINT_CUTS
void print_cycle(cycle *s_cycle)
{
int e;
printf("\n content of cycle: weight = %f, length = %d\n",
s_cycle->weight,s_cycle->length);
for ( e = 0; e < s_cycle->length; e++ )
print_edge(s_cycle->edge_list[e]);
}
#endif
void free_cycle(cycle *s_cycle)
{
free(s_cycle->edge_list);
free(s_cycle);
}
/* initialize_cycle_list: allocate and initialize the cycle list data structure */
cycle_list *initialize_cycle_list(int max_cyc /* maximum number of cycles in the list */)
{
cycle_list *s_cycle_list;
s_cycle_list = reinterpret_cast<cycle_list *> (calloc(1,sizeof(cycle_list)));
if ( s_cycle_list == NULL ) alloc_error(const_cast<char*>("s_cycle_list"));
s_cycle_list->cnum = 0;
s_cycle_list->list = reinterpret_cast<cycle **> (calloc(max_cyc,sizeof(cycle *)));
if ( s_cycle_list->list == NULL ) alloc_error(const_cast<char*>("s_cycle_list->list"));
return(s_cycle_list);
}
/* add_cycle_to_list: add a new cycle to the cycle list data structure
(if not already in the list) */
cycle_list *add_cycle_to_list(
cycle *s_cycle, /* pointer to the cycle to be added to the list */
cycle_list *s_cycle_list /* input cycle list to be updated */
)
{
int c;
if ( ! simple_cycle(s_cycle) ) {
free_cycle(s_cycle);
return(s_cycle_list);
}
for ( c = 0; c < s_cycle_list->cnum; c++ )
if ( same_cycle(s_cycle,s_cycle_list->list[c]) ) {
free_cycle(s_cycle);
return(s_cycle_list);
}
s_cycle_list->list[s_cycle_list->cnum] = s_cycle;
(s_cycle_list->cnum)++;
return(s_cycle_list);
}
void free_cycle_list(cycle_list *s_cycle_list)
{
int c;
for ( c = 0; c < s_cycle_list->cnum; c++ )
free_cycle(s_cycle_list->list[c]);
free(s_cycle_list->list);
free(s_cycle_list);
}
#ifdef PRINT_CUTS
void print_cycle_list(cycle_list *s_cycle_list)
{
int c;
printf("\n content of cycle_list: cnum = %d\n",s_cycle_list->cnum);
for ( c = 0; c < s_cycle_list->cnum; c++ )
print_cycle(s_cycle_list->list[c]);
}
#endif
/* get_shortest_odd_cycle_list: computation of the shortest odd cycles
visiting a certain node in the separation graph, and each other
possible intermediate node, using the auxiliary graph data structure
for the shortest path computation - all the cycles in the list are
different from each other */
cycle_list *get_shortest_odd_cycle_list(
int j, /* first node to be visited by the odd cycle */
separation_graph *s_graph, /* current separation graph */
auxiliary_graph *a_graph /* auxiliary graph for the shortest path computation */
)
{
int source, sink, curr, pred, totedges, k, t, kt;
double weight;
#ifndef CGL_NEW_SHORT
//node *source_ptr, *sink_ptr, *first_ptr;
#else
//cgl_node *source_ptr, *sink_ptr, *first_ptr;
#endif
edge *curr_edge;
short_path_node *forw_arb, *backw_arb;
cycle *s_cycle;
cycle_list *s_cycle_list;
#ifdef TIME
second_(&tsi);
#endif
s_cycle_list = initialize_cycle_list((a_graph->nnodes)-2);
source = AG_TWIN1(j); sink = AG_TWIN2(j);
//source_ptr = &(a_graph->nodes[source]);
//sink_ptr = &(a_graph->nodes[sink]);
//first_ptr = &(a_graph->nodes[0]);
/* compute the shortest path arborescence rooted at source and
the shortest path arborescence rooted at sink (that comes for
free due to symmetry) and store them (the path information is
hidden into aux_graph) */
#ifdef TIME
second_(&ti);
#endif
#ifndef CGL_NEW_SHORT
{
int nNodes = a_graph->nnodes;
int nArcs = a_graph->narcs;
cgl_arc * arcs = new cgl_arc [nArcs];
for (int i=0;i<nArcs;i++) {
arcs[i].length=a_graph->arcs[i].len;
arcs[i].to=a_graph->arcs[i].head->index;
}
cgl_node * nodes = new cgl_node[nNodes+1];
for (int i=0;i<nNodes;i++) {
int iArc = a_graph->nodes[i].first-a_graph->arcs;
nodes[i].firstArc=arcs+iArc;
nodes[i].index=i;
}
int iArc = a_graph->nodes[nNodes].first-a_graph->arcs;
nodes[nNodes].firstArc=arcs+iArc;
cgl_graph graph;
graph.nnodes=nNodes;
graph.narcs=nArcs;
graph.nodes=nodes;
graph.arcs=arcs;
cglShortestPath(&graph,source,ISCALE);
dikbd(a_graph->nnodes,first_ptr,source_ptr,ISCALE);
for ( k = 0; k < a_graph->nnodes; k++ ) {
if ( a_graph->nodes[k].parent != NULL ) {
int distance1 = a_graph->nodes[k].dist;
int distance2 = graph.nodes[k].distanceBack;
assert (distance1==distance2);
} else {
//
printf("null parent %d\n",k);
}
}
}
#else
cglShortestPath(a_graph,source,ISCALE);
#endif
#ifdef TIME
second_(&tf);
path_time += tf - ti;
#endif
forw_arb =
reinterpret_cast<short_path_node *> (calloc(a_graph->nnodes,sizeof(short_path_node)));
if ( forw_arb == NULL ) alloc_error(const_cast<char*>("forw_arb"));
for ( k = 0; k < a_graph->nnodes; k++ ) {
#ifndef CGL_NEW_SHORT
if ( a_graph->nodes[k].parent != NULL ) {
forw_arb[k].dist = a_graph->nodes[k].dist;
forw_arb[k].pred = a_graph->nodes[k].parent->index;
}
#else
if ( a_graph->nodes[k].parentNode >=0 ) {
forw_arb[k].dist = a_graph->nodes[k].distanceBack;
forw_arb[k].pred = a_graph->nodes[k].parentNode;
}
#endif
else {
forw_arb[k].dist = COIN_INT_MAX;
forw_arb[k].pred = NONE;
}
}
backw_arb =
reinterpret_cast<short_path_node *> (calloc(a_graph->nnodes,sizeof(short_path_node)));
if ( backw_arb == NULL ) alloc_error(const_cast<char*>("backw_arb"));
for ( k = 0; k < a_graph->nnodes; k++ ) {
#ifndef CGL_NEW_SHORT
if ( a_graph->nodes[k].parent != NULL ) {
backw_arb[AG_MATE(k)].dist = a_graph->nodes[k].dist;
backw_arb[AG_MATE(k)].pred = AG_MATE(a_graph->nodes[k].parent->index);
}
#else
if ( a_graph->nodes[k].parentNode >=0) {
backw_arb[AG_MATE(k)].dist = a_graph->nodes[k].distanceBack;
backw_arb[AG_MATE(k)].pred = AG_MATE(a_graph->nodes[k].parentNode);
}
#endif
else {
backw_arb[AG_MATE(k)].dist = COIN_INT_MAX;
backw_arb[AG_MATE(k)].pred = NONE;
}
}
#ifdef USELESS
/* compute second the shortest path anti-arborescence rooted at sink
(which coincides with the arborescence since aux_graph is
symmetrical) and store it */
#ifdef TIME
second_(&ti);
#endif
cc = dikbd(a_graph->nnodes,first_ptr,sink_ptr,ISCALE);
#ifdef TIME
second_(&tf);
path_time += tf - ti;
#endif
backw_arb =
(short_path_node *) calloc(a_graph->nnodes,sizeof(short_path_node));
if ( backw_arb == NULL ) alloc_error("backw_arb");
for ( k = 0; k < a_graph->nnodes; k++ ) {
backw_arb[k].dist = a_graph->nodes[k].dist;
backw_arb[k].pred = a_graph->nodes[k].parent->index;
}
#endif
/* consider each possible intermediate node in aux_graph */
for ( k = 0; k < s_graph->nnodes; k++ ) {
if ( k != j ) {
for ( t = 1; t <= 2; t++ ) {
if ( t == 1 ) kt = AG_TWIN1(k);
else kt = AG_TWIN2(k);
weight =
(static_cast<double> (forw_arb[kt].dist + backw_arb[kt].dist)) /
(static_cast<double> (ISCALE));
if ( weight < MAX_CYCLE_WEIGHT + EPS ) {
totedges = 0;
/* count how many edges are in the forward path from source ... */
curr = kt;
do {
if (curr<0) {
totedges=-1;
break;
}
curr = forw_arb[curr].pred; totedges++;
} while ( curr != source );
if (totedges>=0) {
/* ... and in the backward path to sink */
curr = kt;
do {
if (curr<0) {
totedges=-1;
break;
}
curr = backw_arb[curr].pred; totedges++;
} while ( curr != sink );
}
if (totedges>0) {
s_cycle = reinterpret_cast<cycle *> (calloc(1,sizeof(cycle)));
if ( s_cycle == NULL ) alloc_error(const_cast<char*>("s_cycle"));
s_cycle->weight = weight;
s_cycle->length = totedges;
s_cycle->edge_list = reinterpret_cast<edge **> (calloc(totedges,sizeof(edge *)));
if ( s_cycle->edge_list == NULL ) alloc_error(const_cast<char*>("s_cycle->edge_list"));
/* define the set of edges corresponding to the paths in sep_graph */
totedges = 0;
/* forward path from source ... */
curr = kt;
do {
pred = forw_arb[curr].pred;
if ( AG_TYPE(pred,curr) == EVEN )
curr_edge = s_graph->even_adj_list
[SG_EDGE_INDEX(s_graph,SG_ORIG(curr),SG_ORIG(pred))];
else
curr_edge = s_graph->odd_adj_list
[SG_EDGE_INDEX(s_graph,SG_ORIG(curr),SG_ORIG(pred))];
s_cycle->edge_list[totedges] = curr_edge;
curr = pred;
totedges++;
} while ( curr != source );
/* ... and backward path to sink */
curr = kt;
do {
pred = backw_arb[curr].pred;
if ( AG_TYPE(pred,curr) == EVEN )
curr_edge = s_graph->even_adj_list
[SG_EDGE_INDEX(s_graph,SG_ORIG(curr),SG_ORIG(pred))];
else
curr_edge = s_graph->odd_adj_list
[SG_EDGE_INDEX(s_graph,SG_ORIG(curr),SG_ORIG(pred))];
s_cycle->edge_list[totedges] = curr_edge;
curr = pred;
totedges++;
} while ( curr != sink );
/* insert the new cycle in the list */
s_cycle_list = add_cycle_to_list(s_cycle,s_cycle_list);
}
}
}
}
}
free(forw_arb);
free(backw_arb);
#ifdef TIME
second_(&tsf);
cycle_time += tsf - tsi;
#endif
return(s_cycle_list);
}
/* cut management subroutines */
/* initialize_cut_list: allocate and initialize the cut list data structure */
cut_list *initialize_cut_list(int max_cut /* maximum number of cuts in the list */)
{
cut_list *cuts;
cuts = reinterpret_cast<cut_list *> (calloc(1,sizeof(cut_list)));
if ( cuts == NULL ) alloc_error(const_cast<char*>("cuts"));
cuts->cnum = 0;
cuts->list = reinterpret_cast<cut **> (calloc(max_cut,sizeof(cut *)));
return(cuts);
}
#ifdef PRINT_CUTS
void Cgl012Cut::print_cut(cut *v_cut)
{
printf("\n content of cut: n_of_constr = %d, cnzcnt = %d, crhs = %d, csense = %c, violation = %f\n",
v_cut->n_of_constr,v_cut->cnzcnt,v_cut->crhs,v_cut->csense,v_cut->violation);
print_int_vect(const_cast<char*>("cind"),v_cut->cind,v_cut->cnzcnt);
print_int_vect(const_cast<char*>("cval"),v_cut->cval,v_cut->cnzcnt);
if ( v_cut->constr_list != NULL )
print_int_vect(const_cast<char*>("constr_list"),v_cut->constr_list,v_cut->n_of_constr);
if ( v_cut->in_constr_list != NULL )
print_short_int_vect(const_cast<char*>("in_constr_list"),v_cut->in_constr_list,inp_ilp->mr);
;
}
void Cgl012Cut::print_cut_list(cut_list *cuts)
{
int c;
printf("\n content of cut_list: cnum = %d\n",cuts->cnum);
for ( c = 0; c < cuts->cnum; c++ )
print_cut(cuts->list[c]);
}
#endif
void free_cut(cut *v_cut)
{
if ( v_cut->constr_list != NULL ) free(v_cut->constr_list);
if ( v_cut->in_constr_list != NULL ) free(v_cut->in_constr_list);
if ( v_cut->cind != NULL ) free(v_cut->cind);
if ( v_cut->cval != NULL ) free(v_cut->cval);
free(v_cut);
}
void free_cut_list(cut_list *cuts)
{
int c;
for ( c = 0; c < cuts->cnum; c++ )
if ( cuts->list[c] != NULL ) free_cut(cuts->list[c]);
free(cuts->list);
free(cuts);
}
/* get_ori_cut_coef: get the coefficients of a cut, before dividing by 2 and
rounding, starting from the list of the constraints combined to get
the cut */
short int Cgl012Cut::get_ori_cut_coef(
int n_of_constr, /* number of constraints combined */
int *constr_list, /* list of the constraints combined */
int *ccoef, /* cut left hand side coefficients */
int *crhs, /* cut right hand side */
short int only_viol /* flag which tells whether only an inequality of
slack smaller than MAX_SLACK is of interest (TRUE)
otherwise (FALSE) */
)
{
int h, i, begi, gcdi, ofsj;
double tot_slack;
/* fast check of the possible violation of the cut */
if ( only_viol ) {
tot_slack = 0.0;
for ( h = 0; h < n_of_constr; h++ ) {
tot_slack += p_ilp->slack[constr_list[h]];
if ( tot_slack > MAX_SLACK - EPS ) return(FALSE);
}
}
//for ( j = 0; j < inp_ilp->mc; j++ )
//ccoef[j] = 0;
memset(ccoef,0,inp_ilp->mc*sizeof(int));
(*crhs) = 0;
for ( h = 0; h < n_of_constr; h++ ) {
i = constr_list[h];
begi = inp_ilp->mtbeg[i]; gcdi = p_ilp->gcd[i];
if ( inp_ilp->msense[i] != 'G' ) {
if ( gcdi == 1 ) {
for ( ofsj = 0; ofsj < inp_ilp->mtcnt[i]; ofsj++ )
ccoef[inp_ilp->mtind[begi+ofsj]] += inp_ilp->mtval[begi+ofsj];
(*crhs) += inp_ilp->mrhs[i];
}
else {
for ( ofsj = 0; ofsj < inp_ilp->mtcnt[i]; ofsj++ )
ccoef[inp_ilp->mtind[begi+ofsj]] += inp_ilp->mtval[begi+ofsj] / gcdi;
(*crhs) += inp_ilp->mrhs[i] / gcdi;
}
}
else {
if ( gcdi == 1 ) {
for ( ofsj = 0; ofsj < inp_ilp->mtcnt[i]; ofsj++ )
ccoef[inp_ilp->mtind[begi+ofsj]] -= inp_ilp->mtval[begi+ofsj];
(*crhs) -= inp_ilp->mrhs[i];
}
else {
for ( ofsj = 0; ofsj < inp_ilp->mtcnt[i]; ofsj++ )
ccoef[inp_ilp->mtind[begi+ofsj]] -= inp_ilp->mtval[begi+ofsj] / gcdi;
(*crhs) -= inp_ilp->mrhs[i] / gcdi;
}
}
}
return(TRUE);
}
/* best_cut: find the coefficients, the rhs and the violation of the
best possible cut that can be obtained by weakening a given set of
coefficients to even and a rhs to odd, dividing by 2 and rounding */
short int Cgl012Cut::best_cut(
int *ccoef, /* vector of the coefficients */
int *crhs, /* pointer to rhs value */
double *violation, /* violation of the cut */
short int update, /* TRUE/FALSE: if TRUE, the new ccoef and crhs are
given on output */
short int only_viol /* flag which tells whether only an inequality of
slack smaller than MAX_SLACK is of interest (TRUE)
otherwise (FALSE) */
)
{
int j, n_to_weak;
short int original_parity;
double original_slack, best_even_slack, best_odd_slack;
int *vars_to_weak;
info_weak *info_even_weak, *info_odd_weak;
/* choose the best weakening for the variables whose coefficient
is not odd - this hopefully produces a stronger cut than that
associated with the weakened inequalities to define the edges */
vars_to_weak = reinterpret_cast<int *> (calloc(inp_ilp->mc,sizeof(int)));
if ( vars_to_weak == NULL ) alloc_error(const_cast<char*>("vars_to_weak"));
n_to_weak = 0;
original_slack = 0.0;
for ( j = 0; j < inp_ilp->mc; j++ ) {
if ( ccoef[j] != 0 ) {
if ( mod2(ccoef[j]) == ODD ) {
vars_to_weak[n_to_weak] = j;
n_to_weak++;
}
original_slack -= inp_ilp->xstar[j] * static_cast<double> (ccoef[j]);
}
}
original_slack += static_cast<double> (*crhs);
if ( original_slack > MAX_SLACK - EPS ) {
free(vars_to_weak);
return(FALSE);
}
original_parity = mod2(*crhs);
if ( best_weakening(n_to_weak,vars_to_weak,
original_parity,original_slack,
&best_even_slack,&best_odd_slack,
&info_even_weak,&info_odd_weak,
TRUE,only_viol) == ODD ) {
*violation = ( 1.0 - best_odd_slack ) / 2.0;
if ( ! update ) {
/* new ccoef and rhs are not required on output */
free(vars_to_weak);
free_info_weak(info_odd_weak);
return(TRUE);
}
/* update ccoef and crhs according to the best odd weakening */
for ( j = 0; j < n_to_weak; j++ )
if ( info_odd_weak->type[j] == LOWER_BOUND ) {
ccoef[vars_to_weak[j]]--;
*crhs -= inp_ilp->vlb[vars_to_weak[j]];
}
else {
ccoef[vars_to_weak[j]]++;
*crhs += inp_ilp->vub[vars_to_weak[j]];
}
/* compute and check the correctness of the cut coefficients */
for ( j = 0; j < inp_ilp->mc; j++ ) {
if ( mod2(ccoef[j]) == ODD ) {
printf("!!! Error 2 in weakening a cut !!!\n");
exit(0);
}
if ( ccoef[j] != 0 ) ccoef[j] /= 2;
}
if ( mod2(*crhs) == EVEN ) {
printf("!!! Error 1 in weakening a cut !!!\n");
exit(0);
}
*crhs = (*crhs - 1) / 2;
free(vars_to_weak);
free_info_weak(info_odd_weak);
return(TRUE);
}
else {
free(vars_to_weak);
return(FALSE);
}
}
/* define_cut: construct a cut data structure from a vector of
coefficients and a right-hand-side */
cut *Cgl012Cut::define_cut(
int *ccoef, /* coefficients of the cut */
int crhs /* right hand side of the cut */
)
{
int cnzcnt, j;
cut *v_cut;
v_cut = reinterpret_cast<cut *> (calloc(1,sizeof(cut)));
if ( v_cut == NULL ) alloc_error(const_cast<char*>("v_cut"));
v_cut->crhs = crhs;
cnzcnt = 0;
for ( j = 0; j < inp_ilp->mc; j++ )
if ( ccoef[j] != 0 ) cnzcnt++;
v_cut->cnzcnt = cnzcnt;
v_cut->csense = 'L';
v_cut->cind = reinterpret_cast<int *> (calloc(cnzcnt,sizeof(int)));
if ( v_cut->cind == NULL ) alloc_error(const_cast<char*>("v_cut->cind"));
v_cut->cval = reinterpret_cast<int *> (calloc(cnzcnt,sizeof(int)));
if ( v_cut->cval == NULL ) alloc_error(const_cast<char*>("v_cut->cval"));
cnzcnt = 0; v_cut->violation = 0.0;
for ( j = 0; j < inp_ilp->mc; j++ )
if ( ccoef[j] != 0 ) {
v_cut->cind[cnzcnt] = j;
v_cut->cval[cnzcnt] = ccoef[j];
v_cut->violation += inp_ilp->xstar[j] * static_cast<double> (ccoef[j]);
cnzcnt++;
}
v_cut->violation -= static_cast<double> (crhs);
return(v_cut);
}
/* get_cut: extract a hopefully violated cut from an odd cycle of the
separation graph */
cut *Cgl012Cut::get_cut(
cycle *s_cyc /* shortest odd cycles identified in the separation graph */
)
{
int i, e, crhs;
short int ok;
/* short int original_parity; */
double violation;
/* double original_slack, best_even_slack, best_odd_slack; */
int *ccoef /*, *vars_to_weak */ ;
/* info_weak *info_even_weak, *info_odd_weak; */
cut *v_cut;
int ncomb;
int *comb;
short int *flag_comb;
#ifndef CGGGGG
static int iter = 0;
static double gap, maxgap = 0.0;
#endif
#ifdef TIME
second_(&tsi);
cut_ncalls++;
#endif
/* compute the cut obtained by adding-up all the constraints
corresponding to edges in the cycle, in their non-weak form */
#ifdef TIME
second_(&tii);
#endif
ccoef = reinterpret_cast<int *> (calloc(inp_ilp->mc,sizeof(int)));
if ( ccoef == NULL ) alloc_error(const_cast<char*>("ccoef"));
ncomb = 0;
comb = reinterpret_cast<int *> (calloc(inp_ilp->mr,sizeof(int)));
if ( comb == NULL ) alloc_error(const_cast<char*>("comb"));
flag_comb = reinterpret_cast<short int *> (calloc(inp_ilp->mr,sizeof(short int)));
if ( flag_comb == NULL ) alloc_error(const_cast<char*>("flag_comb"));
#if 0
// no need to as calloc used
for ( i = 0; i < inp_ilp->mr; i++ ) flag_comb[i] = OUT;
for ( j = 0; j < inp_ilp->mc; j++ )
ccoef[j] = 0;
#endif
crhs = 0;
for ( e = 0; e < s_cyc->length; e++ ) {
i = (s_cyc->edge_list[e])->constr;
if ( i >= 0 ) {
/* the edge is not associated with a bound constraint */
comb[ncomb] = i; ncomb++; flag_comb[i] = IN;
}
}
ok = get_ori_cut_coef(ncomb,comb,ccoef,&crhs,TRUE);
#ifdef TIME
second_(&tff);
coef_time += tff - tii;
#endif
ok = ok && best_cut(ccoef,&crhs,&violation,TRUE,TRUE);
if ( ! ok ) {
free(ccoef);
free(comb);
free(flag_comb);
#ifdef TIME
second_(&tsf);
cut_time += tsf - tsi;
#endif
return(NULL);
}
v_cut = define_cut(ccoef,crhs);
iter++;
if ( v_cut->violation > violation + EPS ||
v_cut->violation < violation - EPS ) {
//printf("Error in violation check\n");
//printf("v_cut->violation %f violation %f gap %f maxgap (previous) %f\n",
// v_cut->violation,violation,v_cut->violation-violation,maxgap);
//printf("iter %d\n",iter);
//exit(0);
free_cut(v_cut);
free(ccoef);
free(comb);
free(flag_comb);
errorNo=1;
return(NULL);
}
gap = v_cut->violation - violation;
if ( gap < 0.0 ) gap = -gap;
if ( gap > maxgap ) maxgap = gap;
v_cut->n_of_constr = ncomb;
v_cut->constr_list = comb;
v_cut->in_constr_list = flag_comb;
free(ccoef);
#ifdef TIME
second_(&tsf);
cut_time += tsf - tsi;
#endif
return(v_cut);
}
/* cut_score: define the score of a (violated) cut */
double Cgl012Cut::cut_score(
int *ccoef, /* cut left hand side coefficients */
int crhs, /* cut right hand side */
double viol, /* cut violation */
short int only_viol /* flag which tells whether only an inequality of
slack smaller than MAX_SLACK is of interest (TRUE)
otherwise (FALSE) */
)
{
int j, norm;
/* very simple score: violation divided/multiplied by the lhs norm */
if ( only_viol && viol < MIN_VIOLATION ) return(-INF);
norm = 0;
for ( j = 0; j < p_ilp->mc; j++ ) {
if ( ccoef[j] != 0 ) norm += ccoef[j] * ccoef[j];
}
if ( viol > 0.0 ) return (viol / sqrt(static_cast<double> (norm)));
else return (viol * sqrt(static_cast<double> (norm)));
}
/* same_cut: check whether two cuts are identical - not too clever
(assumes the sparse coefficients are sorted by column index) */
short int same_cut(cut *cut1, cut *cut2 /* cuts to be compared */)
{
int j;
if ( cut1->cnzcnt != cut2->cnzcnt ) return(FALSE);
if ( cut1->crhs != cut2->crhs ) return(FALSE);
if ( cut1->csense != cut2->csense ) return(FALSE);
for ( j = 0; j < cut1->cnzcnt; j++ ) {
if ( cut1->cind[j] != cut2->cind[j] ) return(FALSE);
if ( cut1->cval[j] != cut2->cval[j] ) return(FALSE);
}
return(TRUE);
}
/* add_cut_to_list: adds a cut to a list after checking that a copy of
the same cut is not already in the list - no checking is made about
cuts dominating each other or implied by other cuts in the list plus
the constraints of the original problem */
cut_list *add_cut_to_list(
cut *v_cut, /* pointer to the violated cut to be added to the list */
cut_list *cuts /* input cut list to be updated */
)
{
int c;
for ( c = 0; c < cuts->cnum; c++ ) {
if ( same_cut(v_cut,cuts->list[c]) ) {
free_cut(v_cut);
return(cuts);
}
}
cuts->list[cuts->cnum] = v_cut;
cuts->cnum++;
return(cuts);
}
/* getcuts: pick the 0-1/2 cuts in the list and give them on output */
void getcuts(
cut_list *cuts, /* input cut list */
int *cnum, /* number of violated 0-1/2 cuts identified by the procedure */
int *cnzcnt, /* overall number of nonzero's in the cuts */
int **cbeg, /* starting position of each cut in arrays cind and cval */
int **ccnt, /* number of entries of each cut in arrays cind and cval */
int **cind, /* column indices of the nonzero entries of the cuts */
int **cval, /* values of the nonzero entries of the cuts */
int **crhs, /* right hand sides of the cuts */
char **csense /* senses of the cuts: 'L', 'G' or 'E' */
)
{
int i, ofsj, count;
cut *cut_ptr;
/* allocate the memory for the output vectors */
(*cnum) = cuts->cnum;
(*cnzcnt) = 0;
for ( i = 0; i < cuts->cnum; i++ )
(*cnzcnt) += (cuts->list[i])->cnzcnt;
/* if ( (*cbeg) != NULL ) free(*cbeg); */
(*cbeg) = reinterpret_cast<int *> (calloc((*cnum),sizeof(int)));
if ( (*cbeg) == NULL ) alloc_error(const_cast<char*>("*cbeg"));
/* if ( (*ccnt) != NULL ) free(*ccnt); */
(*ccnt) = reinterpret_cast<int *> (calloc((*cnum),sizeof(int)));
if ( (*ccnt) == NULL ) alloc_error(const_cast<char*>("*ccnt"));
/* if ( (*crhs) != NULL ) free(*crhs); */
(*crhs) = reinterpret_cast<int *> (calloc((*cnum),sizeof(int)));
if ( (*crhs) == NULL ) alloc_error(const_cast<char*>("*crhs"));
/* if ( (*csense) != NULL ) free(*csense); */
(*csense) = reinterpret_cast<char *> (calloc((*cnum),sizeof(char)));
if ( (*csense) == NULL ) alloc_error(const_cast<char*>("*csense"));
/* if ( (*cind) != NULL ) free(*cind); */
(*cind) = reinterpret_cast<int *> (calloc((*cnzcnt),sizeof(int)));
if ( (*cind) == NULL ) alloc_error(const_cast<char*>("*cind"));
/* if ( (*cval) != NULL ) free(*cval); */
(*cval) = reinterpret_cast<int *> (calloc((*cnzcnt),sizeof(int)));
if ( (*cval) == NULL ) alloc_error(const_cast<char*>("*cval"));
/* transfer the cuts information into the output data structures */
count = 0;
for ( i = 0; i < cuts->cnum; i++ ) {
cut_ptr = cuts->list[i];
(*cbeg)[i] = count;
(*ccnt)[i] = cut_ptr->cnzcnt;
(*crhs)[i] = cut_ptr->crhs;
(*csense)[i] = cut_ptr->csense;
for ( ofsj = 0; ofsj < cut_ptr->cnzcnt; ofsj++ ) {
(*cind)[count] = cut_ptr->cind[ofsj];
(*cval)[count] = cut_ptr->cval[ofsj];
count++;
}
}
}
/* actual separation subroutines */
/* best_weakening: find the best upper/lower bound weakening of a set
of variables */
int Cgl012Cut::best_weakening(
int n_to_weak, /* number of variables to weaken */
int *vars_to_weak, /* indices of the variables to weaken */
short int original_parity, /* original parity of the constraint to weaken */
double original_slack, /* original slack of the constraint to weaken */
double *best_even_slack, /* best possible slack of a weakened constraint
with even right-hand-side */
double *best_odd_slack, /* best possible slack of a weakened constraint
with odd right-hand-side */
info_weak **info_even_weak, /* weakening information about the best possible
even weakened constraint */
info_weak **info_odd_weak, /* weakening information about the best possible
odd weakened constraint */
short int only_odd, /* flag which tells whether only an odd weakening is of
interest (TRUE) or both weakenings are (FALSE) */
short int only_viol /* flag which tells whether only an inequality of
slack smaller than MAX_SLACK is of interest (TRUE)
otherwise (FALSE) */
)
{
int nweak, cntweak, ofsl, l;
short int flag_even, flag_odd, ok_even, ok_odd;
double best_even_e, best_even_o, best_odd_e, best_odd_o;
short int *type_even_weak, *type_odd_weak,
*switch_even_weak, *switch_odd_weak;
#ifdef TIME
second_(&tii);
#endif
type_even_weak = reinterpret_cast<short int *> (calloc(p_ilp->mc,sizeof(short int)));
if (type_even_weak == NULL ) alloc_error(const_cast<char*>("type_even_weak"));
switch_even_weak = reinterpret_cast<short int *> (calloc(p_ilp->mc,sizeof(short int)));
if (switch_even_weak == NULL ) alloc_error(const_cast<char*>("switch_even_weak"));
type_odd_weak = reinterpret_cast<short int *> (calloc(p_ilp->mc,sizeof(short int)));
if (type_odd_weak == NULL ) alloc_error(const_cast<char*>("type_odd_weak"));
switch_odd_weak = reinterpret_cast<short int *> (calloc(p_ilp->mc,sizeof(short int)));
if (switch_odd_weak == NULL ) alloc_error(const_cast<char*>("switch_odd_weak"));
if ( original_parity == EVEN ) {
(*best_even_slack) = original_slack;
(*best_odd_slack) = INF;
}
else {
(*best_odd_slack) = original_slack;
(*best_even_slack) = INF;
}
nweak = 0;
for ( ofsl = 0; ofsl < n_to_weak; ofsl++ ) {
l = vars_to_weak[nweak];
if ( p_ilp->possible_weak[l] == NONE ) {
free(type_even_weak); free(type_odd_weak);
free(switch_even_weak); free(switch_odd_weak);
#ifdef TIME
second_(&tff);
bw_time += tff - tii;
#endif
return(NONE);
}
else if ( p_ilp->possible_weak[l] == EVEN ) {
/* only even weakening of l is possible */
(*best_even_slack) += p_ilp->loss_even_weak[l];
type_even_weak[nweak] = p_ilp->type_even_weak[l];
switch_even_weak[nweak] = FALSE;
(*best_odd_slack) += p_ilp->loss_even_weak[l];
type_odd_weak[nweak] = p_ilp->type_even_weak[l];
switch_odd_weak[nweak] = FALSE;
}
else if ( p_ilp->possible_weak[l] == ODD ) {
/* only odd weakening of l is possible */
best_even_e = (*best_even_slack);
best_odd_o = (*best_odd_slack);
(*best_even_slack) = best_odd_o + p_ilp->loss_odd_weak[l];
type_even_weak[nweak] = p_ilp->type_odd_weak[l];
switch_even_weak[nweak] = TRUE;
(*best_odd_slack) = best_even_e + p_ilp->loss_odd_weak[l];
type_odd_weak[nweak] = p_ilp->type_odd_weak[l];
switch_odd_weak[nweak] = TRUE;
}
else {
/* both weakenings of l are possible */
best_even_e = (*best_even_slack) + p_ilp->loss_even_weak[l];
best_even_o = (*best_odd_slack) + p_ilp->loss_odd_weak[l];
best_odd_e = (*best_odd_slack) + p_ilp->loss_even_weak[l];
best_odd_o = (*best_even_slack) + p_ilp->loss_odd_weak[l];
if ( best_even_e <= best_even_o ) {
(*best_even_slack) = best_even_e;
type_even_weak[nweak] = p_ilp->type_even_weak[l];
switch_even_weak[nweak] = FALSE;
}
else {
(*best_even_slack) = best_even_o;
type_even_weak[nweak] = p_ilp->type_odd_weak[l];
switch_even_weak[nweak] = TRUE;
}
if ( best_odd_e <= best_odd_o ) {
(*best_odd_slack) = best_odd_e;
type_odd_weak[nweak] = p_ilp->type_even_weak[l];
switch_odd_weak[nweak] = FALSE;
}
else {
(*best_odd_slack) = best_odd_o;
type_odd_weak[nweak] = p_ilp->type_odd_weak[l];
switch_odd_weak[nweak] = TRUE;
}
}
if ( ( only_viol ) &&
( (*best_even_slack) > MAX_SLACK - EPS ) &&
( (*best_odd_slack) > MAX_SLACK - EPS ) ) {
free(type_even_weak); free(type_odd_weak);
free(switch_even_weak); free(switch_odd_weak);
#ifdef TIME
second_(&tff);
bw_time += tff - tii;
#endif
return(NONE);
}
nweak++;
}
/* construct the weakening vectors associated with the best
even and odd pairs (if the associated slack is not too big) */
if ( (! only_odd) &&
( ( (*best_even_slack) <= MAX_SLACK - EPS ) ||
( (! only_viol) && (*best_even_slack) <= INF - EPS ) ) ) {
ok_even = TRUE;
(*info_even_weak) = alloc_info_weak(nweak);
(*info_even_weak)->nweak = nweak;
flag_even = EVEN;
cntweak = nweak;
for ( ofsl = n_to_weak - 1; ofsl >= 0; ofsl-- ) {
cntweak--;
(*info_even_weak)->var[cntweak] = vars_to_weak[ofsl];
if ( flag_even == EVEN ) {
(*info_even_weak)->type[cntweak] = type_even_weak[cntweak];
if ( switch_even_weak[cntweak] ) flag_even = ODD;
}
else {
(*info_even_weak)->type[cntweak] = type_odd_weak[cntweak];
if ( switch_odd_weak[cntweak] ) flag_even = EVEN;
}
}
}
else ok_even = FALSE;
if ( ( (*best_odd_slack) <= MAX_SLACK - EPS ) ||
( (! only_viol) && (*best_odd_slack) <= INF - EPS ) ) {
ok_odd = TRUE;
(*info_odd_weak) = alloc_info_weak(nweak);
(*info_odd_weak)->nweak = nweak;
flag_odd = ODD;
cntweak = nweak;
for ( ofsl = n_to_weak - 1; ofsl >= 0; ofsl-- ) {
cntweak--;
(*info_odd_weak)->var[cntweak] = vars_to_weak[ofsl];
if ( flag_odd == EVEN ) {
(*info_odd_weak)->type[cntweak] = type_even_weak[cntweak];
if ( switch_even_weak[cntweak] ) flag_odd = ODD;
}
else {
(*info_odd_weak)->type[cntweak] = type_odd_weak[cntweak];
if ( switch_odd_weak[cntweak] ) flag_odd = EVEN;
}
}
}
else ok_odd = FALSE;
free(type_even_weak); free(type_odd_weak);
free(switch_even_weak); free(switch_odd_weak);
#ifdef TIME
second_(&tff);
bw_time += tff - tii;
#endif
if ( ok_odd && ok_even ) return(BOTH);
if ( ok_even ) return(EVEN);
if ( ok_odd ) return(ODD);
return(NONE);
}
/* basic_separation: try to identify violated 0-1/2 cuts by using the
original procedure described in Caprara and Fischetti's MP paper */
cut_list *Cgl012Cut::basic_separation()
{
int i, j, k, l, begi, special, ofsj, ofsk, ofsl, n_to_weak, c;
short int parity, original_parity, ok_weak;
double weight, original_slack, best_even_slack, best_odd_slack;
int *vars_to_weak;
info_weak *info_even_weak, *info_odd_weak, *i_weak;
separation_graph *sep_graph;
auxiliary_graph *aux_graph;
cycle_list *short_cycle_list;
cut *violated_cut;
cut_list *out_cuts;
/* construct the separation graph by the standard weakening procedure */
#ifdef PRINT
print_parity_ilp();
#endif
#ifdef TIME
second_(&td);
#endif
#ifdef PRINT_TIME
printf("... time elapsed at the beginning of basic_separation: %f\n",td - tti);
#endif
#ifdef TIME
second_(&td);
#endif
#ifdef PRINT_TIME
printf("... time elapsed before initialize_sep_graph: %f\n",td - tti);
#endif
sep_graph = initialize_sep_graph();
special = p_ilp->mc;
#ifdef TIME
second_(&td);
#endif
#ifdef PRINT_TIME
printf("... time elapsed before weakening: %f\n",td - tti);
#endif
/* edges associated with actual constraints in the ILP */
for ( i = 0; i < p_ilp->mr; i++ ) {
if ( ! p_ilp->row_to_delete[i] ) {
begi = p_ilp->mtbeg[i];
if ( p_ilp->mtcnt[i] == 1 ) {
/* row with one odd entry only: edge j -- special */
weight = p_ilp->slack[i];
if ( weight < MAX_SLACK - EPS ) {
j = p_ilp->mtind[begi];
parity = p_ilp->mrhs[i];
i_weak = alloc_info_weak(0);
sep_graph = update_weight_sep_graph
(j,special,weight,parity,i,i_weak,sep_graph);
}
}
else if ( p_ilp->mtcnt[i] == 2 ) {
/* row with two odd entries only: edge j -- k */
weight = p_ilp->slack[i];
if ( weight < MAX_SLACK - EPS ) {
j = p_ilp->mtind[begi];
k = p_ilp->mtind[begi+1];
parity = p_ilp->mrhs[i];
i_weak = alloc_info_weak(0);
sep_graph = update_weight_sep_graph
(j,k,weight,parity,i,i_weak,sep_graph);
}
}
else {
/* row with three or more odd entries: weakening for all 1's pairs */
for ( ofsj = 0; ofsj < p_ilp->mtcnt[i]; ofsj++ ) {
for ( ofsk = ofsj + 1; ofsk < p_ilp->mtcnt[i]; ofsk++ ) {
/* edge(s) j -- k */
j = p_ilp->mtind[begi+ofsj];
k = p_ilp->mtind[begi+ofsk];
original_slack = p_ilp->slack[i];
original_parity = p_ilp->mrhs[i];
n_to_weak = 0;
vars_to_weak = reinterpret_cast<int *> (calloc(inp_ilp->mc,sizeof(int)));
if ( vars_to_weak == NULL ) alloc_error(const_cast<char*>("vars_to_weak"));
for ( ofsl = 0; ofsl < p_ilp->mtcnt[i]; ofsl++ )
if ( ofsl != ofsj && ofsl != ofsk ) {
l = p_ilp->mtind[begi+ofsl];
vars_to_weak[n_to_weak] = l;
n_to_weak++;
}
ok_weak = best_weakening(n_to_weak,vars_to_weak,
original_parity,original_slack,
&best_even_slack,&best_odd_slack,
&info_even_weak,&info_odd_weak,
FALSE,TRUE);
free(vars_to_weak);
if ( ok_weak == NONE ) goto EXITJK;
if ( ok_weak == BOTH || ok_weak == EVEN ) {
if ( best_even_slack < MAX_SLACK - EPS ) {
weight = best_even_slack; parity = EVEN;
sep_graph = update_weight_sep_graph
(j,k,weight,parity,i,info_even_weak,sep_graph);
}
}
if ( ok_weak == BOTH || ok_weak == ODD ) {
if ( best_odd_slack < MAX_SLACK - EPS ) {
weight = best_odd_slack; parity = ODD;
sep_graph = update_weight_sep_graph
(j,k,weight,parity,i,info_odd_weak,sep_graph);
}
}
EXITJK:; }
}
}
}
}
/* edges associated with the bound constraints (probably useless
but necessary in some cases */
for ( j = 0; j < p_ilp->mc; j++ ) {
if ( ! p_ilp->col_to_delete[j] ) {
weight = p_ilp->xstar[j] - inp_ilp->vlb[j];
if ( weight < MAX_SLACK - EPS ) {
parity = mod2(inp_ilp->vlb[j]);
i_weak = alloc_info_weak(0);
sep_graph = update_weight_sep_graph
(j,special,weight,parity,NONE,i_weak,sep_graph);
}
weight = inp_ilp->vub[j] - p_ilp->xstar[j];
if ( weight < MAX_SLACK - EPS ) {
parity = mod2(inp_ilp->vub[j]);
i_weak = alloc_info_weak(0);
sep_graph = update_weight_sep_graph
(j,special,weight,parity,NONE,i_weak,sep_graph);
}
}
}
#ifdef TIME
second_(&tf);
weak_time += tf - ti;
#endif
/* construct the auxiliary graph for the shortest path computation
and compute the smallest cost odd cycle visiting each node -
this part is strongly dependent on the data structure used by
the shortest path subroutine used */
#ifdef PRINT
print_sep_graph(sep_graph);
#endif
#ifdef TIME
second_(&ti);
#endif
#ifdef TIME
second_(&td);
#endif
#ifdef PRINT_TIME
printf("... time elapsed before define_aux_graph: %f\n",td - tti);
#endif
aux_graph = define_aux_graph(sep_graph);
#ifdef TIME
second_(&tf);
aux_time += tf - ti;
#endif
#ifdef PRINT
print_aux_graph(aux_graph);
#endif
/* exit(1); */
#ifdef TIME
second_(&td);
#endif
#ifdef PRINT_TIME
printf("... time elapsed before cycles and cuts: %f\n",td - tti);
printf("%d nodes on list\n",sep_graph->nnodes);
#endif
out_cuts = initialize_cut_list(MAX_CUTS);
for ( j = 0; j < sep_graph->nnodes; j++ ) {
short_cycle_list = get_shortest_odd_cycle_list(j,sep_graph,aux_graph);
if ( short_cycle_list == NULL ) goto EXIT_NODE;
#ifdef PRINT
print_cycle_list(short_cycle_list);
#endif
for ( c = 0; c < short_cycle_list->cnum; c++ ) {
violated_cut = get_cut(short_cycle_list->list[c]);
if ( violated_cut == NULL ) {
if (!errorNo)
continue;
else
break;
}
#ifdef PRINT
print_cut(violated_cut);
#endif
if ( violated_cut->violation > MIN_VIOLATION + EPS ) {
/* violated 0-1/2 cut found */
out_cuts = add_cut_to_list(violated_cut,out_cuts);
if ( out_cuts->cnum >= MAX_CUTS ) {
free_cycle_list(short_cycle_list);
goto EXIT_CUTS;
}
}
else free_cut(violated_cut);
}
/* remove the current node from the auxiliary graph */
EXIT_NODE:
aux_graph = cancel_node_aux_graph(j,aux_graph);
free_cycle_list(short_cycle_list);
}
EXIT_CUTS:
free_sep_graph(sep_graph);
free_aux_graph(aux_graph);
#ifdef PRINT_CUTS
print_cut_list(out_cuts);
#endif
#ifdef TIME
second_(&td);
#endif
#ifdef PRINT_TIME
printf("... time elapsed at the end of basic_separation: %f\n",td - tti);
#endif
return(out_cuts);
}
/*
012cut: main procedure for 0-1/2 cut separation
first release: Aug 12 1996
last revision: Jun 10 1997
*/
/* static data structures for log information about separation */
#ifndef CGGGGG
static int sep_iter = 0; /* number of the current separation iteration */
//#define POOL
#ifdef POOL
static pool_cut_list *pool = NULL; /* information about the cuts separated
so far, used to decide when they should
be added to the current LP */
#endif
static log_var **vlog = NULL; /* information about the value attained
by the variables in the last iterations,
used to possibly set to 0 some coefficient
> 0 in a cut to be added */
static bool aggr; /* flag saying whether as many cuts as possible are required
from the separation procedure (TRUE) or not (FALSE) */
#endif
/* include the reactive local search heuristic */
//was #include "Cgltabu_012.c"
//start include "Cgltabu_012.c"
#define MAX_TABU_ITER 100
#define NUM_HASH_ENTRIES MAX_CUT_POOL
/* initial length of the tabu list */
#define IN_PROHIB_PERIOD 3
#define MAX_TIME_FACTOR 3
/* data structure for the current local search solution */
typedef struct {
int n_of_constr; /* number of constraints in the current cut */
short int *in_constr_list; /* flag saying whether a given constraint is
in the list of constraints of the cut (IN)
or not (OUT) */
int *non_weak_coef; /* coefficients of the cut before weakening */
int non_weak_rhs; /* coefficient of the rhs before weakening */
double slack_sum; /* sum of the slacks of the constraints in the cut */
double min_weak_loss; /* minimum loss by weakening the non even
coefficients */
int one_norm; /* 1-norm of the lhs, i.e. sum of the absolute values of
the coefficients */
short int ok; /* logical flag telling whether the cut could be weakened
to a 0-1/2 cut or not - if false the two fields below
have no meaning */
int *coef; /* actual coefficients of the cut */
int rhs; /* actual rhs of the cut */
double violation; /* violation of the cut */
} tabu_cut;
/* data structure for the hash table used in memory reaction */
typedef struct h_e {
int n_of_el; /* number of components to be considered */
short int *flag_vect; /* vector of flags for the components */
int last_vis; /* last iteration when this element was visited */
struct h_e *next; /* pointer to the next element in the hash chain */
} hash_element;
typedef hash_element **hash_table;
/* global variables for local search */
int n; /* number of variables in the ILP */
int m; /* number of constraints in the ILP */
int it; /* number of tabu search iterations so far */
tabu_cut *cur_cut; /* information about the current cut in local search */
int *last_moved; /* last iteration when a given constraint was added/
deleted from the list of constraints of the cut */
int last_it_add; /* last iteration when a cut was added to the list */
int last_it_restart; /* last iteration when a restart was performed */
int prohib_period; /* current prohibition period */
int last_prohib_period_mod; /* last iteration where prohibition period was modified */
hash_table hash_tab; /* hash table */
int A; /* parameter A in Battiti and Protasi */
int B; /* parameter B in Battiti and Protasi */
float elapsed_time; /* time elapsed since the beginning of the current
tabu search call */
/* clear_cur_cut: clear the current solution (no constraint in the cut) */
void clear_cur_cut()
{
int i, j;
cur_cut->n_of_constr = 0;
cur_cut->rhs = 0;
cur_cut->non_weak_rhs = 0;
cur_cut->violation = 0.0;
cur_cut->slack_sum = 0.0;
cur_cut->min_weak_loss = 0.0;
cur_cut->one_norm = 0;
for ( j = 0; j < n; j++ ) {
cur_cut->coef[j] = 0;
cur_cut->non_weak_coef[j] = 0;
}
for ( i = 0; i < m; i++ ) {
cur_cut->in_constr_list[i] = OUT;
}
cur_cut->ok = FALSE;
}
/* initialize_cur_cut: allocate the memory for cur_cut */
void initialize_cur_cut()
{
cur_cut = reinterpret_cast<tabu_cut *> (calloc(1,sizeof(tabu_cut)));
if ( cur_cut == NULL ) alloc_error(const_cast<char*>("cur_cut"));
cur_cut->coef = reinterpret_cast<int *> (calloc(n,sizeof(int)));
if ( cur_cut->coef == NULL ) alloc_error(const_cast<char*>("cur_cut->coef"));
cur_cut->non_weak_coef = reinterpret_cast<int *> (calloc(n,sizeof(int)));
if ( cur_cut->non_weak_coef == NULL ) alloc_error(const_cast<char*>("cur_cut->non_weak_coef"));
cur_cut->in_constr_list = reinterpret_cast<short int *> (calloc(m,sizeof(short int)));
if ( cur_cut->in_constr_list == NULL ) alloc_error(const_cast<char*>("cur_cut->in_constr_list"));
clear_cur_cut();
}
/* free_cur_cut: free the memory for cur_cut */
void free_cur_cut()
{
free(cur_cut->coef);
free(cur_cut->non_weak_coef);
free(cur_cut->in_constr_list);
free(cur_cut);
}
#ifdef PRINT_TABU
/* print_cur_cut: display cur_cut on output */
void Cgl012Cut::print_cur_cut()
{
int i, j;
printf("iteration %d prohib_period %d\n",it,prohib_period);
printf("\n content of cur_cut data structure: n_of_constr = %d, ok = %d\n", cur_cut->n_of_constr, cur_cut->ok);
for ( i = 0; i < m; i++ )
if ( cur_cut->in_constr_list[i] == IN )
printf("constr. %d\n",i);
/*
printf(" list of constraints:\n");
for ( i = 0; i < m; i++ )
if ( cur_cut->in_constr_list[i] == IN )
print_constr(i);
print_int_vect("non_weak_coef",cur_cut->non_weak_coef,n);
printf(" non_weak_rhs = %d\n",cur_cut->non_weak_rhs);
print_int_vect("coef",cur_cut->coef,n);
printf(" rhs = %d\n",cur_cut->rhs);
*/
for ( j = 0 /* , viol = - (double) cur_cut->rhs */ ; j < n; j++ )
if ( ( p_ilp->xstar[j] > ZERO || p_ilp->xstar[j] < -ZERO ) && cur_cut->non_weak_coef[j] != 0 ) {
printf("var. %d xstar %f non_weak_coef %d coef %d\n", j, p_ilp->xstar[j], cur_cut->non_weak_coef[j], cur_cut->coef[j]);
/* viol += p_ilp->xstar[j] * cur_cut->coef[j]; */
}
printf("rhs %d viol %f slack_sum %f min_weak_loss %f one_norm %d\n",
cur_cut->rhs, cur_cut->violation, cur_cut->slack_sum,
cur_cut->min_weak_loss, cur_cut->one_norm);
}
#endif
/* same_short_vect: check whether two short int vectors have the same content */
short int same_short_vect(
int n_of_el, /* number of components in the vectors */
short int *vec_1,
short int *vec_2 /* vectors to be checked */
)
{
int i;
for ( i = 0; i < n_of_el; i++ )
if ( vec_1[i] != vec_2[i] ) return(FALSE);
return(TRUE);
}
/* initialize_hash_table: allocate the memory for the hash table */
void initialize_hash_table()
{
int i;
hash_tab = reinterpret_cast<hash_element **> (calloc(NUM_HASH_ENTRIES,sizeof(hash_element *)));
if ( hash_tab == NULL ) alloc_error(const_cast<char*>("hash_tab"));
for ( i = 0; i < NUM_HASH_ENTRIES; i++ ) hash_tab[i] = NULL;
}
/* clear_hash_table: clear the current hash table */
void clear_hash_table()
{
int i;
hash_element *hash_ptr, *hash_el;
for ( i = 0; i < NUM_HASH_ENTRIES; i++ ) {
if ( hash_tab[i] != NULL ) {
hash_ptr = hash_tab[i];
do {
hash_el = hash_ptr->next;
free(hash_ptr->flag_vect);
free(hash_ptr);
hash_ptr = hash_el;
} while ( hash_ptr != NULL );
hash_tab[i] = NULL;
}
}
}
/* free_hash_table: deallocate the memory for the hash table */
void free_hash_table()
{
clear_hash_table();
free(hash_tab);
}
/* hash_addr: compute the hash address associated with the current cut */
int hash_addr(
int n_of_el, /* number of elements to be considered */
short int *flag_vect /* vector of flags for the elements */
)
{
int i, addr;
/* very simple algorithm: just add-up the squared indices of the IN elements */
addr = 0;
for ( i = 0; i < n_of_el; i++ )
if ( flag_vect[i] == IN ) addr += i * i;
return(addr % NUM_HASH_ENTRIES);
}
/* hash_search: search for the current cut in the hash list of all cuts -
if found return TRUE and update the last iteration the cut was found */
short int hash_search(int *cyc_len /* length of the cycle if the current cut is found */)
{
int addr;
hash_element *hash_el;
addr = hash_addr(m,cur_cut->in_constr_list);
hash_el = hash_tab[addr];
while ( hash_el != NULL ) {
if ( same_short_vect(m,cur_cut->in_constr_list,hash_el->flag_vect) ) {
*cyc_len = it - hash_el->last_vis;
hash_el->last_vis = it;
return(TRUE);
}
hash_el = hash_el->next;
}
return(FALSE);
}
/* hash_insert: insert a new cut in the hash list of all cuts */
void hash_insert()
{
int addr, i;
hash_element *hash_el, *hash_ptr;
addr = hash_addr(m,cur_cut->in_constr_list);
hash_el = reinterpret_cast<hash_element *> (calloc(1,sizeof(hash_element)));
if ( hash_el == NULL ) alloc_error(const_cast<char*>("hash_el"));
hash_el->n_of_el = m;
hash_el->last_vis = it;
hash_el->next = NULL;
hash_el->flag_vect = reinterpret_cast<short int *> (calloc(m,sizeof(short int)));
if ( hash_el->flag_vect == NULL ) alloc_error(const_cast<char*>("hash_el->flag_vect"));
for ( i = 0; i < m; i++ )
hash_el->flag_vect[i] = cur_cut->in_constr_list[i];
if ( hash_tab[addr] == NULL )
hash_tab[addr] = hash_el;
else {
hash_ptr = hash_tab[addr];
while ( hash_ptr->next != NULL ) {
#if 0
/* this check can be omitted to save time */
if ( same_short_vect(m,cur_cut->in_constr_list,hash_ptr->flag_vect) ) {
printf("attempt to insert in the hash an already present cut\n");
exit(0);
}
#endif
hash_ptr = hash_ptr->next;
}
hash_ptr->next = hash_el;
}
}
/* increase_prohib_period: implemented as in Battiti and Protasi */
void increase_prohib_period()
{
if ( prohib_period * 1.1 > prohib_period + 1 )
if ( prohib_period * 1.1 < m - 2 ) prohib_period =
static_cast<int> (prohib_period*1.1);
else prohib_period = m - 2;
else
if ( prohib_period + 1 < m - 2 ) prohib_period += 1;
else prohib_period = m - 2;
last_prohib_period_mod = it;
}
/* decrease_prohib_period: implemented as in Battiti and Protasi */
void decrease_prohib_period()
{
if ( prohib_period * 0.9 < prohib_period - 1 )
if ( prohib_period * 0.9 > IN_PROHIB_PERIOD ) prohib_period = static_cast<int> (prohib_period* 0.9);
else prohib_period = IN_PROHIB_PERIOD;
else
if ( prohib_period - 1 > IN_PROHIB_PERIOD ) prohib_period -= 1;
else prohib_period = IN_PROHIB_PERIOD;
last_prohib_period_mod = it;
}
/* allowed: check if moving (adding/deleting) a given constraint
is not a tabu move */
short int allowed(int i /* constraint to be checked */)
{
if ( last_moved[i] < it - prohib_period ) {
if ( cur_cut->in_constr_list[i] == IN ) {
if ( cur_cut->n_of_constr > 1 ) return(TRUE);
else return(FALSE);
}
else {
if ( cur_cut->n_of_constr < m - 1 ) return(TRUE);
else return(FALSE);
}
}
else return(FALSE);
}
/* in_cur_cut: check whether a given constraint is in the list of
constraints defining the current cut */
short int in_cur_cut(int i /* constraint to be checked */)
{
if ( cur_cut->in_constr_list[i] == OUT ) return(FALSE);
else return(TRUE);
}
/* tabu_score: define the score of a potential new cut */
double tabu_score(
int *ccoef, /* cut left hand side coefficients */
int crhs, /* cut right hand side */
double viol, /* cut violation */
double norm /* cut norm - 1-norm is used below */
)
{
/* very simple score: violation divided/multiplied by the lhs 1-norm */
if ( norm == 0 ) norm = 1;
if ( viol > 0.0 ) return (viol / norm);
else return (viol * norm);
}
/* score_by_moving: compute the score of the best cut obtainable from
the current local search solution by inserting/deleting a constraint */
double Cgl012Cut::score_by_moving(
int i, /* constraint to be moved */
short int itype, /* type of move - ADD or DEL */
double thresh /* minimum value of an interesting score */
)
{
#define PENALTY_NON_WEAKABLE -1.0
#define FAST_SCORE_EVAL 1
int j, begi, gcdi, ofsj, ij, crhs, one_norm, support_inter;
short int flag_gt, ok;
double slack_sum, weak_loss, score, viol;
int *new_coef, *ccoef;
begi = inp_ilp->mtbeg[i];
gcdi = p_ilp->gcd[i];
/* fast check - optimistic evaluation of the score */
slack_sum = cur_cut->slack_sum;
if ( itype == ADD ) slack_sum += p_ilp->slack[i] / static_cast<double> (gcdi);
else slack_sum -= p_ilp->slack[i] / static_cast<double> (gcdi);
viol = ( 1.0 - slack_sum ) / 2.0;
score = tabu_score(NULL,0,viol,1.0);
/*
printf("Score estimate 1 %f Threshold %f\n",score,thresh);
*/
if ( score < thresh + ZERO ) {
return (score);
}
/* discard the cuts that have empty support intersection with the
current one */
support_inter = 0;
for ( ofsj = 0, ij = begi; ofsj < inp_ilp->mtcnt[i]; ofsj++, ij++ )
if ( cur_cut->non_weak_coef[inp_ilp->mtind[ij]] !=0 ) support_inter++;
if ( support_inter == 0 ) return(-INF);
/* compute the new cut coefficients and rhs (before weakening) */
new_coef = reinterpret_cast<int *> (calloc(inp_ilp->mtcnt[i],sizeof(int)));
if ( new_coef == NULL ) alloc_error(const_cast<char*>("new_coef"));
if ( ( itype == ADD && inp_ilp->msense[i] != 'G' ) ||
( itype == DEL && inp_ilp->msense[i] == 'G' ) ) flag_gt = 1;
else flag_gt = -1;
if ( flag_gt == 1 ) {
if ( gcdi == 1 ) {
for ( ofsj = 0, ij = begi; ofsj < inp_ilp->mtcnt[i]; ofsj++, ij++ ) {
j = inp_ilp->mtind[ij];
new_coef[ofsj] = cur_cut->non_weak_coef[j] + inp_ilp->mtval[ij];
}
crhs = cur_cut->non_weak_rhs + inp_ilp->mrhs[i];
}
else {
for ( ofsj = 0, ij = begi; ofsj < inp_ilp->mtcnt[i]; ofsj++, ij++ ) {
j = inp_ilp->mtind[ij];
new_coef[ofsj] = cur_cut->non_weak_coef[j] + inp_ilp->mtval[ij] / gcdi;
}
crhs = cur_cut->non_weak_rhs + inp_ilp->mrhs[i] / gcdi;
}
}
else {
if ( gcdi == 1 ) {
for ( ofsj = 0, ij= begi; ofsj < inp_ilp->mtcnt[i]; ofsj++, ij++ ) {
j = inp_ilp->mtind[ij];
new_coef[ofsj] = cur_cut->non_weak_coef[j] - inp_ilp->mtval[ij];
}
crhs = cur_cut->non_weak_rhs - inp_ilp->mrhs[i];
}
else {
for ( ofsj = 0, ij = begi; ofsj < inp_ilp->mtcnt[i]; ofsj++, ij++ ) {
j = inp_ilp->mtind[ij];
new_coef[ofsj] = cur_cut->non_weak_coef[j] - inp_ilp->mtval[ij] / gcdi;
}
crhs = cur_cut->non_weak_rhs - inp_ilp->mrhs[i] / gcdi;
}
}
/* other - relatively fast - check by optimistic evaluation of the
cut score */
weak_loss = cur_cut->min_weak_loss;
one_norm = cur_cut->one_norm;
for ( ofsj = 0, ij = begi; ofsj < inp_ilp->mtcnt[i]; ofsj++, ij++ ) {
j = inp_ilp->mtind[ij];
if ( cur_cut->coef[j] > 0 ) one_norm -= cur_cut->coef[j];
else one_norm += cur_cut->coef[j];
if ( new_coef[ofsj] >= 2 ) one_norm += new_coef[ofsj] / 2;
else one_norm -= new_coef[ofsj] / 2;
if ( mod2(cur_cut->non_weak_coef[j]) == ODD ) {
if ( mod2(new_coef[ofsj]) == EVEN )
weak_loss -= p_ilp->min_loss_by_weak[j];
}
else {
if ( mod2(new_coef[ofsj]) == ODD )
weak_loss += p_ilp->min_loss_by_weak[j];
}
}
viol = ( 1.0 - slack_sum - weak_loss ) / 2.0;
score = tabu_score(NULL,0,viol,static_cast<double> (one_norm));
/*
printf("Score estimate 2 %f Threshold %f\n",score,thresh);
*/
if ( score < thresh + ZERO || FAST_SCORE_EVAL ) {
free(new_coef);
return (score);
}
/* get the actual cut coefficients and the violation of the
best cut obtainable trough weakening */
ccoef = reinterpret_cast<int *> (calloc(n,sizeof(int)));
if ( ccoef == NULL ) alloc_error(const_cast<char*>("ccoef"));
for ( j = 0; j < n; j++ ) ccoef[j] = cur_cut->non_weak_coef[j];
for ( ofsj = 0; ofsj < inp_ilp->mtcnt[i]; ofsj++ ) {
ij = begi + ofsj;
j = inp_ilp->mtind[ij];
ccoef[j] = new_coef[ofsj];
}
ok = best_cut(ccoef,&crhs,&viol,FALSE,FALSE);
if ( ok ) score = tabu_score(ccoef,crhs,viol,static_cast<double> (one_norm));
else {
viol = ( - slack_sum - weak_loss ) / 2.0;
score = PENALTY_NON_WEAKABLE + tabu_score(ccoef,crhs,viol,static_cast<double> (one_norm));
}
/*
printf("Score actual %f Threshold %f\n",score,thresh);
*/
free(new_coef);
free(ccoef);
return(score);
}
/* modify_current: update the current local search solution by inserting/
deleting a constraint */
void Cgl012Cut::modify_current(
int i, /* constraint to be moved */
short int itype /* type of move - ADD or DEL */
)
{
int j, begi, gcdi, ofsj, ij;
short int flag_gt;
if ( itype == ADD ) {
cur_cut->n_of_constr++;
cur_cut->in_constr_list[i] = IN;
}
else {
cur_cut->n_of_constr--;
cur_cut->in_constr_list[i] = OUT;
}
last_moved[i] = it;
/* compute the new cut coefficients and rhs (before weakening) */
if ( ( itype == ADD && inp_ilp->msense[i] != 'G' ) ||
( itype == DEL && inp_ilp->msense[i] == 'G' ) ) flag_gt = 1;
else flag_gt = -1;
begi = inp_ilp->mtbeg[i];
gcdi = p_ilp->gcd[i];
for ( ofsj = 0; ofsj < inp_ilp->mtcnt[i]; ofsj++ ) {
ij = begi + ofsj;
j = inp_ilp->mtind[ij];
/* the '*' and '/' operations can be saved by writing some more code ... */
cur_cut->non_weak_coef[j] += flag_gt * (inp_ilp->mtval[ij] / gcdi);
}
cur_cut->non_weak_rhs += flag_gt * (inp_ilp->mrhs[i] / gcdi);
if ( itype == ADD )
cur_cut->slack_sum += p_ilp->slack[i] / static_cast<double> (gcdi);
else
cur_cut->slack_sum -= p_ilp->slack[i] / static_cast<double> (gcdi);
/* get the best possible cut */
cur_cut->min_weak_loss = 0.0;
for ( j = 0; j < n; j++ ) {
cur_cut->coef[j] = cur_cut->non_weak_coef[j];
if ( mod2(cur_cut->coef[j]) == ODD )
cur_cut->min_weak_loss += p_ilp->min_loss_by_weak[j];
}
cur_cut->rhs = cur_cut->non_weak_rhs;
cur_cut->ok =
best_cut(cur_cut->coef,&cur_cut->rhs,&cur_cut->violation,TRUE,FALSE);
cur_cut->one_norm = 0;
for ( j = 0; j < n; j++ ) {
if ( cur_cut->coef[j] > 0 ) cur_cut->one_norm += cur_cut->coef[j];
else cur_cut->one_norm -= cur_cut->coef[j];
}
}
/* get_current_cut: return a cut data type with the information about
the current cut of the search procedure */
cut *Cgl012Cut::get_current_cut()
{
int i, j, nz;
/*double viol;*/
cut *cut_ptr;
cut_ptr = reinterpret_cast<cut *> (calloc(1,sizeof(cut)));
if ( cut_ptr == NULL ) alloc_error(const_cast<char*>("cut_ptr"));
cut_ptr->crhs = cur_cut->rhs;
cut_ptr->csense = 'L';
/* count the number of nonzeroes in the cut */
for ( j = 0, nz = 0; j < n; j++ ) if ( cur_cut->coef[j] != 0 ) nz++;
cut_ptr->cnzcnt = nz;
cut_ptr->cind = reinterpret_cast<int *> (calloc(nz,sizeof(int)));
if ( cut_ptr->cind == NULL ) alloc_error(const_cast<char*>("cut_ptr->cind"));
cut_ptr->cval = reinterpret_cast<int *> (calloc(nz,sizeof(int)));
if ( cut_ptr->cval == NULL ) alloc_error(const_cast<char*>("cut_ptr->cval"));
nz = 0; /*viol = 0.0;*/
for ( j = 0; j < n; j++ ) {
if ( cur_cut->coef[j] != 0 ) {
cut_ptr->cind[nz] = j;
cut_ptr->cval[nz] = cur_cut->coef[j];
nz++;
/* viol += p_ilp->xstar[j] * (double) cur_cut->coef[j]; */
}
}
/* viol -= (double) cur_cut->rhs; */
/* cut_ptr->violation = viol; */
cut_ptr->violation = cur_cut->violation;
cut_ptr->n_of_constr = 0;
cut_ptr->constr_list = reinterpret_cast<int *> (calloc(inp_ilp->mr,sizeof(int)));
if ( cut_ptr->constr_list == NULL ) alloc_error(const_cast<char*>("cut_ptr->constr_list"));
cut_ptr->in_constr_list = reinterpret_cast<short int *> (calloc(inp_ilp->mr,sizeof(short int)));
if ( cut_ptr->in_constr_list == NULL ) alloc_error(const_cast<char*>("cut_ptr->in_constr_list"));
for ( i = 0; i < m; i++ ) {
if ( cur_cut->in_constr_list[i] == IN ) {
cut_ptr->in_constr_list[i] = IN;
cut_ptr->constr_list[cut_ptr->n_of_constr] = i;
(cut_ptr->n_of_constr)++;
}
else cut_ptr->in_constr_list[i] = OUT;
}
return(cut_ptr);
}
/* best neighbour: find the cut to be added/deleted from the current
solution among those allowed by the tabu rules */
short int Cgl012Cut::best_neighbour(cut_list *out_cuts /* list of the violated cuts found */)
{
int i, ibest;
short int itype, itypebest=-1;
double score, max_score;
cut *new_cut;
/* cycle through all the constraints in your problem ... */
max_score = -INF;
ibest = NONE;
for ( i = 0; i < m; i++ ) {
if ( ! p_ilp->row_to_delete[i] && allowed(i) ) {
if ( in_cur_cut(i) ) {
/* constraint i is in the current set of constraints
(those defining the current cut) */
itype = DEL;
}
else {
/* constraint i is not in the current set of constraints
(those defining the current cut) */
itype = ADD;
}
score = score_by_moving(i,itype,max_score);
if ( score > max_score ) {
/* best cut found in this iteration: store it */
ibest = i;
itypebest = itype;
max_score = score;
}
}
} /* for ( i = 0; i < m; i++ ) */
if ( ibest == NONE ) {
#ifdef PRINT_TABU
printf("No move could be performed by best_neighbour\n");
#endif
return(TRUE);
}
modify_current(ibest,itypebest);
if ( cur_cut->violation > MIN_VIOLATION + EPS ) {
#ifdef PRINT_TABU
printf("... adding the current cut to the output list - it = %d viol = %f\n",it, cur_cut->violation);
#endif
new_cut = get_current_cut();
out_cuts = add_cut_to_list(new_cut,out_cuts);
last_it_add = it;
}
return(FALSE);
}
/* memory_reaction: perform the long term reaction by cheching whether the
current solution has already been visited or the best solution has not
been updated for too many iterations */
void memory_reaction()
{
int cycle_length;
if ( hash_search(&cycle_length) ) {
if ( cycle_length < 2 * ( m - 1 ) ) {
increase_prohib_period();
return;
}
}
else hash_insert();
if ( it - last_prohib_period_mod > B )
decrease_prohib_period();
}
/* add_tight_constraint: initialize the current cut by adding a tight
constraint to it */
void Cgl012Cut::add_tight_constraint()
{
int i, ntight;
double smin=COIN_DBL_MAX;
abort();
int *tight;
ntight = 0;
tight = reinterpret_cast<int *> (calloc(m,sizeof(int)));
if ( tight == NULL ) alloc_error(const_cast<char*>("tight"));
for ( i = 0; i < m; i++ ) {
/* search for the tightest constraint never added to cut */
if ( last_moved[i] < 0 && p_ilp->slack[i] < smin ) {
if ( p_ilp->slack[i] < ZERO ) {
/* tight constraint */
smin = ZERO;
tight[ntight] = i;
ntight++;
}
else {
/* best constraint so far */
smin = p_ilp->slack[i];
tight[0] = i;
ntight = 1;
}
}
}
if ( ntight > 0 ) i = tight[rand() % ntight];
/* if all constraints have already been in cur_cut choose first at random */
else i = rand() % m;
free(tight);
modify_current(i,ADD);
}
/* initialize: initialize the data structures for local search */
void Cgl012Cut::initialize()
{
int i;
m = inp_ilp->mr;
n = inp_ilp->mc;
it = 0;
last_it_add = 0;
last_it_restart = 0;
last_prohib_period_mod = 0;
prohib_period = IN_PROHIB_PERIOD;
initialize_cur_cut();
last_moved = reinterpret_cast<int *> (calloc(m,sizeof(int)));
if ( last_moved == NULL ) alloc_error(const_cast<char*>("last_moved"));
for ( i = 0; i < m; i++ ) {
last_moved[i] = -COIN_INT_MAX;
}
initialize_hash_table();
add_tight_constraint();
A = m;
B = 10 * m;
}
/* restart: perform a restart of the search - IMPORTANT: in the current
implementation vector last_moved is not cleared at restart */
void Cgl012Cut::restart(short int failure /* flag forcing the restart if some trouble occurred */)
{
if ( failure || ( it - last_it_add > A && it - last_it_restart > A ) ) {
/* perform restart */
last_it_restart = it;
prohib_period = IN_PROHIB_PERIOD;
last_prohib_period_mod = it;
clear_hash_table();
clear_cur_cut();
add_tight_constraint();
}
}
/* free_memory: free the memory used by local search */
void free_memory()
{
free_cur_cut();
free(last_moved);
free_hash_table();
}
/* tabu_012: try to identify violated 0-1/2 cuts by a simple tabu search
procedure adapted from that used by Battiti and Protasi for finding
large cliques */
cut_list *Cgl012Cut::tabu_012()
{
short int failure;
cut_list *out_cuts;
out_cuts = initialize_cut_list(MAX_CUTS);
initialize();
it = 0;
do {
memory_reaction();
failure = best_neighbour(out_cuts);
#ifdef PRINT_TABU
print_cur_cut();
#endif
it++;
restart(failure);
}
while ( out_cuts->cnum < MAX_CUTS && it < MAX_TABU_ITER );
free_memory();
#ifdef PRINT_TABU
printf("Number of violated cuts found by Tabu Search %d\n",out_cuts->cnum);
printf("Tabu Search timings: best_neighbour %f score_by_moving %f coefficient %f best_cut %f\n",
time_best_neigh, time_scor_by_mov, time_coef, time_best_cut);
#endif
return(out_cuts);
}
//end include "Cgltabu_012.c"
#ifdef POOL
/* same_pool_cut: check whether two pool cuts are in fact the same cut */
short int same_pool_cut(pool_cut *p_cut1, pool_cut *p_cut2)
{
int c;
if ( p_cut1->n_of_constr != p_cut2->n_of_constr ) return(FALSE);
if ( p_cut1->code != p_cut2->code ) return(FALSE);
/* assumes constr_list is sorted for increasing/decreasing constraint index */
for ( c = 1; c < p_cut1->n_of_constr; c++ )
if ( p_cut1->constr_list[c] != p_cut2->constr_list[c] ) return(FALSE);
return(TRUE);
}
/* free_pool_cut: free the memory for a non-empty pool cut */
void free_pool_cut(pool_cut *p_cut)
{
if ( p_cut == NULL ) return;
if ( p_cut->constr_list != NULL ) free(p_cut->constr_list);
free(p_cut);
}
/* initialize_pool: initialize the pool data structure */
void initialize_pool()
{
pool = (pool_cut_list *) calloc(1,sizeof(pool_cut_list));
if ( pool == NULL ) alloc_error(const_cast<char*>("pool"));
pool->cnum = 0;
pool->list = (pool_cut **) calloc(MAX_CUT_POOL,sizeof(pool_cut *));
if ( pool->list == NULL ) alloc_error(const_cast<char*>("pool->list"));
pool->ncod = (int *) calloc(MAX_CUT_COD,sizeof(int));
if ( pool->ncod == NULL ) alloc_error(const_cast<char*>("pool->ncod"));
}
/* free_pool: free the memory used by the pool */
void free_pool()
{
int c;
if ( pool == NULL ) return;
for ( c = 0; c < pool->cnum; c++ )
free_pool_cut(pool->list[c]);
free(pool);
}
/* clean_pool: remove form the pool the cuts which are inactive since a
large number of iterations */
void clean_pool()
{
int c, d;
/* the pool compression could be implemented more efficiently */
for ( c = 0; c < pool->cnum; c++ ) {
if ( pool->list[c]->it_found > MAX_ITER_POOL ) {
free_pool_cut(pool->list[c]);
pool->list[c] = NULL;
for ( d = c ; d < pool->cnum; d++ )
pool->list[d] = pool->list[d + 1];
(pool->cnum)--;
}
}
}
/* insert_cut_in_pool: add a cut to the pool if there is space */
void insert_cut_in_pool(pool_cut *p_cut)
{
if ( pool->cnum == MAX_CUT_POOL ) {
#ifdef COIN_DEVELOP
printf("Warning: pool is full and separated cuts cannot be added\n");
#endif
return;
}
pool->list[pool->cnum] = p_cut;
(pool->cnum)++;
(pool->ncod[p_cut->code])++;
}
/* cut_is_in_pool: check whether a given cut is already in the pool */
short int cut_is_in_pool(cut *v_cut)
{
int c, i, cod;
short int equal;
short int *flag_v;
pool_cut *p_cut;
/*
print_cut(v_cut);
printf("checking for a cut in the pool ...\n");
*/
cod = hash_addr(inp_ilp->mr,v_cut->in_constr_list);
if ( pool->ncod[cod] == 0 ) return(FALSE);
/* trivial sequential search */
/*
printf("... sequential search needed ...\n");
*/
flag_v = (short int *) calloc(inp_ilp->mr,sizeof(short int));
if ( flag_v == NULL ) alloc_error(const_cast<char*>("flag_v"));
for ( c = 0; c < pool->cnum; c++ ) {
if ( pool->list[c]->code != cod ) continue;
p_cut = pool->list[c];
equal = TRUE;
for ( i = 0; i < inp_ilp->mr; i++ ) flag_v[i] = OUT;
for ( i = 0; i < p_cut->n_of_constr; i++ )
flag_v[p_cut->constr_list[i]] = IN;
for ( i = 0; i < inp_ilp->mr; i++ ) {
if ( v_cut->in_constr_list[i] != flag_v[i] ) {
equal = FALSE;
break;
}
}
if ( equal ) {
free(flag_v);
/*
printf("... cut was in the pool!\n");
*/
return(TRUE);
}
}
return(FALSE);
}
/* good_pool_cut: check whether the current cut is worth adding to the pool */
short int good_pool_cut(cut *v_cut)
{
/* no check performed */
return(TRUE);
}
/* add_cuts_to_pool: add the cuts separated to the pool structure */
void add_cuts_to_pool(cut_list *out_cuts)
{
int i, c;
cut *v_cut;
pool_cut *p_cut;
/* float ti, tf, tti, ttf; */
/* static float tcutis = 0.0, taddcut = 0.0; */
/* second_(&tti); */
if ( pool == NULL ) initialize_pool();
if ( pool->cnum >= MAX_CUT_POOL * CLEAN_THRESH ) clean_pool();
for ( i = 0; i < out_cuts->cnum; i++ ) {
v_cut = out_cuts->list[i];
/* second_(&ti); */
if ( cut_is_in_pool(v_cut) ) {
/* second_(&tf); */
/* tcutis += tf - ti; */
continue;
}
else {
/* second_(&tf); */
/* tcutis += tf - ti; */
}
if ( good_pool_cut(v_cut) ) {
/* add the cut to the pool list */
p_cut = (pool_cut *) calloc(1,sizeof(pool_cut));
if ( p_cut == NULL ) alloc_error(const_cast<char*>("p_cut"));
p_cut->n_of_constr = v_cut->n_of_constr;
p_cut->constr_list = (int *) calloc(v_cut->n_of_constr,sizeof(int));
if ( p_cut->constr_list == NULL ) alloc_error(const_cast<char*>("p_cut->constr_list"));
for ( c = 0; c < v_cut->n_of_constr; c++ )
p_cut->constr_list[c] = v_cut->constr_list[c];
p_cut->code = hash_addr(inp_ilp->mr,v_cut->in_constr_list);
p_cut->n_it_violated = 0;
p_cut->it_found = sep_iter;
insert_cut_in_pool(p_cut);
}
}
/* second_(&ttf); */
/* taddcut += ttf - tti; */
/* printf("add_cuts_to_pool: tcutis %f taddcut %f\n",tcutis,taddcut); */
}
/* interesting_var: decides whether a variable is relevant in the
separation or not */
short int interesting_var(int j /* variable to be evaluated */)
{
/* return ( vlog[j]->n_it_zero < MANY_IT_ZERO ); */
/* if ( aggr ) return (TRUE); */
return( ! p_ilp->col_to_delete[j] );
}
/* get_cuts_from_pool: select from the pool a convenient set of violated
constraints to be added to the current LP */
static double max_score_ever = ZERO; /* maximum score of a violated cut during
the whole cutting plane procedure */
#define MIN_CUT_SCORE ( max_score_ever / MIN_SCORE_RANGE )
#define MAX_CUT_SCORE ( max_score_ever / MAX_SCORE_RANGE )
#define MIN_IT_VIOL 2
cut_list *get_cuts_from_pool(
short int after_sep /* flag telling whether the pool is searched after
a new separation in which case only new cuts are
checked */
)
{
int c, crhs, j, k, l, maxc, n_interest_var;
double viol, score, maxscore, min_cut_score, max_cut_score;
short int ok;
int *interest_var, *ccoef;
short int *added;
pool_cut *p_cut;
select_cut **best_var_cut;
cut *a_cut;
cut_list *add_cuts;
/* float ti, tf, tti, ttf; */
/* static float tgetcut = 0.0, talloc = 0.0, tgetcoef = 0.0, tbestcut = 0.0, tscore = 0.0, tupdbest = 0.0, tadd = 0.0; */
/* second_(&tti); */
if ( pool == NULL ) {
add_cuts = initialize_cut_list(1);
return(add_cuts);
}
/* in the current implementation, the cut with best score with nonzero
coefficient for each "interesting" variable is added to the current LP,
provided the cut satisfies some requirements (violation, depth, etc.) */
/* define the set of the interesting variables */
interest_var = (int *) calloc(p_ilp->mc,sizeof(int));
if ( interest_var == NULL ) alloc_error(const_cast<char*>("interest_var"));
n_interest_var = 0;
for ( j = 0; j < p_ilp->mc; j++ ) {
if ( interesting_var(j) ) {
interest_var[n_interest_var] = j;
n_interest_var++;
}
}
best_var_cut = (select_cut **) calloc(n_interest_var,sizeof(select_cut *));
if ( best_var_cut == NULL ) alloc_error(const_cast<char*>("best_var_cut"));
for ( k = 0; k < n_interest_var; k++ ) {
best_var_cut[k] = (select_cut *) calloc(1,sizeof(select_cut));
if ( best_var_cut[k] == NULL ) alloc_error(const_cast<char*>("best_var_cut[k]"));
best_var_cut[k]->ccoef = (int *) calloc(p_ilp->mc,sizeof(int));
if ( best_var_cut[k]->ccoef == NULL )
alloc_error(const_cast<char*>("best_var_cut[k]->ccoef"));
best_var_cut[k]->score = -INF;
}
/* find the cuts with the best scores and the list of the cuts
violated by the current fractional point */
maxscore = -INF;
ccoef = (int *) calloc(p_ilp->mc,sizeof(int));
if ( ccoef == NULL ) alloc_error(const_cast<char*>("ccoef"));
/* second_(&tf); */
/* talloc += tf - tti; */
min_cut_score = MIN_CUT_SCORE;
max_cut_score = MAX_CUT_SCORE;
for ( c = 0; c < pool->cnum; c++ ) {
p_cut = pool->list[c];
/* if a new separation was made before the call check only new cuts */
if ( after_sep && p_cut->it_found != sep_iter ) continue;
/* determine the actual coefficients of the cut */
/* second_(&ti); */
ok = get_ori_cut_coef(p_cut->n_of_constr,p_cut->constr_list,
ccoef,&crhs,TRUE);
/* second_(&tf); */
/* tgetcoef += tf - ti; */
/* second_(&ti); */
ok = ok && best_cut(ccoef,&crhs,&viol,TRUE,TRUE);
/* second_(&tf); */
/* tbestcut += tf - ti; */
if ( ok && viol > MIN_VIOLATION ) {
(p_cut->n_it_violated)++;
/* second_(&ti); */
score = cut_score(ccoef,crhs,viol,TRUE);
if ( score > maxscore ) {
maxscore = score;
maxc = c;
}
/* second_(&tf); */
/* tscore += tf - ti; */
if ( ! aggr ) {
if ( score < min_cut_score ) continue;
if ( score < max_cut_score && p_cut->n_it_violated < MIN_IT_VIOL )
continue;
}
/* second_(&ti); */
for ( k = 0; k < n_interest_var; k++ ) {
j = interest_var[k];
if ( ccoef[j] != 0 && best_var_cut[k]->score < score ) {
best_var_cut[k]->score = score;
for ( l = 0; l < p_ilp->mc; l++ )
best_var_cut[k]->ccoef[l] = ccoef[l];
best_var_cut[k]->crhs = crhs;
best_var_cut[k]->pool_index = c;
}
}
/* second_(&tf); */
/* tupdbest += tf - ti; */
}
else p_cut->n_it_violated = 0;
}
free(ccoef);
/* printf("maxscore of a cut : %f ever: %f\n",maxscore,max_score_ever); */
if ( maxscore > max_score_ever ) max_score_ever = maxscore;
/* second_(&ti); */
add_cuts = initialize_cut_list(n_interest_var);
added = (short int *) calloc(pool->cnum,sizeof(short int));
if ( added == NULL ) alloc_error(const_cast<char*>("added"));
for ( c = 0; c < pool->cnum; c++ ) added[c] = FALSE;
for ( k = 0; k < n_interest_var; k++ ) {
j = interest_var[k];
if ( ! added[best_var_cut[k]->pool_index] &&
best_var_cut[k]->score >= MIN_CUT_SCORE ) {
/* add the cut to the cut list on output */
a_cut = define_cut(best_var_cut[k]->ccoef,best_var_cut[k]->crhs);
add_cuts = add_cut_to_list(a_cut,add_cuts);
added[best_var_cut[k]->pool_index] = TRUE;
}
free(best_var_cut[k]->ccoef);
free(best_var_cut[k]);
}
free(best_var_cut);
free(interest_var);
free(added);
/* second_(&tf); */
/* tadd += tf - ti; */
/* second_(&ttf); */
/* tgetcut += ttf - tti; */
/* printf("get_cuts_from_pool: talloc %f tgetcoef %f tbestcut %f tscore %f tupdbest %f tadd %f tgetcut %f\n",talloc,tgetcoef,tbestcut,tscore,tupdbest,tadd,tgetcut); */
return(add_cuts);
}
#endif
/* initialize_log_var: initialize the log information for the problem variables */
void Cgl012Cut::initialize_log_var()
{
int j;
if (!vlog) {
if (p_ilp->mc) {
vlog = reinterpret_cast<log_var **> (calloc(p_ilp->mc,sizeof(log_var *)));
if ( vlog == NULL ) alloc_error(const_cast<char*>("vlog"));
for ( j = 0; j < p_ilp->mc; j++ ) {
vlog[j] = reinterpret_cast<log_var *> (calloc(1,sizeof(log_var)));
if ( vlog[j] == NULL ) alloc_error(const_cast<char*>("vlog[j]"));
vlog[j]->n_it_zero = 0;
}
}
} else {
// just initialize counts
for ( j = 0; j < p_ilp->mc; j++ ) {
vlog[j]->n_it_zero = 0;
}
}
}
/* free_log_var */
void Cgl012Cut::free_log_var()
{
if (vlog) {
int j;
for ( j = 0; j < p_ilp->mc; j++ ) free(vlog[j]);
free(vlog);
vlog=NULL;
}
}
/* update_log_var: update the log information for the problem variables */
void Cgl012Cut::update_log_var()
{
int j;
/* so far one counts only the number of consecutive iterations with
0 value for each variable */
if ( vlog == NULL ) initialize_log_var();
for ( j = 0; j < p_ilp->mc; j++ ) {
if ( p_ilp->xstar[j] < ZERO && p_ilp->xstar[j] > - ZERO )
vlog[j]->n_it_zero++;
else vlog[j]->n_it_zero = 0;
}
}
/* the final implementation should use the following additional functions:
init_sep_012_cut: defines the inp_ilp and p_ilp data structures and
initializes pool and vlog
sep_012_cut: updates inp_ilp and p_ilp, performs separation and pool
management
kill_sep_012_cut: frees all the permanent memory (inp_ilp, p_ilp,
pool, vlog, etc.)
*/
int Cgl012Cut::sep_012_cut(
/*
INPUT parameters:
*/
int mr, /* number of rows in the ILP matrix */
int mc, /* number of columns in the ILP matrix */
int mnz, /* number of nonzero's in the ILP matrix */
int *mtbeg, /* starting position of each row in arrays mtind and mtval */
int *mtcnt, /* number of entries of each row in arrays mtind and mtval */
int *mtind, /* column indices of the nonzero entries of the ILP matrix */
int *mtval, /* values of the nonzero entries of the ILP matrix */
int *vlb, /* lower bounds on the variables */
int *vub, /* upper bounds on the variables */
int *mrhs, /* right hand sides of the constraints */
char *msense, /* senses of the constraints: 'L', 'G' or 'E' */
const double *xstar, /* current optimal solution of the LP relaxation */
bool aggressive, /* flag asking whether as many cuts as possible are
required on output (TRUE) or not (FALSE) */
/*
OUTPUT parameters (the memory for the vectors is allocated INTERNALLY
by the procedure: if some memory is already allocated, it is FREED):
*/
int *cnum, /* number of violated 0-1/2 cuts identified by the procedure */
int *cnzcnt, /* overall number of nonzero's in the cuts */
int **cbeg, /* starting position of each cut in arrays cind and cval */
int **ccnt, /* number of entries of each cut in arrays cind and cval */
int **cind, /* column indices of the nonzero entries of the cuts */
int **cval, /* values of the nonzero entries of the cuts */
int **crhs, /* right hand sides of the cuts */
char **csense /* senses of the cuts: 'L', 'G' or 'E' */
/*
NOTE that all the numerical input/output vectors are INTEGER (with
the exception of xstar), since the procedure is intended to work
with pure ILP's, and that the ILP matrix has to be given on input
in ROW format.
*/
)
{
#ifdef TIME
float tbasi, tbasf;
#endif
errorNo=0;
cut_list *out_cuts, *add_cuts;
#ifdef TIME
second_(&tti);
#endif
aggr = aggressive;
/* load the input ILP into an internal data structure */
//ilp_load(mr,mc,mnz,mtbeg,mtcnt,mtind,mtval,
// vlb,vub,mrhs,msense);
inp_ilp->xstar = xstar;
/* construct an internal data structure containing all the information
which can be useful for 0-1/2 cut separation - this may in fact be
done only at the first call of the separation procedure */
#ifdef TIME
second_(&ti);
#endif
get_parity_ilp();
/*
print_double_vect("xstar",p_ilp->xstar,p_ilp->mc);
print_parity_ilp();
*/
#ifdef TIME
second_(&tf);
prep_time += tf - ti;
#endif
if ( p_ilp->mnz == 0 ) {
#ifdef COIN_DEVELOP
printf("Warning: no significant constraint for 0-1/2 cut separation\n");
printf("... end separation\n");
#endif
//free_ilp();
//free_parity_ilp();
#ifdef TIME
second_(&ttf);
total_time += ttf - tti;
#endif
return(FALSE);
}
/* print_double_vect("xstar",p_ilp->xstar,p_ilp->mc); */
sep_iter++;
update_log_var();
#ifdef POOL
/* search for possible violated cuts in the pool */
#ifdef TIME
second_(&tpi);
#endif
add_cuts = get_cuts_from_pool(FALSE);
#ifdef TIME
second_(&tpf);
pool_time += tpf - tpi;
#ifdef PRINT_TIME
printf("... time elapsed at the end of get_cuts_from_pool: %f\n",tpf - tti);
#endif
#endif
if ( add_cuts->cnum == 0 ) {
free_cut_list(add_cuts);
}
else {
/* printf("Violated cuts found in the pool - no separation procedure used\n"); */
goto free_memory;
}
#endif
/* try to identify violated 0-1/2 cuts by using the original procedure
described in Caprara and Fischetti's MP paper */
#ifdef TIME
second_(&tbasi);
#endif
out_cuts = basic_separation();
#ifdef TIME
second_(&tbasf);
tot_basic_sep_time += tbasf - tbasi;
avg_basic_sep_time = tot_basic_sep_time / (float) sep_iter;
#endif
//#define TABU_SEARCH
#ifdef TABU_SEARCH
/* try to identify violated cuts by tabu search if none was found */
if ( out_cuts->cnum == 0 ) {
free_cut_list(out_cuts);
#ifdef TIME
second_(&ttabi);
#endif
out_cuts = tabu_012();
#ifdef TIME
second_(&ttabf);
tabu_time += ttabf - ttabi;
#endif
}
#endif
#ifdef POOL
/* add the cuts separated to the pool */
#ifdef TIME
second_(&tpi);
#endif
add_cuts_to_pool(out_cuts);
free_cut_list(out_cuts);
/* select from the pool a convenient set of violated constraints
to be added to the current LP */
add_cuts = get_cuts_from_pool(TRUE);
#ifdef TIME
second_(&tpf);
pool_time += tpf - tpi;
#ifdef PRINT_TIME
printf("... time elapsed at the end of get_cuts_from_pool: %f\n",tpf - tti);
#endif
#endif
#else
/* give on output the cuts separated */
add_cuts = out_cuts;
#endif
//free_ilp();
//free_parity_ilp();
#ifdef TIME
second_(&ttf);
total_time += ttf - tti;
#endif
if ( add_cuts->cnum > 0 ) {
getcuts(add_cuts,cnum,cnzcnt,cbeg,ccnt,cind,cval,crhs,csense);
/* print_cut_list(add_cuts); */
free_cut_list(add_cuts);
return(TRUE);
}
else {
free_cut_list(add_cuts);
return(FALSE);
}
}
#ifdef TIME
/* print_times: print the timings of the separation procedure */
void print_times()
{
printf("... separation timings \n");
printf("times total: %f prep: %f weak: %f aux: %f path: %f cycle: %f cut: %f (%d calls) bw: %f coef: %f pool: %f tabu: %f\n",
total_time, prep_time, weak_time, aux_time, path_time, cycle_time, cut_time, cut_ncalls, bw_time, coef_time, pool_time, tabu_time);
}
#endif
//-------------------------------------------------------------------
// Default Constructor
//-------------------------------------------------------------------
Cgl012Cut::Cgl012Cut () :
inp_ilp(NULL),
p_ilp(NULL),
iter(0),
gap(0.0),
maxgap(0.0),
errorNo(0),
sep_iter(0),
vlog(NULL),
aggr(true)
{
// nothing to do here
}
//-------------------------------------------------------------------
// Copy constructor
//-------------------------------------------------------------------
Cgl012Cut::Cgl012Cut (const Cgl012Cut & rhs) :
inp_ilp(NULL),
p_ilp(NULL),
iter(rhs.iter),
gap(rhs.gap),
maxgap(rhs.maxgap),
errorNo(rhs.errorNo),
sep_iter(rhs.sep_iter),
vlog(NULL),
aggr(rhs.aggr)
{
if (rhs.p_ilp||rhs.vlog||inp_ilp)
abort();
}
//-------------------------------------------------------------------
// Destructor
//-------------------------------------------------------------------
Cgl012Cut::~Cgl012Cut ()
{
free_log_var();
free_parity_ilp();
free_ilp();
}
//----------------------------------------------------------------
// Assignment operator
//-------------------------------------------------------------------
Cgl012Cut &
Cgl012Cut::operator=(
const Cgl012Cut& rhs)
{
if (this != &rhs) {
if (rhs.p_ilp||rhs.vlog||inp_ilp)
abort();
free_log_var();
free_parity_ilp();
free_ilp();
#if 0
inp_ilp = reinterpret_cast<ilp *> (calloc(1,sizeof(ilp)));
if ( inp_ilp == NULL ) alloc_error(const_cast<char*>("inp_ilp"));
inp_ilp->mr = rhs.inp_ilp->mr; inp_ilp->mc = rhs.inp_ilp->mc;
inp_ilp->mnz = rhs.inp_ilp->mnz;
inp_ilp->mtbeg = rhs.inp_ilp->mtbeg; inp_ilp->mtcnt = rhs.inp_ilp->mtcnt;
inp_ilp->mtind = rhs.inp_ilp->mtind; inp_ilp->mtval = rhs.inp_ilp->mtval;
inp_ilp->vlb = rhs.inp_ilp->vlb; inp_ilp->vub = rhs.inp_ilp->vub;
inp_ilp->mrhs = rhs.inp_ilp->mrhs; inp_ilp->msense = rhs.inp_ilp->msense;
#endif
iter = rhs.iter;
gap = rhs.gap;
maxgap = rhs.maxgap;
errorNo = rhs.errorNo;
sep_iter = rhs.sep_iter;
aggr = rhs.aggr;
}
return *this;
}