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

@@ -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

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])

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);

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;

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);
 }

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);

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 *));

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;

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)
 {

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 */