1/*************************************************************************
2 * Copyright (c) 2011 AT&T Intellectual Property
3 * All rights reserved. This program and the accompanying materials
4 * are made available under the terms of the Eclipse Public License v1.0
5 * which accompanies this distribution, and is available at
6 * http://www.eclipse.org/legal/epl-v10.html
7 *
8 * Contributors: See CVS logs. Details at http://www.graphviz.org/
9 *************************************************************************/
10
11#include "config.h"
12
13#include "types.h"
14#include "globals.h"
15#include "general.h"
16#include "SparseMatrix.h"
17#include "edge_bundling.h"
18#include <time.h>
19#include "clustering.h"
20#include "ink.h"
21#include "agglomerative_bundling.h"
22
23#define SMALL 1.e-10
24
25#ifdef OPENGL
26#include "gl.h"
27extern pedge *edges_global;
28extern int *clusters_global;
29#endif
30
31static real norm(int n, real *x){
32 real res = 0;
33 int i;
34
35 for (i = 0; i < n; i++) res += x[i]*x[i];
36 return sqrt(res);
37
38}
39static real sqr_dist(int dim, real *x, real *y){
40 int i;
41 real res = 0;
42 for (i = 0; i < dim; i++) res += (x[i] - y[i])*(x[i] - y[i]);
43 return res;
44}
45static real dist(int dim, real *x, real *y){
46 return sqrt(sqr_dist(dim,x,y));
47}
48
49pedge pedge_new(int np, int dim, real *x){
50 pedge e;
51
52 e = MALLOC(sizeof(struct pedge_struct));
53 e->npoints = np;
54 e->dim = dim;
55 e->len = np;
56 e->x = MALLOC(dim*(e->len)*sizeof(real));
57 MEMCPY(e->x, x, dim*(e->len)*sizeof(real));
58 e->edge_length = dist(dim, &(x[0*dim]), &(x[(np-1)*dim]));
59 e->wgt = 1.;
60 e->wgts = NULL;
61 return e;
62
63}
64pedge pedge_wgt_new(int np, int dim, real *x, real wgt){
65 pedge e;
66 int i;
67
68 e = MALLOC(sizeof(struct pedge_struct));
69 e->npoints = np;
70 e->dim = dim;
71 e->len = np;
72 e->x = MALLOC(dim*(e->len)*sizeof(real));
73 MEMCPY(e->x, x, dim*(e->len)*sizeof(real));
74 e->edge_length = dist(dim, &(x[0*dim]), &(x[(np-1)*dim]));
75 e->wgt = wgt;
76 e->wgts = MALLOC(sizeof(real)*(np - 1));
77 for (i = 0; i < np - 1; i++) e->wgts[i] = wgt;
78 return e;
79
80}
81void pedge_delete(pedge e){
82 FREE(e->x);
83 FREE(e);
84}
85
86pedge pedge_flip(pedge e){
87 /* flip the polyline so that last point becomes the first, second last the second, etc*/
88 real *y;
89 real *x = e->x;
90 int i, dim = e->dim;
91 int n = e->npoints;
92
93 y = MALLOC(sizeof(real)*e->dim);
94 for (i = 0; i < (e->npoints)/2; i++){
95 MEMCPY(y, &(x[i*dim]), sizeof(real)*dim);
96 MEMCPY(&(x[(n-1-i)*dim]), &(x[i*dim]), sizeof(real)*dim);
97 MEMCPY(&(x[i*dim]), y, sizeof(real)*dim);
98 }
99 FREE(y);
100 return e;
101}
102
103static real edge_compatibility(pedge e1, pedge e2){
104 /* two edges are u1->v1, u2->v2.
105 return 1 if two edges are exactly the same, 0 if they are very different.
106 */
107 real *u1, *v1, *u2, *v2, *u, dist1, dist2, len1, len2;
108 int dim = e1->dim, flipped = FALSE;
109
110 u1 = e1->x;
111 v1 = (e1->x)+(e1->npoints)*dim-dim;
112 u2 = e2->x;
113 v2 = (e2->x)+(e2->npoints)*dim-dim;
114 dist1 = sqr_dist(dim, u1, u2) + sqr_dist(dim, v1, v2);
115 dist2 = sqr_dist(dim, u1, v2) + sqr_dist(dim, v1, u2);
116 if (dist1 > dist2){
117 u = u2;
118 u2 = v2;
119 v2 = u;
120 dist1 = dist2;
121 flipped = TRUE;
122 }
123 len1 = dist(dim, u1, v1);
124 len2 = dist(dim, u2, v2);
125 //dist1 = MAX(0.1, dist1/(len1+len2+dist1));
126 dist1 = MAX(0.1, dist1/(len1+len2+0.0001*dist1));
127 if (flipped){
128 return -1/dist1;
129 } else {
130 return 1/dist1;
131 }
132}
133
134static real edge_compatibility_full(pedge e1, pedge e2){
135 /* two edges are u1->v1, u2->v2.
136 return 1 if two edges are exactly the same, 0 if they are very different.
137 This is based on Holten and van Wijk's paper
138 */
139 real *u1, *v1, *u2, *v2, *u, dist1, dist2, len1, len2, len;
140 real tmp, ca, cp, cs;
141 int dim = e1->dim, flipped = FALSE, i;
142
143 u1 = e1->x;
144 v1 = (e1->x)+(e1->npoints)*dim-dim;
145 u2 = e2->x;
146 v2 = (e2->x)+(e2->npoints)*dim-dim;
147 dist1 = sqr_dist(dim, u1, u2) + sqr_dist(dim, v1, v2);
148 dist2 = sqr_dist(dim, u1, v2) + sqr_dist(dim, v1, u2);
149 if (dist1 > dist2){
150 u = u2;
151 u2 = v2;
152 v2 = u;
153 dist1 = dist2;
154 flipped = TRUE;
155 }
156 len1 = MAX(dist(dim, u1, v1), SMALL);
157 len2 = MAX(dist(dim, u2, v2), SMALL);
158 len = 0.5*(len1+len2);
159
160 /* angle compatibility */
161 ca = 0;
162 for (i = 0; i < dim; i++)
163 ca += (v1[i]-u1[i])*(v2[i]-u2[i]);
164 ca = ABS(ca/(len1*len2));
165 assert(ca > -0.001);
166
167 /* scale compatibility */
168 //cs = 2/(len1/len2+len2/len1);
169 cs = 2/(MAX(len1,len2)/len + len/MIN(len1, len2));
170 assert(cs > -0.001 && cs < 1.001);
171
172 /* position compatibility */
173 cp = 0;
174 for (i = 0; i < dim; i++) {
175 tmp = .5*(v1[i]+u1[i])-.5*(v2[i]+u2[i]);
176 cp += tmp*tmp;
177 }
178 cp = sqrt(cp);
179 cp = len/(len + cp);
180 assert(cp > -0.001 && cp < 1.001);
181
182 /* visibility compatibility */
183
184 //dist1 = MAX(0.1, dist1/(len1+len2+dist1));
185 dist1 = cp*ca*cs;
186 if (flipped){
187 return -dist1;
188 } else {
189 return dist1;
190 }
191}
192
193static void fprint_rgb(FILE* fp, int r, int g, int b, int alpha){
194 fprintf(fp,"#");
195 if (r >= 16) {
196 fprintf(fp,"%2x",r);
197 } else {
198 fprintf(fp,"0%1x",r);
199 }
200 if (g >= 16) {
201 fprintf(fp,"%2x",g);
202 } else {
203 fprintf(fp,"0%1x",g);
204 }
205 if (b >= 16) {
206 fprintf(fp,"%2x",b);
207 } else {
208 fprintf(fp,"0%1x",b);
209 }
210 if (alpha >= 16) {
211 fprintf(fp,"%2x",alpha);
212 } else {
213 fprintf(fp,"0%1x",alpha);
214 }
215
216}
217
218void pedge_export_gv(FILE *fp, int ne, pedge *edges){
219 pedge edge;
220 real *x, t;
221 int i, j, k, kk, dim, sta, sto;
222 real maxwgt = 0, len, len_total, len_total0;
223 int r, g, b;
224 // real tt1[3]={0.25,0.5,0.75};
225 // real tt2[4]={0.2,0.4,0.6,0.8};
226 real tt1[3]={0.15,0.5,0.85};
227 real tt2[4]={0.15,0.4,0.6,0.85};
228 real *tt;
229
230 fprintf(fp,"strict graph{\n");
231 /* points */
232 for (i = 0; i < ne; i++){
233 edge = edges[i];
234 x = edge->x;
235 dim = edge->dim;
236 sta = 0; sto = edge->npoints - 1;
237
238 fprintf(fp, "%d [pos=\"", i);
239 for (k = 0; k < dim; k++) {
240 if (k != 0) fprintf(fp, ",");
241 fprintf(fp, "%f", x[sta*dim+k]);
242 }
243 fprintf(fp, "\"];\n");
244
245 fprintf(fp, "%d [pos=\"", i + ne);
246 for (k = 0; k < dim; k++) {
247 if (k != 0) fprintf(fp, ",");
248 fprintf(fp, "%f", x[sto*dim+k]);
249 }
250 fprintf(fp, "\"];\n");
251
252 }
253
254 /* figure out max number of bundled original edges in a pedge */
255 for (i = 0; i < ne; i++){
256 edge = edges[i];
257 if (edge->wgts){
258 for (j = 0; j < edge->npoints - 1; j++){
259 maxwgt = MAX(maxwgt, edge->wgts[j]);
260 }
261 }
262 }
263
264 /* spline and colors */
265 for (i = 0; i < ne; i++){
266 fprintf(fp,"%d -- %d [pos=\"", i, i + ne);
267 edge = edges[i];
268 x = edge->x;
269 dim = edge->dim;
270 /* splines */
271 for (j = 0; j < edge->npoints; j++){
272 if (j != 0) {
273 int mm = 3;
274 fprintf(fp," ");
275 /* there are ninterval+1 points, add 3*ninterval+2 points, get rid of internal ninternal-1 points,
276 make into 3*ninterval+4 points so that gviz spline rendering can work */
277 if (j == 1 || j == edge->npoints - 1) {
278 /* every interval gets 3 points inmserted except the first and last one */
279 tt = tt2;
280 mm = 4;
281 } else {
282 tt = tt1;
283 }
284 for (kk = 1; kk <= mm; kk++){
285 //t = kk/(real) (mm+1);
286 t = tt[kk-1];
287 for (k = 0; k < dim; k++) {
288 if (k != 0) fprintf(fp,",");
289 fprintf(fp, "%f", (x[(j-1)*dim+k]*(1-t)+x[j*dim+k]*(t)));
290 }
291 fprintf(fp," ");
292 }
293 }
294 if (j == 0 || j == edge->npoints - 1){
295 for (k = 0; k < dim; k++) {
296 if (k != 0) fprintf(fp,",");
297 fprintf(fp, "%f", x[j*dim+k]);
298 }
299 }
300 }
301 /* colors based on how much bundling */
302 if (edge->wgts){
303 fprintf(fp, "\", wgts=\"");
304 for (j = 0; j < edge->npoints - 1; j++){
305 if (j != 0) fprintf(fp,",");
306 fprintf(fp, "%f", edge->wgts[j]);
307 }
308
309 len_total = 0;
310 len_total0 = 0;
311 fprintf(fp, "\", color=\"");
312 for (j = 0; j < edge->npoints - 1; j++){
313 len = 0;
314 for (k = 0; k < dim; k++){
315 len += (edge->x[dim*j+k] - edge->x[dim*(j+1)+k])*(edge->x[dim*j+k] - edge->x[dim*(j+1)+k]);
316 }
317 len = sqrt(len/k);
318 len_total0 += len;
319 }
320 for (j = 0; j < edge->npoints - 1; j++){
321 len = 0;
322 for (k = 0; k < dim; k++){
323 len += (edge->x[dim*j+k] - edge->x[dim*(j+1)+k])*(edge->x[dim*j+k] - edge->x[dim*(j+1)+k]);
324 }
325 len = sqrt(len/k);
326 len_total += len;
327 t = edge->wgts[j]/maxwgt;
328 /* interpotate between red (t = 1) to blue (t = 0) */
329 r = 255*t; g = 0; b = 255*(1-t); b = 255*(1-t);
330 if (j != 0) fprintf(fp,":");
331 fprint_rgb(fp, r, g, b, 85);
332 if (j < edge->npoints - 2) fprintf(fp,";%f",len/len_total0);
333 }
334
335 }
336 fprintf(fp, "\"];\n");
337
338 }
339 fprintf(fp,"}\n");
340}
341
342void pedge_export_mma(FILE *fp, int ne, pedge *edges){
343 pedge edge;
344 real *x;
345 int i, j, k, dim;
346
347 fprintf(fp,"Graphics[{");
348 /* points */
349 fprintf(fp,"{Red, ");
350 for (i = 0; i < ne; i++){
351 if (i != 0) fprintf(fp,",");
352 fprintf(fp,"Point[");
353 edge = edges[i];
354 x = edge->x;
355 dim = edge->dim;
356 fprintf(fp, "{");
357 for (j = 0; j < edge->npoints; j+= edge->npoints - 1){
358 if (j != 0) fprintf(fp,",");
359 fprintf(fp, "{");
360 for (k = 0; k < dim; k++) {
361 if (k != 0) fprintf(fp,",");
362 fprintf(fp, "%f", x[j*dim+k]);
363 }
364 fprintf(fp, "}");
365 }
366 fprintf(fp, "}");
367 fprintf(fp, "]");
368 }
369 fprintf(fp,"},\n{GrayLevel[0.5,0.2], ");
370
371 /* spline */
372 for (i = 0; i < ne; i++){
373 if (i != 0) fprintf(fp,",");
374 fprintf(fp,"Spline[");
375 edge = edges[i];
376 x = edge->x;
377 dim = edge->dim;
378 fprintf(fp, "{");
379 for (j = 0; j < edge->npoints; j++){
380 if (j != 0) fprintf(fp,",");
381 fprintf(fp, "{");
382 for (k = 0; k < dim; k++) {
383 if (k != 0) fprintf(fp,",");
384 fprintf(fp, "%f", x[j*dim+k]);
385 }
386 fprintf(fp, "}");
387 }
388 fprintf(fp, "}");
389 fprintf(fp, ", Bezier]");
390 }
391 fprintf(fp,"}");
392
393 fprintf(fp,"}]\n");
394}
395
396#ifdef DEBUG
397static void pedge_print(char *comments, pedge e){
398 int i, j, dim;
399 dim = e->dim;
400 fprintf(stderr,"%s", comments);
401 for (i = 0; i < e->npoints; i++){
402 if (i > 0) fprintf(stderr,",");
403 fprintf(stderr,"{");
404 for (j = 0; j < dim; j++){
405 if (j > 0) fprintf(stderr,",");
406 fprintf(stderr,"%f",e->x[dim*i+j]);
407 }
408 fprintf(stderr,"}");
409 }
410 fprintf(stderr,"\n");
411}
412#endif
413
414pedge pedge_realloc(pedge e, int n){
415 if (n <= e->npoints) return e;
416 e->x = REALLOC(e->x, e->dim*n*sizeof(real));
417 if (e->wgts) e->wgts = REALLOC(e->wgts, (n-1)*sizeof(real));
418 e->len = n;
419 return e;
420}
421pedge pedge_wgts_realloc(pedge e, int n){
422 /* diff from pedge_alloc: allocate wgts if do not exist and initialize to wgt */
423 int i;
424 if (n <= e->npoints) return e;
425 e->x = REALLOC(e->x, e->dim*n*sizeof(real));
426 if (!(e->wgts)){
427 e->wgts = REALLOC(e->wgts, (n-1)*sizeof(real));
428 for (i = 0; i < e->npoints; i++) e->wgts[i] = e->wgt;
429 } else {
430 e->wgts = REALLOC(e->wgts, (n-1)*sizeof(real));
431 }
432 e->len = n;
433 return e;
434}
435
436
437pedge pedge_double(pedge e){
438 /* double the number of points (more precisely, add a point between two points in the polyline */
439 int npoints = e->npoints, len = e->len, i, dim = e->dim;
440 real *x;
441 int j, ii, ii2, np;
442
443 assert(npoints >= 2);
444 // pedge_print("before doubling edge = ", e);
445 if (npoints*2-1 > len){
446 len = 3*npoints;
447 e->x = REALLOC(e->x, dim*len*sizeof(real));
448 }
449 x = e->x;
450
451 for (i = npoints - 1; i >= 0; i--){
452 ii = 2*i;
453 for (j = 0; j < dim; j++){
454 x[dim*ii + j] = x[dim*i + j];
455 }
456 }
457
458 for (i = 0; i < npoints - 1; i++){
459 ii = 2*i;/* left and right interpolant of a new point */
460 ii2 = 2*(i+1);
461 for (j = 0; j < dim; j++){
462 x[dim*(2*i + 1) + j] = 0.5*(x[dim*ii + j] + x[dim*ii2 + j]);
463 }
464 }
465 e->len = len;
466 np = e->npoints = 2*(e->npoints) - 1;
467 e->edge_length = dist(dim, &(x[0*dim]), &(x[(np-1)*dim]));
468
469 //pedge_print("after doubling edge = ", e);
470
471 return e;
472}
473
474static void edge_tension_force(real *force, pedge e){
475 real *x = e->x;
476 int dim = e->dim;
477 int np = e->npoints;
478 int i, left, right, j;
479 //real dist_left, dist_right;
480 real s;
481
482
483 /* tention force = ((np-1)*||2x-xleft-xright||)/||e||, so the force is norminal and unitless
484 */
485 //s = (np-1)*(np-1)/MAX(SMALL, e->edge_length);
486 s = (np-1)/MAX(SMALL, e->edge_length);
487 for (i = 1; i <= np - 2; i++){
488 left = i - 1;
489 right = i + 1;
490 // dist_left = dist(dim, &(x[i*dim]), &(x[left*dim]));
491 // dist_right = dist(dim, &(x[i*dim]), &(x[right*dim]));
492 for (j = 0; j < dim; j++) force[i*dim + j] += s*(x[left*dim + j] - x[i*dim + j]);
493 for (j = 0; j < dim; j++) force[i*dim + j] += s*(x[right*dim + j] - x[i*dim + j]);
494 }
495}
496
497
498#if 0
499static void edge_tension_force2(real *force, pedge e){
500 /* in this version each node is pulled towards its original position on the line */
501 real *x = e->x;
502 int dim = e->dim;
503 int np = e->npoints;
504 int i, j;
505 //int left, right;
506 //real dist_left, dist_right;
507 real s;
508
509
510 /* tention force = ((np-1)*||2x-xleft-xright||)/||e||, so the force is norminal and unitless
511 */
512 s = .1/MAX(SMALL, e->edge_length);
513 for (i = 1; i <= np - 2; i++){
514 //left = i - 1;
515 // right = i + 1;
516 // dist_left = dist(dim, &(x[i*dim]), &(x[left*dim]));
517 // dist_right = dist(dim, &(x[i*dim]), &(x[right*dim]));
518 for (j = 0; j < dim; j++) force[i*dim + j] += s*((i*x[0*dim + j]+(np-1-i)*x[(np-1)*dim + j])/(np-1) - x[i*dim + j]);
519 }
520}
521#endif
522
523static void edge_attraction_force(real similarity, pedge e1, pedge e2, real *force){
524 /* attrractive force from x2 applied to x1 */
525 real *x1 = e1->x, *x2 = e2->x;
526 int dim = e1->dim;
527 int np = e1->npoints;
528 int i, j;
529 real dist, s, ss;
530 real edge_length = e1->edge_length;
531
532
533 assert(e1->npoints == e2->npoints);
534
535 /* attractive force = 1/d where d = D/||e1|| is the relative distance, D is the distance between e1 and e2.
536 so the force is norminal and unitless
537 */
538 if (similarity > 0){
539 s = edge_length;
540 s = similarity*edge_length;
541 for (i = 1; i <= np - 2; i++){
542 dist = sqr_dist(dim, &(x1[i*dim]), &(x2[i*dim]));
543 if (dist < SMALL) dist = SMALL;
544 ss = s/(dist+0.1*edge_length*sqrt(dist));
545 for (j = 0; j < dim; j++) force[i*dim + j] += ss*(x2[i*dim + j] - x1[i*dim + j]);
546 }
547 } else {/* clip e2 */
548 s = -edge_length;
549 s = -similarity*edge_length;
550 for (i = 1; i <= np - 2; i++){
551 dist = sqr_dist(dim, &(x1[i*dim]), &(x2[(np - 1 - i)*dim]));
552 if (dist < SMALL) dist = SMALL;
553 ss = s/(dist+0.1*edge_length*sqrt(dist));
554 for (j = 0; j < dim; j++) force[i*dim + j] += ss*(x2[(np - 1 - i)*dim + j] - x1[i*dim + j]);
555 }
556 }
557
558}
559
560static pedge* force_directed_edge_bundling(SparseMatrix A, pedge* edges, int maxit, real step0, real K, int open_gl){
561 int i, j, ne = A->n, k;
562 int *ia = A->ia, *ja = A->ja, iter = 0;
563 real *a = (real*) A->a;
564 pedge e1, e2;
565 real *force_t, *force_a;
566 int np = edges[0]->npoints, dim = edges[0]->dim;
567 real *x;
568 real step = step0;
569 real fnorm_a, fnorm_t, edge_length, start;
570
571 if (Verbose > 1)
572 fprintf(stderr, "total interaction pairs = %d out of %d, avg neighbors per edge = %f\n",A->nz, A->m*A->m, A->nz/(real) A->m);
573
574 force_t = MALLOC(sizeof(real)*dim*(np));
575 force_a = MALLOC(sizeof(real)*dim*(np));
576 while (step > 0.001 && iter < maxit){
577 start = clock();
578 iter++;
579 for (i = 0; i < ne; i++){
580 for (j = 0; j < dim*np; j++) {
581 force_t[j] = 0.;
582 force_a[j] = 0.;
583 }
584 e1 = edges[i];
585 x = e1->x;
586 //edge_tension_force2(force_t, e1);
587 edge_tension_force(force_t, e1);
588 for (j = ia[i]; j < ia[i+1]; j++){
589 e2 = edges[ja[j]];
590 edge_attraction_force(a[j], e1, e2, force_a);
591 }
592 fnorm_t = MAX(SMALL, norm(dim*(np - 2), &(force_t[1*dim])));
593 fnorm_a = MAX(SMALL, norm(dim*(np - 2), &(force_a[1*dim])));
594 edge_length = e1->edge_length;
595
596 // fprintf(stderr,"tension force norm = %f att force norm = %f step = %f edge_length = %f\n",fnorm_t, fnorm_a, step, edge_length);
597 for (j = 1; j <= np - 2; j++){
598 for (k = 0; k < dim; k++) {
599 x[j*dim + k] += step*edge_length*(force_t[j*dim+k] + K*force_a[j*dim+k])/(sqrt(fnorm_t*fnorm_t + K*K*fnorm_a*fnorm_a));
600 }
601 }
602
603 /*
604 fprintf(stderr,"edge %d",i);
605 for (j = 0; j < np; j++){
606 if (j != 0) fprintf(stderr,",");
607 fprintf(stderr,"{");
608 for (k = 0; k < dim; k++) {
609 if (k != 0) fprintf(stderr,",");
610 fprintf(stderr,"%f", x[j*dim+k]);
611 }
612 fprintf(stderr,"}");
613 }
614 fprintf(stderr,"\n");
615 */
616
617 }
618 step = step*0.9;
619 if (Verbose > 1)
620 fprintf(stderr, "iter ==== %d cpu = %f npoints = %d\n",iter, ((real) (clock() - start))/CLOCKS_PER_SEC, np - 2);
621
622#ifdef OPENGL
623 if (open_gl){
624 edges_global = edges;
625 drawScene();
626 }
627#endif
628
629 }
630
631 FREE(force_t);
632 FREE(force_a);
633 return edges;
634}
635
636static real absfun(real x){
637 return ABS(x);
638}
639
640
641
642
643static pedge* modularity_ink_bundling(int dim, int ne, SparseMatrix B, pedge* edges, real angle_param, real angle){
644 int *assignment = NULL, flag, nclusters;
645 real modularity;
646 int *clusterp, *clusters;
647 SparseMatrix D, C;
648 point_t meet1, meet2;
649 real ink0, ink1;
650 pedge e;
651 int i, j, jj;
652 int use_value_for_clustering = TRUE;
653
654 //int use_value_for_clustering = FALSE;
655
656 SparseMatrix BB;
657
658 /* B may contain negative entries */
659 BB = SparseMatrix_copy(B);
660 BB = SparseMatrix_apply_fun(BB, absfun);
661 modularity_clustering(BB, TRUE, 0, use_value_for_clustering, &nclusters, &assignment, &modularity, &flag);
662 SparseMatrix_delete(BB);
663
664#ifdef OPENGL
665 clusters_global = assignment;
666#endif
667
668 assert(!flag);
669 if (Verbose > 1) fprintf(stderr, "there are %d clusters, modularity = %f\n",nclusters, modularity);
670
671 C = SparseMatrix_new(1, 1, 1, MATRIX_TYPE_PATTERN, FORMAT_COORD);
672
673 for (i = 0; i < ne; i++){
674 jj = assignment[i];
675 SparseMatrix_coordinate_form_add_entries(C, 1, &jj, &i, NULL);
676 }
677
678 D = SparseMatrix_from_coordinate_format(C);
679 SparseMatrix_delete(C);
680 clusterp = D->ia;
681 clusters = D->ja;
682 for (i = 0; i < nclusters; i++){
683 ink1 = ink(edges, clusterp[i+1] - clusterp[i], &(clusters[clusterp[i]]), &ink0, &meet1, &meet2, angle_param, angle);
684 if (Verbose > 1)
685 fprintf(stderr,"nedges = %d ink0 = %f, ink1 = %f\n",clusterp[i+1] - clusterp[i], ink0, ink1);
686 if (ink1 < ink0){
687 for (j = clusterp[i]; j < clusterp[i+1]; j++){
688 /* make this edge 5 points, insert two meeting points at 1 and 2, make 3 the last point */
689 edges[clusters[j]] = pedge_double(edges[clusters[j]]);
690 e = edges[clusters[j]] = pedge_double(edges[clusters[j]]);
691 e->x[1*dim] = meet1.x;
692 e->x[1*dim+1] = meet1.y;
693 e->x[2*dim] = meet2.x;
694 e->x[2*dim+1] = meet2.y;
695 e->x[3*dim] = e->x[4*dim];
696 e->x[3*dim+1] = e->x[4*dim+1];
697 e->npoints = 4;
698 }
699#ifdef OPENGL
700 edges_global = edges;
701 drawScene();
702#endif
703 }
704 }
705 SparseMatrix_delete(D);
706 return edges;
707}
708
709static SparseMatrix check_compatibility(SparseMatrix A, int ne, pedge *edges, int compatibility_method, real tol){
710 /* go through the links and make sure edges are compatible */
711 SparseMatrix B, C;
712 int *ia, *ja, i, j, jj;
713 real start;
714 real dist;
715
716 B = SparseMatrix_new(1, 1, 1, MATRIX_TYPE_REAL, FORMAT_COORD);
717 ia = A->ia; ja = A->ja;
718 //SparseMatrix_print("A",A);
719 start = clock();
720 for (i = 0; i < ne; i++){
721 for (j = ia[i]; j < ia[i+1]; j++){
722 jj = ja[j];
723 if (i == jj) continue;
724 if (compatibility_method == COMPATIBILITY_DIST){
725 dist = edge_compatibility_full(edges[i], edges[jj]);
726 } else if (compatibility_method == COMPATIBILITY_FULL){
727 dist = edge_compatibility(edges[i], edges[jj]);
728 }
729 /*
730 fprintf(stderr,"edge %d = {",i);
731 pedge_print("", edges[i]);
732 fprintf(stderr,"edge %d = {",jj);
733 pedge_print("", edges[jj]);
734 fprintf(stderr,"compatibility=%f\n",dist);
735 */
736
737 if (ABS(dist) > tol){
738 B = SparseMatrix_coordinate_form_add_entries(B, 1, &i, &jj, &dist);
739 B = SparseMatrix_coordinate_form_add_entries(B, 1, &jj, &i, &dist);
740 }
741 }
742 }
743 C = SparseMatrix_from_coordinate_format(B);
744 //SparseMatrix_print("C",C);
745 SparseMatrix_delete(B);
746 B = C;
747 if (Verbose > 1)
748 fprintf(stderr, "edge compatibilitu time = %f\n",((real) (clock() - start))/CLOCKS_PER_SEC);
749 return B;
750}
751
752pedge* edge_bundling(SparseMatrix A0, int dim, real *x, int maxit_outer, real K, int method, int nneighbor, int compatibility_method,
753 int max_recursion, real angle_param, real angle, int open_gl){
754 /* bundle edges.
755 A: edge graph
756 x: edge i is at {p,q},
757 . where p = x[2*dim*i : 2*dim*i+dim-1]
758 . and q = x[2*dim*i+dim : 2*dim*i+2*dim-1]
759 maxit_outer: max outer iteration for force directed bundling. Every outer iter subdivide each edge segment into 2.
760 K: norminal edge length in force directed bundling
761 method: which method to use.
762 nneighbor: number of neighbors to be used in forming nearest neighbor graph. Used only in agglomerative method
763 compatibility_method: which method to use to calculate compatibility. Used only in force directed.
764 max_recursion: used only in agglomerative method. Specify how many level of recursion to do to bundle bundled edges again
765 open_gl: whether to plot in X.
766
767 */
768 int ne = A0->m;
769 pedge *edges;
770 SparseMatrix A = A0, B = NULL;
771 int i;
772 real tol = 0.001;
773 int k;
774 real step0 = 0.1, start = 0.0;
775 int maxit = 10;
776 int flag;
777
778 assert(A->n == ne);
779 edges = MALLOC(sizeof(pedge)*ne);
780
781 for (i = 0; i < ne; i++){
782 edges[i] = pedge_new(2, dim, &(x[dim*2*i]));
783 }
784
785 A = SparseMatrix_symmetrize(A0, TRUE);
786
787
788
789 if (Verbose) start = clock();
790 if (method == METHOD_INK){
791
792 /* go through the links and make sure edges are compatible */
793 B = check_compatibility(A, ne, edges, compatibility_method, tol);
794
795 edges = modularity_ink_bundling(dim, ne, B, edges, angle_param, angle);
796
797 } else if (method == METHOD_INK_AGGLOMERATE){
798#ifdef HAVE_ANN
799 /* plan: merge a node with its neighbors if doing so improve. Form coarsening graph, repeat until no more ink saving */
800 edges = agglomerative_ink_bundling(dim, A, edges, nneighbor, max_recursion, angle_param, angle, open_gl, &flag);
801 assert(!flag);
802#else
803 agerr (AGERR, "Graphviz built without approximate nearest neighbor library ANN; agglomerative inking not available\n");
804 edges = edges;
805#endif
806 } else if (method == METHOD_FD){/* FD method */
807
808 /* go through the links and make sure edges are compatible */
809 B = check_compatibility(A, ne, edges, compatibility_method, tol);
810
811
812 for (k = 0; k < maxit_outer; k++){
813 for (i = 0; i < ne; i++){
814 edges[i] = pedge_double(edges[i]);
815 }
816 step0 = step0/2;
817 edges = force_directed_edge_bundling(B, edges, maxit, step0, K, open_gl);
818 }
819
820 } else if (method == METHOD_NONE){
821 edges = edges;
822 } else {
823 assert(0);
824 }
825 if (Verbose)
826 fprintf(stderr, "total edge bundling cpu = %f\n",((real) (clock() - start))/CLOCKS_PER_SEC);
827 if (B != A) SparseMatrix_delete(B);
828 if (A != A0) SparseMatrix_delete(A);
829
830 return edges;
831}
832