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   energy_init();
  77
  78   printf("X=%d=0x%x  Y=%d=0x%x  DIM=%d\n",X,X,Y,Y,DIM);
  79
  80   minimiser= gsl_multimin_fminimizer_alloc
  81     (gsl_multimin_fminimizer_nmsimplex, DIM);
  82   if (!minimiser) { perror("alloc minimiser"); exit(-1); }
  83
  84   multimin_function.f= minfunc_f;
  85   multimin_function.n= DIM;
  86   multimin_function.params= 0;
  87
  88   initial_f= fopen(input_file,"rb");  if (!initial_f) diee("fopen initial");
  89   errno= 0; r= fread(&initial_full,sizeof(initial_full),1,initial_f);
  90   if (r!=1) diee("fread");
  91   fclose(initial_f);
  92
  93   pmap_dimensions(&initial_full);
  94   unmap_dimensions(initial_half,&initial_full);
  95   for (i=0; i<DIM; i++) step_size[i]= 0.03;
  96
  97   initial_gsl.size= DIM;
  98   initial_gsl.stride= 1;
  99   initial_gsl.block= 0;
 100   initial_gsl.owner= 0;
 101   step_size_gsl= initial_gsl;
 102
 103   initial_gsl.data= initial_half;
 104   step_size_gsl.data= step_size;
 105
 106   GA( gsl_multimin_fminimizer_set(minimiser, &multimin_function,
 107                                   &initial_gsl, &step_size_gsl) );
 108
 109   for (;;) {
 110     GA( gsl_multimin_fminimizer_iterate(minimiser) );
 111
 112     size= gsl_multimin_fminimizer_size(minimiser);
 113     r= gsl_multimin_test_size(size, stop_epsilon);
 114
 115     if (printing_check(pr_size,155))
 116       printf("size %# e, r=%d\n", size, r);
 117     flushoutput();
 118
 119     if (r==GSL_SUCCESS) break;
 120     assert(r==GSL_CONTINUE);
 121   }
 122   return 0;
 123 }
 124
 125 /*---------- printing rate limit ----------*/
 126
 127 static volatile unsigned print_todo;
 128 static sigset_t print_alarmset;
 129
 130 int printing_check(enum printing_instance which, int indent) {
 131   static int skipped[pr__max];
 132
 133   unsigned bits= 1u << which;
 134   int sk;
 135
 136   if (!(print_todo & bits)) {
 137     skipped[which]++;
 138     return 0;;
 139   }
 140
 141   sigprocmask(SIG_BLOCK,&print_alarmset,0);
 142   print_todo &= ~bits;
 143   sigprocmask(SIG_UNBLOCK,&print_alarmset,0);
 144
 145   sk= skipped[which];
 146   if (sk) printf("%*s[%4d] ",indent,"", sk);
 147   else printf("       ");
 148   skipped[which]= 0;
 149
 150   return 1;
 151 }
 152
 153 static void alarmhandler(int ignored) {
 154   print_todo= ~0u;
 155 }
 156
 157 static void printing_init(void) {
 158   struct sigaction sa;
 159   struct itimerval itv;
 160
 161   sigemptyset(&print_alarmset);
 162   sigaddset(&print_alarmset,SIGALRM);
 163
 164   sa.sa_handler= alarmhandler;
 165   sa.sa_mask= print_alarmset;
 166   sa.sa_flags= SA_RESTART;
 167   if (sigaction(SIGALRM,&sa,0)) diee("sigaction ALRM");
 168
 169   itv.it_interval.tv_sec= 0;
 170   itv.it_interval.tv_usec= 200000;
 171   itv.it_value= itv.it_interval;
 172
 173   if (setitimer(ITIMER_REAL,&itv,0)) diee("setitimer REAL");
 174
 175   raise(SIGALRM);
 176 }