Sure, attached is the simplest program I can get to demonstrate the nan
result (on my computer - perhaps this is hardware-dependent). It took me a
little while to come up with it, as the bug seems somewhat difficult to
trigger. Hopefully this helps.
:) David
On Sat, Mar 21, 2015 at 11:17 PM, Patrick Alken <[email protected]>
wrote:
> Thanks for the report David. Do you have a simple program which
> demonstrates the bug?
>
> Patrick
>
> On 03/21/2015 01:59 AM, David Zaslavsky wrote:
> > I'm fairly sure I've identified a small bug in the VEGAS implementation:
> in
> > lines 280-283 of vegas.c on the master branch,
> >
> > if (var > 0) { wgt = 1.0 / var; }
> >
> > very small but nonzero values of var (like 1e-322) will cause wgt to be
> set
> > to inf, and that in turn causes gsl_monte_vegas_integrate to return nan.
> > The logical fix seems to be to treat any value of var satisfying
> > 1.0/var==inf as zero, i.e.
> >
> > - if (var > 0)
> > + if (var > 0 && gsl_finite(1.0 / var))
> >
> > I don't know if that would break anything else, but I can't imagine how.
> > Thoughts?
> >
> > I can send a proper patch file if desired (it's on my other computer).
> >
> > :) David
> >
>
>
>
#include <assert.h>
#include <stdio.h>
#include <stdlib.h>
#include "gsl/gsl_math.h"
#include "gsl/gsl_monte.h"
#include "gsl/gsl_monte_vegas.h"
#define check_finite(v) assert(gsl_finite(v) && !gsl_isnan(v))
static const double inf = 100;
double hypot_sq(double x, double y) {
return gsl_pow_2(x) + gsl_pow_2(y);
}
double integrand(double* coords, size_t ndims, void* params) {
double k = 0.8;
double k2 = k*k;
double lx = coords[0];
double ly = coords[1];
double lpx = coords[2];
double lpy = coords[3];
double lppx = coords[4];
double lppy = coords[5];
double F1 = exp(-hypot_sq(lppx, lppy));
double F2 = exp(-hypot_sq(k - lpx - lppx, -lpy - lppy));
double F3 = exp(-hypot_sq(k - lx - lppx, -ly - lppy));
double result;
assert(ndims == 6);
result = F1 * F2 * F3;
check_finite(result);
return result;
}
int main(const int argc, const char** argv) {
double res, err;
double min[] = {-inf, -inf, -inf, -inf, -inf, -inf};
double max[] = { inf, inf, inf, inf, inf, inf};
const size_t ndim = 6;
const gsl_rng_type* T;
gsl_rng *r;
T = gsl_rng_default;
r = gsl_rng_alloc(T);
gsl_rng_set(r, 100);
gsl_monte_function function_package;
function_package.f = integrand;
function_package.dim = ndim;
function_package.params = NULL;
gsl_monte_vegas_state *s = gsl_monte_vegas_alloc(ndim);
gsl_monte_vegas_integrate (&function_package, (double*)min, (double*)max, (size_t)ndim, 50000, r, s, &res, &err);
fprintf(stderr, "Expected output to demonstrate bug:\n result = -nan; errbound = 0.000000; chisq/dof = nan\n");
fprintf(stderr, "Actual output:\n result = % .6f; errbound = % .6f; chisq/dof = %.2f\n", res, err, gsl_monte_vegas_chisq(s));
gsl_monte_vegas_free(s);
gsl_rng_free(r);
return 0;
}
