X-Git-Url: http://www.chiark.greenend.org.uk/ucgi/~ian/git?p=moebius2.git;a=blobdiff_plain;f=energy.c;h=b2ee07339b2401b610e728065b423421787f26fb;hp=c3b8754aff63dfa661cfd5dcbbbc41e67ae36885;hb=cd8d998ff14921022cecbfc2475365783f5a8090;hpb=732b811081946ad56d05769de8846b27375b7eb7 diff --git a/energy.c b/energy.c index c3b8754..b2ee073 100644 --- a/energy.c +++ b/energy.c @@ -6,21 +6,24 @@ #include "bgl.h" #include "mgraph.h" +#include +#include + #define BEST_F "best" #define INITIAL_F "initial" static double edgewise_vertex_displacement_cost(const Vertices vertices); +static double noncircular_rim_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 double best_energy= DBL_MAX; 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) { +static double compute_energy(const Vertices vertices) { double vertex_areas[N], energy; compute_vertex_areas(vertices,vertex_areas); @@ -29,11 +32,13 @@ static double compute_energy(Vertices vertices) { COST(1000.0, edgewise_vertex_displacement_cost(vertices)); COST(1.0, graph_layout_cost(vertices,vertex_areas)); - COST(1e6, noncircular_edge_cost(vertices)); + COST(1e3, noncircular_rim_cost(vertices)); printf("| total %# e |", energy); if (energy < best_energy) { - FILE *best; + FILE *best_f; + int r; + printf(" BEST"); best_f= fopen(BEST_F ".new","wb"); if (!best_f) diee("fopen new best"); @@ -53,11 +58,9 @@ static void cost(double *energy, double tweight, double tcost) { *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]) { + int v0,v1,v2, e1,e2, k; + FOR_VERTEX(v0) { double total= 0.0; int count= 0; @@ -73,7 +76,7 @@ static void compute_vertex_areas(const Vertices vertices, double areas[N]) { e2v[k]= vertices[v2][k] - vertices[v0][k]; } xprod(av, e1v, e2v); - total += hypotD1(av); + total += magnD(av); count++; } areas[v0]= total / count; @@ -106,30 +109,23 @@ static void compute_vertex_areas(const Vertices vertices, double areas[N]) { * 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 const double stop_epsilon= 1e-4; static double minfunc_f(const gsl_vector *x, void *params) { assert(x->size == DIM); assert(x->stride == 1); - return compute_energy((Vertices)x->data); + return compute_energy((const double(*)[D3])x->data); } int main(int argc, const char *const *argv) { - struct gsl_multimin_function multimin_function; + gsl_multimin_function multimin_function; double size; - Vertices initial; - FILE *initial; - gsl_vector initial_gsl, *step_size; - int r; + Vertices initial, step_size; + FILE *initial_f; + gsl_vector initial_gsl, step_size_gsl; + int r, v, vx,vy, k; if (argc>1) { fputs("takes no arguments\n",stderr); exit(8); } @@ -148,20 +144,23 @@ int main(int argc, const char *const *argv) { initial_gsl.size= DIM; initial_gsl.stride= 1; - initial_gsl.data= initial; initial_gsl.block= 0; initial_gsl.owner= 0; + step_size_gsl= initial_gsl; + + initial_gsl.data= (double*)initial; + step_size_gsl.data= (double*)step_size; - step_size= gsl_vector_alloc(DIM); if (!step_size) gsldie("alloc step"); - gsl_vector_set_all(step_size, 1e-3); + FOR_VERTEX(v) + K step_size[v][k]= 1e-3; + FOR_RIM_VERTEX(vx,vy,v) + step_size[v][3] *= 0.1; - r= gsl_multimin_fminimizer_set(minimiser, &multimin_function, - &initial_gsl, &step_size); - if (r) { gsldie("fminimizer_set",r); } + GA( gsl_multimin_fminimizer_set(minimiser, &multimin_function, + &initial_gsl, &step_size_gsl) ); for (;;) { - r= gsl_multimin_fminimizer_iterate(minimiser); - if (r) { gsldie("fminimizer_iterate",r); } + GA( gsl_multimin_fminimizer_iterate(minimiser) ); size= gsl_multimin_fminimizer_size(minimiser); r= gsl_multimin_test_size(size, stop_epsilon); @@ -172,6 +171,7 @@ int main(int argc, const char *const *argv) { if (r==GSL_SUCCESS) break; assert(r==GSL_CONTINUE); } + return 0; } /*---------- Edgewise vertex displacement ----------*/ @@ -239,14 +239,14 @@ int main(int argc, const char *const *argv) { */ static double edgewise_vertex_displacement_cost(const Vertices vertices) { - static const l3_epsison= 1e-6; + static const double l3_epsilon= 1e-6; int pi,e,qi,ri,si, k; - double m[D3], mprime[D3], b, d2, l, sigma_bd2_l3; + double m[D3], mprime[D3], b, d2, l, sigma_bd2_l3=0; 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; + ri= EDGE_END2(pi,(e+1)%V6); if (ri<0) continue; + si= EDGE_END2(pi,(e+5)%V6); if (si<0) continue; assert(ri == EDGE_END2(qi,(e+2)%V6)); assert(si == EDGE_END2(qi,(e+4)%V6)); @@ -254,10 +254,32 @@ static double edgewise_vertex_displacement_cost(const Vertices vertices) { 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; + l= hypotD(vertices[ri], vertices[si]); + double l3 = l*l*l + l3_epsilon; sigma_bd2_l3 += b * d2 / l3; } return sigma_bd2_l3; } + +/*---------- noncircular rim cost ----------*/ + +static double noncircular_rim_cost(const Vertices vertices) { + int vy,vx,v; + double cost= 0.0; + + FOR_RIM_VERTEX(vy,vx,v) { + double oncircle[3]; + /* By symmetry, nearest point on circle is the one with + * the same angle subtended at the z axis. */ + oncircle[0]= vertices[v][0]; + oncircle[1]= vertices[v][1]; + oncircle[2]= 0; + double mult= 1.0/ magnD(oncircle); + oncircle[0] *= mult; + oncircle[1] *= mult; + double d2= hypotD2(vertices[v], oncircle); + cost += d2*d2; + } + return cost; +}