Fix a memory leak when branch lengths are too close together. Fixes #59

Have cafe_tree_set_birthdeath return its matrix cache for easier management Show a better error message for an invalid family
hahnlab · Mar 7, 2019 · bce1cef · bce1cef
1 parent b5aa0eb
commit bce1cef
Show file tree

Hide file tree

Showing 10 changed files with 48 additions and 19 deletions.
diff --git a/cafe/cafe_commands.cpp b/cafe/cafe_commands.cpp
@@ -1787,7 +1787,12 @@ void tree_set_branch_lengths(pCafeTree pcafe, std::vector<int> lengths)
 			fprintf(stderr, "ERROR: the branch length of node %d is not changed\n", (int)i);
 		}
 	}
-	if (probability_cache) cafe_tree_set_birthdeath(pcafe, probability_cache->maxFamilysize);
+    if (probability_cache)
+    {
+        auto bd_cache = cafe_tree_set_birthdeath(pcafe, probability_cache->maxFamilysize);
+        free_cache_keep_matrices(bd_cache);
+
+    }
 
 }
 

diff --git a/cafe/cafe_main.c b/cafe/cafe_main.c
@@ -322,8 +322,7 @@ void reset_birthdeath_cache(pCafeTree tree, int k_value, family_size_range* rang
 	{
 		birthdeath_cache_array_free(probability_cache);
 	}
-	probability_cache = birthdeath_cache_init(MAX(range->max, range->root_max), &cache);
-	cafe_tree_set_birthdeath(tree, probability_cache->maxFamilysize);
+    probability_cache = cafe_tree_set_birthdeath(tree, MAX(range->max, range->root_max));
 }
 
 /**************************************************************************

diff --git a/cafe/cafe_shell.c b/cafe/cafe_shell.c
@@ -317,7 +317,7 @@ void cafe_shell_set_branchlength(pCafeParam param, int max_family_size)
 		}
 	}
 
-    cafe_tree_set_birthdeath(param->pcafe, max_family_size);
+    free_cache_keep_matrices(cafe_tree_set_birthdeath(param->pcafe, max_family_size));
 }
 
 

diff --git a/cafe/cafe_tree.c b/cafe/cafe_tree.c
@@ -371,8 +371,9 @@ void node_set_birthdeath_matrix(pCafeNode pcnode, pBirthDeathCacheArray cache, i
 
 }
 
-void add_key(pArrayList arr, double branchlength, double lambda, double mu)
+void add_key(pArrayList arr, double d_branchlength, double lambda, double mu)
 {
+    int branchlength = (int)d_branchlength;
     for (int i = 0; i < arr->size; ++i)
     {
         struct BirthDeathCacheKey *key = (struct BirthDeathCacheKey *)arraylist_get(arr, i);
@@ -444,10 +445,20 @@ void gather_keys(pTree ptree, pTreeNode ptnode, va_list ap1)
     get_keys_from_node((pCafeNode)ptnode, arr, pcafe->k);
 }
 
+void free_cache_keep_matrices(pBirthDeathCacheArray bd_cache)
+{
+    // free the cache without deleting the matrices
+    void** keys = NULL;
+    hash_table_get_keys(bd_cache->table, &keys);
+    free(keys);
+    hash_table_delete(bd_cache->table);
+    memory_free(bd_cache);
+}
+
 /**
 *	Set each node's birthdeath matrix based on its values of branchlength, lambdas, and mus
 **/
-void cafe_tree_set_birthdeath(pCafeTree pcafe, int max_family_size)
+pBirthDeathCacheArray cafe_tree_set_birthdeath(pCafeTree pcafe, int max_family_size)
 {
     pArrayList arr = arraylist_new(40);
     tree_traveral_prefix((pTree)pcafe, gather_keys, arr);
@@ -466,12 +477,9 @@ void cafe_tree_set_birthdeath(pCafeTree pcafe, int max_family_size)
 
 	tree_traveral_prefix((pTree)pcafe, do_node_set_birthdeath, bd_cache);
 
-    // free the cache without deleting the matrices
-    void** keys = NULL;
-    hash_table_get_keys(bd_cache->table, &keys);
-    free(keys);
-    hash_table_delete(bd_cache->table);
-    memory_free(bd_cache);
+    arraylist_free(arr, NULL);
+
+    return bd_cache;
 }
 
 void cafe_tree_node_copy(pTreeNode psrc, pTreeNode pdest)
@@ -513,8 +521,8 @@ pCafeTree cafe_tree_split(pCafeTree pcafe, int idx )
 	{
 		__cafe_tree_copy_new_fill(pcafe,psub);
 	}
-	cafe_tree_set_birthdeath(pcafe, probability_cache->maxFamilysize);
-	cafe_tree_set_birthdeath(psub, probability_cache->maxFamilysize);
+	free_cache_keep_matrices(cafe_tree_set_birthdeath(pcafe, probability_cache->maxFamilysize));
+	free_cache_keep_matrices(cafe_tree_set_birthdeath(psub, probability_cache->maxFamilysize));
 	return psub;
 }
 

