exploit symmetry

[moebius2.git] / minimise.c

diff --git a/minimise.c b/minimise.c
new file mode 100644 (file)
index 0000000..d7cc86d
--- /dev/null
+++ b/minimise.c
@@ -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);
+}