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
grbrweight(Gredge_t * e,ssize_t w)261 ssize_t grbrweight(Gredge_t* e, ssize_t w)
262 {
263 if(w > 0)
264 BREDGE(e)->wght = w;
265 return BREDGE(e)->wght;
266 }
267
268
brdata(Gralgo_t * algo,int type,Grdata_t * data)269 static Grdata_t* brdata(Gralgo_t* algo, int type, Grdata_t* data)
270 {
271 Brnode_t* brn;
272 Bredge_t* bre;
273
274 if(algo != Grbranching)
275 return NIL(Grdata_t*);
276
277 switch(type)
278 { default:
279 return NIL(Grdata_t*);
280
281 case GR_NODE|GR_CLOSING:
282 case GR_EDGE|GR_CLOSING:
283 if(data)
284 free(data);
285 return NIL(Grdata_t*);
286
287 case GR_NODE|GR_OPENING:
288 if(!(brn = (Brnode_t*)calloc(1, sizeof(Brnode_t))) )
289 return NIL(Grdata_t*);
290 return (Grdata_t*)brn;
291
292 case GR_EDGE|GR_OPENING:
293 if(!(bre = (Bredge_t*)calloc(1, sizeof(Bredge_t))) )
294 return NIL(Grdata_t*);
295 return (Grdata_t*)bre;
296 }
297 }
298
299 static Gralgo_t _Grbranching = { brdata, 0 };
300 Gralgo_t* Grbranching = &_Grbranching;
301
302 #ifdef KPVTEST
303 #include <stdio.h>
main(int argc,char ** argv)304 main(int argc, char** argv)
305 {
306 Graph_t *gr;
307 Grnode_t *n, *hn, *tn;
308 Gredge_t *e;
309 int h, t, w;
310 char buf[1024];
311 int type = 0;
312
313 if(argc > 1 && strcmp(argv[1],"-l") == 0)
314 type = -1;
315 if(argc > 1 && strcmp(argv[1],"-r") == 0)
316 type = 1;
317
318 gr = gropen(0);
319 while(fgets(buf,sizeof(buf),stdin))
320 { if(buf[0] == '#')
321 continue;
322 if(sscanf(buf,"%d,%d,%d", &t, &h, &w) != 3)
323 continue;
324 tn = grnode(gr, t, 1);
325 hn = grnode(gr, h, 1);
326 if(type > 0 && t >= h) /* only use edges with t < h */
327 continue;
328 if(type < 0 && t <= h) /* only use edges with t > h */
329 continue;
330 e = gredge(gr, tn, hn, 0, 1);
331 grbrweight(e, w);
332 }
333
334 printf("\nTotal weight=%d\n", grbranching(gr));
335 for(n = (Grnode_t*)dtflatten(gr->nodes); n; n = (Grnode_t*)dtlink(gr->nodes, n) )
336 if((e = n->iedge) )
337 fprintf(stderr, "%d -> %d [%d]\n", e->tail->label, e->head->label, grbrweight(e, 0) );
338 }
339 #endif
340