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 "gsl__config.h"
36 #include <stdlib.h>
37 #include "gsl_blas.h"
38 #include "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 = 0, s_hi = 0, lo = 0;
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 = ds_hi = dlo = gsl_vector_get (y1, 0);
346
347 for (i = 1; i < n; i++)
348 {
349 val = (gsl_vector_get (y1, i));
350 if (val < dlo)
351 {
352 dlo = val;
353 lo = i;
354 }
355 else if (val > dhi)
356 {
357 ds_hi = dhi;
358 s_hi = hi;
359 dhi = val;
360 hi = i;
361 }
362 else if (val > ds_hi)
363 {
364 ds_hi = val;
365 s_hi = i;
366 }
367 }
368
369 /* reflect the highest value */
370
371 val = nmsimplex_move_corner (-1.0, state, hi, xc, f);
372
373 if (gsl_finite(val) && val < gsl_vector_get (y1, lo))
374 {
375
376 /* reflected point becomes lowest point, try expansion */
377
378 val2 = nmsimplex_move_corner (-2.0, state, hi, xc2, f);
379
380 if (gsl_finite(val2) && val2 < gsl_vector_get (y1, lo))
381 {
382 gsl_matrix_set_row (x1, hi, xc2);
383 gsl_vector_set (y1, hi, val2);
384 }
385 else
386 {
387 gsl_matrix_set_row (x1, hi, xc);
388 gsl_vector_set (y1, hi, val);
389 }
390 }
391
392 /* reflection does not improve things enough
393 or
394 we got a non-finite (illegal) function value */
395
396 else if (!gsl_finite(val) || val > gsl_vector_get (y1, s_hi))
397 {
398 if (gsl_finite(val) && val <= gsl_vector_get (y1, hi))
399 {
400
401 /* if trial point is better than highest point, replace
402 highest point */
403
404 gsl_matrix_set_row (x1, hi, xc);
405 gsl_vector_set (y1, hi, val);
406 }
407
408 /* try one dimensional contraction */
409
410 val2 = nmsimplex_move_corner (0.5, state, hi, xc2, f);
411
412 if (gsl_finite(val2) && val2 <= gsl_vector_get (y1, hi))
413 {
414 gsl_matrix_set_row (state->x1, hi, xc2);
415 gsl_vector_set (y1, hi, val2);
416 }
417
418 else
419 {
420
421 /* contract the whole simplex in respect to the best point */
422
423 status = nmsimplex_contract_by_best (state, lo, xc, f);
424 if (status != GSL_SUCCESS)
425 {
426 GSL_ERROR ("nmsimplex_contract_by_best failed", GSL_EFAILED);
427 }
428 }
429 }
430 else
431 {
432
433 /* trial point is better than second highest point.
434 Replace highest point by it */
435
436 gsl_matrix_set_row (x1, hi, xc);
437 gsl_vector_set (y1, hi, val);
438 }
439
440 /* return lowest point of simplex as x */
441
442 lo = gsl_vector_min_index (y1);
443 gsl_matrix_get_row (x, x1, lo);
444 *fval = gsl_vector_get (y1, lo);
445
446 /* Update simplex size */
447
448 *size = nmsimplex_size (state);
449
450 return GSL_SUCCESS;
451 }
452
453 static const gsl_multimin_fminimizer_type nmsimplex_type =
454 { "nmsimplex", /* name */
455 sizeof (nmsimplex_state_t),
456 &nmsimplex_alloc,
457 &nmsimplex_set,
458 &nmsimplex_iterate,
459 &nmsimplex_free
460 };
461
462 const gsl_multimin_fminimizer_type
463 * gsl_multimin_fminimizer_nmsimplex = &nmsimplex_type;
464