]> granicus.if.org Git - graphviz/commitdiff
Add mingle program and library to graphviz
authorEmden R. Gansner <erg@research.att.com>
Wed, 14 Aug 2013 21:07:33 +0000 (17:07 -0400)
committerEmden R. Gansner <erg@research.att.com>
Wed, 14 Aug 2013 21:07:33 +0000 (17:07 -0400)
19 files changed:
cmd/Makefile.am
cmd/mingle/Makefile.am [new file with mode: 0644]
cmd/mingle/mingle.vcproj [new file with mode: 0644]
cmd/mingle/minglemain.c [new file with mode: 0644]
configure.ac
graphviz.sln
lib/Makefile.am
lib/mingle/Makefile.am [new file with mode: 0644]
lib/mingle/agglomerative_bundling.c [new file with mode: 0644]
lib/mingle/agglomerative_bundling.h [new file with mode: 0644]
lib/mingle/edge_bundling.c [new file with mode: 0644]
lib/mingle/edge_bundling.h [new file with mode: 0644]
lib/mingle/ink.c [new file with mode: 0644]
lib/mingle/ink.h [new file with mode: 0644]
lib/mingle/minglelib.vcproj [new file with mode: 0644]
lib/mingle/nearest_neighbor_graph.c [new file with mode: 0644]
lib/mingle/nearest_neighbor_graph.h [new file with mode: 0644]
lib/mingle/nearest_neighbor_graph_ann.cpp [new file with mode: 0644]
lib/mingle/nearest_neighbor_graph_ann.h [new file with mode: 0644]

