• Add a vector pattern_weights of pattern weights to pll_partition_t structure
• Add a pll_set_pattern_weights() function
• Modify pll_compute_root_loglikelihood and pll_compute_edge_loglikelihood() to accommodate pattern_weights

Setup wiki documentation for the current milestone:

• Basics about libpll. Partition structure, how to setup a partition.
• Explain how tip vectors should look like in a typical phylogenetic analysis (given a solid and a degenerate character)
• The two ways of setting tip vectors (set clv manually or pass sequence and map)
• How to parse FASTA files (complete or partial parse)
• How to parse newick files (rooted and unrooted)
• How to setup frequencies, rate heterogeneity, substitution parameters
• Explain what is a rate matrix (R), an instantaneous rate matrix (Q-matrix) and a transition probability matrix (P-matrix). Show all the gory details of creating a p-matrix (normalization at diagonal etc). Give a description of the GTR.
• How to compute the p-matrices with the library (update_prob_matrices) given the branch lengths
• How to compute the CLV vectors (pll_update_partials)
• How to evaluate the log likelihood of a rooted tree (pll_compute_root_loglikelihood)
• How to evaluate the log likelihood of an unrooted tree (pll_compute_edge_loglikelihood)
• How to account for invariant sites
• Document all protein models, give references
• Document the demo examples

Perhaps now it would be a good idea to decide the data types of the following pll_partition_t variables:

My justification for the proposed changes is that int (or unsigned int) could be an insufficient range for large phylogenomic datasets. Changing tips, clv_buffers, sites, prob_matrices and scale_buffers would enable the analysis of datasets comprising of more than 4,294,967,295 sites.
The drawback is that the memory requirements will double.

Also, if the proposed changed are done, it would also make sense to change scale_buffers, invariant to unsigned long ** and unsigned long * respectively, because the new data type (unsigned long) would then be of the same size (8 bytes) as a CLV entry (double), and this would improve vectorization efficiency.

What do you think?

My idea is to test several strategies, although we should discuss them before implementing.

• Option for disabling scaling completely
• Option for scaling always
• Option for scaling 'sometimes'

Remove all warnings when using clang together with the switches -Weverything and -Wno-padding, except from the flex/bison generated code.

Update pll_rtree_traverse and pll_utree_traverse in a way to allow a customized traversal using a callback function. For example, to be able to perform a partial postorder traversal based on the validity of CLV at inner nodes. That is if nodes are correctly oriented, then the subtree rooted at those nodes is not traversed.

Check if:

• tips > 0
• clv_buffers > 0
• states > 0
• sites > 0
• rate_matrices > 0
• prob_matrices > 0
• rate_cats > 0

We should do profiling runs for full tree traversals to determine how much time is spent in P-matrix and CLV updates and if it would be worthwhile to further optimize the p matrix calcs.

Implement two functions. One for reading a complete FASTA file in memory, and another for iteratively getting the reads. Will copy the functions from SALT.

PMatrices are allocated according to the padded number of states:

partition->pmatrix[i] = pll_aligned_alloc(states_padded * states_padded *
                                              rate_cats * sizeof(double),
                                              partition->alignment);

However, when computing the likelihood it does not account for the padding. It works well for 4 and 20 states, but it might be a problem if padding is applied.

pmatrix = partition->pmatrix[matrix_index];
for (i = 0; i < partition->rate_cats; ++i)
{
  for (j = 0; j < partition->states; ++j)
  {
     ....
     pmatrix += states;
  }
}

We should check if this problem happens somewhere else in the code.

I also wondered if branch length optimization can be accelerated, if, instead of initial default branch lengths, we compute a good initial value for each branch separately using parsimony

Implement a pll_update_operations_minimal that takes a list of operations and an integer N and re-assigns parent_clv_index, child1_clv_index and child2_clv_index of all nodes in order to use the minimal amount C of CLVs, which is, in the worst case, log(n)+2.

