version 1.1-5

ChangeLog

@@ -1,5 +1,16 @@
-			CHANGES IN OUTBREAKER VERSION 1.1-3
+			CHANGES IN OUTBREAKER VERSION 1.1-5
+BUG FIXES
+	o genetic likelihood with large numbers of mutation was prone to
+	numerical approximation errors; this is no longer the case.
+
+
+			CHANGES IN OUTBREAKER VERSION 1.1-4
+NEW FEATURES
+	o implemented a local likelihood computation providing substantial
+	improvements in speed and scalability
+
 
+			CHANGES IN OUTBREAKER VERSION 1.1-3
 BUG FIXES
 	o fixed a bug due to the odd conversion of long numbers into
 	characters by read.table - handling specific errors with R 3.1.0
@@ -8,7 +19,6 @@ BUG FIXES
 
 
 			CHANGES IN OUTBREAKER VERSION 1.1-2
-
 BUG FIXES
 	o fixed a bug due to the odd conversion of long numbers into
 	characters by read.table.

DESCRIPTION

@@ -1,17 +1,17 @@
 Package: outbreaker
-Version: 1.1-4
-Date: 2014/12/05
-Title: Bayesian reconstruction of disease outbreaks by combining
-        epidemiologic and genomic data
+Version: 1.1-5
+Date: 2015/02/19
+Title: Bayesian Reconstruction of Disease Outbreaks by Combining
+        Epidemiologic and Genomic Data
 Author: Thibaut Jombart <t.jombart@imperial.ac.uk>, Anne Cori, Joel Hellewell
 Maintainer: Thibaut Jombart <t.jombart@imperial.ac.uk>
-Suggests: EpiEstim
+Suggests: EpiEstim, OutbreakTools
 Depends: R (>= 3.0.0), parallel
 Imports: utils, ape, igraph, adegenet
 URL: http://sites.google.com/site/therepiproject/r-pac/outbreaker
 Description: Bayesian reconstruction of disease outbreaks using epidemiological and genetic information.
 License: GPL (>= 2)
-Packaged: 2014-12-05 16:56:12 UTC; thibaut
 NeedsCompilation: yes
+Packaged: 2015-02-20 17:01:00 UTC; thibaut
 Repository: CRAN
-Date/Publication: 2014-12-05 18:08:09
+Date/Publication: 2015-02-20 23:49:08

MD5

@@ -1,6 +1,6 @@
-6f7601de240d16818c58279443234e82 *ChangeLog
-db70d9d3bb0eed05948e2038464b7a52 *DESCRIPTION
-b1c3a9ccc2b95d7b57297d94b2386013 *NAMESPACE
+601162e3337957db19a32260732ea767 *ChangeLog
+ae163759b1dc2fa2f28cc07a70675356 *DESCRIPTION
+8d070a3fc6a69b0038d89e7cb322283f *NAMESPACE
 cbdf308393d62f46f0d8bc2f5d4470b6 *R/interface.R
 c0d0d9160caaf41e22f0265b788c7e00 *R/mutationRates.R
 9a32af7fc041e4b0647256cd3b1c0d78 *R/plot.R
@@ -11,7 +11,7 @@ dc26678070da5a3942b779aa0670d4f3 *R/selectChains.R
 92d8aa7bf09017257369dcfdf00a82cb *R/zzz.R
 72b5047fee7ef404647d8c939b9ec91b *configure
 a6d485ade20b5042ab8459e99f45c895 *configure.ac
-999579b31cc165da3ac161c33966957e *data/fakeOutbreak.RData
+be94f9222da3b7b700fd9465a911de54 *data/fakeOutbreak.RData
 2cb344c2133da976ccf7ab341e81353b *inst/CITATION
 0d5cd8ac8f915095773450ce55e2f669 *man/Rstat.Rd
 c774351eefcfdd1193bc820ceb2efd8c *man/fakeOutbreak.Rd
@@ -30,15 +30,15 @@ d88c062ac774d6aa5583ea359f503b61 *man/tTree.Rd
 06a461664c434629cb567fc5d25d8156 *src/genclasses.h
 3301c6dc9d874a2f6dee7e6170db2985 *src/init.c
 4975d1b53399a67f5f94358a013d26b7 *src/init.h
