chiark / gitweb /
parallel thing compiles now
[moebius2.git] / energy.c
1 /*
2  * We try to find an optimal triangle grid
3  */
4
5 #include "common.h"
6 #include "minimise.h"
7 #include "mgraph.h"
8 #include "parallel.h"
9
10 double vertex_areas[N], vertex_mean_edge_lengths[N], edge_lengths[N][V6];
11
12 static double best_energy= DBL_MAX;
13
14 static void addcost(double *energy, double tweight, double tcost, int pr);
15
16 /*---------- main energy computation, weights, etc. ----------*/
17
18 typedef double CostComputation(const Vertices vertices, int section);
19
20 typedef struct {
21   double weight;
22   CostComputation *fn;
23 } CostContribution;
24
25 static const CostContribution costs[]= {
26 #define COST(weight, compute) { (weight),(compute) },
27
28 #if XBITS==3
29 #define STOP_EPSILON 1e-6;
30     COST(  3e2,   line_bending_cost)
31     COST(  1e3,   edge_length_variation_cost)
32     COST( 0.2e3,  rim_proximity_cost)
33     COST(  1e8,   noncircular_rim_cost)
34 #endif
35
36 #if XBITS==4
37 #define STOP_EPSILON 1e-5;
38     COST(  3e2,   line_bending_cost)
39     COST(  3e3,   edge_length_variation_cost)
40     COST( 3.8e1,  rim_proximity_cost) // 5e1 is too much
41                                          // 2.5e1 is too little
42     // 0.2e1 grows compared to previous ?
43     // 0.6e0 shrinks compared to previous ?
44     COST(  1e12,   noncircular_rim_cost)
45 #endif
46 };
47
48 #define NCOSTS ((sizeof(costs)/sizeof(costs[0])))
49
50 void energy_init(void) {
51   stop_epsilon= STOP_EPSILON;
52 }
53
54 void compute_energy_separately(const struct Vertices *vs,
55                          int section, void *energies_v, void *totals_v) {
56   double *energies= energies_v;
57   int ci;
58     
59   compute_edge_lengths(vs->a, section);
60   compute_vertex_areas(vs->a, section);
61
62   for (ci=0; ci<NCOSTS; ci++)
63     energies[ci]= costs[ci].fn(vs->a, section);
64 }
65
66 /*---------- energy computation machinery ----------*/
67
68 void compute_energy_combine(const struct Vertices *vertices,
69                          int section, void *energies_v, void *totals_v) {
70   int ci;
71   
72   double *energies= energies_v;
73   double *totals= totals_v;
74
75   for (ci=0; ci<NCOSTS; ci++)
76     totals[ci] += energies[ci];
77 }
78
79 double compute_energy(const struct Vertices *vs) {
80   static int bests_unprinted;
81
82   double totals[NCOSTS], energy;
83   int ci, printing;
84
85   printing= printing_check(pr_cost,0);
86
87   if (printing) printf("%15lld c>e |", evaluations);
88
89   inparallel(vs,
90              compute_energy_separately,
91              compute_energy_combine,
92              sizeof(totals) /* really, size of energies */,
93              totals);
94
95   energy= 0;
96
97   for (ci=0; ci<NCOSTS; ci++)
98     addcost(&energy, costs[ci].weight, totals[ci], printing);
99
100   if (printing) printf("| total %# e |", energy);
101
102   if (energy < best_energy) {
103     FILE *best_f;
104     int r;
105
106     if (printing) {
107       printf(" BEST");
108       if (bests_unprinted) printf(" [%4d]",bests_unprinted);
109       bests_unprinted= 0;
110     } else {
111       bests_unprinted++;
112     }
113
114     best_f= fopen(best_file_tmp,"wb");  if (!best_f) diee("fopen new out");
115     r= fwrite(vs->a,sizeof(vs->a),1,best_f); if (r!=1) diee("fwrite");
116     if (fclose(best_f)) diee("fclose new best");
117     if (rename(best_file_tmp,best_file)) diee("rename install new best");
118
119     best_energy= energy;
120   }
121   if (printing) {
122     putchar('\n');
123     flushoutput();
124   }
125
126   evaluations++;
127   return energy;
128 }
129
130 static void addcost(double *energy, double tweight, double tcost, int pr) {
131   double tenergy= tweight * tcost;
132   if (pr) printf(" %# e x %g > %# e* |", tcost, tweight, tenergy);
133   *energy += tenergy;
134 }
135
136 /*---------- Precomputations ----------*/
137
138 void compute_edge_lengths(const Vertices vertices, int section) {
139   int v1,e,v2;
140
141   FOR_EDGE(v1,e,v2,OUTER)
142     edge_lengths[v1][e]= hypotD(vertices[v1],vertices[v2]);
143 }
144
145 void compute_vertex_areas(const Vertices vertices, int section) {
146   int v0,v1,v2, e1,e2;
147 //  int k;
148
149   FOR_VERTEX(v0,OUTER) {
150     double total= 0.0, edges_total=0;
151     int count= 0;
152
153     FOR_VEDGE(v0,e1,v1) {
154       e2= (e1+1) % V6;
155       v2= EDGE_END2(v0,e2);
156       if (v2<0) continue;
157
158       edges_total += edge_lengths[v0][e1];
159
160 //      double e1v[D3], e2v[D3], av[D3];
161 //      K {
162 //      e1v[k]= vertices[v1][k] - vertices[v0][k];
163 //      e2v[k]= vertices[v2][k] - vertices[v0][k];
164 //      }
165 //      xprod(av, e1v, e2v);
166 //      total += magnD(av);
167
168       count++;
169     }
170     vertex_areas[v0]= total / count;
171     vertex_mean_edge_lengths[v0]= edges_total / count;
172   }
173 }
174
175 /*---------- Edgewise vertex displacement ----------*/
176
177   /*
178    * Definition:
179    *
180    *    At each vertex Q, in each direction e:
181    *
182    *                                  e
183    *                           Q ----->----- R
184    *                      _,-'\__/
185    *                  _,-'       delta
186    *               P '
187    *
188    *                      r
189    *       cost    = delta          (we use r=3)
190    *           Q,e
191    *
192    *
193    * Calculation:
194    *
195    *      Let vector A = PQ
196    *                 B = QR
197    *
198    *                   -1   A . B
199    *      delta =  tan     -------
200    *                      | A x B |
201    *
202    *      which is always in the range 0..pi because the denominator
203    *      is nonnegative.  We add epsilon to |AxB| to avoid division
204    *      by zero.
205    *
206    *                     r
207    *      cost    = delta
208    *          Q,e
209    */
210
211 double line_bending_cost(const Vertices vertices, int section) {
212   static const double axb_epsilon= 1e-6;
213   static const double exponent_r= 3;
214
215   int pi,e,qi,ri, k;
216   double  a[D3], b[D3], axb[D3];
217   double total_cost= 0;
218
219   FOR_EDGE(qi,e,ri, OUTER) {
220     pi= EDGE_END2(qi,(e+3)%V6); if (pi<0) continue;
221
222     K a[k]= -vertices[pi][k] + vertices[qi][k];
223     K b[k]= -vertices[qi][k] + vertices[ri][k];
224
225     xprod(axb,a,b);
226
227     double delta= atan2(magnD(axb) + axb_epsilon, dotprod(a,b));
228     double cost= pow(delta,exponent_r);
229
230     if (!e && !(qi & ~XMASK))
231       cost *= 10;
232
233     total_cost += cost;
234   }
235   return total_cost;
236 }
237
238 /*---------- edge length variation ----------*/
239
240   /*
241    * Definition:
242    *
243    *    See the diagram above.
244    *                                r
245    *       cost    = ( |PQ| - |QR| )
246    *           Q,e
247    */
248
249 double edge_length_variation_cost(const Vertices vertices, int section) {
250   double diff, cost= 0, exponent_r= 2;
251   int q, e,r, eback;
252
253   FOR_EDGE(q,e,r, OUTER) {
254     eback= edge_reverse(q,e);
255     diff= edge_lengths[q][e] - edge_lengths[q][eback];
256     cost += pow(diff,exponent_r);
257   }
258   return cost;
259 }
260
261 /*---------- rim proximity cost ----------*/
262
263 static void find_nearest_oncircle(double oncircle[D3], const double p[D3]) {
264   /* By symmetry, nearest point on circle is the one with
265    * the same angle subtended at the z axis. */
266   oncircle[0]= p[0];
267   oncircle[1]= p[1];
268   oncircle[2]= 0;
269   double mult= 1.0/ magnD(oncircle);
270   oncircle[0] *= mult;
271   oncircle[1] *= mult;
272 }
273
274 double rim_proximity_cost(const Vertices vertices, int section) {
275   double oncircle[3], cost=0;
276   int v;
277
278   FOR_VERTEX(v, OUTER) {
279     int y= v >> YSHIFT;
280     int nominal_edge_distance= y <= Y/2 ? y : Y-1-y;
281     if (nominal_edge_distance==0) continue;
282
283     find_nearest_oncircle(oncircle, vertices[v]);
284
285     cost +=
286       vertex_mean_edge_lengths[v] *
287       (nominal_edge_distance*nominal_edge_distance) /
288       (hypotD2(vertices[v], oncircle) + 1e-6);
289   }
290   return cost;
291 }
292
293 /*---------- noncircular rim cost ----------*/
294
295 double noncircular_rim_cost(const Vertices vertices, int section) {
296   int vy,vx,v;
297   double cost= 0.0;
298   double oncircle[3];
299
300   FOR_RIM_VERTEX(vy,vx,v, OUTER) {
301     find_nearest_oncircle(oncircle, vertices[v]);
302
303     double d2= hypotD2(vertices[v], oncircle);
304     cost += d2*d2;
305   }
306   return cost;
307 }