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