anneal.c

   1 /*
   2  * We try to find an optimal triangle grid
   3  *
   4  * Vertices in strip are numbered as follows:
   5  *
   6  *     ___ L-2 ___ L-1 ___ 0   ___ 1   ___ 2   ___ 3   ___ 4   __
   7  *         W-1     W-1      0       0       0       0       0
   8  *        /  \    /  \    /  \    /  \    /  \    /  \    /  \
   9  *       /    \  /    \  /    \  /    \  /    \  /    \  /    \
  10  *     L-3 ___ L-2 ___ L-1 ___ 0   ___ 1   ___ 2   ___  3  ___ 4
  11  *     W-2     W-2     W-2      1       1       1       1       1
  12  *       \    /  \    /  \    /  \    /  \    /  \    /  \    /
  13  *        \  /    \  /    \  /    \  /    \  /    \  /    \  /
  14  *     ___ L-3 ___ L-2 ___ L-1 ___ 0   ___ 1   ___ 2   ___ 3   __
  15  *         W-3     W-3     W-3      2       2       2       2
  16  *
  17  *       .   .   .   .   .   .   .   .   .   .   .   .   .   .   .
  18  *
  19  *      L-4 ___ L-3 ___ L-2 ___ L-1 ___ 0   ___ 1   ___ 2   ___ 3
  20  *       1       1       1       1      W-2     W-2     W-2     W-2
  21  *        \    /  \    /  \    /  \    /  \    /  \    /  \    /
  22  *         \  /    \  /    \  /    \  /    \  /    \  /    \  /
  23  *      ___ L-4 ___ L-3 ___ L-2 ___ L-1 ___ 0   ___ 1   ___ 2   __
  24  *           0       0       0       0      W-1     W-1     W-1
  25  *
  26  * Node x,y for
  27  *   0 <= x < L     x = distance along    L = length
  28  *   0 <= y < W     y = distance across   W = width
  29  *
  30  * Vertices are in reading order from diagram above ie x varies fastest.
  31  *
  32  * W must be even.  The actual location opposite (0,0) is (L-(W-1)/2,0),
  33  * and likewise opposite (0,W-1) is ((W-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 (L-5,*)..(L-1,*).  The
  47  *       numbering has an actual discontinuity between (L-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 LN2W 4
  54 #define L (1<<LN2W)
  55 #define W 10
  56 #define N (L*W)
  57 #define YMASK (~0u << LN2W)
  58 #define XMASK (~YMASK)
  59
  60   /* Energy has the following parts right now:
  61    */
  62   /*
  63    *  Edgewise expected vertex displacement
  64    *
  65    *
  66    *                Q `-_
  67    *              / |    `-_
  68    *   R' - _ _ _/_ |       `-.
  69    *    .       /   M - - - - - S
  70    *    .      /    |      _,-'
  71    *    .     /     |  _,-'
  72    *    .    /    , P '
  73    *    .   /  ,-'
  74    *    .  /,-'
  75    *    . /'
  76    *     R
  77    *
  78    *
  79    *
  80    *  Find R', the `expected' location of R, by
  81    *  reflecting S in M (the midpoint of QP).
  82    *
  83    *  Let  d = |RR'|
  84    *       m = |PQ|
  85    *       l = |RS|
  86    *
  87    *  Giving energy contribution:
  88    *
  89    *                               2
  90    *                            m d
  91    *    E             =  F   .  ----
  92    *     vd, edge PQ      vd      l
  93    *
  94    *  By symmetry, this calculation gives the same answer
  95    *  with R and S exchanged.  Looking at the projection in
  96    *  the RMS plane:
  97    *
  98    *
  99    *                           S'
 100    *                         ,'
 101    *                       ,'
 102    *    R'               ,'             d" = |SS'| = |RR'| = d
 103    *      `-._         ,'
 104    *          `-._   ,'
 105    *              ` M
 106    *              ,'  `-._
 107    *            ,'        `-._
 108    *          ,'              ` S
 109    *        ,'
 110    *      ,'
 111    *     R
 112    *
 113    *  We choose this value for l (rather than |RM|+|MS|, say, or |RM|)
 114    *  because we want this symmetry and because we're happy to punish
 115    *  bending more than uneveness in the metric.
 116    *
 117    *  In practice to avoid division by zero we'll add epsilon to l and
 118    *  the huge energy ought then to be sufficient for the model to
 119    *  avoid being close to R=S.
 120    */
 121
 122 static int edge_end2(int v1, int e) {
 123   if (!(v1 &
 124
 125   if (e==1 || e==2)
 126     if (v2 < L)
 127
 128
 129 #define FOR_VERTEX(v) \
 130   for ((v)=0; (v)<N; (v)++)
 131
 132 #define FOR_VPEDGE(v,e) \
 133   for ((e)=0; (e)<6; (e)++)
 134
 135 #define EDGE_END2 edge_end2
 136
 137 #define FOR_VEDGE(v1,e,v2) \
 138   FOR_VPEDGE((v1),(e))
 139     if (((v2)= EDGE_END2((v1),(e))) >= 0) ; else
 140
 141 #define FOR_EDGE(v1,e,v2) \
 142   FOR_VERTEX((v1)) \
 143     FOR_VEDGE((v1),(e),(v2))
 144
 145 static double energy_function(const double vertices[N][3]) {
 146   int vertex;
 147
 148   FOR_VERTEX {
 149
 150
 151 static gsl_multimin_fminimizer *minimiser;
 152
 153
 154
 155 {
 156   gsl_multimin_fminimizer_alloc(
 157