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

   1 /* monte/miser.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 /* MISER.  Based on the algorithm described in the following article,
  21
  22    W.H. Press, G.R. Farrar, "Recursive Stratified Sampling for
  23    Multidimensional Monte Carlo Integration", Computers in Physics,
  24    v4 (1990), pp190-195.
  25
  26 */
  27
  28 /* Author: MJB */
  29 /* Modified by Brian Gough 12/2000 */
  30
  31 #include <config.h>
  32 #include <math.h>
  33 #include <stdlib.h>
  34
  35 #include <gsl/gsl_math.h>
  36 #include <gsl/gsl_errno.h>
  37 #include <gsl/gsl_rng.h>
  38 #include <gsl/gsl_monte.h>
  39 #include <gsl/gsl_monte_miser.h>
  40
  41 static int
  42 estimate_corrmc (gsl_monte_function * f,
  43                  const double xl[], const double xu[],
  44                  size_t dim, size_t calls,
  45                  gsl_rng * r,
  46                  gsl_monte_miser_state * state,
  47                  double *result, double *abserr,
  48                  const double xmid[], double sigma_l[], double sigma_r[]);
  49
  50
  51 int
  52 gsl_monte_miser_integrate (gsl_monte_function * f,
  53                            const double xl[], const double xu[],
  54                            size_t dim, size_t calls,
  55                            gsl_rng * r,
  56                            gsl_monte_miser_state * state,
  57                            double *result, double *abserr)
  58 {
  59   size_t n, estimate_calls, calls_l, calls_r;
  60   const size_t min_calls = state->min_calls;
  61   size_t i;
  62   size_t i_bisect;
  63   int found_best;
  64
  65   double res_est = 0, err_est = 0;
  66   double res_r = 0, err_r = 0, res_l = 0, err_l = 0;
  67   double xbi_l, xbi_m, xbi_r, s;
  68
  69   double vol;
  70   double weight_l, weight_r;
  71
  72   double *x = state->x;
  73   double *xmid = state->xmid;
  74   double *sigma_l = state->sigma_l, *sigma_r = state->sigma_r;
  75
  76   if (dim != state->dim)
  77     {
  78       GSL_ERROR ("number of dimensions must match allocated size", GSL_EINVAL);
  79     }
  80
  81   for (i = 0; i < dim; i++)
  82     {
  83       if (xu[i] <= xl[i])
  84         {
  85           GSL_ERROR ("xu must be greater than xl", GSL_EINVAL);
  86         }
  87
  88       if (xu[i] - xl[i] > GSL_DBL_MAX)
  89         {
  90           GSL_ERROR ("Range of integration is too large, please rescale",
  91                      GSL_EINVAL);
  92         }
  93     }
  94
  95   if (state->alpha < 0)
  96     {
  97       GSL_ERROR ("alpha must be non-negative", GSL_EINVAL);
  98     }
  99
 100   /* Compute volume */
 101
 102   vol = 1;
 103
 104   for (i = 0; i < dim; i++)
 105     {
 106       vol *= xu[i] - xl[i];
 107     }
 108
 109   if (calls < state->min_calls_per_bisection)
 110     {
 111       double m = 0.0, q = 0.0;
 112
 113       if (calls < 2)
 114         {
 115           GSL_ERROR ("insufficient calls for subvolume", GSL_EFAILED);
 116         }
 117
 118       for (n = 0; n < calls; n++)
 119         {
 120           /* Choose a random point in the integration region */
 121
 122           for (i = 0; i < dim; i++)
 123             {
 124               x[i] = xl[i] + gsl_rng_uniform_pos (r) * (xu[i] - xl[i]);
 125             }
 126
 127           {
 128             double fval = GSL_MONTE_FN_EVAL (f, x);
 129
 130             /* recurrence for mean and variance */
 131
 132             double d = fval - m;
 133             m += d / (n + 1.0);
 134             q += d * d * (n / (n + 1.0));
 135           }
 136         }
 137
 138       *result = vol * m;
 139
 140       *abserr = vol * sqrt (q / (calls * (calls - 1.0)));
 141
 142       return GSL_SUCCESS;
 143     }
 144
 145   estimate_calls = GSL_MAX (min_calls, calls * (state->estimate_frac));
 146
 147   if (estimate_calls < 4 * dim)
 148     {
 149       GSL_ERROR ("insufficient calls to sample all halfspaces", GSL_ESANITY);
 150     }
 151
 152   /* Flip coins to bisect the integration region with some fuzz */
 153
 154   for (i = 0; i < dim; i++)
 155     {
 156       s = (gsl_rng_uniform (r) - 0.5) >= 0.0 ? state->dither : -state->dither;
 157       state->xmid[i] = (0.5 + s) * xl[i] + (0.5 - s) * xu[i];
 158     }
 159
 160   /* The idea is to chose the direction to bisect based on which will
 161      give the smallest total variance.  We could (and may do so later)
 162      use MC to compute these variances.  But the NR guys simply estimate
 163      the variances by finding the min and max function values
 164      for each half-region for each bisection. */
 165
 166   estimate_corrmc (f, xl, xu, dim, estimate_calls,
 167                    r, state, &res_est, &err_est, xmid, sigma_l, sigma_r);
 168
 169   /* We have now used up some calls for the estimation */
 170
 171   calls -= estimate_calls;
 172
 173   /* Now find direction with the smallest total "variance" */
 174
 175   {
 176     double best_var = GSL_DBL_MAX;
 177     double beta = 2.0 / (1.0 + state->alpha);
 178     found_best = 0;
 179     i_bisect = 0;
 180     weight_l = weight_r = 1.0;
 181
 182     for (i = 0; i < dim; i++)
 183       {
 184         if (sigma_l[i] >= 0 && sigma_r[i] >= 0)
 185           {
 186             /* estimates are okay */
 187             double var = pow (sigma_l[i], beta) + pow (sigma_r[i], beta);
 188
 189             if (var <= best_var)
 190               {
 191                 found_best = 1;
 192                 best_var = var;
 193                 i_bisect = i;
 194                 weight_l = pow (sigma_l[i], beta);
 195                 weight_r = pow (sigma_r[i], beta);
 196               }
 197           }
 198         else
 199           {
 200             if (sigma_l[i] < 0)
 201               {
 202                 GSL_ERROR ("no points in left-half space!", GSL_ESANITY);
 203               }
 204             if (sigma_r[i] < 0)
 205               {
 206                 GSL_ERROR ("no points in right-half space!", GSL_ESANITY);
 207               }
 208           }
 209       }
 210   }
 211
 212   if (!found_best)
 213     {
 214       /* All estimates were the same, so chose a direction at random */
 215
 216       i_bisect = gsl_rng_uniform_int (r, dim);
 217     }
 218
 219   xbi_l = xl[i_bisect];
 220   xbi_m = xmid[i_bisect];
 221   xbi_r = xu[i_bisect];
 222
 223   /* Get the actual fractional sizes of the two "halves", and
 224      distribute the remaining calls among them */
 225
 226   {
 227     double fraction_l = fabs ((xbi_m - xbi_l) / (xbi_r - xbi_l));
 228     double fraction_r = 1 - fraction_l;
 229
 230     double a = fraction_l * weight_l;
 231     double b = fraction_r * weight_r;
 232
 233     calls_l = min_calls + (calls - 2 * min_calls) * a / (a + b);
 234     calls_r = min_calls + (calls - 2 * min_calls) * b / (a + b);
 235   }
 236
 237   /* Compute the integral for the left hand side of the bisection */
 238
 239   /* Due to the recursive nature of the algorithm we must allocate
 240      some new memory for each recursive call */
 241
 242   {
 243     int status;
 244
 245     double *xu_tmp = (double *) malloc (dim * sizeof (double));
 246
 247     if (xu_tmp == 0)
 248       {
 249         GSL_ERROR_VAL ("out of memory for left workspace", GSL_ENOMEM, 0);
 250       }
 251
 252     for (i = 0; i < dim; i++)
 253       {
 254         xu_tmp[i] = xu[i];
 255       }
 256
 257     xu_tmp[i_bisect] = xbi_m;
 258
 259     status = gsl_monte_miser_integrate (f, xl, xu_tmp,
 260                                         dim, calls_l, r, state,
 261                                         &res_l, &err_l);
 262     free (xu_tmp);
 263
 264     if (status != GSL_SUCCESS)
 265       {
 266         return status;
 267       }
 268   }
 269
 270   /* Compute the integral for the right hand side of the bisection */
 271
 272   {
 273     int status;
 274
 275     double *xl_tmp = (double *) malloc (dim * sizeof (double));
 276
 277     if (xl_tmp == 0)
 278       {
 279         GSL_ERROR_VAL ("out of memory for right workspace", GSL_ENOMEM, 0);
 280       }
 281
 282     for (i = 0; i < dim; i++)
 283       {
 284         xl_tmp[i] = xl[i];
 285       }
 286
 287     xl_tmp[i_bisect] = xbi_m;
 288
 289     status = gsl_monte_miser_integrate (f, xl_tmp, xu,
 290                                         dim, calls_r, r, state,
 291                                         &res_r, &err_r);
 292     free (xl_tmp);
 293
 294     if (status != GSL_SUCCESS)
 295       {
 296         return status;
 297       }
 298   }
 299
 300   *result = res_l + res_r;
 301   *abserr = sqrt (err_l * err_l + err_r * err_r);
 302
 303   return GSL_SUCCESS;
 304 }
 305
 306 gsl_monte_miser_state *
 307 gsl_monte_miser_alloc (size_t dim)
 308 {
 309   gsl_monte_miser_state *s =
 310     (gsl_monte_miser_state *) malloc (sizeof (gsl_monte_miser_state));
 311
 312   if (s == 0)
 313     {
 314       GSL_ERROR_VAL ("failed to allocate space for miser state struct",
 315                      GSL_ENOMEM, 0);
 316     }
 317
 318   s->x = (double *) malloc (dim * sizeof (double));
 319
 320   if (s->x == 0)
 321     {
 322       free (s);
 323       GSL_ERROR_VAL ("failed to allocate space for x", GSL_ENOMEM, 0);
 324     }
 325
 326   s->xmid = (double *) malloc (dim * sizeof (double));
 327
 328   if (s->xmid == 0)
 329     {
 330       free (s->x);
 331       free (s);
 332       GSL_ERROR_VAL ("failed to allocate space for xmid", GSL_ENOMEM, 0);
 333     }
 334
 335   s->sigma_l = (double *) malloc (dim * sizeof (double));
 336
 337   if (s->sigma_l == 0)
 338     {
 339       free (s->xmid);
 340       free (s->x);
 341       free (s);
 342       GSL_ERROR_VAL ("failed to allocate space for sigma_l", GSL_ENOMEM, 0);
 343     }
 344
 345   s->sigma_r = (double *) malloc (dim * sizeof (double));
 346
 347   if (s->sigma_r == 0)
 348     {
 349       free (s->sigma_l);
 350       free (s->xmid);
 351       free (s->x);
 352       free (s);
 353       GSL_ERROR_VAL ("failed to allocate space for sigma_r", GSL_ENOMEM, 0);
 354     }
 355
 356   s->fmax_l = (double *) malloc (dim * sizeof (double));
 357
 358   if (s->fmax_l == 0)
 359     {
 360       free (s->sigma_r);
 361       free (s->sigma_l);
 362       free (s->xmid);
 363       free (s->x);
 364       free (s);
 365       GSL_ERROR_VAL ("failed to allocate space for fmax_l", GSL_ENOMEM, 0);
 366     }
 367
 368   s->fmax_r = (double *) malloc (dim * sizeof (double));
 369
 370   if (s->fmax_r == 0)
 371     {
 372       free (s->fmax_l);
 373       free (s->sigma_r);
 374       free (s->sigma_l);
 375       free (s->xmid);
 376       free (s->x);
 377       free (s);
 378       GSL_ERROR_VAL ("failed to allocate space for fmax_r", GSL_ENOMEM, 0);
 379     }
 380
 381   s->fmin_l = (double *) malloc (dim * sizeof (double));
 382
 383   if (s->fmin_l == 0)
 384     {
 385       free (s->fmax_r);
 386       free (s->fmax_l);
 387       free (s->sigma_r);
 388       free (s->sigma_l);
 389       free (s->xmid);
 390       free (s->x);
 391       free (s);
 392       GSL_ERROR_VAL ("failed to allocate space for fmin_l", GSL_ENOMEM, 0);
 393     }
 394
 395   s->fmin_r = (double *) malloc (dim * sizeof (double));
 396
 397   if (s->fmin_r == 0)
 398     {
 399       free (s->fmin_l);
 400       free (s->fmax_r);
 401       free (s->fmax_l);
 402       free (s->sigma_r);
 403       free (s->sigma_l);
 404       free (s->xmid);
 405       free (s->x);
 406       free (s);
 407       GSL_ERROR_VAL ("failed to allocate space for fmin_r", GSL_ENOMEM, 0);
 408     }
 409
 410   s->fsum_l = (double *) malloc (dim * sizeof (double));
 411
 412   if (s->fsum_l == 0)
 413     {
 414       free (s->fmin_r);
 415       free (s->fmin_l);
 416       free (s->fmax_r);
 417       free (s->fmax_l);
 418       free (s->sigma_r);
 419       free (s->sigma_l);
 420       free (s->xmid);
 421       free (s->x);
 422       free (s);
 423       GSL_ERROR_VAL ("failed to allocate space for fsum_l", GSL_ENOMEM, 0);
 424     }
 425
 426   s->fsum_r = (double *) malloc (dim * sizeof (double));
 427
 428   if (s->fsum_r == 0)
 429     {
 430       free (s->fsum_l);
 431       free (s->fmin_r);
 432       free (s->fmin_l);
 433       free (s->fmax_r);
 434       free (s->fmax_l);
 435       free (s->sigma_r);
 436       free (s->sigma_l);
 437       free (s->xmid);
 438       free (s->x);
 439       free (s);
 440       GSL_ERROR_VAL ("failed to allocate space for fsum_r", GSL_ENOMEM, 0);
 441     }
 442
 443   s->fsum2_l = (double *) malloc (dim * sizeof (double));
 444
 445   if (s->fsum2_l == 0)
 446     {
 447       free (s->fsum_r);
 448       free (s->fsum_l);
 449       free (s->fmin_r);
 450       free (s->fmin_l);
 451       free (s->fmax_r);
 452       free (s->fmax_l);
 453       free (s->sigma_r);
 454       free (s->sigma_l);
 455       free (s->xmid);
 456       free (s->x);
 457       free (s);
 458       GSL_ERROR_VAL ("failed to allocate space for fsum2_l", GSL_ENOMEM, 0);
 459     }
 460
 461   s->fsum2_r = (double *) malloc (dim * sizeof (double));
 462
 463   if (s->fsum2_r == 0)
 464     {
 465       free (s->fsum2_l);
 466       free (s->fsum_r);
 467       free (s->fsum_l);
 468       free (s->fmin_r);
 469       free (s->fmin_l);
 470       free (s->fmax_r);
 471       free (s->fmax_l);
 472       free (s->sigma_r);
 473       free (s->sigma_l);
 474       free (s->xmid);
 475       free (s->x);
 476       free (s);
 477       GSL_ERROR_VAL ("failed to allocate space for fsum2_r", GSL_ENOMEM, 0);
 478     }
 479
 480
 481   s->hits_r = (size_t *) malloc (dim * sizeof (size_t));
 482
 483   if (s->hits_r == 0)
 484     {
 485       free (s->fsum2_r);
 486       free (s->fsum2_l);
 487       free (s->fsum_r);
 488       free (s->fsum_l);
 489       free (s->fmin_r);
 490       free (s->fmin_l);
 491       free (s->fmax_r);
 492       free (s->fmax_l);
 493       free (s->sigma_r);
 494       free (s->sigma_l);
 495       free (s->xmid);
 496       free (s->x);
 497       free (s);
 498       GSL_ERROR_VAL ("failed to allocate space for fsum2_r", GSL_ENOMEM, 0);
 499     }
 500
 501   s->hits_l = (size_t *) malloc (dim * sizeof (size_t));
 502
 503   if (s->hits_l == 0)
 504     {
 505       free (s->hits_r);
 506       free (s->fsum2_r);
 507       free (s->fsum2_l);
 508       free (s->fsum_r);
 509       free (s->fsum_l);
 510       free (s->fmin_r);
 511       free (s->fmin_l);
 512       free (s->fmax_r);
 513       free (s->fmax_l);
 514       free (s->sigma_r);
 515       free (s->sigma_l);
 516       free (s->xmid);
 517       free (s->x);
 518       free (s);
 519       GSL_ERROR_VAL ("failed to allocate space for fsum2_r", GSL_ENOMEM, 0);
 520     }
 521
 522   s->dim = dim;
 523
 524   gsl_monte_miser_init (s);
 525
 526   return s;
 527 }
 528
 529 int
 530 gsl_monte_miser_init (gsl_monte_miser_state * s)
 531 {
 532   /* We use 8 points in each halfspace to estimate the variance. There are
 533      2*dim halfspaces. A variance estimate requires a minimum of 2 points. */
 534   s->min_calls = 16 * s->dim;
 535   s->min_calls_per_bisection = 32 * s->min_calls;
 536   s->estimate_frac = 0.1;
 537   s->alpha = 2.0;
 538   s->dither = 0.0;
 539
 540   return GSL_SUCCESS;
 541 }
 542
 543 void
 544 gsl_monte_miser_free (gsl_monte_miser_state * s)
 545 {
 546   free (s->hits_r);
 547   free (s->hits_l);
 548   free (s->fsum2_r);
 549   free (s->fsum2_l);
 550   free (s->fsum_r);
 551   free (s->fsum_l);
 552   free (s->fmin_r);
 553   free (s->fmin_l);
 554   free (s->fmax_r);
 555   free (s->fmax_l);
 556   free (s->sigma_r);
 557   free (s->sigma_l);
 558   free (s->xmid);
 559   free (s->x);
 560   free (s);
 561 }
 562
 563 static int
 564 estimate_corrmc (gsl_monte_function * f,
 565                  const double xl[], const double xu[],
 566                  size_t dim, size_t calls,
 567                  gsl_rng * r,
 568                  gsl_monte_miser_state * state,
 569                  double *result, double *abserr,
 570                  const double xmid[], double sigma_l[], double sigma_r[])
 571 {
 572   size_t i, n;
 573
 574   double *x = state->x;
 575   double *fsum_l = state->fsum_l;
 576   double *fsum_r = state->fsum_r;
 577   double *fsum2_l = state->fsum2_l;
 578   double *fsum2_r = state->fsum2_r;
 579   size_t *hits_l = state->hits_l;
 580   size_t *hits_r = state->hits_r;
 581
 582   double m = 0.0, q = 0.0;
 583   double vol = 1.0;
 584
 585   for (i = 0; i < dim; i++)
 586     {
 587       vol *= xu[i] - xl[i];
 588       hits_l[i] = hits_r[i] = 0;
 589       fsum_l[i] = fsum_r[i] = 0.0;
 590       fsum2_l[i] = fsum2_r[i] = 0.0;
 591       sigma_l[i] = sigma_r[i] = -1;
 592     }
 593
 594   for (n = 0; n < calls; n++)
 595     {
 596       double fval;
 597
 598       unsigned int j = (n/2) % dim;
 599       unsigned int side = (n % 2);
 600
 601       for (i = 0; i < dim; i++)
 602         {
 603           double z = gsl_rng_uniform_pos (r) ;
 604
 605           if (i != j)
 606             {
 607               x[i] = xl[i] + z * (xu[i] - xl[i]);
 608             }
 609           else
 610             {
 611               if (side == 0)
 612                 {
 613                   x[i] = xmid[i] + z * (xu[i] - xmid[i]);
 614                 }
 615               else
 616                 {
 617                   x[i] = xl[i] + z * (xmid[i] - xl[i]);
 618                 }
 619             }
 620         }
 621
 622       fval = GSL_MONTE_FN_EVAL (f, x);
 623
 624       /* recurrence for mean and variance */
 625       {
 626         double d = fval - m;
 627         m += d / (n + 1.0);
 628         q += d * d * (n / (n + 1.0));
 629       }
 630
 631       /* compute the variances on each side of the bisection */
 632       for (i = 0; i < dim; i++)
 633         {
 634           if (x[i] <= xmid[i])
 635             {
 636               fsum_l[i] += fval;
 637               fsum2_l[i] += fval * fval;
 638               hits_l[i]++;
 639             }
 640           else
 641             {
 642               fsum_r[i] += fval;
 643               fsum2_r[i] += fval * fval;
 644               hits_r[i]++;
 645             }
 646         }
 647     }
 648
 649   for (i = 0; i < dim; i++)
 650     {
 651       double fraction_l = (xmid[i] - xl[i]) / (xu[i] - xl[i]);
 652
 653       if (hits_l[i] > 0)
 654         {
 655           fsum_l[i] /= hits_l[i];
 656           sigma_l[i] = sqrt (fsum2_l[i] - fsum_l[i] * fsum_l[i] / hits_l[i]);
 657           sigma_l[i] *= fraction_l * vol / hits_l[i];
 658         }
 659
 660       if (hits_r[i] > 0)
 661         {
 662           fsum_r[i] /= hits_r[i];
 663           sigma_r[i] = sqrt (fsum2_r[i] - fsum_r[i] * fsum_r[i] / hits_r[i]);
 664           sigma_r[i] *= (1 - fraction_l) * vol / hits_r[i];
 665         }
 666     }
 667
 668   *result = vol * m;
 669
 670   if (calls < 2)
 671     {
 672       *abserr = GSL_POSINF;
 673     }
 674   else
 675     {
 676       *abserr = vol * sqrt (q / (calls * (calls - 1.0)));
 677     }
 678
 679   return GSL_SUCCESS;
 680 }
 681