Graphviz  2.35.20130930.0449
stress.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 #include "neato.h"
16 #include "dijkstra.h"
17 #include "bfs.h"
18 #include "pca.h"
19 #include "matrix_ops.h"
20 #include "conjgrad.h"
21 #include "embed_graph.h"
22 #include "kkutils.h"
23 #include "stress.h"
24 #include <math.h>
25 #include <stdlib.h>
26 #include <time.h>
27 
28 
29 #ifndef HAVE_DRAND48
30 extern double drand48(void);
31 #endif
32 
33 #define Dij2 /* If defined, the terms in the stress energy are normalized
34  by d_{ij}^{-2} otherwise, they are normalized by d_{ij}^{-1}
35  */
36 
37 #ifdef NONCORE
38 /* Set 'max_nodes_in_mem' so that
39  * 4*(max_nodes_in_mem^2) is smaller than the available memory (in bytes)
40  * 4 = sizeof(float)
41  */
42 #define max_nodes_in_mem 18000
43 #endif
44 
45  /* relevant when using sparse distance matrix not within subspace */
46 #define smooth_pivots true
47 
48 /* dimensionality of subspace; relevant
49  * when optimizing within subspace)
50  */
51 #define stress_pca_dim 50
52 
53  /* a structure used for storing sparse distance matrix */
54 typedef struct {
55  int nedges;
56  int *edges;
58  boolean free_mem;
59 } dist_data;
60 
61 static double compute_stressf(float **coords, float *lap, int dim, int n, int exp)
62 {
63  /* compute the overall stress */
64 
65  int i, j, l, neighbor, count;
66  double sum, dist, Dij;
67  sum = 0;
68  for (count = 0, i = 0; i < n - 1; i++) {
69  count++; /* skip diagonal entry */
70  for (j = 1; j < n - i; j++, count++) {
71  dist = 0;
72  neighbor = i + j;
73  for (l = 0; l < dim; l++) {
74  dist +=
75  (coords[l][i] - coords[l][neighbor]) * (coords[l][i] -
76  coords[l]
77  [neighbor]);
78  }
79  dist = sqrt(dist);
80  if (exp == 2) {
81 #ifdef Dij2
82  Dij = 1.0 / sqrt(lap[count]);
83  sum += (Dij - dist) * (Dij - dist) * (lap[count]);
84 #else
85  Dij = 1.0 / lap[count];
86  sum += (Dij - dist) * (Dij - dist) * (lap[count]);
87 #endif
88  } else {
89  Dij = 1.0 / lap[count];
90  sum += (Dij - dist) * (Dij - dist) * (lap[count]);
91  }
92  }
93  }
94 
95  return sum;
96 }
97 
98 static double
99 compute_stress1(double **coords, dist_data * distances, int dim, int n, int exp)
100 {
101  /* compute the overall stress */
102 
103  int i, j, l, node;
104  double sum, dist, Dij;
105  sum = 0;
106  if (exp == 2) {
107  for (i = 0; i < n; i++) {
108  for (j = 0; j < distances[i].nedges; j++) {
109  node = distances[i].edges[j];
110  if (node <= i) {
111  continue;
112  }
113  dist = 0;
114  for (l = 0; l < dim; l++) {
115  dist +=
116  (coords[l][i] - coords[l][node]) * (coords[l][i] -
117  coords[l]
118  [node]);
119  }
120  dist = sqrt(dist);
121  Dij = distances[i].edist[j];
122 #ifdef Dij2
123  sum += (Dij - dist) * (Dij - dist) / (Dij * Dij);
124 #else
125  sum += (Dij - dist) * (Dij - dist) / Dij;
126 #endif
127  }
128  }
129  } else {
130  for (i = 0; i < n; i++) {
131  for (j = 0; j < distances[i].nedges; j++) {
132  node = distances[i].edges[j];
133  if (node <= i) {
134  continue;
135  }
136  dist = 0;
137  for (l = 0; l < dim; l++) {
138  dist +=
139  (coords[l][i] - coords[l][node]) * (coords[l][i] -
140  coords[l]
141  [node]);
142  }
143  dist = sqrt(dist);
144  Dij = distances[i].edist[j];
145  sum += (Dij - dist) * (Dij - dist) / Dij;
146  }
147  }
148  }
149 
150  return sum;
151 }
152 
153 /* initLayout:
154  * Initialize node coordinates. If the node already has
155  * a position, use it.
156  * Return true if some node is fixed.
157  */
158 int
159 initLayout(vtx_data * graph, int n, int dim, double **coords,
160  node_t ** nodes)
161 {
162  node_t *np;
163  double *xp;
164  double *yp;
165  double *pt;
166  int i, d;
167  int pinned = 0;
168 
169  xp = coords[0];
170  yp = coords[1];
171  for (i = 0; i < n; i++) {
172  np = nodes[i];
173  if (hasPos(np)) {
174  pt = ND_pos(np);
175  *xp++ = *pt++;
176  *yp++ = *pt++;
177  if (dim > 2) {
178  for (d = 2; d < dim; d++)
179  coords[d][i] = *pt++;
180  }
181  if (isFixed(np))
182  pinned = 1;
183  } else {
184  *xp++ = drand48();
185  *yp++ = drand48();
186  if (dim > 2) {
187  for (d = 2; d < dim; d++)
188  coords[d][i] = drand48();
189  }
190  }
191  }
192 
193  for (d = 0; d < dim; d++)
194  orthog1(n, coords[d]);
195 
196  return pinned;
197 }
198 
199 float *circuitModel(vtx_data * graph, int nG)
200 {
201  int i, j, e, rv, count;
202  float *Dij = N_NEW(nG * (nG + 1) / 2, float);
203  double **Gm;
204  double **Gm_inv;
205 
206  Gm = new_array(nG, nG, 0.0);
207  Gm_inv = new_array(nG, nG, 0.0);
208 
209  /* set non-diagonal entries */
210  if (graph->ewgts) {
211  for (i = 0; i < nG; i++) {
212  for (e = 1; e < graph[i].nedges; e++) {
213  j = graph[i].edges[e];
214  /* conductance is 1/resistance */
215  Gm[i][j] = Gm[j][i] = -1.0 / graph[i].ewgts[e]; /* negate */
216  }
217  }
218  } else {
219  for (i = 0; i < nG; i++) {
220  for (e = 1; e < graph[i].nedges; e++) {
221  j = graph[i].edges[e];
222  /* conductance is 1/resistance */
223  Gm[i][j] = Gm[j][i] = -1.0; /* ewgts are all 1 */
224  }
225  }
226  }
227 
228  rv = solveCircuit(nG, Gm, Gm_inv);
229 
230  if (rv) {
231  float v;
232  count = 0;
233  for (i = 0; i < nG; i++) {
234  for (j = i; j < nG; j++) {
235  if (i == j)
236  v = 0.0;
237  else
238  v = (float) (Gm_inv[i][i] + Gm_inv[j][j] -
239  2.0 * Gm_inv[i][j]);
240  Dij[count++] = v;
241  }
242  }
243  } else {
244  free(Dij);
245  Dij = NULL;
246  }
247  free_array(Gm);
248  free_array(Gm_inv);
249  return Dij;
250 }
251 
252 /* sparse_stress_subspace_majorization_kD:
253  * Optimization of the stress function using sparse distance matrix, within a vector-space
254  * Fastest and least accurate method
255  *
256  * NOTE: We use integral shortest path values here, assuming
257  * this is only to get an initial layout. In general, if edge lengths
258  * are involved, we may end up with 0 length edges.
259  */
260 static int sparse_stress_subspace_majorization_kD(vtx_data * graph, /* Input graph in sparse representation */
261  int n, /* Number of nodes */
262  int nedges_graph, /* Number of edges */
263  double **coords, /* coordinates of nodes (output layout) */
264  int dim, /* dimemsionality of layout */
265  int smart_ini, /* smart initialization */
266  int exp, /* scale exponent */
267  int reweight_graph, /* difference model */
268  int n_iterations, /* max #iterations */
269  int dist_bound, /* neighborhood size in sparse distance matrix */
270  int num_centers /* #pivots in sparse distance matrix */
271  )
272 {
273  int iterations; /* output: number of iteration of the process */
274 
275  double conj_tol = tolerance_cg; /* tolerance of Conjugate Gradient */
276 
277  /*************************************************
278  ** Computation of pivot-based, sparse, subspace-restricted **
279  ** k-D stress minimization by majorization **
280  *************************************************/
281 
282  int i, j, k, node;
283 
284  /*************************************************
285  ** First compute the subspace in which we optimize **
286  ** The subspace is the high-dimensional embedding **
287  *************************************************/
288 
289  int subspace_dim = MIN(stress_pca_dim, n); /* overall dimensionality of subspace */
290  double **subspace = N_GNEW(subspace_dim, double *);
291  double *d_storage = N_GNEW(subspace_dim * n, double);
292  int num_centers_local;
293  DistType **full_coords;
294  /* if i is a pivot than CenterIndex[i] is its index, otherwise CenterIndex[i]= -1 */
295  int *CenterIndex;
296  int *invCenterIndex; /* list the pivot nodes */
297  Queue Q;
298  float *old_weights;
299  /* this matrix stores the distance between each node and each "center" */
300  DistType **Dij;
301  /* this vector stores the distances of each node to the selected "centers" */
302  DistType *dist;
303  DistType max_dist;
304  DistType *storage;
305  int *visited_nodes;
306  dist_data *distances;
307  int available_space;
308  int *storage1 = NULL;
309  DistType *storage2 = NULL;
310  int num_visited_nodes;
311  int num_neighbors;
312  int index;
313  int nedges;
314  DistType *dist_list;
315  vtx_data *lap;
316  int *edges;
317  float *ewgts;
318  double degree;
319  double **directions;
320  float **tmp_mat;
321  float **matrix;
322  double dist_ij;
323  double *b;
324  double *b_restricted;
325  double L_ij;
326  double old_stress, new_stress;
327  boolean converged;
328 
329  for (i = 0; i < subspace_dim; i++) {
330  subspace[i] = d_storage + i * n;
331  }
332 
333  /* compute PHDE: */
334  num_centers_local = MIN(n, MAX(2 * subspace_dim, 50));
335  full_coords = NULL;
336  /* High dimensional embedding */
337  embed_graph(graph, n, num_centers_local, &full_coords, reweight_graph);
338  /* Centering coordinates */
339  center_coordinate(full_coords, n, num_centers_local);
340  /* PCA */
341  PCA_alloc(full_coords, num_centers_local, n, subspace, subspace_dim);
342 
343  free(full_coords[0]);
344  free(full_coords);
345 
346  /*************************************************
347  ** Compute the sparse-shortest-distances matrix 'distances' **
348  *************************************************/
349 
350  CenterIndex = N_GNEW(n, int);
351  for (i = 0; i < n; i++) {
352  CenterIndex[i] = -1;
353  }
354  invCenterIndex = NULL;
355 
356  mkQueue(&Q, n);
357  old_weights = graph[0].ewgts;
358 
359  if (reweight_graph) {
360  /* weight graph to separate high-degree nodes */
361  /* in the future, perform slower Dijkstra-based computation */
362  compute_new_weights(graph, n);
363  }
364 
365  /* compute sparse distance matrix */
366  /* first select 'num_centers' pivots from which we compute distance */
367  /* to all other nodes */
368 
369  Dij = NULL;
370  dist = N_GNEW(n, DistType);
371  if (num_centers == 0) { /* no pivots, skip pivots-to-nodes distance calculation */
372  goto after_pivots_selection;
373  }
374 
375  invCenterIndex = N_GNEW(num_centers, int);
376 
377  storage = N_GNEW(n * num_centers, DistType);
378  Dij = N_GNEW(num_centers, DistType *);
379  for (i = 0; i < num_centers; i++)
380  Dij[i] = storage + i * n;
381 
382  /* select 'num_centers' pivots that are uniformaly spreaded over the graph */
383 
384  /* the first pivots is selected randomly */
385  node = rand() % n;
386  CenterIndex[node] = 0;
387  invCenterIndex[0] = node;
388 
389  if (reweight_graph) {
390  dijkstra(node, graph, n, Dij[0]);
391  } else {
392  bfs(node, graph, n, Dij[0], &Q);
393  }
394 
395  /* find the most distant node from first pivot */
396  max_dist = 0;
397  for (i = 0; i < n; i++) {
398  dist[i] = Dij[0][i];
399  if (dist[i] > max_dist) {
400  node = i;
401  max_dist = dist[i];
402  }
403  }
404  /* select other dim-1 nodes as pivots */
405  for (i = 1; i < num_centers; i++) {
406  CenterIndex[node] = i;
407  invCenterIndex[i] = node;
408  if (reweight_graph) {
409  dijkstra(node, graph, n, Dij[i]);
410  } else {
411  bfs(node, graph, n, Dij[i], &Q);
412  }
413  max_dist = 0;
414  for (j = 0; j < n; j++) {
415  dist[j] = MIN(dist[j], Dij[i][j]);
416  if (dist[j] > max_dist
417  || (dist[j] == max_dist && rand() % (j + 1) == 0)) {
418  node = j;
419  max_dist = dist[j];
420  }
421  }
422  }
423 
424  after_pivots_selection:
425 
426  /* Construct a sparse distance matrix 'distances' */
427 
428  /* initialize dist to -1, important for 'bfs_bounded(..)' */
429  for (i = 0; i < n; i++) {
430  dist[i] = -1;
431  }
432 
433  visited_nodes = N_GNEW(n, int);
434  distances = N_GNEW(n, dist_data);
435  available_space = 0;
436  nedges = 0;
437  for (i = 0; i < n; i++) {
438  if (CenterIndex[i] >= 0) { /* a pivot node */
439  distances[i].edges = N_GNEW(n - 1, int);
440  distances[i].edist = N_GNEW(n - 1, DistType);
441  distances[i].nedges = n - 1;
442  nedges += n - 1;
443  distances[i].free_mem = TRUE;
444  index = CenterIndex[i];
445  for (j = 0; j < i; j++) {
446  distances[i].edges[j] = j;
447  distances[i].edist[j] = Dij[index][j];
448  }
449  for (j = i + 1; j < n; j++) {
450  distances[i].edges[j - 1] = j;
451  distances[i].edist[j - 1] = Dij[index][j];
452  }
453  continue;
454  }
455 
456  /* a non pivot node */
457 
458  if (dist_bound > 0) {
459  if (reweight_graph) {
460  num_visited_nodes =
461  dijkstra_bounded(i, graph, n, dist, dist_bound,
462  visited_nodes);
463  } else {
464  num_visited_nodes =
465  bfs_bounded(i, graph, n, dist, &Q, dist_bound,
466  visited_nodes);
467  }
468  /* filter the pivots out of the visited nodes list, and the self loop: */
469  for (j = 0; j < num_visited_nodes;) {
470  if (CenterIndex[visited_nodes[j]] < 0
471  && visited_nodes[j] != i) {
472  /* not a pivot or self loop */
473  j++;
474  } else {
475  dist[visited_nodes[j]] = -1;
476  visited_nodes[j] = visited_nodes[--num_visited_nodes];
477  }
478  }
479  } else {
480  num_visited_nodes = 0;
481  }
482  num_neighbors = num_visited_nodes + num_centers;
483  if (num_neighbors > available_space) {
484  available_space = (dist_bound + 1) * n;
485  storage1 = N_GNEW(available_space, int);
486  storage2 = N_GNEW(available_space, DistType);
487  distances[i].free_mem = TRUE;
488  } else {
489  distances[i].free_mem = FALSE;
490  }
491  distances[i].edges = storage1;
492  distances[i].edist = storage2;
493  distances[i].nedges = num_neighbors;
494  nedges += num_neighbors;
495  for (j = 0; j < num_visited_nodes; j++) {
496  storage1[j] = visited_nodes[j];
497  storage2[j] = dist[visited_nodes[j]];
498  dist[visited_nodes[j]] = -1;
499  }
500  /* add all pivots: */
501  for (j = num_visited_nodes; j < num_neighbors; j++) {
502  index = j - num_visited_nodes;
503  storage1[j] = invCenterIndex[index];
504  storage2[j] = Dij[index][i];
505  }
506 
507  storage1 += num_neighbors;
508  storage2 += num_neighbors;
509  available_space -= num_neighbors;
510  }
511 
512  free(dist);
513  free(visited_nodes);
514 
515  if (Dij != NULL) {
516  free(Dij[0]);
517  free(Dij);
518  }
519 
520  /*************************************************
521  ** Laplacian computation **
522  *************************************************/
523 
524  lap = N_GNEW(n, vtx_data);
525  edges = N_GNEW(nedges + n, int);
526  ewgts = N_GNEW(nedges + n, float);
527  for (i = 0; i < n; i++) {
528  lap[i].edges = edges;
529  lap[i].ewgts = ewgts;
530  lap[i].nedges = distances[i].nedges + 1; /*add the self loop */
531  dist_list = distances[i].edist - 1; /* '-1' since edist[0] goes for number '1' entry in the lap */
532  degree = 0;
533  if (exp == 2) {
534  for (j = 1; j < lap[i].nedges; j++) {
535  edges[j] = distances[i].edges[j - 1];
536 #ifdef Dij2
537  ewgts[j] = (float) -1.0 / ((float) dist_list[j] * (float) dist_list[j]); /* cast to float to prevent overflow */
538 #else
539  ewgts[j] = -1.0 / (float) dist_list[j];
540 #endif
541  degree -= ewgts[j];
542  }
543  } else {
544  for (j = 1; j < lap[i].nedges; j++) {
545  edges[j] = distances[i].edges[j - 1];
546  ewgts[j] = -1.0 / (float) dist_list[j];
547  degree -= ewgts[j];
548  }
549  }
550  edges[0] = i;
551  ewgts[0] = (float) degree;
552  edges += lap[i].nedges;
553  ewgts += lap[i].nedges;
554  }
555 
556  /*************************************************
557  ** initialize direction vectors **
558  ** to get an intial layout **
559  *************************************************/
560 
561  /* the layout is subspace*directions */
562  directions = N_GNEW(dim, double *);
563  directions[0] = N_GNEW(dim * subspace_dim, double);
564  for (i = 1; i < dim; i++) {
565  directions[i] = directions[0] + i * subspace_dim;
566  }
567 
568  if (smart_ini) {
569  /* smart initialization */
570  for (k = 0; k < dim; k++) {
571  for (i = 0; i < subspace_dim; i++) {
572  directions[k][i] = 0;
573  }
574  }
575  if (dim != 2) {
576  /* use the first vectors in the eigenspace */
577  /* each direction points to its "principal axes" */
578  for (k = 0; k < dim; k++) {
579  directions[k][k] = 1;
580  }
581  } else {
582  /* for the frequent 2-D case we prefer iterative-PCA over PCA */
583  /* Note that we don't want to mix the Lap's eigenspace with the HDE */
584  /* in the computation since they have different scales */
585 
586  directions[0][0] = 1; /* first pca projection vector */
587  if (!iterativePCA_1D(subspace, subspace_dim, n, directions[1])) {
588  for (k = 0; k < subspace_dim; k++) {
589  directions[1][k] = 0;
590  }
591  directions[1][1] = 1;
592  }
593  }
594 
595  } else {
596  /* random initialization */
597  for (k = 0; k < dim; k++) {
598  for (i = 0; i < subspace_dim; i++) {
599  directions[k][i] = (double) (rand()) / RAND_MAX;
600  }
601  }
602  }
603 
604  /* compute initial k-D layout */
605 
606  for (k = 0; k < dim; k++) {
607  right_mult_with_vector_transpose(subspace, n, subspace_dim,
608  directions[k], coords[k]);
609  }
610 
611  /*************************************************
612  ** compute restriction of the laplacian to subspace: **
613  *************************************************/
614 
615  tmp_mat = NULL;
616  matrix = NULL;
617  mult_sparse_dense_mat_transpose(lap, subspace, n, subspace_dim,
618  &tmp_mat);
619  mult_dense_mat(subspace, tmp_mat, subspace_dim, n, subspace_dim,
620  &matrix);
621  free(tmp_mat[0]);
622  free(tmp_mat);
623 
624  /*************************************************
625  ** Layout optimization **
626  *************************************************/
627 
628  b = N_GNEW(n, double);
629  b_restricted = N_GNEW(subspace_dim, double);
630  old_stress = compute_stress1(coords, distances, dim, n, exp);
631  for (converged = FALSE, iterations = 0;
632  iterations < n_iterations && !converged; iterations++) {
633 
634  /* Axis-by-axis optimization: */
635  for (k = 0; k < dim; k++) {
636  /* compute the vector b */
637  /* multiply on-the-fly with distance-based laplacian */
638  /* (for saving storage we don't construct this Lap explicitly) */
639  for (i = 0; i < n; i++) {
640  degree = 0;
641  b[i] = 0;
642  dist_list = distances[i].edist - 1;
643  edges = lap[i].edges;
644  ewgts = lap[i].ewgts;
645  for (j = 1; j < lap[i].nedges; j++) {
646  node = edges[j];
647  dist_ij = distance_kD(coords, dim, i, node);
648  if (dist_ij > 1e-30) { /* skip zero distances */
649  L_ij = -ewgts[j] * dist_list[j] / dist_ij; /* L_ij=w_{ij}*d_{ij}/dist_{ij} */
650  degree -= L_ij;
651  b[i] += L_ij * coords[k][node];
652  }
653  }
654  b[i] += degree * coords[k][i];
655  }
656  right_mult_with_vector_d(subspace, subspace_dim, n, b,
657  b_restricted);
658  if (conjugate_gradient_f(matrix, directions[k], b_restricted,
659  subspace_dim, conj_tol, subspace_dim,
660  FALSE)) {
661  iterations = -1;
662  goto finish0;
663  }
664  right_mult_with_vector_transpose(subspace, n, subspace_dim,
665  directions[k], coords[k]);
666  }
667 
668  if ((converged = (iterations % 2 == 0))) { /* check for convergence each two iterations */
669  new_stress = compute_stress1(coords, distances, dim, n, exp);
670  converged =
671  fabs(new_stress - old_stress) / (new_stress + 1e-10) <
672  Epsilon;
673  old_stress = new_stress;
674  }
675  }
676 finish0:
677  free(b_restricted);
678  free(b);
679 
680  if (reweight_graph) {
681  restore_old_weights(graph, n, old_weights);
682  }
683 
684  for (i = 0; i < n; i++) {
685  if (distances[i].free_mem) {
686  free(distances[i].edges);
687  free(distances[i].edist);
688  }
689  }
690 
691  free(distances);
692  free(lap[0].edges);
693  free(lap[0].ewgts);
694  free(lap);
695  free(CenterIndex);
696  free(invCenterIndex);
697  free(directions[0]);
698  free(directions);
699  if (matrix != NULL) {
700  free(matrix[0]);
701  free(matrix);
702  }
703  free(subspace[0]);
704  free(subspace);
705  freeQueue(&Q);
706 
707  return iterations;
708 }
709 
710 /* compute_weighted_apsp_packed:
711  * Edge lengths can be any float > 0
712  */
713 static float *compute_weighted_apsp_packed(vtx_data * graph, int n)
714 {
715  int i, j, count;
716  float *Dij = N_NEW(n * (n + 1) / 2, float);
717 
718  float *Di = N_NEW(n, float);
719  Queue Q;
720 
721  mkQueue(&Q, n);
722 
723  count = 0;
724  for (i = 0; i < n; i++) {
725  dijkstra_f(i, graph, n, Di);
726  for (j = i; j < n; j++) {
727  Dij[count++] = Di[j];
728  }
729  }
730  free(Di);
731  freeQueue(&Q);
732  return Dij;
733 }
734 
735 
736 /* mdsModel:
737  * Update matrix with actual edge lengths
738  */
739 float *mdsModel(vtx_data * graph, int nG)
740 {
741  int i, j, e;
742  float *Dij;
743  int shift = 0;
744  double delta;
745 
746  if (graph->ewgts == NULL)
747  return 0;
748 
749  /* first, compute shortest paths to fill in non-edges */
750  Dij = compute_weighted_apsp_packed(graph, nG);
751 
752  /* then, replace edge entries will user-supplied len */
753  for (i = 0; i < nG; i++) {
754  shift += i;
755  for (e = 1; e < graph[i].nedges; e++) {
756  j = graph[i].edges[e];
757  if (j < i)
758  continue;
759  delta += abs(Dij[i * nG + j - shift] - graph[i].ewgts[e]);
760  Dij[i * nG + j - shift] = graph[i].ewgts[e];
761  }
762  }
763  if (Verbose) {
764  fprintf(stderr, "mdsModel: delta = %f\n", delta);
765  }
766  return Dij;
767 }
768 
769 /* compute_apsp_packed:
770  * Assumes integral weights > 0.
771  */
772 float *compute_apsp_packed(vtx_data * graph, int n)
773 {
774  int i, j, count;
775  float *Dij = N_NEW(n * (n + 1) / 2, float);
776 
777  DistType *Di = N_NEW(n, DistType);
778  Queue Q;
779 
780  mkQueue(&Q, n);
781 
782  count = 0;
783  for (i = 0; i < n; i++) {
784  bfs(i, graph, n, Di, &Q);
785  for (j = i; j < n; j++) {
786  Dij[count++] = ((float) Di[j]);
787  }
788  }
789  free(Di);
790  freeQueue(&Q);
791  return Dij;
792 }
793 
794 #define max(x,y) ((x)>(y)?(x):(y))
795 
797 {
798  /* compute all-pairs-shortest-path-length while weighting the graph */
799  /* so high-degree nodes are distantly located */
800 
801  float *Dij;
802  int i, j;
803  float *old_weights = graph[0].ewgts;
804  int nedges = 0;
805  float *weights;
806  int *vtx_vec;
807  int deg_i, deg_j, neighbor;
808 
809  for (i = 0; i < n; i++) {
810  nedges += graph[i].nedges;
811  }
812 
813  weights = N_NEW(nedges, float);
814  vtx_vec = N_NEW(n, int);
815  for (i = 0; i < n; i++) {
816  vtx_vec[i] = 0;
817  }
818 
819  if (graph->ewgts) {
820  for (i = 0; i < n; i++) {
821  fill_neighbors_vec_unweighted(graph, i, vtx_vec);
822  deg_i = graph[i].nedges - 1;
823  for (j = 1; j <= deg_i; j++) {
824  neighbor = graph[i].edges[j];
825  deg_j = graph[neighbor].nedges - 1;
826  weights[j] = (float)
827  max((float)
828  (deg_i + deg_j -
829  2 * common_neighbors(graph, i, neighbor,
830  vtx_vec)),
831  graph[i].ewgts[j]);
832  }
833  empty_neighbors_vec(graph, i, vtx_vec);
834  graph[i].ewgts = weights;
835  weights += graph[i].nedges;
836  }
837  Dij = compute_weighted_apsp_packed(graph, n);
838  } else {
839  for (i = 0; i < n; i++) {
840  graph[i].ewgts = weights;
841  fill_neighbors_vec_unweighted(graph, i, vtx_vec);
842  deg_i = graph[i].nedges - 1;
843  for (j = 1; j <= deg_i; j++) {
844  neighbor = graph[i].edges[j];
845  deg_j = graph[neighbor].nedges - 1;
846  weights[j] =
847  ((float) deg_i + deg_j -
848  2 * common_neighbors(graph, i, neighbor, vtx_vec));
849  }
850  empty_neighbors_vec(graph, i, vtx_vec);
851  weights += graph[i].nedges;
852  }
853  Dij = compute_apsp_packed(graph, n);
854  }
855 
856  free(vtx_vec);
857  free(graph[0].ewgts);
858  graph[0].ewgts = NULL;
859  if (old_weights != NULL) {
860  for (i = 0; i < n; i++) {
861  graph[i].ewgts = old_weights;
862  old_weights += graph[i].nedges;
863  }
864  }
865  return Dij;
866 }
867 
868 #ifdef DEBUG
869 static void dumpMatrix(float *Dij, int n)
870 {
871  int i, j, count = 0;
872  for (i = 0; i < n; i++) {
873  for (j = i; j < n; j++) {
874  fprintf(stderr, "%.02f ", Dij[count++]);
875  }
876  fputs("\n", stderr);
877  }
878 }
879 #endif
880 
881 /* Accumulator type for diagonal of Laplacian. Needs to be as large
882  * as possible. Use long double; configure to double if necessary.
883  */
884 #define DegType long double
885 
886 /* stress_majorization_kD_mkernel:
887  * At present, if any nodes have pos set, smart_ini is false.
888  */
889 int stress_majorization_kD_mkernel(vtx_data * graph, /* Input graph in sparse representation */
890  int n, /* Number of nodes */
891  int nedges_graph, /* Number of edges */
892  double **d_coords, /* coordinates of nodes (output layout) */
893  node_t ** nodes, /* original nodes */
894  int dim, /* dimemsionality of layout */
895  int opts, /* options */
896  int model, /* model */
897  int maxi /* max iterations */
898  )
899 {
900  int iterations; /* output: number of iteration of the process */
901 
902  double conj_tol = tolerance_cg; /* tolerance of Conjugate Gradient */
903  float *Dij = NULL;
904  int i, j, k;
905  float **coords = NULL;
906  float *f_storage = NULL;
907  float constant_term;
908  int count;
909  DegType degree;
910  int lap_length;
911  float *lap2 = NULL;
912  DegType *degrees = NULL;
913  int step;
914  float val;
915  double old_stress, new_stress;
916  boolean converged;
917  float **b = NULL;
918  float *tmp_coords = NULL;
919  float *dist_accumulator = NULL;
920  float *lap1 = NULL;
921  int smart_ini = opts & opt_smart_init;
922  int exp = opts & opt_exp_flag;
923  int len;
924  int havePinned; /* some node is pinned */
925 #ifdef ALTERNATIVE_STRESS_CALC
926  double mat_stress;
927 #endif
928 #ifdef NONCORE
929  FILE *fp = NULL;
930 #endif
931 
932 
933  /*************************************************
934  ** Computation of full, dense, unrestricted k-D **
935  ** stress minimization by majorization **
936  *************************************************/
937 
938  /****************************************************
939  ** Compute the all-pairs-shortest-distances matrix **
940  ****************************************************/
941 
942  if (maxi < 0)
943  return 0;
944 
945  if (Verbose)
946  start_timer();
947 
948  if (model == MODEL_SUBSET) {
949  /* weight graph to separate high-degree nodes */
950  /* and perform slower Dijkstra-based computation */
951  if (Verbose)
952  fprintf(stderr, "Calculating subset model");
954  } else if (model == MODEL_CIRCUIT) {
955  Dij = circuitModel(graph, n);
956  if (!Dij) {
957  agerr(AGWARN,
958  "graph is disconnected. Hence, the circuit model\n");
959  agerr(AGPREV,
960  "is undefined. Reverting to the shortest path model.\n");
961  }
962  } else if (model == MODEL_MDS) {
963  if (Verbose)
964  fprintf(stderr, "Calculating MDS model");
965  Dij = mdsModel(graph, n);
966  }
967  if (!Dij) {
968  if (Verbose)
969  fprintf(stderr, "Calculating shortest paths");
970  if (graph->ewgts)
971  Dij = compute_weighted_apsp_packed(graph, n);
972  else
973  Dij = compute_apsp_packed(graph, n);
974  }
975 
976  if (Verbose) {
977  fprintf(stderr, ": %.2f sec\n", elapsed_sec());
978  fprintf(stderr, "Setting initial positions");
979  start_timer();
980  }
981 
982  /**************************
983  ** Layout initialization **
984  **************************/
985 
986  if (smart_ini && (n > 1)) {
987  havePinned = 0;
988  /* optimize layout quickly within subspace */
989  /* perform at most 50 iterations within 30-D subspace to
990  get an estimate */
991  if (sparse_stress_subspace_majorization_kD(graph, n, nedges_graph,
992  d_coords, dim, smart_ini, exp,
993  (model == MODEL_SUBSET), 50,
995  num_pivots_stress) < 0) {
996  iterations = -1;
997  goto finish1;
998  }
999 
1000  for (i = 0; i < dim; i++) {
1001  /* for numerical stability, scale down layout */
1002  double max = 1;
1003  for (j = 0; j < n; j++) {
1004  if (fabs(d_coords[i][j]) > max) {
1005  max = fabs(d_coords[i][j]);
1006  }
1007  }
1008  for (j = 0; j < n; j++) {
1009  d_coords[i][j] /= max;
1010  }
1011  /* add small random noise */
1012  for (j = 0; j < n; j++) {
1013  d_coords[i][j] += 1e-6 * (drand48() - 0.5);
1014  }
1015  orthog1(n, d_coords[i]);
1016  }
1017  } else {
1018  havePinned = initLayout(graph, n, dim, d_coords, nodes);
1019  }
1020  if (Verbose)
1021  fprintf(stderr, ": %.2f sec", elapsed_sec());
1022  if ((n == 1) || (maxi == 0))
1023  return 0;
1024 
1025  if (Verbose) {
1026  fprintf(stderr, ": %.2f sec\n", elapsed_sec());
1027  fprintf(stderr, "Setting up stress function");
1028  start_timer();
1029  }
1030  coords = N_NEW(dim, float *);
1031  f_storage = N_NEW(dim * n, float);
1032  for (i = 0; i < dim; i++) {
1033  coords[i] = f_storage + i * n;
1034  for (j = 0; j < n; j++) {
1035  coords[i][j] = ((float) d_coords[i][j]);
1036  }
1037  }
1038 
1039  /* compute constant term in stress sum */
1040  /* which is \sum_{i<j} w_{ij}d_{ij}^2 */
1041  if (exp) {
1042 #ifdef Dij2
1043  constant_term = ((float) n * (n - 1) / 2);
1044 #else
1045  constant_term = 0;
1046  for (count = 0, i = 0; i < n - 1; i++) {
1047  count++; /* skip self distance */
1048  for (j = 1; j < n - i; j++, count++) {
1049  constant_term += Dij[count];
1050  }
1051  }
1052 #endif
1053  } else {
1054  constant_term = 0;
1055  for (count = 0, i = 0; i < n - 1; i++) {
1056  count++; /* skip self distance */
1057  for (j = 1; j < n - i; j++, count++) {
1058  constant_term += Dij[count];
1059  }
1060  }
1061  }
1062 
1063  /**************************
1064  ** Laplacian computation **
1065  **************************/
1066 
1067  lap_length = n * (n + 1) / 2;
1068  lap2 = Dij;
1069  if (exp == 2) {
1070 #ifdef Dij2
1071  square_vec(lap_length, lap2);
1072 #endif
1073  }
1074  /* compute off-diagonal entries */
1075  invert_vec(lap_length, lap2);
1076 
1077  /* compute diagonal entries */
1078  count = 0;
1079  degrees = N_NEW(n, DegType);
1080  /* set_vector_val(n, 0, degrees); */
1081  memset(degrees, 0, n * sizeof(DegType));
1082  for (i = 0; i < n - 1; i++) {
1083  degree = 0;
1084  count++; /* skip main diag entry */
1085  for (j = 1; j < n - i; j++, count++) {
1086  val = lap2[count];
1087  degree += val;
1088  degrees[i + j] -= val;
1089  }
1090  degrees[i] -= degree;
1091  }
1092  for (step = n, count = 0, i = 0; i < n; i++, count += step, step--) {
1093  lap2[count] = degrees[i];
1094  }
1095 
1096 #ifdef NONCORE
1097  if (n > max_nodes_in_mem) {
1098 #define FILENAME "tmp_Dij$$$.bin"
1099  fp = fopen(FILENAME, "wb");
1100  fwrite(lap2, sizeof(float), lap_length, fp);
1101  fclose(fp);
1102  fp = NULL;
1103  }
1104 #endif
1105 
1106  /*************************
1107  ** Layout optimization **
1108  *************************/
1109 
1110  b = N_NEW(dim, float *);
1111  b[0] = N_NEW(dim * n, float);
1112  for (k = 1; k < dim; k++) {
1113  b[k] = b[0] + k * n;
1114  }
1115 
1116  tmp_coords = N_NEW(n, float);
1117  dist_accumulator = N_NEW(n, float);
1118  lap1 = NULL;
1119 #ifdef NONCORE
1120  if (n <= max_nodes_in_mem) {
1121  lap1 = N_NEW(lap_length, float);
1122  } else {
1123  lap1 = lap2;
1124  fp = fopen(FILENAME, "rb");
1125  fgetpos(fp, &pos);
1126  }
1127 #else
1128  lap1 = N_NEW(lap_length, float);
1129 #endif
1130 
1131 
1132 #ifdef USE_MAXFLOAT
1133  old_stress = MAXFLOAT; /* at least one iteration */
1134 #else
1135  old_stress = MAXDOUBLE; /* at least one iteration */
1136 #endif
1137  if (Verbose) {
1138  fprintf(stderr, ": %.2f sec\n", elapsed_sec());
1139  fprintf(stderr, "Solving model: ");
1140  start_timer();
1141  }
1142 
1143  for (converged = FALSE, iterations = 0;
1144  iterations < maxi && !converged; iterations++) {
1145 
1146  /* First, construct Laplacian of 1/(d_ij*|p_i-p_j|) */
1147  /* set_vector_val(n, 0, degrees); */
1148  memset(degrees, 0, n * sizeof(DegType));
1149  if (exp == 2) {
1150 #ifdef Dij2
1151 #ifdef NONCORE
1152  if (n <= max_nodes_in_mem) {
1153  sqrt_vecf(lap_length, lap2, lap1);
1154  } else {
1155  sqrt_vec(lap_length, lap1);
1156  }
1157 #else
1158  sqrt_vecf(lap_length, lap2, lap1);
1159 #endif
1160 #endif
1161  }
1162  for (count = 0, i = 0; i < n - 1; i++) {
1163  len = n - i - 1;
1164  /* init 'dist_accumulator' with zeros */
1165  set_vector_valf(len, 0, dist_accumulator);
1166 
1167  /* put into 'dist_accumulator' all squared distances between 'i' and 'i'+1,...,'n'-1 */
1168  for (k = 0; k < dim; k++) {
1169  set_vector_valf(len, coords[k][i], tmp_coords);
1170  vectors_mult_additionf(len, tmp_coords, -1,
1171  coords[k] + i + 1);
1172  square_vec(len, tmp_coords);
1173  vectors_additionf(len, tmp_coords, dist_accumulator,
1174  dist_accumulator);
1175  }
1176 
1177  /* convert to 1/d_{ij} */
1178  invert_sqrt_vec(len, dist_accumulator);
1179  /* detect overflows */
1180  for (j = 0; j < len; j++) {
1181  if (dist_accumulator[j] >= MAXFLOAT
1182  || dist_accumulator[j] < 0) {
1183  dist_accumulator[j] = 0;
1184  }
1185  }
1186 
1187  count++; /* save place for the main diagonal entry */
1188  degree = 0;
1189  if (exp == 2) {
1190  for (j = 0; j < len; j++, count++) {
1191 #ifdef Dij2
1192  val = lap1[count] *= dist_accumulator[j];
1193 #else
1194  val = lap1[count] = dist_accumulator[j];
1195 #endif
1196  degree += val;
1197  degrees[i + j + 1] -= val;
1198  }
1199  } else {
1200  for (j = 0; j < len; j++, count++) {
1201  val = lap1[count] = dist_accumulator[j];
1202  degree += val;
1203  degrees[i + j + 1] -= val;
1204  }
1205  }
1206  degrees[i] -= degree;
1207  }
1208  for (step = n, count = 0, i = 0; i < n; i++, count += step, step--) {
1209  lap1[count] = degrees[i];
1210  }
1211 
1212  /* Now compute b[] */
1213  for (k = 0; k < dim; k++) {
1214  /* b[k] := lap1*coords[k] */
1215  right_mult_with_vector_ff(lap1, n, coords[k], b[k]);
1216  }
1217 
1218 
1219  /* compute new stress */
1220  /* remember that the Laplacians are negated, so we subtract instead of add and vice versa */
1221  new_stress = 0;
1222  for (k = 0; k < dim; k++) {
1223  new_stress += vectors_inner_productf(n, coords[k], b[k]);
1224  }
1225  new_stress *= 2;
1226  new_stress += constant_term; /* only after mult by 2 */
1227 #ifdef NONCORE
1228  if (n > max_nodes_in_mem) {
1229  /* restore lap2 from memory */
1230  fsetpos(fp, &pos);
1231  fread(lap2, sizeof(float), lap_length, fp);
1232  }
1233 #endif
1234  for (k = 0; k < dim; k++) {
1235  right_mult_with_vector_ff(lap2, n, coords[k], tmp_coords);
1236  new_stress -= vectors_inner_productf(n, coords[k], tmp_coords);
1237  }
1238 #ifdef ALTERNATIVE_STRESS_CALC
1239  mat_stress = new_stress;
1240  new_stress = compute_stressf(coords, lap2, dim, n);
1241  if (fabs(mat_stress - new_stress) / min(mat_stress, new_stress) >
1242  0.001) {
1243  fprintf(stderr, "Diff stress vals: %lf %lf (iteration #%d)\n",
1244  mat_stress, new_stress, iterations);
1245  }
1246 #endif
1247  /* Invariant: old_stress > 0. In theory, old_stress >= new_stress
1248  * but we use fabs in case of numerical error.
1249  */
1250  {
1251  double diff = old_stress - new_stress;
1252  double change = ABS(diff);
1253  converged = (((change / old_stress) < Epsilon)
1254  || (new_stress < Epsilon));
1255  }
1256  old_stress = new_stress;
1257 
1258  for (k = 0; k < dim; k++) {
1259  node_t *np;
1260  if (havePinned) {
1261  copy_vectorf(n, coords[k], tmp_coords);
1262  if (conjugate_gradient_mkernel(lap2, tmp_coords, b[k], n,
1263  conj_tol, n) < 0) {
1264  iterations = -1;
1265  goto finish1;
1266  }
1267  for (i = 0; i < n; i++) {
1268  np = nodes[i];
1269  if (!isFixed(np))
1270  coords[k][i] = tmp_coords[i];
1271  }
1272  } else {
1273  if (conjugate_gradient_mkernel(lap2, coords[k], b[k], n,
1274  conj_tol, n) < 0) {
1275  iterations = -1;
1276  goto finish1;
1277  }
1278  }
1279  }
1280  if (Verbose && (iterations % 5 == 0)) {
1281  fprintf(stderr, "%.3f ", new_stress);
1282  if ((iterations + 5) % 50 == 0)
1283  fprintf(stderr, "\n");
1284  }
1285  }
1286  if (Verbose) {
1287  fprintf(stderr, "\nfinal e = %f %d iterations %.2f sec\n",
1288  compute_stressf(coords, lap2, dim, n, exp),
1289  iterations, elapsed_sec());
1290  }
1291 
1292  for (i = 0; i < dim; i++) {
1293  for (j = 0; j < n; j++) {
1294  d_coords[i][j] = coords[i][j];
1295  }
1296  }
1297 #ifdef NONCORE
1298  if (fp)
1299  fclose(fp);
1300 #endif
1301 finish1:
1302  free(f_storage);
1303  free(coords);
1304 
1305  free(lap2);
1306  if (b) {
1307  free(b[0]);
1308  free(b);
1309  }
1310  free(tmp_coords);
1311  free(dist_accumulator);
1312  free(degrees);
1313  free(lap1);
1314  return iterations;
1315 }