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, *best_file;
  41 char *best_file_tmp;
  42 long long evaluations;
  43 double stop_epsilon= 1e-20;
  44
  45 static void printing_init(void);
  46
  47 static gsl_multimin_fminimizer *minimiser;
  48 static const char *final_file;
  49 static char *final_file_tmp;
  50
  51 static double minfunc_f(const gsl_vector *x, void *params) {
  52   struct Vertices vs;
  53
  54   assert(x->size == DIM);
  55   assert(x->stride == 1);
  56   map_dimensions(x->data, vs.a);
  57   return compute_energy(&vs);
  58 }
  59
  60 static void badusage(void) {
  61   fputs("usage: minimise <input> [-i<intermediate>] -o<output\n",stderr);
  62   exit(8);
  63 }
  64
  65 static sig_atomic_t quit_requested;
  66 sig_atomic_t quitting_reported_threads;
  67 int quitting_last_iteration;
  68
  69 static void sigint_handler(int ignored) {
  70   quit_requested= 1;
  71 }
  72
  73 int main(int argc, const char *const *argv) {
  74   gsl_multimin_function multimin_function;
  75   double size;
  76   struct Vertices initial_full;
  77   double initial_half[DIM], step_size[DIM];
  78   FILE *initial_f;
  79   gsl_vector initial_gsl, step_size_gsl;
  80   int r, i, size_check_counter= 0;
  81   struct sigaction sa;
  82
  83   if (argc==3) {
  84     input_file= argv[1];
  85     if (strncmp(argv[2],"-o",2)) badusage();   best_file= argv[2]+2;
  86     final_file= 0;
  87   } else if (argc==4) {
  88     input_file= argv[1];
  89     if (strncmp(argv[2],"-i",2)) badusage();   best_file= argv[2]+2;
  90     if (strncmp(argv[3],"-o",2)) badusage();   final_file= argv[3]+2;
  91   } else {
  92     badusage();
  93   }
  94   if (argv[1][0]=='-') badusage();
  95
  96   if (asprintf(&best_file_tmp,"%s.new",best_file) <= 0) diee("asprintf");
  97   if (final_file)
  98     if (asprintf(&final_file_tmp,"%s.new",final_file) <= 0) diee("asprintf");
  99
 100   memset(&sa,0,sizeof(sa));
 101   sa.sa_handler= sigint_handler;
 102   sa.sa_flags= SA_RESETHAND;
 103   r= sigaction(SIGINT,&sa,0);
 104   if (r) diee("sigaction SIGINT");
 105
 106   graph_layout_prepare();
 107   printing_init();
 108   energy_init();
 109
 110   printf("X=%d=0x%x  Y=%d=0x%x  DIM=%d\n",X,X,Y,Y,DIM);
 111
 112   minimiser= gsl_multimin_fminimizer_alloc
 113     (gsl_multimin_fminimizer_nmsimplex, DIM);
 114   if (!minimiser) { perror("alloc minimiser"); exit(-1); }
 115
 116   multimin_function.f= minfunc_f;
 117   multimin_function.n= DIM;
 118   multimin_function.params= 0;
 119
 120   initial_f= fopen(input_file,"rb");  if (!initial_f) diee("fopen initial");
 121   errno= 0; r= fread(&initial_full,sizeof(initial_full),1,initial_f);
 122   if (r!=1) diee("fread");
 123   fclose(initial_f);
 124
 125   pmap_dimensions(&initial_full);
 126   unmap_dimensions(initial_half,&initial_full);
 127   for (i=0; i<DIM; i++) step_size[i]= 0.03;
 128
 129   initial_gsl.size= DIM;
 130   initial_gsl.stride= 1;
 131   initial_gsl.block= 0;
 132   initial_gsl.owner= 0;
 133   step_size_gsl= initial_gsl;
 134
 135   initial_gsl.data= initial_half;
 136   step_size_gsl.data= step_size;
 137
 138   GA( gsl_multimin_fminimizer_set(minimiser, &multimin_function,
 139                                   &initial_gsl, &step_size_gsl) );
 140
 141   for (;;) {
 142     if (quit_requested) {
 143       fprintf(stderr,"SIGINT caught, quitting soon.\n");
 144       quitting_last_iteration= 1;
 145     }
 146
 147     GA( gsl_multimin_fminimizer_iterate(minimiser) );
 148
 149     if (!(size_check_counter++ % DIM)) {
 150       size= gsl_multimin_fminimizer_size(minimiser);
 151       r= gsl_multimin_test_size(size, stop_epsilon);
 152
 153       if (printing_check(pr_size,215))
 154         printf("r=%2d, size %# e\n", r, size);
 155       flushoutput();
 156
 157       if (r==GSL_SUCCESS) break;
 158       assert(r==GSL_CONTINUE);
 159     }
 160
 161     if (quitting_reported_threads)
 162       exit(1);
 163   }
 164
 165   if (final_file) {
 166     if (unlink(final_file_tmp) && errno != ENOENT) diee("unlink final .tmp");
 167     if (link(best_file,final_file_tmp)) diee("link final .tmp");
 168     if (rename(final_file_tmp,final_file)) diee("install final");
 169   }
 170   printf("FINISHED (%lld evaluations)\n",evaluations);
 171
 172   return 0;
 173 }
 174
 175 /*---------- printing rate limit ----------*/
 176
 177 static volatile unsigned print_todo;
 178 static sigset_t print_alarmset;
 179
 180 int printing_check(enum printing_instance which, int indent) {
 181   static int skipped[pr__max];
 182
 183   unsigned bits= 1u << which;
 184   int sk;
 185
 186   if (!(print_todo & bits)) {
 187     skipped[which]++;
 188     return 0;;
 189   }
 190
 191   sigprocmask(SIG_BLOCK,&print_alarmset,0);
 192   print_todo &= ~bits;
 193   sigprocmask(SIG_UNBLOCK,&print_alarmset,0);
 194
 195   printf("%*s",indent,"");
 196   sk= skipped[which];
 197   if (sk) printf("[%4d] ",sk);
 198   else printf("      ");
 199   skipped[which]= 0;
 200
 201   return 1;
 202 }
 203
 204 static void alarmhandler(int ignored) {
 205   print_todo= ~0u;
 206 }
 207
 208 static void printing_init(void) {
 209   struct sigaction sa;
 210   struct itimerval itv;
 211
 212   sigemptyset(&print_alarmset);
 213   sigaddset(&print_alarmset,SIGALRM);
 214
 215   sa.sa_handler= alarmhandler;
 216   sa.sa_mask= print_alarmset;
 217   sa.sa_flags= SA_RESTART;
 218   if (sigaction(SIGALRM,&sa,0)) diee("sigaction ALRM");
 219
 220   itv.it_interval.tv_sec= 0;
 221   itv.it_interval.tv_usec= 200000;
 222   itv.it_value= itv.it_interval;
 223
 224   if (setitimer(ITIMER_REAL,&itv,0)) diee("setitimer REAL");
 225
 226   raise(SIGALRM);
 227 }