chiark / gitweb /
exploit symmetry
authorIan Jackson <ian@davenant.greenend.org.uk>
Sun, 6 Jan 2008 18:11:21 +0000 (18:11 +0000)
committerIan Jackson <ian@davenant.greenend.org.uk>
Sun, 6 Jan 2008 18:11:21 +0000 (18:11 +0000)
Makefile
energy.c
half.c [new file with mode: 0644]
half.h [new file with mode: 0644]
mgraph.h
minimise.c [new file with mode: 0644]
minimise.h
primer.c

index e1a54dd723204884cda7c5f36ae03ac2d2747156..15e4ae0f1eb983d28c66d84fe104baf35b875ca8 100644 (file)
--- a/Makefile
+++ b/Makefile
@@ -14,7 +14,7 @@ o= >$@.new && mv -f $@.new $@
 
 all:           $(TARGETS)
 
-minimise:      energy.o graph.o common.o mgraph.o
+minimise:      energy.o graph.o common.o mgraph.o minimise.o half.o
                $(CXX) $(CXXFLAGS) -o $@ $^ $(LIBGSL)
 
 primer:                primer.o common.o
index 78b8f31ff44be4942efb40cb91f61d21718217b2..ed50bc7e0effa6d235d6c29d02c22aab7feee5b5 100644 (file)
--- a/energy.c
+++ b/energy.c
@@ -6,41 +6,28 @@
 #include "minimise.h"
 #include "mgraph.h"
 
-#include <gsl/gsl_errno.h>
-#include <gsl/gsl_multimin.h>
-
-#include <signal.h>
-#include <sys/time.h>
-
-static const char *input_file, *output_file;
-static char *output_file_tmp;
-
 static void compute_vertex_areas(const Vertices vertices, double areas[N]);
 static double best_energy= DBL_MAX;
 
-enum printing_instance { pr_cost, pr_size, pr__max };
-
 static void addcost(double *energy, double tweight, double tcost, int pr);
 #define COST(weight, compute) addcost(&energy, (weight), (compute), printing)
-static int printing_check(enum printing_instance);
-static void printing_init(void);
 
 /*---------- main energy computation and subroutines ----------*/
 