-e646b5a26c4bd117ef9d914958f6e8b4 *src/likelihood.c
-4a0b6f2e48d93243737e624fd822b598 *src/likelihood.h
+83a270f777d54a124b3681a42bf7103a *src/likelihood.c
+38fcf57840a35802c45118d43ced95bd *src/likelihood.h
 02ba8d225eb3b405cd237e879b0abbba *src/matvec.c
 b5e29a394d577993c7b301e2235c40ec *src/matvec.h
 8325cc315d2881988b72a666ee60d945 *src/mcmc.c
 1aef3fd5435b818318bb0fc3af1d1295 *src/mcmc.h
 7ea50bca494303845e8ba4e62cd465ff *src/moves.c
 d5e38df1932ec27dfa188ca56583868b *src/moves.h
-6849d5b7b70fbb2dcc404db37bbe3d6b *src/outbreaker.c
+17cf15ebb49399c476af5e40b998862a *src/outbreaker.c
 ab6b405475a66ad44d59aa153d341ed3 *src/prior.c
 8a56d478207d395beaf8a9c34758fa63 *src/prior.h
 1aa3691b3ecdeb22e3661e1872d6f918 *src/structures.c

NAMESPACE

@@ -1,8 +1,4 @@
 
-## Export all names
-exportPattern("^[^\\.]")
-
-
 ## Import all packages listed as Imports or Depends
 importFrom(utils, "packageDescription")
 
