Skip to content

Commit

Permalink
Fixed likelihood tests for new data format. Added a function that cre…
Browse files Browse the repository at this point in the history 
…ates random dummy genetic data
  • Loading branch information
@jhellewell14
jhellewell14 committed Nov 12, 2024
1 parent cf08062 commit 09d046f
Show file tree
Hide file tree
Showing 4 changed files with 67 additions and 89 deletions.
21 changes: 21 additions & 0 deletions src/genetic_data.rs
View file
Original file line number Diff line number Diff line change
@@ -1,10 +1,12 @@
use std::os::unix::thread;
use std::thread::current;
use std::collections::HashMap;
use needletail::parse_fastx_file;
use crate::topology::Topology;
use ndarray::s;
use logaddexp::LogAddExp;
use crate::moves::MoveFn;
use rand::{thread_rng, Rng};

const NEGINF: f64 = -f64::INFINITY;
// (A, C, G, T)
Expand Down Expand Up @@ -103,6 +105,25 @@ pub fn create_internal_data(mut gen_data: ndarray::ArrayBase<ndarray::OwnedRepr<

gen_data
}

pub fn create_dummy_gendata(n_bases: usize, topology: &Topology, rate_matrix: &na::Matrix4<f64>) -> ndarray::ArrayBase<ndarray::OwnedRepr<f64>, ndarray::Dim<[usize; 3]>> {

let n_seqs = topology.count_leaves();

let mut gen_data: ndarray::ArrayBase<ndarray::OwnedRepr<f64>, ndarray::Dim<[usize; 3]>> =
ndarray::Array3::from_elem(((2 * n_seqs) + 1, n_bases, 4), 0.0);

let mut rng = thread_rng();

for i in 0..n_seqs {
for j in 0..n_bases {
let k = rng.gen_range(0..4);
gen_data[[i, j, k]] = NEGINF;
}
}

create_internal_data(gen_data, topology, rate_matrix)
}

