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