best so far perhaps

diff --git a/energy.c b/energy.c
index 3459a726351724226f4fe7994c3fdc75d5b4d909..9c351a52b9349b327547b7122172944071a9e0d3 100644 (file)
--- a/energy.c
+++ b/energy.c
@@ -1,180 +1,182 @@
-/*
- * We try to find an optimal triangle grid
- */
-
-#include "common.h"
-#include "minimise.h"
-#include "mgraph.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)
-
-void energy_init(void) {
-}
-
-/*---------- main energy computation and subroutines ----------*/
-
-double compute_energy(const struct Vertices *vs) {
-  static int bests_unprinted;
-
-  double energy;
-  int printing;
-
-  compute_edge_lengths(vs->a);
-  compute_vertex_areas(vs->a);
-  energy= 0;
-
-  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();
-  }
-
-  if (printing) printf("| total %# e |", energy);
-
-  if (energy < best_energy) {
-    FILE *best_f;
-    int r;
-
-    if (printing) {
-      printf(" BEST");
-      if (bests_unprinted) printf(" [%4d]",bests_unprinted);
-      bests_unprinted= 0;
-    } else {
-      bests_unprinted++;
-    }
-
-    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(best_file_tmp,best_file)) diee("rename install new best");
-
-    best_energy= energy;
-  }
-  if (printing) {
-    putchar('\n');
-    flushoutput();
-  }
-
-  evaluations++;
-  return energy;
-}
-
-static void addcost(double *energy, double tweight, double tcost, int pr) {
-  double tenergy= tweight * tcost;
-  if (pr) printf(" %# e x %g > %# e* |", tcost, tweight, tenergy);
-  *energy += tenergy;
-}
-
-/*---------- Precomputations ----------*/
-
-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;
-//  int k;
-
-  FOR_VERTEX(v0) {
-    double total= 0.0, edges_total=0;
-    int count= 0;
-
-    FOR_VEDGE(v0,e1,v1) {
-      e2= (e1+1) % V6;
-      v2= EDGE_END2(v0,e2);
-      if (v2<0) continue;
-
-      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);
-
-      count++;
-    }
-    vertex_areas[v0]= total / count;
-    vertex_mean_edge_lengths[v0]= edges_total / count;
-  }
-}
-
-/*---------- Edgewise vertex displacement ----------*/
-
-  /*
-   * Definition:
-   *
-   *    At each vertex Q, in each direction e:
-   *
-   *                                         e
-   *                           Q ----->----- R
-   *                     _,-'\__/
-   *                         _,-'       delta
-   *              P '
-   *
-   *                      r
-   *       cost    = delta          (we use r=3)
-   *           Q,e
-   *
-   *
-   * Calculation:
-   *
-   *      Let vector A = PQ
-   *                 B = QR
-   *
-   *                          -1   A . B
-   *      delta =  tan     -------
-   *                     | A x B |
-   *
-   *      which is always in the range 0..pi because the denominator
-   *      is nonnegative.  We add epsilon to |AxB| to avoid division
-   *      by zero.
-   *
-   *                     r
-   *      cost    = delta
-   *          Q,e
-   */
-
-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;
-  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;
+ /*
+  * We try to find an optimal triangle grid
+  */
+
+ #include "common.h"
+ #include "minimise.h"
+ #include "mgraph.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)
+
+ void energy_init(void) {
+ }
+
+ /*---------- main energy computation and subroutines ----------*/
+
+ double compute_energy(const struct Vertices *vs) {
+   static int bests_unprinted;
+
+   double energy;
+   int printing;
+
+   compute_edge_lengths(vs->a);
+   compute_vertex_areas(vs->a);
+   energy= 0;
+
+   printing= printing_check(pr_cost,0);
+
+   if (printing) printf("%15lld c>e |", evaluations);
+
+   if (XBITS==3) {
+     COST(  3e3,   line_bending_cost(vs->a));
+     COST(  3e3,   edge_length_variation_cost(vs->a));
+     COST( 0.4e3,  rim_proximity_cost(vs->a));
+     COST(  1e6,   edge_angle_cost(vs->a, 0.5/1.7));
+ //    COST(  1e1,   small_triangles_cost(vs->a));
+     COST(  1e12,   noncircular_rim_cost(vs->a));
+     stop_epsilon= 1e-6;
+   } else if (XBITS==4) {
+     COST(  3e5,   line_bending_cost(vs->a));
+     COST( 10e2,   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, 0.5/1.3));
+     COST(  1e18,   noncircular_rim_cost(vs->a));
+     stop_epsilon= 1e-6;
+   } else {
+     abort();
+   }
+
+   if (printing) printf("| total %# e |", energy);
+
+   if (energy < best_energy) {
+     FILE *best_f;
+     int r;
+
+     if (printing) {
+       printf(" BEST");
+       if (bests_unprinted) printf(" [%4d]",bests_unprinted);
+       bests_unprinted= 0;
+     } else {
+       bests_unprinted++;
+     }
+
+     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(best_file_tmp,best_file)) diee("rename install new best");
+
+     best_energy= energy;
+   }
+   if (printing) {
+     putchar('\n');
+     flushoutput();
+   }
+
+   evaluations++;
+   return energy;
+ }
+
+ static void addcost(double *energy, double tweight, double tcost, int pr) {
+   double tenergy= tweight * tcost;
+   if (pr) printf(" %# e x %g > %# e* |", tcost, tweight, tenergy);
+   *energy += tenergy;
+ }
+
+ /*---------- Precomputations ----------*/
+
+ 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;
+ //  int k;
+
+   FOR_VERTEX(v0) {
+     double total= 0.0, edges_total=0;
+     int count= 0;
+
+     FOR_VEDGE(v0,e1,v1) {
+       e2= (e1+1) % V6;
+       v2= EDGE_END2(v0,e2);
+       if (v2<0) continue;
+
+       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);
+
+       count++;
+     }
+     vertex_areas[v0]= total / count;
+     vertex_mean_edge_lengths[v0]= edges_total / count;
+   }
+ }
+
+ /*---------- Edgewise vertex displacement ----------*/
+
+   /*
+    * Definition:
+    *
+    *    At each vertex Q, in each direction e:
+    *
+    *                                        e
+    *                           Q ----->----- R
+    *                    _,-'\__/
+    *                        _,-'       delta
+    *             P '
+    *
+    *                      r
+    *       cost    = delta          (we use r=3)
+    *           Q,e
+    *
+    *
+    * Calculation:
+    *
+    *      Let vector A = PQ
+    *                 B = QR
+    *
+    *                         -1   A . B
+    *      delta =  tan     -------
+    *                    | A x B |
+    *
+    *      which is always in the range 0..pi because the denominator
+    *      is nonnegative.  We add epsilon to |AxB| to avoid division
+    *      by zero.
+    *
+    *                     r
+    *      cost    = delta
+    *          Q,e
+    */
+
+ double line_bending_cost(const Vertices vertices) {
+   static const double axb_epsilon= 1e-6;
+   static const double exponent_r= 4;
+
+   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;
+
+//if (!(qi&XMASK)) fprintf(stderr,"%02x-%02x-%02x (%d)\n",pi,qi,ri,e);

 
     K a[k]= -vertices[pi][k] + vertices[qi][k];
     K b[k]= -vertices[qi][k] + vertices[ri][k];
 
     K a[k]= -vertices[pi][k] + vertices[qi][k];
     K b[k]= -vertices[qi][k] + vertices[ri][k];


@@ -184,9 +186,6 @@ double line_bending_cost(const Vertices vertices) {
     double delta= atan2(magnD(axb) + axb_epsilon, dotprod(a,b));
     double cost= pow(delta,exponent_r);
 
     double delta= atan2(magnD(axb) + axb_epsilon, dotprod(a,b));
     double cost= pow(delta,exponent_r);
 

-    if (!e && !(qi & ~XMASK))
-      cost *= 10;
-

     total_cost += cost;
   }
   return total_cost;
     total_cost += cost;
   }
   return total_cost;


