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