bodkan / slendr

Population genetic simulations in R 🌍
https://bodkan.net/slendr
Other
54 stars 5 forks source link

Full implementation of the msprime back end for slendr models #76

Closed bodkan closed 2 years ago

bodkan commented 2 years ago

This is a massive update which not only adds an alternative backend for executing slendr models without any change, but it also implements a billion unit tests verifying that for a given model, the standard SLiM backend gives the same results as what the same model produces when run as an msprime simulation.

The unit tests currently test for matching allele frequency spectra for various simple demographic models (constant, contracting, expanding, exponentially contracting and exponentially expanding populations) and some f4 statistics (various models with or without gene flow between lineages) as computed from SLiM and msprime tree sequence files.

Additionally, all of these models are checked in pairs -- slendr allows flexible specification of model time units either in forward time (simulations running from generation N -> generation M where N < M) or backward times (M < N), so each model was expressed in forward time units and backward time units, and the tests make sure that a given model gives the same result (AFS vectors and f-statistics) regardless of in which orientation it is run.

This was a ton of work but I think it was worth it:

The last point was my initial motivation, the first two points are just a bonus, really.

bhaller commented 2 years ago

Awesome! This also serves as a new way of testing SLiM for correctness, too. I appreciate that additional validation; let me know if you find any bugs! :->

bodkan commented 2 years ago

I would be extremely surprised if I found something off with SLiM with these simple tests. After all, that those kinds of demographic models must have been run trillion times across hundreds of publications at this point, making SLiM extremely well battle-tested. :)

That said, these tests will certainly help me sleep a little better! The forward/backward voodoo I'm doing in slendr is a little scary and will need some serious refactoring... some day! For now, I'm glad the basic models work as expected.

(Also, thanks for keeping an eye on my progress here! ❤️)

bhaller commented 2 years ago

I would be extremely surprised if I found something off with SLiM with these simple tests. After all, that those kinds of demographic models must have been run trillion times across hundreds of publications at this point, making SLiM extremely well battle-tested. :)

I do not have your level of confidence. :-O There are a great many different code paths in SLiM; lots of dusty corners where bugs could lurk. Some code paths, and combinations of code paths, do not get hit by SLiM's self-tests and have probably never been hit by any user code either, but might get hit by a test suite like yours. And even the more well-used code paths could always acquire new bugs as a result of changes to the code; I try to touch the SLiM core as little as possible nowadays, for fear of introducing new bugs, but these things can happen. I live in fear of discovering (or introducing!) a fundamental bug in SLiM that invalidates results for published papers. Software validation is an incredibly difficult problem, especially software like SLiM that is (a) stochastic in some sense, and (b) optimized for speed. So I'm always glad to see more angles on its validation, such as what you have done here; and I always urge users to compare to results from other models/software, to think very carefully about what they expect their model to do and make sure it is working as they expect, etc. Bugs can and do happen.

That said, these tests will certainly help me sleep a little better! The forward/backward voodoo I'm doing in slendr is a little scary and will need some serious refactoring... some day! For now, I'm glad the basic models work as expected.

(Also, thanks for keeping an eye on my progress here! ❤️)

:->

bodkan commented 2 years ago

Hmm, yeah, that makes perfect sense. You raise a very good point with those corner cases and the stochasticity.

One idea I had in context of slendr was that it would be fun to simulate tons of random trees of different sizes and topologies (with something like rtree in ape), turn them into slendr models by sampling random N trajectories for different branches (i.e. populations) of the tree (constant, step size changes, exponential growing/shrinking), adding random "gene flow edges" between them... and running those models with SLiM and msprime backends, comparing them to one another using various tree sequence statistics. I'm not sure if this would have any practical scientific use beyond really comprehensive unit testing, but even that on itself would be more than enough.

Actually, this sounds like a nice, self-contained project for an intern interested in learning a bit about simulations in general. I'm making a note of this discussion for future reference when @FerRacimo and I chat about potential project ideas for new students. :)

bhaller commented 2 years ago

Sounds awesome!