Msa refactor (#17)

* Initial PIP branch length optimisation

Initial impl of branch length optimisation under PIP

* Refactoring: tree stores all nodes in one list

Refactoring the tree structure to simplify access/creation:
now all nodes are stored in a single vector instead of leaves and internals separately.
Additional modifications to make all other structures work with the new version.

* Removing get_ prefixes

* PhyloInfo refactor to include builder

Big refactoring for PhyloInfo to create the object using a builder rather than random
functions.
Refactor of the alphabet module (renamed from "sequences") to have own alphabet not relying
on bio::alphabets.

* Refactoring alphabets

Refactoring alphabets:
    made checking whether something is a word faster using a HashSet;
    removed all small letters, only using uppercase now;
    rename nucleotide/aa sets to improve readability and remove mixing
DNA/nucleotide and Protein/aminoacid.

* Alignment structure overhaul

Alignment structure overhaul to include both pairwise and a complete alignment and
to work with a local and global representation. Added tests for all the representations.

* Creating proper alignments from aligned sequences

Creating a proper aligment from supplied aligned sequences, reconstructing the local
alignments per node, saving leaf maps for later, removing gaps from saved leaf sequences.

* Tree refac: adding accessors, fields more private

Minor refactoring of the tree class: made a few fields more private (nodes.
prostorder, preorder) to discourage direct access (for myself), added accessors
for nodes and node ids, and len() and is_empty() methods.

* Fix duplicate code for Node/NodeIdx Display

Removed duplicate code from the Display impl in Node and NodeIdx.
Added explicit Debug to NodeIdx for easier debugging.

* Modify usize::from for NodeIdx to reduce deref

Modified usize::from for both &NodeIdx and NodeIdx to reduce dereferencing.

* Refactor alignment to use builder, class for seqs

Refactored alignment creation into a builder struct, the internal representation
can be built from either a given MSA or from a map of internal nodes.
Now the alignment class also stores the ungapped sequences and, for a current alignment,
the maps of leaf sequences to gaps/characters.
Sequences have a separate struct for storage that has several methods to simplify access.

* Changes throughout to use Sequences and Alignment

Massive refactoring to use Sequences and Alignment structs throughout.
PhyloInfo can only be generated at the moment with aligned sequences,
but leaf encodings are stored for sequences without gaps and are used to compute
conditional probabilities at the tips with leaf maps from the alignment.

* Add test for gaps treated as X for non-gap models

Added a test case checking whether gaps are indeed treated as ambiguous characters when
a non-gap model is given.

* Added missing test data file

* Refactor AlignBuilder to separate file

Refactored AlignmentBuilder to a separate file to clean up code.
Rename alignment_tests.rs to tests.rs to remove repetition in names.

* Sequences to uppercase when reading from file

Sequences are now transformed into uppercase as they are read from file, not
afterwards + doctest for this.
Signature for write_sequences_to_file now is &Pathbuf, not Pathbuf to avoid copying.

* PhyloInfoBuilder to new module + docs

Moved PhyloInfoBuilder to a separate file and added docs + doctests.

* Docs for PhyloInfo and move leaf_encoding creation

Added docs for PhyloInfo struct.
Moved leaf_encoding creation to a more logical place within PhyloInfo
rather than in the builder.

* Alphabet is defined in Sequences and only there

The sequence struct now defines the Alphabet for the whole setup and that is the only
place where the Alphabet exists. The Alphabet is used to compute leaf sequence encodings
to simplify likelihood computation.

* move alphabets_tests -> tests

* Remove unnecessary method

* Remove unnecessary GapHandling enum

Removed the GapHandling enum for good, the evolutionary models themselves
now handle getting the correct encodings for gaps and characters.

* Removed stale files

* Remove vscode settings file

* Remove vscode files from git

Author: junniest

.gitignore

@@ -18,3 +18,6 @@ Cargo.lock
 
 # Mac files
 .DS_Store
+
+# VSCode files
+.vscode/

.vscode/settings.json

@@ -1,8 +1,8 @@
 >A
-aAAAAAA
+aAAAAAAa
 >B
-ccccCCC
+ccccCCCc
 >C
-gGggggg
+gGgggggG
 >D
-ttttttt
+tttttttT

phylo/data/sequences_DNA_small.fasta

@@ -0,0 +1 @@
+(((A:1.0,B:2.0):3.0,(C:3.5,D:0.5):5.25):1.0,(E:4.0,F:5.0):0.5);

phylo/data/tree_diff_branch_lengths_3.newick

@@ -0,0 +1,121 @@
+use std::collections::HashMap;
+
+use anyhow::bail;
+
+use crate::align;
+use crate::alignment::{Alignment, InternalMapping, Mapping, PairwiseAlignment, Sequences};
+use crate::tree::{NodeIdx, NodeIdx::Internal as Int, NodeIdx::Leaf, Tree};
+use crate::Result;
+
+pub struct AlignmentBuilder<'a> {
+    tree: &'a Tree,
+    seqs: Sequences,
+    node_map: InternalMapping,
+}
+
+impl<'a> AlignmentBuilder<'a> {
+    pub fn new(tree: &'a Tree, seqs: Sequences) -> AlignmentBuilder<'a> {
+        AlignmentBuilder {
+            tree,
+            seqs,
+            node_map: InternalMapping::new(),
+        }
+    }
+
+    pub fn msa(mut self, msa: InternalMapping) -> Self {
+        self.node_map = msa;
+        self
+    }
+
+    pub fn build(self) -> Result<Alignment> {
+        if self.node_map.is_empty() {
+            if self.seqs.aligned {
+                self.build_from_seqs()
+            } else {
+                self.build_from_unaligned()
+            }
+        } else {
+            self.build_from_map()
+        }
+    }
+
+    /// This assumes that the tree structure matches the alignment structure and that the sequences are aligned.
+    fn build_from_seqs(self) -> Result<Alignment> {
+        let msa_len = self.seqs.msa_len();
+        let mut stack = HashMap::<NodeIdx, Mapping>::with_capacity(self.tree.len());
+        let mut msa = InternalMapping::with_capacity(self.tree.n);
+        for node in self.tree.postorder.iter() {
+            match node {
+                Int(_) => {
+                    let childs = self.tree.children(node);
+                    let map_x = stack[&childs[0]].clone();
+                    let map_y = stack[&childs[1]].clone();
+                    stack.insert(*node, Self::stack_maps(msa_len, &map_x, &map_y));
+                    msa.insert(*node, Self::clear_common_gaps(msa_len, &map_x, &map_y));
+                }
+                Leaf(_) => {
+                    let seq = self.seqs.get_by_id(self.tree.node_id(node)).seq();
+                    stack.insert(*node, align!(seq).clone());
+                }
+            }
+        }
+        let leaf_maps = stack
+            .iter()
+            .filter_map(|(idx, map)| match idx {
+                Leaf(_) => Some((*idx, map.clone())),
+                _ => None,
+            })
+            .collect();
+        Ok(Alignment {
+            seqs: self.seqs.without_gaps(),
+            leaf_map: leaf_maps,
+            node_map: msa,
+        })
+    }
+
+    fn build_from_unaligned(self) -> Result<Alignment> {
+        // TODO: use parsimony to align the sequences.
+        bail!("Unaligned sequences are not yet supported.")
+    }
+
+    /// This assumes that the tree structure matches the alignment structure.
+    fn build_from_map(self) -> Result<Alignment> {
+        let mut alignment = Alignment {
+            seqs: Sequences::new(Vec::new()),
+            leaf_map: HashMap::new(),
+            node_map: self.node_map,
+        };
+        let leaf_map = alignment.compile_leaf_map(self.tree.root, self.tree)?;
+        alignment.leaf_map = leaf_map;
+        alignment.seqs = self.seqs.without_gaps();
+        Ok(alignment)
+    }
+
+    fn stack_maps(msa_len: usize, map_x: &Mapping, map_y: &Mapping) -> Mapping {
+        let mut map = Vec::with_capacity(msa_len);
+        let mut ind: usize = 0;
+        for (x, y) in map_x.iter().zip(map_y.iter()) {
+            if x.is_none() && y.is_none() {
+                map.push(None);
+            } else {
+                map.push(Some(ind));
+                ind += 1;
+            }
+        }
+        map
+    }
+
+    fn clear_common_gaps(msa_len: usize, map_x: &Mapping, map_y: &Mapping) -> PairwiseAlignment {
+        let mut upd_map_x = Vec::with_capacity(msa_len);
+        let mut upd_map_y = Vec::with_capacity(msa_len);
+        for (x, y) in map_x.iter().zip(map_y.iter()) {
+            if x.is_none() && y.is_none() {
+                continue;
+            } else {
+                upd_map_x.push(*x);
+                upd_map_y.push(*y);
+            }
+        }
+        PairwiseAlignment::new(upd_map_x, upd_map_y)
+    }
+}