33static const double C = 0.2;
40static const double bh = 0.6;
44static const double tol = 0.001;
46static const double cool = 0.90;
72 "NONE",
"STRESS_MAJORIZATION_GRAPH_DIST",
"STRESS_MAJORIZATION_AVG_DIST",
"STRESS_MAJORIZATION_POWER_DIST",
"SPRING",
"TRIANGLE",
"RNG"
76 "NONE",
"NORMAL",
"FAST",
"HYBRID"
80 fprintf (stderr,
"spring_electrical_control:\n");
81 fprintf (stderr,
" repulsive exponent: %.03f\n", ctrl.
p);
82 fprintf(stderr,
" random start %d seed %d\n", (
int)ctrl.
random_start,
84 fprintf (stderr,
" K : %.03f C : %.03f\n", ctrl.
K,
C);
85 fprintf (stderr,
" max levels %d\n", ctrl.
multilevels);
87 fprintf (stderr,
" Barnes-Hutt constant %.03f tolerance %.03f maxiter %d\n",
bh,
tol, ctrl.
maxiter);
88 fprintf(stderr,
" cooling %.03f step size %.03f adaptive %d\n",
cool,
90 fprintf (stderr,
" beautify_leaves %d node weights %d rotation %.03f\n",
92 fprintf (stderr,
" smoothing %s overlap %d initial_scaling %.03f do_shrinking %d\n",
136 if (opt->
work[i] < opt->
work[i+1] && opt->
i > 0){
137 opt->
i =
MAX(0, opt->
i-1);
152 int *ia =
A->ia, *ja =
A->ja, i, j, k;
155 if (ia[
A->m] == 0)
return 1;
156 for (i = 0; i <
A->m; i++){
157 for (j = ia[i]; j < ia[i+1]; j++){
159 for (k = 0; k <
dim; k++){
165 return dist/ia[
A->m];
168static double update_step(
bool adaptive_cooling,
double step,
double Fnorm,
171 if (!adaptive_cooling) {
174 if (Fnorm >= Fnorm0){
176 }
else if (Fnorm > 0.95*Fnorm0){
179 step = 0.99*step/
cool;
185#define node_degree(i) (ia[(i)+1] - ia[(i)])
195 int m =
A->m, i, j, *ia =
A->ia, *ja =
A->ja;
204 for (i = 0; i < m; i++){
205 if (ia[i+1] - ia[i] != 1)
continue;
212 for (j = ia[p]; j < ia[p+1]; j++){
216 ints_append(&leaves, ja[j]);
219 assert(!ints_is_empty(&leaves));
220 dist /= (double)ints_size(&leaves);
222 double ang2 = 2 *
M_PI;
223 const double pad = 0.1;
227 assert(ang2 >= ang1);
229 if (ints_size(&leaves) > 1) step = (ang2 - ang1) / (
double)ints_size(&leaves);
230 for (
size_t k = 0; k < ints_size(&leaves); k++) {
243 double *x,
int *flag) {
248 double p = ctrl->
p, K = ctrl->
K, CRK, maxiter = ctrl->
maxiter, step = ctrl->
step, KP;
250 double *f =
NULL,
dist,
F, Fnorm = 0, Fnorm0;
253 double counts[4], *force =
NULL;
255 clock_t start, end, start0;
256 double qtree_cpu = 0, qtree_new_cpu = 0;
257 double total_cpu = 0;
262 if (!
A || maxiter <= 0)
return;
265 if (n <= 0 ||
dim <= 0)
return;
281 for (i = 0; i <
dim*n; i++) x[i] =
drand();
286 if (p >= 0) ctrl->
p = p = -1;
288 CRK = pow(
C, (2.-p)/3.)/K;
305 qtree_new_cpu += (double)(clock() - start) / CLOCKS_PER_SEC;
317 qtree_cpu += (double)(end - start) / CLOCKS_PER_SEC;
321 for (i = 0; i < n; i++){
323 for (j = ia[i]; j < ia[i+1]; j++){
324 if (ja[j] == i)
continue;
326 for (k = 0; k <
dim; k++){
327 f[k] -= CRK*(x[i*
dim+k] - x[ja[j]*
dim+k])*
dist;
334 for (i = 0; i < n; i++){
337 for (k = 0; k <
dim; k++)
F += f[k]*f[k];
340 if (
F > 0)
for (k = 0; k <
dim; k++) f[k] /=
F;
341 for (k = 0; k <
dim; k++) x[i*
dim+k] += step*f[k];
353 qtree_new_cpu += (double)(end - start) / CLOCKS_PER_SEC;
357 counts[0] + 0.85 * counts[1] + 3.3 * counts[2]);
360 fprintf(stderr,
"\r iter = %d, step = %f Fnorm = %f nz = %d K = %f ",iter, step, Fnorm,
A->nz,K);
364 step =
update_step(adaptive_cooling, step, Fnorm, Fnorm0);
365 }
while (step >
tol && iter < maxiter);
369 fprintf(stderr,
"\n iter = %d, step = %f Fnorm = %f nz = %d K = %f ",iter, step, Fnorm,
A->nz, K);
376 total_cpu += (double)(clock() - start0) / CLOCKS_PER_SEC;
377 if (
Verbose) fprintf(stderr,
"\n time for qtree = %f, qtree_force = %f, total cpu = %f\n",qtree_new_cpu, qtree_cpu, total_cpu);
390 double *x,
int *flag) {
396 double p = ctrl->
p, K = ctrl->
K, CRK, maxiter = ctrl->
maxiter, step = ctrl->
step, KP;
398 double *f =
NULL,
dist,
F, Fnorm = 0, Fnorm0;
403 clock_t start, end, start0, start2;
404 double total_cpu = 0;
408 fprintf(stderr,
"spring_electrical_embedding_slow");
409 if (!
A || maxiter <= 0)
return;
412 if (n <= 0 ||
dim <= 0)
return;
427 for (i = 0; i <
dim*n; i++) x[i] =
drand();
432 if (p >= 0) ctrl->
p = p = -1;
434 CRK = pow(
C, (2.-p)/3.)/K;
438 for (i = 0; i <
dim*n; i++) force[i] = 0;
449 for (i = 0; i < n; i++){
450 for (k = 0; k <
dim; k++) f[k] = 0.;
452 for (j = 0; j < n; j++){
453 if (j == i)
continue;
455 for (k = 0; k <
dim; k++){
456 f[k] += KP*(x[i*
dim+k] - x[j*
dim+k])/pow(
dist, 1.- p);
459 for (k = 0; k <
dim; k++) force[i*
dim+k] += f[k];
464 for (i = 0; i < n; i++){
465 for (k = 0; k <
dim; k++) f[k] = 0.;
467 for (j = ia[i]; j < ia[i+1]; j++){
468 if (ja[j] == i)
continue;
470 for (k = 0; k <
dim; k++){
471 f[k] -= CRK*(x[i*
dim+k] - x[ja[j]*
dim+k])*
dist;
474 for (k = 0; k <
dim; k++) force[i*
dim+k] += f[k];
479 for (i = 0; i < n; i++){
481 for (k = 0; k <
dim; k++) f[k] = force[i*
dim+k];
484 for (k = 0; k <
dim; k++)
F += f[k]*f[k];
488 if (
F > 0)
for (k = 0; k <
dim; k++) f[k] /=
F;
490 for (k = 0; k <
dim; k++) x[i*
dim+k] += step*f[k];
494 step =
update_step(adaptive_cooling, step, Fnorm, Fnorm0);
495 }
while (step >
tol && iter < maxiter);
499 fprintf(stderr,
"iter = %d, step = %f Fnorm = %f nsuper = 0 nz = %d K = %f ",iter, step, Fnorm,
A->nz,K);
506 total_cpu += (double)(clock() - start0) / CLOCKS_PER_SEC;
507 if (
Verbose) fprintf(stderr,
"time for supernode = 0, total cpu = %f\n", total_cpu);
523 double p = ctrl->
p, K = ctrl->
K, CRK, maxiter = ctrl->
maxiter, step = ctrl->
step, KP;
525 double *f =
NULL,
dist,
F, Fnorm = 0, Fnorm0;
529 int nsuper = 0, nsupermax = 10;
530 double *
center =
NULL, *supernode_wgts =
NULL, *distances =
NULL, nsuper_avg, counts = 0, counts_avg = 0;
532 clock_t start, end, start0, start2;
533 double qtree_cpu = 0;
534 double total_cpu = 0;
540 if (!
A || maxiter <= 0)
return;
543 if (n <= 0 ||
dim <= 0)
return;
549 supernode_wgts =
gv_calloc(nsupermax,
sizeof(
double));
550 distances =
gv_calloc(nsupermax,
sizeof(
double));
564 for (i = 0; i <
dim*n; i++) x[i] =
drand();
569 if (p >= 0) ctrl->
p = p = -1;
571 CRK = pow(
C, (2.-p)/3.)/K;
594 for (i = 0; i < n; i++){
595 for (k = 0; k <
dim; k++) f[k] = 0.;
597 for (j = ia[i]; j < ia[i+1]; j++){
598 if (ja[j] == i)
continue;
600 for (k = 0; k <
dim; k++){
601 f[k] -= CRK*(x[i*
dim+k] - x[ja[j]*
dim+k])*
dist;
611 &
center, &supernode_wgts, &distances, &counts);
615 qtree_cpu += (double)(end - start) / CLOCKS_PER_SEC;
617 counts_avg += counts;
618 nsuper_avg += nsuper;
619 for (j = 0; j < nsuper; j++){
621 for (k = 0; k <
dim; k++){
622 f[k] += supernode_wgts[j]*KP*(x[i*
dim+k] -
center[j*
dim+k])/pow(
dist, 1.- p);
626 for (j = 0; j < n; j++){
627 if (j == i)
continue;
629 for (k = 0; k <
dim; k++){
630 f[k] += KP*(x[i*
dim+k] - x[j*
dim+k])/pow(
dist, 1.- p);
637 for (k = 0; k <
dim; k++)
F += f[k]*f[k];
641 if (
F > 0)
for (k = 0; k <
dim; k++) f[k] /=
F;
643 for (k = 0; k <
dim; k++) x[i*
dim+k] += step*f[k];
654 step =
update_step(adaptive_cooling, step, Fnorm, Fnorm0);
655 }
while (step >
tol && iter < maxiter);
660 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);
662 fprintf(stderr,
"iter = %d, step = %f Fnorm = %f nsuper = %d nz = %d K = %f ",iter, step, Fnorm, (
int) nsuper_avg,
A->nz,K);
670 total_cpu += (double)(clock() - start0) / CLOCKS_PER_SEC;
671 if (
Verbose) fprintf(stderr,
"time for supernode = %f, total cpu = %f\n",qtree_cpu, total_cpu);
681 free(supernode_wgts);
687 double *x,
int *flag) {
694 double p = ctrl->
p, K = ctrl->
K, CRK, maxiter = ctrl->
maxiter, step = ctrl->
step, KP;
699 double *f =
NULL,
dist,
F, Fnorm = 0, Fnorm0;
703 int nsuper = 0, nsupermax = 10;
705 int max_qtree_level = 10;
707 if (!
A || maxiter <= 0)
return;
709 if (n <= 0 ||
dim <= 0)
return;
714 supernode_wgts =
gv_calloc(nsupermax,
sizeof(
double));
715 distances =
gv_calloc(nsupermax,
sizeof(
double));
732 for (i = 0; i <
dim*n; i++) x[i] =
drand();
737 if (p >= 0) ctrl->
p = p = -1;
739 CRK = pow(
C, (2.-p)/3.)/K;
745 memcpy(xold, x,
sizeof(
double)*
dim*n);
754 for (i = 0; i < n; i++){
755 for (k = 0; k <
dim; k++) f[k] = 0.;
758 for (j = ia[i]; j < ia[i+1]; j++){
759 if (ja[j] == i)
continue;
761 for (k = 0; k <
dim; k++){
762 f[k] -= CRK*(x[i*
dim+k] - x[ja[j]*
dim+k])*
dist;
766 for (j =
id[i]; j <
id[i+1]; j++){
767 if (jd[j] == i)
continue;
769 for (k = 0; k <
dim; k++){
782 &
center, &supernode_wgts, &distances, &counts);
783 for (j = 0; j < nsuper; j++){
785 for (k = 0; k <
dim; k++){
786 f[k] += supernode_wgts[j]*KP*(x[i*
dim+k] -
center[j*
dim+k])/pow(
dist, 1.- p);
790 for (j = 0; j < n; j++){
791 if (j == i)
continue;
793 for (k = 0; k <
dim; k++){
794 f[k] += KP*(x[i*
dim+k] - x[j*
dim+k])/pow(
dist, 1.- p);
801 for (k = 0; k <
dim; k++)
F += f[k]*f[k];
805 if (
F > 0)
for (k = 0; k <
dim; k++) f[k] /=
F;
807 for (k = 0; k <
dim; k++) x[i*
dim+k] += step*f[k];
813 step =
update_step(adaptive_cooling, step, Fnorm, Fnorm0);
814 }
while (step >
tol && iter < maxiter);
823 free(supernode_wgts);
828 int i, j, k, *ia =
A->ia, *ja =
A->ja, nz;
829 double alpha = 0.5, beta;
832 for (i = 0; i <
A->m; i++){
833 for (k = 0; k <
dim; k++) y[k] = 0;
835 for (j = ia[i]; j < ia[i+1]; j++){
836 if (ja[j] == i)
continue;
838 for (k = 0; k <
dim; k++){
839 y[k] += x[ja[j]*
dim + k];
844 for (k = 0; k <
dim; k++) x[i*
dim+k] =
alpha*x[i*
dim+k] + beta*y[k];
851 int nc, *ia, *ja, i, j, k;
858 for (i = 0; i < nc; i++){
859 for (j = ia[i]+1; j < ia[i+1]; j++){
860 for (k = 0; k <
dim; k++){
868 int m, max = 0, i, *ia =
A->ia, *ja =
A->ja, j, deg;
871 int *mask =
gv_calloc(m + 1,
sizeof(
int));
873 for (i = 0; i < m + 1; i++){
877 for (i = 0; i < m; i++){
879 for (j = ia[i]; j < ia[i+1]; j++){
880 if (i == ja[j])
continue;
884 max =
MAX(max, mask[deg]);
886 if (mask[1] > 0.8*max && mask[1] > 0.3*m) res =
true;
893 double y[4], axis[2],
center[2],
dist, x0, x1;
896 for (i = 0; i <
dim*
dim; i++) y[i] = 0;
898 for (i = 0; i < n; i++){
899 for (k = 0; k <
dim; k++){
903 for (i = 0; i <
dim; i++)
center[i] /= n;
904 for (i = 0; i < n; i++){
905 for (k = 0; k <
dim; k++){
910 for (i = 0; i < n; i++){
911 for (k = 0; k <
dim; k++){
912 for (l = 0; l <
dim; l++){
918 axis[0] = 0; axis[1] = 1;
926 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]);
929 dist = sqrt(1+axis[0]*axis[0]);
930 axis[0] = axis[0]/
dist;
931 axis[1] = axis[1]/
dist;
932 for (i = 0; i < n; i++){
933 x0 = x[
dim*i]*axis[0]+x[
dim*i+1]*axis[1];
934 x1 = -x[
dim*i]*axis[1]+x[
dim*i+1]*axis[0];
943static void rotate(
int n,
int dim,
double *x,
double angle){
945 double axis[2],
center[2], x0, x1;
946 double radian = 3.14159/180;
950 for (i = 0; i < n; i++){
951 for (k = 0; k <
dim; k++){
955 for (i = 0; i <
dim; i++)
center[i] /= n;
956 for (i = 0; i < n; i++){
957 for (k = 0; k <
dim; k++){
961 axis[0] = cos(-angle*radian);
962 axis[1] = sin(-angle*radian);
963 for (i = 0; i < n; i++){
964 x0 = x[
dim*i]*axis[0]+x[
dim*i+1]*axis[1];
965 x1 = -x[
dim*i]*axis[1]+x[
dim*i+1]*axis[0];
980 for (i = 0; i <
A->m; i++) mask[i] = 1;
981 for (i = 0; i < n_edge_label_nodes; i++) {
982 if (edge_label_nodes[i] >= 0 && edge_label_nodes[i] <
A->m) mask[edge_label_nodes[i]] = -1;
985 for (i = 0; i <
A->m; i++) {
986 if (mask[i] >= 0) mask[i] = nnodes++;
990 for (i = 0; i <
A->m; i++){
992 for (k = 0; k <
dim; k++) x[i*
dim+k] = x2[mask[i]*
dim+k];
996 for (i = 0; i < n_edge_label_nodes; i++){
997 ii = edge_label_nodes[i];
998 len =
A->ia[ii+1] -
A->ia[ii];
1000 assert(mask[ii] < 0);
1001 for (k = 0; k <
dim; k++) {
1004 for (j =
A->ia[ii]; j <
A->ia[ii+1]; j++){
1005 for (k = 0; k <
dim; k++){
1006 x[ii *
dim + k] += x[
A->ja[j] *
dim + k];
1009 for (k = 0; k <
dim; k++) {
1018 int i,
id = 0, nz, j, jj, ii;
1019 int *ia =
A->
ia, *ja =
A->ja, *irn =
NULL, *jcn =
NULL;
1024 for (i = 0; i <
A->m; i++) mask[i] = 1;
1026 for (i = 0; i < n_edge_label_nodes; i++){
1027 mask[edge_label_nodes[i]] = -1;
1030 for (i = 0; i <
A->m; i++) {
1031 if (mask[i] > 0) mask[i] =
id++;
1035 for (i = 0; i <
A->m; i++){
1036 if (mask[i] < 0)
continue;
1037 for (j = ia[i]; j < ia[i+1]; j++){
1038 if (mask[ja[j]] >= 0) {
1043 for (jj = ia[ii]; jj < ia[ii+1]; jj++){
1044 if (ja[jj] != i && mask[ja[jj]] >= 0) nz++;
1055 for (i = 0; i <
A->m; i++){
1056 if (mask[i] < 0)
continue;
1057 for (j = ia[i]; j < ia[i+1]; j++){
1058 if (mask[ja[j]] >= 0) {
1060 jcn[nz++] = mask[ja[j]];
1064 for (jj = ia[ii]; jj < ia[ii+1]; jj++){
1065 if (ja[jj] != i && mask[ja[jj]] >= 0) {
1067 jcn[nz++] = mask[ja[jj]];
1084 double *label_sizes,
double *x,
1085 int n_edge_label_nodes,
1086 int *edge_label_nodes,
int *flag) {
1105 if (n <= 0 ||
dim <= 0)
return;
1115 && n_edge_label_nodes > 0){
1146 if (plg) ctrl->
p = -1.8;
1154 fprintf(stderr,
"coarsest level -- %d, n = %d\n",
grid->level,
grid->n);
1156 fprintf(stderr,
"level -- %d, n = %d\n",
grid->level,
grid->n);
1164 fprintf(stderr,
"QUAD_TREE_HYBRID, size larger than %d, switch to fast quadtree",
QUAD_TREE_HYBRID_SIZE);
1186 ctrl->
K = ctrl->
K * 0.75;
1193 fprintf(stderr,
"layout time %f\n", (
double)(clock() - cpu) / CLOCKS_PER_SEC);
1199 if (
Verbose) fprintf(stderr,
"ctrl->overlap=%d\n",ctrl->
overlap);
void print_padding(int n)
Multilevel Multilevel_get_coarsest(Multilevel grid)
void Multilevel_delete(Multilevel grid)
Multilevel Multilevel_new(SparseMatrix A0, const Multilevel_control ctrl)
#define Multilevel_is_coarsest(grid)
#define Multilevel_is_finest(grid)
void QuadTree_get_repulsive_force(QuadTree qt, double *force, double *x, double bh, double p, double KP, double *counts)
QuadTree QuadTree_new_from_point_list(int dim, int n, int max_level, double *coord)
void QuadTree_delete(QuadTree q)
void QuadTree_get_supernodes(QuadTree qt, double bh, double *pt, int nodeid, int *nsuper, int *nsupermax, double **center, double **supernode_wgts, double **distances, double *counts)
SparseMatrix SparseMatrix_symmetrize(SparseMatrix A, bool pattern_symmetric_only)
void SparseMatrix_multiply_dense(SparseMatrix A, const double *v, double *res, int dim)
bool SparseMatrix_is_symmetric(SparseMatrix A, bool test_pattern_symmetry_only)
SparseMatrix SparseMatrix_from_coordinate_arrays(int nz, int m, int n, int *irn, int *jcn, void *val0, int type, size_t sz)
void SparseMatrix_delete(SparseMatrix A)
SparseMatrix SparseMatrix_get_real_adjacency_matrix_symmetrized(SparseMatrix A)
bool SparseMatrix_has_diagonal(SparseMatrix A)
SparseMatrix SparseMatrix_remove_diagonal(SparseMatrix A)
Memory allocation wrappers that exit on failure.
static void * gv_calloc(size_t nmemb, size_t size)
API for compacted arrays of booleans.
static bitarray_t bitarray_new(size_t size_bits)
create an array of the given element length
static bool bitarray_get(bitarray_t self, size_t index)
get the value of the given element
static void bitarray_set(bitarray_t *self, size_t index, bool value)
set or clear the value of the given element
static void bitarray_reset(bitarray_t *self)
free underlying resources and leave a bit array empty
static double dist(int dim, double *x, double *y)
double distance(double *x, int dim, int i, int j)
double distance_cropped(double *x, int dim, int i, int j)
static double len(glCompPoint p)
#define DEFINE_LIST(name, type)
void remove_overlap(int dim, SparseMatrix A, double *x, double *label_sizes, int ntry, double initial_scaling, int edge_labeling_scheme, int n_constr_nodes, int *constr_nodes, SparseMatrix A_constr, bool do_shrinking)
@ ELSCHEME_STRAIGHTLINE_PENALTY2
@ ELSCHEME_STRAIGHTLINE_PENALTY
void post_process_smoothing(int dim, SparseMatrix A, spring_electrical_control ctrl, double *x)
static bool power_law_graph(SparseMatrix A)
void spring_electrical_embedding(int dim, SparseMatrix A0, spring_electrical_control *ctrl, double *x, int *flag)
void spring_electrical_embedding_fast(int dim, SparseMatrix A0, spring_electrical_control *ctrl, double *x, int *flag)
static char * smoothings[]
void spring_electrical_control_print(spring_electrical_control ctrl)
void multilevel_spring_electrical_embedding(int dim, SparseMatrix A0, spring_electrical_control *ctrl, double *label_sizes, double *x, int n_edge_label_nodes, int *edge_label_nodes, int *flag)
static double update_step(bool adaptive_cooling, double step, double Fnorm, double Fnorm0)
spring_electrical_control spring_electrical_control_new(void)
static const int quadtree_size
cut off size above which quadtree approximation is used
static SparseMatrix shorting_edge_label_nodes(SparseMatrix A, int n_edge_label_nodes, int *edge_label_nodes)
static int oned_optimizer_get(const oned_optimizer opt)
static void prolongate(int dim, SparseMatrix A, SparseMatrix P, SparseMatrix R, double *x, double *y, double delta)
static void beautify_leaves(int dim, SparseMatrix A, double *x)
static oned_optimizer oned_optimizer_new(int i)
static void attach_edge_label_coordinates(int dim, SparseMatrix A, int n_edge_label_nodes, int *edge_label_nodes, double *x, double *x2)
double average_edge_length(SparseMatrix A, int dim, double *coord)
static void oned_optimizer_train(oned_optimizer *opt, double work)
void pcp_rotate(int n, int dim, double *x)
static void rotate(int n, int dim, double *x, double angle)
void spring_electrical_spring_embedding(int dim, SparseMatrix A0, SparseMatrix D, spring_electrical_control *ctrl, double *x, int *flag)
static void set_leaves(double *x, int dim, double dist, double ang, int i, int j)
static void spring_electrical_embedding_slow(int dim, SparseMatrix A0, spring_electrical_control *ctrl, double *x, int *flag)
static void interpolate_coord(int dim, SparseMatrix A, double *x)
@ ERROR_NOT_SQUARE_MATRIX
bool random_start
whether to apply SE from a random layout, or from exisiting layout
static point center(point vertex[], size_t n)
graphs, nodes and edges info: Agraphinfo_t, Agnodeinfo_t and Agedgeinfo_t