diff --git a/cafe/lambda.cpp b/cafe/lambda.cpp
@@ -646,11 +646,27 @@ double* cafe_best_lambda_by_fminsearch(pCafeParam param, int lambda_len, int k)
 	return param->input.parameters;
 }
 
+void cafe_tree_string_likelihood(pString pstr, pPhylogenyNode pnode)
+{
+    pCafeNode pcnode = (pCafeNode)pnode;
+    double val = *std::max_element(pcnode->likelihoods, pcnode->likelihoods+50);
+    if (pnode->name) string_fadd(pstr, "%s", pnode->name);
+    string_fadd(pstr, "<%f>", val);
+}
+
 posterior compute_posterior(pCafeFamilyItem pitem, pCafeTree pcafe, const std::vector<double>& prior_rfsize)
 {
     posterior result;
     compute_tree_likelihoods(pcafe);
 
+#if 0
+    // print probabiities in Newick format, for debugging
+    cout << "Family: " << pitem->id;
+    pString pstr = phylogeny_string_newick((pTree)pcafe, cafe_tree_string_likelihood, PS_SKIP_BL);
+    cout << pstr->buf << endl;
+    string_free(pstr);
+#endif
+
     double *likelihood = get_likelihoods(pcafe);		// likelihood of the whole tree = multiplication of likelihood of all nodes
 
     result.max_likelihood = __max(likelihood, pcafe->rfsize);			// this part find root size condition with maxlikelihood for each family			

diff --git a/cafe/likelihood_ratio.cpp b/cafe/likelihood_ratio.cpp
@@ -89,7 +89,7 @@ double __cafe_lhr_get_likelihood_for_diff_lambdas(pCafeParam param, int idx, int
 		memcpy(param->lambda, lambda_cache[t], sizeof(double)*param->num_lambdas);
         cafe_shell_set_lambdas(param, param->lambda);
 		probability_cache = PBDC[t];
-		cafe_tree_set_birthdeath(param->pcafe, probability_cache->maxFamilysize);
+		free_cache_keep_matrices(cafe_tree_set_birthdeath(param->pcafe, probability_cache->maxFamilysize));
 	}
 	int i;
   pCafeFamilyItem pitem = (pCafeFamilyItem)param->pfamily->flist->array[idx];

diff --git a/cafe/pvalue.cpp b/cafe/pvalue.cpp
@@ -28,7 +28,7 @@ void check_cache_and_compute_likelihoods(pCafeTree pTree, int max, pBirthDeathCa
 		{
 			int remaxFamilysize = MAX(range.max, range.root_max);
 			birthdeath_cache_resize(cache, remaxFamilysize);
-			cafe_tree_set_birthdeath(pTree, cache->maxFamilysize);
+			free_cache_keep_matrices(cafe_tree_set_birthdeath(pTree, cache->maxFamilysize));
 		}
 		else
 		{

diff --git a/configure.ac b/configure.ac
@@ -1,4 +1,4 @@
-AC_INIT(cafe, [4.2], [email protected],, hahnlab.github.io/CAFE)
+AC_INIT(cafe, [4.2.1], [email protected],, hahnlab.github.io/CAFE)
 
 m4_include([m4/ax_lib_readline.m4])
 

diff --git a/libtree/family.h b/libtree/family.h
@@ -178,5 +178,6 @@ struct tagCafeParam
 
 extern void thread_run_with_arraylist(int numthreads, void* (*run)(void*), pArrayList pal );
 // cafe tree
-extern void cafe_tree_set_birthdeath(pCafeTree pcafe, int max_family_size);
+extern pBirthDeathCacheArray cafe_tree_set_birthdeath(pCafeTree pcafe, int max_family_size);
+void free_cache_keep_matrices(pBirthDeathCacheArray cache);
 #endif
diff --git a/tests/command_tests.cpp b/tests/command_tests.cpp
@@ -93,7 +93,7 @@ TEST(CommandTests, cafe_cmd_date)
 	char outbuf[10000];
 	globals.param.flog = fmemopen(outbuf, 999, "w");
 	cafe_cmd_date(globals, tokens);
-	STRCMP_CONTAINS("2018", outbuf);	// this will start to fail in 2018
+	STRCMP_CONTAINS("2019", outbuf);	// this will start to fail in 2020
 	fclose(globals.param.flog);
 }
-Original file line number
+Diff line change
@@ Expand Up @@
     		}
     	}
-        cafe_tree_set_birthdeath(param->pcafe, max_family_size);
+        free_cache_keep_matrices(cafe_tree_set_birthdeath(param->pcafe, max_family_size));
     }
@@ Expand Down @@