htswanalysis/MACS/lib/gsl/gsl-1.11/multimin/simplex.c

   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
  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
  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
 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
 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
 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
 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
 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
 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;