index 5b0a1dece1412f94678cf9ff8e67a758db34da78..01c6472eb0438466a75677661c9f4a5ba7e12358 100644 (file)
@@ -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 (file)
index 0000000..2ffa7cf
--- /dev/null
@@ -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 (file)
index 0000000..6a8fea6
--- /dev/null
@@ -0,0 +1,203 @@
+<?xml version="1.0" encoding="Windows-1252"?>\r
+<VisualStudioProject\r
+       ProjectType="Visual C++"\r
+       Version="9.00"\r
+       Name="mingle"\r
+       ProjectGUID="{A6DF0D74-E4D1-4EF4-A642-59B0CF3E74BB}"\r
+       RootNamespace="mingle"\r
+       Keyword="Win32Proj"\r
+       TargetFrameworkVersion="196613"\r
+       >\r
+       <Platforms>\r
+               <Platform\r
+                       Name="Win32"\r
+               />\r
+       </Platforms>\r
+       <ToolFiles>\r
+       </ToolFiles>\r
+       <Configurations>\r
+               <Configuration\r
+                       Name="Debug|Win32"\r
+                       OutputDirectory="Debug"\r
+                       IntermediateDirectory="Debug"\r
+                       ConfigurationType="1"\r
+                       CharacterSet="1"\r
+                       >\r
+                       <Tool\r
+                               Name="VCPreBuildEventTool"\r
+                       />\r
+                       <Tool\r
+                               Name="VCCustomBuildTool"\r
+                       />\r
+                       <Tool\r
+                               Name="VCXMLDataGeneratorTool"\r
+                       />\r
+                       <Tool\r
+                               Name="VCWebServiceProxyGeneratorTool"\r
+                       />\r
+                       <Tool\r
+                               Name="VCMIDLTool"\r
+                       />\r
+                       <Tool\r
+                               Name="VCCLCompilerTool"\r
+                               Optimization="0"\r
+                               AdditionalIncludeDirectories="../../lib/ingraphs;../../lib/common;../../lib/cdt;../../lib/cgraph;../../lib/sparse;../../lib/sfdpgen;../../;../../lib/mingle"\r
+                               PreprocessorDefinitions="WIN32;_DEBUG;_LIB;HAVE_CONFIG_H"\r
+                               MinimalRebuild="true"\r
+                               BasicRuntimeChecks="3"\r
+                               RuntimeLibrary="3"\r
+                               UsePrecompiledHeader="0"\r
+                               WarningLevel="3"\r
+                               DebugInformationFormat="4"\r
+                       />\r
+                       <Tool\r
+                               Name="VCManagedResourceCompilerTool"\r
+                       />\r
+                       <Tool\r
+                               Name="VCResourceCompilerTool"\r
+                       />\r
+                       <Tool\r
+                               Name="VCPreLinkEventTool"\r
+                       />\r
+                       <Tool\r
+                               Name="VCLinkerTool"\r
+                               OutputFile="c:/graphviz-ms/bin/$(ProjectName).exe"\r
+                               LinkIncremental="2"\r
+                               AdditionalLibraryDirectories="&quot;C:\graphviz-ms\bin&quot;;C:\gtk\lib"\r
+                               GenerateDebugInformation="true"\r
+                               SubSystem="1"\r
+                               TargetMachine="1"\r
+                       />\r
+                       <Tool\r
+                               Name="VCALinkTool"\r
+                       />\r
+                       <Tool\r
+                               Name="VCManifestTool"\r
+                       />\r
+                       <Tool\r
+                               Name="VCXDCMakeTool"\r
+                       />\r
+                       <Tool\r
+                               Name="VCBscMakeTool"\r
+                       />\r
+                       <Tool\r
+                               Name="VCFxCopTool"\r
+                       />\r
+                       <Tool\r
+                               Name="VCAppVerifierTool"\r
+                       />\r
+                       <Tool\r
+                               Name="VCPostBuildEventTool"\r
+                       />\r
+               </Configuration>\r
+               <Configuration\r
+                       Name="Release|Win32"\r
+                       OutputDirectory="Release"\r
+                       IntermediateDirectory="Release"\r
+                       ConfigurationType="1"\r
+                       CharacterSet="1"\r
+                       WholeProgramOptimization="1"\r
+                       >\r
+                       <Tool\r
+                               Name="VCPreBuildEventTool"\r
+                       />\r
+                       <Tool\r
+                               Name="VCCustomBuildTool"\r
+                       />\r
+                       <Tool\r
+                               Name="VCXMLDataGeneratorTool"\r
+                       />\r
+                       <Tool\r
+                               Name="VCWebServiceProxyGeneratorTool"\r
+                       />\r
+                       <Tool\r
+                               Name="VCMIDLTool"\r
+                       />\r
+                       <Tool\r
+                               Name="VCCLCompilerTool"\r
+                               Optimization="2"\r
+                               EnableIntrinsicFunctions="true"\r
+                               AdditionalIncludeDirectories="../../lib/ingraphs;../../lib/common;../../lib/cdt;../../lib/cgraph;../../lib/sparse;../../lib/sfdpgen;../../;../../lib/mingle"\r
+                               PreprocessorDefinitions="WIN32;_DEBUG;_LIB;HAVE_CONFIG_H"\r
+                               RuntimeLibrary="2"\r
+                               EnableFunctionLevelLinking="true"\r
+                               UsePrecompiledHeader="0"\r
+                               WarningLevel="3"\r
+                               DebugInformationFormat="3"\r
+                       />\r
+                       <Tool\r
+                               Name="VCManagedResourceCompilerTool"\r
+                       />\r
+                       <Tool\r
+                               Name="VCResourceCompilerTool"\r
+                       />\r
+                       <Tool\r
+                               Name="VCPreLinkEventTool"\r
+                       />\r
+                       <Tool\r
+                               Name="VCLinkerTool"\r
+                               OutputFile="c:/graphviz-ms/bin/$(ProjectName).exe"\r
+                               LinkIncremental="2"\r
+                               AdditionalLibraryDirectories="&quot;C:\graphviz-ms\bin&quot;;C:\gtk\lib"\r
+                               GenerateDebugInformation="true"\r
+                               SubSystem="1"\r
+                               OptimizeReferences="2"\r
+                               EnableCOMDATFolding="2"\r
+                               TargetMachine="1"\r
+                       />\r
+                       <Tool\r
+                               Name="VCALinkTool"\r
+                       />\r
+                       <Tool\r
+                               Name="VCManifestTool"\r
+                       />\r
+                       <Tool\r
+                               Name="VCXDCMakeTool"\r
+                       />\r
+                       <Tool\r
+                               Name="VCBscMakeTool"\r
+                       />\r
+                       <Tool\r
+                               Name="VCFxCopTool"\r
+                       />\r
+                       <Tool\r
+                               Name="VCAppVerifierTool"\r
+                       />\r
+                       <Tool\r
+                               Name="VCPostBuildEventTool"\r
+                       />\r
+               </Configuration>\r
+       </Configurations>\r
+       <References>\r
+       </References>\r
+       <Files>\r
+               <Filter\r
+                       Name="Source Files"\r
+                       Filter="cpp;c;cc;cxx;def;odl;idl;hpj;bat;asm;asmx"\r
+                       UniqueIdentifier="{4FC737F1-C7A5-4376-A066-2A32D752A2FF}"\r
+                       >\r
+                       <File\r
+                               RelativePath=".\minglemain.c"\r
+                               >\r
+                       </File>\r
+               </Filter>\r
+               <Filter\r
+                       Name="Header Files"\r
+                       Filter="h;hpp;hxx;hm;inl;inc;xsd"\r
+                       UniqueIdentifier="{93995380-89BD-4b04-88EB-625FBE52EBFB}"\r
+                       >\r
+               </Filter>\r
+               <Filter\r
+                       Name="Resource Files"\r
+                       Filter="rc;ico;cur;bmp;dlg;rc2;rct;bin;rgs;gif;jpg;jpeg;jpe;resx;tiff;tif;png;wav"\r
+                       UniqueIdentifier="{67DA6AB6-F800-4c08-8B7A-83BB121AAD01}"\r
+                       >\r
+               </Filter>\r
+               <File\r
+                       RelativePath=".\ReadMe.txt"\r
+                       >\r
+               </File>\r
+       </Files>\r
+       <Globals>\r
+       </Globals>\r
+</VisualStudioProject>\r
diff --git a/cmd/mingle/minglemain.c b/cmd/mingle/minglemain.c
new file mode 100644 (file)
index 0000000..3e0dae0
--- /dev/null
@@ -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 <cgraph.h>
+#include <ingraphs.h>
+#ifdef HAVE_GETOPT_H
+#include <getopt.h>
+#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 <options> <file>\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;
+}
index 9b78ffe9680f409e23c1d3900ff282c57eea1b6b..d5b13185eff5ada8ab18bc9a0db78eacb034c5cf 100644 (file)
@@ -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"
index c2088ec690c5eaf7f1533158e193711a2f47919b..357739505e67ba795792d9782619bc2c533060dd 100644 (file)
@@ -1,5 +1,4 @@
 Microsoft Visual Studio Solution File, Format Version 10.00\r
