Scippy

SCIP

Solving Constraint Integer Programs

GomoryHuTree.cpp
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-2018 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 GomoryHuTree.cpp
17  * @brief generator for global cuts in undirected graphs
18  * @author Georg Skorobohatyj
19  * @author Timo Berthold
20  */
21 
22 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
23 
24 
25 #include <stdio.h>
26 #include <assert.h>
27 
28 #include "objscip/objscip.h"
29 #include "GomoryHuTree.h"
30 
31 #define EPS 1.0E-10
32 
33 static GRAPHNODE **active;
34 static long *number;
35 static long max_dist, bound;
37 
38 SCIP_Bool create_graph (int n, int m, GRAPH** gr)
39 {
40  assert( gr != NULL );
41 
42  BMSallocMemory(gr);
43  if( *gr == NULL )
44  return FALSE;
45 
46  BMSallocMemoryArray( &(*gr)->nodes, n );
47  if( (*gr)->nodes == NULL )
48  {
49  BMSfreeMemory(gr);
50  return FALSE;
51  }
52 
53  BMSallocMemoryArray( &(*gr)->edges, m );
54  if( (*gr)->edges == NULL )
55  {
56  BMSfreeMemoryArray(&(*gr)->nodes);
57  BMSfreeMemory(gr);
58  return FALSE;
59  }
60  (*gr)->nuses = 1;
61  (*gr)->nnodes = n;
62  (*gr)->nedges = m/2;
63  (*gr)->nedgesnonzero = m/2;
64  return TRUE;
65 }
66 
67 static
68 void free_graph (GRAPH** gr)
69 {
70  assert(gr != NULL);
71  assert(*gr != NULL);
72  assert((*gr)->nuses == 0);
73  BMSfreeMemory(&(*gr)->nodes);
74  BMSfreeMemory(&(*gr)->edges);
75  BMSfreeMemory(gr);
76 }
77 
78 void capture_graph (GRAPH* gr)
79 {
80  assert(gr != NULL);
81  gr->nuses++;
82 }
83 
84 void release_graph (GRAPH** gr)
85 {
86  assert(gr != NULL);
87  assert(*gr != NULL);
88  (*gr)->nuses--;
89  if( (*gr)->nuses == 0 )
90  free_graph(gr);
91  *gr = NULL;
92 }
93 
95 {
96  active = (GRAPHNODE **) malloc ((n+1L) * sizeof (GRAPHNODE *));
97  /* holds stacks of active nodes arranged by distances */
98  if ( active == (GRAPHNODE **) 0 )
99  {
100  printf ("Unable to allocate memory\n");
101  return FALSE;
102  }
103  number = (long *) malloc ((n+1L) * sizeof (long));
104  /* counts occurences of node distances in set
105  of alive nodes, i.e. nodes not contained in
106  the set of nodes disconnected from the sink */
107  if ( number == (long *) 0 )
108  {
109  printf ("Unable to allocate memory\n");
110  return FALSE;
111  }
112  co_check = TRUE;
113  return TRUE;
114 
115 } /* init_maxflow */
116 
117 
119 {
120  free (active);
121  free (number);
122 
123 } /* fini_maxflow */
124 
125 
126 void global_relabel (GRAPH *gr, GRAPHNODE *tptr)
127 {
128  /* breadth first search to get exact distance labels
129  from sink with reordering of stack of active nodes */
130 
131  GRAPHNODE *front, *rear, *nptr, **ptr;
132  GRAPHEDGE *eptr;
133  long n, level, count, i;
134 
135  n = gr->nnodes;
136  for ( nptr = &(gr->nodes[n-1L]); nptr >= gr->nodes; nptr-- )
137  {
138  nptr->unmarked = TRUE;
139  nptr->stack_link = NULL;
140  nptr->scan_ptr = nptr->first_edge;
141  while( nptr->scan_ptr != NULL && nptr->scan_ptr->cap <= EPS )
142  nptr->scan_ptr = nptr->scan_ptr->next;
143  }
144  tptr->unmarked = FALSE;
145 
146  /* initialize stack of active nodes */
147  for ( ptr = &(active[n]); ptr >= active; ptr-- )
148  *ptr = NULL;
149  for ( i = 0L; i <= n; i++ )
150  number[i] = 0L;
151  max_dist = 0L;
152  count = 1L; /* number of alive nodes */
153  front = tptr;
154  rear = front;
155 
156  bfs_next:
157  level = rear->dist + 1L;
158  eptr = rear->first_edge;
159  while ( eptr != NULL )
160  {
161  nptr = eptr->adjac;
162  if ( nptr->alive && nptr->unmarked && eptr->back->rcap > EPS )
163  {
164  assert(eptr->cap > EPS);
165  assert(eptr->back->cap > EPS);
166  nptr->unmarked = FALSE;
167  nptr->dist = (int) level;
168  ++count;
169  ++number[level];
170  if ( nptr->excess > EPS )
171  {
172  nptr->stack_link = active[level];
173  active[level] = nptr;
174  max_dist = level;
175  }
176  front->bfs_link = nptr;
177  front = nptr;
178  }
179  eptr = eptr->next;
180  }
181  if ( front == rear )
182  goto bfs_ready;
183 
184  rear = rear->bfs_link;
185  goto bfs_next;
186 
187  bfs_ready:
188 
189  if ( count < bound )
190  {
191  /* identify nodes that are marked alive but have
192  not been reached by BFS and mark them as dead */
193  for ( nptr = &(gr->nodes[n-1L]); nptr >= gr->nodes; nptr-- )
194  {
195  if ( nptr->unmarked && nptr->alive )
196  {
197  nptr->dist = (int) n;
198  nptr->alive = FALSE;
199  }
200  }
201  bound = count;
202  }
203 
204 } /* global_relabel */
205 
206 
207 double maxflow (GRAPH *gr, GRAPHNODE *s_ptr, GRAPHNODE *t_ptr)
208 {
209  /* Determines maximum flow and minimum cut between nodes
210  s (= *s_ptr) and t (= *t_ptr) in an undirected graph
211 
212  References:
213  ----------
214  A. Goldberg/ E. Tarjan: "A New Approach to the
215  Maximum Flow Problem", in Proc. 18th ACM Symp.
216  on Theory of Computing, 1986.
217  */
218  GRAPHNODE *aptr, *nptr, *q_front, *q_rear;
219  GRAPHEDGE *eptr;
220  long n, m, m0, level, i, n_discharge;
221  double incre;
222  long dmin;
223  double cap;
224 
225  /* node ids range from 1 to n, node array indices
226  range from 0 to n-1 */
227 
228  n = gr->nnodes;
229  for ( nptr = &(gr->nodes[n-1L]); nptr >= gr->nodes; nptr-- )
230  {
231  nptr->scan_ptr = nptr->first_edge;
232  while( nptr->scan_ptr != NULL && nptr->scan_ptr->cap <= EPS )
233  nptr->scan_ptr = nptr->scan_ptr->next;
234 
235  if ( nptr->scan_ptr == NULL )
236  {
237  fprintf (stderr, "isolated node in input graph\n");
238  return (FALSE);
239  }
240 
241  nptr->excess = 0.0L;
242  nptr->stack_link = NULL;
243  nptr->alive = TRUE;
244  nptr->unmarked = TRUE;
245  }
246 
247  m = gr->nedgesnonzero;
248  m0 = gr->nedges;
249  for ( eptr = &(gr->edges[m-1L]); eptr >= gr->edges; eptr-- )
250  eptr->rcap = eptr->cap;
251  for ( eptr = &(gr->edges[m0+m-1L]); eptr >= &(gr->edges[m0]); eptr-- )
252  eptr->rcap = eptr->cap;
253 
254  for ( i = n; i >= 0L; i-- )
255  {
256  number[i] = 0L;
257  active[i] = NULL;
258  }
259  t_ptr->dist = 0L;
260 
261  /* breadth first search to get exact distances
262  from sink and for test of graph connectivity */
263 
264  t_ptr->unmarked = FALSE;
265  q_front = t_ptr;
266  q_rear = q_front;
267  bfs_next:
268  level = q_rear->dist + 1L;
269  eptr = q_rear->first_edge;
270  while ( eptr != NULL )
271  {
272  assert(eptr->back->cap == eptr->back->rcap); /*lint !e777*/
273  if ( eptr->adjac->unmarked && eptr->back->rcap > EPS )
274  {
275  nptr = eptr->adjac;
276  nptr->unmarked = FALSE;
277  nptr->dist = (int) level;
278  ++number[level];
279  q_front->bfs_link = nptr;
280  q_front = nptr;
281  }
282  eptr = eptr->next;
283  }
284  if ( q_rear == q_front )
285  goto bfs_ready;
286 
287  q_rear = q_rear->bfs_link;
288  goto bfs_next;
289 
290  bfs_ready:
291  if ( co_check )
292  {
293  co_check = FALSE;
294  for ( nptr = &(gr->nodes[n-1]); nptr >= gr->nodes; --nptr )
295  {
296  if ( nptr->unmarked )
297  return (-1.0L);
298  }
299  }
300 
301  s_ptr->dist = (int) n; /* number[0] and number[n] not required */
302  t_ptr->dist = 0L;
303  t_ptr->excess = 1.0L; /* to be subtracted again */
304 
305 
306  /* initial preflow push from source node */
307 
308  max_dist = 0L; /* = max_dist of active nodes */
309  eptr = s_ptr->first_edge;
310  while ( eptr != NULL )
311  {
312  if( eptr->cap > EPS )
313  {
314  nptr = eptr->adjac;
315  cap = eptr->rcap;
316  nptr->excess += cap;
317  s_ptr->excess -= cap;
318  eptr->back->rcap += cap;
319  eptr->rcap = 0.0L;
320 
321  if ( nptr != t_ptr && nptr->excess <= cap + EPS )
322  {
323  /* push node nptr onto stack for nptr->dist,
324  but only once in case of double edges */
325  nptr->stack_link = active[nptr->dist];
326  active[nptr->dist] = nptr;
327  if ( nptr->dist > max_dist )
328  max_dist = nptr->dist;
329  }
330  }
331  eptr = eptr->next;
332  }
333 
334  s_ptr->alive = FALSE;
335  bound = n;
336  n_discharge = 0L;
337 
338  /* main loop */
339 
340  do
341  {
342  /* get maximum distance active node */
343  aptr = active[max_dist];
344  while ( aptr != NULL )
345  {
346  active[max_dist] = aptr->stack_link;
347  eptr = aptr->scan_ptr;
348  assert(eptr != NULL);
349  assert(eptr->back != NULL);
350  assert(eptr->cap > EPS);
351 
352  edge_scan: /* for current active node */
353  assert(eptr != NULL);
354  assert(eptr->back != NULL);
355  assert(eptr->cap > EPS);
356 
357  nptr = eptr->adjac;
358  assert(nptr != NULL);
359  if ( nptr->dist == aptr->dist - 1L && eptr->rcap > EPS )
360  {
361  incre = aptr->excess;
362  if ( incre <= eptr->rcap )
363  {
364  /* perform a non saturating push */
365  eptr->rcap -= incre;
366  eptr->back->rcap += incre;
367  aptr->excess = 0.0L;
368  nptr->excess += incre;
369  if ( nptr->excess <= incre + EPS )
370  {
371  /* push nptr onto active stack */
372  nptr->stack_link = active[nptr->dist];
373  active[nptr->dist] = nptr;
374  }
375  assert(eptr->cap > EPS);
376  aptr->scan_ptr = eptr;
377  goto node_ready;
378  }
379  else
380  {
381  /* perform a saturating push */
382  incre = eptr->rcap;
383  eptr->back->rcap += incre;
384  aptr->excess -= incre;
385  nptr->excess += incre;
386  eptr->rcap = 0.0L;
387  if ( nptr->excess <= incre + EPS )
388  {
389  /* push nptr onto active stack */
390  nptr->stack_link = active[nptr->dist];
391  active[nptr->dist] = nptr;
392  }
393  if ( aptr->excess <= EPS )
394  {
395  assert(eptr->cap > EPS);
396  aptr->scan_ptr = eptr;
397  goto node_ready;
398  }
399  }
400  }
401 
402  /* go to the next non-zero edge */
403  do
404  {
405  eptr = eptr->next;
406  }
407  while( eptr != NULL && eptr->cap <= EPS );
408 
409  if ( eptr == NULL )
410  {
411  /* all incident arcs scanned, but node still
412  has positive excess, check if for all nptr
413  nptr->dist != aptr->dist */
414 
415  if ( number[aptr->dist] == 1L )
416  {
417  /* put all nodes v with dist[v] >= dist[a]
418  into the set of "dead" nodes since they
419  are disconnected from the sink */
420 
421  for ( nptr = &(gr->nodes[n-1L]); nptr >= gr->nodes; nptr-- )
422  {
423  if ( nptr->alive && nptr->dist > aptr->dist )
424  {
425  --number[nptr->dist];
426  active[nptr->dist] = NULL;
427  nptr->alive = FALSE;
428  nptr->dist = (int) n;
429  --bound;
430  }
431  }
432  --number[aptr->dist];
433  active[aptr->dist] = NULL;
434  aptr->alive = FALSE;
435  aptr->dist = (int) n;
436  --bound;
437  goto node_ready;
438  }
439  else
440  {
441  /* determine new label value */
442  dmin = n;
443  aptr->scan_ptr = NULL;
444  eptr = aptr->first_edge;
445  while ( eptr != NULL )
446  {
447  assert(eptr->rcap <= EPS || eptr->cap > EPS);
448  if ( eptr->adjac->dist < dmin && eptr->rcap > EPS )
449  {
450  assert(eptr->cap > EPS);
451  dmin = eptr->adjac->dist;
452  if ( aptr->scan_ptr == NULL )
453  aptr->scan_ptr = eptr;
454  }
455  eptr = eptr->next;
456  }
457 
458  if ( ++dmin < bound )
459  {
460  /* ordinary relabel operation */
461  --number[aptr->dist];
462  aptr->dist = (int) dmin;
463  ++number[dmin];
464  max_dist = dmin;
465  eptr = aptr->scan_ptr;
466  assert(eptr != NULL);
467  assert(eptr->cap > EPS);
468  goto edge_scan;
469  }
470  else
471  {
472  aptr->alive = FALSE;
473  --number[aptr->dist];
474  aptr->dist = (int) n;
475  --bound;
476  goto node_ready;
477  }
478  }
479  }
480  else
481  goto edge_scan;
482 
483 
484  node_ready:
485  ++n_discharge;
486  if ( n_discharge == n )
487  {
488  n_discharge = 0L;
489  global_relabel (gr, t_ptr);
490  }
491  aptr = active[max_dist];
492  } /* aptr != NULL */
493  --max_dist;
494  }
495  while ( max_dist > 0L );
496 
497  return (int) (t_ptr->excess - 1.0L);
498 }
499 
500 SCIP_Bool nodeOnRootPath(GRAPH* gr, int i, int j)
501 {
502  while( i != j && j != 0 )
503  {
504  j = gr->nodes[j].parent->id;
505  }
506  if( i == j )
507  return TRUE;
508  else
509  return FALSE;
510 }
511 
512 // constructs a list of cuts for a TSP relaxation polytope from a Gomory Hu Tree
513 void constructCutList(GRAPH *gr, SCIP_Bool** cuts, int* ncuts, double minviol)
514 {
515  int k = 0;
516  for( int i = 1; i < gr->nnodes; i++ )
517  {
518  if( gr->nodes[i].mincap < 2.0 - minviol )
519  {
520  cuts[k][0] = FALSE;
521  for( int j = 1 ; j < gr->nnodes; j++ )
522  cuts[k][j] = nodeOnRootPath(gr, i, j);
523  k++;
524  }
525  }
526  *ncuts = k;
527 }
528 
529 // if the non-zero-edges of a TSP relaxation induce a non-connected graph,
530 // an according cut is generated, using information from BFS in method maxflow
532 {
533  cuts[0][0] = FALSE;
534  for( int i = 1; i < gr->nnodes; i++ )
535  cuts[0][i]=gr->nodes[i].unmarked;
536 }
537 
538 SCIP_Bool ghc_tree (GRAPH *gr, SCIP_Bool** cuts, int* ncuts, double minviol)
539 {
540  /* Determines Gomory/Hu cut tree for input graph with
541  capacitated edges, the tree structures is represented
542  by parent pointers which are part of the node structure,
543  the capacity of a tree edge is stored at the child node,
544  the root of the cut tree is the first node in the list
545  of graph nodes (&gr->nodes[0]). The implementation is
546  described in [1].
547 
548  References:
549  ----------
550  1) D. Gusfield: "Very Simple Algorithms and Programs for
551  All Pairs Network Flow Analysis", Computer Science Di-
552  vision, University of California, Davis, 1987.
553 
554  2) R.E. Gomory and T.C. Hu: "Multi-Terminal Network Flows",
555  SIAM J. Applied Math. 9 (1961), 551-570.
556 
557  */
558 
559  GRAPHNODE *nptr, *nptr1, *nptrn, *sptr, *tptr;
560  long n;
561  double maxfl;
562 
563  n = gr->nnodes;
564 
565  if ( ! init_maxflow (n) )
566  return FALSE;
567 
568  nptr1 = gr->nodes;
569  nptrn = &(gr->nodes[n-1L]);
570  for ( nptr = nptrn; nptr >= nptr1; nptr-- )
571  nptr->parent = nptr1;
572 
573  for ( sptr = &(gr->nodes[1L]); sptr <= nptrn; sptr++ )
574  {
575  tptr = sptr->parent;
576  maxfl = maxflow(gr, sptr, tptr);
577 
578  //maxfl < 0 <=> graph is not connected => generate cut
579  if ( maxfl < 0L )
580  {
581  constructSingleCut(gr,cuts);
582  *ncuts = 1;
583  fini_maxflow ();
584  return TRUE;
585  }
586 
587  sptr->mincap = maxfl;
588  for ( nptr = &(gr->nodes[1L]); nptr <= nptrn; nptr++ )
589  if ( nptr != sptr && ! nptr->alive && nptr->parent == tptr )
590  nptr->parent = sptr;
591  if ( ! tptr->parent->alive )
592  {
593  sptr->parent = tptr->parent;
594  tptr->parent = sptr;
595  sptr->mincap = tptr->mincap;
596  tptr->mincap = maxfl;
597  }
598  }
599  fini_maxflow ();
600  constructCutList(gr,cuts,ncuts,minviol);
601 
602  return TRUE;
603 } /* ghc_tree */
struct GraphEdge * next
Definition: GomoryHuTree.h:59
#define NULL
Definition: def.h:239
Definition: grph.h:57
SCIP_Bool create_graph(int n, int m, GRAPH **gr)
static long bound
#define FALSE
Definition: def.h:65
#define TRUE
Definition: def.h:64
#define BMSallocMemoryArray(ptr, num)
Definition: memory.h:105
static GRAPHNODE ** active
generator for global cuts in undirected graphs
#define BMSfreeMemory(ptr)
Definition: memory.h:127
double maxflow(GRAPH *gr, GRAPHNODE *s_ptr, GRAPHNODE *t_ptr)
void global_relabel(GRAPH *gr, GRAPHNODE *tptr)
static SCIP_Bool co_check
struct GraphEdge * scan_ptr
Definition: GomoryHuTree.h:44
#define BMSfreeMemoryArray(ptr)
Definition: memory.h:129
#define EPS
double cap
Definition: GomoryHuTree.h:55
GRAPHNODE * adjac
Definition: GomoryHuTree.h:62
struct GraphEdge * back
Definition: GomoryHuTree.h:60
struct GraphNode * parent
Definition: GomoryHuTree.h:48
struct GraphEdge * first_edge
Definition: GomoryHuTree.h:43
C++ wrapper classes for SCIP.
SCIP_Bool unmarked
Definition: GomoryHuTree.h:40
double excess
Definition: GomoryHuTree.h:37
void fini_maxflow()
struct GraphNode * stack_link
Definition: GomoryHuTree.h:47
#define SCIP_Bool
Definition: def.h:62
double rcap
Definition: GomoryHuTree.h:56
SCIP_Bool alive
Definition: GomoryHuTree.h:41
static long max_dist
static long * number
double mincap
Definition: GomoryHuTree.h:38
SCIP_Bool ghc_tree(GRAPH *gr, SCIP_Bool **cuts, int *ncuts, double minviol)
void constructCutList(GRAPH *gr, SCIP_Bool **cuts, int *ncuts, double minviol)
static void free_graph(GRAPH **gr)
#define BMSallocMemory(ptr)
Definition: memory.h:101
void constructSingleCut(GRAPH *gr, SCIP_Bool **cuts)
SCIP_Bool init_maxflow(long n)
int edges
Definition: grph.h:92
struct GraphNode * bfs_link
Definition: GomoryHuTree.h:46
void capture_graph(GRAPH *gr)
void release_graph(GRAPH **gr)
SCIP_Bool nodeOnRootPath(GRAPH *gr, int i, int j)