[Bug-gsl] Bug in gsl-1.10/monte/miser.c

[Brad Bell]

I wish to report what I think is a bug in
   Bug in gsl-1.10/monte/miser.c
An extensive discussion of this bug (and how it was found) can be found at
   http://bugzilla.rfpk.washington.edu/show_bug.cgi?id=585

I have further isolated the problem and am attaching a simple bash 
script (bug.sh) that both demonstrates the problem and its fix. When I 
run this script on my Cygwin system I get the attached output 
(bug.out).  You will notice that the patched version of  miser.c 
estimates the integral while the original one does not.

Running the command
   ./bug.sh
will create a patched version of miser.c called
   ./gsl_miser_bug/miser.c
A description of the problem can be found in the patched version of 
miser.c (search for the text "Begin Patch" in miser.c).

Note. You will need to edit the bug.sh file to specify the location of 
the gsl source code on your machine.
#! /bin/bash
#
# You must replace this with the proper path for your system
GSL_SOURCE=$HOME/trash/gsl-1.10
#
if [ ! -d $GSL_SOURCE ]
then
        echo "You must edit this file and specify the gsl source location."
        echo "You can get a of the gsl source with the following commands:"
        echo "  wget ftp://sources.redhat.com/pub/gsl/gsl-1.10.tar.gz";
        echo "  tar -xzf gsl-1.10.tar.gz"
        exit 1
fi
#
if [ -e gsl_miser_bug ]
then
        rm gsl_miser_bug/bug.*
        rm gsl_miser_bug/miser.*
        if ! rmdir gsl_miser_bug
        then
                echo "cannot remove old version of gsl_miser_bug"
        fi
fi
mkdir gsl_miser_bug
if ! cd gsl_miser_bug
then
        echo "cannot create new gsl_miser_bug directory"
        exit 1
fi
#
#
echo "Create a local patched version of miser.c"
cat << EOF > miser.sed
s%#include <config.h>%// Begin Patch: Brad Bell 22-03-2008 \\
// Comment out this include so that miser.c will compile separately\\
// &\\
// End Patch%
#
# Note the best fix, but this way it is done in one line
s%^\( *\)if \((sigma_l.i. >= 0 && sigma_r.i. >= 0)\)%// Begin Patch: Brad Bell 
22-03-2008 \\
// Protect against the case where both sigma_l and sigma_r are zero.\\
// This case leads to weight_l = 0, weight_r = 0, and\\
// values calls_l and calls_r are conversion of nan to int.\\
\\1if ( \\2 \\
          \&\& (sigma_l[i] != 0 || sigma_r[i] != 0) \\
\\1)\\
// End Patch%
EOF
sed -f miser.sed < $GSL_SOURCE/monte/miser.c > miser.c
#
# Test program (modified version of $GSL_SOURCE/doc/examples/monte.c)
cat << EOF > bug.c
#include <stdlib.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_monte.h>
#include <gsl/gsl_monte_miser.h>

double exact = .01;

double
g (double *k, size_t dim, void *params)
{
  if( k[0] < 0 || k[0] > 1. )
    return 1.;
  else
    return 0.;
}

int
main (void)
{
  double res, err;

  double xl[3] = { 0, 0, 0 };
  double xu[3] = { 1, 1, 1 };

  xl[0] -= exact / 2.;
  xu[0] += exact / 2.;


  const gsl_rng_type *T;
  gsl_rng *r;

  gsl_monte_function G = { &g, 3, 0 };

  size_t calls = 500000;

  gsl_rng_env_setup ();

  T = gsl_rng_default;
  r = gsl_rng_alloc (T);

  {
    gsl_monte_miser_state *s = gsl_monte_miser_alloc (3);
    gsl_monte_miser_integrate (&G, xl, xu, 3, calls, r, s,
                               &res, &err);
    gsl_monte_miser_free (s);

    printf ("result = % .6f\n", res);
    printf ("sigma  = % .6f\n", err);
    printf ("exact  = % .6f\n", exact);
  }

  gsl_rng_free (r);

  return 0;
}
EOF
#
echo
echo "Compile, and execute the test program with patched version of miser."
echo "gcc bug.c miser.c -lgsl -o bug"
gcc bug.c miser.c -lgsl -o bug
echo "./bug"
./bug
#
echo
echo "Compile, and execute the test program with library version of miser."
echo "gcc bug.c -lgsl -o bug"
gcc bug.c -lgsl -o bug
echo "./bug"
./bug

Create a local patched version of miser.c

Compile, and execute the test program with patched version of miser.
gcc bug.c miser.c -lgsl -o bug
Info: resolving _gsl_rng_default by linking to __imp__gsl_rng_default 
(auto-import)
./bug
result =  0.010132
sigma  =  0.000309
exact  =  0.010000

Compile, and execute the test program with library version of miser.
gcc bug.c -lgsl -o bug
Info: resolving _gsl_rng_default by linking to __imp__gsl_rng_default 
(auto-import)
./bug
gsl: /usr/src/gsl-1.10-1/src/gsl-1.10/monte/miser.c:115: ERROR: insufficient 
calls for subvolume
Default GSL error handler invoked.
./bug.sh: line 122:  2776 Aborted                 (core dumped) ./bug