-static double compute_energy(const Vertices vertices) {
+double compute_energy(const struct Vertices *vs) {
   double vertex_areas[N], energy;
   int printing;
 
-  compute_vertex_areas(vertices,vertex_areas);
+  compute_vertex_areas(vs->a, vertex_areas);
   energy= 0;
 
   printing= printing_check(pr_cost);
 
   if (printing) printf("cost > energy |");
 
-  COST(1e2, edgewise_vertex_displacement_cost(vertices));
-  COST(1e2, graph_layout_cost(vertices,vertex_areas));
-//  COST(1e4, noncircular_rim_cost(vertices));
+  COST(1e2, edgewise_vertex_displacement_cost(vs->a));
+  COST(1e2, graph_layout_cost(vs->a,vertex_areas));
+//  COST(1e4, noncircular_rim_cost(vs->a));
 
   if (printing) printf("| total %# e |", energy);
 
@@ -51,7 +38,7 @@ static double compute_energy(const Vertices vertices) {
     if (printing) printf(" BEST");
 
     best_f= fopen(output_file_tmp,"wb");  if (!best_f) diee("fopen new out");
-    r= fwrite(vertices,sizeof(Vertices),1,best_f); if (r!=1) diee("fwrite");
+    r= fwrite(vs->a,sizeof(vs->a),1,best_f); if (r!=1) diee("fwrite");
     if (fclose(best_f)) diee("fclose new best");
     if (rename(output_file_tmp,output_file)) diee("rename install new best");
 
@@ -96,107 +83,6 @@ static void compute_vertex_areas(const Vertices vertices, double areas[N]) {
   }
 }
 
-/*---------- use of GSL ----------*/
-
-  /* We want to do multidimensional minimisation.
-   *
-   * We don't think there are any local minima.  Or at least, if there
-   * are, the local minimum which will be found from the starting
-   * state is the one we want.
-   *
-   * We don't want to try to provide a derivative of the cost
-   * function.  That's too tedious (and anyway the polynomial
-   * approximation to our our cost function sometimes has high degree
-   * in the inputs which means the quadratic model implied by most of
-   * the gradient descent minimisers is not ideal).
-   *
-   * This eliminates most of the algorithms.  Nelder and Mead's
-   * simplex algorithm is still available and we will try that.
-   *
-   * In our application we are searching for the optimal locations of
-   * N actualvertices in D3 (3) dimensions - ie, we are searching for
-   * the optimal metapoint in an N*D3-dimensional space.
-   *
-   * So eg with X=Y=100, the simplex will contain 300 metavertices
-   * each of which is an array of 300 doubles for the actualvertex
-   * coordinates.  Hopefully this won't be too slow ...
-   */
-
-static gsl_multimin_fminimizer *minimiser;
-
-static const double stop_epsilon= 1e-6;
-
-static double minfunc_f(const gsl_vector *x, void *params) {
-  assert(x->size == DIM);
-  assert(x->stride == 1);
-  return compute_energy((const double(*)[D3])x->data);
-}
-
-int main(int argc, const char *const *argv) {
-  gsl_multimin_function multimin_function;
-  double size;
-  Vertices initial, step_size;
-  FILE *initial_f;
-  gsl_vector initial_gsl, step_size_gsl;
-  int r, v, k;
-
-  if (argc!=3 || argv[1][0]=='-' || strncmp(argv[2],"-o",2))
-    { fputs("usage: minimise <input> -o<output\n",stderr); exit(8); }
-
-  input_file= argv[1];
-  output_file= argv[2]+2;
-  if (asprintf(&output_file_tmp,"%s.new",output_file) <= 0) diee("asprintf");
-
-  graph_layout_prepare();
-  printing_init();
-  
-  minimiser= gsl_multimin_fminimizer_alloc
-    (gsl_multimin_fminimizer_nmsimplex, DIM);
-  if (!minimiser) { perror("alloc minimiser"); exit(-1); }
-
-  multimin_function.f= minfunc_f;
-  multimin_function.n= DIM;
-  multimin_function.params= 0;
-
-  initial_f= fopen(input_file,"rb");  if (!initial_f) diee("fopen initial");
-  errno= 0; r= fread(initial,sizeof(initial),1,initial_f);
-  if (r!=1) diee("fread");
-  fclose(initial_f);
-
-  initial_gsl.size= DIM;
-  initial_gsl.stride= 1;
-  initial_gsl.block= 0;
-  initial_gsl.owner= 0;
-  step_size_gsl= initial_gsl;
-
-  initial_gsl.data= &initial[0][0];
-  step_size_gsl.data= &step_size[0][0];
-
-  FOR_VERTEX(v)
-    K step_size[v][k]= 0.03;
-//int vx,vy;
-//  FOR_RIM_VERTEX(vx,vy,v)
-//    step_size[v][3] *= 0.1;
-
-  GA( gsl_multimin_fminimizer_set(minimiser, &multimin_function,
-                                 &initial_gsl, &step_size_gsl) );
-
-  for (;;) {
-    GA( gsl_multimin_fminimizer_iterate(minimiser) );
-
-    size= gsl_multimin_fminimizer_size(minimiser);
-    r= gsl_multimin_test_size(size, stop_epsilon);
-
-    if (printing_check(pr_size))
-      printf("%*s size %# e, r=%d\n", 135,"", size, r);
-    flushoutput();
-
-    if (r==GSL_SUCCESS) break;
-    assert(r==GSL_CONTINUE);
-  }
-  return 0;
-}
-
 /*---------- Edgewise vertex displacement ----------*/
 
   /*
@@ -286,56 +172,3 @@ double noncircular_rim_cost(const Vertices vertices) {
   }
   return cost;
 }
-
-/*---------- printing rate limit ----------*/
-
-static volatile unsigned print_todo;
-static sigset_t print_alarmset;
-
-static int printing_check(enum printing_instance which) {
-  static int skipped[pr__max];
-  
-  unsigned bits= 1u << which;
-  int sk;
-
-  if (!(print_todo & bits)) {
-    skipped[which]++;
-    return 0;;
-  }
-
-  sigprocmask(SIG_BLOCK,&print_alarmset,0);
-  print_todo &= ~bits;
-  sigprocmask(SIG_UNBLOCK,&print_alarmset,0);
-
-  sk= skipped[which];
-  if (sk) printf("[%4d] ",sk);
-  else printf("       ");
-  skipped[which]= 0;
-
-  return 1;
-}
-
-static void alarmhandler(int ignored) {
-  print_todo= ~0u;
-}
-
-static void printing_init(void) {
-  struct sigaction sa;
-  struct itimerval itv;
-
-  sigemptyset(&print_alarmset);
-  sigaddset(&print_alarmset,SIGALRM);
-
-  sa.sa_handler= alarmhandler;
-  sa.sa_mask= print_alarmset;
-  sa.sa_flags= SA_RESTART;
-  if (sigaction(SIGALRM,&sa,0)) diee("sigaction ALRM");
-  
-  itv.it_interval.tv_sec= 0;
-  itv.it_interval.tv_usec= 200000;
-  itv.it_value= itv.it_interval;
-
-  if (setitimer(ITIMER_REAL,&itv,0)) diee("setitimer REAL");
-
-  raise(SIGALRM);
-}
diff --git a/half.c b/half.c
new file mode 100644 (file)
index 0000000..5ac69f8
--- /dev/null
+++ b/half.c
@@ -0,0 +1,157 @@
+/*
+ * actually do mapping from metaspace to realspace
+ */
+
+#include "half.h"
+
+/*
+ * Our metaspace dimensions are, in order:
+ *   - y realspace coordinate of symmetry axis vertices
+ *   - x,y,z realspace coordinates of even-y opposite-of-axis vertices
+ *   - x,y,z realspace coordinates of even-y non-axis vertices
+ *   - x,y,z realspace coordinates of odd-y vertices in canonical order
+ */
+
+#define FOR_MAP_DIMENSION__VERTEX_PAIR(body,mirror)    \
+  k=0; sign=-1; { body } { mirror } i++;               \
+  k=1; sign=+1; { body } { mirror } i++;               \
+  k=2; sign=-1; { body } { mirror } i++;
+
+/*#define FOR_MAP_DIMENSION_VPRIME             \
+  int vprime= N-1-v
+  */
+
+#define FOR_MAP_DIMENSION(body, mirror, degenerate)            \
+  do{                                                          \
+    int i,v,vprime,k,sign; /* public */                                \
+    int fmd_i, fmd_x, fmd_y;                                   \
+    i= 0;                                                      \
+    for (fmd_i=0; fmd_i<HALF_DEGRADED_NVERTICES; fmd_i++) {    \
+      v= fmd_i*X*2;                                            \
+      k=0;          { degenerate }                             \
+      k=1; sign=+1; { body } i++;                              \
+      k=2;          { degenerate }                             \
+    }                                                          \
+    for (fmd_y=0; fmd_y<Y; fmd_y+=2) {                         \
+      v= fmd_y << YSHIFT | (X/2);                              \
+      FOR_MAP_DIMENSION__VERTEX_PAIR(body,)                    \
+    }                                                          \
+    for (fmd_y=0; fmd_y<Y; fmd_y+=2)                           \
+      for (fmd_x=1; fmd_x<X/2; fmd_x++) {                      \
+       v= fmd_y << YSHIFT | fmd_x;                             \
+        vprime= ((Y-1) << YSHIFT | X) - v;                     \
+        FOR_MAP_DIMENSION__VERTEX_PAIR(body,mirror)            \
+      }                                                                \
+    for (fmd_y=1; fmd_y<Y; fmd_y+=2)                           \
+      for (fmd_x=0; fmd_x<X/2; fmd_x++)        {                       \
+       v= fmd_y << YSHIFT | fmd_x;                             \
+        vprime= ((Y-1) << YSHIFT | (X-1)) - v;                 \
+        FOR_MAP_DIMENSION__VERTEX_PAIR(body,mirror)            \
+      }                                                                \
+    assert(i == DIM);                                          \
+  }while(0)
+
+void pmap_dimensions(const struct Vertices *vs) {
+  int ixfrom[N][D3], count[N][D3], bad;
+  int vc,kc;
+
+  bad= 0;
+
+  FOR_VERTEX(vc)
+    FOR_COORD(kc) {
+      count[vc][kc]= 0;
+      ixfrom[vc][kc]= -1;
+    }
+  
+#define CHECK_V(v) assert(v >= 0); assert(v < N);
+#define CHECK_NEAR(a,b)                                        \
+  if (abs(a-b) >= 1e-6) { printf("    BAD!"); bad++; } \
+  putchar('\n');
+
+  FOR_MAP_DIMENSION(
+                   (printf("i=%2d  v=%03x k=%d  %# 20g\n",
+                           i,v,k,
+                           vs->a[v][k]
+                           ));
+                   CHECK_V(v)
+                   ixfrom[v][k]=i;
+                   count[v][k]++;,
+
+//                 FOR_MAP_DIMENSION_VPRIME;
+                   (printf("i=%2d v'=%03x k=%d  %# 20g := %# 20g",
+                           i,vprime,k,
+                           vs->a[vprime][k], sign*vs->a[v][k]
+                           ));
+                   CHECK_V(vprime)
+                   CHECK_NEAR(vs->a[vprime][k], sign*vs->a[v][k])
+                   ixfrom[vprime][k]=-i;
+                   count[vprime][k]++;,
+
+                   (printf("zero  v=%03x k=%d  %# 20g := 0",
+                           v,k,
+                           vs->a[v][k]
+                           ));
+                   CHECK_V(v)
+                   CHECK_NEAR(vs->a[v][k], 0)
+                   count[v][k]++;
+                   );
+
+  FOR_VERTEX(vc)
+    FOR_COORD(kc) {
+      if (count[vc][kc]==1) continue;
+      printf("BAD %03x %d count=%d ixfrom=%d\n",
+            vc,kc, count[vc][kc], ixfrom[vc][kc]);
+      bad++;
+    }
+  assert(!bad);
+}
+
+void map_dimensions(const double metavertex[DIM], Vertices out) {
+  FOR_MAP_DIMENSION(out[v     ][k]=      metavertex[i];,
+
+//                 FOR_MAP_DIMENSION_VPRIME;
+                   out[vprime][k]= sign*metavertex[i];,
+
+                   out[v     ][k]=      0;
+                   );
+}
+
+void unmap_dimensions(double metavertex[DIM], const struct Vertices *in) {
+  FOR_MAP_DIMENSION(
+                   metavertex[i]= in->a[v][k];,
+                   ,
+                   );
+}
+
+
+/*
+ *  OLD NUMBERING
+ *
+ * Our numbering is as follows:
+ *
+ *          :axis of symmetry
+ *               | :
+ *               | :
+ *     ___| 0  ___  1  ___  2  ___  3  ___  4  ___    up to X"-1
+ *        | :                                           where X"=(X/2)
+ *        / :\    /  \    /  \    /  \    /  \
+ *       /| : \  /    \  /    \  /    \  /    \
+ *        |___ 0   ___ 1   ___ 2   ___ 3   ___ 4
+ *        | :  +X"     +X"     +X"     +X"     +X"
+ *       \| : /  \    /  \    /  \    /  \    /
+ *        \ :/    \  /    \  /    \  /    \  /
+ *     ___|2*X"___ 1+  ___ 2+  ___ 3+  __  4+  ___
+ *        | :     2*X"   2*X"    2*X"    2*X"
+ *        / :\    /  \    /  \    /  \    /  \
+ *       /| : \  /    \  /    \  /    \  /    \
+ *        |___ 0+  ___ 1+  ___ 2+  ___ 3+  ___ 4+
+ *        | : 3*X"    3*X"    3*X"    3*X"    3*X"
+ *        | :
+ *
+ * X" is still a power of two so a vertex pair can be named by a
+ * number structured like this:
+ *
+ * vertex number":   0000 | y     | x
+ *                         YBITS   XBITS-1
+ *
+ */
diff --git a/half.h b/half.h
new file mode 100644 (file)
index 0000000..88c2e1a
--- /dev/null
+++ b/half.h
@@ -0,0 +1,38 @@
+/*
+ * During the actual optimisation we constrain it to be symmetrical.
+ * That means we have half as many vertices.  We process the vertices
+ * with X >=0.
+ *
+ * We conventionally suffix all of our variables, macros, constants,
+ * etc., with `h', to represent the double-prime ".
+ *
+ * Vertices with x==0 are special as they lie on the symmetry axis:
+ * all but one of the coordinates are fixed.  We do not offer those
+ * coordinates to our optimisation algorithm.
+ *
+ * We provide to the rest of the program just:
+ *    NH                   number of vertices"
+ *    VH2V(v)              maps a v" to the representative vertex number v
+ *    FOR_VERTEX_H(vh,v)   iterates over v", computing each v
+ *    DIM                  dimensionality of coordinate space
+ *    [un]map_dimensions   maps double[DIM] to Vertices and back
+ */
+
+#ifndef HALF_H
+#define HALF_H
+
+#include "mgraph.h"
+
+#define NH (N/2)
+#define VH2V(v) ((v<<1) & YMASK | v & (XMASK>>1))
+#define FOR_VERTEX_H(vh,v) \
+  for ((vh)=0; (v)= VH2V((vh)), (vh) < NH; (vh)++)
+
+#define HALF_DEGRADED_NVERTICES ((Y+1)/2)
+#define DIM (NH*D3 + HALF_DEGRADED_NVERTICES*(1-D3) + ((Y+1)/2)*D3)
+
+void unmap_dimensions(double metavertex[DIM], const struct Vertices *in);
+void map_dimensions(const double metavertex[DIM], Vertices out);
+void pmap_dimensions(const struct Vertices *vs);
+
+#endif /*HALF_H*/
index eb4bd129fb990bb650d1da0605e780e1a13680af..8636b91bff06efdcbdfcf275123874f514eb9402 100644 (file)
--- a/mgraph.h
+++ b/mgraph.h
@@ -4,34 +4,38 @@
 /*
  * Vertices in strip are numbered as follows:
  *
- *                       |
+ *                          :axis of symmetry
+ *                       | :
+ *                       | :
  *     ___ X-2 ___ X-1 ___| 0  ___  1  ___  2  ___  3  ___  4  __
- *         Y-1     Y-1    |0       0       0       0       0
- *        /  \    /  \    /  \    /  \    /  \    /  \    /  \
- *       /    \  /    \  /|   \  /    \  /    \  /    \  /    \
+ *         Y-1     Y-1    |0:      0       0       0       0
+ *        /  \    /  \    / :\    /  \    /  \    /  \    / !! \
+ *       /    \  /    \  /| : \  /    \  /    \  /    \  /    \
  *     X-3 ___ X-2 ___ X-1|___  0  ___  1  ___  2  ___  3  ___  4
- *     Y-2     Y-2     Y-2|    1       1       1       1       1
- *       \    /  \    /  \|   /  \    /  \    /  \    /  \    /
- *        \  /    \  /    \  /    \  /    \  /    \  /    \  /
+ *     Y-2     Y-2     Y-2| :  1       1       1       1       1
+ *       \    /  \    /  \| : /  \    /  \    /  \    /  \    /
+ *        \  /    \  /    \ :/    \  /    \  /    \  /    \  /
  *     ___ X-2 ___ X-1 ___| 0  ___  1  ___  2  ___  3  __   4 ___
- *         Y-3     Y-3    |2       2       2       2      2
- *        /  \    /  \    /  \    /  \    /  \    /  \    /  \
- *       /    \  /    \  /|   \  /    \  /    \  /    \  /    \
+ *         Y-3     Y-3    |2:      2       2       2      2
+ *        /  \    /  \    / :\    /  \    /  \    /  \    /  \
+ *       /    \  /    \  /| : \  /    \  /    \  /    \  /    \
  *     X-3 ___ X-2 ___ X-1|___  0  ___  1  ___  2  ___  3  ___  4
- *     Y-4     Y-4     Y-4|    3       3       3       3       3
- *                       |
- *       .   .   .   .   .|  .   .   .   .   .   .   .   .   .   .
- *                       |
+ *     Y-4     Y-4     Y-4| :  3       3       3       3       3
+ *                       | :
+ *       .   .   .   .   .| :.   .   .   .   .   .   .   .   .   .
+ *                       | :
  *     ___ X-2 ___ X-1 ___| 0  ___  1  ___  2  ___  3  ___  4 ___
  *         2       2      |Y-3     Y-3     Y-3     Y-3     Y-3
- *        /  \    /  \    /  \    /  \    /  \    /  \    /  \
- *       /    \  /    \  /|   \  /    \  /    \  /    \  /    \
+ *        /  \    /  \    / :\    /  \    /  \    /  \    /  \
+ *       /    \  /    \  /| : \  /    \  /    \  /    \  /    \
  *          __ X-2 ___ X-1|___  0  ___  1  ___  2  ___  3  ___  3  ___  4
- *             1       1  |    Y-2     Y-2     Y-2     Y-2     Y-2     Y-2
- *    /  \    /  \    /  \|   /  \    /  \    /         \    /  \    /
- *   /    \  /    \  /    \  /    \  /    \  /   \  /    \  /
+ *             1       1  | :  Y-2     Y-2     Y-2     Y-2     Y-2     Y-2
+ *    /  \    /  \    /  \| : /  \    /  \    /         \    /  \    /
+ *   /    \  /    \  /    \ :/    \  /    \  /   \  /    \  /
  *  -3 ___ X-2 ___ X-1 ___| 0  ___  1  ___  2  ___  3  ___  4 ___
  *  0      0       0      |Y-1     Y-1     Y-1     Y-1     Y-1
+ *                       | :
+ *                       | :
  *                       |
  *                        ^ join, where there is
  *                           a discontinuity in numbering
@@ -55,6 +59,9 @@
  *               3    1   0
  *                  /  \
  *                4/   5\
+ *
+ * vertex number:   0000 | y     | x
+ *                        YBITS   XBITS
  */
 
 #ifndef MGRAPH_H
 
 #include "common.h"
 
-#define XBITS 3
+#define XBITS 4
 #define X (1<<XBITS)
-#define YBITS 3
+#define YBITS 4
 #define Y ((1<<YBITS) - 1)
 
-/* 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 << YSHIFT)
 
-#define DIM (N*D3)
-
 #define V6 6
 
 #define FOR_VERTEX(v) \
@@ -112,5 +113,6 @@ extern int edge_end2(unsigned v1, int e);
     for ((vx)=0; (v)= (vy)<<YSHIFT | (vx), (vx)<X; (vx)++)
 
 typedef double Vertices[N][D3];
+struct Vertices { Vertices a; };
 
 #endif /*MGRAPH_H*/
diff --git a/minimise.c b/minimise.c
new file mode 100644 (file)
index 0000000..d7cc86d
--- /dev/null
@@ -0,0 +1,175 @@
+/*
+ * use of GSL
+ */
+
+  /* We want to do multidimensional minimisation.
+   *
+   * We don't think there are any local minima.  Or at least, if there
+   * are, the local minimum which will be found from the starting
+   * state is the one we want.
+   *
+   * We don't want to try to provide a derivative of the cost
+   * function.  That's too tedious (and anyway the polynomial
+   * approximation to our our cost function sometimes has high degree
+   * in the inputs which means the quadratic model implied by most of
+   * the gradient descent minimisers is not ideal).
+   *
+   * This eliminates most of the algorithms.  Nelder and Mead's
+   * simplex algorithm is still available and we will try that.
+   *
+   * In our application we are searching for the optimal locations of
+   * N actualvertices in D3 (3) dimensions - ie, we are searching for
+   * the optimal metapoint in an N*D3-dimensional space.
+   *
+   * So eg with X=Y=100, the simplex will contain 300 metavertices
+   * each of which is an array of 300 doubles for the actualvertex
+   * coordinates.  Hopefully this won't be too slow ...
+   */
+
+#include "common.h"
+
+#include <gsl/gsl_errno.h>
+#include <gsl/gsl_multimin.h>
+
+#include <signal.h>
+#include <sys/time.h>
+
+#include "half.h"
+#include "minimise.h"
+
+const char *input_file, *output_file;
+char *output_file_tmp;
+
+static void printing_init(void);
+
+static gsl_multimin_fminimizer *minimiser;
+
+static const double stop_epsilon= 1e-6;
+
+static double minfunc_f(const gsl_vector *x, void *params) {
+  struct Vertices vs;
+  
+  assert(x->size == DIM);
+  assert(x->stride == 1);
+  map_dimensions(x->data, vs.a);
+  return compute_energy(&vs);
+}
+
+int main(int argc, const char *const *argv) {
+  gsl_multimin_function multimin_function;
+  double size;
+  struct Vertices initial_full;
+  double initial_half[DIM], step_size[DIM];
+  FILE *initial_f;
+  gsl_vector initial_gsl, step_size_gsl;
+  int r, i;
+
+  if (argc!=3 || argv[1][0]=='-' || strncmp(argv[2],"-o",2))
+    { fputs("usage: minimise <input> -o<output\n",stderr); exit(8); }
+
+  input_file= argv[1];
+  output_file= argv[2]+2;
+  if (asprintf(&output_file_tmp,"%s.new",output_file) <= 0) diee("asprintf");
+
+  graph_layout_prepare();
+  printing_init();
+
+  printf("X=%d=0x%x  Y=%d=0x%x  DIM=%d\n",X,X,Y,Y,DIM);
+  
+  minimiser= gsl_multimin_fminimizer_alloc
+    (gsl_multimin_fminimizer_nmsimplex, DIM);
+  if (!minimiser) { perror("alloc minimiser"); exit(-1); }
+
+  multimin_function.f= minfunc_f;
+  multimin_function.n= DIM;
+  multimin_function.params= 0;
+
+  initial_f= fopen(input_file,"rb");  if (!initial_f) diee("fopen initial");
+  errno= 0; r= fread(&initial_full,sizeof(initial_full),1,initial_f);
+  if (r!=1) diee("fread");
+  fclose(initial_f);
+
+  pmap_dimensions(&initial_full);
+  unmap_dimensions(initial_half,&initial_full);
+  for (i=0; i<DIM; i++) step_size[i]= 0.03;
+
+  initial_gsl.size= DIM;
+  initial_gsl.stride= 1;
+  initial_gsl.block= 0;
+  initial_gsl.owner= 0;
+  step_size_gsl= initial_gsl;
+
+  initial_gsl.data= initial_half;
+  step_size_gsl.data= step_size;
+
+  GA( gsl_multimin_fminimizer_set(minimiser, &multimin_function,
+                                 &initial_gsl, &step_size_gsl) );
+
+  for (;;) {
+    GA( gsl_multimin_fminimizer_iterate(minimiser) );
+
+    size= gsl_multimin_fminimizer_size(minimiser);
+    r= gsl_multimin_test_size(size, stop_epsilon);
+
+    if (printing_check(pr_size))
+      printf("%*s size %# e, r=%d\n", 135,"", size, r);
+    flushoutput();
+
+    if (r==GSL_SUCCESS) break;
+    assert(r==GSL_CONTINUE);
+  }
+  return 0;
+}
+
+/*---------- printing rate limit ----------*/
+
+static volatile unsigned print_todo;
+static sigset_t print_alarmset;
+
+int printing_check(enum printing_instance which) {
+  static int skipped[pr__max];
+  
+  unsigned bits= 1u << which;
+  int sk;
+
+  if (!(print_todo & bits)) {
+    skipped[which]++;
+    return 0;;
+  }
+
+  sigprocmask(SIG_BLOCK,&print_alarmset,0);
+  print_todo &= ~bits;
+  sigprocmask(SIG_UNBLOCK,&print_alarmset,0);
+
+  sk= skipped[which];
+  if (sk) printf("[%4d] ",sk);
+  else printf("       ");
+  skipped[which]= 0;
+
+  return 1;
+}
+
+static void alarmhandler(int ignored) {
+  print_todo= ~0u;
+}
+
+static void printing_init(void) {
+  struct sigaction sa;
+  struct itimerval itv;
+
+  sigemptyset(&print_alarmset);
+  sigaddset(&print_alarmset,SIGALRM);
+
+  sa.sa_handler= alarmhandler;
+  sa.sa_mask= print_alarmset;
+  sa.sa_flags= SA_RESTART;
+  if (sigaction(SIGALRM,&sa,0)) diee("sigaction ALRM");
+  
+  itv.it_interval.tv_sec= 0;
+  itv.it_interval.tv_usec= 200000;
+  itv.it_value= itv.it_interval;
+
+  if (setitimer(ITIMER_REAL,&itv,0)) diee("setitimer REAL");
+
+  raise(SIGALRM);
+}
index 85dc96f31df978b0b0408bf556a292a0231a55d3..74788d3fbe862b55e3ba4459adbf98dd58dd6685 100644 (file)
@@ -7,10 +7,18 @@
 
 #include "mgraph.h"
 
+double compute_energy(const struct Vertices *vs);
+
 double graph_layout_cost(const Vertices v, const double vertex_areas[N]);
 void graph_layout_prepare();
 
 double noncircular_rim_cost(const Vertices vertices);
 double edgewise_vertex_displacement_cost(const Vertices vertices);
 
+extern const char *input_file, *output_file;
+extern char *output_file_tmp;
+
+enum printing_instance { pr_cost, pr_size, pr__max };
+int printing_check(enum printing_instance);
+
 #endif /*MINIMISE_H*/
index 1f87036119b45ac54d2fc254276031e727290b96..70cbcfeda45883205414c05c6b285df43fbb790c 100644 (file)
--- a/primer.c
+++ b/primer.c
@@ -14,7 +14,7 @@ int main(int argc, const char **argv) {
   if (argc!=1) { fputs("need no args\n",stderr); exit(8); }
 
   printf("%d %d %d %d %d\n%%-%d.%dg\n",
-        DIM, N, X, Y, D3,
+        X*Y, N, X, Y, D3,
         prec+5,prec);
 
   FOR_VERTEX(vi) {
@@ -26,7 +26,9 @@ int main(int argc, const char **argv) {
      * So that corresponds to 0..X (since 0==X in our scheme).
      * Vertices with odd y coordinate are halfway to the next x coordinate.
      */
-    double v= (x*2 + (y&1)) * M_PI / (X*2);
+    double v= (x*2 + (y&1)) * 1.0 / (X*2);
+    v += 0.5;
+    v *= M_PI;
 
     printf("%-*.*g %-*.*g # %03x %2d %2d\n",
           prec+5,prec,u, prec+5,prec,v,