anneal.c

   1 /*
   2  * We try to find an optimal triangle grid
   3  *
   4  * Vertices in strip are numbered as follows:
   5  *
   6  *     ___ X-2 ___ X-1 ___ 0   ___ 1   ___ 2   ___ 3   ___ 4   __
   7  *         Y-1     Y-1      0       0       0       0       0
   8  *        /  \    /  \    /  \    /  \    /  \    /  \    /  \
   9  *       /    \  /    \  /    \  /    \  /    \  /    \  /    \
  10  *     X-3 ___ X-2 ___ X-1 ___ 0   ___ 1   ___ 2   ___  3  ___ 4
  11  *     Y-2     Y-2     Y-2      1       1       1       1       1
  12  *       \    /  \    /  \    /  \    /  \    /  \    /  \    /
  13  *        \  /    \  /    \  /    \  /    \  /    \  /    \  /
  14  *     ___ X-3 ___ X-2 ___ X-1 ___ 0   ___ 1   ___ 2   ___ 3   __
  15  *         Y-3     Y-3     Y-3      2       2       2       2
  16  *
  17  *       .   .   .   .   .   .   .   .   .   .   .   .   .   .   .
  18  *
  19  *      X-4 ___ X-3 ___ X-2 ___ X-1 ___ 0   ___ 1   ___ 2   ___ 3
  20  *       1       1       1       1      Y-2     Y-2     Y-2     Y-2
  21  *        \    /  \    /  \    /  \    /  \    /  \    /  \    /
  22  *         \  /    \  /    \  /    \  /    \  /    \  /    \  /
  23  *      ___ X-4 ___ X-3 ___ X-2 ___ X-1 ___ 0   ___ 1   ___ 2   __
  24  *           0       0       0       0      Y-1     Y-1     Y-1
  25  *
  26  * Node x,y for
  27  *   0 <= x < X     x = distance along
  28  *   0 <= y < Y     y = distance across
  29  *
  30  * Vertices are in reading order from diagram above ie x varies fastest.
  31  *
  32  * Y must be even.  The actual location opposite (0,0) is (X-(Y-1)/2,0),
  33  * and likewise opposite (0,Y-1) is ((Y-1)/2,0).
  34  *
  35  * To label edges, we counte anticlockwise[*] from to-the-right:
  36  *
  37  *                          \2   /1
  38  *                           \  /
  39  *                        ___ 0   __
  40  *                        3    1   0
  41  *                           /  \
  42  *                         4/   5\
  43  *
  44  *   [*] That is, in the direction represented as anticlockwise for
  45  *       the vertices (0,*)..(4,*) in the diagram above; and of course
  46  *       that is clockwise for the vertices (X-5,*)..(X-1,*).  The
  47  *       numbering has an actual discontinuity between (X-1,*) and (0,*).
  48  *
  49  * When we iterate over edges, we iterate first over vertices and then
  50  * over edges 0 to 2, disregarding edges 3 to 5.
  51  */
  52
  53 #define XBITS 4
  54 #define X (1<<XBITS)
  55 #define YBITS 4
  56 #define Y (1<<YBITS)
  57
  58 /* vertex number:   0000 | y     | x
  59  *                        YBITS   XBITS
  60  */
  61
  62 #define N (X*Y)
  63 #define XMASK (X-1)
  64 #define YSHIFT XBITS
  65 #define Y1 (1 << YSHIFT)
  66 #define YMASK (Y-1 << YSHIFT)
  67
  68   /* Energy has the following parts right now:
  69    */
  70   /*
  71    *  Edgewise expected vertex displacement
  72    *
  73    *
  74    *                Q `-_
  75    *              / |    `-_
  76    *   R' - _ _ _/_ |       `-.
  77    *    .       /   M - - - - - S
  78    *    .      /    |      _,-'
  79    *    .     /     |  _,-'
  80    *    .    /    , P '
  81    *    .   /  ,-'
  82    *    .  /,-'
  83    *    . /'
  84    *     R
  85    *
  86    *
  87    *
  88    *  Find R', the `expected' location of R, by
  89    *  reflecting S in M (the midpoint of QP).
  90    *
  91    *  Let 2d = |RR'|
  92    *       b = |PQ|
  93    *       l = |RS|
  94    *
  95    *  Giving energy contribution:
  96    *
  97    *                               2
  98    *                            b d
  99    *    E             =  F   .  ----
 100    *     vd, edge PQ      vd      3
 101    *                             l
 102    *
 103    *  By symmetry, this calculation gives the same answer
 104    *  with R and S exchanged.  Looking at the projection in
 105    *  the RMS plane:
 106    *
 107    *
 108    *                           S'
 109    *                         ,'
 110    *                       ,'
 111    *    R'               ,'             2d" = |SS'| = |RR'| = 2d
 112    *      `-._         ,'
 113    *          `-._   ,'                 By congruent triangles,
 114    *              ` M                   with M' = midpoint of RS,
 115    *              ,'  `-._              |MM'| = |RR'|/2 = d
 116    *            ,'        `-._
 117    *          ,'              ` S       So use
 118    *        ,'       M' _ , - '            d = |MM'|
 119    *      ,'   _ , - '
 120    *     R - '
 121    *
 122    *  We choose this value for l (rather than |RM|+|MS|, say, or |RM|)
 123    *  because we want this symmetry and because we're happy to punish
 124    *  bending more than uneveness in the metric.
 125    *
 126    *  In practice to avoid division by zero we'll add epsilon to l^3
 127    *  and the huge energy ought then to be sufficient for the model to
 128    *  avoid being close to R=S.  */
 129
 130 #define V6 6
 131
 132 static int edge_end2(unsigned v1, int e) {
 133   /* The topology is equivalent to that of a square lattice with only
 134    * half of the diagonals.  Ie, the result of shearing the triangular
 135    * lattice to make the lines of constant x vertical.  This gives
 136    * these six directions:
 137    *
 138    *                2  1
 139    *                | /
 140    *                |/
 141    *             3--*--0
 142    *               /|
 143    *              / |
 144    *             4  5
 145    *
 146    * This also handily makes vertical the numbering discontinuity,
 147    * where the join happens.
 148    */
 149   static const unsigned dx[V6]= {  +1,  +1,   0,  -1,  -1,   0  },
 150                         dy[V6]= {   0, +Y1, +Y1,   0, -Y1, -Y1  };
 151   unsigned x, y;
 152
 153   y= (v1 & YMASK) + dy[e];
 154   if (y & ~YMASK) return -1;
 155
 156   x= (v1 & XMASK) + dx[e];
 157   if (x & ~XMASK) {
 158     y= (Y-1)*Y1 - y;
 159     x &= XMASK;;
 160   }
 161   return x | y;
 162 }
 163
 164 #define D3 3
 165
 166 #define FOR_VERTEX(v) \
 167   for ((v)=0; (v)<N; (v)++)
 168
 169 #define FOR_VPEDGE(v,e) \
 170   for ((e)=0; (e)<V6; (e)++)
 171
 172 #define EDGE_END2 edge_end2
 173
 174 #define FOR_VEDGE(v1,e,v2) \
 175   FOR_VPEDGE((v1),(e))
 176     if (((v2)= EDGE_END2((v1),(e))) < 0) ; else
 177
 178 #define FOR_EDGE(v1,e,v2) \
 179   FOR_VERTEX((v1)) \
 180     FOR_VEDGE((v1),(e),(v2))
 181
 182 #define FOR_COORD(k) \
 183   for ((k)=0; (k)<D3; (k)++)
 184
 185 #define K FOR_COORD(k)
 186
 187 static double hypotD(const double p[], const double q[]) {
 188   int k;
 189   double pq[D3];
 190   gsl_vector v;
 191
 192   K pq[_k]= p[k] - q[k];
 193   v.size= D3;
 194   v.stride= 1;
 195   v.vector.data= pq;
 196   /* owner and block ought not to be used */
 197
 198   return gsl_blas_snrm2(&v);
 199 }
 200
 201 #ifdef FP_FAST_FMA
 202 # define ffsqa(factor,term) fma((factor),(factor),(term))
 203 #else
 204 # define ffsqu(factor,term) ((factor)*(factor)+(term))
 205 #endif
 206
 207 static const l3_epsison= 1e-6;
 208
 209 static double energy_function(const double vertices[N][D3]) {
 210   int pi,e,qi,ri,si, k;
 211   double m[D3], mprime[D3], b, d2, l, sigma_bd2_l;
 212
 213   FOR_EDGE(pi,e,qi) {
 214     ri= EDGE_END2(pi,(e+1)%V6); if (r<0) continue;
 215     si= EDGE_END2(pi,(e+5)%V6); if (s<0) continue;
 216     assert(ri == EDGE_END2(qi,(e+2)%V6));
 217     assert(si == EDGE_END2(qi,(e+4)%V6));
 218
 219     K m[k]= (vertices[pi][k] + vertices[qi][k]) * 0.5;
 220     K mprime[k]= (vertices[ri][k] + vertices[si][k]) * 0.5;
 221     b= hypotD(vertices[pi], vertices[qi]);
 222     d2= 0; K d2= ffsqa(m[k] - mprime[k], d2);
 223     l= hypotD(vertices[ri][k] - vertices[si][k]);
 224     l3 = l*l*l + l3_epsilon;
 225
 226     sigma_bd2_l += b * d2 / l3;
 227   }
 228
 229
 230
 231 static gsl_multimin_fminimizer *minimiser;
 232
 233
 234
 235 {
 236   gsl_multimin_fminimizer_alloc(
 237