Graphviz  2.41.20170921.2350
spring_electrical.c
Go to the documentation of this file.
1 /* $Id$Revision: */
2 /* vim:set shiftwidth=4 ts=8: */
3 
4 /*************************************************************************
5  * Copyright (c) 2011 AT&T Intellectual Property
6  * All rights reserved. This program and the accompanying materials
7  * are made available under the terms of the Eclipse Public License v1.0
8  * which accompanies this distribution, and is available at
9  * http://www.eclipse.org/legal/epl-v10.html
10  *
11  * Contributors: See CVS logs. Details at http://www.graphviz.org/
12  *************************************************************************/
13 
14 #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 
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;
67  ctrl->initial_scaling = -4;
68  ctrl->rotation = 0.;
69  ctrl->edge_labeling_scheme = 0;
70  return ctrl;
71 }
72 
74  FREE(ctrl);
75 }
76 
77 static char* smoothings[] = {
78  "NONE", "STRESS_MAJORIZATION_GRAPH_DIST", "STRESS_MAJORIZATION_AVG_DIST", "STRESS_MAJORIZATION_POWER_DIST", "SPRING", "TRIANGLE", "RNG"
79 };
80 
81 static char* tschemes[] = {
82  "NONE", "NORMAL", "FAST", "HYBRID"
83 };
84 
85 static char* methods[] = {
86  "SPRING_ELECTRICAL", "SPRING_MAXENT", "STRESS_MAXENT", "STRESS_APPROX", "STRESS", "UNIFORM_STRESS", "FULL_STRESS", "NONE"
87 };
88 
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,
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 
108  FREE(opt);
109 }
110 
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 
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 
158  return opt->i;
159 }
160 
161 
163  real dist = 0, d;
164  int *ia = A->ia, *ja = A->ja, i, j, k;
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
181 static 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 
218 void 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 
301 static 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 
319 void 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 }
326 void 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 
334 real 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 
355 int 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 }
366 static void sort_real(int n, real *a){
367  qsort(a, n, sizeof(real), comp_real);
368 }
369 
370 
371 static 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 
376 static 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 
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;
434 ang1 = 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 
451 void 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 
488 void 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);
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 
686 void 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);
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 
948 void 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);
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 
1218 static 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 
1245 static 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 
1258 void 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 
1535 void 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 
1746 void 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 /*
1762 static 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 
1789 void 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 }
1812 static 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 
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 
1858 void 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 
1910 static 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 
1940 static 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 
1985 static 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 
2055 static 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 
2081  if (ctrl->method == METHOD_SPRING_MAXENT){
2083  assert(D0);
2085  } else {
2087  }
2088  } else {
2089  if (ctrl->method == METHOD_SPRING_MAXENT){
2090  assert(D0);
2092  }
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 */
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 
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 &&
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);
2235  Multilevel_delete(grid0);
2236 }
2237 
2238 #ifdef GVIEWER
2239 struct multilevel_spring_electrical_embedding_data {
2240  int dim;
2241  SparseMatrix A;
2242  SparseMatrix D;
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 
2252 void 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 }
2259 void 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
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
SparseMatrix P
Definition: Multilevel.h:27
#define Multilevel_is_finest(grid)
Definition: Multilevel.h:67
#define MAX(a, b)
Definition: agerror.c:17
#define ABS(a)
Definition: arith.h:45
spring_electrical_control spring_electrical_control_new()
void QuadTree_get_repulsive_force(QuadTree qt, real *force, real *x, real bh, real p, real KP, real *counts, int *flag)
Definition: QuadTree.c:317
double xmax
Definition: geometry.c:20
real drand()
Definition: general.c:52
void spring_electrical_control_delete(spring_electrical_control ctrl)
SparseMatrix SparseMatrix_symmetrize_nodiag(SparseMatrix A, int pattern_symmetric_only)
Definition: SparseMatrix.c:163
#define MIN(a, b)
Definition: arith.h:35
oned_optimizer oned_optimizer_new(int i)
Multilevel Multilevel_get_coarsest(Multilevel grid)
Definition: Multilevel.c:1312
#define C
Definition: pack.c:29
#define FREE
Definition: PriorityQueue.c:23
#define assert(x)
Definition: cghdr.h:47
#define node_degree(i)
void spring_electrical_embedding(int dim, SparseMatrix A0, spring_electrical_control ctrl, real *node_weights, real *x, int *flag)
double xmin
Definition: geometry.c:20
SparseMatrix SparseMatrix_remove_diagonal(SparseMatrix A)
double ymax
Definition: geometry.c:20
real distance(real *x, int dim, int i, int j)
void interpolate_coord(int dim, SparseMatrix A, real *x)
void check_int_array_size(int **a, int len, int *lenmax)
void SparseMatrix_multiply_dense(SparseMatrix A, int ATransposed, real *v, int vTransposed, real **res, int res_transposed, int dim)
double ymin
Definition: geometry.c:20
int SparseMatrix_has_diagonal(SparseMatrix A)
int SparseMatrix_is_symmetric(SparseMatrix A, int test_pattern_symmetry_only)
Definition: SparseMatrix.c:178
Multilevel Multilevel_new(SparseMatrix A0, SparseMatrix D0, real *node_weights, Multilevel_control ctrl)
Definition: Multilevel.c:1294
#define Multilevel_is_coarsest(grid)
Definition: Multilevel.h:68
void spring_electrical_embedding_fast(int dim, SparseMatrix A0, spring_electrical_control ctrl, real *node_weights, real *x, int *flag)
void Multilevel_control_delete(Multilevel_control ctrl)
Definition: Multilevel.c:44
void oned_optimizer_delete(oned_optimizer opt)
void spring_electrical_embedding_slow(int dim, SparseMatrix A0, spring_electrical_control ctrl, real *node_weights, real *x, int *flag)
void export_embedding(FILE *fp, int dim, SparseMatrix A, real *x, real *width)
#define PI
#define AUTOP
void print_matrix(real *x, int n, int dim)
void spring_electrical_spring_embedding(int dim, SparseMatrix A0, SparseMatrix D, spring_electrical_control ctrl, real *node_weights, real *x, int *flag)
void print_padding(int n)
Definition: Multilevel.c:1249
#define max(x, y)
Definition: stress.c:794
void force_print(FILE *fp, int n, int dim, real *x, real *force)
void spring_electrical_control_print(spring_electrical_control ctrl)
SparseMatrix SparseMatrix_get_real_adjacency_matrix_symmetrized(SparseMatrix A)
real average_edge_length(SparseMatrix A, int dim, real *coord)
real get_angle(real *x, int dim, int i, int j)
Multilevel_control Multilevel_control_new(int scheme, int mode)
Definition: Multilevel.c:22
#define MINDIST
Definition: circular.c:20
#define MALLOC
Definition: PriorityQueue.c:21
void QuadTree_get_supernodes(QuadTree qt, real bh, real *point, int nodeid, int *nsuper, int *nsupermax, real **center, real **supernode_wgts, real **distances, double *counts, int *flag)
Definition: QuadTree.c:126
void multilevel_spring_electrical_embedding(int dim, SparseMatrix A, SparseMatrix D, spring_electrical_control ctrl, real *node_weights, real *label_sizes, real *x, int n_edge_label_nodes, int *edge_label_nodes, int *flag)
SparseMatrix SparseMatrix_symmetrize(SparseMatrix A, int pattern_symmetric_only)
Definition: SparseMatrix.c:151
void pcp_rotate(int n, int dim, real *x)
void SparseMatrix_delete(SparseMatrix A)
Definition: SparseMatrix.c:411
Definition: grammar.c:79
#define REALLOC
Definition: PriorityQueue.c:22
#define MEMCPY
Definition: PriorityQueue.c:24
if(aagss+aagstacksize-1<=aagssp)
Definition: grammar.c:1332
void post_process_smoothing(int dim, SparseMatrix A, spring_electrical_control ctrl, real *node_weights, real *x, int *flag)
void check_real_array_size(real **a, int len, int *lenmax)
int oned_optimizer_get(oned_optimizer opt)
void QuadTree_delete(QuadTree q)
Definition: QuadTree.c:421
SparseMatrix A
Definition: Multilevel.h:24
#define NULL
Definition: logic.h:39
SparseMatrix SparseMatrix_from_coordinate_arrays(int nz, int m, int n, int *irn, int *jcn, void *val0, int type, size_t sz)
Definition: SparseMatrix.c:957
int power_law_graph(SparseMatrix A)
QuadTree QuadTree_new_from_point_list(int dim, int n, int max_level, real *coord, real *weight)
Definition: QuadTree.c:347
EXTERN unsigned char Verbose
Definition: globals.h:64
void oned_optimizer_train(oned_optimizer opt, real work)
SparseMatrix R
Definition: Multilevel.h:28
void spring_maxent_embedding(int dim, SparseMatrix A0, SparseMatrix D, spring_electrical_control ctrl, real *node_weights, real *x, real rho, int *flag)
#define alpha
Definition: shapes.c:3902
Multilevel prev
Definition: Multilevel.h:31
real * node_weights
Definition: Multilevel.h:29
SparseMatrix D
Definition: Multilevel.h:25
double dist(Site *s, Site *t)
Definition: site.c:41
real distance_cropped(real *x, int dim, int i, int j)
pointf coord(node_t *n)
Definition: utils.c:202
Multilevel next
Definition: Multilevel.h:30
#define FALSE
Definition: cgraph.h:35
int comp_real(const void *x, const void *y)
Definition: legal.c:60
void Multilevel_delete(Multilevel grid)
Definition: Multilevel.c:66
void remove_overlap(int dim, SparseMatrix A, int m, real *x, real *label_sizes, int ntry, real initial_scaling, int do_shrinking, int *flag)
Definition: overlap.c:686
#define TRUE
Definition: cgraph.h:38
#define real
Definition: general.h:34