chiark / gitweb /
reorganisation is complete, need to implement various things, make it compile, etc.
authorIan Jackson <ian@davenant.relativity.greenend.org.uk>
Sat, 29 Dec 2007 20:42:39 +0000 (20:42 +0000)
committerIan Jackson <ian@davenant.relativity.greenend.org.uk>
Sat, 29 Dec 2007 20:42:39 +0000 (20:42 +0000)
anneal.c [deleted file]
bgl.cpp
bgl.h [new file with mode: 0644]
common.c [new file with mode: 0644]
common.h [new file with mode: 0644]
energy.c [new file with mode: 0644]
mgraph.c
mgraph.h

diff --git a/anneal.c b/anneal.c
deleted file mode 100644 (file)
index cb4793a..0000000
--- a/anneal.c
+++ /dev/null
@@ -1,131 +0,0 @@
-/*
- * We try to find an optimal triangle grid
- */
-
-  /* Energy has the following parts right now:
-   */
-  /*
-   *  Edgewise expected vertex displacement
-   *
-   *
-   *                       Q `-_
-   *              / |    `-_
-   *   R' - _ _ _/_ |       `-.
-   *    .       /   M - - - - - S
-   *    .      /    |      _,-'
-   *    .     /     |  _,-'
-   *    .    /    , P '
-   *    .   /  ,-'
-   *    .  /,-'
-   *    . /'
-   *            R
-   *
-   *
-   *
-   *  Find R', the `expected' location of R, by
-   *  reflecting S in M (the midpoint of QP).
-   *
-   *  Let 2d = |RR'|
-   *       b = |PQ|
-   *       l = |RS|
-   *
-   *  Giving energy contribution:
-   *
-   *                               2
-   *                            b d
-   *    E             =  F   .  ----
-   *     vd, edge PQ      vd      3
-   *                             l
-   *
-   *  (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'
-   *                        ,'
-   *                      ,'
-   *           R'               ,'             2d" = |SS'| = |RR'| = 2d
-   *     `-._         ,'
-   *                 `-._   ,'                 By congruent triangles,
-   *              ` M                  with M' = midpoint of RS,
-   *             ,'  `-._              |MM'| = |RR'|/2 = d
-   *           ,'        `-._
-   *         ,'              ` S       So use
-   *       ,'       M' _ , - '            d = |MM'|
-   *     ,'   _ , - '
-   *    R - '
-   *
-   *  We choose this value for l (rather than |RM|+|MS|, say, or |RM|)
-   *  because we want this symmetry and because we're happy to punish
-   *  bending more than uneveness in the metric.
-   *
-   *  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.
-   */
-
-double hypotD(const double p[D3], const double q[D3]) {
-  int k;
-  double pq[D3];
-  gsl_vector v;
-
-  K pq[_k]= p[k] - q[k];
-  v.size= D3;
-  v.stride= 1;
-  v.vector.data= pq;
-  /* owner and block ought not to be used */
-
-  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;
-}
-
-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 fma_fast fma
-#else
-# 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;
-
-static double energy_function(const double vertices[N][D3]) {
-  int pi,e,qi,ri,si, k;
-  double m[D3], mprime[D3], b, d2, l, sigma_bd2_l;
-
-  FOR_EDGE(pi,e,qi) {
-    ri= EDGE_END2(pi,(e+1)%V6); if (r<0) continue;
-    si= EDGE_END2(pi,(e+5)%V6); if (s<0) continue;
-    assert(ri == EDGE_END2(qi,(e+2)%V6));
-    assert(si == EDGE_END2(qi,(e+4)%V6));
-    
-    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= hypotD2(m, mprime);
-    l= hypotD(vertices[ri][k] - vertices[si][k]);
-    l3 = l*l*l + l3_epsilon;
-
-    sigma_bd2_l += b * d2 / l3;
-  }
-
-  
-
-static gsl_multimin_fminimizer *minimiser;
-
-
-
-{
-  gsl_multimin_fminimizer_alloc(
-
diff --git a/bgl.cpp b/bgl.cpp
index eb7e096..fbc5b7d 100644 (file)
--- a/bgl.cpp
+++ b/bgl.cpp
@@ -1,5 +1,14 @@
+/*
+ * Everything that needs the Boost Graph Library and C++ templates etc.
+ * (and what a crazy set of stuff that all is)
+ */
+
+#include <math.h>
+
 extern "C" {
+#include "bgl.h"
 #include "mgraph.h"
+#include "common.h"
 }
 
 /*
@@ -27,9 +36,11 @@ extern "C" {
 #define VMASK (YMASK|XMASK)
 #define ESHIFT (YBITS|XBITS)
 
+class Graph { }; // this is a dummy as our graph has no actual representation
+
 namespace boost {
-  // We make Layout a model of various BGL Graph concepts.
-  // This mainly means that graph_traits<Layout> has lots of stuff.
+  // We make Graph a model of various BGL Graph concepts.
+  // This mainly means that graph_traits<Graph> has lots of stuff.
 
   // First, some definitions used later:
   
@@ -44,7 +55,7 @@ namespace boost {
     OutEdgeIncrable(int v, int e) : f(v | (e << ESHIFT)) { }
   };
  
-  struct graph_traits<Layout> {
+  struct graph_traits<Graph> {
 
     // Concept Graph:
   
@@ -61,37 +72,38 @@ namespace boost {
       forward_iterator_tag> out_edge_iterator;
     typedef int degree_size_type;
     
-    inline int source(int f, const Layout&) { return f&VMASK; }
-    inline int target(int f, const Layout&) { return EDGE_END2(f&VMASK, f>>ESHIFT); }
+    inline int source(int f, const Graph&) { return f&VMASK; }
+    inline int target(int f, const Graph&) { return EDGE_END2(f&VMASK, f>>ESHIFT); }
     inline std::pair<out_edge_iterator,out_edge_iterator>
-    out_edges(int v, const Layout&) {
+    out_edges(int v, const Graph&) {
       return std::make_pair(out_edge_iterator(OutEdgeIncrable(v, VE_MIN(v))),
                            out_edge_iterator(OutEdgeIncrable(v, VE_MAX(v))));
     }
-    out_degree(int v, const Layout&) { return VE_MAX(v) - VE_MIN(v); }
+    inline out_degree(int v, const Graph&) { return VE_MAX(v) - VE_MIN(v); }
 
     // Concept VertexListGraph:
     typedef counting_iterator<int> vertex_iterator;
     typedef unsigned vertices_size_type;
     inline std::pair<vertex_iterator,vertex_iterator>
-    vertices(const Layout&) {
+    vertices(const Graph&) {
       return std::make_pair(vertex_iterator(0), vertex_iterator(N));
     }
-    inline unsigned num_vertices(const Layout&) { return N; }
-
-}
+    inline unsigned num_vertices(const Graph&) { return N; }
+  };
+};
 
 static void single_source_shortest_paths(int v1,
                                         const double edge_weights[/*f*/],
                                         double vertex_distances[/*v*/]) {
-  boost::dijkstra_shortest_paths
-    (g, v1,
+  Graph g;
+
+  boost::dijkstra_shortest_paths(g, v1,
      weight_map(edge_weights).
      vertex_index_map(identity_property_map()).
      distance_map(vertex_distances));
 }
     
-double graph_layout_cost(const Layout *g, const double vertex_areas[N]) {
+double graph_layout_cost(const Vertices v, 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|.
    *
@@ -108,23 +120,25 @@ double graph_layout_cost(const Layout *g, const double vertex_areas[N]) {
    * (In practice we compute d^2+epsilon and use it for the
    *  divisions, to avoid division by zero.)
    */
+  static const d2_epsilon= 1e-6;
+  
   double edge_weights[N*V6], vertex_distances[N], total_cost;
-  int v1, e, f;
+  int v1,v2,e,f;
 
-  FOR_VEDGE_X(v1,e,
+  FOR_VEDGE_X(v1,e,v2,
              f= v1 | e << ESHIFT,
              edge_weights[f]= NaN)
-    edge_weights[f]= hypotD(g.v[v1], g.v[v2]);
+    edge_weights[f]= hypotD(v[v1], v[v2]);
 
   FOR_VERTEX(v1) {
     double a1= vertex_areas[v1];
     single_source_shortest_paths(v1, edge_weights, vertex_distances);
     FOR_VERTEX(v2) {
       double a2= vertex_areas[v2];
-      double d2= hypotD2plus(g->v[v1],g->v[v2], d2_epsilon);
+      double d2= hypotD2plus(v[v1],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);
+      total_cost += a1*a2 * (sd2 - 1) / (d2*d2);
     }
   }
   return total_cost;
diff --git a/bgl.h b/bgl.h
new file mode 100644 (file)
index 0000000..68700b9
--- /dev/null
+++ b/bgl.h
@@ -0,0 +1,10 @@
+/*
+ * Definitions exported and imported by our BGL-using stuff
+ */
+
+#ifndef BGL_H
+#define BGL_H
+
+double graph_layout_cost(const Vertices v, const double vertex_areas[N]);
+
+#endif /*BGL_H*/
diff --git a/common.c b/common.c
new file mode 100644 (file)
index 0000000..146e5d5
--- /dev/null
+++ b/common.c
@@ -0,0 +1,41 @@
+/*
+ * Generally useful stuff.
+ */
+
+#include "common.h"
+
+double hypotD1(const double pq[D3]) {
+  gsl_vector v;
+
+  v.size= D3;
+  v.stride= 1;
+  v.vector.data= pq;
+  /* owner and block ought not to be used */
+
+  return gsl_blas_snrm2(&v);
+}
+
+double hypotD(const double p[D3], const double q[D3]) {
+  int k;
+  double pq[D3];
+
+  K pq[k]= p[k] - q[k];
+  return hypotD1(pq);
+}
+
+double hypotD2(const double p[D3], const double q[D3]) {
+  double d2= 0;
+  K d2= ffsqa(p[k] - q[k], d2);
+  return d2;
+}
+
+double hypotD2plus(const double p[D3], const double q[D3], double d2) {
+  K d2= ffsqa(p[k] - q[k], d2);
+  return d2;
+}
+
+void xprod(double r[D3], const double a[D3], const double b[D3]) {
+  r[0]= a[1]*b[2] - a[2]*b[1]);
+  r[1]= a[2]*b[0] - a[0]*b[2]);
+  r[2]= a[0]*b[1] - a[1]*b[0]);
+}
diff --git a/common.h b/common.h
new file mode 100644 (file)
index 0000000..1a358d1
--- /dev/null
+++ b/common.h
@@ -0,0 +1,23 @@
+/*
+ * Generally useful stuff.
+ */
+
+#ifndef COMMON_H
+#define COMMON_H
+
+#define _GNU_SOURCE
+#include <math.h>
+#include <limits.h>
+
+double hypotD(const double p[D3], const double q[D3]);
+double hypotD2(const double p[D3], const double q[D3]);
+double hypotD2plus(const double p[D3], const double q[D3], double add);
+
+#ifdef FP_FAST_FMA
+# define fma_fast fma
+#else
+# define fma_fast(f1,f2,t) ((f1)*(f2)+(t))
+#endif
+#define ffsqa(factor,term) fma_fast((factor),(factor),(term))
+
+#endif /*COMMON_H*/
diff --git a/energy.c b/energy.c
new file mode 100644 (file)
index 0000000..c3b8754
--- /dev/null
+++ b/energy.c
@@ -0,0 +1,263 @@
+/*
+ * We try to find an optimal triangle grid
+ */
+
+#include "common.h"
+#include "bgl.h"
+#include "mgraph.h"
+
+#define BEST_F "best"
+#define INITIAL_F "initial"
+
+static double edgewise_vertex_displacement_cost(const Vertices vertices);
+
+static void compute_vertex_areas(const Vertices vertices, double areas[N]);
+static double best_energy= DOUBLE_MAX;
+static void flushoutput(void);
+
+static void cost(double *energy, double tweight, double tcost);
+#define COST(weight, compute) cost(&energy, (weight), (compute))
+
+/*---------- main energy computation and subroutines ----------*/
+
+static double compute_energy(Vertices vertices) {
+  double vertex_areas[N], energy;
+
+  compute_vertex_areas(vertices,vertex_areas);
+  energy= 0;
+  printf("cost > energy |");
+
+  COST(1000.0, edgewise_vertex_displacement_cost(vertices));
+  COST(1.0,    graph_layout_cost(vertices,vertex_areas));
+  COST(1e6,    noncircular_edge_cost(vertices));
+  
+  printf("| total %# e |", energy);
+  if (energy < best_energy) {
+    FILE *best;
+    printf(" BEST");
+    
+    best_f= fopen(BEST_F ".new","wb");  if (!best_f) diee("fopen new best");
+    r= fwrite(vertices,sizeof(vertices),1,best_f); if (r!=1) diee("fwrite");
+    if (fclose(best_f)) diee("fclose new best");
+    if (rename(BEST_F ".new", BEST_F)) diee("rename install new best");
+  }
+  putchar('\n');
+  flushoutput();
+
+  return energy;
+}    
+
+static void cost(double *energy, double tweight, double tcost) {
+  double tenergy= tweight * tcost;
+  printf(" %# e > %# e |", tcost, tenergy);
+  *energy += tenergy;
+}
+
+static void flushoutput(void) {
+  if (fflush(stdout) || ferror(stdout)) { perror("stdout"); exit(-1); }
+}
+
+static void compute_vertex_areas(const Vertices vertices, double areas[N]) {
+  FOR_VERTEX(v0) {
+    double total= 0.0;
+    int count= 0;
+    
+    FOR_VEDGE(v0,e1,v1) {
+      e2= (e1+1) % V6;
+      v2= EDGE_END2(v0,e2);
+      if (v2<0) continue;
+      
+      double e1v[D3], e2v[D3], av[D3];
+      K {
+       e1v[k]= vertices[v1][k] - vertices[v0][k];
+       e2v[k]= vertices[v2][k] - vertices[v0][k];
+      }
+      xprod(av, e1v, e2v);
+      total += hypotD1(av);
+      count++;
+    }
+    areas[v0]= total / count;
+  }
+}
+
+/*---------- use of GSL ----------*/
+
+  /* We want to do multidimensional minimisation.
+   *
+   * We don't think there are any local minima.  Or at least, if there
+   * are, the local minimum which will be found from the starting
+   * state is the one we want.
+   *
+   * We don't want to try to provide a derivative of the cost
+   * function.  That's too tedious (and anyway the polynomial
+   * approximation to our our cost function sometimes has high degree
+   * in the inputs which means the quadratic model implied by most of
+   * the gradient descent minimisers is not ideal).
+   *
+   * This eliminates most of the algorithms.  Nelder and Mead's
+   * simplex algorithm is still available and we will try that.
+   *
+   * In our application we are searching for the optimal locations of
+   * N actualvertices in D3 (3) dimensions - ie, we are searching for
+   * the optimal metapoint in an N*D3-dimensional space.
+   * 
+   * So eg with X=Y=100, the simplex will contain 300 metavertices
+   * each of which is an array of 300 doubles for the actualvertex
+   * coordinates.  Hopefully this won't be too slow ...
+   */
+
+static void gsldie(const char *what, int status) {
+  fprintf(stderr,"gsl function failed: %s: %s\n", what, gsl_strerror(status));
+  exit(-1);
+}
+
+static gsl_multimin_fminimizer *minimiser;
+
+static const stop_epsilon= 1e-4;
+
+#define DIM (N*D3)
+
+static double minfunc_f(const gsl_vector *x, void *params) {
+  assert(x->size == DIM);
+  assert(x->stride == 1);
+  return compute_energy((Vertices)x->data);
+}
+
+int main(int argc, const char *const *argv) {
+  struct gsl_multimin_function multimin_function;
+  double size;
+  Vertices initial;
+  FILE *initial;
+  gsl_vector initial_gsl, *step_size;
+  int r;
+  
+  if (argc>1) { fputs("takes no arguments\n",stderr); exit(8); }
+
+  minimiser= gsl_multimin_fminimizer_alloc
+    (gsl_multimin_fminimizer_nmsimplex, DIM);
+  if (!minimiser) { perror("alloc minimiser"); exit(-1); }
+
+  multimin_function.f= minfunc_f;
+  multimin_function.n= DIM;
+  multimin_function.params= 0;
+
+  initial_f= fopen(INITIAL_F,"rb");  if (!initial_f) diee("fopen initial");
+  errno= 0; r= fread(initial,sizeof(initial),1,initial_f);
+  if (r!=1) diee("fread");
+  fclose(initial_f);
+
+  initial_gsl.size= DIM;
+  initial_gsl.stride= 1;
+  initial_gsl.data= initial;
+  initial_gsl.block= 0;
+  initial_gsl.owner= 0;
+
+  step_size= gsl_vector_alloc(DIM);  if (!step_size) gsldie("alloc step");
+  gsl_vector_set_all(step_size, 1e-3);
+
+  r= gsl_multimin_fminimizer_set(minimiser, &multimin_function,
+                                &initial_gsl, &step_size);
+  if (r) { gsldie("fminimizer_set",r); }
+  
+  for (;;) {
+    r= gsl_multimin_fminimizer_iterate(minimiser);
+    if (r) { gsldie("fminimizer_iterate",r); }
+
+    size= gsl_multimin_fminimizer_size(minimiser);
+    r= gsl_multimin_test_size(size, stop_epsilon);
+
+    printf("size %# e, r=%d\n", size, r);
+    flushoutput();
+
+    if (r==GSL_SUCCESS) break;
+    assert(r==GSL_CONTINUE);
+  }
+}
+
+/*---------- Edgewise vertex displacement ----------*/
+
+  /*
+   *  
+   *
+   *
+   *                       Q `-_
+   *              / |    `-_
+   *   R' - _ _ _/_ |       `-.
+   *    .       /   M - - - - - S
+   *    .      /    |      _,-'
+   *    .     /     |  _,-'
+   *    .    /    , P '
+   *    .   /  ,-'
+   *    .  /,-'
+   *    . /'
+   *            R
+   *
+   *
+   *
+   *  Find R', the `expected' location of R, by
+   *  reflecting S in M (the midpoint of QP).
+   *
+   *  Let 2d = |RR'|
+   *       b = |PQ|
+   *       l = |RS|
+   *
+   *  Giving energy contribution:
+   *
+   *                               2
+   *                            b d
+   *    E             =  F   .  ----
+   *     vd, edge PQ      vd      3
+   *                             l
+   *
+   *  (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'
+   *                        ,'
+   *                      ,'
+   *           R'               ,'             2d" = |SS'| = |RR'| = 2d
+   *     `-._         ,'
+   *                 `-._   ,'                 By congruent triangles,
+   *              ` M                  with M' = midpoint of RS,
+   *             ,'  `-._              |MM'| = |RR'|/2 = d
+   *           ,'        `-._
+   *         ,'              ` S       So use
+   *       ,'       M' _ , - '            d = |MM'|
+   *     ,'   _ , - '
+   *    R - '
+   *
+   *  We choose this value for l (rather than |RM|+|MS|, say, or |RM|)
+   *  because we want this symmetry and because we're happy to punish
+   *  bending more than uneveness in the metric.
+   *
+   *  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.
+   */
+
+static double edgewise_vertex_displacement_cost(const Vertices vertices) {
+  static const l3_epsison= 1e-6;
+
+  int pi,e,qi,ri,si, k;
+  double m[D3], mprime[D3], b, d2, l, sigma_bd2_l3;
+
+  FOR_EDGE(pi,e,qi) {
+    ri= EDGE_END2(pi,(e+1)%V6); if (r<0) continue;
+    si= EDGE_END2(pi,(e+5)%V6); if (s<0) continue;
+    assert(ri == EDGE_END2(qi,(e+2)%V6));
+    assert(si == EDGE_END2(qi,(e+4)%V6));
+    
+    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= hypotD2(m, mprime);
+    l= hypotD(vertices[ri][k] - vertices[si][k]);
+    l3 = l*l*l + l3_epsilon;
+
+    sigma_bd2_l3 += b * d2 / l3;
+  }
+  return sigma_bd2_l3;
+}
index 3233890..222bd09 100644 (file)
--- a/mgraph.c
+++ b/mgraph.c
@@ -1,3 +1,7 @@
+/*
+ * Graph topology
+ */
+
 #include "mgraph.h"
 
 static const unsigned dx[V6]= {   0,  +1,  -1,  +1,  -1,   0  },
index ae341fa..0991a7f 100644 (file)
--- a/mgraph.h
+++ b/mgraph.h
@@ -1,4 +1,7 @@
 /*
+ * Graph topology
+ */
+/*
  * Vertices in strip are numbered as follows:
  *
  *     ___ X-2 ___ X-1 ___ 0   ___ 1   ___ 2   ___ 3   ___ 4   __
@@ -80,7 +83,7 @@
 extern int edge_end2(unsigned v1, int e);
 #define EDGE_END2 edge_end2
 
-#define FOR_VEDGE_X(v1,e,init,otherwise)       \
+#define FOR_VEDGE_X(v1,e,v2,init,otherwise)    \
   FOR_VPEDGE((v1),(e))                         \
     if (((v2)= EDGE_END2((v1),(e)),            \
         (init),                                \
@@ -88,11 +91,11 @@ extern int edge_end2(unsigned v1, int e);
 
 #define NOTHING ((void)0)
 
-#define FOR_VEDGE(v1,e) \
-  FOR_VEDGE_X(v1,e,NOTHING,NOTHING)
+#define FOR_VEDGE(v1,e,v2)                     \
+  FOR_VEDGE_X(v1,e,v2,NOTHING,NOTHING)
 
-#define FOR_EDGE(v1,e,v2) \
-  FOR_VERTEX((v1)) \
+#define FOR_EDGE(v1,e,v2)                      \
+  FOR_VERTEX((v1))                             \
     FOR_VEDGE((v1),(e),(v2))
 
 #define FOR_COORD(k) \
@@ -100,11 +103,6 @@ extern int edge_end2(unsigned v1, int e);
 
 #define K FOR_COORD(k)
 
-typedef struct {
-  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]);
+typedef double Vertices[N][D3];
 
 #endif /*MGRAPH_H*/