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  *
  20  *      _ L-4 ___ L-3 ___ L-2 ___ L-1 ___ 0   ___ 1   ___ 2   ___ 3
  21  *         1       1       1       1      W-2     W-2     W-2     W-2
  22  *       /  \    /  \    /  \    /  \    /  \    /  \    /  \    /
  23  *      /    \  /    \  /    \  /    \  /    \  /    \  /    \  /
  24  *    L-5 ___ L-4 ___ L-3 ___ L-2 ___ L-1 ___ 0   ___ 1   ___ 2   __
  25  *     0       0       0       0       0      W-1     W-1     W-1
  26  *
  27  * Node x,y for
  28  *   0 <= x < L     x = distance along    L = length
  29  *   0 <= y < W     y = distance across   W = width
  30  *
  31  * Vertices are in reading order from diagram above ie x varies fastest.
  32  *
  33  * Where necessary, edges are counted anticlockwise from to-the-right:
  34  *
  35  *                          \2   /1
  36  *                           \  /
  37  *                        ___ 0   __
  38  *                        3    1   0
  39  *                           /  \
  40  *                         4/   5\
  41  *
  42  * When we iterate over edges, we iterate first over vertices and then
  43  * over edges 0 to 2, disregarding edges 3 to 5.
  44  */
  45
  46 #define L 10
  47 #define W 10
  48 #define N (L*W)
  49
  50 static double energy_function(const double vertices[N][3]) {
  51   /* Energy has the following parts right now:
  52    */
  53   /*
  54    *  Edgewise expected vertex displacement
  55    *
  56    *
  57    *                Q `-_
  58    *              / |    `-_
  59    *   R' - _ _ _/_ |       `-.
  60    *    .       /   M - - - - - S
  61    *    .      /    |      _,-'
  62    *    .     /     |  _,-'
  63    *    .    /    , P '
  64    *    .   /  ,-'
  65    *    .  /,-'
  66    *    . /'
  67    *     R
  68    *
  69    *
  70    *
  71    *  Find R', the `expected' location of R, by
  72    *  reflecting S in M (the midpoint of QP).
  73    *
  74    *  Let  d = |RR'|
  75    *       m = |PQ|
  76    *       l = |RS|
  77    *
  78    *  Giving energy contribution:
  79    *
  80    *                               2
  81    *                            m d
  82    *    E             =  F   .  ----
  83    *     vd, edge PQ      vd      l
  84    *
  85    *  By symmetry, this calculation gives the same answer
  86    *  with R and S exchanged.  Looking at the projection in
  87    *  the RMS plane:
  88    *
  89    *
  90    *                           S'
  91    *                         ,'
  92    *                       ,'
  93    *    R'               ,'             d" = |SS'| = |RR'| = d
  94    *      `-._         ,'
  95    *          `-._   ,'
  96    *              ` M
  97    *              ,'  `-._
  98    *            ,'        `-._
  99    *          ,'              ` S
 100    *        ,'
 101    *      ,'
 102    *     R
 103    *
 104    *  We choose this value for l (rather than |RM|+|MS|, say, or |RM|)
 105    *  because we want this symmetry and because we're happy to punish
 106    *  bending more than uneveness in the metric.
 107    *
 108    *  In practice to avoid division by zero we'll add epsilon to l and
 109    *  the huge energy ought then to be sufficient for the model to
 110    *  avoid being close to R=S.
 111    */
 112
 113 static gsl_multimin_fminimizer *minimiser;
 114
 115
 116
 117 {
 118   gsl_multimin_fminimizer_alloc(
 119