Skip to content

Commit

Permalink
Fixing clippy suggestions
Browse files Browse the repository at this point in the history
  • Loading branch information
jhellewell14 committed May 22, 2024
1 parent 7350a30 commit 096ecea
Show file tree
Hide file tree
Showing 4 changed files with 69 additions and 71 deletions.
34 changes: 21 additions & 13 deletions src/lib.rs
Original file line number Diff line number Diff line change
@@ -1,14 +1,13 @@
mod mutation;
mod build_tree;
mod hillclimb;
mod likelihoods;
mod mutation;
mod node;
mod build_tree;
mod tests;
mod tree;
mod hillclimb;
mod tree_iterators;
mod tree_to_newick;

use crate::hillclimb::peturb_vector;
use crate::build_tree::*;
use crate::tree::Tree;
extern crate nalgebra as na;
Expand All @@ -17,16 +16,28 @@ use crate::cli::*;
use std::time::Instant;

pub fn main() {
let args = cli_args();
let args = cli_args();

// Define rate matrix
let q: na::Matrix4<f64> = na::Matrix4::new(
-1.0, 1.0 / 3.0, 1.0 / 3.0, 1.0 / 3.0,
1.0 / 3.0, -1.0, 1.0 / 3.0, 1.0 / 3.0,
1.0 / 3.0, 1.0 / 3.0, -1.0, 1.0 / 3.0,
1.0 / 3.0, 1.0 / 3.0, 1.0 / 3.0, -1.0,
-1.0,
1.0 / 3.0,
1.0 / 3.0,
1.0 / 3.0,
1.0 / 3.0,
-1.0,
1.0 / 3.0,
1.0 / 3.0,
1.0 / 3.0,
1.0 / 3.0,
-1.0,
1.0 / 3.0,
1.0 / 3.0,
1.0 / 3.0,
1.0 / 3.0,
-1.0,
);

let mut tr = vector_to_tree(&random_vector(27));
tr.add_genetic_data(&args.alignment);

Expand All @@ -45,7 +56,4 @@ pub fn main() {
// eprintln!("Done in {}ms", end.duration_since(start).as_millis());
// eprintln!("Done in {}ns", end.duration_since(start).as_nanos());
}



}
54 changes: 27 additions & 27 deletions src/likelihoods.rs
Original file line number Diff line number Diff line change
Expand Up @@ -5,22 +5,29 @@ use logaddexp::LogAddExp;
const BF_DEFAULT: [f64; 4] = [0.25, 0.25, 0.25, 0.25];

