before new shortest path graph layout cost function

[moebius2.git] / anneal.c

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
    *
-   *  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'
@@ -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
-   *  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;
@@ -78,10 +80,16 @@ static double hypotD(const double p[], const double q[]) {
   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
-# define ffsqu(factor,term) ((factor)*(factor)+(term))
+# define ffsqa(factor,term) ((factor)*(factor)+(term))
 #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]);
-    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;