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