// Calculates the likelihood at a base given the bases at each child and probability matrices
pub fn base_likelihood(mut1: &Mutation, mut2: &Mutation, p1: &na::Matrix4<f64>, p2: &na::Matrix4<f64>) -> Mutation {

pub fn base_likelihood(
mut1: &Mutation,
mut2: &Mutation,
p1: &na::Matrix4<f64>,
p2: &na::Matrix4<f64>,
) -> Mutation {
child_likelihood(mut1, p1).add(child_likelihood(mut2, p2))

}

// Calculates the likelihood for one child at a base
pub fn child_likelihood(muta: &Mutation, p: &na::Matrix4<f64>) -> Mutation {

let mut outmut: Mutation = Mutation::default();

// This does the correct row multiplications between the existing likelihood values and the probability matrix
for i in 0..=3 {
if let Some(val) = outmut.get_mut(i) {
*val = p.row(i).iter().zip(muta.iter()).map(|(a, b)| a.ln() + b)
.reduce(|a, b| a.ln_add_exp(b)).unwrap();
*val = p
.row(i)
.iter()
.zip(muta.iter())
.map(|(a, b)| a.ln() + b)
.reduce(|a, b| a.ln_add_exp(b))
.unwrap();
}
}

Expand All @@ -34,49 +41,48 @@ pub fn calculate_likelihood(
p1: &na::Matrix4<f64>,
p2: &na::Matrix4<f64>,
) -> Vec<Mutation> {

// Iterate over bases and calculate likelihood
let out: Vec<Mutation> = seq1.iter()
.zip(seq2.iter())
.map(|(b1, b2)| base_likelihood(b1, b2, &p1, &p2))
.collect();
let out: Vec<Mutation> = seq1
.iter()
.zip(seq2.iter())
.map(|(b1, b2)| base_likelihood(b1, b2, p1, p2))
.collect();

out
}

// LogSumExp function that includes base frequency values for final likelihood calculation
pub fn base_freq_logse(muta: &Mutation, bf: [f64; 4]) -> f64 {
muta.iter().zip(bf.iter()).fold(0.0, |tot, (muta, bf)| tot + muta.exp() * bf).ln()
muta.iter()
.zip(bf.iter())
.fold(0.0, |tot, (muta, bf)| tot + muta.exp() * bf)
.ln()
}

impl Tree {

// Updates the genetic likelihood at a given node
pub fn update_node_likelihood(&mut self, index: usize, rate_matrix: &na::Matrix4<f64>) {
if let (Some(ch1), Some(ch2)) = self.get_node(index).unwrap().children {

let p1 = na::Matrix::exp(&(rate_matrix * self.get_branchlength(ch1)));
let p2 = na::Matrix::exp(&(rate_matrix * self.get_branchlength(ch2)));

let seq1 = self.mutation_lists.get(ch1).unwrap();
let seq2 = self.mutation_lists.get(ch2).unwrap();

self.mutation_lists[index] = calculate_likelihood(seq1, seq2, &p1, &p2);
}
}

// Goes through all nodes that have changed and updates genetic likelihood
// Used after tree.update()
pub fn update_likelihood(&mut self, rate_matrix: &na::Matrix4<f64>) {

if self.changes.is_empty() {
return
return;
}

let max_depth: usize = *self.changes.keys().max().unwrap();

for current_depth in (0..=max_depth).rev() {

let mut nodes: Vec<usize> = self.changes.remove(&current_depth).unwrap();
nodes.sort();
nodes.dedup();
Expand All @@ -88,7 +94,6 @@ impl Tree {

// Traverse all nodes at current_depth
for node in nodes {

self.update_node_likelihood(node, rate_matrix);

if current_depth > 0 {
Expand All @@ -106,7 +111,6 @@ impl Tree {
}
}
}

}

// Traverses tree in post-order below given node (except leaves), updating likelihood
Expand All @@ -120,7 +124,6 @@ impl Tree {
for node in nodes {
self.update_node_likelihood(node, rate_matrix);
}

}

// Fetches likelihood value for a tree
Expand All @@ -129,9 +132,6 @@ impl Tree {
.get(self.get_root().unwrap().index)
.unwrap()
.iter()
.fold(0.0, |acc, muta| {
acc + base_freq_logse(muta, BF_DEFAULT)
})
.fold(0.0, |acc, muta| acc + base_freq_logse(muta, BF_DEFAULT))
}

}
}
50 changes: 20 additions & 30 deletions src/mutation.rs
Original file line number Diff line number Diff line change
@@ -1,9 +1,8 @@

#[derive(Debug, Copy, Clone)]
pub struct Mutation(pub f64, pub f64, pub f64, pub f64);

const NEGINF: f64 = -f64::INFINITY;
// (A, C, G, T)
// (A, C, G, T)
const AMUT: Mutation = Mutation(0.0, NEGINF, NEGINF, NEGINF);
const CMUT: Mutation = Mutation(NEGINF, 0.0, NEGINF, NEGINF);
const GMUT: Mutation = Mutation(NEGINF, NEGINF, 0.0, NEGINF);
Expand Down Expand Up @@ -48,7 +47,8 @@ pub fn char_to_mutation(e: &char) -> Mutation {
'V' => VMUT,
'-' => {
// This way of coding gives same answer as other tree programs
Mutation::default()},
Mutation::default()
}
_ => panic!("Unrecognised character: {}", e),
}
}
Expand All @@ -63,49 +63,39 @@ pub fn char_to_mutation(e: &char) -> Mutation {
// }

impl Mutation {

// Converts a Mutation to a vector
// pub fn to_vector(self) -> Vec<f64> {
// vec![self.0, self.1, self.2, self.3]
// }

// Adds two Mutations together
pub fn add(self, r: Mutation) -> Mutation {
Mutation(
self.0 + r.0,
self.1 + r.1,
self.2 + r.2,
self.3 + r.3,
)
Mutation(self.0 + r.0, self.1 + r.1, self.2 + r.2, self.3 + r.3)
}

pub fn get(self, i: usize) -> Option<f64> {
match i {
0 => {Some(self.0)},
1 => {Some(self.1)},
2 => {Some(self.2)},
3 => {Some(self.3)},
_ => {None},
}
match i {
0 => Some(self.0),
1 => Some(self.1),
2 => Some(self.2),
3 => Some(self.3),
_ => None,
}
}

pub fn get_mut<'a>(&'a mut self, i: usize) -> Option<&mut f64> {
pub fn get_mut(&mut self, i: usize) -> Option<&mut f64> {
match i {
0 => {Some(&mut self.0)},
1 => {Some(&mut self.1)},
2 => {Some(&mut self.2)},
3 => {Some(&mut self.3)},
_ => {None},
0 => Some(&mut self.0),
1 => Some(&mut self.1),
2 => Some(&mut self.2),
3 => Some(&mut self.3),
_ => None,
}
}

pub fn iter(self) -> MutationIter {
MutationIter {
ind: 0,
muta: self
}
MutationIter { ind: 0, muta: self }
}

}

#[derive(Debug)]
Expand All @@ -114,7 +104,7 @@ pub struct MutationIter {
muta: Mutation,
}

impl<'a> Iterator for MutationIter {
impl Iterator for MutationIter {
type Item = f64;

fn next(&mut self) -> Option<Self::Item> {
Expand All @@ -128,4 +118,4 @@ impl Default for Mutation {
fn default() -> Self {
Mutation(0.0, 0.0, 0.0, 0.0)
}
}
}
2 changes: 1 addition & 1 deletion src/phylo2vec.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -509,4 +509,4 @@ std::unique_ptr<std::vector<int>> doToVector(std::string& newick, int num_leaves
}

return std::make_unique<std::vector<int>>(converted_v);
}
}

0 comments on commit 096ecea

Please sign in to comment.