/* 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, &params) ; params.alpha = 1.5 ; gsl_monte_miser_params_set(s, &params) ; 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, &params) ; params.alpha = 1.5 ; gsl_monte_miser_params_set(s, &params) ; 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, &params) ; params.alpha = 2 ; params.iterations = 3 ; gsl_monte_vegas_params_set(s, &params) ; 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);

}