molpopgen / fwdpp

fwdpp is a C++ template library for implementing efficient forward-time population genetic simulations
http://fwdpp.readthedocs.io
GNU General Public License v3.0
27 stars 11 forks source link

Likely bug in "sugar" code. #17

Closed molpopgen closed 9 years ago

molpopgen commented 9 years ago

This is a bit informal, but here's the scoop:

While working on the next release, I noticed that Hudson and Kaplan's Rmin (http://www.ncbi.nlm.nih.gov/pmc/articles/PMC1462432/) was not coming out right for neutral simulations with recombination.

I checked previous versions of the library, and 0.2.9 "is fine" but 0.3.0 shows this problem. 0.3.0 is also when the example programs switched from using the original interface to using the new "sugar" interface, which streamlines setting up simulations, etc.

I have confirmed that compiling diploid_ind.cc from 0.2.9 (again, this doesn't use the "sugar" code to implement a simulation) against fwdpp's dev branch gives the correct distribution of Rmin. This observation puts the bug somewhere in the sugar code, so now I have to track that down.

molpopgen commented 9 years ago

The bug is in fwdpp/sugar/GSLrng_t.hpp, which defines a "smart pointer" wrapper to hold a gsl_rng *. It appears that my implementation is, in fact, not all that smart.

The workaround is to not use this type in simulations right now, and instead initialized the gsl_rng objects manually, which is quite easy. For example:

gsl_rng * r = gsl_rng_alloc(gsl_rng_mt19937); gsl_rng_set(r,seed);

molpopgen commented 9 years ago

This issue was due to the GSLrng_t object inadvertently copying its underlying pointer, and doing so incorrectly. The commits are here: https://github.com/molpopgen/fwdpp/commit/050e8b7c5e49dadb30d0b3f4b64ef80f8db0c394

The required API change is that the implicit conversion from GSLrng_t to gsl_rng * has been removed from the fwdpp sugar layer. Instead, the member function get() has been added, which wraps the underlying unique_ptr's get().

molpopgen commented 9 years ago

Closed with release of 0.3.2