From: Ian Jackson <ian@davenant.relativity.greenend.org.uk>
Date: Sat, 29 Dec 2007 18:07:10 +0000 (+0000)
Subject: done graph layout cost function core, now reorganise
X-Git-Url: https://www.chiark.greenend.org.uk/ucgi/~ian/git?a=commitdiff_plain;h=1050281a7d39f56431d9f36a6af4b19774419b75;p=moebius2.git

done graph layout cost function core, now reorganise
---

diff --git a/anneal.c b/anneal.c
index fd725f1..cb4793a 100644
--- a/anneal.c
+++ b/anneal.c
@@ -86,11 +86,17 @@ double hypotD2(const double p[D3], const double q[D3]) {
   return d2;
 }
 
+double hypotD2plus(const double p[D3], const double q[D3], double d2) {
+  K d2= ffsqa(p[k] - q[k], d2);
+  return d2;
+}
+
 #ifdef FP_FAST_FMA
-# define ffsqa(factor,term) fma((factor),(factor),(term))
+# define fma_fast fma
 #else
-# define ffsqa(factor,term) ((factor)*(factor)+(term))
+# define fma_fast(f1,f2,t) ((f1)*(f2)+(t))
 #endif
+#define ffsqa(factor,term) fma_fast((factor),(factor),(term))
 
 static const l3_epsison= 1e-6;
 
diff --git a/bgl.cpp b/bgl.cpp
index 775fb12..eb7e096 100644
--- a/bgl.cpp
+++ b/bgl.cpp
@@ -91,18 +91,24 @@ static void single_source_shortest_paths(int v1,
      distance_map(vertex_distances));
 }
     
-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.
+double graph_layout_cost(const Layout *g, const double vertex_areas[N]) {
+  /* For each (vi,vj) computes shortest path s_ij = |vi..vj|
+   * along edges, and actual distance d_ij = |vi-vj|.
+   *
+   * We will also use the `vertex areas': for each vertex vi the
+   * vertex area a_vi is the mean area of the incident triangles.
+   * This is computed elsewhere.
    *
    * Energy contribution is proportional to
    *
-   *     pij - dij   
-   *     ---------  .  Sigma                                 length
-   *           2        e member (edges(vi) union edges(vj)        e
-   *        dij
+   *               -4          2
+   *    a  a   .  d   . [ (s/d)  - 1 ]
+   *     vi vj
+   *
+   * (In practice we compute d^2+epsilon and use it for the
+   *  divisions, to avoid division by zero.)
    */
-  double edge_weights[N*V6], vertex_distances[N];
+  double edge_weights[N*V6], vertex_distances[N], total_cost;
   int v1, e, f;
 
   FOR_VEDGE_X(v1,e,
@@ -111,23 +117,15 @@ double graph_layout_energy(const Layout *g) {
     edge_weights[f]= hypotD(g.v[v1], g.v[v2]);
 
   FOR_VERTEX(v1) {
+    double a1= vertex_areas[v1];
     single_source_shortest_paths(v1, edge_weights, vertex_distances);
     FOR_VERTEX(v2) {
-      
-	     
-      
-    
-    { ) {
-    
-       
- 0);
-    
-			    /* weight_map(). ? */
-			    /* 
-
-			    
-			    predecessor_map().
-			    distance_map()
-
-
-
+      double a2= vertex_areas[v2];
+      double d2= hypotD2plus(g->v[v1],g->v[v2], d2_epsilon);
+      double sd= vertex_distances[v2] / d2;
+      double sd2= sd*sd;
+      total_cost= fma_fast(a1*a2, (sd2 - 1.0)/(d2*d2), total_cost);
+    }
+  }
+  return total_cost;
+}