X-Git-Url: http://www.chiark.greenend.org.uk/ucgi/~ian/git?p=moebius2.git;a=blobdiff_plain;f=energy.c;h=dd9206f4b7a5cb80c3178cfc9f35ae603982158f;hp=4723e0ce408a5510d29644432c8cb096f897123a;hb=a0f36c8849d41b4be54f00e3f9e19ce79a67abae;hpb=13f94397661f67256a68849da13d2029880647c2 diff --git a/energy.c b/energy.c index 4723e0c..dd9206f 100644 --- a/energy.c +++ b/energy.c @@ -13,9 +13,14 @@ 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) +void energy_init(void) { +} + /*---------- main energy computation and subroutines ----------*/ double compute_energy(const struct Vertices *vs) { + static int bests_unprinted; + double energy; int printing; @@ -25,13 +30,12 @@ double compute_energy(const struct Vertices *vs) { printing= printing_check(pr_cost,0); - if (printing) printf("cost > energy |"); + if (printing) printf("%15lld c>e |", evaluations); - COST(3e2, edgewise_vertex_displacement_cost(vs->a)); - COST(1e3, edge_length_variation_cost(vs->a)); - COST(0.2e3, rim_proximity_cost(vs->a)); -// COST(1e2, graph_layout_cost(vs->a)); - COST(1e8, noncircular_rim_cost(vs->a)); + COST( 3e2, line_bending_cost(vs->a)); + COST( 1e3, edge_length_variation_cost(vs->a)); + COST( 0.2e3, rim_proximity_cost(vs->a)); + COST( 1e8, noncircular_rim_cost(vs->a)); if (printing) printf("| total %# e |", energy); @@ -39,12 +43,18 @@ double compute_energy(const struct Vertices *vs) { FILE *best_f; int r; - if (printing) printf(" BEST"); + if (printing) { + printf(" BEST"); + if (bests_unprinted) printf(" [%4d]",bests_unprinted); + bests_unprinted= 0; + } else { + bests_unprinted++; + } - best_f= fopen(output_file_tmp,"wb"); if (!best_f) diee("fopen new out"); + best_f= fopen(best_file_tmp,"wb"); if (!best_f) diee("fopen new out"); 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"); + if (rename(best_file_tmp,best_file)) diee("rename install new best"); best_energy= energy; } @@ -53,12 +63,13 @@ double compute_energy(const struct Vertices *vs) { flushoutput(); } + evaluations++; return energy; } static void addcost(double *energy, double tweight, double tcost, int pr) { double tenergy= tweight * tcost; - if (pr) printf(" %# e > %# e* |", tcost, tenergy); + if (pr) printf(" %# e x %g > %# e* |", tcost, tweight, tenergy); *energy += tenergy; } @@ -66,13 +77,14 @@ static void addcost(double *energy, double tweight, double tcost, int pr) { void compute_edge_lengths(const Vertices vertices) { int v1,e,v2; - + FOR_EDGE(v1,e,v2) edge_lengths[v1][e]= hypotD(vertices[v1],vertices[v2]); } void compute_vertex_areas(const Vertices vertices) { - int v0,v1,v2, e1,e2, k; + int v0,v1,v2, e1,e2; +// int k; FOR_VERTEX(v0) { double total= 0.0, edges_total=0; @@ -85,14 +97,14 @@ void compute_vertex_areas(const Vertices vertices) { edges_total += edge_lengths[v0][e1]; - 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 += magnD(av); - +// 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 += magnD(av); + count++; } vertex_areas[v0]= total / count; @@ -103,62 +115,57 @@ void compute_vertex_areas(const Vertices vertices) { /*---------- Edgewise vertex displacement ----------*/ /* + * Definition: * + * At each vertex Q, in each direction e: * + * e + * Q ----->----- R + * _,-'\__/ + * _,-' delta + * P ' * - * Q `-_ - * / | `-_ - * / | `-. - * / M - - - - - S - * / ' | _,-' - * / ' | _,-' - * / ' , P ' - * / ',-' - * /,-' - * /' - * R - * - * Let delta = 180deg - angle RMS + * r + * cost = delta (we use r=3) + * Q,e * - * Let l = |PQ| - * d = |RS| * - * Giving energy contribution: + * Calculation: * - * 3 - * l delta - * E = F . -------- - * vd, edge PQ vd d + * Let vector A = PQ + * B = QR * + * -1 A . B + * delta = tan ------- + * | A x B | * - * (The dimensions of this are those of F_vd.) + * which is always in the range 0..pi because the denominator + * is nonnegative. We add epsilon to |AxB| to avoid division + * by zero. * - * We calculate delta as atan2(|AxB|, A.B) - * where A = PQ, B = QR - * - * In practice to avoid division by zero we'll add epsilon to d and - * |AxB| and the huge energy ought then to be sufficient for the - * model to avoid being close to R=S. + * r + * cost = delta + * Q,e */ -double edgewise_vertex_displacement_cost(const Vertices vertices) { +double line_bending_cost(const Vertices vertices) { static const double axb_epsilon= 1e-6; + static const double exponent_r= 3; - int pi,e,qi,ri, k; //,si - double a[D3], b[D3], axb[D3]; //m[D3], + int pi,e,qi,ri, k; + double a[D3], b[D3], axb[D3]; double total_cost= 0; FOR_EDGE(qi,e,ri) { pi= EDGE_END2(qi,(e+3)%V6); if (pi<0) continue; -// K m[k]= (vertices[pi][k] + vertices[qi][k]) * 0.5; K a[k]= -vertices[pi][k] + vertices[qi][k]; K b[k]= -vertices[qi][k] + vertices[ri][k]; xprod(axb,a,b); - + double delta= atan2(magnD(axb) + axb_epsilon, dotprod(a,b)); - double cost= pow(delta,3); + double cost= pow(delta,exponent_r); if (!e && !(qi & YMASK)) cost *= 10; @@ -170,17 +177,26 @@ double edgewise_vertex_displacement_cost(const Vertices vertices) { /*---------- edge length variation ----------*/ + /* + * Definition: + * + * See the diagram above. + * r + * cost = ( |PQ| - |QR| ) + * Q,e + */ + double edge_length_variation_cost(const Vertices vertices) { - double diff, cost= 0; - int v0, efwd,vfwd, eback; + double diff, cost= 0, exponent_r= 2; + int q, e,r, eback; - FOR_EDGE(v0,efwd,vfwd) { - eback= edge_reverse(v0,efwd); - diff= edge_lengths[v0][efwd] - edge_lengths[v0][eback]; - cost += diff*diff; + FOR_EDGE(q,e,r) { + eback= edge_reverse(q,e); + diff= edge_lengths[q][e] - edge_lengths[q][eback]; + cost += pow(diff,exponent_r); } return cost; -} +} /*---------- rim proximity cost ----------*/ @@ -223,7 +239,7 @@ double noncircular_rim_cost(const Vertices vertices) { FOR_RIM_VERTEX(vy,vx,v) { find_nearest_oncircle(oncircle, vertices[v]); - + double d2= hypotD2(vertices[v], oncircle); cost += d2*d2; }