X-Git-Url: http://www.chiark.greenend.org.uk/ucgi/~ian/git?a=blobdiff_plain;ds=sidebyside;f=anneal.c;h=ab521f4700da30b1b9fa0d5989ffa6d6d89665f5;hb=7904a232f94d57fb0884725481d0c616f02c2520;hp=e13c0b0cb8faaa4d03515148a5f32cf035c1ecd0;hpb=e38ba6e5733db1faab4c146fbd9b3114e10b29f1;p=moebius2.git

diff --git a/anneal.c b/anneal.c
index e13c0b0..ab521f4 100644
--- a/anneal.c
+++ b/anneal.c
@@ -3,34 +3,36 @@
  *
  * Vertices in strip are numbered as follows:
  *
- *     ___ L-2 ___ L-1 ___ 0   ___ 1   ___ 2   ___ 3   ___ 4   __
- *         W-1     W-1      0       0       0       0       0
+ *     ___ X-2 ___ X-1 ___ 0   ___ 1   ___ 2   ___ 3   ___ 4   __
+ *         Y-1     Y-1      0       0       0       0       0
  *        /  \    /  \    /  \    /  \    /  \    /  \    /  \
  *       /    \  /    \  /    \  /    \  /    \  /    \  /    \
- *     L-3 ___ L-2 ___ L-1 ___ 0   ___ 1   ___ 2   ___  3  ___ 4
- *     W-2     W-2     W-2      1       1       1       1       1
+ *     X-3 ___ X-2 ___ X-1 ___ 0   ___ 1   ___ 2   ___  3  ___ 4
+ *     Y-2     Y-2     Y-2      1       1       1       1       1
  *       \    /  \    /  \    /  \    /  \    /  \    /  \    /
  *        \  /    \  /    \  /    \  /    \  /    \  /    \  /
- *     ___ L-3 ___ L-2 ___ L-1 ___ 0   ___ 1   ___ 2   ___ 3   __
- *         W-3     W-3     W-3      2       2       2       2
- *        /  \    /  \    /  \    /  \    /  \    /  \    /  \
- *       /    \  /    \  /    \  /    \  /    \  /    \  /    \
- *    .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .
+ *     ___ X-3 ___ X-2 ___ X-1 ___ 0   ___ 1   ___ 2   ___ 3   __
+ *         Y-3     Y-3     Y-3      2       2       2       2
+ *
+ *       .   .   .   .   .   .   .   .   .   .   .   .   .   .   .
  *
- *      _ L-4 ___ L-3 ___ L-2 ___ L-1 ___ 0   ___ 1   ___ 2   ___ 3
- *         1       1       1       1      W-2     W-2     W-2     W-2
- *       /  \    /  \    /  \    /  \    /  \    /  \    /  \    /
- *      /    \  /    \  /    \  /    \  /    \  /    \  /    \  /
- *    L-5 ___ L-4 ___ L-3 ___ L-2 ___ L-1 ___ 0   ___ 1   ___ 2   __
- *     0       0       0       0       0      W-1     W-1     W-1
+ *      X-4 ___ X-3 ___ X-2 ___ X-1 ___ 0   ___ 1   ___ 2   ___ 3
+ *       1       1       1       1      Y-2     Y-2     Y-2     Y-2
+ *        \    /  \    /  \    /  \    /  \    /  \    /  \    /
+ *         \  /    \  /    \  /    \  /    \  /    \  /    \  /
+ *      ___ X-4 ___ X-3 ___ X-2 ___ X-1 ___ 0   ___ 1   ___ 2   __
+ *           0       0       0       0      Y-1     Y-1     Y-1
  *
  * Node x,y for
- *   0 <= x < L     x = distance along    L = length
- *   0 <= y < W     y = distance across   W = width
+ *   0 <= x < X     x = distance along
+ *   0 <= y < Y     y = distance across
  *
  * Vertices are in reading order from diagram above ie x varies fastest.
  *
- * Where necessary, edges are counted anticlockwise from to-the-right:
+ * Y must be even.  The actual location opposite (0,0) is (X-(Y-1)/2,0),
+ * and likewise opposite (0,Y-1) is ((Y-1)/2,0).
+ *
+ * To label edges, we counte anticlockwise[*] from to-the-right:
  *
  *			    \2   /1
  *			     \  /
