chiark / gitweb /
exploit symmetry
[moebius2.git] / 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))
115       printf("%*s size %# e, r=%d\n", 135,"", 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) {
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("[%4d] ",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 }