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#include "config.h"
15
16#include "SparseMatrix.h"
17#include "spring_electrical.h"
18#include "QuadTree.h"
19#include "Multilevel.h"
20#include "post_process.h"
21#include "overlap.h"
22#include "types.h"
23#include "memory.h"
24#include "arith.h"
25#include "logic.h"
26#include "math.h"
27#include "globals.h"
28#include <string.h>
29#include <time.h>
30
31#define PI M_PI
32
33spring_electrical_control spring_electrical_control_new(){
34 spring_electrical_control ctrl;
35 ctrl = MALLOC(sizeof(struct spring_electrical_control_struct));
36 ctrl->p = AUTOP;/*a negativve number default to -1. repulsive force = dist^p */
37 ctrl->q = 1;/*a positive number default to 1. Only apply to maxent.
38 attractive force = dist^q. Stress energy = (||x_i-x_j||-d_ij)^{q+1} */
39 ctrl->random_start = TRUE;/* whether to apply SE from a random layout, or from exisiting layout */
40 ctrl->K = -1;/* the natural distance. If K < 0, K will be set to the average distance of an edge */
41 ctrl->C = 0.2;/* another parameter. f_a(i,j) = C*dist(i,j)^2/K * d_ij, f_r(i,j) = K^(3-p)/dist(i,j)^(-p). By default C = 0.2. */
42 ctrl->multilevels = FALSE;/* if <=1, single level */
43
44 //ctrl->multilevel_coarsen_scheme = COARSEN_INDEPENDENT_EDGE_SET;
45 //ctrl->multilevel_coarsen_mode = COARSEN_MODE_GENTLE;
46
47 ctrl->multilevel_coarsen_scheme = COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_EDGE_PERNODE_SUPERNODES_FIRST; /* pass on to Multilevel_control->coarsen_scheme */
48 ctrl->multilevel_coarsen_mode = COARSEN_MODE_FORCEFUL;/*alternative: COARSEN_MODE_GENTLE. pass on to Multilevel_control->coarsen_mode */
49
50
51 ctrl->quadtree_size = 45;/* cut off size above which quadtree approximation is used */
52 ctrl->max_qtree_level = 10;/* max level of quadtree */
53 ctrl->bh = 0.6;/* Barnes-Hutt constant, if width(snode)/dist[i,snode] < bh, treat snode as a supernode.*/
54 ctrl->tol = 0.001;/* minimum different between two subsequence config before terminating. ||x-xold||_infinity < tol/K */
55 ctrl->maxiter = 500;
56 ctrl->cool = 0.90;/* default 0.9 */
57 ctrl->step = 0.1;
58 ctrl->adaptive_cooling = TRUE;
59 ctrl->random_seed = 123;
60 ctrl->beautify_leaves = FALSE;
61 ctrl->use_node_weights = FALSE;
62 ctrl->smoothing = SMOOTHING_NONE;
63 ctrl->overlap = 0;
64 ctrl->do_shrinking = 1;
65 ctrl->tscheme = QUAD_TREE_HYBRID;
66 ctrl->method = METHOD_SPRING_ELECTRICAL;
67 ctrl->initial_scaling = -4;
68 ctrl->rotation = 0.;
69 ctrl->edge_labeling_scheme = 0;
70 return ctrl;
71}
72
73void spring_electrical_control_delete(spring_electrical_control ctrl){
74 FREE(ctrl);
75}
76
77static char* smoothings[] = {
78 "NONE", "STRESS_MAJORIZATION_GRAPH_DIST", "STRESS_MAJORIZATION_AVG_DIST", "STRESS_MAJORIZATION_POWER_DIST", "SPRING", "TRIANGLE", "RNG"
79};
80
81static char* tschemes[] = {
82 "NONE", "NORMAL", "FAST", "HYBRID"
83};
84
85static char* methods[] = {
86 "SPRING_ELECTRICAL", "SPRING_MAXENT", "STRESS_MAXENT", "STRESS_APPROX", "STRESS", "UNIFORM_STRESS", "FULL_STRESS", "NONE"
87};
88
89void spring_electrical_control_print(spring_electrical_control ctrl){
90 fprintf (stderr, "spring_electrical_control:\n");
91 fprintf (stderr, " repulsive and attractive exponents: %.03f %.03f\n", ctrl->p, ctrl->q);
92 fprintf (stderr, " random start %d seed %d\n", ctrl->random_start, ctrl->random_seed);
93 fprintf (stderr, " K : %.03f C : %.03f\n", ctrl->K, ctrl->C);
94 fprintf (stderr, " max levels %d coarsen_scheme %d coarsen_node %d\n", ctrl->multilevels,
95 ctrl->multilevel_coarsen_scheme,ctrl->multilevel_coarsen_mode);
96 fprintf (stderr, " quadtree size %d max_level %d\n", ctrl->quadtree_size, ctrl->max_qtree_level);
97 fprintf (stderr, " Barnes-Hutt constant %.03f tolerance %.03f maxiter %d\n", ctrl->bh, ctrl->tol, ctrl->maxiter);
98 fprintf (stderr, " cooling %.03f step size %.03f adaptive %d\n", ctrl->cool, ctrl->step, ctrl->adaptive_cooling);
99 fprintf (stderr, " beautify_leaves %d node weights %d rotation %.03f\n",
100 ctrl->beautify_leaves, ctrl->use_node_weights, ctrl->rotation);
101 fprintf (stderr, " smoothing %s overlap %d initial_scaling %.03f do_shrinking %d\n",
102 smoothings[ctrl->smoothing], ctrl->overlap, ctrl->initial_scaling, ctrl->do_shrinking);
103 fprintf (stderr, " octree scheme %s method %s\n", tschemes[ctrl->tscheme], methods[ctrl->method]);
104 fprintf (stderr, " edge_labeling_scheme %d\n", ctrl->edge_labeling_scheme);
105}
106
107void oned_optimizer_delete(oned_optimizer opt){
108 FREE(opt);
109}
110
111oned_optimizer oned_optimizer_new(int i){
112 oned_optimizer opt;
113 opt = MALLOC(sizeof(struct oned_optimizer_struct));
114 opt->i = i;
115 opt->direction = OPT_INIT;
116 return opt;
117}
118
119void oned_optimizer_train(oned_optimizer opt, real work){
120 int i = opt->i;
121
122 assert(i >= 0);
123 opt->work[i] = work;
124 if (opt->direction == OPT_INIT){
125 if (opt->i == MAX_I){
126 opt->direction = OPT_DOWN;
127 opt->i = opt->i - 1;
128 } else {
129 opt->direction = OPT_UP;
130 opt->i = MIN(MAX_I, opt->i + 1);
131 }
132 } else if (opt->direction == OPT_UP){
133 /* fprintf(stderr, "{current_level, prev_level} = {%d,%d}, {work, work_prev} = {%f,%f}",i,i-1,opt->work[i], opt->work[i-1]);*/
134 assert(i >= 1);
135 if (opt->work[i] < opt->work[i-1] && opt->i < MAX_I){
136 /* fprintf(stderr, "keep going up to level %d\n",opt->i+1);*/
137 opt->i = MIN(MAX_I, opt->i + 1);
138 } else {
139 /* fprintf(stderr, "going down to level %d\n",opt->i-1);*/
140 (opt->i)--;
141 opt->direction = OPT_DOWN;
142 }
143 } else {
144 assert(i < MAX_I);
145 /* fprintf(stderr, "{current_level, prev_level} = {%d,%d}, {work, work_prev} = {%f,%f}",i,i+1,opt->work[i], opt->work[i+1]);*/
146 if (opt->work[i] < opt->work[i+1] && opt->i > 0){
147 /* fprintf(stderr, "keep going down to level %d\n",opt->i-1);*/
148 opt->i = MAX(0, opt->i-1);
149 } else {
150 /* fprintf(stderr, "keep up to level %d\n",opt->i+1);*/
151 (opt->i)++;
152 opt->direction = OPT_UP;
153 }
154 }
155}
156
157int oned_optimizer_get(oned_optimizer opt){
158 return opt->i;
159}
160
161
162real average_edge_length(SparseMatrix A, int dim, real *coord){
163 real dist = 0, d;
164 int *ia = A->ia, *ja = A->ja, i, j, k;
165 assert(SparseMatrix_is_symmetric(A, TRUE));
166
167 if (ia[A->m] == 0) return 1;
168 for (i = 0; i < A->m; i++){
169 for (j = ia[i]; j < ia[i+1]; j++){
170 d = 0;
171 for (k = 0; k < dim; k++){
172 d += (coord[dim*i+k] - coord[dim*ja[j]])*(coord[dim*i+k] - coord[dim*ja[j]]);
173 }
174 dist += sqrt(d);
175 }
176 }
177 return dist/ia[A->m];
178}
179
180#ifdef ENERGY
181static real spring_electrical_energy(int dim, SparseMatrix A, real *x, real p, real CRK, real KP){
182 /* 1. Grad[||x-y||^k,x] = k||x-y||^(k-1)*0.5*(x-y)/||x-y|| = k/2*||x-y||^(k-2) (x-y)
183 which should equal to -force (force = -gradient),
184 hence energy for force ||x-y||^m (y-x) is ||x-y||^(m+2)*2/(m+2) where m != 2
185 2. Grad[Log[||x-y||],x] = 1/||x-y||*0.5*(x-y)/||x-y|| = 0.5*(x-y)/||x-y||^2,
186 hence the energy to give force ||x-y||^-2 (x-y) is -2*Log[||x-y||]
187
188 */
189 int i, j, k, *ia = A->ia, *ja = A->ja, n = A->m;
190 real energy = 0, dist;
191
192 for (i = 0; i < n; i++){
193 /* attractive force C^((2-p)/3) ||x_i-x_j||/K * (x_j - x_i) */
194 for (j = ia[i]; j < ia[i+1]; j++){
195 if (ja[j] == i) continue;
196 dist = distance(x, dim, i, ja[j]);
197 energy += CRK*pow(dist, 3.)*2./3.;
198 }
199
200 /* repulsive force K^(1 - p)/||x_i-x_j||^(1 - p) (x_i - x_j) */
201 for (j = 0; j < n; j++){
202 if (j == i) continue;
203 dist = distance_cropped(x, dim, i, j);
204 for (k = 0; k < dim; k++){
205 if (p == -1){
206 energy += -KP*2*log(dist);
207 } else {
208 energy += -KP*pow(dist,p+1)*2/(p+1);
209 }
210 }
211 }
212 }
213 return energy;
214}
215
216#endif
217
218void export_embedding(FILE *fp, int dim, SparseMatrix A, real *x, real *width){
219 int i, j, k, *ia=A->ia, *ja = A->ja;
220 int ne = 0;
221 real xsize, ysize, xmin, xmax, ymin, ymax;
222
223 xmax = xmin = x[0];
224 ymax = ymin = x[1];
225 for (i = 0; i < A->m; i++){
226 xmax = MAX(xmax, x[i*dim]);
227 xmin = MIN(xmin, x[i*dim]);
228 ymax = MAX(ymax, x[i*dim+1]);
229 ymin = MIN(ymin, x[i*dim+1]);
230 }
231 xsize = xmax-xmin;
232 ysize = ymax-ymin;
233 xsize = MAX(xsize, ysize);
234
235 if (dim == 2){
236 fprintf(fp,"Graphics[{GrayLevel[0.5],Line[{");
237 } else {
238 fprintf(fp,"Graphics3D[{GrayLevel[0.5],Line[{");
239 }
240 for (i = 0; i < A->m; i++){
241 for (j = ia[i]; j < ia[i+1]; j++){
242 if (ja[j] == i) continue;
243 ne++;
244 if (ne > 1) fprintf(fp, ",");
245 fprintf(fp, "{{");
246 for (k = 0; k < dim; k++) {
247 if (k > 0) fprintf(fp,",");
248 fprintf(fp, "%f",x[i*dim+k]);
249 }
250 fprintf(fp, "},{");
251 for (k = 0; k < dim; k++) {
252 if (k > 0) fprintf(fp,",");
253 fprintf(fp, "%f",x[ja[j]*dim+k]);
254 }
255 fprintf(fp, "}}");
256 }
257 }
258
259 fprintf(fp,"}],Hue[%f]",/*drand()*/1.);
260
261 if (width && dim == 2){
262 for (i = 0; i < A->m; i++){
263 if (i >= 0) fprintf(fp,",");
264 fprintf(fp,"(*width={%f,%f}, x = {%f,%f}*){GrayLevel[.5,.5],Rectangle[{%f,%f},{%f,%f}]}", width[i*dim], width[i*dim+1], x[i*dim], x[i*dim + 1],
265 x[i*dim] - width[i*dim], x[i*dim+1] - width[i*dim+1],
266 x[i*dim] + width[i*dim], x[i*dim+1] + width[i*dim+1]);
267 }
268 }
269
270 if (A->m < 100){
271 for (i = 0; i < A->m; i++){
272 if (i >= 0) fprintf(fp,",");
273 fprintf(fp,"Text[%d,{",i+1);
274 for (k = 0; k < dim; k++) {
275 if (k > 0) fprintf(fp,",");
276 fprintf(fp, "%f",x[i*dim+k]);
277 }
278 fprintf(fp,"}]");
279 }
280 } else if (A->m < 500000){
281 fprintf(fp, ", Point[{");
282 for (i = 0; i < A->m; i++){
283 if (i > 0) fprintf(fp,",");
284 fprintf(fp,"{");
285 for (k = 0; k < dim; k++) {
286 if (k > 0) fprintf(fp,",");
287 fprintf(fp, "%f",x[i*dim+k]);
288 }
289 fprintf(fp,"}");
290 }
291 fprintf(fp, "}]");
292 } else {
293 fprintf(fp,"{}");
294 }
295
296
297 fprintf(fp,"},ImageSize->%f]\n", 2*xsize/2);
298
299}
300
301static real update_step(int adaptive_cooling, real step, real Fnorm, real Fnorm0, real cool){
302
303 if (!adaptive_cooling) {
304 return cool*step;
305 }
306 if (Fnorm >= Fnorm0){
307 step = cool*step;
308 } else if (Fnorm > 0.95*Fnorm0){
309 // step = step;
310 } else {
311 step = 0.99*step/cool;
312 }
313 return step;
314}
315
316
317#define node_degree(i) (ia[(i)+1] - ia[(i)])
318
319void check_real_array_size(real **a, int len, int *lenmax){
320 if (len >= *lenmax){
321 *lenmax = len + MAX((int) 0.2*len, 10);
322 *a = REALLOC(*a, sizeof(real)*(*lenmax));
323 }
324
325}
326void check_int_array_size(int **a, int len, int *lenmax){
327 if (len >= *lenmax){
328 *lenmax = len + MAX((int) 0.2*len, 10);
329 *a = REALLOC(*a, sizeof(int)*(*lenmax));
330 }
331
332}
333
334real get_angle(real *x, int dim, int i, int j){
335 /* between [0, 2Pi)*/
336 int k;
337 real y[2], res;
338 real eps = 0.00001;
339 for (k = 0; k < 2; k++){
340 y[k] = x[j*dim+k] - x[i*dim+k];
341 }
342 if (ABS(y[0]) <= ABS(y[1])*eps){
343 if (y[1] > 0) return 0.5*PI;
344 return 1.5*PI;
345 }
346 res = atan(y[1]/y[0]);
347 if (y[0] > 0){
348 if (y[1] < 0) res = 2*PI+res;
349 } else if (y[0] < 0){
350 res = res + PI;
351 }
352 return res;
353}
354
355int comp_real(const void *x, const void *y){
356 real *xx = (real*) x;
357 real *yy = (real*) y;
358
359 if (*xx > *yy){
360 return 1;
361 } else if (*xx < *yy){
362 return -1;
363 }
364 return 0;
365}
366static void sort_real(int n, real *a){
367 qsort(a, n, sizeof(real), comp_real);
368}
369
370
371static void set_leaves(real *x, int dim, real dist, real ang, int i, int j){
372 x[dim*j] = cos(ang)*dist + x[dim*i];
373 x[dim*j+1] = sin(ang)*dist + x[dim*i+1];
374}
375
376static void beautify_leaves(int dim, SparseMatrix A, real *x){
377 int m = A->m, i, j, *ia = A->ia, *ja = A->ja, k;
378 int *checked, p;
379 real dist;
380 int nleaves, nleaves_max = 10;
381 real *angles, maxang, ang1 = 0, ang2 = 0, pad, step;
382 int *leaves, nangles_max = 10, nangles;
383
384 assert(!SparseMatrix_has_diagonal(A));
385
386 checked = MALLOC(sizeof(int)*m);
387 angles = MALLOC(sizeof(real)*nangles_max);
388 leaves = MALLOC(sizeof(int)*nleaves_max);
389
390
391 for (i = 0; i < m; i++) checked[i] = FALSE;
392
393 for (i = 0; i < m; i++){
394 if (ia[i+1] - ia[i] != 1) continue;
395 if (checked[i]) continue;
396 p = ja[ia[i]];
397 if (!checked[p]){
398 checked[p] = TRUE;
399 dist = 0; nleaves = 0; nangles = 0;
400 for (j = ia[p]; j < ia[p+1]; j++){
401 if (node_degree(ja[j]) == 1){
402 checked[ja[j]] = TRUE;
403 check_int_array_size(&leaves, nleaves, &nleaves_max);
404 dist += distance(x, dim, p, ja[j]);
405 leaves[nleaves] = ja[j];
406 nleaves++;
407 } else {
408 check_real_array_size(&angles, nangles, &nangles_max);
409 angles[nangles++] = get_angle(x, dim, p, ja[j]);
410 }
411 }
412 assert(nleaves > 0);
413 dist /= nleaves;
414 if (nangles > 0){
415 sort_real(nangles, angles);
416 maxang = 0;
417 for (k = 0; k < nangles - 1; k++){
418 if (angles[k+1] - angles[k] > maxang){
419 maxang = angles[k+1] - angles[k];
420 ang1 = angles[k]; ang2 = angles[k+1];
421 }
422 }
423 if (2*PI + angles[0] - angles[nangles - 1] > maxang){
424 maxang = 2*PI + angles[0] - angles[nangles - 1];
425 ang1 = angles[nangles - 1];
426 ang2 = 2*PI + angles[0];
427 }
428 } else {
429 ang1 = 0; ang2 = 2*PI; maxang = 2*PI;
430 }
431 pad = MAX(maxang - PI*0.166667*(nleaves-1), 0)*0.5;
432 ang1 += pad*0.95;
433 ang2 -= pad*0.95;
434ang1 = 0; ang2 = 2*PI; maxang = 2*PI;
435 assert(ang2 >= ang1);
436 step = 0.;
437 if (nleaves > 1) step = (ang2 - ang1)/(nleaves - 1);
438 for (i = 0; i < nleaves; i++) {
439 set_leaves(x, dim, dist, ang1, p, leaves[i]);
440 ang1 += step;
441 }
442 }
443 }
444
445
446 FREE(checked);
447 FREE(angles);
448 FREE(leaves);
449}
450
451void force_print(FILE *fp, int n, int dim, real *x, real *force){
452 int i, k;
453
454 fprintf(fp,"Graphics[{");
455 for (i = 0; i < n; i++){
456 if (i > 0) fprintf(fp, ",");
457 fprintf(fp, "Arrow[{{");
458 for (k = 0; k < dim; k++){
459 if (k > 0) fprintf(fp, ",");
460 fprintf(fp, "%f",x[i*dim+k]);
461 }
462 fprintf(fp, "},{");
463 for (k = 0; k < dim; k++){
464 if (k > 0) fprintf(fp, ",");
465 fprintf(fp, "%f",x[i*dim+k]+0.5*force[i*dim+k]);
466 }
467 fprintf(fp, "}}]");
468 }
469 fprintf(fp,",");
470 for (i = 0; i < n; i++){
471 if (i > 0) fprintf(fp, ",");
472 fprintf(fp, "Tooltip[Point[{");
473 for (k = 0; k < dim; k++){
474 if (k > 0) fprintf(fp, ",");
475 fprintf(fp, "%f",x[i*dim+k]);
476 }
477 fprintf(fp, "}],%d]",i);
478 }
479
480
481
482
483 fprintf(fp,"}]\n");
484
485}
486
487
488void spring_electrical_embedding_fast(int dim, SparseMatrix A0, spring_electrical_control ctrl, real *node_weights, real *x, int *flag){
489 /* x is a point to a 1D array, x[i*dim+j] gives the coordinate of the i-th node at dimension j. */
490 SparseMatrix A = A0;
491 int m, n;
492 int i, j, k;
493 real p = ctrl->p, K = ctrl->K, C = ctrl->C, CRK, tol = ctrl->tol, maxiter = ctrl->maxiter, cool = ctrl->cool, step = ctrl->step, KP;
494 int *ia = NULL, *ja = NULL;
495 real *xold = NULL;
496 real *f = NULL, dist, F, Fnorm = 0, Fnorm0;
497 int iter = 0;
498 int adaptive_cooling = ctrl->adaptive_cooling;
499 QuadTree qt = NULL;
500 real counts[4], *force = NULL;
501#ifdef TIME
502 clock_t start, end, start0;
503 real qtree_cpu = 0, qtree_cpu0 = 0, qtree_new_cpu = 0, qtree_new_cpu0 = 0;
504 real total_cpu = 0;
505 start0 = clock();
506#endif
507 int max_qtree_level = ctrl->max_qtree_level;
508 oned_optimizer qtree_level_optimizer = NULL;
509
510 if (!A || maxiter <= 0) return;
511
512 m = A->m, n = A->n;
513 if (n <= 0 || dim <= 0) return;
514
515 qtree_level_optimizer = oned_optimizer_new(max_qtree_level);
516
517 *flag = 0;
518 if (m != n) {
519 *flag = ERROR_NOT_SQUARE_MATRIX;
520 goto RETURN;
521 }
522 assert(A->format == FORMAT_CSR);
523 A = SparseMatrix_symmetrize(A, TRUE);
524 ia = A->ia;
525 ja = A->ja;
526
527 if (ctrl->random_start){
528 srand(ctrl->random_seed);
529 for (i = 0; i < dim*n; i++) x[i] = drand();
530 }
531 if (K < 0){
532 ctrl->K = K = average_edge_length(A, dim, x);
533 }
534 if (C < 0) ctrl->C = C = 0.2;
535 if (p >= 0) ctrl->p = p = -1;
536 KP = pow(K, 1 - p);
537 CRK = pow(C, (2.-p)/3.)/K;
538
539 xold = MALLOC(sizeof(real)*dim*n);
540 force = MALLOC(sizeof(real)*dim*n);
541
542 do {
543#ifdef TIME
544 //start2 = clock();
545#endif
546
547#ifdef GVIEWER
548 if (Gviewer){
549 char *lab;
550 lab = MALLOC(sizeof(char)*100);
551 sprintf(lab,"sfdp, iter=%d", iter);
552 gviewer_set_label(lab);
553 gviewer_reset_graph_coord(A, dim, x);
554 drawScene();
555 gviewer_dump_current_frame();
556 //if ((adaptive_cooling && iter%100 == 0) || (!adaptive_cooling && iter%10 == 0)) gviewer_dump_current_frame();
557 FREE(lab);
558 }
559#endif
560
561 iter++;
562 xold = MEMCPY(xold, x, sizeof(real)*dim*n);
563 Fnorm0 = Fnorm;
564 Fnorm = 0.;
565
566 max_qtree_level = oned_optimizer_get(qtree_level_optimizer);
567
568#ifdef TIME
569 start = clock();
570#endif
571 if (ctrl->use_node_weights){
572 qt = QuadTree_new_from_point_list(dim, n, max_qtree_level, x, node_weights);
573 } else {
574 qt = QuadTree_new_from_point_list(dim, n, max_qtree_level, x, NULL);
575 }
576
577#ifdef TIME
578 qtree_new_cpu += ((real) (clock() - start))/CLOCKS_PER_SEC;
579#endif
580
581 /* repulsive force */
582#ifdef TIME
583 start = clock();
584#endif
585
586 QuadTree_get_repulsive_force(qt, force, x, ctrl->bh, p, KP, counts, flag);
587
588 assert(!(*flag));
589
590#ifdef TIME
591 end = clock();
592 qtree_cpu += ((real) (end - start)) / CLOCKS_PER_SEC;
593#endif
594
595 /* attractive force C^((2-p)/3) ||x_i-x_j||/K * (x_j - x_i) */
596 for (i = 0; i < n; i++){
597 f = &(force[i*dim]);
598 for (j = ia[i]; j < ia[i+1]; j++){
599 if (ja[j] == i) continue;
600 dist = distance(x, dim, i, ja[j]);
601 for (k = 0; k < dim; k++){
602 f[k] -= CRK*(x[i*dim+k] - x[ja[j]*dim+k])*dist;
603 }
604 }
605 }
606
607
608 /* move */
609 for (i = 0; i < n; i++){
610 f = &(force[i*dim]);
611 F = 0.;
612 for (k = 0; k < dim; k++) F += f[k]*f[k];
613 F = sqrt(F);
614 Fnorm += F;
615 if (F > 0) for (k = 0; k < dim; k++) f[k] /= F;
616 for (k = 0; k < dim; k++) x[i*dim+k] += step*f[k];
617 }/* done vertex i */
618
619
620
621 if (qt) {
622#ifdef TIME
623 start = clock();
624#endif
625 QuadTree_delete(qt);
626#ifdef TIME
627 end = clock();
628 qtree_new_cpu += ((real) (end - start)) / CLOCKS_PER_SEC;
629#endif
630
631#ifdef TIME
632 qtree_cpu0 = qtree_cpu - qtree_cpu0;
633 qtree_new_cpu0 = qtree_new_cpu - qtree_new_cpu0;
634 /* if (Verbose) fprintf(stderr, "\r iter=%d cpu=%.2f, quadtree=%.2f quad_force=%.2f other=%.2f counts={%.2f,%.2f,%.2f} step=%f Fnorm=%f nz=%d K=%f qtree_lev = %d",
635 iter, ((real) (clock() - start2)) / CLOCKS_PER_SEC, qtree_new_cpu0,
636 qtree_cpu0,((real) (clock() - start2))/CLOCKS_PER_SEC - qtree_cpu0 - qtree_new_cpu0,
637 counts[0], counts[1], counts[2],
638 step, Fnorm, A->nz,K,max_qtree_level);
639 */
640 qtree_cpu0 = qtree_cpu;
641 qtree_new_cpu0 = qtree_new_cpu;
642#endif
643 oned_optimizer_train(qtree_level_optimizer, counts[0]+0.85*counts[1]+3.3*counts[2]);
644 } else {
645 if (Verbose) {
646 fprintf(stderr, "\r iter = %d, step = %f Fnorm = %f nz = %d K = %f ",iter, step, Fnorm, A->nz,K);
647#ifdef ENERGY
648 fprintf(stderr, "energy = %f\n",spring_electrical_energy(dim, A, x, p, CRK, KP));
649#endif
650 }
651 }
652
653 step = update_step(adaptive_cooling, step, Fnorm, Fnorm0, cool);
654 } while (step > tol && iter < maxiter);
655
656#ifdef DEBUG_PRINT
657 if (Verbose && 0) fputs("\n", stderr);
658#endif
659
660
661#ifdef DEBUG_PRINT
662 if (Verbose) {
663 fprintf(stderr, "\n iter = %d, step = %f Fnorm = %f nz = %d K = %f ",iter, step, Fnorm, A->nz, K);
664 }
665#endif
666
667 if (ctrl->beautify_leaves) beautify_leaves(dim, A, x);
668
669#ifdef TIME
670 total_cpu += ((real) (clock() - start0)) / CLOCKS_PER_SEC;
671 if (Verbose) fprintf(stderr, "\n time for qtree = %f, qtree_force = %f, total cpu = %f\n",qtree_new_cpu, qtree_cpu, total_cpu);
672#endif
673
674
675 RETURN:
676 oned_optimizer_delete(qtree_level_optimizer);
677 ctrl->max_qtree_level = max_qtree_level;
678
679 if (xold) FREE(xold);
680 if (A != A0) SparseMatrix_delete(A);
681 if (force) FREE(force);
682
683}
684
685
686void spring_electrical_embedding_slow(int dim, SparseMatrix A0, spring_electrical_control ctrl, real *node_weights, real *x, int *flag){
687 /* a version that does vertex moves in one go, instead of one at a time, use for debugging the fast version. Quadtree is not used. */
688 /* x is a point to a 1D array, x[i*dim+j] gives the coordinate of the i-th node at dimension j. */
689 SparseMatrix A = A0;
690 int m, n;
691 int i, j, k;
692 real p = ctrl->p, K = ctrl->K, C = ctrl->C, CRK, tol = ctrl->tol, maxiter = ctrl->maxiter, cool = ctrl->cool, step = ctrl->step, KP;
693 int *ia = NULL, *ja = NULL;
694 real *xold = NULL;
695 real *f = NULL, dist, F, Fnorm = 0, Fnorm0;
696 int iter = 0;
697 int adaptive_cooling = ctrl->adaptive_cooling;
698 QuadTree qt = NULL;
699 int USE_QT = FALSE;
700 int nsuper = 0, nsupermax = 10;
701 real *center = NULL, *supernode_wgts = NULL, *distances = NULL, nsuper_avg, counts = 0, counts_avg = 0;
702 real *force;
703#ifdef TIME
704 clock_t start, end, start0, start2;
705 real qtree_cpu = 0, qtree_cpu0 = 0;
706 real total_cpu = 0;
707 start0 = clock();
708#endif
709 int max_qtree_level = ctrl->max_qtree_level;
710 oned_optimizer qtree_level_optimizer = NULL;
711
712 fprintf(stderr,"spring_electrical_embedding_slow");
713 if (!A || maxiter <= 0) return;
714
715 m = A->m, n = A->n;
716 if (n <= 0 || dim <= 0) return;
717 force = MALLOC(sizeof(real)*n*dim);
718
719 if (n >= ctrl->quadtree_size) {
720 USE_QT = TRUE;
721 qtree_level_optimizer = oned_optimizer_new(max_qtree_level);
722 center = MALLOC(sizeof(real)*nsupermax*dim);
723 supernode_wgts = MALLOC(sizeof(real)*nsupermax);
724 distances = MALLOC(sizeof(real)*nsupermax);
725 }
726 USE_QT = FALSE;
727 *flag = 0;
728 if (m != n) {
729 *flag = ERROR_NOT_SQUARE_MATRIX;
730 goto RETURN;
731 }
732 assert(A->format == FORMAT_CSR);
733 A = SparseMatrix_symmetrize(A, TRUE);
734 ia = A->ia;
735 ja = A->ja;
736
737 if (ctrl->random_start){
738 srand(ctrl->random_seed);
739 for (i = 0; i < dim*n; i++) x[i] = drand();
740 }
741 if (K < 0){
742 ctrl->K = K = average_edge_length(A, dim, x);
743 }
744 if (C < 0) ctrl->C = C = 0.2;
745 if (p >= 0) ctrl->p = p = -1;
746 KP = pow(K, 1 - p);
747 CRK = pow(C, (2.-p)/3.)/K;
748
749#ifdef DEBUG_0
750 {
751 FILE *f;
752 char fname[10000];
753 strcpy(fname,"/tmp/graph_layout_0_");
754 sprintf(&(fname[strlen(fname)]), "%d",n);
755 f = fopen(fname,"w");
756 export_embedding(f, dim, A, x, NULL);
757 fclose(f);
758 }
759#endif
760
761 f = MALLOC(sizeof(real)*dim);
762 xold = MALLOC(sizeof(real)*dim*n);
763 do {
764 for (i = 0; i < dim*n; i++) force[i] = 0;
765
766 iter++;
767 xold = MEMCPY(xold, x, sizeof(real)*dim*n);
768 Fnorm0 = Fnorm;
769 Fnorm = 0.;
770 nsuper_avg = 0;
771
772 if (USE_QT) {
773 max_qtree_level = oned_optimizer_get(qtree_level_optimizer);
774 if (ctrl->use_node_weights){
775 qt = QuadTree_new_from_point_list(dim, n, max_qtree_level, x, node_weights);
776 } else {
777 qt = QuadTree_new_from_point_list(dim, n, max_qtree_level, x, NULL);
778 }
779 }
780#ifdef TIME
781 start2 = clock();
782#endif
783
784
785 for (i = 0; i < n; i++){
786 for (k = 0; k < dim; k++) f[k] = 0.;
787 /* repulsive force K^(1 - p)/||x_i-x_j||^(1 - p) (x_i - x_j) */
788 if (USE_QT){
789#ifdef TIME
790 start = clock();
791#endif
792 QuadTree_get_supernodes(qt, ctrl->bh, &(x[dim*i]), i, &nsuper, &nsupermax,
793 &center, &supernode_wgts, &distances, &counts, flag);
794#ifdef TIME
795 end = clock();
796 qtree_cpu += ((real) (end - start)) / CLOCKS_PER_SEC;
797#endif
798 counts_avg += counts;
799 nsuper_avg += nsuper;
800 if (*flag) goto RETURN;
801 for (j = 0; j < nsuper; j++){
802 dist = MAX(distances[j], MINDIST);
803 for (k = 0; k < dim; k++){
804 if (p == -1){
805 f[k] += supernode_wgts[j]*KP*(x[i*dim+k] - center[j*dim+k])/(dist*dist);
806 } else {
807 f[k] += supernode_wgts[j]*KP*(x[i*dim+k] - center[j*dim+k])/pow(dist, 1.- p);
808 }
809 }
810 }
811 } else {
812 if (ctrl->use_node_weights && node_weights){
813 for (j = 0; j < n; j++){
814 if (j == i) continue;
815 dist = distance_cropped(x, dim, i, j);
816 for (k = 0; k < dim; k++){
817 if (p == -1){
818 f[k] += node_weights[j]*KP*(x[i*dim+k] - x[j*dim+k])/(dist*dist);
819 } else {
820 f[k] += node_weights[j]*KP*(x[i*dim+k] - x[j*dim+k])/pow(dist, 1.- p);
821 }
822 }
823 }
824 } else {
825 for (j = 0; j < n; j++){
826 if (j == i) continue;
827 dist = distance_cropped(x, dim, i, j);
828 for (k = 0; k < dim; k++){
829 if (p == -1){
830 f[k] += KP*(x[i*dim+k] - x[j*dim+k])/(dist*dist);
831 } else {
832 f[k] += KP*(x[i*dim+k] - x[j*dim+k])/pow(dist, 1.- p);
833 }
834 }
835 }
836 }
837 }
838 for (k = 0; k < dim; k++) force[i*dim+k] += f[k];
839 }
840
841
842
843 for (i = 0; i < n; i++){
844 for (k = 0; k < dim; k++) f[k] = 0.;
845 /* attractive force C^((2-p)/3) ||x_i-x_j||/K * (x_j - x_i) */
846 for (j = ia[i]; j < ia[i+1]; j++){
847 if (ja[j] == i) continue;
848 dist = distance(x, dim, i, ja[j]);
849 for (k = 0; k < dim; k++){
850 f[k] -= CRK*(x[i*dim+k] - x[ja[j]*dim+k])*dist;
851 }
852 }
853 for (k = 0; k < dim; k++) force[i*dim+k] += f[k];
854 }
855
856
857
858 for (i = 0; i < n; i++){
859 /* normalize force */
860 for (k = 0; k < dim; k++) f[k] = force[i*dim+k];
861
862 F = 0.;
863 for (k = 0; k < dim; k++) F += f[k]*f[k];
864 F = sqrt(F);
865 Fnorm += F;
866
867 if (F > 0) for (k = 0; k < dim; k++) f[k] /= F;
868
869 for (k = 0; k < dim; k++) x[i*dim+k] += step*f[k];
870
871 }/* done vertex i */
872
873 if (qt) {
874 QuadTree_delete(qt);
875 nsuper_avg /= n;
876 counts_avg /= n;
877#ifdef TIME
878 qtree_cpu0 = qtree_cpu - qtree_cpu0;
879 if (Verbose && 0) fprintf(stderr, "\n cpu this outer iter = %f, quadtree time = %f other time = %f\n",((real) (clock() - start2)) / CLOCKS_PER_SEC, qtree_cpu0,((real) (clock() - start2))/CLOCKS_PER_SEC - qtree_cpu0);
880 qtree_cpu0 = qtree_cpu;
881#endif
882 if (Verbose && 0) fprintf(stderr, "nsuper_avg=%f, counts_avg = %f 2*nsuper+counts=%f\n",nsuper_avg,counts_avg, 2*nsuper_avg+counts_avg);
883 oned_optimizer_train(qtree_level_optimizer, 5*nsuper_avg + counts_avg);
884 }
885
886#ifdef ENERGY
887 if (Verbose) {
888 fprintf(stderr, "\r iter = %d, step = %f Fnorm = %f nsuper = %d nz = %d K = %f ",iter, step, Fnorm, (int) nsuper_avg,A->nz,K);
889 fprintf(stderr, "energy = %f\n",spring_electrical_energy(dim, A, x, p, CRK, KP));
890 }
891#endif
892
893
894 step = update_step(adaptive_cooling, step, Fnorm, Fnorm0, cool);
895 } while (step > tol && iter < maxiter);
896
897#ifdef DEBUG_PRINT
898 if (Verbose && 0) fputs("\n", stderr);
899#endif
900
901#ifdef DEBUG_PRINT_0
902 {
903 FILE *f;
904 char fname[10000];
905 strcpy(fname,"/tmp/graph_layout");
906 sprintf(&(fname[strlen(fname)]), "%d",n);
907 f = fopen(fname,"w");
908 export_embedding(f, dim, A, x, NULL);
909 fclose(f);
910 }
911#endif
912
913
914#ifdef DEBUG_PRINT
915 if (Verbose) {
916 if (USE_QT){
917 fprintf(stderr, "iter = %d, step = %f Fnorm = %f qt_level = %d nsuper = %d nz = %d K = %f ",iter, step, Fnorm, max_qtree_level, (int) nsuper_avg,A->nz,K);
918 } else {
919 fprintf(stderr, "iter = %d, step = %f Fnorm = %f nsuper = %d nz = %d K = %f ",iter, step, Fnorm, (int) nsuper_avg,A->nz,K);
920 }
921 }
922#endif
923
924 if (ctrl->beautify_leaves) beautify_leaves(dim, A, x);
925
926#ifdef TIME
927 total_cpu += ((real) (clock() - start0)) / CLOCKS_PER_SEC;
928 if (Verbose) fprintf(stderr, "time for supernode = %f, total cpu = %f\n",qtree_cpu, total_cpu);
929#endif
930
931 RETURN:
932 if (USE_QT) {
933 oned_optimizer_delete(qtree_level_optimizer);
934 ctrl->max_qtree_level = max_qtree_level;
935 }
936 if (xold) FREE(xold);
937 if (A != A0) SparseMatrix_delete(A);
938 if (f) FREE(f);
939 if (center) FREE(center);
940 if (supernode_wgts) FREE(supernode_wgts);
941 if (distances) FREE(distances);
942 FREE(force);
943
944}
945
946
947
948void spring_electrical_embedding(int dim, SparseMatrix A0, spring_electrical_control ctrl, real *node_weights, real *x, int *flag){
949 /* x is a point to a 1D array, x[i*dim+j] gives the coordinate of the i-th node at dimension j. */
950 SparseMatrix A = A0;
951 int m, n;
952 int i, j, k;
953 real p = ctrl->p, K = ctrl->K, C = ctrl->C, CRK, tol = ctrl->tol, maxiter = ctrl->maxiter, cool = ctrl->cool, step = ctrl->step, KP;
954 int *ia = NULL, *ja = NULL;
955 real *xold = NULL;
956 real *f = NULL, dist, F, Fnorm = 0, Fnorm0;
957 int iter = 0;
958 int adaptive_cooling = ctrl->adaptive_cooling;
959 QuadTree qt = NULL;
960 int USE_QT = FALSE;
961 int nsuper = 0, nsupermax = 10;
962 real *center = NULL, *supernode_wgts = NULL, *distances = NULL, nsuper_avg, counts = 0, counts_avg = 0;
963#ifdef TIME
964 clock_t start, end, start0, start2;
965 real qtree_cpu = 0, qtree_cpu0 = 0;
966 real total_cpu = 0;
967 start0 = clock();
968#endif
969 int max_qtree_level = ctrl->max_qtree_level;
970 oned_optimizer qtree_level_optimizer = NULL;
971
972 if (!A || maxiter <= 0) return;
973
974 m = A->m, n = A->n;
975 if (n <= 0 || dim <= 0) return;
976
977 if (n >= ctrl->quadtree_size) {
978 USE_QT = TRUE;
979 qtree_level_optimizer = oned_optimizer_new(max_qtree_level);
980 center = MALLOC(sizeof(real)*nsupermax*dim);
981 supernode_wgts = MALLOC(sizeof(real)*nsupermax);
982 distances = MALLOC(sizeof(real)*nsupermax);
983 }
984 *flag = 0;
985 if (m != n) {
986 *flag = ERROR_NOT_SQUARE_MATRIX;
987 goto RETURN;
988 }
989 assert(A->format == FORMAT_CSR);
990 A = SparseMatrix_symmetrize(A, TRUE);
991 ia = A->ia;
992 ja = A->ja;
993
994 if (ctrl->random_start){
995 srand(ctrl->random_seed);
996 for (i = 0; i < dim*n; i++) x[i] = drand();
997 }
998 if (K < 0){
999 ctrl->K = K = average_edge_length(A, dim, x);
1000 }
1001 if (C < 0) ctrl->C = C = 0.2;
1002 if (p >= 0) ctrl->p = p = -1;
1003 KP = pow(K, 1 - p);
1004 CRK = pow(C, (2.-p)/3.)/K;
1005
1006#ifdef DEBUG_0
1007 {
1008 FILE *f;
1009 char fname[10000];
1010 strcpy(fname,"/tmp/graph_layout_0_");
1011 sprintf(&(fname[strlen(fname)]), "%d",n);
1012 f = fopen(fname,"w");
1013 export_embedding(f, dim, A, x, NULL);
1014 fclose(f);
1015 }
1016#endif
1017
1018 f = MALLOC(sizeof(real)*dim);
1019 xold = MALLOC(sizeof(real)*dim*n);
1020 do {
1021
1022 //#define VIS_MULTILEVEL
1023#ifdef VIS_MULTILEVEL
1024 {
1025 FILE *f;
1026 char fname[10000];
1027 static int count = 0;
1028 sprintf(fname, "/tmp/multilevel_%d",count++);
1029 f = fopen(fname,"w");
1030 export_embedding(f, dim, A, x, NULL);
1031 fclose(f);
1032 }
1033#endif
1034#ifdef GVIEWER
1035 if (Gviewer){
1036 char *lab;
1037 lab = MALLOC(sizeof(char)*100);
1038 sprintf(lab,"sfdp, adaptive_cooling = %d iter=%d", adaptive_cooling, iter);
1039 gviewer_set_label(lab);
1040 gviewer_reset_graph_coord(A, dim, x);
1041 drawScene();
1042 gviewer_dump_current_frame();
1043 //if ((adaptive_cooling && iter%100 == 0) || (!adaptive_cooling && iter%10 == 0)) gviewer_dump_current_frame();
1044 FREE(lab);
1045 }
1046#endif
1047
1048 iter++;
1049 xold = MEMCPY(xold, x, sizeof(real)*dim*n);
1050 Fnorm0 = Fnorm;
1051 Fnorm = 0.;
1052 nsuper_avg = 0;
1053 counts_avg = 0;
1054
1055 if (USE_QT) {
1056
1057 max_qtree_level = oned_optimizer_get(qtree_level_optimizer);
1058 if (ctrl->use_node_weights){
1059 qt = QuadTree_new_from_point_list(dim, n, max_qtree_level, x, node_weights);
1060 } else {
1061 qt = QuadTree_new_from_point_list(dim, n, max_qtree_level, x, NULL);
1062 }
1063
1064
1065 }
1066#ifdef TIME
1067 start2 = clock();
1068#endif
1069
1070 for (i = 0; i < n; i++){
1071 for (k = 0; k < dim; k++) f[k] = 0.;
1072 /* attractive force C^((2-p)/3) ||x_i-x_j||/K * (x_j - x_i) */
1073 for (j = ia[i]; j < ia[i+1]; j++){
1074 if (ja[j] == i) continue;
1075 dist = distance(x, dim, i, ja[j]);
1076 for (k = 0; k < dim; k++){
1077 f[k] -= CRK*(x[i*dim+k] - x[ja[j]*dim+k])*dist;
1078 }
1079 }
1080
1081 /* repulsive force K^(1 - p)/||x_i-x_j||^(1 - p) (x_i - x_j) */
1082 if (USE_QT){
1083#ifdef TIME
1084 start = clock();
1085#endif
1086 QuadTree_get_supernodes(qt, ctrl->bh, &(x[dim*i]), i, &nsuper, &nsupermax,
1087 &center, &supernode_wgts, &distances, &counts, flag);
1088
1089#ifdef TIME
1090 end = clock();
1091 qtree_cpu += ((real) (end - start)) / CLOCKS_PER_SEC;
1092#endif
1093 counts_avg += counts;
1094 nsuper_avg += nsuper;
1095 if (*flag) goto RETURN;
1096 for (j = 0; j < nsuper; j++){
1097 dist = MAX(distances[j], MINDIST);
1098 for (k = 0; k < dim; k++){
1099 if (p == -1){
1100 f[k] += supernode_wgts[j]*KP*(x[i*dim+k] - center[j*dim+k])/(dist*dist);
1101 } else {
1102 f[k] += supernode_wgts[j]*KP*(x[i*dim+k] - center[j*dim+k])/pow(dist, 1.- p);
1103 }
1104 }
1105 }
1106 } else {
1107 if (ctrl->use_node_weights && node_weights){
1108 for (j = 0; j < n; j++){
1109 if (j == i) continue;
1110 dist = distance_cropped(x, dim, i, j);
1111 for (k = 0; k < dim; k++){
1112 if (p == -1){
1113 f[k] += node_weights[j]*KP*(x[i*dim+k] - x[j*dim+k])/(dist*dist);
1114 } else {
1115 f[k] += node_weights[j]*KP*(x[i*dim+k] - x[j*dim+k])/pow(dist, 1.- p);
1116 }
1117 }
1118 }
1119 } else {
1120 for (j = 0; j < n; j++){
1121 if (j == i) continue;
1122 dist = distance_cropped(x, dim, i, j);
1123 for (k = 0; k < dim; k++){
1124 if (p == -1){
1125 f[k] += KP*(x[i*dim+k] - x[j*dim+k])/(dist*dist);
1126 } else {
1127 f[k] += KP*(x[i*dim+k] - x[j*dim+k])/pow(dist, 1.- p);
1128 }
1129 }
1130 }
1131 }
1132 }
1133
1134 /* normalize force */
1135 F = 0.;
1136 for (k = 0; k < dim; k++) F += f[k]*f[k];
1137 F = sqrt(F);
1138 Fnorm += F;
1139
1140 if (F > 0) for (k = 0; k < dim; k++) f[k] /= F;
1141
1142 for (k = 0; k < dim; k++) x[i*dim+k] += step*f[k];
1143
1144 }/* done vertex i */
1145
1146 if (qt) {
1147 QuadTree_delete(qt);
1148 nsuper_avg /= n;
1149 counts_avg /= n;
1150#ifdef TIME
1151 qtree_cpu0 = qtree_cpu - qtree_cpu0;
1152 if (Verbose && 0) fprintf(stderr, "\n cpu this outer iter = %f, quadtree time = %f other time = %f\n",((real) (clock() - start2)) / CLOCKS_PER_SEC, qtree_cpu0,((real) (clock() - start2))/CLOCKS_PER_SEC - qtree_cpu0);
1153 qtree_cpu0 = qtree_cpu;
1154#endif
1155 if (Verbose & 0) fprintf(stderr, "nsuper_avg=%f, counts_avg = %f 2*nsuper+counts=%f\n",nsuper_avg,counts_avg, 2*nsuper_avg+counts_avg);
1156 oned_optimizer_train(qtree_level_optimizer, 5*nsuper_avg + counts_avg);
1157 }
1158
1159#ifdef ENERGY
1160 if (Verbose) {
1161 fprintf(stderr, "\r iter = %d, step = %f Fnorm = %f nsuper = %d nz = %d K = %f ",iter, step, Fnorm, (int) nsuper_avg,A->nz,K);
1162 fprintf(stderr, "energy = %f\n",spring_electrical_energy(dim, A, x, p, CRK, KP));
1163 }
1164#endif
1165
1166
1167 step = update_step(adaptive_cooling, step, Fnorm, Fnorm0, cool);
1168 } while (step > tol && iter < maxiter);
1169
1170#ifdef DEBUG_PRINT
1171 if (Verbose && 0) fputs("\n", stderr);
1172#endif
1173
1174#ifdef DEBUG_PRINT_0
1175 {
1176 FILE *f;
1177 char fname[10000];
1178 strcpy(fname,"/tmp/graph_layout");
1179 sprintf(&(fname[strlen(fname)]), "%d",n);
1180 f = fopen(fname,"w");
1181 export_embedding(f, dim, A, x, NULL);
1182 fclose(f);
1183 }
1184#endif
1185
1186
1187#ifdef DEBUG_PRINT
1188 if (Verbose) {
1189 if (USE_QT){
1190 fprintf(stderr, "iter = %d, step = %f Fnorm = %f qt_level = %d nsuper = %d nz = %d K = %f ",iter, step, Fnorm, max_qtree_level, (int) nsuper_avg,A->nz,K);
1191 } else {
1192 fprintf(stderr, "iter = %d, step = %f Fnorm = %f nsuper = %d nz = %d K = %f ",iter, step, Fnorm, (int) nsuper_avg,A->nz,K);
1193 }
1194 }
1195#endif
1196
1197 if (ctrl->beautify_leaves) beautify_leaves(dim, A, x);
1198
1199#ifdef TIME
1200 total_cpu += ((real) (clock() - start0)) / CLOCKS_PER_SEC;
1201 if (Verbose) fprintf(stderr, "time for supernode = %f, total cpu = %f\n",qtree_cpu, total_cpu);
1202#endif
1203
1204 RETURN:
1205 if (USE_QT) {
1206 oned_optimizer_delete(qtree_level_optimizer);
1207 ctrl->max_qtree_level = max_qtree_level;
1208 }
1209 if (xold) FREE(xold);
1210 if (A != A0) SparseMatrix_delete(A);
1211 if (f) FREE(f);
1212 if (center) FREE(center);
1213 if (supernode_wgts) FREE(supernode_wgts);
1214 if (distances) FREE(distances);
1215
1216}
1217
1218static void scale_coord(int n, int dim, real *x, int *id, int *jd, real *d, real dj){
1219 int i, j, k;
1220 real w_ij, dist, s = 0, stop = 0, sbot = 0., nz = 0;
1221
1222 if (dj == 0.) return;
1223 for (i = 0; i < n; i++){
1224 for (j = id[i]; j < id[i+1]; j++){
1225 if (jd[j] == i) continue;
1226 dist = distance_cropped(x, dim, i, jd[j]);
1227 if (d){
1228 dj = d[j];
1229 }
1230 assert(dj > 0);
1231 w_ij = 1./(dj*dj);
1232 /* spring force */
1233 for (k = 0; k < dim; k++){
1234 stop += w_ij*dj*dist;
1235 sbot += w_ij*dist*dist;
1236 }
1237 s += dist; nz++;
1238 }
1239 }
1240 s = stop/sbot;
1241 for (i = 0; i < n*dim; i++) x[i] *= s;
1242 fprintf(stderr,"scaling factor = %f\n",s);
1243}
1244
1245static real dmean_get(int n, int *id, int *jd, real* d){
1246 real dmean = 0;
1247 int i, j;
1248
1249 if (!d) return 1.;
1250 for (i = 0; i < n; i++){
1251 for (j = id[i]; j < id[i+1]; j++){
1252 dmean += d[j];
1253 }
1254 }
1255 return dmean/((real) id[n]);
1256}
1257
1258void spring_maxent_embedding(int dim, SparseMatrix A0, SparseMatrix D, spring_electrical_control ctrl, real *node_weights, real *x, real rho, int *flag){
1259 /* x is a point to a 1D array, x[i*dim+j] gives the coordinate of the i-th node at dimension j.
1260
1261 Minimize \Sum_{(i,j)\in E} w_ij (||x_i-x_j||-d_ij)^2 - \rho \Sum_{(i,j)\NotIn E} Log ||x_i-x_j||
1262
1263 or
1264
1265 Minimize \Sum_{(i,j)\in E} w_ij (||x_i-x_j||-d_ij)^2 - \rho \Sum_{(i,j)\NotIn E} ||x_i-x_j||^p
1266
1267 The derivatives are
1268
1269 d E/d x_i = \Sum_{(i,j)\in E} w_ij (||x_i-x_j||-d_ij) (x_i-x_j)/||x_i-x_j|| - \rho \Sum_{(i,j)\NotIn E} (x_i-x_j)/||x_i-x_j||^2
1270
1271 or
1272
1273 d E/d x_i = \Sum_{(i,j)\in E} w_ij (||x_i-x_j||-d_ij) (x_i-x_j)/||x_i-x_j|| - \rho \Sum_{(i,j)\NotIn E} ||x_i-x_j||^(p-2) (x_i-x_j)
1274
1275 if D == NULL, unit weight assumed
1276
1277 */
1278 SparseMatrix A = A0;
1279 int m, n;
1280 int i, j, k;
1281 real p = ctrl->p, C = ctrl->C, tol = ctrl->tol, maxiter = ctrl->maxiter, cool = ctrl->cool, step = ctrl->step, w_ij, dj = 1.;
1282 int *ia = NULL, *ja = NULL;
1283 int *id = NULL, *jd = NULL;
1284 real *d, dmean;
1285 real *xold = NULL;
1286 real *f = NULL, dist, F, Fnorm = 0, Fnorm0;
1287 int iter = 0;
1288 int adaptive_cooling = ctrl->adaptive_cooling;
1289 QuadTree qt = NULL;
1290 int USE_QT = FALSE;
1291 int nsuper = 0, nsupermax = 10;
1292 real *center = NULL, *supernode_wgts = NULL, *distances = NULL, nsuper_avg, counts = 0;
1293 int max_qtree_level = 10;
1294#ifdef DEBUG
1295 double stress = 0;
1296#endif
1297
1298 if (!A || maxiter <= 0) return;
1299 m = A->m, n = A->n;
1300 if (n <= 0 || dim <= 0) return;
1301
1302 if (ctrl->tscheme != QUAD_TREE_NONE && n >= ctrl->quadtree_size) {
1303 USE_QT = TRUE;
1304 center = MALLOC(sizeof(real)*nsupermax*dim);
1305 supernode_wgts = MALLOC(sizeof(real)*nsupermax);
1306 distances = MALLOC(sizeof(real)*nsupermax);
1307 }
1308
1309 *flag = 0;
1310 if (m != n) {
1311 *flag = ERROR_NOT_SQUARE_MATRIX;
1312 goto RETURN;
1313 }
1314
1315
1316 assert(A->format == FORMAT_CSR);
1317 A = SparseMatrix_symmetrize(A, TRUE);
1318 ia = A->ia;
1319 ja = A->ja;
1320 if (D){
1321 id = D->ia;
1322 jd = D->ja;
1323 d = (real*) D->a;
1324 } else {
1325 id = ia; jd = ja; d = NULL;
1326 }
1327 if (rho < 0) {
1328 dmean = dmean_get(n, id, jd, d);
1329 rho = rho*(id[n]/((((real) n)*((real) n)) - id[n]))/pow(dmean, p+1);
1330 fprintf(stderr,"dmean = %f, rho = %f\n",dmean, rho);
1331 }
1332
1333 if (ctrl->random_start){
1334 fprintf(stderr, "send random coordinates\n");
1335 srand(ctrl->random_seed);
1336 for (i = 0; i < dim*n; i++) x[i] = drand();
1337 /* rescale x to give minimum stress:
1338 Min \Sum_{(i,j)\in E} w_ij (s ||x_i-x_j||-d_ij)^2
1339 thus
1340 s = (\Sum_{(ij)\in E} w_ij d_ij ||x_i-x_j||)/(\Sum_{(i,j)\in E} w_ij ||x_i-x_j||^2)
1341 */
1342
1343 }
1344 scale_coord(n, dim, x, id, jd, d, dj);
1345
1346
1347
1348 if (C < 0) ctrl->C = C = 0.2;
1349 if (p >= 0) ctrl->p = p = -1;
1350
1351#ifdef DEBUG_0
1352 {
1353 FILE *f;
1354 char fname[10000];
1355 strcpy(fname,"/tmp/graph_layout_0_");
1356 sprintf(&(fname[strlen(fname)]), "%d",n);
1357 f = fopen(fname,"w");
1358 export_embedding(f, dim, A, x, NULL);
1359 fclose(f);
1360 }
1361#endif
1362
1363 f = MALLOC(sizeof(real)*dim);
1364 xold = MALLOC(sizeof(real)*dim*n);
1365 do {
1366 iter++;
1367 xold = MEMCPY(xold, x, sizeof(real)*dim*n);
1368 Fnorm0 = Fnorm;
1369 Fnorm = 0.;
1370 nsuper_avg = 0;
1371#ifdef DEBUG
1372 stress = 0;
1373#endif
1374
1375 if (USE_QT) {
1376 if (ctrl->use_node_weights){
1377 qt = QuadTree_new_from_point_list(dim, n, max_qtree_level, x, node_weights);
1378 } else {
1379 qt = QuadTree_new_from_point_list(dim, n, max_qtree_level, x, NULL);
1380 }
1381 }
1382
1383 /*
1384 . d E/d x_i = \Sum_{(i,j)\in E} w_ij (||x_i-x_j||-d_ij) (x_i-x_j)/||x_i-x_j|| - \rho \Sum_{(i,j)\NotIn E} (x_i-x_j)/||x_i-x_j||^2
1385 or
1386 . d E/d x_i = \Sum_{(i,j)\in E} w_ij (||x_i-x_j||-d_ij) (x_i-x_j)/||x_i-x_j|| - \rho \Sum_{(i,j)\NotIn E} ||x_i-x_j||^(p-2) (x_i-x_j)
1387 */
1388 for (i = 0; i < n; i++){
1389 for (k = 0; k < dim; k++) f[k] = 0.;
1390
1391 /* spring (attractive or repulsive) force w_ij (||x_i-x_j||-d_ij) (x_i-x_j)/||x_i-x_j|| */
1392 for (j = id[i]; j < id[i+1]; j++){
1393 if (jd[j] == i) continue;
1394 dist = distance_cropped(x, dim, i, jd[j]);
1395 if (d){
1396 dj = d[j];
1397 }
1398 assert(dj > 0);
1399 /* spring force */
1400 if (ctrl->q == 2){
1401 w_ij = 1./(dj*dj*dj);
1402 for (k = 0; k < dim; k++){
1403 f[k] += -w_ij*(x[i*dim+k] - x[jd[j]*dim+k])*(dist - dj)*(dist - dj)/dist;
1404 }
1405 } else if (ctrl->q == 1){/* square stress force */
1406 w_ij = 1./(dj*dj);
1407 for (k = 0; k < dim; k++){
1408 f[k] += -w_ij*(x[i*dim+k] - x[jd[j]*dim+k])*(dist - dj)/dist;
1409 }
1410 } else {
1411 w_ij = 1./pow(dj, ctrl->q + 1);
1412 for (k = 0; k < dim; k++){
1413 f[k] += -w_ij*(x[i*dim+k] - x[jd[j]*dim+k])*pow(dist - dj, ctrl->q)/dist;
1414 }
1415 }
1416
1417#ifdef DEBUG
1418 w_ij = 1./(dj*dj);
1419 for (k = 0; k < dim; k++){
1420 stress += (dist - dj)*(dist - dj)*w_ij;
1421 }
1422#endif
1423
1424
1425 /* discount repulsive force between neighboring vertices which will be applied next, that way there is no
1426 repulsive forces between neighboring vertices */
1427 if (ctrl->use_node_weights && node_weights){
1428 for (k = 0; k < dim; k++){
1429 if (p == -1){
1430 f[k] -= rho*node_weights[j]*(x[i*dim+k] - x[jd[j]*dim+k])/(dist*dist);
1431 } else {
1432 f[k] -= rho*node_weights[j]*(x[i*dim+k] - x[jd[j]*dim+k])/pow(dist, 1.- p);
1433 }
1434 }
1435 } else {
1436 for (k = 0; k < dim; k++){
1437 if (p == -1){
1438 f[k] -= rho*(x[i*dim+k] - x[jd[j]*dim+k])/(dist*dist);
1439 } else {
1440 f[k] -= rho*(x[i*dim+k] - x[jd[j]*dim+k])/pow(dist, 1.- p);
1441 }
1442 }
1443
1444 }
1445
1446 }
1447
1448 /* repulsive force ||x_i-x_j||^(1 - p) (x_i - x_j) */
1449 if (USE_QT){
1450 QuadTree_get_supernodes(qt, ctrl->bh, &(x[dim*i]), i, &nsuper, &nsupermax,
1451 &center, &supernode_wgts, &distances, &counts, flag);
1452 nsuper_avg += nsuper;
1453 if (*flag) goto RETURN;
1454 for (j = 0; j < nsuper; j++){
1455 dist = MAX(distances[j], MINDIST);
1456 for (k = 0; k < dim; k++){
1457 if (p == -1){
1458 f[k] += rho*supernode_wgts[j]*(x[i*dim+k] - center[j*dim+k])/(dist*dist);
1459 } else {
1460 f[k] += rho*supernode_wgts[j]*(x[i*dim+k] - center[j*dim+k])/pow(dist, 1.- p);
1461 }
1462 }
1463 }
1464 } else {
1465 if (ctrl->use_node_weights && node_weights){
1466 for (j = 0; j < n; j++){
1467 if (j == i) continue;
1468 dist = distance_cropped(x, dim, i, j);
1469 for (k = 0; k < dim; k++){
1470 if (p == -1){
1471 f[k] += rho*node_weights[j]*(x[i*dim+k] - x[j*dim+k])/(dist*dist);
1472 } else {
1473 f[k] += rho*node_weights[j]*(x[i*dim+k] - x[j*dim+k])/pow(dist, 1.- p);
1474 }
1475 }
1476 }
1477 } else {
1478 for (j = 0; j < n; j++){
1479 if (j == i) continue;
1480 dist = distance_cropped(x, dim, i, j);
1481 for (k = 0; k < dim; k++){
1482 if (p == -1){
1483 f[k] += rho*(x[i*dim+k] - x[j*dim+k])/(dist*dist);
1484 } else {
1485 f[k] += rho*(x[i*dim+k] - x[j*dim+k])/pow(dist, 1.- p);
1486 }
1487 }
1488 }
1489 }
1490 }
1491
1492 /* normalize force */
1493 F = 0.;
1494 for (k = 0; k < dim; k++) F += f[k]*f[k];
1495 F = sqrt(F);
1496 Fnorm += F;
1497
1498 if (F > 0) for (k = 0; k < dim; k++) f[k] /= F;
1499
1500 for (k = 0; k < dim; k++) x[i*dim+k] += step*f[k];
1501
1502 }/* done vertex i */
1503
1504 if (qt) QuadTree_delete(qt);
1505 nsuper_avg /= n;
1506#ifdef DEBUG_PRINT
1507 stress /= (double) A->nz;
1508 if (Verbose) {
1509 fprintf(stderr, "\r iter = %d, step = %f Fnorm = %f nsuper = %d nz = %d stress = %f ",iter, step, Fnorm, (int) nsuper_avg,A->nz, stress);
1510 }
1511#endif
1512
1513 step = update_step(adaptive_cooling, step, Fnorm, Fnorm0, cool);
1514 } while (step > tol && iter < maxiter);
1515
1516#ifdef DEBUG_PRINT
1517 if (Verbose) fputs("\n", stderr);
1518#endif
1519
1520
1521 if (ctrl->beautify_leaves) beautify_leaves(dim, A, x);
1522
1523 RETURN:
1524 if (xold) FREE(xold);
1525 if (A != A0) SparseMatrix_delete(A);
1526 if (f) FREE(f);
1527 if (center) FREE(center);
1528 if (supernode_wgts) FREE(supernode_wgts);
1529 if (distances) FREE(distances);
1530
1531}
1532
1533
1534
1535void spring_electrical_spring_embedding(int dim, SparseMatrix A0, SparseMatrix D, spring_electrical_control ctrl, real *node_weights, real *x, int *flag){
1536 /* x is a point to a 1D array, x[i*dim+j] gives the coordinate of the i-th node at dimension j. Same as the spring-electrical except we also
1537 introduce force due to spring length
1538 */
1539 SparseMatrix A = A0;
1540 int m, n;
1541 int i, j, k;
1542 real p = ctrl->p, K = ctrl->K, C = ctrl->C, CRK, tol = ctrl->tol, maxiter = ctrl->maxiter, cool = ctrl->cool, step = ctrl->step, KP;
1543 int *ia = NULL, *ja = NULL;
1544 int *id = NULL, *jd = NULL;
1545 real *d;
1546 real *xold = NULL;
1547 real *f = NULL, dist, F, Fnorm = 0, Fnorm0;
1548 int iter = 0;
1549 int adaptive_cooling = ctrl->adaptive_cooling;
1550 QuadTree qt = NULL;
1551 int USE_QT = FALSE;
1552 int nsuper = 0, nsupermax = 10;
1553 real *center = NULL, *supernode_wgts = NULL, *distances = NULL, nsuper_avg, counts = 0;
1554 int max_qtree_level = 10;
1555
1556 if (!A || maxiter <= 0) return;
1557 m = A->m, n = A->n;
1558 if (n <= 0 || dim <= 0) return;
1559
1560 if (n >= ctrl->quadtree_size) {
1561 USE_QT = TRUE;
1562 center = MALLOC(sizeof(real)*nsupermax*dim);
1563 supernode_wgts = MALLOC(sizeof(real)*nsupermax);
1564 distances = MALLOC(sizeof(real)*nsupermax);
1565 }
1566 *flag = 0;
1567 if (m != n) {
1568 *flag = ERROR_NOT_SQUARE_MATRIX;
1569 goto RETURN;
1570 }
1571 assert(A->format == FORMAT_CSR);
1572 A = SparseMatrix_symmetrize(A, TRUE);
1573 ia = A->ia;
1574 ja = A->ja;
1575 id = D->ia;
1576 jd = D->ja;
1577 d = (real*) D->a;
1578
1579 if (ctrl->random_start){
1580 srand(ctrl->random_seed);
1581 for (i = 0; i < dim*n; i++) x[i] = drand();
1582 }
1583 if (K < 0){
1584 ctrl->K = K = average_edge_length(A, dim, x);
1585 }
1586 if (C < 0) ctrl->C = C = 0.2;
1587 if (p >= 0) ctrl->p = p = -1;
1588 KP = pow(K, 1 - p);
1589 CRK = pow(C, (2.-p)/3.)/K;
1590
1591#ifdef DEBUG_0
1592 {
1593 FILE *f;
1594 char fname[10000];
1595 strcpy(fname,"/tmp/graph_layout_0_");
1596 sprintf(&(fname[strlen(fname)]), "%d",n);
1597 f = fopen(fname,"w");
1598 export_embedding(f, dim, A, x, NULL);
1599 fclose(f);
1600 }
1601#endif
1602
1603 f = MALLOC(sizeof(real)*dim);
1604 xold = MALLOC(sizeof(real)*dim*n);
1605 do {
1606 iter++;
1607 xold = MEMCPY(xold, x, sizeof(real)*dim*n);
1608 Fnorm0 = Fnorm;
1609 Fnorm = 0.;
1610 nsuper_avg = 0;
1611
1612 if (USE_QT) {
1613 if (ctrl->use_node_weights){
1614 qt = QuadTree_new_from_point_list(dim, n, max_qtree_level, x, node_weights);
1615 } else {
1616 qt = QuadTree_new_from_point_list(dim, n, max_qtree_level, x, NULL);
1617 }
1618 }
1619
1620 for (i = 0; i < n; i++){
1621 for (k = 0; k < dim; k++) f[k] = 0.;
1622 /* attractive force C^((2-p)/3) ||x_i-x_j||/K * (x_j - x_i) */
1623
1624 for (j = ia[i]; j < ia[i+1]; j++){
1625 if (ja[j] == i) continue;
1626 dist = distance(x, dim, i, ja[j]);
1627 for (k = 0; k < dim; k++){
1628 f[k] -= CRK*(x[i*dim+k] - x[ja[j]*dim+k])*dist;
1629 }
1630 }
1631
1632 for (j = id[i]; j < id[i+1]; j++){
1633 if (jd[j] == i) continue;
1634 dist = distance_cropped(x, dim, i, jd[j]);
1635 for (k = 0; k < dim; k++){
1636 if (dist < d[j]){
1637 f[k] += 0.2*CRK*(x[i*dim+k] - x[jd[j]*dim+k])*(dist - d[j])*(dist - d[j])/dist;
1638 } else {
1639 f[k] -= 0.2*CRK*(x[i*dim+k] - x[jd[j]*dim+k])*(dist - d[j])*(dist - d[j])/dist;
1640 }
1641 /* f[k] -= 0.2*CRK*(x[i*dim+k] - x[jd[j]*dim+k])*(dist - d[j]);*/
1642 }
1643 }
1644
1645 /* repulsive force K^(1 - p)/||x_i-x_j||^(1 - p) (x_i - x_j) */
1646 if (USE_QT){
1647 QuadTree_get_supernodes(qt, ctrl->bh, &(x[dim*i]), i, &nsuper, &nsupermax,
1648 &center, &supernode_wgts, &distances, &counts, flag);
1649 nsuper_avg += nsuper;
1650 if (*flag) goto RETURN;
1651 for (j = 0; j < nsuper; j++){
1652 dist = MAX(distances[j], MINDIST);
1653 for (k = 0; k < dim; k++){
1654 if (p == -1){
1655 f[k] += supernode_wgts[j]*KP*(x[i*dim+k] - center[j*dim+k])/(dist*dist);
1656 } else {
1657 f[k] += supernode_wgts[j]*KP*(x[i*dim+k] - center[j*dim+k])/pow(dist, 1.- p);
1658 }
1659 }
1660 }
1661 } else {
1662 if (ctrl->use_node_weights && node_weights){
1663 for (j = 0; j < n; j++){
1664 if (j == i) continue;
1665 dist = distance_cropped(x, dim, i, j);
1666 for (k = 0; k < dim; k++){
1667 if (p == -1){
1668 f[k] += node_weights[j]*KP*(x[i*dim+k] - x[j*dim+k])/(dist*dist);
1669 } else {
1670 f[k] += node_weights[j]*KP*(x[i*dim+k] - x[j*dim+k])/pow(dist, 1.- p);
1671 }
1672 }
1673 }
1674 } else {
1675 for (j = 0; j < n; j++){
1676 if (j == i) continue;
1677 dist = distance_cropped(x, dim, i, j);
1678 for (k = 0; k < dim; k++){
1679 if (p == -1){
1680 f[k] += KP*(x[i*dim+k] - x[j*dim+k])/(dist*dist);
1681 } else {
1682 f[k] += KP*(x[i*dim+k] - x[j*dim+k])/pow(dist, 1.- p);
1683 }
1684 }
1685 }
1686 }
1687 }
1688
1689 /* normalize force */
1690 F = 0.;
1691 for (k = 0; k < dim; k++) F += f[k]*f[k];
1692 F = sqrt(F);
1693 Fnorm += F;
1694
1695 if (F > 0) for (k = 0; k < dim; k++) f[k] /= F;
1696
1697 for (k = 0; k < dim; k++) x[i*dim+k] += step*f[k];
1698
1699 }/* done vertex i */
1700
1701 if (qt) QuadTree_delete(qt);
1702 nsuper_avg /= n;
1703#ifdef DEBUG_PRINT
1704 if (Verbose && 0) {
1705 fprintf(stderr, "\r iter = %d, step = %f Fnorm = %f nsuper = %d nz = %d K = %f ",iter, step, Fnorm, (int) nsuper_avg,A->nz,K);
1706#ifdef ENERGY
1707 fprintf(stderr, "energy = %f\n",spring_electrical_energy(dim, A, x, p, CRK, KP));
1708#endif
1709 }
1710#endif
1711
1712 step = update_step(adaptive_cooling, step, Fnorm, Fnorm0, cool);
1713 } while (step > tol && iter < maxiter);
1714
1715#ifdef DEBUG_PRINT
1716 if (Verbose && 0) fputs("\n", stderr);
1717#endif
1718
1719#ifdef DEBUG_PRINT_0
1720 {
1721 FILE *f;
1722 char fname[10000];
1723 strcpy(fname,"/tmp/graph_layout");
1724 sprintf(&(fname[strlen(fname)]), "%d",n);
1725 f = fopen(fname,"w");
1726 export_embedding(f, dim, A, x, NULL);
1727 fclose(f);
1728 }
1729#endif
1730
1731 if (ctrl->beautify_leaves) beautify_leaves(dim, A, x);
1732
1733 RETURN:
1734 if (xold) FREE(xold);
1735 if (A != A0) SparseMatrix_delete(A);
1736 if (f) FREE(f);
1737 if (center) FREE(center);
1738 if (supernode_wgts) FREE(supernode_wgts);
1739 if (distances) FREE(distances);
1740
1741}
1742
1743
1744
1745
1746void print_matrix(real *x, int n, int dim){
1747 int i, k;
1748 printf("{");
1749 for (i = 0; i < n; i++){
1750 if (i != 0) printf(",");
1751 printf("{");
1752 for (k = 0; k < dim; k++) {
1753 if (k != 0) printf(",");
1754 printf("%f",x[i*dim+k]);
1755 }
1756 printf("}");
1757 }
1758 printf("}\n");
1759}
1760
1761/*
1762static void interpolate2(int dim, SparseMatrix A, real *x){
1763 int i, j, k, *ia = A->ia, *ja = A->ja, nz, m = A->m;
1764 real alpha = 0.5, beta, *y;
1765
1766 y = MALLOC(sizeof(real)*dim*m);
1767 for (k = 0; k < dim*m; k++) y[k] = 0;
1768 for (i = 0; i < m; i++){
1769 nz = 0;
1770 for (j = ia[i]; j < ia[i+1]; j++){
1771 if (ja[j] == i) continue;
1772 nz++;
1773 for (k = 0; k < dim; k++){
1774 y[i*dim+k] += x[ja[j]*dim + k];
1775 }
1776 }
1777 if (nz > 0){
1778 beta = (1-alpha)/nz;
1779 for (k = 0; k < dim; k++) y[i*dim+k] = alpha*x[i*dim+k] + beta*y[i*dim+k];
1780 }
1781 }
1782 for (k = 0; k < dim*m; k++) x[k] = y[k];
1783
1784 FREE(y);
1785}
1786
1787*/
1788
1789void interpolate_coord(int dim, SparseMatrix A, real *x){
1790 int i, j, k, *ia = A->ia, *ja = A->ja, nz;
1791 real alpha = 0.5, beta, *y;
1792
1793 y = MALLOC(sizeof(real)*dim);
1794 for (i = 0; i < A->m; i++){
1795 for (k = 0; k < dim; k++) y[k] = 0;
1796 nz = 0;
1797 for (j = ia[i]; j < ia[i+1]; j++){
1798 if (ja[j] == i) continue;
1799 nz++;
1800 for (k = 0; k < dim; k++){
1801 y[k] += x[ja[j]*dim + k];
1802 }
1803 }
1804 if (nz > 0){
1805 beta = (1-alpha)/nz;
1806 for (k = 0; k < dim; k++) x[i*dim+k] = alpha*x[i*dim+k] + beta*y[k];
1807 }
1808 }
1809
1810 FREE(y);
1811}
1812static void prolongate(int dim, SparseMatrix A, SparseMatrix P, SparseMatrix R, real *x, real *y, int coarsen_scheme_used, real delta){
1813 int nc, *ia, *ja, i, j, k;
1814 SparseMatrix_multiply_dense(P, FALSE, x, FALSE, &y, FALSE, dim);
1815
1816 /* xu yao rao dong */
1817 if (coarsen_scheme_used > EDGE_BASED_STA && coarsen_scheme_used < EDGE_BASED_STO){
1818 interpolate_coord(dim, A, y);
1819 nc = R->m;
1820 ia = R->ia;
1821 ja = R->ja;
1822 for (i = 0; i < nc; i++){
1823 for (j = ia[i]+1; j < ia[i+1]; j++){
1824 for (k = 0; k < dim; k++){
1825 y[ja[j]*dim + k] += delta*(drand() - 0.5);
1826 }
1827 }
1828 }
1829 }
1830}
1831
1832
1833
1834int power_law_graph(SparseMatrix A){
1835 int *mask, m, max = 0, i, *ia = A->ia, *ja = A->ja, j, deg;
1836 int res = FALSE;
1837 m = A->m;
1838 mask = MALLOC(sizeof(int)*(m+1));
1839
1840 for (i = 0; i < m + 1; i++){
1841 mask[i] = 0;
1842 }
1843
1844 for (i = 0; i < m; i++){
1845 deg = 0;
1846 for (j = ia[i]; j < ia[i+1]; j++){
1847 if (i == ja[j]) continue;
1848 deg++;
1849 }
1850 mask[deg]++;
1851 max = MAX(max, mask[deg]);
1852 }
1853 if (mask[1] > 0.8*max && mask[1] > 0.3*m) res = TRUE;
1854 FREE(mask);
1855 return res;
1856}
1857
1858void pcp_rotate(int n, int dim, real *x){
1859 int i, k,l;
1860 real y[4], axis[2], center[2], dist, x0, x1;
1861
1862 assert(dim == 2);
1863 for (i = 0; i < dim*dim; i++) y[i] = 0;
1864 for (i = 0; i < dim; i++) center[i] = 0;
1865 for (i = 0; i < n; i++){
1866 for (k = 0; k < dim; k++){
1867 center[k] += x[i*dim+k];
1868 }
1869 }
1870 for (i = 0; i < dim; i++) center[i] /= n;
1871 for (i = 0; i < n; i++){
1872 for (k = 0; k < dim; k++){
1873 x[dim*i+k] = x[dim*i+k] - center[k];
1874 }
1875 }
1876
1877 for (i = 0; i < n; i++){
1878 for (k = 0; k < dim; k++){
1879 for (l = 0; l < dim; l++){
1880 y[dim*k+l] += x[i*dim+k]*x[i*dim+l];
1881 }
1882 }
1883 }
1884 if (y[1] == 0) {
1885 axis[0] = 0; axis[1] = 1;
1886 } else {
1887 /* Eigensystem[{{x0, x1}, {x1, x3}}] =
1888 {{(x0 + x3 - Sqrt[x0^2 + 4*x1^2 - 2*x0*x3 + x3^2])/2,
1889 (x0 + x3 + Sqrt[x0^2 + 4*x1^2 - 2*x0*x3 + x3^2])/2},
1890 {{-(-x0 + x3 + Sqrt[x0^2 + 4*x1^2 - 2*x0*x3 + x3^2])/(2*x1), 1},
1891 {-(-x0 + x3 - Sqrt[x0^2 + 4*x1^2 - 2*x0*x3 + x3^2])/(2*x1), 1}}}
1892 */
1893 axis[0] = -(-y[0] + y[3] - sqrt(y[0]*y[0]+4*y[1]*y[1]-2*y[0]*y[3]+y[3]*y[3]))/(2*y[1]);
1894 axis[1] = 1;
1895 }
1896 dist = sqrt(1+axis[0]*axis[0]);
1897 axis[0] = axis[0]/dist;
1898 axis[1] = axis[1]/dist;
1899 for (i = 0; i < n; i++){
1900 x0 = x[dim*i]*axis[0]+x[dim*i+1]*axis[1];
1901 x1 = -x[dim*i]*axis[1]+x[dim*i+1]*axis[0];
1902 x[dim*i] = x0;
1903 x[dim*i + 1] = x1;
1904
1905 }
1906
1907
1908}
1909
1910static void rotate(int n, int dim, real *x, real angle){
1911 int i, k;
1912 real axis[2], center[2], x0, x1;
1913 real radian = 3.14159/180;
1914
1915 assert(dim == 2);
1916 for (i = 0; i < dim; i++) center[i] = 0;
1917 for (i = 0; i < n; i++){
1918 for (k = 0; k < dim; k++){
1919 center[k] += x[i*dim+k];
1920 }
1921 }
1922 for (i = 0; i < dim; i++) center[i] /= n;
1923 for (i = 0; i < n; i++){
1924 for (k = 0; k < dim; k++){
1925 x[dim*i+k] = x[dim*i+k] - center[k];
1926 }
1927 }
1928 axis[0] = cos(-angle*radian);
1929 axis[1] = sin(-angle*radian);
1930 for (i = 0; i < n; i++){
1931 x0 = x[dim*i]*axis[0]+x[dim*i+1]*axis[1];
1932 x1 = -x[dim*i]*axis[1]+x[dim*i+1]*axis[0];
1933 x[dim*i] = x0;
1934 x[dim*i + 1] = x1;
1935 }
1936
1937
1938}
1939
1940static void attach_edge_label_coordinates(int dim, SparseMatrix A, int n_edge_label_nodes, int *edge_label_nodes, real *x, real *x2){
1941 int *mask;
1942 int i, ii, j, k;
1943 int nnodes = 0;
1944 real len;
1945
1946 mask = MALLOC(sizeof(int)*A->m);
1947
1948 for (i = 0; i < A->m; i++) mask[i] = 1;
1949 for (i = 0; i < n_edge_label_nodes; i++) {
1950 if (edge_label_nodes[i] >= 0 && edge_label_nodes[i] < A->m) mask[edge_label_nodes[i]] = -1;
1951 }
1952
1953 for (i = 0; i < A->m; i++) {
1954 if (mask[i] >= 0) mask[i] = nnodes++;
1955 }
1956
1957
1958 for (i = 0; i < A->m; i++){
1959 if (mask[i] >= 0){
1960 for (k = 0; k < dim; k++) x[i*dim+k] = x2[mask[i]*dim+k];
1961 }
1962 }
1963
1964 for (i = 0; i < n_edge_label_nodes; i++){
1965 ii = edge_label_nodes[i];
1966 len = A->ia[ii+1] - A->ia[ii];
1967 assert(len >= 2); /* should just be 2 */
1968 assert(mask[ii] < 0);
1969 for (k = 0; k < dim; k++) {
1970 x[ii*dim+k] = 0;
1971 }
1972 for (j = A->ia[ii]; j < A->ia[ii+1]; j++){
1973 for (k = 0; k < dim; k++){
1974 x[ii*dim+k] += x[(A->ja[j])*dim+k];
1975 }
1976 }
1977 for (k = 0; k < dim; k++) {
1978 x[ii*dim+k] /= len;
1979 }
1980 }
1981
1982 FREE(mask);
1983}
1984
1985static SparseMatrix shorting_edge_label_nodes(SparseMatrix A, int n_edge_label_nodes, int *edge_label_nodes){
1986 int *mask;
1987 int i, id = 0, nz, j, jj, ii;
1988 int *ia = A->ia, *ja = A->ja, *irn = NULL, *jcn = NULL;
1989 SparseMatrix B;
1990
1991 mask = MALLOC(sizeof(int)*A->m);
1992
1993 for (i = 0; i < A->m; i++) mask[i] = 1;
1994
1995 for (i = 0; i < n_edge_label_nodes; i++){
1996 mask[edge_label_nodes[i]] = -1;
1997 }
1998
1999 for (i = 0; i < A->m; i++) {
2000 if (mask[i] > 0) mask[i] = id++;
2001 }
2002
2003 nz = 0;
2004 for (i = 0; i < A->m; i++){
2005 if (mask[i] < 0) continue;
2006 for (j = ia[i]; j < ia[i+1]; j++){
2007 if (mask[ja[j]] >= 0) {
2008 nz++;
2009 continue;
2010 }
2011 ii = ja[j];
2012 for (jj = ia[ii]; jj < ia[ii+1]; jj++){
2013 if (ja[jj] != i && mask[ja[jj]] >= 0) nz++;
2014 }
2015 }
2016 }
2017
2018 if (nz > 0) {
2019 irn = MALLOC(sizeof(int)*nz);
2020 jcn = MALLOC(sizeof(int)*nz);
2021 }
2022
2023 nz = 0;
2024 for (i = 0; i < A->m; i++){
2025 if (mask[i] < 0) continue;
2026 for (j = ia[i]; j < ia[i+1]; j++){
2027 if (mask[ja[j]] >= 0) {
2028 irn[nz] = mask[i];
2029 jcn[nz++] = mask[ja[j]];
2030 continue;
2031 }
2032 ii = ja[j];
2033 for (jj = ia[ii]; jj < ia[ii+1]; jj++){
2034 if (ja[jj] != i && mask[ja[jj]] >= 0) {
2035 irn[nz] = mask[i];
2036 jcn[nz++] = mask[ja[jj]];
2037 if (mask[i] == 68 || mask[ja[jj]] == 68){
2038 fprintf(stderr, "%d %d\n",mask[i], mask[ja[jj]]);
2039 mask[i] = mask[i];
2040 }
2041 }
2042 }
2043 }
2044 }
2045
2046 B = SparseMatrix_from_coordinate_arrays(nz, id, id, irn, jcn, NULL, MATRIX_TYPE_PATTERN, sizeof(real));
2047
2048 FREE(irn);
2049 FREE(jcn);
2050 FREE(mask);
2051 return B;
2052
2053}
2054
2055static void multilevel_spring_electrical_embedding_core(int dim, SparseMatrix A0, SparseMatrix D0, spring_electrical_control ctrl, real *node_weights, real *label_sizes,
2056 real *x, int n_edge_label_nodes, int *edge_label_nodes, int *flag){
2057
2058
2059 Multilevel_control mctrl = NULL;
2060 int n, plg, coarsen_scheme_used;
2061 SparseMatrix A = A0, D = D0, P = NULL;
2062 Multilevel grid, grid0;
2063 real *xc = NULL, *xf = NULL;
2064 struct spring_electrical_control_struct ctrl0;
2065#ifdef TIME
2066 clock_t cpu;
2067#endif
2068
2069 ctrl0 = *ctrl;
2070
2071#ifdef TIME
2072 cpu = clock();
2073#endif
2074
2075 *flag = 0;
2076 if (!A) return;
2077 n = A->n;
2078 if (n <= 0 || dim <= 0) return;
2079
2080 if (!SparseMatrix_is_symmetric(A, FALSE) || A->type != MATRIX_TYPE_REAL){
2081 if (ctrl->method == METHOD_SPRING_MAXENT){
2082 A = SparseMatrix_symmetrize_nodiag(A, FALSE);
2083 assert(D0);
2084 D = SparseMatrix_symmetrize_nodiag(D, FALSE);
2085 } else {
2086 A = SparseMatrix_get_real_adjacency_matrix_symmetrized(A);
2087 }
2088 } else {
2089 if (ctrl->method == METHOD_SPRING_MAXENT){
2090 assert(D0);
2091 D = SparseMatrix_remove_diagonal(D);
2092 }
2093 A = SparseMatrix_remove_diagonal(A);
2094 }
2095
2096 /* we first generate a layout discarding (shorting) the edge labels nodes, then assign the edge label nodes at the average of their neighbors */
2097 if ((ctrl->edge_labeling_scheme == ELSCHEME_STRAIGHTLINE_PENALTY || ctrl->edge_labeling_scheme == ELSCHEME_STRAIGHTLINE_PENALTY2)
2098 && n_edge_label_nodes > 0){
2099 SparseMatrix A2;
2100
2101 real *x2 = MALLOC(sizeof(real)*(A->m)*dim);
2102 A2 = shorting_edge_label_nodes(A, n_edge_label_nodes, edge_label_nodes);
2103 multilevel_spring_electrical_embedding(dim, A2, NULL, ctrl, NULL, NULL, x2, 0, NULL, flag);
2104
2105 assert(!(*flag));
2106 attach_edge_label_coordinates(dim, A, n_edge_label_nodes, edge_label_nodes, x, x2);
2107 remove_overlap(dim, A, x, label_sizes, ctrl->overlap, ctrl->initial_scaling,
2108 ctrl->edge_labeling_scheme, n_edge_label_nodes, edge_label_nodes, A, ctrl->do_shrinking, flag);
2109 SparseMatrix_delete(A2);
2110 FREE(x2);
2111 if (A != A0) SparseMatrix_delete(A);
2112
2113 return;
2114 }
2115
2116 mctrl = Multilevel_control_new(ctrl->multilevel_coarsen_scheme, ctrl->multilevel_coarsen_mode);
2117 mctrl->maxlevel = ctrl->multilevels;
2118 grid0 = Multilevel_new(A, D, node_weights, mctrl);
2119
2120 grid = Multilevel_get_coarsest(grid0);
2121 if (Multilevel_is_finest(grid)){
2122 xc = x;
2123 } else {
2124 xc = MALLOC(sizeof(real)*grid->n*dim);
2125 }
2126
2127 plg = power_law_graph(A);
2128 if (ctrl->p == AUTOP){
2129 ctrl->p = -1;
2130 if (plg) ctrl->p = -1.8;
2131 }
2132
2133 do {
2134#ifdef DEBUG_PRINT
2135 if (Verbose) {
2136 print_padding(grid->level);
2137 if (Multilevel_is_coarsest(grid)){
2138 fprintf(stderr, "coarsest level -- %d, n = %d\n", grid->level, grid->n);
2139 } else {
2140 fprintf(stderr, "level -- %d, n = %d\n", grid->level, grid->n);
2141 }
2142 }
2143#endif
2144 if (ctrl->method == METHOD_SPRING_ELECTRICAL){
2145 if (ctrl->tscheme == QUAD_TREE_NONE){
2146 spring_electrical_embedding_slow(dim, grid->A, ctrl, grid->node_weights, xc, flag);
2147 } else if (ctrl->tscheme == QUAD_TREE_FAST || (ctrl->tscheme == QUAD_TREE_HYBRID && grid->A->m > QUAD_TREE_HYBRID_SIZE)){
2148 if (ctrl->tscheme == QUAD_TREE_HYBRID && grid->A->m > 10 && Verbose){
2149 fprintf(stderr, "QUAD_TREE_HYBRID, size larger than %d, switch to fast quadtree", QUAD_TREE_HYBRID_SIZE);
2150 }
2151 spring_electrical_embedding_fast(dim, grid->A, ctrl, grid->node_weights, xc, flag);
2152 } else if (ctrl->tscheme == QUAD_TREE_NORMAL){
2153 spring_electrical_embedding(dim, grid->A, ctrl, grid->node_weights, xc, flag);
2154 } else {
2155 spring_electrical_embedding(dim, grid->A, ctrl, grid->node_weights, xc, flag);
2156 }
2157 } else if (ctrl->method == METHOD_SPRING_MAXENT){
2158 double rho = 0.05;
2159
2160 ctrl->step = 1;
2161 ctrl->adaptive_cooling = TRUE;
2162 if (Multilevel_is_coarsest(grid)){
2163 ctrl->maxiter=500;
2164 rho = 0.5;
2165 } else {
2166 ctrl->maxiter=100;
2167 }
2168
2169 if (Multilevel_is_finest(grid)) {/* gradually reduce influence of entropy */
2170 spring_maxent_embedding(dim, grid->A, grid->D, ctrl, grid->node_weights, xc, rho, flag);
2171 ctrl->random_start = FALSE;
2172 ctrl->step = .05;
2173 ctrl->adaptive_cooling = FALSE;
2174 spring_maxent_embedding(dim, grid->A, grid->D, ctrl, grid->node_weights, xc, rho/2, flag);
2175 spring_maxent_embedding(dim, grid->A, grid->D, ctrl, grid->node_weights, xc, rho/8, flag);
2176 spring_maxent_embedding(dim, grid->A, grid->D, ctrl, grid->node_weights, xc, rho/32, flag);
2177 } else {
2178 spring_maxent_embedding(dim, grid->A, grid->D, ctrl, grid->node_weights, xc, rho, flag);
2179 }
2180 } else {
2181 assert(0);
2182 }
2183 if (Multilevel_is_finest(grid)) break;
2184 if (*flag) {
2185 FREE(xc);
2186 goto RETURN;
2187 }
2188 P = grid->P;
2189 coarsen_scheme_used = grid->coarsen_scheme_used;
2190 grid = grid->prev;
2191 if (Multilevel_is_finest(grid)){
2192 xf = x;
2193 } else {
2194 xf = MALLOC(sizeof(real)*grid->n*dim);
2195 }
2196 prolongate(dim, grid->A, P, grid->R, xc, xf, coarsen_scheme_used, (ctrl->K)*0.001);
2197 FREE(xc);
2198 xc = xf;
2199 ctrl->random_start = FALSE;
2200 ctrl->K = ctrl->K * 0.75;
2201 ctrl->adaptive_cooling = FALSE;
2202 if (grid->next->coarsen_scheme_used > VERTEX_BASED_STA &&
2203 grid->next->coarsen_scheme_used < VERTEX_BASED_STO){
2204 ctrl->step = 1;
2205 } else {
2206 ctrl->step = .1;
2207 }
2208 } while (grid);
2209
2210#ifdef TIME
2211 if (Verbose)
2212 fprintf(stderr, "layout time %f\n",((real) (clock() - cpu)) / CLOCKS_PER_SEC);
2213 cpu = clock();
2214#endif
2215
2216 post_process_smoothing(dim, A, ctrl, node_weights, x, flag);
2217
2218 if (Verbose) fprintf(stderr, "ctrl->overlap=%d\n",ctrl->overlap);
2219
2220 /* rotation has to be done before overlap removal, since rotation could induce overlaps */
2221 if (dim == 2){
2222 pcp_rotate(n, dim, x);
2223 }
2224 if (ctrl->rotation != 0) rotate(n, dim, x, ctrl->rotation);
2225
2226
2227 remove_overlap(dim, A, x, label_sizes, ctrl->overlap, ctrl->initial_scaling,
2228 ctrl->edge_labeling_scheme, n_edge_label_nodes, edge_label_nodes, A, ctrl->do_shrinking, flag);
2229
2230 RETURN:
2231 *ctrl = ctrl0;
2232 if (A != A0) SparseMatrix_delete(A);
2233 if (D && D != D0) SparseMatrix_delete(D);
2234 Multilevel_control_delete(mctrl);
2235 Multilevel_delete(grid0);
2236}
2237
2238#ifdef GVIEWER
2239struct multilevel_spring_electrical_embedding_data {
2240 int dim;
2241 SparseMatrix A;
2242 SparseMatrix D;
2243 spring_electrical_control ctrl;
2244 real *node_weights;
2245 real *label_sizes;
2246 real *x;
2247 int n_edge_label_nodes;
2248 int *edge_label_nodes;
2249 int *flag;
2250};
2251
2252void multilevel_spring_electrical_embedding_gv(void* data){
2253 struct multilevel_spring_electrical_embedding_data* d;
2254
2255 d = (struct multilevel_spring_electrical_embedding_data*) data;
2256 multilevel_spring_electrical_embedding_core(d->dim, d->A, d->D, d->ctrl, d->node_weights, d->label_sizes, d->x, d->n_edge_label_nodes, d->edge_label_nodes, d->flag);
2257 gviewer_reset_graph_coord(d->A, d->dim, d->x);/* A inside spring_electrical gets deleted */
2258}
2259void multilevel_spring_electrical_embedding(int dim, SparseMatrix A, SparseMatrix D, spring_electrical_control ctrl, real *node_weights, real *label_sizes,
2260 real *x, int n_edge_label_nodes, int *edge_label_nodes, int *flag){
2261 struct multilevel_spring_electrical_embedding_data data = {dim, A, D, ctrl, node_weights, label_sizes, x, n_edge_label_nodes, edge_label_nodes, flag};
2262
2263 int argcc = 1;
2264 char **argvv;
2265
2266 if (!Gviewer) return multilevel_spring_electrical_embedding_core(dim, A, D, ctrl, node_weights, label_sizes, x, n_edge_label_nodes, edge_label_nodes, flag);
2267
2268 argcc = 1;
2269 argvv = malloc(sizeof(char*)*argcc);
2270 argvv[0] = malloc(sizeof(char));
2271 argvv[0][0] = '1';
2272
2273 gviewer_set_edge_color_scheme(COLOR_SCHEME_NO);
2274 //gviewer_set_edge_color_scheme(COLOR_SCHEME_MEDIAN_AS_GREEN);
2275 gviewer_toggle_bgcolor();
2276 //gviewer_toggle_vertex();
2277 //gviewer_init(&argcc, argvv, 0.01, 20, 60, 2*1010, 2*770, A, dim, x, &(data), multilevel_spring_electrical_embedding_gv);
2278 gviewer_init(&argcc, argvv, 0.01, 20, 60, 320, 320, A, dim, x, &(data), multilevel_spring_electrical_embedding_gv);
2279 free(argvv);
2280
2281}
2282#else
2283void multilevel_spring_electrical_embedding(int dim, SparseMatrix A, SparseMatrix D, spring_electrical_control ctrl, real *node_weights, real *label_sizes,
2284 real *x, int n_edge_label_nodes, int *edge_label_nodes, int *flag){
2285 multilevel_spring_electrical_embedding_core(dim, A, D, ctrl, node_weights, label_sizes, x, n_edge_label_nodes, edge_label_nodes, flag);
2286}
2287#endif
2288