@@ -39,75 +41,189 @@
  *			     /  \
  *			   4/   5\
  *
+ *   [*] That is, in the direction represented as anticlockwise for
+ *       the vertices (0,*)..(4,*) in the diagram above; and of course
+ *       that is clockwise for the vertices (X-5,*)..(X-1,*).  The
+ *       numbering has an actual discontinuity between (X-1,*) and (0,*).
+ *
  * When we iterate over edges, we iterate first over vertices and then
  * over edges 0 to 2, disregarding edges 3 to 5.
  */
 
-#define L 10
-#define W 10
-#define N (L*W)
+#define XBITS 4
+#define X (1<<XBITS)
+#define YBITS 4
+#define Y (1<<YBITS)
+
+/* vertex number:   0000 | y     | x
+ *                        YBITS   XBITS
+ */
+
+#define N (X*Y)
+#define XMASK (X-1)
+#define YSHIFT XBITS
+#define Y1 (1 << YSHIFT)
+#define YMASK (Y-1 << YSHIFT)
 
-static double energy_function(const double vertices[N][3]) {
   /* Energy has the following parts right now:
-   *  - worst curvature radius
-   *  - total angle subtended
    */
   /*
+   *  Edgewise expected vertex displacement
    *
-   *	            Q
-   *              / | `-.
-   *             /  |    `-.
-   *            /   P ------ S
-   *           /.-'
-   *          /'
-   *         R
-   */
-  /*
-   * Curvature radius:
    *
-   *  First flatten into plane _|_ PQ:
+   *   	            Q `-_
+   *              / |    `-_
+   *   R' - _ _ _/_ |       `-.
+   *    .       /   M - - - - - S
+   *    .      /    |      _,-'
+   *    .     /     |  _,-'
+   *    .    /    , P '
+   *    .   /  ,-'
+   *    .  /,-'
+   *    . /'
+   *   	 R
    *
-   *   R' = (R-P) x (Q-P)/|Q-P|
-   *   S' = (S-P) x (Q-P)/|Q-P|
-   *   Q' = (Q-P) x (Q-P)/|Q-P| = 0 = P'
    *
    *
-   *	      R'.
-   *   	       	 \.
-   *		   \.
-   *                 0
-   *                 |
-   *                 |
-   *   		     |
-   *		     S'
+   *  Find R', the `expected' location of R, by
+   *  reflecting S in M (the midpoint of QP).
    *
-   *  Now radius formula (wikipedia `Circle'), P1=R', P2=0, P3=S':
+   *  Let 2d = |RR'|
+   *       b = |PQ|
+   *       l = |RS|
    *
-   *          |R'| * |S'| * | R' - S' |
-   *     r =  -------------------------
-   *               2 * | R' x S' |
+   *  Giving energy contribution:
    *
-   *  Energy
-   *                            2
-   *     E  = C  min           r
-   *      r    r   j in edges   j
+   *                               2
+   *                            b d
+   *    E             =  F   .  ----
+   *     vd, edge PQ      vd      l
    *
-   */
-  /*
-   * Angle subtended:
+   *  By symmetry, this calculation gives the same answer
+   *  with R and S exchanged.  Looking at the projection in
+   *  the RMS plane:
    *
-   *                -1   | R' . S' |
-   *    gamma =  cos     -----------
-   *                     |R'| * |S'|
    *
-   *    l = | P - Q |
+   *	       	       	       S'
+   *	   	    	     ,'
+   *	       	       	   ,'
+   *   	R'           	 ,'    	       	2d" = |SS'| = |RR'| = 2d
+   *	  `-._         ,'
+   *   	      `-._   ,'	      		By congruent triangles,
+   *              ` M	      		with M' = midpoint of RS,
+   *	       	  ,'  `-._     	       	|MM'| = |RR'|/2 = d
+   *	       	,'     	  `-._
+   *	      ,'       	      ` S       So use
+   *	    ,'	     M' _ , - '            d = |MM'|
+   *	  ,'   _ , - '
+   *	 R - '
    *
-   *  Energy
+   *  We choose this value for l (rather than |RM|+|MS|, say, or |RM|)
+   *  because we want this symmetry and because we're happy to punish
+   *  bending more than uneveness in the metric.
    *
-   *    E      = C      sum           gamma  * l
-   *     gamma    gamma   j in edges       j    j
+   *  In practice to avoid division by zero we'll add epsilon to l and
+   *  the huge energy ought then to be sufficient for the model to
+   *  avoid being close to R=S.
+   */
+
+#define V6 6
+
+static int edge_end2(unsigned v1, int e) {
+  /* The topology is equivalent to that of a square lattice with only
+   * half of the diagonals.  Ie, the result of shearing the triangular
+   * lattice to make the lines of constant x vertical.  This gives
+   * these six directions:
+   *
+   *                2  1
+   *                | /
+   *                |/
+   *             3--*--0
+   * 		   /|
+   * 		  / |
+   * 		 4  5
    *
+   * This also handily makes vertical the numbering discontinuity,
+   * where the join happens.
    */
+  static const unsigned dx[V6]= {  +1,  +1,   0,  -1,  -1,   0  },
+                        dy[V6]= {   0, +Y1, +Y1,   0, -Y1, -Y1  };
+  unsigned x, y;
+
+  y= (v1 & YMASK) + dy[e];
+  if (y & ~YMASK) return -1;
+
+  x= (v1 & XMASK) + dx[e];
+  if (x & ~XMASK) {
+    y= (Y-1)*Y1 - y;
+    x &= XMASK;;
+  }
+  return x | y;
+}
+
+#define D3 3
+
+#define FOR_VERTEX(v) \
+  for ((v)=0; (v)<N; (v)++)
+
+#define FOR_VPEDGE(v,e) \
+  for ((e)=0; (e)<V6; (e)++)
+
+#define EDGE_END2 edge_end2
+
+#define FOR_VEDGE(v1,e,v2) \
+  FOR_VPEDGE((v1),(e))
+    if (((v2)= EDGE_END2((v1),(e))) < 0) ; else
+
+#define FOR_EDGE(v1,e,v2) \
+  FOR_VERTEX((v1)) \
+    FOR_VEDGE((v1),(e),(v2))
+
+#define FOR_COORD(k) \
+  for ((k)=0; (k)<D3; (k)++)
+
+#define K FOR_COORD(k)
+
+static double hypotD(const double p[], const double q[]) {
+  int k;
+  double pq[D3];
+  gsl_vector v;
+
+  K pq[_k]= p[k] - q[k];
+  v.size= D3;
+  v.stride= 1;
+  v.vector.data= pq;
+  /* owner and block ought not to be used */
+
+  return gsl_blas_snrm2(&v);
+}
+
+#ifdef FP_FAST_FMA
+# define ffsqa(factor,term) fma((factor),(factor),(term))
+#else
+# define ffsqu(factor,term) ((factor)*(factor)+(term))
+#endif
+
+static double energy_function(const double vertices[N][D3]) {
+  int pi,e,qi,ri,si, k;
+  double m[D3], mprime[D3], b, d2, l, sigma_bd2_l;
+
+  FOR_EDGE(pi,e,qi) {
+    ri= EDGE_END2(pi,(e+1)%V6); if (r<0) continue;
+    si= EDGE_END2(pi,(e+5)%V6); if (s<0) continue;
+    assert(ri == EDGE_END2(qi,(e+2)%V6));
+    assert(si == EDGE_END2(qi,(e+4)%V6));
+    
+    K m[k]= (vertices[pi][k] + vertices[qi][k]) * 0.5;
+    K mprime[k]= (vertices[ri][k] + vertices[si][k]) * 0.5;
+    b= hypotD(vertices[pi], vertices[qi]);
+    d2= 0; K d2= ffsqa(m[k] - mprime[k], d2);
+    l= hypotD(vertices[ri][k] - vertices[si][k]);
+
+    sigma_bd2_l += b * d2 / l;
+  }
+
+  
 
 static gsl_multimin_fminimizer *minimiser;
 
@@ -115,3 +231,4 @@ static gsl_multimin_fminimizer *minimiser;
 
 {
   gsl_multimin_fminimizer_alloc(
+