chiark / gitweb /
set threads_started
[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, *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 int main(int argc, const char *const *argv) {
66   gsl_multimin_function multimin_function;
67   double size;
68   struct Vertices initial_full;
69   double initial_half[DIM], step_size[DIM];
70   FILE *initial_f;
71   gsl_vector initial_gsl, step_size_gsl;
72   int r, i;
73
74   if (argc==3) {
75     input_file= argv[1];
76     if (strncmp(argv[2],"-o",2)) badusage();   best_file= argv[2]+2;
77     final_file= 0;
78   } else if (argc==4) {
79     input_file= argv[1];
80     if (strncmp(argv[2],"-i",2)) badusage();   best_file= argv[2]+2;
81     if (strncmp(argv[3],"-o",2)) badusage();   final_file= argv[3]+2;
82   } else {
83     badusage();
84   }
85   if (argv[1][0]=='-') badusage();
86
87   if (asprintf(&best_file_tmp,"%s.new",best_file) <= 0) diee("asprintf");
88   if (final_file)
89     if (asprintf(&final_file_tmp,"%s.new",final_file) <= 0) diee("asprintf");
90
91   graph_layout_prepare();
92   printing_init();
93   energy_init();
94
95   printf("X=%d=0x%x  Y=%d=0x%x  DIM=%d\n",X,X,Y,Y,DIM);
96   
97   minimiser= gsl_multimin_fminimizer_alloc
98     (gsl_multimin_fminimizer_nmsimplex, DIM);
99   if (!minimiser) { perror("alloc minimiser"); exit(-1); }
100
101   multimin_function.f= minfunc_f;
102   multimin_function.n= DIM;
103   multimin_function.params= 0;
104
105   initial_f= fopen(input_file,"rb");  if (!initial_f) diee("fopen initial");
106   errno= 0; r= fread(&initial_full,sizeof(initial_full),1,initial_f);
107   if (r!=1) diee("fread");
108   fclose(initial_f);
109
110   pmap_dimensions(&initial_full);
111   unmap_dimensions(initial_half,&initial_full);
112   for (i=0; i<DIM; i++) step_size[i]= 0.03;
113
114   initial_gsl.size= DIM;
115   initial_gsl.stride= 1;
116   initial_gsl.block= 0;
117   initial_gsl.owner= 0;
118   step_size_gsl= initial_gsl;
119
120   initial_gsl.data= initial_half;
121   step_size_gsl.data= step_size;
122
123   GA( gsl_multimin_fminimizer_set(minimiser, &multimin_function,
124                                   &initial_gsl, &step_size_gsl) );
125
126   for (;;) {
127     GA( gsl_multimin_fminimizer_iterate(minimiser) );
128
129     size= gsl_multimin_fminimizer_size(minimiser);
130     r= gsl_multimin_test_size(size, stop_epsilon);
131
132     if (printing_check(pr_size,215))
133       printf("r=%2d, size %# e\n", r, size);
134     flushoutput();
135
136     if (r==GSL_SUCCESS) break;
137     assert(r==GSL_CONTINUE);
138   }
139
140   if (final_file) {
141     if (unlink(final_file_tmp) && errno != ENOENT) diee("unlink final .tmp");
142     if (link(best_file,final_file_tmp)) diee("link final .tmp");
143     if (rename(final_file_tmp,final_file)) diee("install final");
144   }
145   printf("FINISHED (%lld evaluations)\n",evaluations);
146   
147   return 0;
148 }
149
150 /*---------- printing rate limit ----------*/
151
152 static volatile unsigned print_todo;
153 static sigset_t print_alarmset;
154
155 int printing_check(enum printing_instance which, int indent) {
156   static int skipped[pr__max];
157   
158   unsigned bits= 1u << which;
159   int sk;
160
161   if (!(print_todo & bits)) {
162     skipped[which]++;
163     return 0;;
164   }
165
166   sigprocmask(SIG_BLOCK,&print_alarmset,0);
167   print_todo &= ~bits;
168   sigprocmask(SIG_UNBLOCK,&print_alarmset,0);
169
170   printf("%*s",indent,"");
171   sk= skipped[which];
172   if (sk) printf("[%4d] ",sk);
173   else printf("      ");
174   skipped[which]= 0;
175
176   return 1;
177 }
178
179 static void alarmhandler(int ignored) {
180   print_todo= ~0u;
181 }
182
183 static void printing_init(void) {
184   struct sigaction sa;
185   struct itimerval itv;
186
187   sigemptyset(&print_alarmset);
188   sigaddset(&print_alarmset,SIGALRM);
189
190   sa.sa_handler= alarmhandler;
191   sa.sa_mask= print_alarmset;
192   sa.sa_flags= SA_RESTART;
193   if (sigaction(SIGALRM,&sa,0)) diee("sigaction ALRM");
194   
195   itv.it_interval.tv_sec= 0;
196   itv.it_interval.tv_usec= 200000;
197   itv.it_value= itv.it_interval;
198
199   if (setitimer(ITIMER_REAL,&itv,0)) diee("setitimer REAL");
200
201   raise(SIGALRM);
202 }