thibautjombart / outbreaker

Disease outbreak reconstruction from epidemiological and genetic data
2 stars 1 forks source link

Error...[in: structures.c->gentime_dens] Trying to get density for 150 time units (max: 150). Exiting. #1

Closed pcdjohnson closed 8 years ago

pcdjohnson commented 9 years ago

Hi,

While running outbreaker on a large data set of 2026 cases with complete sequence data, I got this type of error: Error in outbreaker( ... [in: structures.c->gentime_dens] Trying to get density for 150 time units (max: 150). Exiting. I reproduced this error in a simple analysis of 2000 simulated dates with no sequence data. The console output from the full session is pasted below, followed by the sessionInfo() output.

Thanks for any help, Paul

PS this issue was originally posted here: https://groups.google.com/forum/?hl=en#!topic/r-epi/l4lOnp-eUKk

R version 3.2.1 (2015-06-18) -- "World-Famous Astronaut"
[snip]
> library(outbreaker)
Loading required package: parallel

   === outbreaker 1.1-5 is loaded ===

Information & Documentation -> check the R-epi project: 
http://sites.google.com/site/therepiproject

Questions -> ask the R-epi forum: 
r-epi@googlegroups.com
> set.seed(4321)
> dates <- cumsum(rpois(2000, 0.1))
> res <-  
+   outbreaker(
+     dates = dates,
+     w.dens = c(0, rep(1/14, 14)), 
+     f.dens = c(0, rep(1/7, 7)),
+     init.tree = "random",
+     n.iter = 1e5, burnin = 5e4, seed = 1234)
Error in outbreaker(dates = dates, w.dens = c(0, rep(1/14, 14)), f.dens = c(0,  : 

[in: structures.c->gentime_dens]
Trying to get density for 150 time units (max: 150). Exiting.

> sessionInfo()
R version 3.2.1 (2015-06-18)
Platform: x86_64-apple-darwin13.4.0 (64-bit)
Running under: OS X 10.10.5 (Yosemite)

locale:
[1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8

attached base packages:
[1] parallel  stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
[1] outbreaker_1.1-5

loaded via a namespace (and not attached):
 [1] Rcpp_0.12.0      spdep_0.5-88     plyr_1.8.3       LearnBayes_2.15 
 [5] tools_3.2.1      boot_1.3-17      digest_0.6.8     gtable_0.1.2    
 [9] nlme_3.1-122     lattice_0.20-33  Matrix_1.2-2     igraph_1.0.1    
[13] shiny_0.12.2     DBI_0.3.1        proto_0.3-10     coda_0.17-1     
[17] dplyr_0.4.2      stringr_1.0.0    ade4_1.7-2       grid_3.2.1      
[21] R6_2.1.1         sp_1.1-1         ggplot2_1.0.1    reshape2_1.4.1  
[25] seqinr_3.1-3     deldir_0.1-9     magrittr_1.5     scales_0.3.0    
[29] htmltools_0.2.6  MASS_7.3-43      splines_3.2.1    adegenet_2.0.0  
[33] assertthat_0.1   mime_0.3         ape_3.3          colorspace_1.2-6
[37] xtable_1.7-4     httpuv_1.3.3     stringi_0.5-5    munsell_0.4.2   
> 
thibautjombart commented 9 years ago

Hi there, sorry for the late reply. This was due to long time intervals generated by the initial random tree. Fixed as of 83b4d2f7d4abe420eef836d9a6ce43a89c1ee030 Behaviour of the function is now a bit more clever. The generation time distribution is automatically completed if needed by an exponential tail summing to a low value (1e-4) to cover the temporal range, so that even initial star trees should not return a -Inf temporal log likelihood. Your example now runs for the first few hundred iterations. I'll wait for you to confirm to close the issue.

pcdjohnson commented 9 years ago

Great, thanks very much -- that seems to have fixed the problem. I haven't had the patience to get the simple simulated example above to complete over 1e5 iterations, but a short run using find.import.n = 5, sample.every = 50, n.iter = 1e4, burnin = 5e3 did work. I also ran it on the the real data (~2000 real cases with complete simulated sequence data + spatial data) that generated the error in the first place. It worked very nicely.