1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2 /* */
3 /* This file is part of the program and library */
4 /* SCIP --- Solving Constraint Integer Programs */
5 /* */
6 /* Copyright (C) 2002-2021 Konrad-Zuse-Zentrum */
7 /* fuer Informationstechnik Berlin */
8 /* */
9 /* SCIP is distributed under the terms of the ZIB Academic License. */
10 /* */
11 /* You should have received a copy of the ZIB Academic License */
12 /* along with SCIP; see the file COPYING. If not visit scipopt.org. */
13 /* */
14 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
15 /**@file sepa_subtour.c
16 * @brief If there exists a transition forward along the cycle, then the state that the transition originates from can
17 * be reached only after another ncluster - 1 transitions. Therefore cycles with a number of transitions smaller than
18 * that can be separated.
19 * @author Leon Eifler
20 */
21
22 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
23
24 #include <assert.h>
25 #include <string.h>
26
27 #include "sepa_subtour.h"
28
29 #include "probdata_cyc.h"
30 #include "scip/cons_linear.h"
31 #include "scip/pub_misc.h"
32
33 #define SEPA_NAME "subtour"
34 #define SEPA_DESC "separator that elininates subtours of length smaller than |NCluster|"
35 #define SEPA_PRIORITY 1000
36 #define SEPA_FREQ 5
37 #define SEPA_MAXBOUNDDIST 0.0
38 #define SEPA_USESSUBSCIP FALSE /**< does the separator use a secondary SCIP instance? */
39 #define SEPA_DELAY FALSE /**< should separation method be delayed, if other separators found cuts? */
40 #define MAXCUTS 2000
41 #define MAXROUNDS 15
42
43 #ifdef SCIP_DEBUG
44 /** Print a cycle to the command line. For debugging purposes */
45 static
printCycle(SCIP * scip,int * cycle,int cyclelength,int nstates)46 void printCycle(
47 SCIP* scip, /**< SCIP data structure */
48 int* cycle, /**< The cycle to be printed */
49 int cyclelength, /**< The length of the cycle */
50 int nstates /**< The number of states */
51 )
52 {
53 int i;
54
55 SCIPinfoMessage(scip, NULL, "cycle_l%d_c: %d", cyclelength, cycle[0]);
56 for( i = 0; i < cyclelength; ++i )
57 {
58 SCIPinfoMessage(scip, NULL, " -> %d", cycle[i+1]);
59 }
60 SCIPinfoMessage(scip, NULL, "\n");
61 }
62 #endif
63
64 /** get distance of longest path between two states with exactly n arcs from the matrix */
65 static
getDist(SCIP_Real *** adjacencymatrix,int n,int state1,int state2)66 SCIP_Real getDist(
67 SCIP_Real*** adjacencymatrix, /**< the adjacency-matrices of all paths with 1,...,|Clutster| arcs */
68 int n, /**< length */
69 int state1, /**< starting state */
70 int state2 /**< end state */
71 )
72 {
73 assert(adjacencymatrix[n] != NULL);
74 assert(adjacencymatrix[n][state1] != NULL);
75
76 return adjacencymatrix[n][state1][state2];
77 }
78
79 /** After finding a violation, construct and add all violated subtour cuts to scip */
80 static
addSubtourCuts(SCIP * scip,SCIP_SEPA * sepa,SCIP_Real *** adjacencymatrix,SCIP_DIGRAPH * adjacencygraph,int ** iscontracted,int cyclelength,SCIP_RESULT * result,int * ncuts)81 SCIP_RETCODE addSubtourCuts(
82 SCIP* scip, /**< SCIP data structure. */
83 SCIP_SEPA* sepa, /**< the subtour separator */
84 SCIP_Real*** adjacencymatrix, /**< the adjacency-matrices of all paths with 1,...,|Clutster| arcs */
85 SCIP_DIGRAPH* adjacencygraph, /**< the directed edge-graph */
86 int** iscontracted, /**< information of intermediate contraction-nodes for contracted arcs */
87 int cyclelength, /**< the length of the subtours to add */
88 SCIP_RESULT* result, /**< pointer to store the result of separation */
89 int* ncuts /**< pointer to store number of cuts */
90 )
91 {
92 SCIP_VAR**** edgevars;
93 char cutname[SCIP_MAXSTRLEN];
94 SCIP_ROW* cut;
95 int** subtours;
96 int* insubtour;
97 int* successors;
98 int nsuccessors;
99 int nstates;
100 int currentnode;
101 int successor;
102 int intermediate;
103 int anchor;
104 int ncontractions;
105 int liftabley;
106 int liftablez;
107 int greater;
108 int smaller;
109 int c;
110 int k;
111 int l;
112 SCIP_Bool isduplicate;
113
114 edgevars = SCIPcycGetEdgevars(scip);
115 nstates = SCIPdigraphGetNNodes(adjacencygraph);
116
117 SCIP_CALL( SCIPallocBlockMemoryArray(scip, &insubtour, nstates) );
118 SCIP_CALL( SCIPallocBlockMemoryArray(scip, &subtours, nstates) );
119
120 for( k = 0; k < nstates; ++k )
121 {
122 SCIP_CALL( SCIPallocClearBlockMemoryArray(scip, &subtours[k], cyclelength + 1) ); /*lint !e866, !e776*/
123 insubtour[k] = -1;
124 }
125
126 /* for each state, check if a subtour inequality is violated */
127 for( anchor = 0; anchor < nstates; ++anchor )
128 {
129 /* while reconstructing the subtour, count the number of contractions */
130 ncontractions = 0;
131
132 /* a cycle inequality is violated if the following is true */
133 if( SCIPisGT(scip, getDist(adjacencymatrix, cyclelength - 1, anchor, anchor), cyclelength - 1.0) )
134 {
135 subtours[anchor][0] = anchor;
136 if( insubtour[anchor] == -1 )
137 insubtour[anchor] = anchor;
138
139 /* traverse the cycle */
140 for( k = 0; k < cyclelength -1; ++k )
141 {
142 currentnode = subtours[anchor][k];
143
144 assert(0 <= currentnode && currentnode < nstates);
145
146 successors = SCIPdigraphGetSuccessors(adjacencygraph, currentnode);
147 nsuccessors = SCIPdigraphGetNSuccessors(adjacencygraph, currentnode);
148
149 /* find the next state along the subtour */
150 for( l = 0; l < nsuccessors; l++ )
151 {
152 successor = successors[l];
153
154 assert(0 <= successor && successor < nstates);
155
156 /* check if this successor of the current node is the one in the cycle. If so add it. */
157 if( SCIPisEQ(scip, getDist(adjacencymatrix, 0, currentnode, successor)
158 + getDist(adjacencymatrix, cyclelength - (k + 2), successor, anchor),
159 getDist(adjacencymatrix, cyclelength - (k + 1), currentnode, anchor)) )
160 {
161 subtours[anchor][k + 1] = successor;
162 insubtour[successor] = anchor;
163
164 if( iscontracted[currentnode][successor] != -1 )
165 ncontractions++;
166
167 break;
168 }
169 }
170 }
171
172 /* start and endnode are always the same in a cycle */
173 subtours[anchor][cyclelength] = anchor;
174
175 /* check last arc for a contraction */
176 if( iscontracted[subtours[anchor][cyclelength - 1]][anchor] != -1 )
177 ncontractions++;
178
179 isduplicate = FALSE;
180
181 /* if this anchor is already in another subtour, we check if the subtour is the same, since we don't want to
182 * add duplicates
183 */
184 if( insubtour[anchor] != anchor )
185 {
186 c = 0;
187 isduplicate = TRUE;
188
189 while( subtours[insubtour[anchor]][c] != anchor )
190 c++;
191
192 for( k = 0; k < cyclelength && isduplicate; ++k )
193 {
194 if( subtours[insubtour[anchor]][(k + c) % cyclelength] != subtours[anchor][k] )
195 isduplicate = FALSE;
196 }
197 }
198
199 if( isduplicate )
200 continue;
201
202 /* set the amount of y and z variables that we can still lift into the inequality */
203 liftabley = cyclelength - 1;
204 liftablez = SCIPcycGetNCluster(scip) - cyclelength - 1;
205
206 /* Now build the cut and add the subtour inequality */
207 (void)SCIPsnprintf(cutname, SCIP_MAXSTRLEN, "subtour_%d_length_%d_contracted_%d", anchor,
208 cyclelength, ncontractions );
209 SCIP_CALL( SCIPcreateEmptyRowSepa(scip, &cut,sepa, cutname, -SCIPinfinity(scip),
210 cyclelength + ncontractions - 1.0, FALSE, FALSE, TRUE) );
211
212 SCIP_CALL( SCIPcacheRowExtensions(scip, cut) );
213
214 for( k = 0; k < cyclelength; ++k )
215 {
216 currentnode = subtours[anchor][k];
217 successor = subtours[anchor][k+1];
218 intermediate = iscontracted[currentnode][successor];
219
220 if( intermediate != -1 )
221 {
222 SCIP_CALL( SCIPaddVarToRow(scip, cut, getEdgevar(edgevars, currentnode, intermediate, 1), 1.0) );
223 SCIP_CALL( SCIPaddVarToRow(scip, cut,
224 getEdgevar(edgevars, MAX(intermediate, successor), MIN(intermediate, successor), 0), 1.0) );
225
226 greater = intermediate > currentnode ? intermediate : currentnode;
227 smaller = intermediate < currentnode ? intermediate : currentnode;
228
229 if( liftabley > 0 && SCIPvarGetLPSol(getEdgevar(edgevars, greater, smaller, 0)) > 0 )
230 {
231 SCIP_CALL( SCIPaddVarToRow(scip, cut, getEdgevar(edgevars, greater, smaller, 0), 1.0) );
232 liftabley--;
233 }
234 if( liftablez > 0 && SCIPvarGetLPSol(getEdgevar(edgevars, intermediate, successor, 1)) > 0 )
235 {
236 SCIP_CALL( SCIPaddVarToRow(scip, cut, getEdgevar(edgevars, intermediate, successor, 1), 1.0) );
237 liftablez--;
238 }
239 }
240 else
241 {
242 SCIP_CALL( SCIPaddVarToRow(scip, cut, getEdgevar(edgevars, currentnode, successor, 1), 1.0) );
243 if( SCIPvarGetLPSol(getEdgevar(edgevars, MAX(currentnode, successor), MIN(currentnode, successor), 0))
244 > 0 && liftabley > 0 )
245 {
246 SCIP_CALL( SCIPaddVarToRow(scip, cut,
247 getEdgevar(edgevars, MAX(currentnode, successor), MIN(currentnode, successor), 0), 1.0) );
248 liftabley--;
249 }
250 }
251 }
252
253 SCIP_CALL( SCIPflushRowExtensions(scip, cut) );
254 SCIP_CALL( SCIPaddPoolCut(scip, cut) );
255
256 /* print for debugging purposes */
257 SCIPdebug( SCIP_CALL( SCIPprintRow(scip, cut, NULL) ) );
258
259 /* release data and increment cut counter */
260 SCIP_CALL( SCIPreleaseRow(scip, &cut) );
261
262 *result = SCIP_SEPARATED;
263 (*ncuts)++;
264 }
265 }
266
267 for( k = 0; k < nstates; ++k )
268 {
269 SCIPfreeBlockMemoryArray(scip, &(subtours[k]), cyclelength + 1);
270 }
271 SCIPfreeBlockMemoryArray(scip, &subtours, nstates);
272 SCIPfreeBlockMemoryArray(scip, &insubtour, nstates);
273
274 return SCIP_OKAY;
275 }
276
277 /** Detect if path inequalities are violated and if so, add them to scip */
278 static
addPathCuts(SCIP * scip,SCIP_SEPA * sepa,SCIP_Real *** adjacencymatrix,SCIP_DIGRAPH * adjacencygraph,int ** iscontracted,int pathlength,SCIP_RESULT * result,int * ncuts)279 SCIP_RETCODE addPathCuts(
280 SCIP* scip, /**< SCIP data structure. */
281 SCIP_SEPA* sepa, /**< the subtour separator */
282 SCIP_Real*** adjacencymatrix, /**< the adjacency-matrix of all paths with 1,...,|Clutster| arcs */
283 SCIP_DIGRAPH* adjacencygraph, /**< the directed edge-graph */
284 int** iscontracted, /**< information of intermediate contraction-nodes for contracted arcs */
285 int pathlength, /**< the length of the subtours to add */
286 SCIP_RESULT* result, /**< pointer to store the result of separation */
287 int* ncuts /**< pointer to store number of cuts */
288 )
289 {
290 SCIP_VAR**** edgevars;
291 char cutname[SCIP_MAXSTRLEN];
292 SCIP_ROW* cut;
293 int* path;
294 int nstates;
295 int currentnode;
296 int successor;
297 int* successors;
298 int nsuccessors;
299 int intermediate;
300 int start;
301 int end;
302 int ncontractions;
303 int k;
304 int i;
305 int j;
306 int nz;
307 int ny;
308
309 edgevars = SCIPcycGetEdgevars(scip);
310 nstates = SCIPdigraphGetNNodes(adjacencygraph);
311
312 SCIP_CALL( SCIPallocMemoryArray(scip, &path, pathlength + 1) );
313
314 for( start = 0; start < nstates; ++start )
315 {
316 path[0] = start;
317
318 for( j = 0; j < SCIPdigraphGetNSuccessors(adjacencygraph, start); ++j )
319 {
320 ncontractions = 0;
321
322 end = SCIPdigraphGetSuccessors(adjacencygraph, start)[j];
323 path[pathlength] = end;
324
325 /* check if path-inequality is violated */
326 if( SCIPisGT(scip, getDist(adjacencymatrix, pathlength - 1, start, end)
327 + getDist(adjacencymatrix, 0, start, end), (SCIP_Real) pathlength) )
328 {
329 /*reconstruct the path */
330 for( k = 0; k < pathlength - 1; ++k )
331 {
332 currentnode = path[k];
333
334 assert(0 <= currentnode && currentnode < nstates);
335
336 successors = SCIPdigraphGetSuccessors(adjacencygraph, currentnode);
337 nsuccessors = SCIPdigraphGetNSuccessors(adjacencygraph, currentnode);
338
339 for( i = 0; i < nsuccessors; ++i )
340 {
341 successor = successors[i];
342
343 assert(0 <= successor && successor < nstates);
344
345 if( SCIPisEQ(scip, getDist(adjacencymatrix, 0, currentnode, successor)
346 + getDist(adjacencymatrix, pathlength - (k + 2), successor, end),
347 getDist(adjacencymatrix, pathlength - (k + 1), currentnode, end)) )
348 {
349 path[k + 1] = successor;
350
351 if( iscontracted[currentnode][successor] != -1 )
352 ncontractions++;
353
354 break;
355 }
356 }
357 }
358
359 /* check the last arc along the path and the direct arc from start to end for contractions */
360 if( iscontracted[path[pathlength - 1]][end] != -1 )
361 ncontractions++;
362
363 if( iscontracted[start][end] != -1 )
364 ncontractions++;
365
366 nz = pathlength;
367 ny = 0;
368
369 /* construct the corresponding inequality and add it to scip */
370 (void)SCIPsnprintf(cutname, SCIP_MAXSTRLEN, "path_%d_%d_length_%d_contracted_%d",
371 start, end, pathlength, ncontractions );
372 SCIP_CALL( SCIPcreateEmptyRowSepa(scip, &cut,sepa, cutname, -SCIPinfinity(scip),
373 (SCIP_Real) pathlength + ncontractions, FALSE, FALSE, TRUE) );
374
375 SCIP_CALL( SCIPcacheRowExtensions(scip, cut) );
376
377 for( k = 0; k < pathlength; ++k )
378 {
379 currentnode = path[k];
380 successor = path[k+1];
381 intermediate = iscontracted[currentnode][successor];
382
383 if( intermediate != -1 )
384 {
385 SCIP_CALL( SCIPaddVarToRow(scip, cut, getEdgevar(edgevars, currentnode, intermediate, 1), 1.0) );
386 SCIP_CALL( SCIPaddVarToRow(scip, cut,
387 getEdgevar(edgevars, MAX(intermediate, successor), MIN(intermediate, successor), 0), 1.0) );
388
389 if( nz < SCIPcycGetNCluster(scip)
390 && SCIPisPositive(scip, SCIPvarGetLPSol(getEdgevar(edgevars, intermediate, successor, 1))) )
391 {
392 SCIP_CALL( SCIPaddVarToRow(scip, cut, getEdgevar(edgevars, intermediate, successor, 1), 1.0) );
393 nz++;
394 }
395
396 if( ny < pathlength - 2 && SCIPisPositive(scip, SCIPvarGetLPSol(
397 getEdgevar(edgevars, MAX(currentnode, intermediate), MIN(currentnode, intermediate), 0))) )
398 {
399 SCIP_CALL( SCIPaddVarToRow(scip, cut,
400 getEdgevar(edgevars, MAX(currentnode, intermediate), MIN(currentnode, intermediate), 0), 1.0) );
401 ny++;
402 }
403 }
404 else
405 {
406 SCIP_CALL( SCIPaddVarToRow(scip, cut, getEdgevar(edgevars, currentnode, successor, 1), 1.0) );
407
408 if( ny < pathlength - 2 && SCIPisPositive(scip, SCIPvarGetLPSol(
409 getEdgevar(edgevars, MAX(currentnode, successor), MIN(currentnode, successor), 0))) )
410 {
411 SCIP_CALL( SCIPaddVarToRow(scip, cut,
412 getEdgevar(edgevars, MAX(currentnode, successor), MIN(currentnode, successor), 0), 1.0) );
413 ny++;
414 }
415 }
416 }
417
418 /* add the direct arc from start to end */
419 intermediate = iscontracted[start][end];
420
421 if( iscontracted[start][end] != -1 )
422 {
423 SCIP_CALL( SCIPaddVarToRow(scip, cut, getEdgevar(edgevars, start, intermediate, 1), 1.0) );
424 SCIP_CALL( SCIPaddVarToRow(scip, cut, getEdgevar(edgevars,
425 MAX(intermediate, end), MIN(intermediate, end), 0), 1.0) );
426 }
427 else
428 {
429 assert( NULL != getEdgevar(edgevars, start, end, 1));
430
431 SCIP_CALL( SCIPaddVarToRow(scip, cut, getEdgevar(edgevars, start, end, 1), 1.0) );
432 }
433
434 SCIP_CALL( SCIPflushRowExtensions(scip, cut) );
435
436 /* print row if in debug mode */
437 SCIPdebug( SCIPprintRow(scip, cut, NULL) );
438
439 /* if an arc appears twice then the path inequality should not be used */
440 if( SCIPisEQ(scip, SCIPgetRowMaxCoef(scip, cut), 1.0) )
441 {
442 SCIP_CALL( SCIPaddPoolCut(scip, cut) );
443 *result = SCIP_SEPARATED;
444 (*ncuts)++;
445 }
446
447 SCIP_CALL( SCIPreleaseRow(scip, &cut) );
448 }
449 }
450 }
451
452 SCIPfreeMemoryArray(scip, &path);
453
454 return SCIP_OKAY;
455 }
456
457 /** detect if path inequalities are violated and if so, add them to scip */
458 static
addTourCuts(SCIP * scip,SCIP_SEPA * sepa,SCIP_Real *** adjacencymatrix,SCIP_DIGRAPH * adjacencygraph,int ** iscontracted,int tourlength,SCIP_RESULT * result,int * ncuts)459 SCIP_RETCODE addTourCuts(
460 SCIP* scip, /**< SCIP data structure. */
461 SCIP_SEPA* sepa, /**< the subtour separator */
462 SCIP_Real*** adjacencymatrix, /**< the adjacency-matrix of all paths with 1,...,|Clutster| arcs */
463 SCIP_DIGRAPH* adjacencygraph, /**< the directed edge-graph */
464 int** iscontracted, /**< information of intermediate contraction-nodes for contracted arcs */
465 int tourlength, /**< the length of the subtours to add */
466 SCIP_RESULT* result, /**< pointer to store the result of separation */
467 int* ncuts /**< pointer to store number of cuts */
468 )
469 {
470 SCIP_VAR**** edgevars;
471 char cutname[SCIP_MAXSTRLEN];
472 SCIP_ROW* cut;
473 int* tour;
474 int* successors;
475 int* succerssorsstart;
476 int nsuccessorsstart;
477 int nsuccessors;
478 int nstates;
479 int currentnode;
480 int successor;
481 int intermediate;
482 int start;
483 int end;
484 int ncontractions;
485 int k;
486 int i;
487 int j;
488
489 edgevars = SCIPcycGetEdgevars(scip);
490 nstates = SCIPdigraphGetNNodes(adjacencygraph);
491
492 SCIP_CALL( SCIPallocMemoryArray(scip, &tour, tourlength + 1) );
493
494 for( start = 0; start < nstates; ++start )
495 {
496 tour[0] = start;
497 succerssorsstart = SCIPdigraphGetSuccessors(adjacencygraph, start);
498 nsuccessorsstart = SCIPdigraphGetNSuccessors(adjacencygraph, start);
499
500 for( j = 0; j < nsuccessorsstart; ++j )
501 {
502 ncontractions = 0;
503
504 end = succerssorsstart[j];
505 tour[tourlength] = end;
506
507 /* check if tour-inequality is violated */
508 if( SCIPisGT(scip, getDist(adjacencymatrix, tourlength - 1, start, end)
509 - getDist(adjacencymatrix, 0, end, start), (SCIP_Real) tourlength - 1) )
510 {
511 /*reconstruct the tour */
512 for( k = 0; k < tourlength - 1; ++k )
513 {
514 currentnode = tour[k];
515 successors = SCIPdigraphGetSuccessors(adjacencygraph, currentnode);
516 nsuccessors = SCIPdigraphGetNSuccessors(adjacencygraph, currentnode);
517
518 for( i = 0; i < nsuccessors; ++i )
519 {
520 successor = successors[i];
521
522 if( SCIPisEQ(scip, getDist(adjacencymatrix, 0, currentnode, successor)
523 + getDist(adjacencymatrix, tourlength - (k + 2), successor, end)
524 , getDist(adjacencymatrix, tourlength - (k + 1), currentnode, end)) )
525 {
526 tour[k + 1] = successor;
527
528 if( iscontracted[currentnode][successor] != -1 )
529 ncontractions++;
530 break;
531 }
532 }
533 }
534
535 /* check the last arc along the tour and the direct arc from start to end for contractions */
536 if( iscontracted[tour[tourlength - 1]][end] != -1 )
537 ncontractions++;
538 if( iscontracted[end][start] != -1 )
539 ncontractions++;
540
541 /* construct the corresponding inequality and add it to scip */
542 (void)SCIPsnprintf(cutname, SCIP_MAXSTRLEN, "tour_%d_%d_length_%d_contracted_%d",
543 start, end, tourlength, ncontractions );
544 SCIP_CALL( SCIPcreateEmptyRowSepa(scip, &cut,sepa, cutname, -SCIPinfinity(scip),
545 (SCIP_Real) tourlength + ncontractions - 1, FALSE, FALSE, TRUE) );
546
547 SCIP_CALL( SCIPcacheRowExtensions(scip, cut) );
548
549 for( k = 0; k < tourlength; ++k )
550 {
551 currentnode = tour[k];
552 successor = tour[k+1];
553 intermediate = iscontracted[currentnode][successor];
554
555 if( intermediate != -1 )
556 {
557 SCIP_CALL( SCIPaddVarToRow(scip, cut, getEdgevar(edgevars, currentnode, intermediate, 1), 1.0) );
558 SCIP_CALL( SCIPaddVarToRow(scip, cut,
559 getEdgevar(edgevars, MAX(intermediate, successor), MIN(intermediate, successor), 0), 1.0) );
560 }
561 else
562 {
563 SCIP_CALL( SCIPaddVarToRow(scip, cut, getEdgevar(edgevars, currentnode, successor, 1), 1.0) );
564 }
565 }
566
567 /* add the direct arc from start to end */
568 intermediate = iscontracted[end][start];
569 if( iscontracted[end][start] != -1 )
570 {
571 SCIP_CALL( SCIPaddVarToRow(scip, cut, getEdgevar(edgevars, end, intermediate, 1), -1.0) );
572 SCIP_CALL( SCIPaddVarToRow(scip, cut,
573 getEdgevar(edgevars, MAX(intermediate, start), MIN(intermediate, start), 0), 1.0) );
574 }
575 else
576 {
577 SCIP_CALL( SCIPaddVarToRow(scip, cut, getEdgevar(edgevars, end, start, 1), -1.0) );
578 }
579
580 SCIP_CALL( SCIPflushRowExtensions(scip, cut) );
581
582 /* print row if in debug mode */
583 SCIPdebug( SCIPprintRow(scip, cut, NULL) );
584
585 /* if an arc appears twice then the tour inequality should not be used */
586 if( SCIPisEQ(scip, SCIPgetRowMaxCoef(scip, cut), 1.0) )
587 {
588 SCIP_CALL( SCIPaddPoolCut(scip, cut) );
589 *result = SCIP_SEPARATED;
590 (*ncuts)++;
591 }
592
593 SCIP_CALL( SCIPreleaseRow(scip, &cut) );
594 }
595 }
596 }
597
598 SCIPfreeMemoryArray(scip, &tour);
599
600 return SCIP_OKAY;
601 }
602
603 /** compute the next matrix with the weight off all the longest paths with exactly narcs and store it in
604 * adjacencymatrix[narcs - 1]. For this, simply compute
605 * \f{align*}{ d^{k}(currentnode,successor) = max_{l=1,\ldots,n} \{d^{k-1}(currentnode,l) + d^1(l,successor) \} \f}.
606 */
607 static
computeNextAdjacency(SCIP * scip,SCIP_Real *** adjacencymatrix,SCIP_DIGRAPH * adjacencygraph,int narcs)608 SCIP_Bool computeNextAdjacency
609 (
610 SCIP* scip, /**< SCIP data structure */
611 SCIP_Real*** adjacencymatrix, /**< the max-distance matrices for all number of arcs less than narcs. */
612 SCIP_DIGRAPH* adjacencygraph, /**< the directed edge-graph */
613 int narcs /**< the current number of arcs in the paths */
614 )
615 {
616 int* intermediates;
617 int nintermediates;
618 int currentnode;
619 int intermediate;
620 int successor;
621 int l;
622 int nnodes;
623 SCIP_Bool foundviolation;
624
625 foundviolation = FALSE;
626 nnodes = SCIPdigraphGetNNodes(adjacencygraph);
627
628 for( currentnode = 0; currentnode < nnodes; ++currentnode )
629 {
630 intermediates = SCIPdigraphGetSuccessors(adjacencygraph, currentnode);
631 nintermediates = SCIPdigraphGetNSuccessors(adjacencygraph, currentnode);
632
633 for( l = 0; l < nintermediates; ++l )
634 {
635 intermediate = intermediates[l];
636
637 assert(0 <= intermediate && intermediate < nnodes);
638
639 for( successor = 0; successor < nnodes; ++successor )
640 {
641 if( SCIPisPositive(scip, getDist(adjacencymatrix, 0, currentnode, intermediate))
642 && SCIPisPositive(scip, getDist(adjacencymatrix, narcs - 2, intermediate, successor)) )
643 {
644 if( SCIPisGT(scip, getDist(adjacencymatrix, 0, currentnode, intermediate)
645 + getDist(adjacencymatrix, narcs - 2, intermediate, successor),
646 getDist(adjacencymatrix, narcs - 1, currentnode, successor)) )
647 {
648 adjacencymatrix[narcs - 1][currentnode][successor] = getDist(adjacencymatrix, 0, currentnode, intermediate)
649 + getDist(adjacencymatrix, narcs - 2, intermediate, successor);
650 }
651 }
652 }
653 }
654 }
655
656 /* check if we have found a violated subtour constraint */
657 for( currentnode = 0; currentnode < nnodes; ++currentnode )
658 {
659 if( SCIPisGT(scip, getDist(adjacencymatrix, narcs - 1, currentnode, currentnode), narcs - 1.0) )
660 foundviolation = TRUE;
661 }
662 return foundviolation;
663 }
664
665 /** copy method for separator plugins (called when SCIP copies plugins) */
666 static
SCIP_DECL_SEPACOPY(sepaCopySubtour)667 SCIP_DECL_SEPACOPY(sepaCopySubtour)
668 { /*lint --e{715}*/
669 assert(scip != NULL);
670 assert(sepa != NULL);
671 assert(strcmp(SCIPsepaGetName(sepa), SEPA_NAME) == 0);
672
673 /* call inclusion method of constraint handler */
674 SCIP_CALL( SCIPincludeSepaSubtour(scip) );
675
676 return SCIP_OKAY;
677 }
678
679 /** LP solution separation method of separator */
680 static
SCIP_DECL_SEPAEXECLP(sepaExeclpSubtour)681 SCIP_DECL_SEPAEXECLP(sepaExeclpSubtour)
682 { /*lint --e{715}*/
683 SCIP_VAR**** edgevars;
684 SCIP_Real*** adjacencymatrix;
685 SCIP_DIGRAPH* adjacencygraph;
686 SCIP_DIGRAPH* edgegraph;
687 int** iscontracted;
688 SCIP_Bool violation;
689 int* successors1;
690 int* successors2;
691 int nsuccessors1;
692 int nsuccessors2;
693 int ncuts;
694 int nstates;
695 int ncluster;
696 int cyclelength;
697 int rounds;
698 int i;
699 int j;
700 int k;
701 int state1;
702 int state2;
703 int state3;
704
705 /* get problem information */
706 rounds = SCIPsepaGetNCallsAtNode(sepa);
707 ncluster = SCIPcycGetNCluster(scip);
708 edgevars = SCIPcycGetEdgevars(scip);
709 nstates = SCIPcycGetNBins(scip);
710 edgegraph = SCIPcycGetEdgeGraph(scip);
711 ncuts = 0;
712
713 if( rounds >= MAXROUNDS )
714 {
715 *result = SCIP_DIDNOTRUN;
716 return SCIP_OKAY;
717 }
718
719 assert(nstates > 0);
720 assert(ncluster > 0 && ncluster < nstates);
721 assert(NULL != edgevars);
722 assert(NULL != edgegraph);
723
724 /* allocate memory */
725 SCIP_CALL( SCIPallocBlockMemoryArray(scip, &adjacencymatrix, ncluster) );
726
727 for( k = 0; k < ncluster; ++k )
728 {
729 SCIP_CALL( SCIPallocBlockMemoryArray(scip, &adjacencymatrix[k], nstates) ); /*lint !e866*/
730
731 for( j = 0; j < nstates; ++j )
732 {
733 SCIP_CALL( SCIPallocClearBlockMemoryArray(scip, &adjacencymatrix[k][j], nstates) ); /*lint !e866*/
734 }
735 }
736
737 /* create Digraph from the current LP-Solution */
738 SCIP_CALL( SCIPcreateDigraph(scip, &adjacencygraph, nstates) );
739 SCIP_CALL( SCIPallocBlockMemoryArray(scip, &iscontracted, nstates) );
740
741
742 /* get the values of the lp-solution */
743 for( i = 0; i < nstates; ++i )
744 {
745 SCIP_CALL( SCIPallocBlockMemoryArray(scip, &iscontracted[i], nstates) );
746
747 for( j = 0; j < nstates; ++j )
748 {
749 iscontracted[i][j] = -1;
750
751 if( edgevars[i] != NULL && edgevars[i][j] != NULL && getEdgevar(edgevars, i, j, 1) != NULL )
752 adjacencymatrix[0][i][j] = SCIPvarGetLPSol(getEdgevar(edgevars, i, j, 1));
753 }
754 }
755
756 /* contract the adjacency matrix if it is better to take z_{ij} + y_{jk} rather than z_{ik} directly,
757 * this stores j at position (i,k)
758 */
759 for( i = 0; i < nstates; ++i )
760 {
761 state1 = i;
762
763 assert( edgevars[state1] != NULL);
764
765 successors1 = SCIPdigraphGetSuccessors(edgegraph, state1);
766 nsuccessors1 = SCIPdigraphGetNSuccessors(edgegraph, state1);
767
768 for( j = 0; j < nsuccessors1; ++j )
769 {
770 state2 = successors1[j];
771
772 assert( edgevars[state2] != NULL);
773
774 successors2 = SCIPdigraphGetSuccessors(edgegraph, state2);
775 nsuccessors2 = SCIPdigraphGetNSuccessors(edgegraph, state2);
776
777 for( k = 0 ; k < nsuccessors2; ++k )
778 {
779 state3 = successors2[k];
780
781 if( edgevars[state1][state2] == NULL || edgevars[state2][state3] == NULL || edgevars[state1][state3] == NULL )
782 continue;
783
784 if( SCIPisLT( scip, getDist(adjacencymatrix, 0, state1, state3),
785 SCIPvarGetLPSol(getEdgevar(edgevars, state1, state2, 1))
786 + SCIPvarGetLPSol(getEdgevar(edgevars, MAX(state2, state3), MIN(state2, state3), 0)) - 1) )
787 {
788 adjacencymatrix[0][state1][state3] = SCIPvarGetLPSol(getEdgevar(edgevars, state1, state2, 1))
789 + SCIPvarGetLPSol(getEdgevar(edgevars, MAX(state2, state3), MIN(state2, state3), 0)) - 1;
790
791 iscontracted[state1][state3] = state2;
792 }
793 }
794 }
795 }
796
797 /* save the contracted matrix as a digraph to be able to reuse it quicker */
798 for( i = 0; i < nstates; ++i )
799 {
800 for( j = 0; j < nstates; ++j )
801 {
802 if( !SCIPisZero(scip, getDist(adjacencymatrix, 0, i, j)) )
803 {
804 SCIP_CALL( SCIPdigraphAddArc(adjacencygraph, i , j, NULL) );
805 }
806 }
807 }
808
809 /* a cyclelength of one does not make sense as there are no loops */
810 cyclelength = 2;
811 *result = SCIP_DIDNOTFIND;
812
813 /* Iterate until we have found a sufficient number of cuts or until we have checked all possible violations */
814 while( cyclelength < ncluster )
815 {
816 /* Compute the next adjacency matrix */
817 violation = computeNextAdjacency(scip, adjacencymatrix, adjacencygraph, cyclelength);
818
819 /* if we found a violation separate it */
820 if( violation )
821 {
822 SCIP_CALL( addSubtourCuts(scip, sepa, adjacencymatrix, adjacencygraph, iscontracted, cyclelength,
823 result, &ncuts) );
824 }
825
826 /* check if any path-inequalities are violated and sepatare them */
827 SCIP_CALL( addPathCuts(scip, sepa, adjacencymatrix, adjacencygraph, iscontracted, cyclelength, result, &ncuts) );
828
829 if( cyclelength == ncluster - 1 )
830 {
831 SCIP_CALL( addTourCuts(scip, sepa, adjacencymatrix, adjacencygraph, iscontracted, cyclelength,
832 result, &ncuts) );
833 }
834
835 /* stop if we added maximal number of cuts */
836 if( ncuts >= MAXCUTS )
837 break;
838
839 cyclelength++;
840 }
841
842 SCIPdigraphFreeComponents(adjacencygraph);
843 SCIPdigraphFree(&adjacencygraph);
844
845 /* free allocated memory */
846 for( i = 0; i < nstates; ++i )
847 {
848 SCIPfreeBlockMemoryArray(scip, &iscontracted[i], nstates);
849 }
850 SCIPfreeBlockMemoryArray(scip, &iscontracted, nstates);
851
852 for( i = 0; i < ncluster; ++i )
853 {
854 for( j = 0; j < nstates; ++j )
855 {
856 SCIPfreeBlockMemoryArray(scip, &adjacencymatrix[i][j], nstates);
857 }
858 SCIPfreeBlockMemoryArray(scip, &adjacencymatrix[i], nstates);
859 }
860 SCIPfreeBlockMemoryArray(scip, &adjacencymatrix, ncluster);
861
862 return SCIP_OKAY;
863 }
864
865
866 /** creates the Subtour separator and includes it in SCIP */
SCIPincludeSepaSubtour(SCIP * scip)867 SCIP_RETCODE SCIPincludeSepaSubtour(
868 SCIP* scip /**< SCIP data structure */
869 )
870 {
871 SCIP_SEPA* sepa;
872
873 /* include separator */
874
875 SCIP_CALL( SCIPincludeSepaBasic(scip, &sepa, SEPA_NAME, SEPA_DESC, SEPA_PRIORITY, SEPA_FREQ, SEPA_MAXBOUNDDIST,
876 SEPA_USESSUBSCIP, SEPA_DELAY,
877 sepaExeclpSubtour, NULL,
878 NULL) );
879
880 assert(sepa != NULL);
881
882 /* set non fundamental callbacks via setter functions */
883 SCIP_CALL( SCIPsetSepaCopy(scip, sepa, sepaCopySubtour) );
884
885
886 return SCIP_OKAY;
887 }
888