-#
 # Visual Studio 2008\r
 Project("{8BC9CEB8-8B4A-11D0-8D11-00A0C91BC942}") = "Pathplan", "lib\pathplan\Pathplan.vcproj", "{BD347753-A09D-48B4-8752-F1D8D9CF235D}"\r
        ProjectSection(ProjectDependencies) = postProject\r
@@ -262,6 +261,20 @@ Project("{8BC9CEB8-8B4A-11D0-8D11-00A0C91BC942}") = "gvedit", "cmd\gvedit\gvedit
                {15229511-9F6C-48A5-9194-660CA6492563} = {15229511-9F6C-48A5-9194-660CA6492563}\r
        EndProjectSection\r
 EndProject\r
+Project("{8BC9CEB8-8B4A-11D0-8D11-00A0C91BC942}") = "minglelib", "lib\mingle\minglelib.vcproj", "{70575BD2-A532-41B8-9A17-B9876E992A84}"\r
+EndProject\r
+Project("{8BC9CEB8-8B4A-11D0-8D11-00A0C91BC942}") = "mingle", "cmd\mingle\mingle.vcproj", "{A6DF0D74-E4D1-4EF4-A642-59B0CF3E74BB}"\r
+       ProjectSection(ProjectDependencies) = postProject\r
+               {C0663A08-F276-4DD6-B17C-E501EE066F7C} = {C0663A08-F276-4DD6-B17C-E501EE066F7C}\r
+               {15229511-9F6C-48A5-9194-660CA6492563} = {15229511-9F6C-48A5-9194-660CA6492563}\r
+               {D6CEB142-BF8E-471C-AE16-4300F2D5DEDA} = {D6CEB142-BF8E-471C-AE16-4300F2D5DEDA}\r
+               {B76BCE8C-63CC-4A99-88B5-D621D563E699} = {B76BCE8C-63CC-4A99-88B5-D621D563E699}\r
+               {C5CEB09E-79AF-4091-ACCF-D28EC3D7D90F} = {C5CEB09E-79AF-4091-ACCF-D28EC3D7D90F}\r
+               {443EB1A7-C634-4292-9F2D-C752BB2BF40F} = {443EB1A7-C634-4292-9F2D-C752BB2BF40F}\r
+               {70575BD2-A532-41B8-9A17-B9876E992A84} = {70575BD2-A532-41B8-9A17-B9876E992A84}\r
+               {D6FD0DE5-5305-458E-8CA5-FCA4B8E05B04} = {D6FD0DE5-5305-458E-8CA5-FCA4B8E05B04}\r
+       EndProjectSection\r
+EndProject\r
 Global\r
        GlobalSection(SolutionConfigurationPlatforms) = preSolution\r
                Debug|Win32 = Debug|Win32\r
