X-Git-Url: http://www.chiark.greenend.org.uk/ucgi/~ian/git?p=moebius2.git;a=blobdiff_plain;f=energy.c;h=308aebfa3f890aea945bafc0aa38fc1feb18f9c8;hp=3459a726351724226f4fe7994c3fdc75d5b4d909;hb=6df67fc7061e6fae41fec8c05b9e49e4e393c733;hpb=2f8c87aa65ec594687eac88d4311475e289d8e51 diff --git a/energy.c b/energy.c index 3459a72..308aebf 100644 --- a/energy.c +++ b/energy.c @@ -5,54 +5,129 @@ #include "common.h" #include "minimise.h" #include "mgraph.h" +#include "parallel.h" double vertex_areas[N], vertex_mean_edge_lengths[N], edge_lengths[N][V6]; static double best_energy= DBL_MAX; static void addcost(double *energy, double tweight, double tcost, int pr); -#define COST(weight, compute) addcost(&energy, (weight), (compute), printing) + +/*---------- main energy computation, weights, etc. ----------*/ + +typedef double CostComputation(const Vertices vertices, int section); +typedef void PreComputation(const Vertices vertices, int section); + +typedef struct { + double weight; + CostComputation *fn; +} CostContribution; + +#define NPRECOMPS ((sizeof(precomps)/sizeof(precomps[0]))) +#define NCOSTS ((sizeof(costs)/sizeof(costs[0]))) +#define COST(weight, compute) { (weight),(compute) }, + +static PreComputation *const precomps[]= { + compute_edge_lengths, + compute_vertex_areas +}; + +static const CostContribution costs[]= { + +#if XBITS==3 +#define STOP_EPSILON 1e-6 + COST( 3e3, line_bending_cost) + COST( 3e3, edge_length_variation_cost) + COST( 0.4e3, rim_proximity_cost) + COST( 1e6, edge_angle_cost) + #define EDGE_ANGLE_COST_CIRCCIRCRAT (0.5/1.7) +// COST( 1e1, small_triangles_cost) + COST( 1e12, noncircular_rim_cost) +#endif + +#if XBITS==4 +#define STOP_EPSILON 1.2e-4 + COST( 3e5, line_bending_cost) + COST( 10e3, edge_length_variation_cost) + COST( 9.0e3, rim_proximity_cost) // 5e1 is too much + // 2.5e1 is too little + // 0.2e1 grows compared to previous ? + // 0.6e0 shrinks compared to previous ? + +// COST( 1e12, edge_angle_cost) + #define EDGE_ANGLE_COST_CIRCCIRCRAT (0.5/1.3) + COST( 1e18, noncircular_rim_cost) +#endif + +#if XBITS>=5 /* nonsense follows but never mind */ +#define STOP_EPSILON 1e-6 + COST( 3e5, line_bending_cost) + COST( 10e2, edge_length_variation_cost) + COST( 9.0e1, rim_proximity_cost) // 5e1 is too much + // 2.5e1 is too little + // 0.2e1 grows compared to previous ? + // 0.6e0 shrinks compared to previous ? + + COST( 1e12, edge_angle_cost) + #define EDGE_ANGLE_COST_CIRCCIRCRAT (0.5/1.3) + COST( 1e18, noncircular_rim_cost) +#endif + +}; + +const double edge_angle_cost_circcircrat= EDGE_ANGLE_COST_CIRCCIRCRAT; void energy_init(void) { + stop_epsilon= STOP_EPSILON; +} + +/*---------- energy computation machinery ----------*/ + +void compute_energy_separately(const struct Vertices *vs, + int section, void *energies_v, void *totals_v) { + double *energies= energies_v; + int ci; + + for (ci=0; cia, section); + inparallel_barrier(); + } + for (ci=0; cia, section); } -/*---------- main energy computation and subroutines ----------*/ +void compute_energy_combine(const struct Vertices *vertices, + int section, void *energies_v, void *totals_v) { + int ci; + double *energies= energies_v; + double *totals= totals_v; + + for (ci=0; cia); - compute_vertex_areas(vs->a); - energy= 0; + double totals[NCOSTS], energy; + int ci, printing; printing= printing_check(pr_cost,0); if (printing) printf("%15lld c>e |", evaluations); - if (XBITS==3) { - COST( 3e2, line_bending_cost(vs->a)); - COST( 1e3, edge_length_variation_cost(vs->a)); - COST( 0.4e3, rim_proximity_cost(vs->a)); - COST( 1e6, edge_angle_cost(vs->a)); -// COST( 1e1, small_triangles_cost(vs->a)); - COST( 1e12, noncircular_rim_cost(vs->a)); - stop_epsilon= 1e-6; - } else if (XBITS==4) { - COST( 3e2, line_bending_cost(vs->a)); - COST( 3e3, edge_length_variation_cost(vs->a)); - COST( 9.0e1, rim_proximity_cost(vs->a)); // 5e1 is too much - // 2.5e1 is too little - // 0.2e1 grows compared to previous ? - // 0.6e0 shrinks compared to previous ? - COST( 1e12, edge_angle_cost(vs->a)); - COST( 1e12, noncircular_rim_cost(vs->a)); - stop_epsilon= 1e-5; - } else { - abort(); - } + for (ci=0; ci> YSHIFT; int nominal_edge_distance= y <= Y/2 ? y : Y-1-y; if (nominal_edge_distance==0) continue; @@ -249,12 +323,12 @@ double rim_proximity_cost(const Vertices vertices) { /*---------- noncircular rim cost ----------*/ -double noncircular_rim_cost(const Vertices vertices) { +double noncircular_rim_cost(const Vertices vertices, int section) { int vy,vx,v; double cost= 0.0; double oncircle[3]; - FOR_RIM_VERTEX(vy,vx,v) { + FOR_RIM_VERTEX(vy,vx,v, OUTER) { find_nearest_oncircle(oncircle, vertices[v]); double d2= hypotD2(vertices[v], oncircle); @@ -290,16 +364,16 @@ double noncircular_rim_cost(const Vertices vertices) { * vd, edge PQ vd */ -double edge_angle_cost(const Vertices vertices) { +double edge_angle_cost(const Vertices vertices, int section) { double pq1[D3], rp[D3], ps[D3], rp_2d[D3], ps_2d[D3], rs_2d[D3]; double a,b,c,s,r; - const double minradius= 0.5; + const double minradius_base= 0.2; int pi,e,qi,ri,si, k; // double our_epsilon=1e-6; double total_cost= 0; - FOR_EDGE(pi,e,qi) { + FOR_EDGE(pi,e,qi, OUTER) { // if (!(RIM_VERTEX_P(pi) || RIM_VERTEX_P(qi))) continue; si= EDGE_END2(pi,(e+V6-1)%V6); if (si<0) continue; @@ -322,8 +396,8 @@ double edge_angle_cost(const Vertices vertices) { c= magnD(rs_2d); s= 0.5*(a+b+c); r= a*b*c / sqrt((a+b+c)*(a-b+c)*(b-c+a)*(c-a+b) + 1e-6); -//fprintf(stderr,"triangle a=%g b=%g c=%g -> s=%g r=%g\n",a,b,c,s,r); + double minradius= minradius_base + edge_angle_cost_circcircrat*(a+b); double deficit= minradius - r; if (deficit < 0) continue; double cost= deficit*deficit; @@ -359,14 +433,14 @@ double edge_angle_cost(const Vertices vertices) { * vd, edge PQ vd */ -double small_triangles_cost(const Vertices vertices) { +double small_triangles_cost(const Vertices vertices, int section) { double pq[D3], ps[D3]; double x[D3]; int pi,e,qi,si, k; // double our_epsilon=1e-6; double total_cost= 0; - FOR_EDGE(pi,e,qi) { + FOR_EDGE(pi,e,qi, OUTER) { // if (!(RIM_VERTEX_P(pi) || RIM_VERTEX_P(qi))) continue; si= EDGE_END2(pi,(e+V6-1)%V6); if (si<0) continue;