|
Packit |
67cb25 |
/* min/test.c
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Brian Gough
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This program is free software; you can redistribute it and/or modify
|
|
Packit |
67cb25 |
* it under the terms of the GNU General Public License as published by
|
|
Packit |
67cb25 |
* the Free Software Foundation; either version 3 of the License, or (at
|
|
Packit |
67cb25 |
* your option) any later version.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This program is distributed in the hope that it will be useful, but
|
|
Packit |
67cb25 |
* WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
Packit |
67cb25 |
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
|
Packit |
67cb25 |
* General Public License for more details.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* You should have received a copy of the GNU General Public License
|
|
Packit |
67cb25 |
* along with this program; if not, write to the Free Software
|
|
Packit |
67cb25 |
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
#include <config.h>
|
|
Packit |
67cb25 |
#include <stdlib.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_math.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_min.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_errno.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_test.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_ieee_utils.h>
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include "test.h"
|
|
Packit |
67cb25 |
#include "min.h"
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* stopping parameters */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
const double EPSABS = 0.001 ;
|
|
Packit |
67cb25 |
const double EPSREL = 0.001 ;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
const unsigned int MAX_ITERATIONS = 100;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
void my_error_handler (const char *reason, const char *file,
|
|
Packit |
67cb25 |
int line, int err);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#define WITHIN_TOL(a, b, epsrel, epsabs) \
|
|
Packit |
67cb25 |
(fabs((a) - (b)) < (epsrel) * GSL_MIN(fabs(a), fabs(b)) + (epsabs))
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
main (void)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_function F_cos, F_func1, F_func2, F_func3, F_func4;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
const gsl_min_fminimizer_type * fminimizer[4] ;
|
|
Packit |
67cb25 |
const gsl_min_fminimizer_type ** T;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_ieee_env_setup ();
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
fminimizer[0] = gsl_min_fminimizer_goldensection;
|
|
Packit |
67cb25 |
fminimizer[1] = gsl_min_fminimizer_brent;
|
|
Packit |
67cb25 |
fminimizer[2] = gsl_min_fminimizer_quad_golden;
|
|
Packit |
67cb25 |
fminimizer[3] = 0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
F_cos = create_function (f_cos) ;
|
|
Packit |
67cb25 |
F_func1 = create_function (func1) ;
|
|
Packit |
67cb25 |
F_func2 = create_function (func2) ;
|
|
Packit |
67cb25 |
F_func3 = create_function (func3) ;
|
|
Packit |
67cb25 |
F_func4 = create_function (func4) ;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_set_error_handler (&my_error_handler);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (T = fminimizer ; *T != 0 ; T++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
test_f (*T, "cos(x) [0 (3) 6]", &F_cos, 0.0, 3.0, 6.0, M_PI);
|
|
Packit |
67cb25 |
test_f (*T, "x^4 - 1 [-3 (-1) 17]", &F_func1, -3.0, -1.0, 17.0, 0.0);
|
|
Packit |
67cb25 |
test_f (*T, "sqrt(|x|) [-2 (-1) 1.5]", &F_func2, -2.0, -1.0, 1.5, 0.0);
|
|
Packit |
67cb25 |
test_f (*T, "func3(x) [-2 (3) 4]", &F_func3, -2.0, 3.0, 4.0, 1.0);
|
|
Packit |
67cb25 |
test_f (*T, "func4(x) [0 (0.782) 1]", &F_func4, 0, 0.782, 1.0, 0.8);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
test_f_e (*T, "invalid range check [4, 0]", &F_cos, 4.0, 3.0, 0.0, M_PI);
|
|
Packit |
67cb25 |
test_f_e (*T, "invalid range check [1, 1]", &F_cos, 1.0, 1.0, 1.0, M_PI);
|
|
Packit |
67cb25 |
test_f_e (*T, "invalid range check [-1, 1]", &F_cos, -1.0, 0.0, 1.0, M_PI);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
test_bracket("cos(x) [1,2]",&F_cos,1.0,2.0,15);
|
|
Packit |
67cb25 |
test_bracket("sqrt(|x|) [-1,0]",&F_func2,-1.0,0.0,15);
|
|
Packit |
67cb25 |
test_bracket("sqrt(|x|) [-1,-0.6]",&F_func2,-1.0,-0.6,15);
|
|
Packit |
67cb25 |
test_bracket("sqrt(|x|) [-1,1]",&F_func2,-1.0,1.0,15);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
exit (gsl_test_summary ());
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
void
|
|
Packit |
67cb25 |
test_f (const gsl_min_fminimizer_type * T,
|
|
Packit |
67cb25 |
const char * description, gsl_function *f,
|
|
Packit |
67cb25 |
double lower_bound, double middle, double upper_bound,
|
|
Packit |
67cb25 |
double correct_minimum)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int status;
|
|
Packit |
67cb25 |
size_t iterations = 0;
|
|
Packit |
67cb25 |
double m, a, b;
|
|
Packit |
67cb25 |
double x_lower, x_upper;
|
|
Packit |
67cb25 |
gsl_min_fminimizer * s;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
x_lower = lower_bound;
|
|
Packit |
67cb25 |
x_upper = upper_bound;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
s = gsl_min_fminimizer_alloc (T) ;
|
|
Packit |
67cb25 |
gsl_min_fminimizer_set (s, f, middle, x_lower, x_upper) ;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
do
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
iterations++ ;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
status = gsl_min_fminimizer_iterate (s);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
m = gsl_min_fminimizer_x_minimum(s);
|
|
Packit |
67cb25 |
a = gsl_min_fminimizer_x_lower(s);
|
|
Packit |
67cb25 |
b = gsl_min_fminimizer_x_upper(s);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#ifdef DEBUG
|
|
Packit |
67cb25 |
printf("%.12f %.18f %.12f %.18f %.12f %.18f status=%d\n",
|
|
Packit |
67cb25 |
a, GSL_FN_EVAL(f, a), m, GSL_FN_EVAL(f, m), b, GSL_FN_EVAL(f, b), status);
|
|
Packit |
67cb25 |
#endif
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (a > b)
|
|
Packit |
67cb25 |
gsl_test (GSL_FAILURE, "interval is invalid (%g,%g)", a, b);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (m < a || m > b)
|
|
Packit |
67cb25 |
gsl_test (GSL_FAILURE, "m lies outside interval %g (%g,%g)", m, a, b);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (status) break ;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
status = gsl_min_test_interval (a, b, EPSABS, EPSREL);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
while (status == GSL_CONTINUE && iterations < MAX_ITERATIONS);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_test (status, "%s, %s (%g obs vs %g expected) ",
|
|
Packit |
67cb25 |
gsl_min_fminimizer_name(s), description,
|
|
Packit |
67cb25 |
gsl_min_fminimizer_x_minimum(s), correct_minimum);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* check the validity of the returned result */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (!WITHIN_TOL (m, correct_minimum, EPSREL, EPSABS))
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_test (GSL_FAILURE, "incorrect precision (%g obs vs %g expected)",
|
|
Packit |
67cb25 |
m, correct_minimum);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_min_fminimizer_free (s);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
void
|
|
Packit |
67cb25 |
test_f_e (const gsl_min_fminimizer_type * T,
|
|
Packit |
67cb25 |
const char * description, gsl_function *f,
|
|
Packit |
67cb25 |
double lower_bound, double middle, double upper_bound,
|
|
Packit |
67cb25 |
double correct_minimum)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int status;
|
|
Packit |
67cb25 |
size_t iterations = 0;
|
|
Packit |
67cb25 |
double x_lower, x_upper;
|
|
Packit |
67cb25 |
double a, b;
|
|
Packit |
67cb25 |
gsl_min_fminimizer * s;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
x_lower = lower_bound;
|
|
Packit |
67cb25 |
x_upper = upper_bound;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
s = gsl_min_fminimizer_alloc (T) ;
|
|
Packit |
67cb25 |
status = gsl_min_fminimizer_set (s, f, middle, x_lower, x_upper) ;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (status != GSL_SUCCESS)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_min_fminimizer_free (s) ;
|
|
Packit |
67cb25 |
gsl_test (status == GSL_SUCCESS, "%s, %s", T->name, description);
|
|
Packit |
67cb25 |
return ;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
do
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
iterations++ ;
|
|
Packit |
67cb25 |
gsl_min_fminimizer_iterate (s);
|
|
Packit |
67cb25 |
a = gsl_min_fminimizer_x_lower(s);
|
|
Packit |
67cb25 |
b = gsl_min_fminimizer_x_upper(s);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
status = gsl_min_test_interval (a, b, EPSABS, EPSREL);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
while (status == GSL_CONTINUE && iterations < MAX_ITERATIONS);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_test (!status, "%s, %s", gsl_min_fminimizer_name(s), description,
|
|
Packit |
67cb25 |
gsl_min_fminimizer_x_minimum(s) - correct_minimum);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_min_fminimizer_free (s);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
void
|
|
Packit |
67cb25 |
my_error_handler (const char *reason, const char *file, int line, int err)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (0)
|
|
Packit |
67cb25 |
printf ("(caught [%s:%d: %s (%d)])\n", file, line, reason, err);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
test_bracket (const char * description,gsl_function *f,double lower_bound,
|
|
Packit |
67cb25 |
double upper_bound, unsigned int max)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int status;
|
|
Packit |
67cb25 |
double x_lower, x_upper;
|
|
Packit |
67cb25 |
double f_upper,f_lower,f_minimum;
|
|
Packit |
67cb25 |
double x_minimum;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
x_lower=lower_bound;
|
|
Packit |
67cb25 |
x_upper=upper_bound;
|
|
Packit |
67cb25 |
SAFE_FUNC_CALL (f,x_lower,&f_lower);
|
|
Packit |
67cb25 |
SAFE_FUNC_CALL (f,x_upper,&f_upper);
|
|
Packit |
67cb25 |
status=gsl_min_find_bracket(f,&x_minimum,&f_minimum,&x_lower,&f_lower,&x_upper,&f_upper,max);
|
|
Packit |
67cb25 |
gsl_test (status,"%s, interval: [%g,%g], values: (%g,%g), minimum at: %g, value: %g",
|
|
Packit |
67cb25 |
description,x_lower,x_upper,f_lower,f_upper,x_minimum,f_minimum);
|
|
Packit |
67cb25 |
return status;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|