1 // $Id$
2 // Copyright (C) 2010, International Business Machines
3 // Corporation and others. All Rights Reserved.
4 // This code is licensed under the terms of the Eclipse Public License (EPL).
5 /** @file 012cut.c Definition file for C coded 0-1/2 separator */
6 #include "CoinFinite.hpp"
7 #include "CoinTime.hpp"
8 #include "Cgl012cut.hpp"
9 #include "CglZeroHalf.hpp"
10 static int MAX_CUTS = 10000000;
11 //#define PRINT_TABU
12 //#define PRINT_CUTS
13 //#define PRINT_TIME
14 //#define TIME
15
16 /* #define TIME */
17 #undef TIME
18
19 #define TRUE 1
20 #define FALSE 0
21
22 #define ODD 1
23 #define EVEN 0
24
25 #define NONE -1
26 #define BOTH 2
27
28 #define IN 1
29 #define OUT 0
30 #define ADD 1
31 #define DEL 0
32
33 #define LOWER_BOUND 0
34 #define UPPER_BOUND 1
35
36 #define EPS 0.0001 /* small tolerance */
37 //#define EPS 0.000001 /* small tolerance */
38 #define ZERO 0.000001 /* estimated accuracy for doubles */
39 //#define ZERO 0.0001 /* estimated accuracy for doubles */
40 #define INF 1000000000.0
41 #define IINF 1000000000
42
43 #define MAX_SLACK 1.0
44 #define MAX_LOSS 1.0
45 #define MAX_CYCLE_WEIGHT 1.0
46 #define MIN_VIOLATION 0.001
47 #define MIN_SCORE_RANGE 10.0
48 #define MAX_SCORE_RANGE ZERO /* 1.0 */
49
50 #define ISCALE 10000
51
52 //#define MAX_CUTS 10
53
54 #define MAX_CUT_POOL 10000
55 #define MAX_CUT_COD 10000
56 #define MAX_ITER_POOL 100
57 #define CLEAN_THRESH 0.9
58 #define MANY_IT_ZERO 10
59
60 #define mod2(I) ( I % 2 == 0 ? 0 : 1 )
61
62
63 #ifdef TIME
64
65 static float tot_basic_sep_time = 0.0; /* total time spent for basic
66 separation */
67 static float avg_basic_sep_time; /* average time per iteration spent for basic
68 separation */
69 static float total_time = 0.0; /* total time spent in the separation */
70 static float prep_time = 0.0; /* time spent for the definition of the
71 parity ILP data structure */
72 static float weak_time = 0.0; /* time spent for the construction of the
73 separation graph by weakening */
74 static float aux_time = 0.0; /* time spent for the definition of the
75 auxiliary graph */
76 static float path_time = 0.0; /* time spent in the computation of the
77 shortest paths */
78 static float cycle_time = 0.0; /* time spent in the determination of the
79 shortest cycles */
80 static float cut_time = 0.0; /* time spent in the determination of the
81 violated cuts */
82 static float bw_time = 0.0; /* time spent in best_weakening */
83 static float coef_time = 0.0; /* time spent in the initial computation
84 of coef in get_cut */
85 static float pool_time = 0.0; /* time spent for the addition and
86 extraction of cuts from the pool */
87 static int cut_ncalls = 0; /* number of calls to get_cut */
88 static float tabu_time = 0.0; /* time spent within tabu search */
89 float ti, tf, td, tti, ttf, tsi, tsf, tii, tff, tpi, tpf, ttabi, ttabf;
90 void
second_(float * t)91 second_(float *t) {*t=CoinCpuTime();}
92 #endif
93
94 /* #endif */
95
96 /* global data structures */
97
98 #define CGGGGG
99 #ifndef CGGGGG
100 static ilp *inp_ilp; /* input ILP data structure */
101 static parity_ilp *p_ilp; /* parity ILP data structure */
102 #endif
103
104 #ifdef PRINT_CUTS
105 /* utility subroutines */
106
print_int_vect(char * s,int * v,int n)107 void print_int_vect(char *s,int *v,int n)
108 {
109 int i;
110
111 printf("integer vector %s:",s);
112 for ( i = 0; i < n; i++ ) printf(" %d",v[i]);
113 printf("\n");
114 }
115
print_short_int_vect(char * s,short int * v,int n)116 void print_short_int_vect(char *s,short int *v,int n)
117 {
118 int i;
119
120 printf("short integer vector %s:",s);
121 for ( i = 0; i < n; i++ ) printf(" %d",v[i]);
122 printf("\n");
123 }
124
print_double_vect(char * s,double * v,int n)125 void print_double_vect(char *s,double *v,int n)
126 {
127 int i;
128
129 printf("double vector %s:",s);
130 for ( i = 0; i < n; i++ ) printf(" %f",v[i]);
131 printf("\n");
132 }
133 #endif
alloc_error(char * s)134 void alloc_error(char *s)
135 {
136 printf("\n Warning: Not enough memory to allocate %s\n",s);
137 printf("\n Cannot proceed with 0-1/2 cut separation\n");
138 exit(FALSE);
139 }
140
141 /* double2int: compute the integer equivalent of a double */
142
double2int(double x)143 int double2int(double x)
144 {
145 if ( x > IINF ) return (IINF);
146 if ( x < - IINF ) return (- IINF);
147 if ( x < ZERO && x > - ZERO ) return(0);
148 if ( x > 0.0 ) return(static_cast<int> (x + ZERO));
149 return(static_cast<int> (x - ZERO));
150 }
151
152 /* gcd: compute the greatest common divisor of two integers */
153
gcd(int a,int b)154 int gcd(int a,int b)
155 {
156 int c;
157
158 if ( a < 0 ) a = - a;
159 if ( b < 0 ) b = - b;
160 if ( a < b ) { c = a; a = b; b = c; }
161 while ( b != 0 ) {
162 c = a % b; a = b; b = c;
163 }
164 return(a);
165 }
166
167 /* ILP data structures subroutines */
168
169 /* ilp_load: load the input ILP into an internal data structure */
170
ilp_load(int mr,int mc,int mnz,int * mtbeg,int * mtcnt,int * mtind,int * mtval,int * vlb,int * vub,int * mrhs,char * msense)171 void Cgl012Cut::ilp_load(
172 int mr, /* number of rows in the ILP matrix */
173 int mc, /* number of columns in the ILP matrix */
174 int mnz, /* number of nonzero's in the ILP matrix */
175 int *mtbeg, /* starting position of each row in arrays mtind and mtval */
176 int *mtcnt, /* number of entries of each row in arrays mtind and mtval */
177 int *mtind, /* column indices of the nonzero entries of the ILP matrix */
178 int *mtval, /* values of the nonzero entries of the ILP matrix */
179 int *vlb, /* lower bounds on the variables */
180 int *vub, /* upper bounds on the variables */
181 int *mrhs, /* right hand sides of the constraints */
182 char *msense /* senses of the constraints: 'L', 'G' or 'E' */
183 )
184 {
185 inp_ilp = reinterpret_cast<ilp *> (calloc(1,sizeof(ilp)));
186 if ( inp_ilp == NULL ) alloc_error(const_cast<char*>("inp_ilp"));
187
188 inp_ilp->mr = mr; inp_ilp->mc = mc; inp_ilp->mnz = mnz;
189 inp_ilp->mtbeg = mtbeg; inp_ilp->mtcnt = mtcnt;
190 inp_ilp->mtind = mtind; inp_ilp->mtval = mtval;
191 inp_ilp->vlb = vlb; inp_ilp->vub = vub;
192 inp_ilp->mrhs = mrhs; inp_ilp->msense = msense;
193 }
free_ilp()194 void Cgl012Cut::free_ilp()
195 {
196 free(inp_ilp);
197 inp_ilp=NULL;
198 }
199
200 #ifdef PRINT_CUTS
print_constr(int i)201 void Cgl012Cut::print_constr(int i /* constraint to be printed */)
202 {
203
204 printf("\n content of constraint %d: nzcnt = %d, rhs = %d, sense = %c, slack = %f\n",
205 i, inp_ilp->mtcnt[i], inp_ilp->mrhs[i], inp_ilp->msense[i], p_ilp->slack[i]);
206 print_int_vect(const_cast<char*>("ind"),inp_ilp->mtind + inp_ilp->mtbeg[i],inp_ilp->mtcnt[i]);
207 print_int_vect(const_cast<char*>("val"),inp_ilp->mtval + inp_ilp->mtbeg[i],inp_ilp->mtcnt[i]);
208 }
209 #endif
210
211 /* alloc_parity_ilp: allocate the memory for the parity ILP data structure */
212
alloc_parity_ilp(int mr,int mc,int mnz)213 void Cgl012Cut::alloc_parity_ilp(
214 int mr, /* number of rows in the ILP matrix */
215 int mc, /* number of columns in the ILP matrix */
216 int mnz /* number of nonzero's in the ILP matrix */
217 )
218 {
219 p_ilp = reinterpret_cast<parity_ilp *> (calloc(1,sizeof(parity_ilp)));
220 if ( p_ilp == NULL ) alloc_error(const_cast<char*>("p_ilp"));
221
222 p_ilp->mtbeg = reinterpret_cast<int *> (calloc(mr,sizeof(int)));
223 if ( p_ilp->mtbeg == NULL ) alloc_error(const_cast<char*>("p_ilp->mtbeg"));
224 p_ilp->mtcnt = reinterpret_cast<int *> (calloc(mr,sizeof(int)));
225 if ( p_ilp->mtcnt == NULL ) alloc_error(const_cast<char*>("p_ilp->mtcnt"));
226 p_ilp->mtind = reinterpret_cast<int *> (calloc(mnz,sizeof(int)));
227 if ( p_ilp->mtind == NULL ) alloc_error(const_cast<char*>("p_ilp->mtind"));
228 p_ilp->mrhs = reinterpret_cast<short int *> (calloc(mr,sizeof(short int)));
229 if ( p_ilp->mrhs== NULL ) alloc_error(const_cast<char*>("p_ilp->mrhs"));
230 p_ilp->xstar = reinterpret_cast<double *> (calloc(mc,sizeof(double)));
231 if ( p_ilp->xstar== NULL ) alloc_error(const_cast<char*>("p_ilp->xstar"));
232 p_ilp->slack = reinterpret_cast<double *> (calloc(mr,sizeof(double)));
233 if ( p_ilp->slack == NULL ) alloc_error(const_cast<char*>("p_ilp->slack"));
234 p_ilp->row_to_delete = reinterpret_cast<short int *> (calloc(mr,sizeof(short int)));
235 if ( p_ilp->row_to_delete == NULL ) alloc_error(const_cast<char*>("p_ilp->row_to_delete"));
236 p_ilp->col_to_delete = reinterpret_cast<short int *> (calloc(mc,sizeof(short int)));
237 if ( p_ilp->col_to_delete == NULL ) alloc_error(const_cast<char*>("p_ilp->col_to_delete"));
238 p_ilp->gcd = reinterpret_cast<int *> (calloc(mr,sizeof(int)));
239 if ( p_ilp->gcd == NULL ) alloc_error(const_cast<char*>("p_ilp->gcd"));
240 p_ilp->possible_weak = reinterpret_cast<short int *> (calloc(mc,sizeof(short int)));
241 if ( p_ilp->possible_weak == NULL ) alloc_error(const_cast<char*>("p_ilp->possible_weak"));
242 p_ilp->type_even_weak = reinterpret_cast<short int *> (calloc(mc,sizeof(short int)));
243 if ( p_ilp->type_even_weak == NULL ) alloc_error(const_cast<char*>("p_ilp->type_even_weak"));
244 p_ilp->type_odd_weak = reinterpret_cast<short int *> (calloc(mc,sizeof(short int)));
245 if ( p_ilp->type_odd_weak == NULL ) alloc_error(const_cast<char*>("p_ilp->type_odd_weak"));
246 p_ilp->loss_even_weak = reinterpret_cast<double *> (calloc(mc,sizeof(double)));
247 if ( p_ilp->loss_even_weak == NULL ) alloc_error(const_cast<char*>("p_ilp->loss_even_weak"));
248 p_ilp->loss_odd_weak = reinterpret_cast<double *> (calloc(mc,sizeof(double)));
249 if ( p_ilp->loss_odd_weak == NULL ) alloc_error(const_cast<char*>("p_ilp->loss_odd_weak"));
250 p_ilp->min_loss_by_weak = reinterpret_cast<double *> (calloc(mc,sizeof(double)));
251 if ( p_ilp->min_loss_by_weak == NULL ) alloc_error(const_cast<char*>("p_ilp->min_loss_by_weak"));
252 p_ilp->mr=mr;
253 p_ilp->mc=mc;
254 p_ilp->mnz=mnz;
255 }
256
257 #ifdef PRINT_CUTS
print_parity_ilp()258 void Cgl012Cut::print_parity_ilp()
259 {
260 printf("\n content of parity_ilp data structure: mc = %d, mr = %d, mnz = %d\n",
261 p_ilp->mc,p_ilp->mr,p_ilp->mnz);
262 print_int_vect(const_cast<char*>("mtbeg"),p_ilp->mtbeg,p_ilp->mr);
263 print_int_vect(const_cast<char*>("mtcnt"),p_ilp->mtcnt,p_ilp->mr);
264 print_int_vect(const_cast<char*>("mtind"),p_ilp->mtind,p_ilp->mnz);
265 print_short_int_vect(const_cast<char*>("mrhs"),p_ilp->mrhs,p_ilp->mr);
266 print_double_vect(const_cast<char*>("xstar"),p_ilp->xstar,p_ilp->mc);
267 print_double_vect(const_cast<char*>("slack"),p_ilp->slack,p_ilp->mr);
268 print_short_int_vect(const_cast<char*>("row_to_delete"),p_ilp->row_to_delete,p_ilp->mr);
269 print_short_int_vect(const_cast<char*>("col_to_delete"),p_ilp->col_to_delete,p_ilp->mc);
270 print_int_vect(const_cast<char*>("gcd"),p_ilp->gcd,p_ilp->mr);
271 print_short_int_vect(const_cast<char*>("possible_weak"),p_ilp->possible_weak,p_ilp->mc);
272 print_short_int_vect(const_cast<char*>("type_even_weak"),p_ilp->type_even_weak,p_ilp->mc);
273 print_short_int_vect(const_cast<char*>("type_odd_weak"),p_ilp->type_odd_weak,p_ilp->mc);
274 print_double_vect(const_cast<char*>("loss_even_weak"),p_ilp->loss_even_weak,p_ilp->mc);
275 print_double_vect(const_cast<char*>("loss_odd_weak"),p_ilp->loss_odd_weak,p_ilp->mc);
276 print_double_vect(const_cast<char*>("min_loss_by_weak"),p_ilp->min_loss_by_weak,p_ilp->mc);
277 }
278 #endif
279
free_parity_ilp()280 void Cgl012Cut::free_parity_ilp()
281 {
282 if (p_ilp) {
283 free(p_ilp->mtbeg);
284 free(p_ilp->mtcnt);
285 free(p_ilp->mtind);
286 free(p_ilp->mrhs);
287 free(p_ilp->xstar);
288 free(p_ilp->slack);
289 free(p_ilp->row_to_delete);
290 free(p_ilp->col_to_delete);
291 free(p_ilp->gcd);
292 free(p_ilp->possible_weak);
293 free(p_ilp->type_even_weak);
294 free(p_ilp->type_odd_weak);
295 free(p_ilp->loss_even_weak);
296 free(p_ilp->loss_odd_weak);
297 free(p_ilp->min_loss_by_weak);
298 free(p_ilp);
299 p_ilp=NULL;
300 }
301 }
302
303 /* alloc_info_weak: allocate memory for the weakening info data structure */
304
alloc_info_weak(int nweak)305 info_weak *alloc_info_weak(int nweak /* number of variables to be weakened */)
306 {
307 info_weak *i_weak;
308
309 i_weak = reinterpret_cast<info_weak *> (calloc(1,sizeof(info_weak)));
310 if ( i_weak == NULL ) alloc_error(const_cast<char*>("i_weak"));
311 if ( nweak > 0 ) {
312 i_weak->var = reinterpret_cast<int *> (calloc(nweak,sizeof(int)));
313 if ( i_weak->var == NULL ) alloc_error(const_cast<char*>("i_weak->var"));
314 i_weak->type = reinterpret_cast<short int *> (calloc(nweak,sizeof(short int)));
315 if ( i_weak->type == NULL ) alloc_error(const_cast<char*>("i_weak->type"));
316 }
317
318 return(i_weak);
319 }
320
321 #ifdef PRINT_CUTS
print_info_weak(info_weak * i_weak)322 void print_info_weak(info_weak *i_weak)
323 {
324 printf("\n content of info_weak: nweak = %d\n",i_weak->nweak);
325 if ( i_weak->nweak > 0 ) {
326 print_int_vect(const_cast<char*>("var"),i_weak->var,i_weak->nweak);
327 print_short_int_vect(const_cast<char*>("type"),i_weak->type,i_weak->nweak);
328 }
329 }
330 #endif
331
free_info_weak(info_weak * i_weak)332 void free_info_weak(info_weak *i_weak)
333 {
334 if ( i_weak->nweak > 0 ) {
335 free(i_weak->var);
336 free(i_weak->type);
337 }
338 free(i_weak);
339 }
340
341 /* get_parity_ilp: construct an internal data structure containing all the
342 information which can be useful for 0-1/2 cut separation */
343
get_parity_ilp()344 void Cgl012Cut::get_parity_ilp()
345 {
346 int i, j, h, ij, aij, cnti, cnttot, begi, begh, ofsj, gcdi, ubj, lbj;
347 double slacki, xstarj, loss_upper, loss_lower;
348 short int parity_col_removed, equalih;
349
350 /* allocate the memory for the parity ILP data structure */
351
352 //alloc_parity_ilp(inp_ilp->mr,inp_ilp->mc,inp_ilp->mnz);
353
354 p_ilp->mr = inp_ilp->mr;
355 p_ilp->mc = inp_ilp->mc;
356
357 /* mark the variables equal to their lower/upper bound */
358
359 parity_col_removed = 0;
360 for ( j = 0; j < inp_ilp->mc; j++ ) {
361 xstarj = p_ilp->xstar[j] = inp_ilp->xstar[j];
362 ubj = inp_ilp->vub[j]; lbj = inp_ilp->vlb[j];
363 if ( xstarj > static_cast<double> (ubj - ZERO) ) {
364 /* variable at its upper bound */
365 p_ilp->col_to_delete[j] = TRUE;
366 if ( mod2(ubj) == ODD ) {
367 p_ilp->possible_weak[j] = ODD;
368 p_ilp->type_odd_weak[j] = UPPER_BOUND;
369 p_ilp->loss_odd_weak[j] = 0.0;
370 if ( parity_col_removed == EVEN ) parity_col_removed = ODD;
371 else parity_col_removed = EVEN;
372 }
373 else {
374 p_ilp->possible_weak[j] = EVEN;
375 p_ilp->type_even_weak[j] = UPPER_BOUND;
376 p_ilp->loss_even_weak[j] = 0.0;
377 }
378 p_ilp->min_loss_by_weak[j] = 0.0;
379 }
380 else if ( xstarj < static_cast<double> (lbj) + ZERO ) {
381 /* variable at its lower bound */
382 p_ilp->col_to_delete[j] = TRUE;
383 if ( mod2(lbj) == ODD ) {
384 p_ilp->possible_weak[j] = ODD;
385 p_ilp->type_odd_weak[j] = LOWER_BOUND;
386 p_ilp->loss_odd_weak[j] = 0.0;
387 p_ilp->min_loss_by_weak[j] = 0.0;
388 if ( parity_col_removed == EVEN ) parity_col_removed = ODD;
389 else parity_col_removed = EVEN;
390 }
391 else {
392 p_ilp->possible_weak[j] = EVEN;
393 p_ilp->type_even_weak[j] = LOWER_BOUND;
394 p_ilp->loss_even_weak[j] = 0.0;
395 p_ilp->min_loss_by_weak[j] = 0.0;
396 }
397 p_ilp->min_loss_by_weak[j] = 0.0;
398 }
399 else {
400 /* variable neither at its lower nor at its upper bound */
401 p_ilp->col_to_delete[j] = FALSE;
402 loss_upper = static_cast<double> (ubj) - xstarj;
403 loss_lower = xstarj - static_cast<double> (lbj);
404 if ( ( loss_upper > MAX_LOSS ) && ( loss_lower > MAX_LOSS ) )
405 /* no weakening for the variable */
406 p_ilp->possible_weak[j] = NONE;
407 else if ( loss_upper > MAX_LOSS ) {
408 /* lower weakening only */
409 if ( mod2(lbj) == EVEN ) {
410 p_ilp->possible_weak[j] = EVEN;
411 p_ilp->type_even_weak[j] = LOWER_BOUND;
412 p_ilp->loss_even_weak[j] = loss_lower;
413 }
414 else {
415 p_ilp->possible_weak[j] = ODD;
416 p_ilp->type_odd_weak[j] = LOWER_BOUND;
417 p_ilp->loss_odd_weak[j] = loss_lower;
418 }
419 }
420 else if ( loss_lower > MAX_LOSS ) {
421 /* upper weakening only */
422 if ( mod2(ubj) == EVEN ) {
423 p_ilp->possible_weak[j] = EVEN;
424 p_ilp->type_even_weak[j] = UPPER_BOUND;
425 p_ilp->loss_even_weak[j] = loss_upper;
426 }
427 else {
428 p_ilp->possible_weak[j] = ODD;
429 p_ilp->type_odd_weak[j] = UPPER_BOUND;
430 p_ilp->loss_odd_weak[j] = loss_upper;
431 }
432 }
433 else if ( mod2(ubj) == mod2(lbj) ) {
434 /* lower and upper bound have the same parity:
435 choose the best weakening */
436 if ( mod2(ubj) == EVEN ) {
437 p_ilp->possible_weak[j] = EVEN;
438 if ( loss_lower <= loss_upper ) {
439 p_ilp->type_even_weak[j] = LOWER_BOUND;
440 p_ilp->loss_even_weak[j] = loss_lower;
441 }
442 else {
443 p_ilp->type_even_weak[j] = UPPER_BOUND;
444 p_ilp->loss_even_weak[j] = loss_upper;
445 }
446 }
447 else {
448 p_ilp->possible_weak[j] = ODD;
449 if ( loss_lower <= loss_upper ) {
450 p_ilp->type_odd_weak[j] = LOWER_BOUND;
451 p_ilp->loss_odd_weak[j] = loss_lower;
452 }
453 else {
454 p_ilp->type_odd_weak[j] = UPPER_BOUND;
455 p_ilp->loss_odd_weak[j] = loss_upper;
456 }
457 }
458 }
459 else {
460 /* lower and upper bound have different parities:
461 consider both weakenings */
462 p_ilp->possible_weak[j] = BOTH;
463 if ( mod2(ubj) == EVEN ) {
464 p_ilp->type_even_weak[j] = UPPER_BOUND;
465 p_ilp->loss_even_weak[j] = loss_upper;
466 p_ilp->type_odd_weak[j] = LOWER_BOUND;
467 p_ilp->loss_odd_weak[j] = loss_lower;
468 }
469 else {
470 p_ilp->type_even_weak[j] = LOWER_BOUND;
471 p_ilp->loss_even_weak[j] = loss_lower;
472 p_ilp->type_odd_weak[j] = UPPER_BOUND;
473 p_ilp->loss_odd_weak[j] = loss_upper;
474 }
475 }
476 if ( loss_upper > loss_lower ) p_ilp->min_loss_by_weak[j] = loss_lower;
477 else p_ilp->min_loss_by_weak[j] = loss_upper;
478 }
479 }
480
481 /* scan the constraints and delete those which are trivially useless
482 in the 0-1/2 cut separation */
483
484 cnttot = 0;
485 for ( i = 0; i < inp_ilp->mr; i++ ) {
486 begi = inp_ilp->mtbeg[i];
487
488 /* compute the row slack and the GCD of the entries of the row */
489
490 slacki = static_cast<double> (inp_ilp->mrhs[i]);
491 gcdi = inp_ilp->mrhs[i];
492 for ( ofsj = 0; ofsj < inp_ilp->mtcnt[i]; ofsj++ ) {
493 ij = begi + ofsj;
494 j = inp_ilp->mtind[ij];
495 aij = inp_ilp->mtval[ij];
496 slacki -= static_cast<double> (aij ) * ( inp_ilp->xstar[j] );
497 gcdi = gcd(gcdi,aij);
498 }
499 if ( inp_ilp->msense[i] == 'G' ) slacki = -slacki;
500 if ( slacki < -ZERO || ( inp_ilp->msense[i] == 'E' && slacki > ZERO ) ) {
501 #ifdef COIN_DEVELOP
502 printf("\n Warning: constraint %d in the model is violated:\n",i);
503 printf("\n 0-1/2 cut separation is not possible\n");
504 printf("\nnumber of nonzero's %d\n",inp_ilp->mtcnt[i]);
505 printf("nonzero's (col,coef,xstar) ");
506 for (ofsj=0;ofsj<inp_ilp->mtcnt[i];ofsj++) printf("(%d,%d,%f) ",
507 inp_ilp->mtind[begi+ofsj], inp_ilp->mtval[begi+ofsj],
508 inp_ilp->xstar[inp_ilp->mtind[begi+ofsj]]);
509 printf("\n");
510 printf("sense %c and rhs %d and slack %.5e\n",inp_ilp->msense[i],inp_ilp->mrhs[i], slacki);
511 #endif
512 //exit(0);
513 slacki = INF;
514 }
515 p_ilp->slack[i] = slacki;
516
517 /* mark the rows with slack greater than the maximum allowed */
518
519 if ( slacki > MAX_SLACK - EPS ) p_ilp->row_to_delete[i] = TRUE;
520 else p_ilp->row_to_delete[i] = FALSE;
521
522 /* store the odd entries in the (possibly scaled) row i */
523
524 //if ( gcdi != 1 )
525 //printf("Warning: constraint %d with nonprime coefficients\n",i);
526 p_ilp->gcd[i] = gcdi;
527 p_ilp->mrhs[i] = mod2(( inp_ilp->mrhs[i] / gcdi ));
528 p_ilp->mtbeg[i] = cnttot;
529 cnti = 0;
530 for ( ofsj = 0; ofsj < inp_ilp->mtcnt[i]; ofsj++ ) {
531 ij = begi + ofsj;
532 j = inp_ilp->mtind[ij];
533 aij = mod2(( inp_ilp->mtval[ij] / gcdi ));
534 if ( aij == ODD ) {
535 if ( ! p_ilp->col_to_delete[j] ) {
536 p_ilp->mtind[cnttot] = j;
537 cnti++; cnttot++;
538 }
539 else if ( p_ilp->possible_weak[j] == ODD ) {
540 if ( p_ilp->mrhs[i] == EVEN ) p_ilp->mrhs[i] = ODD;
541 else p_ilp->mrhs[i] = EVEN;
542 }
543 }
544 }
545 p_ilp->mtcnt[i] = cnti;
546 if ( cnti == 0 ) /* (scaled) row with even entries only */
547 p_ilp->row_to_delete[i] = TRUE;
548 else {
549 if ( cnti == 1 && slacki < EPS ) { /* the row could be deleted */
550 #ifdef PRINT
551 printf("get_parity_ilp: row %d could be deleted since it\n",i);
552 printf("has only one odd entry and is tight, but it is not ...\n");
553 #endif
554 }
555 }
556 }
557 p_ilp->mnz = cnttot;
558
559 #ifdef REDUCTION
560
561 /* remove identical rows in the parity matrix */
562 /* very trivial implementation */
563
564 for ( i = 0; i < p_ilp->mr; i++ )
565 for ( h = i+1; h < p_ilp->mr; h++ )
566 if ( ( p_ilp->mrhs[i] == p_ilp->mrhs[h] ) &&
567 ( p_ilp->mtcnt[i] == p_ilp->mtcnt[h] ) &&
568 ( ! p_ilp->row_to_delete[i] ) &&
569 ( ! p_ilp->row_to_delete[h] ) ) {
570 begi = p_ilp->mtbeg[i]; begh = p_ilp->mtbeg[h];
571 equalih = TRUE;
572 for ( ofsj = 0; ofsj < p_ilp->mtcnt[i]; ofsj++ )
573 /* the check assumes the indexes of the columns associated
574 with each row are ordered in p_ilp->mtind[] ... */
575 if ( p_ilp->mtind[begi+ofsj] != p_ilp->mtind[begh+ofsj] ) {
576 equalih = FALSE;
577 break;
578 }
579 if ( equalih ) {
580 if ( p_ilp->slack[h] > p_ilp->slack[i] )
581 p_ilp->row_to_delete[h] = TRUE;
582 else
583 p_ilp->row_to_delete[i] = TRUE;
584 }
585 }
586
587 /* check for the existence of separate connected components in the
588 parity matrix row intersection graph */
589 /* not implemented so far - if ever, the availability of the parity
590 matrix in column form also would be really convenient */
591
592 #endif
593
594 }
595
596 /* separation graph subroutines */
597
598 #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 ) )
599
600 /* initialize_sep_graph: allocate and initialize the data structure
601 to contain the information associated with a separation graph */
602
initialize_sep_graph()603 separation_graph *Cgl012Cut::initialize_sep_graph()
604 {
605 int maxnodes, maxedges, nnodes, j, jk;
606 int *nodes, *ind;
607 separation_graph *s_graph;
608
609 s_graph = reinterpret_cast<separation_graph *> (calloc(1,sizeof(separation_graph)));
610 if ( s_graph == NULL ) alloc_error(const_cast<char*>("s_graph"));
611
612 maxnodes = p_ilp->mc + 1;
613 nnodes = 0;
614 nodes = reinterpret_cast<int *> (calloc(maxnodes,sizeof(int)));
615 if ( nodes == NULL ) alloc_error(const_cast<char*>("nodes"));
616 ind = reinterpret_cast<int *> (calloc(maxnodes,sizeof(int)));
617 if ( ind == NULL ) alloc_error(const_cast<char*>("ind"));
618 for ( j = 0; j < p_ilp->mc; j++ )
619 if ( ! p_ilp->col_to_delete[j] ) {
620 /* variable not removed from the separation problem */
621 nodes[nnodes] = j;
622 ind[j] = nnodes;
623 nnodes++;
624 }
625
626 /* take into account the special node */
627 nodes[nnodes] = maxnodes - 1;
628 ind[maxnodes-1] = nnodes;
629 nnodes++;
630
631 s_graph->nnodes = nnodes;
632 s_graph->nedges = 0;
633 s_graph->nodes = reinterpret_cast<int *> (malloc(nnodes*sizeof(int)));
634 if ( s_graph->nodes == NULL ) alloc_error(const_cast<char*>("s_graph->nodes"));
635 for ( j = 0; j < nnodes; j++ ) s_graph->nodes[j] = nodes[j];
636 free(nodes);
637 s_graph->ind = reinterpret_cast<int *> (malloc(maxnodes*sizeof(int)));
638 if ( s_graph->ind == NULL ) alloc_error(const_cast<char*>("s_graph->ind"));
639 for ( j = 0; j < maxnodes; j++ ) s_graph->ind[j] = ind[j];
640 free(ind);
641 maxedges = (nnodes * (nnodes - 1)) / 2;
642 s_graph->even_adj_list = reinterpret_cast<edge **> (malloc(maxedges*sizeof(edge *)));
643 if ( s_graph->even_adj_list == NULL ) alloc_error(const_cast<char*>("s_graph->even_adj_list"));
644 s_graph->odd_adj_list = reinterpret_cast<edge **> (malloc(maxedges*sizeof(edge *)));
645 if ( s_graph->odd_adj_list == NULL ) alloc_error(const_cast<char*>("s_graph->odd_adj_list"));
646 for ( jk = 0; jk < maxedges; jk++ )
647 s_graph->even_adj_list[jk] = s_graph->odd_adj_list[jk] = NULL;
648
649 return(s_graph);
650 }
651
652 /* update_weight_sep_graph: consider a new edge obtained from the
653 (weakened) parity ILP and (possibly) add it to the separation graph */
654
update_weight_sep_graph(int j,int k,double weight,short int parity,int i,info_weak * i_weak,separation_graph * s_graph)655 separation_graph *update_weight_sep_graph(
656 int j, int k, /* endpoints of the new edge */
657 double weight, /* weight of the new edge */
658 short int parity, /* parity of the new edge */
659 int i, /* constraint associated with the new edge */
660 info_weak *i_weak, /* information associated with the weakening */
661 separation_graph *s_graph /* separation graph to be updated */
662 )
663 {
664 int indj, indk, indjk;
665 edge *old_edge, *new_edge;
666
667 indj = s_graph->ind[j]; indk = s_graph->ind[k];
668 indjk = SG_EDGE_INDEX(s_graph,indj,indk);
669 if ( parity == EVEN ) old_edge = s_graph->even_adj_list[indjk];
670 else old_edge = s_graph->odd_adj_list[indjk];
671 if ( old_edge == NULL ) {
672 /* edge is not in the graph */
673 new_edge = reinterpret_cast<edge *> (calloc(1,sizeof(edge)));
674 if ( new_edge == NULL ) alloc_error(const_cast<char*>("new_edge"));
675 new_edge->endpoint1 = indj; new_edge->endpoint2 = indk;
676 new_edge->weight = weight; new_edge->parity = parity;
677 new_edge->constr = i; new_edge->weak = i_weak;
678 (s_graph->nedges)++;
679 if ( parity == EVEN ) s_graph->even_adj_list[indjk] = new_edge;
680 else s_graph->odd_adj_list[indjk] = new_edge;
681 }
682 else {
683 /* edge is already in the graph */
684 if ( old_edge->weight > weight ) {
685 /* replace the old edge */
686 old_edge->weight = weight; old_edge->constr = i;
687 free_info_weak(old_edge->weak);
688 old_edge->weak = i_weak;
689 }
690 else {
691 /* keep the old edge */
692 free_info_weak(i_weak);
693 }
694 }
695
696 return(s_graph);
697 }
698
699 #ifdef PRINT_CUTS
print_edge(edge * e)700 void print_edge(edge *e)
701 {
702 printf("\n content of edge: endpoint1 = %d, endpoint2 = %d, weight = %f, parity = %d, constr = %d\n",
703 e->endpoint1,e->endpoint2,e->weight,e->parity,e->constr);
704 print_info_weak(e->weak);
705 }
706 #endif
707
free_edge(edge * e)708 void free_edge(edge *e)
709 {
710 if ( e->weak != NULL )
711 free_info_weak(e->weak);
712 free(e);
713 }
714
715 #ifdef PRINT_CUTS
print_sep_graph(separation_graph * s_graph)716 void print_sep_graph(separation_graph *s_graph)
717 {
718 int nnodes, maxedges, jk;
719
720 nnodes = s_graph->nnodes;
721 maxedges = (nnodes * (nnodes - 1)) / 2;
722 printf("\n content of separation_graph: nnodes = %d, nedges = %d\n",
723 nnodes, s_graph->nedges);
724 print_int_vect(const_cast<char*>("nodes"),s_graph->nodes,nnodes);
725 print_int_vect(const_cast<char*>("ind"),s_graph->ind,nnodes);
726 for ( jk = 0; jk < maxedges; jk++ ) {
727 if ( s_graph->even_adj_list[jk] != NULL )
728 print_edge(s_graph->even_adj_list[jk]);
729 if ( s_graph->odd_adj_list[jk] != NULL )
730 print_edge(s_graph->odd_adj_list[jk]);
731 }
732 }
733 #endif
734
free_sep_graph(separation_graph * s_graph)735 void free_sep_graph(separation_graph *s_graph)
736 {
737 int nnodes, maxedges, jk;
738
739 nnodes = s_graph->nnodes;
740 maxedges = (nnodes * (nnodes - 1)) / 2;
741 for ( jk = 0; jk < maxedges; jk++ ) {
742 if ( s_graph->even_adj_list[jk] != NULL )
743 free_edge(s_graph->even_adj_list[jk]);
744 if ( s_graph->odd_adj_list[jk] != NULL )
745 free_edge(s_graph->odd_adj_list[jk]);
746 }
747 free(s_graph->nodes);
748 free(s_graph->ind);
749 free(s_graph->even_adj_list);
750 free(s_graph->odd_adj_list);
751 free(s_graph);
752 }
753
754 /* auxiliary graph subroutines - depend on the shortest path code used */
755
756 #ifndef CGL_NEW_SHORT
757 // will error if we get here
758 #include "Cgldikbd.c"
759 #endif
760
761 #define AG_TWIN1(J) 2 * J
762 #define AG_TWIN2(J) 2 * J + 1
763 #define AG_MATE(J) 2 * static_cast<int> ( J / 2 ) + ( J % 2 == EVEN ? 1 : 0 )
764 #define AG_TYPE(J,K) ( (J % 2) == (K % 2) ? EVEN : ODD )
765 #define SG_ORIG(J) static_cast<int> (J / 2)
766
767 /* define_aux_graph: construct the auxiliary graph for the shortest
768 path computation - the data structure is based on that used by
769 Cherkassky, Goldberg and Radzik's shortest path codes */
770
define_aux_graph(separation_graph * s_graph)771 auxiliary_graph *define_aux_graph(separation_graph *s_graph /* input separation graph */)
772 {
773 int j, k, indjk, auxj1, auxj2, auxk1, auxk2, noutj, totoutj, narcs;
774 edge *s_edge;
775 auxiliary_graph *a_graph;
776
777 a_graph = reinterpret_cast<auxiliary_graph *> (calloc(1,sizeof(auxiliary_graph)));
778 if ( a_graph == NULL ) alloc_error(const_cast<char*>("a_graph"));
779
780 a_graph->nnodes = 2 * s_graph->nnodes;
781 a_graph->narcs = 4 * s_graph->nedges;
782
783 #ifndef CGL_NEW_SHORT
784 a_graph->nodes = reinterpret_cast<node *> (calloc((a_graph->nnodes + 1),sizeof(node)));
785 #else
786 a_graph->nodes = reinterpret_cast<cgl_node *> (calloc((a_graph->nnodes + 1),sizeof(cgl_node)));
787 #endif
788 if ( a_graph->nodes == NULL ) alloc_error(const_cast<char*>("a_graph->nodes"));
789 #ifndef CGL_NEW_SHORT
790 a_graph->arcs = reinterpret_cast<arc *> (calloc(((a_graph->narcs) + 1),sizeof(arc)));
791 #else
792 a_graph->arcs = reinterpret_cast<cgl_arc *> (calloc(((a_graph->narcs) + 1),sizeof(cgl_arc)));
793 #endif
794 if ( a_graph->arcs == NULL ) alloc_error(const_cast<char*>("a_graph->arcs"));
795
796 narcs = 0;
797 for ( j = 0; j < s_graph->nnodes; j++ ) {
798 /* count the number of edges incident with j in the separation graph */
799 totoutj = 0;
800 for ( k = 0; k < s_graph->nnodes; k++ )
801 if ( k != j ) {
802 indjk = SG_EDGE_INDEX(s_graph,j,k);
803 if ( s_graph->even_adj_list[indjk] != NULL ) totoutj++;
804 if ( s_graph->odd_adj_list[indjk] != NULL ) totoutj++;
805 }
806 auxj1 = AG_TWIN1(j); auxj2 = AG_TWIN2(j);
807 a_graph->nodes[auxj1].index = auxj1;
808 a_graph->nodes[auxj2].index = auxj2;
809 #ifndef CGL_NEW_SHORT
810 a_graph->nodes[auxj1].first = &(a_graph->arcs[narcs]);
811 a_graph->nodes[auxj2].first = &(a_graph->arcs[narcs+totoutj]);
812 #else
813 a_graph->nodes[auxj1].firstArc = &(a_graph->arcs[narcs]);
814 a_graph->nodes[auxj2].firstArc = &(a_graph->arcs[narcs+totoutj]);
815 #endif
816 /* add the edges as arcs outgoing from j to the auxiliary graph */
817 noutj = 0;
818 for ( k = 0; k < s_graph->nnodes; k++ ) {
819 if ( k != j ) {
820 auxk1 = AG_TWIN1(k); auxk2 = AG_TWIN2(k);
821 indjk = SG_EDGE_INDEX(s_graph,j,k);
822 s_edge = s_graph->even_adj_list[indjk];
823 if ( s_edge != NULL ) {
824 /* there is an even edge between j and k */
825 #ifndef CGL_NEW_SHORT
826 a_graph->arcs[narcs].len = a_graph->arcs[narcs+totoutj].len =
827 (int) (s_edge->weight * ISCALE);
828 a_graph->arcs[narcs].head = &(a_graph->nodes[auxk1]);
829 a_graph->arcs[narcs+totoutj].head = &(a_graph->nodes[auxk2]);
830 #else
831 a_graph->arcs[narcs].length = a_graph->arcs[narcs+totoutj].length =
832 static_cast<int> (s_edge->weight * ISCALE);
833 a_graph->arcs[narcs].to = auxk1;
834 a_graph->arcs[narcs+totoutj].to = auxk2;
835 #endif
836 narcs++; noutj++;
837 }
838 s_edge = s_graph->odd_adj_list[indjk];
839 if ( s_edge != NULL ) {
840 /* there is an odd edge between j and k */
841 #ifndef CGL_NEW_SHORT
842 a_graph->arcs[narcs].len = a_graph->arcs[narcs+totoutj].len =
843 (int) (s_edge->weight * ISCALE);
844 a_graph->arcs[narcs].head = &(a_graph->nodes[auxk2]);
845 a_graph->arcs[narcs+totoutj].head = &(a_graph->nodes[auxk1]);
846 #else
847 a_graph->arcs[narcs].length = a_graph->arcs[narcs+totoutj].length =
848 static_cast<int> (s_edge->weight * ISCALE);
849 a_graph->arcs[narcs].to = auxk2;
850 a_graph->arcs[narcs+totoutj].to = auxk1;
851 #endif
852 /* this looks really useless - to be removed ...
853 if ( noutj == 0 ) {
854 a_graph->nodes[auxj1].first = &(a_graph->arcs[narcs]);
855 a_graph->nodes[auxj2].first = &(a_graph->arcs[narcs+totoutj]);
856 }
857 ... */
858 narcs++; noutj++;
859 }
860 }
861 }
862 narcs += totoutj;
863 }
864 #ifndef CGL_NEW_SHORT
865 a_graph->nodes[a_graph->nnodes].first = &(a_graph->arcs[narcs]);
866 #else
867 a_graph->nodes[a_graph->nnodes].firstArc = &(a_graph->arcs[narcs]);
868 #endif
869
870 return(a_graph);
871 }
872
873 /* cancel_node_aux_graph: remove the node j in the separation graph
874 from the auxiliary graph - all the outgoing arc lengths are set to
875 a large value */
876
cancel_node_aux_graph(int j,auxiliary_graph * a_graph)877 auxiliary_graph *cancel_node_aux_graph(
878 int j, /* index of the node in the separation graph */
879 auxiliary_graph *a_graph /* auxiliary graph to be updated */
880 )
881 {
882 int auxj1, auxj2;
883 #ifndef CGL_NEW_SHORT
884 arc *arc_ptr;
885 #else
886 cgl_arc *arc_ptr;
887 #endif
888
889 auxj1 = AG_TWIN1(j); auxj2 = AG_TWIN2(j);
890 #ifndef CGL_NEW_SHORT
891 for ( arc_ptr = a_graph->nodes[auxj1].first;
892 arc_ptr < a_graph->nodes[auxj1+1].first;
893 arc_ptr++ )
894 (*arc_ptr).len = ISCALE;
895 for ( arc_ptr = a_graph->nodes[auxj2].first;
896 arc_ptr < a_graph->nodes[auxj2+1].first;
897 arc_ptr++ )
898 (*arc_ptr).len = ISCALE;
899 #else
900 for ( arc_ptr = a_graph->nodes[auxj1].firstArc;
901 arc_ptr < a_graph->nodes[auxj1+1].firstArc;
902 arc_ptr++ )
903 (*arc_ptr).length = ISCALE;
904 for ( arc_ptr = a_graph->nodes[auxj2].firstArc;
905 arc_ptr < a_graph->nodes[auxj2+1].firstArc;
906 arc_ptr++ )
907 (*arc_ptr).length = ISCALE;
908 #endif
909
910 return(a_graph);
911 }
912
913 #ifdef PRINT_CUTS
914 #ifndef CGL_NEW_SHORT
print_node(node * n)915 void print_node(node *n)
916 {
917 printf("\n content of node (addr = %d): first = %d, dist = %d, parent = %d, next = %d, prev = %d, status = %d\n",
918 (int)n,(int)(*n).first,(*n).dist,(int)(*n).parent,(int)(*n).next,(int)(*n).prev,
919 (*n).status);
920 }
921
print_arc(arc * a)922 void print_arc(arc *a)
923 {
924 printf("\n content of arc (addr = %d): len = %d, head = %d\n",
925 (int)a,(*a).len,(int)(*a).head);
926 }
927
print_node_vect(char * s,node * v,int n)928 void print_node_vect(char *s,node *v,int n)
929 {
930 int i;
931
932 printf("node vector %s:",s);
933 for ( i = 0; i < n; i++ ) print_node(&v[i]);
934 printf("\n");
935 }
936
print_arc_vect(char * s,arc * v,int n)937 void print_arc_vect(char *s,arc *v,int n)
938 {
939 int i;
940
941 printf("arc vector %s:",s);
942 for ( i = 0; i < n; i++ ) print_arc(&v[i]);
943 printf("\n");
944 }
945 #else
print_node(cgl_node * n)946 void print_node(cgl_node *n)
947 {
948 printf("\n content of node (addr = %p): first = %p, dist = %d, parent = %p\n",
949 reinterpret_cast<void *>(n),
950 reinterpret_cast<void *>((*n).firstArc),
951 (*n).distanceBack,
952 static_cast<int>((*n).parentNode));
953 }
954
print_arc(cgl_arc * a)955 void print_arc(cgl_arc *a)
956 {
957 printf("\n content of arc (addr = %p): len = %d, head = %d\n",
958 reinterpret_cast<void *>(a),(*a).length,static_cast<int>((*a).to));
959 }
960
print_node_vect(char * s,cgl_node * v,int n)961 void print_node_vect(char *s,cgl_node *v,int n)
962 {
963 int i;
964
965 printf("node vector %s:",s);
966 for ( i = 0; i < n; i++ ) print_node(&v[i]);
967 printf("\n");
968 }
969
print_arc_vect(char * s,cgl_arc * v,int n)970 void print_arc_vect(char *s,cgl_arc *v,int n)
971 {
972 int i;
973
974 printf("arc vector %s:",s);
975 for ( i = 0; i < n; i++ ) print_arc(&v[i]);
976 printf("\n");
977 }
978 #endif
979
print_aux_graph(auxiliary_graph * a_graph)980 void print_aux_graph(auxiliary_graph *a_graph)
981 {
982 printf("\n content of auxiliary graph: nnodes = %d, narcs = %d\n",
983 a_graph->nnodes,a_graph->narcs);
984 print_node_vect(const_cast<char*>("nodes"),a_graph->nodes,a_graph->nnodes);
985 print_arc_vect(const_cast<char*>("nodes"),a_graph->arcs,a_graph->narcs);
986 }
987 #endif
988
free_aux_graph(auxiliary_graph * a_graph)989 void free_aux_graph(auxiliary_graph *a_graph)
990 {
991 free(a_graph->nodes);
992 free(a_graph->arcs);
993 free(a_graph);
994 }
995
996 /* odd cycles management subroutines */
997
998 /* simple_cycle: check whether a given cycle is simple
999 (and therefore may correspond to a non-dominated ineq.) */
1000
simple_cycle(cycle * s_cyc)1001 short int simple_cycle(cycle *s_cyc /* cycle to be checked */)
1002 {
1003 int i, e, maxnodes;
1004 int *cnt;
1005
1006 maxnodes = 0;
1007 for ( e = 0; e < s_cyc->length; e++ ) {
1008 if (!s_cyc->edge_list[e]) {
1009 // bad
1010 maxnodes=-1;
1011 abort();//break;
1012 }
1013 i = s_cyc->edge_list[e]->endpoint1;
1014 if ( i > maxnodes ) maxnodes = i;
1015 i = s_cyc->edge_list[e]->endpoint2;
1016 if ( i > maxnodes ) maxnodes = i;
1017 }
1018 if (maxnodes<0)
1019 return FALSE;
1020 cnt = reinterpret_cast<int *> (calloc(maxnodes+1,sizeof(int)));
1021 if ( cnt == NULL ) alloc_error(const_cast<char*>("cnt"));
1022 //for ( i = 0; i <= maxnodes; i++ ) cnt[i] = 0;
1023
1024 for ( e = 0; e < s_cyc->length; e++ ) {
1025 i = s_cyc->edge_list[e]->endpoint1;
1026 cnt[i]++;
1027 if ( cnt[i] > 2 ) {
1028 free(cnt);
1029 return(FALSE);
1030 }
1031 i = s_cyc->edge_list[e]->endpoint2;
1032 cnt[i]++;
1033 if ( cnt[i] > 2 ) {
1034 free(cnt);
1035 return(FALSE);
1036 }
1037 }
1038
1039 free(cnt);
1040 return(TRUE);
1041 }
1042
1043 /* same_cycle: check whether two cycles are identical
1044 (assumes the first nodes of the cycles coincide) */
1045
same_cycle(cycle * s_cyc1,cycle * s_cyc2)1046 short int same_cycle(cycle *s_cyc1, cycle *s_cyc2 /* cycles to be compared */)
1047 {
1048 int e, eb;
1049 short int same;
1050
1051 if ( s_cyc1->length != s_cyc2->length ) return(FALSE);
1052 /* check the cycles in the same direction ... */
1053 same = TRUE;
1054 for ( e = 0; e < s_cyc1->length; e++ ) {
1055 if ( s_cyc1->edge_list[e] != s_cyc2->edge_list[e] ) {
1056 same = FALSE;
1057 break;
1058 }
1059 }
1060 if ( same ) return(TRUE);
1061 /* ... and in reverse direction */
1062 same = TRUE;
1063 for ( e = 0, eb = s_cyc2->length - 1; e < s_cyc1->length; e++, eb-- ) {
1064 if ( s_cyc1->edge_list[e] != s_cyc2->edge_list[eb] ) {
1065 same = FALSE;
1066 break;
1067 }
1068 }
1069 if ( same ) return(TRUE);
1070 return(FALSE);
1071 }
1072
1073 #ifdef PRINT_CUTS
print_cycle(cycle * s_cycle)1074 void print_cycle(cycle *s_cycle)
1075 {
1076 int e;
1077
1078 printf("\n content of cycle: weight = %f, length = %d\n",
1079 s_cycle->weight,s_cycle->length);
1080 for ( e = 0; e < s_cycle->length; e++ )
1081 print_edge(s_cycle->edge_list[e]);
1082 }
1083 #endif
1084
free_cycle(cycle * s_cycle)1085 void free_cycle(cycle *s_cycle)
1086 {
1087 free(s_cycle->edge_list);
1088 free(s_cycle);
1089 }
1090
1091 /* initialize_cycle_list: allocate and initialize the cycle list data structure */
1092
initialize_cycle_list(int max_cyc)1093 cycle_list *initialize_cycle_list(int max_cyc /* maximum number of cycles in the list */)
1094 {
1095 cycle_list *s_cycle_list;
1096
1097 s_cycle_list = reinterpret_cast<cycle_list *> (calloc(1,sizeof(cycle_list)));
1098 if ( s_cycle_list == NULL ) alloc_error(const_cast<char*>("s_cycle_list"));
1099 s_cycle_list->cnum = 0;
1100 s_cycle_list->list = reinterpret_cast<cycle **> (calloc(max_cyc,sizeof(cycle *)));
1101 if ( s_cycle_list->list == NULL ) alloc_error(const_cast<char*>("s_cycle_list->list"));
1102 return(s_cycle_list);
1103 }
1104
1105 /* add_cycle_to_list: add a new cycle to the cycle list data structure
1106 (if not already in the list) */
1107
add_cycle_to_list(cycle * s_cycle,cycle_list * s_cycle_list)1108 cycle_list *add_cycle_to_list(
1109 cycle *s_cycle, /* pointer to the cycle to be added to the list */
1110 cycle_list *s_cycle_list /* input cycle list to be updated */
1111 )
1112 {
1113 int c;
1114
1115 if ( ! simple_cycle(s_cycle) ) {
1116 free_cycle(s_cycle);
1117 return(s_cycle_list);
1118 }
1119
1120 for ( c = 0; c < s_cycle_list->cnum; c++ )
1121 if ( same_cycle(s_cycle,s_cycle_list->list[c]) ) {
1122 free_cycle(s_cycle);
1123 return(s_cycle_list);
1124 }
1125
1126 s_cycle_list->list[s_cycle_list->cnum] = s_cycle;
1127 (s_cycle_list->cnum)++;
1128
1129 return(s_cycle_list);
1130 }
1131
free_cycle_list(cycle_list * s_cycle_list)1132 void free_cycle_list(cycle_list *s_cycle_list)
1133 {
1134 int c;
1135
1136 for ( c = 0; c < s_cycle_list->cnum; c++ )
1137 free_cycle(s_cycle_list->list[c]);
1138 free(s_cycle_list->list);
1139 free(s_cycle_list);
1140 }
1141
1142 #ifdef PRINT_CUTS
print_cycle_list(cycle_list * s_cycle_list)1143 void print_cycle_list(cycle_list *s_cycle_list)
1144 {
1145 int c;
1146
1147 printf("\n content of cycle_list: cnum = %d\n",s_cycle_list->cnum);
1148 for ( c = 0; c < s_cycle_list->cnum; c++ )
1149 print_cycle(s_cycle_list->list[c]);
1150 }
1151 #endif
1152
1153 /* get_shortest_odd_cycle_list: computation of the shortest odd cycles
1154 visiting a certain node in the separation graph, and each other
1155 possible intermediate node, using the auxiliary graph data structure
1156 for the shortest path computation - all the cycles in the list are
1157 different from each other */
1158
get_shortest_odd_cycle_list(int j,separation_graph * s_graph,auxiliary_graph * a_graph)1159 cycle_list *get_shortest_odd_cycle_list(
1160 int j, /* first node to be visited by the odd cycle */
1161 separation_graph *s_graph, /* current separation graph */
1162 auxiliary_graph *a_graph /* auxiliary graph for the shortest path computation */
1163 )
1164 {
1165 int source, sink, curr, pred, totedges, k, t, kt;
1166 double weight;
1167 #ifndef CGL_NEW_SHORT
1168 //node *source_ptr, *sink_ptr, *first_ptr;
1169 #else
1170 //cgl_node *source_ptr, *sink_ptr, *first_ptr;
1171 #endif
1172 edge *curr_edge;
1173 short_path_node *forw_arb, *backw_arb;
1174 cycle *s_cycle;
1175 cycle_list *s_cycle_list;
1176
1177 #ifdef TIME
1178 second_(&tsi);
1179 #endif
1180
1181 s_cycle_list = initialize_cycle_list((a_graph->nnodes)-2);
1182
1183 source = AG_TWIN1(j); sink = AG_TWIN2(j);
1184 //source_ptr = &(a_graph->nodes[source]);
1185 //sink_ptr = &(a_graph->nodes[sink]);
1186 //first_ptr = &(a_graph->nodes[0]);
1187
1188 /* compute the shortest path arborescence rooted at source and
1189 the shortest path arborescence rooted at sink (that comes for
1190 free due to symmetry) and store them (the path information is
1191 hidden into aux_graph) */
1192
1193 #ifdef TIME
1194 second_(&ti);
1195 #endif
1196 #ifndef CGL_NEW_SHORT
1197 {
1198 int nNodes = a_graph->nnodes;
1199 int nArcs = a_graph->narcs;
1200 cgl_arc * arcs = new cgl_arc [nArcs];
1201 for (int i=0;i<nArcs;i++) {
1202 arcs[i].length=a_graph->arcs[i].len;
1203 arcs[i].to=a_graph->arcs[i].head->index;
1204 }
1205 cgl_node * nodes = new cgl_node[nNodes+1];
1206 for (int i=0;i<nNodes;i++) {
1207 int iArc = a_graph->nodes[i].first-a_graph->arcs;
1208 nodes[i].firstArc=arcs+iArc;
1209 nodes[i].index=i;
1210 }
1211 int iArc = a_graph->nodes[nNodes].first-a_graph->arcs;
1212 nodes[nNodes].firstArc=arcs+iArc;
1213 cgl_graph graph;
1214 graph.nnodes=nNodes;
1215 graph.narcs=nArcs;
1216 graph.nodes=nodes;
1217 graph.arcs=arcs;
1218 cglShortestPath(&graph,source,ISCALE);
1219 dikbd(a_graph->nnodes,first_ptr,source_ptr,ISCALE);
1220 for ( k = 0; k < a_graph->nnodes; k++ ) {
1221 if ( a_graph->nodes[k].parent != NULL ) {
1222 int distance1 = a_graph->nodes[k].dist;
1223 int distance2 = graph.nodes[k].distanceBack;
1224 assert (distance1==distance2);
1225 } else {
1226 //
1227 printf("null parent %d\n",k);
1228 }
1229 }
1230 }
1231 #else
1232 cglShortestPath(a_graph,source,ISCALE);
1233 #endif
1234 #ifdef TIME
1235 second_(&tf);
1236 path_time += tf - ti;
1237 #endif
1238 forw_arb =
1239 reinterpret_cast<short_path_node *> (calloc(a_graph->nnodes,sizeof(short_path_node)));
1240 if ( forw_arb == NULL ) alloc_error(const_cast<char*>("forw_arb"));
1241 for ( k = 0; k < a_graph->nnodes; k++ ) {
1242 #ifndef CGL_NEW_SHORT
1243 if ( a_graph->nodes[k].parent != NULL ) {
1244 forw_arb[k].dist = a_graph->nodes[k].dist;
1245 forw_arb[k].pred = a_graph->nodes[k].parent->index;
1246 }
1247 #else
1248 if ( a_graph->nodes[k].parentNode >=0 ) {
1249 forw_arb[k].dist = a_graph->nodes[k].distanceBack;
1250 forw_arb[k].pred = a_graph->nodes[k].parentNode;
1251 }
1252 #endif
1253 else {
1254 forw_arb[k].dist = COIN_INT_MAX;
1255 forw_arb[k].pred = NONE;
1256 }
1257 }
1258 backw_arb =
1259 reinterpret_cast<short_path_node *> (calloc(a_graph->nnodes,sizeof(short_path_node)));
1260 if ( backw_arb == NULL ) alloc_error(const_cast<char*>("backw_arb"));
1261 for ( k = 0; k < a_graph->nnodes; k++ ) {
1262 #ifndef CGL_NEW_SHORT
1263 if ( a_graph->nodes[k].parent != NULL ) {
1264 backw_arb[AG_MATE(k)].dist = a_graph->nodes[k].dist;
1265 backw_arb[AG_MATE(k)].pred = AG_MATE(a_graph->nodes[k].parent->index);
1266 }
1267 #else
1268 if ( a_graph->nodes[k].parentNode >=0) {
1269 backw_arb[AG_MATE(k)].dist = a_graph->nodes[k].distanceBack;
1270 backw_arb[AG_MATE(k)].pred = AG_MATE(a_graph->nodes[k].parentNode);
1271 }
1272 #endif
1273 else {
1274 backw_arb[AG_MATE(k)].dist = COIN_INT_MAX;
1275 backw_arb[AG_MATE(k)].pred = NONE;
1276 }
1277 }
1278
1279 #ifdef USELESS
1280 /* compute second the shortest path anti-arborescence rooted at sink
1281 (which coincides with the arborescence since aux_graph is
1282 symmetrical) and store it */
1283
1284 #ifdef TIME
1285 second_(&ti);
1286 #endif
1287 cc = dikbd(a_graph->nnodes,first_ptr,sink_ptr,ISCALE);
1288 #ifdef TIME
1289 second_(&tf);
1290 path_time += tf - ti;
1291 #endif
1292 backw_arb =
1293 (short_path_node *) calloc(a_graph->nnodes,sizeof(short_path_node));
1294 if ( backw_arb == NULL ) alloc_error("backw_arb");
1295 for ( k = 0; k < a_graph->nnodes; k++ ) {
1296 backw_arb[k].dist = a_graph->nodes[k].dist;
1297 backw_arb[k].pred = a_graph->nodes[k].parent->index;
1298 }
1299 #endif
1300
1301 /* consider each possible intermediate node in aux_graph */
1302
1303 for ( k = 0; k < s_graph->nnodes; k++ ) {
1304 if ( k != j ) {
1305 for ( t = 1; t <= 2; t++ ) {
1306 if ( t == 1 ) kt = AG_TWIN1(k);
1307 else kt = AG_TWIN2(k);
1308 weight =
1309 (static_cast<double> (forw_arb[kt].dist + backw_arb[kt].dist)) /
1310 (static_cast<double> (ISCALE));
1311 if ( weight < MAX_CYCLE_WEIGHT + EPS ) {
1312 totedges = 0;
1313 /* count how many edges are in the forward path from source ... */
1314 curr = kt;
1315 do {
1316 if (curr<0) {
1317 totedges=-1;
1318 break;
1319 }
1320 curr = forw_arb[curr].pred; totedges++;
1321 } while ( curr != source );
1322 if (totedges>=0) {
1323 /* ... and in the backward path to sink */
1324 curr = kt;
1325 do {
1326 if (curr<0) {
1327 totedges=-1;
1328 break;
1329 }
1330 curr = backw_arb[curr].pred; totedges++;
1331 } while ( curr != sink );
1332 }
1333 if (totedges>0) {
1334 s_cycle = reinterpret_cast<cycle *> (calloc(1,sizeof(cycle)));
1335 if ( s_cycle == NULL ) alloc_error(const_cast<char*>("s_cycle"));
1336 s_cycle->weight = weight;
1337 s_cycle->length = totedges;
1338 s_cycle->edge_list = reinterpret_cast<edge **> (calloc(totedges,sizeof(edge *)));
1339 if ( s_cycle->edge_list == NULL ) alloc_error(const_cast<char*>("s_cycle->edge_list"));
1340 /* define the set of edges corresponding to the paths in sep_graph */
1341 totedges = 0;
1342 /* forward path from source ... */
1343 curr = kt;
1344 do {
1345 pred = forw_arb[curr].pred;
1346 if ( AG_TYPE(pred,curr) == EVEN )
1347 curr_edge = s_graph->even_adj_list
1348 [SG_EDGE_INDEX(s_graph,SG_ORIG(curr),SG_ORIG(pred))];
1349 else
1350 curr_edge = s_graph->odd_adj_list
1351 [SG_EDGE_INDEX(s_graph,SG_ORIG(curr),SG_ORIG(pred))];
1352 s_cycle->edge_list[totedges] = curr_edge;
1353 curr = pred;
1354 totedges++;
1355 } while ( curr != source );
1356 /* ... and backward path to sink */
1357 curr = kt;
1358 do {
1359 pred = backw_arb[curr].pred;
1360 if ( AG_TYPE(pred,curr) == EVEN )
1361 curr_edge = s_graph->even_adj_list
1362 [SG_EDGE_INDEX(s_graph,SG_ORIG(curr),SG_ORIG(pred))];
1363 else
1364 curr_edge = s_graph->odd_adj_list
1365 [SG_EDGE_INDEX(s_graph,SG_ORIG(curr),SG_ORIG(pred))];
1366 s_cycle->edge_list[totedges] = curr_edge;
1367 curr = pred;
1368 totedges++;
1369 } while ( curr != sink );
1370 /* insert the new cycle in the list */
1371 s_cycle_list = add_cycle_to_list(s_cycle,s_cycle_list);
1372 }
1373 }
1374 }
1375 }
1376 }
1377 free(forw_arb);
1378 free(backw_arb);
1379
1380 #ifdef TIME
1381 second_(&tsf);
1382 cycle_time += tsf - tsi;
1383 #endif
1384
1385 return(s_cycle_list);
1386 }
1387
1388 /* cut management subroutines */
1389
1390 /* initialize_cut_list: allocate and initialize the cut list data structure */
1391
initialize_cut_list(int max_cut)1392 cut_list *initialize_cut_list(int max_cut /* maximum number of cuts in the list */)
1393 {
1394 cut_list *cuts;
1395
1396 cuts = reinterpret_cast<cut_list *> (calloc(1,sizeof(cut_list)));
1397 if ( cuts == NULL ) alloc_error(const_cast<char*>("cuts"));
1398 cuts->cnum = 0;
1399 cuts->list = reinterpret_cast<cut **> (calloc(max_cut,sizeof(cut *)));
1400
1401 return(cuts);
1402 }
1403
1404 #ifdef PRINT_CUTS
print_cut(cut * v_cut)1405 void Cgl012Cut::print_cut(cut *v_cut)
1406 {
1407 printf("\n content of cut: n_of_constr = %d, cnzcnt = %d, crhs = %d, csense = %c, violation = %f\n",
1408 v_cut->n_of_constr,v_cut->cnzcnt,v_cut->crhs,v_cut->csense,v_cut->violation);
1409 print_int_vect(const_cast<char*>("cind"),v_cut->cind,v_cut->cnzcnt);
1410 print_int_vect(const_cast<char*>("cval"),v_cut->cval,v_cut->cnzcnt);
1411 if ( v_cut->constr_list != NULL )
1412 print_int_vect(const_cast<char*>("constr_list"),v_cut->constr_list,v_cut->n_of_constr);
1413 if ( v_cut->in_constr_list != NULL )
1414 print_short_int_vect(const_cast<char*>("in_constr_list"),v_cut->in_constr_list,inp_ilp->mr);
1415 ;
1416 }
1417
print_cut_list(cut_list * cuts)1418 void Cgl012Cut::print_cut_list(cut_list *cuts)
1419 {
1420 int c;
1421
1422 printf("\n content of cut_list: cnum = %d\n",cuts->cnum);
1423 for ( c = 0; c < cuts->cnum; c++ )
1424 print_cut(cuts->list[c]);
1425 }
1426 #endif
1427
free_cut(cut * v_cut)1428 void free_cut(cut *v_cut)
1429 {
1430 if ( v_cut->constr_list != NULL ) free(v_cut->constr_list);
1431 if ( v_cut->in_constr_list != NULL ) free(v_cut->in_constr_list);
1432 if ( v_cut->cind != NULL ) free(v_cut->cind);
1433 if ( v_cut->cval != NULL ) free(v_cut->cval);
1434 free(v_cut);
1435 }
1436
free_cut_list(cut_list * cuts)1437 void free_cut_list(cut_list *cuts)
1438 {
1439 int c;
1440
1441 for ( c = 0; c < cuts->cnum; c++ )
1442 if ( cuts->list[c] != NULL ) free_cut(cuts->list[c]);
1443 free(cuts->list);
1444 free(cuts);
1445 }
1446
1447 /* get_ori_cut_coef: get the coefficients of a cut, before dividing by 2 and
1448 rounding, starting from the list of the constraints combined to get
1449 the cut */
1450
get_ori_cut_coef(int n_of_constr,int * constr_list,int * ccoef,int * crhs,short int only_viol)1451 short int Cgl012Cut::get_ori_cut_coef(
1452 int n_of_constr, /* number of constraints combined */
1453 int *constr_list, /* list of the constraints combined */
1454 int *ccoef, /* cut left hand side coefficients */
1455 int *crhs, /* cut right hand side */
1456 short int only_viol /* flag which tells whether only an inequality of
1457 slack smaller than MAX_SLACK is of interest (TRUE)
1458 otherwise (FALSE) */
1459 )
1460 {
1461 int h, i, begi, gcdi, ofsj;
1462 double tot_slack;
1463
1464 /* fast check of the possible violation of the cut */
1465 if ( only_viol ) {
1466 tot_slack = 0.0;
1467 for ( h = 0; h < n_of_constr; h++ ) {
1468 tot_slack += p_ilp->slack[constr_list[h]];
1469 if ( tot_slack > MAX_SLACK - EPS ) return(FALSE);
1470 }
1471 }
1472
1473 //for ( j = 0; j < inp_ilp->mc; j++ )
1474 //ccoef[j] = 0;
1475 memset(ccoef,0,inp_ilp->mc*sizeof(int));
1476 (*crhs) = 0;
1477
1478 for ( h = 0; h < n_of_constr; h++ ) {
1479 i = constr_list[h];
1480 begi = inp_ilp->mtbeg[i]; gcdi = p_ilp->gcd[i];
1481 if ( inp_ilp->msense[i] != 'G' ) {
1482 if ( gcdi == 1 ) {
1483 for ( ofsj = 0; ofsj < inp_ilp->mtcnt[i]; ofsj++ )
1484 ccoef[inp_ilp->mtind[begi+ofsj]] += inp_ilp->mtval[begi+ofsj];
1485 (*crhs) += inp_ilp->mrhs[i];
1486 }
1487 else {
1488 for ( ofsj = 0; ofsj < inp_ilp->mtcnt[i]; ofsj++ )
1489 ccoef[inp_ilp->mtind[begi+ofsj]] += inp_ilp->mtval[begi+ofsj] / gcdi;
1490 (*crhs) += inp_ilp->mrhs[i] / gcdi;
1491 }
1492 }
1493 else {
1494 if ( gcdi == 1 ) {
1495 for ( ofsj = 0; ofsj < inp_ilp->mtcnt[i]; ofsj++ )
1496 ccoef[inp_ilp->mtind[begi+ofsj]] -= inp_ilp->mtval[begi+ofsj];
1497 (*crhs) -= inp_ilp->mrhs[i];
1498 }
1499 else {
1500 for ( ofsj = 0; ofsj < inp_ilp->mtcnt[i]; ofsj++ )
1501 ccoef[inp_ilp->mtind[begi+ofsj]] -= inp_ilp->mtval[begi+ofsj] / gcdi;
1502 (*crhs) -= inp_ilp->mrhs[i] / gcdi;
1503 }
1504 }
1505 }
1506
1507 return(TRUE);
1508 }
1509
1510 /* best_cut: find the coefficients, the rhs and the violation of the
1511 best possible cut that can be obtained by weakening a given set of
1512 coefficients to even and a rhs to odd, dividing by 2 and rounding */
1513
best_cut(int * ccoef,int * crhs,double * violation,short int update,short int only_viol)1514 short int Cgl012Cut::best_cut(
1515 int *ccoef, /* vector of the coefficients */
1516 int *crhs, /* pointer to rhs value */
1517 double *violation, /* violation of the cut */
1518 short int update, /* TRUE/FALSE: if TRUE, the new ccoef and crhs are
1519 given on output */
1520 short int only_viol /* flag which tells whether only an inequality of
1521 slack smaller than MAX_SLACK is of interest (TRUE)
1522 otherwise (FALSE) */
1523 )
1524 {
1525 int j, n_to_weak;
1526 short int original_parity;
1527 double original_slack, best_even_slack, best_odd_slack;
1528 int *vars_to_weak;
1529 info_weak *info_even_weak, *info_odd_weak;
1530
1531 /* choose the best weakening for the variables whose coefficient
1532 is not odd - this hopefully produces a stronger cut than that
1533 associated with the weakened inequalities to define the edges */
1534
1535 vars_to_weak = reinterpret_cast<int *> (calloc(inp_ilp->mc,sizeof(int)));
1536 if ( vars_to_weak == NULL ) alloc_error(const_cast<char*>("vars_to_weak"));
1537
1538 n_to_weak = 0;
1539 original_slack = 0.0;
1540 for ( j = 0; j < inp_ilp->mc; j++ ) {
1541 if ( ccoef[j] != 0 ) {
1542 if ( mod2(ccoef[j]) == ODD ) {
1543 vars_to_weak[n_to_weak] = j;
1544 n_to_weak++;
1545 }
1546 original_slack -= inp_ilp->xstar[j] * static_cast<double> (ccoef[j]);
1547 }
1548 }
1549 original_slack += static_cast<double> (*crhs);
1550 if ( original_slack > MAX_SLACK - EPS ) {
1551 free(vars_to_weak);
1552 return(FALSE);
1553 }
1554 original_parity = mod2(*crhs);
1555
1556 if ( best_weakening(n_to_weak,vars_to_weak,
1557 original_parity,original_slack,
1558 &best_even_slack,&best_odd_slack,
1559 &info_even_weak,&info_odd_weak,
1560 TRUE,only_viol) == ODD ) {
1561 *violation = ( 1.0 - best_odd_slack ) / 2.0;
1562 if ( ! update ) {
1563 /* new ccoef and rhs are not required on output */
1564 free(vars_to_weak);
1565 free_info_weak(info_odd_weak);
1566 return(TRUE);
1567 }
1568 /* update ccoef and crhs according to the best odd weakening */
1569 for ( j = 0; j < n_to_weak; j++ )
1570 if ( info_odd_weak->type[j] == LOWER_BOUND ) {
1571 ccoef[vars_to_weak[j]]--;
1572 *crhs -= inp_ilp->vlb[vars_to_weak[j]];
1573 }
1574 else {
1575 ccoef[vars_to_weak[j]]++;
1576 *crhs += inp_ilp->vub[vars_to_weak[j]];
1577 }
1578 /* compute and check the correctness of the cut coefficients */
1579 for ( j = 0; j < inp_ilp->mc; j++ ) {
1580 if ( mod2(ccoef[j]) == ODD ) {
1581 printf("!!! Error 2 in weakening a cut !!!\n");
1582 exit(0);
1583 }
1584 if ( ccoef[j] != 0 ) ccoef[j] /= 2;
1585 }
1586 if ( mod2(*crhs) == EVEN ) {
1587 printf("!!! Error 1 in weakening a cut !!!\n");
1588 exit(0);
1589 }
1590 *crhs = (*crhs - 1) / 2;
1591 free(vars_to_weak);
1592 free_info_weak(info_odd_weak);
1593 return(TRUE);
1594 }
1595 else {
1596 free(vars_to_weak);
1597 return(FALSE);
1598 }
1599 }
1600
1601 /* define_cut: construct a cut data structure from a vector of
1602 coefficients and a right-hand-side */
1603
define_cut(int * ccoef,int crhs)1604 cut *Cgl012Cut::define_cut(
1605 int *ccoef, /* coefficients of the cut */
1606 int crhs /* right hand side of the cut */
1607 )
1608 {
1609 int cnzcnt, j;
1610 cut *v_cut;
1611
1612 v_cut = reinterpret_cast<cut *> (calloc(1,sizeof(cut)));
1613 if ( v_cut == NULL ) alloc_error(const_cast<char*>("v_cut"));
1614 v_cut->crhs = crhs;
1615 cnzcnt = 0;
1616 for ( j = 0; j < inp_ilp->mc; j++ )
1617 if ( ccoef[j] != 0 ) cnzcnt++;
1618 v_cut->cnzcnt = cnzcnt;
1619 v_cut->csense = 'L';
1620 v_cut->cind = reinterpret_cast<int *> (calloc(cnzcnt,sizeof(int)));
1621 if ( v_cut->cind == NULL ) alloc_error(const_cast<char*>("v_cut->cind"));
1622 v_cut->cval = reinterpret_cast<int *> (calloc(cnzcnt,sizeof(int)));
1623 if ( v_cut->cval == NULL ) alloc_error(const_cast<char*>("v_cut->cval"));
1624 cnzcnt = 0; v_cut->violation = 0.0;
1625 for ( j = 0; j < inp_ilp->mc; j++ )
1626 if ( ccoef[j] != 0 ) {
1627 v_cut->cind[cnzcnt] = j;
1628 v_cut->cval[cnzcnt] = ccoef[j];
1629 v_cut->violation += inp_ilp->xstar[j] * static_cast<double> (ccoef[j]);
1630 cnzcnt++;
1631 }
1632 v_cut->violation -= static_cast<double> (crhs);
1633 return(v_cut);
1634 }
1635
1636 /* get_cut: extract a hopefully violated cut from an odd cycle of the
1637 separation graph */
1638
get_cut(cycle * s_cyc)1639 cut *Cgl012Cut::get_cut(
1640 cycle *s_cyc /* shortest odd cycles identified in the separation graph */
1641 )
1642 {
1643 int i, e, crhs;
1644 short int ok;
1645 /* short int original_parity; */
1646 double violation;
1647 /* double original_slack, best_even_slack, best_odd_slack; */
1648 int *ccoef /*, *vars_to_weak */ ;
1649 /* info_weak *info_even_weak, *info_odd_weak; */
1650 cut *v_cut;
1651 int ncomb;
1652 int *comb;
1653 short int *flag_comb;
1654 #ifndef CGGGGG
1655 static int iter = 0;
1656 static double gap, maxgap = 0.0;
1657 #endif
1658
1659 #ifdef TIME
1660 second_(&tsi);
1661 cut_ncalls++;
1662 #endif
1663
1664 /* compute the cut obtained by adding-up all the constraints
1665 corresponding to edges in the cycle, in their non-weak form */
1666
1667 #ifdef TIME
1668 second_(&tii);
1669 #endif
1670
1671 ccoef = reinterpret_cast<int *> (calloc(inp_ilp->mc,sizeof(int)));
1672 if ( ccoef == NULL ) alloc_error(const_cast<char*>("ccoef"));
1673 ncomb = 0;
1674 comb = reinterpret_cast<int *> (calloc(inp_ilp->mr,sizeof(int)));
1675 if ( comb == NULL ) alloc_error(const_cast<char*>("comb"));
1676 flag_comb = reinterpret_cast<short int *> (calloc(inp_ilp->mr,sizeof(short int)));
1677 if ( flag_comb == NULL ) alloc_error(const_cast<char*>("flag_comb"));
1678 #if 0
1679 // no need to as calloc used
1680 for ( i = 0; i < inp_ilp->mr; i++ ) flag_comb[i] = OUT;
1681
1682 for ( j = 0; j < inp_ilp->mc; j++ )
1683 ccoef[j] = 0;
1684 #endif
1685 crhs = 0;
1686 for ( e = 0; e < s_cyc->length; e++ ) {
1687 i = (s_cyc->edge_list[e])->constr;
1688 if ( i >= 0 && flag_comb[i] != IN) {
1689 /* the edge is not associated with a bound constraint */
1690 assert (ncomb<inp_ilp->mr);
1691 comb[ncomb] = i; ncomb++; flag_comb[i] = IN;
1692 }
1693 }
1694 ok = get_ori_cut_coef(ncomb,comb,ccoef,&crhs,TRUE);
1695
1696 #ifdef TIME
1697 second_(&tff);
1698 coef_time += tff - tii;
1699 #endif
1700
1701 ok = ok && best_cut(ccoef,&crhs,&violation,TRUE,TRUE);
1702 if ( ! ok ) {
1703 free(ccoef);
1704 free(comb);
1705 free(flag_comb);
1706 #ifdef TIME
1707 second_(&tsf);
1708 cut_time += tsf - tsi;
1709 #endif
1710 return(NULL);
1711 }
1712
1713 v_cut = define_cut(ccoef,crhs);
1714 iter++;
1715 if ( v_cut->violation > violation + EPS ||
1716 v_cut->violation < violation - EPS ) {
1717 //printf("Error in violation check\n");
1718 //printf("v_cut->violation %f violation %f gap %f maxgap (previous) %f\n",
1719 // v_cut->violation,violation,v_cut->violation-violation,maxgap);
1720 //printf("iter %d\n",iter);
1721 //exit(0);
1722 free_cut(v_cut);
1723 free(ccoef);
1724 free(comb);
1725 free(flag_comb);
1726 errorNo=1;
1727 return(NULL);
1728 }
1729 gap = v_cut->violation - violation;
1730 if ( gap < 0.0 ) gap = -gap;
1731 if ( gap > maxgap ) maxgap = gap;
1732 v_cut->n_of_constr = ncomb;
1733 v_cut->constr_list = comb;
1734 v_cut->in_constr_list = flag_comb;
1735
1736 free(ccoef);
1737
1738 #ifdef TIME
1739 second_(&tsf);
1740 cut_time += tsf - tsi;
1741 #endif
1742
1743 return(v_cut);
1744 }
1745
1746 /* cut_score: define the score of a (violated) cut */
1747
cut_score(int * ccoef,int crhs,double viol,short int only_viol)1748 double Cgl012Cut::cut_score(
1749 int *ccoef, /* cut left hand side coefficients */
1750 int crhs, /* cut right hand side */
1751 double viol, /* cut violation */
1752 short int only_viol /* flag which tells whether only an inequality of
1753 slack smaller than MAX_SLACK is of interest (TRUE)
1754 otherwise (FALSE) */
1755 )
1756 {
1757 int j, norm;
1758
1759 /* very simple score: violation divided/multiplied by the lhs norm */
1760 if ( only_viol && viol < MIN_VIOLATION ) return(-INF);
1761 norm = 0;
1762 for ( j = 0; j < p_ilp->mc; j++ ) {
1763 if ( ccoef[j] != 0 ) norm += ccoef[j] * ccoef[j];
1764 }
1765 if ( viol > 0.0 ) return (viol / sqrt(static_cast<double> (norm)));
1766 else return (viol * sqrt(static_cast<double> (norm)));
1767 }
1768
1769 /* same_cut: check whether two cuts are identical - not too clever
1770 (assumes the sparse coefficients are sorted by column index) */
1771
same_cut(cut * cut1,cut * cut2)1772 short int same_cut(cut *cut1, cut *cut2 /* cuts to be compared */)
1773 {
1774 int j;
1775
1776 if ( cut1->cnzcnt != cut2->cnzcnt ) return(FALSE);
1777 if ( cut1->crhs != cut2->crhs ) return(FALSE);
1778 if ( cut1->csense != cut2->csense ) return(FALSE);
1779 for ( j = 0; j < cut1->cnzcnt; j++ ) {
1780 if ( cut1->cind[j] != cut2->cind[j] ) return(FALSE);
1781 if ( cut1->cval[j] != cut2->cval[j] ) return(FALSE);
1782 }
1783 return(TRUE);
1784 }
1785
1786 /* add_cut_to_list: adds a cut to a list after checking that a copy of
1787 the same cut is not already in the list - no checking is made about
1788 cuts dominating each other or implied by other cuts in the list plus
1789 the constraints of the original problem */
1790
add_cut_to_list(cut * v_cut,cut_list * cuts)1791 cut_list *add_cut_to_list(
1792 cut *v_cut, /* pointer to the violated cut to be added to the list */
1793 cut_list *cuts /* input cut list to be updated */
1794 )
1795 {
1796 int c;
1797
1798 for ( c = 0; c < cuts->cnum; c++ ) {
1799 if ( same_cut(v_cut,cuts->list[c]) ) {
1800 free_cut(v_cut);
1801 return(cuts);
1802 }
1803 }
1804 cuts->list[cuts->cnum] = v_cut;
1805 cuts->cnum++;
1806 return(cuts);
1807 }
1808
1809 /* getcuts: pick the 0-1/2 cuts in the list and give them on output */
1810
getcuts(cut_list * cuts,int * cnum,int * cnzcnt,int ** cbeg,int ** ccnt,int ** cind,int ** cval,int ** crhs,char ** csense)1811 void getcuts(
1812 cut_list *cuts, /* input cut list */
1813 int *cnum, /* number of violated 0-1/2 cuts identified by the procedure */
1814 int *cnzcnt, /* overall number of nonzero's in the cuts */
1815 int **cbeg, /* starting position of each cut in arrays cind and cval */
1816 int **ccnt, /* number of entries of each cut in arrays cind and cval */
1817 int **cind, /* column indices of the nonzero entries of the cuts */
1818 int **cval, /* values of the nonzero entries of the cuts */
1819 int **crhs, /* right hand sides of the cuts */
1820 char **csense /* senses of the cuts: 'L', 'G' or 'E' */
1821 )
1822 {
1823 int i, ofsj, count;
1824 cut *cut_ptr;
1825
1826 /* allocate the memory for the output vectors */
1827
1828 (*cnum) = cuts->cnum;
1829 (*cnzcnt) = 0;
1830 for ( i = 0; i < cuts->cnum; i++ )
1831 (*cnzcnt) += (cuts->list[i])->cnzcnt;
1832 /* if ( (*cbeg) != NULL ) free(*cbeg); */
1833 (*cbeg) = reinterpret_cast<int *> (calloc((*cnum),sizeof(int)));
1834 if ( (*cbeg) == NULL ) alloc_error(const_cast<char*>("*cbeg"));
1835 /* if ( (*ccnt) != NULL ) free(*ccnt); */
1836 (*ccnt) = reinterpret_cast<int *> (calloc((*cnum),sizeof(int)));
1837 if ( (*ccnt) == NULL ) alloc_error(const_cast<char*>("*ccnt"));
1838 /* if ( (*crhs) != NULL ) free(*crhs); */
1839 (*crhs) = reinterpret_cast<int *> (calloc((*cnum),sizeof(int)));
1840 if ( (*crhs) == NULL ) alloc_error(const_cast<char*>("*crhs"));
1841 /* if ( (*csense) != NULL ) free(*csense); */
1842 (*csense) = reinterpret_cast<char *> (calloc((*cnum),sizeof(char)));
1843 if ( (*csense) == NULL ) alloc_error(const_cast<char*>("*csense"));
1844 /* if ( (*cind) != NULL ) free(*cind); */
1845 (*cind) = reinterpret_cast<int *> (calloc((*cnzcnt),sizeof(int)));
1846 if ( (*cind) == NULL ) alloc_error(const_cast<char*>("*cind"));
1847 /* if ( (*cval) != NULL ) free(*cval); */
1848 (*cval) = reinterpret_cast<int *> (calloc((*cnzcnt),sizeof(int)));
1849 if ( (*cval) == NULL ) alloc_error(const_cast<char*>("*cval"));
1850
1851 /* transfer the cuts information into the output data structures */
1852
1853 count = 0;
1854 for ( i = 0; i < cuts->cnum; i++ ) {
1855 cut_ptr = cuts->list[i];
1856 (*cbeg)[i] = count;
1857 (*ccnt)[i] = cut_ptr->cnzcnt;
1858 (*crhs)[i] = cut_ptr->crhs;
1859 (*csense)[i] = cut_ptr->csense;
1860 for ( ofsj = 0; ofsj < cut_ptr->cnzcnt; ofsj++ ) {
1861 (*cind)[count] = cut_ptr->cind[ofsj];
1862 (*cval)[count] = cut_ptr->cval[ofsj];
1863 count++;
1864 }
1865 }
1866 }
1867
1868 /* actual separation subroutines */
1869
1870 /* best_weakening: find the best upper/lower bound weakening of a set
1871 of variables */
1872
best_weakening(int n_to_weak,int * vars_to_weak,short int original_parity,double original_slack,double * best_even_slack,double * best_odd_slack,info_weak ** info_even_weak,info_weak ** info_odd_weak,short int only_odd,short int only_viol)1873 int Cgl012Cut::best_weakening(
1874 int n_to_weak, /* number of variables to weaken */
1875 int *vars_to_weak, /* indices of the variables to weaken */
1876 short int original_parity, /* original parity of the constraint to weaken */
1877 double original_slack, /* original slack of the constraint to weaken */
1878 double *best_even_slack, /* best possible slack of a weakened constraint
1879 with even right-hand-side */
1880 double *best_odd_slack, /* best possible slack of a weakened constraint
1881 with odd right-hand-side */
1882 info_weak **info_even_weak, /* weakening information about the best possible
1883 even weakened constraint */
1884 info_weak **info_odd_weak, /* weakening information about the best possible
1885 odd weakened constraint */
1886 short int only_odd, /* flag which tells whether only an odd weakening is of
1887 interest (TRUE) or both weakenings are (FALSE) */
1888 short int only_viol /* flag which tells whether only an inequality of
1889 slack smaller than MAX_SLACK is of interest (TRUE)
1890 otherwise (FALSE) */
1891 )
1892 {
1893 int nweak, cntweak, ofsl, l;
1894 short int flag_even, flag_odd, ok_even, ok_odd;
1895 double best_even_e, best_even_o, best_odd_e, best_odd_o;
1896 short int *type_even_weak, *type_odd_weak,
1897 *switch_even_weak, *switch_odd_weak;
1898
1899 #ifdef TIME
1900 second_(&tii);
1901 #endif
1902
1903 type_even_weak = reinterpret_cast<short int *> (calloc(p_ilp->mc,sizeof(short int)));
1904 if (type_even_weak == NULL ) alloc_error(const_cast<char*>("type_even_weak"));
1905 switch_even_weak = reinterpret_cast<short int *> (calloc(p_ilp->mc,sizeof(short int)));
1906 if (switch_even_weak == NULL ) alloc_error(const_cast<char*>("switch_even_weak"));
1907 type_odd_weak = reinterpret_cast<short int *> (calloc(p_ilp->mc,sizeof(short int)));
1908 if (type_odd_weak == NULL ) alloc_error(const_cast<char*>("type_odd_weak"));
1909 switch_odd_weak = reinterpret_cast<short int *> (calloc(p_ilp->mc,sizeof(short int)));
1910 if (switch_odd_weak == NULL ) alloc_error(const_cast<char*>("switch_odd_weak"));
1911
1912 if ( original_parity == EVEN ) {
1913 (*best_even_slack) = original_slack;
1914 (*best_odd_slack) = INF;
1915 }
1916 else {
1917 (*best_odd_slack) = original_slack;
1918 (*best_even_slack) = INF;
1919 }
1920 nweak = 0;
1921 for ( ofsl = 0; ofsl < n_to_weak; ofsl++ ) {
1922 l = vars_to_weak[nweak];
1923 if ( p_ilp->possible_weak[l] == NONE ) {
1924 free(type_even_weak); free(type_odd_weak);
1925 free(switch_even_weak); free(switch_odd_weak);
1926 #ifdef TIME
1927 second_(&tff);
1928 bw_time += tff - tii;
1929 #endif
1930 return(NONE);
1931 }
1932 else if ( p_ilp->possible_weak[l] == EVEN ) {
1933 /* only even weakening of l is possible */
1934 (*best_even_slack) += p_ilp->loss_even_weak[l];
1935 type_even_weak[nweak] = p_ilp->type_even_weak[l];
1936 switch_even_weak[nweak] = FALSE;
1937 (*best_odd_slack) += p_ilp->loss_even_weak[l];
1938 type_odd_weak[nweak] = p_ilp->type_even_weak[l];
1939 switch_odd_weak[nweak] = FALSE;
1940 }
1941 else if ( p_ilp->possible_weak[l] == ODD ) {
1942 /* only odd weakening of l is possible */
1943 best_even_e = (*best_even_slack);
1944 best_odd_o = (*best_odd_slack);
1945 (*best_even_slack) = best_odd_o + p_ilp->loss_odd_weak[l];
1946 type_even_weak[nweak] = p_ilp->type_odd_weak[l];
1947 switch_even_weak[nweak] = TRUE;
1948 (*best_odd_slack) = best_even_e + p_ilp->loss_odd_weak[l];
1949 type_odd_weak[nweak] = p_ilp->type_odd_weak[l];
1950 switch_odd_weak[nweak] = TRUE;
1951 }
1952 else {
1953 /* both weakenings of l are possible */
1954 best_even_e = (*best_even_slack) + p_ilp->loss_even_weak[l];
1955 best_even_o = (*best_odd_slack) + p_ilp->loss_odd_weak[l];
1956 best_odd_e = (*best_odd_slack) + p_ilp->loss_even_weak[l];
1957 best_odd_o = (*best_even_slack) + p_ilp->loss_odd_weak[l];
1958 if ( best_even_e <= best_even_o ) {
1959 (*best_even_slack) = best_even_e;
1960 type_even_weak[nweak] = p_ilp->type_even_weak[l];
1961 switch_even_weak[nweak] = FALSE;
1962 }
1963 else {
1964 (*best_even_slack) = best_even_o;
1965 type_even_weak[nweak] = p_ilp->type_odd_weak[l];
1966 switch_even_weak[nweak] = TRUE;
1967 }
1968 if ( best_odd_e <= best_odd_o ) {
1969 (*best_odd_slack) = best_odd_e;
1970 type_odd_weak[nweak] = p_ilp->type_even_weak[l];
1971 switch_odd_weak[nweak] = FALSE;
1972 }
1973 else {
1974 (*best_odd_slack) = best_odd_o;
1975 type_odd_weak[nweak] = p_ilp->type_odd_weak[l];
1976 switch_odd_weak[nweak] = TRUE;
1977 }
1978 }
1979 if ( ( only_viol ) &&
1980 ( (*best_even_slack) > MAX_SLACK - EPS ) &&
1981 ( (*best_odd_slack) > MAX_SLACK - EPS ) ) {
1982 free(type_even_weak); free(type_odd_weak);
1983 free(switch_even_weak); free(switch_odd_weak);
1984 #ifdef TIME
1985 second_(&tff);
1986 bw_time += tff - tii;
1987 #endif
1988 return(NONE);
1989 }
1990 nweak++;
1991 }
1992
1993 /* construct the weakening vectors associated with the best
1994 even and odd pairs (if the associated slack is not too big) */
1995
1996 if ( (! only_odd) &&
1997 ( ( (*best_even_slack) <= MAX_SLACK - EPS ) ||
1998 ( (! only_viol) && (*best_even_slack) <= INF - EPS ) ) ) {
1999 ok_even = TRUE;
2000 (*info_even_weak) = alloc_info_weak(nweak);
2001 (*info_even_weak)->nweak = nweak;
2002 flag_even = EVEN;
2003 cntweak = nweak;
2004 for ( ofsl = n_to_weak - 1; ofsl >= 0; ofsl-- ) {
2005 cntweak--;
2006 (*info_even_weak)->var[cntweak] = vars_to_weak[ofsl];
2007 if ( flag_even == EVEN ) {
2008 (*info_even_weak)->type[cntweak] = type_even_weak[cntweak];
2009 if ( switch_even_weak[cntweak] ) flag_even = ODD;
2010 }
2011 else {
2012 (*info_even_weak)->type[cntweak] = type_odd_weak[cntweak];
2013 if ( switch_odd_weak[cntweak] ) flag_even = EVEN;
2014 }
2015 }
2016 }
2017 else ok_even = FALSE;
2018
2019 if ( ( (*best_odd_slack) <= MAX_SLACK - EPS ) ||
2020 ( (! only_viol) && (*best_odd_slack) <= INF - EPS ) ) {
2021 ok_odd = TRUE;
2022 (*info_odd_weak) = alloc_info_weak(nweak);
2023 (*info_odd_weak)->nweak = nweak;
2024 flag_odd = ODD;
2025 cntweak = nweak;
2026 for ( ofsl = n_to_weak - 1; ofsl >= 0; ofsl-- ) {
2027 cntweak--;
2028 (*info_odd_weak)->var[cntweak] = vars_to_weak[ofsl];
2029 if ( flag_odd == EVEN ) {
2030 (*info_odd_weak)->type[cntweak] = type_even_weak[cntweak];
2031 if ( switch_even_weak[cntweak] ) flag_odd = ODD;
2032 }
2033 else {
2034 (*info_odd_weak)->type[cntweak] = type_odd_weak[cntweak];
2035 if ( switch_odd_weak[cntweak] ) flag_odd = EVEN;
2036 }
2037 }
2038 }
2039 else ok_odd = FALSE;
2040
2041 free(type_even_weak); free(type_odd_weak);
2042 free(switch_even_weak); free(switch_odd_weak);
2043 #ifdef TIME
2044 second_(&tff);
2045 bw_time += tff - tii;
2046 #endif
2047
2048 if ( ok_odd && ok_even ) return(BOTH);
2049 if ( ok_even ) return(EVEN);
2050 if ( ok_odd ) return(ODD);
2051 return(NONE);
2052 }
2053
2054 /* basic_separation: try to identify violated 0-1/2 cuts by using the
2055 original procedure described in Caprara and Fischetti's MP paper */
2056
basic_separation()2057 cut_list *Cgl012Cut::basic_separation()
2058 {
2059 int i, j, k, l, begi, special, ofsj, ofsk, ofsl, n_to_weak, c;
2060 short int parity, original_parity, ok_weak;
2061 double weight, original_slack, best_even_slack, best_odd_slack;
2062 int *vars_to_weak;
2063 info_weak *info_even_weak, *info_odd_weak, *i_weak;
2064 separation_graph *sep_graph;
2065 auxiliary_graph *aux_graph;
2066 cycle_list *short_cycle_list;
2067 cut *violated_cut;
2068 cut_list *out_cuts;
2069
2070 /* construct the separation graph by the standard weakening procedure */
2071
2072 #ifdef PRINT
2073 print_parity_ilp();
2074 #endif
2075
2076 #ifdef TIME
2077 second_(&td);
2078 #endif
2079 #ifdef PRINT_TIME
2080 printf("... time elapsed at the beginning of basic_separation: %f\n",td - tti);
2081 #endif
2082 #ifdef TIME
2083 second_(&td);
2084 #endif
2085 #ifdef PRINT_TIME
2086 printf("... time elapsed before initialize_sep_graph: %f\n",td - tti);
2087 #endif
2088 sep_graph = initialize_sep_graph();
2089 special = p_ilp->mc;
2090 #ifdef TIME
2091 second_(&td);
2092 #endif
2093 #ifdef PRINT_TIME
2094 printf("... time elapsed before weakening: %f\n",td - tti);
2095 #endif
2096
2097 /* edges associated with actual constraints in the ILP */
2098
2099 for ( i = 0; i < p_ilp->mr; i++ ) {
2100 if ( ! p_ilp->row_to_delete[i] ) {
2101 begi = p_ilp->mtbeg[i];
2102 if ( p_ilp->mtcnt[i] == 1 ) {
2103 /* row with one odd entry only: edge j -- special */
2104 weight = p_ilp->slack[i];
2105 if ( weight < MAX_SLACK - EPS ) {
2106 j = p_ilp->mtind[begi];
2107 parity = p_ilp->mrhs[i];
2108 i_weak = alloc_info_weak(0);
2109 sep_graph = update_weight_sep_graph
2110 (j,special,weight,parity,i,i_weak,sep_graph);
2111 }
2112 }
2113 else if ( p_ilp->mtcnt[i] == 2 ) {
2114 /* row with two odd entries only: edge j -- k */
2115 weight = p_ilp->slack[i];
2116 if ( weight < MAX_SLACK - EPS ) {
2117 j = p_ilp->mtind[begi];
2118 k = p_ilp->mtind[begi+1];
2119 parity = p_ilp->mrhs[i];
2120 i_weak = alloc_info_weak(0);
2121 sep_graph = update_weight_sep_graph
2122 (j,k,weight,parity,i,i_weak,sep_graph);
2123 }
2124 }
2125 else {
2126 /* row with three or more odd entries: weakening for all 1's pairs */
2127 for ( ofsj = 0; ofsj < p_ilp->mtcnt[i]; ofsj++ ) {
2128 for ( ofsk = ofsj + 1; ofsk < p_ilp->mtcnt[i]; ofsk++ ) {
2129 /* edge(s) j -- k */
2130 j = p_ilp->mtind[begi+ofsj];
2131 k = p_ilp->mtind[begi+ofsk];
2132 original_slack = p_ilp->slack[i];
2133 original_parity = p_ilp->mrhs[i];
2134 n_to_weak = 0;
2135 vars_to_weak = reinterpret_cast<int *> (calloc(inp_ilp->mc,sizeof(int)));
2136 if ( vars_to_weak == NULL ) alloc_error(const_cast<char*>("vars_to_weak"));
2137 for ( ofsl = 0; ofsl < p_ilp->mtcnt[i]; ofsl++ )
2138 if ( ofsl != ofsj && ofsl != ofsk ) {
2139 l = p_ilp->mtind[begi+ofsl];
2140 vars_to_weak[n_to_weak] = l;
2141 n_to_weak++;
2142 }
2143 ok_weak = best_weakening(n_to_weak,vars_to_weak,
2144 original_parity,original_slack,
2145 &best_even_slack,&best_odd_slack,
2146 &info_even_weak,&info_odd_weak,
2147 FALSE,TRUE);
2148 free(vars_to_weak);
2149 if ( ok_weak == NONE ) goto EXITJK;
2150 if ( ok_weak == BOTH || ok_weak == EVEN ) {
2151 if ( best_even_slack < MAX_SLACK - EPS ) {
2152 weight = best_even_slack; parity = EVEN;
2153 sep_graph = update_weight_sep_graph
2154 (j,k,weight,parity,i,info_even_weak,sep_graph);
2155 }
2156 }
2157 if ( ok_weak == BOTH || ok_weak == ODD ) {
2158 if ( best_odd_slack < MAX_SLACK - EPS ) {
2159 weight = best_odd_slack; parity = ODD;
2160 sep_graph = update_weight_sep_graph
2161 (j,k,weight,parity,i,info_odd_weak,sep_graph);
2162 }
2163 }
2164 EXITJK:; }
2165 }
2166 }
2167 }
2168 }
2169
2170 /* edges associated with the bound constraints (probably useless
2171 but necessary in some cases */
2172
2173 for ( j = 0; j < p_ilp->mc; j++ ) {
2174 if ( ! p_ilp->col_to_delete[j] ) {
2175 weight = p_ilp->xstar[j] - inp_ilp->vlb[j];
2176 if ( weight < MAX_SLACK - EPS ) {
2177 parity = mod2(inp_ilp->vlb[j]);
2178 i_weak = alloc_info_weak(0);
2179 sep_graph = update_weight_sep_graph
2180 (j,special,weight,parity,NONE,i_weak,sep_graph);
2181 }
2182 weight = inp_ilp->vub[j] - p_ilp->xstar[j];
2183 if ( weight < MAX_SLACK - EPS ) {
2184 parity = mod2(inp_ilp->vub[j]);
2185 i_weak = alloc_info_weak(0);
2186 sep_graph = update_weight_sep_graph
2187 (j,special,weight,parity,NONE,i_weak,sep_graph);
2188 }
2189 }
2190 }
2191 #ifdef TIME
2192 second_(&tf);
2193 weak_time += tf - ti;
2194 #endif
2195
2196 /* construct the auxiliary graph for the shortest path computation
2197 and compute the smallest cost odd cycle visiting each node -
2198 this part is strongly dependent on the data structure used by
2199 the shortest path subroutine used */
2200
2201 #ifdef PRINT
2202 print_sep_graph(sep_graph);
2203 #endif
2204
2205 #ifdef TIME
2206 second_(&ti);
2207 #endif
2208 #ifdef TIME
2209 second_(&td);
2210 #endif
2211 #ifdef PRINT_TIME
2212 printf("... time elapsed before define_aux_graph: %f\n",td - tti);
2213 #endif
2214 aux_graph = define_aux_graph(sep_graph);
2215 #ifdef TIME
2216 second_(&tf);
2217 aux_time += tf - ti;
2218 #endif
2219 #ifdef PRINT
2220 print_aux_graph(aux_graph);
2221 #endif
2222
2223 /* exit(1); */
2224
2225 #ifdef TIME
2226 second_(&td);
2227 #endif
2228 #ifdef PRINT_TIME
2229 printf("... time elapsed before cycles and cuts: %f\n",td - tti);
2230 printf("%d nodes on list\n",sep_graph->nnodes);
2231 #endif
2232 out_cuts = initialize_cut_list(MAX_CUTS);
2233 for ( j = 0; j < sep_graph->nnodes; j++ ) {
2234 short_cycle_list = get_shortest_odd_cycle_list(j,sep_graph,aux_graph);
2235 if ( short_cycle_list == NULL ) goto EXIT_NODE;
2236 #ifdef PRINT
2237 print_cycle_list(short_cycle_list);
2238 #endif
2239 for ( c = 0; c < short_cycle_list->cnum; c++ ) {
2240 violated_cut = get_cut(short_cycle_list->list[c]);
2241 if ( violated_cut == NULL ) {
2242 if (!errorNo)
2243 continue;
2244 else
2245 break;
2246 }
2247 #ifdef PRINT
2248 print_cut(violated_cut);
2249 #endif
2250 if ( violated_cut->violation > MIN_VIOLATION + EPS ) {
2251 /* violated 0-1/2 cut found */
2252 out_cuts = add_cut_to_list(violated_cut,out_cuts);
2253 if ( out_cuts->cnum >= MAX_CUTS ) {
2254 free_cycle_list(short_cycle_list);
2255 goto EXIT_CUTS;
2256 }
2257 }
2258 else free_cut(violated_cut);
2259 }
2260 /* remove the current node from the auxiliary graph */
2261 EXIT_NODE:
2262 aux_graph = cancel_node_aux_graph(j,aux_graph);
2263 free_cycle_list(short_cycle_list);
2264 }
2265
2266 EXIT_CUTS:
2267 free_sep_graph(sep_graph);
2268 free_aux_graph(aux_graph);
2269
2270 #ifdef PRINT_CUTS
2271 print_cut_list(out_cuts);
2272 #endif
2273 #ifdef TIME
2274 second_(&td);
2275 #endif
2276 #ifdef PRINT_TIME
2277 printf("... time elapsed at the end of basic_separation: %f\n",td - tti);
2278 #endif
2279
2280 return(out_cuts);
2281 }
2282
2283 /*
2284 012cut: main procedure for 0-1/2 cut separation
2285 first release: Aug 12 1996
2286 last revision: Jun 10 1997
2287 */
2288
2289 /* static data structures for log information about separation */
2290
2291 #ifndef CGGGGG
2292 static int sep_iter = 0; /* number of the current separation iteration */
2293 //#define POOL
2294 #ifdef POOL
2295 static pool_cut_list *pool = NULL; /* information about the cuts separated
2296 so far, used to decide when they should
2297 be added to the current LP */
2298 #endif
2299 static log_var **vlog = NULL; /* information about the value attained
2300 by the variables in the last iterations,
2301 used to possibly set to 0 some coefficient
2302 > 0 in a cut to be added */
2303 static bool aggr; /* flag saying whether as many cuts as possible are required
2304 from the separation procedure (TRUE) or not (FALSE) */
2305 #endif
2306 /* include the reactive local search heuristic */
2307
2308 //was #include "Cgltabu_012.c"
2309 //start include "Cgltabu_012.c"
2310
2311 #define MAX_TABU_ITER 100
2312 #define NUM_HASH_ENTRIES MAX_CUT_POOL
2313 /* initial length of the tabu list */
2314 #define IN_PROHIB_PERIOD 3
2315 #define MAX_TIME_FACTOR 3
2316
2317 /* data structure for the current local search solution */
2318
2319 typedef struct {
2320 int n_of_constr; /* number of constraints in the current cut */
2321 short int *in_constr_list; /* flag saying whether a given constraint is
2322 in the list of constraints of the cut (IN)
2323 or not (OUT) */
2324 int *non_weak_coef; /* coefficients of the cut before weakening */
2325 int non_weak_rhs; /* coefficient of the rhs before weakening */
2326 double slack_sum; /* sum of the slacks of the constraints in the cut */
2327 double min_weak_loss; /* minimum loss by weakening the non even
2328 coefficients */
2329 int one_norm; /* 1-norm of the lhs, i.e. sum of the absolute values of
2330 the coefficients */
2331 short int ok; /* logical flag telling whether the cut could be weakened
2332 to a 0-1/2 cut or not - if false the two fields below
2333 have no meaning */
2334 int *coef; /* actual coefficients of the cut */
2335 int rhs; /* actual rhs of the cut */
2336 double violation; /* violation of the cut */
2337 } tabu_cut;
2338
2339 /* data structure for the hash table used in memory reaction */
2340
2341 typedef struct h_e {
2342 int n_of_el; /* number of components to be considered */
2343 short int *flag_vect; /* vector of flags for the components */
2344 int last_vis; /* last iteration when this element was visited */
2345 struct h_e *next; /* pointer to the next element in the hash chain */
2346 } hash_element;
2347
2348 typedef hash_element **hash_table;
2349
2350 /* global variables for local search */
2351
2352 static int n; /* number of variables in the ILP */
2353 static int m; /* number of constraints in the ILP */
2354 static int it; /* number of tabu search iterations so far */
2355 static tabu_cut *cur_cut; /* information about the current cut in local search */
2356 static int *last_moved; /* last iteration when a given constraint was added/
2357 deleted from the list of constraints of the cut */
2358 static int last_it_add; /* last iteration when a cut was added to the list */
2359 static int last_it_restart; /* last iteration when a restart was performed */
2360 static int prohib_period; /* current prohibition period */
2361 static int last_prohib_period_mod; /* last iteration where prohibition period was modified */
2362 static hash_table hash_tab; /* hash table */
2363 static int A; /* parameter A in Battiti and Protasi */
2364 static int B; /* parameter B in Battiti and Protasi */
2365 static float elapsed_time; /* time elapsed since the beginning of the current
2366 tabu search call */
2367
2368 /* clear_cur_cut: clear the current solution (no constraint in the cut) */
2369
clear_cur_cut()2370 void clear_cur_cut()
2371 {
2372 int i, j;
2373
2374 cur_cut->n_of_constr = 0;
2375 cur_cut->rhs = 0;
2376 cur_cut->non_weak_rhs = 0;
2377 cur_cut->violation = 0.0;
2378 cur_cut->slack_sum = 0.0;
2379 cur_cut->min_weak_loss = 0.0;
2380 cur_cut->one_norm = 0;
2381 for ( j = 0; j < n; j++ ) {
2382 cur_cut->coef[j] = 0;
2383 cur_cut->non_weak_coef[j] = 0;
2384 }
2385 for ( i = 0; i < m; i++ ) {
2386 cur_cut->in_constr_list[i] = OUT;
2387 }
2388 cur_cut->ok = FALSE;
2389 }
2390
2391 /* initialize_cur_cut: allocate the memory for cur_cut */
2392
initialize_cur_cut()2393 void initialize_cur_cut()
2394 {
2395 cur_cut = reinterpret_cast<tabu_cut *> (calloc(1,sizeof(tabu_cut)));
2396 if ( cur_cut == NULL ) alloc_error(const_cast<char*>("cur_cut"));
2397 cur_cut->coef = reinterpret_cast<int *> (calloc(n,sizeof(int)));
2398 if ( cur_cut->coef == NULL ) alloc_error(const_cast<char*>("cur_cut->coef"));
2399 cur_cut->non_weak_coef = reinterpret_cast<int *> (calloc(n,sizeof(int)));
2400 if ( cur_cut->non_weak_coef == NULL ) alloc_error(const_cast<char*>("cur_cut->non_weak_coef"));
2401 cur_cut->in_constr_list = reinterpret_cast<short int *> (calloc(m,sizeof(short int)));
2402 if ( cur_cut->in_constr_list == NULL ) alloc_error(const_cast<char*>("cur_cut->in_constr_list"));
2403 clear_cur_cut();
2404 }
2405
2406 /* free_cur_cut: free the memory for cur_cut */
2407
free_cur_cut()2408 void free_cur_cut()
2409 {
2410 free(cur_cut->coef);
2411 free(cur_cut->non_weak_coef);
2412 free(cur_cut->in_constr_list);
2413 free(cur_cut);
2414 }
2415
2416 #ifdef PRINT_TABU
2417 /* print_cur_cut: display cur_cut on output */
2418
print_cur_cut()2419 void Cgl012Cut::print_cur_cut()
2420 {
2421 int i, j;
2422
2423 printf("iteration %d prohib_period %d\n",it,prohib_period);
2424 printf("\n content of cur_cut data structure: n_of_constr = %d, ok = %d\n", cur_cut->n_of_constr, cur_cut->ok);
2425 for ( i = 0; i < m; i++ )
2426 if ( cur_cut->in_constr_list[i] == IN )
2427 printf("constr. %d\n",i);
2428 /*
2429 printf(" list of constraints:\n");
2430 for ( i = 0; i < m; i++ )
2431 if ( cur_cut->in_constr_list[i] == IN )
2432 print_constr(i);
2433 print_int_vect("non_weak_coef",cur_cut->non_weak_coef,n);
2434 printf(" non_weak_rhs = %d\n",cur_cut->non_weak_rhs);
2435 print_int_vect("coef",cur_cut->coef,n);
2436 printf(" rhs = %d\n",cur_cut->rhs);
2437 */
2438 for ( j = 0 /* , viol = - (double) cur_cut->rhs */ ; j < n; j++ )
2439 if ( ( p_ilp->xstar[j] > ZERO || p_ilp->xstar[j] < -ZERO ) && cur_cut->non_weak_coef[j] != 0 ) {
2440 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]);
2441 /* viol += p_ilp->xstar[j] * cur_cut->coef[j]; */
2442 }
2443 printf("rhs %d viol %f slack_sum %f min_weak_loss %f one_norm %d\n",
2444 cur_cut->rhs, cur_cut->violation, cur_cut->slack_sum,
2445 cur_cut->min_weak_loss, cur_cut->one_norm);
2446 }
2447 #endif
2448 /* same_short_vect: check whether two short int vectors have the same content */
2449
same_short_vect(int n_of_el,short int * vec_1,short int * vec_2)2450 short int same_short_vect(
2451 int n_of_el, /* number of components in the vectors */
2452 short int *vec_1,
2453 short int *vec_2 /* vectors to be checked */
2454 )
2455 {
2456 int i;
2457 for ( i = 0; i < n_of_el; i++ )
2458 if ( vec_1[i] != vec_2[i] ) return(FALSE);
2459 return(TRUE);
2460 }
2461
2462 /* initialize_hash_table: allocate the memory for the hash table */
2463
initialize_hash_table()2464 void initialize_hash_table()
2465 {
2466 int i;
2467 hash_tab = reinterpret_cast<hash_element **> (calloc(NUM_HASH_ENTRIES,sizeof(hash_element *)));
2468 if ( hash_tab == NULL ) alloc_error(const_cast<char*>("hash_tab"));
2469 for ( i = 0; i < NUM_HASH_ENTRIES; i++ ) hash_tab[i] = NULL;
2470 }
2471
2472 /* clear_hash_table: clear the current hash table */
2473
clear_hash_table()2474 void clear_hash_table()
2475 {
2476 int i;
2477 hash_element *hash_ptr, *hash_el;
2478
2479 for ( i = 0; i < NUM_HASH_ENTRIES; i++ ) {
2480 if ( hash_tab[i] != NULL ) {
2481 hash_ptr = hash_tab[i];
2482 do {
2483 hash_el = hash_ptr->next;
2484 free(hash_ptr->flag_vect);
2485 free(hash_ptr);
2486 hash_ptr = hash_el;
2487 } while ( hash_ptr != NULL );
2488 hash_tab[i] = NULL;
2489 }
2490 }
2491 }
2492
2493 /* free_hash_table: deallocate the memory for the hash table */
2494
free_hash_table()2495 void free_hash_table()
2496 {
2497 clear_hash_table();
2498 free(hash_tab);
2499 }
2500
2501 /* hash_addr: compute the hash address associated with the current cut */
2502
hash_addr(int n_of_el,short int * flag_vect)2503 int hash_addr(
2504 int n_of_el, /* number of elements to be considered */
2505 short int *flag_vect /* vector of flags for the elements */
2506 )
2507 {
2508 int i, addr;
2509
2510 /* very simple algorithm: just add-up the squared indices of the IN elements */
2511 addr = 0;
2512 for ( i = 0; i < n_of_el; i++ )
2513 if ( flag_vect[i] == IN ) addr += i * i;
2514 return(addr % NUM_HASH_ENTRIES);
2515 }
2516
2517 /* hash_search: search for the current cut in the hash list of all cuts -
2518 if found return TRUE and update the last iteration the cut was found */
2519
hash_search(int * cyc_len)2520 short int hash_search(int *cyc_len /* length of the cycle if the current cut is found */)
2521 {
2522 int addr;
2523 hash_element *hash_el;
2524
2525 addr = hash_addr(m,cur_cut->in_constr_list);
2526 hash_el = hash_tab[addr];
2527 while ( hash_el != NULL ) {
2528 if ( same_short_vect(m,cur_cut->in_constr_list,hash_el->flag_vect) ) {
2529 *cyc_len = it - hash_el->last_vis;
2530 hash_el->last_vis = it;
2531 return(TRUE);
2532 }
2533 hash_el = hash_el->next;
2534 }
2535 return(FALSE);
2536 }
2537
2538 /* hash_insert: insert a new cut in the hash list of all cuts */
2539
hash_insert()2540 void hash_insert()
2541 {
2542 int addr, i;
2543 hash_element *hash_el, *hash_ptr;
2544
2545 addr = hash_addr(m,cur_cut->in_constr_list);
2546 hash_el = reinterpret_cast<hash_element *> (calloc(1,sizeof(hash_element)));
2547 if ( hash_el == NULL ) alloc_error(const_cast<char*>("hash_el"));
2548 hash_el->n_of_el = m;
2549 hash_el->last_vis = it;
2550 hash_el->next = NULL;
2551 hash_el->flag_vect = reinterpret_cast<short int *> (calloc(m,sizeof(short int)));
2552 if ( hash_el->flag_vect == NULL ) alloc_error(const_cast<char*>("hash_el->flag_vect"));
2553 for ( i = 0; i < m; i++ )
2554 hash_el->flag_vect[i] = cur_cut->in_constr_list[i];
2555 if ( hash_tab[addr] == NULL )
2556 hash_tab[addr] = hash_el;
2557 else {
2558 hash_ptr = hash_tab[addr];
2559 while ( hash_ptr->next != NULL ) {
2560 #if 0
2561 /* this check can be omitted to save time */
2562 if ( same_short_vect(m,cur_cut->in_constr_list,hash_ptr->flag_vect) ) {
2563 printf("attempt to insert in the hash an already present cut\n");
2564 exit(0);
2565 }
2566 #endif
2567 hash_ptr = hash_ptr->next;
2568 }
2569 hash_ptr->next = hash_el;
2570 }
2571 }
2572
2573 /* increase_prohib_period: implemented as in Battiti and Protasi */
2574
increase_prohib_period()2575 void increase_prohib_period()
2576 {
2577 if ( prohib_period * 1.1 > prohib_period + 1 )
2578 if ( prohib_period * 1.1 < m - 2 ) prohib_period =
2579 static_cast<int> (prohib_period*1.1);
2580 else prohib_period = m - 2;
2581 else
2582 if ( prohib_period + 1 < m - 2 ) prohib_period += 1;
2583 else prohib_period = m - 2;
2584 last_prohib_period_mod = it;
2585 }
2586
2587 /* decrease_prohib_period: implemented as in Battiti and Protasi */
2588
decrease_prohib_period()2589 void decrease_prohib_period()
2590 {
2591 if ( prohib_period * 0.9 < prohib_period - 1 )
2592 if ( prohib_period * 0.9 > IN_PROHIB_PERIOD ) prohib_period = static_cast<int> (prohib_period* 0.9);
2593 else prohib_period = IN_PROHIB_PERIOD;
2594 else
2595 if ( prohib_period - 1 > IN_PROHIB_PERIOD ) prohib_period -= 1;
2596 else prohib_period = IN_PROHIB_PERIOD;
2597 last_prohib_period_mod = it;
2598 }
2599
2600 /* allowed: check if moving (adding/deleting) a given constraint
2601 is not a tabu move */
2602
allowed(int i)2603 short int allowed(int i /* constraint to be checked */)
2604 {
2605 if ( last_moved[i] < it - prohib_period ) {
2606 if ( cur_cut->in_constr_list[i] == IN ) {
2607 if ( cur_cut->n_of_constr > 1 ) return(TRUE);
2608 else return(FALSE);
2609 }
2610 else {
2611 if ( cur_cut->n_of_constr < m - 1 ) return(TRUE);
2612 else return(FALSE);
2613 }
2614 }
2615 else return(FALSE);
2616 }
2617
2618 /* in_cur_cut: check whether a given constraint is in the list of
2619 constraints defining the current cut */
2620
in_cur_cut(int i)2621 short int in_cur_cut(int i /* constraint to be checked */)
2622 {
2623 if ( cur_cut->in_constr_list[i] == OUT ) return(FALSE);
2624 else return(TRUE);
2625 }
2626
2627 /* tabu_score: define the score of a potential new cut */
2628
tabu_score(int * ccoef,int crhs,double viol,double norm)2629 double tabu_score(
2630 int *ccoef, /* cut left hand side coefficients */
2631 int crhs, /* cut right hand side */
2632 double viol, /* cut violation */
2633 double norm /* cut norm - 1-norm is used below */
2634 )
2635 {
2636 /* very simple score: violation divided/multiplied by the lhs 1-norm */
2637
2638 if ( norm == 0 ) norm = 1;
2639 if ( viol > 0.0 ) return (viol / norm);
2640 else return (viol * norm);
2641 }
2642
2643 /* score_by_moving: compute the score of the best cut obtainable from
2644 the current local search solution by inserting/deleting a constraint */
2645
score_by_moving(int i,short int itype,double thresh)2646 double Cgl012Cut::score_by_moving(
2647 int i, /* constraint to be moved */
2648 short int itype, /* type of move - ADD or DEL */
2649 double thresh /* minimum value of an interesting score */
2650 )
2651 {
2652 #define PENALTY_NON_WEAKABLE -1.0
2653 #define FAST_SCORE_EVAL 1
2654 int j, begi, gcdi, ofsj, ij, crhs, one_norm, support_inter;
2655 short int flag_gt, ok;
2656 double slack_sum, weak_loss, score, viol;
2657 int *new_coef, *ccoef;
2658 begi = inp_ilp->mtbeg[i];
2659 gcdi = p_ilp->gcd[i];
2660
2661 /* fast check - optimistic evaluation of the score */
2662
2663 slack_sum = cur_cut->slack_sum;
2664
2665 if ( itype == ADD ) slack_sum += p_ilp->slack[i] / static_cast<double> (gcdi);
2666 else slack_sum -= p_ilp->slack[i] / static_cast<double> (gcdi);
2667
2668 viol = ( 1.0 - slack_sum ) / 2.0;
2669
2670 score = tabu_score(NULL,0,viol,1.0);
2671 /*
2672 printf("Score estimate 1 %f Threshold %f\n",score,thresh);
2673 */
2674 if ( score < thresh + ZERO ) {
2675 return (score);
2676 }
2677
2678 /* discard the cuts that have empty support intersection with the
2679 current one */
2680
2681 support_inter = 0;
2682
2683 for ( ofsj = 0, ij = begi; ofsj < inp_ilp->mtcnt[i]; ofsj++, ij++ )
2684 if ( cur_cut->non_weak_coef[inp_ilp->mtind[ij]] !=0 ) support_inter++;
2685
2686 if ( support_inter == 0 ) return(-INF);
2687
2688 /* compute the new cut coefficients and rhs (before weakening) */
2689
2690 new_coef = reinterpret_cast<int *> (calloc(inp_ilp->mtcnt[i],sizeof(int)));
2691 if ( new_coef == NULL ) alloc_error(const_cast<char*>("new_coef"));
2692
2693 if ( ( itype == ADD && inp_ilp->msense[i] != 'G' ) ||
2694 ( itype == DEL && inp_ilp->msense[i] == 'G' ) ) flag_gt = 1;
2695 else flag_gt = -1;
2696 if ( flag_gt == 1 ) {
2697 if ( gcdi == 1 ) {
2698 for ( ofsj = 0, ij = begi; ofsj < inp_ilp->mtcnt[i]; ofsj++, ij++ ) {
2699 j = inp_ilp->mtind[ij];
2700 new_coef[ofsj] = cur_cut->non_weak_coef[j] + inp_ilp->mtval[ij];
2701 }
2702 crhs = cur_cut->non_weak_rhs + inp_ilp->mrhs[i];
2703 }
2704 else {
2705 for ( ofsj = 0, ij = begi; ofsj < inp_ilp->mtcnt[i]; ofsj++, ij++ ) {
2706 j = inp_ilp->mtind[ij];
2707 new_coef[ofsj] = cur_cut->non_weak_coef[j] + inp_ilp->mtval[ij] / gcdi;
2708 }
2709 crhs = cur_cut->non_weak_rhs + inp_ilp->mrhs[i] / gcdi;
2710 }
2711 }
2712 else {
2713 if ( gcdi == 1 ) {
2714 for ( ofsj = 0, ij= begi; ofsj < inp_ilp->mtcnt[i]; ofsj++, ij++ ) {
2715 j = inp_ilp->mtind[ij];
2716 new_coef[ofsj] = cur_cut->non_weak_coef[j] - inp_ilp->mtval[ij];
2717 }
2718 crhs = cur_cut->non_weak_rhs - inp_ilp->mrhs[i];
2719 }
2720 else {
2721 for ( ofsj = 0, ij = begi; ofsj < inp_ilp->mtcnt[i]; ofsj++, ij++ ) {
2722 j = inp_ilp->mtind[ij];
2723 new_coef[ofsj] = cur_cut->non_weak_coef[j] - inp_ilp->mtval[ij] / gcdi;
2724 }
2725 crhs = cur_cut->non_weak_rhs - inp_ilp->mrhs[i] / gcdi;
2726 }
2727 }
2728
2729 /* other - relatively fast - check by optimistic evaluation of the
2730 cut score */
2731
2732 weak_loss = cur_cut->min_weak_loss;
2733 one_norm = cur_cut->one_norm;
2734 for ( ofsj = 0, ij = begi; ofsj < inp_ilp->mtcnt[i]; ofsj++, ij++ ) {
2735 j = inp_ilp->mtind[ij];
2736 if ( cur_cut->coef[j] > 0 ) one_norm -= cur_cut->coef[j];
2737 else one_norm += cur_cut->coef[j];
2738 if ( new_coef[ofsj] >= 2 ) one_norm += new_coef[ofsj] / 2;
2739 else one_norm -= new_coef[ofsj] / 2;
2740 if ( mod2(cur_cut->non_weak_coef[j]) == ODD ) {
2741 if ( mod2(new_coef[ofsj]) == EVEN )
2742 weak_loss -= p_ilp->min_loss_by_weak[j];
2743 }
2744 else {
2745 if ( mod2(new_coef[ofsj]) == ODD )
2746 weak_loss += p_ilp->min_loss_by_weak[j];
2747 }
2748 }
2749
2750 viol = ( 1.0 - slack_sum - weak_loss ) / 2.0;
2751
2752 score = tabu_score(NULL,0,viol,static_cast<double> (one_norm));
2753 /*
2754 printf("Score estimate 2 %f Threshold %f\n",score,thresh);
2755 */
2756 if ( score < thresh + ZERO || FAST_SCORE_EVAL ) {
2757
2758 free(new_coef);
2759 return (score);
2760 }
2761 /* get the actual cut coefficients and the violation of the
2762 best cut obtainable trough weakening */
2763
2764 ccoef = reinterpret_cast<int *> (calloc(n,sizeof(int)));
2765 if ( ccoef == NULL ) alloc_error(const_cast<char*>("ccoef"));
2766 for ( j = 0; j < n; j++ ) ccoef[j] = cur_cut->non_weak_coef[j];
2767 for ( ofsj = 0; ofsj < inp_ilp->mtcnt[i]; ofsj++ ) {
2768 ij = begi + ofsj;
2769 j = inp_ilp->mtind[ij];
2770 ccoef[j] = new_coef[ofsj];
2771 }
2772
2773 ok = best_cut(ccoef,&crhs,&viol,FALSE,FALSE);
2774 if ( ok ) score = tabu_score(ccoef,crhs,viol,static_cast<double> (one_norm));
2775 else {
2776 viol = ( - slack_sum - weak_loss ) / 2.0;
2777 score = PENALTY_NON_WEAKABLE + tabu_score(ccoef,crhs,viol,static_cast<double> (one_norm));
2778 }
2779
2780 /*
2781 printf("Score actual %f Threshold %f\n",score,thresh);
2782 */
2783
2784 free(new_coef);
2785 free(ccoef);
2786 return(score);
2787 }
2788
2789 /* modify_current: update the current local search solution by inserting/
2790 deleting a constraint */
2791
modify_current(int i,short int itype)2792 void Cgl012Cut::modify_current(
2793 int i, /* constraint to be moved */
2794 short int itype /* type of move - ADD or DEL */
2795 )
2796 {
2797 int j, begi, gcdi, ofsj, ij;
2798 short int flag_gt;
2799
2800 if ( itype == ADD ) {
2801 cur_cut->n_of_constr++;
2802 cur_cut->in_constr_list[i] = IN;
2803 }
2804 else {
2805 cur_cut->n_of_constr--;
2806 cur_cut->in_constr_list[i] = OUT;
2807 }
2808 last_moved[i] = it;
2809
2810 /* compute the new cut coefficients and rhs (before weakening) */
2811
2812 if ( ( itype == ADD && inp_ilp->msense[i] != 'G' ) ||
2813 ( itype == DEL && inp_ilp->msense[i] == 'G' ) ) flag_gt = 1;
2814 else flag_gt = -1;
2815 begi = inp_ilp->mtbeg[i];
2816 gcdi = p_ilp->gcd[i];
2817 for ( ofsj = 0; ofsj < inp_ilp->mtcnt[i]; ofsj++ ) {
2818 ij = begi + ofsj;
2819 j = inp_ilp->mtind[ij];
2820 /* the '*' and '/' operations can be saved by writing some more code ... */
2821 cur_cut->non_weak_coef[j] += flag_gt * (inp_ilp->mtval[ij] / gcdi);
2822 }
2823 cur_cut->non_weak_rhs += flag_gt * (inp_ilp->mrhs[i] / gcdi);
2824
2825 if ( itype == ADD )
2826 cur_cut->slack_sum += p_ilp->slack[i] / static_cast<double> (gcdi);
2827 else
2828 cur_cut->slack_sum -= p_ilp->slack[i] / static_cast<double> (gcdi);
2829
2830 /* get the best possible cut */
2831
2832 cur_cut->min_weak_loss = 0.0;
2833 for ( j = 0; j < n; j++ ) {
2834 cur_cut->coef[j] = cur_cut->non_weak_coef[j];
2835 if ( mod2(cur_cut->coef[j]) == ODD )
2836 cur_cut->min_weak_loss += p_ilp->min_loss_by_weak[j];
2837 }
2838 cur_cut->rhs = cur_cut->non_weak_rhs;
2839 cur_cut->ok =
2840 best_cut(cur_cut->coef,&cur_cut->rhs,&cur_cut->violation,TRUE,FALSE);
2841 cur_cut->one_norm = 0;
2842 for ( j = 0; j < n; j++ ) {
2843 if ( cur_cut->coef[j] > 0 ) cur_cut->one_norm += cur_cut->coef[j];
2844 else cur_cut->one_norm -= cur_cut->coef[j];
2845 }
2846 }
2847
2848 /* get_current_cut: return a cut data type with the information about
2849 the current cut of the search procedure */
2850
get_current_cut()2851 cut *Cgl012Cut::get_current_cut()
2852 {
2853 int i, j, nz;
2854 /*double viol;*/
2855 cut *cut_ptr;
2856
2857 cut_ptr = reinterpret_cast<cut *> (calloc(1,sizeof(cut)));
2858 if ( cut_ptr == NULL ) alloc_error(const_cast<char*>("cut_ptr"));
2859 cut_ptr->crhs = cur_cut->rhs;
2860 cut_ptr->csense = 'L';
2861 /* count the number of nonzeroes in the cut */
2862 for ( j = 0, nz = 0; j < n; j++ ) if ( cur_cut->coef[j] != 0 ) nz++;
2863 cut_ptr->cnzcnt = nz;
2864 cut_ptr->cind = reinterpret_cast<int *> (calloc(nz,sizeof(int)));
2865 if ( cut_ptr->cind == NULL ) alloc_error(const_cast<char*>("cut_ptr->cind"));
2866 cut_ptr->cval = reinterpret_cast<int *> (calloc(nz,sizeof(int)));
2867 if ( cut_ptr->cval == NULL ) alloc_error(const_cast<char*>("cut_ptr->cval"));
2868 nz = 0; /*viol = 0.0;*/
2869 for ( j = 0; j < n; j++ ) {
2870 if ( cur_cut->coef[j] != 0 ) {
2871 cut_ptr->cind[nz] = j;
2872 cut_ptr->cval[nz] = cur_cut->coef[j];
2873 nz++;
2874 /* viol += p_ilp->xstar[j] * (double) cur_cut->coef[j]; */
2875 }
2876 }
2877 /* viol -= (double) cur_cut->rhs; */
2878 /* cut_ptr->violation = viol; */
2879 cut_ptr->violation = cur_cut->violation;
2880 cut_ptr->n_of_constr = 0;
2881 cut_ptr->constr_list = reinterpret_cast<int *> (calloc(inp_ilp->mr,sizeof(int)));
2882 if ( cut_ptr->constr_list == NULL ) alloc_error(const_cast<char*>("cut_ptr->constr_list"));
2883 cut_ptr->in_constr_list = reinterpret_cast<short int *> (calloc(inp_ilp->mr,sizeof(short int)));
2884 if ( cut_ptr->in_constr_list == NULL ) alloc_error(const_cast<char*>("cut_ptr->in_constr_list"));
2885 for ( i = 0; i < m; i++ ) {
2886 if ( cur_cut->in_constr_list[i] == IN ) {
2887 cut_ptr->in_constr_list[i] = IN;
2888 cut_ptr->constr_list[cut_ptr->n_of_constr] = i;
2889 (cut_ptr->n_of_constr)++;
2890 }
2891 else cut_ptr->in_constr_list[i] = OUT;
2892 }
2893 return(cut_ptr);
2894 }
2895
2896 /* best neighbour: find the cut to be added/deleted from the current
2897 solution among those allowed by the tabu rules */
2898
best_neighbour(cut_list * out_cuts)2899 short int Cgl012Cut::best_neighbour(cut_list *out_cuts /* list of the violated cuts found */)
2900 {
2901 int i, ibest;
2902 short int itype, itypebest=-1;
2903 double score, max_score;
2904 cut *new_cut;
2905
2906 /* cycle through all the constraints in your problem ... */
2907
2908 max_score = -INF;
2909 ibest = NONE;
2910 for ( i = 0; i < m; i++ ) {
2911 if ( ! p_ilp->row_to_delete[i] && allowed(i) ) {
2912 if ( in_cur_cut(i) ) {
2913 /* constraint i is in the current set of constraints
2914 (those defining the current cut) */
2915 itype = DEL;
2916 }
2917 else {
2918 /* constraint i is not in the current set of constraints
2919 (those defining the current cut) */
2920 itype = ADD;
2921 }
2922 score = score_by_moving(i,itype,max_score);
2923 if ( score > max_score ) {
2924 /* best cut found in this iteration: store it */
2925 ibest = i;
2926 itypebest = itype;
2927 max_score = score;
2928 }
2929 }
2930 } /* for ( i = 0; i < m; i++ ) */
2931
2932 if ( ibest == NONE ) {
2933 #ifdef PRINT_TABU
2934 printf("No move could be performed by best_neighbour\n");
2935 #endif
2936 return(TRUE);
2937 }
2938 modify_current(ibest,itypebest);
2939 if ( cur_cut->violation > MIN_VIOLATION + EPS ) {
2940 #ifdef PRINT_TABU
2941 printf("... adding the current cut to the output list - it = %d viol = %f\n",it, cur_cut->violation);
2942 #endif
2943 new_cut = get_current_cut();
2944 out_cuts = add_cut_to_list(new_cut,out_cuts);
2945 last_it_add = it;
2946 }
2947 return(FALSE);
2948 }
2949
2950 /* memory_reaction: perform the long term reaction by cheching whether the
2951 current solution has already been visited or the best solution has not
2952 been updated for too many iterations */
2953
memory_reaction()2954 void memory_reaction()
2955 {
2956 int cycle_length;
2957
2958 if ( hash_search(&cycle_length) ) {
2959 if ( cycle_length < 2 * ( m - 1 ) ) {
2960 increase_prohib_period();
2961 return;
2962 }
2963 }
2964 else hash_insert();
2965 if ( it - last_prohib_period_mod > B )
2966 decrease_prohib_period();
2967 }
2968
2969 /* add_tight_constraint: initialize the current cut by adding a tight
2970 constraint to it */
2971
add_tight_constraint()2972 void Cgl012Cut::add_tight_constraint()
2973 {
2974 int i, ntight;
2975 double smin=COIN_DBL_MAX;
2976 abort();
2977 int *tight;
2978
2979 ntight = 0;
2980 tight = reinterpret_cast<int *> (calloc(m,sizeof(int)));
2981 if ( tight == NULL ) alloc_error(const_cast<char*>("tight"));
2982 for ( i = 0; i < m; i++ ) {
2983 /* search for the tightest constraint never added to cut */
2984 if ( last_moved[i] < 0 && p_ilp->slack[i] < smin ) {
2985 if ( p_ilp->slack[i] < ZERO ) {
2986 /* tight constraint */
2987 smin = ZERO;
2988 tight[ntight] = i;
2989 ntight++;
2990 }
2991 else {
2992 /* best constraint so far */
2993 smin = p_ilp->slack[i];
2994 tight[0] = i;
2995 ntight = 1;
2996 }
2997 }
2998 }
2999 if ( ntight > 0 ) i = tight[rand() % ntight];
3000 /* if all constraints have already been in cur_cut choose first at random */
3001 else i = rand() % m;
3002 free(tight);
3003 modify_current(i,ADD);
3004 }
3005
3006 /* initialize: initialize the data structures for local search */
3007
initialize()3008 void Cgl012Cut::initialize()
3009 {
3010 int i;
3011
3012 m = inp_ilp->mr;
3013 n = inp_ilp->mc;
3014 it = 0;
3015 last_it_add = 0;
3016 last_it_restart = 0;
3017 last_prohib_period_mod = 0;
3018 prohib_period = IN_PROHIB_PERIOD;
3019 initialize_cur_cut();
3020 last_moved = reinterpret_cast<int *> (calloc(m,sizeof(int)));
3021 if ( last_moved == NULL ) alloc_error(const_cast<char*>("last_moved"));
3022 for ( i = 0; i < m; i++ ) {
3023 last_moved[i] = -COIN_INT_MAX;
3024 }
3025 initialize_hash_table();
3026 add_tight_constraint();
3027 A = m;
3028 B = 10 * m;
3029 }
3030
3031 /* restart: perform a restart of the search - IMPORTANT: in the current
3032 implementation vector last_moved is not cleared at restart */
3033
restart(short int failure)3034 void Cgl012Cut::restart(short int failure /* flag forcing the restart if some trouble occurred */)
3035 {
3036 if ( failure || ( it - last_it_add > A && it - last_it_restart > A ) ) {
3037 /* perform restart */
3038 last_it_restart = it;
3039 prohib_period = IN_PROHIB_PERIOD;
3040 last_prohib_period_mod = it;
3041 clear_hash_table();
3042 clear_cur_cut();
3043 add_tight_constraint();
3044 }
3045 }
3046
3047 /* free_memory: free the memory used by local search */
3048
free_memory()3049 void free_memory()
3050 {
3051 free_cur_cut();
3052 free(last_moved);
3053 free_hash_table();
3054 }
3055
3056 /* tabu_012: try to identify violated 0-1/2 cuts by a simple tabu search
3057 procedure adapted from that used by Battiti and Protasi for finding
3058 large cliques */
3059
tabu_012()3060 cut_list *Cgl012Cut::tabu_012()
3061 {
3062 short int failure;
3063 cut_list *out_cuts;
3064
3065 out_cuts = initialize_cut_list(MAX_CUTS);
3066 initialize();
3067
3068 it = 0;
3069 do {
3070 memory_reaction();
3071 failure = best_neighbour(out_cuts);
3072 #ifdef PRINT_TABU
3073 print_cur_cut();
3074 #endif
3075 it++;
3076 restart(failure);
3077
3078 }
3079 while ( out_cuts->cnum < MAX_CUTS && it < MAX_TABU_ITER );
3080 free_memory();
3081 #ifdef PRINT_TABU
3082 printf("Number of violated cuts found by Tabu Search %d\n",out_cuts->cnum);
3083 printf("Tabu Search timings: best_neighbour %f score_by_moving %f coefficient %f best_cut %f\n",
3084 time_best_neigh, time_scor_by_mov, time_coef, time_best_cut);
3085 #endif
3086 return(out_cuts);
3087 }
3088 //end include "Cgltabu_012.c"
3089 #ifdef POOL
3090
3091 /* same_pool_cut: check whether two pool cuts are in fact the same cut */
3092
same_pool_cut(pool_cut * p_cut1,pool_cut * p_cut2)3093 short int same_pool_cut(pool_cut *p_cut1, pool_cut *p_cut2)
3094 {
3095 int c;
3096
3097 if ( p_cut1->n_of_constr != p_cut2->n_of_constr ) return(FALSE);
3098 if ( p_cut1->code != p_cut2->code ) return(FALSE);
3099 /* assumes constr_list is sorted for increasing/decreasing constraint index */
3100 for ( c = 1; c < p_cut1->n_of_constr; c++ )
3101 if ( p_cut1->constr_list[c] != p_cut2->constr_list[c] ) return(FALSE);
3102 return(TRUE);
3103 }
3104
3105 /* free_pool_cut: free the memory for a non-empty pool cut */
3106
free_pool_cut(pool_cut * p_cut)3107 void free_pool_cut(pool_cut *p_cut)
3108 {
3109 if ( p_cut == NULL ) return;
3110 if ( p_cut->constr_list != NULL ) free(p_cut->constr_list);
3111 free(p_cut);
3112 }
3113
3114 /* initialize_pool: initialize the pool data structure */
3115
initialize_pool()3116 void initialize_pool()
3117 {
3118 pool = (pool_cut_list *) calloc(1,sizeof(pool_cut_list));
3119 if ( pool == NULL ) alloc_error(const_cast<char*>("pool"));
3120 pool->cnum = 0;
3121 pool->list = (pool_cut **) calloc(MAX_CUT_POOL,sizeof(pool_cut *));
3122 if ( pool->list == NULL ) alloc_error(const_cast<char*>("pool->list"));
3123 pool->ncod = (int *) calloc(MAX_CUT_COD,sizeof(int));
3124 if ( pool->ncod == NULL ) alloc_error(const_cast<char*>("pool->ncod"));
3125 }
3126
3127 /* free_pool: free the memory used by the pool */
3128
free_pool()3129 void free_pool()
3130 {
3131 int c;
3132
3133 if ( pool == NULL ) return;
3134 for ( c = 0; c < pool->cnum; c++ )
3135 free_pool_cut(pool->list[c]);
3136 free(pool);
3137 }
3138
3139 /* clean_pool: remove form the pool the cuts which are inactive since a
3140 large number of iterations */
3141
clean_pool()3142 void clean_pool()
3143 {
3144 int c, d;
3145
3146 /* the pool compression could be implemented more efficiently */
3147
3148 for ( c = 0; c < pool->cnum; c++ ) {
3149 if ( pool->list[c]->it_found > MAX_ITER_POOL ) {
3150 free_pool_cut(pool->list[c]);
3151 pool->list[c] = NULL;
3152 for ( d = c ; d < pool->cnum; d++ )
3153 pool->list[d] = pool->list[d + 1];
3154 (pool->cnum)--;
3155 }
3156 }
3157 }
3158
3159 /* insert_cut_in_pool: add a cut to the pool if there is space */
3160
insert_cut_in_pool(pool_cut * p_cut)3161 void insert_cut_in_pool(pool_cut *p_cut)
3162 {
3163
3164 if ( pool->cnum == MAX_CUT_POOL ) {
3165 #ifdef COIN_DEVELOP
3166 printf("Warning: pool is full and separated cuts cannot be added\n");
3167 #endif
3168 return;
3169 }
3170 pool->list[pool->cnum] = p_cut;
3171 (pool->cnum)++;
3172 (pool->ncod[p_cut->code])++;
3173 }
3174
3175 /* cut_is_in_pool: check whether a given cut is already in the pool */
3176
cut_is_in_pool(cut * v_cut)3177 short int cut_is_in_pool(cut *v_cut)
3178 {
3179 int c, i, cod;
3180 short int equal;
3181 short int *flag_v;
3182 pool_cut *p_cut;
3183
3184 /*
3185 print_cut(v_cut);
3186 printf("checking for a cut in the pool ...\n");
3187 */
3188 cod = hash_addr(inp_ilp->mr,v_cut->in_constr_list);
3189 if ( pool->ncod[cod] == 0 ) return(FALSE);
3190 /* trivial sequential search */
3191 /*
3192 printf("... sequential search needed ...\n");
3193 */
3194 flag_v = (short int *) calloc(inp_ilp->mr,sizeof(short int));
3195 if ( flag_v == NULL ) alloc_error(const_cast<char*>("flag_v"));
3196 for ( c = 0; c < pool->cnum; c++ ) {
3197 if ( pool->list[c]->code != cod ) continue;
3198 p_cut = pool->list[c];
3199 equal = TRUE;
3200 for ( i = 0; i < inp_ilp->mr; i++ ) flag_v[i] = OUT;
3201 for ( i = 0; i < p_cut->n_of_constr; i++ )
3202 flag_v[p_cut->constr_list[i]] = IN;
3203 for ( i = 0; i < inp_ilp->mr; i++ ) {
3204 if ( v_cut->in_constr_list[i] != flag_v[i] ) {
3205 equal = FALSE;
3206 break;
3207 }
3208 }
3209 if ( equal ) {
3210 free(flag_v);
3211 /*
3212 printf("... cut was in the pool!\n");
3213 */
3214 return(TRUE);
3215 }
3216 }
3217 return(FALSE);
3218 }
3219
3220 /* good_pool_cut: check whether the current cut is worth adding to the pool */
3221
good_pool_cut(cut * v_cut)3222 short int good_pool_cut(cut *v_cut)
3223 {
3224 /* no check performed */
3225 return(TRUE);
3226 }
3227
3228 /* add_cuts_to_pool: add the cuts separated to the pool structure */
3229
add_cuts_to_pool(cut_list * out_cuts)3230 void add_cuts_to_pool(cut_list *out_cuts)
3231 {
3232 int i, c;
3233 cut *v_cut;
3234 pool_cut *p_cut;
3235 /* float ti, tf, tti, ttf; */
3236 /* static float tcutis = 0.0, taddcut = 0.0; */
3237
3238 /* second_(&tti); */
3239
3240 if ( pool == NULL ) initialize_pool();
3241 if ( pool->cnum >= MAX_CUT_POOL * CLEAN_THRESH ) clean_pool();
3242
3243 for ( i = 0; i < out_cuts->cnum; i++ ) {
3244 v_cut = out_cuts->list[i];
3245 /* second_(&ti); */
3246 if ( cut_is_in_pool(v_cut) ) {
3247 /* second_(&tf); */
3248 /* tcutis += tf - ti; */
3249 continue;
3250 }
3251 else {
3252 /* second_(&tf); */
3253 /* tcutis += tf - ti; */
3254 }
3255 if ( good_pool_cut(v_cut) ) {
3256 /* add the cut to the pool list */
3257 p_cut = (pool_cut *) calloc(1,sizeof(pool_cut));
3258 if ( p_cut == NULL ) alloc_error(const_cast<char*>("p_cut"));
3259 p_cut->n_of_constr = v_cut->n_of_constr;
3260 p_cut->constr_list = (int *) calloc(v_cut->n_of_constr,sizeof(int));
3261 if ( p_cut->constr_list == NULL ) alloc_error(const_cast<char*>("p_cut->constr_list"));
3262 for ( c = 0; c < v_cut->n_of_constr; c++ )
3263 p_cut->constr_list[c] = v_cut->constr_list[c];
3264 p_cut->code = hash_addr(inp_ilp->mr,v_cut->in_constr_list);
3265 p_cut->n_it_violated = 0;
3266 p_cut->it_found = sep_iter;
3267 insert_cut_in_pool(p_cut);
3268 }
3269 }
3270
3271 /* second_(&ttf); */
3272 /* taddcut += ttf - tti; */
3273 /* printf("add_cuts_to_pool: tcutis %f taddcut %f\n",tcutis,taddcut); */
3274
3275 }
3276
3277 /* interesting_var: decides whether a variable is relevant in the
3278 separation or not */
3279
interesting_var(int j)3280 short int interesting_var(int j /* variable to be evaluated */)
3281 {
3282 /* return ( vlog[j]->n_it_zero < MANY_IT_ZERO ); */
3283 /* if ( aggr ) return (TRUE); */
3284 return( ! p_ilp->col_to_delete[j] );
3285 }
3286
3287 /* get_cuts_from_pool: select from the pool a convenient set of violated
3288 constraints to be added to the current LP */
3289
3290 static double max_score_ever = ZERO; /* maximum score of a violated cut during
3291 the whole cutting plane procedure */
3292 #define MIN_CUT_SCORE ( max_score_ever / MIN_SCORE_RANGE )
3293 #define MAX_CUT_SCORE ( max_score_ever / MAX_SCORE_RANGE )
3294 #define MIN_IT_VIOL 2
3295
get_cuts_from_pool(short int after_sep)3296 cut_list *get_cuts_from_pool(
3297 short int after_sep /* flag telling whether the pool is searched after
3298 a new separation in which case only new cuts are
3299 checked */
3300 )
3301 {
3302 int c, crhs, j, k, l, maxc, n_interest_var;
3303 double viol, score, maxscore, min_cut_score, max_cut_score;
3304 short int ok;
3305 int *interest_var, *ccoef;
3306 short int *added;
3307 pool_cut *p_cut;
3308 select_cut **best_var_cut;
3309 cut *a_cut;
3310 cut_list *add_cuts;
3311 /* float ti, tf, tti, ttf; */
3312 /* static float tgetcut = 0.0, talloc = 0.0, tgetcoef = 0.0, tbestcut = 0.0, tscore = 0.0, tupdbest = 0.0, tadd = 0.0; */
3313
3314 /* second_(&tti); */
3315
3316 if ( pool == NULL ) {
3317 add_cuts = initialize_cut_list(1);
3318 return(add_cuts);
3319 }
3320
3321 /* in the current implementation, the cut with best score with nonzero
3322 coefficient for each "interesting" variable is added to the current LP,
3323 provided the cut satisfies some requirements (violation, depth, etc.) */
3324
3325 /* define the set of the interesting variables */
3326 interest_var = (int *) calloc(p_ilp->mc,sizeof(int));
3327 if ( interest_var == NULL ) alloc_error(const_cast<char*>("interest_var"));
3328 n_interest_var = 0;
3329 for ( j = 0; j < p_ilp->mc; j++ ) {
3330 if ( interesting_var(j) ) {
3331 interest_var[n_interest_var] = j;
3332 n_interest_var++;
3333 }
3334 }
3335 best_var_cut = (select_cut **) calloc(n_interest_var,sizeof(select_cut *));
3336 if ( best_var_cut == NULL ) alloc_error(const_cast<char*>("best_var_cut"));
3337 for ( k = 0; k < n_interest_var; k++ ) {
3338 best_var_cut[k] = (select_cut *) calloc(1,sizeof(select_cut));
3339 if ( best_var_cut[k] == NULL ) alloc_error(const_cast<char*>("best_var_cut[k]"));
3340 best_var_cut[k]->ccoef = (int *) calloc(p_ilp->mc,sizeof(int));
3341 if ( best_var_cut[k]->ccoef == NULL )
3342 alloc_error(const_cast<char*>("best_var_cut[k]->ccoef"));
3343 best_var_cut[k]->score = -INF;
3344 }
3345
3346 /* find the cuts with the best scores and the list of the cuts
3347 violated by the current fractional point */
3348 maxscore = -INF;
3349 ccoef = (int *) calloc(p_ilp->mc,sizeof(int));
3350 if ( ccoef == NULL ) alloc_error(const_cast<char*>("ccoef"));
3351
3352 /* second_(&tf); */
3353 /* talloc += tf - tti; */
3354
3355 min_cut_score = MIN_CUT_SCORE;
3356 max_cut_score = MAX_CUT_SCORE;
3357
3358 for ( c = 0; c < pool->cnum; c++ ) {
3359 p_cut = pool->list[c];
3360 /* if a new separation was made before the call check only new cuts */
3361 if ( after_sep && p_cut->it_found != sep_iter ) continue;
3362 /* determine the actual coefficients of the cut */
3363 /* second_(&ti); */
3364 ok = get_ori_cut_coef(p_cut->n_of_constr,p_cut->constr_list,
3365 ccoef,&crhs,TRUE);
3366 /* second_(&tf); */
3367 /* tgetcoef += tf - ti; */
3368 /* second_(&ti); */
3369 ok = ok && best_cut(ccoef,&crhs,&viol,TRUE,TRUE);
3370 /* second_(&tf); */
3371 /* tbestcut += tf - ti; */
3372 if ( ok && viol > MIN_VIOLATION ) {
3373 (p_cut->n_it_violated)++;
3374 /* second_(&ti); */
3375 score = cut_score(ccoef,crhs,viol,TRUE);
3376 if ( score > maxscore ) {
3377 maxscore = score;
3378 maxc = c;
3379 }
3380 /* second_(&tf); */
3381 /* tscore += tf - ti; */
3382 if ( ! aggr ) {
3383 if ( score < min_cut_score ) continue;
3384 if ( score < max_cut_score && p_cut->n_it_violated < MIN_IT_VIOL )
3385 continue;
3386 }
3387 /* second_(&ti); */
3388 for ( k = 0; k < n_interest_var; k++ ) {
3389 j = interest_var[k];
3390 if ( ccoef[j] != 0 && best_var_cut[k]->score < score ) {
3391 best_var_cut[k]->score = score;
3392 for ( l = 0; l < p_ilp->mc; l++ )
3393 best_var_cut[k]->ccoef[l] = ccoef[l];
3394 best_var_cut[k]->crhs = crhs;
3395 best_var_cut[k]->pool_index = c;
3396 }
3397 }
3398 /* second_(&tf); */
3399 /* tupdbest += tf - ti; */
3400 }
3401 else p_cut->n_it_violated = 0;
3402 }
3403 free(ccoef);
3404
3405 /* printf("maxscore of a cut : %f ever: %f\n",maxscore,max_score_ever); */
3406 if ( maxscore > max_score_ever ) max_score_ever = maxscore;
3407
3408 /* second_(&ti); */
3409
3410 add_cuts = initialize_cut_list(n_interest_var);
3411 added = (short int *) calloc(pool->cnum,sizeof(short int));
3412 if ( added == NULL ) alloc_error(const_cast<char*>("added"));
3413 for ( c = 0; c < pool->cnum; c++ ) added[c] = FALSE;
3414
3415 for ( k = 0; k < n_interest_var; k++ ) {
3416 j = interest_var[k];
3417 if ( ! added[best_var_cut[k]->pool_index] &&
3418 best_var_cut[k]->score >= MIN_CUT_SCORE ) {
3419 /* add the cut to the cut list on output */
3420 a_cut = define_cut(best_var_cut[k]->ccoef,best_var_cut[k]->crhs);
3421 add_cuts = add_cut_to_list(a_cut,add_cuts);
3422 added[best_var_cut[k]->pool_index] = TRUE;
3423 }
3424 free(best_var_cut[k]->ccoef);
3425 free(best_var_cut[k]);
3426 }
3427 free(best_var_cut);
3428 free(interest_var);
3429 free(added);
3430
3431 /* second_(&tf); */
3432 /* tadd += tf - ti; */
3433
3434 /* second_(&ttf); */
3435 /* tgetcut += ttf - tti; */
3436 /* 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); */
3437
3438 return(add_cuts);
3439 }
3440 #endif
3441 /* initialize_log_var: initialize the log information for the problem variables */
3442
initialize_log_var()3443 void Cgl012Cut::initialize_log_var()
3444 {
3445 int j;
3446 if (!vlog) {
3447 if (p_ilp->mc) {
3448 vlog = reinterpret_cast<log_var **> (calloc(p_ilp->mc,sizeof(log_var *)));
3449 if ( vlog == NULL ) alloc_error(const_cast<char*>("vlog"));
3450 for ( j = 0; j < p_ilp->mc; j++ ) {
3451 vlog[j] = reinterpret_cast<log_var *> (calloc(1,sizeof(log_var)));
3452 if ( vlog[j] == NULL ) alloc_error(const_cast<char*>("vlog[j]"));
3453 vlog[j]->n_it_zero = 0;
3454 }
3455 }
3456 } else {
3457 // just initialize counts
3458 for ( j = 0; j < p_ilp->mc; j++ ) {
3459 vlog[j]->n_it_zero = 0;
3460 }
3461 }
3462 }
3463
3464 /* free_log_var */
3465
free_log_var()3466 void Cgl012Cut::free_log_var()
3467 {
3468 if (vlog) {
3469 int j;
3470 for ( j = 0; j < p_ilp->mc; j++ ) free(vlog[j]);
3471 free(vlog);
3472 vlog=NULL;
3473 }
3474 }
3475
3476 /* update_log_var: update the log information for the problem variables */
3477
update_log_var()3478 void Cgl012Cut::update_log_var()
3479 {
3480 int j;
3481
3482 /* so far one counts only the number of consecutive iterations with
3483 0 value for each variable */
3484
3485 if ( vlog == NULL ) initialize_log_var();
3486
3487 for ( j = 0; j < p_ilp->mc; j++ ) {
3488 if ( p_ilp->xstar[j] < ZERO && p_ilp->xstar[j] > - ZERO )
3489 vlog[j]->n_it_zero++;
3490 else vlog[j]->n_it_zero = 0;
3491 }
3492 }
3493
3494 /* the final implementation should use the following additional functions:
3495 init_sep_012_cut: defines the inp_ilp and p_ilp data structures and
3496 initializes pool and vlog
3497 sep_012_cut: updates inp_ilp and p_ilp, performs separation and pool
3498 management
3499 kill_sep_012_cut: frees all the permanent memory (inp_ilp, p_ilp,
3500 pool, vlog, etc.)
3501 */
3502
sep_012_cut(int mr,int mc,int mnz,int * mtbeg,int * mtcnt,int * mtind,int * mtval,int * vlb,int * vub,int * mrhs,char * msense,const double * xstar,bool aggressive,int * cnum,int * cnzcnt,int ** cbeg,int ** ccnt,int ** cind,int ** cval,int ** crhs,char ** csense)3503 int Cgl012Cut::sep_012_cut(
3504 /*
3505 INPUT parameters:
3506 */
3507 int mr, /* number of rows in the ILP matrix */
3508 int mc, /* number of columns in the ILP matrix */
3509 int mnz, /* number of nonzero's in the ILP matrix */
3510 int *mtbeg, /* starting position of each row in arrays mtind and mtval */
3511 int *mtcnt, /* number of entries of each row in arrays mtind and mtval */
3512 int *mtind, /* column indices of the nonzero entries of the ILP matrix */
3513 int *mtval, /* values of the nonzero entries of the ILP matrix */
3514 int *vlb, /* lower bounds on the variables */
3515 int *vub, /* upper bounds on the variables */
3516 int *mrhs, /* right hand sides of the constraints */
3517 char *msense, /* senses of the constraints: 'L', 'G' or 'E' */
3518 const double *xstar, /* current optimal solution of the LP relaxation */
3519 bool aggressive, /* flag asking whether as many cuts as possible are
3520 required on output (TRUE) or not (FALSE) */
3521 /*
3522 OUTPUT parameters (the memory for the vectors is allocated INTERNALLY
3523 by the procedure: if some memory is already allocated, it is FREED):
3524 */
3525 int *cnum, /* number of violated 0-1/2 cuts identified by the procedure */
3526 int *cnzcnt, /* overall number of nonzero's in the cuts */
3527 int **cbeg, /* starting position of each cut in arrays cind and cval */
3528 int **ccnt, /* number of entries of each cut in arrays cind and cval */
3529 int **cind, /* column indices of the nonzero entries of the cuts */
3530 int **cval, /* values of the nonzero entries of the cuts */
3531 int **crhs, /* right hand sides of the cuts */
3532 char **csense /* senses of the cuts: 'L', 'G' or 'E' */
3533 /*
3534 NOTE that all the numerical input/output vectors are INTEGER (with
3535 the exception of xstar), since the procedure is intended to work
3536 with pure ILP's, and that the ILP matrix has to be given on input
3537 in ROW format.
3538 */
3539 )
3540 {
3541 #ifdef TIME
3542 float tbasi, tbasf;
3543 #endif
3544 errorNo=0;
3545 cut_list *out_cuts, *add_cuts;
3546 #ifdef TIME
3547 second_(&tti);
3548 #endif
3549
3550 aggr = aggressive;
3551
3552 /* load the input ILP into an internal data structure */
3553
3554 //ilp_load(mr,mc,mnz,mtbeg,mtcnt,mtind,mtval,
3555 // vlb,vub,mrhs,msense);
3556 if (!inp_ilp)
3557 return FALSE;
3558 inp_ilp->xstar = xstar;
3559
3560
3561 /* construct an internal data structure containing all the information
3562 which can be useful for 0-1/2 cut separation - this may in fact be
3563 done only at the first call of the separation procedure */
3564
3565 #ifdef TIME
3566 second_(&ti);
3567 #endif
3568
3569 get_parity_ilp();
3570
3571 /*
3572 print_double_vect("xstar",p_ilp->xstar,p_ilp->mc);
3573 print_parity_ilp();
3574 */
3575
3576 #ifdef TIME
3577 second_(&tf);
3578 prep_time += tf - ti;
3579 #endif
3580
3581 if ( p_ilp->mnz == 0 ) {
3582 #ifdef COIN_DEVELOP
3583 printf("Warning: no significant constraint for 0-1/2 cut separation\n");
3584 printf("... end separation\n");
3585 #endif
3586 //free_ilp();
3587 //free_parity_ilp();
3588 #ifdef TIME
3589 second_(&ttf);
3590 total_time += ttf - tti;
3591 #endif
3592 return(FALSE);
3593 }
3594
3595 /* print_double_vect("xstar",p_ilp->xstar,p_ilp->mc); */
3596
3597 sep_iter++;
3598 update_log_var();
3599
3600 #ifdef POOL
3601
3602 /* search for possible violated cuts in the pool */
3603
3604 #ifdef TIME
3605 second_(&tpi);
3606 #endif
3607 add_cuts = get_cuts_from_pool(FALSE);
3608 #ifdef TIME
3609 second_(&tpf);
3610 pool_time += tpf - tpi;
3611 #ifdef PRINT_TIME
3612 printf("... time elapsed at the end of get_cuts_from_pool: %f\n",tpf - tti);
3613 #endif
3614 #endif
3615
3616 if ( add_cuts->cnum == 0 ) {
3617 free_cut_list(add_cuts);
3618 }
3619 else {
3620 /* printf("Violated cuts found in the pool - no separation procedure used\n"); */
3621 goto free_memory;
3622 }
3623
3624 #endif
3625
3626 /* try to identify violated 0-1/2 cuts by using the original procedure
3627 described in Caprara and Fischetti's MP paper */
3628
3629 #ifdef TIME
3630 second_(&tbasi);
3631 #endif
3632
3633 out_cuts = basic_separation();
3634
3635 #ifdef TIME
3636 second_(&tbasf);
3637 tot_basic_sep_time += tbasf - tbasi;
3638 avg_basic_sep_time = tot_basic_sep_time / (float) sep_iter;
3639 #endif
3640
3641 //#define TABU_SEARCH
3642 #ifdef TABU_SEARCH
3643
3644 /* try to identify violated cuts by tabu search if none was found */
3645
3646 if ( out_cuts->cnum == 0 ) {
3647 free_cut_list(out_cuts);
3648
3649 #ifdef TIME
3650 second_(&ttabi);
3651 #endif
3652
3653 out_cuts = tabu_012();
3654
3655 #ifdef TIME
3656 second_(&ttabf);
3657 tabu_time += ttabf - ttabi;
3658 #endif
3659
3660 }
3661
3662 #endif
3663
3664 #ifdef POOL
3665
3666 /* add the cuts separated to the pool */
3667
3668 #ifdef TIME
3669 second_(&tpi);
3670 #endif
3671 add_cuts_to_pool(out_cuts);
3672 free_cut_list(out_cuts);
3673
3674 /* select from the pool a convenient set of violated constraints
3675 to be added to the current LP */
3676
3677 add_cuts = get_cuts_from_pool(TRUE);
3678 #ifdef TIME
3679 second_(&tpf);
3680 pool_time += tpf - tpi;
3681 #ifdef PRINT_TIME
3682 printf("... time elapsed at the end of get_cuts_from_pool: %f\n",tpf - tti);
3683 #endif
3684 #endif
3685
3686 #else
3687
3688 /* give on output the cuts separated */
3689
3690 add_cuts = out_cuts;
3691
3692 #endif
3693
3694 //free_ilp();
3695 //free_parity_ilp();
3696
3697 #ifdef TIME
3698 second_(&ttf);
3699 total_time += ttf - tti;
3700 #endif
3701
3702 if ( add_cuts->cnum > 0 ) {
3703 getcuts(add_cuts,cnum,cnzcnt,cbeg,ccnt,cind,cval,crhs,csense);
3704 /* print_cut_list(add_cuts); */
3705 free_cut_list(add_cuts);
3706 return(TRUE);
3707 }
3708 else {
3709 free_cut_list(add_cuts);
3710 return(FALSE);
3711 }
3712 }
3713
3714 #ifdef TIME
3715 /* print_times: print the timings of the separation procedure */
3716
print_times()3717 void print_times()
3718 {
3719 printf("... separation timings \n");
3720 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",
3721 total_time, prep_time, weak_time, aux_time, path_time, cycle_time, cut_time, cut_ncalls, bw_time, coef_time, pool_time, tabu_time);
3722 }
3723 #endif
3724 //-------------------------------------------------------------------
3725 // Default Constructor
3726 //-------------------------------------------------------------------
Cgl012Cut()3727 Cgl012Cut::Cgl012Cut () :
3728 inp_ilp(NULL),
3729 p_ilp(NULL),
3730 iter(0),
3731 gap(0.0),
3732 maxgap(0.0),
3733 errorNo(0),
3734 sep_iter(0),
3735 vlog(NULL),
3736 aggr(true)
3737 {
3738 // nothing to do here
3739 }
3740 //-------------------------------------------------------------------
3741 // Copy constructor
3742 //-------------------------------------------------------------------
Cgl012Cut(const Cgl012Cut & rhs)3743 Cgl012Cut::Cgl012Cut (const Cgl012Cut & rhs) :
3744 inp_ilp(NULL),
3745 p_ilp(NULL),
3746 iter(rhs.iter),
3747 gap(rhs.gap),
3748 maxgap(rhs.maxgap),
3749 errorNo(rhs.errorNo),
3750 sep_iter(rhs.sep_iter),
3751 vlog(NULL),
3752 aggr(rhs.aggr)
3753 {
3754 if (rhs.p_ilp||rhs.vlog||inp_ilp)
3755 abort();
3756 }
3757
3758 //-------------------------------------------------------------------
3759 // Destructor
3760 //-------------------------------------------------------------------
~Cgl012Cut()3761 Cgl012Cut::~Cgl012Cut ()
3762 {
3763 free_log_var();
3764 free_parity_ilp();
3765 free_ilp();
3766 }
3767
3768 //----------------------------------------------------------------
3769 // Assignment operator
3770 //-------------------------------------------------------------------
3771 Cgl012Cut &
operator =(const Cgl012Cut & rhs)3772 Cgl012Cut::operator=(
3773 const Cgl012Cut& rhs)
3774 {
3775 if (this != &rhs) {
3776 if (rhs.p_ilp||rhs.vlog||inp_ilp)
3777 abort();
3778 free_log_var();
3779 free_parity_ilp();
3780 free_ilp();
3781 #if 0
3782 inp_ilp = reinterpret_cast<ilp *> (calloc(1,sizeof(ilp)));
3783 if ( inp_ilp == NULL ) alloc_error(const_cast<char*>("inp_ilp"));
3784
3785 inp_ilp->mr = rhs.inp_ilp->mr; inp_ilp->mc = rhs.inp_ilp->mc;
3786 inp_ilp->mnz = rhs.inp_ilp->mnz;
3787 inp_ilp->mtbeg = rhs.inp_ilp->mtbeg; inp_ilp->mtcnt = rhs.inp_ilp->mtcnt;
3788 inp_ilp->mtind = rhs.inp_ilp->mtind; inp_ilp->mtval = rhs.inp_ilp->mtval;
3789 inp_ilp->vlb = rhs.inp_ilp->vlb; inp_ilp->vub = rhs.inp_ilp->vub;
3790 inp_ilp->mrhs = rhs.inp_ilp->mrhs; inp_ilp->msense = rhs.inp_ilp->msense;
3791 #endif
3792 iter = rhs.iter;
3793 gap = rhs.gap;
3794 maxgap = rhs.maxgap;
3795 errorNo = rhs.errorNo;
3796 sep_iter = rhs.sep_iter;
3797 aggr = rhs.aggr;
3798 }
3799 return *this;
3800 }
3801