From 732b811081946ad56d05769de8846b27375b7eb7 Mon Sep 17 00:00:00 2001 From: Ian Jackson Date: Sat, 29 Dec 2007 20:42:39 +0000 Subject: [PATCH] reorganisation is complete, need to implement various things, make it compile, etc. --- anneal.c | 131 --------------------------- bgl.cpp | 52 +++++++---- bgl.h | 10 +++ common.c | 41 +++++++++ common.h | 23 +++++ energy.c | 263 +++++++++++++++++++++++++++++++++++++++++++++++++++++++ mgraph.c | 4 + mgraph.h | 20 ++--- 8 files changed, 383 insertions(+), 161 deletions(-) delete mode 100644 anneal.c create mode 100644 bgl.h create mode 100644 common.c create mode 100644 common.h create mode 100644 energy.c diff --git a/anneal.c b/anneal.c deleted file mode 100644 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 --- 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 + 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 has lots of stuff. + // We make Graph a model of various BGL Graph concepts. + // This mainly means that graph_traits 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 { + struct graph_traits { // 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_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 vertex_iterator; typedef unsigned vertices_size_type; inline std::pair - 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 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 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 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 +#include + +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 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; +} diff --git a/mgraph.c b/mgraph.c index 3233890..222bd09 100644 --- 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 }, diff --git a/mgraph.h b/mgraph.h index ae341fa..0991a7f 100644 --- a/mgraph.h +++ b/mgraph.h @@ -1,3 +1,6 @@ +/* + * Graph topology + */ /* * Vertices in strip are numbered as follows: * @@ -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*/ -- 2.30.2