From bc3b5b8813e120b28b3148d8c7412fa9937f6cb1 Mon Sep 17 00:00:00 2001 From: Ian Jackson Date: Sun, 6 Jan 2008 18:11:21 +0000 Subject: [PATCH] exploit symmetry --- Makefile | 2 +- energy.c | 179 ++--------------------------------------------------- half.c | 157 ++++++++++++++++++++++++++++++++++++++++++++++ half.h | 38 ++++++++++++ mgraph.h | 56 +++++++++-------- minimise.c | 175 +++++++++++++++++++++++++++++++++++++++++++++++++++ minimise.h | 8 +++ primer.c | 6 +- 8 files changed, 418 insertions(+), 203 deletions(-) create mode 100644 half.c create mode 100644 half.h create mode 100644 minimise.c diff --git a/Makefile b/Makefile index e1a54dd..15e4ae0 100644 --- a/Makefile +++ b/Makefile @@ -14,7 +14,7 @@ o= >$@.new && mv -f $@.new $@ all: $(TARGETS) -minimise: energy.o graph.o common.o mgraph.o +minimise: energy.o graph.o common.o mgraph.o minimise.o half.o $(CXX) $(CXXFLAGS) -o $@ $^ $(LIBGSL) primer: primer.o common.o diff --git a/energy.c b/energy.c index 78b8f31..ed50bc7 100644 --- a/energy.c +++ b/energy.c @@ -6,41 +6,28 @@ #include "minimise.h" #include "mgraph.h" -#include -#include - -#include -#include - -static const char *input_file, *output_file; -static char *output_file_tmp; - static void compute_vertex_areas(const Vertices vertices, double areas[N]); static double best_energy= DBL_MAX; -enum printing_instance { pr_cost, pr_size, pr__max }; - static void addcost(double *energy, double tweight, double tcost, int pr); #define COST(weight, compute) addcost(&energy, (weight), (compute), printing) -static int printing_check(enum printing_instance); -static void printing_init(void); /*---------- main energy computation and subroutines ----------*/ -static double compute_energy(const Vertices vertices) { +double compute_energy(const struct Vertices *vs) { double vertex_areas[N], energy; int printing; - compute_vertex_areas(vertices,vertex_areas); + compute_vertex_areas(vs->a, vertex_areas); energy= 0; printing= printing_check(pr_cost); if (printing) printf("cost > energy |"); - COST(1e2, edgewise_vertex_displacement_cost(vertices)); - COST(1e2, graph_layout_cost(vertices,vertex_areas)); -// COST(1e4, noncircular_rim_cost(vertices)); + COST(1e2, edgewise_vertex_displacement_cost(vs->a)); + COST(1e2, graph_layout_cost(vs->a,vertex_areas)); +// COST(1e4, noncircular_rim_cost(vs->a)); if (printing) printf("| total %# e |", energy); @@ -51,7 +38,7 @@ static double compute_energy(const Vertices vertices) { if (printing) printf(" BEST"); best_f= fopen(output_file_tmp,"wb"); if (!best_f) diee("fopen new out"); - r= fwrite(vertices,sizeof(Vertices),1,best_f); if (r!=1) diee("fwrite"); + r= fwrite(vs->a,sizeof(vs->a),1,best_f); if (r!=1) diee("fwrite"); if (fclose(best_f)) diee("fclose new best"); if (rename(output_file_tmp,output_file)) diee("rename install new best"); @@ -96,107 +83,6 @@ static void compute_vertex_areas(const Vertices vertices, double areas[N]) { } } -/*---------- 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 gsl_multimin_fminimizer *minimiser; - -static const double stop_epsilon= 1e-6; - -static double minfunc_f(const gsl_vector *x, void *params) { - assert(x->size == DIM); - assert(x->stride == 1); - return compute_energy((const double(*)[D3])x->data); -} - -int main(int argc, const char *const *argv) { - gsl_multimin_function multimin_function; - double size; - Vertices initial, step_size; - FILE *initial_f; - gsl_vector initial_gsl, step_size_gsl; - int r, v, k; - - if (argc!=3 || argv[1][0]=='-' || strncmp(argv[2],"-o",2)) - { fputs("usage: minimise -o= 0); assert(v < N); +#define CHECK_NEAR(a,b) \ + if (abs(a-b) >= 1e-6) { printf(" BAD!"); bad++; } \ + putchar('\n'); + + FOR_MAP_DIMENSION( + (printf("i=%2d v=%03x k=%d %# 20g\n", + i,v,k, + vs->a[v][k] + )); + CHECK_V(v) + ixfrom[v][k]=i; + count[v][k]++;, + +// FOR_MAP_DIMENSION_VPRIME; + (printf("i=%2d v'=%03x k=%d %# 20g := %# 20g", + i,vprime,k, + vs->a[vprime][k], sign*vs->a[v][k] + )); + CHECK_V(vprime) + CHECK_NEAR(vs->a[vprime][k], sign*vs->a[v][k]) + ixfrom[vprime][k]=-i; + count[vprime][k]++;, + + (printf("zero v=%03x k=%d %# 20g := 0", + v,k, + vs->a[v][k] + )); + CHECK_V(v) + CHECK_NEAR(vs->a[v][k], 0) + count[v][k]++; + ); + + FOR_VERTEX(vc) + FOR_COORD(kc) { + if (count[vc][kc]==1) continue; + printf("BAD %03x %d count=%d ixfrom=%d\n", + vc,kc, count[vc][kc], ixfrom[vc][kc]); + bad++; + } + assert(!bad); +} + +void map_dimensions(const double metavertex[DIM], Vertices out) { + FOR_MAP_DIMENSION(out[v ][k]= metavertex[i];, + +// FOR_MAP_DIMENSION_VPRIME; + out[vprime][k]= sign*metavertex[i];, + + out[v ][k]= 0; + ); +} + +void unmap_dimensions(double metavertex[DIM], const struct Vertices *in) { + FOR_MAP_DIMENSION( + metavertex[i]= in->a[v][k];, + , + ); +} + + +/* + * OLD NUMBERING + * + * Our numbering is as follows: + * + * :axis of symmetry + * | : + * | : + * ___| 0 ___ 1 ___ 2 ___ 3 ___ 4 ___ up to X"-1 + * | : where X"=(X/2) + * / :\ / \ / \ / \ / \ + * /| : \ / \ / \ / \ / \ + * |___ 0 ___ 1 ___ 2 ___ 3 ___ 4 + * | : +X" +X" +X" +X" +X" + * \| : / \ / \ / \ / \ / + * \ :/ \ / \ / \ / \ / + * ___|2*X"___ 1+ ___ 2+ ___ 3+ __ 4+ ___ + * | : 2*X" 2*X" 2*X" 2*X" + * / :\ / \ / \ / \ / \ + * /| : \ / \ / \ / \ / \ + * |___ 0+ ___ 1+ ___ 2+ ___ 3+ ___ 4+ + * | : 3*X" 3*X" 3*X" 3*X" 3*X" + * | : + * + * X" is still a power of two so a vertex pair can be named by a + * number structured like this: + * + * vertex number": 0000 | y | x + * YBITS XBITS-1 + * + */ diff --git a/half.h b/half.h new file mode 100644 index 0000000..88c2e1a --- /dev/null +++ b/half.h @@ -0,0 +1,38 @@ +/* + * During the actual optimisation we constrain it to be symmetrical. + * That means we have half as many vertices. We process the vertices + * with X >=0. + * + * We conventionally suffix all of our variables, macros, constants, + * etc., with `h', to represent the double-prime ". + * + * Vertices with x==0 are special as they lie on the symmetry axis: + * all but one of the coordinates are fixed. We do not offer those + * coordinates to our optimisation algorithm. + * + * We provide to the rest of the program just: + * NH number of vertices" + * VH2V(v) maps a v" to the representative vertex number v + * FOR_VERTEX_H(vh,v) iterates over v", computing each v + * DIM dimensionality of coordinate space + * [un]map_dimensions maps double[DIM] to Vertices and back + */ + +#ifndef HALF_H +#define HALF_H + +#include "mgraph.h" + +#define NH (N/2) +#define VH2V(v) ((v<<1) & YMASK | v & (XMASK>>1)) +#define FOR_VERTEX_H(vh,v) \ + for ((vh)=0; (v)= VH2V((vh)), (vh) < NH; (vh)++) + +#define HALF_DEGRADED_NVERTICES ((Y+1)/2) +#define DIM (NH*D3 + HALF_DEGRADED_NVERTICES*(1-D3) + ((Y+1)/2)*D3) + +void unmap_dimensions(double metavertex[DIM], const struct Vertices *in); +void map_dimensions(const double metavertex[DIM], Vertices out); +void pmap_dimensions(const struct Vertices *vs); + +#endif /*HALF_H*/ diff --git a/mgraph.h b/mgraph.h index eb4bd12..8636b91 100644 --- a/mgraph.h +++ b/mgraph.h @@ -4,34 +4,38 @@ /* * Vertices in strip are numbered as follows: * - * | + * :axis of symmetry + * | : + * | : * ___ X-2 ___ X-1 ___| 0 ___ 1 ___ 2 ___ 3 ___ 4 __ - * Y-1 Y-1 |0 0 0 0 0 - * / \ / \ / \ / \ / \ / \ / \ - * / \ / \ /| \ / \ / \ / \ / \ + * Y-1 Y-1 |0: 0 0 0 0 + * / \ / \ / :\ / \ / \ / \ / !! \ + * / \ / \ /| : \ / \ / \ / \ / \ * X-3 ___ X-2 ___ X-1|___ 0 ___ 1 ___ 2 ___ 3 ___ 4 - * Y-2 Y-2 Y-2| 1 1 1 1 1 - * \ / \ / \| / \ / \ / \ / \ / - * \ / \ / \ / \ / \ / \ / \ / + * Y-2 Y-2 Y-2| : 1 1 1 1 1 + * \ / \ / \| : / \ / \ / \ / \ / + * \ / \ / \ :/ \ / \ / \ / \ / * ___ X-2 ___ X-1 ___| 0 ___ 1 ___ 2 ___ 3 __ 4 ___ - * Y-3 Y-3 |2 2 2 2 2 - * / \ / \ / \ / \ / \ / \ / \ - * / \ / \ /| \ / \ / \ / \ / \ + * Y-3 Y-3 |2: 2 2 2 2 + * / \ / \ / :\ / \ / \ / \ / \ + * / \ / \ /| : \ / \ / \ / \ / \ * X-3 ___ X-2 ___ X-1|___ 0 ___ 1 ___ 2 ___ 3 ___ 4 - * Y-4 Y-4 Y-4| 3 3 3 3 3 - * | - * . . . . .| . . . . . . . . . . - * | + * Y-4 Y-4 Y-4| : 3 3 3 3 3 + * | : + * . . . . .| :. . . . . . . . . . + * | : * ___ X-2 ___ X-1 ___| 0 ___ 1 ___ 2 ___ 3 ___ 4 ___ * 2 2 |Y-3 Y-3 Y-3 Y-3 Y-3 - * / \ / \ / \ / \ / \ / \ / \ - * / \ / \ /| \ / \ / \ / \ / \ + * / \ / \ / :\ / \ / \ / \ / \ + * / \ / \ /| : \ / \ / \ / \ / \ * __ X-2 ___ X-1|___ 0 ___ 1 ___ 2 ___ 3 ___ 3 ___ 4 - * 1 1 | Y-2 Y-2 Y-2 Y-2 Y-2 Y-2 - * / \ / \ / \| / \ / \ / \ / \ / - * / \ / \ / \ / \ / \ / \ / \ / + * 1 1 | : Y-2 Y-2 Y-2 Y-2 Y-2 Y-2 + * / \ / \ / \| : / \ / \ / \ / \ / + * / \ / \ / \ :/ \ / \ / \ / \ / * -3 ___ X-2 ___ X-1 ___| 0 ___ 1 ___ 2 ___ 3 ___ 4 ___ * 0 0 0 |Y-1 Y-1 Y-1 Y-1 Y-1 + * | : + * | : * | * ^ join, where there is * a discontinuity in numbering @@ -55,6 +59,9 @@ * 3 1 0 * / \ * 4/ 5\ + * + * vertex number: 0000 | y | x + * YBITS XBITS */ #ifndef MGRAPH_H @@ -62,23 +69,17 @@ #include "common.h" -#define XBITS 3 +#define XBITS 4 #define X (1< +#include + +#include +#include + +#include "half.h" +#include "minimise.h" + +const char *input_file, *output_file; +char *output_file_tmp; + +static void printing_init(void); + +static gsl_multimin_fminimizer *minimiser; + +static const double stop_epsilon= 1e-6; + +static double minfunc_f(const gsl_vector *x, void *params) { + struct Vertices vs; + + assert(x->size == DIM); + assert(x->stride == 1); + map_dimensions(x->data, vs.a); + return compute_energy(&vs); +} + +int main(int argc, const char *const *argv) { + gsl_multimin_function multimin_function; + double size; + struct Vertices initial_full; + double initial_half[DIM], step_size[DIM]; + FILE *initial_f; + gsl_vector initial_gsl, step_size_gsl; + int r, i; + + if (argc!=3 || argv[1][0]=='-' || strncmp(argv[2],"-o",2)) + { fputs("usage: minimise -o