pub fn child_likelihood_i(i: usize, ll: ndarray::ArrayBase<ndarray::ViewRepr<&f64>, ndarray::Dim<[usize; 1]>>, p: &na::Matrix4<f64>) -> f64 {
p
Expand Down
16 changes: 1 addition & 15 deletions src/lib.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -34,27 +34,13 @@ pub fn main() {
let p = &rate_matrix::GTR::default();
let mut gen_data = create_genetic_data(&args.alignment, &t, &p.get_matrix());


// println!("{:?}", slice_gen_data(1, &gen_data));

// tr.initialise_likelihood();
// println!("tree likelihood {}", tr.get_tree_likelihood());
println!("{:?}", likelihood(&t, &gen_data));
println!("{:?}", t.get_newick());
println!("{:?}", t.tree_vec);

// let new_v = random_vector(27);
// let mut t2: Topology = Topology::from_vec(&new_v);

// let mv = ExactMove{target_vector: new_v};
// t.apply_move(mv, hillclimb_accept, &mut gen_data, &p.get_matrix());
// println!("updated ll: {:?}, initialised ll: {:?}", likelihood(&t, &gen_data), likelihood(&t2, &gen_data));



if !args.no_optimise {
let start = Instant::now();
for i in 0..10 {
for i in 0..1 {
println!{"Step {}", i};
// let new_v = random_vector(27);
// let mv = ExactMove{target_vector: new_v};
Expand Down
4 changes: 4 additions & 0 deletions src/moves.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -143,3 +143,7 @@ impl MoveFn for NearestNeighbour {
pub fn hillclimb_accept(old_ll: &f64, new_ll: &f64) -> bool {
new_ll.gt(old_ll)
}

pub fn always_accept(old_ll: &f64, new_ll: &f64) -> bool {
true
}
115 changes: 41 additions & 74 deletions src/tests.rs
View file
Original file line number Diff line number Diff line change
@@ -1,8 +1,13 @@
#[cfg(test)]
mod tests {
use crate::newick_to_vec::*;
use crate::rate_matrix::GTR;
use predicates::constant::always;

use crate::{likelihood, newick_to_vec::*};
use crate::rate_matrix::{RateMatrix, GTR};
use crate::topology::Topology;
use crate::moves::ExactMove;
use crate::create_dummy_gendata;
use crate::always_accept;

#[test]
fn check_topology_build_manual() {
Expand Down Expand Up @@ -79,84 +84,46 @@ mod tests {
assert_eq!(top.nodes[6].get_parent(), None);
}

// #[test]
// fn update_tree() {
// let mut tree_1 = vector_to_tree(&vec![0, 0, 1, 0], &GTR::default());

// let vecs: Vec<Vec<usize>> = vec![vec![0, 0, 0, 0], vec![0, 0, 1, 0], vec![0, 0, 1, 2], vec![0, 0, 1, 1]];

// for vec in vecs {
// let tree_2 = vector_to_tree(&vec, &GTR::default());
// tree_1.update(&vec);

// for i in 0..=tree_1.tree_vec.len() {
// assert_eq!(
// tree_1.get_node(i].get_parent(),
// tree_2.get_node(i].get_parent());
// assert_eq!(
// tree_1.get_node(i).unwrap().index,
// tree_2.get_node(i).unwrap().index
// );
// }
// }
#[test]
fn update_tree() {
let p = GTR::default();
let mut t_1 = Topology::from_vec(&vec![0, 0, 1, 0]);

let mut gen_data = create_dummy_gendata(2, &t_1, &p.get_matrix());

let vecs: Vec<Vec<usize>> = vec![vec![0, 0, 0, 0], vec![0, 0, 1, 0], vec![0, 0, 1, 2], vec![0, 0, 1, 1]];

for vec in vecs {
let t_2 = Topology::from_vec(&vec);
let mv = ExactMove{target_vector: vec};
t_1.apply_move(mv, always_accept, &mut gen_data, &p.get_matrix());

for i in 0..t_1.nodes.len() {
assert_eq!(t_1.nodes[i].get_parent(), t_2.nodes[i].get_parent());
assert_eq!(t_1.nodes[i].get_id(), t_2.nodes[i].get_id());
};
}

// }

// #[test]
// fn likelihood_internal_consistency_check() {
// // let q: na::Matrix4<f64> = na::Matrix4::new(
// // -3.0, 1.0, 1.0, 1.0, 1.0, -3.0, 1.0, 1.0, 1.0, 1.0, -3.0, 1.0, 1.0, 1.0, 1.0, -3.0,
// // );

// let mut tr = vector_to_tree(&vec![0, 0, 0, 0], &GTR::default());

// let genetic_data = vec![
// vec![
// Mutation(1.0, 0.0, 0.0, 0.0),
// Mutation(1.0, 0.0, 0.0, 0.0),
// ],
// vec![
// Mutation(0.0, 1.0, 0.0, 0.0),
// Mutation(1.0, 0.0, 0.0, 0.0),
// ],
// vec![
// Mutation(0.0, 0.0, 1.0, 0.0),
// Mutation(1.0, 0.0, 0.0, 0.0),
// ],
// vec![
// Mutation(1.0, 0.0, 0.0, 0.0),
// Mutation(0.0, 0.0, 0.0, 1.0),
// ],
// vec![
// Mutation(0.0, 1.0, 0.0, 0.0),
// Mutation(0.0, 0.0, 0.0, 1.0),
// ],
// vec![
// Mutation(0.0, 0.0, 1.0, 0.0),
// Mutation(0.0, 0.0, 0.0, 1.0),
// ],
// vec![
// Mutation(0.0, 1.0, 0.0, 0.0),
// Mutation(1.0, 0.0, 0.0, 0.0),
// ],
// ];

// tr.mutation_lists = genetic_data;

// tr.initialise_likelihood();
}

// let old_likelihood = tr.get_tree_likelihood();
#[test]
fn likelihood_internal_consistency_check() {
let p = GTR::default();
let mut t = Topology::from_vec(&vec![0, 0, 0, 0]);
let mut gen_data = create_dummy_gendata(5, &t, &p.get_matrix());

// tr.update(&vec![0, 0, 0, 1]);
// tr.update_likelihood();
let old_likelihood = likelihood(&t, &gen_data);

// tr.update(&vec![0, 0, 0, 0]);
// tr.update_likelihood();
let mv = ExactMove{target_vector: vec![0, 0, 0, 1]};
t.apply_move(mv, always_accept, &mut gen_data, &p.get_matrix());

// let new_likelihood = tr.get_tree_likelihood();
let mv = ExactMove{target_vector: vec![0, 0, 0, 0]};
t.apply_move(mv, always_accept, &mut gen_data, &p.get_matrix());

// assert_eq!(old_likelihood, new_likelihood);
// }
let new_likelihood = likelihood(&t, &gen_data);

assert_eq!(old_likelihood, new_likelihood);
}

#[test]
fn manual_parent_check () {
Expand Down

0 comments on commit 09d046f

Please sign in to comment.