From: erg Date: Mon, 28 Feb 2011 16:19:43 +0000 (+0000) Subject: Add new features to sfdp X-Git-Tag: LAST_LIBGRAPH~32^2~979 X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=44e27a0185c39a6e4f196c61c50d7cad8e9f15aa;p=graphviz Add new features to sfdp --- diff --git a/lib/sfdpgen/LinkedList.c b/lib/sfdpgen/LinkedList.c index b18888247..2d4c62198 100644 --- a/lib/sfdpgen/LinkedList.c +++ b/lib/sfdpgen/LinkedList.c @@ -28,6 +28,14 @@ SingleLinkedList SingleLinkedList_new(void *data){ return head; } +SingleLinkedList SingleLinkedList_new_int(int i){ + int *data; + data = malloc(sizeof(int)); + data[0] = i; + return SingleLinkedList_new((void*) data); +} + + void SingleLinkedList_delete(SingleLinkedList head, void (*linklist_deallocator)(void*)){ SingleLinkedList next; @@ -48,6 +56,13 @@ SingleLinkedList SingleLinkedList_prepend(SingleLinkedList l, void *data){ return head; } +SingleLinkedList SingleLinkedList_prepend_int(SingleLinkedList l, int i){ + int *data; + data = malloc(sizeof(int)); + data[0] = i; + return SingleLinkedList_prepend(l, (void*) data); +} + void* SingleLinkedList_get_data(SingleLinkedList l){ return l->data; } diff --git a/lib/sfdpgen/LinkedList.h b/lib/sfdpgen/LinkedList.h index 52abcf929..b36683f25 100644 --- a/lib/sfdpgen/LinkedList.h +++ b/lib/sfdpgen/LinkedList.h @@ -22,8 +22,10 @@ struct SingleLinkedList_struct { }; SingleLinkedList SingleLinkedList_new(void *data); +SingleLinkedList SingleLinkedList_new_int(int i); void SingleLinkedList_delete(SingleLinkedList head, void (*linklist_deallocator)(void*)); SingleLinkedList SingleLinkedList_prepend(SingleLinkedList l, void *data); +SingleLinkedList SingleLinkedList_prepend_int(SingleLinkedList l, int i); void* SingleLinkedList_get_data(SingleLinkedList l); diff --git a/lib/sfdpgen/Makefile.am b/lib/sfdpgen/Makefile.am index 4d8e0f277..6b8a36cbe 100644 --- a/lib/sfdpgen/Makefile.am +++ b/lib/sfdpgen/Makefile.am @@ -21,6 +21,7 @@ AM_CPPFLAGS = \ noinst_HEADERS = sfdpinternal.h spring_electrical.h \ LinkedList.h sparse_solve.h post_process.h \ + stress_model.h uniform_stress.h \ QuadTree.h Multilevel.h sfdp.h PriorityQueue.h if WITH_SFDP @@ -29,6 +30,7 @@ endif libsfdpgen_C_la_SOURCES = sfdpinit.c spring_electrical.c \ LinkedList.c sparse_solve.c post_process.c \ + stress_model.c uniform_stress.c \ QuadTree.c Multilevel.c PriorityQueue.c EXTRA_DIST = Makefile.old sfdp.vcproj diff --git a/lib/sfdpgen/Multilevel.c b/lib/sfdpgen/Multilevel.c index 28041efc8..c30d52895 100644 --- a/lib/sfdpgen/Multilevel.c +++ b/lib/sfdpgen/Multilevel.c @@ -18,40 +18,8 @@ #include "assert.h" #include "arith.h" -#define MALLOC gmalloc -#define REALLOC grealloc -#define FREE free -#define MEMCPY memcpy - -static int irand(int n){ - /* 0, 1, ..., n-1 */ - assert(n > 1); - /*return (int) MIN(floor(drand()*n),n-1);*/ - return rand()%n; -} -static int *random_permutation(int n) -{ - int *p; - int i, j, pp, len; - if (n <= 0) - return NULL; - p = N_GNEW(n, int); - for (i = 0; i < n; i++) - p[i] = i; - - len = n; - while (len > 1) { - j = irand(len); - pp = p[len - 1]; - p[len - 1] = p[j]; - p[j] = pp; - len--; - } - return p; -} - -Multilevel_control Multilevel_control_new(){ +Multilevel_control Multilevel_control_new(int scheme, int mode){ Multilevel_control ctrl; ctrl = N_GNEW(1,struct Multilevel_control_struct); @@ -59,9 +27,17 @@ Multilevel_control Multilevel_control_new(){ ctrl->min_coarsen_factor = 0.75; ctrl->maxlevel = 1<<30; ctrl->randomize = TRUE; - ctrl->coarsen_scheme = COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_CLUSTER_PERNODE_LEAVES_FIRST; - ctrl->coarsen_scheme = COARSEN_HYBRID; - /* ctrl->coarsen_scheme = COARSEN_INDEPENDENT_VERTEX_SET_RS;*/ + /* now set in spring_electrical_control_new(), as well as by command line argument -c + ctrl->coarsen_scheme = COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_CLUSTER_PERNODE_LEAVES_FIRST; + ctrl->coarsen_scheme = COARSEN_INDEPENDENT_VERTEX_SET_RS; + ctrl->coarsen_scheme = COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_EDGE_PERNODE; + ctrl->coarsen_scheme = COARSEN_HYBRID; + ctrl->coarsen_scheme = COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_EDGE_PERNODE_SUPERNODES_FIRST; + ctrl->coarsen_mode = COARSEN_MODE_FORCEFUL; or COARSEN_MODE_GENTLE; + */ + + ctrl->coarsen_scheme = scheme; + ctrl->coarsen_mode = mode; return ctrl; } @@ -69,7 +45,7 @@ void Multilevel_control_delete(Multilevel_control ctrl){ FREE(ctrl); } -static Multilevel Multilevel_init(SparseMatrix A, real *node_weights){ +static Multilevel Multilevel_init(SparseMatrix A, SparseMatrix D, real *node_weights){ Multilevel grid; if (!A) return NULL; assert(A->m == A->n); @@ -77,6 +53,7 @@ static Multilevel Multilevel_init(SparseMatrix A, real *node_weights){ grid->level = 0; grid->n = A->n; grid->A = A; + grid->D = D; grid->P = NULL; grid->R = NULL; grid->node_weights = node_weights; @@ -90,9 +67,13 @@ void Multilevel_delete(Multilevel grid){ if (!grid) return; if (grid->A){ if (grid->level == 0) { - if (grid->delete_top_level_A) SparseMatrix_delete(grid->A); + if (grid->delete_top_level_A) { + SparseMatrix_delete(grid->A); + if (grid->D) SparseMatrix_delete(grid->D); + } } else { SparseMatrix_delete(grid->A); + if (grid->D) SparseMatrix_delete(grid->D); } } SparseMatrix_delete(grid->P); @@ -831,9 +812,143 @@ static void maximal_independent_edge_set_heavest_edge_pernode_scaled(SparseMatri } } -static void Multilevel_coarsen(SparseMatrix A, SparseMatrix *cA, real *node_wgt, real **cnode_wgt, - SparseMatrix *P, SparseMatrix *R, Multilevel_control ctrl, int *coarsen_scheme_used){ - int *matching = NULL, nmatch, nc, nzc, n, i; +SparseMatrix DistanceMatrix_restrict_cluster(int ncluster, int *clusterp, int *cluster, SparseMatrix P, SparseMatrix R, SparseMatrix D){ +#if 0 + /* this construct a distance matrix of a coarse graph, for a coarsen give by merging all nodes in eahc cluster */ + SparseMatrix cD = NULL; + int i, j, nzc; + int **irn, **jcn; + real **val; + int n = D->m; + int *assignment = NULL; + int nz; + int *id = D->ia, jd = D->ja; + int *mask = NULL; + int *nnodes, *mask; + real *d = NULL; + + + assert(D->m == D->n); + if (!D) return NULL; + if (D->a && D->type == MATRIX_TYPE_REAL) d = (real*) D->val; + + irn = N_GNEW(ncluster,int*); + jcn = N_GNEW(ncluster,int*); + val = N_GNEW(ncluster,real*); + assignment = N_GNEW(n,int); + nz = N_GNEW(ncluster,int); + mask = N_GNEW(n,int); + nnodes = N_GNEW(ncluster,int); + + + /* find ncluster-subgrahs induced by the ncluster -clusters, find the diameter of each, + then use the radius as the distance from the supernode to the rest of the "world" + */ + for (i = 0; i < ncluster; i++) nz[i] = 0; + for (i = 0; i < ncluster; i++){ + for (j = clusterp[i]; j < clusterp[i+1]; j++){ + assert(clusterp[i+1] > clusterp[i]); + assignment[cluster[j]] = i; + } + } + + for (i = 0; i < n; i++){/* figure out how many entries per submatrix */ + ic = asignment[i]; + for (j = id[i]; j < id[i+1]; j++){ + if (i != jd[j] && ic == assignment[jd[j]]) { + nz[ic]++; + } + } + } + for (i = 0; i < ncluster; i++) { + irn[i] = N_GNEW(nz[i],int); + jcn[i] = N_GNEW(nz[i],int); + val[i] = N_GNEW(nz[i],int); + val[i] = NULL; + } + + + for (i = 0; i < ncluster; i++) nz[i] = 0;/* get subgraphs */ + for (i = 0; i < n; i++) mask[i] = -1; + for (i = 0; i < ncluster; i++) nnodes[i] = -1; + for (i = 0; i < n; i++){ + ic = asignment[i]; + ii = mask[i]; + if (ii < 0){ + mask[i] = ii = nnodes[ic]; + nnodes[ic]++; + } + for (j = id[i]; j < id[i+1]; j++){ + jc = assignment[jd[j]]; + if (i != jd[j] && ic == jc) { + jj = mask[jd[j]]; + if (jj < 0){ + mask[jd[j]] = jj = nnodes[jc]; + nnodes[jc]++; + } + irn[ic][nz[ic]] = ii; + jcn[ic][nz[ic]] = jj; + if (d) val[ic][nz[ic]] = d[j]; + } + } + } + + for (i = 0; i < ncluster; i++){/* form subgraphs */ + SparseMatrix A; + A = SparseMatrix_from_coordinate_arrays(nz[nz[i]], nnodes[i], nnodes[i], irn[i], jcn[i], (void*) val[i], MATRIX_TYPE_REAL); + + SparseMatrix_delete(A); + } + + + for (i = 0; i < ncluster; i++){ + for (j = clusterp[i]; j < clusterp[i+1]; j++){ + assert(clusterp[i+1] > clusterp[i]); + irn[nzc] = cluster[j]; + jcn[nzc] = i; + val[nzc++] = 1.; + } + } + assert(nzc == n); + cD = SparseMatrix_multiply3(R, D, P); + + SparseMatrix_set_symmetric(cD); + SparseMatrix_set_pattern_symmetric(cD); + cD = SparseMatrix_remove_diagonal(cD); + + FREE(nz); + FREE(assignment); + for (i = 0; i < ncluster; i++){ + FREE(irn[i]); + FREE(jcn[i]); + FREE(val[i]); + } + FREE(irn); FREE(jcn); FREE(val); + FREE(mask); + FREE(nnodes); + + return cD; +#endif + return NULL; +} + +SparseMatrix DistanceMatrix_restrict_matching(int *matching, SparseMatrix D){ + if (!D) return NULL; + assert(0);/* not yet implemented! */ + return NULL; +} + +SparseMatrix DistanceMatrix_restrict_filtering(int *mask, int is_C, int is_F, SparseMatrix D){ + /* max independent vtx set based coarsening. Coarsen nodes has mask >= is_C. Fine nodes == is_F. */ + if (!D) return NULL; + assert(0);/* not yet implemented! */ + return NULL; +} + +static void Multilevel_coarsen_internal(SparseMatrix A, SparseMatrix *cA, SparseMatrix D, SparseMatrix *cD, + real *node_wgt, real **cnode_wgt, + SparseMatrix *P, SparseMatrix *R, Multilevel_control ctrl, int *coarsen_scheme_used){ + int *matching = NULL, nmatch = 0, nc, nzc, n, i; int *irn = NULL, *jcn = NULL, *ia = NULL, *ja = NULL; real *val = NULL; SparseMatrix B = NULL; @@ -842,6 +957,7 @@ static void Multilevel_coarsen(SparseMatrix A, SparseMatrix *cA, real *node_wgt, assert(A->m == A->n); *cA = NULL; + *cD = NULL; *P = NULL; *R = NULL; n = A->m; @@ -855,7 +971,7 @@ static void Multilevel_coarsen(SparseMatrix A, SparseMatrix *cA, real *node_wgt, fprintf(stderr, "hybrid scheme, try COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_EDGE_PERNODE_LEAVES_FIRST first\n"); #endif *coarsen_scheme_used = ctrl->coarsen_scheme = COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_EDGE_PERNODE_LEAVES_FIRST; - Multilevel_coarsen(A, cA, node_wgt, cnode_wgt, P, R, ctrl, coarsen_scheme_used); + Multilevel_coarsen_internal(A, cA, D, cD, node_wgt, cnode_wgt, P, R, ctrl, coarsen_scheme_used); if (!(*cA)) { #ifdef DEBUG_PRINT @@ -863,7 +979,7 @@ static void Multilevel_coarsen(SparseMatrix A, SparseMatrix *cA, real *node_wgt, fprintf(stderr, "switching to COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_EDGE_PERNODE_SUPERNODES_FIRST\n"); #endif *coarsen_scheme_used = ctrl->coarsen_scheme = COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_EDGE_PERNODE_SUPERNODES_FIRST; - Multilevel_coarsen(A, cA, node_wgt, cnode_wgt, P, R, ctrl, coarsen_scheme_used); + Multilevel_coarsen_internal(A, cA, D, cD, node_wgt, cnode_wgt, P, R, ctrl, coarsen_scheme_used); } if (!(*cA)) { @@ -872,7 +988,7 @@ static void Multilevel_coarsen(SparseMatrix A, SparseMatrix *cA, real *node_wgt, fprintf(stderr, "switching to COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_CLUSTER_PERNODE_LEAVES_FIRST\n"); #endif *coarsen_scheme_used = ctrl->coarsen_scheme = COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_CLUSTER_PERNODE_LEAVES_FIRST; - Multilevel_coarsen(A, cA, node_wgt, cnode_wgt, P, R, ctrl, coarsen_scheme_used); + Multilevel_coarsen_internal(A, cA, D, cD, node_wgt, cnode_wgt, P, R, ctrl, coarsen_scheme_used); } if (!(*cA)) { @@ -881,7 +997,7 @@ static void Multilevel_coarsen(SparseMatrix A, SparseMatrix *cA, real *node_wgt, fprintf(stderr, "switching to COARSEN_INDEPENDENT_VERTEX_SET\n"); #endif *coarsen_scheme_used = ctrl->coarsen_scheme = COARSEN_INDEPENDENT_VERTEX_SET; - Multilevel_coarsen(A, cA, node_wgt, cnode_wgt, P, R, ctrl, coarsen_scheme_used); + Multilevel_coarsen_internal(A, cA, D, cD, node_wgt, cnode_wgt, P, R, ctrl, coarsen_scheme_used); } @@ -891,7 +1007,7 @@ static void Multilevel_coarsen(SparseMatrix A, SparseMatrix *cA, real *node_wgt, fprintf(stderr, "switching to COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_EDGE_PERNODE\n"); #endif *coarsen_scheme_used = ctrl->coarsen_scheme = COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_EDGE_PERNODE; - Multilevel_coarsen(A, cA, node_wgt, cnode_wgt, P, R, ctrl, coarsen_scheme_used); + Multilevel_coarsen_internal(A, cA, D, cD, node_wgt, cnode_wgt, P, R, ctrl, coarsen_scheme_used); } ctrl->coarsen_scheme = COARSEN_HYBRID; break; @@ -907,7 +1023,7 @@ static void Multilevel_coarsen(SparseMatrix A, SparseMatrix *cA, real *node_wgt, } assert(ncluster <= n); nc = ncluster; - if (nc > ctrl->min_coarsen_factor*n || nc < ctrl->minsize) { + if ((ctrl->coarsen_mode == COARSEN_MODE_GENTLE && nc > ctrl->min_coarsen_factor*n) || nc < ctrl->minsize) { #ifdef DEBUG_PRINT if (Verbose) fprintf(stderr, "nc = %d, nf = %d, minsz = %d, coarsen_factor = %f coarsening stops\n",nc, n, ctrl->minsize, ctrl->min_coarsen_factor); @@ -929,15 +1045,26 @@ static void Multilevel_coarsen(SparseMatrix A, SparseMatrix *cA, real *node_wgt, assert(nzc == n); *P = SparseMatrix_from_coordinate_arrays(nzc, n, nc, irn, jcn, (void *) val, MATRIX_TYPE_REAL); *R = SparseMatrix_transpose(*P); - B = SparseMatrix_multiply(*R, A); - if (!B) goto RETURN; - *cA = SparseMatrix_multiply(B, *P); + + *cD = DistanceMatrix_restrict_cluster(ncluster, clusterp, cluster, *P, *R, D); + + *cA = SparseMatrix_multiply3(*R, A, *P); + + /* + B = SparseMatrix_multiply(*R, A); + if (!B) goto RETURN; + *cA = SparseMatrix_multiply(B, *P); + */ if (!*cA) goto RETURN; + SparseMatrix_multiply_vector(*R, node_wgt, cnode_wgt, FALSE); *R = SparseMatrix_divide_row_by_degree(*R); SparseMatrix_set_symmetric(*cA); SparseMatrix_set_pattern_symmetric(*cA); *cA = SparseMatrix_remove_diagonal(*cA); + + + break; case COARSEN_INDEPENDENT_EDGE_SET: maximal_independent_edge_set(A, ctrl->randomize, &matching, &nmatch); @@ -948,7 +1075,7 @@ static void Multilevel_coarsen(SparseMatrix A, SparseMatrix *cA, real *node_wgt, if (ctrl->coarsen_scheme == COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_EDGE_PERNODE_DEGREE_SCALED) maximal_independent_edge_set_heavest_edge_pernode_scaled(A, ctrl->randomize, &matching, &nmatch); nc = nmatch; - if (nc > ctrl->min_coarsen_factor*n || nc < ctrl->minsize) { + if ((ctrl->coarsen_mode == COARSEN_MODE_GENTLE && nc > ctrl->min_coarsen_factor*n) || nc < ctrl->minsize) { #ifdef DEBUG_PRINT if (Verbose) fprintf(stderr, "nc = %d, nf = %d, minsz = %d, coarsen_factor = %f coarsening stops\n",nc, n, ctrl->minsize, ctrl->min_coarsen_factor); @@ -982,15 +1109,23 @@ static void Multilevel_coarsen(SparseMatrix A, SparseMatrix *cA, real *node_wgt, assert(nzc == n); *P = SparseMatrix_from_coordinate_arrays(nzc, n, nc, irn, jcn, (void *) val, MATRIX_TYPE_REAL); *R = SparseMatrix_transpose(*P); - B = SparseMatrix_multiply(*R, A); - if (!B) goto RETURN; - *cA = SparseMatrix_multiply(B, *P); + *cA = SparseMatrix_multiply3(*R, A, *P); + /* + B = SparseMatrix_multiply(*R, A); + if (!B) goto RETURN; + *cA = SparseMatrix_multiply(B, *P); + */ if (!*cA) goto RETURN; SparseMatrix_multiply_vector(*R, node_wgt, cnode_wgt, FALSE); *R = SparseMatrix_divide_row_by_degree(*R); SparseMatrix_set_symmetric(*cA); SparseMatrix_set_pattern_symmetric(*cA); *cA = SparseMatrix_remove_diagonal(*cA); + + + *cD = DistanceMatrix_restrict_matching(matching, D); + *cD=NULL; + break; case COARSEN_INDEPENDENT_VERTEX_SET: case COARSEN_INDEPENDENT_VERTEX_SET_RS: @@ -1002,7 +1137,7 @@ static void Multilevel_coarsen(SparseMatrix A, SparseMatrix *cA, real *node_wgt, ia = A->ia; ja = A->ja; nc = nvset; - if (nc > ctrl->min_coarsen_factor*n || nc < ctrl->minsize) { + if ((ctrl->coarsen_mode == COARSEN_MODE_GENTLE && nc > ctrl->min_coarsen_factor*n) || nc < ctrl->minsize) { #ifdef DEBUG_PRINT if (Verbose) fprintf(stderr, "nc = %d, nf = %d, minsz = %d, coarsen_factor = %f coarsening stops\n",nc, n, ctrl->minsize, ctrl->min_coarsen_factor); @@ -1039,14 +1174,14 @@ static void Multilevel_coarsen(SparseMatrix A, SparseMatrix *cA, real *node_wgt, *P = SparseMatrix_from_coordinate_arrays(nzc, n, nc, irn, jcn, (void *) val, MATRIX_TYPE_REAL); *R = SparseMatrix_transpose(*P); - B = SparseMatrix_multiply(*R, A); - if (!B) goto RETURN; - *cA = SparseMatrix_multiply(B, *P); + *cA = SparseMatrix_multiply3(*R, A, *P); if (!*cA) goto RETURN; SparseMatrix_multiply_vector(*R, node_wgt, cnode_wgt, FALSE); SparseMatrix_set_symmetric(*cA); SparseMatrix_set_pattern_symmetric(*cA); *cA = SparseMatrix_remove_diagonal(*cA); + + *cD = DistanceMatrix_restrict_filtering(vset, MAX_IND_VTX_SET_C, MAX_IND_VTX_SET_F, D); break; default: goto RETURN; @@ -1060,6 +1195,54 @@ static void Multilevel_coarsen(SparseMatrix A, SparseMatrix *cA, real *node_wgt, if (B) SparseMatrix_delete(B); } +void Multilevel_coarsen(SparseMatrix A, SparseMatrix *cA, SparseMatrix D, SparseMatrix *cD, real *node_wgt, real **cnode_wgt, + SparseMatrix *P, SparseMatrix *R, Multilevel_control ctrl, int *coarsen_scheme_used){ + SparseMatrix cA0 = A, cD0 = NULL, P0 = NULL, R0 = NULL, M; + real *cnode_wgt0 = NULL; + int nc = 0, n; + + *P = NULL; *R = NULL; *cA = NULL; *cnode_wgt = NULL, *cD = NULL; + + n = A->n; + + do {/* this loop force a sufficient reduction */ + node_wgt = cnode_wgt0; + Multilevel_coarsen_internal(A, &cA0, D, &cD0, node_wgt, &cnode_wgt0, &P0, &R0, ctrl, coarsen_scheme_used); + if (!cA0) return; + nc = cA0->n; +#ifdef DEBUG_PRINT + if (Verbose) fprintf(stderr,"nc=%d n = %d\n",nc,n); +#endif + if (*P){ + assert(*R); + M = SparseMatrix_multiply(*P, P0); + SparseMatrix_delete(*P); + SparseMatrix_delete(P0); + *P = M; + M = SparseMatrix_multiply(R0, *R); + SparseMatrix_delete(*R); + SparseMatrix_delete(R0); + *R = M; + } else { + *P = P0; + *R = R0; + } + + if (*cA) SparseMatrix_delete(*cA); + *cA = cA0; + if (*cD) SparseMatrix_delete(*cD); + *cD = cD0; + + if (*cnode_wgt) FREE(*cnode_wgt); + *cnode_wgt = cnode_wgt0; + A = cA0; + D = cD0; + node_wgt = cnode_wgt0; + cnode_wgt0 = NULL; + } while (nc > ctrl->min_coarsen_factor*n && ctrl->coarsen_mode == COARSEN_MODE_FORCEFUL); + +} + void print_padding(int n){ int i; for (i = 0; i < n; i++) fputs (" ", stderr); @@ -1068,7 +1251,7 @@ static Multilevel Multilevel_establish(Multilevel grid, Multilevel_control ctrl) Multilevel cgrid; int coarsen_scheme_used; real *cnode_weights = NULL; - SparseMatrix P, R, A, cA; + SparseMatrix P, R, A, cA, D, cD; #ifdef DEBUG_PRINT if (Verbose) { @@ -1077,6 +1260,7 @@ static Multilevel Multilevel_establish(Multilevel grid, Multilevel_control ctrl) } #endif A = grid->A; + D = grid->D; if (grid->level >= ctrl->maxlevel - 1) { #ifdef DEBUG_PRINT if (Verbose) { @@ -1086,15 +1270,16 @@ static Multilevel Multilevel_establish(Multilevel grid, Multilevel_control ctrl) #endif return grid; } - Multilevel_coarsen(A, &cA, grid->node_weights, &cnode_weights, &P, &R, ctrl, &coarsen_scheme_used); + Multilevel_coarsen(A, &cA, D, &cD, grid->node_weights, &cnode_weights, &P, &R, ctrl, &coarsen_scheme_used); if (!cA) return grid; - cgrid = Multilevel_init(cA, cnode_weights); + cgrid = Multilevel_init(cA, cD, cnode_weights); grid->next = cgrid; cgrid->coarsen_scheme_used = coarsen_scheme_used; cgrid->level = grid->level + 1; cgrid->n = cA->m; cgrid->A = cA; + cgrid->D = cD; cgrid->P = P; grid->R = R; cgrid->prev = grid; @@ -1103,14 +1288,18 @@ static Multilevel Multilevel_establish(Multilevel grid, Multilevel_control ctrl) } -Multilevel Multilevel_new(SparseMatrix A0, real *node_weights, Multilevel_control ctrl){ +Multilevel Multilevel_new(SparseMatrix A0, SparseMatrix D0, real *node_weights, Multilevel_control ctrl){ + /* A: the weighting matrix. D: the distance matrix, could be NULL. If not null, the two matrices must have the same sparsity pattern */ Multilevel grid; - SparseMatrix A = A0; + SparseMatrix A = A0, D = D0; if (!SparseMatrix_is_symmetric(A, FALSE) || A->type != MATRIX_TYPE_REAL){ A = SparseMatrix_get_real_adjacency_matrix_symmetrized(A); } - grid = Multilevel_init(A, node_weights); + if (D && (!SparseMatrix_is_symmetric(D, FALSE) || D->type != MATRIX_TYPE_REAL)){ + D = SparseMatrix_symmetrize_nodiag(D, FALSE); + } + grid = Multilevel_init(A, D, node_weights); grid = Multilevel_establish(grid, ctrl); if (A != A0) grid->delete_top_level_A = TRUE;/* be sure to clean up later */ return grid; diff --git a/lib/sfdpgen/Multilevel.h b/lib/sfdpgen/Multilevel.h index 64314284c..1f981dce3 100644 --- a/lib/sfdpgen/Multilevel.h +++ b/lib/sfdpgen/Multilevel.h @@ -21,7 +21,9 @@ typedef struct Multilevel_struct *Multilevel; struct Multilevel_struct { int level;/* 0, 1, ... */ int n; - SparseMatrix A; + SparseMatrix A;/* the weighting matrix */ + SparseMatrix D;/* the distance matrix. A and D should have same pattern, + but different entry values. For spring-electrical method, D = NULL. */ SparseMatrix P; SparseMatrix R; real *node_weights; @@ -37,6 +39,7 @@ enum {MAX_CLUSTER_SIZE = 4}; enum {EDGE_BASED_STA, COARSEN_INDEPENDENT_EDGE_SET, COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_EDGE_PERNODE, COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_EDGE_PERNODE_LEAVES_FIRST, COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_EDGE_PERNODE_SUPERNODES_FIRST, COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_EDGE_PERNODE_DEGREE_SCALED, COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_CLUSTER_PERNODE_LEAVES_FIRST, EDGE_BASED_STO, VERTEX_BASED_STA, COARSEN_INDEPENDENT_VERTEX_SET, COARSEN_INDEPENDENT_VERTEX_SET_RS, VERTEX_BASED_STO, COARSEN_HYBRID}; +enum {COARSEN_MODE_GENTLE, COARSEN_MODE_FORCEFUL}; struct Multilevel_control_struct { int minsize; @@ -44,17 +47,18 @@ struct Multilevel_control_struct { int maxlevel; int randomize; int coarsen_scheme; + int coarsen_mode; }; typedef struct Multilevel_control_struct *Multilevel_control; -Multilevel_control Multilevel_control_new(void); +Multilevel_control Multilevel_control_new(int scheme, int mode); void Multilevel_control_delete(Multilevel_control ctrl); void Multilevel_delete(Multilevel grid); -Multilevel Multilevel_new(SparseMatrix A, real *node_weights, Multilevel_control ctrl); +Multilevel Multilevel_new(SparseMatrix A, SparseMatrix D, real *node_weights, Multilevel_control ctrl); Multilevel Multilevel_get_coarsest(Multilevel grid); @@ -63,4 +67,6 @@ void print_padding(int n); #define Multilevel_is_finest(grid) (!((grid)->prev)) #define Multilevel_is_coarsest(grid) (!((grid)->next)) +void Multilevel_coarsen(SparseMatrix A, SparseMatrix *cA, SparseMatrix D, SparseMatrix *cD, real *node_wgt, real **cnode_wgt, + SparseMatrix *P, SparseMatrix *R, Multilevel_control ctrl, int *coarsen_scheme_used); #endif diff --git a/lib/sfdpgen/QuadTree.c b/lib/sfdpgen/QuadTree.c index 58564189d..cde02b0ec 100644 --- a/lib/sfdpgen/QuadTree.c +++ b/lib/sfdpgen/QuadTree.c @@ -11,17 +11,12 @@ * Contributors: See CVS logs. Details at http://www.graphviz.org/ *************************************************************************/ -#include -typedef double real; +#include "general.h" #include "geom.h" #include "arith.h" -#include "QuadTree.h" -#include "memory.h" #include "math.h" -#define MALLOC gmalloc -#define REALLOC grealloc -#define FREE free -#define MEMCPY memcpy +#include "LinkedList.h" +#include "QuadTree.h" extern real distance_cropped(real *x, int dim, int i, int j); @@ -34,20 +29,13 @@ struct node_data_struct { typedef struct node_data_struct *node_data; -real point_distance(real *p1, real *p2, int dim){ - int i; - real dist; - dist = 0; - for (i = 0; i < dim; i++) dist += (p1[i] - p2[i])*(p1[i] - p2[i]); - return sqrt(dist); -} static node_data node_data_new(int dim, real weight, real *coord, int id){ node_data nd; int i; - nd = N_GNEW(1,struct node_data_struct); + nd = MALLOC(sizeof(struct node_data_struct)); nd->node_weight = weight; - nd->coord = N_GNEW(dim,real); + nd->coord = MALLOC(sizeof(real)*dim); nd->id = id; for (i = 0; i < dim; i++) nd->coord[i] = coord[i]; nd->data = NULL; @@ -73,7 +61,7 @@ real* node_data_get_coord(void *d){ int node_data_get_id(void *d){ node_data nd = (node_data) d; - return (int)nd->id; + return nd->id; } #define node_data_get_data(d) (((node_data) (d))->data) @@ -83,9 +71,9 @@ void check_or_realloc_arrays(int dim, int *nsuper, int *nsupermax, real **center if (*nsuper >= *nsupermax) { *nsupermax = *nsuper + MAX(10, (int) 0.2*(*nsuper)); - *center = RALLOC((*nsupermax)*dim,*center,real); - *supernode_wgts = RALLOC((*nsupermax),*supernode_wgts,real); - *distances = RALLOC((*nsupermax),*distances,real); + *center = REALLOC(*center, sizeof(real)*(*nsupermax)*dim); + *supernode_wgts = REALLOC(*supernode_wgts, sizeof(real)*(*nsupermax)); + *distances = REALLOC(*distances, sizeof(real)*(*nsupermax)); } } @@ -145,9 +133,9 @@ void QuadTree_get_supernodes(QuadTree qt, real bh, real *point, int nodeid, int *flag = 0; *nsupermax = 10; - if (!*center) *center = N_GNEW((*nsupermax)*dim,real); - if (!*supernode_wgts) *supernode_wgts = N_GNEW((*nsupermax),real); - if (!*distances) *distances = N_GNEW((*nsupermax),real); + if (!*center) *center = MALLOC(sizeof(real)*(*nsupermax)*dim); + if (!*supernode_wgts) *supernode_wgts = MALLOC(sizeof(real)*(*nsupermax)); + if (!*distances) *distances = MALLOC(sizeof(real)*(*nsupermax)); QuadTree_get_supernodes_internal(qt, bh, point, nodeid, nsuper, nsupermax, center, supernode_wgts, distances, counts, flag); } @@ -167,7 +155,7 @@ static real *get_or_alloc_force_qt(QuadTree qt, int dim){ int i; real *force = (real*) qt->data; if (!force){ - qt->data = N_GNEW(dim,real); + qt->data = MALLOC(sizeof(real)*dim); force = (real*) qt->data; for (i = 0; i < dim; i++) force[i] = 0.; } @@ -326,7 +314,7 @@ static void QuadTree_repulsive_force_accumulate(QuadTree qt, real *force, real * } -void QuadTree_get_repulvice_force(QuadTree qt, real *force, real *x, real bh, real p, real KP, real *counts, int *flag){ +void QuadTree_get_repulsive_force(QuadTree qt, real *force, real *x, real bh, real p, real KP, real *counts, int *flag){ /* get repulsice force by a more efficient algortihm: we consider two cells, if they are well separated, we calculate the overall repulsive force on the cell level, if not well separated, we divide one of the cell. If both cells are at the leaf level, we calcuaulate repulsicve force among individual nodes. Finally @@ -365,9 +353,9 @@ QuadTree QuadTree_new_from_point_list(int dim, int n, int max_level, real *coord QuadTree qt = NULL; int i, k; - xmin = N_GNEW(dim,real); - xmax = N_GNEW(dim,real); - center = N_GNEW(dim,real); + xmin = MALLOC(sizeof(real)*dim); + xmax = MALLOC(sizeof(real)*dim); + center = MALLOC(sizeof(real)*dim); if (!xmin || !xmax || !center) return NULL; for (i = 0; i < dim; i++) xmin[i] = coord[i]; @@ -408,10 +396,10 @@ QuadTree QuadTree_new_from_point_list(int dim, int n, int max_level, real *coord QuadTree QuadTree_new(int dim, real *center, real width, int max_level){ QuadTree q; int i; - q = N_GNEW(1,struct QuadTree_struct); + q = MALLOC(sizeof(struct QuadTree_struct)); q->dim = dim; q->n = 0; - q->center = N_GNEW(dim,real); + q->center = MALLOC(sizeof(real)*dim); for (i = 0; i < dim; i++) q->center[i] = center[i]; assert(width > 0); q->width = width; @@ -489,14 +477,20 @@ static QuadTree QuadTree_add_internal(QuadTree q, real *coord, real weight, int /* Make sure that coord is within bounding box */ for (i = 0; i < q->dim; i++) { - if (coord[i] < q->center[i] - q->width || coord[i] > q->center[i] + q->width) return NULL; + if (coord[i] < q->center[i] - q->width - 1.e5*MACHINEACC || coord[i] > q->center[i] + q->width + 1.e5*MACHINEACC) { +#ifdef DEBUG_PRINT + fprintf(stderr,"(q->center[i] - q->width) - coord[i] =%g, coord[i]-(q->center[i] + q->width) = %g\n", + (q->center[i] - q->width) - coord[i], coord[i]-(q->center[i] + q->width)); +#endif + return NULL; + } } if (q->n == 0){ /* if this node is empty right now */ q->n = 1; q->total_weight = weight; - q->average = N_GNEW(dim,real); + q->average = MALLOC(sizeof(real)*dim); for (i = 0; i < q->dim; i++) q->average[i] = coord[i]; nd = node_data_new(q->dim, weight, coord, id); assert(!(q->l)); @@ -506,7 +500,7 @@ static QuadTree QuadTree_add_internal(QuadTree q, real *coord, real weight, int q->total_weight += weight; for (i = 0; i < q->dim; i++) q->average[i] = ((q->average[i])*q->n + coord[i])/(q->n + 1); if (!q->qts){ - q->qts = N_GNEW((1<qts = MALLOC(sizeof(QuadTree)*(1<qts[i] = NULL; }/* done adding new quadtree, now add points to them */ @@ -667,3 +661,71 @@ void QuadTree_print(FILE *fp, QuadTree q){ fprintf(fp, "}, PlotRange -> All]\n"); } } + + + + +static void QuadTree_get_nearest_internal(QuadTree qt, real *x, real *y, real *min, int *imin, int tentative, int *flag){ + /* get the narest point years to {x[0], ..., x[dim]} and store in y.*/ + SingleLinkedList l; + real *coord, dist; + int dim, i, iq = -1; + real qmin; + real *point = x; + + *flag = 0; + if (!qt) return; + dim = qt->dim; + l = qt->l; + if (l){ + while (l){ + coord = node_data_get_coord(SingleLinkedList_get_data(l)); + dist = point_distance(point, coord, dim); + if(*min < 0 || dist < *min) { + *min = dist; + *imin = node_data_get_id(SingleLinkedList_get_data(l)); + for (i = 0; i < dim; i++) y[i] = coord[i]; + } + l = SingleLinkedList_get_next(l); + } + } + + if (qt->qts){ + dist = point_distance(qt->center, point, dim); + if (*min >= 0 && (dist - sqrt((real) dim) * qt->width > *min)){ + return; + } else { + if (tentative){/* quick first approximation*/ + qmin = -1; + for (i = 0; i < 1<qts[i]){ + dist = point_distance(qt->qts[i]->average, point, dim); + if (dist < qmin || qmin < 0){ + qmin = dist; iq = i; + } + } + } + assert(iq >= 0); + QuadTree_get_nearest_internal(qt->qts[iq], x, y, min, imin, tentative, flag); + } else { + for (i = 0; i < 1<qts[i], x, y, min, imin, tentative, flag); + } + } + } + } + +} + + +void QuadTree_get_nearest(QuadTree qt, real *x, real *ymin, int *imin, real *min, int *flag){ + + *flag = 0; + *min = -1; + + QuadTree_get_nearest_internal(qt, x, ymin, min, imin, TRUE, flag); + + QuadTree_get_nearest_internal(qt, x, ymin, min, imin, FALSE, flag); + + +} diff --git a/lib/sfdpgen/QuadTree.h b/lib/sfdpgen/QuadTree.h index e71bf1536..bbd904859 100644 --- a/lib/sfdpgen/QuadTree.h +++ b/lib/sfdpgen/QuadTree.h @@ -15,7 +15,7 @@ #define QUAD_TREE_H #include "LinkedList.h" -#include "sfdpinternal.h" +/* #include "sfdpinternal.h" */ #include typedef struct QuadTree_struct *QuadTree; @@ -54,5 +54,10 @@ real point_distance(real *p1, real *p2, int dim); void QuadTree_get_supernodes(QuadTree qt, real bh, real *point, int nodeid, int *nsuper, int *nsupermax, real **center, real **supernode_wgts, real **distances, real *counts, int *flag); -void QuadTree_get_repulvice_force(QuadTree qt, real *force, real *x, real bh, real p, real KP, real *counts, int *flag); +void QuadTree_get_repulsive_force(QuadTree qt, real *force, real *x, real bh, real p, real KP, real *counts, int *flag); + +/* find the nearest point and put in ymin, index in imin and distance in min */ +void QuadTree_get_nearest(QuadTree qt, real *x, real *ymin, int *imin, real *min, int *flag); + + #endif diff --git a/lib/sfdpgen/post_process.c b/lib/sfdpgen/post_process.c index 3d9fd16c7..0beeb3875 100644 --- a/lib/sfdpgen/post_process.c +++ b/lib/sfdpgen/post_process.c @@ -25,22 +25,13 @@ #include "sparse_solve.h" #include "post_process.h" +#include "overlap.h" #include "spring_electrical.h" #include "call_tri.h" #include "sfdpinternal.h" -#define MALLOC gmalloc -#define REALLOC grealloc -#define FREE free -#define MEMCPY memcpy - - - - #define node_degree(i) (ia[(i)+1] - ia[(i)]) - - SparseMatrix ideal_distance_matrix(SparseMatrix A, int dim, real *x){ /* find the ideal distance between edges, either 1, or |N[i] \Union N[j]| - |N[i] \Intersection N[j]| */ @@ -109,8 +100,9 @@ SparseMatrix ideal_distance_matrix(SparseMatrix A, int dim, real *x){ return D; } -StressMajorizationSmoother StressMajorizationSmoother_new(SparseMatrix A, int dim, real lambda0, real *x, +StressMajorizationSmoother StressMajorizationSmoother2_new(SparseMatrix A, int dim, real lambda0, real *x, int ideal_dist_scheme){ + /* use up to dist 2 neighbor */ StressMajorizationSmoother sm; int i, j, k, l, m = A->m, *ia = A->ia, *ja = A->ja, *iw, *jw, *id, *jd; int *mask, nz; @@ -124,6 +116,10 @@ StressMajorizationSmoother StressMajorizationSmoother_new(SparseMatrix A, int di dd = (real*) ID->a; sm = N_GNEW(1,struct StressMajorizationSmoother_struct); + sm->scaling = 1.; + sm->data = NULL; + sm->scheme = SM_SCHEME_NORMAL; + lambda = sm->lambda = N_GNEW(m,real); for (i = 0; i < m; i++) sm->lambda[i] = lambda0; mask = N_GNEW(m,int); @@ -223,6 +219,7 @@ StressMajorizationSmoother StressMajorizationSmoother_new(SparseMatrix A, int di } } + /* distance 2 neighbors */ for (j = ia[i]; j < ia[i+1]; j++){ k = ja[j]; for (l = ia[k]; l < ia[k+1]; l++){ @@ -280,6 +277,7 @@ StressMajorizationSmoother StressMajorizationSmoother_new(SparseMatrix A, int di s = stop/sbot; for (i = 0; i < nz; i++) d[i] *= s; + sm->scaling = s; sm->Lw->nz = nz; sm->Lwd->nz = nz; @@ -288,6 +286,113 @@ StressMajorizationSmoother StressMajorizationSmoother_new(SparseMatrix A, int di SparseMatrix_delete(ID); return sm; } + +StressMajorizationSmoother SparseStressMajorizationSmoother_new(SparseMatrix A, int dim, real lambda0, real *x, int weighting_scheme){ + /* solve a stress model to achieve the ideal distance among a sparse set of edges recorded in A. + A must be a real matrix. + */ + StressMajorizationSmoother sm; + int i, j, k, m = A->m, *ia, *ja, *iw, *jw, *id, *jd; + int nz; + real *d, *w, *lambda; + real diag_d, diag_w, *a, dist, s = 0, stop = 0, sbot = 0; + real xdot = 0; + + assert(SparseMatrix_is_symmetric(A, FALSE) && A->type == MATRIX_TYPE_REAL); + + /* if x is all zero, make it random */ + for (i = 0; i < m*dim; i++) xdot += x[i]*x[i]; + if (xdot == 0){ + for (i = 0; i < m*dim; i++) x[i] = 72*drand(); + } + + ia = A->ia; + ja = A->ja; + a = (real*) A->a; + + + sm = MALLOC(sizeof(struct StressMajorizationSmoother_struct)); + sm->scaling = 1.; + sm->data = NULL; + sm->scheme = SM_SCHEME_NORMAL; + + lambda = sm->lambda = MALLOC(sizeof(real)*m); + for (i = 0; i < m; i++) sm->lambda[i] = lambda0; + + nz = A->nz; + + sm->Lw = SparseMatrix_new(m, m, nz + m, MATRIX_TYPE_REAL, FORMAT_CSR); + sm->Lwd = SparseMatrix_new(m, m, nz + m, MATRIX_TYPE_REAL, FORMAT_CSR); + if (!(sm->Lw) || !(sm->Lwd)) { + StressMajorizationSmoother_delete(sm); + return NULL; + } + + iw = sm->Lw->ia; jw = sm->Lw->ja; + id = sm->Lwd->ia; jd = sm->Lwd->ja; + w = (real*) sm->Lw->a; d = (real*) sm->Lwd->a; + iw[0] = id[0] = 0; + + nz = 0; + for (i = 0; i < m; i++){ + diag_d = diag_w = 0; + for (j = ia[i]; j < ia[i+1]; j++){ + k = ja[j]; + if (k != i){ + + jw[nz] = k; + dist = a[j]; + switch (weighting_scheme){ + case WEIGHTING_SCHEME_SQR_DIST: + if (dist*dist == 0){ + w[nz] = -100000; + } else { + w[nz] = -1/(dist*dist); + } + break; + case WEIGHTING_SCHEME_NONE: + w[nz] = -1; + break; + default: + assert(0); + return NULL; + } + diag_w += w[nz]; + jd[nz] = k; + d[nz] = w[nz]*dist; + + stop += d[nz]*distance(x,dim,i,k); + sbot += d[nz]*dist; + diag_d += d[nz]; + + nz++; + } + } + + jw[nz] = i; + lambda[i] *= (-diag_w);/* alternatively don't do that then we have a constant penalty term scaled by lambda0 */ + w[nz] = -diag_w + lambda[i]; + + jd[nz] = i; + d[nz] = -diag_d; + nz++; + + iw[i+1] = nz; + id[i+1] = nz; + } + s = stop/sbot; + if (s == 0) { + return NULL; + } + for (i = 0; i < nz; i++) d[i] *= s; + + + sm->scaling = s; + sm->Lw->nz = nz; + sm->Lwd->nz = nz; + + return sm; +} static real total_distance(int m, int dim, real* x, real* y){ real total = 0, dist = 0; @@ -306,10 +411,172 @@ static real total_distance(int m, int dim, real* x, real* y){ -void StressMajorizationSmoother_smooth(StressMajorizationSmoother sm, int dim, real *x, int maxit_sm) { +void SparseStressMajorizationSmoother_delete(SparseStressMajorizationSmoother sm){ + return StressMajorizationSmoother_delete(sm); +} + + +real SparseStressMajorizationSmoother_smooth(SparseStressMajorizationSmoother sm, int dim, real *x, int maxit_sm, real tol){ + + return StressMajorizationSmoother_smooth(sm, dim, x, maxit_sm, tol); + + +} +static void get_edge_label_matrix(relative_position_constraints data, int m, int dim, real *x, SparseMatrix *LL, real **rhs){ + int edge_labeling_scheme = data->edge_labeling_scheme; + int n_constr_nodes = data->n_constr_nodes; + int *constr_nodes = data->constr_nodes; + SparseMatrix A_constr = data->A_constr; + int *ia = A_constr->ia, *ja = A_constr->ja, ii, jj, nz, l, ll, i, j; + int *irn = data->irn, *jcn = data->jcn; + real *val = data->val, dist, kk, k; + real *x00 = NULL; + SparseMatrix Lc = NULL; + real constr_penalty = data->constr_penalty; + + if (edge_labeling_scheme == ELSCHEME_PENALTY || edge_labeling_scheme == ELSCHEME_STRAIGHTLINE_PENALTY){ + /* for an node with two neighbors j--i--k, and assume i needs to be between j and k, then the contribution to P is + . i j k + i 1 -0.5 -0.5 + j -0.5 0.25 0.25 + k -0.5 0.25 0.25 + in general, if i has m neighbors j, k, ..., then + p_ii = 1 + p_jj = 1/m^2 + p_ij = -1/m + p jk = 1/m^2 + . i j k + i 1 -1/m -1/m ... + j -1/m 1/m^2 1/m^2 ... + k -1/m 1/m^2 1/m^2 ... + */ + if (!irn){ + assert((!jcn) && (!val)); + nz = 0; + for (i = 0; i < n_constr_nodes; i++){ + ii = constr_nodes[i]; + k = ia[ii+1] - ia[ii];/*usually k = 2 */ + nz += (k+1)*(k+1); + + } + irn = data->irn = MALLOC(sizeof(int)*nz); + jcn = data->jcn = MALLOC(sizeof(int)*nz); + val = data->val = MALLOC(sizeof(double)*nz); + } + nz = 0; + for (i = 0; i < n_constr_nodes; i++){ + ii = constr_nodes[i]; + jj = ja[ia[ii]]; ll = ja[ia[ii] + 1]; + if (jj == ll) continue; /* do not do loops */ + dist = distance_cropped(x, dim, jj, ll); + dist *= dist; + + k = ia[ii+1] - ia[ii];/* usually k = 2 */ + kk = k*k; + irn[nz] = ii; jcn[nz] = ii; val[nz++] = constr_penalty/(dist); + k = constr_penalty/(k*dist); kk = constr_penalty/(kk*dist); + for (j = ia[ii]; j < ia[ii+1]; j++){ + irn[nz] = ii; jcn[nz] = ja[j]; val[nz++] = -k; + } + for (j = ia[ii]; j < ia[ii+1]; j++){ + jj = ja[j]; + irn[nz] = jj; jcn[nz] = ii; val[nz++] = -k; + for (l = ia[ii]; l < ia[ii+1]; l++){ + ll = ja[l]; + irn[nz] = jj; jcn[nz] = ll; val[nz++] = kk; + } + } + } + Lc = SparseMatrix_from_coordinate_arrays(nz, m, m, irn, jcn, val, MATRIX_TYPE_REAL); + } else if (edge_labeling_scheme == ELSCHEME_PENALTY2 || edge_labeling_scheme == ELSCHEME_STRAIGHTLINE_PENALTY2){ + /* for an node with two neighbors j--i--k, and assume i needs to be between the old position of j and k, then the contribution to P is + 1/d_jk, and to the right hand side: {0,...,average_position_of_i's neighbor if i is an edge node,...} + */ + if (!irn){ + assert((!jcn) && (!val)); + nz = n_constr_nodes; + irn = data->irn = MALLOC(sizeof(int)*nz); + jcn = data->jcn = MALLOC(sizeof(int)*nz); + val = data->val = MALLOC(sizeof(double)*nz); + } + x00 = MALLOC(sizeof(real)*m*dim); + for (i = 0; i < m*dim; i++) x00[i] = 0; + nz = 0; + for (i = 0; i < n_constr_nodes; i++){ + ii = constr_nodes[i]; + jj = ja[ia[ii]]; ll = ja[ia[ii] + 1]; + dist = distance_cropped(x, dim, jj, ll); + irn[nz] = ii; jcn[nz] = ii; val[nz++] = constr_penalty/(dist); + for (j = ia[ii]; j < ia[ii+1]; j++){ + jj = ja[j]; + for (l = 0; l < dim; l++){ + x00[ii*dim+l] += x[jj*dim+l]; + } + } + for (l = 0; l < dim; l++) { + x00[ii*dim+l] *= constr_penalty/(dist)/(ia[ii+1] - ia[ii]); + } + } + Lc = SparseMatrix_from_coordinate_arrays(nz, m, m, irn, jcn, val, MATRIX_TYPE_REAL); + + } + *LL = Lc; + *rhs = x00; +} + +real get_stress(int m, int dim, int *iw, int *jw, real *w, real *d, real *x, real scaling, void *data, int weighted){ + int i, j; + real res = 0., dist; + /* we use the fact that d_ij = w_ij*dist(i,j). Also, d_ij and x are scalinged by *scaling, so divide by it to get actual unscaled streee. */ + for (i = 0; i < m; i++){ + for (j = iw[i]; j < iw[i+1]; j++){ + if (i == jw[j]) { + continue; + } + dist = d[j]/w[j];/* both negative*/ + if (weighted){ + res += -w[j]*(dist - distance(x, dim, i, jw[j]))*(dist - distance(x, dim, i, jw[j])); + } else { + res += (dist - distance(x, dim, i, jw[j]))*(dist - distance(x, dim, i, jw[j])); + } + } + } + return res/scaling/scaling; + +} + +static void uniform_stress_augment_rhs(int m, int dim, real *x, real *y, real alpha, real M){ + int i, j, k; + real dist, distij; + for (i = 0; i < m; i++){ + for (j = i+1; j < m; j++){ + dist = distance_cropped(x, dim, i, j); + for (k = 0; k < dim; k++){ + distij = (x[i*dim+k] - x[j*dim+k])/dist; + y[i*dim+k] += alpha*M*distij; + y[j*dim+k] += alpha*M*(-distij); + } + } + } +} + +static real uniform_stress_solve(SparseMatrix Lw, real alpha, int dim, real *x0, real *rhs, real tol, int maxit, int *flag){ + Operator Ax; + Operator Precon; + + Ax = Operator_uniform_stress_matmul(Lw, alpha); + Precon = Operator_uniform_stress_diag_precon_new(Lw, alpha); + + return cg(Ax, Precon, Lw->m, dim, x0, rhs, tol, maxit, flag); + +} + + +real StressMajorizationSmoother_smooth(StressMajorizationSmoother sm, int dim, real *x, int maxit_sm, real tol) { SparseMatrix Lw = sm->Lw, Lwd = sm->Lwd, Lwdd = NULL; - int i, j, m, *id, *jd, idiag, flag = 0, iter = 0; - real *dd, *d, *y = NULL, *x0 = NULL, diag, diff = 1, tol = 0.001, *lambda = sm->lambda, maxit, res; + int i, j, m, *id, *jd, *iw, *jw, idiag, flag = 0, iter = 0; + real *w, *dd, *d, *y = NULL, *x0 = NULL, *x00 = NULL, diag, diff = 1, *lambda = sm->lambda, maxit, res, alpha = 0., M = 0.; + SparseMatrix Lc = NULL; Lwdd = SparseMatrix_copy(Lwd); m = Lw->m; @@ -323,8 +590,20 @@ void StressMajorizationSmoother_smooth(StressMajorizationSmoother sm, int dim, r id = Lwd->ia; jd = Lwd->ja; d = (real*) Lwd->a; dd = (real*) Lwdd->a; + w = (real*) Lw->a; + iw = Lw->ia; jw = Lw->ja; + + /* for the additional matrix L due to the position constraints */ + if (sm->scheme == SM_SCHEME_NORMAL_ELABEL){ + get_edge_label_matrix(sm->data, m, dim, x, &Lc, &x00); + if (Lc) Lw = SparseMatrix_add(Lw, Lc); + } else if (sm->scheme == SM_SCHEME_UNIFORM_STRESS){ + alpha = ((real*) (sm->data))[0]; + M = ((real*) (sm->data))[1]; + } while (iter++ < maxit_sm && diff > tol){ + for (i = 0; i < m; i++){ idiag = -1; diag = 0.; @@ -339,21 +618,45 @@ void StressMajorizationSmoother_smooth(StressMajorizationSmoother sm, int dim, r assert(idiag >= 0); dd[idiag] = -diag; } + /* solve (Lw+lambda*I) x = Lwdd y + lambda x0 */ SparseMatrix_multiply_dense(Lwdd, FALSE, x, FALSE, &y, FALSE, dim); - for (i = 0; i < m; i++){ - for (j = 0; j < dim; j++){ - y[i*dim+j] += lambda[i]*x0[i*dim+j]; + + if (lambda){/* is there a penalty term? */ + for (i = 0; i < m; i++){ + for (j = 0; j < dim; j++){ + y[i*dim+j] += lambda[i]*x0[i*dim+j]; + } } } - maxit = sqrt((double) m); - res = SparseMatrix_solve(Lw, dim, x, y, 0.01, (int)maxit, SOLVE_METHOD_CG, &flag); + if (sm->scheme == SM_SCHEME_NORMAL_ELABEL){ + for (i = 0; i < m; i++){ + for (j = 0; j < dim; j++){ + y[i*dim+j] += x00[i*dim+j]; + } + } + } else if (sm->scheme == SM_SCHEME_UNIFORM_STRESS){/* this part can be done more efficiently using octree approximation */ + uniform_stress_augment_rhs(m, dim, x, y, alpha, M); + } + +#ifdef DEBUG_PRINT + if (Verbose) fprintf(stderr, "stress1 = %f\n",get_stress(m, dim, iw, jw, w, d, x, sm->scaling, sm->data, 0)); +#endif + maxit = sqrt((double) m); + if (sm->scheme == SM_SCHEME_UNIFORM_STRESS){ + res = uniform_stress_solve(Lw, alpha, dim, x, y, 0.01, maxit, &flag); + } else { + res = SparseMatrix_solve(Lw, dim, x, y, 0.01, maxit, SOLVE_METHOD_CG, &flag); + } if (flag) goto RETURN; +#ifdef DEBUG_PRINT + if (Verbose) fprintf(stderr, "stress2 = %f\n",get_stress(m, dim, iw, jw, w, d, y, sm->scaling, sm->data, 0)); +#endif - diff = total_distance(m, dim, x, y)/m; + diff = total_distance(m, dim, x, y)/sqrt(vector_product(m*dim, x, x)); #ifdef DEBUG_PRINT if (Verbose){ fprintf(stderr, "Outer iter = %d, res = %g Stress Majorization diff = %g\n",iter, res, diff); @@ -368,8 +671,16 @@ void StressMajorizationSmoother_smooth(StressMajorizationSmoother sm, int dim, r RETURN: SparseMatrix_delete(Lwdd); - if (!x0) FREE(x0); - if (!y) FREE(y); + if (Lc) { + SparseMatrix_delete(Lc); + SparseMatrix_delete(Lw); + } + + if (x0) FREE(x0); + if (y) FREE(y); + if (x00) FREE(x00); + return diff; + } void StressMajorizationSmoother_delete(StressMajorizationSmoother sm){ @@ -377,6 +688,8 @@ void StressMajorizationSmoother_delete(StressMajorizationSmoother sm){ if (sm->Lw) SparseMatrix_delete(sm->Lw); if (sm->Lwd) SparseMatrix_delete(sm->Lwd); if (sm->lambda) FREE(sm->lambda); + if (sm->data) sm->data_deallocator(sm->data); + FREE(sm); } @@ -404,6 +717,10 @@ TriangleSmoother TriangleSmoother_new(SparseMatrix A, int dim, real lambda0, rea } sm = N_GNEW(1,struct TriangleSmoother_struct); + sm->scaling = 1; + sm->data = NULL; + sm->scheme = SM_SCHEME_NORMAL; + lambda = sm->lambda = N_GNEW(m,real); for (i = 0; i < m; i++) sm->lambda[i] = lambda0; @@ -417,12 +734,6 @@ TriangleSmoother TriangleSmoother_new(SparseMatrix A, int dim, real lambda0, rea B = SparseMatrix_copy(A); } - /* - {FILE *fp; - fp = fopen("/tmp/111","w"); - export_embedding(fp, dim, B, x, NULL); - fclose(fp);} - */ sm->Lw = SparseMatrix_add(A, B); @@ -472,6 +783,7 @@ TriangleSmoother TriangleSmoother_new(SparseMatrix A, int dim, real lambda0, rea s = stop/sbot; for (i = 0; i < iw[m]; i++) d[i] *= s; + sm->scaling = s; FREE(avg_dist); @@ -486,7 +798,7 @@ void TriangleSmoother_delete(TriangleSmoother sm){ void TriangleSmoother_smooth(TriangleSmoother sm, int dim, real *x){ - StressMajorizationSmoother_smooth(sm, dim, x, 50); + StressMajorizationSmoother_smooth(sm, dim, x, 50, 0.001); } @@ -660,8 +972,8 @@ void post_process_smoothing(int dim, SparseMatrix A, spring_electrical_control c } for (k = 0; k < 1; k++){ - sm = StressMajorizationSmoother_new(A, dim, 0.05, x, dist_scheme); - StressMajorizationSmoother_smooth(sm, dim, x, 50); + sm = StressMajorizationSmoother2_new(A, dim, 0.05, x, dist_scheme); + StressMajorizationSmoother_smooth(sm, dim, x, 50, 0.001); StressMajorizationSmoother_delete(sm); } break; diff --git a/lib/sfdpgen/post_process.h b/lib/sfdpgen/post_process.h index 745b1b744..970148f40 100644 --- a/lib/sfdpgen/post_process.h +++ b/lib/sfdpgen/post_process.h @@ -16,10 +16,16 @@ #include "spring_electrical.h" +enum {SM_SCHEME_NORMAL, SM_SCHEME_NORMAL_ELABEL, SM_SCHEME_UNIFORM_STRESS}; + struct StressMajorizationSmoother_struct { SparseMatrix Lw; SparseMatrix Lwd; real* lambda; + void (*data_deallocator)(void*); + void *data; + int scheme; + real scaling;/* scaling applied to the distance. need to divide coordinate at the end of the stress majorization process */ }; typedef struct StressMajorizationSmoother_struct *StressMajorizationSmoother; @@ -27,9 +33,9 @@ typedef struct StressMajorizationSmoother_struct *StressMajorizationSmoother; void StressMajorizationSmoother_delete(StressMajorizationSmoother sm); enum {IDEAL_GRAPH_DIST, IDEAL_AVG_DIST, IDEAL_POWER_DIST}; -StressMajorizationSmoother StressMajorizationSmoother_new(SparseMatrix A, int dim, real lambda, real *x, int ideal_dist_scheme); +StressMajorizationSmoother StressMajorizationSmoother2_new(SparseMatrix A, int dim, real lambda, real *x, int ideal_dist_scheme); -void StressMajorizationSmoother_smooth(StressMajorizationSmoother sm, int dim, real *x, int maxit); +real StressMajorizationSmoother_smooth(StressMajorizationSmoother sm, int dim, real *x, int maxit, real tol); /*-------------------- triangle/neirhborhood graph based smoother ------------------- */ typedef StressMajorizationSmoother TriangleSmoother; @@ -60,4 +66,21 @@ void SpringSmoother_smooth(SpringSmoother sm, SparseMatrix A, real *node_weights /*------------------------------------------------------------------*/ void post_process_smoothing(int dim, SparseMatrix A, spring_electrical_control ctrl, real *node_weights, real *x, int *flag); + +/*-------------------- sparse stress majorizationp ------------------- */ +typedef StressMajorizationSmoother SparseStressMajorizationSmoother; + +#define SparseStressMajorizationSmoother_struct StressMajorizationSmoother_struct + +void SparseStressMajorizationSmoother_delete(SparseStressMajorizationSmoother sm); + +enum {WEIGHTING_SCHEME_NONE, WEIGHTING_SCHEME_SQR_DIST}; +SparseStressMajorizationSmoother SparseStressMajorizationSmoother_new(SparseMatrix A, int dim, real lambda, real *x, int weighting_scheme); + +real SparseStressMajorizationSmoother_smooth(SparseStressMajorizationSmoother sm, int dim, real *x, int maxit_sm, real tol); + +real get_stress(int m, int dim, int *iw, int *jw, real *w, real *d, real *x, real scaling, void *data, int weighted); +/*--------------------------------------------------------------*/ + #endif + diff --git a/lib/sfdpgen/sfdpinit.c b/lib/sfdpgen/sfdpinit.c index 77c26e2e1..959fff49b 100644 --- a/lib/sfdpgen/sfdpinit.c +++ b/lib/sfdpgen/sfdpinit.c @@ -30,10 +30,8 @@ #include #include #include - -#ifdef DEBUG -double _statistics[10]; -#endif +#include +#include static void sfdp_init_edge(edge_t * e) { @@ -65,7 +63,6 @@ static void sfdp_init_node_edge(graph_t * g) static void sfdp_init_graph(Agraph_t * g) { int outdim; - int temp; setEdgeType(g, ET_LINE); outdim = late_int(g, agfindgraphattr(g, "dimen"), 2, 2); @@ -86,17 +83,13 @@ static real *getPos(Agraph_t * g, spring_electrical_control ctrl) if (agfindnodeattr(g, "pos") == NULL) return pos; - ctrl->posSet = N_NEW(agnnodes(g), char); for (n = agfstnode(g); n; n = agnxtnode(g, n)) { i = ND_id(n); if (hasPos(n)) { for (ix = 0; ix < Ndim; ix++) { pos[i * Ndim + ix] = ND_pos(n)[ix]; } - ctrl->posSet[i] = 1; } - else - ctrl->posSet[i] = 0; } return pos; @@ -109,15 +102,40 @@ static void sfdpLayout(graph_t * g, spring_electrical_control ctrl, real *pos; Agnode_t *n; int flag, i; + int n_edge_label_nodes = 0, *edge_label_nodes = NULL; + SparseMatrix D = NULL; + SparseMatrix A; + + if (ctrl->method == METHOD_SPRING_MAXENT) /* maxent can work with distance matrix */ + A = makeMatrix(g, Ndim, &D); + else + A = makeMatrix(g, Ndim, NULL); - SparseMatrix A = makeMatrix(g, Ndim); - if (ctrl->overlap >= 0) - sizes = getSizes(g, pad); + if (ctrl->overlap >= 0) { + if (ctrl->edge_labeling_scheme > 0) + sizes = getSizes(g, pad, &n_edge_label_nodes, &edge_label_nodes); + else + sizes = getSizes(g, pad, NULL, NULL); + } else sizes = NULL; pos = getPos(g, ctrl); - multilevel_spring_electrical_embedding(Ndim, A, ctrl, NULL, sizes, pos, &flag); + switch (ctrl->method) { + case METHOD_SPRING_ELECTRICAL: + case METHOD_SPRING_MAXENT: + multilevel_spring_electrical_embedding(Ndim, A, D, ctrl, NULL, sizes, pos, n_edge_label_nodes, edge_label_nodes, &flag); + break; + case METHOD_UNIFORM_STRESS: + uniform_stress(Ndim, A, pos, &flag); + break; + case METHOD_STRESS:{ + int maxit = 1000; + real tol = 0.001; + stress_model(Ndim, A, &pos, maxit, tol, &flag); + } + break; + } for (n = agfstnode(g); n; n = agnxtnode(g, n)) { real *npos = pos + (Ndim * ND_id(n)); @@ -128,8 +146,47 @@ static void sfdpLayout(graph_t * g, spring_electrical_control ctrl, free(sizes); free(pos); - free(ctrl->posSet); SparseMatrix_delete (A); + if (D) SparseMatrix_delete (D); + if (edge_label_nodes) FREE(edge_label_nodes); +} + +static int +late_mode (graph_t* g, Agsym_t* sym, int dflt) +{ + char* s; + int v; + int rv; + + if (!sym) return dflt; +#ifdef WITH_CGRAPH + s = agxget (g, sym); +#else + s = agxget (g, sym->index); +#endif + if (isdigit(*s)) { + if ((v = atoi (s)) <= METHOD_UNIFORM_STRESS) + rv = v; + else + rv = dflt; + } + else if (isalpha(*s)) { + if (!strcasecmp(s, "spring")) + rv = METHOD_SPRING_ELECTRICAL; + else if (!strcasecmp(s, "maxent")) + rv = METHOD_SPRING_MAXENT; + else if (!strcasecmp(s, "stress")) + rv = METHOD_STRESS; + else if (!strcasecmp(s, "uniform")) + rv = METHOD_UNIFORM_STRESS; + else { + agerr (AGWARN, "Unknown value \"%s\" for mode attribute\n", s); + rv = dflt; + } + } + else + rv = dflt; + return rv; } static int @@ -242,13 +299,20 @@ tuneControl (graph_t* g, spring_electrical_control ctrl) ctrl->multilevels = late_int(g, agfindgraphattr(g, "levels"), INT_MAX, 0); ctrl->smoothing = late_smooth(g, agfindgraphattr(g, "smoothing"), SMOOTHING_NONE); ctrl->tscheme = late_quadtree_scheme(g, agfindgraphattr(g, "quadtree"), QUAD_TREE_NORMAL); + ctrl->method = late_mode(g, agfindgraphattr(g, "quadtree"), METHOD_SPRING_ELECTRICAL); + ctrl->rotation = late_double(g, agfindgraphattr(g, "rotation"), 0.0, -MAXDOUBLE); + ctrl->edge_labeling_scheme = late_int(g, agfindgraphattr(g, "label_scheme"), 0, 0); + if (ctrl->edge_labeling_scheme > 4) { + agerr (AGWARN, "label_scheme = %d > 4 : ignoring\n", ctrl->edge_labeling_scheme); + ctrl->edge_labeling_scheme = 0; + } } void sfdp_layout(graph_t * g) { int doAdjust; adjust_data am; - sfdp_init_graph(g); + sfdp_init_graph(g); doAdjust = (Ndim == 2); if (agnnodes(g)) { diff --git a/lib/sfdpgen/sfdpinternal.h b/lib/sfdpgen/sfdpinternal.h index 2340b2111..ada05e255 100644 --- a/lib/sfdpgen/sfdpinternal.h +++ b/lib/sfdpgen/sfdpinternal.h @@ -16,9 +16,5 @@ #include -#ifdef DEBUG -extern double _statistics[10]; -#endif - #endif diff --git a/lib/sfdpgen/sparse_solve.c b/lib/sfdpgen/sparse_solve.c index 61ae7add6..372d5eead 100644 --- a/lib/sfdpgen/sparse_solve.c +++ b/lib/sfdpgen/sparse_solve.c @@ -23,45 +23,49 @@ #include "globals.h" #define DEBUG_PRINT -#define MALLOC gmalloc -#define REALLOC grealloc -#define FREE free -#define MEMCPY memcpy - -static real *vector_subtract_to(int n, real * x, real * y) -{ - /* y = x-y */ - int i; - for (i = 0; i < n; i++) - y[i] = x[i] - y[i]; - return y; -} -static real *vector_saxpy(int n, real * x, real * y, real beta) -{ - /* y = x+beta*y */ - int i; - for (i = 0; i < n; i++) - y[i] = x[i] + beta * y[i]; - return y; +struct uniform_stress_matmul_data{ + real alpha; + SparseMatrix A; +}; + + +void Operator_uniform_stress_matmul_delete(Operator o){ + FREE(o->data); } -static real *vector_saxpy2(int n, real * x, real * y, real beta) -{ - /* x = x+beta*y */ - int i; - for (i = 0; i < n; i++) - x[i] = x[i] + beta * y[i]; - return x; +real *Operator_uniform_stress_matmul_apply(Operator o, real *x, real *y){ + struct uniform_stress_matmul_data *d = (struct uniform_stress_matmul_data*) (o->data);; + SparseMatrix A = d->A; + real alpha = d->alpha; + real xsum = 0.; + int m = A->m, i; + + SparseMatrix_multiply_vector(A, x, &y, FALSE); + + /* alpha*V*x */ + for (i = 0; i < m; i++) xsum += x[i]; + + for (i = 0; i < m; i++) y[i] += alpha*(m*x[i] - xsum); + + return y; } -static real vector_product(int n, real *x, real *y){ - real res = 0; - int i; - for (i = 0; i < n; i++) res += x[i]*y[i]; - return res; + + +Operator Operator_uniform_stress_matmul(SparseMatrix A, real alpha){ + Operator o; + struct uniform_stress_matmul_data *d; + + o = MALLOC(sizeof(struct Operator_struct)); + o->data = d = MALLOC(sizeof(struct uniform_stress_matmul_data)); + d->alpha = alpha; + d->A = A; + o->Operator_apply = Operator_uniform_stress_matmul_apply; + return o; } + real *Operator_matmul_apply(Operator o, real *x, real *y){ SparseMatrix A = (SparseMatrix) o->data; SparseMatrix_multiply_vector(A, x, &y, FALSE); @@ -79,7 +83,7 @@ Operator Operator_matmul_new(SparseMatrix A){ void Operator_matmul_delete(Operator o){ - + if (o) FREE(o); } @@ -92,6 +96,36 @@ real* Operator_diag_precon_apply(Operator o, real *x, real *y){ return y; } + +Operator Operator_uniform_stress_diag_precon_new(SparseMatrix A, real alpha){ + Operator o; + real *diag; + int i, j, m = A->m, *ia = A->ia, *ja = A->ja; + real *a = (real*) A->a; + + assert(A->type == MATRIX_TYPE_REAL); + + assert(a); + + o = MALLOC(sizeof(struct Operator_struct)); + o->data = MALLOC(sizeof(real)*(m + 1)); + diag = (real*) o->data; + + diag[0] = m; + diag++; + for (i = 0; i < m; i++){ + diag[i] = 1./(m-1); + for (j = ia[i]; j < ia[i+1]; j++){ + if (i == ja[j] && ABS(a[j]) > 0) diag[i] = 1./((m-1)*alpha+a[j]); + } + } + + o->Operator_apply = Operator_diag_precon_apply; + + return o; +} + + Operator Operator_diag_precon_new(SparseMatrix A){ Operator o; real *diag; @@ -121,7 +155,8 @@ Operator Operator_diag_precon_new(SparseMatrix A){ } void Operator_diag_precon_delete(Operator o){ - FREE(o->data); + if (o->data) FREE(o->data); + if (o) FREE(o); } static real conjugate_gradient(Operator A, Operator precon, int n, real *x, real *rhs, real tol, int maxit, int *flag){ @@ -176,6 +211,7 @@ static real conjugate_gradient(Operator A, Operator precon, int n, real *x, real rho_old = rho; } + FREE(z); FREE(r); FREE(p); FREE(q); #ifdef DEBUG _statistics[0] += iter - 1; #endif @@ -188,36 +224,41 @@ static real conjugate_gradient(Operator A, Operator precon, int n, real *x, real return res; } +real cg(Operator Ax, Operator precond, int n, int dim, real *x0, real *rhs, real tol, int maxit, int *flag){ + real *x, *b, res = 0; + int k, i; + x = N_GNEW(n, real); + b = N_GNEW(n, real); + for (k = 0; k < dim; k++){ + for (i = 0; i < n; i++) { + x[i] = x0[i*dim+k]; + b[i] = rhs[i*dim+k]; + } + + res += conjugate_gradient(Ax, precond, n, x, b, tol, maxit, flag); + for (i = 0; i < n; i++) { + rhs[i*dim+k] = x[i]; + } + } + FREE(x); + FREE(b); + return res; + +} + real SparseMatrix_solve(SparseMatrix A, int dim, real *x0, real *rhs, real tol, int maxit, int method, int *flag){ Operator Ax, precond; - real *x, *b, res = 0; - int n = A->m, k, i; - + int n = A->m; + real res = 0; *flag = 0; switch (method){ case SOLVE_METHOD_CG: Ax = Operator_matmul_new(A); precond = Operator_diag_precon_new(A); - - x = N_GNEW(n,real); - b = N_GNEW(n,real); - for (k = 0; k < dim; k++){ - for (i = 0; i < n; i++) { - x[i] = x0[i*dim+k]; - b[i] = rhs[i*dim+k]; - } - - res += conjugate_gradient(Ax, precond, n, x, b, tol, maxit, flag); - for (i = 0; i < n; i++) { - rhs[i*dim+k] = x[i]; - } - } + res = cg(Ax, precond, n, dim, x0, rhs, tol, maxit, flag); Operator_matmul_delete(Ax); Operator_diag_precon_delete(precond); - FREE(x); - FREE(b); - break; default: assert(0); diff --git a/lib/sfdpgen/sparse_solve.h b/lib/sfdpgen/sparse_solve.h index ede73f3ae..c0d21360a 100644 --- a/lib/sfdpgen/sparse_solve.h +++ b/lib/sfdpgen/sparse_solve.h @@ -27,7 +27,13 @@ struct Operator_struct { real* (*Operator_apply)(Operator o, real *in, real *out); }; +real cg(Operator Ax, Operator precond, int n, int dim, real *x0, real *rhs, real tol, int maxit, int *flag); + real SparseMatrix_solve(SparseMatrix A, int dim, real *x0, real *rhs, real tol, int maxit, int method, int *flag); +Operator Operator_uniform_stress_matmul(SparseMatrix A, real alpha); + +Operator Operator_uniform_stress_diag_precon_new(SparseMatrix A, real alpha); + #endif diff --git a/lib/sfdpgen/spring_electrical.c b/lib/sfdpgen/spring_electrical.c index dbbcc80a4..62a21a1b7 100644 --- a/lib/sfdpgen/spring_electrical.c +++ b/lib/sfdpgen/spring_electrical.c @@ -1,3 +1,4 @@ +/* $Id$Revision: */ /* vim:set shiftwidth=4 ts=8: */ /************************************************************************* @@ -10,8 +11,7 @@ * Contributors: See CVS logs. Details at http://www.graphviz.org/ *************************************************************************/ - - +#include "SparseMatrix.h" #include "spring_electrical.h" #include "QuadTree.h" #include "Multilevel.h" @@ -26,23 +26,24 @@ #include #include -#define drand() (rand()/(real) RAND_MAX) - -#define DEBUG_PRINT -#define MALLOC gmalloc -#define REALLOC grealloc -#define FREE free -#define MEMCPY memcpy - - spring_electrical_control spring_electrical_control_new(){ spring_electrical_control ctrl; - ctrl = N_GNEW(1,struct spring_electrical_control_struct); + ctrl = MALLOC(sizeof(struct spring_electrical_control_struct)); ctrl->p = AUTOP;/*a negativve number default to -1. repulsive force = dist^p */ + ctrl->q = 1;/*a positive number default to 1. Only apply to maxent. + attractive force = dist^q. Stress energy = (||x_i-x_j||-d_ij)^{q+1} */ ctrl->random_start = TRUE;/* whether to apply SE from a random layout, or from exisiting layout */ ctrl->K = -1;/* the natural distance. If K < 0, K will be set to the average distance of an edge */ 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. */ ctrl->multilevels = FALSE;/* if <=1, single level */ + + //ctrl->multilevel_coarsen_scheme = COARSEN_INDEPENDENT_EDGE_SET; + //ctrl->multilevel_coarsen_mode = COARSEN_MODE_GENTLE; + + ctrl->multilevel_coarsen_scheme = COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_EDGE_PERNODE_SUPERNODES_FIRST; /* pass on to Multilevel_control->coarsen_scheme */ + ctrl->multilevel_coarsen_mode = COARSEN_MODE_FORCEFUL;/*alternative: COARSEN_MODE_GENTLE. pass on to Multilevel_control->coarsen_mode */ + + ctrl->quadtree_size = 45;/* cut off size above which quadtree approximation is used */ ctrl->max_qtree_level = 10;/* max level of quadtree */ ctrl->bh = 0.6;/* Barnes-Hutt constant, if width(snode)/dist[i,snode] < bh, treat snode as a supernode.*/ @@ -56,10 +57,11 @@ spring_electrical_control spring_electrical_control_new(){ ctrl->use_node_weights = FALSE; ctrl->smoothing = SMOOTHING_NONE; ctrl->overlap = 0; - ctrl->tscheme = QUAD_TREE_NORMAL; - ctrl->posSet = NULL; + ctrl->tscheme = QUAD_TREE_HYBRID; + ctrl->method = METHOD_SPRING_ELECTRICAL; ctrl->initial_scaling = -4; - + ctrl->rotation = 0.; + ctrl->edge_labeling_scheme = 0; return ctrl; } @@ -67,6 +69,8 @@ void spring_electrical_control_delete(spring_electrical_control ctrl){ FREE(ctrl); } + + enum {MAX_I = 20, OPT_UP = 1, OPT_DOWN = -1, OPT_INIT = 0}; struct oned_optimizer_struct{ int i; @@ -83,7 +87,7 @@ static void oned_optimizer_delete(oned_optimizer opt){ static oned_optimizer oned_optimizer_new(int i){ oned_optimizer opt; - opt = N_GNEW(1,struct oned_optimizer_struct); + opt = MALLOC(sizeof(struct oned_optimizer_struct)); opt->i = i; opt->direction = OPT_INIT; return opt; @@ -150,49 +154,6 @@ real average_edge_length(SparseMatrix A, int dim, real *coord){ return dist/ia[A->m]; } -static real -initPos (SparseMatrix A, int n, int dim, real* x, spring_electrical_control ctrl) -{ - int i, j; - char* posSet; - - if (ctrl->random_start){ - srand(ctrl->random_seed); - if ((posSet = ctrl->posSet)) { - for (i = 0; i < n; i++) { - if (posSet[i]) continue; - for (j = 0; j < dim; j++) { - x[dim*i + j] = drand(); - } - } - } - else - for (i = 0; i < dim*n; i++) x[i] = drand(); - } - - if (ctrl->K < 0){ - ctrl->K = average_edge_length(A, dim, x); - } - - return ctrl->K; -} - -real distance_cropped(real *x, int dim, int i, int j){ - int k; - real dist = 0.; - for (k = 0; k < dim; k++) dist += (x[i*dim+k] - x[j*dim + k])*(x[i*dim+k] - x[j*dim + k]); - dist = sqrt(dist); - return MAX(dist, MINDIST); -} - -real distance(real *x, int dim, int i, int j){ - int k; - real dist = 0.; - for (k = 0; k < dim; k++) dist += (x[i*dim+k] - x[j*dim + k])*(x[i*dim+k] - x[j*dim + k]); - dist = sqrt(dist); - return dist; -} - #ifdef ENERGY static real spring_electrical_energy(int dim, SparseMatrix A, real *x, real p, real CRK, real KP){ /* 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) @@ -234,6 +195,20 @@ static real spring_electrical_energy(int dim, SparseMatrix A, real *x, real p, r void export_embedding(FILE *fp, int dim, SparseMatrix A, real *x, real *width){ int i, j, k, *ia=A->ia, *ja = A->ja; int ne = 0; + real xsize, ysize, xmin, xmax, ymin, ymax; + + xmax = xmin = x[0]; + ymax = ymin = x[1]; + for (i = 0; i < A->m; i++){ + xmax = MAX(xmax, x[i*dim]); + xmin = MIN(xmin, x[i*dim]); + ymax = MAX(ymax, x[i*dim+1]); + ymin = MIN(ymin, x[i*dim+1]); + } + xsize = xmax-xmin; + ysize = ymax-ymin; + xsize = MAX(xsize, ysize); + if (dim == 2){ fprintf(fp,"Graphics[{GrayLevel[0.5],Line[{"); } else { @@ -258,11 +233,19 @@ void export_embedding(FILE *fp, int dim, SparseMatrix A, real *x, real *width){ } } - fprintf(fp,"}],Hue[%f],",/*drand()*/1.); + fprintf(fp,"}],Hue[%f]",/*drand()*/1.); + if (width && dim == 2){ + for (i = 0; i < A->m; i++){ + if (i >= 0) fprintf(fp,","); + fprintf(fp,"(*%f,%f*){GrayLevel[.5,.5],Rectangle[{%f,%f},{%f,%f}]}", width[i*dim], width[i*dim+1],x[i*dim] - width[i*dim] + 5, x[i*dim+1] - width[i*dim+1] + 5, + x[i*dim] + width[i*dim] - 5, x[i*dim+1] + width[i*dim+1] - 5); + } + } + if (A->m < 100){ for (i = 0; i < A->m; i++){ - if (i > 0) fprintf(fp,","); + if (i >= 0) fprintf(fp,","); fprintf(fp,"Text[%d,{",i+1); for (k = 0; k < dim; k++) { if (k > 0) fprintf(fp,","); @@ -271,7 +254,7 @@ void export_embedding(FILE *fp, int dim, SparseMatrix A, real *x, real *width){ fprintf(fp,"}]"); } } else if (A->m < 500000){ - fprintf(fp, "Point[{"); + fprintf(fp, ", Point[{"); for (i = 0; i < A->m; i++){ if (i > 0) fprintf(fp,","); fprintf(fp,"{"); @@ -286,16 +269,8 @@ void export_embedding(FILE *fp, int dim, SparseMatrix A, real *x, real *width){ fprintf(fp,"{}"); } - if (width && dim == 2){ - fprintf(fp, ","); - for (i = 0; i < A->m; i++){ - if (i > 0) fprintf(fp,","); - fprintf(fp,"(*%f,%f*){GrayLevel[.5,.5],Rectangle[{%f,%f},{%f,%f}]}", width[i*dim], width[i*dim+1],x[i*dim] - width[i*dim], x[i*dim+1] - width[i*dim+1], - x[i*dim] + width[i*dim], x[i*dim+1] + width[i*dim+1]); - } - } - fprintf(fp,"},ImageSize->600]\n"); + fprintf(fp,"},ImageSize->%f]\n", 2*xsize/2); } @@ -320,14 +295,14 @@ static real update_step(int adaptive_cooling, real step, real Fnorm, real Fnorm0 void check_real_array_size(real **a, int len, int *lenmax){ if (len >= *lenmax){ *lenmax = len + MAX((int) 0.2*len, 10); - *a = RALLOC((*lenmax),*a,real); + *a = REALLOC(*a, sizeof(real)*(*lenmax)); } } void check_int_array_size(int **a, int len, int *lenmax){ if (len >= *lenmax){ *lenmax = len + MAX((int) 0.2*len, 10); - *a = RALLOC((*lenmax),*a,int); + *a = REALLOC(*a, sizeof(int)*(*lenmax)); } } @@ -341,14 +316,14 @@ real get_angle(real *x, int dim, int i, int j){ y[k] = x[j*dim+k] - x[i*dim+k]; } if (ABS(y[0]) <= ABS(y[1])*eps){ - if (y[1] > 0) return 0.5*M_PI; - return 1.5*M_PI; + if (y[1] > 0) return 0.5*PI; + return 1.5*PI; } res = atan(y[1]/y[0]); if (y[0] > 0){ - if (y[1] < 0) res = 2*M_PI+res; + if (y[1] < 0) res = 2*PI+res; } else if (y[0] < 0){ - res = res + M_PI; + res = res + PI; } return res; } @@ -384,9 +359,9 @@ static void beautify_leaves(int dim, SparseMatrix A, real *x){ assert(!SparseMatrix_has_diagonal(A)); - checked = N_GNEW(m,int); - angles = N_GNEW(nangles_max,real); - leaves = N_GNEW(nleaves_max,int); + checked = MALLOC(sizeof(int)*m); + angles = MALLOC(sizeof(real)*nangles_max); + leaves = MALLOC(sizeof(int)*nleaves_max); for (i = 0; i < m; i++) checked[i] = FALSE; @@ -421,15 +396,15 @@ static void beautify_leaves(int dim, SparseMatrix A, real *x){ ang1 = angles[k]; ang2 = angles[k+1]; } } - if (2*M_PI + angles[0] - angles[nangles - 1] > maxang){ - maxang = 2*M_PI + angles[0] - angles[nangles - 1]; + if (2*PI + angles[0] - angles[nangles - 1] > maxang){ + maxang = 2*PI + angles[0] - angles[nangles - 1]; ang1 = angles[nangles - 1]; - ang2 = 2*M_PI + angles[0]; + ang2 = 2*PI + angles[0]; } } else { - ang1 = 0; ang2 = 2*M_PI; maxang = 2*M_PI; + ang1 = 0; ang2 = 2*PI; maxang = 2*PI; } - pad = MAX(maxang - M_PI*0.166667*(nleaves-1), 0)*0.5; + pad = MAX(maxang - PI*0.166667*(nleaves-1), 0)*0.5; ang1 += pad*0.95; ang2 -= pad*0.95; assert(ang2 >= ang1); @@ -524,15 +499,20 @@ void spring_electrical_embedding_fast(int dim, SparseMatrix A0, spring_electrica ia = A->ia; ja = A->ja; - K = initPos (A, n, dim, x, ctrl); - + if (ctrl->random_start){ + srand(ctrl->random_seed); + for (i = 0; i < dim*n; i++) x[i] = drand(); + } + if (K < 0){ + ctrl->K = K = average_edge_length(A, dim, x); + } if (C < 0) ctrl->C = C = 0.2; if (p >= 0) ctrl->p = p = -1; KP = pow(K, 1 - p); CRK = pow(C, (2.-p)/3.)/K; - xold = N_GNEW(dim*n,real); - force = N_GNEW(dim*n,real); + xold = MALLOC(sizeof(real)*dim*n); + force = MALLOC(sizeof(real)*dim*n); do { #ifdef TIME @@ -564,7 +544,7 @@ void spring_electrical_embedding_fast(int dim, SparseMatrix A0, spring_electrica start = clock(); #endif - QuadTree_get_repulvice_force(qt, force, x, ctrl->bh, p, KP, counts, flag); + QuadTree_get_repulsive_force(qt, force, x, ctrl->bh, p, KP, counts, flag); assert(!(*flag)); @@ -612,11 +592,12 @@ void spring_electrical_embedding_fast(int dim, SparseMatrix A0, spring_electrica #ifdef TIME qtree_cpu0 = qtree_cpu - qtree_cpu0; qtree_new_cpu0 = qtree_new_cpu - qtree_new_cpu0; - 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", + /* 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", iter, ((real) (clock() - start2)) / CLOCKS_PER_SEC, qtree_new_cpu0, qtree_cpu0,((real) (clock() - start2))/CLOCKS_PER_SEC - qtree_cpu0 - qtree_new_cpu0, counts[0], counts[1], counts[2], step, Fnorm, A->nz,K,max_qtree_level); + */ qtree_cpu0 = qtree_cpu; qtree_new_cpu0 = qtree_new_cpu; #endif @@ -696,14 +677,14 @@ void spring_electrical_embedding_slow(int dim, SparseMatrix A0, spring_electrica m = A->m, n = A->n; if (n <= 0 || dim <= 0) return; - force = N_GNEW(n*dim,real); + force = MALLOC(sizeof(real)*n*dim); if (n >= ctrl->quadtree_size) { USE_QT = TRUE; qtree_level_optimizer = oned_optimizer_new(max_qtree_level); - center = N_GNEW(nsupermax*dim,real); - supernode_wgts = N_GNEW(nsupermax,real); - distances = N_GNEW(nsupermax,real); + center = MALLOC(sizeof(real)*nsupermax*dim); + supernode_wgts = MALLOC(sizeof(real)*nsupermax); + distances = MALLOC(sizeof(real)*nsupermax); } USE_QT = FALSE; *flag = 0; @@ -716,8 +697,13 @@ void spring_electrical_embedding_slow(int dim, SparseMatrix A0, spring_electrica ia = A->ia; ja = A->ja; - K = initPos (A, n, dim, x, ctrl); - + if (ctrl->random_start){ + srand(ctrl->random_seed); + for (i = 0; i < dim*n; i++) x[i] = drand(); + } + if (K < 0){ + ctrl->K = K = average_edge_length(A, dim, x); + } if (C < 0) ctrl->C = C = 0.2; if (p >= 0) ctrl->p = p = -1; KP = pow(K, 1 - p); @@ -735,8 +721,8 @@ void spring_electrical_embedding_slow(int dim, SparseMatrix A0, spring_electrica } #endif - f = N_GNEW(dim,real); - xold = N_GNEW(dim*n,real); + f = MALLOC(sizeof(real)*dim); + xold = MALLOC(sizeof(real)*dim*n); do { for (i = 0; i < dim*n; i++) force[i] = 0; @@ -954,9 +940,9 @@ void spring_electrical_embedding(int dim, SparseMatrix A0, spring_electrical_con if (n >= ctrl->quadtree_size) { USE_QT = TRUE; qtree_level_optimizer = oned_optimizer_new(max_qtree_level); - center = N_GNEW(nsupermax*dim,real); - supernode_wgts = N_GNEW(nsupermax,real); - distances = N_GNEW(nsupermax,real); + center = MALLOC(sizeof(real)*nsupermax*dim); + supernode_wgts = MALLOC(sizeof(real)*nsupermax); + distances = MALLOC(sizeof(real)*nsupermax); } *flag = 0; if (m != n) { @@ -968,8 +954,13 @@ void spring_electrical_embedding(int dim, SparseMatrix A0, spring_electrical_con ia = A->ia; ja = A->ja; - K = initPos (A, n, dim, x, ctrl); - + if (ctrl->random_start){ + srand(ctrl->random_seed); + for (i = 0; i < dim*n; i++) x[i] = drand(); + } + if (K < 0){ + ctrl->K = K = average_edge_length(A, dim, x); + } if (C < 0) ctrl->C = C = 0.2; if (p >= 0) ctrl->p = p = -1; KP = pow(K, 1 - p); @@ -987,8 +978,8 @@ void spring_electrical_embedding(int dim, SparseMatrix A0, spring_electrical_con } #endif - f = N_GNEW(dim,real); - xold = N_GNEW(dim*n,real); + f = MALLOC(sizeof(real)*dim); + xold = MALLOC(sizeof(real)*dim*n); do { iter++; xold = MEMCPY(xold, x, sizeof(real)*dim*n); @@ -1160,6 +1151,320 @@ void spring_electrical_embedding(int dim, SparseMatrix A0, spring_electrical_con } +static void scale_coord(int n, int dim, real *x, int *id, int *jd, real *d, real dj){ + int i, j, k; + real w_ij, dist, s = 0, stop = 0, sbot = 0., nz = 0; + + if (dj == 0.) return; + for (i = 0; i < n; i++){ + for (j = id[i]; j < id[i+1]; j++){ + if (jd[j] == i) continue; + dist = distance_cropped(x, dim, i, jd[j]); + if (d){ + dj = d[j]; + } + assert(dj > 0); + w_ij = 1./(dj*dj); + /* spring force */ + for (k = 0; k < dim; k++){ + stop += w_ij*dj*dist; + sbot += w_ij*dist*dist; + } + s += dist; nz++; + } + } + s = stop/sbot; + for (i = 0; i < n*dim; i++) x[i] *= s; + fprintf(stderr,"scaling factor = %f\n",s); +} + +static real dmean_get(int n, int *id, int *jd, real* d){ + real dmean = 0; + int i, j; + + if (!d) return 1.; + for (i = 0; i < n; i++){ + for (j = id[i]; j < id[i+1]; j++){ + dmean += d[j]; + } + } + return dmean/((real) id[n]); +} + +void spring_maxent_embedding(int dim, SparseMatrix A0, SparseMatrix D, spring_electrical_control ctrl, real *node_weights, real *x, real rho, int *flag){ + /* x is a point to a 1D array, x[i*dim+j] gives the coordinate of the i-th node at dimension j. + + 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|| + + or + + 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 + + The derivatives are + + 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 + + or + + 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) + + if D == NULL, unit weight assumed + + */ + SparseMatrix A = A0; + int m, n; + int i, j, k; + real p = ctrl->p, C = ctrl->C, tol = ctrl->tol, maxiter = ctrl->maxiter, cool = ctrl->cool, step = ctrl->step, w_ij, dj = 1.; + int *ia = NULL, *ja = NULL; + int *id = NULL, *jd = NULL; + real *d, dmean; + real *xold = NULL; + real *f = NULL, dist, F, Fnorm = 0, Fnorm0; + int iter = 0; + int adaptive_cooling = ctrl->adaptive_cooling; + QuadTree qt = NULL; + int USE_QT = FALSE; + int nsuper = 0, nsupermax = 10; + real *center = NULL, *supernode_wgts = NULL, *distances = NULL, nsuper_avg, counts = 0; + int max_qtree_level = 10; +#ifdef DEBUG + double stress = 0; +#endif + + if (!A) return; + m = A->m, n = A->n; + if (n <= 0 || dim <= 0) return; + + if (ctrl->tscheme != QUAD_TREE_NONE && n >= ctrl->quadtree_size) { + USE_QT = TRUE; + center = MALLOC(sizeof(real)*nsupermax*dim); + supernode_wgts = MALLOC(sizeof(real)*nsupermax); + distances = MALLOC(sizeof(real)*nsupermax); + } + + *flag = 0; + if (m != n) { + *flag = ERROR_NOT_SQUARE_MATRIX; + goto RETURN; + } + + + assert(A->format == FORMAT_CSR); + A = SparseMatrix_symmetrize(A, TRUE); + ia = A->ia; + ja = A->ja; + if (D){ + id = D->ia; + jd = D->ja; + d = (real*) D->a; + } else { + id = ia; jd = ja; d = NULL; + } + if (rho < 0) { + dmean = dmean_get(n, id, jd, d); + rho = rho*(id[n]/((((real) n)*((real) n)) - id[n]))/pow(dmean, p+1); + fprintf(stderr,"dmean = %f, rho = %f\n",dmean, rho); + } + + if (ctrl->random_start){ + fprintf(stderr, "send random coordinates\n"); + srand(ctrl->random_seed); + for (i = 0; i < dim*n; i++) x[i] = drand(); + /* rescale x to give minimum stress: + Min \Sum_{(i,j)\in E} w_ij (s ||x_i-x_j||-d_ij)^2 + thus + 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) + */ + + } + scale_coord(n, dim, x, id, jd, d, dj); + + + + if (C < 0) ctrl->C = C = 0.2; + if (p >= 0) ctrl->p = p = -1; + +#ifdef DEBUG_0 + { + FILE *f; + char fname[10000]; + strcpy(fname,"/tmp/graph_layout_0_"); + sprintf(&(fname[strlen(fname)]), "%d",n); + f = fopen(fname,"w"); + export_embedding(f, dim, A, x, NULL); + fclose(f); + } +#endif + + f = MALLOC(sizeof(real)*dim); + xold = MALLOC(sizeof(real)*dim*n); + do { + iter++; + xold = MEMCPY(xold, x, sizeof(real)*dim*n); + Fnorm0 = Fnorm; + Fnorm = 0.; + nsuper_avg =- 0; +#ifdef DEBUG + stress = 0; +#endif + + if (USE_QT) { + if (ctrl->use_node_weights){ + qt = QuadTree_new_from_point_list(dim, n, max_qtree_level, x, node_weights); + } else { + qt = QuadTree_new_from_point_list(dim, n, max_qtree_level, x, NULL); + } + } + + /* + . 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 + or + . 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) + */ + for (i = 0; i < n; i++){ + for (k = 0; k < dim; k++) f[k] = 0.; + + /* spring (attractive or repulsive) force w_ij (||x_i-x_j||-d_ij) (x_i-x_j)/||x_i-x_j|| */ + for (j = id[i]; j < id[i+1]; j++){ + if (jd[j] == i) continue; + dist = distance_cropped(x, dim, i, jd[j]); + if (d){ + dj = d[j]; + } + assert(dj > 0); + /* spring force */ + if (ctrl->q == 2){ + w_ij = 1./(dj*dj*dj); + for (k = 0; k < dim; k++){ + f[k] += -w_ij*(x[i*dim+k] - x[jd[j]*dim+k])*(dist - dj)*(dist - dj)/dist; + } + } else if (ctrl->q == 1){/* square stress force */ + w_ij = 1./(dj*dj); + for (k = 0; k < dim; k++){ + f[k] += -w_ij*(x[i*dim+k] - x[jd[j]*dim+k])*(dist - dj)/dist; + } + } else { + w_ij = 1./pow(dj, ctrl->q + 1); + for (k = 0; k < dim; k++){ + f[k] += -w_ij*(x[i*dim+k] - x[jd[j]*dim+k])*pow(dist - dj, ctrl->q)/dist; + } + } + +#ifdef DEBUG + w_ij = 1./(dj*dj); + for (k = 0; k < dim; k++){ + stress += (dist - dj)*(dist - dj)*w_ij; + } +#endif + + + /* discount repulsive force between neighboring vertices which will be applied next, that way there is no + repulsive forces between neighboring vertices */ + if (ctrl->use_node_weights && node_weights){ + for (k = 0; k < dim; k++){ + if (p == -1){ + f[k] -= rho*node_weights[j]*(x[i*dim+k] - x[jd[j]*dim+k])/(dist*dist); + } else { + f[k] -= rho*node_weights[j]*(x[i*dim+k] - x[jd[j]*dim+k])/pow(dist, 1.- p); + } + } + } else { + for (k = 0; k < dim; k++){ + if (p == -1){ + f[k] -= rho*(x[i*dim+k] - x[jd[j]*dim+k])/(dist*dist); + } else { + f[k] -= rho*(x[i*dim+k] - x[jd[j]*dim+k])/pow(dist, 1.- p); + } + } + + } + + } + + /* repulsive force ||x_i-x_j||^(1 - p) (x_i - x_j) */ + if (USE_QT){ + QuadTree_get_supernodes(qt, ctrl->bh, &(x[dim*i]), i, &nsuper, &nsupermax, + ¢er, &supernode_wgts, &distances, &counts, flag); + nsuper_avg += nsuper; + if (*flag) goto RETURN; + for (j = 0; j < nsuper; j++){ + dist = MAX(distances[j], MINDIST); + for (k = 0; k < dim; k++){ + if (p == -1){ + f[k] += rho*supernode_wgts[j]*(x[i*dim+k] - center[j*dim+k])/(dist*dist); + } else { + f[k] += rho*supernode_wgts[j]*(x[i*dim+k] - center[j*dim+k])/pow(dist, 1.- p); + } + } + } + } else { + if (ctrl->use_node_weights && node_weights){ + for (j = 0; j < n; j++){ + if (j == i) continue; + dist = distance_cropped(x, dim, i, j); + for (k = 0; k < dim; k++){ + if (p == -1){ + f[k] += rho*node_weights[j]*(x[i*dim+k] - x[j*dim+k])/(dist*dist); + } else { + f[k] += rho*node_weights[j]*(x[i*dim+k] - x[j*dim+k])/pow(dist, 1.- p); + } + } + } + } else { + for (j = 0; j < n; j++){ + if (j == i) continue; + dist = distance_cropped(x, dim, i, j); + for (k = 0; k < dim; k++){ + if (p == -1){ + f[k] += rho*(x[i*dim+k] - x[j*dim+k])/(dist*dist); + } else { + f[k] += rho*(x[i*dim+k] - x[j*dim+k])/pow(dist, 1.- p); + } + } + } + } + } + + /* normalize force */ + F = 0.; + for (k = 0; k < dim; k++) F += f[k]*f[k]; + F = sqrt(F); + Fnorm += F; + + if (F > 0) for (k = 0; k < dim; k++) f[k] /= F; + + for (k = 0; k < dim; k++) x[i*dim+k] += step*f[k]; + + }/* done vertex i */ + + if (qt) QuadTree_delete(qt); + nsuper_avg /= n; +#ifdef DEBUG_PRINT + stress /= (double) A->nz; + if (Verbose) { + fprintf(stderr, "\r iter = %d, step = %f Fnorm = %f nsuper = %d nz = %d stress = %f ",iter, step, Fnorm, (int) nsuper_avg,A->nz, stress); + } +#endif + + step = update_step(adaptive_cooling, step, Fnorm, Fnorm0, cool); + } while (step > tol && iter < maxiter); + +#ifdef DEBUG_PRINT + if (Verbose) fputs("\n", stderr); +#endif + + + if (ctrl->beautify_leaves) beautify_leaves(dim, A, x); + + RETURN: + if (xold) FREE(xold); + if (A != A0) SparseMatrix_delete(A); + if (f) FREE(f); + if (center) FREE(center); + if (supernode_wgts) FREE(supernode_wgts); + if (distances) FREE(distances); + +} @@ -1190,9 +1495,9 @@ void spring_electrical_spring_embedding(int dim, SparseMatrix A0, SparseMatrix D if (n >= ctrl->quadtree_size) { USE_QT = TRUE; - center = N_GNEW(nsupermax*dim,real); - supernode_wgts = N_GNEW(nsupermax,real); - distances = N_GNEW(nsupermax,real); + center = MALLOC(sizeof(real)*nsupermax*dim); + supernode_wgts = MALLOC(sizeof(real)*nsupermax); + distances = MALLOC(sizeof(real)*nsupermax); } *flag = 0; if (m != n) { @@ -1207,8 +1512,13 @@ void spring_electrical_spring_embedding(int dim, SparseMatrix A0, SparseMatrix D jd = D->ja; d = (real*) D->a; - K = initPos (A, n, dim, x, ctrl); - + if (ctrl->random_start){ + srand(ctrl->random_seed); + for (i = 0; i < dim*n; i++) x[i] = drand(); + } + if (K < 0){ + ctrl->K = K = average_edge_length(A, dim, x); + } if (C < 0) ctrl->C = C = 0.2; if (p >= 0) ctrl->p = p = -1; KP = pow(K, 1 - p); @@ -1226,8 +1536,8 @@ void spring_electrical_spring_embedding(int dim, SparseMatrix A0, SparseMatrix D } #endif - f = N_GNEW(dim,real); - xold = N_GNEW(dim*n,real); + f = MALLOC(sizeof(real)*dim); + xold = MALLOC(sizeof(real)*dim*n); do { iter++; xold = MEMCPY(xold, x, sizeof(real)*dim*n); @@ -1389,7 +1699,7 @@ static void interpolate2(int dim, SparseMatrix A, real *x){ int i, j, k, *ia = A->ia, *ja = A->ja, nz, m = A->m; real alpha = 0.5, beta, *y; - y = N_GNEW(dim*m,real); + y = MALLOC(sizeof(real)*dim*m); for (k = 0; k < dim*m; k++) y[k] = 0; for (i = 0; i < m; i++){ nz = 0; @@ -1412,11 +1722,11 @@ static void interpolate2(int dim, SparseMatrix A, real *x){ */ -static void interpolate0(int dim, SparseMatrix A, real *x){ +static void interpolate(int dim, SparseMatrix A, real *x){ int i, j, k, *ia = A->ia, *ja = A->ja, nz; real alpha = 0.5, beta, *y; - y = N_GNEW(dim,real); + y = MALLOC(sizeof(real)*dim); for (i = 0; i < A->m; i++){ for (k = 0; k < dim; k++) y[k] = 0; nz = 0; @@ -1441,7 +1751,7 @@ static void prolongate(int dim, SparseMatrix A, SparseMatrix P, SparseMatrix R, /* xu yao rao dong */ if (coarsen_scheme_used > EDGE_BASED_STA && coarsen_scheme_used < EDGE_BASED_STO){ - interpolate0(dim, A, y); + interpolate(dim, A, y); nc = R->m; ia = R->ia; ja = R->ja; @@ -1461,7 +1771,7 @@ int power_law_graph(SparseMatrix A){ int *mask, m, max = 0, i, *ia = A->ia, *ja = A->ja, j, deg; int res = FALSE; m = A->m; - mask = N_GNEW((m+1),int); + mask = MALLOC(sizeof(int)*(m+1)); for (i = 0; i < m + 1; i++){ mask[i] = 0; @@ -1532,21 +1842,170 @@ void pcp_rotate(int n, int dim, real *x){ } + +static void rotate(int n, int dim, real *x, real angle){ + int i, k; + real y[4], axis[2], center[2], x0, x1; + real radian = 3.14159/180; + assert(dim == 2); + for (i = 0; i < dim*dim; i++) y[i] = 0; + for (i = 0; i < dim; i++) center[i] = 0; + for (i = 0; i < n; i++){ + for (k = 0; k < dim; k++){ + center[k] += x[i*dim+k]; + } + } + for (i = 0; i < dim; i++) center[i] /= n; + for (i = 0; i < n; i++){ + for (k = 0; k < dim; k++){ + x[dim*i+k] = x[dim*i+k] - center[k]; + } + } + axis[0] = cos(-angle*radian); + axis[1] = sin(-angle*radian); + for (i = 0; i < n; i++){ + x0 = x[dim*i]*axis[0]+x[dim*i+1]*axis[1]; + x1 = -x[dim*i]*axis[1]+x[dim*i+1]*axis[0]; + x[dim*i] = x0; + x[dim*i + 1] = x1; + } + + +} + +static void attach_edge_label_coordinates(int dim, SparseMatrix A, int n_edge_label_nodes, int *edge_label_nodes, real *x, real *x2){ + int *mask; + int i, ii, j, k; + int nnodes = 0; + real len; + + mask = MALLOC(sizeof(int)*A->m); + + for (i = 0; i < A->m; i++) mask[i] = 1; + for (i = 0; i < n_edge_label_nodes; i++) { + if (edge_label_nodes[i] >= 0 && edge_label_nodes[i] < A->m) mask[edge_label_nodes[i]] = -1; + } + + for (i = 0; i < A->m; i++) { + if (mask[i] >= 0) mask[i] = nnodes++; + } + + + for (i = 0; i < A->m; i++){ + if (mask[i] >= 0){ + for (k = 0; k < dim; k++) x[i*dim+k] = x2[mask[i]*dim+k]; + } + } + + for (i = 0; i < n_edge_label_nodes; i++){ + ii = edge_label_nodes[i]; + len = A->ia[ii+1] - A->ia[ii]; + assert(len >= 2); /* should just be 2 */ + assert(mask[ii] < 0); + for (k = 0; k < dim; k++) { + x[ii*dim+k] = 0; + } + for (j = A->ia[ii]; j < A->ia[ii+1]; j++){ + for (k = 0; k < dim; k++){ + x[ii*dim+k] += x[(A->ja[j])*dim+k]; + } + } + for (k = 0; k < dim; k++) { + x[ii*dim+k] /= len; + } + } + + FREE(mask); +} + +static SparseMatrix shorting_edge_label_nodes(SparseMatrix A, int n_edge_label_nodes, int *edge_label_nodes){ + int *mask; + int i, id = 0, nz, j, jj, ii; + int *ia = A->ia, *ja = A->ja, *irn = NULL, *jcn = NULL; + SparseMatrix B; + + mask = MALLOC(sizeof(int)*A->m); + + for (i = 0; i < A->m; i++) mask[i] = 1; + + for (i = 0; i < n_edge_label_nodes; i++){ + mask[edge_label_nodes[i]] = -1; + } -void multilevel_spring_electrical_embedding(int dim, SparseMatrix A0, spring_electrical_control ctrl, real *node_weights, real *label_sizes, - real *x, int *flag){ + for (i = 0; i < A->m; i++) { + if (mask[i] > 0) mask[i] = id++; + } + + nz = 0; + for (i = 0; i < A->m; i++){ + if (mask[i] < 0) continue; + for (j = ia[i]; j < ia[i+1]; j++){ + if (mask[ja[j]] >= 0) { + nz++; + continue; + } + ii = ja[j]; + for (jj = ia[ii]; jj < ia[ii+1]; jj++){ + if (ja[jj] != i && mask[ja[jj]] >= 0) nz++; + } + } + } + + if (nz > 0) { + irn = MALLOC(sizeof(int)*nz); + jcn = MALLOC(sizeof(int)*nz); + } + + nz = 0; + for (i = 0; i < A->m; i++){ + if (mask[i] < 0) continue; + for (j = ia[i]; j < ia[i+1]; j++){ + if (mask[ja[j]] >= 0) { + irn[nz] = mask[i]; + jcn[nz++] = mask[ja[j]]; + if (mask[i] == 68 || mask[ja[j]] == 68){ + fprintf(stderr, "xxx %d %d\n",mask[i], mask[ja[j]]); + mask[i] = mask[i]; + } + continue; + } + ii = ja[j]; + for (jj = ia[ii]; jj < ia[ii+1]; jj++){ + if (ja[jj] != i && mask[ja[jj]] >= 0) { + irn[nz] = mask[i]; + jcn[nz++] = mask[ja[jj]]; + if (mask[i] == 68 || mask[ja[jj]] == 68){ + fprintf(stderr, "%d %d\n",mask[i], mask[ja[jj]]); + mask[i] = mask[i]; + } + } + } + } + } + + B = SparseMatrix_from_coordinate_arrays(nz, id, id, irn, jcn, NULL, MATRIX_TYPE_PATTERN); + + FREE(irn); + FREE(jcn); + FREE(mask); + return B; + +} + +void multilevel_spring_electrical_embedding(int dim, SparseMatrix A0, SparseMatrix D0, spring_electrical_control ctrl, real *node_weights, real *label_sizes, + real *x, int n_edge_label_nodes, int *edge_label_nodes, int *flag){ Multilevel_control mctrl = NULL; int n, plg, coarsen_scheme_used; - SparseMatrix A = A0, P = NULL; + SparseMatrix A = A0, D = D0, P = NULL; Multilevel grid, grid0; real *xc = NULL, *xf = NULL; + struct spring_electrical_control_struct ctrl0; #ifdef TIME clock_t cpu; #endif - struct spring_electrical_control_struct ctrl0; ctrl0 = *ctrl; @@ -1560,20 +2019,50 @@ void multilevel_spring_electrical_embedding(int dim, SparseMatrix A0, spring_ele if (n <= 0 || dim <= 0) return; if (!SparseMatrix_is_symmetric(A, FALSE) || A->type != MATRIX_TYPE_REAL){ - A = SparseMatrix_get_real_adjacency_matrix_symmetrized(A); + if (ctrl->method == METHOD_SPRING_MAXENT){ + A = SparseMatrix_symmetrize_nodiag(A, FALSE); + assert(D0); + D = SparseMatrix_symmetrize_nodiag(D, FALSE); + } else { + A = SparseMatrix_get_real_adjacency_matrix_symmetrized(A); + } } else { + if (ctrl->method == METHOD_SPRING_MAXENT){ + assert(D0); + D = SparseMatrix_remove_diagonal(D); + } A = SparseMatrix_remove_diagonal(A); } - mctrl = Multilevel_control_new(); + /* we first generate a layout discarding (shorting) the edge labels nodes, then assign the edge label nodes at the average of their neighbors */ + if ((ctrl->edge_labeling_scheme == ELSCHEME_STRAIGHTLINE_PENALTY || ctrl->edge_labeling_scheme == ELSCHEME_STRAIGHTLINE_PENALTY2) + && n_edge_label_nodes > 0){ + SparseMatrix A2; + + real *x2 = MALLOC(sizeof(real)*(A->m)*dim); + A2 = shorting_edge_label_nodes(A, n_edge_label_nodes, edge_label_nodes); + multilevel_spring_electrical_embedding(dim, A2, NULL, ctrl, NULL, NULL, x2, 0, NULL, flag); + + assert(!(*flag)); + attach_edge_label_coordinates(dim, A, n_edge_label_nodes, edge_label_nodes, x, x2); + remove_overlap(dim, A, x, label_sizes, ctrl->overlap, ctrl->initial_scaling, + ctrl->edge_labeling_scheme, n_edge_label_nodes, edge_label_nodes, A, flag); + SparseMatrix_delete(A2); + FREE(x2); + if (A != A0) SparseMatrix_delete(A); + + return; + } + + mctrl = Multilevel_control_new(ctrl->multilevel_coarsen_scheme, ctrl->multilevel_coarsen_mode); mctrl->maxlevel = ctrl->multilevels; - grid0 = Multilevel_new(A, node_weights, mctrl); + grid0 = Multilevel_new(A, D, node_weights, mctrl); grid = Multilevel_get_coarsest(grid0); if (Multilevel_is_finest(grid)){ xc = x; } else { - xc = N_GNEW(grid->n*dim,real); + xc = MALLOC(sizeof(real)*grid->n*dim); } plg = power_law_graph(A); @@ -1593,12 +2082,42 @@ void multilevel_spring_electrical_embedding(int dim, SparseMatrix A0, spring_ele } } #endif - if (ctrl->tscheme == QUAD_TREE_NONE){ - spring_electrical_embedding_slow(dim, grid->A, ctrl, grid->node_weights, xc, flag); - } else if (ctrl->tscheme == QUAD_TREE_NORMAL){ - spring_electrical_embedding(dim, grid->A, ctrl, grid->node_weights, xc, flag); - } else if (ctrl->tscheme == QUAD_TREE_FAST){ - spring_electrical_embedding_fast(dim, grid->A, ctrl, grid->node_weights, xc, flag); + if (ctrl->method == METHOD_SPRING_ELECTRICAL){ + if (ctrl->tscheme == QUAD_TREE_NONE){ + spring_electrical_embedding_slow(dim, grid->A, ctrl, grid->node_weights, xc, flag); + } else if (ctrl->tscheme == QUAD_TREE_FAST || (ctrl->tscheme == QUAD_TREE_HYBRID && grid->A->m > QUAD_TREE_HYBRID_SIZE)){ + if (ctrl->tscheme == QUAD_TREE_HYBRID && grid->A->m > 10 && Verbose){ + fprintf(stderr, "QUAD_TREE_HYBRID, size larger than %d, switch to fast quadtree", QUAD_TREE_HYBRID_SIZE); + } + spring_electrical_embedding_fast(dim, grid->A, ctrl, grid->node_weights, xc, flag); + } else if (ctrl->tscheme == QUAD_TREE_NORMAL){ + spring_electrical_embedding(dim, grid->A, ctrl, grid->node_weights, xc, flag); + } else { + spring_electrical_embedding(dim, grid->A, ctrl, grid->node_weights, xc, flag); + } + } else if (ctrl->method == METHOD_SPRING_MAXENT){ + double rho = 0.05; + + ctrl->step = 1; + ctrl->adaptive_cooling = TRUE; + if (Multilevel_is_coarsest(grid)){ + ctrl->maxiter=500; + rho = 0.5; + } else { + ctrl->maxiter=100; + } + + if (Multilevel_is_finest(grid)) {/* gradually reduce influence of entropy */ + spring_maxent_embedding(dim, grid->A, grid->D, ctrl, grid->node_weights, xc, rho, flag); + ctrl->random_start = FALSE; + ctrl->step = .05; + ctrl->adaptive_cooling = FALSE; + spring_maxent_embedding(dim, grid->A, grid->D, ctrl, grid->node_weights, xc, rho/2, flag); + spring_maxent_embedding(dim, grid->A, grid->D, ctrl, grid->node_weights, xc, rho/8, flag); + spring_maxent_embedding(dim, grid->A, grid->D, ctrl, grid->node_weights, xc, rho/32, flag); + } else { + spring_maxent_embedding(dim, grid->A, grid->D, ctrl, grid->node_weights, xc, rho, flag); + } } else { assert(0); } @@ -1613,7 +2132,7 @@ void multilevel_spring_electrical_embedding(int dim, SparseMatrix A0, spring_ele if (Multilevel_is_finest(grid)){ xf = x; } else { - xf = N_GNEW(grid->n*dim,real); + xf = MALLOC(sizeof(real)*grid->n*dim); } prolongate(dim, grid->A, P, grid->R, xc, xf, coarsen_scheme_used, (ctrl->K)*0.001); FREE(xc); @@ -1640,17 +2159,17 @@ void multilevel_spring_electrical_embedding(int dim, SparseMatrix A0, spring_ele if (dim == 2){ pcp_rotate(n, dim, x); } + if (ctrl->rotation != 0) rotate(n, dim, x, ctrl->rotation); - if (Verbose) - fprintf(stderr, "sfdp: overlap=%d scaling %.02f\n", - ctrl->overlap, ctrl->initial_scaling); + if (Verbose) fprintf(stderr, "ctrl->overlap=%d\n",ctrl->overlap); - if (ctrl->overlap >= 0) - remove_overlap(dim, A, A->m, x, label_sizes, ctrl->overlap, ctrl->initial_scaling, flag); + remove_overlap(dim, A, x, label_sizes, ctrl->overlap, ctrl->initial_scaling, + ctrl->edge_labeling_scheme, n_edge_label_nodes, edge_label_nodes, A, flag); RETURN: *ctrl = ctrl0; if (A != A0) SparseMatrix_delete(A); + if (D && D != D0) SparseMatrix_delete(D); Multilevel_control_delete(mctrl); Multilevel_delete(grid0); } diff --git a/lib/sfdpgen/spring_electrical.h b/lib/sfdpgen/spring_electrical.h index 62ba3ab9e..68289b805 100644 --- a/lib/sfdpgen/spring_electrical.h +++ b/lib/sfdpgen/spring_electrical.h @@ -21,18 +21,23 @@ enum {ERROR_NOT_SQUARE_MATRIX = -100}; /* a flag to indicate that p should be set auto */ #define AUTOP -1.0001234 -#define MINDIST 1.e-15 - enum {SMOOTHING_NONE, SMOOTHING_STRESS_MAJORIZATION_GRAPH_DIST, SMOOTHING_STRESS_MAJORIZATION_AVG_DIST, SMOOTHING_STRESS_MAJORIZATION_POWER_DIST, SMOOTHING_SPRING, SMOOTHING_TRIANGLE, SMOOTHING_RNG}; -enum {QUAD_TREE_NONE = 0, QUAD_TREE_NORMAL, QUAD_TREE_FAST}; +enum {QUAD_TREE_HYBRID_SIZE = 10000}; + +enum {QUAD_TREE_NONE = 0, QUAD_TREE_NORMAL, QUAD_TREE_FAST, QUAD_TREE_HYBRID}; + +enum {METHOD_SPRING_ELECTRICAL = 0, METHOD_SPRING_MAXENT, METHOD_STRESS, METHOD_UNIFORM_STRESS}; struct spring_electrical_control_struct { - real p;/*a negativve number default to -1. repulsive force = dist^p */ + real p;/*a negativve real number default to -1. repulsive force = dist^p */ + real q;/*a positive real number default to 2. attractive force = dist^q */ int random_start;/* whether to apply SE from a random layout, or from exisiting layout */ real K;/* the natural distance. If K < 0, K will be set to the average distance of an edge */ real C;/* 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. */ int multilevels;/* if <=1, single level */ + int multilevel_coarsen_scheme;/* pass on to Multilevel_control->coarsen_scheme */ + int multilevel_coarsen_mode;/* pass on to Multilevel_control->coarsen_mode */ int quadtree_size;/* cut off size above which quadtree approximation is used */ int max_qtree_level;/* max level of quadtree */ real bh;/* Barnes-Hutt constant, if width(snode)/dist[i,snode] < bh, treat snode as a supernode. default 0.2*/ @@ -47,11 +52,16 @@ struct spring_electrical_control_struct { int smoothing; int overlap; int tscheme; /* octree scheme. 0 (no octree), 1 (normal), 2 (fast) */ - char* posSet; + int method;/* spring_electical, spring_maxent */ real initial_scaling;/* how to scale the layout of the graph before passing to overlap removal algorithm. positive values are absolute in points, negative values are relative to average label size. */ + real rotation;/* degree of rotation */ + int edge_labeling_scheme; /* specifying whether to treat node of the form |edgelabel|* as a special node representing an edge label. + 0 (no action, default), 1 (penalty based method to make that kind of node close to the center of its neighbor), + 1 (penalty based method to make that kind of node close to the old center of its neighbor), + 3 (two step process of overlap removal and straightening) */ }; typedef struct spring_electrical_control_struct *spring_electrical_control; @@ -61,14 +71,13 @@ spring_electrical_control spring_electrical_control_new(void); void spring_electrical_embedding(int dim, SparseMatrix A0, spring_electrical_control ctrl, real *node_weights, real *x, int *flag); void spring_electrical_embedding_fast(int dim, SparseMatrix A0, spring_electrical_control ctrl, real *node_weights, real *x, int *flag); -void multilevel_spring_electrical_embedding(int dim, SparseMatrix A0, spring_electrical_control ctrl0, real *node_weights, real *label_sizes, real *x, int *flag); +void multilevel_spring_electrical_embedding(int dim, SparseMatrix A0, 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); void export_embedding(FILE *fp, int dim, SparseMatrix A, real *x, real *width); void spring_electrical_control_delete(spring_electrical_control ctrl); void print_matrix(real *x, int n, int dim); -real distance(real *x, int dim, int i, int j); -real distance_cropped(real *x, int dim, int i, int j); real average_edge_length(SparseMatrix A, int dim, real *coord); void spring_electrical_spring_embedding(int dim, SparseMatrix A, SparseMatrix D, spring_electrical_control ctrl, real *node_weights, real *x, int *flag); diff --git a/lib/sfdpgen/stress_model.c b/lib/sfdpgen/stress_model.c new file mode 100644 index 000000000..f3bec14b3 --- /dev/null +++ b/lib/sfdpgen/stress_model.c @@ -0,0 +1,49 @@ +#include "general.h" +#include "SparseMatrix.h" +#include "spring_electrical.h" +#include "post_process.h" +#include "stress_model.h" + +void stress_model(int dim, SparseMatrix B, real **x, int maxit_sm, real tol, int *flag){ + SparseStressMajorizationSmoother sm; + real lambda = 0; + /*int maxit_sm = 1000, i; tol = 0.001*/ + int i; + SparseMatrix A = B; + + if (!SparseMatrix_is_symmetric(A, FALSE) || A->type != MATRIX_TYPE_REAL){ + if (A->type == MATRIX_TYPE_REAL){ + A = SparseMatrix_symmetrize(A, FALSE); + A = SparseMatrix_remove_diagonal(A); + } else { + A = SparseMatrix_get_real_adjacency_matrix_symmetrized(A); + } + } + A = SparseMatrix_remove_diagonal(A); + + *flag = 0; + int m = A->m; + if (!x) { + *x = MALLOC(sizeof(real)*m*dim); + srand(123); + for (i = 0; i < dim*m; i++) (*x)[i] = drand(); + } + + sm = SparseStressMajorizationSmoother_new(A, dim, lambda, *x, WEIGHTING_SCHEME_NONE);/* do not under weight the long distances */ + + if (!sm) { + *flag = -1; + goto RETURN; + } + + + SparseStressMajorizationSmoother_smooth(sm, dim, *x, maxit_sm, 0.001); + for (i = 0; i < dim*m; i++) { + (*x)[i] /= sm->scaling; + } + SparseStressMajorizationSmoother_delete(sm); + + RETURN: + if (A != B) SparseMatrix_delete(A); + +} diff --git a/lib/sfdpgen/stress_model.h b/lib/sfdpgen/stress_model.h new file mode 100644 index 000000000..4e8c61276 --- /dev/null +++ b/lib/sfdpgen/stress_model.h @@ -0,0 +1,6 @@ +#ifndef STRESS_MODEL_H +#define STRESS_MODEL_H + +void stress_model(int dim, SparseMatrix A, real **x, int maxit, real tol, int *flag); + +#endif diff --git a/lib/sfdpgen/uniform_stress.c b/lib/sfdpgen/uniform_stress.c new file mode 100644 index 000000000..a0a72592a --- /dev/null +++ b/lib/sfdpgen/uniform_stress.c @@ -0,0 +1,165 @@ +#include "general.h" +#include "SparseMatrix.h" +#include "spring_electrical.h" +#include "post_process.h" +#include "sparse_solve.h" +#include +#include "uniform_stress.h" + + +UniformStressSmoother UniformStressSmoother_new(int dim, SparseMatrix A, real *x, real alpha, real M, int *flag){ + UniformStressSmoother sm; + int i, j, k, m = A->m, *ia = A->ia, *ja = A->ja, *iw, *jw, *id, *jd; + int nz; + real *d, *w, *a = (real*) A->a; + real diag_d, diag_w, dist, epsilon = 0.01; + + assert(SparseMatrix_is_symmetric(A, FALSE)); + + sm = MALLOC(sizeof(struct StressMajorizationSmoother_struct)); + sm->data = NULL; + sm->scheme = SM_SCHEME_UNIFORM_STRESS; + sm->lambda = NULL; + sm->data = MALLOC(sizeof(real)*2); + ((real*) sm->data)[0] = alpha; + ((real*) sm->data)[1] = M; + sm->data_deallocator = FREE; + + /* Lw and Lwd have diagonals */ + sm->Lw = SparseMatrix_new(m, m, A->nz + m, MATRIX_TYPE_REAL, FORMAT_CSR); + sm->Lwd = SparseMatrix_new(m, m, A->nz + m, MATRIX_TYPE_REAL, FORMAT_CSR); + iw = sm->Lw->ia; jw = sm->Lw->ja; + id = sm->Lwd->ia; jd = sm->Lwd->ja; + w = (real*) sm->Lw->a; d = (real*) sm->Lwd->a; + + if (!(sm->Lw) || !(sm->Lwd)) { + StressMajorizationSmoother_delete(sm); + return NULL; + } + + iw = sm->Lw->ia; jw = sm->Lw->ja; + id = sm->Lwd->ia; jd = sm->Lwd->ja; + w = (real*) sm->Lw->a; d = (real*) sm->Lwd->a; + iw[0] = id[0] = 0; + + nz = 0; + for (i = 0; i < m; i++){ + diag_d = diag_w = 0; + for (j = ia[i]; j < ia[i+1]; j++){ + k = ja[j]; + if (k != i){ + dist = MAX(ABS(a[j]), epsilon); + jd[nz] = jw[nz] = k; + w[nz] = -1/(dist*dist); + w[nz] = -1.; + d[nz] = w[nz]*dist; + diag_w += w[nz]; + diag_d += d[nz]; + nz++; + } + } + jd[nz] = jw[nz] = i; + w[nz] = -diag_w; + d[nz] = -diag_d; + nz++; + + iw[i+1] = nz; + id[i+1] = nz; + + } + + sm->Lw->nz = nz; + sm->Lwd->nz = nz; + + return sm; +} + + +void UniformStressSmoother_delete(UniformStressSmoother sm){ + + StressMajorizationSmoother_delete(sm); + +} + +real UniformStressSmoother_smooth(UniformStressSmoother sm, int dim, real *x, int maxit_sm) { + + return StressMajorizationSmoother_smooth(sm, dim, x, maxit_sm, 0.001); + +} + +SparseMatrix get_distance_matrix(SparseMatrix A, real scaling){ + /* get a distance matrix from a graph, at the moment we just symmetrize the matrix. At the moment if the matrix is not real, + we just assume distance of 1 among edges. Then we apply scaling to the entire matrix */ + SparseMatrix B; + real *val; + int i; + + if (A->type == MATRIX_TYPE_REAL){ + B = SparseMatrix_symmetrize(A, FALSE); + } else { + B = SparseMatrix_get_real_adjacency_matrix_symmetrized(A); + } + val = (real*) B->a; + if (scaling != 1) for (i = 0; i < B->nz; i++) val[i] *= scaling; + return B; +} + +extern void scale_to_box(real xmin, real ymin, real xmax, real ymax, int n, int dim, real *x); + +void uniform_stress(int dim, SparseMatrix A, real *x, int *flag){ + UniformStressSmoother sm; + real lambda0 = 10.1, M = 100, scaling = 1.; + real res; + int maxit = 300, samepoint = TRUE, i, k, n = A->m; + SparseMatrix B = NULL; + + *flag = 0; + + /* just set random initial for now */ + for (i = 0; i < dim*n; i++) { + x[i] = M*drand(); + } + + /* make sure x is not all at the same point */ + for (i = 1; i < n; i++){ + for (k = 0; k < dim; k++) { + if (ABS(x[0*dim+k] - x[i*dim+k]) > MACHINEACC){ + samepoint = FALSE; + i = n; + break; + } + } + } + + if (samepoint){ + srand(1); +#ifdef DEBUG_PRINT + fprintf(stderr,"input coordinates to uniform_stress are the same, use random coordinates as initial input"); +#endif + for (i = 0; i < dim*n; i++) x[i] = M*drand(); + } + + B = get_distance_matrix(A, scaling); + assert(SparseMatrix_is_symmetric(B, FALSE)); + + sm = UniformStressSmoother_new(dim, B, x, 1000000*lambda0, M, flag); + res = UniformStressSmoother_smooth(sm, dim, x, maxit); + UniformStressSmoother_delete(sm); + + sm = UniformStressSmoother_new(dim, B, x, 10000*lambda0, M, flag); + res = UniformStressSmoother_smooth(sm, dim, x, maxit); + UniformStressSmoother_delete(sm); + + sm = UniformStressSmoother_new(dim, B, x, 100*lambda0, M, flag); + res = UniformStressSmoother_smooth(sm, dim, x, maxit); + UniformStressSmoother_delete(sm); + + sm = UniformStressSmoother_new(dim, B, x, lambda0, M, flag); + res = UniformStressSmoother_smooth(sm, dim, x, maxit); + UniformStressSmoother_delete(sm); + + scale_to_box(0,0,7*70,10*70,A->m,dim,x);; + + SparseMatrix_delete(B); + +} diff --git a/lib/sfdpgen/uniform_stress.h b/lib/sfdpgen/uniform_stress.h new file mode 100644 index 000000000..40114d168 --- /dev/null +++ b/lib/sfdpgen/uniform_stress.h @@ -0,0 +1,18 @@ +/* +#ifndef UniformStress_H +#define UniformStress_H +*/ + +typedef StressMajorizationSmoother UniformStressSmoother; + +#define UniformStressSmoother_struct StressMajorizationSmoother_struct + +void UniformStressSmoother_delete(UniformStressSmoother sm); + +UniformStressSmoother UniformStressSmoother_new(int dim, SparseMatrix A, real *x, real alpha, real M, int *flag); + +void uniform_stress(int dim, SparseMatrix A, real *x, int *flag); + +/* +#endif +*/