Another useful function would be: Given a tree, and not requiring a partition instance, compute the minimal number C, such that this amount may be passed later on to the pll_partition_create function, in case we are working with very large datasets.

to reflect the modified rates/freqs.

Implement a function that given the alpha shape parameter of the Gamma distirbution, generates the specified number of discretized rate categories.

I believe that implementing separate, sepcifically optimized functions for the tips should yield huge performance improvements, keep in mind that CLVs above the tips make up for about 50% of all CLV updates

Several issues:

• instead of the source code there is a binary executable
• missing documentation (overall description) as a README.md file in the example dir.

Offer core LLPLL functions that do not require passing a partition data type as argument, but just C arrays and variables, this will facilitate integration with ExaML.

The Makefile doesn't work on Mac OS X because -soname is not a valid compiler option. I think, it should be substituted with -install_name on OS X (if that's possible).

Extend/adapt partition datatype such that it can store the overall length of the partition and also the number of sites it does computations on, for instance, when the partition has been split up across one or more processors.

Just a few comments on the vectorization if you want to make it generic:

1. Keep in mind that for a truly generic implementation there might be data-types that are not multiples of the vector width, e.g., a 7-state model (this exists for RNA secondary structure models!)
2. Thus, you may have to use uload for vector words and also need to add additional loops outside the vectorized loops for handling those indexes that don't fit into a vector word.
3. Alternatively, you can of course use padding, but that may be more complicated to implement

It would be nice to have test cases for Win{32,64} and Mac.

For example, something like ./runtest.py win64, which will use a cross-compiler (i.e. mingw) to build a PE file.

If the shape of the Gamma distribution is low enough, it might happen that the likelihood scores of different categories require different scaling. That means that at some point it is very likely that some per-category likelihoods go under precision. We can keep one individual scaler for each site and category (thus this requires 'c' times more memory for the scalers, where 'c' is the number of different categories), and the likelihood computation is modified as follows:

s_{i,j} is the scaler for site 'i' and category 'j'.
K is the likelihood scaling threshold

Apart from the additional memory required, which is not excessive compared to the CLVs, we need to check whether each term with scaling involved in the sum (i.e., s_{i,j} > s_{i,min}) may or may not go out of precision. If we assume that K is fixed to 2^256, then a scaling >= 2 is enough for discarding the term. However, we might want to implement this in a more generic way and decide which would be the highest supported difference in the scalers. For example, we could change the scaling threshold at some point for testing purposes.

A segmentation fault occurs when scaling is disabled for the parent node, i.e. operations[i].parent_scaler_index = PLL_SCALE_BUFFER_NONE, in which case the scaler array is set to NULL. If, however, scaling is required, the function update_partials() will attempt to fill the scaler array.

Depending on the starting node we pass, pll_traverse_utree will set the tip CLV indices in the operations structure on a visit-first order (post-order traversal). Therefore, the indices may (and will) not correspond to the indices used in pll_set_tip_states or pll_set_tip_clv when setting tip CLVs using sequences.

Proposed solution

Add a CLV index variable in the pll_utree_t structure for tip nodes, which must be filled, for instance when setting tip CLVs from a tree traversal. The function pll_traverse_utree will then uses those CLV indices, instead of (wrongly) enumerating on a first-visit basis.

A CLV update function, that only updates between indices i and j, where 0 <= i < j < num_sites

Motivation:
Suppose you're calculating the inner node CLV between two leafs with sequences

-  -  -  -  -  A  T  A  G  C  T  -  -  -  -  -  - 
-  -  -  -  -  A  T  T  G  C  T  A  -  -  -  -  - 
0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16

then ideally we would like to skip the computation for the outer parts where neither sequence has valid data.

Compare performance with PLL, right now, looking at the code (nicely written by the way) I can still not believe that the performance differences are so small ....

I am just working on developing a PThreads-version of ExaML.

In this context I observed the following:

