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;
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;
}
};
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);
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
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
#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);
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;
}
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);
grid->level = 0;
grid->n = A->n;
grid->A = A;
+ grid->D = D;
grid->P = NULL;
grid->R = NULL;
grid->node_weights = node_weights;
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);
}
}
-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;
assert(A->m == A->n);
*cA = NULL;
+ *cD = NULL;
*P = NULL;
*R = NULL;
n = A->m;
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
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)) {
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)) {
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);
}
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;
}
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);
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);
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);
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:
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);
*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;
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);
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) {
}
#endif
A = grid->A;
+ D = grid->D;
if (grid->level >= ctrl->maxlevel - 1) {
#ifdef DEBUG_PRINT
if (Verbose) {
#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;
}
-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;
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;
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;
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);
#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
* Contributors: See CVS logs. Details at http://www.graphviz.org/
*************************************************************************/
-#include <assert.h>
-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);
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;
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)
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));
}
}
*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);
}
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.;
}
}
-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
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];
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;
/* 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));
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<<dim),QuadTree);
+ q->qts = MALLOC(sizeof(QuadTree)*(1<<dim));
for (i = 0; i < 1<<dim; i++) q->qts[i] = NULL;
}/* done adding new quadtree, now add points to them */
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<<dim; i++){
+ if (qt->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<<dim; i++){
+ QuadTree_get_nearest_internal(qt->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);
+
+
+}
#define QUAD_TREE_H
#include "LinkedList.h"
-#include "sfdpinternal.h"
+/* #include "sfdpinternal.h" */
#include <stdio.h>
typedef struct QuadTree_struct *QuadTree;
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
#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]|
*/
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;
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);
}
}
+ /* distance 2 neighbors */
for (j = ia[i]; j < ia[i+1]; j++){
k = ja[j];
for (l = ia[k]; l < ia[k+1]; l++){
s = stop/sbot;
for (i = 0; i < nz; i++) d[i] *= s;
+ sm->scaling = s;
sm->Lw->nz = nz;
sm->Lwd->nz = nz;
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;
-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;
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.;
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);
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){
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);
}
}
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;
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);
s = stop/sbot;
for (i = 0; i < iw[m]; i++) d[i] *= s;
+ sm->scaling = s;
FREE(avg_dist);
void TriangleSmoother_smooth(TriangleSmoother sm, int dim, real *x){
- StressMajorizationSmoother_smooth(sm, dim, x, 50);
+ StressMajorizationSmoother_smooth(sm, dim, x, 50, 0.001);
}
}
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;
#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;
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;
/*------------------------------------------------------------------*/
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
+
#include <ctype.h>
#include <spring_electrical.h>
#include <overlap.h>
-
-#ifdef DEBUG
-double _statistics[10];
-#endif
+#include <uniform_stress.h>
+#include <stress_model.h>
static void sfdp_init_edge(edge_t * e)
{
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);
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;
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));
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
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)) {
#include <sfdp.h>
-#ifdef DEBUG
-extern double _statistics[10];
-#endif
-
#endif
#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);
void Operator_matmul_delete(Operator o){
-
+ if (o) FREE(o);
}
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;
}
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){
rho_old = rho;
}
+ FREE(z); FREE(r); FREE(p); FREE(q);
#ifdef DEBUG
_statistics[0] += iter - 1;
#endif
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);
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
+/* $Id$Revision: */
/* vim:set shiftwidth=4 ts=8: */
/*************************************************************************
* Contributors: See CVS logs. Details at http://www.graphviz.org/
*************************************************************************/
-
-
+#include "SparseMatrix.h"
#include "spring_electrical.h"
#include "QuadTree.h"
#include "Multilevel.h"
#include <string.h>
#include <time.h>
-#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.*/
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;
}
FREE(ctrl);
}
+
+
enum {MAX_I = 20, OPT_UP = 1, OPT_DOWN = -1, OPT_INIT = 0};
struct oned_optimizer_struct{
int i;
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;
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)
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 {
}
}
- 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,",");
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,"{");
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);
}
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));
}
}
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;
}
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;
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);
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
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));
#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
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;
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);
}
#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;
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) {
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);
}
#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);
}
+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);
+
+}
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) {
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);
}
#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);
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;
*/
-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;
/* 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;
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;
}
+
+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;
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);
}
}
#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);
}
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);
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);
}
/* 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*/
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;
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);
--- /dev/null
+#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);
+
+}
--- /dev/null
+#ifndef STRESS_MODEL_H
+#define STRESS_MODEL_H
+
+void stress_model(int dim, SparseMatrix A, real **x, int maxit, real tol, int *flag);
+
+#endif
--- /dev/null
+#include "general.h"
+#include "SparseMatrix.h"
+#include "spring_electrical.h"
+#include "post_process.h"
+#include "sparse_solve.h"
+#include <time.h>
+#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);
+
+}
--- /dev/null
+/*
+#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
+*/
projects / graphviz / commitdiff