Skip to content

Commit

Permalink
added wrapper tree structure, updated examples, minor fixes in unit t…
Browse files Browse the repository at this point in the history 
…ests and examples - issue #126
  • Loading branch information
flouris authored and flouris committed Mar 30, 2017
1 parent 5558192 commit 61b78f7
Show file tree
Hide file tree
Showing 19 changed files with 944 additions and 642 deletions.
140 changes: 76 additions & 64 deletions examples/lg4/lg4.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,42 +29,49 @@

static void fatal(const char * format, ...) __attribute__ ((noreturn));

static void * xmalloc(size_t size)
{
void * t;
t = malloc(size);
if (!t)
fatal("Unable to allocate enough memory.");

return t;
}

static char * xstrdup(const char * s)
{
size_t len = strlen(s);
char * p = (char *)xmalloc(len+1);
return strcpy(p,s);
}

/* a callback function for performing a full traversal */
static int cb_full_traversal(pll_utree_t * node)
static int cb_full_traversal(pll_unode_t * node)
{
return 1;
}

static void set_missing_branch_length_recursive(pll_utree_t * tree,
double length)
/* branch lengths not present in the newick file get a value of 0.000001 */
static void set_missing_branch_length(pll_utree_t * tree, double length)
{
if (tree)
{
/* set branch length to default if not set */
if (!tree->length)
tree->length = length;

if (tree->next)
{
if (!tree->next->length)
tree->next->length = length;
unsigned int i;

if (!tree->next->next->length)
tree->next->next->length = length;
for (i = 0; i < tree->tip_count; ++i)
if (!tree->nodes[i]->length)
tree->nodes[i]->length = length;

set_missing_branch_length_recursive(tree->next->back, length);
set_missing_branch_length_recursive(tree->next->next->back, length);
}
for (i = tree->tip_count; i < tree->tip_count + tree->inner_count; ++i)
{
if (!tree->nodes[i]->length)
tree->nodes[i]->length = length;
if (!tree->nodes[i]->next->length)
tree->nodes[i]->next->length = length;
if (!tree->nodes[i]->next->next->length)
tree->nodes[i]->next->next->length = length;
}
}

/* branch lengths not present in the newick file get a value of 0.000001 */
static void set_missing_branch_length(pll_utree_t * tree, double length)
{
set_missing_branch_length_recursive(tree, length);
set_missing_branch_length_recursive(tree->back, length);
}

static void fatal(const char * format, ...)
{
va_list argptr;
Expand All @@ -84,7 +91,7 @@ int main(int argc, char * argv[])
double * branch_lengths;
pll_partition_t * partition;
pll_operation_t * operations;
pll_utree_t ** travbuffer;
pll_unode_t ** travbuffer;

/* we accept only two arguments - the newick tree (unrooted binary) and the
alignment in the form of FASTA reads */
Expand All @@ -93,13 +100,16 @@ int main(int argc, char * argv[])

/* parse the unrooted binary tree in newick format, and store the number
of tip nodes in tip_nodes_count */
pll_utree_t * tree = pll_utree_parse_newick(argv[1], &tip_nodes_count);
pll_utree_t * tree = pll_utree_parse_newick(argv[1]);
if (!tree)
fatal("%s", pll_errmsg);

tip_nodes_count = tree->tip_count;

/* fix all missing branch lengths (i.e. those that did not appear in the
newick) to 0.000001 */
set_missing_branch_length(tree, 0.000001);


/* compute and show node count information */
inner_nodes_count = tip_nodes_count - 2;
nodes_count = inner_nodes_count + tip_nodes_count;
Expand All @@ -117,24 +127,23 @@ int main(int argc, char * argv[])
char * newick = pll_utree_export_newick(tree);
printf("%s\n", newick);
free(newick);
// */

/* obtain an array of pointers to tip nodes */
pll_utree_t ** tipnodes = (pll_utree_t **)calloc(tip_nodes_count,
sizeof(pll_utree_t *));
pll_utree_query_tipnodes(tree, tipnodes);
*/

/* create a libc hash table of size tip_nodes_count */
hcreate(tip_nodes_count);

/* populate a libc hash table with tree tip labels */
unsigned int * data = (unsigned int *)malloc(tip_nodes_count *
sizeof(unsigned int));
unsigned int * data = (unsigned int *)xmalloc(tip_nodes_count *
sizeof(unsigned int));
for (i = 0; i < tip_nodes_count; ++i)
{
data[i] = tipnodes[i]->clv_index;
data[i] = tree->nodes[i]->clv_index;
ENTRY entry;
entry.key = tipnodes[i]->label;
#ifdef __APPLE__
entry.key = xstrdup(tree->nodes[i]->label);
#else
entry.key = tree->nodes[i]->label;
#endif
entry.data = (void *)(data+i);
hsearch(entry, ENTER);
}
Expand Down Expand Up @@ -229,7 +238,6 @@ int main(int argc, char * argv[])

/* we no longer need these two arrays (keys and values of hash table... */
free(data);
free(tipnodes);

/* ...neither the sequences and the headers as they are already
present in the form of probabilities in the tip CLVs */
Expand All @@ -243,17 +251,21 @@ int main(int argc, char * argv[])

/* allocate a buffer for storing pointers to nodes of the tree in postorder
traversal */
travbuffer = (pll_utree_t **)malloc(nodes_count * sizeof(pll_utree_t *));
travbuffer = (pll_unode_t **)xmalloc(nodes_count * sizeof(pll_unode_t *));

branch_lengths = (double *)malloc(branch_count * sizeof(double));
matrix_indices = (unsigned int *)malloc(branch_count * sizeof(unsigned int));
operations = (pll_operation_t *)malloc(inner_nodes_count *
sizeof(pll_operation_t));
branch_lengths = (double *)xmalloc(branch_count * sizeof(double));
matrix_indices = (unsigned int *)xmalloc(branch_count * sizeof(unsigned int));
operations = (pll_operation_t *)xmalloc(inner_nodes_count *
sizeof(pll_operation_t));

/* compute a partial traversal starting from the randomly selected
inner node */
/* compute a partial postorder traversal starting from the inner node that
was the the root of the parsed 'unrooted' binary tree */

pll_unode_t * root = tree->nodes[tip_nodes_count+inner_nodes_count-1];
unsigned int traversal_size;
if (!pll_utree_traverse(tree,

if (!pll_utree_traverse(root,
PLL_TREE_TRAVERSE_POSTORDER,
cb_full_traversal,
travbuffer,
&traversal_size))
Expand Down Expand Up @@ -323,11 +335,11 @@ int main(int argc, char * argv[])
frequency vector is to be used */

double logl = pll_compute_edge_loglikelihood(partition,
tree->clv_index,
tree->scaler_index,
tree->back->clv_index,
tree->back->scaler_index,
tree->pmatrix_index,
root->clv_index,
root->scaler_index,
root->back->clv_index,
root->back->scaler_index,
root->pmatrix_index,
params_indices,
NULL);

Expand All @@ -351,8 +363,8 @@ int main(int argc, char * argv[])
double weights[4] = { 0.209224645, 0.224707726, 0.277599198, 0.288468431 };
double rates[4] = { 0.498991136, 0.563680734, 0.808264032, 1.887769458 };

pll_set_category_rates (partition, rates);
pll_set_category_weights (partition, weights);
pll_set_category_rates(partition, rates);
pll_set_category_weights(partition, weights);

/* update transition probability matrices */
pll_update_prob_matrices(partition,
Expand All @@ -374,11 +386,11 @@ int main(int argc, char * argv[])
pll_update_partials(partition, operations, ops_count);

logl = pll_compute_edge_loglikelihood(partition,
tree->clv_index,
tree->scaler_index,
tree->back->clv_index,
tree->back->scaler_index,
tree->pmatrix_index,
root->clv_index,
root->scaler_index,
root->back->clv_index,
root->back->scaler_index,
root->pmatrix_index,
params_indices,
NULL);

Expand All @@ -391,11 +403,11 @@ int main(int argc, char * argv[])
printf("\nLog-L (LG4X): %f\n", logl);

logl = pll_compute_edge_loglikelihood(partition,
tree->back->clv_index,
tree->back->scaler_index,
tree->clv_index,
tree->scaler_index,
tree->pmatrix_index,
root->back->clv_index,
root->back->scaler_index,
root->clv_index,
root->scaler_index,
root->pmatrix_index,
params_indices,
NULL);

Expand Down
20 changes: 17 additions & 3 deletions examples/load-utree/load-utree.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -42,18 +42,26 @@ pll_utree_t * load_tree_unrooted(const char * filename,

if (!(rtree = pll_rtree_parse_newick(filename, tip_count)))
{
if (!(utree = pll_utree_parse_newick(filename, tip_count)))
if (!(utree = pll_utree_parse_newick(filename)))
{
fprintf(stderr, "%s\n", pll_errmsg);
return NULL;
}
else
{
*tip_count = utree->tip_count;
}
}
else
{
utree = pll_rtree_unroot(rtree);

pll_unode_t * root = utree->nodes[utree->tip_count+utree->inner_count-1];

/* optional step if using default PLL clv/pmatrix index assignments */
pll_utree_reset_template_indices(utree, *tip_count);
pll_utree_reset_template_indices(root, *tip_count);

pll_rtree_destroy(rtree,NULL);
}

return utree;
Expand All @@ -70,7 +78,13 @@ int main(int argc, char * argv[])
if (!utree)
fatal("Tree must be a rooted or unrooted binary.");

char * newick = pll_utree_export_newick(utree);
/* select a random inner node */
long int r = random() % utree->inner_count;
pll_unode_t * root = utree->nodes[utree->tip_count + r];

/* export the tree to newick format with the selected inner node as the root
of the unrooted binary tree */
char * newick = pll_utree_export_newick(root);

printf("%s\n", newick);

Expand Down
Loading

0 comments on commit 61b78f7

Please sign in to comment.