1 /* multimin/simplex2.c
2 *
3 * Copyright (C) 2007, 2008, 2009 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 - Optimisations added by Brian Gough <bjg@network-theory.co.uk>
27 + use BLAS for frequently-called functions
28 + keep track of the center to avoid unnecessary computation
29 + compute size as RMS value, allowing linear update on each step
30 instead of recomputing from all N+1 vectors.
31 */
32
33 /* The Simplex method of Nelder and Mead, also known as the polytope
34 search alogorithm. Ref: Nelder, J.A., Mead, R., Computer Journal 7
35 (1965) pp. 308-313.
36
37 This implementation uses n+1 corner points in the simplex.
38 */
39
40 #include <config.h>
41 #include <stdlib.h>
42 #include <gsl/gsl_blas.h>
43 #include <gsl/gsl_multimin.h>
44 #include <gsl/gsl_matrix_double.h>
45
46 typedef struct
47 {
48 gsl_matrix *x1; /* simplex corner points */
49 gsl_vector *y1; /* function value at corner points */
50 gsl_vector *ws1; /* workspace 1 for algorithm */
51 gsl_vector *ws2; /* workspace 2 for algorithm */
52 gsl_vector *center; /* center of all points */
53 gsl_vector *delta; /* current step */
54 gsl_vector *xmc; /* x - center (workspace) */
55 double S2;
56 unsigned long count;
57 }
58 nmsimplex_state_t;
59
60 static int
61 compute_center (const nmsimplex_state_t * state, gsl_vector * center);
62 static double
63 compute_size (nmsimplex_state_t * state, const gsl_vector * center);
64
65 static double
try_corner_move(const double coeff,const nmsimplex_state_t * state,size_t corner,gsl_vector * xc,const gsl_multimin_function * f)66 try_corner_move (const double coeff,
67 const nmsimplex_state_t * state,
68 size_t corner,
69 gsl_vector * xc, const gsl_multimin_function * f)
70 {
71 /* moves a simplex corner scaled by coeff (negative value represents
72 mirroring by the middle point of the "other" corner points)
73 and gives new corner in xc and function value at xc as a
74 return value
75 */
76
77 gsl_matrix *x1 = state->x1;
78 const size_t P = x1->size1;
79 double newval;
80
81 /* xc = (1-coeff)*(P/(P-1)) * center(all) + ((P*coeff-1)/(P-1))*x_corner */
82 {
83 double alpha = (1 - coeff) * P / (P - 1.0);
84 double beta = (P * coeff - 1.0) / (P - 1.0);
85 gsl_vector_const_view row = gsl_matrix_const_row (x1, corner);
86
87 gsl_vector_memcpy (xc, state->center);
88 gsl_blas_dscal (alpha, xc);
89 gsl_blas_daxpy (beta, &row.vector, xc);
90 }
91
92 newval = GSL_MULTIMIN_FN_EVAL (f, xc);
93
94 return newval;
95 }
96
97
98 static void
update_point(nmsimplex_state_t * state,size_t i,const gsl_vector * x,double val)99 update_point (nmsimplex_state_t * state, size_t i,
100 const gsl_vector * x, double val)
101 {
102 gsl_vector_const_view x_orig = gsl_matrix_const_row (state->x1, i);
103 const size_t P = state->x1->size1;
104
105 /* Compute delta = x - x_orig */
106 gsl_vector_memcpy (state->delta, x);
107 gsl_blas_daxpy (-1.0, &x_orig.vector, state->delta);
108
109 /* Compute xmc = x_orig - c */
110 gsl_vector_memcpy (state->xmc, &x_orig.vector);
111 gsl_blas_daxpy (-1.0, state->center, state->xmc);
112
113 /* Update size: S2' = S2 + (2/P) * (x_orig - c).delta + (P-1)*(delta/P)^2 */
114 {
115 double d = gsl_blas_dnrm2 (state->delta);
116 double xmcd;
117 gsl_blas_ddot (state->xmc, state->delta, &xmcd);
118 state->S2 += (2.0 / P) * xmcd + ((P - 1.0) / P) * (d * d / P);
119 }
120
121 /* Update center: c' = c + (x - x_orig) / P */
122
123 {
124 double alpha = 1.0 / P;
125 gsl_blas_daxpy (-alpha, &x_orig.vector, state->center);
126 gsl_blas_daxpy (alpha, x, state->center);
127 }
128
129 gsl_matrix_set_row (state->x1, i, x);
130 gsl_vector_set (state->y1, i, val);
131 }
132
133 static int
contract_by_best(nmsimplex_state_t * state,size_t best,gsl_vector * xc,gsl_multimin_function * f)134 contract_by_best (nmsimplex_state_t * state, size_t best,
135 gsl_vector * xc, gsl_multimin_function * f)
136 {
137
138 /* Function contracts the simplex in respect to best valued
139 corner. That is, all corners besides the best corner are moved.
140 (This function is rarely called in practice, since it is the last
141 choice, hence not optimised - BJG) */
142
143 /* the xc vector is simply work space here */
144
145 gsl_matrix *x1 = state->x1;
146 gsl_vector *y1 = state->y1;
147
148 size_t i, j;
149 double newval;
150
151 int status = GSL_SUCCESS;
152
153 for (i = 0; i < x1->size1; i++)
154 {
155 if (i != best)
156 {
157 for (j = 0; j < x1->size2; j++)
158 {
159 newval = 0.5 * (gsl_matrix_get (x1, i, j)
160 + gsl_matrix_get (x1, best, j));
161 gsl_matrix_set (x1, i, j, newval);
162 }
163
164 /* evaluate function in the new point */
165
166 gsl_matrix_get_row (xc, x1, i);
167 newval = GSL_MULTIMIN_FN_EVAL (f, xc);
168 gsl_vector_set (y1, i, newval);
169
170 /* notify caller that we found at least one bad function value.
171 we finish the contraction (and do not abort) to allow the user
172 to handle the situation */
173
174 if (!gsl_finite (newval))
175 {
176 status = GSL_EBADFUNC;
177 }
178 }
179 }
180
181 /* We need to update the centre and size as well */
182 compute_center (state, state->center);
183 compute_size (state, state->center);
184
185 return status;
186 }
187
188 static int
compute_center(const nmsimplex_state_t * state,gsl_vector * center)189 compute_center (const nmsimplex_state_t * state, gsl_vector * center)
190 {
191 /* calculates the center of the simplex and stores in center */
192
193 gsl_matrix *x1 = state->x1;
194 const size_t P = x1->size1;
195 size_t i;
196
197 gsl_vector_set_zero (center);
198
199 for (i = 0; i < P; i++)
200 {
201 gsl_vector_const_view row = gsl_matrix_const_row (x1, i);
202 gsl_blas_daxpy (1.0, &row.vector, center);
203 }
204
205 {
206 const double alpha = 1.0 / P;
207 gsl_blas_dscal (alpha, center);
208 }
209
210 return GSL_SUCCESS;
211 }
212
213 static double
compute_size(nmsimplex_state_t * state,const gsl_vector * center)214 compute_size (nmsimplex_state_t * state, const gsl_vector * center)
215 {
216 /* calculates simplex size as rms sum of length of vectors
217 from simplex center to corner points:
218
219 sqrt( sum ( || y - y_middlepoint ||^2 ) / n )
220 */
221
222 gsl_vector *s = state->ws1;
223 gsl_matrix *x1 = state->x1;
224 const size_t P = x1->size1;
225 size_t i;
226
227 double ss = 0.0;
228
229 for (i = 0; i < P; i++)
230 {
231 double t;
232 gsl_matrix_get_row (s, x1, i);
233 gsl_blas_daxpy (-1.0, center, s);
234 t = gsl_blas_dnrm2 (s);
235 ss += t * t;
236 }
237
238 /* Store squared size in the state */
239 state->S2 = (ss / P);
240
241 return sqrt (ss / P);
242 }
243
244 static int
nmsimplex_alloc(void * vstate,size_t n)245 nmsimplex_alloc (void *vstate, size_t n)
246 {
247 nmsimplex_state_t *state = (nmsimplex_state_t *) vstate;
248
249 if (n == 0)
250 {
251 GSL_ERROR ("invalid number of parameters specified", GSL_EINVAL);
252 }
253
254 state->x1 = gsl_matrix_alloc (n + 1, n);
255
256 if (state->x1 == NULL)
257 {
258 GSL_ERROR ("failed to allocate space for x1", GSL_ENOMEM);
259 }
260
261 state->y1 = gsl_vector_alloc (n + 1);
262
263 if (state->y1 == NULL)
264 {
265 gsl_matrix_free (state->x1);
266 GSL_ERROR ("failed to allocate space for y", GSL_ENOMEM);
267 }
268
269 state->ws1 = gsl_vector_alloc (n);
270
271 if (state->ws1 == NULL)
272 {
273 gsl_matrix_free (state->x1);
274 gsl_vector_free (state->y1);
275 GSL_ERROR ("failed to allocate space for ws1", GSL_ENOMEM);
276 }
277
278 state->ws2 = gsl_vector_alloc (n);
279
280 if (state->ws2 == NULL)
281 {
282 gsl_matrix_free (state->x1);
283 gsl_vector_free (state->y1);
284 gsl_vector_free (state->ws1);
285 GSL_ERROR ("failed to allocate space for ws2", GSL_ENOMEM);
286 }
287
288 state->center = gsl_vector_alloc (n);
289
290 if (state->center == NULL)
291 {
292 gsl_matrix_free (state->x1);
293 gsl_vector_free (state->y1);
294 gsl_vector_free (state->ws1);
295 gsl_vector_free (state->ws2);
296 GSL_ERROR ("failed to allocate space for center", GSL_ENOMEM);
297 }
298
299 state->delta = gsl_vector_alloc (n);
300
301 if (state->delta == NULL)
302 {
303 gsl_matrix_free (state->x1);
304 gsl_vector_free (state->y1);
305 gsl_vector_free (state->ws1);
306 gsl_vector_free (state->ws2);
307 gsl_vector_free (state->center);
308 GSL_ERROR ("failed to allocate space for delta", GSL_ENOMEM);
309 }
310
311 state->xmc = gsl_vector_alloc (n);
312
313 if (state->xmc == NULL)
314 {
315 gsl_matrix_free (state->x1);
316 gsl_vector_free (state->y1);
317 gsl_vector_free (state->ws1);
318 gsl_vector_free (state->ws2);
319 gsl_vector_free (state->center);
320 gsl_vector_free (state->delta);
321 GSL_ERROR ("failed to allocate space for xmc", GSL_ENOMEM);
322 }
323
324 state->count = 0;
325
326 return GSL_SUCCESS;
327 }
328
329 static void
nmsimplex_free(void * vstate)330 nmsimplex_free (void *vstate)
331 {
332 nmsimplex_state_t *state = (nmsimplex_state_t *) vstate;
333
334 gsl_matrix_free (state->x1);
335 gsl_vector_free (state->y1);
336 gsl_vector_free (state->ws1);
337 gsl_vector_free (state->ws2);
338 gsl_vector_free (state->center);
339 gsl_vector_free (state->delta);
340 gsl_vector_free (state->xmc);
341 }
342
343 static int
nmsimplex_set(void * vstate,gsl_multimin_function * f,const gsl_vector * x,double * size,const gsl_vector * step_size)344 nmsimplex_set (void *vstate, gsl_multimin_function * f,
345 const gsl_vector * x,
346 double *size, const gsl_vector * step_size)
347 {
348 int status;
349 size_t i;
350 double val;
351
352 nmsimplex_state_t *state = (nmsimplex_state_t *) vstate;
353
354 gsl_vector *xtemp = state->ws1;
355
356 if (xtemp->size != x->size)
357 {
358 GSL_ERROR ("incompatible size of x", GSL_EINVAL);
359 }
360
361 if (xtemp->size != step_size->size)
362 {
363 GSL_ERROR ("incompatible size of step_size", GSL_EINVAL);
364 }
365
366 /* first point is the original x0 */
367
368 val = GSL_MULTIMIN_FN_EVAL (f, x);
369
370 if (!gsl_finite (val))
371 {
372 GSL_ERROR ("non-finite function value encountered", GSL_EBADFUNC);
373 }
374
375 gsl_matrix_set_row (state->x1, 0, x);
376 gsl_vector_set (state->y1, 0, val);
377
378 /* following points are initialized to x0 + step_size */
379
380 for (i = 0; i < x->size; i++)
381 {
382 status = gsl_vector_memcpy (xtemp, x);
383
384 if (status != 0)
385 {
386 GSL_ERROR ("vector memcopy failed", GSL_EFAILED);
387 }
388
389 {
390 double xi = gsl_vector_get (x, i);
391 double si = gsl_vector_get (step_size, i);
392
393 gsl_vector_set (xtemp, i, xi + si);
394 val = GSL_MULTIMIN_FN_EVAL (f, xtemp);
395 }
396
397 if (!gsl_finite (val))
398 {
399 GSL_ERROR ("non-finite function value encountered", GSL_EBADFUNC);
400 }
401
402 gsl_matrix_set_row (state->x1, i + 1, xtemp);
403 gsl_vector_set (state->y1, i + 1, val);
404 }
405
406 compute_center (state, state->center);
407
408 /* Initialize simplex size */
409 *size = compute_size (state, state->center);
410
411 state->count++;
412
413 return GSL_SUCCESS;
414 }
415
416 static int
nmsimplex_iterate(void * vstate,gsl_multimin_function * f,gsl_vector * x,double * size,double * fval)417 nmsimplex_iterate (void *vstate, gsl_multimin_function * f,
418 gsl_vector * x, double *size, double *fval)
419 {
420
421 /* Simplex iteration tries to minimize function f value */
422 /* Includes corrections from Ivo Alxneit <ivo.alxneit@psi.ch> */
423
424 nmsimplex_state_t *state = (nmsimplex_state_t *) vstate;
425
426 /* xc and xc2 vectors store tried corner point coordinates */
427
428 gsl_vector *xc = state->ws1;
429 gsl_vector *xc2 = state->ws2;
430 gsl_vector *y1 = state->y1;
431 gsl_matrix *x1 = state->x1;
432
433 const size_t n = y1->size;
434 size_t i;
435 size_t hi, s_hi, lo;
436 double dhi, ds_hi, dlo;
437 int status;
438 double val, val2;
439
440 if (xc->size != x->size)
441 {
442 GSL_ERROR ("incompatible size of x", GSL_EINVAL);
443 }
444
445 /* get index of highest, second highest and lowest point */
446
447 dhi = dlo = gsl_vector_get (y1, 0);
448 hi = 0;
449 lo = 0;
450
451 ds_hi = gsl_vector_get (y1, 1);
452 s_hi = 1;
453
454 for (i = 1; i < n; i++)
455 {
456 val = (gsl_vector_get (y1, i));
457 if (val < dlo)
458 {
459 dlo = val;
460 lo = i;
461 }
462 else if (val > dhi)
463 {
464 ds_hi = dhi;
465 s_hi = hi;
466 dhi = val;
467 hi = i;
468 }
469 else if (val > ds_hi)
470 {
471 ds_hi = val;
472 s_hi = i;
473 }
474 }
475
476 /* try reflecting the highest value point */
477
478 val = try_corner_move (-1.0, state, hi, xc, f);
479
480 if (gsl_finite (val) && val < gsl_vector_get (y1, lo))
481 {
482 /* reflected point is lowest, try expansion */
483
484 val2 = try_corner_move (-2.0, state, hi, xc2, f);
485
486 if (gsl_finite (val2) && val2 < gsl_vector_get (y1, lo))
487 {
488 update_point (state, hi, xc2, val2);
489 }
490 else
491 {
492 update_point (state, hi, xc, val);
493 }
494 }
495 else if (!gsl_finite (val) || val > gsl_vector_get (y1, s_hi))
496 {
497 /* reflection does not improve things enough, or we got a
498 non-finite function value */
499
500 if (gsl_finite (val) && val <= gsl_vector_get (y1, hi))
501 {
502 /* if trial point is better than highest point, replace
503 highest point */
504
505 update_point (state, hi, xc, val);
506 }
507
508 /* try one-dimensional contraction */
509
510 val2 = try_corner_move (0.5, state, hi, xc2, f);
511
512 if (gsl_finite (val2) && val2 <= gsl_vector_get (y1, hi))
513 {
514 update_point (state, hi, xc2, val2);
515 }
516 else
517 {
518 /* contract the whole simplex about the best point */
519
520 status = contract_by_best (state, lo, xc, f);
521
522 if (status != GSL_SUCCESS)
523 {
524 GSL_ERROR ("contraction failed", GSL_EFAILED);
525 }
526 }
527 }
528 else
529 {
530 /* trial point is better than second highest point. Replace
531 highest point by it */
532
533 update_point (state, hi, xc, val);
534 }
535
536 /* return lowest point of simplex as x */
537
538 lo = gsl_vector_min_index (y1);
539 gsl_matrix_get_row (x, x1, lo);
540 *fval = gsl_vector_get (y1, lo);
541
542 /* Update simplex size */
543
544 {
545 double S2 = state->S2;
546
547 if (S2 > 0)
548 {
549 *size = sqrt (S2);
550 }
551 else
552 {
553 /* recompute if accumulated error has made size invalid */
554 *size = compute_size (state, state->center);
555 }
556 }
557
558 return GSL_SUCCESS;
559 }
560
561 static const gsl_multimin_fminimizer_type nmsimplex_type =
562 { "nmsimplex2", /* name */
563 sizeof (nmsimplex_state_t),
564 &nmsimplex_alloc,
565 &nmsimplex_set,
566 &nmsimplex_iterate,
567 &nmsimplex_free
568 };
569
570 const gsl_multimin_fminimizer_type
571 * gsl_multimin_fminimizer_nmsimplex2 = &nmsimplex_type;
572
573
574 static inline double
ran_unif(unsigned long * seed)575 ran_unif (unsigned long *seed)
576 {
577 unsigned long s = *seed;
578 *seed = (s * 69069 + 1) & 0xffffffffUL;
579 return (*seed) / 4294967296.0;
580 }
581
582 static int
nmsimplex_set_rand(void * vstate,gsl_multimin_function * f,const gsl_vector * x,double * size,const gsl_vector * step_size)583 nmsimplex_set_rand (void *vstate, gsl_multimin_function * f,
584 const gsl_vector * x,
585 double *size, const gsl_vector * step_size)
586 {
587 size_t i, j;
588 double val;
589
590 nmsimplex_state_t *state = (nmsimplex_state_t *) vstate;
591
592 gsl_vector *xtemp = state->ws1;
593
594 if (xtemp->size != x->size)
595 {
596 GSL_ERROR ("incompatible size of x", GSL_EINVAL);
597 }
598
599 if (xtemp->size != step_size->size)
600 {
601 GSL_ERROR ("incompatible size of step_size", GSL_EINVAL);
602 }
603
604 /* first point is the original x0 */
605
606 val = GSL_MULTIMIN_FN_EVAL (f, x);
607
608 if (!gsl_finite (val))
609 {
610 GSL_ERROR ("non-finite function value encountered", GSL_EBADFUNC);
611 }
612
613 gsl_matrix_set_row (state->x1, 0, x);
614 gsl_vector_set (state->y1, 0, val);
615
616 {
617 gsl_matrix_view m =
618 gsl_matrix_submatrix (state->x1, 1, 0, x->size, x->size);
619
620 /* generate a random orthornomal basis */
621 unsigned long seed = state->count ^ 0x12345678;
622
623 ran_unif (&seed); /* warm it up */
624
625 gsl_matrix_set_identity (&m.matrix);
626
627 /* start with random reflections */
628 for (i = 0; i < x->size; i++)
629 {
630 double s = ran_unif (&seed);
631 if (s > 0.5) gsl_matrix_set (&m.matrix, i, i, -1.0);
632 }
633
634 /* apply random rotations */
635 for (i = 0; i < x->size; i++)
636 {
637 for (j = i + 1; j < x->size; j++)
638 {
639 /* rotate columns i and j by a random angle */
640 double angle = 2.0 * M_PI * ran_unif (&seed);
641 double c = cos (angle), s = sin (angle);
642 gsl_vector_view c_i = gsl_matrix_column (&m.matrix, i);
643 gsl_vector_view c_j = gsl_matrix_column (&m.matrix, j);
644 gsl_blas_drot (&c_i.vector, &c_j.vector, c, s);
645 }
646 }
647
648 /* scale the orthonormal basis by the user-supplied step_size in
649 each dimension, and use as an offset from the central point x */
650
651 for (i = 0; i < x->size; i++)
652 {
653 double x_i = gsl_vector_get (x, i);
654 double s_i = gsl_vector_get (step_size, i);
655 gsl_vector_view c_i = gsl_matrix_column (&m.matrix, i);
656
657 for (j = 0; j < x->size; j++)
658 {
659 double x_ij = gsl_vector_get (&c_i.vector, j);
660 gsl_vector_set (&c_i.vector, j, x_i + s_i * x_ij);
661 }
662 }
663
664 /* compute the function values at each offset point */
665
666 for (i = 0; i < x->size; i++)
667 {
668 gsl_vector_view r_i = gsl_matrix_row (&m.matrix, i);
669
670 val = GSL_MULTIMIN_FN_EVAL (f, &r_i.vector);
671
672 if (!gsl_finite (val))
673 {
674 GSL_ERROR ("non-finite function value encountered", GSL_EBADFUNC);
675 }
676
677 gsl_vector_set (state->y1, i + 1, val);
678 }
679 }
680
681 compute_center (state, state->center);
682
683 /* Initialize simplex size */
684 *size = compute_size (state, state->center);
685
686 state->count++;
687
688 return GSL_SUCCESS;
689 }
690
691 static const gsl_multimin_fminimizer_type nmsimplex2rand_type =
692 { "nmsimplex2rand", /* name */
693 sizeof (nmsimplex_state_t),
694 &nmsimplex_alloc,
695 &nmsimplex_set_rand,
696 &nmsimplex_iterate,
697 &nmsimplex_free
698 };
699
700 const gsl_multimin_fminimizer_type
701 * gsl_multimin_fminimizer_nmsimplex2rand = &nmsimplex2rand_type;
702