1 /* multimin/simplex.c
2 *
3 * Copyright (C) 2007 Brian Gough
4 * Copyright (C) 2002 Tuomo Keskitalo, Ivo Alxneit
5 *
6 * This program is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation; either version 3 of the License, or (at
9 * your option) any later version.
10 *
11 * This program is distributed in the hope that it will be useful, but
12 * WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 * General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with this program; if not, write to the Free Software
18 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
19 */
20
21 /*
22 - Originally written by Tuomo Keskitalo <tuomo.keskitalo@iki.fi>
23 - Corrections to nmsimplex_iterate and other functions
24 by Ivo Alxneit <ivo.alxneit@psi.ch>
25 - Additional help by Brian Gough <bjg@network-theory.co.uk>
26 */
27
28 /* The Simplex method of Nelder and Mead,
29 also known as the polytope search alogorithm. Ref:
30 Nelder, J.A., Mead, R., Computer Journal 7 (1965) pp. 308-313.
31
32 This implementation uses n+1 corner points in the simplex.
33 */
34
35 #include <config.h>
36 #include <stdlib.h>
37 #include <gsl/gsl_blas.h>
38 #include <gsl/gsl_multimin.h>
39
40 typedef struct
41 {
42 gsl_matrix *x1; /* simplex corner points */
43 gsl_vector *y1; /* function value at corner points */
44 gsl_vector *ws1; /* workspace 1 for algorithm */
45 gsl_vector *ws2; /* workspace 2 for algorithm */
46 }
47 nmsimplex_state_t;
48
49 static double
nmsimplex_move_corner(const double coeff,const nmsimplex_state_t * state,size_t corner,gsl_vector * xc,const gsl_multimin_function * f)50 nmsimplex_move_corner (const double coeff, const nmsimplex_state_t * state,
51 size_t corner, gsl_vector * xc,
52 const gsl_multimin_function * f)
53 {
54 /* moves a simplex corner scaled by coeff (negative value represents
55 mirroring by the middle point of the "other" corner points)
56 and gives new corner in xc and function value at xc as a
57 return value
58 */
59
60 gsl_matrix *x1 = state->x1;
61
62 size_t i, j;
63 double newval, mp;
64
65 for (j = 0; j < x1->size2; j++)
66 {
67 mp = 0.0;
68 for (i = 0; i < x1->size1; i++)
69 {
70 if (i != corner)
71 {
72 mp += (gsl_matrix_get (x1, i, j));
73 }
74 }
75 mp /= (double) (x1->size1 - 1);
76 newval = mp - coeff * (mp - gsl_matrix_get (x1, corner, j));
77 gsl_vector_set (xc, j, newval);
78 }
79
80 newval = GSL_MULTIMIN_FN_EVAL (f, xc);
81
82 return newval;
83 }
84
85 static int
nmsimplex_contract_by_best(nmsimplex_state_t * state,size_t best,gsl_vector * xc,gsl_multimin_function * f)86 nmsimplex_contract_by_best (nmsimplex_state_t * state, size_t best,
87 gsl_vector * xc, gsl_multimin_function * f)
88 {
89
90 /* Function contracts the simplex in respect to
91 best valued corner. That is, all corners besides the
92 best corner are moved. */
93
94 /* the xc vector is simply work space here */
95
96 gsl_matrix *x1 = state->x1;
97 gsl_vector *y1 = state->y1;
98
99 size_t i, j;
100 double newval;
101
102 int status = GSL_SUCCESS;
103
104 for (i = 0; i < x1->size1; i++)
105 {
106 if (i != best)
107 {
108 for (j = 0; j < x1->size2; j++)
109 {
110 newval = 0.5 * (gsl_matrix_get (x1, i, j)
111 + gsl_matrix_get (x1, best, j));
112 gsl_matrix_set (x1, i, j, newval);
113 }
114
115 /* evaluate function in the new point */
116
117 gsl_matrix_get_row (xc, x1, i);
118 newval = GSL_MULTIMIN_FN_EVAL (f, xc);
119 gsl_vector_set (y1, i, newval);
120
121 /* notify caller that we found at least one bad function value.
122 we finish the contraction (and do not abort) to allow the user
123 to handle the situation */
124
125 if(!gsl_finite(newval))
126 {
127 status = GSL_EBADFUNC;
128 }
129 }
130 }
131
132 return status;
133 }
134
135 static int
nmsimplex_calc_center(const nmsimplex_state_t * state,gsl_vector * mp)136 nmsimplex_calc_center (const nmsimplex_state_t * state, gsl_vector * mp)
137 {
138 /* calculates the center of the simplex to mp */
139
140 gsl_matrix *x1 = state->x1;
141
142 size_t i, j;
143 double val;
144
145 for (j = 0; j < x1->size2; j++)
146 {
147 val = 0.0;
148 for (i = 0; i < x1->size1; i++)
149 {
150 val += gsl_matrix_get (x1, i, j);
151 }
152 val /= x1->size1;
153 gsl_vector_set (mp, j, val);
154 }
155
156 return GSL_SUCCESS;
157 }
158
159 static double
nmsimplex_size(nmsimplex_state_t * state)160 nmsimplex_size (nmsimplex_state_t * state)
161 {
162 /* calculates simplex size as average sum of length of vectors
163 from simplex center to corner points:
164
165 ( sum ( || y - y_middlepoint || ) ) / n
166 */
167
168 gsl_vector *s = state->ws1;
169 gsl_vector *mp = state->ws2;
170
171 gsl_matrix *x1 = state->x1;
172 size_t i;
173
174 double ss = 0.0;
175
176 /* Calculate middle point */
177 nmsimplex_calc_center (state, mp);
178
179 for (i = 0; i < x1->size1; i++)
180 {
181 gsl_matrix_get_row (s, x1, i);
182 gsl_blas_daxpy (-1.0, mp, s);
183 ss += gsl_blas_dnrm2 (s);
184 }
185
186 return ss / (double) (x1->size1);
187 }
188
189 static int
nmsimplex_alloc(void * vstate,size_t n)190 nmsimplex_alloc (void *vstate, size_t n)
191 {
192 nmsimplex_state_t *state = (nmsimplex_state_t *) vstate;
193
194 if (n == 0)
195 {
196 GSL_ERROR("invalid number of parameters specified", GSL_EINVAL);
197 }
198
199 state->x1 = gsl_matrix_alloc (n + 1, n);
200
201 if (state->x1 == NULL)
202 {
203 GSL_ERROR ("failed to allocate space for x1", GSL_ENOMEM);
204 }
205
206 state->y1 = gsl_vector_alloc (n + 1);
207
208 if (state->y1 == NULL)
209 {
210 gsl_matrix_free(state->x1);
211 GSL_ERROR ("failed to allocate space for y", GSL_ENOMEM);
212 }
213
214 state->ws1 = gsl_vector_alloc (n);
215
216 if (state->ws1 == NULL)
217 {
218 gsl_matrix_free(state->x1);
219 gsl_vector_free(state->y1);
220 GSL_ERROR ("failed to allocate space for ws1", GSL_ENOMEM);
221 }
222
223 state->ws2 = gsl_vector_alloc (n);
224
225 if (state->ws2 == NULL)
226 {
227 gsl_matrix_free(state->x1);
228 gsl_vector_free(state->y1);
229 gsl_vector_free(state->ws1);
230 GSL_ERROR ("failed to allocate space for ws2", GSL_ENOMEM);
231 }
232
233 return GSL_SUCCESS;
234 }
235
236 static int
nmsimplex_set(void * vstate,gsl_multimin_function * f,const gsl_vector * x,double * size,const gsl_vector * step_size)237 nmsimplex_set (void *vstate, gsl_multimin_function * f,
238 const gsl_vector * x,
239 double *size, const gsl_vector * step_size)
240 {
241 int status;
242 size_t i;
243 double val;
244
245 nmsimplex_state_t *state = (nmsimplex_state_t *) vstate;
246
247 gsl_vector *xtemp = state->ws1;
248
249 if (xtemp->size != x->size)
250 {
251 GSL_ERROR("incompatible size of x", GSL_EINVAL);
252 }
253
254 if (xtemp->size != step_size->size)
255 {
256 GSL_ERROR("incompatible size of step_size", GSL_EINVAL);
257 }
258
259 /* first point is the original x0 */
260
261 val = GSL_MULTIMIN_FN_EVAL (f, x);
262
263 if (!gsl_finite(val))
264 {
265 GSL_ERROR("non-finite function value encountered", GSL_EBADFUNC);
266 }
267
268 gsl_matrix_set_row (state->x1, 0, x);
269 gsl_vector_set (state->y1, 0, val);
270
271 /* following points are initialized to x0 + step_size */
272
273 for (i = 0; i < x->size; i++)
274 {
275 status = gsl_vector_memcpy (xtemp, x);
276
277 if (status != 0)
278 {
279 GSL_ERROR ("vector memcopy failed", GSL_EFAILED);
280 }
281
282 val = gsl_vector_get (xtemp, i) + gsl_vector_get (step_size, i);
283 gsl_vector_set (xtemp, i, val);
284 val = GSL_MULTIMIN_FN_EVAL (f, xtemp);
285
286 if (!gsl_finite(val))
287 {
288 GSL_ERROR("non-finite function value encountered", GSL_EBADFUNC);
289 }
290
291 gsl_matrix_set_row (state->x1, i + 1, xtemp);
292 gsl_vector_set (state->y1, i + 1, val);
293 }
294
295 /* Initialize simplex size */
296
297 *size = nmsimplex_size (state);
298
299 return GSL_SUCCESS;
300 }
301
302 static void
nmsimplex_free(void * vstate)303 nmsimplex_free (void *vstate)
304 {
305 nmsimplex_state_t *state = (nmsimplex_state_t *) vstate;
306
307 gsl_matrix_free (state->x1);
308 gsl_vector_free (state->y1);
309 gsl_vector_free (state->ws1);
310 gsl_vector_free (state->ws2);
311 }
312
313 static int
nmsimplex_iterate(void * vstate,gsl_multimin_function * f,gsl_vector * x,double * size,double * fval)314 nmsimplex_iterate (void *vstate, gsl_multimin_function * f,
315 gsl_vector * x, double *size, double *fval)
316 {
317
318 /* Simplex iteration tries to minimize function f value */
319 /* Includes corrections from Ivo Alxneit <ivo.alxneit@psi.ch> */
320
321 nmsimplex_state_t *state = (nmsimplex_state_t *) vstate;
322
323 /* xc and xc2 vectors store tried corner point coordinates */
324
325 gsl_vector *xc = state->ws1;
326 gsl_vector *xc2 = state->ws2;
327 gsl_vector *y1 = state->y1;
328 gsl_matrix *x1 = state->x1;
329
330 size_t n = y1->size;
331 size_t i;
332 size_t hi, s_hi, lo;
333 double dhi, ds_hi, dlo;
334 int status;
335 double val, val2;
336
337
338 if (xc->size != x->size)
339 {
340 GSL_ERROR("incompatible size of x", GSL_EINVAL);
341 }
342
343 /* get index of highest, second highest and lowest point */
344
345 dhi = dlo = gsl_vector_get (y1, 0);
346 hi = 0; lo = 0;
347
348 ds_hi = gsl_vector_get(y1, 1);
349 s_hi = 1;
350
351 for (i = 1; i < n; i++)
352 {
353 val = (gsl_vector_get (y1, i));
354 if (val < dlo)
355 {
356 dlo = val;
357 lo = i;
358 }
359 else if (val > dhi)
360 {
361 ds_hi = dhi;
362 s_hi = hi;
363 dhi = val;
364 hi = i;
365 }
366 else if (val > ds_hi)
367 {
368 ds_hi = val;
369 s_hi = i;
370 }
371 }
372
373 /* reflect the highest value */
374
375 val = nmsimplex_move_corner (-1.0, state, hi, xc, f);
376
377 if (gsl_finite(val) && val < gsl_vector_get (y1, lo))
378 {
379
380 /* reflected point becomes lowest point, try expansion */
381
382 val2 = nmsimplex_move_corner (-2.0, state, hi, xc2, f);
383
384 if (gsl_finite(val2) && val2 < gsl_vector_get (y1, lo))
385 {
386 gsl_matrix_set_row (x1, hi, xc2);
387 gsl_vector_set (y1, hi, val2);
388 }
389 else
390 {
391 gsl_matrix_set_row (x1, hi, xc);
392 gsl_vector_set (y1, hi, val);
393 }
394 }
395
396 /* reflection does not improve things enough
397 or
398 we got a non-finite (illegal) function value */
399
400 else if (!gsl_finite(val) || val > gsl_vector_get (y1, s_hi))
401 {
402 if (gsl_finite(val) && val <= gsl_vector_get (y1, hi))
403 {
404
405 /* if trial point is better than highest point, replace
406 highest point */
407
408 gsl_matrix_set_row (x1, hi, xc);
409 gsl_vector_set (y1, hi, val);
410 }
411
412 /* try one dimensional contraction */
413
414 val2 = nmsimplex_move_corner (0.5, state, hi, xc2, f);
415
416 if (gsl_finite(val2) && val2 <= gsl_vector_get (y1, hi))
417 {
418 gsl_matrix_set_row (state->x1, hi, xc2);
419 gsl_vector_set (y1, hi, val2);
420 }
421
422 else
423 {
424
425 /* contract the whole simplex in respect to the best point */
426
427 status = nmsimplex_contract_by_best (state, lo, xc, f);
428 if (status != GSL_SUCCESS)
429 {
430 GSL_ERROR ("nmsimplex_contract_by_best failed", GSL_EFAILED);
431 }
432 }
433 }
434 else
435 {
436
437 /* trial point is better than second highest point.
438 Replace highest point by it */
439
440 gsl_matrix_set_row (x1, hi, xc);
441 gsl_vector_set (y1, hi, val);
442 }
443
444 /* return lowest point of simplex as x */
445
446 lo = gsl_vector_min_index (y1);
447 gsl_matrix_get_row (x, x1, lo);
448 *fval = gsl_vector_get (y1, lo);
449
450 /* Update simplex size */
451
452 *size = nmsimplex_size (state);
453
454 return GSL_SUCCESS;
455 }
456
457 static const gsl_multimin_fminimizer_type nmsimplex_type =
458 { "nmsimplex", /* name */
459 sizeof (nmsimplex_state_t),
460 &nmsimplex_alloc,
461 &nmsimplex_set,
462 &nmsimplex_iterate,
463 &nmsimplex_free
464 };
465
466 const gsl_multimin_fminimizer_type
467 * gsl_multimin_fminimizer_nmsimplex = &nmsimplex_type;
468