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

   1 /* monte/plain.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 /* Plain Monte-Carlo. */
  21
  22 /* Author: MJB */
  23
  24 #include <config.h>
  25 #include <math.h>
  26 #include <gsl/gsl_math.h>
  27 #include <gsl/gsl_rng.h>
  28 #include <gsl/gsl_monte_plain.h>
  29
  30 int
  31 gsl_monte_plain_integrate (const gsl_monte_function * f,
  32                            const double xl[], const double xu[],
  33                            const size_t dim,
  34                            const size_t calls,
  35                            gsl_rng * r,
  36                            gsl_monte_plain_state * state,
  37                            double *result, double *abserr)
  38 {
  39   double vol, m = 0, q = 0;
  40   double *x = state->x;
  41   size_t n, i;
  42
  43   if (dim != state->dim)
  44     {
  45       GSL_ERROR ("number of dimensions must match allocated size", GSL_EINVAL);
  46     }
  47
  48   for (i = 0; i < dim; i++)
  49     {
  50       if (xu[i] <= xl[i])
  51         {
  52           GSL_ERROR ("xu must be greater than xl", GSL_EINVAL);
  53         }
  54
  55       if (xu[i] - xl[i] > GSL_DBL_MAX)
  56         {
  57           GSL_ERROR ("Range of integration is too large, please rescale",
  58                      GSL_EINVAL);
  59         }
  60     }
  61
  62   /* Compute the volume of the region */
  63
  64   vol = 1;
  65
  66   for (i = 0; i < dim; i++)
  67     {
  68       vol *= xu[i] - xl[i];
  69     }
  70
  71   for (n = 0; n < calls; n++)
  72     {
  73       /* Choose a random point in the integration region */
  74
  75       for (i = 0; i < dim; i++)
  76         {
  77           x[i] = xl[i] + gsl_rng_uniform_pos (r) * (xu[i] - xl[i]);
  78         }
  79
  80       {
  81         double fval = GSL_MONTE_FN_EVAL (f, x);
  82
  83         /* recurrence for mean and variance */
  84
  85         double d = fval - m;
  86         m += d / (n + 1.0);
  87         q += d * d * (n / (n + 1.0));
  88       }
  89     }
  90
  91   *result = vol * m;
  92
  93   if (calls < 2)
  94     {
  95       *abserr = GSL_POSINF;
  96     }
  97   else
  98     {
  99       *abserr = vol * sqrt (q / (calls * (calls - 1.0)));
 100     }
 101
 102   return GSL_SUCCESS;
 103 }
 104
 105 gsl_monte_plain_state *
 106 gsl_monte_plain_alloc (size_t dim)
 107 {
 108   gsl_monte_plain_state *s =
 109     (gsl_monte_plain_state *) malloc (sizeof (gsl_monte_plain_state));
 110
 111   if (s == 0)
 112     {
 113       GSL_ERROR_VAL ("failed to allocate space for state struct",
 114                      GSL_ENOMEM, 0);
 115     }
 116
 117   s->x = (double *) malloc (dim * sizeof (double));
 118
 119   if (s->x == 0)
 120     {
 121       free (s);
 122       GSL_ERROR_VAL ("failed to allocate space for working vector",
 123                      GSL_ENOMEM, 0);
 124     }
 125
 126   s->dim = dim;
 127
 128   return s;
 129 }
 130
 131 /* Set some default values and whatever */
 132
 133 int
 134 gsl_monte_plain_init (gsl_monte_plain_state * s)
 135 {
 136   size_t i;
 137
 138   for (i = 0; i < s->dim; i++)
 139     {
 140       s->x[i] = 0.0;
 141     }
 142
 143   return GSL_SUCCESS;
 144 }
 145
 146 void
 147 gsl_monte_plain_free (gsl_monte_plain_state * s)
 148 {
 149   free (s->x);
 150   free (s);
 151 }