@@ -508,6 +521,14 @@ Global
                {33E97266-F213-3639-BCAB-2F3F95E15B16}.Debug|Win32.Build.0 = Debug|Win32\r
                {33E97266-F213-3639-BCAB-2F3F95E15B16}.Release|Win32.ActiveCfg = Release|Win32\r
                {33E97266-F213-3639-BCAB-2F3F95E15B16}.Release|Win32.Build.0 = Release|Win32\r
+               {70575BD2-A532-41B8-9A17-B9876E992A84}.Debug|Win32.ActiveCfg = Debug|Win32\r
+               {70575BD2-A532-41B8-9A17-B9876E992A84}.Debug|Win32.Build.0 = Debug|Win32\r
+               {70575BD2-A532-41B8-9A17-B9876E992A84}.Release|Win32.ActiveCfg = Release|Win32\r
+               {70575BD2-A532-41B8-9A17-B9876E992A84}.Release|Win32.Build.0 = Release|Win32\r
+               {A6DF0D74-E4D1-4EF4-A642-59B0CF3E74BB}.Debug|Win32.ActiveCfg = Debug|Win32\r
+               {A6DF0D74-E4D1-4EF4-A642-59B0CF3E74BB}.Debug|Win32.Build.0 = Debug|Win32\r
+               {A6DF0D74-E4D1-4EF4-A642-59B0CF3E74BB}.Release|Win32.ActiveCfg = Release|Win32\r
+               {A6DF0D74-E4D1-4EF4-A642-59B0CF3E74BB}.Release|Win32.Build.0 = Release|Win32\r
        EndGlobalSection\r
        GlobalSection(SolutionProperties) = preSolution\r
                HideSolutionNode = FALSE\r
