1/* $Id$ $Revision$ */
2/* vim:set shiftwidth=4 ts=8: */
3
4/*************************************************************************
5 * Copyright (c) 2011 AT&T Intellectual Property
6 * All rights reserved. This program and the accompanying materials
7 * are made available under the terms of the Eclipse Public License v1.0
8 * which accompanies this distribution, and is available at
9 * http://www.eclipse.org/legal/epl-v10.html
10 *
11 * Contributors: See CVS logs. Details at http://www.graphviz.org/
12 *************************************************************************/
13
14
15#include "config.h"
16
17#include "neato.h"
18#include "stress.h"
19#include <time.h>
20#ifndef _WIN32
21#include <unistd.h>
22#endif
23
24static double Epsilon2;
25
26
27double fpow32(double x)
28{
29 x = sqrt(x);
30 return x * x * x;
31}
32
33double distvec(double *p0, double *p1, double *vec)
34{
35 int k;
36 double dist = 0.0;
37
38 for (k = 0; k < Ndim; k++) {
39 vec[k] = p0[k] - p1[k];
40 dist += (vec[k] * vec[k]);
41 }
42 dist = sqrt(dist);
43 return dist;
44}
45
46double **new_array(int m, int n, double ival)
47{
48 double **rv;
49 double *mem;
50 int i, j;
51
52 rv = N_NEW(m, double *);
53 mem = N_NEW(m * n, double);
54 for (i = 0; i < m; i++) {
55 rv[i] = mem;
56 mem = mem + n;
57 for (j = 0; j < n; j++)
58 rv[i][j] = ival;
59 }
60 return rv;
61}
62
63void free_array(double **rv)
64{
65 if (rv) {
66 free(rv[0]);
67 free(rv);
68 }
69}
70
71
72static double ***new_3array(int m, int n, int p, double ival)
73{
74 double ***rv;
75 int i, j, k;
76
77 rv = N_NEW(m + 1, double **);
78 for (i = 0; i < m; i++) {
79 rv[i] = N_NEW(n + 1, double *);
80 for (j = 0; j < n; j++) {
81 rv[i][j] = N_NEW(p, double);
82 for (k = 0; k < p; k++)
83 rv[i][j][k] = ival;
84 }
85 rv[i][j] = NULL; /* NULL terminate so we can clean up */
86 }
87 rv[i] = NULL;
88 return rv;
89}
90
91static void free_3array(double ***rv)
92{
93 int i, j;
94
95 if (rv) {
96 for (i = 0; rv[i]; i++) {
97 for (j = 0; rv[i][j]; j++)
98 free(rv[i][j]);
99 free(rv[i]);
100 }
101 free(rv);
102 }
103}
104
105
106/* lenattr:
107 * Return 1 if attribute not defined
108 * Return 2 if attribute string bad
109 */
110static int lenattr(edge_t* e, Agsym_t* index, double* val)
111{
112 char* s;
113
114 if (index == NULL)
115 return 1;
116
117 s = agxget(e, index);
118 if (*s == '\0') return 1;
119
120 if ((sscanf(s, "%lf", val) < 1) || (*val < 0) || ((*val == 0) && !Nop)) {
121 agerr(AGWARN, "bad edge len \"%s\"", s);
122 return 2;
123 }
124 else
125 return 0;
126}
127
128
129/* degreeKind;
130 * Returns degree of n ignoring loops and multiedges.
131 * Returns 0, 1 or many (2)
132 * For case of 1, returns other endpoint of edge.
133 */
134static int degreeKind(graph_t * g, node_t * n, node_t ** op)
135{
136 edge_t *ep;
137 int deg = 0;
138 node_t *other = NULL;
139
140 for (ep = agfstedge(g, n); ep; ep = agnxtedge(g, ep, n)) {
141 if (aghead(ep) == agtail(ep))
142 continue; /* ignore loops */
143 if (deg == 1) {
144 if (((agtail(ep) == n) && (aghead(ep) == other)) || /* ignore multiedge */
145 ((agtail(ep) == other) && (aghead(ep) == n)))
146 continue;
147 return 2;
148 } else { /* deg == 0 */
149 if (agtail(ep) == n)
150 other = aghead(ep);
151 else
152 other = agtail(ep);
153 *op = other;
154 deg++;
155 }
156 }
157 return deg;
158}
159
160
161/* prune:
162 * np is at end of a chain. If its degree is 0, remove it.
163 * If its degree is 1, remove it and recurse.
164 * If its degree > 1, stop.
165 * The node next is the next node to be visited during iteration.
166 * If it is equal to a node being deleted, set it to next one.
167 * Delete from root graph, since G may be a connected component subgraph.
168 * Return next.
169 */
170static node_t *prune(graph_t * G, node_t * np, node_t * next)
171{
172 node_t *other;
173 int deg;
174
175 while (np) {
176 deg = degreeKind(G, np, &other);
177 if (deg == 0) {
178 if (next == np)
179 next = agnxtnode(G, np);
180 agdelete(G->root, np);
181 np = 0;
182 } else if (deg == 1) {
183 if (next == np)
184 next = agnxtnode(G, np);
185 agdelete(G->root, np);
186 np = other;
187 } else
188 np = 0;
189
190 }
191 return next;
192}
193
194static double setEdgeLen(graph_t * G, node_t * np, Agsym_t* lenx, double dfltlen)
195{
196 edge_t *ep;
197 double total_len = 0.0;
198 double len;
199 int err;
200
201 for (ep = agfstout(G, np); ep; ep = agnxtout(G, ep)) {
202 if ((err = lenattr(ep, lenx, &len))) {
203 if (err == 2) agerr(AGPREV, " in %s - setting to %.02f\n", agnameof(G), dfltlen);
204 len = dfltlen;
205 }
206 ED_dist(ep) = len;
207 total_len += len;
208 }
209 return total_len;
210}
211
212/* scan_graph_mode:
213 * Prepare the graph and data structures depending on the layout mode.
214 * If Reduce is true, eliminate singletons and trees. Since G may be a
215 * subgraph, we remove the nodes from the root graph.
216 * Return the number of nodes in the reduced graph.
217 */
218int scan_graph_mode(graph_t * G, int mode)
219{
220 int i, nV, nE, deg;
221 char *str;
222 node_t *np, *xp, *other;
223 double total_len = 0.0;
224 double dfltlen = 1.0;
225 Agsym_t* lenx;
226
227 if (Verbose)
228 fprintf(stderr, "Scanning graph %s, %d nodes\n", agnameof(G),
229 agnnodes(G));
230
231
232 /* Eliminate singletons and trees */
233 if (Reduce) {
234 for (np = agfstnode(G); np; np = xp) {
235 xp = agnxtnode(G, np);
236 deg = degreeKind(G, np, &other);
237 if (deg == 0) { /* singleton node */
238 agdelete(G->root, np);
239 } else if (deg == 1) {
240 agdelete(G->root, np);
241 xp = prune(G, other, xp);
242 }
243 }
244 }
245
246 nV = agnnodes(G);
247 nE = agnedges(G);
248
249 lenx = agattr(G, AGEDGE, "len", 0);
250 if (mode == MODE_KK) {
251 Epsilon = .0001 * nV;
252 getdouble(G, "epsilon", &Epsilon);
253 if ((str = agget(G->root, "Damping")))
254 Damping = atof(str);
255 else
256 Damping = .99;
257 GD_neato_nlist(G) = N_NEW(nV + 1, node_t *);
258 for (i = 0, np = agfstnode(G); np; np = agnxtnode(G, np)) {
259 GD_neato_nlist(G)[i] = np;
260 ND_id(np) = i++;
261 ND_heapindex(np) = -1;
262 total_len += setEdgeLen(G, np, lenx, dfltlen);
263 }
264 } else {
265 Epsilon = DFLT_TOLERANCE;
266 getdouble(G, "epsilon", &Epsilon);
267 for (i = 0, np = agfstnode(G); np; np = agnxtnode(G, np)) {
268 ND_id(np) = i++;
269 total_len += setEdgeLen(G, np, lenx, dfltlen);
270 }
271 }
272
273 str = agget(G, "defaultdist");
274 if (str && str[0])
275 Initial_dist = MAX(Epsilon, atof(str));
276 else
277 Initial_dist = total_len / (nE > 0 ? nE : 1) * sqrt(nV) + 1;
278
279 if (!Nop && (mode == MODE_KK)) {
280 GD_dist(G) = new_array(nV, nV, Initial_dist);
281 GD_spring(G) = new_array(nV, nV, 1.0);
282 GD_sum_t(G) = new_array(nV, Ndim, 1.0);
283 GD_t(G) = new_3array(nV, nV, Ndim, 0.0);
284 }
285
286 return nV;
287}
288
289int scan_graph(graph_t * g)
290{
291 return scan_graph_mode(g, MODE_KK);
292}
293
294void free_scan_graph(graph_t * g)
295{
296 free(GD_neato_nlist(g));
297 if (!Nop) {
298 free_array(GD_dist(g));
299 free_array(GD_spring(g));
300 free_array(GD_sum_t(g));
301 free_3array(GD_t(g));
302 GD_t(g) = NULL;
303 }
304}
305
306void jitter_d(node_t * np, int nG, int n)
307{
308 int k;
309 for (k = n; k < Ndim; k++)
310 ND_pos(np)[k] = nG * drand48();
311}
312
313void jitter3d(node_t * np, int nG)
314{
315 jitter_d(np, nG, 2);
316}
317
318void randompos(node_t * np, int nG)
319{
320 ND_pos(np)[0] = nG * drand48();
321 ND_pos(np)[1] = nG * drand48();
322 if (Ndim > 2)
323 jitter3d(np, nG);
324}
325
326void initial_positions(graph_t * G, int nG)
327{
328 int init, i;
329 node_t *np;
330 static int once = 0;
331
332 if (Verbose)
333 fprintf(stderr, "Setting initial positions\n");
334
335 init = checkStart(G, nG, INIT_RANDOM);
336 if (init == INIT_REGULAR)
337 return;
338 if ((init == INIT_SELF) && (once == 0)) {
339 agerr(AGWARN, "start=%s not supported with mode=self - ignored\n");
340 once = 1;
341 }
342
343 for (i = 0; (np = GD_neato_nlist(G)[i]); i++) {
344 if (hasPos(np))
345 continue;
346 randompos(np, 1);
347 }
348}
349
350void diffeq_model(graph_t * G, int nG)
351{
352 int i, j, k;
353 double dist, **D, **K, del[MAXDIM], f;
354 node_t *vi, *vj;
355 edge_t *e;
356
357 if (Verbose) {
358 fprintf(stderr, "Setting up spring model: ");
359 start_timer();
360 }
361 /* init springs */
362 K = GD_spring(G);
363 D = GD_dist(G);
364 for (i = 0; i < nG; i++) {
365 for (j = 0; j < i; j++) {
366 f = Spring_coeff / (D[i][j] * D[i][j]);
367 if ((e = agfindedge(G, GD_neato_nlist(G)[i], GD_neato_nlist(G)[j])))
368 f = f * ED_factor(e);
369 K[i][j] = K[j][i] = f;
370 }
371 }
372
373 /* init differential equation solver */
374 for (i = 0; i < nG; i++)
375 for (k = 0; k < Ndim; k++)
376 GD_sum_t(G)[i][k] = 0.0;
377
378 for (i = 0; (vi = GD_neato_nlist(G)[i]); i++) {
379 for (j = 0; j < nG; j++) {
380 if (i == j)
381 continue;
382 vj = GD_neato_nlist(G)[j];
383 dist = distvec(ND_pos(vi), ND_pos(vj), del);
384 for (k = 0; k < Ndim; k++) {
385 GD_t(G)[i][j][k] =
386 GD_spring(G)[i][j] * (del[k] -
387 GD_dist(G)[i][j] * del[k] /
388 dist);
389 GD_sum_t(G)[i][k] += GD_t(G)[i][j][k];
390 }
391 }
392 }
393 if (Verbose) {
394 fprintf(stderr, "%.2f sec\n", elapsed_sec());
395 }
396}
397
398/* total_e:
399 * Return 2*energy of system.
400 */
401static double total_e(graph_t * G, int nG)
402{
403 int i, j, d;
404 double e = 0.0; /* 2*energy */
405 double t0; /* distance squared */
406 double t1;
407 node_t *ip, *jp;
408
409 for (i = 0; i < nG - 1; i++) {
410 ip = GD_neato_nlist(G)[i];
411 for (j = i + 1; j < nG; j++) {
412 jp = GD_neato_nlist(G)[j];
413 for (t0 = 0.0, d = 0; d < Ndim; d++) {
414 t1 = (ND_pos(ip)[d] - ND_pos(jp)[d]);
415 t0 += t1 * t1;
416 }
417 e = e + GD_spring(G)[i][j] *
418 (t0 + GD_dist(G)[i][j] * GD_dist(G)[i][j]
419 - 2.0 * GD_dist(G)[i][j] * sqrt(t0));
420 }
421 }
422 return e;
423}
424
425void solve_model(graph_t * G, int nG)
426{
427 node_t *np;
428
429 Epsilon2 = Epsilon * Epsilon;
430
431 while ((np = choose_node(G, nG))) {
432 move_node(G, nG, np);
433 }
434 if (Verbose) {
435 fprintf(stderr, "\nfinal e = %f", total_e(G, nG));
436 fprintf(stderr, " %d%s iterations %.2f sec\n",
437 GD_move(G), (GD_move(G) == MaxIter ? "!" : ""),
438 elapsed_sec());
439 }
440 if (GD_move(G) == MaxIter)
441 agerr(AGWARN, "Max. iterations (%d) reached on graph %s\n",
442 MaxIter, agnameof(G));
443}
444
445void update_arrays(graph_t * G, int nG, int i)
446{
447 int j, k;
448 double del[MAXDIM], dist, old;
449 node_t *vi, *vj;
450
451 vi = GD_neato_nlist(G)[i];
452 for (k = 0; k < Ndim; k++)
453 GD_sum_t(G)[i][k] = 0.0;
454 for (j = 0; j < nG; j++) {
455 if (i == j)
456 continue;
457 vj = GD_neato_nlist(G)[j];
458 dist = distvec(ND_pos(vi), ND_pos(vj), del);
459 for (k = 0; k < Ndim; k++) {
460 old = GD_t(G)[i][j][k];
461 GD_t(G)[i][j][k] =
462 GD_spring(G)[i][j] * (del[k] -
463 GD_dist(G)[i][j] * del[k] / dist);
464 GD_sum_t(G)[i][k] += GD_t(G)[i][j][k];
465 old = GD_t(G)[j][i][k];
466 GD_t(G)[j][i][k] = -GD_t(G)[i][j][k];
467 GD_sum_t(G)[j][k] += (GD_t(G)[j][i][k] - old);
468 }
469 }
470}
471
472#define Msub(i,j) M[(i)*Ndim+(j)]
473void D2E(graph_t * G, int nG, int n, double *M)
474{
475 int i, l, k;
476 node_t *vi, *vn;
477 double scale, sq, t[MAXDIM];
478 double **K = GD_spring(G);
479 double **D = GD_dist(G);
480
481 vn = GD_neato_nlist(G)[n];
482 for (l = 0; l < Ndim; l++)
483 for (k = 0; k < Ndim; k++)
484 Msub(l, k) = 0.0;
485 for (i = 0; i < nG; i++) {
486 if (n == i)
487 continue;
488 vi = GD_neato_nlist(G)[i];
489 sq = 0.0;
490 for (k = 0; k < Ndim; k++) {
491 t[k] = ND_pos(vn)[k] - ND_pos(vi)[k];
492 sq += (t[k] * t[k]);
493 }
494 scale = 1 / fpow32(sq);
495 for (k = 0; k < Ndim; k++) {
496 for (l = 0; l < k; l++)
497 Msub(l, k) += K[n][i] * D[n][i] * t[k] * t[l] * scale;
498 Msub(k, k) +=
499 K[n][i] * (1.0 - D[n][i] * (sq - (t[k] * t[k])) * scale);
500 }
501 }
502 for (k = 1; k < Ndim; k++)
503 for (l = 0; l < k; l++)
504 Msub(k, l) = Msub(l, k);
505}
506
507void final_energy(graph_t * G, int nG)
508{
509 fprintf(stderr, "iterations = %d final e = %f\n", GD_move(G),
510 total_e(G, nG));
511}
512
513node_t *choose_node(graph_t * G, int nG)
514{
515 int i, k;
516 double m, max;
517 node_t *choice, *np;
518 static int cnt = 0;
519#if 0
520 double e;
521 static double save_e = MAXDOUBLE;
522#endif
523
524 cnt++;
525 if (GD_move(G) >= MaxIter)
526 return NULL;
527#if 0
528 if ((cnt % 100) == 0) {
529 e = total_e(G, nG);
530 if (e - save_e > 0)
531 return NULL;
532 save_e = e;
533 }
534#endif
535 max = 0.0;
536 choice = NULL;
537 for (i = 0; i < nG; i++) {
538 np = GD_neato_nlist(G)[i];
539 if (ND_pinned(np) > P_SET)
540 continue;
541 for (m = 0.0, k = 0; k < Ndim; k++)
542 m += (GD_sum_t(G)[i][k] * GD_sum_t(G)[i][k]);
543 /* could set the color=energy of the node here */
544 if (m > max) {
545 choice = np;
546 max = m;
547 }
548 }
549 if (max < Epsilon2)
550 choice = NULL;
551 else {
552 if (Verbose && (cnt % 100 == 0)) {
553 fprintf(stderr, "%.3f ", sqrt(max));
554 if (cnt % 1000 == 0)
555 fprintf(stderr, "\n");
556 }
557#if 0
558 e = total_e(G, nG);
559 if (fabs((e - save_e) / save_e) < 1e-5) {
560 choice = NULL;
561 }
562#endif
563 }
564 return choice;
565}
566
567void move_node(graph_t * G, int nG, node_t * n)
568{
569 int i, m;
570 static double *a, b[MAXDIM], c[MAXDIM];
571
572 m = ND_id(n);
573 a = ALLOC(Ndim * Ndim, a, double);
574 D2E(G, nG, m, a);
575 for (i = 0; i < Ndim; i++)
576 c[i] = -GD_sum_t(G)[m][i];
577 solve(a, b, c, Ndim);
578 for (i = 0; i < Ndim; i++) {
579 b[i] = (Damping + 2 * (1 - Damping) * drand48()) * b[i];
580 ND_pos(n)[i] += b[i];
581 }
582 GD_move(G)++;
583 update_arrays(G, nG, m);
584 if (test_toggle()) {
585 double sum = 0;
586 for (i = 0; i < Ndim; i++) {
587 sum += fabs(b[i]);
588 } /* Why not squared? */
589 sum = sqrt(sum);
590 fprintf(stderr, "%s %.3f\n", agnameof(n), sum);
591 }
592}
593
594static node_t **Heap;
595static int Heapsize;
596static node_t *Src;
597
598void heapup(node_t * v)
599{
600 int i, par;
601 node_t *u;
602
603 for (i = ND_heapindex(v); i > 0; i = par) {
604 par = (i - 1) / 2;
605 u = Heap[par];
606 if (ND_dist(u) <= ND_dist(v))
607 break;
608 Heap[par] = v;
609 ND_heapindex(v) = par;
610 Heap[i] = u;
611 ND_heapindex(u) = i;
612 }
613}
614
615void heapdown(node_t * v)
616{
617 int i, left, right, c;
618 node_t *u;
619
620 i = ND_heapindex(v);
621 while ((left = 2 * i + 1) < Heapsize) {
622 right = left + 1;
623 if ((right < Heapsize)
624 && (ND_dist(Heap[right]) < ND_dist(Heap[left])))
625 c = right;
626 else
627 c = left;
628 u = Heap[c];
629 if (ND_dist(v) <= ND_dist(u))
630 break;
631 Heap[c] = v;
632 ND_heapindex(v) = c;
633 Heap[i] = u;
634 ND_heapindex(u) = i;
635 i = c;
636 }
637}
638
639void neato_enqueue(node_t * v)
640{
641 int i;
642
643 assert(ND_heapindex(v) < 0);
644 i = Heapsize++;
645 ND_heapindex(v) = i;
646 Heap[i] = v;
647 if (i > 0)
648 heapup(v);
649}
650
651node_t *neato_dequeue(void)
652{
653 int i;
654 node_t *rv, *v;
655
656 if (Heapsize == 0)
657 return NULL;
658 rv = Heap[0];
659 i = --Heapsize;
660 v = Heap[i];
661 Heap[0] = v;
662 ND_heapindex(v) = 0;
663 if (i > 1)
664 heapdown(v);
665 ND_heapindex(rv) = -1;
666 return rv;
667}
668
669void shortest_path(graph_t * G, int nG)
670{
671 node_t *v;
672
673 Heap = N_NEW(nG + 1, node_t *);
674 if (Verbose) {
675 fprintf(stderr, "Calculating shortest paths: ");
676 start_timer();
677 }
678 for (v = agfstnode(G); v; v = agnxtnode(G, v))
679 s1(G, v);
680 if (Verbose) {
681 fprintf(stderr, "%.2f sec\n", elapsed_sec());
682 }
683 free(Heap);
684}
685
686void s1(graph_t * G, node_t * node)
687{
688 node_t *v, *u;
689 edge_t *e;
690 int t;
691 double f;
692
693 for (t = 0; (v = GD_neato_nlist(G)[t]); t++)
694 ND_dist(v) = Initial_dist;
695 Src = node;
696 ND_dist(Src) = 0;
697 ND_hops(Src) = 0;
698 neato_enqueue(Src);
699
700 while ((v = neato_dequeue())) {
701 if (v != Src)
702 make_spring(G, Src, v, ND_dist(v));
703 for (e = agfstedge(G, v); e; e = agnxtedge(G, e, v)) {
704 if ((u = agtail(e)) == v)
705 u = aghead(e);
706 f = ND_dist(v) + ED_dist(e);
707 if (ND_dist(u) > f) {
708 ND_dist(u) = f;
709 if (ND_heapindex(u) >= 0)
710 heapup(u);
711 else {
712 ND_hops(u) = ND_hops(v) + 1;
713 neato_enqueue(u);
714 }
715 }
716 }
717 }
718}
719
720void make_spring(graph_t * G, node_t * u, node_t * v, double f)
721{
722 int i, j;
723
724 i = ND_id(u);
725 j = ND_id(v);
726 GD_dist(G)[i][j] = GD_dist(G)[j][i] = f;
727}
728