minimise.c

   1 /*
   2  * use of GSL
   3  */
   4
   5   /* We want to do multidimensional minimisation.
   6    *
   7    * We don't think there are any local minima.  Or at least, if there
   8    * are, the local minimum which will be found from the starting
   9    * state is the one we want.
  10    *
  11    * We don't want to try to provide a derivative of the cost
  12    * function.  That's too tedious (and anyway the polynomial
  13    * approximation to our our cost function sometimes has high degree
  14    * in the inputs which means the quadratic model implied by most of
  15    * the gradient descent minimisers is not ideal).
  16    *
  17    * This eliminates most of the algorithms.  Nelder and Mead's
  18    * simplex algorithm is still available and we will try that.
  19    *
  20    * In our application we are searching for the optimal locations of
  21    * N actualvertices in D3 (3) dimensions - ie, we are searching for
  22    * the optimal metapoint in an N*D3-dimensional space.
  23    *
  24    * So eg with X=Y=100, the simplex will contain 300 metavertices
  25    * each of which is an array of 300 doubles for the actualvertex
  26    * coordinates.  Hopefully this won't be too slow ...
  27    */
  28
  29 #include "common.h"
  30
  31 #include <gsl/gsl_errno.h>
  32 #include <gsl/gsl_multimin.h>
  33
  34 #include <signal.h>
  35 #include <sys/time.h>
  36
  37 #include "half.h"
  38 #include "minimise.h"
  39
  40 const char *input_file, *output_file;
  41 char *output_file_tmp;
  42
  43 static void printing_init(void);
  44
  45 static gsl_multimin_fminimizer *minimiser;
  46
  47 static const double stop_epsilon= 1e-6;
  48
  49 static double minfunc_f(const gsl_vector *x, void *params) {
  50   struct Vertices vs;
  51
  52   assert(x->size == DIM);
  53   assert(x->stride == 1);
  54   map_dimensions(x->data, vs.a);
  55   return compute_energy(&vs);
  56 }
  57
  58 int main(int argc, const char *const *argv) {
  59   gsl_multimin_function multimin_function;
  60   double size;
  61   struct Vertices initial_full;
  62   double initial_half[DIM], step_size[DIM];
  63   FILE *initial_f;
  64   gsl_vector initial_gsl, step_size_gsl;
  65   int r, i;
  66
  67   if (argc!=3 || argv[1][0]=='-' || strncmp(argv[2],"-o",2))
  68     { fputs("usage: minimise <input> -o<output\n",stderr); exit(8); }
  69
  70   input_file= argv[1];
  71   output_file= argv[2]+2;
  72   if (asprintf(&output_file_tmp,"%s.new",output_file) <= 0) diee("asprintf");
  73
  74   graph_layout_prepare();
  75   printing_init();
  76
  77   printf("X=%d=0x%x  Y=%d=0x%x  DIM=%d\n",X,X,Y,Y,DIM);
  78
  79   minimiser= gsl_multimin_fminimizer_alloc
  80     (gsl_multimin_fminimizer_nmsimplex, DIM);
  81   if (!minimiser) { perror("alloc minimiser"); exit(-1); }
  82
  83   multimin_function.f= minfunc_f;
  84   multimin_function.n= DIM;
  85   multimin_function.params= 0;
  86
  87   initial_f= fopen(input_file,"rb");  if (!initial_f) diee("fopen initial");
  88   errno= 0; r= fread(&initial_full,sizeof(initial_full),1,initial_f);
  89   if (r!=1) diee("fread");
  90   fclose(initial_f);
  91
  92   pmap_dimensions(&initial_full);
  93   unmap_dimensions(initial_half,&initial_full);
  94   for (i=0; i<DIM; i++) step_size[i]= 0.03;
  95
  96   initial_gsl.size= DIM;
  97   initial_gsl.stride= 1;
  98   initial_gsl.block= 0;
  99   initial_gsl.owner= 0;
 100   step_size_gsl= initial_gsl;
 101
 102   initial_gsl.data= initial_half;
 103   step_size_gsl.data= step_size;
 104
 105   GA( gsl_multimin_fminimizer_set(minimiser, &multimin_function,
 106                                   &initial_gsl, &step_size_gsl) );
 107
 108   for (;;) {
 109     GA( gsl_multimin_fminimizer_iterate(minimiser) );
 110
 111     size= gsl_multimin_fminimizer_size(minimiser);
 112     r= gsl_multimin_test_size(size, stop_epsilon);
 113
 114     if (printing_check(pr_size,155))
 115       printf("size %# e, r=%d\n", size, r);
 116     flushoutput();
 117
 118     if (r==GSL_SUCCESS) break;
 119     assert(r==GSL_CONTINUE);
 120   }
 121   return 0;
 122 }
 123
 124 /*---------- printing rate limit ----------*/
 125
 126 static volatile unsigned print_todo;
 127 static sigset_t print_alarmset;
 128
 129 int printing_check(enum printing_instance which, int indent) {
 130   static int skipped[pr__max];
 131
 132   unsigned bits= 1u << which;
 133   int sk;
 134
 135   if (!(print_todo & bits)) {
 136     skipped[which]++;
 137     return 0;;
 138   }
 139
 140   sigprocmask(SIG_BLOCK,&print_alarmset,0);
 141   print_todo &= ~bits;
 142   sigprocmask(SIG_UNBLOCK,&print_alarmset,0);
 143
 144   sk= skipped[which];
 145   if (sk) printf("%*s[%4d] ",indent,"", sk);
 146   else printf("       ");
 147   skipped[which]= 0;
 148
 149   return 1;
 150 }
 151
 152 static void alarmhandler(int ignored) {
 153   print_todo= ~0u;
 154 }
 155
 156 static void printing_init(void) {
 157   struct sigaction sa;
 158   struct itimerval itv;
 159
 160   sigemptyset(&print_alarmset);
 161   sigaddset(&print_alarmset,SIGALRM);
 162
 163   sa.sa_handler= alarmhandler;
 164   sa.sa_mask= print_alarmset;
 165   sa.sa_flags= SA_RESTART;
 166   if (sigaction(SIGALRM,&sa,0)) diee("sigaction ALRM");
 167
 168   itv.it_interval.tv_sec= 0;
 169   itv.it_interval.tv_usec= 200000;
 170   itv.it_value= itv.it_interval;
 171
 172   if (setitimer(ITIMER_REAL,&itv,0)) diee("setitimer REAL");
 173
 174   raise(SIGALRM);
 175 }