index 792a1342471cf8ae78a9ab443907b6ee481243ac..6a44bf7cf361c16d86b025d4717c738d1547f984 100644 (file)
@@ -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 (file)
index 0000000..a92a4a5
--- /dev/null
@@ -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 (file)
index 0000000..61de39f
--- /dev/null
@@ -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, &current_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 (file)
index 0000000..0cbfd8e
--- /dev/null
@@ -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 (file)
index 0000000..fa114b5
--- /dev/null
@@ -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 <time.h>
+#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 (file)
index 0000000..7d6018d
--- /dev/null
@@ -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 <SparseMatrix.h>
+
+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 (file)
index 0000000..9bee033
--- /dev/null
@@ -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 <math.h>
+#include <stdlib.h>
+#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; i<npoints; i++) {
+    diff_x = points[i].x-meeting.x;
+    diff_y = points[i].y-meeting.y;
+    len = sqrt(diff_x*diff_x+diff_y*diff_y);
+    sum += len;
+    cos_theta = (diff_x0*diff_x + diff_y0*diff_y)/MAX((len*len0),0.00001);
+    cos_max = MAX(cos_max, cos_theta);
+  }
+
+  // distance of single line from 'meeting' to 'end'
+  return sum*(cos_max + angle_param);/* straight line gives angle_param - 1, turning angle of 180 degree gives angle_param + 1 */
+}
+double sumLengths(point_t* points, int npoints, point_t end, point_t meeting) 
+{
+  int i;
+  double sum = 0;
+  double diff_x, diff_y;
+
+    // distance form each of 'points' till 'meeting'
+  for (i=0; i<npoints; i++) {
+        diff_x = points[i].x-meeting.x;
+        diff_y = points[i].y-meeting.y;
+        sum += sqrt(diff_x*diff_x+diff_y*diff_y);
+  }
+    // distance of single line from 'meeting' to 'end'
+  diff_x = end.x-meeting.x;
+  diff_y = end.y-meeting.y;
+  sum += sqrt(diff_x*diff_x+diff_y*diff_y);
+  return sum;
+}
+
+/* bestInk:
+ */
+double bestInk(point_t* points, int npoints, point_t begin, point_t end, double prec, point_t *meet, real angle_param)
+{
+  point_t first, second, third, fourth, diff, meeting;
+  double value1, value2, value3, value4;
+
+  first = begin;
+  fourth = end;
+
+  do {
+    diff = subPoint(fourth,first);
+    second = addPoint(first,scalePoint(diff,1.0/3.0));
+    third = addPoint(first,scalePoint(diff,2.0/3.0));
+
+    if (angle_param < 1){
+      value1 = sumLengths(points, npoints, end, first);
+      value2 = sumLengths(points, npoints, end, second);
+      value3 = sumLengths(points, npoints, end, third);
+      value4 = sumLengths(points, npoints, end, fourth);
+    } else {
+      value1 = sumLengths_avoid_bad_angle(points, npoints, end, first, angle_param);
+      value2 = sumLengths_avoid_bad_angle(points, npoints, end, second, angle_param);
+      value3 = sumLengths_avoid_bad_angle(points, npoints, end, third, angle_param);
+      value4 = sumLengths_avoid_bad_angle(points, npoints, end, fourth, angle_param);
+    }
+
+    if (value1<value2) {
+      if (value1<value3) {
+        if (value1<value4) {
+          // first is smallest
+          fourth = second;
+        }
+        else {
+          // fourth is smallest
+          first = third;
+        }
+      }
+      else {
+        if (value3<value4) {
+          // third is smallest
+          first = second;
+        }
+        else {
+          // fourth is smallest
+          first = third;
+        }
+      }
+    }
+    else {
+      if (value2<value3) {
+        if (value2<value4) {
+          // second is smallest
+          fourth = third;
+        }
+        else {
+          // fourth is smallest
+          first = third;
+        }
+      }
+      else {
+        if (value3<value4) {
+          // third is smallest
+          first = second;
+        }
+        else {
+          // fourth is smallest
+          first = third;
+        }
+      }
+    }
+  } while (fabs(value1-value4)/(MIN(value1,value4)+1e-10)>prec && 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 (file)
index 0000000..09be9f5
--- /dev/null
@@ -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 <edge_bundling.h>
+
+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 (file)
index 0000000..3132ffe
--- /dev/null
@@ -0,0 +1,213 @@
+<?xml version="1.0" encoding="Windows-1252"?>\r
+<VisualStudioProject\r
+       ProjectType="Visual C++"\r
+       Version="9.00"\r
+       Name="minglelib"\r
+       ProjectGUID="{70575BD2-A532-41B8-9A17-B9876E992A84}"\r
+       RootNamespace="minglelib"\r
+       Keyword="Win32Proj"\r
+       TargetFrameworkVersion="196613"\r
+       >\r
+       <Platforms>\r
+               <Platform\r
+                       Name="Win32"\r
+               />\r
+       </Platforms>\r
+       <ToolFiles>\r
+       </ToolFiles>\r
+       <Configurations>\r
+               <Configuration\r
+                       Name="Debug|Win32"\r
+                       OutputDirectory="Debug"\r
+                       IntermediateDirectory="Debug"\r
+                       ConfigurationType="4"\r
+                       CharacterSet="1"\r
+                       >\r
+                       <Tool\r
+                               Name="VCPreBuildEventTool"\r
+                       />\r
+                       <Tool\r
+                               Name="VCCustomBuildTool"\r
+                       />\r
+                       <Tool\r
+                               Name="VCXMLDataGeneratorTool"\r
+                       />\r
+                       <Tool\r
+                               Name="VCWebServiceProxyGeneratorTool"\r
+                       />\r
+                       <Tool\r
+                               Name="VCMIDLTool"\r
+                       />\r
+                       <Tool\r
+                               Name="VCCLCompilerTool"\r
+                               Optimization="0"\r
+                               AdditionalIncludeDirectories="&quot;$(SolutionDir)/lib/sparse&quot;;&quot;$(SolutionDir)/lib/sfpdpgen&quot;;&quot;$(SolutionDir)/lib/common&quot;;&quot;$(SolutionDir)/lib/gvc&quot;;&quot;$(SolutionDir)/lib/cdt&quot;;&quot;$(SolutionDir)/lib/pathplan&quot;;&quot;$(SolutionDir)/lib/cgraph&quot;;&quot;$(SolutionDir)&quot;;&quot;C:\graphviz-ms\graphviz2\lib\mingle&quot;"\r
+                               PreprocessorDefinitions="WIN32;_DEBUG;_LIB;HAVE_CONFIG_H"\r
+                               MinimalRebuild="true"\r
+                               BasicRuntimeChecks="3"\r
+                               RuntimeLibrary="3"\r
+                               UsePrecompiledHeader="0"\r
+                               WarningLevel="3"\r
+                               DebugInformationFormat="4"\r
+                       />\r
+                       <Tool\r
+                               Name="VCManagedResourceCompilerTool"\r
+                       />\r
+                       <Tool\r
+                               Name="VCResourceCompilerTool"\r
+                       />\r
+                       <Tool\r
+                               Name="VCPreLinkEventTool"\r
+                       />\r
+                       <Tool\r
+                               Name="VCLibrarianTool"\r
+                       />\r
+                       <Tool\r
+                               Name="VCALinkTool"\r
+                       />\r
+                       <Tool\r
+                               Name="VCXDCMakeTool"\r
+                       />\r
+                       <Tool\r
+                               Name="VCBscMakeTool"\r
+                       />\r
+                       <Tool\r
+                               Name="VCFxCopTool"\r
+                       />\r
+                       <Tool\r
+                               Name="VCPostBuildEventTool"\r
+                       />\r
+               </Configuration>\r
+               <Configuration\r
+                       Name="Release|Win32"\r
+                       OutputDirectory="Release"\r
+                       IntermediateDirectory="Release"\r
+                       ConfigurationType="4"\r
+                       CharacterSet="1"\r
+                       WholeProgramOptimization="1"\r
+                       >\r
+                       <Tool\r
+                               Name="VCPreBuildEventTool"\r
+                       />\r
+                       <Tool\r
+                               Name="VCCustomBuildTool"\r
+                       />\r
+                       <Tool\r
+                               Name="VCXMLDataGeneratorTool"\r
+                       />\r
+                       <Tool\r
+                               Name="VCWebServiceProxyGeneratorTool"\r
+                       />\r
+                       <Tool\r
+                               Name="VCMIDLTool"\r
+                       />\r
+                       <Tool\r
+                               Name="VCCLCompilerTool"\r
+                               Optimization="2"\r
+                               EnableIntrinsicFunctions="true"\r
+                               AdditionalIncludeDirectories="&quot;$(SolutionDir)/lib/sparse&quot;;&quot;$(SolutionDir)/lib/sfpdpgen&quot;;&quot;$(SolutionDir)/lib/common&quot;;&quot;$(SolutionDir)/lib/gvc&quot;;&quot;$(SolutionDir)/lib/cdt&quot;;&quot;$(SolutionDir)/lib/pathplan&quot;;&quot;$(SolutionDir)/lib/cgraph&quot;;&quot;$(SolutionDir)&quot;;&quot;C:\graphviz-ms\graphviz2\lib\mingle&quot;"\r
+                               PreprocessorDefinitions="WIN32;_DEBUG;_LIB;HAVE_CONFIG_H"\r
+                               RuntimeLibrary="2"\r
+                               EnableFunctionLevelLinking="true"\r
+                               UsePrecompiledHeader="0"\r
+                               WarningLevel="3"\r
+                               DebugInformationFormat="3"\r
+                       />\r
+                       <Tool\r
+                               Name="VCManagedResourceCompilerTool"\r
+                       />\r
+                       <Tool\r
+                               Name="VCResourceCompilerTool"\r
+                       />\r
+                       <Tool\r
+                               Name="VCPreLinkEventTool"\r
+                       />\r
+                       <Tool\r
+                               Name="VCLibrarianTool"\r
+                       />\r
+                       <Tool\r
+                               Name="VCALinkTool"\r
+                       />\r
+                       <Tool\r
+                               Name="VCXDCMakeTool"\r
+                       />\r
+                       <Tool\r
+                               Name="VCBscMakeTool"\r
+                       />\r
+                       <Tool\r
+                               Name="VCFxCopTool"\r
+                       />\r
+                       <Tool\r
+                               Name="VCPostBuildEventTool"\r
+                       />\r
+               </Configuration>\r
+       </Configurations>\r
+       <References>\r
+       </References>\r
+       <Files>\r
+               <Filter\r
+                       Name="Source Files"\r
+                       Filter="cpp;c;cc;cxx;def;odl;idl;hpj;bat;asm;asmx"\r
+                       UniqueIdentifier="{4FC737F1-C7A5-4376-A066-2A32D752A2FF}"\r
+                       >\r
+                       <File\r
+                               RelativePath=".\agglomerative_bundling.c"\r
+                               >\r
+                       </File>\r
+                       <File\r
+                               RelativePath=".\edge_bundling.c"\r
+                               >\r
+                       </File>\r
+                       <File\r
+                               RelativePath=".\ink.c"\r
+                               >\r
+                       </File>\r
+                       <File\r
+                               RelativePath=".\nearest_neighbor_graph.c"\r
+                               >\r
+                       </File>\r
+                       <File\r
+                               RelativePath=".\nearest_neighbor_graph_ann.cpp"\r
+                               >\r
+                       </File>\r
+               </Filter>\r
+               <Filter\r
+                       Name="Header Files"\r
+                       Filter="h;hpp;hxx;hm;inl;inc;xsd"\r
+                       UniqueIdentifier="{93995380-89BD-4b04-88EB-625FBE52EBFB}"\r
+                       >\r
+                       <File\r
+                               RelativePath=".\agglomerative_bundling.h"\r
+                               >\r
+                       </File>\r
+                       <File\r
+                               RelativePath=".\edge_bundling.h"\r
+                               >\r
+                       </File>\r
+                       <File\r
+                               RelativePath=".\ink.h"\r
+                               >\r
+                       </File>\r
+                       <File\r
+                               RelativePath=".\nearest_neighbor_graph.h"\r
+                               >\r
+                       </File>\r
+                       <File\r
+                               RelativePath=".\nearest_neighbor_graph_ann.h"\r
+                               >\r
+                       </File>\r
+               </Filter>\r
+               <Filter\r
+                       Name="Resource Files"\r
+                       Filter="rc;ico;cur;bmp;dlg;rc2;rct;bin;rgs;gif;jpg;jpeg;jpe;resx;tiff;tif;png;wav"\r
+                       UniqueIdentifier="{67DA6AB6-F800-4c08-8B7A-83BB121AAD01}"\r
+                       >\r
+               </Filter>\r
+               <File\r
+                       RelativePath=".\ReadMe.txt"\r
+                       >\r
+               </File>\r
+       </Files>\r
+       <Globals>\r
+       </Globals>\r
+</VisualStudioProject>\r
diff --git a/lib/mingle/nearest_neighbor_graph.c b/lib/mingle/nearest_neighbor_graph.c
new file mode 100644 (file)
index 0000000..02707c3
--- /dev/null
@@ -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 (file)
index 0000000..3a042eb
--- /dev/null
@@ -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 (file)
index 0000000..02db85d
--- /dev/null
@@ -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 <cstdlib>                                             // C standard library
+#include <cstdio>                                              // C I/O (for sscanf)
+#include <cstring>                                             // string manipulation
+#include <fstream>                                             // file I/O
+#include <ANN/ANN.h>                                   // 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 (file)
index 0000000..f30ee5b
--- /dev/null
@@ -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 */