Scippy

SCIP

Solving Constraint Integer Programs

sepa_subtour.c
Go to the documentation of this file.
1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2 /* */
3 /* This file is part of the program and library */
4 /* SCIP --- Solving Constraint Integer Programs */
5 /* */
6 /* Copyright (C) 2002-2020 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
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]);x
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
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
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
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
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
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
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 */
675 
676  return SCIP_OKAY;
677 }
678 
679 /** LP solution separation method of separator */
680 static
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 */
868  SCIP* scip /**< SCIP data structure */
869 )
870 {
871  SCIP_SEPA* sepa;
872 
873  /* include separator */
874 
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 }
enum SCIP_Result SCIP_RESULT
Definition: type_result.h:52
SCIP_Bool SCIPisEQ(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
#define SCIPfreeBlockMemoryArray(scip, ptr, num)
Definition: scip_mem.h:97
static SCIP_Bool computeNextAdjacency(SCIP *scip, SCIP_Real ***adjacencymatrix, SCIP_DIGRAPH *adjacencygraph, int narcs)
Definition: sepa_subtour.c:609
SCIP_Bool SCIPisPositive(SCIP *scip, SCIP_Real val)
#define narcs
Definition: gastrans.c:68
#define SCIPallocBlockMemoryArray(scip, ptr, num)
Definition: scip_mem.h:80
SCIP_RETCODE SCIPcreateDigraph(SCIP *scip, SCIP_DIGRAPH **digraph, int nnodes)
int SCIPdigraphGetNNodes(SCIP_DIGRAPH *digraph)
Definition: misc.c:7648
Separate Subtours-Elimination inequalities in Cycle-Clustering Applications.
#define SCIPfreeMemoryArray(scip, ptr)
Definition: scip_mem.h:69
SCIP_VAR * getEdgevar(SCIP_VAR ****edgevars, int state1, int state2, int direction)
#define SCIP_MAXSTRLEN
Definition: def.h:273
#define SCIPallocClearBlockMemoryArray(scip, ptr, num)
Definition: scip_mem.h:84
SCIP_EXPORT SCIP_Real SCIPvarGetLPSol(SCIP_VAR *var)
Definition: var.c:18041
#define SCIPallocMemoryArray(scip, ptr, num)
Definition: scip_mem.h:53
SCIP_RETCODE SCIPprintRow(SCIP *scip, SCIP_ROW *row, FILE *file)
Definition: scip_lp.c:2152
SCIP_VAR **** SCIPcycGetEdgevars(SCIP *scip)
#define FALSE
Definition: def.h:73
#define TRUE
Definition: def.h:72
#define SCIPdebug(x)
Definition: pub_message.h:84
enum SCIP_Retcode SCIP_RETCODE
Definition: type_retcode.h:54
int SCIPdigraphGetNSuccessors(SCIP_DIGRAPH *digraph, int node)
Definition: misc.c:7706
void SCIPdigraphFree(SCIP_DIGRAPH **digraph)
Definition: misc.c:7470
SCIP_RETCODE SCIPsetSepaCopy(SCIP *scip, SCIP_SEPA *sepa, SCIP_DECL_SEPACOPY((*sepacopy)))
Definition: scip_sepa.c:142
SCIP_RETCODE SCIPflushRowExtensions(SCIP *scip, SCIP_ROW *row)
Definition: scip_lp.c:1604
SCIP_VAR ** x
Definition: circlepacking.c:54
static SCIP_RETCODE addTourCuts(SCIP *scip, SCIP_SEPA *sepa, SCIP_Real ***adjacencymatrix, SCIP_DIGRAPH *adjacencygraph, int **iscontracted, int tourlength, SCIP_RESULT *result, int *ncuts)
Definition: sepa_subtour.c:459
static SCIP_DECL_SEPAEXECLP(sepaExeclpSubtour)
Definition: sepa_subtour.c:681
static SCIP_Real getDist(SCIP_Real ***adjacencymatrix, int n, int state1, int state2)
Definition: sepa_subtour.c:66
SCIP_RETCODE SCIPcacheRowExtensions(SCIP *scip, SCIP_ROW *row)
Definition: scip_lp.c:1581
#define SEPA_FREQ
Definition: sepa_subtour.c:36
SCIP_RETCODE SCIPreleaseRow(SCIP *scip, SCIP_ROW **row)
Definition: scip_lp.c:1508
#define SEPA_PRIORITY
Definition: sepa_subtour.c:35
SCIP_Bool SCIPisZero(SCIP *scip, SCIP_Real val)
#define NULL
Definition: lpi_spx1.cpp:155
#define SCIP_CALL(x)
Definition: def.h:364
SCIP_EXPORT const char * SCIPsepaGetName(SCIP_SEPA *sepa)
Definition: sepa.c:697
#define SEPA_MAXBOUNDDIST
Definition: sepa_subtour.c:37
int SCIPcycGetNBins(SCIP *scip)
static SCIP_DECL_SEPACOPY(sepaCopySubtour)
Definition: sepa_subtour.c:667
SCIP_Real SCIPinfinity(SCIP *scip)
public data structures and miscellaneous methods
#define SCIP_Bool
Definition: def.h:70
#define MAXCUTS
Definition: sepa_subtour.c:40
SCIP_RETCODE SCIPincludeSepaBasic(SCIP *scip, SCIP_SEPA **sepa, const char *name, const char *desc, int priority, int freq, SCIP_Real maxbounddist, SCIP_Bool usessubscip, SCIP_Bool delay, SCIP_DECL_SEPAEXECLP((*sepaexeclp)), SCIP_DECL_SEPAEXECSOL((*sepaexecsol)), SCIP_SEPADATA *sepadata)
Definition: scip_sepa.c:100
int * SCIPdigraphGetSuccessors(SCIP_DIGRAPH *digraph, int node)
Definition: misc.c:7721
static SCIP_RETCODE addPathCuts(SCIP *scip, SCIP_SEPA *sepa, SCIP_Real ***adjacencymatrix, SCIP_DIGRAPH *adjacencygraph, int **iscontracted, int pathlength, SCIP_RESULT *result, int *ncuts)
Definition: sepa_subtour.c:279
#define MAX(x, y)
Definition: tclique_def.h:83
problem data for cycle clustering problem
Constraint handler for linear constraints in their most general form, .
#define SEPA_NAME
Definition: sepa_subtour.c:33
#define SEPA_DESC
Definition: sepa_subtour.c:34
int SCIPcycGetNCluster(SCIP *scip)
#define SEPA_USESSUBSCIP
Definition: sepa_subtour.c:38
SCIP_RETCODE SCIPaddVarToRow(SCIP *scip, SCIP_ROW *row, SCIP_VAR *var, SCIP_Real val)
Definition: scip_lp.c:1641
int SCIPsnprintf(char *t, int len, const char *s,...)
Definition: misc.c:10604
SCIP_DIGRAPH * SCIPcycGetEdgeGraph(SCIP *scip)
#define SCIP_Real
Definition: def.h:163
SCIP_Bool SCIPisLT(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
SCIP_Bool SCIPisGT(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
SCIP_RETCODE SCIPaddPoolCut(SCIP *scip, SCIP_ROW *row)
Definition: scip_cut.c:332
#define nnodes
Definition: gastrans.c:65
SCIP_Real SCIPgetRowMaxCoef(SCIP *scip, SCIP_ROW *row)
Definition: scip_lp.c:1862
void SCIPinfoMessage(SCIP *scip, FILE *file, const char *formatstr,...)
Definition: scip_message.c:199
SCIP_EXPORT int SCIPsepaGetNCallsAtNode(SCIP_SEPA *sepa)
Definition: sepa.c:824
static SCIP_RETCODE addSubtourCuts(SCIP *scip, SCIP_SEPA *sepa, SCIP_Real ***adjacencymatrix, SCIP_DIGRAPH *adjacencygraph, int **iscontracted, int cyclelength, SCIP_RESULT *result, int *ncuts)
Definition: sepa_subtour.c:81
#define MAXROUNDS
Definition: sepa_subtour.c:41
void SCIPdigraphFreeComponents(SCIP_DIGRAPH *digraph)
Definition: misc.c:8419
SCIP_RETCODE SCIPincludeSepaSubtour(SCIP *scip)
Definition: sepa_subtour.c:867
#define SEPA_DELAY
Definition: sepa_subtour.c:39
SCIP_RETCODE SCIPcreateEmptyRowSepa(SCIP *scip, SCIP_ROW **row, SCIP_SEPA *sepa, const char *name, SCIP_Real lhs, SCIP_Real rhs, SCIP_Bool local, SCIP_Bool modifiable, SCIP_Bool removable)
Definition: scip_lp.c:1399
SCIP_RETCODE SCIPdigraphAddArc(SCIP_DIGRAPH *digraph, int startnode, int endnode, void *data)
Definition: misc.c:7564