| /* monte/test.c |
| * |
| * Copyright (C) 1996, 1997, 1998, 1999, 2000 Michael Booth |
| * |
| * This program is free software; you can redistribute it and/or modify |
| * it under the terms of the GNU General Public License as published by |
| * the Free Software Foundation; either version 2 of the License, or (at |
| * your option) any later version. |
| * |
| * This program is distributed in the hope that it will be useful, but |
| * WITHOUT ANY WARRANTY; without even the implied warranty of |
| * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
| * General Public License for more details. |
| * |
| * You should have received a copy of the GNU General Public License |
| * along with this program; if not, write to the Free Software |
| * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. |
| */ |
| |
| #include <config.h> |
| #include <stdlib.h> |
| #include <math.h> |
| #include <stdio.h> |
| |
| #include <gsl/gsl_math.h> |
| #include <gsl/gsl_errno.h> |
| #include <gsl/gsl_rng.h> |
| #include <gsl/gsl_test.h> |
| #include <gsl/gsl_ieee_utils.h> |
| |
| #include <gsl/gsl_rng.h> |
| #include <gsl/gsl_monte_plain.h> |
| #include <gsl/gsl_monte_miser.h> |
| #include <gsl/gsl_monte_vegas.h> |
| |
| #define CONSTANT |
| #define PRODUCT |
| #define GAUSSIAN |
| #define DBLGAUSSIAN |
| #define TSUDA |
| |
| #define PLAIN |
| #define MISER |
| #define VEGAS |
| |
| double xl[11] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }; |
| double xu[11] = { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 }; |
| double xu2[11] = { 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2 }; |
| double xu3[2] = { GSL_DBL_MAX, GSL_DBL_MAX }; |
| |
| double fconst (double x[], size_t d, void *params); |
| double f0 (double x[], size_t d, void *params); |
| double f1 (double x[], size_t d, void *params); |
| double f2 (double x[], size_t d, void *params); |
| double f3 (double x[], size_t d, void *params); |
| |
| void my_error_handler (const char *reason, const char *file, |
| int line, int err); |
| |
| struct problem { |
| gsl_monte_function * f; |
| double * xl; |
| double * xu; |
| size_t dim; |
| size_t calls; |
| double expected_result; |
| double expected_error; |
| char * description; |
| } ; |
| |
| gsl_monte_function |
| make_function (double (*f)(double *, size_t, void *), size_t d, void * p); |
| |
| gsl_monte_function |
| make_function (double (*f)(double *, size_t, void *), size_t d, void * p) |
| { |
| gsl_monte_function f_new; |
| |
| f_new.f = f; |
| f_new.dim = d; |
| f_new.params = p; |
| |
| return f_new; |
| } |
| |
| |
| void |
| add (struct problem * problems, int * n, |
| gsl_monte_function * f, double xl[], double xu[], size_t dim, size_t calls, |
| double result, double err, char * description); |
| |
| void |
| add (struct problem * problems, int * n, |
| gsl_monte_function * f, double xl[], double xu[], size_t dim, size_t calls, |
| double result, double err, char * description) |
| { |
| int i = *n; |
| problems[i].f = f; |
| problems[i].xl = xl; |
| problems[i].xu = xu; |
| problems[i].dim = dim; |
| problems[i].calls = calls; |
| problems[i].expected_result = result; |
| problems[i].expected_error = err; |
| problems[i].description = description; |
| (*n)++; |
| } |
| |
| #define TRIALS 10 |
| |
| int |
| main (void) |
| { |
| double result[TRIALS], error[TRIALS]; |
| double a = 0.1; |
| double c = (1.0 + sqrt (10.0)) / 9.0; |
| |
| gsl_monte_function Fc = make_function(&fconst, 0, 0); |
| gsl_monte_function F0 = make_function(&f0, 0, &a); |
| gsl_monte_function F1 = make_function(&f1, 0, &a); |
| gsl_monte_function F2 = make_function(&f2, 0, &a); |
| gsl_monte_function F3 = make_function(&f3, 0, &c); |
| |
| /* The relationship between the variance of the function itself, the |
| error on the integral and the number of calls is, |
| |
| sigma = sqrt(variance/N) |
| |
| where the variance is the <(f - <f>)^2> where <.> denotes the |
| volume average (integral over the integration region divided by |
| the volume) */ |
| |
| int n = 0; |
| struct problem * I; |
| struct problem problems[256]; |
| |
| #ifdef CONSTANT |
| /* variance(Fc) = 0 */ |
| |
| add(problems,&n, &Fc, xl, xu, 1, 1000, 1.0, 0.0, "constant, 1d"); |
| add(problems,&n, &Fc, xl, xu, 2, 1000, 1.0, 0.0, "constant, 2d"); |
| add(problems,&n, &Fc, xl, xu, 3, 1000, 1.0, 0.0, "constant, 3d"); |
| add(problems,&n, &Fc, xl, xu, 4, 1000, 1.0, 0.0, "constant, 4d"); |
| add(problems,&n, &Fc, xl, xu, 5, 1000, 1.0, 0.0, "constant, 5d"); |
| add(problems,&n, &Fc, xl, xu, 6, 1000, 1.0, 0.0, "constant, 6d"); |
| add(problems,&n, &Fc, xl, xu, 7, 1000, 1.0, 0.0, "constant, 7d"); |
| add(problems,&n, &Fc, xl, xu, 8, 1000, 1.0, 0.0, "constant, 8d"); |
| add(problems,&n, &Fc, xl, xu, 9, 1000, 1.0, 0.0, "constant, 9d"); |
| add(problems,&n, &Fc, xl, xu, 10, 1000, 1.0, 0.0, "constant, 10d"); |
| #endif |
| |
| #ifdef PRODUCT |
| /* variance(F0) = (4/3)^d - 1 */ |
| |
| add(problems,&n, &F0, xl, xu, 1, 3333, 1.0, 0.01, "product, 1d" ); |
| add(problems,&n, &F0, xl, xu, 2, 7777, 1.0, 0.01, "product, 2d" ); |
| add(problems,&n, &F0, xl, xu, 3, 13703, 1.0, 0.01, "product, 3d" ); |
| add(problems,&n, &F0, xl, xu, 4, 21604, 1.0, 0.01, "product, 4d" ); |
| add(problems,&n, &F0, xl, xu, 5, 32139, 1.0, 0.01, "product, 5d" ); |
| add(problems,&n, &F0, xl, xu, 6, 46186, 1.0, 0.01, "product, 6d" ); |
| add(problems,&n, &F0, xl, xu, 7, 64915, 1.0, 0.01, "product, 7d" ); |
| add(problems,&n, &F0, xl, xu, 8, 89887, 1.0, 0.01, "product, 8d" ); |
| add(problems,&n, &F0, xl, xu, 9, 123182, 1.0, 0.01, "product, 9d" ); |
| add(problems,&n, &F0, xl, xu, 10, 167577, 1.0, 0.01, "product, 10d" ); |
| #endif |
| |
| #ifdef GAUSSIAN |
| /* variance(F1) = (1/(a sqrt(2 pi)))^d - 1 */ |
| |
| add(problems,&n, &F1, xl, xu, 1, 298, 1.0, 0.1, "gaussian, 1d" ); |
| add(problems,&n, &F1, xl, xu, 2, 1492, 1.0, 0.1, "gaussian, 2d" ); |
| add(problems,&n, &F1, xl, xu, 3, 6249, 1.0, 0.1, "gaussian, 3d" ); |
| add(problems,&n, &F1, xl, xu, 4, 25230, 1.0, 0.1, "gaussian, 4d" ); |
| add(problems,&n, &F1, xl, xu, 5, 100953, 1.0, 0.1, "gaussian, 5d" ); |
| add(problems,&n, &F1, xl, xu, 6, 44782, 1.0, 0.3, "gaussian, 6d" ); |
| add(problems,&n, &F1, xl, xu, 7, 178690, 1.0, 0.3, "gaussian, 7d" ); |
| add(problems,&n, &F1, xl, xu, 8, 712904, 1.0, 0.3, "gaussian, 8d" ); |
| add(problems,&n, &F1, xl, xu, 9, 2844109, 1.0, 0.3, "gaussian, 9d" ); |
| add(problems,&n, &F1, xl, xu, 10, 11346390, 1.0, 0.3, "gaussian, 10d" ); |
| #endif |
| |
| #ifdef DBLGAUSSIAN |
| /* variance(F2) = 0.5 * (1/(a sqrt(2 pi)))^d - 1 */ |
| |
| add(problems,&n, &F2, xl, xu, 1, 9947, 1.0, 0.01, "double gaussian, 1d" ); |
| add(problems,&n, &F2, xl, xu, 2, 695, 1.0, 0.1, "double gaussian, 2d" ); |
| add(problems,&n, &F2, xl, xu, 3, 3074, 1.0, 0.1, "double gaussian, 3d" ); |
| add(problems,&n, &F2, xl, xu, 4, 12565, 1.0, 0.1, "double gaussian, 4d" ); |
| add(problems,&n, &F2, xl, xu, 5, 50426, 1.0, 0.1, "double gaussian, 5d" ); |
| add(problems,&n, &F2, xl, xu, 6, 201472, 1.0, 0.1, "double gaussian, 6d" ); |
| add(problems,&n, &F2, xl, xu, 7, 804056, 1.0, 0.1, "double gaussian, 7d" ); |
| add(problems,&n, &F2, xl, xu, 8, 356446, 1.0, 0.3, "double gaussian, 8d" ); |
| add(problems,&n, &F2, xl, xu, 9, 1422049, 1.0, 0.3, "double gaussian, 9d" ); |
| add(problems,&n, &F2, xl, xu, 10, 5673189, 1.0, 0.3, "double gaussian, 10d" ); |
| #endif |
| |
| #ifdef TSUDA |
| /* variance(F3) = ((c^2 + c + 1/3)/(c(c+1)))^d - 1 */ |
| |
| add(problems,&n, &F3, xl, xu, 1, 4928, 1.0, 0.01, "tsuda function, 1d" ); |
| add(problems,&n, &F3, xl, xu, 2, 12285, 1.0, 0.01, "tsuda function, 2d" ); |
| add(problems,&n, &F3, xl, xu, 3, 23268, 1.0, 0.01, "tsuda function, 3d" ); |
| add(problems,&n, &F3, xl, xu, 4, 39664, 1.0, 0.01, "tsuda function, 4d" ); |
| add(problems,&n, &F3, xl, xu, 5, 64141, 1.0, 0.01, "tsuda function, 5d" ); |
| add(problems,&n, &F3, xl, xu, 6, 100680, 1.0, 0.01, "tsuda function, 6d" ); |
| add(problems,&n, &F3, xl, xu, 7, 155227, 1.0, 0.01, "tsuda function, 7d" ); |
| add(problems,&n, &F3, xl, xu, 8, 236657, 1.0, 0.01, "tsuda function, 8d" ); |
| add(problems,&n, &F3, xl, xu, 9, 358219, 1.0, 0.01, "tsuda function, 9d" ); |
| add(problems,&n, &F3, xl, xu, 10, 539690, 1.0, 0.01, "tsuda function, 10d" ); |
| #endif |
| |
| add(problems,&n, 0, 0, 0, 0, 0, 0, 0, 0 ); |
| |
| |
| /* gsl_set_error_handler (&my_error_handler); */ |
| gsl_ieee_env_setup (); |
| gsl_rng_env_setup (); |
| |
| #ifdef A |
| printf ("testing allocation/input checks\n"); |
| |
| status = gsl_monte_plain_validate (s, xl, xu3, 1, 1); |
| gsl_test (status != 0, "error if limits too large"); |
| status = gsl_monte_plain_validate (s, xl, xu, 0, 10); |
| gsl_test (status != 0, "error if num_dim = 0"); |
| status = gsl_monte_plain_validate (s, xl, xu, 1, 0); |
| gsl_test (status != 0, "error if calls = 0"); |
| status = gsl_monte_plain_validate (s, xu, xl, 1, 10); |
| gsl_test (status != 0, "error if xu < xl"); |
| #endif |
| |
| #ifdef PLAIN |
| #define NAME "plain" |
| #define MONTE_STATE gsl_monte_plain_state |
| #define MONTE_ALLOC gsl_monte_plain_alloc |
| #define MONTE_INTEGRATE gsl_monte_plain_integrate |
| #define MONTE_FREE gsl_monte_plain_free |
| #define MONTE_SPEEDUP 1 |
| #define MONTE_ERROR_TEST(err,expected) gsl_test_factor(err,expected, 5.0, NAME ", %s, abserr[%d]", I->description, i) |
| #include "test_main.c" |
| #undef NAME |
| #undef MONTE_STATE |
| #undef MONTE_ALLOC |
| #undef MONTE_INTEGRATE |
| #undef MONTE_FREE |
| #undef MONTE_ERROR_TEST |
| #undef MONTE_SPEEDUP |
| #endif |
| |
| #ifdef MISER |
| #define NAME "miser" |
| #define MONTE_STATE gsl_monte_miser_state |
| #define MONTE_ALLOC gsl_monte_miser_alloc |
| #define MONTE_INTEGRATE gsl_monte_miser_integrate |
| #define MONTE_FREE gsl_monte_miser_free |
| #define MONTE_SPEEDUP 2 |
| #define MONTE_ERROR_TEST(err,expected) gsl_test(err > 5.0 * expected, NAME ", %s, abserr[%d] (obs %g vs plain %g)", I->description, i, err, expected) |
| #include "test_main.c" |
| #undef NAME |
| #undef MONTE_STATE |
| #undef MONTE_ALLOC |
| #undef MONTE_INTEGRATE |
| #undef MONTE_FREE |
| #undef MONTE_ERROR_TEST |
| #undef MONTE_SPEEDUP |
| #endif |
| |
| #ifdef VEGAS |
| #define NAME "vegas" |
| #define MONTE_STATE gsl_monte_vegas_state |
| #define MONTE_ALLOC gsl_monte_vegas_alloc |
| #define MONTE_INTEGRATE(f,xl,xu,dim,calls,r,s,res,err) { gsl_monte_vegas_integrate(f,xl,xu,dim,calls,r,s,res,err) ; if (s->chisq < 0.5 || s->chisq > 2) gsl_monte_vegas_integrate(f,xl,xu,dim,calls,r,s,res,err); } |
| #define MONTE_FREE gsl_monte_vegas_free |
| #define MONTE_SPEEDUP 3 |
| #define MONTE_ERROR_TEST(err,expected) gsl_test(err > 3.0 * (expected == 0 ? 1.0/(I->calls/MONTE_SPEEDUP) : expected), NAME ", %s, abserr[%d] (obs %g vs exp %g)", I->description, i, err, expected) |
| #include "test_main.c" |
| #undef NAME |
| #undef MONTE_STATE |
| #undef MONTE_ALLOC |
| #undef MONTE_INTEGRATE |
| #undef MONTE_FREE |
| #undef MONTE_ERROR_TEST |
| #undef MONTE_SPEEDUP |
| #endif |
| |
| exit (gsl_test_summary ()); |
| } |
| |
| /* Simple constant function */ |
| double |
| fconst (double x[], size_t num_dim, void *params) |
| { |
| return 1; |
| } |
| |
| /* Simple product function */ |
| double |
| f0 (double x[], size_t num_dim, void *params) |
| { |
| double prod = 1.0; |
| unsigned int i; |
| |
| for (i = 0; i < num_dim; ++i) |
| { |
| prod *= 2.0 * x[i]; |
| } |
| |
| return prod; |
| } |
| |
| /* Gaussian centered at 1/2. */ |
| |
| double |
| f1 (double x[], size_t num_dim, void *params) |
| { |
| double a = *(double *)params; |
| double sum = 0.; |
| |
| unsigned int i; |
| for (i = 0; i < num_dim; i++) |
| { |
| double dx = x[i] - 0.5; |
| sum += dx * dx; |
| } |
| return (pow (M_2_SQRTPI / (2. * a), (double) num_dim) * |
| exp (-sum / (a * a))); |
| } |
| |
| /* double gaussian */ |
| double |
| f2 (double x[], size_t num_dim, void *params) |
| { |
| double a = *(double *)params; |
| double sum1 = 0.; |
| double sum2 = 0.; |
| |
| unsigned int i; |
| for (i = 0; i < num_dim; i++) |
| { |
| double dx1 = x[i] - 1. / 3.; |
| double dx2 = x[i] - 2. / 3.; |
| sum1 += dx1 * dx1; |
| sum2 += dx2 * dx2; |
| } |
| return 0.5 * pow (M_2_SQRTPI / (2. * a), num_dim) |
| * (exp (-sum1 / (a * a)) + exp (-sum2 / (a * a))); |
| } |
| |
| /* Tsuda's example */ |
| double |
| f3 (double x[], size_t num_dim, void *params) |
| { |
| double c = *(double *)params; |
| |
| double prod = 1.; |
| |
| unsigned int i; |
| |
| for (i = 0; i < num_dim; i++) |
| { |
| prod *= c / (c + 1) * pow((c + 1) / (c + x[i]), 2.0); |
| } |
| |
| return prod; |
| } |
| |
| |
| void |
| my_error_handler (const char *reason, const char *file, int line, int err) |
| { |
| if (0) |
| printf ("(caught [%s:%d: %s (%d)])\n", file, line, reason, err); |
| } |