@@ -18,5 +14,16 @@ importFrom(adegenet, "bluepal", "redpal", "greenpal", "seasun", "spectral", "fun
 
 importFrom(parallel, "detectCores", "makeCluster", "clusterEvalQ", "clusterExport", "parLapply", "stopCluster")
 
+
+## Export all names
+exportPattern("^[^\\.]")
+
+## register S3 methods
+S3method(labels, simOutbreak)
+S3method(plot, simOutbreak)
+S3method(print, simOutbreak)
+S3method(plot, tTree)
+
+
 ## Load DLL
 useDynLib(outbreaker)

src/likelihood.c

@@ -126,10 +126,23 @@ double gsl_ran_poisson_pdf_fixed(unsigned int k, double mu){
 /* Formula: p = mu^nbmut x (1-mu)^{(kappa x nbnucl) - nbmut} */
 double proba_mut(int nbmut, int nbnucl, int kappa, double mu){
     double out=0.0;
+    /* old version, numerical approx problems for nbmut high */
     out = gsl_sf_pow_int(mu, nbmut) * gsl_sf_pow_int(1.0-mu, kappa*nbnucl - nbmut);
     return out;
 }
 
+/* Compute log-proba of ... mutations */
+/* - binomial density without the binomial coefficient - */
+/* Assume that no site can mutate twice over kappa generations */
+/* there are kappa*nbnucl 'draws' (possible mutation occurences) */
+/* Formula: p = mu^nbmut x (1-mu)^{(kappa x nbnucl) - nbmut} */
+double log_proba_mut(int nbmut, int nbnucl, int kappa, double mu){
+    double out=0.0;
+    /* note: integers are automatically promoted to double here */
+    out = nbmut * log(mu) + (kappa*nbnucl - nbmut) * log(1.0-mu);
+    return out;
+}
+
 
 
 
@@ -240,7 +253,7 @@ double loglikelihood_gen_i(int i, data *dat, dna_dist *dnaInfo, param *par, gsl_
 	if(mutation1_ij(i, ances, dat, dnaInfo)<0 || par->kappa_temp<0 || par->mu1<0 || par->mu1 >1){ /* fool proof */
 	  out += NEARMINUSINF;
 	} else {
-	  out += log(proba_mut(mutation1_ij(i, ances, dat, dnaInfo), com_nucl_ij(i, ances, dat, dnaInfo), par->kappa_temp, par->mu1));
+	  out += log_proba_mut(mutation1_ij(i, ances, dat, dnaInfo), com_nucl_ij(i, ances, dat, dnaInfo), par->kappa_temp, par->mu1);
 	}
       }
       break;
@@ -252,10 +265,10 @@ double loglikelihood_gen_i(int i, data *dat, dna_dist *dnaInfo, param *par, gsl_
 	  out += NEARMINUSINF;
 	} else {
 	  /* transitions */
-	  out += log(proba_mut(mutation1_ij(i, ances, dat, dnaInfo), com_nucl_ij(i, ances, dat, dnaInfo), par->kappa_temp, par->mu1));
+	  out += log_proba_mut(mutation1_ij(i, ances, dat, dnaInfo), com_nucl_ij(i, ances, dat, dnaInfo), par->kappa_temp, par->mu1);
 
 	  /* transversions */
-	  out += log(proba_mut(mutation2_ij(i, ances, dat, dnaInfo), com_nucl_ij(i, ances, dat, dnaInfo), par->kappa_temp,  par->gamma * par->mu1));
+	  out += log_proba_mut(mutation2_ij(i, ances, dat, dnaInfo), com_nucl_ij(i, ances, dat, dnaInfo), par->kappa_temp,  par->gamma * par->mu1);
 	}
       }
       break;

src/likelihood.h

@@ -20,6 +20,9 @@ double gsl_ran_poisson_pdf_fixed(unsigned int k, double mu);
 
 double proba_mut(int nbmut, int nbnucl, int kappa, double mu);
 
+double log_proba_mut(int nbmut, int nbnucl, int kappa, double mu);
+
+
 /*
   ====================
   LIKELIHOOD FUNCTIONS

src/outbreaker.c

@@ -91,18 +91,18 @@ void R_outbreaker(unsigned char *DNAbinInput, int *Tcollec, int *n, int *nSeq, i
 
     /* COMPUTE PRIORS */
     double logPrior = logprior_all(par);
-    Rprintf("\nPrior value (log): %.10f\n", logPrior);/* fflush(stdout); */
+    /* Rprintf("\nPrior value (log): %.10f\n", logPrior);/\* fflush(stdout); *\/ */
 
    /* COMPUTE LIKELIHOOD */
     double logLike = loglikelihood_all(dat, dnaInfo, spatialInfo, gen, par, rng);
-    Rprintf("\n\n = Initial Log-likelihood value: %f\n", logLike);
+    /* Rprintf("\n\n = Initial Log-likelihood value: %f\n", logLike); */
 
     /* COMPUTE POSTERIOR */
     double logPost = logposterior_all(dat, dnaInfo, spatialInfo, gen, par, rng);
-    Rprintf("\nLog-posterior value: %.10f\n", logPost);/* fflush(stdout); */
+    /* Rprintf("\nLog-posterior value: %.10f\n", logPost);/\* fflush(stdout); *\/ */
 
     /* ALLOCATE AND INITIALIZE MCMC PARAMETERS */
-    Rprintf("\nBefore check init LL\n");/* fflush(stdout); */
+    /* Rprintf("\nBefore check init LL\n");/\* fflush(stdout); *\/ */
 
     mcmcPar = alloc_mcmc_param(N);
     init_mcmc_param(mcmcPar, par, dat, (bool) *moveMut, moveAlpha, moveKappa, (bool) *moveTinf, 
@@ -116,8 +116,8 @@ void R_outbreaker(unsigned char *DNAbinInput, int *Tcollec, int *n, int *nSeq, i
       warning("\n\n!WARNING! Initial state of the chain has a likelihood of zero. The chain may never converge. Please consider using a different initial tree.\n");
     }
 
-    Rprintf("\nAfter check init LL\n");/* fflush(stdout); */
-    Rprintf("\nBefore MCMC\n");/* fflush(stdout); */
+    /* Rprintf("\nAfter check init LL\n");/\* fflush(stdout); *\/ */
+    /* Rprintf("\nBefore MCMC\n");/\* fflush(stdout); *\/ */
 
     /* RUN MCMC */
     mcmc(*nIter, *outputEvery, *resFileName, *tuneFileName, *tuneEvery,