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
30extern 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 */
54typedef struct {
55 int nedges;
56 int *edges;
57 DistType *edist;
58 boolean free_mem;
59} dist_data;
60
61static 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
98static double
99compute_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 */
158int
159initLayout(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
199float *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 */
260static 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 spread 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 initial 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 }
676finish0:
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 */
713static 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 */
739float *mdsModel(vtx_data * graph, int nG)
740{
741 int i, j, e;
742 float *Dij;
743 int shift = 0;
744 double delta = 0.0;
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 += fabsf(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 */
772float *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
796float *compute_apsp_artifical_weights_packed(vtx_data * graph, int n)
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#if DEBUG > 1
869static 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 */
889int 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");
953 Dij = compute_apsp_artifical_weights_packed(graph, n);
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,
994 neighborhood_radius_subspace,
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
1301finish1:
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}
1316