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

   1 /* monte/test.c
   2  *
   3  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Michael Booth
   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
  20 #include <config.h>
  21 #include <stdlib.h>
  22 #include <math.h>
  23 #include <stdio.h>
  24
  25 #include <gsl/gsl_math.h>
  26 #include <gsl/gsl_errno.h>
  27 #include <gsl/gsl_rng.h>
  28 #include <gsl/gsl_test.h>
  29 #include <gsl/gsl_ieee_utils.h>
  30
  31 #include <gsl/gsl_rng.h>
  32 #include <gsl/gsl_monte_plain.h>
  33 #include <gsl/gsl_monte_miser.h>
  34 #include <gsl/gsl_monte_vegas.h>
  35
  36 #define CONSTANT
  37 #define PRODUCT
  38 #define GAUSSIAN
  39 #define DBLGAUSSIAN
  40 #define TSUDA
  41
  42 #define PLAIN
  43 #define MISER
  44 #define VEGAS
  45
  46 double xl[11]  = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
  47 double xu[11]  = { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 };
  48 double xu2[11] = { 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2 };
  49 double xu3[2]  = { GSL_DBL_MAX, GSL_DBL_MAX };
  50
  51 double fconst (double x[], size_t d, void *params);
  52 double f0 (double x[], size_t d, void *params);
  53 double f1 (double x[], size_t d, void *params);
  54 double f2 (double x[], size_t d, void *params);
  55 double f3 (double x[], size_t d, void *params);
  56
  57 void my_error_handler (const char *reason, const char *file,
  58                        int line, int err);
  59
  60 struct problem {
  61   gsl_monte_function * f;
  62   double * xl;
  63   double * xu;
  64   size_t dim;
  65   size_t calls;
  66   double expected_result;
  67   double expected_error;
  68   char * description;
  69 } ;
  70
  71 gsl_monte_function
  72 make_function (double (*f)(double *, size_t, void *), size_t d, void * p);
  73
  74 gsl_monte_function
  75 make_function (double (*f)(double *, size_t, void *), size_t d, void * p)
  76 {
  77   gsl_monte_function f_new;
  78
  79   f_new.f = f;
  80   f_new.dim = d;
  81   f_new.params = p;
  82
  83   return f_new;
  84 }
  85
  86
  87 void
  88 add (struct problem * problems, int * n,
  89      gsl_monte_function * f, double xl[], double xu[], size_t dim, size_t calls,
  90      double result, double err, char * description);
  91
  92 void
  93 add (struct problem * problems, int * n,
  94      gsl_monte_function * f, double xl[], double xu[], size_t dim, size_t calls,
  95      double result, double err, char * description)
  96 {
  97   int i = *n;
  98   problems[i].f = f;
  99   problems[i].xl = xl;
 100   problems[i].xu = xu;
 101   problems[i].dim = dim;
 102   problems[i].calls = calls;
 103   problems[i].expected_result = result;
 104   problems[i].expected_error = err;
 105   problems[i].description = description;
 106   (*n)++;
 107 }
 108
 109 #define TRIALS 10
 110
 111 int
 112 main (void)
 113 {
 114   double result[TRIALS], error[TRIALS];
 115   double a = 0.1;
 116   double c = (1.0 + sqrt (10.0)) / 9.0;
 117
 118   gsl_monte_function Fc = make_function(&fconst, 0, 0);
 119   gsl_monte_function F0 = make_function(&f0, 0, &a);
 120   gsl_monte_function F1 = make_function(&f1, 0, &a);
 121   gsl_monte_function F2 = make_function(&f2, 0, &a);
 122   gsl_monte_function F3 = make_function(&f3, 0, &c);
 123
 124   /* The relationship between the variance of the function itself, the
 125      error on the integral and the number of calls is,
 126
 127      sigma = sqrt(variance/N)
 128
 129      where the variance is the <(f - <f>)^2> where <.> denotes the
 130      volume average (integral over the integration region divided by
 131      the volume) */
 132
 133   int n = 0;
 134   struct problem * I;
 135   struct problem problems[256];
 136
 137 #ifdef CONSTANT
 138     /* variance(Fc) = 0 */
 139
 140     add(problems,&n, &Fc, xl, xu,  1, 1000, 1.0, 0.0, "constant, 1d");
 141     add(problems,&n, &Fc, xl, xu,  2, 1000, 1.0, 0.0, "constant, 2d");
 142     add(problems,&n, &Fc, xl, xu,  3, 1000, 1.0, 0.0, "constant, 3d");
 143     add(problems,&n, &Fc, xl, xu,  4, 1000, 1.0, 0.0, "constant, 4d");
 144     add(problems,&n, &Fc, xl, xu,  5, 1000, 1.0, 0.0, "constant, 5d");
 145     add(problems,&n, &Fc, xl, xu,  6, 1000, 1.0, 0.0, "constant, 6d");
 146     add(problems,&n, &Fc, xl, xu,  7, 1000, 1.0, 0.0, "constant, 7d");
 147     add(problems,&n, &Fc, xl, xu,  8, 1000, 1.0, 0.0, "constant, 8d");
 148     add(problems,&n, &Fc, xl, xu,  9, 1000, 1.0, 0.0, "constant, 9d");
 149     add(problems,&n, &Fc, xl, xu, 10, 1000, 1.0, 0.0, "constant, 10d");
 150 #endif
 151
 152 #ifdef PRODUCT
 153     /* variance(F0) = (4/3)^d - 1 */
 154
 155     add(problems,&n, &F0, xl, xu,  1, 3333,   1.0, 0.01, "product, 1d" );
 156     add(problems,&n, &F0, xl, xu,  2, 7777,   1.0, 0.01, "product, 2d" );
 157     add(problems,&n, &F0, xl, xu,  3, 13703,  1.0, 0.01, "product, 3d" );
 158     add(problems,&n, &F0, xl, xu,  4, 21604,  1.0, 0.01, "product, 4d" );
 159     add(problems,&n, &F0, xl, xu,  5, 32139,  1.0, 0.01, "product, 5d" );
 160     add(problems,&n, &F0, xl, xu,  6, 46186,  1.0, 0.01, "product, 6d" );
 161     add(problems,&n, &F0, xl, xu,  7, 64915,  1.0, 0.01, "product, 7d" );
 162     add(problems,&n, &F0, xl, xu,  8, 89887,  1.0, 0.01, "product, 8d" );
 163     add(problems,&n, &F0, xl, xu,  9, 123182, 1.0, 0.01, "product, 9d" );
 164     add(problems,&n, &F0, xl, xu, 10, 167577, 1.0, 0.01, "product, 10d" );
 165 #endif
 166
 167 #ifdef GAUSSIAN
 168     /* variance(F1) = (1/(a sqrt(2 pi)))^d - 1 */
 169
 170     add(problems,&n, &F1, xl, xu,  1, 298,      1.0, 0.1, "gaussian, 1d" );
 171     add(problems,&n, &F1, xl, xu,  2, 1492,     1.0, 0.1, "gaussian, 2d" );
 172     add(problems,&n, &F1, xl, xu,  3, 6249,     1.0, 0.1, "gaussian, 3d" );
 173     add(problems,&n, &F1, xl, xu,  4, 25230,    1.0, 0.1, "gaussian, 4d" );
 174     add(problems,&n, &F1, xl, xu,  5, 100953,   1.0, 0.1, "gaussian, 5d" );
 175     add(problems,&n, &F1, xl, xu,  6, 44782,    1.0, 0.3, "gaussian, 6d" );
 176     add(problems,&n, &F1, xl, xu,  7, 178690,   1.0, 0.3, "gaussian, 7d" );
 177     add(problems,&n, &F1, xl, xu,  8, 712904,   1.0, 0.3, "gaussian, 8d" );
 178     add(problems,&n, &F1, xl, xu,  9, 2844109,  1.0, 0.3, "gaussian, 9d" );
 179     add(problems,&n, &F1, xl, xu, 10, 11346390, 1.0, 0.3, "gaussian, 10d" );
 180 #endif
 181
 182 #ifdef DBLGAUSSIAN
 183     /* variance(F2) = 0.5 * (1/(a sqrt(2 pi)))^d - 1 */
 184
 185     add(problems,&n, &F2, xl, xu,  1, 9947,  1.0, 0.01, "double gaussian, 1d" );
 186     add(problems,&n, &F2, xl, xu,  2, 695,   1.0, 0.1, "double gaussian, 2d" );
 187     add(problems,&n, &F2, xl, xu,  3, 3074,  1.0, 0.1, "double gaussian, 3d" );
 188     add(problems,&n, &F2, xl, xu,  4, 12565, 1.0, 0.1, "double gaussian, 4d" );
 189     add(problems,&n, &F2, xl, xu,  5, 50426, 1.0, 0.1, "double gaussian, 5d" );
 190     add(problems,&n, &F2, xl, xu,  6, 201472, 1.0, 0.1, "double gaussian, 6d" );
 191     add(problems,&n, &F2, xl, xu,  7, 804056, 1.0, 0.1, "double gaussian, 7d" );
 192     add(problems,&n, &F2, xl, xu,  8, 356446, 1.0, 0.3, "double gaussian, 8d" );
 193     add(problems,&n, &F2, xl, xu,  9, 1422049, 1.0, 0.3, "double gaussian, 9d" );
 194     add(problems,&n, &F2, xl, xu, 10, 5673189, 1.0, 0.3, "double gaussian, 10d" );
 195 #endif
 196
 197 #ifdef TSUDA
 198     /* variance(F3) = ((c^2 + c + 1/3)/(c(c+1)))^d - 1 */
 199
 200     add(problems,&n, &F3, xl, xu,  1, 4928,   1.0, 0.01, "tsuda function, 1d" );
 201     add(problems,&n, &F3, xl, xu,  2, 12285,  1.0, 0.01, "tsuda function, 2d" );
 202     add(problems,&n, &F3, xl, xu,  3, 23268,  1.0, 0.01, "tsuda function, 3d" );
 203     add(problems,&n, &F3, xl, xu,  4, 39664,  1.0, 0.01, "tsuda function, 4d" );
 204     add(problems,&n, &F3, xl, xu,  5, 64141,  1.0, 0.01, "tsuda function, 5d" );
 205     add(problems,&n, &F3, xl, xu,  6, 100680, 1.0, 0.01, "tsuda function, 6d" );
 206     add(problems,&n, &F3, xl, xu,  7, 155227, 1.0, 0.01, "tsuda function, 7d" );
 207     add(problems,&n, &F3, xl, xu,  8, 236657, 1.0, 0.01, "tsuda function, 8d" );
 208     add(problems,&n, &F3, xl, xu,  9, 358219, 1.0, 0.01, "tsuda function, 9d" );
 209     add(problems,&n, &F3, xl, xu, 10, 539690, 1.0, 0.01, "tsuda function, 10d" );
 210 #endif
 211
 212     add(problems,&n,   0,  0,  0, 0,    0,   0,   0, 0  );
 213
 214
 215   /* gsl_set_error_handler (&my_error_handler); */
 216   gsl_ieee_env_setup ();
 217   gsl_rng_env_setup ();
 218
 219 #ifdef A
 220   printf ("testing allocation/input checks\n");
 221
 222   status = gsl_monte_plain_validate (s, xl, xu3, 1, 1);
 223   gsl_test (status != 0, "error if limits too large");
 224   status = gsl_monte_plain_validate (s, xl, xu, 0, 10);
 225   gsl_test (status != 0, "error if num_dim = 0");
 226   status = gsl_monte_plain_validate (s, xl, xu, 1, 0);
 227   gsl_test (status != 0, "error if calls = 0");
 228   status = gsl_monte_plain_validate (s, xu, xl, 1, 10);
 229   gsl_test (status != 0, "error if xu < xl");
 230 #endif
 231
 232 #ifdef PLAIN
 233 #define NAME "plain"
 234 #define MONTE_STATE gsl_monte_plain_state
 235 #define MONTE_ALLOC gsl_monte_plain_alloc
 236 #define MONTE_INTEGRATE gsl_monte_plain_integrate
 237 #define MONTE_FREE gsl_monte_plain_free
 238 #define MONTE_SPEEDUP 1
 239 #define MONTE_ERROR_TEST(err,expected) gsl_test_factor(err,expected, 5.0, NAME ", %s, abserr[%d]", I->description, i)
 240 #include "test_main.c"
 241 #undef NAME
 242 #undef MONTE_STATE
 243 #undef MONTE_ALLOC
 244 #undef MONTE_INTEGRATE
 245 #undef MONTE_FREE
 246 #undef MONTE_ERROR_TEST
 247 #undef MONTE_SPEEDUP
 248 #endif
 249
 250 #ifdef MISER
 251 #define NAME "miser"
 252 #define MONTE_STATE gsl_monte_miser_state
 253 #define MONTE_ALLOC gsl_monte_miser_alloc
 254 #define MONTE_INTEGRATE gsl_monte_miser_integrate
 255 #define MONTE_FREE gsl_monte_miser_free
 256 #define MONTE_SPEEDUP 2
 257 #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)
 258 #include "test_main.c"
 259 #undef NAME
 260 #undef MONTE_STATE
 261 #undef MONTE_ALLOC
 262 #undef MONTE_INTEGRATE
 263 #undef MONTE_FREE
 264 #undef MONTE_ERROR_TEST
 265 #undef MONTE_SPEEDUP
 266 #endif
 267
 268 #ifdef VEGAS
 269 #define NAME "vegas"
 270 #define MONTE_STATE gsl_monte_vegas_state
 271 #define MONTE_ALLOC gsl_monte_vegas_alloc
 272 #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); }
 273 #define MONTE_FREE gsl_monte_vegas_free
 274 #define MONTE_SPEEDUP 3
 275 #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)
 276 #include "test_main.c"
 277 #undef NAME
 278 #undef MONTE_STATE
 279 #undef MONTE_ALLOC
 280 #undef MONTE_INTEGRATE
 281 #undef MONTE_FREE
 282 #undef MONTE_ERROR_TEST
 283 #undef MONTE_SPEEDUP
 284 #endif
 285
 286   exit (gsl_test_summary ());
 287 }
 288
 289 /* Simple constant function */
 290 double
 291 fconst (double x[], size_t num_dim, void *params)
 292 {
 293   return 1;
 294 }
 295
 296 /* Simple product function */
 297 double
 298 f0 (double x[], size_t num_dim, void *params)
 299 {
 300   double prod = 1.0;
 301   unsigned int i;
 302
 303   for (i = 0; i < num_dim; ++i)
 304     {
 305       prod *= 2.0 * x[i];
 306     }
 307
 308   return prod;
 309 }
 310
 311 /* Gaussian centered at 1/2. */
 312
 313 double
 314 f1 (double x[], size_t num_dim, void *params)
 315 {
 316   double a = *(double *)params;
 317   double sum = 0.;
 318
 319   unsigned int i;
 320   for (i = 0; i < num_dim; i++)
 321     {
 322       double dx = x[i] - 0.5;
 323       sum += dx * dx;
 324     }
 325   return (pow (M_2_SQRTPI / (2. * a), (double) num_dim) *
 326           exp (-sum / (a * a)));
 327 }
 328
 329 /* double gaussian */
 330 double
 331 f2 (double x[], size_t num_dim, void *params)
 332 {
 333   double a = *(double *)params;
 334   double sum1 = 0.;
 335   double sum2 = 0.;
 336
 337   unsigned int i;
 338   for (i = 0; i < num_dim; i++)
 339     {
 340       double dx1 = x[i] - 1. / 3.;
 341       double dx2 = x[i] - 2. / 3.;
 342       sum1 += dx1 * dx1;
 343       sum2 += dx2 * dx2;
 344     }
 345   return 0.5 * pow (M_2_SQRTPI / (2. * a), num_dim)
 346     * (exp (-sum1 / (a * a)) + exp (-sum2 / (a * a)));
 347 }
 348
 349 /* Tsuda's example */
 350 double
 351 f3 (double x[], size_t num_dim, void *params)
 352 {
 353   double c = *(double *)params;
 354
 355   double prod = 1.;
 356
 357   unsigned int i;
 358
 359   for (i = 0; i < num_dim; i++)
 360     {
 361       prod *= c / (c + 1) * pow((c + 1) / (c + x[i]), 2.0);
 362     }
 363
 364   return prod;
 365 }
 366
 367
 368 void
 369 my_error_handler (const char *reason, const char *file, int line, int err)
 370 {
 371   if (0)
 372     printf ("(caught [%s:%d: %s (%d)])\n", file, line, reason, err);
 373 }