role-model / roleR

R package implementing the RoLE model
https://role-model.github.io/roleR
GNU General Public License v3.0
1 stars 2 forks source link

Stochastic bug with msprime "Attempt to sample lineage in a population with size=0" #132

Closed ajrominger closed 1 year ago

ajrominger commented 1 year ago

running this code:

p <- roleParams(dispersal_prob = 0.5, alpha = 5000, individuals_local = 100,
                speciation_local = 0, mutation_rate = 1e-06, num_basepairs = 500,
                niter = 15000, niterTimestep = 5000)

m <- runRole(roleModel(p))

It works >50% of the time, but occasionally throws this error:

Error: msprime._msprime.InputError: Input error in initialise: Attempt to sample lineage in a population with size=0

Traceback shows

10.
stop(structure(list(message = "msprime._msprime.InputError: Input error in initialise: Attempt to sample lineage in a population with size=0\n",
call = NULL, cppstack = structure(list(file = "", line = -1L,
stack = c("1 reticulate.so 0x0000000106d917de _ZN4Rcpp9exceptionC2EPKcb + 222",
"2 reticulate.so 0x0000000106d99665 _ZN4Rcpp4stopERKNSt3__112basic_stringIcNS0_11char_traitsIcEENS0_9allocatorIcEEEE + 53", ...
9.
__init__ at ancestry.py#1318
8.
_parse_sim_ancestry at ancestry.py#1010
7.
sim_ancestry at ancestry.py#1207
6.
py_msprime_simulate at role_msprime.py#104
5.
py_msprime_simulate(J_m, J, curtime, metaTree, metaAbund, localAbund,
spAbundHarmMean, localTDiv, alpha, sequence_length, mu, verbose = FALSE) at sim_seqs.R#47
4.
sim_seqs(m) at runRole.R#66
3.
.local(x)
2.
runRole(roleModel(p)) at runRole.R#29
1.
runRole(roleModel(p))
ajrominger commented 1 year ago

and the frequency of the error increased with individuals_local such that when individuals_local = 1000 the error happens 100% of the time

ajrominger commented 1 year ago

caused by #127. solved thanks to @isaacovercast but when asked how he would like to be credited he says "don't worry about it" yet again closing issues without touching a keyboard

ajrominger commented 1 year ago

even after the fix in b669ff7 there is still a problem here. Seems like for some reason on rare occasion the harmonic mean will be 0 but the species abundance will be 1

ajrominger commented 1 year ago

by pure and unbelievable luck i found an off-by-one bug in the for loop that updates harmonic means. the bug made it so the last element of the harmonic means vector never gets updated.

it was

for(int s = 0; s < d.nTipsP(0) - 1; s++)

I made it

for(int s = 0; s < d.nTipsP(0); s++)

no more msprime errors, but we should look for more of these off-by-oners

ajrominger commented 1 year ago

fixed by ee4fdfc