From 1b847b5abf50d2ab0701523b90a20306d0ee5528 Mon Sep 17 00:00:00 2001 From: ellson Date: Wed, 3 Jun 2009 01:10:59 +0000 Subject: [PATCH] try again to complete merge into HEAD --- cmd/tools/matrix_market.c | 550 +++++----- cmd/tools/matrix_market.h | 25 +- cmd/tools/mmio.c | 528 +++++----- cmd/tools/mmio.h | 48 +- lib/sfdpgen/LinkedList.c | 178 ++++ lib/sfdpgen/LinkedList.h | 55 + lib/sfdpgen/Makefile.old | 37 + lib/sfdpgen/Multilevel.c | 1129 +++++++++++++++++++++ lib/sfdpgen/Multilevel.h | 69 ++ lib/sfdpgen/QuadTree.c | 673 +++++++++++++ lib/sfdpgen/QuadTree.h | 61 ++ lib/sfdpgen/layout_similarity.c | 198 ++++ lib/sfdpgen/post_process.c | 690 +++++++++++++ lib/sfdpgen/post_process.h | 66 ++ lib/sfdpgen/sfdp.h | 31 + lib/sfdpgen/sfdpinit.c | 331 ++++++ lib/sfdpgen/sfdpinternal.h | 27 + lib/sfdpgen/sparse_solve.c | 232 +++++ lib/sfdpgen/sparse_solve.h | 36 + lib/sfdpgen/spring_electrical.c | 1660 +++++++++++++++++++++++++++++++ lib/sfdpgen/spring_electrical.h | 80 ++ 21 files changed, 6101 insertions(+), 603 deletions(-) create mode 100644 lib/sfdpgen/LinkedList.c create mode 100644 lib/sfdpgen/LinkedList.h create mode 100644 lib/sfdpgen/Makefile.old create mode 100644 lib/sfdpgen/Multilevel.c create mode 100644 lib/sfdpgen/Multilevel.h create mode 100644 lib/sfdpgen/QuadTree.c create mode 100644 lib/sfdpgen/QuadTree.h create mode 100644 lib/sfdpgen/layout_similarity.c create mode 100644 lib/sfdpgen/post_process.c create mode 100644 lib/sfdpgen/post_process.h create mode 100644 lib/sfdpgen/sfdp.h create mode 100644 lib/sfdpgen/sfdpinit.c create mode 100644 lib/sfdpgen/sfdpinternal.h create mode 100644 lib/sfdpgen/sparse_solve.c create mode 100644 lib/sfdpgen/sparse_solve.h create mode 100644 lib/sfdpgen/spring_electrical.c create mode 100644 lib/sfdpgen/spring_electrical.h diff --git a/cmd/tools/matrix_market.c b/cmd/tools/matrix_market.c index d96dc3804..a8c956164 100644 --- a/cmd/tools/matrix_market.c +++ b/cmd/tools/matrix_market.c @@ -1,326 +1,304 @@ -/* $Id$Revision: */ -/* vim:set shiftwidth=4 ts=8: */ - -/********************************************************** -* This software is part of the graphviz package * -* http://www.graphviz.org/ * -* * -* Copyright (c) 1994-2004 AT&T Corp. * -* and is licensed under the * -* Common Public License, Version 1.0 * -* by AT&T Corp. * -* * -* Information and Software Systems Research * -* AT&T Research, Florham Park NJ * -**********************************************************/ #include "SparseMatrix.h" #include "mmio.h" #include "matrix_market.h" -#include "render.h" #define MALLOC gmalloc #define REALLOC grealloc #define FREE free #define MEMCPY memcpy -int mm_get_type(MM_typecode typecode) -{ - if (mm_is_complex(typecode)) { - return MATRIX_TYPE_COMPLEX; - } else if (mm_is_real(typecode)) { - return MATRIX_TYPE_REAL; - } else if (mm_is_integer(typecode)) { - return MATRIX_TYPE_INTEGER; - } else if (mm_is_pattern(typecode)) { - return MATRIX_TYPE_PATTERN; - } - return MATRIX_TYPE_UNKNOWN; +int mm_get_type(MM_typecode typecode){ + if (mm_is_complex(typecode)){ + return MATRIX_TYPE_COMPLEX; + } else if (mm_is_real(typecode)) { + return MATRIX_TYPE_REAL; + } else if (mm_is_integer(typecode)) { + return MATRIX_TYPE_INTEGER; + } else if (mm_is_pattern(typecode)){ + return MATRIX_TYPE_PATTERN; + } + return MATRIX_TYPE_UNKNOWN; } -static void set_mm_typecode(int type, MM_typecode * typecode) -{ - switch (type) { - case MATRIX_TYPE_COMPLEX: - mm_set_complex(typecode); - break; - case MATRIX_TYPE_REAL: - mm_set_real(typecode); - break; - case MATRIX_TYPE_INTEGER: - mm_set_integer(typecode); - break; - case MATRIX_TYPE_PATTERN: - mm_set_pattern(typecode); - break; - default: - break; - } +static void set_mm_typecode(int type, MM_typecode *typecode){ + switch (type){ + case MATRIX_TYPE_COMPLEX: + mm_set_complex(typecode); + break; + case MATRIX_TYPE_REAL: + mm_set_real(typecode); + break; + case MATRIX_TYPE_INTEGER: + mm_set_integer(typecode); + break; + case MATRIX_TYPE_PATTERN: + mm_set_pattern(typecode); + break; + default: + break; + } } -SparseMatrix SparseMatrix_import_matrix_market(FILE * f, int format) -{ - int ret_code, type; - MM_typecode matcode; - real *val = NULL, *v; - int *vali = NULL, i, m, n, *I = NULL, *J = NULL, nz; - void *vp = NULL; - SparseMatrix A = NULL; - int nzold; - int c; +SparseMatrix SparseMatrix_import_matrix_market(FILE *f, int format){ + int ret_code, type; + MM_typecode matcode; + real *val = NULL, *v; + int *vali = NULL, i, m, n, *I = NULL, *J = NULL, nz; + void *vp = NULL; + SparseMatrix A = NULL; + int nzold; + int c; - if ((c = fgetc(f)) != '%') { - ungetc(c, f); - return NULL; - } + if ((c=fgetc(f)) != '%'){ ungetc(c, f); - if (mm_read_banner(f, &matcode) != 0) { + return NULL; + } + ungetc(c, f); + if (mm_read_banner(f, &matcode) != 0) { #ifdef DEBUG - printf("Could not process Matrix Market banner.\n"); + printf("Could not process Matrix Market banner.\n"); #endif - return NULL; - } - - - /* This is how one can screen matrix types if their application */ - /* only supports a subset of the Matrix Market data types. */ - - if (!mm_is_matrix(matcode) || !mm_is_sparse(matcode)) { - assert(0); - /* - printf("Sorry, this application does not support dense matrix"); - printf("Market Market type: [%s]\n", mm_typecode_to_str(matcode)); - */ - return NULL; - } - - /* find out size of sparse matrix .... */ - if ((ret_code = mm_read_mtx_crd_size(f, &m, &n, &nz)) != 0) { - assert(0); - return NULL; - } - /* reseve memory for matrices */ + return NULL; + } + + + /* This is how one can screen matrix types if their application */ + /* only supports a subset of the Matrix Market data types. */ + + if (!mm_is_matrix(matcode) || + !mm_is_sparse(matcode) ){ + assert(0); + /* + printf("Sorry, this application does not support dense matrix"); + printf("Market Market type: [%s]\n", mm_typecode_to_str(matcode)); + */ + return NULL; + } + + /* find out size of sparse matrix .... */ + if ((ret_code = mm_read_mtx_crd_size(f, &m, &n, &nz)) !=0){ + assert(0); + return NULL; + } + /* reseve memory for matrices */ + + I = MALLOC(nz * sizeof(int)); + J = MALLOC(nz * sizeof(int)); + + - I = MALLOC(nz * sizeof(int)); - J = MALLOC(nz * sizeof(int)); + switch (format){ + case FORMAT_CSC: + assert(0);/* not supported yet */ + break; + case FORMAT_CSR: + case FORMAT_COORD: - - - switch (format) { - case FORMAT_CSC: - assert(0); /* not supported yet */ - break; - case FORMAT_CSR: - case FORMAT_COORD: - - /* NOTE: when reading in doubles, ANSI C requires the use of the "l" */ - /* specifier as in "%lg", "%lf", "%le", otherwise errors will occur */ - /* (ANSI C X3.159-1989, Sec. 4.9.6.2, p. 136 lines 13-15) */ - type = mm_get_type(matcode); - switch (type) { - case MATRIX_TYPE_REAL: - val = (real *) malloc(nz * sizeof(real)); - for (i = 0; i < nz; i++) { - fscanf(f, "%d %d %lg\n", &I[i], &J[i], &val[i]); - I[i]--; /* adjust from 1-based to 0-based */ - J[i]--; - } - if (mm_is_symmetric(matcode)) { - I = REALLOC(I, 2 * sizeof(int) * nz); - J = REALLOC(J, 2 * sizeof(int) * nz); - val = REALLOC(val, 2 * sizeof(real) * nz); - nzold = nz; - for (i = 0; i < nzold; i++) { - if (I[i] != J[i]) { - I[nz] = J[i]; - J[nz] = I[i]; - val[nz++] = val[i]; - } - } - } else if (mm_is_skew(matcode)) { - I = REALLOC(I, 2 * sizeof(int) * nz); - J = REALLOC(J, 2 * sizeof(int) * nz); - val = REALLOC(val, 2 * sizeof(real) * nz); - vp = (void *) val; - nzold = nz; - for (i = 0; i < nzold; i++) { - assert(I[i] != J[i]); /* skew symm has no diag */ - I[nz] = J[i]; - J[nz] = I[i]; - val[nz++] = -val[i]; - } - } else { - assert(!mm_is_hermitian(matcode)); - } - vp = (void *) val; - break; - case MATRIX_TYPE_INTEGER: - vali = (int *) malloc(nz * sizeof(int)); - for (i = 0; i < nz; i++) { - fscanf(f, "%d %d %d\n", &I[i], &J[i], &vali[i]); - I[i]--; /* adjust from 1-based to 0-based */ - J[i]--; - } - if (mm_is_symmetric(matcode)) { - I = REALLOC(I, 2 * sizeof(int) * nz); - J = REALLOC(J, 2 * sizeof(int) * nz); - vali = REALLOC(vali, 2 * sizeof(int) * nz); - nzold = nz; - for (i = 0; i < nzold; i++) { - if (I[i] != J[i]) { - I[nz] = J[i]; - J[nz] = I[i]; - vali[nz++] = vali[i]; - } - } - } else if (mm_is_skew(matcode)) { - I = REALLOC(I, 2 * sizeof(int) * nz); - J = REALLOC(J, 2 * sizeof(int) * nz); - vali = REALLOC(vali, 2 * sizeof(int) * nz); - vp = (void *) val; - nzold = nz; - for (i = 0; i < nzold; i++) { - assert(I[i] != J[i]); /* skew symm has no diag */ - I[nz] = J[i]; - J[nz] = I[i]; - vali[nz++] = -vali[i]; - } - } else { - assert(!mm_is_hermitian(matcode)); - } - vp = (void *) vali; - break; - case MATRIX_TYPE_PATTERN: - for (i = 0; i < nz; i++) { - fscanf(f, "%d %d\n", &I[i], &J[i]); - I[i]--; /* adjust from 1-based to 0-based */ - J[i]--; - } - if (mm_is_symmetric(matcode) || mm_is_skew(matcode)) { - I = REALLOC(I, 2 * sizeof(int) * nz); - J = REALLOC(J, 2 * sizeof(int) * nz); - nzold = nz; - for (i = 0; i < nzold; i++) { - if (I[i] != J[i]) { - I[nz] = J[i]; - J[nz++] = I[i]; - } - } - } else { - assert(!mm_is_hermitian(matcode)); - } - break; - case MATRIX_TYPE_COMPLEX: - val = (real *) malloc(2 * nz * sizeof(real)); - v = val; - for (i = 0; i < nz; i++) { - fscanf(f, "%d %d %lg %lg\n", &I[i], &J[i], &v[0], &v[1]); - v += 2; - I[i]--; /* adjust from 1-based to 0-based */ - J[i]--; - } - if (mm_is_symmetric(matcode)) { - I = REALLOC(I, 2 * sizeof(int) * nz); - J = REALLOC(J, 2 * sizeof(int) * nz); - val = REALLOC(val, 4 * sizeof(real) * nz); - nzold = nz; - for (i = 0; i < nzold; i++) { - if (I[i] != J[i]) { - I[nz] = J[i]; - J[nz] = I[i]; - val[2 * nz] = val[2 * i]; - val[2 * nz + 1] = val[2 * i + 1]; - nz++; - } - } - } else if (mm_is_skew(matcode)) { - I = REALLOC(I, 2 * sizeof(int) * nz); - J = REALLOC(J, 2 * sizeof(int) * nz); - val = REALLOC(val, 4 * sizeof(real) * nz); - vp = (void *) val; - nzold = nz; - for (i = 0; i < nzold; i++) { - assert(I[i] != J[i]); /* skew symm has no diag */ - I[nz] = J[i]; - J[nz] = I[i]; - val[2 * nz] = -val[2 * i]; - val[2 * nz + 1] = -val[2 * i + 1]; - nz++; - - } - } else if (mm_is_hermitian(matcode)) { - I = REALLOC(I, 2 * sizeof(int) * nz); - J = REALLOC(J, 2 * sizeof(int) * nz); - val = REALLOC(val, 4 * sizeof(real) * nz); - vp = (void *) val; - nzold = nz; - for (i = 0; i < nzold; i++) { - if (I[i] != J[i]) { - I[nz] = J[i]; - J[nz] = I[i]; - val[2 * nz] = val[2 * i]; - val[2 * nz + 1] = -val[2 * i + 1]; - nz++; - } - } - } - vp = (void *) val; - break; - default: - return 0; + /* NOTE: when reading in doubles, ANSI C requires the use of the "l" */ + /* specifier as in "%lg", "%lf", "%le", otherwise errors will occur */ + /* (ANSI C X3.159-1989, Sec. 4.9.6.2, p. 136 lines 13-15) */ + type = mm_get_type(matcode); + switch (type) { + case MATRIX_TYPE_REAL: + val = (real *) malloc(nz * sizeof(real)); + for (i=0; itype, &matcode); + mm_initialize_typecode(&matcode); + mm_set_matrix(&matcode); + mm_set_sparse(&matcode); + mm_set_general(&matcode); + set_mm_typecode(A->type, &matcode); - mm_write_banner(file, matcode); - mm_write_comment(file, comment); + mm_write_banner(file, matcode); + mm_write_comment(file, comment); - SparseMatrix_export(file, A); + SparseMatrix_export(file, A); -} +} diff --git a/cmd/tools/matrix_market.h b/cmd/tools/matrix_market.h index 6bfcd9f71..50b9f8297 100644 --- a/cmd/tools/matrix_market.h +++ b/cmd/tools/matrix_market.h @@ -1,26 +1,5 @@ -/* $Id$Revision: */ -/* vim:set shiftwidth=4 ts=8: */ - -/********************************************************** -* This software is part of the graphviz package * -* http://www.graphviz.org/ * -* * -* Copyright (c) 1994-2004 AT&T Corp. * -* and is licensed under the * -* Common Public License, Version 1.0 * -* by AT&T Corp. * -* * -* Information and Software Systems Research * -* AT&T Research, Florham Park NJ * -**********************************************************/ -#ifndef MATRIX_MARKET_H -#define MATRIX_MARKET_H - #include "mmio.h" #include "SparseMatrix.h" int mm_get_type(MM_typecode typecode); -void SparseMatrix_export_matrix_market(FILE * file, SparseMatrix A, - char *comment); -SparseMatrix SparseMatrix_import_matrix_market(FILE * f, int format); - -#endif +void SparseMatrix_export_matrix_market(FILE *file, SparseMatrix A, char *comment); +SparseMatrix SparseMatrix_import_matrix_market(FILE *f, int format); diff --git a/cmd/tools/mmio.c b/cmd/tools/mmio.c index 9a45bd01f..2e4cb8738 100644 --- a/cmd/tools/mmio.c +++ b/cmd/tools/mmio.c @@ -1,18 +1,3 @@ -/* $Id$Revision: */ -/* vim:set shiftwidth=4 ts=8: */ - -/********************************************************** -* This software is part of the graphviz package * -* http://www.graphviz.org/ * -* * -* Copyright (c) 1994-2004 AT&T Corp. * -* and is licensed under the * -* Common Public License, Version 1.0 * -* by AT&T Corp. * -* * -* Information and Software Systems Research * -* AT&T Research, Florham Park NJ * -**********************************************************/ /* * Matrix Market I/O library for ANSI C * @@ -25,12 +10,13 @@ #include #include #include +/*#include */ #include #include "mmio.h" -int mm_read_unsymmetric_sparse(const char *fname, int *M_, int *N_, - int *nz_, double **val_, int **I_, int **J_) +int mm_read_unsymmetric_sparse(const char *fname, int *M_, int *N_, int *nz_, + double **val_, int **I_, int **J_) { FILE *f; MM_typecode matcode; @@ -38,166 +24,170 @@ int mm_read_unsymmetric_sparse(const char *fname, int *M_, int *N_, int i; double *val; int *I, *J; - + if ((f = fopen(fname, "r")) == NULL) - return -1; - - - if (mm_read_banner(f, &matcode) != 0) { - fprintf(stderr, - "mm_read_unsymetric: Could not process Matrix Market banner "); - fprintf(stderr, " in file [%s]\n", fname); - return -1; + return -1; + + + if (mm_read_banner(f, &matcode) != 0) + { + fprintf(stderr, "mm_read_unsymetric: Could not process Matrix Market banner "); + fprintf(stderr, " in file [%s]\n", fname); + return -1; } - - - - if (!(mm_is_real(matcode) && mm_is_matrix(matcode) && - mm_is_sparse(matcode))) { - fprintf(stderr, "Sorry, this application does not support "); - fprintf(stderr, "Market Market type: [%s]\n", - mm_typecode_to_str(matcode)); - return -1; + + + + if ( !(mm_is_real(matcode) && mm_is_matrix(matcode) && + mm_is_sparse(matcode))) + { + fprintf(stderr, "Sorry, this application does not support "); + fprintf(stderr, "Market Market type: [%s]\n", + mm_typecode_to_str(matcode)); + return -1; } - + /* find out size of sparse matrix: M, N, nz .... */ - - if (mm_read_mtx_crd_size(f, &M, &N, &nz) != 0) { - fprintf(stderr, - "read_unsymmetric_sparse(): could not parse matrix size.\n"); - return -1; + + if (mm_read_mtx_crd_size(f, &M, &N, &nz) !=0) + { + fprintf(stderr, "read_unsymmetric_sparse(): could not parse matrix size.\n"); + return -1; } - + *M_ = M; *N_ = N; *nz_ = nz; - + /* reseve memory for matrices */ - + I = (int *) malloc(nz * sizeof(int)); J = (int *) malloc(nz * sizeof(int)); val = (double *) malloc(nz * sizeof(double)); - + *val_ = val; *I_ = I; *J_ = J; - + /* NOTE: when reading in doubles, ANSI C requires the use of the "l" */ /* specifier as in "%lg", "%lf", "%le", otherwise errors will occur */ /* (ANSI C X3.159-1989, Sec. 4.9.6.2, p. 136 lines 13-15) */ - - for (i = 0; i < nz; i++) { - fscanf(f, "%d %d %lg\n", &I[i], &J[i], &val[i]); - I[i]--; /* adjust from 1-based to 0-based */ - J[i]--; + + for (i=0; idata = data; + head->next = NULL; + return head; +} + +void SingleLinkedList_delete(SingleLinkedList head, void (*linklist_deallocator)(void*)){ + SingleLinkedList next; + + if (!head) return; + do { + next = head->next; + if (head->data) linklist_deallocator(head->data); + if (head) FREE(head); + head = next; + } while (head); + +} + + +SingleLinkedList SingleLinkedList_prepend(SingleLinkedList l, void *data){ + SingleLinkedList head = SingleLinkedList_new(data); + head->next = l; + return head; +} + +void* SingleLinkedList_get_data(SingleLinkedList l){ + return l->data; +} + +SingleLinkedList SingleLinkedList_get_next(SingleLinkedList l){ + return l->next; +} +void SingleLinkedList_print(SingleLinkedList head, void (*linkedlist_print)(void*)){ + + if (!head) return; + do { + if (head->data) linkedlist_print(head->data); + head = head->next; + } while (head); + +} + + +DoubleLinkedList DoubleLinkedList_new(void *data){ + DoubleLinkedList head; + head = N_GNEW(1,struct DoubleLinkedList_struct); + head->data = data; + head->next = NULL; + head->prev = NULL; + return head; +} + +void DoubleLinkedList_delete(DoubleLinkedList head, void (*linklist_deallocator)(void*)){ + DoubleLinkedList next; + + if (!head) return; + do { + next = head->next; + if (head->data) linklist_deallocator(head->data); + if (head) FREE(head); + head = next; + } while (head); + +} + + +DoubleLinkedList DoubleLinkedList_prepend(DoubleLinkedList l, void *data){ + DoubleLinkedList head = DoubleLinkedList_new(data); + if (l){ + head->next = l; + l->prev = head; + } + return head; +} + +void* DoubleLinkedList_get_data(DoubleLinkedList l){ + return l->data; +} + +DoubleLinkedList DoubleLinkedList_get_next(DoubleLinkedList l){ + return l->next; +} + +void DoubleLinkedList_print(DoubleLinkedList head, void (*linkedlist_print)(void*)){ + + if (!head) return; + do { + if (head->data) linkedlist_print(head->data); + head = head->next; + } while (head); + +} + +void DoubleLinkedList_delete_element(DoubleLinkedList l, void (*linklist_deallocator)(void*), DoubleLinkedList *head){ + /* delete an entry in the chain of linked list. If the head changes due to this (if l is the first element in the list), update */ + DoubleLinkedList next, prev; + + if (l){ + next = l->next; + prev = l->prev; + + if (l->data) linklist_deallocator(l->data); + FREE(l); + l = NULL; + + if (next) next->prev = prev; + if (prev) prev->next = next; + if (!prev) *head = next; + } +} + + +/* +static void print_int(void *d){ + int *i = (int*) d; + printf("%d\n",*i); +} + +main(){ + DoubleLinkedList l, ll; + + int i, *j; + + for (;;){ + j = malloc(sizeof(int)); + j[0] = -1; + l = DoubleLinkedList_new((void*) j); + + for (i = 0; i < 10; i++){ + j = malloc(sizeof(int)); + j[0] = i; + l = DoubleLinkedList_prepend(l, (void*) j); + + } + DoubleLinkedList_print(l, print_int); + + ll = DoubleLinkedList_get_next(l); + DoubleLinkedList_delete_element(ll, free, &l); + printf("after delete 8\n"); + DoubleLinkedList_print(l, print_int); + + DoubleLinkedList_delete_element(l, free, &l); + printf("after delete first elemnt\n"); + DoubleLinkedList_print(l, print_int); + + DoubleLinkedList_delete(l, free); + } +} + +*/ diff --git a/lib/sfdpgen/LinkedList.h b/lib/sfdpgen/LinkedList.h new file mode 100644 index 000000000..86c5aa18c --- /dev/null +++ b/lib/sfdpgen/LinkedList.h @@ -0,0 +1,55 @@ +/* $Id$ $Revision$ */ +/* vim:set shiftwidth=4 ts=8: */ + +/********************************************************** +* This software is part of the graphviz package * +* http://www.graphviz.org/ * +* * +* Copyright (c) 1994-2004 AT&T Corp. * +* and is licensed under the * +* Common Public License, Version 1.0 * +* by AT&T Corp. * +* * +* Information and Software Systems Research * +* AT&T Research, Florham Park NJ * +**********************************************************/ + +#ifndef LINKED_LIST_H +#define LINKED_LIST_H + +typedef struct SingleLinkedList_struct* SingleLinkedList; + +struct SingleLinkedList_struct { + void *data; + SingleLinkedList next; +}; + +SingleLinkedList SingleLinkedList_new(void *data); +void SingleLinkedList_delete(SingleLinkedList head, void (*linklist_deallocator)(void*)); +SingleLinkedList SingleLinkedList_prepend(SingleLinkedList l, void *data); + +void* SingleLinkedList_get_data(SingleLinkedList l); + +SingleLinkedList SingleLinkedList_get_next(SingleLinkedList l); + +void SingleLinkedList_print(SingleLinkedList head, void (*linkedlist_print)(void*)); + + +typedef struct DoubleLinkedList_struct* DoubleLinkedList; + +struct DoubleLinkedList_struct { + void *data; + DoubleLinkedList next; + DoubleLinkedList prev; +}; + +DoubleLinkedList DoubleLinkedList_new(void *data); +void DoubleLinkedList_delete(DoubleLinkedList head, void (*linklist_deallocator)(void*)); +DoubleLinkedList DoubleLinkedList_prepend(DoubleLinkedList l, void *data); + +void* DoubleLinkedList_get_data(DoubleLinkedList l); + +DoubleLinkedList DoubleLinkedList_get_next(DoubleLinkedList l); + +void DoubleLinkedList_delete_element(DoubleLinkedList l, void (*linklist_deallocator)(void*), DoubleLinkedList *head); +#endif diff --git a/lib/sfdpgen/Makefile.old b/lib/sfdpgen/Makefile.old new file mode 100644 index 000000000..50571e479 --- /dev/null +++ b/lib/sfdpgen/Makefile.old @@ -0,0 +1,37 @@ +all: libsfdpgen.a +ROOT=../.. +include $(ROOT)/Config.mk +include $(ROOT)/makearch/$(ARCH) + +INCS = -I. -I$(ROOT) \ + -I../common \ + -I../neatogen \ + -I../pack \ + -I../gvc \ + -I../pathplan \ + -I../sparse \ + -I../rbtree \ + -I../graph \ + -I../cdt \ + -I../gd + +DEFINES = -DDEBUG -DHAVE_CONFIG_H + +OBJS = sfdpinit.o spring_electrical.o LinkedList.o \ + sparse_solve.o post_process.o \ + QuadTree.o Multilevel.o PriorityQueue.o + +libsfdpgen.a : $(OBJS) + $(RM) libsfdpgen.a + $(AR) cr libsfdpgen.a $(OBJS) + $(RANLIB) libsfdpgen.a + +install: libsfdpgen.a + $(MKPATH) $(LIBDIR) + $(INSTALL) libsfdpgen.a $(LIBDIR) + +clean: + $(RM) *.o core *.a + +distclean: clean + $(RM) lib*.so.* diff --git a/lib/sfdpgen/Multilevel.c b/lib/sfdpgen/Multilevel.c new file mode 100644 index 000000000..8af20bcd5 --- /dev/null +++ b/lib/sfdpgen/Multilevel.c @@ -0,0 +1,1129 @@ +/* $Id$ $Revision$ */ +/* vim:set shiftwidth=4 ts=8: */ + +/********************************************************** +* This software is part of the graphviz package * +* http://www.graphviz.org/ * +* * +* Copyright (c) 1994-2004 AT&T Corp. * +* and is licensed under the * +* Common Public License, Version 1.0 * +* by AT&T Corp. * +* * +* Information and Software Systems Research * +* AT&T Research, Florham Park NJ * +**********************************************************/ + +#include "Multilevel.h" +#include "PriorityQueue.h" +#include "memory.h" +#include "logic.h" +#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 ctrl; + + ctrl = N_GNEW(1,struct Multilevel_control_struct); + ctrl->minsize = 4; + 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;*/ + return ctrl; +} + +void Multilevel_control_delete(Multilevel_control ctrl){ + FREE(ctrl); +} + +static Multilevel Multilevel_init(SparseMatrix A, real *node_weights){ + Multilevel grid; + if (!A) return NULL; + assert(A->m == A->n); + grid = N_GNEW(1,struct Multilevel_struct); + grid->level = 0; + grid->n = A->n; + grid->A = A; + grid->P = NULL; + grid->R = NULL; + grid->node_weights = node_weights; + grid->next = NULL; + grid->prev = NULL; + grid->delete_top_level_A = FALSE; + return grid; +} + +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); + } else { + SparseMatrix_delete(grid->A); + } + } + SparseMatrix_delete(grid->P); + SparseMatrix_delete(grid->R); + if (grid->node_weights && grid->level > 0) FREE(grid->node_weights); + Multilevel_delete(grid->next); + FREE(grid); +} + +static void maximal_independent_vertex_set(SparseMatrix A, int randomize, int **vset, int *nvset, int *nzc){ + int i, ii, j, *ia, *ja, m, n, *p = NULL; + assert(A); + assert(SparseMatrix_known_strucural_symmetric(A)); + ia = A->ia; + ja = A->ja; + m = A->m; + n = A->n; + assert(n == m); + *vset = N_GNEW(m,int); + for (i = 0; i < m; i++) (*vset)[i] = MAX_IND_VTX_SET_U; + *nvset = 0; + *nzc = 0; + + if (!randomize){ + for (i = 0; i < m; i++){ + if ((*vset)[i] == MAX_IND_VTX_SET_U){ + (*vset)[i] = (*nvset)++; + for (j = ia[i]; j < ia[i+1]; j++){ + if (i == ja[j]) continue; + (*vset)[ja[j]] = MAX_IND_VTX_SET_F; + (*nzc)++; + } + } + } + } else { + p = random_permutation(m); + for (ii = 0; ii < m; ii++){ + i = p[ii]; + if ((*vset)[i] == MAX_IND_VTX_SET_U){ + (*vset)[i] = (*nvset)++; + for (j = ia[i]; j < ia[i+1]; j++){ + if (i == ja[j]) continue; + (*vset)[ja[j]] = MAX_IND_VTX_SET_F; + (*nzc)++; + } + } + } + FREE(p); + } + (*nzc) += *nvset; +} + + +static void maximal_independent_vertex_set_RS(SparseMatrix A, int randomize, int **vset, int *nvset, int *nzc){ + /* The Ruge-Stuben coarsening scheme. Initially all vertices are in the U set (with marker MAX_IND_VTX_SET_U), + with gain equal to their degree. Select vertex with highest gain into a C set (with + marker >= MAX_IND_VTX_SET_C), and their neighbors j in F set (with marker MAX_IND_VTX_SET_F). The neighbors of + j that are in the U set get their gains incremented by 1. So overall + gain[k] = |{neighbor of k in U set}|+2*|{neighbors of k in F set}|. + nzc is the number of entries in the restriction matrix + */ + int i, jj, ii, *p = NULL, j, k, *ia, *ja, m, n, gain, removed, nf = 0; + PriorityQueue q; + assert(A); + assert(SparseMatrix_known_strucural_symmetric(A)); + + ia = A->ia; + ja = A->ja; + m = A->m; + n = A->n; + assert(n == m); + *vset = N_GNEW(m,int); + for (i = 0; i < m; i++) { + (*vset)[i] = MAX_IND_VTX_SET_U; + } + *nvset = 0; + *nzc = 0; + + q = PriorityQueue_new(m, 2*(m-1)); + + if (!randomize){ + for (i = 0; i < m; i++) + PriorityQueue_push(q, i, ia[i+1] - ia[i]); + } else { + p = random_permutation(m); + for (ii = 0; ii < m; ii++){ + i = p[ii]; + PriorityQueue_push(q, i, ia[i+1] - ia[i]); + } + FREE(p); + } + + while (PriorityQueue_pop(q, &i, &gain)){ + assert((*vset)[i] == MAX_IND_VTX_SET_U); + (*vset)[i] = (*nvset)++; + for (j = ia[i]; j < ia[i+1]; j++){ + jj = ja[j]; + assert((*vset)[jj] == MAX_IND_VTX_SET_U || (*vset)[jj] == MAX_IND_VTX_SET_F); + if (i == jj) continue; + + if ((*vset)[jj] == MAX_IND_VTX_SET_U){ + removed = PriorityQueue_remove(q, jj); + assert(removed); + (*vset)[jj] = MAX_IND_VTX_SET_F; + nf++; + + for (k = ia[jj]; k < ia[jj+1]; k++){ + if (jj == ja[k]) continue; + if ((*vset)[ja[k]] == MAX_IND_VTX_SET_U){ + gain = PriorityQueue_get_gain(q, ja[k]); + assert(gain >= 0); + PriorityQueue_push(q, ja[k], gain + 1); + } + } + } + (*nzc)++; + } + } + (*nzc) += *nvset; + PriorityQueue_delete(q); + +} + + + +static void maximal_independent_edge_set(SparseMatrix A, int randomize, int **matching, int *nmatch){ + int i, ii, j, *ia, *ja, m, n, *p = NULL; + assert(A); + assert(SparseMatrix_known_strucural_symmetric(A)); + ia = A->ia; + ja = A->ja; + m = A->m; + n = A->n; + assert(n == m); + *matching = N_GNEW(m,int); + for (i = 0; i < m; i++) (*matching)[i] = i; + *nmatch = n; + + if (!randomize){ + for (i = 0; i < m; i++){ + for (j = ia[i]; j < ia[i+1]; j++){ + if (i == ja[j]) continue; + if ((*matching)[ja[j]] == ja[j] && (*matching)[i] == i){ + (*matching)[ja[j]] = i; + (*matching)[i] = ja[j]; + (*nmatch)--; + } + } + } + } else { + p = random_permutation(m); + for (ii = 0; ii < m; ii++){ + i = p[ii]; + for (j = ia[i]; j < ia[i+1]; j++){ + if (i == ja[j]) continue; + if ((*matching)[ja[j]] == ja[j] && (*matching)[i] == i){ + (*matching)[ja[j]] = i; + (*matching)[i] = ja[j]; + (*nmatch)--; + } + } + } + FREE(p); + } +} + + + +static void maximal_independent_edge_set_heavest_edge_pernode(SparseMatrix A, int randomize, int **matching, int *nmatch){ + int i, ii, j, *ia, *ja, m, n, *p = NULL; + real *a, amax = 0; + int first = TRUE, jamax = 0; + + assert(A); + assert(SparseMatrix_known_strucural_symmetric(A)); + ia = A->ia; + ja = A->ja; + m = A->m; + n = A->n; + assert(n == m); + *matching = N_GNEW(m,int); + for (i = 0; i < m; i++) (*matching)[i] = i; + *nmatch = n; + + assert(SparseMatrix_is_symmetric(A, FALSE)); + assert(A->type == MATRIX_TYPE_REAL); + + a = (real*) A->a; + if (!randomize){ + for (i = 0; i < m; i++){ + first = TRUE; + for (j = ia[i]; j < ia[i+1]; j++){ + if (i == ja[j]) continue; + if ((*matching)[ja[j]] == ja[j] && (*matching)[i] == i){ + if (first) { + amax = a[j]; + jamax = ja[j]; + first = FALSE; + } else { + if (a[j] > amax){ + amax = a[j]; + jamax = ja[j]; + } + } + } + } + if (!first){ + (*matching)[jamax] = i; + (*matching)[i] = jamax; + (*nmatch)--; + } + } + } else { + p = random_permutation(m); + for (ii = 0; ii < m; ii++){ + i = p[ii]; + if ((*matching)[i] != i) continue; + first = TRUE; + for (j = ia[i]; j < ia[i+1]; j++){ + if (i == ja[j]) continue; + if ((*matching)[ja[j]] == ja[j] && (*matching)[i] == i){ + if (first) { + amax = a[j]; + jamax = ja[j]; + first = FALSE; + } else { + if (a[j] > amax){ + amax = a[j]; + jamax = ja[j]; + } + } + } + } + if (!first){ + (*matching)[jamax] = i; + (*matching)[i] = jamax; + (*nmatch)--; + } + } + FREE(p); + } +} + + + + + +#define node_degree(i) (ia[(i)+1] - ia[(i)]) + +static void maximal_independent_edge_set_heavest_edge_pernode_leaves_first(SparseMatrix A, int randomize, int **cluster, int **clusterp, int *ncluster){ + int i, ii, j, *ia, *ja, m, n, *p = NULL, q; + real *a, amax = 0; + int first = TRUE, jamax = 0; + int *matched, nz, ncmax = 0, nz0, nzz,k ; + enum {UNMATCHED = -2, MATCHED = -1}; + + assert(A); + assert(SparseMatrix_known_strucural_symmetric(A)); + ia = A->ia; + ja = A->ja; + m = A->m; + n = A->n; + assert(n == m); + *cluster = N_GNEW(m,int); + *clusterp = N_GNEW((m+1),int); + matched = N_GNEW(m,int); + + for (i = 0; i < m; i++) matched[i] = i; + + assert(SparseMatrix_is_symmetric(A, FALSE)); + assert(A->type == MATRIX_TYPE_REAL); + + *ncluster = 0; + (*clusterp)[0] = 0; + nz = 0; + a = (real*) A->a; + if (!randomize){ + for (i = 0; i < m; i++){ + if (matched[i] == MATCHED || node_degree(i) != 1) continue; + q = ja[ia[i]]; + assert(matched[q] != MATCHED); + matched[q] = MATCHED; + (*cluster)[nz++] = q; + for (j = ia[q]; j < ia[q+1]; j++){ + if (q == ja[j]) continue; + if (node_degree(ja[j]) == 1){ + matched[ja[j]] = MATCHED; + (*cluster)[nz++] = ja[j]; + } + } + ncmax = MAX(ncmax, nz - (*clusterp)[*ncluster]); + nz0 = (*clusterp)[*ncluster]; + if (nz - nz0 <= MAX_CLUSTER_SIZE){ + (*clusterp)[++(*ncluster)] = nz; + } else { + (*clusterp)[++(*ncluster)] = ++nz0; + nzz = nz0; + for (k = nz0; k < nz && nzz < nz; k++){ + nzz += MAX_CLUSTER_SIZE - 1; + nzz = MIN(nz, nzz); + (*clusterp)[++(*ncluster)] = nzz; + } + } + + } + #ifdef DEBUG_print + if (Verbose) + fprintf(stderr, "%d leaves and parents for %d clusters, largest cluster = %d\n",nz, *ncluster, ncmax); +#endif + for (i = 0; i < m; i++){ + first = TRUE; + if (matched[i] == MATCHED) continue; + for (j = ia[i]; j < ia[i+1]; j++){ + if (i == ja[j]) continue; + if (matched[ja[j]] != MATCHED && matched[i] != MATCHED){ + if (first) { + amax = a[j]; + jamax = ja[j]; + first = FALSE; + } else { + if (a[j] > amax){ + amax = a[j]; + jamax = ja[j]; + } + } + } + } + if (!first){ + matched[jamax] = MATCHED; + matched[i] = MATCHED; + (*cluster)[nz++] = i; + (*cluster)[nz++] = jamax; + (*clusterp)[++(*ncluster)] = nz; + } + } + + /* dan yi dian, wu ban */ + for (i = 0; i < m; i++){ + if (matched[i] == i){ + (*cluster)[nz++] = i; + (*clusterp)[++(*ncluster)] = nz; + } + } + assert(nz == n); + + } else { + p = random_permutation(m); + for (ii = 0; ii < m; ii++){ + i = p[ii]; + if (matched[i] == MATCHED || node_degree(i) != 1) continue; + q = ja[ia[i]]; + assert(matched[q] != MATCHED); + matched[q] = MATCHED; + (*cluster)[nz++] = q; + for (j = ia[q]; j < ia[q+1]; j++){ + if (q == ja[j]) continue; + if (node_degree(ja[j]) == 1){ + matched[ja[j]] = MATCHED; + (*cluster)[nz++] = ja[j]; + } + } + ncmax = MAX(ncmax, nz - (*clusterp)[*ncluster]); + nz0 = (*clusterp)[*ncluster]; + if (nz - nz0 <= MAX_CLUSTER_SIZE){ + (*clusterp)[++(*ncluster)] = nz; + } else { + (*clusterp)[++(*ncluster)] = ++nz0; + nzz = nz0; + for (k = nz0; k < nz && nzz < nz; k++){ + nzz += MAX_CLUSTER_SIZE - 1; + nzz = MIN(nz, nzz); + (*clusterp)[++(*ncluster)] = nzz; + } + } + } + + #ifdef DEBUG_print + if (Verbose) + fprintf(stderr, "%d leaves and parents for %d clusters, largest cluster = %d\n",nz, *ncluster, ncmax); +#endif + for (ii = 0; ii < m; ii++){ + i = p[ii]; + first = TRUE; + if (matched[i] == MATCHED) continue; + for (j = ia[i]; j < ia[i+1]; j++){ + if (i == ja[j]) continue; + if (matched[ja[j]] != MATCHED && matched[i] != MATCHED){ + if (first) { + amax = a[j]; + jamax = ja[j]; + first = FALSE; + } else { + if (a[j] > amax){ + amax = a[j]; + jamax = ja[j]; + } + } + } + } + if (!first){ + matched[jamax] = MATCHED; + matched[i] = MATCHED; + (*cluster)[nz++] = i; + (*cluster)[nz++] = jamax; + (*clusterp)[++(*ncluster)] = nz; + } + } + + /* dan yi dian, wu ban */ + for (i = 0; i < m; i++){ + if (matched[i] == i){ + (*cluster)[nz++] = i; + (*clusterp)[++(*ncluster)] = nz; + } + } + + FREE(p); + } + + FREE(matched); +} + + + +static void maximal_independent_edge_set_heavest_edge_pernode_supernodes_first(SparseMatrix A, int randomize, int **cluster, int **clusterp, int *ncluster){ + int i, ii, j, *ia, *ja, m, n, *p = NULL; + real *a, amax = 0; + int first = TRUE, jamax = 0; + int *matched, nz, nz0; + enum {UNMATCHED = -2, MATCHED = -1}; + int nsuper, *super = NULL, *superp = NULL; + + assert(A); + assert(SparseMatrix_known_strucural_symmetric(A)); + ia = A->ia; + ja = A->ja; + m = A->m; + n = A->n; + assert(n == m); + *cluster = N_GNEW(m,int); + *clusterp = N_GNEW((m+1),int); + matched = N_GNEW(m,int); + + for (i = 0; i < m; i++) matched[i] = i; + + assert(SparseMatrix_is_symmetric(A, FALSE)); + assert(A->type == MATRIX_TYPE_REAL); + + SparseMatrix_decompose_to_supervariables(A, &nsuper, &super, &superp); + + *ncluster = 0; + (*clusterp)[0] = 0; + nz = 0; + a = (real*) A->a; + + for (i = 0; i < nsuper; i++){ + if (superp[i+1] - superp[i] <= 1) continue; + nz0 = (*clusterp)[*ncluster]; + for (j = superp[i]; j < superp[i+1]; j++){ + matched[super[j]] = MATCHED; + (*cluster)[nz++] = super[j]; + if (nz - nz0 >= MAX_CLUSTER_SIZE){ + (*clusterp)[++(*ncluster)] = nz; + nz0 = nz; + } + } + if (nz > nz0) (*clusterp)[++(*ncluster)] = nz; + } + + if (!randomize){ + for (i = 0; i < m; i++){ + first = TRUE; + if (matched[i] == MATCHED) continue; + for (j = ia[i]; j < ia[i+1]; j++){ + if (i == ja[j]) continue; + if (matched[ja[j]] != MATCHED && matched[i] != MATCHED){ + if (first) { + amax = a[j]; + jamax = ja[j]; + first = FALSE; + } else { + if (a[j] > amax){ + amax = a[j]; + jamax = ja[j]; + } + } + } + } + if (!first){ + matched[jamax] = MATCHED; + matched[i] = MATCHED; + (*cluster)[nz++] = i; + (*cluster)[nz++] = jamax; + (*clusterp)[++(*ncluster)] = nz; + } + } + + /* dan yi dian, wu ban */ + for (i = 0; i < m; i++){ + if (matched[i] == i){ + (*cluster)[nz++] = i; + (*clusterp)[++(*ncluster)] = nz; + } + } + assert(nz == n); + + } else { + p = random_permutation(m); + for (ii = 0; ii < m; ii++){ + i = p[ii]; + first = TRUE; + if (matched[i] == MATCHED) continue; + for (j = ia[i]; j < ia[i+1]; j++){ + if (i == ja[j]) continue; + if (matched[ja[j]] != MATCHED && matched[i] != MATCHED){ + if (first) { + amax = a[j]; + jamax = ja[j]; + first = FALSE; + } else { + if (a[j] > amax){ + amax = a[j]; + jamax = ja[j]; + } + } + } + } + if (!first){ + matched[jamax] = MATCHED; + matched[i] = MATCHED; + (*cluster)[nz++] = i; + (*cluster)[nz++] = jamax; + (*clusterp)[++(*ncluster)] = nz; + } + } + + /* dan yi dian, wu ban */ + for (i = 0; i < m; i++){ + if (matched[i] == i){ + (*cluster)[nz++] = i; + (*clusterp)[++(*ncluster)] = nz; + } + } + FREE(p); + + } + + FREE(super); + + FREE(superp); + + FREE(matched); +} + +static int scomp(const void *s1, const void *s2){ + real *ss1, *ss2; + ss1 = (real*) s1; + ss2 = (real*) s2; + + if ((ss1)[1] > (ss2)[1]){ + return -1; + } else if ((ss1)[1] < (ss2)[1]){ + return 1; + } + return 0; +} + +static void maximal_independent_edge_set_heavest_cluster_pernode_leaves_first(SparseMatrix A, int csize, + int randomize, int **cluster, int **clusterp, int *ncluster){ + int i, ii, j, *ia, *ja, m, n, *p = NULL, q, iv; + real *a; + int *matched, nz, nz0, nzz,k, nv; + enum {UNMATCHED = -2, MATCHED = -1}; + real *vlist; + + assert(A); + assert(SparseMatrix_known_strucural_symmetric(A)); + ia = A->ia; + ja = A->ja; + m = A->m; + n = A->n; + assert(n == m); + *cluster = N_GNEW(m,int); + *clusterp = N_GNEW((m+1),int); + matched = N_GNEW(m,int); + vlist = N_GNEW(2*m,real); + + for (i = 0; i < m; i++) matched[i] = i; + + assert(SparseMatrix_is_symmetric(A, FALSE)); + assert(A->type == MATRIX_TYPE_REAL); + + *ncluster = 0; + (*clusterp)[0] = 0; + nz = 0; + a = (real*) A->a; + + p = random_permutation(m); + for (ii = 0; ii < m; ii++){ + i = p[ii]; + if (matched[i] == MATCHED || node_degree(i) != 1) continue; + q = ja[ia[i]]; + assert(matched[q] != MATCHED); + matched[q] = MATCHED; + (*cluster)[nz++] = q; + for (j = ia[q]; j < ia[q+1]; j++){ + if (q == ja[j]) continue; + if (node_degree(ja[j]) == 1){ + matched[ja[j]] = MATCHED; + (*cluster)[nz++] = ja[j]; + } + } + nz0 = (*clusterp)[*ncluster]; + if (nz - nz0 <= MAX_CLUSTER_SIZE){ + (*clusterp)[++(*ncluster)] = nz; + } else { + (*clusterp)[++(*ncluster)] = ++nz0; + nzz = nz0; + for (k = nz0; k < nz && nzz < nz; k++){ + nzz += MAX_CLUSTER_SIZE - 1; + nzz = MIN(nz, nzz); + (*clusterp)[++(*ncluster)] = nzz; + } + } + } + + for (ii = 0; ii < m; ii++){ + i = p[ii]; + if (matched[i] == MATCHED) continue; + nv = 0; + for (j = ia[i]; j < ia[i+1]; j++){ + if (i == ja[j]) continue; + if (matched[ja[j]] != MATCHED && matched[i] != MATCHED){ + vlist[2*nv] = ja[j]; + vlist[2*nv+1] = a[j]; + nv++; + } + } + if (nv > 0){ + qsort(vlist, nv, sizeof(real)*2, scomp); + for (j = 0; j < MIN(csize - 1, nv); j++){ + iv = (int) vlist[2*j]; + matched[iv] = MATCHED; + (*cluster)[nz++] = iv; + } + matched[i] = MATCHED; + (*cluster)[nz++] = i; + (*clusterp)[++(*ncluster)] = nz; + } + } + + /* dan yi dian, wu ban */ + for (i = 0; i < m; i++){ + if (matched[i] == i){ + (*cluster)[nz++] = i; + (*clusterp)[++(*ncluster)] = nz; + } + } + FREE(p); + + + FREE(matched); +} +static void maximal_independent_edge_set_heavest_edge_pernode_scaled(SparseMatrix A, int randomize, int **matching, int *nmatch){ + int i, ii, j, *ia, *ja, m, n, *p = NULL; + real *a, amax = 0; + int first = TRUE, jamax = 0; + + assert(A); + assert(SparseMatrix_known_strucural_symmetric(A)); + ia = A->ia; + ja = A->ja; + m = A->m; + n = A->n; + assert(n == m); + *matching = N_GNEW(m,int); + for (i = 0; i < m; i++) (*matching)[i] = i; + *nmatch = n; + + assert(SparseMatrix_is_symmetric(A, FALSE)); + assert(A->type == MATRIX_TYPE_REAL); + + a = (real*) A->a; + if (!randomize){ + for (i = 0; i < m; i++){ + first = TRUE; + for (j = ia[i]; j < ia[i+1]; j++){ + if (i == ja[j]) continue; + if ((*matching)[ja[j]] == ja[j] && (*matching)[i] == i){ + if (first) { + amax = a[j]/(ia[i+1]-ia[i])/(ia[ja[j]+1]-ia[ja[j]]); + jamax = ja[j]; + first = FALSE; + } else { + if (a[j]/(ia[i+1]-ia[i])/(ia[ja[j]+1]-ia[ja[j]]) > amax){ + amax = a[j]/(ia[i+1]-ia[i])/(ia[ja[j]+1]-ia[ja[j]]); + jamax = ja[j]; + } + } + } + } + if (!first){ + (*matching)[jamax] = i; + (*matching)[i] = jamax; + (*nmatch)--; + } + } + } else { + p = random_permutation(m); + for (ii = 0; ii < m; ii++){ + i = p[ii]; + if ((*matching)[i] != i) continue; + first = TRUE; + for (j = ia[i]; j < ia[i+1]; j++){ + if (i == ja[j]) continue; + if ((*matching)[ja[j]] == ja[j] && (*matching)[i] == i){ + if (first) { + amax = a[j]/(ia[i+1]-ia[i])/(ia[ja[j]+1]-ia[ja[j]]); + jamax = ja[j]; + first = FALSE; + } else { + if (a[j]/(ia[i+1]-ia[i])/(ia[ja[j]+1]-ia[ja[j]]) > amax){ + amax = a[j]/(ia[i+1]-ia[i])/(ia[ja[j]+1]-ia[ja[j]]); + jamax = ja[j]; + } + } + } + } + if (!first){ + (*matching)[jamax] = i; + (*matching)[i] = jamax; + (*nmatch)--; + } + } + FREE(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; + int *irn = NULL, *jcn = NULL, *ia = NULL, *ja = NULL; + real *val = NULL; + SparseMatrix B = NULL; + int *vset = NULL, nvset, ncov, j; + int *cluster, *clusterp, ncluster; + + assert(A->m == A->n); + *cA = NULL; + *P = NULL; + *R = NULL; + n = A->m; + + *coarsen_scheme_used = ctrl->coarsen_scheme; + + switch (ctrl->coarsen_scheme){ + case COARSEN_HYBRID: +#ifdef DEBUG_PRINT + if (Verbose) + 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); + + if (!(*cA)) { +#ifdef DEBUG_PRINT + if (Verbose) + 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); + } + + if (!(*cA)) { +#ifdef DEBUG_PRINT + if (Verbose) + 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); + } + + if (!(*cA)) { +#ifdef DEBUG_PRINT + if (Verbose) + 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); + } + + + if (!(*cA)) { +#ifdef DEBUG_PRINT + if (Verbose) + 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); + } + ctrl->coarsen_scheme = COARSEN_HYBRID; + break; + case COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_EDGE_PERNODE_SUPERNODES_FIRST: + case COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_CLUSTER_PERNODE_LEAVES_FIRST: + case COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_EDGE_PERNODE_LEAVES_FIRST: + if (ctrl->coarsen_scheme == COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_EDGE_PERNODE_LEAVES_FIRST) { + maximal_independent_edge_set_heavest_edge_pernode_leaves_first(A, ctrl->randomize, &cluster, &clusterp, &ncluster); + } else if (ctrl->coarsen_scheme == COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_EDGE_PERNODE_SUPERNODES_FIRST) { + maximal_independent_edge_set_heavest_edge_pernode_supernodes_first(A, ctrl->randomize, &cluster, &clusterp, &ncluster); + } else { + maximal_independent_edge_set_heavest_cluster_pernode_leaves_first(A, 4, ctrl->randomize, &cluster, &clusterp, &ncluster); + } + assert(ncluster <= n); + nc = ncluster; + if (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); +#endif + goto RETURN; + } + irn = N_GNEW(n,int); + jcn = N_GNEW(n,int); + val = N_GNEW(n,real); + nzc = 0; + 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); + *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); + 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); + case COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_EDGE_PERNODE: + if (ctrl->coarsen_scheme == COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_EDGE_PERNODE) + maximal_independent_edge_set_heavest_edge_pernode(A, ctrl->randomize, &matching, &nmatch); + case COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_EDGE_PERNODE_DEGREE_SCALED: + 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) { +#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); +#endif + goto RETURN; + } + irn = N_GNEW(n,int); + jcn = N_GNEW(n,int); + val = N_GNEW(n,real); + nzc = 0; nc = 0; + for (i = 0; i < n; i++){ + if (matching[i] >= 0){ + if (matching[i] == i){ + irn[nzc] = i; + jcn[nzc] = nc; + val[nzc++] = 1.; + } else { + irn[nzc] = i; + jcn[nzc] = nc; + val[nzc++] = 1; + irn[nzc] = matching[i]; + jcn[nzc] = nc; + val[nzc++] = 1; + matching[matching[i]] = -1; + } + nc++; + matching[i] = -1; + } + } + assert(nc == nmatch); + 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); + 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_VERTEX_SET: + case COARSEN_INDEPENDENT_VERTEX_SET_RS: + if (ctrl->coarsen_scheme == COARSEN_INDEPENDENT_VERTEX_SET){ + maximal_independent_vertex_set(A, ctrl->randomize, &vset, &nvset, &nzc); + } else { + maximal_independent_vertex_set_RS(A, ctrl->randomize, &vset, &nvset, &nzc); + } + ia = A->ia; + ja = A->ja; + nc = nvset; + if (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); +#endif + goto RETURN; + } + irn = N_GNEW(nzc,int); + jcn = N_GNEW(nzc,int); + val = N_GNEW(nzc,real); + nzc = 0; + for (i = 0; i < n; i++){ + if (vset[i] == MAX_IND_VTX_SET_F){ + ncov = 0; + for (j = ia[i]; j < ia[i+1]; j++){ + if (vset[ja[j]] >= MAX_IND_VTX_SET_C){ + ncov++; + } + } + assert(ncov > 0); + for (j = ia[i]; j < ia[i+1]; j++){ + if (vset[ja[j]] >= MAX_IND_VTX_SET_C){ + irn[nzc] = i; + jcn[nzc] = vset[ja[j]]; + val[nzc++] = 1./(double) ncov; + } + } + } else { + assert(vset[i] >= MAX_IND_VTX_SET_C); + irn[nzc] = i; + jcn[nzc] = vset[i]; + val[nzc++] = 1.; + } + } + + *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); + 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); + break; + default: + goto RETURN; + } + RETURN: + if (matching) FREE(matching); + if (vset) FREE(vset); + if (irn) FREE(irn); + if (jcn) FREE(jcn); + if (val) FREE(val); + if (B) SparseMatrix_delete(B); +} + +void print_padding(int n){ + int i; + for (i = 0; i < n; i++) fputs (" ", stderr); +} +static Multilevel Multilevel_establish(Multilevel grid, Multilevel_control ctrl){ + Multilevel cgrid; + int coarsen_scheme_used; + real *cnode_weights = NULL; + SparseMatrix P, R, A, cA; + +#ifdef DEBUG_PRINT + if (Verbose) { + print_padding(grid->level); + fprintf(stderr, "level -- %d, n = %d, nz = %d nz/n = %f\n", grid->level, grid->n, grid->A->nz, grid->A->nz/(double) grid->n); + } +#endif + A = grid->A; + if (grid->level >= ctrl->maxlevel - 1) { +#ifdef DEBUG_PRINT + if (Verbose) { + print_padding(grid->level); + fprintf(stderr, " maxlevel reached, coarsening stops\n"); + } +#endif + return grid; + } + Multilevel_coarsen(A, &cA, grid->node_weights, &cnode_weights, &P, &R, ctrl, &coarsen_scheme_used); + if (!cA) return grid; + + cgrid = Multilevel_init(cA, 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->P = P; + grid->R = R; + cgrid->prev = grid; + cgrid = Multilevel_establish(cgrid, ctrl); + return grid; + +} + +Multilevel Multilevel_new(SparseMatrix A0, real *node_weights, Multilevel_control ctrl){ + Multilevel grid; + SparseMatrix A = A0; + + 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); + grid = Multilevel_establish(grid, ctrl); + if (A != A0) grid->delete_top_level_A = TRUE;/* be sure to clean up later */ + return grid; +} + + +Multilevel Multilevel_get_coarsest(Multilevel grid){ + while (grid->next){ + grid = grid->next; + } + return grid; +} + diff --git a/lib/sfdpgen/Multilevel.h b/lib/sfdpgen/Multilevel.h new file mode 100644 index 000000000..f1469632c --- /dev/null +++ b/lib/sfdpgen/Multilevel.h @@ -0,0 +1,69 @@ +/* $Id$ $Revision$ */ +/* vim:set shiftwidth=4 ts=8: */ + +/********************************************************** +* This software is part of the graphviz package * +* http://www.graphviz.org/ * +* * +* Copyright (c) 1994-2004 AT&T Corp. * +* and is licensed under the * +* Common Public License, Version 1.0 * +* by AT&T Corp. * +* * +* Information and Software Systems Research * +* AT&T Research, Florham Park NJ * +**********************************************************/ + +#ifndef MULTILEVEL_H +#define MULTILEVEL_H + +#include "SparseMatrix.h" + +typedef struct Multilevel_struct *Multilevel; + +struct Multilevel_struct { + int level;/* 0, 1, ... */ + int n; + SparseMatrix A; + SparseMatrix P; + SparseMatrix R; + real *node_weights; + Multilevel next; + Multilevel prev; + int delete_top_level_A; + int coarsen_scheme_used;/* to get from previous level to here */ +}; + +enum {MAX_IND_VTX_SET_U = -100, MAX_IND_VTX_SET_F = -1, MAX_IND_VTX_SET_C = 0}; + +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}; + + +struct Multilevel_control_struct { + int minsize; + real min_coarsen_factor; + int maxlevel; + int randomize; + int coarsen_scheme; +}; + +typedef struct Multilevel_control_struct *Multilevel_control; + +Multilevel_control Multilevel_control_new(void); + +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_get_coarsest(Multilevel grid); + +void print_padding(int n); + +#define Multilevel_is_finest(grid) (!((grid)->prev)) +#define Multilevel_is_coarsest(grid) (!((grid)->next)) + +#endif diff --git a/lib/sfdpgen/QuadTree.c b/lib/sfdpgen/QuadTree.c new file mode 100644 index 000000000..cc5fe6e02 --- /dev/null +++ b/lib/sfdpgen/QuadTree.c @@ -0,0 +1,673 @@ +/* $Id$ $Revision$ */ +/* vim:set shiftwidth=4 ts=8: */ + +/********************************************************** +* This software is part of the graphviz package * +* http://www.graphviz.org/ * +* * +* Copyright (c) 1994-2004 AT&T Corp. * +* and is licensed under the * +* Common Public License, Version 1.0 * +* by AT&T Corp. * +* * +* Information and Software Systems Research * +* AT&T Research, Florham Park NJ * +**********************************************************/ + +#include + +typedef double real; + +#include "QuadTree.h" +#include "memory.h" +#include "arith.h" +#include "math.h" +#define MALLOC gmalloc +#define REALLOC grealloc +#define FREE free +#define MEMCPY memcpy + +extern real distance_cropped(real *x, int dim, int i, int j); + +struct node_data_struct { + real node_weight; + real *coord; + real id; + void *data; +}; + +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->node_weight = weight; + nd->coord = N_GNEW(dim,real); + nd->id = id; + for (i = 0; i < dim; i++) nd->coord[i] = coord[i]; + nd->data = NULL; + return nd; +} + +void node_data_delete(void *d){ + node_data nd = (node_data) d; + FREE(nd->coord); + /*delete outside if (nd->data) FREE(nd->data);*/ + FREE(nd); +} + +real node_data_get_weight(void *d){ + node_data nd = (node_data) d; + return nd->node_weight; +} + +real* node_data_get_coord(void *d){ + node_data nd = (node_data) d; + return nd->coord; +} + +int node_data_get_id(void *d){ + node_data nd = (node_data) d; + return nd->id; +} + +#define node_data_get_data(d) (((node_data) (d))->data) + + +void check_or_realloc_arrays(int dim, int *nsuper, int *nsupermax, real **center, real **supernode_wgts, real **distances){ + + 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); + } +} + +void QuadTree_get_supernodes_internal(QuadTree qt, real bh, real *point, int nodeid, int *nsuper, int *nsupermax, real **center, real **supernode_wgts, real **distances, real *counts, int *flag){ + SingleLinkedList l; + real *coord, dist; + int dim, i; + + (*counts)++; + + if (!qt) return; + dim = qt->dim; + l = qt->l; + if (l){ + while (l){ + check_or_realloc_arrays(dim, nsuper, nsupermax, center, supernode_wgts, distances); + if (node_data_get_id(SingleLinkedList_get_data(l)) != nodeid){ + coord = node_data_get_coord(SingleLinkedList_get_data(l)); + for (i = 0; i < dim; i++){ + (*center)[dim*(*nsuper)+i] = coord[i]; + } + (*supernode_wgts)[*nsuper] = node_data_get_weight(SingleLinkedList_get_data(l)); + (*distances)[*nsuper] = point_distance(point, coord, dim); + (*nsuper)++; + } + l = SingleLinkedList_get_next(l); + } + } + + if (qt->qts){ + dist = point_distance(qt->center, point, dim); + if (qt->width < bh*dist){ + check_or_realloc_arrays(dim, nsuper, nsupermax, center, supernode_wgts, distances); + for (i = 0; i < dim; i++){ + (*center)[dim*(*nsuper)+i] = qt->average[i]; + } + (*supernode_wgts)[*nsuper] = qt->total_weight; + (*distances)[*nsuper] = point_distance(qt->average, point, dim); + (*nsuper)++; + } else { + for (i = 0; i < 1<qts[i], bh, point, nodeid, nsuper, nsupermax, center, + supernode_wgts, distances, counts, flag); + } + } + } + +} + +void QuadTree_get_supernodes(QuadTree qt, real bh, real *point, int nodeid, int *nsuper, + int *nsupermax, real **center, real **supernode_wgts, real **distances, double *counts, int *flag){ + int dim = qt->dim; + + (*counts) = 0; + + *nsuper = 0; + + *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); + QuadTree_get_supernodes_internal(qt, bh, point, nodeid, nsuper, nsupermax, center, supernode_wgts, distances, counts, flag); + +} + + +static real *get_or_assign_node_force(real *force, int i, SingleLinkedList l, int dim){ + + real *f = (real*) node_data_get_data(SingleLinkedList_get_data(l)); + + if (!f){ + node_data_get_data(SingleLinkedList_get_data(l)) = &(force[i*dim]); + f = (real*) node_data_get_data(SingleLinkedList_get_data(l)); + } + return f; +} +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); + force = (real*) qt->data; + for (i = 0; i < dim; i++) force[i] = 0.; + } + return force; +} + +static void QuadTree_repulsive_force_interact(QuadTree qt1, QuadTree qt2, real *x, real *force, real bh, real p, real KP, real *counts){ + /* calculate the all to all reopulsive force and accumulate on each node of the quadtree if an interaction is possible. + force[i*dim+j], j=1,...,dim is teh force on node i + */ + SingleLinkedList l1, l2; + real *x1, *x2, dist, wgt1, wgt2, f, *f1, *f2, w1, w2; + int dim, i, j, i1, i2, k; + QuadTree qt11, qt12; + + if (!qt1 || !qt2) return; + assert(qt1->n > 0 && qt2->n > 0); + dim = qt1->dim; + + l1 = qt1->l; + l2 = qt2->l; + + /* far enough, calculate repulsive force */ + dist = point_distance(qt1->average, qt2->average, dim); + if (qt1->width + qt2->width < bh*dist){ + counts[0]++; + x1 = qt1->average; + w1 = qt1->total_weight; + f1 = get_or_alloc_force_qt(qt1, dim); + x2 = qt2->average; + w2 = qt2->total_weight; + f2 = get_or_alloc_force_qt(qt2, dim); + assert(dist > 0); + for (k = 0; k < dim; k++){ + if (p == -1){ + f = w1*w2*KP*(x1[k] - x2[k])/(dist*dist); + } else { + f = w1*w2*KP*(x1[k] - x2[k])/pow(dist, 1.- p); + } + f1[k] += f; + f2[k] -= f; + } + return; + } + + + /* both at leaves, calculate repulsive force */ + if (l1 && l2){ + while (l1){ + x1 = node_data_get_coord(SingleLinkedList_get_data(l1)); + wgt1 = node_data_get_weight(SingleLinkedList_get_data(l1)); + i1 = node_data_get_id(SingleLinkedList_get_data(l1)); + f1 = get_or_assign_node_force(force, i1, l1, dim); + l2 = qt2->l; + while (l2){ + x2 = node_data_get_coord(SingleLinkedList_get_data(l2)); + wgt2 = node_data_get_weight(SingleLinkedList_get_data(l2)); + i2 = node_data_get_id(SingleLinkedList_get_data(l2)); + f2 = get_or_assign_node_force(force, i2, l2, dim); + if ((qt1 == qt2 && i2 < i1) || i1 == i2) { + l2 = SingleLinkedList_get_next(l2); + continue; + } + counts[1]++; + dist = distance_cropped(x, dim, i1, i2); + for (k = 0; k < dim; k++){ + if (p == -1){ + f = wgt1*wgt2*KP*(x1[k] - x2[k])/(dist*dist); + } else { + f = wgt1*wgt2*KP*(x1[k] - x2[k])/pow(dist, 1.- p); + } + f1[k] += f; + f2[k] -= f; + } + l2 = SingleLinkedList_get_next(l2); + } + l1 = SingleLinkedList_get_next(l1); + } + return; + } + + + /* identical, split one */ + if (qt1 == qt2){ + for (i = 0; i < 1<qts[i]; + for (j = i; j < 1<qts[j]; + QuadTree_repulsive_force_interact(qt11, qt12, x, force, bh, p, KP, counts); + } + } + } else { + /* split the one with bigger box, or one not at the last level */ + if (qt1->width > qt2->width && !l1){ + for (i = 0; i < 1<qts[i]; + QuadTree_repulsive_force_interact(qt11, qt2, x, force, bh, p, KP, counts); + } + } else if (qt2->width > qt1->width && !l2){ + for (i = 0; i < 1<qts[i]; + QuadTree_repulsive_force_interact(qt11, qt1, x, force, bh, p, KP, counts); + } + } else if (!l1){/* pick one that is not at the last level */ + for (i = 0; i < 1<qts[i]; + QuadTree_repulsive_force_interact(qt11, qt2, x, force, bh, p, KP, counts); + } + } else if (!l2){ + for (i = 0; i < 1<qts[i]; + QuadTree_repulsive_force_interact(qt11, qt1, x, force, bh, p, KP, counts); + } + } else { + assert(0); /* can be both at the leaf level since that should be catched at the beginning of this func. */ + } + } +} + +static void QuadTree_repulsive_force_accumulate(QuadTree qt, real *force, real *counts){ + /* push down forces on cells into the node level */ + real wgt, wgt2; + real *f, *f2; + SingleLinkedList l = qt->l; + int i, k, dim; + QuadTree qt2; + + dim = qt->dim; + wgt = qt->total_weight; + f = get_or_alloc_force_qt(qt, dim); + assert(wgt > 0); + counts[2]++; + + if (l){ + while (l){ + i = node_data_get_id(SingleLinkedList_get_data(l)); + f2 = get_or_assign_node_force(force, i, l, dim); + wgt2 = node_data_get_weight(SingleLinkedList_get_data(l)); + wgt2 = wgt2/wgt; + for (k = 0; k < dim; k++) f2[k] += wgt2*f[k]; + l = SingleLinkedList_get_next(l); + } + return; + } + + for (i = 0; i < 1<qts[i]; + if (!qt2) continue; + assert(qt2->n > 0); + f2 = get_or_alloc_force_qt(qt2, dim); + wgt2 = qt2->total_weight; + wgt2 = wgt2/wgt; + for (k = 0; k < dim; k++) f2[k] += wgt2*f[k]; + QuadTree_repulsive_force_accumulate(qt2, force, counts); + } + +} + +void QuadTree_get_repulvice_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 + we accumulate forces at the cell levels to teh node level + qt: the quadtree + x: current coordinates for node i is x[i*dim+j], j = 0, ..., dim-1 + force: the repulsice force, an array of length dim*nnodes, the force for node i is at force[i*dim+j], j = 0, ..., dim - 1 + bh: Barnes-Hut coefficient. If width_cell1+width_cell2 < bh*dist_between_cells, we treat each cell as a supernode. + p: the repulsive force power + KP: pow(K, 1 - p) + counts: array of size 4. + . counts[0]: number of cell-cell interaction + . counts[1]: number of cell-node interaction + . counts[2]: number of total cells in the quadtree + . Al normalized by dividing by number of nodes + */ + int n = qt->n, dim = qt->dim, i; + + for (i = 0; i < 4; i++) counts[i] = 0; + + *flag = 0; + + for (i = 0; i < dim*n; i++) force[i] = 0; + + QuadTree_repulsive_force_interact(qt, qt, x, force, bh, p, KP, counts); + QuadTree_repulsive_force_accumulate(qt, force, counts); + for (i = 0; i < 4; i++) counts[i] /= n; + +} +QuadTree QuadTree_new_from_point_list(int dim, int n, int max_level, real *coord, real *weight){ + /* form a new QuadTree data structure from a list of coordinates of n points + coord: of length n*dim, point i sits at [i*dim, i*dim+dim - 1] + weight: node weight of lentgth n. If NULL, unit weight assumed. + */ + real *xmin, *xmax, *center, width; + QuadTree qt = NULL; + int i, k; + + xmin = N_GNEW(dim,real); + xmax = N_GNEW(dim,real); + center = N_GNEW(dim,real); + if (!xmin || !xmax || !center) return NULL; + + for (i = 0; i < dim; i++) xmin[i] = coord[i]; + for (i = 0; i < dim; i++) xmax[i] = coord[i]; + + for (i = 1; i < n; i++){ + for (k = 0; k < dim; k++){ + xmin[k] = MIN(xmin[k], coord[i*dim+k]); + xmax[k] = MAX(xmax[k], coord[i*dim+k]); + } + } + + width = xmax[0] - xmin[0]; + for (i = 0; i < dim; i++) { + center[i] = (xmin[i] + xmax[i])*0.5; + width = MAX(width, xmax[i] - xmin[i]); + } + width *= 0.52; + qt = QuadTree_new(dim, center, width, max_level); + + if (weight){ + for (i = 0; i < n; i++){ + qt = QuadTree_add(qt, &(coord[i*dim]), weight[i], i); + } + } else { + for (i = 0; i < n; i++){ + qt = QuadTree_add(qt, &(coord[i*dim]), 1, i); + } + } + + + FREE(xmin); + FREE(xmax); + FREE(center); + return qt; +} + +QuadTree QuadTree_new(int dim, real *center, real width, int max_level){ + QuadTree q; + int i; + q = N_GNEW(1,struct QuadTree_struct); + q->dim = dim; + q->n = 0; + q->center = N_GNEW(dim,real); + for (i = 0; i < dim; i++) q->center[i] = center[i]; + assert(width > 0); + q->width = width; + q->total_weight = 0; + q->average = NULL; + q->qts = NULL; + q->l = NULL; + q->max_level = max_level; + q->data = NULL; + return q; +} + +void QuadTree_delete(QuadTree q){ + int i, dim; + if (!q) return; + dim = q->dim; + FREE(q->center); + FREE(q->average); + if (q->data) FREE(q->data); + if (q->qts){ + for (i = 0; i < 1<qts[i]); + } + FREE(q->qts); + } + SingleLinkedList_delete(q->l, node_data_delete); + FREE(q); +} + +static int QuadTree_get_quadrant(int dim, real *center, real *coord){ + /* find the quadrant that a point of coordinates coord is going into with reference to the center. + if coord - center == {+,-,+,+} = {1,0,1,1}, then it will sit in the i-quadrant where + i's binary representation is 1011 (that is, decimal 11). + */ + int d = 0, i; + + for (i = dim - 1; i >= 0; i--){ + if (coord[i] - center[i] < 0){ + d = 2*d; + } else { + d = 2*d+1; + } + } + return d; +} + +static QuadTree QuadTree_new_in_quadrant(int dim, real *center, real width, int max_level, int i){ + /* a new quadtree in quadrant i of the original cell. The original cell is centered at 'center". + The new cell have width "width". + */ + QuadTree qt; + int k; + + qt = QuadTree_new(dim, center, width, max_level); + center = qt->center;/* right now this has the center for the parent */ + for (k = 0; k < dim; k++){/* decompose child id into binary, if {1,0}, say, then + add {width/2, -width/2} to the parents' center + to get the child's center. */ + if (i%2 == 0){ + center[k] -= width; + } else { + center[k] += width; + } + i = (i - i%2)/2; + } + return qt; + +} +static QuadTree QuadTree_add_internal(QuadTree q, real *coord, real weight, int id, int level){ + int i, dim = q->dim, ii; + node_data nd = NULL; + + int max_level = q->max_level; + int idd; + + /* 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 (q->n == 0){ + /* if this node is empty right now */ + q->n = 1; + q->total_weight = weight; + q->average = N_GNEW(dim,real); + 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->l = SingleLinkedList_new(nd); + } else if (level < max_level){ + /* otherwise open up into 2^dim quadtrees unless the level is too high */ + 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[i] = NULL; + }/* done adding new quadtree, now add points to them */ + + /* insert the old node (if exist) and the current node into the appropriate child quadtree */ + ii = QuadTree_get_quadrant(dim, q->center, coord); + assert(ii < 1<= 0); + if (q->qts[ii] == NULL) q->qts[ii] = QuadTree_new_in_quadrant(q->dim, q->center, (q->width)/2, max_level, ii); + + q->qts[ii] = QuadTree_add_internal(q->qts[ii], coord, weight, id, level + 1); + assert(q->qts[ii]); + + if (q->l){ + idd = node_data_get_id(SingleLinkedList_get_data(q->l)); + assert(q->n == 1); + coord = node_data_get_coord(SingleLinkedList_get_data(q->l)); + weight = node_data_get_weight(SingleLinkedList_get_data(q->l)); + ii = QuadTree_get_quadrant(dim, q->center, coord); + assert(ii < 1<= 0); + + if (q->qts[ii] == NULL) q->qts[ii] = QuadTree_new_in_quadrant(q->dim, q->center, (q->width)/2, max_level, ii); + + q->qts[ii] = QuadTree_add_internal(q->qts[ii], coord, weight, idd, level + 1); + assert(q->qts[ii]); + + /* delete the old node data on parent */ + SingleLinkedList_delete(q->l, node_data_delete); + q->l = NULL; + } + + (q->n)++; + } else { + assert(!(q->qts)); + /* level is too high, append data in the linked list */ + (q->n)++; + q->total_weight += weight; + for (i = 0; i < q->dim; i++) q->average[i] = ((q->average[i])*q->n + coord[i])/(q->n + 1); + nd = node_data_new(q->dim, weight, coord, id); + assert(q->l); + q->l = SingleLinkedList_prepend(q->l, nd); + } + return q; +} + + +QuadTree QuadTree_add(QuadTree q, real *coord, real weight, int id){ + if (!q) return q; + return QuadTree_add_internal(q, coord, weight, id, 0); + +} + +static void draw_polygon(FILE *fp, int dim, real *center, real width){ + /* pliot the enclosing square */ + if (dim < 2 || dim > 3) return; + fprintf(fp,"(*in c*){Line[{"); + + if (dim == 2){ + fprintf(fp, "{%f, %f}", center[0] + width, center[1] + width); + fprintf(fp, ",{%f, %f}", center[0] - width, center[1] + width); + fprintf(fp, ",{%f, %f}", center[0] - width, center[1] - width); + fprintf(fp, ",{%f, %f}", center[0] + width, center[1] - width); + fprintf(fp, ",{%f, %f}", center[0] + width, center[1] + width); + } else if (dim == 3){ + /* top */ + fprintf(fp,"{"); + fprintf(fp, "{%f, %f, %f}", center[0] + width, center[1] + width, center[2] + width); + fprintf(fp, ",{%f, %f, %f}", center[0] - width, center[1] + width, center[2] + width); + fprintf(fp, ",{%f, %f, %f}", center[0] - width, center[1] - width, center[2] + width); + fprintf(fp, ",{%f, %f, %f}", center[0] + width, center[1] - width, center[2] + width); + fprintf(fp, ",{%f, %f, %f}", center[0] + width, center[1] + width, center[2] + width); + fprintf(fp,"},"); + /* bot */ + fprintf(fp,"{"); + fprintf(fp, "{%f, %f, %f}", center[0] + width, center[1] + width, center[2] - width); + fprintf(fp, ",{%f, %f, %f}", center[0] - width, center[1] + width, center[2] - width); + fprintf(fp, ",{%f, %f, %f}", center[0] - width, center[1] - width, center[2] - width); + fprintf(fp, ",{%f, %f, %f}", center[0] + width, center[1] - width, center[2] - width); + fprintf(fp, ",{%f, %f, %f}", center[0] + width, center[1] + width, center[2] - width); + fprintf(fp,"},"); + /* for sides */ + fprintf(fp,"{"); + fprintf(fp, "{%f, %f, %f}", center[0] + width, center[1] + width, center[2] - width); + fprintf(fp, ",{%f, %f, %f}", center[0] + width, center[1] + width, center[2] + width); + fprintf(fp,"},"); + + fprintf(fp,"{"); + fprintf(fp, "{%f, %f, %f}", center[0] - width, center[1] + width, center[2] - width); + fprintf(fp, ",{%f, %f, %f}", center[0] - width, center[1] + width, center[2] + width); + fprintf(fp,"},"); + + fprintf(fp,"{"); + fprintf(fp, "{%f, %f, %f}", center[0] + width, center[1] - width, center[2] - width); + fprintf(fp, ",{%f, %f, %f}", center[0] + width, center[1] - width, center[2] + width); + fprintf(fp,"},"); + + fprintf(fp,"{"); + fprintf(fp, "{%f, %f, %f}", center[0] - width, center[1] - width, center[2] - width); + fprintf(fp, ",{%f, %f, %f}", center[0] - width, center[1] - width, center[2] + width); + fprintf(fp,"}"); + } + + + fprintf(fp, "}]}(*end C*)"); + + +} +static void QuadTree_print_internal(FILE *fp, QuadTree q, int level){ + /* dump a quad tree in Mathematica format. */ + SingleLinkedList l, l0; + real *coord; + int i, dim; + + if (!q) return; + + draw_polygon(fp, q->dim, q->center, q->width); + dim = q->dim; + + l0 = l = q->l; + if (l){ + printf(",(*a*) {Red,"); + while (l){ + if (l != l0) printf(","); + coord = node_data_get_coord(SingleLinkedList_get_data(l)); + fprintf(fp, "(*node %d*) Point[{", node_data_get_id(SingleLinkedList_get_data(l))); + for (i = 0; i < dim; i++){ + if (i != 0) printf(","); + fprintf(fp, "%f",coord[i]); + } + fprintf(fp, "}]"); + l = SingleLinkedList_get_next(l); + } + fprintf(fp, "}"); + } + + if (q->qts){ + for (i = 0; i < 1<qts[i], level + 1); + fprintf(fp, "}"); + } + } + + +} + +void QuadTree_print(FILE *fp, QuadTree q){ + if (!fp) return; + if (q->dim == 2){ + fprintf(fp, "Graphics[{"); + } else if (q->dim == 3){ + fprintf(fp, "Graphics3D[{"); + } else { + return; + } + QuadTree_print_internal(fp, q, 0); + if (q->dim == 2){ + fprintf(fp, "}, PlotRange -> All, Frame -> True, FrameTicks -> True]\n"); + } else { + fprintf(fp, "}, PlotRange -> All]\n"); + } +} diff --git a/lib/sfdpgen/QuadTree.h b/lib/sfdpgen/QuadTree.h new file mode 100644 index 000000000..651488aaf --- /dev/null +++ b/lib/sfdpgen/QuadTree.h @@ -0,0 +1,61 @@ +/* $Id$ $Revision$ */ +/* vim:set shiftwidth=4 ts=8: */ + +/********************************************************** +* This software is part of the graphviz package * +* http://www.graphviz.org/ * +* * +* Copyright (c) 1994-2004 AT&T Corp. * +* and is licensed under the * +* Common Public License, Version 1.0 * +* by AT&T Corp. * +* * +* Information and Software Systems Research * +* AT&T Research, Florham Park NJ * +**********************************************************/ + +#ifndef QUAD_TREE_H +#define QUAD_TREE_H + +#include "LinkedList.h" +#include "sfdpinternal.h" +#include + +typedef struct QuadTree_struct *QuadTree; + +struct QuadTree_struct { + /* a data structure containing coordinates of n items, their average is in "average". + The current level is a square or cube of width "width", which is subdivided into + 2^dim QuadTrees qts. At the last level, all coordinates are stored in a single linked list l. + total_weight is the combined weights of the nodes */ + int n;/* number of items */ + real total_weight; + int dim; + real *center;/* center of the bounding box, array of dimension dim. Allocated inside */ + real width;/* center +/- width gives the lower/upper bound, so really width is the + "radius" */ + real *average;/* the average coordinates. Array of length dim. Allocated inside */ + QuadTree *qts;/* subtree . If dim = 2, there are 4, dim = 3 gives 8 */ + SingleLinkedList l; + int max_level; + void *data; +}; + + +QuadTree QuadTree_new(int dim, real *center, real width, int max_level); + +void QuadTree_delete(QuadTree q); + +QuadTree QuadTree_add(QuadTree q, real *coord, real weight, int id);/* coord is copied in */ + +void QuadTree_print(FILE *fp, QuadTree q); + +QuadTree QuadTree_new_from_point_list(int dim, int n, int max_level, real *coord, real *weight); + +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); +#endif diff --git a/lib/sfdpgen/layout_similarity.c b/lib/sfdpgen/layout_similarity.c new file mode 100644 index 000000000..d708ca6bc --- /dev/null +++ b/lib/sfdpgen/layout_similarity.c @@ -0,0 +1,198 @@ +/* $Id$ $Revision$ */ +#include +#include +#include +#include "general.h" +#include "SparseMatrix.h" +#include "call_tri.h" + +static double boundingbox_area(int n, double *x){ + double xmax, xmin, ymax, ymin; + int i; + xmax = xmin = x[0]; + ymax = ymin = x[1]; + for (i = 0; i < n; i++){ + xmax = MAX(xmax, x[2*i]); + xmin = MIN(xmin, x[2*i]); + ymax = MAX(ymax, x[2*i+1]); + ymin = MIN(ymin, x[2*i+1]); + } + return (xmax-xmin)*(ymax-ymin); + +} + + +#define MINDIST 0.000000001 + +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; +} + +static real get_r(int n, real *x, real *y, real theta, real *r, real *p, real *q){ + /* find the best scaling factor r* such that + f[theta, r] = Sum[(x[i,1] + p - r(Sin[theta] y[i,1] + Cos[theta] y[i,2]))^2+(x[i,2] + q - r(Cos[theta] y[i,1] -Sin[theta] y[i,2]))^2, i = 1, n] is minimized. Return the value f[theta, r*], + also return r. + */ + real top = 0, bot = 0, c, s, f = 0, yy, yx, xxavg = 0, xyavg = 0, yxavg = 0, yyavg = 0; + int i; + + for (i = 0; i < n; i++){ + xxavg += x[2*i]; + xyavg += x[2*i+1]; + } + xxavg /= n; + xyavg /= n; + + for (i = 0; i < n; i++){ + s = sin(theta); + c = cos(theta); + yx = c*y[2*i] - s*y[2*i+1]; + yy = s*y[2*i] + c*y[2*i+1]; + top += (x[2*i] - xxavg)*yx + (x[2*i+1] - xyavg)*yy; + bot += yx*yx+yy*yy; + yxavg += yx; + yyavg += yy; + } + *r = top/(bot - yxavg*yxavg/n - yyavg*yyavg/n); + + *p = yxavg/n*(*r)-xxavg; + *q = yyavg/n*(*r)-xyavg; + + for (i = 0; i < n; i++){ + s = sin(theta); + c = cos(theta); + yx = c*y[2*i] - s*y[2*i+1]; + yy = s*y[2*i] + c*y[2*i+1]; + f += (x[2*i] + *p - (*r)*yx)*(x[2*i] + *p - (*r)*yx)+(x[2*i+1] + *q - (*r)*yy)*(x[2*i+1] + *q - (*r)*yy); + } + return f; +} + +static real dispacement(int n, real *x, real *y, real *rmin, real *thetamin, real *pmin, real *qmin){ + /* find the distance between two sets of cordinates by finding a suitable + rotation and scaling to minimize the distance */ + real dist = 0, theta, normx = 0, r, d, p, q; + int i; + + *thetamin = 0.; + dist = get_r(n, x, y, 0., &r, &p, &q); + *rmin = r; *pmin = p; *qmin = q; + + for (i = 0; i < 180; i++){ + theta = 3.1415626*i/180; + d = get_r(n, x, y, theta, &r, &p, &q); + if (d < dist){ + dist = d; + *rmin = r; *pmin = p; *qmin = q; + *thetamin = theta; + } + } + + + for (i = 0; i < n; i++) normx += x[2*i]*x[2*i]+x[2*i+1]*x[2*i+1]; + return sqrt(dist/normx); +} + +int main(int argc, char *argv[]) +{ + + + char *infile; + + SparseMatrix B = NULL; + int dim = 2, n = 0, i, *ia, *ja, j, k, nz = 0; + real *x, *y, mean, dev, ratio, *z, r, theta, p, q, xx; + + FILE *fp; + + + if (argc < 2){ + fprintf(stderr,"Usage: similarity graph layout1 layout2\n"); + } + + infile = argv[1]; + fp = fopen(infile,"r"); + while (fscanf(fp,"%lf %lf\n",&xx, &xx) == 2) n++; + fclose(fp); + + infile = argv[1]; + fp = fopen(infile,"r"); + x = N_GNEW(n*2,real); + for (i = 0; i < n; i++){ + fscanf(fp,"%lf %lf\n",&(x[2*i]), &(x[2*i+1])); + } + fclose(fp); + + + infile = argv[2]; + fp = fopen(infile,"r"); + y = N_GNEW(n*2,real); + nz = 0; + for (i = 0; i < n; i++){ + if (fscanf(fp,"%lf %lf\n",&(y[2*i]), &(y[2*i+1])) != 2) goto ERROR; + } + fclose(fp); + + B = call_tri(n, 2, x); + + z = N_GNEW(B->nz,real); + ia = B->ia; ja = B->ja; + nz = 0; + for (i = 0; i < n; i++){ + for (j = ia[i]; j < ia[i+1]; j++){ + k = ja[j]; + if (k == i) continue; + z[nz++] = distance(y,dim,i,k)/distance_cropped(x,dim,i,k); + } + } + + /* mean */ + mean = 0; + for (i = 0; i < nz; i++) mean += z[i]; + mean /= nz; + + /* deviation*/ + dev = 0; + for (i = 0; i < nz; i++) { + dev += (z[i]-mean)*(z[i]-mean); + } + dev /= nz; + dev = sqrt(dev); + + /* bounding box area */ + ratio = boundingbox_area(n, y); + /*/MAX(boundingbox_area(n, x), 0.001);*/ + + + fprintf(stderr, "mean = %f std = %f disimilarity = %f area = %f badness = %f displacement = %f\n", + mean, dev, dev/mean, ratio, (dev/mean+1)*ratio, dispacement(n, x, y, &r, &theta, &p, &q)); + fprintf(stderr, "theta = %f scaling = %f, shift = {%f, %f}\n",theta, 1/r, p, q); + + + printf("%.2f %.2f %.2f\n",dev/mean, dispacement(n, x, y, &r, &theta, &p, &q),ratio/1000000.); + + SparseMatrix_delete(B); + FREE(x); + FREE(y); + FREE(z); + + + return 0; + + ERROR: + printf("- - -\n"); + return 0; +} + diff --git a/lib/sfdpgen/post_process.c b/lib/sfdpgen/post_process.c new file mode 100644 index 000000000..ec52d7050 --- /dev/null +++ b/lib/sfdpgen/post_process.c @@ -0,0 +1,690 @@ +/* $Id$ $Revision$ */ +/* vim:set shiftwidth=4 ts=8: */ + +/********************************************************** +* This software is part of the graphviz package * +* http://www.graphviz.org/ * +* * +* Copyright (c) 1994-2004 AT&T Corp. * +* and is licensed under the * +* Common Public License, Version 1.0 * +* by AT&T Corp. * +* * +* Information and Software Systems Research * +* AT&T Research, Florham Park NJ * +**********************************************************/ + +#ifdef HAVE_CONFIG_H +#include "config.h" +#endif + +#include +#include +#include + +#include "types.h" +#include "memory.h" +#include "globals.h" + +#include "sparse_solve.h" +#include "post_process.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]| + */ + SparseMatrix D; + int *ia, *ja, i, j, k, l, nz; + real *d; + int *mask = NULL; + real len, di, sum, sumd; + + assert(SparseMatrix_is_symmetric(A, FALSE)); + + D = SparseMatrix_copy(A); + ia = D->ia; + ja = D->ja; + if (D->type != MATRIX_TYPE_REAL){ + FREE(D->a); + D->type = MATRIX_TYPE_REAL; + D->a = N_GNEW(D->nz,real); + } + d = (real*) D->a; + + mask = N_GNEW(D->m,int); + for (i = 0; i < D->m; i++) mask[i] = -1; + + for (i = 0; i < D->m; i++){ + di = node_degree(i); + mask[i] = i; + for (j = ia[i]; j < ia[i+1]; j++){ + if (i == ja[j]) continue; + mask[ja[j]] = i; + } + for (j = ia[i]; j < ia[i+1]; j++){ + k = ja[j]; + if (i == k) continue; + len = di + node_degree(k); + for (l = ia[k]; l < ia[k+1]; l++){ + if (mask[ja[l]] == i) len--; + } + d[j] = len; + assert(len > 0); + } + + } + + sum = 0; sumd = 0; + nz = 0; + for (i = 0; i < D->m; i++){ + for (j = ia[i]; j < ia[i+1]; j++){ + if (i == ja[j]) continue; + nz++; + sum += distance(x, dim, i, ja[j]); + sumd += d[j]; + } + } + sum /= nz; sumd /= nz; + sum = sum/sumd; + + for (i = 0; i < D->m; i++){ + for (j = ia[i]; j < ia[i+1]; j++){ + if (i == ja[j]) continue; + d[j] = sum*d[j]; + } + } + + + return D; +} + +StressMajorizationSmoother StressMajorizationSmoother_new(SparseMatrix A, int dim, real lambda0, real *x, + int ideal_dist_scheme){ + StressMajorizationSmoother sm; + int i, j, k, l, m = A->m, *ia = A->ia, *ja = A->ja, *iw, *jw, *id, *jd; + int *mask, nz; + real *d, *w, *lambda; + real *avg_dist, diag_d, diag_w, *dd, dist, s = 0, stop = 0, sbot = 0; + SparseMatrix ID; + + assert(SparseMatrix_is_symmetric(A, FALSE)); + + ID = ideal_distance_matrix(A, dim, x); + dd = (real*) ID->a; + + sm = N_GNEW(1,struct StressMajorizationSmoother_struct); + lambda = sm->lambda = N_GNEW(m,real); + for (i = 0; i < m; i++) sm->lambda[i] = lambda0; + mask = N_GNEW(m,int); + + avg_dist = N_GNEW(m,real); + + for (i = 0; i < m ;i++){ + avg_dist[i] = 0; + nz = 0; + for (j = ia[i]; j < ia[i+1]; j++){ + if (i == ja[j]) continue; + avg_dist[i] += distance(x, dim, i, ja[j]); + nz++; + } + assert(nz > 0); + avg_dist[i] /= nz; + } + + + for (i = 0; i < m; i++) mask[i] = -1; + + nz = 0; + for (i = 0; i < m; i++){ + mask[i] = i; + for (j = ia[i]; j < ia[i+1]; j++){ + k = ja[j]; + if (mask[k] != i){ + mask[k] = i; + nz++; + } + } + for (j = ia[i]; j < ia[i+1]; j++){ + k = ja[j]; + for (l = ia[k]; l < ia[k+1]; l++){ + if (mask[ja[l]] != i){ + mask[ja[l]] = i; + 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++){ + mask[i] = i+m; + diag_d = diag_w = 0; + for (j = ia[i]; j < ia[i+1]; j++){ + k = ja[j]; + if (mask[k] != i+m){ + mask[k] = i+m; + + jw[nz] = k; + if (ideal_dist_scheme == IDEAL_GRAPH_DIST){ + dist = 1; + } else if (ideal_dist_scheme == IDEAL_AVG_DIST){ + dist = (avg_dist[i] + avg_dist[k])*0.5; + } else if (ideal_dist_scheme == IDEAL_POWER_DIST){ + dist = pow(distance_cropped(x,dim,i,k),.4); + } else { + fprintf(stderr,"ideal_dist_scheme value wrong"); + assert(0); + exit(1); + } + + /* + w[nz] = -1./(ia[i+1]-ia[i]+ia[ja[j]+1]-ia[ja[j]]); + w[nz] = -2./(avg_dist[i]+avg_dist[k]);*/ + /* w[nz] = -1.;*//* use unit weight for now, later can try 1/(deg(i)+deg(k)) */ + w[nz] = -1/(dist*dist); + + diag_w += w[nz]; + + jd[nz] = k; + /* + d[nz] = w[nz]*distance(x,dim,i,k); + d[nz] = w[nz]*dd[j]; + d[nz] = w[nz]*(avg_dist[i] + avg_dist[k])*0.5; + */ + d[nz] = w[nz]*dist; + stop += d[nz]*distance(x,dim,i,k); + sbot += d[nz]*dist; + diag_d += d[nz]; + + nz++; + } + } + + for (j = ia[i]; j < ia[i+1]; j++){ + k = ja[j]; + for (l = ia[k]; l < ia[k+1]; l++){ + if (mask[ja[l]] != i+m){ + mask[ja[l]] = i+m; + + if (ideal_dist_scheme == IDEAL_GRAPH_DIST){ + dist = 2; + } else if (ideal_dist_scheme == IDEAL_AVG_DIST){ + dist = (avg_dist[i] + 2*avg_dist[k] + avg_dist[ja[l]])*0.5; + } else if (ideal_dist_scheme == IDEAL_POWER_DIST){ + dist = pow(distance_cropped(x,dim,i,ja[l]),.4); + } else { + fprintf(stderr,"ideal_dist_scheme value wrong"); + assert(0); + exit(1); + } + + jw[nz] = ja[l]; + /* + w[nz] = -1/(ia[i+1]-ia[i]+ia[ja[l]+1]-ia[ja[l]]); + w[nz] = -2/(avg_dist[i] + 2*avg_dist[k] + avg_dist[ja[l]]);*/ + /* w[nz] = -1.;*//* use unit weight for now, later can try 1/(deg(i)+deg(k)) */ + + w[nz] = -1/(dist*dist); + + diag_w += w[nz]; + + jd[nz] = ja[l]; + /* + d[nz] = w[nz]*(distance(x,dim,i,k)+distance(x,dim,k,ja[l])); + d[nz] = w[nz]*(dd[j]+dd[l]); + d[nz] = w[nz]*(avg_dist[i] + 2*avg_dist[k] + avg_dist[ja[l]])*0.5; + */ + d[nz] = w[nz]*dist; + stop += d[nz]*distance(x,dim,ja[l],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; + for (i = 0; i < nz; i++) d[i] *= s; + + sm->Lw->nz = nz; + sm->Lwd->nz = nz; + + FREE(mask); + FREE(avg_dist); + SparseMatrix_delete(ID); + return sm; +} + +static real total_distance(int m, int dim, real* x, real* y){ + real total = 0, dist = 0; + int i, j; + + for (i = 0; i < m; i++){ + dist = 0.; + for (j = 0; j < dim; j++){ + dist += (y[i*dim+j] - x[i*dim+j])*(y[i*dim+j] - x[i*dim+j]); + } + total += sqrt(dist); + } + return total; + +} + + + +void StressMajorizationSmoother_smooth(StressMajorizationSmoother sm, int dim, real *x, int maxit_sm) { + 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; + + Lwdd = SparseMatrix_copy(Lwd); + m = Lw->m; + x0 = N_GNEW(dim*m,real); + if (!x0) goto RETURN; + + x0 = MEMCPY(x0, x, sizeof(real)*dim*m); + y = N_GNEW(dim*m,real); + if (!y) goto RETURN; + + id = Lwd->ia; jd = Lwd->ja; + d = (real*) Lwd->a; + dd = (real*) Lwdd->a; + + while (iter++ < maxit_sm && diff > tol){ + for (i = 0; i < m; i++){ + idiag = -1; + diag = 0.; + for (j = id[i]; j < id[i+1]; j++){ + if (i == jd[j]) { + idiag = j; + continue; + } + dd[j] = d[j]/distance_cropped(x, dim, i, jd[j]); + diag += dd[j]; + } + 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]; + } + } + + maxit = sqrt((double) m); + res = SparseMatrix_solve(Lw, dim, x, y, 0.01, maxit, SOLVE_METHOD_CG, &flag); + + if (flag) goto RETURN; + + diff = total_distance(m, dim, x, y)/m; +#ifdef DEBUG_PRINT + if (Verbose){ + fprintf(stderr, "Outer iter = %d, res = %g Stress Majorization diff = %g\n",iter, res, diff); + } +#endif + MEMCPY(x, y, sizeof(real)*m*dim); + } + +#ifdef DEBUG + _statistics[1] += iter-1; +#endif + + RETURN: + SparseMatrix_delete(Lwdd); + if (!x0) FREE(x0); + if (!y) FREE(y); +} + +void StressMajorizationSmoother_delete(StressMajorizationSmoother sm){ + if (!sm) return; + if (sm->Lw) SparseMatrix_delete(sm->Lw); + if (sm->Lwd) SparseMatrix_delete(sm->Lwd); + if (sm->lambda) FREE(sm->lambda); +} + + +TriangleSmoother TriangleSmoother_new(SparseMatrix A, int dim, real lambda0, real *x, int use_triangularization){ + TriangleSmoother sm; + int i, j, k, m = A->m, *ia = A->ia, *ja = A->ja, *iw, *jw, *id, *jd, jdiag, nz; + SparseMatrix B; + real *avg_dist, *lambda, *d, *w, diag_d, diag_w, dist; + real s = 0, stop = 0, sbot = 0; + + assert(SparseMatrix_is_symmetric(A, FALSE)); + + avg_dist = N_GNEW(m,real); + + for (i = 0; i < m ;i++){ + avg_dist[i] = 0; + nz = 0; + for (j = ia[i]; j < ia[i+1]; j++){ + if (i == ja[j]) continue; + avg_dist[i] += distance(x, dim, i, ja[j]); + nz++; + } + assert(nz > 0); + avg_dist[i] /= nz; + } + + sm = N_GNEW(1,struct TriangleSmoother_struct); + lambda = sm->lambda = N_GNEW(m,real); + for (i = 0; i < m; i++) sm->lambda[i] = lambda0; + + if (m > 2){ + if (use_triangularization){ + B= call_tri(m, dim, x); + } else { + B= call_tri2(m, dim, x); + } + } else { + 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); + + SparseMatrix_delete(B); + sm->Lwd = SparseMatrix_copy(sm->Lw); + if (!(sm->Lw) || !(sm->Lwd)) { + TriangleSmoother_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; + + for (i = 0; i < m; i++){ + diag_d = diag_w = 0; + jdiag = -1; + for (j = iw[i]; j < iw[i+1]; j++){ + k = jw[j]; + if (k == i){ + jdiag = j; + continue; + } + /* w[j] = -1./(ia[i+1]-ia[i]+ia[ja[j]+1]-ia[ja[j]]); + w[j] = -2./(avg_dist[i]+avg_dist[k]); + w[j] = -1.*/;/* use unit weight for now, later can try 1/(deg(i)+deg(k)) */ + dist = pow(distance_cropped(x,dim,i,k),0.6); + w[j] = 1/(dist*dist); + diag_w += w[j]; + + /* d[j] = w[j]*distance(x,dim,i,k); + d[j] = w[j]*(avg_dist[i] + avg_dist[k])*0.5;*/ + d[j] = w[j]*dist; + stop += d[j]*distance(x,dim,i,k); + sbot += d[j]*dist; + diag_d += d[j]; + + } + + lambda[i] *= (-diag_w);/* alternatively don't do that then we have a constant penalty term scaled by lambda0 */ + + assert(jdiag >= 0); + w[jdiag] = -diag_w + lambda[i]; + d[jdiag] = -diag_d; + } + + s = stop/sbot; + for (i = 0; i < iw[m]; i++) d[i] *= s; + + FREE(avg_dist); + + return sm; +} + +void TriangleSmoother_delete(TriangleSmoother sm){ + + StressMajorizationSmoother_delete(sm); + +} + +void TriangleSmoother_smooth(TriangleSmoother sm, int dim, real *x){ + + StressMajorizationSmoother_smooth(sm, dim, x, 50); +} + + + + +/* ================================ spring and spring-electrical based smoother ================ */ +SpringSmoother SpringSmoother_new(SparseMatrix A, int dim, spring_electrical_control ctrl, real *x){ + SpringSmoother sm; + int i, j, k, l, m = A->m, *ia = A->ia, *ja = A->ja, *id, *jd; + int *mask, nz; + real *d, *dd; + real *avg_dist; + SparseMatrix ID = NULL; + + assert(SparseMatrix_is_symmetric(A, FALSE)); + + ID = ideal_distance_matrix(A, dim, x); + dd = (real*) ID->a; + + sm = N_GNEW(1,struct SpringSmoother_struct); + mask = N_GNEW(m,int); + + avg_dist = N_GNEW(m,real); + + for (i = 0; i < m ;i++){ + avg_dist[i] = 0; + nz = 0; + for (j = ia[i]; j < ia[i+1]; j++){ + if (i == ja[j]) continue; + avg_dist[i] += distance(x, dim, i, ja[j]); + nz++; + } + assert(nz > 0); + avg_dist[i] /= nz; + } + + + for (i = 0; i < m; i++) mask[i] = -1; + + nz = 0; + for (i = 0; i < m; i++){ + mask[i] = i; + for (j = ia[i]; j < ia[i+1]; j++){ + k = ja[j]; + if (mask[k] != i){ + mask[k] = i; + nz++; + } + } + for (j = ia[i]; j < ia[i+1]; j++){ + k = ja[j]; + for (l = ia[k]; l < ia[k+1]; l++){ + if (mask[ja[l]] != i){ + mask[ja[l]] = i; + nz++; + } + } + } + } + + sm->D = SparseMatrix_new(m, m, nz, MATRIX_TYPE_REAL, FORMAT_CSR); + if (!(sm->D)){ + SpringSmoother_delete(sm); + return NULL; + } + + id = sm->D->ia; jd = sm->D->ja; + d = (real*) sm->D->a; + id[0] = 0; + + nz = 0; + for (i = 0; i < m; i++){ + mask[i] = i+m; + for (j = ia[i]; j < ia[i+1]; j++){ + k = ja[j]; + if (mask[k] != i+m){ + mask[k] = i+m; + jd[nz] = k; + d[nz] = (avg_dist[i] + avg_dist[k])*0.5; + d[nz] = dd[j]; + nz++; + } + } + + for (j = ia[i]; j < ia[i+1]; j++){ + k = ja[j]; + for (l = ia[k]; l < ia[k+1]; l++){ + if (mask[ja[l]] != i+m){ + mask[ja[l]] = i+m; + jd[nz] = ja[l]; + d[nz] = (avg_dist[i] + 2*avg_dist[k] + avg_dist[ja[l]])*0.5; + d[nz] = dd[j]+dd[l]; + nz++; + } + } + } + id[i+1] = nz; + } + sm->D->nz = nz; + sm->ctrl = spring_electrical_control_new(); + *(sm->ctrl) = *ctrl; + sm->ctrl->random_start = FALSE; + sm->ctrl->multilevels = 1; + sm->ctrl->step /= 2; + sm->ctrl->maxiter = 20; + + FREE(mask); + FREE(avg_dist); + SparseMatrix_delete(ID); + + return sm; +} + + +void SpringSmoother_delete(SpringSmoother sm){ + if (!sm) return; + if (sm->D) SparseMatrix_delete(sm->D); + if (sm->ctrl) spring_electrical_control_delete(sm->ctrl); +} + + + + +void SpringSmoother_smooth(SpringSmoother sm, SparseMatrix A, real *node_weights, int dim, real *x){ + int flag = 0; + + spring_electrical_spring_embedding(dim, A, sm->D, sm->ctrl, node_weights, x, &flag); + assert(!flag); + +} + +/*=============================== end of spring and spring-electrical based smoother =========== */ + +void post_process_smoothing(int dim, SparseMatrix A, spring_electrical_control ctrl, real *node_weights, real *x, int *flag){ +#ifdef TIME + clock_t cpu; +#endif + +#ifdef TIME + cpu = clock(); +#endif + *flag = 0; + + switch (ctrl->smoothing){ + case SMOOTHING_RNG: + case SMOOTHING_TRIANGLE:{ + TriangleSmoother sm; + + if (ctrl->smoothing == SMOOTHING_RNG){ + sm = TriangleSmoother_new(A, dim, 0, x, FALSE); + } else { + sm = TriangleSmoother_new(A, dim, 0, x, TRUE); + } + TriangleSmoother_smooth(sm, dim, x); + TriangleSmoother_delete(sm); + break; + } + case SMOOTHING_STRESS_MAJORIZATION_GRAPH_DIST: + case SMOOTHING_STRESS_MAJORIZATION_POWER_DIST: + case SMOOTHING_STRESS_MAJORIZATION_AVG_DIST: + { + StressMajorizationSmoother sm; + int k, dist_scheme = IDEAL_AVG_DIST; + + if (ctrl->smoothing == SMOOTHING_STRESS_MAJORIZATION_GRAPH_DIST){ + dist_scheme = IDEAL_GRAPH_DIST; + } else if (ctrl->smoothing == SMOOTHING_STRESS_MAJORIZATION_AVG_DIST){ + dist_scheme = IDEAL_AVG_DIST; + } else if (ctrl->smoothing == SMOOTHING_STRESS_MAJORIZATION_POWER_DIST){ + dist_scheme = IDEAL_POWER_DIST; + } + + for (k = 0; k < 1; k++){ + sm = StressMajorizationSmoother_new(A, dim, 0.05, x, dist_scheme); + StressMajorizationSmoother_smooth(sm, dim, x, 50); + StressMajorizationSmoother_delete(sm); + } + break; + } + case SMOOTHING_SPRING:{ + SpringSmoother sm; + int k; + + for (k = 0; k < 1; k++){ + sm = SpringSmoother_new(A, dim, ctrl, x); + SpringSmoother_smooth(sm, A, node_weights, dim, x); + SpringSmoother_delete(sm); + } + + break; + } + + }/* end switch between smoothing methods */ + +#ifdef TIME + if (Verbose) fprintf(stderr, "post processing %f\n",((real) (clock() - cpu)) / CLOCKS_PER_SEC); +#endif +} diff --git a/lib/sfdpgen/post_process.h b/lib/sfdpgen/post_process.h new file mode 100644 index 000000000..d9b2b499e --- /dev/null +++ b/lib/sfdpgen/post_process.h @@ -0,0 +1,66 @@ +/* $Id$ $Revision$ */ +/* vim:set shiftwidth=4 ts=8: */ + +/********************************************************** +* This software is part of the graphviz package * +* http://www.graphviz.org/ * +* * +* Copyright (c) 1994-2004 AT&T Corp. * +* and is licensed under the * +* Common Public License, Version 1.0 * +* by AT&T Corp. * +* * +* Information and Software Systems Research * +* AT&T Research, Florham Park NJ * +**********************************************************/ + +#ifndef POST_PROCESS_H +#define POST_PROCESS_H + +#include "spring_electrical.h" + +struct StressMajorizationSmoother_struct { + SparseMatrix Lw; + SparseMatrix Lwd; + real* lambda; +}; + +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); + +void StressMajorizationSmoother_smooth(StressMajorizationSmoother sm, int dim, real *x, int maxit); +/*-------------------- triangle/neirhborhood graph based smoother ------------------- */ +typedef StressMajorizationSmoother TriangleSmoother; + +#define TriangleSmoother_struct StressMajorizationSmoother_struct + +void TriangleSmoother_delete(TriangleSmoother sm); + +TriangleSmoother TriangleSmoother_new(SparseMatrix A, int dim, real lambda, real *x, int use_triangularization); + +void TriangleSmoother_smooth(TriangleSmoother sm, int dim, real *x); + + + +/*------------------ spring and spring-electrical based smoother */ + +struct SpringSmoother_struct { + SparseMatrix D; + spring_electrical_control ctrl; +}; + +typedef struct SpringSmoother_struct *SpringSmoother; + +SpringSmoother SpringSmoother_new(SparseMatrix A, int dim, spring_electrical_control ctrl, real *x); + +void SpringSmoother_delete(SpringSmoother sm); + +void SpringSmoother_smooth(SpringSmoother sm, SparseMatrix A, real *node_weights, int dim, real *x); +/*------------------------------------------------------------------*/ + +void post_process_smoothing(int dim, SparseMatrix A, spring_electrical_control ctrl, real *node_weights, real *x, int *flag); +#endif diff --git a/lib/sfdpgen/sfdp.h b/lib/sfdpgen/sfdp.h new file mode 100644 index 000000000..39a688322 --- /dev/null +++ b/lib/sfdpgen/sfdp.h @@ -0,0 +1,31 @@ +/* $Id$ $Revision$ */ +/* vim:set shiftwidth=4 ts=8: */ + +/********************************************************** +* This software is part of the graphviz package * +* http://www.graphviz.org/ * +* * +* Copyright (c) 1994-2004 AT&T Corp. * +* and is licensed under the * +* Common Public License, Version 1.0 * +* by AT&T Corp. * +* * +* Information and Software Systems Research * +* AT&T Research, Florham Park NJ * +**********************************************************/ + +#ifndef SFDP_H +#define SFDP_H + +#ifdef HAVE_CONFIG_H +#include "config.h" +#endif + +#include "render.h" + +void sfdp_layout (graph_t * g); +void sfdp_cleanup (graph_t * g); +int fdpAdjust (graph_t * g); + +#endif + diff --git a/lib/sfdpgen/sfdpinit.c b/lib/sfdpgen/sfdpinit.c new file mode 100644 index 000000000..f46b36a32 --- /dev/null +++ b/lib/sfdpgen/sfdpinit.c @@ -0,0 +1,331 @@ +/* $Id$ $Revision$ */ +/* vim:set shiftwidth=4 ts=8: */ + +/********************************************************** +* This software is part of the graphviz package * +* http://www.graphviz.org/ * +* * +* Copyright (c) 1994-2004 AT&T Corp. * +* and is licensed under the * +* Common Public License, Version 1.0 * +* by AT&T Corp. * +* * +* Information and Software Systems Research * +* AT&T Research, Florham Park NJ * +**********************************************************/ + + +#ifdef HAVE_CONFIG_H +#include "config.h" +#endif +#ifdef HAVE_LIMITS_H +#include +#else +#ifdef HAVE_VALUES_H +#include +#endif +#endif +#include +#include +#include +#include +#include +#include +#include +#include + +#ifdef DEBUG +double _statistics[10]; +#endif + +static void sfdp_init_node(node_t * n) +{ + common_init_node(n); + ND_pos(n) = N_NEW(GD_ndim(n->graph), double); + gv_nodesize(n, GD_flip(n->graph)); +} + +static void sfdp_init_edge(edge_t * e) +{ + common_init_edge(e); +} + +static void sfdp_init_node_edge(graph_t * g) +{ + node_t *n; + edge_t *e; +#if 0 + int nnodes = agnnodes(g); + attrsym_t *N_pos = agfindnodeattr(g, "pos"); +#endif + + for (n = agfstnode(g); n; n = agnxtnode(g, n)) { + sfdp_init_node(n); +#if 0 + FIX so that user positions works with multiscale + user_pos(N_pos, NULL, n, nnodes); +#endif + } + for (n = agfstnode(g); n; n = agnxtnode(g, n)) { + for (e = agfstout(g, n); e; e = agnxtout(g, e)) + sfdp_init_edge(e); + } +} + +static void sfdp_init_graph(Agraph_t * g) +{ + int outdim; + + setEdgeType(g, ET_LINE); + outdim = late_int(g, agfindgraphattr(g, "dimen"), 2, 2); + GD_ndim(agroot(g)) = late_int(g, agfindgraphattr(g, "dim"), outdim, 2); + Ndim = GD_ndim(agroot(g)) = MIN(GD_ndim(agroot(g)), MAXDIM); + GD_odim(agroot(g)) = MIN(outdim, Ndim); + sfdp_init_node_edge(g); +} + +/* getPos: + */ +static real *getPos(Agraph_t * g, spring_electrical_control ctrl) +{ + Agnode_t *n; + real *pos = N_NEW(Ndim * agnnodes(g), real); + int ix, i; + + 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; +} + +static void sfdpLayout(graph_t * g, spring_electrical_control ctrl, + pointf pad) +{ + real *sizes; + real *pos; + Agnode_t *n; + int flag, i; + + SparseMatrix A = makeMatrix(g, Ndim); + if (ctrl->overlap >= 0) + sizes = getSizes(g, pad); + else + sizes = NULL; + pos = getPos(g, ctrl); + + multilevel_spring_electrical_embedding(Ndim, A, ctrl, NULL, sizes, pos, &flag); + + for (n = agfstnode(g); n; n = agnxtnode(g, n)) { + real *npos = pos + (Ndim * ND_id(n)); + for (i = 0; i < Ndim; i++) { + ND_pos(n)[i] = npos[i]; + } + } + + free(sizes); + free(pos); + free(ctrl->posSet); + SparseMatrix_delete (A); +} + +static int +late_smooth (graph_t* g, Agsym_t* sym, int dflt) +{ + char* s; + int v; + int rv; + + if (!sym) return dflt; + s = agxget (g, sym->index); + if (isdigit(*s)) { +#if (HAVE_GTS || HAVE_TRIANGLE) + if ((v = atoi (s)) <= SMOOTHING_RNG) +#else + if ((v = atoi (s)) <= SMOOTHING_SPRING) +#endif + rv = v; + else + rv = dflt; + } + else if (isalpha(*s)) { + if (!strcasecmp(s, "avg_dist")) + rv = SMOOTHING_STRESS_MAJORIZATION_AVG_DIST; + else if (!strcasecmp(s, "graph_dist")) + rv = SMOOTHING_STRESS_MAJORIZATION_GRAPH_DIST; + else if (!strcasecmp(s, "none")) + rv = SMOOTHING_NONE; + else if (!strcasecmp(s, "power_dist")) + rv = SMOOTHING_STRESS_MAJORIZATION_POWER_DIST; +#if (HAVE_GTS || HAVE_TRIANGLE) + else if (!strcasecmp(s, "rng")) + rv = SMOOTHING_RNG; +#endif + else if (!strcasecmp(s, "spring")) + rv = SMOOTHING_SPRING; +#if (HAVE_GTS || HAVE_TRIANGLE) + else if (!strcasecmp(s, "triangle")) + rv = SMOOTHING_TRIANGLE; +#endif + else + rv = dflt; + } + else + rv = dflt; + return rv; +} + +static int +late_quadtree_scheme (graph_t* g, Agsym_t* sym, int dflt) +{ + char* s; + int v; + int rv; + + if (!sym) return dflt; + s = agxget (g, sym->index); + if (isdigit(*s)) { + if ((v = atoi (s)) <= QUAD_TREE_FAST && v >= QUAD_TREE_NONE){ + rv = v; + } else { + rv = dflt; + } + } else if (isalpha(*s)) { + if (!strcasecmp(s, "none") || !strcasecmp(s, "false") ){ + rv = QUAD_TREE_NONE; + } else if (!strcasecmp(s, "normal") || !strcasecmp(s, "true") || !strcasecmp(s, "yes")){ + rv = QUAD_TREE_NORMAL; + } else if (!strcasecmp(s, "fast")){ + rv = QUAD_TREE_FAST; + } else { + rv = dflt; + } + } else { + rv = dflt; + } + return rv; +} + + +/* tuneControl: + * Use user values to reset control + * + * Possible parameters: + * ctrl->use_node_weights + */ +static void +tuneControl (graph_t* g, spring_electrical_control ctrl) +{ + long seed; + int init; + + seed = ctrl->random_seed; + init = setSeed (g, INIT_RANDOM, &seed); + if (init != INIT_RANDOM) { + agerr(AGWARN, "sfdp only supports start=random\n"); + } + ctrl->random_seed = seed; + + ctrl->K = late_double(g, agfindgraphattr(g, "K"), -1.0, 0.0); + ctrl->p = -1.0*late_double(g, agfindgraphattr(g, "repulsiveforce"), -AUTOP, 0.0); + 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); +} + +void sfdp_layout(graph_t * g) +{ + int doAdjust; + sfdp_init_graph(g); + doAdjust = (Ndim == 2); + adjust_data am; + + if (agnnodes(g)) { + Agraph_t **ccs; + Agraph_t *sg; + int ncc; + int i; + expand_t sep; + pointf pad; + spring_electrical_control ctrl = spring_electrical_control_new(); + + tuneControl (g, ctrl); + graphAdjustMode(g, &am, "prism0"); + + if ((am.mode == AM_PRISM) && doAdjust) { + doAdjust = 0; /* overlap removal done in sfpd */ + ctrl->overlap = am.value; + ctrl->initial_scaling = am.scaling; + sep = sepFactor(g); + if (sep.doAdd) { + pad.x = PS2INCH(sep.x); + pad.y = PS2INCH(sep.y); + } else { + pad.x = PS2INCH(DFLT_MARGIN); + pad.y = PS2INCH(DFLT_MARGIN); + } + } + else { + /* Turn off overlap removal in sfdp if prism not used */ + ctrl->overlap = -1; + } + + ccs = ccomps(g, &ncc, 0); + if (ncc == 1) { + sfdpLayout(g, ctrl, pad); + if (doAdjust) removeOverlapWith(g, &am); + spline_edges(g); + } else { + pack_info pinfo; + getPackInfo(g, l_node, CL_OFFSET, &pinfo); + pinfo.doSplines = 1; + + for (i = 0; i < ncc; i++) { + sg = ccs[i]; + nodeInduce(sg); + sfdpLayout(sg, ctrl, pad); + if (doAdjust) removeOverlapWith(sg, &am); + } + spline_edges(g); + packSubgraphs(ncc, ccs, g, &pinfo); + } + for (i = 0; i < ncc; i++) { + agdelete(g, ccs[i]); + } + free(ccs); + spring_electrical_control_delete(ctrl); + } + + dotneato_postprocess(g); +} + +static void sfdp_cleanup_graph(graph_t * g) +{ +} + +void sfdp_cleanup(graph_t * g) +{ + node_t *n; + edge_t *e; + + for (n = agfstnode(g); n; n = agnxtnode(g, n)) { + for (e = agfstout(g, n); e; e = agnxtout(g, e)) { + gv_cleanup_edge(e); + } + gv_cleanup_node(n); + } + sfdp_cleanup_graph(g); +} + diff --git a/lib/sfdpgen/sfdpinternal.h b/lib/sfdpgen/sfdpinternal.h new file mode 100644 index 000000000..d2fecd431 --- /dev/null +++ b/lib/sfdpgen/sfdpinternal.h @@ -0,0 +1,27 @@ +/* $Id$ $Revision$ */ +/* vim:set shiftwidth=4 ts=8: */ + +/********************************************************** +* This software is part of the graphviz package * +* http://www.graphviz.org/ * +* * +* Copyright (c) 1994-2004 AT&T Corp. * +* and is licensed under the * +* Common Public License, Version 1.0 * +* by AT&T Corp. * +* * +* Information and Software Systems Research * +* AT&T Research, Florham Park NJ * +**********************************************************/ + +#ifndef SFDPINTERNAL_H +#define SFDPINTERNAL_H + +#include + +#ifdef DEBUG +extern double _statistics[10]; +#endif + +#endif + diff --git a/lib/sfdpgen/sparse_solve.c b/lib/sfdpgen/sparse_solve.c new file mode 100644 index 000000000..445743d08 --- /dev/null +++ b/lib/sfdpgen/sparse_solve.c @@ -0,0 +1,232 @@ +/* $Id$ $Revision$ */ +/* vim:set shiftwidth=4 ts=8: */ + +/********************************************************** +* This software is part of the graphviz package * +* http://www.graphviz.org/ * +* * +* Copyright (c) 1994-2004 AT&T Corp. * +* and is licensed under the * +* Common Public License, Version 1.0 * +* by AT&T Corp. * +* * +* Information and Software Systems Research * +* AT&T Research, Florham Park NJ * +**********************************************************/ + +#include +#include +#include "sparse_solve.h" +#include "sfdpinternal.h" +#include "memory.h" +#include "logic.h" +#include "math.h" +#include "arith.h" +#include "types.h" +#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; +} + +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; +} + +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; +} + +real *Operator_matmul_apply(Operator o, real *x, real *y){ + SparseMatrix A = (SparseMatrix) o->data; + SparseMatrix_multiply_vector(A, x, &y, FALSE); + return y; +} + +Operator Operator_matmul_new(SparseMatrix A){ + Operator o; + + o = N_GNEW(1,struct Operator_struct); + o->data = (void*) A; + o->Operator_apply = Operator_matmul_apply; + return o; +} + + +void Operator_matmul_delete(Operator o){ + +} + + +real* Operator_diag_precon_apply(Operator o, real *x, real *y){ + int i, m; + real *diag = (real*) o->data; + m = (int) diag[0]; + diag++; + for (i = 0; i < m; i++) y[i] = x[i]*diag[i]; + return y; +} + +Operator Operator_diag_precon_new(SparseMatrix A){ + 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 = N_GNEW(1,struct Operator_struct); + o->data = N_GNEW((A->m + 1),real); + diag = (real*) o->data; + + diag[0] = m; + diag++; + for (i = 0; i < m; i++){ + diag[i] = 1.; + for (j = ia[i]; j < ia[i+1]; j++){ + if (i == ja[j] && ABS(a[j]) > 0) diag[i] = 1./a[j]; + } + } + + o->Operator_apply = Operator_diag_precon_apply; + + return o; +} + +void Operator_diag_precon_delete(Operator o){ + FREE(o->data); +} + +static real conjugate_gradient(Operator A, Operator precon, int n, real *x, real *rhs, real tol, int maxit, int *flag){ + real *z, *r, *p, *q, res = 10*tol, alpha; + real rho = 1.0e20, rho_old = 1, res0, beta; + real* (*Ax)(Operator o, real *in, real *out) = A->Operator_apply; + real* (*Minvx)(Operator o, real *in, real *out) = precon->Operator_apply; + int iter = 0; + + z = N_GNEW(n,real); + r = N_GNEW(n,real); + p = N_GNEW(n,real); + q = N_GNEW(n,real); + + r = Ax(A, x, r); + r = vector_subtract_to(n, rhs, r); + + res0 = res = sqrt(vector_product(n, r, r))/n; +#ifdef DEBUG_PRINT + if (Verbose && 0){ + fprintf(stderr, " cg iter = %d, residual = %g\n", iter, res); + } +#endif + + while ((iter++) < maxit && res > tol*res0){ + z = Minvx(precon, r, z); + rho = vector_product(n, r, z); + + if (iter > 1){ + beta = rho/rho_old; + p = vector_saxpy(n, z, p, beta); + } else { + MEMCPY(p, z, sizeof(real)*n); + } + + q = Ax(A, p, q); + + alpha = rho/vector_product(n, p, q); + + x = vector_saxpy2(n, x, p, alpha); + r = vector_saxpy2(n, r, q, -alpha); + + res = sqrt(vector_product(n, r, r))/n; + +#ifdef DEBUG_PRINT + if (Verbose && 0){ + fprintf(stderr, " cg iter = %d, residual = %g\n", iter, res); + } +#endif + + + + rho_old = rho; + } +#ifdef DEBUG + _statistics[0] += iter - 1; +#endif + +#ifdef DEBUG_PRINT + if (Verbose && 0){ + fprintf(stderr, " cg iter = %d, residual = %g\n", iter, res); + } +#endif + 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; + + *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]; + } + } + Operator_matmul_delete(Ax); + Operator_diag_precon_delete(precond); + FREE(x); + FREE(b); + + break; + default: + assert(0); + break; + + } + return res; +} + diff --git a/lib/sfdpgen/sparse_solve.h b/lib/sfdpgen/sparse_solve.h new file mode 100644 index 000000000..c18d1ad2f --- /dev/null +++ b/lib/sfdpgen/sparse_solve.h @@ -0,0 +1,36 @@ +/* $Id$ $Revision$ */ +/* vim:set shiftwidth=4 ts=8: */ + +/********************************************************** +* This software is part of the graphviz package * +* http://www.graphviz.org/ * +* * +* Copyright (c) 1994-2004 AT&T Corp. * +* and is licensed under the * +* Common Public License, Version 1.0 * +* by AT&T Corp. * +* * +* Information and Software Systems Research * +* AT&T Research, Florham Park NJ * +**********************************************************/ + + +#ifndef SPARSE_SOLVER_H +#define SPARSE_SOLVER_H + +#include "SparseMatrix.h" + +enum {SOLVE_METHOD_CG}; + + +typedef struct Operator_struct *Operator; + +struct Operator_struct { + void *data; + real* (*Operator_apply)(Operator o, real *in, real *out); +}; + +real SparseMatrix_solve(SparseMatrix A, int dim, real *x0, real *rhs, real tol, int maxit, int method, int *flag); + +#endif + diff --git a/lib/sfdpgen/spring_electrical.c b/lib/sfdpgen/spring_electrical.c new file mode 100644 index 000000000..51f1d24b0 --- /dev/null +++ b/lib/sfdpgen/spring_electrical.c @@ -0,0 +1,1660 @@ +/* vim:set shiftwidth=4 ts=8: */ + +/********************************************************** +* This software is part of the graphviz package * +* http://www.graphviz.org/ * +* * +* Copyright (c) 1994-2004 AT&T Corp. * +* and is licensed under the * +* Common Public License, Version 1.0 * +* by AT&T Corp. * +* * +* Information and Software Systems Research * +* AT&T Research, Florham Park NJ * +**********************************************************/ + + + +#include "spring_electrical.h" +#include "QuadTree.h" +#include "Multilevel.h" +#include "post_process.h" +#include "overlap.h" +#include "types.h" +#include "memory.h" +#include "arith.h" +#include "logic.h" +#include "math.h" +#include "globals.h" +#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->p = AUTOP;/*a negativve number default to -1. repulsive force = dist^p */ + 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->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->tol = 0.001;/* minimum different between two subsequence config before terminating. ||x-xold||_infinity < tol/K */ + ctrl->maxiter = 500; + ctrl->cool = 0.90;/* default 0.9 */ + ctrl->step = 0.1; + ctrl->adaptive_cooling = TRUE; + ctrl->random_seed = 123; + ctrl->beautify_leaves = FALSE; + ctrl->use_node_weights = FALSE; + ctrl->smoothing = SMOOTHING_NONE; + ctrl->overlap = 0; + ctrl->tscheme = QUAD_TREE_NORMAL; + ctrl->posSet = NULL; + ctrl->initial_scaling = -4; + + return ctrl; +} + +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; + real work[MAX_I+1]; + int direction; +}; + +typedef struct oned_optimizer_struct *oned_optimizer; + + +static void oned_optimizer_delete(oned_optimizer opt){ + FREE(opt); +} + +static oned_optimizer oned_optimizer_new(int i){ + oned_optimizer opt; + opt = N_GNEW(1,struct oned_optimizer_struct); + opt->i = i; + opt->direction = OPT_INIT; + return opt; +} + +static void oned_optimizer_train(oned_optimizer opt, real work){ + int i = opt->i; + + assert(i >= 0); + opt->work[i] = work; + if (opt->direction == OPT_INIT){ + if (opt->i == MAX_I){ + opt->direction = OPT_DOWN; + opt->i = opt->i - 1; + } else { + opt->direction = OPT_UP; + opt->i = MIN(MAX_I, opt->i + 1); + } + } else if (opt->direction == OPT_UP){ + /* fprintf(stderr, "{current_level, prev_level} = {%d,%d}, {work, work_prev} = {%f,%f}",i,i-1,opt->work[i], opt->work[i-1]);*/ + assert(i >= 1); + if (opt->work[i] < opt->work[i-1] && opt->i < MAX_I){ + /* fprintf(stderr, "keep going up to level %d\n",opt->i+1);*/ + opt->i = MIN(MAX_I, opt->i + 1); + } else { + /* fprintf(stderr, "going down to level %d\n",opt->i-1);*/ + (opt->i)--; + opt->direction = OPT_DOWN; + } + } else { + assert(i < MAX_I); + /* fprintf(stderr, "{current_level, prev_level} = {%d,%d}, {work, work_prev} = {%f,%f}",i,i+1,opt->work[i], opt->work[i+1]);*/ + if (opt->work[i] < opt->work[i+1] && opt->i > 0){ + /* fprintf(stderr, "keep going down to level %d\n",opt->i-1);*/ + opt->i = MAX(0, opt->i-1); + } else { + /* fprintf(stderr, "keep up to level %d\n",opt->i+1);*/ + (opt->i)++; + opt->direction = OPT_UP; + } + } +} + +static int oned_optimizer_get(oned_optimizer opt){ + return opt->i; +} + + +real average_edge_length(SparseMatrix A, int dim, real *coord){ + real dist = 0, d; + int *ia = A->ia, *ja = A->ja, i, j, k; + assert(SparseMatrix_is_symmetric(A, TRUE)); + + if (ia[A->m] == 0) return 1; + for (i = 0; i < A->m; i++){ + for (j = ia[i]; j < ia[i+1]; j++){ + d = 0; + for (k = 0; k < dim; k++){ + d += (coord[dim*i+k] - coord[dim*ja[j]])*(coord[dim*i+k] - coord[dim*ja[j]]); + } + dist += sqrt(d); + } + } + 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) + which should equal to -force (force = -gradient), + hence energy for force ||x-y||^m (y-x) is ||x-y||^(m+2)*2/(m+2) where m != 2 + 2. Grad[Log[||x-y||],x] = 1/||x-y||*0.5*(x-y)/||x-y|| = 0.5*(x-y)/||x-y||^2, + hence the energy to give force ||x-y||^-2 (x-y) is -2*Log[||x-y||] + + */ + int i, j, k, *ia = A->ia, *ja = A->ja, n = A->m; + real energy = 0, dist; + + for (i = 0; i < n; i++){ + /* attractive force C^((2-p)/3) ||x_i-x_j||/K * (x_j - x_i) */ + for (j = ia[i]; j < ia[i+1]; j++){ + if (ja[j] == i) continue; + dist = distance(x, dim, i, ja[j]); + energy += CRK*pow(dist, 3.)*2./3.; + } + + /* repulsive force K^(1 - p)/||x_i-x_j||^(1 - p) (x_i - x_j) */ + 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){ + energy += -KP*2*log(dist); + } else { + energy += -KP*pow(dist,p+1)*2/(p+1); + } + } + } + } + return energy; +} + +#endif + +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; + if (dim == 2){ + fprintf(fp,"Graphics[{GrayLevel[0.5],Line[{"); + } else { + fprintf(fp,"Graphics3D[{GrayLevel[0.5],Line[{"); + } + for (i = 0; i < A->m; i++){ + for (j = ia[i]; j < ia[i+1]; j++){ + if (ja[j] == i) continue; + ne++; + if (ne > 1) fprintf(fp, ","); + fprintf(fp, "{{"); + for (k = 0; k < dim; k++) { + if (k > 0) fprintf(fp,","); + fprintf(fp, "%f",x[i*dim+k]); + } + fprintf(fp, "},{"); + for (k = 0; k < dim; k++) { + if (k > 0) fprintf(fp,","); + fprintf(fp, "%f",x[ja[j]*dim+k]); + } + fprintf(fp, "}}"); + } + } + + fprintf(fp,"}],Hue[%f],",/*drand()*/1.); + + if (A->m < 100){ + for (i = 0; i < A->m; i++){ + if (i > 0) fprintf(fp,","); + fprintf(fp,"Text[%d,{",i+1); + for (k = 0; k < dim; k++) { + if (k > 0) fprintf(fp,","); + fprintf(fp, "%f",x[i*dim+k]); + } + fprintf(fp,"}]"); + } + } else if (A->m < 500000){ + fprintf(fp, "Point[{"); + for (i = 0; i < A->m; i++){ + if (i > 0) fprintf(fp,","); + fprintf(fp,"{"); + for (k = 0; k < dim; k++) { + if (k > 0) fprintf(fp,","); + fprintf(fp, "%f",x[i*dim+k]); + } + fprintf(fp,"}"); + } + fprintf(fp, "}]"); + } else { + 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"); + +} + +static real update_step(int adaptive_cooling, real step, real Fnorm, real Fnorm0, real cool){ + + if (!adaptive_cooling) { + return cool*step; + } + if (Fnorm >= Fnorm0){ + step = cool*step; + } else if (Fnorm > 0.95*Fnorm0){ + step = step; + } else { + step = 0.99*step/cool; + } + return step; +} + + +#define node_degree(i) (ia[(i)+1] - ia[(i)]) + +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); + } + +} +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); + } + +} + +real get_angle(real *x, int dim, int i, int j){ + /* between [0, 2Pi)*/ + int k; + real y[2], res; + real eps = 0.00001; + for (k = 0; k < 2; k++){ + 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; + } + res = atan(y[1]/y[0]); + if (y[0] > 0){ + if (y[1] < 0) res = 2*M_PI+res; + } else if (y[0] < 0){ + res = res + M_PI; + } + return res; +} + +int comp_real(const void *x, const void *y){ + real *xx = (real*) x; + real *yy = (real*) y; + + if (*xx > *yy){ + return 1; + } else if (*xx < *yy){ + return -1; + } + return 0; +} +static void sort_real(int n, real *a){ + qsort(a, n, sizeof(real), comp_real); +} + + +static void set_leaves(real *x, int dim, real dist, real ang, int i, int j){ + x[dim*j] = cos(ang)*dist + x[dim*i]; + x[dim*j+1] = sin(ang)*dist + x[dim*i+1]; +} + +static void beautify_leaves(int dim, SparseMatrix A, real *x){ + int m = A->m, i, j, *ia = A->ia, *ja = A->ja, k; + int *checked, p; + real dist; + int nleaves, nleaves_max = 10; + real *angles, maxang, ang1 = 0, ang2 = 0, pad, step; + int *leaves, nangles_max = 10, nangles; + + assert(!SparseMatrix_has_diagonal(A)); + + checked = N_GNEW(m,int); + angles = N_GNEW(nangles_max,real); + leaves = N_GNEW(nleaves_max,int); + + + for (i = 0; i < m; i++) checked[i] = FALSE; + + for (i = 0; i < m; i++){ + if (ia[i+1] - ia[i] != 1) continue; + if (checked[i]) continue; + p = ja[ia[i]]; + if (!checked[p]){ + checked[p] = TRUE; + dist = 0; nleaves = 0; nangles = 0; + for (j = ia[p]; j < ia[p+1]; j++){ + if (node_degree(ja[j]) == 1){ + checked[ja[j]] = TRUE; + check_int_array_size(&leaves, nleaves, &nleaves_max); + dist += distance(x, dim, p, ja[j]); + leaves[nleaves] = ja[j]; + nleaves++; + } else { + check_real_array_size(&angles, nangles, &nangles_max); + angles[nangles++] = get_angle(x, dim, p, ja[j]); + } + } + assert(nleaves > 0); + dist /= nleaves; + if (nangles > 0){ + sort_real(nangles, angles); + maxang = 0; + for (k = 0; k < nangles - 1; k++){ + if (angles[k+1] - angles[k] > maxang){ + maxang = angles[k+1] - angles[k]; + 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]; + ang1 = angles[nangles - 1]; + ang2 = 2*M_PI + angles[0]; + } + } else { + ang1 = 0; ang2 = 2*M_PI; maxang = 2*M_PI; + } + pad = MAX(maxang - M_PI*0.166667*(nleaves-1), 0)*0.5; + ang1 += pad*0.95; + ang2 -= pad*0.95; + assert(ang2 >= ang1); + step = 0.; + if (nleaves > 1) step = (ang2 - ang1)/(nleaves - 1); + for (i = 0; i < nleaves; i++) { + set_leaves(x, dim, dist, ang1, p, leaves[i]); + ang1 += step; + } + } + } + + + FREE(checked); + FREE(angles); + FREE(leaves); +} + +void force_print(FILE *fp, int n, int dim, real *x, real *force){ + int i, k; + + fprintf(fp,"Graphics[{"); + for (i = 0; i < n; i++){ + if (i > 0) fprintf(fp, ","); + fprintf(fp, "Arrow[{{"); + for (k = 0; k < dim; k++){ + if (k > 0) fprintf(fp, ","); + fprintf(fp, "%f",x[i*dim+k]); + } + fprintf(fp, "},{"); + for (k = 0; k < dim; k++){ + if (k > 0) fprintf(fp, ","); + fprintf(fp, "%f",x[i*dim+k]+0.5*force[i*dim+k]); + } + fprintf(fp, "}}]"); + } + fprintf(fp,","); + for (i = 0; i < n; i++){ + if (i > 0) fprintf(fp, ","); + fprintf(fp, "Tooltip[Point[{"); + for (k = 0; k < dim; k++){ + if (k > 0) fprintf(fp, ","); + fprintf(fp, "%f",x[i*dim+k]); + } + fprintf(fp, "}],%d]",i); + } + + + + + fprintf(fp,"}]\n"); + +} + + +void spring_electrical_embedding_fast(int dim, SparseMatrix A0, spring_electrical_control ctrl, real *node_weights, real *x, 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. */ + SparseMatrix A = A0; + int m, n; + int i, j, k; + real p = ctrl->p, K = ctrl->K, C = ctrl->C, CRK, tol = ctrl->tol, maxiter = ctrl->maxiter, cool = ctrl->cool, step = ctrl->step, KP; + int *ia = NULL, *ja = NULL; + real *xold = NULL; + real *f = NULL, dist, F, Fnorm = 0, Fnorm0; + int iter = 0; + int adaptive_cooling = ctrl->adaptive_cooling; + QuadTree qt = NULL; + real counts[4], *force = NULL; +#ifdef TIME + clock_t start, end, start0, start2; + real qtree_cpu = 0, qtree_cpu0 = 0, qtree_new_cpu = 0, qtree_new_cpu0 = 0; + real total_cpu = 0; + start0 = clock(); +#endif + int max_qtree_level = ctrl->max_qtree_level; + oned_optimizer qtree_level_optimizer = NULL; + + if (!A) return; + + m = A->m, n = A->n; + if (n <= 0 || dim <= 0) return; + + qtree_level_optimizer = oned_optimizer_new(max_qtree_level); + + *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; + + K = initPos (A, n, dim, x, ctrl); + + 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); + + do { +#ifdef TIME + start2 = clock(); +#endif + + iter++; + xold = MEMCPY(xold, x, sizeof(real)*dim*n); + Fnorm0 = Fnorm; + Fnorm = 0.; + + max_qtree_level = oned_optimizer_get(qtree_level_optimizer); + +#ifdef TIME + start = clock(); +#endif + 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); + } + +#ifdef TIME + qtree_new_cpu += ((real) (clock() - start))/CLOCKS_PER_SEC; +#endif + + /* repulsive force */ +#ifdef TIME + start = clock(); +#endif + + QuadTree_get_repulvice_force(qt, force, x, ctrl->bh, p, KP, counts, flag); + + assert(!(*flag)); + +#ifdef TIME + end = clock(); + qtree_cpu += ((real) (end - start)) / CLOCKS_PER_SEC; +#endif + + /* attractive force C^((2-p)/3) ||x_i-x_j||/K * (x_j - x_i) */ + for (i = 0; i < n; i++){ + f = &(force[i*dim]); + for (j = ia[i]; j < ia[i+1]; j++){ + if (ja[j] == i) continue; + dist = distance(x, dim, i, ja[j]); + for (k = 0; k < dim; k++){ + f[k] -= CRK*(x[i*dim+k] - x[ja[j]*dim+k])*dist; + } + } + } + + + /* move */ + for (i = 0; i < n; i++){ + f = &(force[i*dim]); + 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) { +#ifdef TIME + start = clock(); +#endif + QuadTree_delete(qt); +#ifdef TIME + end = clock(); + qtree_new_cpu += ((real) (end - start)) / CLOCKS_PER_SEC; +#endif + +#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", + 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 + oned_optimizer_train(qtree_level_optimizer, counts[0]+0.85*counts[1]+3.3*counts[2]); + } else { +#ifdef DEBUG_PRINT + if (Verbose) { + fprintf(stderr, "\r iter = %d, step = %f Fnorm = %f nz = %d K = %f ",iter, step, Fnorm, A->nz,K); +#ifdef ENERGY + fprintf(stderr, "energy = %f\n",spring_electrical_energy(dim, A, x, p, CRK, KP)); +#endif + } +#endif + } + + step = update_step(adaptive_cooling, step, Fnorm, Fnorm0, cool); + } while (step > tol && iter < maxiter); + +#ifdef DEBUG_PRINT + if (Verbose && 0) fputs("\n", stderr); +#endif + + +#ifdef DEBUG_PRINT + if (Verbose) { + fprintf(stderr, "\n iter = %d, step = %f Fnorm = %f nz = %d K = %f ",iter, step, Fnorm, A->nz, K); + } +#endif + + if (ctrl->beautify_leaves) beautify_leaves(dim, A, x); + +#ifdef TIME + total_cpu += ((real) (clock() - start0)) / CLOCKS_PER_SEC; + if (Verbose) fprintf(stderr, "\n time for qtree = %f, qtree_force = %f, total cpu = %f\n",qtree_new_cpu, qtree_cpu, total_cpu); +#endif + + + RETURN: + oned_optimizer_delete(qtree_level_optimizer); + ctrl->max_qtree_level = max_qtree_level; + + if (xold) FREE(xold); + if (A != A0) SparseMatrix_delete(A); + if (force) FREE(force); + +} + + +void spring_electrical_embedding_slow(int dim, SparseMatrix A0, spring_electrical_control ctrl, real *node_weights, real *x, int *flag){ + /* a version that does vertex moves in one go, instead of one at a time, use for debugging the fast version. Quadtree is not used. */ + /* x is a point to a 1D array, x[i*dim+j] gives the coordinate of the i-th node at dimension j. */ + SparseMatrix A = A0; + int m, n; + int i, j, k; + real p = ctrl->p, K = ctrl->K, C = ctrl->C, CRK, tol = ctrl->tol, maxiter = ctrl->maxiter, cool = ctrl->cool, step = ctrl->step, KP; + int *ia = NULL, *ja = NULL; + 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, counts_avg = 0; + real *force; +#ifdef TIME + clock_t start, end, start0, start2; + real qtree_cpu = 0, qtree_cpu0 = 0; + real total_cpu = 0; + start0 = clock(); +#endif + int max_qtree_level = ctrl->max_qtree_level; + oned_optimizer qtree_level_optimizer = NULL; + + fprintf(stderr,"spring_electrical_embedding_slow"); + if (!A) return; + + m = A->m, n = A->n; + if (n <= 0 || dim <= 0) return; + force = N_GNEW(n*dim,real); + + 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); + } + USE_QT = FALSE; + *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; + + K = initPos (A, n, dim, x, ctrl); + + 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; + +#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 = N_GNEW(dim,real); + xold = N_GNEW(dim*n,real); + do { + for (i = 0; i < dim*n; i++) force[i] = 0; + + iter++; + xold = MEMCPY(xold, x, sizeof(real)*dim*n); + Fnorm0 = Fnorm; + Fnorm = 0.; + nsuper_avg = 0; + + if (USE_QT) { + max_qtree_level = oned_optimizer_get(qtree_level_optimizer); + 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); + } + } +#ifdef TIME + start2 = clock(); +#endif + + + for (i = 0; i < n; i++){ + for (k = 0; k < dim; k++) f[k] = 0.; + /* repulsive force K^(1 - p)/||x_i-x_j||^(1 - p) (x_i - x_j) */ + if (USE_QT){ +#ifdef TIME + start = clock(); +#endif + QuadTree_get_supernodes(qt, ctrl->bh, &(x[dim*i]), i, &nsuper, &nsupermax, + ¢er, &supernode_wgts, &distances, &counts, flag); +#ifdef TIME + end = clock(); + qtree_cpu += ((real) (end - start)) / CLOCKS_PER_SEC; +#endif + counts_avg += counts; + 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] += supernode_wgts[j]*KP*(x[i*dim+k] - center[j*dim+k])/(dist*dist); + } else { + f[k] += supernode_wgts[j]*KP*(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] += node_weights[j]*KP*(x[i*dim+k] - x[j*dim+k])/(dist*dist); + } else { + f[k] += node_weights[j]*KP*(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] += KP*(x[i*dim+k] - x[j*dim+k])/(dist*dist); + } else { + f[k] += KP*(x[i*dim+k] - x[j*dim+k])/pow(dist, 1.- p); + } + } + } + } + } + for (k = 0; k < dim; k++) force[i*dim+k] += f[k]; + } + + + + for (i = 0; i < n; i++){ + for (k = 0; k < dim; k++) f[k] = 0.; + /* attractive force C^((2-p)/3) ||x_i-x_j||/K * (x_j - x_i) */ + for (j = ia[i]; j < ia[i+1]; j++){ + if (ja[j] == i) continue; + dist = distance(x, dim, i, ja[j]); + for (k = 0; k < dim; k++){ + f[k] -= CRK*(x[i*dim+k] - x[ja[j]*dim+k])*dist; + } + } + for (k = 0; k < dim; k++) force[i*dim+k] += f[k]; + } + + + + for (i = 0; i < n; i++){ + /* normalize force */ + for (k = 0; k < dim; k++) f[k] = force[i*dim+k]; + + 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; + counts_avg /= n; +#ifdef TIME + qtree_cpu0 = qtree_cpu - qtree_cpu0; + if (Verbose && 0) fprintf(stderr, "\n cpu this outer iter = %f, quadtree time = %f other time = %f\n",((real) (clock() - start2)) / CLOCKS_PER_SEC, qtree_cpu0,((real) (clock() - start2))/CLOCKS_PER_SEC - qtree_cpu0); + qtree_cpu0 = qtree_cpu; +#endif + if (Verbose && 0) fprintf(stderr, "nsuper_avg=%f, counts_avg = %f 2*nsuper+counts=%f\n",nsuper_avg,counts_avg, 2*nsuper_avg+counts_avg); + oned_optimizer_train(qtree_level_optimizer, 5*nsuper_avg + counts_avg); + } + +#ifdef ENERGY + if (Verbose) { + fprintf(stderr, "\r iter = %d, step = %f Fnorm = %f nsuper = %d nz = %d K = %f ",iter, step, Fnorm, (int) nsuper_avg,A->nz,K); + fprintf(stderr, "energy = %f\n",spring_electrical_energy(dim, A, x, p, CRK, KP)); + } +#endif + + + step = update_step(adaptive_cooling, step, Fnorm, Fnorm0, cool); + } while (step > tol && iter < maxiter); + +#ifdef DEBUG_PRINT + if (Verbose && 0) fputs("\n", stderr); +#endif + +#ifdef DEBUG_PRINT_0 + { + FILE *f; + char fname[10000]; + strcpy(fname,"/tmp/graph_layout"); + sprintf(&(fname[strlen(fname)]), "%d",n); + f = fopen(fname,"w"); + export_embedding(f, dim, A, x, NULL); + fclose(f); + } +#endif + + +#ifdef DEBUG_PRINT + if (Verbose) { + if (USE_QT){ + fprintf(stderr, "iter = %d, step = %f Fnorm = %f qt_level = %d nsuper = %d nz = %d K = %f ",iter, step, Fnorm, max_qtree_level, (int) nsuper_avg,A->nz,K); + } else { + fprintf(stderr, "iter = %d, step = %f Fnorm = %f nsuper = %d nz = %d K = %f ",iter, step, Fnorm, (int) nsuper_avg,A->nz,K); + } + } +#endif + + if (ctrl->beautify_leaves) beautify_leaves(dim, A, x); + +#ifdef TIME + total_cpu += ((real) (clock() - start0)) / CLOCKS_PER_SEC; + if (Verbose) fprintf(stderr, "time for supernode = %f, total cpu = %f\n",qtree_cpu, total_cpu); +#endif + + RETURN: + if (USE_QT) { + oned_optimizer_delete(qtree_level_optimizer); + ctrl->max_qtree_level = max_qtree_level; + } + 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); + FREE(force); + +} + + + +void spring_electrical_embedding(int dim, SparseMatrix A0, spring_electrical_control ctrl, real *node_weights, real *x, 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. */ + SparseMatrix A = A0; + int m, n; + int i, j, k; + real p = ctrl->p, K = ctrl->K, C = ctrl->C, CRK, tol = ctrl->tol, maxiter = ctrl->maxiter, cool = ctrl->cool, step = ctrl->step, KP; + int *ia = NULL, *ja = NULL; + 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, counts_avg = 0; +#ifdef TIME + clock_t start, end, start0, start2; + real qtree_cpu = 0, qtree_cpu0 = 0; + real total_cpu = 0; + start0 = clock(); +#endif + int max_qtree_level = ctrl->max_qtree_level; + oned_optimizer qtree_level_optimizer = NULL; + + if (!A) return; + + m = A->m, n = A->n; + if (n <= 0 || dim <= 0) return; + + 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); + } + *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; + + K = initPos (A, n, dim, x, ctrl); + + 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; + +#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 = N_GNEW(dim,real); + xold = N_GNEW(dim*n,real); + do { + iter++; + xold = MEMCPY(xold, x, sizeof(real)*dim*n); + Fnorm0 = Fnorm; + Fnorm = 0.; + nsuper_avg = 0; + counts_avg = 0; + + if (USE_QT) { + + max_qtree_level = oned_optimizer_get(qtree_level_optimizer); + 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); + } + + + } +#ifdef TIME + start2 = clock(); +#endif + + for (i = 0; i < n; i++){ + for (k = 0; k < dim; k++) f[k] = 0.; + /* attractive force C^((2-p)/3) ||x_i-x_j||/K * (x_j - x_i) */ + for (j = ia[i]; j < ia[i+1]; j++){ + if (ja[j] == i) continue; + dist = distance(x, dim, i, ja[j]); + for (k = 0; k < dim; k++){ + f[k] -= CRK*(x[i*dim+k] - x[ja[j]*dim+k])*dist; + } + } + + /* repulsive force K^(1 - p)/||x_i-x_j||^(1 - p) (x_i - x_j) */ + if (USE_QT){ +#ifdef TIME + start = clock(); +#endif + QuadTree_get_supernodes(qt, ctrl->bh, &(x[dim*i]), i, &nsuper, &nsupermax, + ¢er, &supernode_wgts, &distances, &counts, flag); + +#ifdef TIME + end = clock(); + qtree_cpu += ((real) (end - start)) / CLOCKS_PER_SEC; +#endif + counts_avg += counts; + 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] += supernode_wgts[j]*KP*(x[i*dim+k] - center[j*dim+k])/(dist*dist); + } else { + f[k] += supernode_wgts[j]*KP*(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] += node_weights[j]*KP*(x[i*dim+k] - x[j*dim+k])/(dist*dist); + } else { + f[k] += node_weights[j]*KP*(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] += KP*(x[i*dim+k] - x[j*dim+k])/(dist*dist); + } else { + f[k] += KP*(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; + counts_avg /= n; +#ifdef TIME + qtree_cpu0 = qtree_cpu - qtree_cpu0; + if (Verbose && 0) fprintf(stderr, "\n cpu this outer iter = %f, quadtree time = %f other time = %f\n",((real) (clock() - start2)) / CLOCKS_PER_SEC, qtree_cpu0,((real) (clock() - start2))/CLOCKS_PER_SEC - qtree_cpu0); + qtree_cpu0 = qtree_cpu; +#endif + if (Verbose & 0) fprintf(stderr, "nsuper_avg=%f, counts_avg = %f 2*nsuper+counts=%f\n",nsuper_avg,counts_avg, 2*nsuper_avg+counts_avg); + oned_optimizer_train(qtree_level_optimizer, 5*nsuper_avg + counts_avg); + } + +#ifdef ENERGY + if (Verbose) { + fprintf(stderr, "\r iter = %d, step = %f Fnorm = %f nsuper = %d nz = %d K = %f ",iter, step, Fnorm, (int) nsuper_avg,A->nz,K); + fprintf(stderr, "energy = %f\n",spring_electrical_energy(dim, A, x, p, CRK, KP)); + } +#endif + + + step = update_step(adaptive_cooling, step, Fnorm, Fnorm0, cool); + } while (step > tol && iter < maxiter); + +#ifdef DEBUG_PRINT + if (Verbose && 0) fputs("\n", stderr); +#endif + +#ifdef DEBUG_PRINT_0 + { + FILE *f; + char fname[10000]; + strcpy(fname,"/tmp/graph_layout"); + sprintf(&(fname[strlen(fname)]), "%d",n); + f = fopen(fname,"w"); + export_embedding(f, dim, A, x, NULL); + fclose(f); + } +#endif + + +#ifdef DEBUG_PRINT + if (Verbose) { + if (USE_QT){ + fprintf(stderr, "iter = %d, step = %f Fnorm = %f qt_level = %d nsuper = %d nz = %d K = %f ",iter, step, Fnorm, max_qtree_level, (int) nsuper_avg,A->nz,K); + } else { + fprintf(stderr, "iter = %d, step = %f Fnorm = %f nsuper = %d nz = %d K = %f ",iter, step, Fnorm, (int) nsuper_avg,A->nz,K); + } + } +#endif + + if (ctrl->beautify_leaves) beautify_leaves(dim, A, x); + +#ifdef TIME + total_cpu += ((real) (clock() - start0)) / CLOCKS_PER_SEC; + if (Verbose) fprintf(stderr, "time for supernode = %f, total cpu = %f\n",qtree_cpu, total_cpu); +#endif + + RETURN: + if (USE_QT) { + oned_optimizer_delete(qtree_level_optimizer); + ctrl->max_qtree_level = max_qtree_level; + } + 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); + +} + + + + +void spring_electrical_spring_embedding(int dim, SparseMatrix A0, SparseMatrix D, spring_electrical_control ctrl, real *node_weights, real *x, 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. Same as the spring-electrical except we also + introduce force due to spring length + */ + SparseMatrix A = A0; + int m, n; + int i, j, k; + real p = ctrl->p, K = ctrl->K, C = ctrl->C, CRK, tol = ctrl->tol, maxiter = ctrl->maxiter, cool = ctrl->cool, step = ctrl->step, KP; + int *ia = NULL, *ja = NULL; + int *id = NULL, *jd = NULL; + real *d; + 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; + + if (!A) return; + m = A->m, n = A->n; + if (n <= 0 || dim <= 0) return; + + 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); + } + *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; + id = D->ia; + jd = D->ja; + d = (real*) D->a; + + K = initPos (A, n, dim, x, ctrl); + + 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; + +#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 = N_GNEW(dim,real); + xold = N_GNEW(dim*n,real); + do { + iter++; + xold = MEMCPY(xold, x, sizeof(real)*dim*n); + Fnorm0 = Fnorm; + Fnorm = 0.; + nsuper_avg =- 0; + + 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); + } + } + + for (i = 0; i < n; i++){ + for (k = 0; k < dim; k++) f[k] = 0.; + /* attractive force C^((2-p)/3) ||x_i-x_j||/K * (x_j - x_i) */ + + for (j = ia[i]; j < ia[i+1]; j++){ + if (ja[j] == i) continue; + dist = distance(x, dim, i, ja[j]); + for (k = 0; k < dim; k++){ + f[k] -= CRK*(x[i*dim+k] - x[ja[j]*dim+k])*dist; + } + } + + for (j = id[i]; j < id[i+1]; j++){ + if (jd[j] == i) continue; + dist = distance_cropped(x, dim, i, jd[j]); + for (k = 0; k < dim; k++){ + if (dist < d[j]){ + f[k] += 0.2*CRK*(x[i*dim+k] - x[jd[j]*dim+k])*(dist - d[j])*(dist - d[j])/dist; + } else { + f[k] -= 0.2*CRK*(x[i*dim+k] - x[jd[j]*dim+k])*(dist - d[j])*(dist - d[j])/dist; + } + /* f[k] -= 0.2*CRK*(x[i*dim+k] - x[jd[j]*dim+k])*(dist - d[j]);*/ + } + } + + /* repulsive force K^(1 - p)/||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] += supernode_wgts[j]*KP*(x[i*dim+k] - center[j*dim+k])/(dist*dist); + } else { + f[k] += supernode_wgts[j]*KP*(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] += node_weights[j]*KP*(x[i*dim+k] - x[j*dim+k])/(dist*dist); + } else { + f[k] += node_weights[j]*KP*(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] += KP*(x[i*dim+k] - x[j*dim+k])/(dist*dist); + } else { + f[k] += KP*(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 + if (Verbose && 0) { + fprintf(stderr, "\r iter = %d, step = %f Fnorm = %f nsuper = %d nz = %d K = %f ",iter, step, Fnorm, (int) nsuper_avg,A->nz,K); +#ifdef ENERGY + fprintf(stderr, "energy = %f\n",spring_electrical_energy(dim, A, x, p, CRK, KP)); +#endif + } +#endif + + step = update_step(adaptive_cooling, step, Fnorm, Fnorm0, cool); + } while (step > tol && iter < maxiter); + +#ifdef DEBUG_PRINT + if (Verbose && 0) fputs("\n", stderr); +#endif + +#ifdef DEBUG_PRINT_0 + { + FILE *f; + char fname[10000]; + strcpy(fname,"/tmp/graph_layout"); + sprintf(&(fname[strlen(fname)]), "%d",n); + f = fopen(fname,"w"); + export_embedding(f, dim, A, x, NULL); + fclose(f); + } +#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); + +} + + + + +void print_matrix(real *x, int n, int dim){ + int i, k; + printf("{"); + for (i = 0; i < n; i++){ + if (i != 0) printf(","); + printf("{"); + for (k = 0; k < dim; k++) { + if (k != 0) printf(","); + printf("%f",x[i*dim+k]); + } + printf("}"); + } + printf("}\n"); +} + +/* +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); + for (k = 0; k < dim*m; k++) y[k] = 0; + for (i = 0; i < m; i++){ + nz = 0; + for (j = ia[i]; j < ia[i+1]; j++){ + if (ja[j] == i) continue; + nz++; + for (k = 0; k < dim; k++){ + y[i*dim+k] += x[ja[j]*dim + k]; + } + } + if (nz > 0){ + beta = (1-alpha)/nz; + for (k = 0; k < dim; k++) y[i*dim+k] = alpha*x[i*dim+k] + beta*y[i*dim+k]; + } + } + for (k = 0; k < dim*m; k++) x[k] = y[k]; + + FREE(y); +} + +*/ + +static void interpolate0(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); + for (i = 0; i < A->m; i++){ + for (k = 0; k < dim; k++) y[k] = 0; + nz = 0; + for (j = ia[i]; j < ia[i+1]; j++){ + if (ja[j] == i) continue; + nz++; + for (k = 0; k < dim; k++){ + y[k] += x[ja[j]*dim + k]; + } + } + if (nz > 0){ + beta = (1-alpha)/nz; + for (k = 0; k < dim; k++) x[i*dim+k] = alpha*x[i*dim+k] + beta*y[k]; + } + } + + FREE(y); +} +static void prolongate(int dim, SparseMatrix A, SparseMatrix P, SparseMatrix R, real *x, real *y, int coarsen_scheme_used, real delta){ + int nc, *ia, *ja, i, j, k; + SparseMatrix_multiply_dense(P, FALSE, x, FALSE, &y, FALSE, dim); + + /* xu yao rao dong */ + if (coarsen_scheme_used > EDGE_BASED_STA && coarsen_scheme_used < EDGE_BASED_STO){ + interpolate0(dim, A, y); + nc = R->m; + ia = R->ia; + ja = R->ja; + for (i = 0; i < nc; i++){ + for (j = ia[i]+1; j < ia[i+1]; j++){ + for (k = 0; k < dim; k++){ + y[ja[j]*dim + k] += delta*(drand() - 0.5); + } + } + } + } +} + + + +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); + + for (i = 0; i < m + 1; i++){ + mask[i] = 0; + } + + for (i = 0; i < m; i++){ + deg = 0; + for (j = ia[i]; j < ia[i+1]; j++){ + if (i == ja[j]) continue; + deg++; + } + mask[deg]++; + max = MAX(max, mask[deg]); + } + if (mask[1] > 0.8*max && mask[1] > 0.3*m) res = TRUE; + FREE(mask); + return res; +} + +void pcp_rotate(int n, int dim, real *x){ + int i, k,l; + real y[4], axis[2], center[2], dist, x0, x1; + + 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]; + } + } + + for (i = 0; i < n; i++){ + for (k = 0; k < dim; k++){ + for (l = 0; l < dim; l++){ + y[dim*k+l] += x[i*dim+k]*x[i*dim+l]; + } + } + } + if (y[1] == 0) { + axis[0] = 0; axis[1] = 1; + } else { + /* Eigensystem[{{x0, x1}, {x1, x3}}] = + {{(x0 + x3 - Sqrt[x0^2 + 4*x1^2 - 2*x0*x3 + x3^2])/2, + (x0 + x3 + Sqrt[x0^2 + 4*x1^2 - 2*x0*x3 + x3^2])/2}, + {{-(-x0 + x3 + Sqrt[x0^2 + 4*x1^2 - 2*x0*x3 + x3^2])/(2*x1), 1}, + {-(-x0 + x3 - Sqrt[x0^2 + 4*x1^2 - 2*x0*x3 + x3^2])/(2*x1), 1}}} + */ + axis[0] = -(-y[0] + y[3] - sqrt(y[0]*y[0]+4*y[1]*y[1]-2*y[0]*y[3]+y[3]*y[3]))/(2*y[1]); + axis[1] = 1; + } + dist = sqrt(1+axis[0]*axis[0]); + axis[0] = axis[0]/dist; + axis[1] = axis[1]/dist; + 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; + + } + + +} + + +void multilevel_spring_electrical_embedding(int dim, SparseMatrix A0, spring_electrical_control ctrl, real *node_weights, real *label_sizes, + real *x, int *flag){ + + + Multilevel_control mctrl = NULL; + int n, plg, coarsen_scheme_used; + SparseMatrix A = A0, P = NULL; + Multilevel grid, grid0; + real *xc = NULL, *xf = NULL; +#ifdef TIME + clock_t cpu; +#endif + struct spring_electrical_control_struct ctrl0; + + ctrl0 = *ctrl; + +#ifdef TIME + cpu = clock(); +#endif + + *flag = 0; + if (!A) return; + n = A->n; + 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); + } else { + A = SparseMatrix_remove_diagonal(A); + } + + mctrl = Multilevel_control_new(); + mctrl->maxlevel = ctrl->multilevels; + grid0 = Multilevel_new(A, node_weights, mctrl); + + grid = Multilevel_get_coarsest(grid0); + if (Multilevel_is_finest(grid)){ + xc = x; + } else { + xc = N_GNEW(grid->n*dim,real); + } + + plg = power_law_graph(A); + if (ctrl->p == AUTOP){ + ctrl->p = -1; + if (plg) ctrl->p = -1.8; + } + + do { +#ifdef DEBUG_PRINT + if (Verbose) { + print_padding(grid->level); + if (Multilevel_is_coarsest(grid)){ + fprintf(stderr, "coarsest level -- %d, n = %d\n", grid->level, grid->n); + } else { + fprintf(stderr, "level -- %d, n = %d\n", grid->level, grid->n); + } + } +#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); + } else { + assert(0); + } + if (Multilevel_is_finest(grid)) break; + if (*flag) { + FREE(xc); + goto RETURN; + } + P = grid->P; + coarsen_scheme_used = grid->coarsen_scheme_used; + grid = grid->prev; + if (Multilevel_is_finest(grid)){ + xf = x; + } else { + xf = N_GNEW(grid->n*dim,real); + } + prolongate(dim, grid->A, P, grid->R, xc, xf, coarsen_scheme_used, (ctrl->K)*0.001); + FREE(xc); + xc = xf; + ctrl->random_start = FALSE; + ctrl->K = ctrl->K * 0.75; + ctrl->adaptive_cooling = FALSE; + if (grid->next->coarsen_scheme_used > VERTEX_BASED_STA && + grid->next->coarsen_scheme_used < VERTEX_BASED_STO){ + ctrl->step = 1; + } else { + ctrl->step = .1; + } + } while (grid); + +#ifdef TIME + if (Verbose) + fprintf(stderr, "layout time %f\n",((real) (clock() - cpu)) / CLOCKS_PER_SEC); + cpu = clock(); +#endif + + post_process_smoothing(dim, A, ctrl, node_weights, x, flag); + + if (dim == 2){ + pcp_rotate(n, dim, x); + } + + if (Verbose) + fprintf(stderr, "sfdp: overlap=%d scaling %.02f\n", + ctrl->overlap, ctrl->initial_scaling); + + if (ctrl->overlap >= 0) + remove_overlap(dim, A, A->m, x, label_sizes, ctrl->overlap, ctrl->initial_scaling, flag); + + RETURN: + *ctrl = ctrl0; + if (A != A0) SparseMatrix_delete(A); + Multilevel_control_delete(mctrl); + Multilevel_delete(grid0); +} + diff --git a/lib/sfdpgen/spring_electrical.h b/lib/sfdpgen/spring_electrical.h new file mode 100644 index 000000000..bccf1658c --- /dev/null +++ b/lib/sfdpgen/spring_electrical.h @@ -0,0 +1,80 @@ +/* $Id$ $Revision$ */ +/* vim:set shiftwidth=4 ts=8: */ + +/********************************************************** +* This software is part of the graphviz package * +* http://www.graphviz.org/ * +* * +* Copyright (c) 1994-2004 AT&T Corp. * +* and is licensed under the * +* Common Public License, Version 1.0 * +* by AT&T Corp. * +* * +* Information and Software Systems Research * +* AT&T Research, Florham Park NJ * +**********************************************************/ + +#ifndef SPRING_ELECTRICAL_H +#define SPRING_ELECTRICAL_H + +#include + +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}; + +struct spring_electrical_control_struct { + real p;/*a negativve number default to -1. repulsive force = dist^p */ + 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 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*/ + real tol;/* minimum different between two subsequence config before terminating. ||x-xold|| < tol */ + int maxiter; + real cool;/* default 0.9 */ + real step;/* initial step size */ + int adaptive_cooling; + int random_seed; + int beautify_leaves; + int use_node_weights; + int smoothing; + int overlap; + int tscheme; /* octree scheme. 0 (no octree), 1 (normal), 2 (fast) */ + char* posSet; + 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. + */ +}; + +typedef struct spring_electrical_control_struct *spring_electrical_control; + +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 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); +void force_print(FILE *fp, int n, int dim, real *x, real *force); + +#endif -- 2.40.0