1. We mainly need to implement a function analogous to the MPI allreduce operation in PThreads.
2. Thinking of auto-tuned likelihood implementations, we may also have to consider aut-tuned barrier implementations. The reason is that barrier performance and types required will probably be very different if we just have one partition or several partitions that are mostly distributed monolithically to threads.

This could be a nice technical paper in the context of a master's thesis.

Modify libpll such that compilation with mingw is possible for Win64

• We should discuss how to proceed on this

1. extend pll_set_tip_states() with a map parameter that will specify the conversion of characters to states. Currently the maps map_nt and map_aa are implicitely used if states=4 and states=20 respectively.
2. add a function pll_set_tip_clv() which can be used instead of pll_set_tip_states for setting the tip clv manually.

• Evaluating the log-likelihood of a rooted tree
• Evaluating the log-likelihood of an unrooted tree
• Perhaps one of them with DNA data, the other with AA
• Show invariant sites usage
• Include parsing of fasta/newick files
• Show both methods for setting the tip CLVs.

for testing purposes

a) Parallel Newton-Raphson optimization

Instead of iteratively optimize and apply the branch lengths, optimize each of them and store the optimized value (proposal) in a buffer. Once all branches have been optimized, applied them and perform a full traversal.

The advantage of this approach is that branches can be optimized in parallel and in an arbitrary order, leading always to the same branch lengths (reproducible result regardless how they are computed). There is no need to update CLVs/scalers between branches, so the optimization of the set of branches in a single core is faster for each iteration. However, the convergence of the branch length optimization process will probably take a longer time, and perhaps it is prone to get stuck in local optima easier, or to end with a lower global likelihood.

b) Multiple-branch BFGS

Optimize simultaneously small subsets of branches (e.g., 5 to 10) using BFGS. Instead of approximating the gradients, they can be calculated using the first derivative.

There are missing amino acid characters in pll_map_fasta

E,I,L,...

For the 2000 taxon dataset it seems that for some models we are not getting the correct log-likelihood. Specifically:

DAYHOFF

DCMUT

CPREV

VT

BLOSUM62

MTZOA

PMB

HIVB

JTTDCMUT

Please add a wiki documentation section with the steps on how to run the testing framework. It's always a painful thing to make it run, especially with the increasing number of header files.

Or is the library intended to only be compiled with asserts enabled?

The really anal approach is to write:

assert(root != NULL);

as:

assert(root != NULL);
if (root == NULL)
{
  *some code here*
}

IMHO, this is too anal. Just thought I'd ask...

The following API function names will be modified

Add an array of integers to the partition structure which will hold the frequencies of each site. It will be initialized to all ones at the pll_partition_create() call, and then the user will be able to set it by calling pll_set_site_weights().

The prototype of the function should be:

PLL_EXPORT void pll_set_site_weights(pll_partition_t * partition,
                                     const unsigned int * weights);

Apart from minor modifications, what is needed is:

• Append a map_bin in maps.c

1. Create a simple rooted tree structure.
2. Copy the flex/bison parser from fastdate.

Hi guys,

Check the latest RAxML commit (optimizeModel.c) where I fixed a numerical issue in BFGS:

stamatak/standard-RAxML@9a02452

You need to be very careful about specifying the lower bounds for the variables you wish to optimize, otherwise you'll spend a day debugging at least. The lower bounds should not be smaller than the error margin!

Alexis

Three things

• missing documentation (overall description) as a README.md file in the example dir.
• keep consistent directory structures (examples dont have their own src).
• currently broken with the TIP-TIP addition.

Need to pass a map to the pll_create_partition. However, the map is decided after the partition is created. I guess it's a simple fix to move the switch deciding the map before the pll_create_partition.

But i'm not well oriented in this example, I guess you can do it quickly @ddarriba

Set the matrices for the popular models and append an interface function for setting the parameters in the partition.

Multiply Q with square root of the diagonal matrix of frequencies (pi) and its inverse, i.e.

sqrt(pi) Q (sqrt(pi))^-1

1. Create raxm-like tree structure
2. Modify the rooted tree flex/bison parser