From df2782ba954212f19c649495d953b3039d70e609 Mon Sep 17 00:00:00 2001 From: flouris Date: Tue, 27 Sep 2016 14:47:54 +0200 Subject: [PATCH] fixed floating point exception when constructing random trees, fixed malloc issue which caused unititalized data when cloning unrooted trees in MCMC, fixed sample size in AICc --- ChangeLog.md | 11 +++++++++++ configure.ac | 2 +- src/aic.c | 12 ++++++------ src/dp.c | 4 ++-- src/likelihood.c | 1 + src/mptp.h | 1 + src/random.c | 7 ++++++- src/rtree.c | 2 +- src/util.c | 22 ++++++++++++++++++++++ src/utree.c | 14 +++++++------- 10 files changed, 58 insertions(+), 18 deletions(-) create mode 100644 ChangeLog.md diff --git a/ChangeLog.md b/ChangeLog.md new file mode 100644 index 0000000..7e579d0 --- /dev/null +++ b/ChangeLog.md @@ -0,0 +1,11 @@ +# Change Log +All notable changes to `mptp` will be documented in this file. +This project adheres to [Semantic Versioning](http://semver.org/). + +## [0.2.0] - 2016-09-27 +### Fixed + - Floating point exception error when constructing random trees caused from + division by zero + - Allocation with malloc caused uninitialized variables when converting unrooted + tree to rooted for the MCMC method + - Sample size for the the AIC with a correction for finite sample sizes diff --git a/configure.ac b/configure.ac index 64c9213..8c70600 100644 --- a/configure.ac +++ b/configure.ac @@ -2,7 +2,7 @@ # Process this file with autoconf to produce a configure script. AC_PREREQ([2.63]) -AC_INIT([mptp], [0.1.1], [Tomas.Flouri@h-its.org]) +AC_INIT([mptp], [0.2.0], [Tomas.Flouri@h-its.org]) AM_INIT_AUTOMAKE([subdir-objects]) AC_LANG([C]) AC_CONFIG_SRCDIR([src/mptp.c]) diff --git a/src/aic.c b/src/aic.c index a17f906..836bfbb 100644 --- a/src/aic.c +++ b/src/aic.c @@ -705,7 +705,7 @@ void aic_mcmc(rtree_t * tree, double max_logl_aic = (method == PTP_METHOD_MULTI) ? vec[best_index].score_multi : vec[best_index].score_single; - double max_aic = aic(max_logl_aic, species_count, tree->leaves); + double max_aic = aic(max_logl_aic, species_count, tree->leaves+2); long coal_edge_count = 0; @@ -814,7 +814,7 @@ void aic_mcmc(rtree_t * tree, if (opt_mcmc_burnin == 1) { //densities[species_count].logl += logl; - densities[species_count].logl += -aic(logl, species_count, tree->leaves); + densities[species_count].logl += -aic(logl, species_count, tree->leaves+2); } if (opt_mcmc_sample == 1) @@ -909,8 +909,8 @@ void aic_mcmc(rtree_t * tree, *mcmc_min_logl = new_logl; - double aic_new_logl = -aic(new_logl, species_count+1, tree->leaves); - double aic_logl = -aic(logl, species_count, tree->leaves); + double aic_new_logl = -aic(new_logl, species_count+1, tree->leaves+2); + double aic_logl = -aic(logl, species_count, tree->leaves+2); /* Hastings ratio */ double a = exp(aic_new_logl - aic_logl) * (old_crnodes_count / new_snodes_count); @@ -1049,8 +1049,8 @@ void aic_mcmc(rtree_t * tree, else if (new_logl < *mcmc_min_logl) *mcmc_min_logl = new_logl; - double aic_new_logl = -aic(new_logl, species_count-1, tree->leaves); - double aic_logl = -aic(logl, species_count, tree->leaves); + double aic_new_logl = -aic(new_logl, species_count-1, tree->leaves+2); + double aic_logl = -aic(logl, species_count, tree->leaves+2); /* Hastings ratio */ double a = exp(aic_new_logl - aic_logl) * (old_snodes_count / new_crnodes_count); diff --git a/src/dp.c b/src/dp.c index c62aa94..74aae3c 100644 --- a/src/dp.c +++ b/src/dp.c @@ -187,12 +187,12 @@ void dp_ptp(rtree_t * tree, long method) if (method == PTP_METHOD_MULTI) { max = vec[0].score_multi; - double min_aic_score = aic(vec[0].score_multi, vec[0].species_count, tree->leaves); + double min_aic_score = aic(vec[0].score_multi, vec[0].species_count, tree->leaves+2); for (i = 1; i < tree->edge_count; i++) { if (vec[i].filled) { - double aic_score = aic(vec[i].score_multi, vec[i].species_count, tree->leaves); + double aic_score = aic(vec[i].score_multi, vec[i].species_count, tree->leaves+2); if (aic_score < min_aic_score) { min_aic_score = aic_score; diff --git a/src/likelihood.c b/src/likelihood.c index 0d15d5f..eb9a07b 100644 --- a/src/likelihood.c +++ b/src/likelihood.c @@ -50,5 +50,6 @@ int lrt(double nullmodel_logl, double ptp_logl, unsigned int df, double * pvalue double aic(double logl, long k, long n) { if (k > 1) k++; + return -2*logl + 2*k + (double)(2*k*(k + 1)) / (double)(n-k-1); } diff --git a/src/mptp.h b/src/mptp.h index dfa9768..efd3e3c 100644 --- a/src/mptp.h +++ b/src/mptp.h @@ -251,6 +251,7 @@ void progress_init(const char * prompt, unsigned long size); void progress_update(unsigned int progress); void progress_done(void); void * xmalloc(size_t size); +void * xcalloc(size_t nmemb, size_t size); void * xrealloc(void *ptr, size_t size); char * xstrchrnul(char *s, int c); char * xstrdup(const char * s); diff --git a/src/random.c b/src/random.c index 498b56c..96a75cd 100644 --- a/src/random.c +++ b/src/random.c @@ -86,9 +86,14 @@ double random_delimitation(rtree_t * root, max_species = root->max_species_count; g_rstate = rstate; + printf("Number of tips: %d\n", root->leaves); + assert(0); rand_long = nrand48(rstate); - species_count = (rand_long % root->max_species_count) + 1; + if (!root->max_species_count) + species_count = (rand_long % root->leaves) + 1; + else + species_count = (rand_long % root->max_species_count) + 1; rtree_t ** inner_node_list = (rtree_t **)xmalloc((size_t)species_count * sizeof(rtree_t *)); diff --git a/src/rtree.c b/src/rtree.c index c1aaeed..195a2c9 100644 --- a/src/rtree.c +++ b/src/rtree.c @@ -412,7 +412,7 @@ rtree_t * rtree_clone(rtree_t * node, rtree_t * parent) if (!node) return NULL; /* clone node */ - rtree_t * clone = (rtree_t *)xmalloc(sizeof(rtree_t)); + rtree_t * clone = (rtree_t *)xcalloc(1,sizeof(rtree_t)); memcpy(clone,node,sizeof(rtree_t)); clone->parent = parent; clone->data = NULL; diff --git a/src/util.c b/src/util.c index 4c94393..b195b51 100644 --- a/src/util.c +++ b/src/util.c @@ -69,6 +69,7 @@ void progress_done() fprintf(stderr, " \r%s %.0f%%\n", progress_prompt, 100.0); } +#if 0 void * xmalloc(size_t size) { const size_t alignment = 16; @@ -81,6 +82,27 @@ void * xmalloc(size_t size) return t; } +#else +void * xmalloc(size_t size) +{ + void * t; + t = malloc(size); + if (!t) + fatal("Unable to allocate enough memory."); + + return t; +} +#endif + +void * xcalloc(size_t nmemb, size_t size) +{ + void * t; + t = calloc(nmemb,size); + if (!t) + fatal("Unable to allocate enough memory."); + + return t; +} void * xrealloc(void *ptr, size_t size) { diff --git a/src/utree.c b/src/utree.c index 6572277..6fb6191 100644 --- a/src/utree.c +++ b/src/utree.c @@ -354,7 +354,7 @@ int utree_query_innernodes(utree_t * root, static rtree_t * utree_rtree(utree_t * unode) { - rtree_t * rnode = (rtree_t *)xmalloc(sizeof(rtree_t)); + rtree_t * rnode = (rtree_t *)xcalloc(1,sizeof(rtree_t)); rnode->event = EVENT_COALESCENT; @@ -390,7 +390,7 @@ utree_t * utree_longest_branchtip(utree_t * node, unsigned int tip_count) utree_t * outgroup = NULL; /* query tip nodes */ - utree_t ** tip_nodes_list = (utree_t **)xmalloc((size_t)tip_count * sizeof(utree_t *)); + utree_t ** tip_nodes_list = (utree_t **)xcalloc(1,(size_t)tip_count * sizeof(utree_t *)); utree_query_tipnodes(node, tip_nodes_list); for (i = 0; i < tip_count; ++i) @@ -413,7 +413,7 @@ rtree_t * utree_crop(utree_t * lca) if (!lca->back->next) return NULL; - rtree_t * root = (rtree_t *)xmalloc(sizeof(rtree_t)); + rtree_t * root = (rtree_t *)xcalloc(1,sizeof(rtree_t)); /* clone the two subtrees */ root->left = utree_rtree(lca->back->next->back); @@ -435,7 +435,7 @@ rtree_t * utree_crop(utree_t * lca) rtree_t * utree_convert_rtree(utree_t * outgroup) { - rtree_t * root = (rtree_t *)xmalloc(sizeof(rtree_t)); + rtree_t * root = (rtree_t *)xcalloc(1,sizeof(rtree_t)); root->left = utree_rtree(outgroup); root->right = utree_rtree(outgroup->back); @@ -477,11 +477,11 @@ static utree_t ** utree_tipstring_nodes(utree_t * root, if (tipstring[i] == ',') commas_count++; - utree_t ** node_list = (utree_t **)xmalloc((size_t)utree_tip_count * + utree_t ** node_list = (utree_t **)xcalloc(1,(size_t)utree_tip_count * sizeof(utree_t *)); utree_query_tipnodes(root, node_list); - utree_t ** out_node_list = (utree_t **)xmalloc((commas_count+1) * + utree_t ** out_node_list = (utree_t **)xcalloc(1,(commas_count+1) * sizeof(utree_t *)); /* create a hashtable of tip labels */ @@ -547,7 +547,7 @@ static utree_t * utree_lca(utree_t ** tip_nodes, utree_t ** path; /* allocate a path */ - path = (utree_t **)xmalloc((size_t)utree_tip_count * + path = (utree_t **)xcalloc(1,(size_t)utree_tip_count * sizeof(utree_t **)); /* mark all tip nodes */