Scippy

SCIP

Solving Constraint Integer Programs

mincut.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-2022 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 scip.zib.de. */
13 /* */
14 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
15 
16 /**@file mincut.c
17  * @brief Minimum cut routines for Steiner problems
18  * @author Daniel Rehfeldt
19  *
20  * This file encompasses minimum cut routines for Steiner tree problems.
21  *
22  */
23 
24 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
25 
26 
27 //#define SCIP_DEBUG
28 
29 
30 #include "mincut.h"
31 #include "probdata_stp.h"
32 #include "misc_stp.h"
33 #include "stpvector.h"
34 #include "reduce.h"
35 #include "portab.h"
36 
37 #define Q_NULL -1 /* NULL element of queue/list */
38 #define ADDCUTSTOPOOL FALSE
39 #define TERMSEPA_SPARSE_MAXRATIO 4
40 #define TERMSEPA_MAXCUTSIZE_DEFAULT 5
41 #define TERMSEPA_MAXNCUTS_DEFAULT 100
42 #define TERMSEPA_NROOTCANDS 3
43 
44 /* *
45 #define FLOW_FACTOR 100000
46 #define CREEP_VALUE 1 this is the original value todo check what is better
47 */
48 
49 #define FLOW_FACTOR 1000000
50 #define CREEP_VALUE 10
51 
52 /** minimum cut helper */
53 typedef struct minimum_cut_helper
54 {
55  const SCIP_Real* xval;
57  int* edges_capa;
58  int* terms;
59  int* csr_start;
60  int* rootcut;
65  int* terms_minsepasize; /**< size of smallest separator for given terminal (non-defined for Steiner nodes) */
66  int* terms_mincompsize; /**< size of smallest component for given terminal (non-defined for Steiner nodes) */
67  STP_Bool* edges_isRemoved; /**< only used for LP cuts */
71  int root;
74  SCIP_RANDNUMGEN* randnumgen; /**< random number generator or NULL */
75  SCIP_Bool isLpcut; /**< cut for LP? */
76 } MINCUT;
77 
78 
79 
80 /** single separator info */
81 typedef struct terminal_separator
82 {
83  int sinkterm;
85 } TSEPA;
86 
87 
88 /** storage */
90 {
91  int* nsepas; /**< of size maxsepasize + 1 */
92  int* currsepa_n; /**< of size maxsepasize + 1 */
98  int root;
99  int maxnsepas; /**< maximum number of separators to compute */
100  int maxsepasize; /**< maximum size of separator */
101 };
102 
103 /*
104  * Local methods
105  */
106 
107 //#define STP_MAXFLOW_WRITE
108 
109 #ifdef STP_MAXFLOW_WRITE
110 /* writes flow problem in extended dimacs format */
111 static
112 void writeFlowProb(
113  const GRAPH* g, /**< graph data structure */
114  const int* capa, /**< edge capacity */
115  int sinkterm /**< sink */
116  )
117 {
118  FILE *fptr;
119  const int nnodes = graph_get_nNodes(g);
120  const int nedges = graph_get_nEdges(g);
121 
122  fptr = fopen("flow", "w+");
123  assert(fptr != NULL);
124 
125  fprintf(fptr, "p max %d %d \n", nnodes, nedges);
126  fprintf(fptr, "n %d s \n", g->source + 1);
127  fprintf(fptr, "n %d t \n", sinkterm + 1);
128 
129  for( int k = 0; k < nnodes; k++ )
130  {
131  for( int e = g->outbeg[k]; e != EAT_LAST; e = g->oeat[e] )
132  {
133  fprintf(fptr, "a %d %d %d \n", k + 1, g->head[e] + 1, capa[e]);
134  }
135  }
136 
137  fprintf(fptr, "x\n");
138 
139  fclose(fptr);
140 }
141 #endif
142 
143 
144 
145 #ifdef SCIP_DEBUG
146 static inline
147 void debugPrintCutNodes(
148  const GRAPH* g, /**< the graph */
149  const MINCUT* mincut /**< minimum cut */
150 )
151 {
152  const int nnodes = graph_get_nNodes(g);
153  const int* const nodes_wakeState = mincut->nodes_wakeState;
154 
155  printf("root component nodes: \n");
156 
157  for( int i = 0; i < nnodes; i++ )
158  {
159  if( nodes_wakeState[i] )
160  graph_knot_printInfo(g, i);
161  }
162 
163  printf("remaining nodes: \n");
164 
165  for( int i = 0; i < nnodes; i++ )
166  {
167  if( !nodes_wakeState[i] )
168  graph_knot_printInfo(g, i);
169  }
170 }
171 
172 static inline
173 void debugPrintCutEdges(
174  const GRAPH* g, /**< the graph */
175  const MINCUT* mincut /**< minimum cut */
176 )
177 {
178  const int nnodes = graph_get_nNodes(g);
179  const int* const nodes_wakeState = mincut->nodes_wakeState;
180 
181  printf("cut edges: \n");
182 
183  for( int i = 0; i < nnodes; i++ )
184  {
185  if( !nodes_wakeState[i] )
186  continue;
187 
188  for( int e = g->outbeg[i]; e != EAT_LAST; e = g->oeat[e] )
189  {
190  const int head = g->head[e];
191 
192  if( nodes_wakeState[head] == 0 )
193  {
194  graph_edge_printInfo(g, e);
195  }
196  }
197  }
198 }
199 
200 
201 static inline
202 void debugPrintCsrCutEdges(
203  const GRAPH* g, /**< the graph */
204  const MINCUT* mincut /**< minimum cut */
205 )
206 {
207  const int* const edges_capa = mincut->edges_capa;
208  const int* const nodes_wakeState = mincut->nodes_wakeState;
209  const int* const csr_start = mincut->csr_start;
210  const int* const csr_headarr = mincut->csr_headarr;
211  const int nnodes_extended = mincut->termsepa_nnodes;
212 
213  assert(!mincut->isLpcut);
214 
215  printf("CSR cut edges: \n");
216 
217  for( int i = 0; i < nnodes_extended; i++ )
218  {
219  const int start = csr_start[i];
220  const int end = csr_start[i + 1];
221 
222  if( !nodes_wakeState[i] )
223  continue;
224 
225  for( int j = start; j != end; j++ )
226  {
227  const int head = csr_headarr[j];
228 
229  if( nodes_wakeState[head] == 0 )
230  {
231  printf("%d->%d, c=%d \n", i, head, edges_capa[j]);
232  }
233  }
234  }
235 }
236 
237 
238 
239 /** prints extended graph */
240 static inline
241 void debugPrintCsr(
242  const GRAPH* g, /**< the graph */
243  const MINCUT* mincut /**< minimum cut */
244 )
245 {
246  const int* const residual = g->mincut_r;
247  const int* const csr_start = mincut->csr_start;
248  const int* const csr_headarr = mincut->csr_headarr;
249  const int* const csr_edgeflipped = mincut->csr_edgeflipped;
250  const int nnodes_extended = mincut->termsepa_nnodes;
251 
252  assert(!mincut->isLpcut);
253 
254  for( int i = 0; i < nnodes_extended; i++ )
255  {
256  const int start = csr_start[i];
257  const int end = csr_start[i + 1];
258 
259  for( int j = start; j != end; j++ )
260  {
261  const int head = csr_headarr[j];
262 
263  printf("edge %d: %d->%d, r=%d (fliped=%d) \n", j, i, head, residual[j], csr_edgeflipped[j]);
264  }
265  }
266 }
267 #endif
268 
269 
270 #ifndef NDEBUG
271 /** valid flip-edges? */
272 static inline
274  const GRAPH* g, /**< the graph */
275  const MINCUT* mincut /**< minimum cut */
276 )
277 {
278  const int* const residual = g->mincut_r;
279  const int* const csr_start = mincut->csr_start;
280  const int* const csr_headarr = mincut->csr_headarr;
281  const int* const csr_edgeflipped = mincut->csr_edgeflipped;
282  const int nnodes_extended = mincut->termsepa_nnodes;
283 
284  assert(!mincut->isLpcut);
285 
286  for( int i = 0; i < nnodes_extended; i++ )
287  {
288  const int start = csr_start[i];
289  const int end = csr_start[i + 1];
290 
291  for( int j = start; j != end; j++ )
292  {
293  assert(csr_edgeflipped[j] >= 0);
294 
295  if( csr_headarr[csr_edgeflipped[j]] != i )
296  {
297  return FALSE;
298  }
299 
300  assert(residual[j] > 0 || residual[csr_edgeflipped[j]] > 0);
301  }
302  }
303 
304  return TRUE;
305 }
306 
307 
308 /** checks cut */
309 static
311  SCIP* scip, /**< SCIP data structure */
312  const GRAPH* g, /**< the graph */
313  int ncutterms, /**< number of cut terminals */
314  const int* cutterms, /**< stores cut terminals */
315  int sinkterm, /**< sink terminal */
316  const MINCUT* mincut /**< minimum cut */
317 )
318 {
319  SCIP_Bool isCorrect = TRUE;
320  const int* const nodes_wakeState = mincut->nodes_wakeState;
321  STP_Vectype(int) stack = NULL;
322  SCIP_Bool* nodes_isBlocked;
323  SCIP_Bool* nodes_isVisited;
324  const int nnodes = graph_get_nNodes(g);
325  int nvisted = 1;
326  int nsinknodes = 0;
327 
328  assert(ncutterms > 0);
329  assert(nodes_wakeState[sinkterm] == 0);
330 
331  SCIP_CALL_ABORT( SCIPallocMemoryArray(scip, &nodes_isBlocked, nnodes) );
332  SCIP_CALL_ABORT( SCIPallocMemoryArray(scip, &nodes_isVisited, nnodes) );
333 
334  for( int i = 0; i < nnodes; i++ )
335  {
336  nodes_isBlocked[i] = FALSE;
337  nodes_isVisited[i] = FALSE;
338  }
339 
340  for( int i = 0; i < ncutterms; i++ )
341  {
342  assert(graph_knot_isInRange(g, cutterms[i]));
343  assert(!nodes_isBlocked[cutterms[i]]);
344 
345  nodes_isBlocked[cutterms[i]] = TRUE;
346  }
347 
348  assert(!nodes_isBlocked[mincut->root]);
349 
350  for( int i = 0; i < nnodes; i++ )
351  {
352  if( nodes_wakeState[i] == 0 )
353  {
354  nsinknodes++;
355  }
356  }
357 
358  /* DFS from sink terminal */
359  StpVecReserve(scip, stack, nnodes);
360  StpVecPushBack(scip, stack, sinkterm);
361  nodes_isVisited[sinkterm] = TRUE;
362 
363  while( StpVecGetSize(stack) )
364  {
365  const int k = stack[StpVecGetSize(stack) - 1];
366  StpVecPopBack(stack);
367 
368  assert(nodes_isVisited[k]);
369 
370  for( int e = g->outbeg[k]; e != EAT_LAST; e = g->oeat[e] )
371  {
372  const int head = g->head[e];
373  if( !nodes_isVisited[head] && !nodes_isBlocked[head] )
374  {
375  nvisted++;
376  nodes_isVisited[head] = TRUE;
377  StpVecPushBack(scip, stack, head);
378  }
379  }
380  }
381 
382  StpVecFree(scip, stack);
383 
384  if( nvisted > nsinknodes )
385  {
386  SCIPerrorMessage("nvisted=%d nsinknodes=%d\n", nvisted, nsinknodes);
387  isCorrect = FALSE;
388  }
389 
390  if( nodes_isVisited[mincut->root] )
391  {
392  SCIPerrorMessage("root is in cut! \n");
393  isCorrect = FALSE;
394  }
395 
396 
397  SCIPfreeMemoryArray(scip, &nodes_isVisited);
398  SCIPfreeMemoryArray(scip, &nodes_isBlocked);
399 
400 
401  return isCorrect;
402 }
403 #endif
404 
405 
406 /** gets infinity */
407 static
409  const GRAPH* g, /**< the graph */
410  const MINCUT* mincut /**< minimum cut */
411 )
412 {
413  return g->terms;
414 }
415 
416 
417 /** updates; takes root candidate if greater than current */
418 static
420  SCIP* scip, /**< SCIP data structure */
421  int rootcand, /**< new root candidate */
422  const GRAPH* g, /**< the graph */
423  int* root, /**< incumbent (IN/OUT) */
424  int* rootcompsize /**< size of component (IN/OUT) */
425 )
426 {
427  STP_Vectype(int) queue = NULL;
428  SCIP_Bool* nodes_isVisited;
429  const int nnodes = graph_get_nNodes(g);
430  int compsize = 0;
431 
432  assert(graph_knot_isInRange(g, rootcand) && Is_term(g->term[rootcand]));
433 
434  SCIP_CALL_ABORT( SCIPallocCleanBufferArray(scip, &nodes_isVisited, nnodes) );
435  StpVecReserve(scip, queue, nnodes);
436  StpVecPushBack(scip, queue, rootcand);
437  nodes_isVisited[rootcand] = TRUE;
438 
439  /* BFS loop stopping at roots */
440  for( int i = 0; i < StpVecGetSize(queue); i++ )
441  {
442  const int node = queue[i];
443  assert(nodes_isVisited[node]);
444 
445  compsize++;
446 
447  if( Is_term(g->term[node]) && node != rootcand )
448  continue;
449 
450  for( int e = g->outbeg[node]; e >= 0; e = g->oeat[e] )
451  {
452  const int head = g->head[e];
453  if( !nodes_isVisited[head] )
454  {
455  nodes_isVisited[head] = TRUE;
456  StpVecPushBack(scip, queue, head);
457  }
458  }
459  }
460 
461  for( int i = 0; i < StpVecGetSize(queue); i++ )
462  {
463  const int node = queue[i];
464  assert(nodes_isVisited[node]);
465  nodes_isVisited[node] = FALSE;
466  }
467 
468  StpVecFree(scip, queue);
469 
470 #ifndef NDEBUG
471  for( int i = 0; i < nnodes; i++ )
472  {
473  assert(nodes_isVisited[i] == FALSE);
474  }
475 #endif
476 
477  if( compsize > *rootcompsize )
478  {
479  *rootcompsize = compsize;
480  *root = rootcand;
481  }
482 
483  SCIPfreeCleanBufferArray(scip, &nodes_isVisited);
484 }
485 
486 
487 /** finds a root terminal */
488 static
490  SCIP* scip, /**< SCIP data structure */
491  const GRAPH* g, /**< the graph */
492  const MINCUT* mincut /**< minimum cut */
493 )
494 {
495  int source = -1;
496 
497  /* NOTE: hack to allow for stable unit tests */
498  if( !mincut->randnumgen )
499  {
500  source = g->source;
501  }
502  else
503  {
504  int* terms;
505  int sourcecompsize = 0;
506  const int ntries = MIN(TERMSEPA_NROOTCANDS, g->terms);
507 
508  SCIP_CALL_ABORT( SCIPallocMemoryArray(scip, &terms, g->terms) );
509  graph_getTerms(g, terms);
510 
511  for( int i = 0; i < ntries; i++ )
512  {
513  const int pos = SCIPrandomGetInt(mincut->randnumgen, 0, g->terms - 1);
514  const int cand = terms[pos];
515 
516  updateTerminalSource(scip, cand, g, &source, &sourcecompsize );
517  }
518 
519  assert(sourcecompsize > 0);
520 
521  SCIPfreeMemoryArray(scip, &terms);
522  }
523 
524  assert(source >= 0);
525  assert(Is_term(g->term[source]));
526 
527  return source;
528 }
529 
530 
531 /** collect cut nodes */
532 static
534  const GRAPH* g, /**< the graph */
535  const TERMSEPAS* termsepas, /**< terminal separator storage */
536  const MINCUT* mincut, /**< minimum cut */
537  int sinkterm, /**< sink terminal */
538  int* cutterms, /**< terminals */
539  int* ncutterms, /**< number of terminals */
540  SCIP_Bool* cutIsGood /**< is cut nodes */
541 )
542 {
543  const int* const edges_capa = mincut->edges_capa;
544  const int* const nodes_wakeState = mincut->nodes_wakeState;
545  const int* const csr_start = mincut->csr_start;
546  const int* const csr_headarr = mincut->csr_headarr;
547  const int nnodes_extended = mincut->termsepa_nnodes;
548  const int capa_inf = termsepaGetCapaInf(g, mincut);
549  const int maxsepasize = termsepas->maxsepasize;
550  SCIP_Bool isGood = TRUE;
551  int n = 0;
552 
553  for( int i = 0; i < nnodes_extended && isGood; i++ )
554  {
555  const int start = csr_start[i];
556  const int end = csr_start[i + 1];
557 
558  if( !nodes_wakeState[i] )
559  continue;
560 
561  for( int j = start; j != end; j++ )
562  {
563  const int head = csr_headarr[j];
564 
565  if( nodes_wakeState[head] == 0 && edges_capa[j] > 0 )
566  {
567  assert(n <= maxsepasize);
568  assert(edges_capa[j] == 1 || edges_capa[j] == capa_inf);
569 
570  if( edges_capa[j] == capa_inf || n == maxsepasize )
571  {
572  isGood = FALSE;
573  break;
574  }
575 
576  cutterms[n++] = i;
577  assert(Is_term(g->term[i]));
578  }
579  }
580  }
581 
582  *cutIsGood = isGood;
583  *ncutterms = n;
584 }
585 
586 
587 
588 /** helper; traverses and does additional optional stuff */
589 static
591  SCIP* scip, /**< SCIP data structure */
592  const GRAPH* g, /**< the graph */
593  SCIP_Bool removeTerms, /**< remove terminals reachable from sink via non-terminal paths? */
594  SCIP_Bool updateVisitedTerms, /**< do update? */
595  int ncutterms, /**< size of the separator */
596  const int* cutterms, /**< the separator nodes (all terminals) */
597  int sinkterm, /**< sink terminal of current cut */
598  MINCUT* mincut, /**< minimum cut */
599  int* ncompnodes /**< number of nodes in component (OUT) */
600 )
601 {
602  int* RESTRICT termcands = mincut->terms;
603  STP_Vectype(int) queue = NULL;
604  SCIP_Bool* nodes_isVisited;
605  const int nnodes = graph_get_nNodes(g);
606 
607  assert(ncutterms > 0);
608  assert(mincut->nodes_wakeState[sinkterm] == 0);
609 
610  SCIP_CALL( SCIPallocCleanBufferArray(scip, &nodes_isVisited, nnodes) );
611 
612  for( int i = 0; i < ncutterms; i++ )
613  {
614  assert(graph_knot_isInRange(g, cutterms[i]));
615  assert(!nodes_isVisited[cutterms[i]]);
616 
617  nodes_isVisited[cutterms[i]] = TRUE;
618  }
619 
620  assert(!nodes_isVisited[mincut->root]);
621 
622  /* BFS from sink terminal */
623  StpVecReserve(scip, queue, nnodes);
624  StpVecPushBack(scip, queue, sinkterm);
625  nodes_isVisited[sinkterm] = TRUE;
626 
627  /* BFS loop */
628  for( int i = 0; i < StpVecGetSize(queue); i++ )
629  {
630  const int node = queue[i];
631  assert(nodes_isVisited[node]);
632 
633  if( removeTerms && Is_term(g->term[node]) && node != sinkterm )
634  {
635  const int ntermcands = mincut->ntermcands;
636 
637  assert(mincut->nodes_wakeState[node] == 0);
638 
639  /* todo more efficiently */
640  for( int t = 0; t < ntermcands; t++ )
641  {
642  if( node == termcands[t] )
643  {
644  // printf("removing terminal %d \n", node);
645 
646  SWAP_INTS(termcands[mincut->ntermcands - 1], termcands[t]);
647  mincut->ntermcands--;
648  break;
649  }
650  }
651 
652  continue;
653  }
654 
655  for( int e = g->outbeg[node]; e >= 0; e = g->oeat[e] )
656  {
657  const int head = g->head[e];
658  if( !nodes_isVisited[head] )
659  {
660  nodes_isVisited[head] = TRUE;
661  StpVecPushBack(scip, queue, head);
662  }
663  }
664  }
665 
666  if( ncompnodes )
667  *ncompnodes = StpVecGetSize(queue);
668 
669  for( int i = 0; i < StpVecGetSize(queue); i++ )
670  {
671  const int node = queue[i];
672  assert(nodes_isVisited[node]);
673  nodes_isVisited[node] = FALSE;
674 
675  if( updateVisitedTerms && Is_term(g->term[node]) && node != sinkterm )
676  {
677  assert(mincut->terms_minsepasize && mincut->terms_mincompsize);
678 
679  if( mincut->terms_minsepasize[node] > ncutterms )
680  mincut->terms_minsepasize[node] = ncutterms;
681 
682  if( mincut->terms_mincompsize[node] > *ncompnodes )
683  mincut->terms_mincompsize[node] = *ncompnodes;
684  }
685  }
686 
687  for( int i = 0; i < ncutterms; i++ )
688  {
689  assert(nodes_isVisited[cutterms[i]]);
690  nodes_isVisited[cutterms[i]] = FALSE;
691  }
692 
693  StpVecFree(scip, queue);
694 
695 #ifndef NDEBUG
696  for( int i = 0; i < nnodes; i++ )
697  {
698  assert(nodes_isVisited[i] == FALSE);
699  }
700 #endif
701 
702  SCIPfreeCleanBufferArray(scip, &nodes_isVisited);
703 
704  return SCIP_OKAY;
705 }
706 
707 
708 /** removes terminals in current cut */
709 static
711  SCIP* scip, /**< SCIP data structure */
712  const GRAPH* g, /**< the graph */
713  int ncutterms, /**< size of the separator */
714  const int* cutterms, /**< the separator nodes (all terminals) */
715  int sinkterm, /**< sink terminal of current cut */
716  MINCUT* mincut /**< minimum cut */
717 )
718 {
719  const SCIP_Bool removeTerms = TRUE;
720  const SCIP_Bool updateVisitedTerms = FALSE;
721 
722  SCIP_CALL( termsepaTraverseSinkComp(scip, g, removeTerms, updateVisitedTerms, ncutterms, cutterms, sinkterm, mincut, NULL) );
723 
724  return SCIP_OKAY;
725 }
726 
727 
728 /** gets number of nodes */
729 static
731  SCIP* scip, /**< SCIP data structure */
732  const GRAPH* g, /**< the graph */
733  int ncutterms, /**< size of the separator */
734  const int* cutterms, /**< the separator nodes (all terminals) */
735  int sinkterm, /**< sink terminal of current cut */
736  MINCUT* mincut, /**< minimum cut */
737  int* ncompnodes /**< number of nodes in component */
738 )
739 {
740  const SCIP_Bool removeTerms = FALSE;
741  const SCIP_Bool updateVisitedTerms = TRUE;
742 
743  SCIP_CALL( termsepaTraverseSinkComp(scip, g, removeTerms, updateVisitedTerms, ncutterms, cutterms, sinkterm, mincut, ncompnodes) );
744 
745  return SCIP_OKAY;
746 }
747 
748 
749 /** stores cut
750  * NOTE: this methods is call once the cut vertices are already stored in the CSR array */
751 static
753  SCIP* scip, /**< SCIP data structure */
754  const GRAPH* g, /**< the graph */
755  int sinkterm, /**< sink terminal */
756  MINCUT* mincut, /**< minimum cut */
757  int ncutterms, /**< number of cut terminals */
758  const int* cutterms, /**< cut terminals */
759  TERMSEPAS* termsepas, /**< terminal separator storage */
760  SCIP_Bool* success /**< success? */
761 )
762 {
763  TSEPA* sepas = termsepas->sepas;
764  int nsinknodes = 0;
765  const int nsepas_all = termsepas->nsepas_all;
766 
767  assert(0 <= nsepas_all && nsepas_all < termsepas->maxnsepas);
768  assert(0 <= ncutterms && ncutterms <= termsepas->maxsepasize);
769  assert(termsepas->nsepaterms_csr + ncutterms <= termsepas->maxnsepas * termsepas->maxsepasize);
770 
771  *success = TRUE;
772 
773  SCIP_CALL( termsepaGetCompNnodes(scip, g, ncutterms, cutterms, sinkterm, mincut, &nsinknodes) );
774 
775  /* NOTE: if the sink is already contained in a separated component and does not make anything smaller, we don't want to take it */
776  if( ncutterms >= mincut->terms_minsepasize[sinkterm] && nsinknodes >= mincut->terms_mincompsize[sinkterm] )
777  {
778  SCIPdebugMessage("discarding already covered component because of SEPA size: %d>=%d, COMP size: %d>=%d \n", ncutterms
779  , mincut->terms_minsepasize[sinkterm], nsinknodes, mincut->terms_mincompsize[sinkterm]);
780 
781  *success = FALSE;
782  return SCIP_OKAY;
783  }
784 
785  sepas[nsepas_all].sinkterm = sinkterm;
786  sepas[nsepas_all].nsinknodes = nsinknodes;
787 
788 #ifndef NDEBUG
789  {
790  const int* const nodes_wakeState = mincut->nodes_wakeState;
791  int nsinknodes_dbg = 0;
792  const int nnodes = graph_get_nNodes(g);
793  for( int i = 0; i < nnodes; i++ )
794  {
795  if( nodes_wakeState[i] == 0 )
796  nsinknodes_dbg++;
797  }
798 
799  assert(nsinknodes_dbg >= nsinknodes);
800  }
801 #endif
802 
803  termsepas->sepastarts_csr[nsepas_all + 1] = termsepas->sepastarts_csr[nsepas_all] + ncutterms;
804  termsepas->nsepaterms_csr += ncutterms;
805  termsepas->nsepas_all++;
806  termsepas->nsepas[ncutterms]++;
807 
808  assert(termsepas->sepastarts_csr[nsepas_all + 1] == termsepas->nsepaterms_csr);
809 
810 #ifdef SCIP_DEBUG
811  if( ncutterms <= TERMSEPA_MAXCUTSIZE )
812  {
813  printf("terminal cut of size %d for sink %d \n", ncutterms, sinkterm);
814  printf("nsinknodes=%d \n", nsinknodes);
815 
816  for( int i = 0; i < ncutterms; i++ )
817  {
818  printf("%d \n", cutterms[i]);
819  }
820  }
821 #endif
822 
823  return SCIP_OKAY;
824 }
825 
826 
827 /** tries to add cut */
828 static
830  SCIP* scip, /**< SCIP data structure */
831  const GRAPH* g, /**< the graph */
832  int sinkterm, /**< sink terminal */
833  MINCUT* mincut, /**< minimum cut */
834  TERMSEPAS* termsepas /**< terminal separator storage */
835 )
836 {
837  int* cutterms;
838  int ncutterms;
839  SCIP_Bool isGoodCut;
840 
841  assert(!mincut->isLpcut);
842  assert(graph_knot_isInRange(g, sinkterm) && Is_term(g->term[sinkterm]));
843  assert(mincut->nodes_wakeState[sinkterm] == 0);
844 
845  cutterms = &(termsepas->sepaterms_csr[termsepas->nsepaterms_csr]);
846  termsepaCollectCutNodes(g, termsepas, mincut, sinkterm, cutterms, &ncutterms, &isGoodCut);
847 
848  if( isGoodCut )
849  {
850  assert(termsepaCutIsCorrect(scip, g, ncutterms, cutterms, sinkterm, mincut));
851 
852  SCIP_CALL( termsepaStoreCutFinalize(scip, g, sinkterm, mincut, ncutterms, cutterms, termsepas, &isGoodCut) );
853 
854  if( isGoodCut )
855  {
856  SCIP_CALL( termsepaRemoveCutTerminals(scip, g, ncutterms, cutterms, sinkterm, mincut) );
857  }
858  }
859 
860  return SCIP_OKAY;
861 }
862 
863 /** initializes */
864 static
866  const GRAPH* g, /**< the graph */
867  MINCUT* mincut /**< minimum cut */
868 )
869 {
870  int* RESTRICT csr_edgeDefaultToCsr = mincut->csr_edgeDefaultToCsr;
871  const int nedges = graph_get_nEdges(g);
872 
873  assert(nedges >= 2);
874 
875  for( int e = 0; e < nedges; e++ )
876  csr_edgeDefaultToCsr[e] = -1;
877 
878  return SCIP_OKAY;
879 }
880 
881 
882 /** creates copies of potential separator terminals */
883 static
885  const GRAPH* g, /**< the graph */
886  MINCUT* mincut /**< minimum cut */
887 )
888 {
889  int* RESTRICT excess = g->mincut_e;
890  int* RESTRICT nodes_termToCopy = mincut->termsepa_termToCopy;
891  int* RESTRICT terms = mincut->terms;
892  const int* const nodes_wakeState = mincut->nodes_wakeState;
893  const int root = mincut->root;
894  const int nnodes_org = graph_get_nNodes(g);
895  int nsepaterms = 0;
896  int ntermcands = 0;
897 
898  assert(nodes_wakeState && nodes_termToCopy);
899 
900  for( int k = 0; k < nnodes_org; k++ )
901  {
902  SCIP_Bool isSeparator;
903  assert(nodes_termToCopy[k] == -1);
904 
905  if( !Is_term(g->term[k]) || k == root )
906  {
907  continue;
908  }
909 
910  isSeparator = FALSE;
911 
912  for( int e = g->outbeg[k]; e != EAT_LAST; e = g->oeat[e] )
913  {
914  const int head = g->head[e];
915 
916  /* head not in root cut? */
917  if( nodes_wakeState[head] == 0 )
918  {
919  isSeparator = TRUE;
920  break;
921  }
922  }
923 
924  if( nodes_wakeState[k] == 0 )
925  {
926  SCIPdebugMessage("add terminal candidate %d \n", k);
927  terms[ntermcands++] = k;
928  }
929 
930  if( isSeparator )
931  {
932  nodes_termToCopy[k] = nnodes_org + nsepaterms;
933  nsepaterms++;
934 
935  assert(excess[nodes_termToCopy[k]] == 0);
936 
937  if( nodes_wakeState[k] != 0 )
938  {
939  /* NOTE: 1 is the default value for separator edges...so we basically push everything out of k */
940  excess[nodes_termToCopy[k]] = 1;
941  }
942 
943  SCIPdebugMessage("adding separator terminal %d (copyindex=%d) \n", k, nodes_termToCopy[k]);
944 
945  }
946  }
947 
948  // todo check on large test set ... does not seem to help so far
949 #ifdef SCIP_DISABLED
950  if( mincut->randnumgen )
951  {
952  SCIPrandomPermuteIntArray(mincut->randnumgen, terms, 0, ntermcands);
953  printf("PERMUTE \n\n \n");
954  }
955 #endif
956 
957  assert(nnodes_org + nsepaterms <= mincut->termsepa_nnodes);
958 
959  mincut->termsepa_nnodes = nnodes_org + nsepaterms;
960  mincut->ntermcands = ntermcands;
961 }
962 
963 
964 /** adds CSR edges */
965 static
967  const GRAPH* g, /**< the graph */
968  MINCUT* mincut /**< minimum cut */
969 )
970 {
971  int* RESTRICT edges_capa = mincut->edges_capa;
972  int* RESTRICT residual = g->mincut_r;
973  int* RESTRICT edgecurr = g->mincut_numb;
974  int* RESTRICT csr_edgeDefaultToCsr = mincut->csr_edgeDefaultToCsr;
975  int* RESTRICT csr_start = mincut->csr_start;
976  int* RESTRICT csr_headarr = mincut->csr_headarr;
977  int* RESTRICT nodes_wakeState = mincut->nodes_wakeState;
978  const int* const nodes_termToCopy = mincut->termsepa_termToCopy;
979  const int capa_infinity = termsepaGetCapaInf(g, mincut);
980  const int capa_one = 1;
981  const int nnodes_org = graph_get_nNodes(g);
982  int csr_nedges = 0;
983 
984  assert(residual && edgecurr && edges_capa);
985 
986  /* add edges from original nodes */
987  for( int k = 0; k < nnodes_org; k++ )
988  {
989  const SCIP_Bool kIsSepaTerm = nodes_termToCopy[k] >= 0;
990  edgecurr[k] = -1;
991  csr_start[k] = csr_nedges;
992 
993  /* add edge from terminal to copy if existent */
994  if( kIsSepaTerm )
995  {
996  const int kCopy = nodes_termToCopy[k];
997 
998  assert(Is_term(g->term[k]) && k != mincut->root);
999  assert(nnodes_org <= kCopy && kCopy < mincut->termsepa_nnodes);
1000 
1001  residual[csr_nedges] = capa_one;
1002  csr_headarr[csr_nedges++] = kCopy;
1003  }
1004 
1005  /* non-dormant node or potential separator? */
1006  if( kIsSepaTerm || nodes_wakeState[k] == 0 )
1007  {
1008  assert(k != mincut->root);
1009  edgecurr[k] = csr_nedges;
1010 
1011  /* NOTE: we go two times over the incident edges so that the edges to
1012  * the terminal copies are at the start */
1013  for( int e = g->outbeg[k]; e != EAT_LAST; e = g->oeat[e] )
1014  {
1015  const int head = g->head[e];
1016 
1017  /* is head a separator terminal? */
1018  if( nodes_termToCopy[head] >= 0 )
1019  {
1020  const int headCopy = nodes_termToCopy[head];
1021 
1022  assert(Is_term(g->term[head]) && head != mincut->root);
1023  assert(nnodes_org <= headCopy && headCopy < mincut->termsepa_nnodes);
1024 
1025  residual[csr_nedges] = 0;
1026  csr_headarr[csr_nedges++] = headCopy;
1027  }
1028  }
1029 
1030  for( int e = g->outbeg[k]; e != EAT_LAST; e = g->oeat[e] )
1031  {
1032  const int head = g->head[e];
1033  const SCIP_Bool headIsSepaTerm = nodes_termToCopy[head] >= 0;
1034 
1035  if( kIsSepaTerm && headIsSepaTerm )
1036  continue;
1037 
1038  if( headIsSepaTerm || nodes_wakeState[head] == 0 )
1039  {
1040  csr_edgeDefaultToCsr[e] = csr_nedges;
1041  residual[csr_nedges] = kIsSepaTerm ? 0 : capa_infinity;
1042  csr_headarr[csr_nedges++] = head;
1043  }
1044  }
1045 
1046  /* unreachable node? */
1047  if( edgecurr[k] == csr_nedges )
1048  {
1049  assert(g->grad[k] == 0);
1050  nodes_wakeState[k] = 1;
1051  }
1052  }
1053  }
1054 
1055  /* add edges from copy terminals */
1056  for( int k = 0; k < nnodes_org; k++ )
1057  {
1058  const SCIP_Bool kIsSepaTerm = nodes_termToCopy[k] >= 0;
1059 
1060  if( kIsSepaTerm )
1061  {
1062  const int kCopy = nodes_termToCopy[k];
1063 
1064  assert(nodes_wakeState[kCopy] == 0);
1065  assert(Is_term(g->term[k]) && k != mincut->root);
1066  assert(nnodes_org <= kCopy && kCopy < mincut->termsepa_nnodes);
1067 
1068  edgecurr[kCopy] = csr_nedges;
1069  csr_start[kCopy] = csr_nedges;
1070 
1071  /* edge from copy to k */
1072  residual[csr_nedges] = 0;
1073  csr_headarr[csr_nedges++] = k;
1074 
1075  // should not be necessary todo remove
1076 #ifdef SCIP_DISABLED
1077  for( int e = g->outbeg[k]; e != EAT_LAST; e = g->oeat[e] )
1078  {
1079  const int head = g->head[e];
1080  /* is head a separator terminal? */
1081  if( nodes_termToCopy[head] >= 0 )
1082  {
1083  const int headCopy = nodes_termToCopy[head];
1084 
1085  assert(Is_term(g->term[head]) && head != root);
1086  assert(nnodes_org <= headCopy && headCopy < mincut->termsepa_nnodes);
1087 
1088  residual[csr_nedges] = 0;
1089  csr_headarr[csr_nedges++] = headCopy;
1090  }
1091  }
1092 #endif
1093 
1094  for( int e = g->outbeg[k]; e != EAT_LAST; e = g->oeat[e] )
1095  {
1096  const int head = g->head[e];
1097  const SCIP_Bool headIsSepaTerm = nodes_termToCopy[head] >= 0;
1098 
1099  if( headIsSepaTerm || nodes_wakeState[head] == 0 )
1100  {
1101  residual[csr_nedges] = capa_infinity;
1102  csr_headarr[csr_nedges++] = head;
1103  }
1104  }
1105  }
1106  }
1107 
1108  assert(csr_nedges <= mincut->termsepa_nedges);
1109 
1110  mincut->termsepa_nedges = csr_nedges;
1111  csr_start[mincut->termsepa_nnodes] = csr_nedges;
1112 
1113  BMScopyMemoryArray(edges_capa, residual, csr_nedges);
1114 
1115  mincut->csr_nedges = csr_nedges;
1116 }
1117 
1118 
1119 /** adds CSR reverse edges */
1120 static
1122  const GRAPH* g, /**< the graph */
1123  MINCUT* mincut /**< minimum cut */
1124 )
1125 {
1126  int* RESTRICT csr_edgeflipped = mincut->csr_edgeflipped;
1127  const int* const csr_edgeDefaultToCsr = mincut->csr_edgeDefaultToCsr;
1128  const int* const csr_start = mincut->csr_start;
1129  const int* const csr_headarr = mincut->csr_headarr;
1130  const int nnodes_org = graph_get_nNodes(g);
1131  const int nedges_org = graph_get_nEdges(g);
1132  const int nnodes_extended = mincut->termsepa_nnodes;
1133 
1134  /* go over all terminal copies */
1135  for( int copyterm = nnodes_org; copyterm < nnodes_extended; copyterm++ )
1136  {
1137  const int copyedges_start = csr_start[copyterm];
1138  const int copyedges_end = csr_start[copyterm + 1];
1139 
1140  assert(copyedges_start < copyedges_end);
1141 
1142  for( int copyedge = copyedges_start; copyedge != copyedges_end; copyedge++ )
1143  {
1144  int antiedge;
1145  const int copyneighbor = csr_headarr[copyedge];
1146  const int neighboredges_start = csr_start[copyneighbor];
1147  const int neighboredges_end = csr_start[copyneighbor + 1];
1148 
1149  assert(neighboredges_start < neighboredges_end);
1150 
1151  /* NOTE: copy-edges are at the beginning, so should be half-way efficient */
1152  for( antiedge = neighboredges_start; antiedge != neighboredges_end; antiedge++ )
1153  {
1154  if( csr_headarr[antiedge] == copyterm )
1155  {
1156  break;
1157  }
1158  }
1159 
1160  assert(antiedge < neighboredges_end);
1161  assert(csr_edgeflipped[antiedge] == -1 || csr_edgeflipped[antiedge] == copyedge);
1162  assert(csr_edgeflipped[copyedge] == -1 || csr_edgeflipped[copyedge] == antiedge);
1163 
1164  csr_edgeflipped[antiedge] = copyedge;
1165  csr_edgeflipped[copyedge] = antiedge;
1166  }
1167  }
1168 
1169  for( int e = 0; e < nedges_org; e++ )
1170  {
1171  if( csr_edgeDefaultToCsr[e] >= 0 )
1172  {
1173  const int csr_pos = csr_edgeDefaultToCsr[e];
1174  assert(csr_edgeflipped[csr_pos] == -1);
1175  csr_edgeflipped[csr_pos] = csr_edgeDefaultToCsr[flipedge(e)];
1176  }
1177  }
1178 
1179  assert(csrFlipedgesAreValid(g, mincut));
1180 
1181 }
1182 
1183 
1184 /** gets maximum number of nodes for extended terminal separation graph */
1185 static
1187  const GRAPH* g /**< the graph */
1188 )
1189 {
1190  const int nnodes = graph_get_nNodes(g);
1191  const int nterms = graph_get_nTerms(g);
1192 
1193  assert(nterms >= 2);
1194 
1195  return nnodes + nterms - 1;
1196 }
1197 
1198 
1199 /** gets maximum number of edges for extended terminal separation graph */
1200 static
1202  int root, /**< root of enlarged graph */
1203  const GRAPH* g /**< the graph */
1204 )
1205 {
1206  const int nnodes = graph_get_nNodes(g);
1207  const int nedges = graph_get_nEdges(g);
1208  int nedges_enlarged = nedges;
1209 
1210  assert(graph_knot_isInRange(g, root));
1211  assert(Is_term(g->term[root]));
1212  assert(nedges >= 2);
1213 
1214  for( int i = 0; i < nnodes; i++ )
1215  {
1216  if( !Is_term(g->term[i]) || i == root )
1217  continue;
1218 
1219  /* NOTE: the + 2 is for the two arcs between the terminal and its copy */
1220  nedges_enlarged += 2 * g->grad[i] + 2;
1221  }
1222 
1223  return nedges_enlarged;
1224 }
1225 
1226 
1227 /** builds and marks root component (nodes reachable from root via non-terminal paths) */
1228 static
1230  const GRAPH* g, /**< the graph */
1231  MINCUT* mincut /**< minimum cut */
1232 )
1233 {
1234  int* RESTRICT nodes_wakeState = mincut->nodes_wakeState;
1235  int* RESTRICT rootcut = mincut->rootcut;
1236  const int root = mincut->root;
1237  int rootcutsize = 0;
1238 
1239  nodes_wakeState[root] = 1;
1240  rootcut[rootcutsize++] = root;
1241 
1242  /* BFS loop, ignoring terminals */
1243  for( int i = 0; i < rootcutsize; i++ )
1244  {
1245  const int k = rootcut[i];
1246 
1247  assert(rootcutsize <= g->knots);
1248  assert(k < g->knots);
1249 
1250  if( Is_term(g->term[k]) && k != root )
1251  continue;
1252 
1253  /* traverse outgoing arcs */
1254  for( int e = g->outbeg[k]; e != EAT_LAST; e = g->oeat[e] )
1255  {
1256  const int head = g->head[e];
1257 
1258  /* head not been added to root cut yet? */
1259  if( nodes_wakeState[head] == 0 )
1260  {
1261  nodes_wakeState[head] = 1;
1262  SCIPdebugMessage("add to root cut: %d \n", head);
1263 
1264  rootcut[rootcutsize++] = head;
1265  }
1266  }
1267  }
1268 
1269  mincut->rootcutsize = rootcutsize;
1270 }
1271 
1272 
1273 /** builds CSR representation of enlarged graph; also build terminal candidates */
1274 static
1276  const GRAPH* g, /**< the graph */
1277  MINCUT* mincut /**< minimum cut */
1278 )
1279 {
1280  SCIP_CALL( termsepaCsrInit(g, mincut) );
1281 
1282  /* create copies for potential separator terminals */
1283  termsepaCsrAddTermCopies(g, mincut);
1284 
1285  termsepaCsrAddEdges(g, mincut);
1286 
1287  termsepaCsrAddReverseEdges(g, mincut);
1288 
1289 #ifdef SCIP_DEBUG
1290  //debugPrintCsr(g, mincut);
1291 #endif
1292 
1293  return SCIP_OKAY;
1294 }
1295 
1296 
1297 /** initializes for LP */
1298 static
1300  SCIP* scip, /**< SCIP data structure */
1301  const GRAPH* g, /**< the graph */
1302  MINCUT* mincut /**< minimum cut */
1303 )
1304 {
1305  SCIP_VAR** vars = SCIPprobdataGetVars(scip);
1306  int* nodes_wakeState;
1307  STP_Bool* RESTRICT edges_isRemoved;
1308  const int nedges = graph_get_nEdges(g);
1309  const int nnodes = graph_get_nNodes(g);
1310 
1311  assert(mincut && scip);
1312  assert(mincut->isLpcut);
1313  assert(!mincut->edges_isRemoved);
1314 
1315  mincut->xval = SCIPprobdataGetXval(scip, NULL);
1316  assert(mincut->xval);
1317 
1318  SCIP_CALL( SCIPallocBufferArray(scip, &(mincut->edges_capa), nedges) );
1319  SCIP_CALL( SCIPallocBufferArray(scip, &(mincut->nodes_wakeState), nnodes) );
1320  SCIP_CALL( SCIPallocBufferArray(scip, &(mincut->terms), g->terms) );
1321  SCIP_CALL( SCIPallocBufferArray(scip, &(mincut->csr_edgeDefaultToCsr), nedges) );
1322  SCIP_CALL( SCIPallocBufferArray(scip, &(mincut->csr_headarr), nedges) );
1323  SCIP_CALL( SCIPallocBufferArray(scip, &(mincut->csr_edgeflipped), nedges) );
1324  SCIP_CALL( SCIPallocBufferArray(scip, &(mincut->csr_start), nnodes + 1) );
1325  SCIP_CALL( SCIPallocBufferArray(scip, &(mincut->rootcut), nnodes + 1) );
1326  SCIP_CALL( SCIPallocBufferArray(scip, &(mincut->edges_isRemoved), nedges) );
1327 
1328  edges_isRemoved = mincut->edges_isRemoved;
1329  for( int i = 0; i < nedges; i++ )
1330  {
1331  edges_isRemoved[i] = (SCIPvarGetUbGlobal(vars[i]) < 0.5);
1332  }
1333 
1334  nodes_wakeState = mincut->nodes_wakeState;
1335  for( int k = 0; k < nnodes; k++ )
1336  {
1337  nodes_wakeState[k] = 0;
1338  }
1339 
1340 #ifndef NDEBUG
1341  for( int i = 0; i < nedges; i++ )
1342  {
1343  mincut->csr_edgeflipped[i] = -1;
1344  }
1345 #endif
1346 
1347  return SCIP_OKAY;
1348 }
1349 
1350 
1351 /** initializes for terminal separation */
1352 static
1354  SCIP* scip, /**< SCIP data structure */
1355  GRAPH* g, /**< the graph */
1356  MINCUT* mincut /**< minimum cut */
1357 )
1358 {
1359  int* RESTRICT nodes_termToCopy;
1360  int* RESTRICT nodes_wakeState;
1361  int* RESTRICT terms_sepasize;
1362  int* RESTRICT terms_mincompsize;
1363  int nnodes_enlarged;
1364  int nedges_enlarged;
1365  const int nedges = graph_get_nEdges(g);
1366  const int nnodes = graph_get_nNodes(g);
1367 
1368  assert(mincut && scip);
1369  assert(!mincut->isLpcut);
1370  assert(!mincut->edges_isRemoved);
1371 
1372  mincut->root = termsepaFindTerminalSource(scip, g, mincut);
1373  nnodes_enlarged = termsepaCsrGetMaxNnodes(g);
1374  nedges_enlarged = termsepaCsrGetMaxNedges(mincut->root, g);
1375 
1376  SCIPdebugMessage("selected source %d \n", mincut->root);
1377 
1378  SCIP_CALL( SCIPallocBufferArray(scip, &(mincut->edges_capa), nedges_enlarged) );
1379  SCIP_CALL( SCIPallocBufferArray(scip, &(mincut->nodes_wakeState), nnodes_enlarged) );
1380  SCIP_CALL( SCIPallocBufferArray(scip, &(mincut->terms), g->terms) );
1381  SCIP_CALL( SCIPallocBufferArray(scip, &(mincut->csr_edgeDefaultToCsr), nedges) );
1382  SCIP_CALL( SCIPallocBufferArray(scip, &(mincut->csr_headarr), nedges_enlarged) );
1383  SCIP_CALL( SCIPallocBufferArray(scip, &(mincut->csr_edgeflipped), nedges_enlarged) );
1384  SCIP_CALL( SCIPallocBufferArray(scip, &(mincut->csr_start), nnodes_enlarged + 1) );
1385  SCIP_CALL( SCIPallocBufferArray(scip, &(mincut->rootcut), nnodes + 1) );
1386  SCIP_CALL( SCIPallocBufferArray(scip, &(mincut->termsepa_termToCopy), nnodes) );
1387  SCIP_CALL( SCIPallocBufferArray(scip, &(mincut->terms_minsepasize), nnodes) );
1388  SCIP_CALL( SCIPallocBufferArray(scip, &(mincut->terms_mincompsize), nnodes) );
1389 
1390  nodes_termToCopy = mincut->termsepa_termToCopy;
1391  terms_sepasize = mincut->terms_minsepasize;
1392  terms_mincompsize = mincut->terms_mincompsize;
1393  for( int i = 0; i < nnodes; i++ )
1394  {
1395  nodes_termToCopy[i] = -1;
1396  terms_sepasize[i] = nnodes;
1397  terms_mincompsize[i] = nnodes;
1398  }
1399 
1400  mincut->termsepa_nnodes = nnodes_enlarged;
1401  mincut->termsepa_nedges = nedges_enlarged;
1402 
1403  SCIP_CALL( graph_mincut_reInit(scip, nnodes_enlarged, nedges_enlarged, g) );
1404 
1405  nodes_wakeState = mincut->nodes_wakeState;
1406  for( int k = 0; k < nnodes_enlarged; k++ )
1407  {
1408  nodes_wakeState[k] = 0;
1409  }
1410 
1411 #ifndef NDEBUG
1412  for( int i = 0; i < nedges_enlarged; i++ )
1413  {
1414  mincut->csr_edgeflipped[i] = -1;
1415  mincut->edges_capa[i] = -1;
1416  }
1417 #endif
1418 
1419  return SCIP_OKAY;
1420 }
1421 
1422 
1423 /** initializes */
1424 static
1426  SCIP* scip, /**< SCIP data structure */
1427  SCIP_RANDNUMGEN* randnumgen, /**< random number generator or NULL */
1428  SCIP_Bool isLpcut, /**< for LP cut? */
1429  GRAPH* g, /**< the graph */
1430  MINCUT** mincut /**< minimum cut */
1431 )
1432 {
1433  MINCUT* mcut;
1434 
1435  assert(scip);
1436  assert(graph_get_nNodes(g) > 0 && graph_get_nEdges(g) > 0);
1437 
1438  SCIP_CALL( SCIPallocMemory(scip, mincut) );
1439  mcut = *mincut;
1440 
1441  mcut->xval = NULL;
1442  mcut->edges_isRemoved = NULL;
1443  mcut->isLpcut = FALSE;
1444  mcut->ntermcands = -1;
1445  mcut->rootcutsize = -1;
1446  mcut->csr_nedges = -1;
1447  mcut->root = g->source;
1448  mcut->termsepa_termToCopy = NULL;
1449  mcut->terms_minsepasize = NULL;
1450  mcut->terms_mincompsize = NULL;
1451  mcut->termsepa_nnodes = -1;
1452  mcut->termsepa_nedges = -1;
1453  mcut->randnumgen = randnumgen;
1454  // mcut->termsepa_inEdges = NULL;
1455  // mcut->termsepa_inNeighbors = NULL;
1456  // mcut->termsepa_inEdgesStart = NULL;
1457 
1458  if( isLpcut )
1459  {
1460  mcut->isLpcut = TRUE;
1461  SCIP_CALL( mincutInitForLp(scip, g, mcut) );
1462  }
1463  else
1464  {
1465  mcut->isLpcut = FALSE;
1466  SCIP_CALL( mincutInitForTermSepa(scip, g, mcut) );
1467  }
1468 
1469 
1470  return SCIP_OKAY;
1471 }
1472 
1473 
1474 /** prepares for LP */
1475 static
1477  SCIP* scip, /**< SCIP data structure */
1478  const GRAPH* g, /**< the graph */
1479  MINCUT* mincut /**< minimum cut */
1480 )
1481 {
1482  SCIP_VAR** vars = SCIPprobdataGetVars(scip);
1483  int* RESTRICT nodes_wakeState = mincut->nodes_wakeState;
1484  int* RESTRICT edges_capa = mincut->edges_capa;
1485  int* RESTRICT terms = mincut->terms;
1486  int* RESTRICT csr_start = mincut->csr_start;
1487  int* RESTRICT rootcut = mincut->rootcut;
1488  int* RESTRICT csr_edgeDefaultToCsr = mincut->csr_edgeDefaultToCsr;
1489  int* RESTRICT csr_headarr = mincut->csr_headarr;
1490  int* RESTRICT csr_edgeflipped = mincut->csr_edgeflipped;
1491  const STP_Bool* const edges_isRemoved = mincut->edges_isRemoved;
1492  const SCIP_Real* const xval = mincut->xval;
1493  const int nnodes = graph_get_nNodes(g);
1494  const int nedges = graph_get_nEdges(g);
1495  const int root = g->source;
1496  int rootcutsize;
1497  int csr_nedges;
1498  int ntermcands;
1499  const SCIP_Bool intree = (SCIPgetDepth(scip) > 0);
1500  int* RESTRICT excess = g->mincut_e;
1501  int* RESTRICT residual = g->mincut_r;
1502  int* RESTRICT edgecurr = g->mincut_numb;
1503 
1504  assert(residual && edgecurr);
1505  assert(xval && vars && edges_isRemoved);
1506  assert(mincut->isLpcut);
1507 
1508 #ifdef SCIP_DISABLED
1509  for( int e = 0; e < nedges; e += 2 )
1510  {
1511  const int erev = e + 1;
1512  SCIP_Bool eIsRemoved;
1513 
1514  if( intree )
1515  eIsRemoved = SCIPvarGetUbLocal(vars[e]) < 0.5 && SCIPvarGetUbLocal(vars[erev]) < 0.5;
1516  else
1517  eIsRemoved = edges_isRemoved[e] && edges_isRemoved[erev];
1518 
1519  assert(eIsRemoved || !edges_isRemoved[e] || !edges_isRemoved[erev]);
1520 
1521  if( eIsRemoved )
1522  {
1523  edges_capa[e] = 0;
1524  edges_capa[erev] = 0;
1525  // residual[e] = 0;
1526  // residual[erev] = 0;
1527 
1528  csr_headarr[e] = 1;
1529  csr_headarr[erev] = 1;
1530  }
1531  else
1532  {
1533  edges_capa[e] = (int)(xval[e] * FLOW_FACTOR + 0.5) + CREEP_VALUE;
1534  edges_capa[erev] = (int)(xval[erev] * FLOW_FACTOR + 0.5) + CREEP_VALUE;
1535  // residual[e] = edges_capa[e];
1536  // residual[erev] = edges_capa[erev];
1537 
1538  /* NOTE: here we misuse csr_headarr to mark the free edges for the BFS */
1539  csr_headarr[e] = SCIPisFeasGE(scip, xval[e], 1.0) ? 0 : 1;
1540  csr_headarr[erev] = SCIPisFeasGE(scip, xval[erev], 1.0) ? 0 : 1;
1541  }
1542  }
1543 #endif
1544 
1545  for( int e = 0; e < nedges; e++ )
1546  {
1547  SCIP_Bool eIsRemoved = edges_isRemoved[e];
1548 
1549  if( intree && !eIsRemoved )
1550  {
1551  // todo remove...does not really work together with non-local cuts!
1552  eIsRemoved = (SCIPvarGetUbLocal(vars[e]) < 0.5 && SCIPvarGetUbLocal(vars[flipedge(e)]) < 0.5);
1553  }
1554 
1555  if( eIsRemoved )
1556  {
1557  edges_capa[e] = 0;
1558  csr_headarr[e] = 1;
1559  }
1560  else
1561  {
1562  edges_capa[e] = (int)(xval[e] * FLOW_FACTOR + 0.5) + CREEP_VALUE;
1563 
1564  /* NOTE: here we misuse csr_headarr to mark the free edges for the BFS */
1565  csr_headarr[e] = SCIPisFeasGE(scip, xval[e], 1.0) ? 0 : 1;
1566  }
1567  }
1568 
1569  /*
1570  * BFS along 0 edges from the root
1571  * */
1572 
1573  nodes_wakeState[root] = 1;
1574  rootcutsize = 0;
1575  rootcut[rootcutsize++] = root;
1576 
1577  /* BFS loop */
1578  for( int i = 0; i < rootcutsize; i++ )
1579  {
1580  const int k = rootcut[i];
1581 
1582  assert(rootcutsize <= nnodes);
1583  assert(k < nnodes);
1584 
1585  /* traverse outgoing arcs */
1586  for( int e = g->outbeg[k]; e != EAT_LAST; e = g->oeat[e] )
1587  {
1588  const int head = g->head[e];
1589 
1590  /* head not been added to root cut yet? */
1591  if( nodes_wakeState[head] == 0 )
1592  {
1593  if( csr_headarr[e] == 0 )
1594  {
1595  nodes_wakeState[head] = 1;
1596  rootcut[rootcutsize++] = head;
1597  }
1598  else
1599  {
1600  /* push as much as possible out of perpetually dormant nodes (possibly to other dormant nodes) */
1601  assert(nodes_wakeState[head] == 0);
1602  excess[head] += edges_capa[e];
1603  }
1604  }
1605  }
1606  }
1607 
1608  for( int e = 0; e < nedges; e++ )
1609  csr_edgeDefaultToCsr[e] = -1;
1610 
1611  csr_nedges = 0;
1612  ntermcands = 0;
1613 
1614  /* fill auxiliary adjacent vertex/edges arrays and get useable terms */
1615  for( int k = 0; k < nnodes; k++ )
1616  {
1617  csr_start[k] = csr_nedges;
1618 
1619  /* non-dormant node? */
1620  if( nodes_wakeState[k] == 0 )
1621  {
1622  edgecurr[k] = csr_nedges;
1623  for( int e = g->outbeg[k]; e != EAT_LAST; e = g->oeat[e] )
1624  {
1625  const int head = g->head[e];
1626  if( nodes_wakeState[head] == 0 )
1627  {
1628  if( edges_capa[e] == 0 && edges_capa[flipedge(e)] == 0 )
1629  continue;
1630 
1631  csr_edgeDefaultToCsr[e] = csr_nedges;
1632  residual[csr_nedges] = edges_capa[e];
1633  csr_headarr[csr_nedges++] = head;
1634  }
1635  }
1636 
1637  /* unreachable node? */
1638  if( edgecurr[k] == csr_nedges )
1639  {
1640  nodes_wakeState[k] = 1;
1641  }
1642  else if( Is_term(g->term[k]) )
1643  {
1644  terms[ntermcands++] = k;
1645  }
1646  }
1647  else
1648  {
1649  edgecurr[k] = -1;
1650  }
1651  }
1652 
1653  csr_start[nnodes] = csr_nedges;
1654 
1655  mincut->ntermcands = ntermcands;
1656  mincut->rootcutsize = rootcutsize;
1657  mincut->csr_nedges = csr_nedges;
1658 
1659  /* initialize edgeflipped */
1660  for( int e = 0; e < nedges; e++ )
1661  {
1662  if( csr_edgeDefaultToCsr[e] >= 0 )
1663  {
1664  const int csr_pos = csr_edgeDefaultToCsr[e];
1665  csr_edgeflipped[csr_pos] = csr_edgeDefaultToCsr[flipedge(e)];
1666  }
1667  }
1668 
1669  SCIPdebugMessage("Cut Pretest: %d eliminations\n", g->terms - ntermcands - 1);
1670 
1671  return SCIP_OKAY;
1672 }
1673 
1674 
1675 
1676 /** prepares for terminal separator comptation */
1677 static
1679  SCIP* scip, /**< SCIP data structure */
1680  GRAPH* g, /**< the graph */
1681  MINCUT* mincut /**< minimum cut */
1682 )
1683 {
1684  assert(!mincut->isLpcut);
1685 
1686  termsepaBuildRootcomp(g, mincut);
1687 
1688  /* build enlarged graph (and also terminal cut candidates) */
1689  SCIP_CALL( termsepaBuildCsr(g, mincut) );
1690 
1691  return SCIP_OKAY;
1692 }
1693 
1694 
1695 /** returns next promising terminal for computing cut */
1696 static
1698  const GRAPH* g, /**< graph data structure */
1699  SCIP_Bool firstrun, /**< first run? */
1700  MINCUT* mincut /**< minimum cut */
1701  )
1702 {
1703  int* RESTRICT termcands = mincut->terms;
1704  const int* const w = mincut->nodes_wakeState;
1705  int k;
1706  int t;
1707  int mindist = g->knots + 1;
1708  const int ntermcands = mincut->ntermcands;
1709 
1710  assert(termcands != NULL);
1711  assert(ntermcands > 0);
1712 
1713  if( firstrun )
1714  {
1715  if( mincut->randnumgen && ntermcands > 1 )
1716  {
1717  const int pos = SCIPrandomGetInt(mincut->randnumgen, 0, ntermcands - 1);
1718  assert(0 <= pos && pos <= ntermcands - 1);
1719 
1720  SWAP_INTS(termcands[ntermcands - 1], termcands[pos]);
1721  }
1722 
1723  assert(w[termcands[ntermcands - 1]] == 0);
1724  return termcands[ntermcands - 1];
1725  }
1726 
1727  k = -1;
1728 
1729  for( int i = 0; i < ntermcands; i++ )
1730  {
1731  assert(w[termcands[i]] >= 0);
1732 
1733  if( w[termcands[i]] == 0 )
1734  {
1735  assert(g->mincut_dist[termcands[i]] < g->knots + 1);
1736 
1737  if( g->mincut_dist[termcands[i]] < mindist )
1738  {
1739  k = i;
1740  mindist = g->mincut_dist[termcands[i]];
1741  }
1742  }
1743  }
1744 
1745  if( k == -1 )
1746  {
1747  int wmax = 0;
1748 
1749  for( int i = 0; i < ntermcands; i++ )
1750  {
1751  if( w[termcands[i]] > wmax )
1752  {
1753  k = i;
1754  wmax = w[termcands[i]];
1755  mindist = g->mincut_dist[termcands[i]];
1756  }
1757  else if( w[termcands[i]] == wmax && g->mincut_dist[termcands[i]] < mindist )
1758  {
1759  assert(wmax != 0);
1760 
1761  k = i;
1762  mindist = g->mincut_dist[termcands[i]];
1763  }
1764  }
1765  }
1766 
1767  assert(k >= 0);
1768  assert(k < ntermcands);
1769 
1770  t = termcands[k];
1771  termcands[k] = termcands[ntermcands - 1];
1772 
1773  return t;
1774 }
1775 
1776 
1777 /** executes */
1778 static
1780  const GRAPH* g, /**< the graph */
1781  int sinkterm, /**< sink terminal */
1782  SCIP_Bool wasRerun, /**< not the first run? */
1783  MINCUT* mincut /**< minimum cut */
1784 )
1785 {
1786  const int* capa;
1787  int nnodes;
1788  const int root = mincut->root;
1789 
1790  if( mincut->isLpcut )
1791  {
1792  capa = mincut->edges_capa;
1793  nnodes = graph_get_nNodes(g);
1794  assert(g->source == root);
1795  }
1796  else
1797  {
1798  /* NOTE: in this case mincut->edges_capa refers to the extended CSR graph,
1799  * and cannot be handled by graph_mincut_exec....anyway only needed for debug checks*/
1800  capa = NULL;
1801  nnodes = mincut->termsepa_nnodes;
1802  }
1803 
1804  graph_mincut_exec(g, root, sinkterm, nnodes, mincut->csr_nedges, mincut->rootcutsize,
1805  mincut->rootcut, capa, mincut->nodes_wakeState,
1806  mincut->csr_start, mincut->csr_edgeflipped, mincut->csr_headarr, wasRerun);
1807 }
1808 
1809 
1810 /** frees */
1811 static
1813  SCIP* scip, /**< SCIP data structure */
1814  MINCUT** mincut /**< minimum cut */
1815  )
1816 {
1817  MINCUT* mcut;
1818 
1819  assert(scip && mincut);
1820  assert(*mincut);
1821  mcut = *mincut;
1822 
1826  SCIPfreeBufferArrayNull(scip, &(mcut->edges_isRemoved));
1827  SCIPfreeBufferArray(scip, &(mcut->rootcut));
1828  SCIPfreeBufferArray(scip, &(mcut->csr_start));
1829  SCIPfreeBufferArray(scip, &(mcut->csr_edgeflipped));
1830  SCIPfreeBufferArray(scip, &(mcut->csr_headarr));
1831  SCIPfreeBufferArray(scip, &(mcut->csr_edgeDefaultToCsr));
1832  SCIPfreeBufferArray(scip, &(mcut->terms));
1833  SCIPfreeBufferArray(scip, &(mcut->nodes_wakeState));
1834  SCIPfreeBufferArrayNull(scip, &(mcut->edges_capa));
1835 
1836  SCIPfreeMemory(scip, mincut);
1837 }
1838 
1839 
1840 #ifdef SCIP_DISABLED
1841 /** checks */
1842 static
1843 SCIP_RETCODE lpcutRemoveReachableTerms(
1844  SCIP* scip, /**< SCIP data structure */
1845  const GRAPH* g, /**< the graph */
1846  int sinkterm,
1847  MINCUT* mincut /**< minimum cut */
1848 )
1849 {
1850  SCIP_Bool* nodes_isVisited;
1851  STP_Vectype(int) queue = NULL;
1852  const int* const nodes_wakeState = mincut->nodes_wakeState;
1853  const SCIP_Real* const xval = mincut->xval;
1854  const int nnodes = graph_get_nNodes(g);
1855 
1856  assert(xval);
1857  assert(mincut->isLpcut);
1858 
1859  SCIP_CALL( SCIPallocBufferArray(scip, &nodes_isVisited, nnodes) );
1860 
1861  for( int i = 0; i < nnodes; i++ )
1862  nodes_isVisited[i] = FALSE;
1863 
1864  printf("check for removables from %d \n", sinkterm);
1865 
1866  /* BFS from sink terminal */
1867  StpVecReserve(scip, queue, 8);
1868  StpVecPushBack(scip, queue, sinkterm);
1869  nodes_isVisited[sinkterm] = TRUE;
1870 
1871  /* BFS loop */
1872  for( int i = 0; i < StpVecGetSize(queue); i++ )
1873  {
1874  const int node = queue[i];
1875 
1876  assert(StpVecGetSize(queue) <= nnodes);
1877  assert(graph_knot_isInRange(g, node));
1878 
1879  /* traverse outgoing arcs */
1880  for( int e = g->outbeg[node]; e >= 0; e = g->oeat[e] )
1881  {
1882  const int head = g->head[e];
1883 
1884  if( nodes_wakeState[head] == 0 && !nodes_isVisited[head] )
1885  {
1886  if( SCIPisFeasGE(scip, xval[e], 1.0) )
1887  {
1888  nodes_isVisited[head] = TRUE;
1889  StpVecPushBack(scip, queue, head);
1890 
1891  if( Is_term(g->term[head]) )
1892  {
1893  printf("found terminal %d \n", head);
1894 
1895  const int ntermcands = mincut->ntermcands;
1896 
1897  assert(mincut->nodes_wakeState[node] == 0);
1898 
1899  /* todo more efficiently */
1900  for( int t = 0; t < ntermcands; t++ )
1901  {
1902  if( head == mincut->terms[t] )
1903  {
1904  printf("removing terminal %d \n", head);
1905 
1906  SWAP_INTS(mincut->terms[mincut->ntermcands - 1], mincut->terms[t]);
1907  mincut->ntermcands--;
1908  break;
1909  }
1910  }
1911 
1912  }
1913  }
1914  }
1915  }
1916  }
1917 
1918  StpVecFree(scip, queue);
1919 
1920  SCIPfreeBufferArray(scip, &nodes_isVisited);
1921 
1922  return SCIP_OKAY;
1923 }
1924 #endif
1925 
1926 
1927 /** add a cut */
1928 static
1930  SCIP* scip, /**< SCIP data structure */
1931  SCIP_CONSHDLR* conshdlr, /**< constraint handler */
1932  const GRAPH* g, /**< graph data structure */
1933  const int* nodes_inRootComp, /**< for each node: in root component? */
1934  const STP_Bool* edges_isRemoved, /**< for each edge: removed? */
1935  const SCIP_Real* xvals, /**< edge values */
1936  int* capa, /**< edges capacities (scaled) */
1937  const int updatecapa, /**< update capacities? */
1938  SCIP_Bool local, /**< is the cut local? */
1939  SCIP_Bool* success /**< pointer to store whether add cut be added */
1940  )
1941 {
1942  SCIP_ROW* row;
1943  SCIP_VAR** vars = SCIPprobdataGetVars(scip);
1944  SCIP_Real sum = 0.0;
1945  const int nedges = graph_get_nEdges(g);
1946  const int* const gtail = g->tail;
1947  const int* const ghead = g->head;
1948 
1949  assert(g);
1950  assert(g->knots > 0);
1951  assert(xvals);
1952  assert(!updatecapa);
1953 
1954  (*success) = FALSE;
1955 
1956  SCIP_CALL( SCIPcreateEmptyRowConshdlr(scip, &row, conshdlr, "twocut", 1.0, SCIPinfinity(scip), local, FALSE, TRUE) );
1957  SCIP_CALL( SCIPcacheRowExtensions(scip, row) );
1958 
1959  assert(nodes_inRootComp[g->source]);
1960 
1961  for( int i = 0; i < nedges; i++ )
1962  {
1963  if( !nodes_inRootComp[ghead[i]] && nodes_inRootComp[gtail[i]] )
1964  {
1965 #ifdef STP_USE_ADVANCED_FLOW
1966  if( updatecapa )
1967  {
1968  const SCIP_Bool inccapa = (capa[e] < FLOW_FACTOR);
1969 
1970  capa[e] = FLOW_FACTOR;
1971 
1972  if( !inccapa )
1973  {
1974  SCIP_CALL(SCIPflushRowExtensions(scip, row));
1975  SCIP_CALL(SCIPreleaseRow(scip, &row));
1976  return SCIP_OKAY;
1977  }
1978  }
1979 #endif
1980 
1981  if( edges_isRemoved[i] )
1982  {
1983  assert(EQ(xvals[i], 0.0));
1984  continue;
1985  }
1986 
1987  sum += xvals[i];
1988 
1989  if( SCIPisFeasGE(scip, sum, 1.0) )
1990  {
1991  SCIP_CALL( SCIPflushRowExtensions(scip, row) );
1992  SCIP_CALL( SCIPreleaseRow(scip, &row) );
1993  // printf("Discard cut! %.18f >= 1.0 \n", sum);
1994 
1995  return SCIP_OKAY;
1996  }
1997 
1998  SCIP_CALL( SCIPaddVarToRow(scip, row, vars[i], 1.0) );
1999  }
2000  }
2001 
2002  SCIP_CALL( SCIPflushRowExtensions(scip, row) );
2003 
2004  /* checks whether cut is sufficiently violated */
2005  if( SCIPisCutEfficacious(scip, NULL, row) )
2006  {
2007  SCIP_Bool infeasible;
2008 
2009  SCIPdebug( SCIP_CALL( SCIPprintRow(scip, row, NULL) ) );
2010 
2011  SCIP_CALL( SCIPaddRow(scip, row, FALSE, &infeasible) );
2012  assert(!infeasible);
2013 
2014 #if ADDCUTSTOPOOL
2015  /* if at root node, add cut to pool */
2016  if( !infeasible )
2017  SCIP_CALL( SCIPaddPoolCut(scip, row) );
2018 #endif
2019  (*success) = TRUE;
2020  }
2021 
2022  SCIP_CALL( SCIPreleaseRow(scip, &row) );
2023 
2024  return SCIP_OKAY;
2025 }
2026 
2027 
2028 /** sets (integer) capacity for edges*/
2029 static
2031  const GRAPH* g, /**< graph data structure */
2032  const SCIP_Real* xval, /**< edge values */
2033  SCIP_Bool creep_flow, /**< creep flow? */
2034  SCIP_Bool flip, /**< reverse the flow? */
2035  int* RESTRICT capa /**< edges capacities (scaled) */
2036  )
2037 {
2038  int k;
2039  int krev;
2040  int nedges = g->edges;
2041 
2042  assert(0 && "currently not supported!");
2043  assert(g != NULL);
2044  assert(xval != NULL);
2045 
2046  for( k = 0; k < nedges; k += 2 )
2047  {
2048  krev = k + 1;
2049  if( !flip )
2050  {
2051  capa[k] = (int) (xval[k] * FLOW_FACTOR + 0.5);
2052  capa[krev] = (int) (xval[krev] * FLOW_FACTOR + 0.5);
2053  }
2054  else
2055  {
2056  capa[k] = (int) (xval[krev] * FLOW_FACTOR + 0.5);
2057  capa[krev] = (int) (xval[k] * FLOW_FACTOR + 0.5);
2058  }
2059 
2060  if( creep_flow )
2061  {
2062  capa[k] += CREEP_VALUE;
2063  capa[krev] += CREEP_VALUE;
2064  }
2065  }
2066 }
2067 
2068 /*
2069  * Interface methods
2070  */
2071 
2072 
2073 /** initializes */
2075  SCIP* scip, /**< SCIP */
2076  const GRAPH* g, /**< graph */
2077  int maxnsepas, /**< maximum number of separators to compute */
2078  int maxsepasize, /**< maximum size of separator */
2079  TERMSEPAS** termsepas /**< to initialize */
2080 )
2081 {
2082  TERMSEPAS* tsepas;
2083  assert(scip && g);
2084  assert(maxnsepas >= 1);
2085  assert(maxsepasize >= 2);
2086 
2087  SCIP_CALL( SCIPallocMemory(scip, termsepas) );
2088  tsepas = *termsepas;
2089 
2090  SCIP_CALL( SCIPallocMemoryArray(scip, &(tsepas->nsepas), maxsepasize + 1) );
2091  SCIP_CALL( SCIPallocMemoryArray(scip, &(tsepas->currsepa_n), maxsepasize + 1) );
2092  SCIP_CALL( SCIPallocMemoryArray(scip, &(tsepas->sepas), maxnsepas) );
2093  SCIP_CALL( SCIPallocMemoryArray(scip, &(tsepas->sepastarts_csr), maxnsepas + 1) );
2094  SCIP_CALL( SCIPallocMemoryArray(scip, &(tsepas->sepaterms_csr), maxnsepas * maxsepasize + 1) );
2095  tsepas->maxnsepas = maxnsepas;
2096  tsepas->maxsepasize = maxsepasize;
2097  tsepas->nsepas_all = 0;
2098  tsepas->nsepaterms_csr = 0;
2099  tsepas->sepastarts_csr[0] = 0;
2100  tsepas->root = -1;
2101 
2102  for( int i = 0; i <= maxsepasize; i++ )
2103  {
2104  tsepas->nsepas[i] = 0;
2105  tsepas->currsepa_n[i] = -1;
2106  }
2107 
2108  return SCIP_OKAY;
2109 }
2110 
2111 
2112 /** frees */
2114  SCIP* scip, /**< SCIP */
2115  TERMSEPAS** termsepas /**< to free */
2116 )
2117 {
2118  TERMSEPAS* tsepas;
2119  assert(scip && termsepas);
2120 
2121  tsepas = *termsepas;
2122  assert(tsepas);
2123 
2124  SCIPfreeMemoryArray(scip, &(tsepas->sepaterms_csr));
2125  SCIPfreeMemoryArray(scip, &(tsepas->sepastarts_csr));
2126  SCIPfreeMemoryArray(scip, &(tsepas->sepas));
2127  SCIPfreeMemoryArray(scip, &(tsepas->currsepa_n));
2128  SCIPfreeMemoryArray(scip, &(tsepas->nsepas));
2129 
2130  SCIPfreeMemory(scip, termsepas);
2131 }
2132 
2133 
2134 /** returns number of all separators */
2136  const TERMSEPAS* termsepas /**< terminal separators */
2137 )
2138 {
2139  assert(termsepas);
2140  assert(termsepas->nsepas_all >= 0);
2141 
2142  return termsepas->nsepas_all;
2143 }
2144 
2145 
2146 /** returns number of separators per given size */
2148  const TERMSEPAS* termsepas, /**< terminal separators */
2149  int sepasize /**< size */
2150 )
2151 {
2152  assert(termsepas);
2153  assert(sepasize >= 1 && sepasize <= termsepas->maxsepasize);
2154  assert(termsepas->nsepas[sepasize] >= 0);
2155  assert(termsepas->nsepas[sepasize] <= termsepas->maxnsepas);
2156 
2157  return termsepas->nsepas[sepasize];
2158 }
2159 
2160 
2161 /** Returns first separator of given size. Returns NULL if none is available. */
2163  int sepasize, /**< size */
2164  TERMSEPAS* termsepas, /**< terminal separators */
2165  int* sinkterm, /**< the sink */
2166  int* nsinknodes /**< number of sink-side nodes */
2167 )
2168 {
2169  assert(termsepas && sinkterm && nsinknodes);
2170  assert(sepasize >= 1 && sepasize <= termsepas->maxsepasize);
2171 
2172  termsepas->currsepa_n[sepasize] = -1;
2173 
2174  return mincut_termsepasGetNext(sepasize, termsepas, sinkterm, nsinknodes);
2175 }
2176 
2177 
2178 /** Returns next separator of given size. Returns NULL if none is available. */
2180  int sepasize, /**< size */
2181  TERMSEPAS* termsepas, /**< terminal separators */
2182  int* sinkterm, /**< the sink */
2183  int* nsinknodes /**< number of sink-side nodes */
2184 )
2185 {
2186  assert(termsepas && sinkterm && nsinknodes);
2187  assert(sepasize >= 1 && sepasize <= termsepas->maxsepasize);
2188  assert(termsepas->currsepa_n[sepasize] >= -1);
2189 
2190  *sinkterm = -1;
2191  *nsinknodes = -1;
2192 
2193  if( termsepas->currsepa_n[sepasize] >= termsepas->nsepas_all )
2194  {
2195  assert(termsepas->currsepa_n[sepasize] == termsepas->nsepas_all);
2196  }
2197  else
2198  {
2199  const int* const starts = termsepas->sepastarts_csr;
2200  const int* const terms = termsepas->sepaterms_csr;
2201  const int startpos = termsepas->currsepa_n[sepasize] + 1;
2202  int s;
2203 
2204  assert(0 <= startpos && startpos <= termsepas->nsepas_all);
2205 
2206  for( s = startpos; s < termsepas->nsepas_all; s++ )
2207  {
2208  const int size = starts[s + 1] - starts[s];
2209  assert(sepasize >= 1 && sepasize <= termsepas->maxsepasize);
2210 
2211  if( size == sepasize )
2212  break;
2213  }
2214 
2215  termsepas->currsepa_n[sepasize] = s;
2216 
2217  if( s < termsepas->nsepas_all )
2218  {
2219  TSEPA tsepa = termsepas->sepas[s];
2220  *sinkterm = tsepa.sinkterm;
2221  *nsinknodes = tsepa.nsinknodes;
2222 
2223  return &(terms[starts[s]]);
2224  }
2225  }
2226 
2227  return NULL;
2228 }
2229 
2230 
2231 /** gets terminal separator source */
2233  const TERMSEPAS* termsepas /**< terminal separators */
2234 )
2235 {
2236  assert(termsepas);
2237 
2238  return termsepas->root;
2239 }
2240 
2241 
2242 /** is it promising to look for terminal separators? */
2244  const GRAPH* g /**< graph data structure */
2245  )
2246 {
2247  int nnodes;
2248  int nedges;
2249 
2250  assert(g);
2251 
2252  if( g->terms < 4 )
2253  {
2254  return FALSE;
2255  }
2256 
2257  graph_get_nVET(g, &nnodes, &nedges, NULL);
2258 
2259  if( nedges == 0 )
2260  {
2261  return FALSE;
2262  }
2263 
2264  assert(nnodes > 0 && nedges >= 2);
2265  assert(nedges % 2 == 0);
2266 
2267  nedges /= 2;
2268 
2269  return ((nedges / nnodes) <= TERMSEPA_SPARSE_MAXRATIO);
2270 }
2271 
2272 
2273 /** searches for (small) terminal separators */
2275  SCIP* scip, /**< SCIP data structure */
2276  SCIP_RANDNUMGEN* randnumgen, /**< random number generator or NULL */
2277  GRAPH* g, /**< graph data structure */
2278  TERMSEPAS* termsepas /**< terminal separator storage */
2279  )
2280 {
2281  MINCUT* mincut;
2282  int* nodes_wakeState;
2283  SCIP_Bool wasRerun;
2284 
2285  assert(scip && g && termsepas);
2286  assert(mincut_termsepasGetNall(termsepas) == 0);
2287 
2288  if( g->terms < 3 )
2289  {
2290  SCIPdebugMessage("not enough terminals (%d), aborting terminal separator check \n", g->terms);
2291  return SCIP_OKAY;
2292  }
2293 
2294  SCIP_CALL( reduce_unconnected(scip, g) );
2295 
2296  SCIP_CALL( mincutInit(scip, randnumgen, FALSE, g, &mincut) );
2297  termsepas->root = mincut->root;
2298 
2299  /* sets excess, g->mincut_head, g->mincut_head_inact */
2301 
2302  SCIP_CALL( mincutPrepareForTermSepa(scip, g, mincut) );
2303  nodes_wakeState = mincut->nodes_wakeState;
2304  wasRerun = FALSE;
2305 
2306  assert(nodes_wakeState);
2307  assert(termsepas->nsepas_all == 0);
2308 
2309 #ifdef SCIP_DEBUG
2311  SCIPdebugMessage("number of terminal separation candidates: %d \n", mincut->ntermcands );
2312 #endif
2313 
2314  // todo probably want to bound the maximum number of iterations!
2315  while( mincut->ntermcands > 0 && termsepas->nsepas_all < termsepas->maxnsepas )
2316  {
2317  int sinkterm;
2318 
2319  if( ((unsigned) mincut->ntermcands) % 32 == 0 && SCIPisStopped(scip) )
2320  break;
2321 
2322  /* look for non-reachable terminal */
2323  sinkterm = mincutGetNextSinkTerm(g, !wasRerun, mincut);
2324  mincut->ntermcands--;
2325 
2326  SCIPdebugMessage("computing cut for sink terminal %d \n", sinkterm);
2327 
2328  assert(Is_term(g->term[sinkterm]) && mincut->root != sinkterm);
2329 
2330  /* non-trivial cut? */
2331  if( nodes_wakeState[sinkterm] != 1 )
2332  {
2333  mincutExec(g, sinkterm, wasRerun, mincut);
2334  assert(nodes_wakeState[mincut->root] != 0);
2335 
2336  SCIP_CALL( termsepaStoreCutTry(scip, g, sinkterm, mincut, termsepas) );
2337 #ifdef SCIP_DEBUG
2338  debugPrintCsrCutEdges(g, mincut);
2339 #endif
2340  }
2341  else
2342  {
2343  SCIPdebugMessage("cut is trivial \n");
2344 
2345  assert(wasRerun);
2346  }
2347 
2348  wasRerun = TRUE;
2349  }
2350 
2351  SCIPdebugMessage("number of separators: %d \n", termsepas->nsepas_all);
2352 
2353  mincutFree(scip, &mincut);
2354 
2355  return SCIP_OKAY;
2356 }
2357 
2358 
2359 /** separates Steiner cuts for LP */
2361  SCIP* scip, /**< SCIP data structure */
2362  SCIP_CONSHDLR* conshdlr, /**< constraint handler */
2363  SCIP_RANDNUMGEN* randnumgen, /**< random number generator or NULL */
2364  const int* termorg, /**< original terminals or NULL */
2365  GRAPH* g, /**< graph data structure */
2366  int maxncuts, /**< maximal number of cuts */
2367  int* ncuts /**< pointer to store number of cuts */
2368  )
2369 {
2370  /* we do not longer support any other flow, since they slow everything down and are of
2371  * little use anyway todo remove user parameter */
2372 #ifdef SCIP_DISABLED
2373  const SCIP_Bool nested_cut = conshdlrdata->nestedcut;
2374  const SCIP_Bool back_cut = conshdlrdata->backcut;
2375  const SCIP_Bool creep_flow = conshdlrdata->creepflow;
2376  const SCIP_Bool disjunct_cut = conshdlrdata->disjunctcut;
2377 #endif
2378  const SCIP_Bool nested_cut = FALSE;
2379  const SCIP_Bool creep_flow = TRUE;
2380  const SCIP_Bool disjunct_cut = FALSE;
2381  MINCUT* mincut;
2382  const SCIP_Real* xval;
2383  int* nodes_wakeState;
2385  const int nnodes = graph_get_nNodes(g);
2386  int cutcount;
2387  SCIP_Bool wasRerun;
2388 
2389 #ifdef STP_MAXFLOW_TIME
2390  clock_t startt, endt;
2391  double cpu_time_used;
2392  startt = clock();
2393 #endif
2394 
2395  assert(conshdlr && ncuts);
2396  assert(creep_flow == TRUE);
2397  assert(nested_cut == FALSE);
2398  assert(disjunct_cut == FALSE);
2399 
2400  SCIP_CALL( mincutInit(scip, randnumgen, TRUE, g, &mincut) );
2401 
2402  /* sets excess, g->mincut_head, g->mincut_head_inact */
2404 
2405 #ifdef STP_MAXFLOW_TIME
2406  endt = clock();
2407  cpu_time_used = ((double) (endt - startt)) / CLOCKS_PER_SEC;
2408  startt = clock();
2409 #endif
2410 
2411  SCIP_CALL( mincutPrepareForLp(scip, g, mincut) );
2412 
2413  xval = mincut->xval;
2414  nodes_wakeState = mincut->nodes_wakeState;
2415  edges_isRemoved = mincut->edges_isRemoved;
2416  cutcount = 0;
2417  wasRerun = FALSE;
2418 
2419  assert(xval && nodes_wakeState && edges_isRemoved);
2420 
2421  while( mincut->ntermcands > 0 )
2422  {
2423  int sinkterm;
2424  SCIP_Bool addedcut = FALSE;
2425 
2426  if( ((unsigned) mincut->ntermcands) % 32 == 0 && SCIPisStopped(scip) )
2427  break;
2428 
2429  /* look for non-reachable terminal */
2430  sinkterm = mincutGetNextSinkTerm(g, !wasRerun, mincut);
2431  mincut->ntermcands--;
2432 
2433  assert(Is_term(g->term[sinkterm]) && g->source != sinkterm);
2434 
2435  if( nested_cut && !disjunct_cut )
2436  lpcutSetEdgeCapacity(g, xval, creep_flow, 0, mincut->edges_capa);
2437 
2438  do
2439  {
2440  /* declare cuts on branched-on (artificial) terminals as local */
2441  const SCIP_Bool localcut = (termorg != NULL && termorg[sinkterm] != g->term[sinkterm]);
2442 
2443 #ifdef STP_MAXFLOW_WRITE
2444  writeFlowProb(g, edges_capa, sinkterm);
2445 #endif
2446 
2447  /* non-trivial cut? */
2448  if( nodes_wakeState[sinkterm] != 1 )
2449  {
2450  mincutExec(g, sinkterm, wasRerun, mincut);
2451  assert(nodes_wakeState[g->source] != 0);
2452 
2453  SCIP_CALL( lpcutAdd(scip, conshdlr, g, nodes_wakeState, edges_isRemoved, xval,
2454  mincut->edges_capa, nested_cut || disjunct_cut, localcut, &addedcut) );
2455  }
2456  else
2457  {
2458  SCIP_Real flowsum = 0.0;
2459  assert(wasRerun);
2460 
2461  for( int e = g->inpbeg[sinkterm]; e != EAT_LAST; e = g->ieat[e] )
2462  flowsum += xval[e];
2463 
2464  if( SCIPisFeasGE(scip, flowsum, 1.0) )
2465  continue;
2466 
2467  for( int k = 0; k < nnodes; k++ )
2468  g->mark[k] = TRUE;
2469 
2470  g->mark[sinkterm] = FALSE;
2471 
2472  SCIP_CALL( lpcutAdd(scip, conshdlr, g, g->mark, edges_isRemoved, xval,
2473  mincut->edges_capa, nested_cut || disjunct_cut, localcut, &addedcut) );
2474  }
2475 
2476  wasRerun = TRUE;
2477 
2478  if( addedcut )
2479  {
2480  cutcount++;
2481  (*ncuts)++;
2482 
2483  if( *ncuts >= maxncuts )
2484  goto TERMINATE;
2485  }
2486  else
2487  {
2488  break;
2489  }
2490  }
2491  while( nested_cut ); /* Nested Cut is CONSTANT ! */
2492  } /* while terms > 0 */
2493 
2494 #ifdef STP_MAXFLOW_TIME
2495  endt = clock();
2496  cpu_time_used = ((double) (endt - startt)) / CLOCKS_PER_SEC;
2497 #endif
2498 
2499 #ifdef SCIP_DISABLED
2500  /*
2501  * back cuts currently not supported
2502  * */
2503  /* back cuts enabled? */
2504  if( back_cut )
2505  {
2506  for( k = 0; k < nnodes; k++ )
2507  nodes_wakeState[k] = 0;
2508 
2509  if( !nested_cut || disjunct_cut )
2510  lpcutSetEdgeCapacity(g, creep_flow, 1, edges_capa, xval);
2511 
2512  ntermcands = tsave;
2513 
2514  while( ntermcands > 0 )
2515  {
2516  /* look for reachable terminal */
2517  i = mincutGetNextSinkTerm(g, ntermcands, terms, nodes_wakeState, TRUE);
2518 
2519  ntermcands--;
2520 
2521  assert(g->terms[i] == 0);
2522  assert(g->source != i);
2523 
2524  if( nested_cut && !disjunct_cut )
2525  lpcutSetEdgeCapacity(g, creep_flow, 1, edges_capa, xval);
2526 
2527  wasRerun = FALSE;
2528 
2529  do
2530  {
2531  graph_mincut_exec(g, i, g->source, nedges, edges_capa, nodes_wakeState, csr_start, csr_edgeDefaultToCsr, csr_headarr, wasRerun);
2532 
2533  wasRerun = TRUE;
2534 
2535  for( k = 0; k < nnodes; k++ )
2536  {
2537  g->mark[k] = (nodes_wakeState[k] != 0) ? 0 : 1; // todo not the other way around??
2538  nodes_wakeState[k] = 0;
2539  }
2540 
2541  SCIP_CALL( lpcutAdd(scip, conshdlr, g, xval, edges_capa, nested_cut || disjunct_cut, ncuts, &addedcut) );
2542  if( addedcut )
2543  {
2544  cutcount++;
2545 
2546  if( *ncuts >= maxncuts )
2547  goto TERMINATE;
2548  }
2549  else
2550  break;
2551 #ifdef SCIP_DISABLED
2552  if (nested_cut || disjunct_cut)
2553  for(k = p->beg[p->rcnt - 1]; k < p->nzcnt; k++)
2554  edges_capa[p->ind[k] % nedges
2555  + (((p->ind[k] % nedges) % 2)
2556  ? -1 : 1)] = FLOW_FACTOR;
2557 #endif
2558  }
2559  while( nested_cut ); /* Nested Cut is CONSTANT todo why not only one round? seems to make no sense whatsoever */
2560 
2561  wasRerun = FALSE;
2562  }
2563  }
2564 #endif
2565  TERMINATE:
2566  mincutFree(scip, &mincut);
2567 
2568  SCIPdebugMessage("2-cut Separator: %d Inequalities added\n", cutcount);
2569 
2570  return SCIP_OKAY;
2571 }
#define STP_Vectype(type)
Definition: stpvector.h:44
int graph_get_nEdges(const GRAPH *)
Definition: graph_stats.c:548
static volatile int nterms
Definition: interrupt.c:38
int *RESTRICT mincut_e
Definition: graphdefs.h:250
SCIP_Bool isLpcut
Definition: mincut.c:75
#define StpVecGetSize(vec)
Definition: stpvector.h:139
#define FLOW_FACTOR
Definition: mincut.c:49
#define CREEP_VALUE
Definition: mincut.c:50
SCIP_RETCODE SCIPcacheRowExtensions(SCIP *scip, SCIP_ROW *row)
Definition: scip_lp.c:1626
int *RESTRICT head
Definition: graphdefs.h:224
int source
Definition: graphdefs.h:195
#define Is_term(a)
Definition: graphdefs.h:102
static void termsepaCsrAddEdges(const GRAPH *g, MINCUT *mincut)
Definition: mincut.c:966
SCIP_RETCODE SCIPflushRowExtensions(SCIP *scip, SCIP_ROW *row)
Definition: scip_lp.c:1649
#define SCIPfreeMemoryArray(scip, ptr)
Definition: scip_mem.h:71
int terms
Definition: graphdefs.h:192
SCIP_RETCODE SCIPaddVarToRow(SCIP *scip, SCIP_ROW *row, SCIP_VAR *var, SCIP_Real val)
Definition: scip_lp.c:1686
SCIP_VAR ** SCIPprobdataGetVars(SCIP *scip)
static void updateTerminalSource(SCIP *scip, int rootcand, const GRAPH *g, int *root, int *rootcompsize)
Definition: mincut.c:419
#define SCIPallocMemoryArray(scip, ptr, num)
Definition: scip_mem.h:55
int * csr_headarr
Definition: mincut.c:62
void graph_knot_printInfo(const GRAPH *, int)
Definition: graph_stats.c:151
SCIP_Bool SCIPisFeasGE(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
SCIP_Bool mincut_findTerminalSeparatorsIsPromising(const GRAPH *g)
Definition: mincut.c:2243
#define FALSE
Definition: def.h:87
int mincut_termsepasGetSource(const TERMSEPAS *termsepas)
Definition: mincut.c:2232
int *RESTRICT inpbeg
Definition: graphdefs.h:202
static SCIP_RETCODE lpcutAdd(SCIP *scip, SCIP_CONSHDLR *conshdlr, const GRAPH *g, const int *nodes_inRootComp, const STP_Bool *edges_isRemoved, const SCIP_Real *xvals, int *capa, const int updatecapa, SCIP_Bool local, SCIP_Bool *success)
Definition: mincut.c:1929
SCIP_Real SCIPinfinity(SCIP *scip)
Problem data for stp problem.
#define TRUE
Definition: def.h:86
#define SCIPdebug(x)
Definition: pub_message.h:84
#define SWAP_INTS(int1, int2)
Definition: misc_stp.h:34
enum SCIP_Retcode SCIP_RETCODE
Definition: type_retcode.h:54
static void termsepaCsrAddReverseEdges(const GRAPH *g, MINCUT *mincut)
Definition: mincut.c:1121
static void mincutExec(const GRAPH *g, int sinkterm, SCIP_Bool wasRerun, MINCUT *mincut)
Definition: mincut.c:1779
#define StpVecPopBack(vec)
Definition: stpvector.h:178
const SCIP_Real * xval
Definition: mincut.c:55
#define EAT_LAST
Definition: graphdefs.h:58
#define StpVecReserve(scip, vec, _size_)
Definition: stpvector.h:182
#define SCIPdebugMessage
Definition: pub_message.h:87
int SCIPrandomGetInt(SCIP_RANDNUMGEN *randnumgen, int minrandval, int maxrandval)
Definition: misc.c:10003
SCIP_RETCODE graph_mincut_reInit(SCIP *, int, int, GRAPH *)
Definition: graph_mcut.c:1120
static void termsepaCollectCutNodes(const GRAPH *g, const TERMSEPAS *termsepas, const MINCUT *mincut, int sinkterm, int *cutterms, int *ncutterms, SCIP_Bool *cutIsGood)
Definition: mincut.c:533
#define SCIPfreeBufferArray(scip, ptr)
Definition: scip_mem.h:127
static SCIP_RETCODE termsepaStoreCutTry(SCIP *scip, const GRAPH *g, int sinkterm, MINCUT *mincut, TERMSEPAS *termsepas)
Definition: mincut.c:829
int mincut_termsepasGetNall(const TERMSEPAS *termsepas)
Definition: mincut.c:2135
header only, simple implementation of an STL like vector
#define TERMSEPA_SPARSE_MAXRATIO
Definition: mincut.c:39
int *RESTRICT mark
Definition: graphdefs.h:198
static SCIP_RETCODE mincutInitForTermSepa(SCIP *scip, GRAPH *g, MINCUT *mincut)
Definition: mincut.c:1353
int termsepa_nnodes
Definition: mincut.c:72
int termsepa_nedges
Definition: mincut.c:73
#define SCIPallocCleanBufferArray(scip, ptr, num)
Definition: scip_mem.h:133
SCIP_Real SCIPvarGetUbGlobal(SCIP_VAR *var)
Definition: var.c:17920
static int mincutGetNextSinkTerm(const GRAPH *g, SCIP_Bool firstrun, MINCUT *mincut)
Definition: mincut.c:1697
SCIP_VAR * w
Definition: circlepacking.c:58
int *RESTRICT oeat
Definition: graphdefs.h:231
#define TERMSEPA_NROOTCANDS
Definition: mincut.c:42
SCIP_Bool SCIPisCutEfficacious(SCIP *scip, SCIP_SOL *sol, SCIP_ROW *cut)
Definition: scip_cut.c:108
STP_Bool * edges_isRemoved
Definition: mincut.c:67
static SCIP_RETCODE mincutPrepareForLp(SCIP *scip, const GRAPH *g, MINCUT *mincut)
Definition: mincut.c:1476
int *RESTRICT mincut_dist
Definition: graphdefs.h:243
miscellaneous methods used for solving Steiner problems
int * edges_capa
Definition: mincut.c:57
int graph_get_nTerms(const GRAPH *)
Definition: graph_stats.c:560
#define SCIPerrorMessage
Definition: pub_message.h:55
void graph_get_nVET(const GRAPH *, int *, int *, int *)
Definition: graph_stats.c:571
void graph_mincut_setDefaultVals(GRAPH *)
Definition: graph_mcut.c:1081
static SCIP_RETCODE mincutInitForLp(SCIP *scip, const GRAPH *g, MINCUT *mincut)
Definition: mincut.c:1299
SCIP_RETCODE mincut_findTerminalSeparators(SCIP *scip, SCIP_RANDNUMGEN *randnumgen, GRAPH *g, TERMSEPAS *termsepas)
Definition: mincut.c:2274
static SCIP_Bool termsepaCutIsCorrect(SCIP *scip, const GRAPH *g, int ncutterms, const int *cutterms, int sinkterm, const MINCUT *mincut)
Definition: mincut.c:310
#define SCIPfreeBufferArrayNull(scip, ptr)
Definition: scip_mem.h:128
static SCIP_RETCODE termsepaCsrInit(const GRAPH *g, MINCUT *mincut)
Definition: mincut.c:865
int *RESTRICT grad
Definition: graphdefs.h:201
SCIP_RETCODE reduce_unconnected(SCIP *, GRAPH *)
struct terminal_separator TSEPA
static int termsepaCsrGetMaxNnodes(const GRAPH *g)
Definition: mincut.c:1186
#define NULL
Definition: lpi_spx1.cpp:155
#define EQ(a, b)
Definition: portab.h:79
int knots
Definition: graphdefs.h:190
#define SCIP_CALL(x)
Definition: def.h:384
static SCIP_RETCODE mincutPrepareForTermSepa(SCIP *scip, GRAPH *g, MINCUT *mincut)
Definition: mincut.c:1678
int * termsepa_termToCopy
Definition: mincut.c:64
SCIP_RETCODE SCIPaddRow(SCIP *scip, SCIP_ROW *row, SCIP_Bool forcecut, SCIP_Bool *infeasible)
Definition: scip_cut.c:241
static SCIP_RETCODE mincutInit(SCIP *scip, SCIP_RANDNUMGEN *randnumgen, SCIP_Bool isLpcut, GRAPH *g, MINCUT **mincut)
Definition: mincut.c:1425
unsigned char STP_Bool
Definition: portab.h:34
static void lpcutSetEdgeCapacity(const GRAPH *g, const SCIP_Real *xval, SCIP_Bool creep_flow, SCIP_Bool flip, int *RESTRICT capa)
Definition: mincut.c:2030
#define SCIPallocBufferArray(scip, ptr, num)
Definition: scip_mem.h:115
#define SCIP_Bool
Definition: def.h:84
int *RESTRICT ieat
Definition: graphdefs.h:230
int * terms_mincompsize
Definition: mincut.c:66
int * csr_edgeflipped
Definition: mincut.c:63
int *RESTRICT tail
Definition: graphdefs.h:223
int SCIPgetDepth(SCIP *scip)
Definition: scip_tree.c:661
int mincut_termsepasGetN(const TERMSEPAS *termsepas, int sepasize)
Definition: mincut.c:2147
void mincut_termsepasFree(SCIP *scip, TERMSEPAS **termsepas)
Definition: mincut.c:2113
void SCIPrandomPermuteIntArray(SCIP_RANDNUMGEN *randnumgen, int *array, int begin, int end)
Definition: misc.c:10044
#define flipedge(edge)
Definition: graphdefs.h:84
SCIP_RETCODE SCIPaddPoolCut(SCIP *scip, SCIP_ROW *row)
Definition: scip_cut.c:352
static SCIP_Bool csrFlipedgesAreValid(const GRAPH *g, const MINCUT *mincut)
Definition: mincut.c:273
void graph_mincut_exec(const GRAPH *, const int, const int, const int, const int, const int, const int *, const int *, int *RESTRICT, const int *, const int *, const int *, const SCIP_Bool)
int *RESTRICT term
Definition: graphdefs.h:196
#define BMScopyMemoryArray(ptr, source, num)
Definition: memory.h:127
void graph_edge_printInfo(const GRAPH *, int)
Definition: graph_stats.c:133
struct minimum_cut_helper MINCUT
static int termsepaFindTerminalSource(SCIP *scip, const GRAPH *g, const MINCUT *mincut)
Definition: mincut.c:489
static void mincutFree(SCIP *scip, MINCUT **mincut)
Definition: mincut.c:1812
static void termsepaCsrAddTermCopies(const GRAPH *g, MINCUT *mincut)
Definition: mincut.c:884
Portable definitions.
SCIP_RETCODE mincut_termsepasInit(SCIP *scip, const GRAPH *g, int maxnsepas, int maxsepasize, TERMSEPAS **termsepas)
Definition: mincut.c:2074
int * terms_minsepasize
Definition: mincut.c:65
int *RESTRICT mincut_numb
Definition: graphdefs.h:246
#define SCIPfreeMemory(scip, ptr)
Definition: scip_mem.h:69
SCIP_RETCODE SCIPreleaseRow(SCIP *scip, SCIP_ROW **row)
Definition: scip_lp.c:1553
const int * mincut_termsepasGetFirst(int sepasize, TERMSEPAS *termsepas, int *sinkterm, int *nsinknodes)
Definition: mincut.c:2162
SCIP_RANDNUMGEN * randnumgen
Definition: mincut.c:74
static SCIP_RETCODE termsepaBuildCsr(const GRAPH *g, MINCUT *mincut)
Definition: mincut.c:1275
static int termsepaCsrGetMaxNedges(int root, const GRAPH *g)
Definition: mincut.c:1201
const int * mincut_termsepasGetNext(int sepasize, TERMSEPAS *termsepas, int *sinkterm, int *nsinknodes)
Definition: mincut.c:2179
#define SCIP_Real
Definition: def.h:177
#define SCIPfreeCleanBufferArray(scip, ptr)
Definition: scip_mem.h:137
SCIP_Bool SCIPisStopped(SCIP *scip)
Definition: scip_general.c:694
#define StpVecFree(scip, vec)
Definition: stpvector.h:149
SCIP_Real * SCIPprobdataGetXval(SCIP *scip, SCIP_SOL *sol)
int *RESTRICT outbeg
Definition: graphdefs.h:204
SCIP_RETCODE SCIPprintRow(SCIP *scip, SCIP_ROW *row, FILE *file)
Definition: scip_lp.c:2197
int edges
Definition: graphdefs.h:219
int graph_get_nNodes(const GRAPH *)
Definition: graph_stats.c:536
int * csr_start
Definition: mincut.c:59
#define SCIPallocMemory(scip, ptr)
Definition: scip_mem.h:51
SCIP_Bool graph_knot_isInRange(const GRAPH *, int)
Definition: graph_node.c:52
void graph_getTerms(const GRAPH *, int *)
Definition: graph_base.c:567
int *RESTRICT mincut_r
Definition: graphdefs.h:252
SCIP_Real SCIPvarGetUbLocal(SCIP_VAR *var)
Definition: var.c:17976
#define nnodes
Definition: gastrans.c:65
#define StpVecPushBack(scip, vec, value)
Definition: stpvector.h:163
int * csr_edgeDefaultToCsr
Definition: mincut.c:61
#define SCIP_CALL_ABORT(x)
Definition: def.h:363
SCIP_RETCODE SCIPcreateEmptyRowConshdlr(SCIP *scip, SCIP_ROW **row, SCIP_CONSHDLR *conshdlr, const char *name, SCIP_Real lhs, SCIP_Real rhs, SCIP_Bool local, SCIP_Bool modifiable, SCIP_Bool removable)
Definition: scip_lp.c:1382
Minimum cut routines for Steiner problems.
int * nodes_wakeState
Definition: mincut.c:56
includes various reduction methods for Steiner tree problems
static void termsepaBuildRootcomp(const GRAPH *g, MINCUT *mincut)
Definition: mincut.c:1229
static SCIP_RETCODE termsepaGetCompNnodes(SCIP *scip, const GRAPH *g, int ncutterms, const int *cutterms, int sinkterm, MINCUT *mincut, int *ncompnodes)
Definition: mincut.c:730
void graph_printInfoReduced(const GRAPH *)
Definition: graph_stats.c:375
static int termsepaGetCapaInf(const GRAPH *g, const MINCUT *mincut)
Definition: mincut.c:408
SCIP_RETCODE mincut_separateLp(SCIP *scip, SCIP_CONSHDLR *conshdlr, SCIP_RANDNUMGEN *randnumgen, const int *termorg, GRAPH *g, int maxncuts, int *ncuts)
Definition: mincut.c:2360
static SCIP_RETCODE termsepaStoreCutFinalize(SCIP *scip, const GRAPH *g, int sinkterm, MINCUT *mincut, int ncutterms, const int *cutterms, TERMSEPAS *termsepas, SCIP_Bool *success)
Definition: mincut.c:752
static SCIP_RETCODE termsepaTraverseSinkComp(SCIP *scip, const GRAPH *g, SCIP_Bool removeTerms, SCIP_Bool updateVisitedTerms, int ncutterms, const int *cutterms, int sinkterm, MINCUT *mincut, int *ncompnodes)
Definition: mincut.c:590
static SCIP_RETCODE termsepaRemoveCutTerminals(SCIP *scip, const GRAPH *g, int ncutterms, const int *cutterms, int sinkterm, MINCUT *mincut)
Definition: mincut.c:710