@@ -290,10 +289,10 @@ double noncircular_rim_cost(const Vertices vertices) {
    *     vd, edge PQ      vd
    */
 
    *     vd, edge PQ      vd
    */
 

-double edge_angle_cost(const Vertices vertices) {
+double edge_angle_cost(const Vertices vertices, double circcircrat) {

   double pq1[D3], rp[D3], ps[D3], rp_2d[D3], ps_2d[D3], rs_2d[D3];
   double a,b,c,s,r;
   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;
 
   int pi,e,qi,ri,si, k;
 //  double our_epsilon=1e-6;


@@ -322,8 +321,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);
     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 + circcircrat*(a+b);

     double deficit= minradius - r;
     if (deficit < 0) continue;
     double cost= deficit*deficit;
     double deficit= minradius - r;
     if (deficit < 0) continue;
     double cost= deficit*deficit;

diff --git a/minimise.h b/minimise.h
index 60efb3e53a07d83985cb277e927d400897d4c988..a8a2288f5c14b38a1993ab84b4133070ce317fbc 100644 (file)
--- a/minimise.h
+++ b/minimise.h
@@ -21,7 +21,7 @@ double line_bending_cost(const Vertices vertices);
 double noncircular_rim_cost(const Vertices vertices);
 double edge_length_variation_cost(const Vertices vertices);
 double rim_proximity_cost(const Vertices vertices);
 double noncircular_rim_cost(const Vertices vertices);
 double edge_length_variation_cost(const Vertices vertices);
 double rim_proximity_cost(const Vertices vertices);

-double edge_angle_cost(const Vertices vertices);
+double edge_angle_cost(const Vertices vertices, double circcircrat);

 double small_triangles_cost(const Vertices vertices);
 
 extern const char *input_file, *best_file;
 double small_triangles_cost(const Vertices vertices);
 
 extern const char *input_file, *best_file;