From 0cd264937e210404c56c195b9b74749e41ffa332 Mon Sep 17 00:00:00 2001 From: "Emden R. Gansner" Date: Wed, 14 Aug 2013 17:07:33 -0400 Subject: [PATCH] Add mingle program and library to graphviz --- cmd/Makefile.am | 2 +- cmd/mingle/Makefile.am | 40 ++ cmd/mingle/mingle.vcproj | 203 ++++++ cmd/mingle/minglemain.c | 198 ++++++ configure.ac | 25 + graphviz.sln | 23 +- lib/Makefile.am | 2 +- lib/mingle/Makefile.am | 21 + lib/mingle/agglomerative_bundling.c | 773 ++++++++++++++++++++ lib/mingle/agglomerative_bundling.h | 34 + lib/mingle/edge_bundling.c | 831 ++++++++++++++++++++++ lib/mingle/edge_bundling.h | 45 ++ lib/mingle/ink.c | 357 ++++++++++ lib/mingle/ink.h | 37 + lib/mingle/minglelib.vcproj | 213 ++++++ lib/mingle/nearest_neighbor_graph.c | 56 ++ lib/mingle/nearest_neighbor_graph.h | 16 + lib/mingle/nearest_neighbor_graph_ann.cpp | 182 +++++ lib/mingle/nearest_neighbor_graph_ann.h | 16 + 19 files changed, 3071 insertions(+), 3 deletions(-) create mode 100644 cmd/mingle/Makefile.am create mode 100644 cmd/mingle/mingle.vcproj create mode 100644 cmd/mingle/minglemain.c create mode 100644 lib/mingle/Makefile.am create mode 100644 lib/mingle/agglomerative_bundling.c create mode 100644 lib/mingle/agglomerative_bundling.h create mode 100644 lib/mingle/edge_bundling.c create mode 100644 lib/mingle/edge_bundling.h create mode 100644 lib/mingle/ink.c create mode 100644 lib/mingle/ink.h create mode 100644 lib/mingle/minglelib.vcproj create mode 100644 lib/mingle/nearest_neighbor_graph.c create mode 100644 lib/mingle/nearest_neighbor_graph.h create mode 100644 lib/mingle/nearest_neighbor_graph_ann.cpp create mode 100644 lib/mingle/nearest_neighbor_graph_ann.h diff --git a/cmd/Makefile.am b/cmd/Makefile.am index 5b0a1dece..01c6472eb 100644 --- a/cmd/Makefile.am +++ b/cmd/Makefile.am @@ -1,6 +1,6 @@ # $Id$ $Revision$ ## Process this file with automake to produce Makefile.in -SUBDIRS = dot tools gvpr lefty lneato dotty smyrna gvmap gvedit +SUBDIRS = dot tools gvpr lefty lneato dotty smyrna gvmap gvedit mingle EXTRA_DIST = Makefile.old diff --git a/cmd/mingle/Makefile.am b/cmd/mingle/Makefile.am new file mode 100644 index 000000000..2ffa7cf8e --- /dev/null +++ b/cmd/mingle/Makefile.am @@ -0,0 +1,40 @@ +# $Id$ $Revision$ +## Process this file with automake to produce Makefile.in + +pdfdir = $(pkgdatadir)/doc/pdf + +AM_CPPFLAGS = \ + -I$(top_srcdir) \ + -I$(top_srcdir)/lib/sparse \ + -I$(top_srcdir)/lib/sfdpgen \ + -I$(top_srcdir)/lib/mingle \ + -I$(top_srcdir)/lib/ingraphs \ + -I$(top_srcdir)/lib/cgraph \ + -I$(top_srcdir)/lib/cdt + +bin_PROGRAMS = mingle + +#man_MANS = mingle.1 +#pdf_DATA = mingle.1.pdf +man_MANS = +pdf_DATA = + +mingle_SOURCES = minglemain.c +mingle_CPPFLAGS = $(AM_CPPFLAGS) +mingle_LDADD = \ + $(top_builddir)/lib/mingle/libmingle_C.la \ + $(top_builddir)/lib/ingraphs/libingraphs_C.la \ + $(top_builddir)/lib/sfdpgen/libsfdpgen_C.la \ + $(top_builddir)/lib/neatogen/libneatogen_C.la \ + $(top_builddir)/lib/sparse/libsparse_C.la \ + $(top_builddir)/lib/rbtree/librbtree_C.la \ + $(top_builddir)/lib/common/libcommon_C.la \ + $(top_builddir)/lib/cgraph/libcgraph.la \ + $(ANN_LIBS) -lstdc++ -lm + +mingle.1.pdf: $(srcdir)/mingle.1 + - @GROFF@ -Tps -man -t $(srcdir)/mingle.1 | @PS2PDF@ - - >mingle.1.pdf + +EXTRA_DIST = $(man_MANS) $(pdf_DATA) mingle.vcproj + +DISTCLEANFILES = $(pdf_DATA) diff --git a/cmd/mingle/mingle.vcproj b/cmd/mingle/mingle.vcproj new file mode 100644 index 000000000..6a8fea602 --- /dev/null +++ b/cmd/mingle/mingle.vcproj @@ -0,0 +1,203 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/cmd/mingle/minglemain.c b/cmd/mingle/minglemain.c new file mode 100644 index 000000000..3e0dae064 --- /dev/null +++ b/cmd/mingle/minglemain.c @@ -0,0 +1,198 @@ + +/************************************************************************* + * Copyright (c) 2011 AT&T Intellectual Property + * All rights reserved. This program and the accompanying materials + * are made available under the terms of the Eclipse Public License v1.0 + * which accompanies this distribution, and is available at + * http://www.eclipse.org/legal/epl-v10.html + * + *************************************************************************/ + +#ifdef HAVE_CONFIG_H +#include "config.h" +#endif + +#ifdef WIN32 /*dependencies*/ + #pragma comment( lib, "cgraph.lib" ) + #pragma comment( lib, "ingraphs.lib" ) +#endif + +#include +#include +#ifdef HAVE_GETOPT_H +#include +#else +#include "compat_getopt.h" +#endif + +#include "DotIO.h" +#include "edge_bundling.h" +#include "nearest_neighbor_graph.h" + +typedef struct { + int outer_iter; + int method; + int compatibility_method; + real K; + char* fmt; + int nneighbors; + int max_recursion; + real angle_param; + real angle; +} opts_t; + +static char *fname; +static FILE *outfile; + +static char* use_msg = +"Usage: mingle \n\ + -a t - max. turning angle [0-180] (40)\n\ + -c i - compatability measure; 0 : distance, 1: full (default)\n\ + -i iter: number of outer iterations/subdivisions (4)\n\ + -k k - number of neighbors in the nearest neighbor graph of edges (10)\n\ + -K k - the force constant\n\ + -m method - method used. 0 (force directed), 1 (agglomerative ink saving, default), 2 (cluster+ink saving)\n\ + -o fname - write output to file fname (stdout)\n\ + -p t - balance for avoiding sharp angles\n\ + The larger the t, the more sharp angles are allowed\n\ + -r R - max. recursion level with agglomerative ink saving method (100)\n\ + -T t - output format\n\ + -v - verbose\n"; + +static void +usage (int eval) +{ + fputs (use_msg, stderr); + if (eval >= 0) exit (eval); +} + +static void init(int argc, char *argv[], opts_t* opts) +{ + unsigned int c; + char* CmdName = argv[0]; + + opts->outer_iter = 4; +#ifdef HAVE_ANN + opts->method = METHOD_INK_AGGLOMERATE; +#else + opts->method = METHOD_FD; +#endif + opts->compatibility_method = COMPATIBILITY_FULL; + opts->K = -1; + opts->fmt = "gv"; + opts->nneighbors = 10; + opts->max_recursion = 100; + opts->angle_param = -1; + opts->angle = 40/180*M_PI; + + while ((c = getopt(argc, argv, ":a:c:i:k:K:m:o:p:r:T:v")) != -1) { + switch (c) { + case 'v': + Verbose = 1; + break; + case ':': + fprintf(stderr, "%s: option -%c missing argument\n", CmdName, optopt); + usage(1); + break; + case '?': + if (optopt == '?') + usage(0); + else + fprintf(stderr, "%s: option -%c unrecognized - ignored\n", CmdName, optopt); + break; + } + } + argv += optind; + argc -= optind; + + if (argc > 0) + Files = argv; + outfile = stdout; +} + +static int +bundle (Agraph_t* g, opts_t* opts) +{ + real *x = NULL; + real *label_sizes = NULL; + int n_edge_label_nodes; + int dim = 2; + SparseMatrix A; + SparseMatrix B; + pedge* edges; + real *xx, eps = 0.; + int nz = 0; + int *ia, *ja, i, j; + int rv = 0; + + initDotIO(g); + A = SparseMatrix_import_dot(g, dim, &label_sizes, &x, &n_edge_label_nodes, NULL, FORMAT_CSR, NULL); + if (!A){ + fprintf (stderr, "Error: could not convert graph %s (%s) into matrix\n", agnameof(g), fname); + return 1; + } + + A = SparseMatrix_symmetrize(A, TRUE); + ia = A->ia; ja = A->ja; + nz = A->nz; + xx = MALLOC(sizeof(real)*nz*4); + nz = 0; + dim = 4; + for (i = 0; i < A->m; i++){ + for (j = ia[i]; j < ia[i+1]; j++){ + if (ja[j] > i){ + xx[nz*dim] = x[i*2]; + xx[nz*dim+1] = x[i*2+1]; + xx[nz*dim+2] = x[ja[j]*2]; + xx[nz*dim+3] = x[ja[j]*2+1]; + nz++; + } + } + } + if (Verbose) + fprintf(stderr,"n = %d nz = %d\n",A->m, nz); + + B = nearest_neighbor_graph(nz, MIN(opts->nneighbors, nz), dim, xx, eps); + + SparseMatrix_delete(A); + A = B; + FREE(x); + x = xx; + + dim = 2; + + edges = edge_bundling(A, 2, x, opts->outer_iter, opts->K, opts->method, opts->nneighbors, opts->compatibility_method, opts->max_recursion, opts->angle_param, opts->angle, 0); + + pedge_export_gv(stdout, A->m, edges); + return rv; +} + +static Agraph_t *gread(FILE * fp) +{ + return agread(fp, (Agdisc_t *) 0); +} + +int main(int argc, char *argv[]) +{ + Agraph_t *g; + Agraph_t *prev = NULL; + ingraph_state ig; + int rv = 0; + opts_t opts; + + init(argc, argv, &opts); + newIngraph(&ig, Files, gread); + + while ((g = nextGraph(&ig)) != 0) { + if (prev) + agclose(prev); + prev = g; + fname = fileName(&ig); + if (Verbose) + fprintf(stderr, "Process graph %s in file %s\n", agnameof(g), + fname); + rv |= bundle (g, &opts); + } + + return rv; +} diff --git a/configure.ac b/configure.ac index 9b78ffe96..d5b13185e 100644 --- a/configure.ac +++ b/configure.ac @@ -2474,6 +2474,28 @@ else fi AM_CONDITIONAL(WITH_GTS, [test "x$use_gts" == "xYes"]) +dnl ----------------------------------- +dnl INCLUDES and LIBS for ANN. + +AC_ARG_WITH(ann, + [AS_HELP_STRING([--with-ann=yes],[ANN library])], + [], [with_ann=yes]) + +if test "x$with_ann" != "xyes"; then + use_ann="No (disabled)" +else + PKG_CHECK_MODULES(ANN, [ann],[ + use_ann="Yes" + AC_DEFINE_UNQUOTED(HAVE_ANN,1, + [Define if you have the ann library]) + AC_SUBST([ANN_CFLAGS]) + AC_SUBST([ANN_LIBS]) + ],[ + use_ann="No (ANN library not available)" + ]) +fi +AM_CONDITIONAL(WITH_ANN, [test "x$use_ann" == "xYes"]) + dnl ----------------------------------- dnl INCLUDES and LIBS for GLADE. @@ -3114,6 +3136,7 @@ AC_CONFIG_FILES(Makefile lib/neatogen/Makefile lib/fdpgen/Makefile lib/sparse/Makefile + lib/mingle/Makefile lib/label/Makefile lib/sfdpgen/Makefile lib/osage/Makefile @@ -3179,6 +3202,7 @@ AC_CONFIG_FILES(Makefile cmd/dotty/Makefile cmd/smyrna/Makefile cmd/gvmap/Makefile + cmd/mingle/Makefile cmd/gvedit/Makefile cmd/gvedit/gvedit.pro cmd/gvedit/ui/Makefile @@ -3254,6 +3278,7 @@ echo " expat: $use_expat" echo " fontconfig: $use_fontconfig" echo " freetype: $use_freetype" echo " glut: $use_glut" +echo " ann: $use_ann" echo " gts: $use_gts" echo " ipsepcola: $use_ipsepcola" echo " ltdl: $use_ltdl" diff --git a/graphviz.sln b/graphviz.sln index c2088ec69..357739505 100644 --- a/graphviz.sln +++ b/graphviz.sln @@ -1,5 +1,4 @@ Microsoft Visual Studio Solution File, Format Version 10.00 -# # Visual Studio 2008 Project("{8BC9CEB8-8B4A-11D0-8D11-00A0C91BC942}") = "Pathplan", "lib\pathplan\Pathplan.vcproj", "{BD347753-A09D-48B4-8752-F1D8D9CF235D}" ProjectSection(ProjectDependencies) = postProject @@ -262,6 +261,20 @@ Project("{8BC9CEB8-8B4A-11D0-8D11-00A0C91BC942}") = "gvedit", "cmd\gvedit\gvedit {15229511-9F6C-48A5-9194-660CA6492563} = {15229511-9F6C-48A5-9194-660CA6492563} EndProjectSection EndProject +Project("{8BC9CEB8-8B4A-11D0-8D11-00A0C91BC942}") = "minglelib", "lib\mingle\minglelib.vcproj", "{70575BD2-A532-41B8-9A17-B9876E992A84}" +EndProject +Project("{8BC9CEB8-8B4A-11D0-8D11-00A0C91BC942}") = "mingle", "cmd\mingle\mingle.vcproj", "{A6DF0D74-E4D1-4EF4-A642-59B0CF3E74BB}" + ProjectSection(ProjectDependencies) = postProject + {C0663A08-F276-4DD6-B17C-E501EE066F7C} = {C0663A08-F276-4DD6-B17C-E501EE066F7C} + {15229511-9F6C-48A5-9194-660CA6492563} = {15229511-9F6C-48A5-9194-660CA6492563} + {D6CEB142-BF8E-471C-AE16-4300F2D5DEDA} = {D6CEB142-BF8E-471C-AE16-4300F2D5DEDA} + {B76BCE8C-63CC-4A99-88B5-D621D563E699} = {B76BCE8C-63CC-4A99-88B5-D621D563E699} + {C5CEB09E-79AF-4091-ACCF-D28EC3D7D90F} = {C5CEB09E-79AF-4091-ACCF-D28EC3D7D90F} + {443EB1A7-C634-4292-9F2D-C752BB2BF40F} = {443EB1A7-C634-4292-9F2D-C752BB2BF40F} + {70575BD2-A532-41B8-9A17-B9876E992A84} = {70575BD2-A532-41B8-9A17-B9876E992A84} + {D6FD0DE5-5305-458E-8CA5-FCA4B8E05B04} = {D6FD0DE5-5305-458E-8CA5-FCA4B8E05B04} + EndProjectSection +EndProject Global GlobalSection(SolutionConfigurationPlatforms) = preSolution Debug|Win32 = Debug|Win32 @@ -508,6 +521,14 @@ Global {33E97266-F213-3639-BCAB-2F3F95E15B16}.Debug|Win32.Build.0 = Debug|Win32 {33E97266-F213-3639-BCAB-2F3F95E15B16}.Release|Win32.ActiveCfg = Release|Win32 {33E97266-F213-3639-BCAB-2F3F95E15B16}.Release|Win32.Build.0 = Release|Win32 + {70575BD2-A532-41B8-9A17-B9876E992A84}.Debug|Win32.ActiveCfg = Debug|Win32 + {70575BD2-A532-41B8-9A17-B9876E992A84}.Debug|Win32.Build.0 = Debug|Win32 + {70575BD2-A532-41B8-9A17-B9876E992A84}.Release|Win32.ActiveCfg = Release|Win32 + {70575BD2-A532-41B8-9A17-B9876E992A84}.Release|Win32.Build.0 = Release|Win32 + {A6DF0D74-E4D1-4EF4-A642-59B0CF3E74BB}.Debug|Win32.ActiveCfg = Debug|Win32 + {A6DF0D74-E4D1-4EF4-A642-59B0CF3E74BB}.Debug|Win32.Build.0 = Debug|Win32 + {A6DF0D74-E4D1-4EF4-A642-59B0CF3E74BB}.Release|Win32.ActiveCfg = Release|Win32 + {A6DF0D74-E4D1-4EF4-A642-59B0CF3E74BB}.Release|Win32.Build.0 = Release|Win32 EndGlobalSection GlobalSection(SolutionProperties) = preSolution HideSolutionNode = FALSE diff --git a/lib/Makefile.am b/lib/Makefile.am index 792a13424..6a44bf7cf 100644 --- a/lib/Makefile.am +++ b/lib/Makefile.am @@ -3,7 +3,7 @@ SUBDIRS = cdt graph cgraph pathplan sfio vmalloc ast \ vpsc rbtree ortho sparse patchwork expr common \ - pack xdot label gvc ingraphs topfish glcomp \ + pack xdot label gvc ingraphs topfish glcomp mingle \ circogen dotgen dotgen2 fdpgen neatogen twopigen sfdpgen osage gvpr EXTRA_DIST = Makefile.old gvc.vcproj gvc.def diff --git a/lib/mingle/Makefile.am b/lib/mingle/Makefile.am new file mode 100644 index 000000000..a92a4a57d --- /dev/null +++ b/lib/mingle/Makefile.am @@ -0,0 +1,21 @@ +# $Id$Revision: +## Process this file with automake to produce Makefile.in + +AM_CPPFLAGS = \ + -I$(top_srcdir) \ + -I$(top_srcdir)/lib/common \ + -I$(top_srcdir)/lib/sparse \ + -I$(top_srcdir)/lib/gvc \ + -I$(top_srcdir)/lib/pathplan \ + -I$(top_srcdir)/lib/sfdpgen \ + -I$(top_srcdir)/lib/cgraph \ + -I$(top_srcdir)/lib/cdt $(ANN_CFLAGS) + +noinst_HEADERS = edge_bundling.h ink.h agglomerative_bundling.h nearest_neighbor_graph.h nearest_neighbor_graph_ann.h + +noinst_LTLIBRARIES = libmingle_C.la + +libmingle_C_la_SOURCES = edge_bundling.c ink.c agglomerative_bundling.c nearest_neighbor_graph.c nearest_neighbor_graph_ann.cpp + +EXTRA_DIST = minglelib.vcproj + diff --git a/lib/mingle/agglomerative_bundling.c b/lib/mingle/agglomerative_bundling.c new file mode 100644 index 000000000..61de39f3a --- /dev/null +++ b/lib/mingle/agglomerative_bundling.c @@ -0,0 +1,773 @@ +/************************************************************************* + * Copyright (c) 2011 AT&T Intellectual Property + * All rights reserved. This program and the accompanying materials + * are made available under the terms of the Eclipse Public License v1.0 + * which accompanies this distribution, and is available at + * http://www.eclipse.org/legal/epl-v10.html + * + * Contributors: See CVS logs. Details at http://www.graphviz.org/ + *************************************************************************/ + +#include "types.h" +#include "globals.h" +#include "general.h" +#include "time.h" +#include "SparseMatrix.h" +#include "vector.h" +#include "edge_bundling.h" +#include "ink.h" +#include "agglomerative_bundling.h" +#include "nearest_neighbor_graph.h" + +#if OPENGL +#include "gl.h" +extern pedge *edges_global; +extern int nedges_global; +#endif + +enum {DEBUG=0}; + +static Agglomerative_Ink_Bundling Agglomerative_Ink_Bundling_init(SparseMatrix A, pedge *edges, int level){ + Agglomerative_Ink_Bundling grid; + int n = A->n, i; + + assert(SparseMatrix_is_symmetric(A, TRUE)); + + if (!A) return NULL; + assert(A->m == n); + grid = MALLOC(sizeof(struct Agglomerative_Ink_Bundling_struct)); + grid->level = level; + grid->n = n; + grid->A = A; + grid->P = NULL; + grid->R0 = NULL; + grid->R = NULL; + grid->next = NULL; + grid->prev = NULL; + grid->inks = MALLOC(sizeof(real)*(A->m)); + grid->edges = edges; + grid->delete_top_level_A = 0; + grid->total_ink = -1; + if (level == 0){ + real total_ink = 0; + for (i = 0; i < n; i++) { + (grid->inks)[i] = ink1(edges[i]); + total_ink += (grid->inks)[i]; + } + grid->total_ink = total_ink; + } + return grid; +} + +static void Agglomerative_Ink_Bundling_delete(Agglomerative_Ink_Bundling 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); + /* on level 0, R0 = NULL, on level 1, R0 = R */ + if (grid->level > 1) SparseMatrix_delete(grid->R0); + SparseMatrix_delete(grid->R); + FREE(grid->inks); + + Agglomerative_Ink_Bundling_delete(grid->next); + FREE(grid); +} + +static Agglomerative_Ink_Bundling Agglomerative_Ink_Bundling_establish(Agglomerative_Ink_Bundling grid, int *pick, real angle_param, real angle){ + /* pick is a work array of dimension n, with n the total number of original edges */ + int *matching; + SparseMatrix A = grid->A; + int n = grid->n, level = grid->level, nc = 0; + int *ia = A->ia, *ja = A->ja; + // real *a; + int i, j, k, jj, jc, jmax, ni, nj, npicks; + int *mask; + pedge *edges = grid->edges; + real *inks = grid->inks, *cinks, inki, inkj; + real gain, maxgain, minink, total_gain = 0; + int *ip = NULL, *jp = NULL, ie; + Vector *cedges;/* a table listing the content of bundled edges in the coarsen grid. + cedges[i] contain the list of origonal edges that make up the bundle i in the next level */ + real ink0, ink1, grand_total_ink = 0, grand_total_gain = 0; + point_t meet1, meet2; + + if (Verbose) fprintf(stderr,"level ===================== %d, n = %d\n",grid->level, n); + cedges = MALLOC(sizeof(Vector)*n); + cinks = MALLOC(sizeof(real)*n); + for (i = 0; i < n; i++) cedges[i] = Vector_new(1, sizeof(int), NULL); + + if (grid->level > 0){ + ip = grid->R0->ia; + jp = grid->R0->ja; + } + + matching = MALLOC(sizeof(int)*n); + mask = MALLOC(sizeof(real)*n); + for (i = 0; i < n; i++) mask[i] = -1; + + assert(n == A->n); + for (i = 0; i < n; i++) matching[i] = UNMATCHED; + + // a = (real*) A->a; + for (i = 0; i < n; i++){ + if (matching[i] != UNMATCHED) continue; + + /* find the best matching in ink saving */ + maxgain = 0; + minink = -1; + jmax = -1; + for (j = ia[i]; j < ia[i+1]; j++){ + jj = ja[j]; + if (jj == i) continue; + + /* ink saving of merging i and j */ + if ((jc=matching[jj]) == UNMATCHED){ + /* neither i nor jj are matched */ + inki = inks[i]; inkj = inks[jj]; + if (ip && jp){/* not the first level */ + ni = (ip[i+1] - ip[i]);/* number of edges represented by i */ + nj = (ip[jj+1] - ip[jj]);/* number of edges represented by jj */ + MEMCPY(pick, &(jp[ip[i]]), sizeof(int)*ni); + MEMCPY(pick+ni, &(jp[ip[jj]]), sizeof(int)*nj); + } else {/* first level */ + pick[0] = i; pick[1] = jj; + ni = nj = 1; + } + if (Verbose && DEBUG) fprintf(stderr, "ink(%d)=%f, ink(%d)=%f", i, inki, jj, inkj); + } else { + /* j is already matched. Its content is on cedges[jc] */ + inki = inks[i]; inkj = cinks[jc]; + if (Verbose && DEBUG) fprintf(stderr, "ink(%d)=%f, ink(%d->%d)=%f", i, inki, jj, jc, inkj); + if (ip) { + ni = (ip[i+1] - ip[i]);/* number of edges represented by i */ + MEMCPY(pick, &(jp[ip[i]]), sizeof(int)*ni); + } else { + ni = 1; pick[0] = i; + } + nj = Vector_get_length(cedges[jc]); + npicks = ni; + for (k = 0; k < nj; k++) { + pick[npicks++] = *((int*) Vector_get(cedges[jc], k)); + } + } + + npicks = ni + nj; + ink1 = ink(edges, npicks, pick, &ink0, &meet1, &meet2, angle_param, angle); + if (Verbose && DEBUG) { + fprintf(stderr,", if merging {"); + for (k = 0; k < npicks; k++) fprintf(stderr,"%d,", pick[k]); + fprintf(stderr,"}, "); + fprintf(stderr, " ink0=%f, ink1=%f", inki+inkj, ink1); + } + + gain = inki + inkj - ink1; + if (Verbose && DEBUG) fprintf(stderr, " gain=%f", gain); + if (gain > maxgain){ + maxgain = gain; + minink = ink1; + jmax = jj; + if (Verbose && DEBUG) fprintf(stderr, "maxgain=%f", maxgain); + } + if (Verbose && DEBUG) fprintf(stderr, "\n"); + + + + } + + + /* now merge i and jmax */ + if (maxgain > 0){ + /* a good bundling of i and another edge jmax is found */ + total_gain += maxgain; + jc = matching[jmax]; + if (jc == UNMATCHED){/* i and j both unmatched. Add j in the table first */ + if (Verbose && DEBUG) printf("maxgain=%f, merge %d with best edge: %d to form coarsen edge %d. Ink=%f\n",maxgain, i, jmax, nc, minink); + matching[i] = matching[jmax] = nc; + if (ip){ + for (k = ip[jmax]; k < ip[jmax+1]; k++) { + ie = jp[k]; + Vector_add(cedges[nc], (void*) (&ie)); + } + } else { + Vector_add(cedges[nc], (void*) (&jmax)); + } + jc = nc; + nc++; + } else {/*j is already matched */ + if (Verbose && DEBUG) printf("maxgain=%f, merge %d with existing cluster %d\n",maxgain, i, jc); + matching[i] = jc; + grand_total_ink -= cinks[jc];/* ink of cluster jc was already added, and will be added again as part of a larger cluster with i, so dicount */ + } + } else {/*i can not match/bundle successfully */ + if (Verbose && DEBUG) fprintf(stderr, "no gain in bundling node %d\n",i); + assert(maxgain <= 0); + matching[i] = nc; + jc = nc; + minink = inks[i]; + nc++; + } + + + /* add i to the appropriate table */ + if (ip){ + for (k = ip[i]; k < ip[i+1]; k++) { + ie = jp[k]; + Vector_add(cedges[jc], (void*) (&ie)); + } + } else { + Vector_add(cedges[jc], (void*) (&i)); + } + cinks[jc] = minink; + grand_total_ink += minink; + grand_total_gain += maxgain; + + if (Verbose && DEBUG){ + fprintf(stderr," coarse edge[%d]={",jc); + for (k = 0; k < Vector_get_length(cedges[jc]); k++) { + fprintf(stderr,"%d,", *((int*) Vector_get(cedges[jc], k))); + } + fprintf(stderr,"}, grand_total_gain=%f\n",grand_total_gain); + } + + } + + if (nc >= 1 && total_gain > 0){ + /* now set up restriction and prolongation operator */ + SparseMatrix P, R, R1, R0, B, cA; + real one = 1.; + Agglomerative_Ink_Bundling cgrid; + + R1 = SparseMatrix_new(nc, n, 1, MATRIX_TYPE_REAL, FORMAT_COORD); + for (i = 0; i < n; i++){ + jj = matching[i]; + SparseMatrix_coordinate_form_add_entries(R1, 1, &jj, &i, &one); + } + R = SparseMatrix_from_coordinate_format(R1); + SparseMatrix_delete(R1); + P = SparseMatrix_transpose(R); + B = SparseMatrix_multiply(R, A); + if (!B) goto RETURN; + cA = SparseMatrix_multiply(B, P); + if (!cA) goto RETURN; + SparseMatrix_delete(B); + grid->P = P; + grid->R = R; + + level++; + cgrid = Agglomerative_Ink_Bundling_init(cA, edges, level); + + + /* set up R0!!! */ + if (grid->R0){ + R0 = SparseMatrix_multiply(R, grid->R0); + } else { + assert(grid->level == 0); + R0 = R; + } + cgrid->R0 = R0; + cgrid->inks = cinks; + cgrid->total_ink = grand_total_ink; + + if (Verbose) fprintf(stderr,"level %d->%d, edges %d -> %d, ink %f->%f , gain = %f, or %f\n", grid->level, cgrid->level, grid->n, + cgrid->n, grid->total_ink, grand_total_ink, grid->total_ink - grand_total_ink, grand_total_gain); + assert(ABS(grid->total_ink - cgrid->total_ink - grand_total_gain) <= 0.0001*grid->total_ink); + + cgrid = Agglomerative_Ink_Bundling_establish(cgrid, pick, angle_param, angle); + grid->next = cgrid; + cgrid->prev = grid; + + } else { + if (Verbose) fprintf(stderr,"no more improvement, orig ink = %f, gain = %f, stop and final bundling found\n", grand_total_ink, grand_total_gain); + /* no more improvement, stop and final bundling found */ + for (i = 0; i < n; i++) matching[i] = i; + } + + RETURN: + FREE(matching); + for (i = 0; i < n; i++) Vector_delete(cedges[i]); + FREE(cedges); + FREE(mask); + return grid; +} + + + +static Agglomerative_Ink_Bundling Agglomerative_Ink_Bundling_aggressive_establish(Agglomerative_Ink_Bundling grid, int *pick, real angle_param, real angle){ + /* this does a true single-linkage clustering: find the edge that gives the best saving, merge, find again. + As oppose to: find the edge of node i that gives the best ink saving, merge, then do the same for node i+1, etc etc. + Right now it is implemented as a quick hack to check whether it is worth while: it saves about 3% extra ink on airlines: from + 59% to 62% + */ + /* pick is a work array of dimension n, with n the total number of original edges */ + int *matching; + SparseMatrix A = grid->A; + int n = grid->n, level = grid->level, nc = 0; + int *ia = A->ia, *ja = A->ja; + // real *a; + int i, j, k, jj, jc = -1, jmax, imax, ni, nj, npicks; + int *mask; + pedge *edges = grid->edges; + real *inks = grid->inks, *cinks, inki, inkj; + real gain, maxgain, minink, total_gain = 0; + int *ip = NULL, *jp = NULL, ie; + Vector *cedges;/* a table listing the content of bundled edges in the coarsen grid. + cedges[i] contain the list of origonal edges that make up the bundle i in the next level */ + real ink0, ink1, grand_total_ink = 0, grand_total_gain = 0; + point_t meet1, meet2; + + if (Verbose) fprintf(stderr,"level ===================== %d, n = %d\n",grid->level, n); + cedges = MALLOC(sizeof(Vector)*n); + cinks = MALLOC(sizeof(real)*n); + for (i = 0; i < n; i++) cedges[i] = Vector_new(1, sizeof(int), NULL); + + if (grid->level > 0){ + ip = grid->R0->ia; + jp = grid->R0->ja; + } + + matching = MALLOC(sizeof(int)*n); + mask = MALLOC(sizeof(real)*n); + for (i = 0; i < n; i++) mask[i] = -1; + + assert(n == A->n); + for (i = 0; i < n; i++) matching[i] = UNMATCHED; + + //a = (real*) A->a; + + do { + maxgain = 0; + imax = -1; + jmax = -1; + minink = -1; + for (i = 0; i < n; i++){ + if (matching[i] != UNMATCHED) continue; + + /* find the best matching in ink saving */ + for (j = ia[i]; j < ia[i+1]; j++){ + jj = ja[j]; + if (jj == i) continue; + + /* ink saving of merging i and j */ + if ((jc=matching[jj]) == UNMATCHED){ + /* neither i nor jj are matched */ + inki = inks[i]; inkj = inks[jj]; + if (ip && jp){/* not the first level */ + ni = (ip[i+1] - ip[i]);/* number of edges represented by i */ + nj = (ip[jj+1] - ip[jj]);/* number of edges represented by jj */ + MEMCPY(pick, &(jp[ip[i]]), sizeof(int)*ni); + MEMCPY(pick+ni, &(jp[ip[jj]]), sizeof(int)*nj); + } else {/* first level */ + pick[0] = i; pick[1] = jj; + ni = nj = 1; + } + if (Verbose && DEBUG) fprintf(stderr, "ink(%d)=%f, ink(%d)=%f", i, inki, jj, inkj); + } else { + /* j is already matched. Its content is on cedges[jc] */ + inki = inks[i]; inkj = cinks[jc]; + if (Verbose && DEBUG) fprintf(stderr, "ink(%d)=%f, ink(%d->%d)=%f", i, inki, jj, jc, inkj); + if (ip) { + ni = (ip[i+1] - ip[i]);/* number of edges represented by i */ + MEMCPY(pick, &(jp[ip[i]]), sizeof(int)*ni); + } else { + ni = 1; pick[0] = i; + } + nj = Vector_get_length(cedges[jc]); + npicks = ni; + for (k = 0; k < nj; k++) { + pick[npicks++] = *((int*) Vector_get(cedges[jc], k)); + } + } + + npicks = ni + nj; + ink1 = ink(edges, npicks, pick, &ink0, &meet1, &meet2, angle_param, angle); + if (Verbose && DEBUG) { + fprintf(stderr,", if merging {"); + for (k = 0; k < npicks; k++) fprintf(stderr,"%d,", pick[k]); + fprintf(stderr,"}, "); + fprintf(stderr, " ink0=%f, ink1=%f", inki+inkj, ink1); + } + + gain = inki + inkj - ink1; + if (Verbose && DEBUG) fprintf(stderr, " gain=%f", gain); + if (gain > maxgain){ + maxgain = gain; + minink = ink1; + jmax = jj; + imax = i; + if (Verbose && DEBUG) fprintf(stderr, "maxgain=%f", maxgain); + } + if (Verbose && DEBUG) fprintf(stderr, "\n"); + + + + } + } + + /* now merge i and jmax */ + if (maxgain > 0){ + /* a good bundling of i and another edge jmax is found */ + total_gain += maxgain; + jc = matching[jmax]; + if (jc == UNMATCHED){/* i and j both unmatched. Add j in the table first */ + if (Verbose && DEBUG) printf("maxgain=%f, merge %d with best edge: %d to form coarsen edge %d. Ink=%f\n",maxgain, imax, jmax, nc, minink); + matching[imax] = matching[jmax] = nc; + if (ip){ + for (k = ip[jmax]; k < ip[jmax+1]; k++) { + ie = jp[k]; + Vector_add(cedges[nc], (void*) (&ie)); + } + } else { + Vector_add(cedges[nc], (void*) (&jmax)); + } + jc = nc; + nc++; + } else {/*j is already matched */ + if (Verbose && DEBUG) printf("maxgain=%f, merge %d with existing cluster %d\n",maxgain, i, jc); + matching[imax] = jc; + grand_total_ink -= cinks[jc];/* ink of cluster jc was already added, and will be added again as part of a larger cluster with i, so dicount */ + } + /* add i to the appropriate table */ + if (ip){ + for (k = ip[imax]; k < ip[imax+1]; k++) { + ie = jp[k]; + Vector_add(cedges[jc], (void*) (&ie)); + } + } else { + Vector_add(cedges[jc], (void*) (&imax)); + } + cinks[jc] = minink; + grand_total_ink += minink; + grand_total_gain += maxgain; + } else {/*i can not match/bundle successfully */ + if (Verbose && DEBUG) fprintf(stderr, "no gain in bundling node %d\n",i); + for (i = 0; i < n; i++){ + assert(maxgain <= 0); + if (matching[i] == UNMATCHED){ + imax = i; + matching[imax] = nc; + jc = nc; + minink = inks[imax]; + nc++; + + if (ip){ + for (k = ip[imax]; k < ip[imax+1]; k++) { + ie = jp[k]; + Vector_add(cedges[jc], (void*) (&ie)); + } + } else { + Vector_add(cedges[jc], (void*) (&imax)); + } + cinks[jc] = minink; + grand_total_ink += minink; + grand_total_gain += maxgain; + + + + + } + } + } + + } while (maxgain > 0); + + + if (Verbose && DEBUG){ + fprintf(stderr," coarse edge[%d]={",jc); + for (k = 0; k < Vector_get_length(cedges[jc]); k++) { + fprintf(stderr,"%d,", *((int*) Vector_get(cedges[jc], k))); + } + fprintf(stderr,"}\n"); + } + + if (nc >= 1 && total_gain > 0){ + /* now set up restriction and prolongation operator */ + SparseMatrix P, R, R1, R0, B, cA; + real one = 1.; + Agglomerative_Ink_Bundling cgrid; + + R1 = SparseMatrix_new(nc, n, 1, MATRIX_TYPE_REAL, FORMAT_COORD); + for (i = 0; i < n; i++){ + jj = matching[i]; + SparseMatrix_coordinate_form_add_entries(R1, 1, &jj, &i, &one); + } + R = SparseMatrix_from_coordinate_format(R1); + SparseMatrix_delete(R1); + P = SparseMatrix_transpose(R); + B = SparseMatrix_multiply(R, A); + if (!B) goto RETURN; + cA = SparseMatrix_multiply(B, P); + if (!cA) goto RETURN; + SparseMatrix_delete(B); + grid->P = P; + grid->R = R; + + level++; + cgrid = Agglomerative_Ink_Bundling_init(cA, edges, level); + + + /* set up R0!!! */ + if (grid->R0){ + R0 = SparseMatrix_multiply(R, grid->R0); + } else { + assert(grid->level == 0); + R0 = R; + } + cgrid->R0 = R0; + cgrid->inks = cinks; + cgrid->total_ink = grand_total_ink; + + if (Verbose) fprintf(stderr,"level %d->%d, edges %d -> %d, ink %f->%f , gain = %f, or %f\n", grid->level, cgrid->level, grid->n, + cgrid->n, grid->total_ink, grand_total_ink, grid->total_ink - grand_total_ink, grand_total_gain); + assert(ABS(grid->total_ink - cgrid->total_ink - grand_total_gain) <= 0.0001*grid->total_ink); + + cgrid = Agglomerative_Ink_Bundling_aggressive_establish(cgrid, pick, angle_param, angle); + grid->next = cgrid; + cgrid->prev = grid; + + } else { + if (Verbose) fprintf(stderr,"no more improvement, orig ink = %f, gain = %f, stop and final bundling found\n", grand_total_ink, grand_total_gain); + /* no more improvement, stop and final bundling found */ + for (i = 0; i < n; i++) matching[i] = i; + } + + RETURN: + FREE(matching); + for (i = 0; i < n; i++) Vector_delete(cedges[i]); + FREE(mask); + return grid; +} + +static Agglomerative_Ink_Bundling Agglomerative_Ink_Bundling_new(SparseMatrix A0, pedge *edges, real angle_param, real angle){ + /* give a link of edges and their nearest neighbor graph, return a multilevel of edge bundling based on ink saving */ + Agglomerative_Ink_Bundling grid; + int *pick; + SparseMatrix A = A0; + + if (!SparseMatrix_is_symmetric(A, FALSE) || A->type != MATRIX_TYPE_REAL){ + A = SparseMatrix_get_real_adjacency_matrix_symmetrized(A); + } + grid = Agglomerative_Ink_Bundling_init(A, edges, 0); + + pick = MALLOC(sizeof(int)*A0->m); + + //grid = Agglomerative_Ink_Bundling_aggressive_establish(grid, pick, angle_param, angle); + grid = Agglomerative_Ink_Bundling_establish(grid, pick, angle_param, angle); + FREE(pick); + + if (A != A0) grid->delete_top_level_A = TRUE;/* be sure to clean up later */ + + return grid; +} + +static pedge* agglomerative_ink_bundling_internal(int dim, SparseMatrix A, pedge* edges, int nneighbors, int *recurse_level, int MAX_RECURSE_LEVEL, real angle_param, real angle, int open_gl, real *current_ink, real *ink00, int *flag){ + + int i, j, jj, k; + int *ia, *ja; + int *pick; + Agglomerative_Ink_Bundling grid, cgrid; + SparseMatrix R; + real ink0, ink1; + point_t meet1, meet2; + pedge e; + real TOL = 0.0001, wgt_all; + clock_t start; + + (*recurse_level)++; + if (Verbose) fprintf(stderr, "agglomerative_ink_bundling_internal, recurse level ------- %d\n",*recurse_level); + + assert(A->m == A->n); + + *flag = 0; + + start = clock(); + grid = Agglomerative_Ink_Bundling_new(A, edges, angle_param, angle); + if (Verbose) + fprintf(stderr, "CPU for agglomerative bundling %f\n", ((real) (clock() - start))/CLOCKS_PER_SEC); + ink0 = grid->total_ink; + + /* find coarsest */ + cgrid = grid; + while (cgrid->next){ + cgrid = cgrid->next; + } + ink1 = cgrid->total_ink; + + if (*current_ink < 0){ + *current_ink = *ink00 = ink0; + if (Verbose) + fprintf(stderr,"initial total ink = %f\n",*current_ink); + } + if (ink1 < ink0){ + *current_ink -= ink0 - ink1; + } + + if (Verbose) fprintf(stderr,"ink: %f->%f, edges: %d->%d, current ink = %f, percentage gain over original = %f\n", ink0, ink1, grid->n, cgrid->n, *current_ink, (ink0-ink1)/(*ink00)); + + /* if no meaningful improvement (0.0001%), out, else rebundle the middle section */ + if ((ink0-ink1)/(*ink00) < 0.000001 || *recurse_level > MAX_RECURSE_LEVEL) { + /* project bundles up */ + R = cgrid->R0; + if (R){ + ia = R->ia; + ja = R->ja; + for (i = 0; i < R->m; i++){ + pick = &(ja[ia[i]]); + + if (Verbose && DEBUG) fprintf(stderr,"calling ink2...\n"); + ink1 = ink(edges, ia[i+1]-ia[i], pick, &ink0, &meet1, &meet2, angle_param, angle); + if (Verbose && DEBUG) fprintf(stderr,"finish calling ink2...\n"); + assert(ABS(ink1 - cgrid->inks[i])<=MAX(TOL, TOL*ink1) && ink1 - ink0 <= TOL); + wgt_all = 0.; + if (ia[i+1]-ia[i] > 1){ + for (j = ia[i]; j < ia[i+1]; j++){ + /* make this edge 4 points, insert two meeting points at 1 and 2, make 3 the last point */ + jj = ja[j]; + edges[jj] = pedge_double(edges[jj]);/* has to call pedge_double twice: from 2 points to 3 points to 5 points. The last point not used, may be + improved later */ + e = edges[jj] = pedge_double(edges[jj]); + + e->wgts = REALLOC(e->wgts, sizeof(real)*4); + e->x[1*dim] = meet1.x; + e->x[1*dim+1] = meet1.y; + e->x[2*dim] = meet2.x; + e->x[2*dim+1] = meet2.y; + e->x[3*dim] = e->x[4*dim]; + e->x[3*dim+1] = e->x[4*dim+1]; + e->npoints = 4; + for (k = 0; k < 3; k++) e->wgts[k] = e->wgt; + wgt_all += e->wgt; + + } + for (j = ia[i]; j < ia[i+1]; j++){ + e = edges[ja[j]]; + e->wgts[1] = wgt_all; + } + } + + } + } + } else { + pedge *mid_edges, midedge;/* middle section of edges that will be bundled again */ + real *xx; + int ne, npp, l; + SparseMatrix A_mid; + real eps = 0., wgt, total_wgt = 0; + + /* make new edges using meet1 and meet2. + + call Agglomerative_Ink_Bundling_new + + inherit new edges to old edges + */ + + R = cgrid->R0; + assert(R && cgrid != grid);/* if ink improved, we should have gone at leat 1 level down! */ + ia = R->ia; + ja = R->ja; + ne = R->m; + mid_edges = MALLOC(sizeof(pedge)*ne); + xx = MALLOC(sizeof(real)*4*ne); + for (i = 0; i < R->m; i++){ + pick = &(ja[ia[i]]); + wgt = 0.; + for (j = ia[i]; j < ia[i+1]; j++) wgt += edges[j]->wgt; + total_wgt += wgt; + if (Verbose && DEBUG) fprintf(stderr,"calling ink3...\n"); + ink1 = ink(edges, ia[i+1]-ia[i], pick, &ink0, &meet1, &meet2, angle_param, angle); + if (Verbose && DEBUG) fprintf(stderr,"done calling ink3...\n"); + assert(ABS(ink1 - cgrid->inks[i])<=MAX(TOL, TOL*ink1) && ink1 - ink0 <= TOL); + xx[i*4 + 0] = meet1.x; + xx[i*4 + 1] = meet1.y; + xx[i*4 + 2] = meet2.x; + xx[i*4 + 3] = meet2.y; + mid_edges[i] = pedge_wgt_new(2, dim, &(xx[i*4]), wgt); + //mid_edges[i] = pedge_wgt_new(2, dim, &(xx[i*4]), 1.); + + } + + A_mid = nearest_neighbor_graph(ne, MIN(nneighbors, ne), 4, xx, eps); + + agglomerative_ink_bundling_internal(dim, A_mid, mid_edges, nneighbors, recurse_level, MAX_RECURSE_LEVEL, angle_param, angle, open_gl, current_ink, ink00, flag); + SparseMatrix_delete(A_mid); + FREE(xx); + + /* patching edges with the new mid-section */ + for (i = 0; i < R->m; i++){ + pick = &(ja[ia[i]]); + midedge = mid_edges[i]; + npp = midedge->npoints + 2; + for (j = ia[i]; j < ia[i+1]; j++){ + jj = ja[j]; + e = edges[jj] = pedge_wgts_realloc(edges[jj], npp); + + assert(e->npoints = 2); + for (l = 0; l < dim; l++){/* move the second point to the last */ + e->x[(npp - 1)*dim+l] = e->x[1*dim+l]; + } + + for (k = 0; k < midedge->npoints; k++){ + for (l = 0; l < dim; l++){ + e->x[(k+1)*dim+l] = midedge->x[k*dim+l]; + } + if (k < midedge->npoints - 1){ + if (midedge->wgts){ + e->wgts[(k+1)] = midedge->wgts[k]; + } else { + e->wgts[(k+1)] = midedge->wgt; + } + } + } + e->wgts[npp - 2] = e->wgts[0];/* the last interval take from the 1st interval */ + + + e->npoints = npp; + } +#ifdef OPENGL + if (open_gl & FALSE){ + nedges_global = grid->n; + edges_global = edges; + drawScene(); + // waitie(5./R->m); + } +#endif + } + + for (i = 0; i < ne; i++) pedge_delete(mid_edges[i]); + + } + + +#ifdef OPENGL + if (open_gl){ + nedges_global = grid->n; + edges_global = edges; + drawScene(); + // waitie(5./R->m); + } +#endif + + Agglomerative_Ink_Bundling_delete(grid); + + return edges; +} + + +pedge* agglomerative_ink_bundling(int dim, SparseMatrix A, pedge* edges, int nneighbor, int MAX_RECURSE_LEVEL, real angle_param, real angle, int open_gl, int *flag){ + int recurse_level = 0; + real current_ink = -1, ink0; + pedge *edges2; + + ink_count = 0; + edges2 = agglomerative_ink_bundling_internal(dim, A, edges, nneighbor, &recurse_level, MAX_RECURSE_LEVEL, angle_param, angle, open_gl, ¤t_ink, &ink0, flag); + + + if (Verbose) + fprintf(stderr,"initial total ink = %f, final total ink = %f, inksaving = %f percent, total ink_calc = %f, avg ink_calc per edge = %f\n", ink0, current_ink, (ink0-current_ink)/ink0, ink_count, ink_count/(real) A->m); + return edges2; +} + diff --git a/lib/mingle/agglomerative_bundling.h b/lib/mingle/agglomerative_bundling.h new file mode 100644 index 000000000..0cbfd8e33 --- /dev/null +++ b/lib/mingle/agglomerative_bundling.h @@ -0,0 +1,34 @@ +/************************************************************************* + * Copyright (c) 2011 AT&T Intellectual Property + * All rights reserved. This program and the accompanying materials + * are made available under the terms of the Eclipse Public License v1.0 + * which accompanies this distribution, and is available at + * http://www.eclipse.org/legal/epl-v10.html + * + * Contributors: See CVS logs. Details at http://www.graphviz.org/ + *************************************************************************/ + +#ifndef AGGLOMERATIVE_BUNDLING_H +#define AGGLOMERATIVE_BUNDLING_H + +typedef struct Agglomerative_Ink_Bundling_struct *Agglomerative_Ink_Bundling; + +struct Agglomerative_Ink_Bundling_struct { + int level;/* 0, 1, ... */ + int n; + SparseMatrix A; /* n x n matrix, where n is the number of edges/bundles in this level */ + SparseMatrix P; /* prolongation matrix from level + 1 to level */ + SparseMatrix R0;/* this is basically R[level - 1].R[level - 2]...R[0], which gives the map of bundling i to the original edges: first row of R0 gives + the nodes on the finest grid corresponding to the coarsest node 1, etc */ + SparseMatrix R;/* striction mtrix from level to level + 1*/ + Agglomerative_Ink_Bundling next; + Agglomerative_Ink_Bundling prev; + real *inks; /* amount of ink needed to draw this edge/bundle. Dimension n. */ + real total_ink; /* amount of ink needed to draw this edge/bundle. Dimension n. */ + pedge* edges; /* the original edge info. This does not vary level to level and is of dimenion n0, where n0 is the number of original edges */ + int delete_top_level_A;/*whether the top level matrix should be deleted on garbage collecting the grid */ +}; + +pedge* agglomerative_ink_bundling(int dim, SparseMatrix A, pedge* edges, int nneighbor, int max_recursion, real angle_param, real angle, int open_gl, int *flag); + +#endif /* AGGLOMERATIVE_BUNDLING_H */ diff --git a/lib/mingle/edge_bundling.c b/lib/mingle/edge_bundling.c new file mode 100644 index 000000000..fa114b504 --- /dev/null +++ b/lib/mingle/edge_bundling.c @@ -0,0 +1,831 @@ +/************************************************************************* + * Copyright (c) 2011 AT&T Intellectual Property + * All rights reserved. This program and the accompanying materials + * are made available under the terms of the Eclipse Public License v1.0 + * which accompanies this distribution, and is available at + * http://www.eclipse.org/legal/epl-v10.html + * + * Contributors: See CVS logs. Details at http://www.graphviz.org/ + *************************************************************************/ + +#ifdef HAVE_CONFIG_H +#include "config.h" +#endif + +#include "types.h" +#include "globals.h" +#include "general.h" +#include "SparseMatrix.h" +#include "edge_bundling.h" +#include +#include "clustering.h" +#include "ink.h" +#include "agglomerative_bundling.h" + +#define SMALL 1.e-10 + +#ifdef OPENGL +#include "gl.h" +extern pedge *edges_global; +extern int *clusters_global; +#endif + +static real norm(int n, real *x){ + real res = 0; + int i; + + for (i = 0; i < n; i++) res += x[i]*x[i]; + return sqrt(res); + +} +static real sqr_dist(int dim, real *x, real *y){ + int i; + real res = 0; + for (i = 0; i < dim; i++) res += (x[i] - y[i])*(x[i] - y[i]); + return res; +} +static real dist(int dim, real *x, real *y){ + return sqrt(sqr_dist(dim,x,y)); +} + +pedge pedge_new(int np, int dim, real *x){ + pedge e; + + e = MALLOC(sizeof(struct pedge_struct)); + e->npoints = np; + e->dim = dim; + e->len = np; + e->x = MALLOC(dim*(e->len)*sizeof(real)); + MEMCPY(e->x, x, dim*(e->len)*sizeof(real)); + e->edge_length = dist(dim, &(x[0*dim]), &(x[(np-1)*dim])); + e->wgt = 1.; + e->wgts = NULL; + return e; + +} +pedge pedge_wgt_new(int np, int dim, real *x, real wgt){ + pedge e; + int i; + + e = MALLOC(sizeof(struct pedge_struct)); + e->npoints = np; + e->dim = dim; + e->len = np; + e->x = MALLOC(dim*(e->len)*sizeof(real)); + MEMCPY(e->x, x, dim*(e->len)*sizeof(real)); + e->edge_length = dist(dim, &(x[0*dim]), &(x[(np-1)*dim])); + e->wgt = wgt; + e->wgts = MALLOC(sizeof(real)*(np - 1)); + for (i = 0; i < np - 1; i++) e->wgts[i] = wgt; + return e; + +} +void pedge_delete(pedge e){ + FREE(e->x); + FREE(e); +} + +pedge pedge_flip(pedge e){ + /* flip the polyline so that last point becomes the first, second last the second, etc*/ + real *y; + real *x = e->x; + int i, dim = e->dim; + int n = e->npoints; + + y = MALLOC(sizeof(real)*e->dim); + for (i = 0; i < (e->npoints)/2; i++){ + MEMCPY(y, &(x[i*dim]), sizeof(real)*dim); + MEMCPY(&(x[(n-1-i)*dim]), &(x[i*dim]), sizeof(real)*dim); + MEMCPY(&(x[i*dim]), y, sizeof(real)*dim); + } + FREE(y); + return e; +} + +static real edge_compatibility(pedge e1, pedge e2){ + /* two edges are u1->v1, u2->v2. + return 1 if two edges are exactly the same, 0 if they are very different. + */ + real *u1, *v1, *u2, *v2, *u, dist1, dist2, len1, len2; + int dim = e1->dim, flipped = FALSE; + + u1 = e1->x; + v1 = (e1->x)+(e1->npoints)*dim-dim; + u2 = e2->x; + v2 = (e2->x)+(e2->npoints)*dim-dim; + dist1 = sqr_dist(dim, u1, u2) + sqr_dist(dim, v1, v2); + dist2 = sqr_dist(dim, u1, v2) + sqr_dist(dim, v1, u2); + if (dist1 > dist2){ + u = u2; + u2 = v2; + v2 = u; + dist1 = dist2; + flipped = TRUE; + } + len1 = dist(dim, u1, v1); + len2 = dist(dim, u2, v2); + //dist1 = MAX(0.1, dist1/(len1+len2+dist1)); + dist1 = MAX(0.1, dist1/(len1+len2+0.0001*dist1)); + if (flipped){ + return -1/dist1; + } else { + return 1/dist1; + } +} + +static real edge_compatibility_full(pedge e1, pedge e2){ + /* two edges are u1->v1, u2->v2. + return 1 if two edges are exactly the same, 0 if they are very different. + This is based on Holten and van Wijk's paper + */ + real *u1, *v1, *u2, *v2, *u, dist1, dist2, len1, len2, len; + real tmp, ca, cp, cs; + int dim = e1->dim, flipped = FALSE, i; + + u1 = e1->x; + v1 = (e1->x)+(e1->npoints)*dim-dim; + u2 = e2->x; + v2 = (e2->x)+(e2->npoints)*dim-dim; + dist1 = sqr_dist(dim, u1, u2) + sqr_dist(dim, v1, v2); + dist2 = sqr_dist(dim, u1, v2) + sqr_dist(dim, v1, u2); + if (dist1 > dist2){ + u = u2; + u2 = v2; + v2 = u; + dist1 = dist2; + flipped = TRUE; + } + len1 = MAX(dist(dim, u1, v1), SMALL); + len2 = MAX(dist(dim, u2, v2), SMALL); + len = 0.5*(len1+len2); + + /* angle compatability */ + ca = 0; + for (i = 0; i < dim; i++) + ca += (v1[i]-u1[i])*(v2[i]-u2[i]); + ca = ABS(ca/(len1*len2)); + assert(ca > -0.001); + + /* scale compatability */ + //cs = 2/(len1/len2+len2/len1); + cs = 2/(MAX(len1,len2)/len + len/MIN(len1, len2)); + assert(cs > -0.001 && cs < 1.001); + + /* position compatability */ + cp = 0; + for (i = 0; i < dim; i++) { + tmp = .5*(v1[i]+u1[i])-.5*(v2[i]+u2[i]); + cp += tmp*tmp; + } + cp = sqrt(cp); + cp = len/(len + cp); + assert(cp > -0.001 && cp < 1.001); + + /* visibility compatability */ + + //dist1 = MAX(0.1, dist1/(len1+len2+dist1)); + dist1 = cp*ca*cs; + if (flipped){ + return -dist1; + } else { + return dist1; + } +} + +static void fprint_rgb(FILE* fp, int r, int g, int b, int alpha){ + fprintf(fp,"#"); + if (r >= 16) { + fprintf(fp,"%2x",r); + } else { + fprintf(fp,"0%1x",r); + } + if (g >= 16) { + fprintf(fp,"%2x",g); + } else { + fprintf(fp,"0%1x",g); + } + if (b >= 16) { + fprintf(fp,"%2x",b); + } else { + fprintf(fp,"0%1x",b); + } + if (alpha >= 16) { + fprintf(fp,"%2x",alpha); + } else { + fprintf(fp,"0%1x",alpha); + } + +} + +void pedge_export_gv(FILE *fp, int ne, pedge *edges){ + pedge edge; + real *x, t; + int i, j, k, kk, dim, sta, sto; + real maxwgt = 0, len, len_total, len_total0; + int r, g, b; + // real tt1[3]={0.25,0.5,0.75}; + // real tt2[4]={0.2,0.4,0.6,0.8}; + real tt1[3]={0.15,0.5,0.85}; + real tt2[4]={0.15,0.4,0.6,0.85}; + real *tt; + + fprintf(fp,"strict graph{\n"); + /* points */ + for (i = 0; i < ne; i++){ + edge = edges[i]; + x = edge->x; + dim = edge->dim; + sta = 0; sto = edge->npoints - 1; + + fprintf(fp, "%d [pos=\"", i); + for (k = 0; k < dim; k++) { + if (k != 0) fprintf(fp, ","); + fprintf(fp, "%f", x[sta*dim+k]); + } + fprintf(fp, "\"];\n"); + + fprintf(fp, "%d [pos=\"", i + ne); + for (k = 0; k < dim; k++) { + if (k != 0) fprintf(fp, ","); + fprintf(fp, "%f", x[sto*dim+k]); + } + fprintf(fp, "\"];\n"); + + } + + /* figure out max number of bundled origional edges in a pedge */ + for (i = 0; i < ne; i++){ + edge = edges[i]; + if (edge->wgts){ + for (j = 0; j < edge->npoints - 1; j++){ + maxwgt = MAX(maxwgt, edge->wgts[j]); + } + } + } + + /* spline and colors */ + for (i = 0; i < ne; i++){ + fprintf(fp,"%d -- %d [pos=\"", i, i + ne); + edge = edges[i]; + x = edge->x; + dim = edge->dim; + /* splines */ + for (j = 0; j < edge->npoints; j++){ + if (j != 0) { + int mm = 3; + fprintf(fp," "); + /* there are ninterval+1 points, add 3*ninterval+2 points, get rid of internal ninternal-1 points, + make into 3*ninterval+4 points so that gviz spline rendering can work */ + if (j == 1 || j == edge->npoints - 1) { + /* every interval gets 3 points inmserted except the first and last one */ + tt = tt2; + mm = 4; + } else { + tt = tt1; + } + for (kk = 1; kk <= mm; kk++){ + //t = kk/(real) (mm+1); + t = tt[kk-1]; + for (k = 0; k < dim; k++) { + if (k != 0) fprintf(fp,","); + fprintf(fp, "%f", (x[(j-1)*dim+k]*(1-t)+x[j*dim+k]*(t))); + } + fprintf(fp," "); + } + } + if (j == 0 || j == edge->npoints - 1){ + for (k = 0; k < dim; k++) { + if (k != 0) fprintf(fp,","); + fprintf(fp, "%f", x[j*dim+k]); + } + } + } + /* colors based on how much bundling */ + if (edge->wgts){ + fprintf(fp, "\", wgts=\""); + for (j = 0; j < edge->npoints - 1; j++){ + if (j != 0) fprintf(fp,","); + fprintf(fp, "%f", edge->wgts[j]); + } + + len_total = 0; + len_total0 = 0; + fprintf(fp, "\", color=\""); + for (j = 0; j < edge->npoints - 1; j++){ + len = 0; + for (k = 0; k < dim; k++){ + len += (edge->x[dim*j+k] - edge->x[dim*(j+1)+k])*(edge->x[dim*j+k] - edge->x[dim*(j+1)+k]); + } + len = sqrt(len/k); + len_total0 += len; + } + for (j = 0; j < edge->npoints - 1; j++){ + len = 0; + for (k = 0; k < dim; k++){ + len += (edge->x[dim*j+k] - edge->x[dim*(j+1)+k])*(edge->x[dim*j+k] - edge->x[dim*(j+1)+k]); + } + len = sqrt(len/k); + len_total += len; + t = edge->wgts[j]/maxwgt; + /* interpotate between red (t = 1) to blue (t = 0) */ + r = 255*t; g = 0; b = 255*(1-t); b = 255*(1-t); + if (j != 0) fprintf(fp,":"); + fprint_rgb(fp, r, g, b, 85); + fprintf(fp,",%f",len_total/len_total0); + } + + } + fprintf(fp, "\"];\n"); + + } + fprintf(fp,"}\n"); +} + +void pedge_export_mma(FILE *fp, int ne, pedge *edges){ + pedge edge; + real *x; + int i, j, k, dim; + + fprintf(fp,"Graphics[{"); + /* points */ + fprintf(fp,"{Red, "); + for (i = 0; i < ne; i++){ + if (i != 0) fprintf(fp,","); + fprintf(fp,"Point["); + edge = edges[i]; + x = edge->x; + dim = edge->dim; + fprintf(fp, "{"); + for (j = 0; j < edge->npoints; j+= edge->npoints - 1){ + if (j != 0) fprintf(fp,","); + fprintf(fp, "{"); + for (k = 0; k < dim; k++) { + if (k != 0) fprintf(fp,","); + fprintf(fp, "%f", x[j*dim+k]); + } + fprintf(fp, "}"); + } + fprintf(fp, "}"); + fprintf(fp, "]"); + } + fprintf(fp,"},\n{GrayLevel[0.5,0.2], "); + + /* spline */ + for (i = 0; i < ne; i++){ + if (i != 0) fprintf(fp,","); + fprintf(fp,"Spline["); + edge = edges[i]; + x = edge->x; + dim = edge->dim; + fprintf(fp, "{"); + for (j = 0; j < edge->npoints; j++){ + if (j != 0) fprintf(fp,","); + fprintf(fp, "{"); + for (k = 0; k < dim; k++) { + if (k != 0) fprintf(fp,","); + fprintf(fp, "%f", x[j*dim+k]); + } + fprintf(fp, "}"); + } + fprintf(fp, "}"); + fprintf(fp, ", Bezier]"); + } + fprintf(fp,"}"); + + fprintf(fp,"}]\n"); +} + +void pedge_print(char *comments, pedge e){ + int i, j, dim; + dim = e->dim; + fprintf(stderr,"%s", comments); + for (i = 0; i < e->npoints; i++){ + if (i > 0) fprintf(stderr,","); + fprintf(stderr,"{"); + for (j = 0; j < dim; j++){ + if (j > 0) fprintf(stderr,","); + fprintf(stderr,"%f",e->x[dim*i+j]); + } + fprintf(stderr,"}"); + } + fprintf(stderr,"\n"); +} + +pedge pedge_realloc(pedge e, int n){ + if (n <= e->npoints) return e; + e->x = REALLOC(e->x, e->dim*n*sizeof(real)); + if (e->wgts) e->wgts = REALLOC(e->wgts, (n-1)*sizeof(real)); + e->len = n; + return e; +} +pedge pedge_wgts_realloc(pedge e, int n){ + /* diff from pedge_alloc: allocate wgts if do not exist and initialize to wgt */ + int i; + if (n <= e->npoints) return e; + e->x = REALLOC(e->x, e->dim*n*sizeof(real)); + if (!(e->wgts)){ + e->wgts = REALLOC(e->wgts, (n-1)*sizeof(real)); + for (i = 0; i < e->npoints; i++) e->wgts[i] = e->wgt; + } else { + e->wgts = REALLOC(e->wgts, (n-1)*sizeof(real)); + } + e->len = n; + return e; +} + + +pedge pedge_double(pedge e){ + /* double the number of points (more precisely, add a point between two points in the polyline */ + int npoints = e->npoints, len = e->len, i, dim = e->dim; + real *x; + int j, ii, ii2, np; + + assert(npoints >= 2); + // pedge_print("before doubling edge = ", e); + if (npoints*2-1 > len){ + len = 3*npoints; + e->x = REALLOC(e->x, dim*len*sizeof(real)); + } + x = e->x; + + for (i = npoints - 1; i >= 0; i--){ + ii = 2*i; + for (j = 0; j < dim; j++){ + x[dim*ii + j] = x[dim*i + j]; + } + } + + for (i = 0; i < npoints - 1; i++){ + ii = 2*i;/* left and right interpolant of a new point */ + ii2 = 2*(i+1); + for (j = 0; j < dim; j++){ + x[dim*(2*i + 1) + j] = 0.5*(x[dim*ii + j] + x[dim*ii2 + j]); + } + } + e->len = len; + np = e->npoints = 2*(e->npoints) - 1; + e->edge_length = dist(dim, &(x[0*dim]), &(x[(np-1)*dim])); + + //pedge_print("after doubling edge = ", e); + + return e; +} + +static void edge_tension_force(real *force, pedge e){ + real *x = e->x; + int dim = e->dim; + int np = e->npoints; + int i, left, right, j; + //real dist_left, dist_right; + real s; + + + /* tention force = ((np-1)*||2x-xleft-xright||)/||e||, so the force is norminal and unitless + */ + //s = (np-1)*(np-1)/MAX(SMALL, e->edge_length); + s = (np-1)/MAX(SMALL, e->edge_length); + for (i = 1; i <= np - 2; i++){ + left = i - 1; + right = i + 1; + // dist_left = dist(dim, &(x[i*dim]), &(x[left*dim])); + // dist_right = dist(dim, &(x[i*dim]), &(x[right*dim])); + for (j = 0; j < dim; j++) force[i*dim + j] += s*(x[left*dim + j] - x[i*dim + j]); + for (j = 0; j < dim; j++) force[i*dim + j] += s*(x[right*dim + j] - x[i*dim + j]); + } +} + + +#if 0 +static void edge_tension_force2(real *force, pedge e){ + /* in this version each node is pulled towards its original position on the line */ + real *x = e->x; + int dim = e->dim; + int np = e->npoints; + int i, j; + //int left, right; + //real dist_left, dist_right; + real s; + + + /* tention force = ((np-1)*||2x-xleft-xright||)/||e||, so the force is norminal and unitless + */ + s = .1/MAX(SMALL, e->edge_length); + for (i = 1; i <= np - 2; i++){ + //left = i - 1; + // right = i + 1; + // dist_left = dist(dim, &(x[i*dim]), &(x[left*dim])); + // dist_right = dist(dim, &(x[i*dim]), &(x[right*dim])); + for (j = 0; j < dim; j++) force[i*dim + j] += s*((i*x[0*dim + j]+(np-1-i)*x[(np-1)*dim + j])/(np-1) - x[i*dim + j]); + } +} +#endif + +static void edge_attraction_force(real similarity, pedge e1, pedge e2, real *force){ + /* attrractive force from x2 applied to x1 */ + real *x1 = e1->x, *x2 = e2->x; + int dim = e1->dim; + int np = e1->npoints; + int i, j; + real dist, s, ss; + real edge_length = e1->edge_length; + + + assert(e1->npoints == e2->npoints); + + /* attractive force = 1/d where d = D/||e1|| is the relative distance, D is the distance between e1 and e2. + so the force is norminal and unitless + */ + if (similarity > 0){ + s = edge_length; + s = similarity*edge_length; + for (i = 1; i <= np - 2; i++){ + dist = sqr_dist(dim, &(x1[i*dim]), &(x2[i*dim])); + if (dist < SMALL) dist = SMALL; + ss = s/(dist+0.1*edge_length*sqrt(dist)); + for (j = 0; j < dim; j++) force[i*dim + j] += ss*(x2[i*dim + j] - x1[i*dim + j]); + } + } else {/* clip e2 */ + s = -edge_length; + s = -similarity*edge_length; + for (i = 1; i <= np - 2; i++){ + dist = sqr_dist(dim, &(x1[i*dim]), &(x2[(np - 1 - i)*dim])); + if (dist < SMALL) dist = SMALL; + ss = s/(dist+0.1*edge_length*sqrt(dist)); + for (j = 0; j < dim; j++) force[i*dim + j] += ss*(x2[(np - 1 - i)*dim + j] - x1[i*dim + j]); + } + } + +} + +static pedge* force_directed_edge_bundling(SparseMatrix A, pedge* edges, int maxit, real step0, real K, int open_gl){ + int i, j, ne = A->n, k; + int *ia = A->ia, *ja = A->ja, iter = 0; + real *a = (real*) A->a; + pedge e1, e2; + real *force_t, *force_a; + int np = edges[0]->npoints, dim = edges[0]->dim; + real *x; + real step = step0; + real fnorm_a, fnorm_t, edge_length, start; + + if (Verbose) + fprintf(stderr, "total interaction pairs = %d out of %d, avg neighbors per edge = %f\n",A->nz, A->m*A->m, A->nz/(real) A->m); + + force_t = MALLOC(sizeof(real)*dim*(np)); + force_a = MALLOC(sizeof(real)*dim*(np)); + while (step > 0.001 && iter < maxit){ + start = clock(); + iter++; + for (i = 0; i < ne; i++){ + for (j = 0; j < dim*np; j++) { + force_t[j] = 0.; + force_a[j] = 0.; + } + e1 = edges[i]; + x = e1->x; + //edge_tension_force2(force_t, e1); + edge_tension_force(force_t, e1); + for (j = ia[i]; j < ia[i+1]; j++){ + e2 = edges[ja[j]]; + edge_attraction_force(a[j], e1, e2, force_a); + } + fnorm_t = MAX(SMALL, norm(dim*(np - 2), &(force_t[1*dim]))); + fnorm_a = MAX(SMALL, norm(dim*(np - 2), &(force_a[1*dim]))); + edge_length = e1->edge_length; + + // fprintf(stderr,"tension force norm = %f att force norm = %f step = %f edge_length = %f\n",fnorm_t, fnorm_a, step, edge_length); + for (j = 1; j <= np - 2; j++){ + for (k = 0; k < dim; k++) { + x[j*dim + k] += step*edge_length*(force_t[j*dim+k] + K*force_a[j*dim+k])/(sqrt(fnorm_t*fnorm_t + K*K*fnorm_a*fnorm_a)); + } + } + + /* + fprintf(stderr,"edge %d",i); + for (j = 0; j < np; j++){ + if (j != 0) fprintf(stderr,","); + fprintf(stderr,"{"); + for (k = 0; k < dim; k++) { + if (k != 0) fprintf(stderr,","); + fprintf(stderr,"%f", x[j*dim+k]); + } + fprintf(stderr,"}"); + } + fprintf(stderr,"\n"); + */ + + } + step = step*0.9; + fprintf(stderr, "iter ==== %d cpu = %f npoints = %d\n",iter, ((real) (clock() - start))/CLOCKS_PER_SEC, np - 2); + +#ifdef OPENGL + if (open_gl){ + edges_global = edges; + drawScene(); + } +#endif + + } + + FREE(force_t); + FREE(force_a); + return edges; +} + +static real absfun(real x){ + return ABS(x); +} + + + + +pedge* modularity_ink_bundling(int dim, int ne, SparseMatrix B, pedge* edges, real angle_param, real angle){ + int *assignment = NULL, flag, nclusters; + real modularity; + int *clusterp, *clusters; + SparseMatrix D, C; + point_t meet1, meet2; + real ink0, ink1; + pedge e; + int i, j, jj; + int use_value_for_clustering = TRUE; + + //int use_value_for_clustering = FALSE; + + SparseMatrix BB; + + /* B may contain negative entries */ + BB = SparseMatrix_copy(B); + BB = SparseMatrix_apply_fun(BB, absfun); + modularity_clustering(BB, TRUE, 0, use_value_for_clustering, &nclusters, &assignment, &modularity, &flag); + SparseMatrix_delete(BB); + +#ifdef OPENGL + clusters_global = assignment; +#endif + + assert(!flag); + if (Verbose) fprintf(stderr, "there are %d clusters, modularity = %f\n",nclusters, modularity); + + C = SparseMatrix_new(1, 1, 1, MATRIX_TYPE_PATTERN, FORMAT_COORD); + + for (i = 0; i < ne; i++){ + jj = assignment[i]; + SparseMatrix_coordinate_form_add_entries(C, 1, &jj, &i, NULL); + } + + D = SparseMatrix_from_coordinate_format(C); + SparseMatrix_delete(C); + clusterp = D->ia; + clusters = D->ja; + for (i = 0; i < nclusters; i++){ + ink1 = ink(edges, clusterp[i+1] - clusterp[i], &(clusters[clusterp[i]]), &ink0, &meet1, &meet2, angle_param, angle); + fprintf(stderr,"nedges = %d ink0 = %f, ink1 = %f\n",clusterp[i+1] - clusterp[i], ink0, ink1); + if (ink1 < ink0){ + for (j = clusterp[i]; j < clusterp[i+1]; j++){ + /* make this edge 5 points, insert two meeting points at 1 and 2, make 3 the last point */ + edges[clusters[j]] = pedge_double(edges[clusters[j]]); + e = edges[clusters[j]] = pedge_double(edges[clusters[j]]); + e->x[1*dim] = meet1.x; + e->x[1*dim+1] = meet1.y; + e->x[2*dim] = meet2.x; + e->x[2*dim+1] = meet2.y; + e->x[3*dim] = e->x[4*dim]; + e->x[3*dim+1] = e->x[4*dim+1]; + e->npoints = 4; + } +#ifdef OPENGL + edges_global = edges; + drawScene(); +#endif + } + } + SparseMatrix_delete(D); + return edges; +} + +static SparseMatrix check_compatibility(SparseMatrix A, int ne, pedge *edges, int compatibility_method, real tol){ + /* go through the links and make sure edges are compatable */ + SparseMatrix B, C; + int *ia, *ja, i, j, jj; + real start; + real dist; + + B = SparseMatrix_new(1, 1, 1, MATRIX_TYPE_REAL, FORMAT_COORD); + ia = A->ia; ja = A->ja; + //SparseMatrix_print("A",A); + start = clock(); + for (i = 0; i < ne; i++){ + for (j = ia[i]; j < ia[i+1]; j++){ + jj = ja[j]; + if (i == jj) continue; + if (compatibility_method == COMPATIBILITY_DIST){ + dist = edge_compatibility_full(edges[i], edges[jj]); + } else if (compatibility_method == COMPATIBILITY_FULL){ + dist = edge_compatibility(edges[i], edges[jj]); + } + /* + fprintf(stderr,"edge %d = {",i); + pedge_print("", edges[i]); + fprintf(stderr,"edge %d = {",jj); + pedge_print("", edges[jj]); + fprintf(stderr,"compatibility=%f\n",dist); + */ + + if (ABS(dist) > tol){ + B = SparseMatrix_coordinate_form_add_entries(B, 1, &i, &jj, &dist); + B = SparseMatrix_coordinate_form_add_entries(B, 1, &jj, &i, &dist); + } + } + } + C = SparseMatrix_from_coordinate_format(B); + //SparseMatrix_print("C",C); + SparseMatrix_delete(B); + B = C; + fprintf(stderr, "edge compatibilitu time = %f\n",((real) (clock() - start))/CLOCKS_PER_SEC); + return B; +} + +pedge* edge_bundling(SparseMatrix A0, int dim, real *x, int maxit_outer, real K, int method, int nneighbor, int compatibility_method, + int max_recursion, real angle_param, real angle, int open_gl){ + /* bundle edges. + A: edge graph + x: edge i is at {p,q}, + . where p = x[2*dim*i : 2*dim*i+dim-1] + . and q = x[2*dim*i+dim : 2*dim*i+2*dim-1] + maxit_outer: max outer iteration for force directed bundling. Every outer iter subdivide each edge segment into 2. + K: norminal edge length in force directed bundling + method: which method to use. + nneighbor: number of neighbors to be used in forming nearest neighbor graph. Used only in agglomerative method + compatibility_method: which method to use to calculate compatibility. Used only in force directed. + max_recursion: used only in agglomerative method. Specify how many level of recursion to do to bundle bundled edges again + open_gl: whether to plot in X. + + */ + int ne = A0->m; + pedge *edges; + SparseMatrix A = A0, B = NULL; + int i; + real tol = 0.001; + int k; + real step0 = 0.1, start; + int maxit = 10; + int flag; + + assert(A->n == ne); + edges = MALLOC(sizeof(pedge)*ne); + + for (i = 0; i < ne; i++){ + edges[i] = pedge_new(2, dim, &(x[dim*2*i])); + } + + A = SparseMatrix_symmetrize(A0, TRUE); + + + + if (Verbose) start = clock(); +#ifndef HAVE_ANN + method = METHOD_INK; +#endif + if (method == METHOD_INK){ + + /* go through the links and make sure edges are compatable */ + B = check_compatibility(A, ne, edges, compatibility_method, tol); + + edges = modularity_ink_bundling(dim, ne, B, edges, angle_param, angle); + + } else if (method == METHOD_INK_AGGLOMERATE){ +#ifdef HAVE_ANN + /* plan: merge a node with its neighbors if doing so improve. Form coarsening graph, repeat until no more ink saving */ + edges = agglomerative_ink_bundling(dim, A, edges, nneighbor, max_recursion, angle_param, angle, open_gl, &flag); + assert(!flag); +#else + agerr (AGERR, "Graphviz built without approximate nearest neighbor library ANN; agglomerative inking not available\n"); + edges = edges; +#endif + } else if (method == METHOD_FD){/* FD method */ + + /* go through the links and make sure edges are compatable */ + B = check_compatibility(A, ne, edges, compatibility_method, tol); + + + for (k = 0; k < maxit_outer; k++){ + for (i = 0; i < ne; i++){ + edges[i] = pedge_double(edges[i]); + } + step0 = step0/2; + edges = force_directed_edge_bundling(B, edges, maxit, step0, K, open_gl); + } + + } else if (method == METHOD_NONE){ + edges = edges; + } else { + assert(0); + } + if (Verbose) + fprintf(stderr, "total edge bundling cpu = %f\n",((real) (clock() - start))/CLOCKS_PER_SEC); + if (B != A) SparseMatrix_delete(B); + if (A != A0) SparseMatrix_delete(A); + + return edges; +} diff --git a/lib/mingle/edge_bundling.h b/lib/mingle/edge_bundling.h new file mode 100644 index 000000000..7d6018d91 --- /dev/null +++ b/lib/mingle/edge_bundling.h @@ -0,0 +1,45 @@ +/************************************************************************* + * Copyright (c) 2011 AT&T Intellectual Property + * All rights reserved. This program and the accompanying materials + * are made available under the terms of the Eclipse Public License v1.0 + * which accompanies this distribution, and is available at + * http://www.eclipse.org/legal/epl-v10.html + * + * Contributors: See CVS logs. Details at http://www.graphviz.org/ + *************************************************************************/ +#ifndef EDGE_BUNDLING_H +#define EDGE_BUNDLING_H + +#include + +struct pedge_struct { + real wgt; /* weight, telling how many original edges this edge represent. If this edge consists of multiple sections of different weights then this is a lower bound. This only applied for agglomerative bundling */ + int npoints;/* number of poly points */ + int len;/* length of arra x. len >= npoints */ + int dim;/* dim >= 2. Point i is stored from x[i*dim]*/ + real edge_length; + real *x;/* coordinates of the npoints poly points. Dimension dim*npoints */ + real *wgts;/* number of original edges each section represnet. Dimension npoint - 1. This only applied for agglomerative bundling Null for other methods */ +}; + +typedef struct pedge_struct* pedge; + +pedge* edge_bundling(SparseMatrix A, int dim, real *x, int maxit_outer, real K, int method, int nneighbor, int compatibility_method, int max_recursion, real angle_param, real angle, int open_gl); +void pedge_delete(pedge e); +pedge pedge_realloc(pedge e, int np); +pedge pedge_wgts_realloc(pedge e, int n); +void pedge_export_mma(FILE *fp, int ne, pedge *edges); +void pedge_export_gv(FILE *fp, int ne, pedge *edges); +enum {METHOD_NONE = -1, METHOD_FD, METHOD_INK_AGGLOMERATE, METHOD_INK}; +enum {COMPATIBILITY_DIST = 0, COMPATIBILITY_FULL}; +pedge pedge_new(int np, int dim, real *x); +pedge pedge_wgt_new(int np, int dim, real *x, real wgt); +pedge pedge_double(pedge e); + +/* flip the polyline so that last point becomes the first, second last the second, etc*/ +pedge pedge_flip(pedge e); + + +#endif /* EDGE_BUNDLING_H */ + + diff --git a/lib/mingle/ink.c b/lib/mingle/ink.c new file mode 100644 index 000000000..9bee033f6 --- /dev/null +++ b/lib/mingle/ink.c @@ -0,0 +1,357 @@ +/************************************************************************* + * Copyright (c) 2011 AT&T Intellectual Property + * All rights reserved. This program and the accompanying materials + * are made available under the terms of the Eclipse Public License v1.0 + * which accompanies this distribution, and is available at + * http://www.eclipse.org/legal/epl-v10.html + * + * Contributors: See CVS logs. Details at http://www.graphviz.org/ + *************************************************************************/ + +#include +#include +#include "types.h" +#include "globals.h" +#include "general.h" +#include "SparseMatrix.h" +#include "edge_bundling.h" +#include "ink.h" + +double ink_count; + +point_t addPoint (point_t a, point_t b) +{ + a.x += b.x; + a.y += b.y; + return a; +} + +point_t subPoint (point_t a, point_t b) +{ + a.x -= b.x; + a.y -= b.y; + return a; +} + +point_t scalePoint (point_t a, double d) +{ + a.x *= d; + a.y *= d; + return a; +} + +double dotPoint(point_t a, point_t b){ + return a.x*b.x + a.y*b.y; +} + + +point_t Origin; + +/* sumLengths: + */ +double sumLengths_avoid_bad_angle(point_t* points, int npoints, point_t end, point_t meeting, real angle_param) +{ + /* avoid sharp turns, we want cos_theta to be as close to -1 as possible */ + int i; + double len0, len, sum = 0; + double diff_x, diff_y, diff_x0, diff_y0; + double cos_theta, cos_max = -10; + + diff_x0 = end.x-meeting.x; + diff_y0 = end.y-meeting.y; + len0 = sum = sqrt(diff_x0*diff_x0+diff_y0*diff_y0); + + // distance form each of 'points' till 'meeting' + for (i=0; iprec && dotPoint(diff, diff) > 1.e-20); + + meeting = scalePoint(addPoint(first,fourth),0.5); + *meet = meeting; + + return sumLengths(points, npoints, end, meeting); + +} + +static double project_to_line(point_t pt, point_t left, point_t right, real angle){ + /* pt + ^ ^ + . \ \ + . \ \ + d . a\ \ + . \ \ + . \ \ + . c \ alpha \ b + .<------left:0 ----------------------------> right:1. Find the projection of pt on the left--right line. If the turning angle is small, + | | + |<-------f--------- + we should get a negative number. Let a := left->pt, b := left->right, then we are calculating: + c = |a| cos(a,b)/|b| b + d = a - c + f = -ctan(alpha)*|d|/|b| b + and the project of alpha degree on the left->right line is + c-f = |a| cos(a,b)/|b| b - -ctan(alpha)*|d|/|b| b + = (|a| a.b/(|a||b|) + ctan(alpha)|a-c|)/|b| b + = (a.b/|b| + ctan(alpha)|a-c|)/|b| b + the dimentionless projection is: + a.b/|b|^2 + ctan(alpha)|a-c|/|b| + = a.b/|b|^2 + ctan(alpha)|d|/|b| + */ + + + point_t b, a; + real bnorm, dnorm; + real alpha, ccord; + + if (angle <=0 || angle >= M_PI) return 2;/* return outside of the interval which should be handled as a sign of infeasible turning angle */ + alpha = angle; + + assert(alpha > 0 && alpha < M_PI); + b = subPoint(right, left); + a = subPoint(pt, left); + bnorm = MAX(1.e-10, dotPoint(b, b)); + ccord = dotPoint(b, a)/bnorm; + dnorm = dotPoint(a,a)/bnorm - ccord*ccord; + if (alpha == M_PI/2){ + return ccord; + } else { + // assert(dnorm >= MIN(-1.e-5, -1.e-5*bnorm)); + return ccord + sqrt(MAX(0, dnorm))/tan(alpha); + } +} + + + + + + + + + +/* ink: + * Compute minimal ink used the input edges are bundled. + * Assumes tails all occur on one side and heads on the other. + */ +double ink(pedge* edges, int numEdges, int *pick, double *ink0, point_t *meet1, point_t *meet2, real angle_param, real angle) +{ + int i; + point_t begin, end, mid, diff; + pedge e; + real *x; + point_t* sources = N_NEW(numEdges, point_t); + point_t* targets = N_NEW(numEdges, point_t); + double inkUsed; + //double eps = 1.0e-2; + double eps = 1.0e-2; + double cend = 0, cbegin = 0; + double wgt = 0; + + // fprintf(stderr,"in ink code ========\n"); + ink_count += numEdges; + + *ink0 = 0; + + /* canonicalize so that edges 1,2,3 and 3,2,1 gives the same optimal ink */ + if (pick) vector_sort_int(numEdges, pick, TRUE); + + begin = end = Origin; + for (i = 0; i < numEdges; i++) { + if (pick) { + e = edges[pick[i]]; + } else { + e = edges[i]; + } + x = e->x; + sources[i].x = x[0]; + sources[i].y = x[1]; + targets[i].x = x[e->dim*e->npoints - e->dim]; + targets[i].y = x[e->dim*e->npoints - e->dim + 1]; + (*ink0) += sqrt((sources[i].x - targets[i].x)*(sources[i].x - targets[i].x) + (sources[i].y - targets[i].y)*(sources[i].y - targets[i].y)); + begin = addPoint (begin, scalePoint(sources[i], e->wgt)); + end = addPoint (end, scalePoint(targets[i], e->wgt)); + wgt += e->wgt; + //fprintf(stderr,"source={%f,%f}, target = {%f,%f}\n",sources[i].x, sources[i].y, + //targets[i].x, targets[i].y); + } + + begin = scalePoint (begin, 1.0/wgt); + end = scalePoint (end, 1.0/wgt); + + + if (numEdges == 1){ + *meet1 = begin; + *meet2 = end; + //fprintf(stderr,"ink used = %f\n",*ink0); + free (sources); + free (targets); + return *ink0; + } + + /* shift the begin and end point to avoid sharp turns */ + for (i = 0; i < numEdges; i++) { + if (pick) { + e = edges[pick[i]]; + } else { + e = edges[i]; + } + x = e->x; + sources[i].x = x[0]; + sources[i].y = x[1]; + targets[i].x = x[e->dim*e->npoints - e->dim]; + targets[i].y = x[e->dim*e->npoints - e->dim + 1]; + /* begin(1) ----------- mid(0) */ + if (i == 0){ + cbegin = project_to_line(sources[i], begin, end, angle); + cend = project_to_line(targets[i], end, begin, angle); + } else { + cbegin = MAX(cbegin, project_to_line(sources[i], begin, end, angle)); + cend = MAX(cend, project_to_line(targets[i], end, begin, angle)); + } + } + + if (angle > 0 && angle < M_PI){ + if (cbegin + cend > 1 || cbegin > 1 || cend > 1){ + /* no point can be found that satisfies the angular constraints, so we give up and set ink to a large value */ + if (Verbose && 0) fprintf(stderr,"no point satisfying any angle constraints can be found. cbeg=%f cend=%f\n",cbegin,cend); + inkUsed = 1000*(*ink0); + *meet1 = *meet2 = mid; + free (sources); + free (targets); + return inkUsed; + } + /* make sure the turning angle is no more than alpha degree */ + cbegin = MAX(0, cbegin);/* make sure the new adjusted point is with in [begin,end] internal */ + diff = subPoint(end, begin); + begin = addPoint(begin, scalePoint(diff, cbegin)); + + cend = MAX(0, cend);/* make sure the new adjusted point is with in [end,begin] internal */ + end = subPoint(end, scalePoint(diff, cend)); + } + mid = scalePoint (addPoint(begin,end),0.5); + + inkUsed = (bestInk (sources, numEdges, begin, mid, eps, meet1, angle_param) + + bestInk (targets, numEdges, end, mid, eps, meet2, angle_param)); + //fprintf(stderr,"beg={%f,%f}, meet1={%f,%f}, meet2={%f,%f}, mid={%f,%f}, end={%f,%f}\n",begin.x, begin.y, meet1->x, meet1->y, meet2->x, meet2->y, + //mid.x, mid.y, end.x, end.y); + + //fprintf(stderr,"ink used = %f\n",inkUsed); + // fprintf(stderr,"{cb,ce}={%f, %f} end={%f,%f}, meet={%f,%f}, mid={%f,%f}\n",cbegin, cend, end.x, end.y, meet2->x, meet2->y, mid.x, mid.y); + free (sources); + free (targets); + return inkUsed; +} + + +double ink1(pedge e){ + + + real *x, xx, yy; + + real ink0 = 0; + + x = e->x; + xx = x[0] - x[e->dim*e->npoints - e->dim]; + yy = x[1] - x[e->dim*e->npoints - e->dim + 1]; + ink0 += sqrt(xx*xx + yy*yy); + return ink0; +} diff --git a/lib/mingle/ink.h b/lib/mingle/ink.h new file mode 100644 index 000000000..09be9f52d --- /dev/null +++ b/lib/mingle/ink.h @@ -0,0 +1,37 @@ +/************************************************************************* + * Copyright (c) 2011 AT&T Intellectual Property + * All rights reserved. This program and the accompanying materials + * are made available under the terms of the Eclipse Public License v1.0 + * which accompanies this distribution, and is available at + * http://www.eclipse.org/legal/epl-v10.html + * + * Contributors: See CVS logs. Details at http://www.graphviz.org/ + *************************************************************************/ + +#ifndef INK_H +#define INK_H + +#include + +typedef struct { + double x, y; +} point_t; + +/* given a list of edges, find the best ink bundling by making them meet at 2 points + \ / + -meet1 ---------- meet2 - + / \ + edges: list of edges + numEdges: number of edges + pick: if not NULL, consider edges pick[0], pick[1], ..., pick[numedges-1], + . othetwise consider edges 0, 1, ..., numEdge-1 + ink0: ink needed if no bundling + meet1, meet2: meeting point + return: best ink needed if bundled. +*/ +double ink(pedge* edges, int numEdges, int *pick, double *ink0, point_t *meet1, point_t *meet2, real angle_param, real angle); +double ink1(pedge e); + +extern double ink_count; + +#endif /* INK_H */ diff --git a/lib/mingle/minglelib.vcproj b/lib/mingle/minglelib.vcproj new file mode 100644 index 000000000..3132ffe69 --- /dev/null +++ b/lib/mingle/minglelib.vcproj @@ -0,0 +1,213 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/lib/mingle/nearest_neighbor_graph.c b/lib/mingle/nearest_neighbor_graph.c new file mode 100644 index 000000000..02707c3d0 --- /dev/null +++ b/lib/mingle/nearest_neighbor_graph.c @@ -0,0 +1,56 @@ +/************************************************************************* + * Copyright (c) 2011 AT&T Intellectual Property + * All rights reserved. This program and the accompanying materials + * are made available under the terms of the Eclipse Public License v1.0 + * which accompanies this distribution, and is available at + * http://www.eclipse.org/legal/epl-v10.html + * + * Contributors: See CVS logs. Details at http://www.graphviz.org/ + *************************************************************************/ + +#ifdef HAVE_CONFIG_H +#include "config.h" +#endif + +#include "general.h" +#include "SparseMatrix.h" +#include "nearest_neighbor_graph_ann.h" +#include "nearest_neighbor_graph.h" + +SparseMatrix nearest_neighbor_graph(int nPts, int num_neigbors, int dim, double *x, double eps){ + /* Gives a nearest neighbor graph of a list of dim-dimendional points. The result is a sparse matrix + of nPts x nPts, with num_neigbors entries per row. + + nPts: number of points + num_neigbors: number of neighbors needed + dim: dimension + eps: error tolerance + x: nPts*dim vector. The i-th point is x[i*dim : i*dim + dim - 1] + + */ + int *irn = NULL, *jcn = NULL, nz; + real *val = NULL; + SparseMatrix A; + int k = num_neigbors; + + /* need to *2 as we do two sweeps of neighbors, so could have repeats */ + irn = MALLOC(sizeof(int)*nPts*k*2); + jcn = MALLOC(sizeof(int)*nPts*k*2); + val = MALLOC(sizeof(double)*nPts*k*2); + +#ifdef HAVE_ANN + nearest_neighbor_graph_ann(nPts, dim, num_neigbors, eps, x, &nz, &irn, &jcn, &val); + + A = SparseMatrix_from_coordinate_arrays(nz, nPts, nPts, irn, jcn, (void *) val, MATRIX_TYPE_REAL, sizeof(real)); +#else + A = NULL; +#endif + + FREE(irn); + FREE(jcn); + FREE(val); + + return A; + + +} diff --git a/lib/mingle/nearest_neighbor_graph.h b/lib/mingle/nearest_neighbor_graph.h new file mode 100644 index 000000000..3a042eb0d --- /dev/null +++ b/lib/mingle/nearest_neighbor_graph.h @@ -0,0 +1,16 @@ +/************************************************************************* + * Copyright (c) 2011 AT&T Intellectual Property + * All rights reserved. This program and the accompanying materials + * are made available under the terms of the Eclipse Public License v1.0 + * which accompanies this distribution, and is available at + * http://www.eclipse.org/legal/epl-v10.html + * + * Contributors: See CVS logs. Details at http://www.graphviz.org/ + *************************************************************************/ + +#ifndef NEAREST_NEIGHBOR_GRAPH_H +#define NEAREST_NEIGHBOR_GRAPH_H + +SparseMatrix nearest_neighbor_graph(int nPts, int num_neigbors, int dim, double *x, double eps); + +#endif /* NEAREST_NEIGHBOR_GRAPH_H */ diff --git a/lib/mingle/nearest_neighbor_graph_ann.cpp b/lib/mingle/nearest_neighbor_graph_ann.cpp new file mode 100644 index 000000000..02db85def --- /dev/null +++ b/lib/mingle/nearest_neighbor_graph_ann.cpp @@ -0,0 +1,182 @@ +/************************************************************************* + * Copyright (c) 2011 AT&T Intellectual Property + * All rights reserved. This program and the accompanying materials + * are made available under the terms of the Eclipse Public License v1.0 + * which accompanies this distribution, and is available at + * http://www.eclipse.org/legal/epl-v10.html + * + * Contributors: See CVS logs. Details at http://www.graphviz.org/ + *************************************************************************/ + +#ifdef HAVE_CONFIG_H +#include "config.h" +#endif + +#ifdef HAVE_ANN +//---------------------------------------------------------------------- +// File: get_nearest_neighb_graph.cpp +//---------------------------------------------------------------------- + +#include // C standard library +#include // C I/O (for sscanf) +#include // string manipulation +#include // file I/O +#include // ANN declarations + +using namespace std; // make std:: accessible +int dim = 4; // dimension + + +static void printPt(ostream &out, ANNpoint p) // print point +{ + out << "" << p[0]; + for (int i = 1; i < dim; i++) { + out << "," << p[i]; + } + out << ""; +} + +static void sortPtsX(int n, ANNpointArray pts){ + /* sort so that edges always go from left to right in x-doordinate */ + ANNpoint p; + ANNcoord x, y; + int i, j; + for (i = 0; i < n; i++){ + for (j = 0; j < dim; j++){ + p = pts[i]; + if (p[0] < p[2] || (p[0] == p[2] && p[1] < p[3])) continue; + x = p[0]; y = p[1]; + p[0] = p[2]; + p[1] = p[3]; + p[2] = x; + p[3] = y; + } + } +} + +static void sortPtsY(int n, ANNpointArray pts){ + /* sort so that edges always go from left to right in x-doordinate */ + ANNpoint p; + ANNcoord x, y; + int i, j; + for (i = 0; i < n; i++){ + for (j = 0; j < dim; j++){ + p = pts[i]; + if (p[1] < p[3] || (p[1] == p[3] && p[0] < p[2])) continue; + x = p[0]; y = p[1]; + p[0] = p[2]; + p[1] = p[3]; + p[2] = x; + p[3] = y; + } + } +} + + +extern "C" void nearest_neighbor_graph_ann(int nPts, int dim, int k, double eps, double *x, int *nz0, int **irn0, int **jcn0, double **val0); + +void nearest_neighbor_graph_ann(int nPts, int dim, int k, double eps, double *x, int *nz0, int **irn0, int **jcn0, double **val0){ + + /* Gives a nearest neighbor graph is a list of dim-dimendional points. The connectivity is in irn/jcn, and the distance in val. + + nPts: number of points + dim: dimension + k: number of neighbors needed + eps: error tolerance + x: nPts*dim vector. The i-th point is x[i*dim : i*dim + dim - 1] + nz: number of entries in the connectivity matrix irn/jcn/val + irn, jcn: the connectivity + val: the distance + + note that there could be repeates + */ + + ANNpointArray dataPts; // data points + ANNidxArray nnIdx; // near neighbor indices + ANNdistArray dists; // near neighbor distances + ANNkd_tree* kdTree; // search structure + + double *xx; + int *irn, *jcn; + double *val; + int nz; + + irn = *irn0; + jcn = *jcn0; + val = *val0; + + + dataPts = annAllocPts(nPts, dim); // allocate data points + nnIdx = new ANNidx[k]; // allocate near neigh indices + dists = new ANNdist[k]; // allocate near neighbor dists + + for (int i = 0; i < nPts; i++){ + xx = dataPts[i]; + for (int j = 0; j < dim; j++) xx[j] = x[i*dim + j]; + } + + //========= graph when sort based on x ======== + nz = 0; + sortPtsX(nPts, dataPts); + kdTree = new ANNkd_tree( // build search structure + dataPts, // the data points + nPts, // number of points + dim); // dimension of space + for (int ip = 0; ip < nPts; ip++){ + kdTree->annkSearch( // search + dataPts[ip], // query point + k, // number of near neighbors + nnIdx, // nearest neighbors (returned) + dists, // distance (returned) + eps); // error bound + + for (int i = 0; i < k; i++) { // print summary + if (nnIdx[i] == ip) continue; + val[nz] = dists[i]; + irn[nz] = ip; + jcn[nz++] = nnIdx[i]; + //cout << ip << "--" << nnIdx[i] << " [len = " << dists[i]<< ", weight = \"" << 1./(dists[i]) << "\", wt = \"" << 1./(dists[i]) << "\"]\n"; + //printPt(cout, dataPts[ip]); + //cout << "--"; + //printPt(cout, dataPts[nnIdx[i]]); + //cout << "\n"; + } + } + + + //========= graph when sort based on y ======== + sortPtsY(nPts, dataPts); + kdTree = new ANNkd_tree( // build search structure + dataPts, // the data points + nPts, // number of points + dim); // dimension of space + for (int ip = 0; ip < nPts; ip++){ + kdTree->annkSearch( // search + dataPts[ip], // query point + k, // number of near neighbors + nnIdx, // nearest neighbors (returned) + dists, // distance (returned) + eps); // error bound + + for (int i = 0; i < k; i++) { // print summary + if (nnIdx[i] == ip) continue; + val[nz] = dists[i]; + irn[nz] = ip; + jcn[nz++] = nnIdx[i]; + // cout << ip << "--" << nnIdx[i] << " [len = " << dists[i]<< ", weight = \"" << 1./(dists[i]) << "\", wt = \"" << 1./(dists[i]) << "\"]\n"; + } + } + + delete [] nnIdx; // clean things up + delete [] dists; + delete kdTree; + + *nz0 = nz; + + annClose(); // done with ANN + + + +} + +#endif /* HAVE_ANN */ diff --git a/lib/mingle/nearest_neighbor_graph_ann.h b/lib/mingle/nearest_neighbor_graph_ann.h new file mode 100644 index 000000000..f30ee5bc4 --- /dev/null +++ b/lib/mingle/nearest_neighbor_graph_ann.h @@ -0,0 +1,16 @@ +/************************************************************************* + * Copyright (c) 2011 AT&T Intellectual Property + * All rights reserved. This program and the accompanying materials + * are made available under the terms of the Eclipse Public License v1.0 + * which accompanies this distribution, and is available at + * http://www.eclipse.org/legal/epl-v10.html + * + * Contributors: See CVS logs. Details at http://www.graphviz.org/ + *************************************************************************/ + +#ifndef NEAREST_NEIGHBOR_GRAPH_ANN_H +#define NEAREST_NEIGHBOR_GRAPH_ANN_H + +void nearest_neighbor_graph_ann(int nPts, int dim, int k, double eps, double *x, int *nz0, int **irn0, int **jcn0, double **val0); + +#endif /* NEAREST_NEIGHBOR_GRAPH_ANN_H */ -- 2.40.0