1 #include "grhdr.h"
2
3 /******************************************************************************
4 Chu-Liu-Edmonds algorithm to compute a maximum branching for a
5 directed graph with edge weights. The algorithm alters the edge
6 lists to contain just the branching edges.
7
8 Written by Kiem-Phong Vo (02/01/2006)
9 ******************************************************************************/
10
11 typedef struct _brcycl_s
12 { Gredge_t* cycl; /* the cycle of strong components */
13 Gredge_t* emin; /* the minimum weight edge of cycle */
14 Gredge_t* entr; /* the edge that came in to the cycle */
15 } Brcycl_t;
16
17 typedef struct _brnode_s
18 { Grdata_t data;
19 ssize_t mark; /* mark 1 if already searched */
20 } Brnode_t;
21
22 typedef struct _bredge_s
23 { Grdata_t data;
24 ssize_t wght; /* the original edge weight */
25 ssize_t wadj; /* adjusted weight during construction */
26 Grnode_t* root; /* strong component before collapsing */
27 Gredge_t* edge; /* branching edge that this represents */
28 } Bredge_t;
29
30 #define BRNODE(n) ((Brnode_t*)grdtnode((n), Grbranching) )
31 #define BREDGE(e) ((Bredge_t*)grdtedge((e), Grbranching) )
32
33 /* compute the representative node in the shadowed union structure */
34 #define GRLINK(n) ((n)->link == (n) ? (n) : grlink(n))
grlink(Grnode_t * n)35 static Grnode_t* grlink(Grnode_t* n)
36 { while(n->link != n)
37 n = n->link;
38 return n;
39 }
40
41 #ifdef DEBUG
42 static int Fd = 2;
predge(Gredge_t * ed)43 static int predge(Gredge_t* ed)
44 {
45 if(!ed)
46 { PRINT(Fd,"Null edge\n",0);
47 return 0;
48 }
49 PRINT(Fd,"%d", (int)ed->tail->label); PRINT(Fd,"(%d) -> ", (int)grfind(ed->tail)->label );
50 PRINT(Fd,"%d", (int)ed->head->label); PRINT(Fd,"(%d), ", (int)grfind(ed->head)->label );
51 PRINT(Fd,"root=%d, ", BREDGE(ed)->root ? (int)BREDGE(ed)->root->label : -1 );
52 PRINT(Fd,"wght=%d, ", (int)BREDGE(ed)->wght); PRINT(Fd,"wadj=%d\n", (int)BREDGE(ed)->wadj);
53 return 0;
54 }
prlink(Gredge_t * e)55 static int prlink(Gredge_t* e)
56 { for(; e; e = e->link)
57 { PRINT(Fd, "\t", 0); predge(e); }
58 PRINT(Fd,"\n",0);
59 }
prinext(Gredge_t * e)60 static int prinext(Gredge_t* e)
61 { for(; e; e = e->inext)
62 { PRINT(Fd, "\t", 0); predge(e); }
63 PRINT(Fd,"\n",0);
64 }
pronext(Gredge_t * e)65 static int pronext(Gredge_t* e)
66 { for(; e; e = e->onext)
67 { PRINT(Fd, "\t", 0); predge(e); }
68 PRINT(Fd,"\n",0);
69 }
prnode(Grnode_t * nd)70 static int prnode(Grnode_t* nd)
71 {
72 if(!nd)
73 { PRINT(Fd,"Null node\n",0);
74 return 0;
75 }
76 PRINT(Fd,"node = %d", (int)nd->label); PRINT(Fd,"(%d), ", grfind(nd)->label);
77 PRINT(Fd,"link = %d, ", nd->link ? (int)nd->link->label : -1);
78 PRINT(Fd,"mark = %d\n", (int)BRNODE(nd)->mark);
79 return 0;
80 }
prid(Graph_t * gr,ssize_t id)81 static int prid(Graph_t* gr, ssize_t id)
82 {
83 return prnode(grnode(gr, (void*)id, 0));
84 }
prcycl(Brcycl_t * cl)85 static int prcycl(Brcycl_t* cl)
86 {
87 Gredge_t *e;
88 PRINT(Fd,"Emin: ",0); predge(cl->emin);
89 PRINT(Fd,"Entr: ",0); predge(cl->entr);
90 prlink(cl->cycl);
91 return 0;
92 }
prgraph(Graph_t * gr)93 static int prgraph(Graph_t* gr)
94 { Grnode_t *nd;
95 Gredge_t *ed;
96 for(nd = (Grnode_t*)dtflatten(gr->nodes); nd; nd = (Grnode_t*)dtlink(gr->nodes,nd) )
97 { if(nd != grfind(nd))
98 continue;
99 prnode(nd);
100 for(ed = nd->iedge; ed; ed = ed->inext)
101 { PRINT(Fd,"\t",0); predge(ed); }
102 }
103 }
104 #endif
105
grbranching(Graph_t * gr)106 ssize_t grbranching(Graph_t* gr)
107 {
108 Grnode_t *n, *nc;
109 Gredge_t *e, *ec, *ep, *en, *path, *emin;
110 Brcycl_t *clist, *cl, *endcl;
111 ssize_t w;
112
113 for(w = 0, n = (Grnode_t*)dtflatten(gr->nodes); n; w += 1, n = (Grnode_t*)dtlink(gr->nodes,n) )
114 { n->link = n->fold = n; /* union structures: link kept as-is, fold does path-compression */
115 n->oedge = NIL(Gredge_t*); /* wipe the out-edges */
116 BRNODE(n)->mark = 0;
117
118 /* compute heaviest weight ec and move it to front */
119 for(ep = en = NIL(Gredge_t*), ec = e = n->iedge; e; en = e, e = e->inext )
120 { e->link = NIL(Gredge_t*);
121 BREDGE(e)->edge = e; /* at start, representing self */
122 if((BREDGE(e)->wadj = BREDGE(e)->wght) > BREDGE(ec)->wadj )
123 { ec = e; ep = en; }
124 }
125 if(ep)
126 { ep->inext = ec->inext; ec->inext = n->iedge; n->iedge = ec; }
127 }
128 if(w == 0)
129 return 0;
130
131 /* space to keep cycle structures */
132 if(!(clist = cl = (Brcycl_t*)calloc(w+1,sizeof(Brcycl_t))) )
133 return -1;
134 endcl = cl+w+1;
135
136 /* search and collapse cycles */
137 for(n = (Grnode_t*)dtflatten(gr->nodes); n; n = (Grnode_t*)dtlink(gr->nodes,n) )
138 { nc = grfind(n);
139 if(BRNODE(nc)->mark) /* already searched */
140 continue;
141
142 for(path = NIL(Gredge_t*); nc; ) /* depth-first search */
143 { if(!BRNODE(nc)->mark) /* not searched yet */
144 { BRNODE(nc)->mark = 1;
145
146 if(!(ec = nc->iedge) ) /* this path cannot be a cycle */
147 break;
148 else /* continue to build the path */
149 { ec->link = path; path = ec;
150 nc = grfind(ec->tail);
151 continue; /* dfs recursion */
152 }
153 }
154
155 /* potential cycle, check that out and also compute min edge */
156 for(emin = NIL(Gredge_t*), ec = path; ec; ec = ec->link)
157 { if(!emin || BREDGE(ec)->wadj < BREDGE(emin)->wadj )
158 emin = ec;
159 if(grfind(ec->head) == nc) /* end of cycle */
160 break;
161 }
162 if(!ec) /* run up against an old path */
163 break;
164
165 if(cl >= endcl) /* hmm, this should not happen */
166 { /**/ASSERT(cl < endcl);
167 free(clist);
168 return -1;
169 }
170
171 cl->cycl = path;
172 path = ec->link; /* will restart search from here */
173 ec->link = NIL(Gredge_t*);
174
175 /* make list of incoming edges and adjust their weights */
176 en = NIL(Gredge_t*);
177 for(ec = cl->cycl; ec; ec = ec->link)
178 { BREDGE(ec)->root = grfind(ec->head); /* save node union structure */
179
180 w = BREDGE(ec)->wadj - BREDGE(emin)->wadj;
181 for(ep = NIL(Gredge_t*), e = grfind(ec->head)->iedge; e; e = e->inext)
182 { BREDGE(e)->wadj -= w;
183 if(!e->inext) /* last of list */
184 ep = e;
185 }
186 ep->inext = en; en = grfind(ec->head)->iedge; /* catenate lists */
187 grfind(ec->head)->iedge = NIL(Gredge_t*);
188 }
189
190 /* collapsing cycle onto nc */
191 for(ec = cl->cycl; ec; ec = ec->link)
192 { if(grfind(ec->head) == nc)
193 continue;
194 grfind(ec->head)->link = nc; /* union history kept as-is */
195 grfind(ec->head)->fold = nc; /* union with path-compression */
196 }
197 nc->fold = nc->link = nc;
198
199 /* make new edge list, keep heaviest edge in front */
200 for(e = en; e; e = en)
201 { en = e->inext;
202 if(BREDGE(e)->wadj <= 0 || grfind(e->head) == grfind(e->tail) )
203 continue;
204 if(!nc->iedge || BREDGE(nc->iedge)->wadj <= BREDGE(e)->wadj )
205 { e->inext = nc->iedge; nc->iedge = e; }
206 else { e->inext = nc->iedge->inext; nc->iedge->inext = e; }
207 }
208
209 cl->emin = emin;
210 cl->entr = nc->iedge;
211 cl += 1;
212 BRNODE(nc)->mark = 0; /* force nc to be searched again */
213 }
214 }
215
216 /* move the remaining branching edges to their real nodes */
217 path = NIL(Gredge_t*);
218 for(n = (Grnode_t*)dtflatten(gr->nodes); n; n = (Grnode_t*)dtlink(gr->nodes,n) )
219 { if(!(e = n->iedge) )
220 continue;
221 e->link = path; path = e;
222 }
223 for(ec = path; ec; ec = ec->link)
224 ec->head->iedge = ec;
225
226 /* unroll collapsed cycles in reverse order to construct branching */
227 for(cl -= 1; cl >= clist; --cl )
228 { /* restore the union structure to just before this cycle collapsed */
229 for(ec = cl->cycl; ec; ec = ec->link)
230 BREDGE(ec)->root->link = BREDGE(ec)->root;
231
232 if((en = cl->entr) )
233 en = BREDGE(en)->edge;
234 for(ec = cl->cycl; ec; ec = ec->link)
235 { if(en && GRLINK(ec->head) == GRLINK(en->head))
236 BREDGE(ec)->edge = en;
237 else ec->head->iedge = ec;
238 }
239 if(!en) /* isolated cycle, remove minimum edge */
240 { cl->emin->head->iedge = NIL(Gredge_t*);
241 if(BREDGE(cl->emin)->edge == cl->emin)
242 BREDGE(cl->emin)->edge = NIL(Gredge_t*);
243 }
244 }
245
246 w = 0; /* construct the external branching representation */
247 for(n = (Grnode_t*)dtflatten(gr->nodes); n; n = (Grnode_t*)dtlink(gr->nodes,n) )
248 { if(!(e = n->iedge) )
249 continue;
250 e->inext = NIL(Gredge_t*);
251 e->onext = e->tail->oedge; e->tail->oedge = e;
252 w += BREDGE(e)->wght;
253 }
254
255 free(clist);
256
257 return w;
258 }
259
260 /* sort edges in reverse order of weights */
gredgesort(Gredge_t * list)261 static Gredge_t* gredgesort(Gredge_t* list)
262 {
263 Gredge_t *link, *equl, *less, *more;
264 ssize_t w, wght;
265
266 if(!list)
267 return NIL(Gredge_t*);
268
269 equl = list; list = list->link; equl->link = NIL(Gredge_t*);
270 wght = BREDGE(equl)->wght; /* partition list by this weight */
271 more = less = NIL(Gredge_t*);
272 for(; list; list = link)
273 { link = list->link;
274
275 if((w = BREDGE(list)->wght) > wght)
276 { list->link = more;
277 more = list;
278 }
279 else if(w == wght)
280 { list->link = equl;
281 equl = list;
282 }
283 else /* if(w < wght) */
284 { list->link = less;
285 less = list;
286 }
287 }
288
289 /* recurse and sort the sublists */
290 if(more && more->link)
291 more = gredgesort(more);
292 if(less && less->link)
293 less = gredgesort(less);
294
295 if((list = more) ) /* heaviest ones go first */
296 { for(link = list; link->link; )
297 link = link->link;
298 link->link = equl; /* link to the equals */
299 }
300 else list = equl; /* no heavier ones than this */
301
302 for(link = equl; link->link; )
303 link = link->link;
304 link->link = less; /* lighter ones go last */
305
306 return list;
307 }
308
309 /* Greedy approximation of a branching */
grbrgreedy(Graph_t * gr)310 ssize_t grbrgreedy(Graph_t* gr)
311 {
312 Grnode_t *n;
313 Gredge_t *e, *list;
314 ssize_t wght;
315
316 list = NIL(Gredge_t*); /* link all edges into a big list */
317 for(n = (Grnode_t*)dtflatten(gr->nodes); n; n = (Grnode_t*)dtlink(gr->nodes,n) )
318 { for(e = n->oedge; e; e = e->onext)
319 { e->link = list;
320 list = e;
321 }
322 n->oedge = n->iedge = NIL(Gredge_t*); /* wipe out edge lists */
323 }
324 list = gredgesort(list); /* sort in reverse order by weights */
325
326 wght = 0;
327 for(e = list; e; e = e->link)
328 { if(e->head == e->tail) /* self-loop cannot be a branching edge */
329 continue;
330 if(e->head->iedge) /* node already got an incoming edge */
331 continue;
332
333 for(n = e->tail;; ) /* check cycle */
334 { if(!n->iedge) /* no cycle */
335 break;
336
337 if((n = n->iedge->tail) == e->head) /* causing a cycle */
338 { n = NIL(Grnode_t*);
339 break;
340 }
341 }
342
343 if(n) /* ok to add to the branching */
344 { e->inext = NIL(Gredge_t*); e->head->iedge = e;
345 e->onext = e->tail->oedge; e->tail->oedge = e;
346 wght += BREDGE(e)->wght;
347 }
348 }
349
350 return wght;
351 }
352
353
354 /* set or query weights of an edge */
grbrweight(Gredge_t * e,ssize_t w)355 ssize_t grbrweight(Gredge_t* e, ssize_t w)
356 {
357 if(w > 0)
358 BREDGE(e)->wght = w;
359 return BREDGE(e)->wght;
360 }
361
362
brdata(Gralgo_t * algo,int type,Grdata_t * data)363 static Grdata_t* brdata(Gralgo_t* algo, int type, Grdata_t* data)
364 {
365 Brnode_t* brn;
366 Bredge_t* bre;
367
368 if(algo != Grbranching)
369 return NIL(Grdata_t*);
370
371 switch(type)
372 { default:
373 return NIL(Grdata_t*);
374
375 case GR_NODE|GR_CLOSING:
376 case GR_EDGE|GR_CLOSING:
377 if(data)
378 free(data);
379 return NIL(Grdata_t*);
380
381 case GR_NODE|GR_OPENING:
382 if(!(brn = (Brnode_t*)calloc(1, sizeof(Brnode_t))) )
383 return NIL(Grdata_t*);
384 return (Grdata_t*)brn;
385
386 case GR_EDGE|GR_OPENING:
387 if(!(bre = (Bredge_t*)calloc(1, sizeof(Bredge_t))) )
388 return NIL(Grdata_t*);
389 return (Grdata_t*)bre;
390 }
391 }
392
393 static Gralgo_t _Grbranching = { brdata, 0 };
394 Gralgo_t* Grbranching = &_Grbranching;
395
396 #ifdef KPVTEST
397 #include <stdio.h>
main(int argc,char ** argv)398 main(int argc, char** argv)
399 {
400 Graph_t *gr;
401 Grnode_t *n, *hn, *tn;
402 Gredge_t *e;
403 int h, t, w;
404 char buf[1024];
405 int type = 0;
406
407 if(argc > 1 && strcmp(argv[1],"-l") == 0)
408 { type = -1;
409 argc--; argv++;
410 }
411 if(argc > 1 && strcmp(argv[1],"-r") == 0)
412 { type = 1;
413 argc--; argv++;
414 }
415
416 gr = gropen(0,0);
417 while(fgets(buf,sizeof(buf),stdin))
418 { if(buf[0] == '#')
419 continue;
420 if(sscanf(buf,"%d,%d,%d", &t, &h, &w) != 3)
421 continue;
422 tn = grnode(gr, (Void_t*)(t), 1);
423 hn = grnode(gr, (Void_t*)(h), 1);
424 if(type > 0 && t >= h) /* only use edges with t < h */
425 continue;
426 if(type < 0 && t <= h) /* only use edges with t > h */
427 continue;
428 e = gredge(gr, tn, hn, 0, 1);
429 grbrweight(e, w);
430 }
431
432 if(argc > 1 && strcmp(argv[1],"-g") == 0)
433 printf("\nTotal weight=%d\n", grbrgreedy(gr));
434 else printf("\nTotal weight=%d\n", grbranching(gr));
435 for(n = (Grnode_t*)dtflatten(gr->nodes); n; n = (Grnode_t*)dtlink(gr->nodes, n) )
436 if((e = n->iedge) )
437 fprintf(stderr, "%d -> %d [%d]\n", e->tail->label, e->head->label, grbrweight(e, 0) );
438 }
439 #endif
440