before new shortest path graph layout cost function

diff --git a/anneal.c b/anneal.c
index 32c06d99cef7faef2a5cd47a8d85a809d459b130..fd725f1844d577488a2e9d154fc03c155f9bf913 100644 (file)
--- a/anneal.c
+++ b/anneal.c
@@ -37,9 +37,10 @@
    *     vd, edge PQ      vd      3
    *                             l
    *
    *     vd, edge PQ      vd      3
    *                             l
    *

-   *  By symmetry, this calculation gives the same answer
-   *  with R and S exchanged.  Looking at the projection in
-   *  the RMS plane:
+   *  (The dimensions of this are those of F_vd.)
+   *
+   *  By symmetry, this calculation gives the same answer with R and S
+   *  exchanged.  Looking at the projection in the RMS plane:

    *
    *
    *                          S'
    *
    *
    *                          S'


@@ -62,9 +63,10 @@
    *
    *  In practice to avoid division by zero we'll add epsilon to l^3
    *  and the huge energy ought then to be sufficient for the model to
    *
    *  In practice to avoid division by zero we'll add epsilon to l^3
    *  and the huge energy ought then to be sufficient for the model to

-   *  avoid being close to R=S.  */
+   *  avoid being close to R=S.
+   */

 
 

-static double hypotD(const double p[], const double q[]) {
+double hypotD(const double p[D3], const double q[D3]) {

   int k;
   double pq[D3];
   gsl_vector v;
   int k;
   double pq[D3];
   gsl_vector v;


@@ -78,10 +80,16 @@ static double hypotD(const double p[], const double q[]) {
   return gsl_blas_snrm2(&v);
 }
 
   return gsl_blas_snrm2(&v);
 }
 

+double hypotD2(const double p[D3], const double q[D3]) {
+  double d2= 0;
+  K d2= ffsqa(p[k] - q[k], d2);
+  return d2;
+}
+

 #ifdef FP_FAST_FMA
 # define ffsqa(factor,term) fma((factor),(factor),(term))
 #else
 #ifdef FP_FAST_FMA
 # define ffsqa(factor,term) fma((factor),(factor),(term))
 #else

-# define ffsqu(factor,term) ((factor)*(factor)+(term))
+# define ffsqa(factor,term) ((factor)*(factor)+(term))

 #endif
 
 static const l3_epsison= 1e-6;
 #endif
 
 static const l3_epsison= 1e-6;


@@ -99,7 +107,7 @@ static double energy_function(const double vertices[N][D3]) {
     K m[k]= (vertices[pi][k] + vertices[qi][k]) * 0.5;
     K mprime[k]= (vertices[ri][k] + vertices[si][k]) * 0.5;
     b= hypotD(vertices[pi], vertices[qi]);
     K m[k]= (vertices[pi][k] + vertices[qi][k]) * 0.5;
     K mprime[k]= (vertices[ri][k] + vertices[si][k]) * 0.5;
     b= hypotD(vertices[pi], vertices[qi]);

-    d2= 0; K d2= ffsqa(m[k] - mprime[k], d2);
+    d2= hypotD2(m, mprime);

     l= hypotD(vertices[ri][k] - vertices[si][k]);
     l3 = l*l*l + l3_epsilon;
 
     l= hypotD(vertices[ri][k] - vertices[si][k]);
     l3 = l*l*l + l3_epsilon;
 

diff --git a/bgl.cpp b/bgl.cpp
index 774be825a1cc53bc017884e94876ddd49c134984..775fb12141c5d806bc3d98f36f9db2600eec7564 100644 (file)
--- a/bgl.cpp
+++ b/bgl.cpp
@@ -81,28 +81,44 @@ namespace boost {
 
 }
 
 
 }
 

-struct VertexIndexMap;
-
-namespace boost {
-  struct property_traits<VertexIndexMap> {
-    // Concept Readable Property Map:
-    typedef int value_type, reference, key_type;
-    category 
-class Boost
-  }};
-
-void single_source_shortest_paths(int v1,
-                                 const double edge_weights[/*f*/],
-                                 ) {
+static void single_source_shortest_paths(int v1,
+                                        const double edge_weights[/*f*/],
+                                        double vertex_distances[/*v*/]) {

   boost::dijkstra_shortest_paths
     (g, v1,
      weight_map(edge_weights).
      vertex_index_map(identity_property_map()).
   boost::dijkstra_shortest_paths
     (g, v1,
      weight_map(edge_weights).
      vertex_index_map(identity_property_map()).

-     
+     distance_map(vertex_distances));
+}

     
     

-void all_pairs_shortest_paths(const Layout *g) {
-  
+double graph_layout_energy(const Layout *g) {
+  /* For each (vi,vj) computes shortest path pij = vi..vj along edges,
+   * and actual distance dij = vi-vj.
+   *
+   * Energy contribution is proportional to
+   *
+   *     pij - dij   
+   *     ---------  .  Sigma                                 length
+   *           2        e member (edges(vi) union edges(vj)        e
+   *        dij
+   */
+  double edge_weights[N*V6], vertex_distances[N];
+  int v1, e, f;
+
+  FOR_VEDGE_X(v1,e,
+             f= v1 | e << ESHIFT,
+             edge_weights[f]= NaN)
+    edge_weights[f]= hypotD(g.v[v1], g.v[v2]);
+

   FOR_VERTEX(v1) {
   FOR_VERTEX(v1) {

+    single_source_shortest_paths(v1, edge_weights, vertex_distances);
+    FOR_VERTEX(v2) {
+      
+            
+      
+    
+    { ) {
+    

        
  0);
     
        
  0);
     

diff --git a/mgraph.h b/mgraph.h
index 9f02e88c418438c682321d6d30ea916c0ac0ca9e..ae341fae532786f45a2b53c779802eeb47ab0f25 100644 (file)
--- a/mgraph.h
+++ b/mgraph.h
@@ -80,9 +80,16 @@
 extern int edge_end2(unsigned v1, int e);
 #define EDGE_END2 edge_end2
 
 extern int edge_end2(unsigned v1, int e);
 #define EDGE_END2 edge_end2
 

+#define FOR_VEDGE_X(v1,e,init,otherwise)       \
+  FOR_VPEDGE((v1),(e))                         \
+    if (((v2)= EDGE_END2((v1),(e)),            \
+        (init),                                \
+        (v2)) < 0) { otherwise } else
+
+#define NOTHING ((void)0)
+

 #define FOR_VEDGE(v1,e) \
 #define FOR_VEDGE(v1,e) \

-  FOR_VPEDGE((v1),(e))
-    if (((v2)= EDGE_END2((v1),(e))) < 0) ; else
+  FOR_VEDGE_X(v1,e,NOTHING,NOTHING)

 
 #define FOR_EDGE(v1,e,v2) \
   FOR_VERTEX((v1)) \
 
 #define FOR_EDGE(v1,e,v2) \
   FOR_VERTEX((v1)) \


@@ -97,4 +104,7 @@ typedef struct {
   double v[N][D3];
 } Layout;
 
   double v[N][D3];
 } Layout;
 

+double hypotD(const double p[D3], const double q[D3]);
+double hypotD2(const double p[D3], const double q[D3]);
+

 #endif /*MGRAPH_H*/
 #endif /*MGRAPH_H*/