/* monte/test.c
*
* Copyright (C) 1996, 1997, 1998, 1999, 2000 Michael Booth
* Copyright (C) 2009 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 3 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 STEP
#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 fstep (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 Fs = make_function(&fstep, 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 STEP
/* variance(Fs) = 0.4/sqrt(1000) */
add(problems,&n, &Fs, xl, xu, 1, 100000, 0.8, 1.264e-3, "step, 1d");
add(problems,&n, &Fs, xl, xu, 2, 100000, 0.8, 1.264e-3, "step, 2d");
add(problems,&n, &Fs, xl, xu, 3, 100000, 0.8, 1.264e-3, "step, 3d");
add(problems,&n, &Fs, xl, xu, 4, 100000, 0.8, 1.264e-3, "step, 4d");
add(problems,&n, &Fs, xl, xu, 5, 100000, 0.8, 1.264e-3, "step, 5d");
add(problems,&n, &Fs, xl, xu, 6, 100000, 0.8, 1.264e-3, "step, 6d");
add(problems,&n, &Fs, xl, xu, 7, 100000, 0.8, 1.264e-3, "step, 7d");
add(problems,&n, &Fs, xl, xu, 8, 100000, 0.8, 1.264e-3, "step, 8d");
add(problems,&n, &Fs, xl, xu, 9, 100000, 0.8, 1.264e-3, "step, 9d");
add(problems,&n, &Fs, xl, xu, 10, 100000, 0.8, 1.264e-3, "step, 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 MISER
#define NAME "miser(params)"
#define MONTE_STATE gsl_monte_miser_state
#define MONTE_ALLOC gsl_monte_miser_alloc
#define MONTE_PARAMS gsl_monte_miser_params
#define MONTE_INTEGRATE(f,xl,xu,dim,calls,r,s,res,err) { gsl_monte_miser_params_get(s, ¶ms) ; params.alpha = 1.5 ; gsl_monte_miser_params_set(s, ¶ms) ; gsl_monte_miser_integrate(f,xl,xu,dim,calls,r,s,res,err); }
#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_PARAMS
#undef MONTE_INTEGRATE
#undef MONTE_FREE
#undef MONTE_ERROR_TEST
#undef MONTE_SPEEDUP
#endif
#ifdef MISER
#define NAME "miser(params)"
#define MONTE_STATE gsl_monte_miser_state
#define MONTE_ALLOC gsl_monte_miser_alloc
#define MONTE_PARAMS gsl_monte_miser_params
#define MONTE_INTEGRATE(f,xl,xu,dim,calls,r,s,res,err) { gsl_monte_miser_params_get(s, ¶ms) ; params.alpha = 1.5 ; gsl_monte_miser_params_set(s, ¶ms) ; gsl_monte_miser_integrate(f,xl,xu,dim,calls,r,s,res,err); }
#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_PARAMS
#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) ; }
#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) ; gsl_test(gsl_monte_vegas_chisq(s) < 0, NAME " returns valid chisq (%g)", gsl_monte_vegas_chisq(s))
#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(warm)"
#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) ; 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); gsl_test(gsl_monte_vegas_chisq(s) < 0, NAME " returns valid chisq (%g)", gsl_monte_vegas_chisq(s))
#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(params)"
#define MONTE_STATE gsl_monte_vegas_state
#define MONTE_ALLOC gsl_monte_vegas_alloc
#define MONTE_PARAMS gsl_monte_vegas_params
#define MONTE_INTEGRATE(f,xl,xu,dim,calls,r,s,res,err) { gsl_monte_vegas_params_get(s, ¶ms) ; params.alpha = 2 ; params.iterations = 3 ; gsl_monte_vegas_params_set(s, ¶ms) ; 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); gsl_test(gsl_monte_vegas_chisq(s) < 0, NAME " returns valid chisq (%g)", gsl_monte_vegas_chisq(s))
#include "test_main.c"
#undef NAME
#undef MONTE_STATE
#undef MONTE_ALLOC
#undef MONTE_PARAMS
#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;
}
/* Step-type (pulse) function */
double
fstep (double x[], size_t num_dim, void *params)
{
return (x[0] > 0.1 && x[0] < 0.9) ? 1 : 0;
}
/* 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);
}