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

   1 /* min/test.c
   2  *
   3  * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Brian Gough
   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 #include <config.h>
  20 #include <stdlib.h>
  21 #include <gsl/gsl_math.h>
  22 #include <gsl/gsl_min.h>
  23 #include <gsl/gsl_errno.h>
  24 #include <gsl/gsl_test.h>
  25 #include <gsl/gsl_ieee_utils.h>
  26
  27 #include "test.h"
  28 #include "min.h"
  29
  30 /* stopping parameters */
  31
  32 const double EPSABS = 0.001 ;
  33 const double EPSREL = 0.001 ;
  34
  35 const unsigned int MAX_ITERATIONS = 100;
  36
  37 void my_error_handler (const char *reason, const char *file,
  38                        int line, int err);
  39
  40 #define WITHIN_TOL(a, b, epsrel, epsabs) \
  41  (fabs((a) - (b)) < (epsrel) * GSL_MIN(fabs(a), fabs(b)) + (epsabs))
  42
  43 int
  44 main (void)
  45 {
  46   gsl_function F_cos, F_func1, F_func2, F_func3, F_func4;
  47
  48   const gsl_min_fminimizer_type * fminimizer[4] ;
  49   const gsl_min_fminimizer_type ** T;
  50
  51   gsl_ieee_env_setup ();
  52
  53   fminimizer[0] = gsl_min_fminimizer_goldensection;
  54   fminimizer[1] = gsl_min_fminimizer_brent;
  55   fminimizer[2] = 0;
  56
  57   F_cos = create_function (f_cos) ;
  58   F_func1 = create_function (func1) ;
  59   F_func2 = create_function (func2) ;
  60   F_func3 = create_function (func3) ;
  61   F_func4 = create_function (func4) ;
  62
  63   gsl_set_error_handler (&my_error_handler);
  64
  65   for (T = fminimizer ; *T != 0 ; T++)
  66     {
  67       test_f (*T, "cos(x) [0 (3) 6]", &F_cos, 0.0, 3.0, 6.0, M_PI);
  68       test_f (*T, "x^4 - 1 [-3 (-1) 17]", &F_func1, -3.0, -1.0, 17.0, 0.0);
  69       test_f (*T, "sqrt(|x|) [-2 (-1) 1.5]", &F_func2, -2.0, -1.0, 1.5, 0.0);
  70       test_f (*T, "func3(x) [-2 (3) 4]", &F_func3, -2.0, 3.0, 4.0, 1.0);
  71       test_f (*T, "func4(x) [0 (0.782) 1]", &F_func4, 0, 0.782, 1.0, 0.8);
  72
  73       test_f_e (*T, "invalid range check [4, 0]", &F_cos, 4.0, 3.0, 0.0, M_PI);
  74       test_f_e (*T, "invalid range check [1, 1]", &F_cos, 1.0, 1.0, 1.0, M_PI);
  75       test_f_e (*T, "invalid range check [-1, 1]", &F_cos, -1.0, 0.0, 1.0, M_PI);
  76     }
  77   test_bracket("cos(x) [1,2]",&F_cos,1.0,2.0,15);
  78   test_bracket("sqrt(|x|) [-1,0]",&F_func2,-1.0,0.0,15);
  79   test_bracket("sqrt(|x|) [-1,-0.6]",&F_func2,-1.0,-0.6,15);
  80   test_bracket("sqrt(|x|) [-1,1]",&F_func2,-1.0,1.0,15);
  81
  82   exit (gsl_test_summary ());
  83 }
  84
  85 void
  86 test_f (const gsl_min_fminimizer_type * T,
  87         const char * description, gsl_function *f,
  88         double lower_bound, double middle, double upper_bound,
  89         double correct_minimum)
  90 {
  91   int status;
  92   size_t iterations = 0;
  93   double m, a, b;
  94   double x_lower, x_upper;
  95   gsl_min_fminimizer * s;
  96
  97   x_lower = lower_bound;
  98   x_upper = upper_bound;
  99
 100   s = gsl_min_fminimizer_alloc (T) ;
 101   gsl_min_fminimizer_set (s, f, middle, x_lower, x_upper) ;
 102
 103   do
 104     {
 105       iterations++ ;
 106
 107       status = gsl_min_fminimizer_iterate (s);
 108
 109       m = gsl_min_fminimizer_x_minimum(s);
 110       a = gsl_min_fminimizer_x_lower(s);
 111       b = gsl_min_fminimizer_x_upper(s);
 112
 113 #ifdef DEBUG
 114       printf("%.12f %.18f %.12f %.18f %.12f %.18f status=%d\n",
 115              a, GSL_FN_EVAL(f, a), m, GSL_FN_EVAL(f, m), b, GSL_FN_EVAL(f, b), status);
 116 #endif
 117
 118       if (a > b)
 119         gsl_test (GSL_FAILURE, "interval is invalid (%g,%g)", a, b);
 120
 121       if (m < a || m > b)
 122         gsl_test (GSL_FAILURE, "m lies outside interval %g (%g,%g)", m, a, b);
 123
 124       if (status) break ;
 125
 126       status = gsl_min_test_interval (a, b, EPSABS, EPSREL);
 127     }
 128   while (status == GSL_CONTINUE && iterations < MAX_ITERATIONS);
 129
 130   gsl_test (status, "%s, %s (%g obs vs %g expected) ",
 131             gsl_min_fminimizer_name(s), description,
 132             gsl_min_fminimizer_x_minimum(s), correct_minimum);
 133
 134   /* check the validity of the returned result */
 135
 136   if (!WITHIN_TOL (m, correct_minimum, EPSREL, EPSABS))
 137     {
 138       gsl_test (GSL_FAILURE, "incorrect precision (%g obs vs %g expected)",
 139                 m, correct_minimum);
 140     }
 141
 142   gsl_min_fminimizer_free (s);
 143
 144 }
 145
 146 void
 147 test_f_e (const gsl_min_fminimizer_type * T,
 148           const char * description, gsl_function *f,
 149           double lower_bound, double middle, double upper_bound,
 150           double correct_minimum)
 151 {
 152   int status;
 153   size_t iterations = 0;
 154   double x_lower, x_upper;
 155   double a, b;
 156   gsl_min_fminimizer * s;
 157
 158   x_lower = lower_bound;
 159   x_upper = upper_bound;
 160
 161   s = gsl_min_fminimizer_alloc (T) ;
 162   status = gsl_min_fminimizer_set (s, f, middle, x_lower, x_upper) ;
 163
 164   if (status != GSL_SUCCESS)
 165     {
 166       gsl_min_fminimizer_free (s) ;
 167       gsl_test (status == GSL_SUCCESS, "%s, %s", T->name, description);
 168       return ;
 169     }
 170
 171   do
 172     {
 173       iterations++ ;
 174       gsl_min_fminimizer_iterate (s);
 175       a = gsl_min_fminimizer_x_lower(s);
 176       b = gsl_min_fminimizer_x_upper(s);
 177
 178       status = gsl_min_test_interval (a, b, EPSABS, EPSREL);
 179     }
 180   while (status == GSL_CONTINUE && iterations < MAX_ITERATIONS);
 181
 182   gsl_test (!status, "%s, %s", gsl_min_fminimizer_name(s), description,
 183             gsl_min_fminimizer_x_minimum(s) - correct_minimum);
 184
 185   gsl_min_fminimizer_free (s);
 186 }
 187
 188 void
 189 my_error_handler (const char *reason, const char *file, int line, int err)
 190 {
 191   if (0)
 192     printf ("(caught [%s:%d: %s (%d)])\n", file, line, reason, err);
 193 }
 194
 195 int
 196 test_bracket (const char * description,gsl_function *f,double lower_bound,
 197               double upper_bound, unsigned int max)
 198 {
 199   int status;
 200   double x_lower, x_upper;
 201   double f_upper,f_lower,f_minimum;
 202   double x_minimum;
 203
 204   x_lower=lower_bound;
 205   x_upper=upper_bound;
 206   SAFE_FUNC_CALL (f,x_lower,&f_lower);
 207   SAFE_FUNC_CALL (f,x_upper,&f_upper);
 208   status=gsl_min_find_bracket(f,&x_minimum,&f_minimum,&x_lower,&f_lower,&x_upper,&f_upper,max);
 209   gsl_test (status,"%s, interval: [%g,%g], values: (%g,%g), minimum at: %g, value: %g",
 210             description,x_lower,x_upper,f_lower,f_upper,x_minimum,f_minimum);
 211   return status;
 212 }
 213
 214
 215