|
Packit |
67cb25 |
/* monte/miser.c
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Copyright (C) 1996, 1997, 1998, 1999, 2000 Michael Booth
|
|
Packit |
67cb25 |
* Copyright (C) 2009 Brian Gough
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This program is free software; you can redistribute it and/or modify
|
|
Packit |
67cb25 |
* it under the terms of the GNU General Public License as published by
|
|
Packit |
67cb25 |
* the Free Software Foundation; either version 3 of the License, or (at
|
|
Packit |
67cb25 |
* your option) any later version.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This program is distributed in the hope that it will be useful, but
|
|
Packit |
67cb25 |
* WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
Packit |
67cb25 |
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
|
Packit |
67cb25 |
* General Public License for more details.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* You should have received a copy of the GNU General Public License
|
|
Packit |
67cb25 |
* along with this program; if not, write to the Free Software
|
|
Packit |
67cb25 |
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* MISER. Based on the algorithm described in the following article,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
W.H. Press, G.R. Farrar, "Recursive Stratified Sampling for
|
|
Packit |
67cb25 |
Multidimensional Monte Carlo Integration", Computers in Physics,
|
|
Packit |
67cb25 |
v4 (1990), pp190-195.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Author: MJB */
|
|
Packit |
67cb25 |
/* Modified by Brian Gough 12/2000 */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include <config.h>
|
|
Packit |
67cb25 |
#include <math.h>
|
|
Packit |
67cb25 |
#include <stdlib.h>
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include <gsl/gsl_math.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_errno.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_rng.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_monte.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_monte_miser.h>
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static int
|
|
Packit |
67cb25 |
estimate_corrmc (gsl_monte_function * f,
|
|
Packit |
67cb25 |
const double xl[], const double xu[],
|
|
Packit |
67cb25 |
size_t dim, size_t calls,
|
|
Packit |
67cb25 |
gsl_rng * r,
|
|
Packit |
67cb25 |
gsl_monte_miser_state * state,
|
|
Packit |
67cb25 |
double *result, double *abserr,
|
|
Packit |
67cb25 |
const double xmid[], double sigma_l[], double sigma_r[]);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_monte_miser_integrate (gsl_monte_function * f,
|
|
Packit |
67cb25 |
const double xl[], const double xu[],
|
|
Packit |
67cb25 |
size_t dim, size_t calls,
|
|
Packit |
67cb25 |
gsl_rng * r,
|
|
Packit |
67cb25 |
gsl_monte_miser_state * state,
|
|
Packit |
67cb25 |
double *result, double *abserr)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
size_t n, estimate_calls, calls_l, calls_r;
|
|
Packit |
67cb25 |
const size_t min_calls = state->min_calls;
|
|
Packit |
67cb25 |
size_t i;
|
|
Packit |
67cb25 |
size_t i_bisect;
|
|
Packit |
67cb25 |
int found_best;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double res_est = 0, err_est = 0;
|
|
Packit |
67cb25 |
double res_r = 0, err_r = 0, res_l = 0, err_l = 0;
|
|
Packit |
67cb25 |
double xbi_l, xbi_m, xbi_r, s;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double vol;
|
|
Packit |
67cb25 |
double weight_l, weight_r;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double *x = state->x;
|
|
Packit |
67cb25 |
double *xmid = state->xmid;
|
|
Packit |
67cb25 |
double *sigma_l = state->sigma_l, *sigma_r = state->sigma_r;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (dim != state->dim)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("number of dimensions must match allocated size", GSL_EINVAL);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < dim; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (xu[i] <= xl[i])
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("xu must be greater than xl", GSL_EINVAL);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (xu[i] - xl[i] > GSL_DBL_MAX)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("Range of integration is too large, please rescale",
|
|
Packit |
67cb25 |
GSL_EINVAL);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (state->alpha < 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("alpha must be non-negative", GSL_EINVAL);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Compute volume */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
vol = 1;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < dim; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
vol *= xu[i] - xl[i];
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (calls < state->min_calls_per_bisection)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double m = 0.0, q = 0.0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (calls < 2)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("insufficient calls for subvolume", GSL_EFAILED);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (n = 0; n < calls; n++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* Choose a random point in the integration region */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < dim; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
x[i] = xl[i] + gsl_rng_uniform_pos (r) * (xu[i] - xl[i]);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double fval = GSL_MONTE_FN_EVAL (f, x);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* recurrence for mean and variance */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double d = fval - m;
|
|
Packit |
67cb25 |
m += d / (n + 1.0);
|
|
Packit |
67cb25 |
q += d * d * (n / (n + 1.0));
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
*result = vol * m;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
*abserr = vol * sqrt (q / (calls * (calls - 1.0)));
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
estimate_calls = GSL_MAX (min_calls, calls * (state->estimate_frac));
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (estimate_calls < 4 * dim)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("insufficient calls to sample all halfspaces", GSL_ESANITY);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Flip coins to bisect the integration region with some fuzz */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < dim; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
s = (gsl_rng_uniform (r) - 0.5) >= 0.0 ? state->dither : -state->dither;
|
|
Packit |
67cb25 |
state->xmid[i] = (0.5 + s) * xl[i] + (0.5 - s) * xu[i];
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* The idea is to chose the direction to bisect based on which will
|
|
Packit |
67cb25 |
give the smallest total variance. We could (and may do so later)
|
|
Packit |
67cb25 |
use MC to compute these variances. But the NR guys simply estimate
|
|
Packit |
67cb25 |
the variances by finding the min and max function values
|
|
Packit |
67cb25 |
for each half-region for each bisection. */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
estimate_corrmc (f, xl, xu, dim, estimate_calls,
|
|
Packit |
67cb25 |
r, state, &res_est, &err_est, xmid, sigma_l, sigma_r);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* We have now used up some calls for the estimation */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
calls -= estimate_calls;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Now find direction with the smallest total "variance" */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double best_var = GSL_DBL_MAX;
|
|
Packit |
67cb25 |
double beta = 2.0 / (1.0 + state->alpha);
|
|
Packit |
67cb25 |
found_best = 0;
|
|
Packit |
67cb25 |
i_bisect = 0;
|
|
Packit |
67cb25 |
weight_l = weight_r = 1.0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < dim; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (sigma_l[i] >= 0 && sigma_r[i] >= 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* estimates are okay */
|
|
Packit |
67cb25 |
double var = pow (sigma_l[i], beta) + pow (sigma_r[i], beta);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (var <= best_var)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
found_best = 1;
|
|
Packit |
67cb25 |
best_var = var;
|
|
Packit |
67cb25 |
i_bisect = i;
|
|
Packit |
67cb25 |
weight_l = pow (sigma_l[i], beta);
|
|
Packit |
67cb25 |
weight_r = pow (sigma_r[i], beta);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (weight_l == 0 && weight_r == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
weight_l = 1;
|
|
Packit |
67cb25 |
weight_r = 1;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (sigma_l[i] < 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("no points in left-half space!", GSL_ESANITY);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
if (sigma_r[i] < 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("no points in right-half space!", GSL_ESANITY);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (!found_best)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* All estimates were the same, so chose a direction at random */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
i_bisect = gsl_rng_uniform_int (r, dim);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
xbi_l = xl[i_bisect];
|
|
Packit |
67cb25 |
xbi_m = xmid[i_bisect];
|
|
Packit |
67cb25 |
xbi_r = xu[i_bisect];
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Get the actual fractional sizes of the two "halves", and
|
|
Packit |
67cb25 |
distribute the remaining calls among them */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double fraction_l = fabs ((xbi_m - xbi_l) / (xbi_r - xbi_l));
|
|
Packit |
67cb25 |
double fraction_r = 1 - fraction_l;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double a = fraction_l * weight_l;
|
|
Packit |
67cb25 |
double b = fraction_r * weight_r;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
calls_l = min_calls + (calls - 2 * min_calls) * a / (a + b);
|
|
Packit |
67cb25 |
calls_r = min_calls + (calls - 2 * min_calls) * b / (a + b);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Compute the integral for the left hand side of the bisection */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Due to the recursive nature of the algorithm we must allocate
|
|
Packit |
67cb25 |
some new memory for each recursive call */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int status;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double *xu_tmp = (double *) malloc (dim * sizeof (double));
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (xu_tmp == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR_VAL ("out of memory for left workspace", GSL_ENOMEM, 0);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < dim; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
xu_tmp[i] = xu[i];
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
xu_tmp[i_bisect] = xbi_m;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
status = gsl_monte_miser_integrate (f, xl, xu_tmp,
|
|
Packit |
67cb25 |
dim, calls_l, r, state,
|
|
Packit |
67cb25 |
&res_l, &err_l);
|
|
Packit |
67cb25 |
free (xu_tmp);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (status != GSL_SUCCESS)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
return status;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Compute the integral for the right hand side of the bisection */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int status;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double *xl_tmp = (double *) malloc (dim * sizeof (double));
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (xl_tmp == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR_VAL ("out of memory for right workspace", GSL_ENOMEM, 0);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < dim; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
xl_tmp[i] = xl[i];
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
xl_tmp[i_bisect] = xbi_m;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
status = gsl_monte_miser_integrate (f, xl_tmp, xu,
|
|
Packit |
67cb25 |
dim, calls_r, r, state,
|
|
Packit |
67cb25 |
&res_r, &err_r);
|
|
Packit |
67cb25 |
free (xl_tmp);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (status != GSL_SUCCESS)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
return status;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
*result = res_l + res_r;
|
|
Packit |
67cb25 |
*abserr = sqrt (err_l * err_l + err_r * err_r);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_monte_miser_state *
|
|
Packit |
67cb25 |
gsl_monte_miser_alloc (size_t dim)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_monte_miser_state *s =
|
|
Packit |
67cb25 |
(gsl_monte_miser_state *) malloc (sizeof (gsl_monte_miser_state));
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (s == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR_VAL ("failed to allocate space for miser state struct",
|
|
Packit |
67cb25 |
GSL_ENOMEM, 0);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
s->x = (double *) malloc (dim * sizeof (double));
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (s->x == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
free (s);
|
|
Packit |
67cb25 |
GSL_ERROR_VAL ("failed to allocate space for x", GSL_ENOMEM, 0);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
s->xmid = (double *) malloc (dim * sizeof (double));
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (s->xmid == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
free (s->x);
|
|
Packit |
67cb25 |
free (s);
|
|
Packit |
67cb25 |
GSL_ERROR_VAL ("failed to allocate space for xmid", GSL_ENOMEM, 0);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
s->sigma_l = (double *) malloc (dim * sizeof (double));
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (s->sigma_l == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
free (s->xmid);
|
|
Packit |
67cb25 |
free (s->x);
|
|
Packit |
67cb25 |
free (s);
|
|
Packit |
67cb25 |
GSL_ERROR_VAL ("failed to allocate space for sigma_l", GSL_ENOMEM, 0);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
s->sigma_r = (double *) malloc (dim * sizeof (double));
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (s->sigma_r == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
free (s->sigma_l);
|
|
Packit |
67cb25 |
free (s->xmid);
|
|
Packit |
67cb25 |
free (s->x);
|
|
Packit |
67cb25 |
free (s);
|
|
Packit |
67cb25 |
GSL_ERROR_VAL ("failed to allocate space for sigma_r", GSL_ENOMEM, 0);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
s->fmax_l = (double *) malloc (dim * sizeof (double));
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (s->fmax_l == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
free (s->sigma_r);
|
|
Packit |
67cb25 |
free (s->sigma_l);
|
|
Packit |
67cb25 |
free (s->xmid);
|
|
Packit |
67cb25 |
free (s->x);
|
|
Packit |
67cb25 |
free (s);
|
|
Packit |
67cb25 |
GSL_ERROR_VAL ("failed to allocate space for fmax_l", GSL_ENOMEM, 0);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
s->fmax_r = (double *) malloc (dim * sizeof (double));
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (s->fmax_r == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
free (s->fmax_l);
|
|
Packit |
67cb25 |
free (s->sigma_r);
|
|
Packit |
67cb25 |
free (s->sigma_l);
|
|
Packit |
67cb25 |
free (s->xmid);
|
|
Packit |
67cb25 |
free (s->x);
|
|
Packit |
67cb25 |
free (s);
|
|
Packit |
67cb25 |
GSL_ERROR_VAL ("failed to allocate space for fmax_r", GSL_ENOMEM, 0);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
s->fmin_l = (double *) malloc (dim * sizeof (double));
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (s->fmin_l == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
free (s->fmax_r);
|
|
Packit |
67cb25 |
free (s->fmax_l);
|
|
Packit |
67cb25 |
free (s->sigma_r);
|
|
Packit |
67cb25 |
free (s->sigma_l);
|
|
Packit |
67cb25 |
free (s->xmid);
|
|
Packit |
67cb25 |
free (s->x);
|
|
Packit |
67cb25 |
free (s);
|
|
Packit |
67cb25 |
GSL_ERROR_VAL ("failed to allocate space for fmin_l", GSL_ENOMEM, 0);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
s->fmin_r = (double *) malloc (dim * sizeof (double));
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (s->fmin_r == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
free (s->fmin_l);
|
|
Packit |
67cb25 |
free (s->fmax_r);
|
|
Packit |
67cb25 |
free (s->fmax_l);
|
|
Packit |
67cb25 |
free (s->sigma_r);
|
|
Packit |
67cb25 |
free (s->sigma_l);
|
|
Packit |
67cb25 |
free (s->xmid);
|
|
Packit |
67cb25 |
free (s->x);
|
|
Packit |
67cb25 |
free (s);
|
|
Packit |
67cb25 |
GSL_ERROR_VAL ("failed to allocate space for fmin_r", GSL_ENOMEM, 0);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
s->fsum_l = (double *) malloc (dim * sizeof (double));
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (s->fsum_l == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
free (s->fmin_r);
|
|
Packit |
67cb25 |
free (s->fmin_l);
|
|
Packit |
67cb25 |
free (s->fmax_r);
|
|
Packit |
67cb25 |
free (s->fmax_l);
|
|
Packit |
67cb25 |
free (s->sigma_r);
|
|
Packit |
67cb25 |
free (s->sigma_l);
|
|
Packit |
67cb25 |
free (s->xmid);
|
|
Packit |
67cb25 |
free (s->x);
|
|
Packit |
67cb25 |
free (s);
|
|
Packit |
67cb25 |
GSL_ERROR_VAL ("failed to allocate space for fsum_l", GSL_ENOMEM, 0);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
s->fsum_r = (double *) malloc (dim * sizeof (double));
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (s->fsum_r == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
free (s->fsum_l);
|
|
Packit |
67cb25 |
free (s->fmin_r);
|
|
Packit |
67cb25 |
free (s->fmin_l);
|
|
Packit |
67cb25 |
free (s->fmax_r);
|
|
Packit |
67cb25 |
free (s->fmax_l);
|
|
Packit |
67cb25 |
free (s->sigma_r);
|
|
Packit |
67cb25 |
free (s->sigma_l);
|
|
Packit |
67cb25 |
free (s->xmid);
|
|
Packit |
67cb25 |
free (s->x);
|
|
Packit |
67cb25 |
free (s);
|
|
Packit |
67cb25 |
GSL_ERROR_VAL ("failed to allocate space for fsum_r", GSL_ENOMEM, 0);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
s->fsum2_l = (double *) malloc (dim * sizeof (double));
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (s->fsum2_l == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
free (s->fsum_r);
|
|
Packit |
67cb25 |
free (s->fsum_l);
|
|
Packit |
67cb25 |
free (s->fmin_r);
|
|
Packit |
67cb25 |
free (s->fmin_l);
|
|
Packit |
67cb25 |
free (s->fmax_r);
|
|
Packit |
67cb25 |
free (s->fmax_l);
|
|
Packit |
67cb25 |
free (s->sigma_r);
|
|
Packit |
67cb25 |
free (s->sigma_l);
|
|
Packit |
67cb25 |
free (s->xmid);
|
|
Packit |
67cb25 |
free (s->x);
|
|
Packit |
67cb25 |
free (s);
|
|
Packit |
67cb25 |
GSL_ERROR_VAL ("failed to allocate space for fsum2_l", GSL_ENOMEM, 0);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
s->fsum2_r = (double *) malloc (dim * sizeof (double));
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (s->fsum2_r == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
free (s->fsum2_l);
|
|
Packit |
67cb25 |
free (s->fsum_r);
|
|
Packit |
67cb25 |
free (s->fsum_l);
|
|
Packit |
67cb25 |
free (s->fmin_r);
|
|
Packit |
67cb25 |
free (s->fmin_l);
|
|
Packit |
67cb25 |
free (s->fmax_r);
|
|
Packit |
67cb25 |
free (s->fmax_l);
|
|
Packit |
67cb25 |
free (s->sigma_r);
|
|
Packit |
67cb25 |
free (s->sigma_l);
|
|
Packit |
67cb25 |
free (s->xmid);
|
|
Packit |
67cb25 |
free (s->x);
|
|
Packit |
67cb25 |
free (s);
|
|
Packit |
67cb25 |
GSL_ERROR_VAL ("failed to allocate space for fsum2_r", GSL_ENOMEM, 0);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
s->hits_r = (size_t *) malloc (dim * sizeof (size_t));
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (s->hits_r == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
free (s->fsum2_r);
|
|
Packit |
67cb25 |
free (s->fsum2_l);
|
|
Packit |
67cb25 |
free (s->fsum_r);
|
|
Packit |
67cb25 |
free (s->fsum_l);
|
|
Packit |
67cb25 |
free (s->fmin_r);
|
|
Packit |
67cb25 |
free (s->fmin_l);
|
|
Packit |
67cb25 |
free (s->fmax_r);
|
|
Packit |
67cb25 |
free (s->fmax_l);
|
|
Packit |
67cb25 |
free (s->sigma_r);
|
|
Packit |
67cb25 |
free (s->sigma_l);
|
|
Packit |
67cb25 |
free (s->xmid);
|
|
Packit |
67cb25 |
free (s->x);
|
|
Packit |
67cb25 |
free (s);
|
|
Packit |
67cb25 |
GSL_ERROR_VAL ("failed to allocate space for fsum2_r", GSL_ENOMEM, 0);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
s->hits_l = (size_t *) malloc (dim * sizeof (size_t));
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (s->hits_l == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
free (s->hits_r);
|
|
Packit |
67cb25 |
free (s->fsum2_r);
|
|
Packit |
67cb25 |
free (s->fsum2_l);
|
|
Packit |
67cb25 |
free (s->fsum_r);
|
|
Packit |
67cb25 |
free (s->fsum_l);
|
|
Packit |
67cb25 |
free (s->fmin_r);
|
|
Packit |
67cb25 |
free (s->fmin_l);
|
|
Packit |
67cb25 |
free (s->fmax_r);
|
|
Packit |
67cb25 |
free (s->fmax_l);
|
|
Packit |
67cb25 |
free (s->sigma_r);
|
|
Packit |
67cb25 |
free (s->sigma_l);
|
|
Packit |
67cb25 |
free (s->xmid);
|
|
Packit |
67cb25 |
free (s->x);
|
|
Packit |
67cb25 |
free (s);
|
|
Packit |
67cb25 |
GSL_ERROR_VAL ("failed to allocate space for fsum2_r", GSL_ENOMEM, 0);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
s->dim = dim;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_monte_miser_init (s);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return s;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_monte_miser_init (gsl_monte_miser_state * s)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* We use 8 points in each halfspace to estimate the variance. There are
|
|
Packit |
67cb25 |
2*dim halfspaces. A variance estimate requires a minimum of 2 points. */
|
|
Packit |
67cb25 |
s->min_calls = 16 * s->dim;
|
|
Packit |
67cb25 |
s->min_calls_per_bisection = 32 * s->min_calls;
|
|
Packit |
67cb25 |
s->estimate_frac = 0.1;
|
|
Packit |
67cb25 |
s->alpha = 2.0;
|
|
Packit |
67cb25 |
s->dither = 0.0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
void
|
|
Packit |
67cb25 |
gsl_monte_miser_free (gsl_monte_miser_state * s)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
RETURN_IF_NULL (s);
|
|
Packit |
67cb25 |
free (s->hits_r);
|
|
Packit |
67cb25 |
free (s->hits_l);
|
|
Packit |
67cb25 |
free (s->fsum2_r);
|
|
Packit |
67cb25 |
free (s->fsum2_l);
|
|
Packit |
67cb25 |
free (s->fsum_r);
|
|
Packit |
67cb25 |
free (s->fsum_l);
|
|
Packit |
67cb25 |
free (s->fmin_r);
|
|
Packit |
67cb25 |
free (s->fmin_l);
|
|
Packit |
67cb25 |
free (s->fmax_r);
|
|
Packit |
67cb25 |
free (s->fmax_l);
|
|
Packit |
67cb25 |
free (s->sigma_r);
|
|
Packit |
67cb25 |
free (s->sigma_l);
|
|
Packit |
67cb25 |
free (s->xmid);
|
|
Packit |
67cb25 |
free (s->x);
|
|
Packit |
67cb25 |
free (s);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
void
|
|
Packit |
67cb25 |
gsl_monte_miser_params_get (const gsl_monte_miser_state * s, gsl_monte_miser_params * p)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
p->estimate_frac = s->estimate_frac;
|
|
Packit |
67cb25 |
p->min_calls = s->min_calls;
|
|
Packit |
67cb25 |
p->min_calls_per_bisection = s->min_calls_per_bisection;
|
|
Packit |
67cb25 |
p->alpha = s->alpha;
|
|
Packit |
67cb25 |
p->dither = s->dither;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
void
|
|
Packit |
67cb25 |
gsl_monte_miser_params_set (gsl_monte_miser_state * s, const gsl_monte_miser_params * p)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
s->estimate_frac = p->estimate_frac;
|
|
Packit |
67cb25 |
s->min_calls = p->min_calls;
|
|
Packit |
67cb25 |
s->min_calls_per_bisection = p->min_calls_per_bisection;
|
|
Packit |
67cb25 |
s->alpha = p->alpha;
|
|
Packit |
67cb25 |
s->dither = p->dither;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static int
|
|
Packit |
67cb25 |
estimate_corrmc (gsl_monte_function * f,
|
|
Packit |
67cb25 |
const double xl[], const double xu[],
|
|
Packit |
67cb25 |
size_t dim, size_t calls,
|
|
Packit |
67cb25 |
gsl_rng * r,
|
|
Packit |
67cb25 |
gsl_monte_miser_state * state,
|
|
Packit |
67cb25 |
double *result, double *abserr,
|
|
Packit |
67cb25 |
const double xmid[], double sigma_l[], double sigma_r[])
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
size_t i, n;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double *x = state->x;
|
|
Packit |
67cb25 |
double *fsum_l = state->fsum_l;
|
|
Packit |
67cb25 |
double *fsum_r = state->fsum_r;
|
|
Packit |
67cb25 |
double *fsum2_l = state->fsum2_l;
|
|
Packit |
67cb25 |
double *fsum2_r = state->fsum2_r;
|
|
Packit |
67cb25 |
size_t *hits_l = state->hits_l;
|
|
Packit |
67cb25 |
size_t *hits_r = state->hits_r;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double m = 0.0, q = 0.0;
|
|
Packit |
67cb25 |
double vol = 1.0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < dim; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
vol *= xu[i] - xl[i];
|
|
Packit |
67cb25 |
hits_l[i] = hits_r[i] = 0;
|
|
Packit |
67cb25 |
fsum_l[i] = fsum_r[i] = 0.0;
|
|
Packit |
67cb25 |
fsum2_l[i] = fsum2_r[i] = 0.0;
|
|
Packit |
67cb25 |
sigma_l[i] = sigma_r[i] = -1;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (n = 0; n < calls; n++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double fval;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
unsigned int j = (n/2) % dim;
|
|
Packit |
67cb25 |
unsigned int side = (n % 2);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < dim; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double z = gsl_rng_uniform_pos (r) ;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (i != j)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
x[i] = xl[i] + z * (xu[i] - xl[i]);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (side == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
x[i] = xmid[i] + z * (xu[i] - xmid[i]);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
x[i] = xl[i] + z * (xmid[i] - xl[i]);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
fval = GSL_MONTE_FN_EVAL (f, x);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* recurrence for mean and variance */
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double d = fval - m;
|
|
Packit |
67cb25 |
m += d / (n + 1.0);
|
|
Packit |
67cb25 |
q += d * d * (n / (n + 1.0));
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute the variances on each side of the bisection */
|
|
Packit |
67cb25 |
for (i = 0; i < dim; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (x[i] <= xmid[i])
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
fsum_l[i] += fval;
|
|
Packit |
67cb25 |
fsum2_l[i] += fval * fval;
|
|
Packit |
67cb25 |
hits_l[i]++;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
fsum_r[i] += fval;
|
|
Packit |
67cb25 |
fsum2_r[i] += fval * fval;
|
|
Packit |
67cb25 |
hits_r[i]++;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < dim; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double fraction_l = (xmid[i] - xl[i]) / (xu[i] - xl[i]);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (hits_l[i] > 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
fsum_l[i] /= hits_l[i];
|
|
Packit |
67cb25 |
sigma_l[i] = sqrt (fsum2_l[i] - fsum_l[i] * fsum_l[i] / hits_l[i]);
|
|
Packit |
67cb25 |
sigma_l[i] *= fraction_l * vol / hits_l[i];
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (hits_r[i] > 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
fsum_r[i] /= hits_r[i];
|
|
Packit |
67cb25 |
sigma_r[i] = sqrt (fsum2_r[i] - fsum_r[i] * fsum_r[i] / hits_r[i]);
|
|
Packit |
67cb25 |
sigma_r[i] *= (1 - fraction_l) * vol / hits_r[i];
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
*result = vol * m;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (calls < 2)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
*abserr = GSL_POSINF;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
*abserr = vol * sqrt (q / (calls * (calls - 1.0)));
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|