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

   1 /* multimin/test.c
   2  *
   3  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Fabrice Rossi
   4  *
   5  * This program is free software; you can redistribute it and/or modify
   6  * it under the terms of the GNU General Public License as published by
   7  * the Free Software Foundation; either version 3 of the License, or (at
   8  * your option) any later version.
   9  *
  10  * This program is distributed in the hope that it will be useful, but
  11  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13  * General Public License for more details.
  14  *
  15  * You should have received a copy of the GNU General Public License
  16  * along with this program; if not, write to the Free Software
  17  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  18  */
  19
  20 /* Modified by Tuomo Keskitalo to add Nelder Mead Simplex test suite */
  21
  22 #include <config.h>
  23 #include <stdlib.h>
  24 #include <gsl/gsl_test.h>
  25 #include <gsl/gsl_blas.h>
  26 #include <gsl/gsl_multimin.h>
  27 #include <gsl/gsl_ieee_utils.h>
  28
  29 #include "test_funcs.h"
  30
  31 unsigned int fcount, gcount;
  32
  33 int
  34 test_fdf(const char * desc, gsl_multimin_function_fdf *f,
  35          initpt_function initpt, const gsl_multimin_fdfminimizer_type *T);
  36
  37 int
  38 test_f(const char * desc, gsl_multimin_function *f, initpt_function initpt);
  39
  40 int
  41 main (void)
  42 {
  43   const gsl_multimin_fdfminimizer_type *fdfminimizers[6];
  44   const gsl_multimin_fdfminimizer_type ** T;
  45
  46   gsl_ieee_env_setup ();
  47
  48   fdfminimizers[0] = gsl_multimin_fdfminimizer_steepest_descent;
  49   fdfminimizers[1] = gsl_multimin_fdfminimizer_conjugate_pr;
  50   fdfminimizers[2] = gsl_multimin_fdfminimizer_conjugate_fr;
  51   fdfminimizers[3] = gsl_multimin_fdfminimizer_vector_bfgs;
  52   fdfminimizers[4] = gsl_multimin_fdfminimizer_vector_bfgs2;
  53   fdfminimizers[5] = 0;
  54
  55   T = fdfminimizers;
  56
  57   while (*T != 0)
  58     {
  59       test_fdf("Roth", &roth, roth_initpt,*T);
  60       test_fdf("Wood", &wood, wood_initpt,*T);
  61       test_fdf("Rosenbrock", &rosenbrock, rosenbrock_initpt,*T);
  62       T++;
  63     }
  64
  65   test_f("Roth", &roth_fmin, roth_initpt);
  66   test_f("Wood", &wood_fmin, wood_initpt);
  67   test_f("Rosenbrock", &rosenbrock_fmin, rosenbrock_initpt);
  68
  69   T = fdfminimizers;
  70
  71   while (*T != 0)
  72     {
  73       test_fdf("NRoth", &Nroth, roth_initpt,*T);
  74       test_fdf("NWood", &Nwood, wood_initpt,*T);
  75       test_fdf("NRosenbrock", &Nrosenbrock, rosenbrock_initpt,*T);
  76       T++;
  77     }
  78
  79   exit (gsl_test_summary());
  80 }
  81
  82 int
  83 test_fdf(const char * desc,
  84          gsl_multimin_function_fdf *f,
  85          initpt_function initpt,
  86          const gsl_multimin_fdfminimizer_type *T)
  87 {
  88   int status;
  89   size_t iter = 0;
  90   double step_size;
  91
  92   gsl_vector *x = gsl_vector_alloc (f->n);
  93
  94   gsl_multimin_fdfminimizer *s;
  95   fcount = 0; gcount = 0;
  96
  97   (*initpt) (x);
  98
  99   step_size = 0.1 * gsl_blas_dnrm2 (x);
 100
 101   s = gsl_multimin_fdfminimizer_alloc(T, f->n);
 102
 103   gsl_multimin_fdfminimizer_set (s, f, x, step_size, 0.1);
 104
 105 #ifdef DEBUG
 106   printf("x "); gsl_vector_fprintf (stdout, s->x, "%g");
 107   printf("g "); gsl_vector_fprintf (stdout, s->gradient, "%g");
 108 #endif
 109
 110   do
 111     {
 112       iter++;
 113       status = gsl_multimin_fdfminimizer_iterate(s);
 114
 115 #ifdef DEBUG
 116       printf("%i: \n",iter);
 117       printf("x "); gsl_vector_fprintf (stdout, s->x, "%g");
 118       printf("g "); gsl_vector_fprintf (stdout, s->gradient, "%g");
 119       printf("f(x) %g\n",s->f);
 120       printf("dx %g\n",gsl_blas_dnrm2(s->dx));
 121       printf("\n");
 122 #endif
 123
 124       status = gsl_multimin_test_gradient(s->gradient,1e-3);
 125     }
 126   while (iter < 5000 && status == GSL_CONTINUE);
 127
 128   status |= (fabs(s->f) > 1e-5);
 129
 130   gsl_test(status, "%s, on %s: %i iters (fn+g=%d+%d), f(x)=%g",
 131            gsl_multimin_fdfminimizer_name(s),desc, iter, fcount, gcount, s->f);
 132
 133   gsl_multimin_fdfminimizer_free(s);
 134   gsl_vector_free(x);
 135
 136   return status;
 137 }
 138
 139 int
 140 test_f(const char * desc, gsl_multimin_function *f, initpt_function initpt)
 141 {
 142   /* currently this function tests only nmsimplex */
 143
 144   int status;
 145   size_t i, iter = 0;
 146   const gsl_multimin_fminimizer_type *T = gsl_multimin_fminimizer_nmsimplex;
 147
 148   gsl_vector *x = gsl_vector_alloc (f->n);
 149
 150   gsl_vector *step_size = gsl_vector_alloc (f->n);
 151
 152   gsl_multimin_fminimizer *s;
 153
 154   fcount = 0; gcount = 0;
 155   (*initpt) (x);
 156
 157   for (i = 0; i < f->n; i++)
 158     gsl_vector_set (step_size, i, 1);
 159
 160   s = gsl_multimin_fminimizer_alloc(T, f->n);
 161
 162   gsl_multimin_fminimizer_set (s, f, x, step_size);
 163
 164 #ifdef DEBUG
 165   printf("x "); gsl_vector_fprintf (stdout, s->x, "%g");
 166 #endif
 167
 168   do
 169     {
 170       iter++;
 171       status = gsl_multimin_fminimizer_iterate(s);
 172
 173 #ifdef DEBUG
 174       printf("%i: \n",iter);
 175       printf("x "); gsl_vector_fprintf (stdout, s->x, "%g");
 176       printf("f(x) %g\n", gsl_multimin_fminimizer_minimum (s));
 177       printf("size: %g\n", gsl_multimin_fminimizer_size (s));
 178       printf("\n");
 179 #endif
 180
 181       status = gsl_multimin_test_size (gsl_multimin_fminimizer_size (s),
 182                                        1e-3);
 183     }
 184   while (iter < 5000 && status == GSL_CONTINUE);
 185
 186   status |= (fabs(s->fval) > 1e-5);
 187
 188   gsl_test(status, "%s, on %s: %d iter (fn=%d), f(x)=%g",
 189            gsl_multimin_fminimizer_name(s),desc, iter, fcount, s->fval);
 190
 191   gsl_multimin_fminimizer_free(s);
 192   gsl_vector_free(x);
 193   gsl_vector_free(step_size);
 194
 195   return status;
 196 }