Renaming and tidying gen_list into mutation.rs

src/graveyard.rs

@@ -147,3 +147,92 @@
 // pub fn get_tips(&self) -> Vec<&Node> {
 //     self.nodes.iter().filter(|n| n.tip).collect()
 // }
+
+
+// pub fn create_dummy_genetic_data(n_leaves: usize, n_mutations: usize, sequence_length: usize) -> Vec<Vec<Mutation>> {
+//     let mut output: Vec<Vec<Mutation>> = Vec::new();
+//     let mut rng = rand::thread_rng();
+
+//     for i in 0..n_leaves {
+//         let mut temp: Vec<Mutation> = Vec::new();
+//         for j in 0..n_mutations {
+//             let mut mutation = Mutation(rng.gen_range(1..sequence_length), 0.0, 0.0, 0.0, 0.0);
+//             match rng.gen_range(1..=4) {
+//                 1 => {mutation.1 = 1.0},
+//                 2 => {mutation.2 = 1.0},
+//                 3 => {mutation.3 = 1.0},
+//                 4 => {mutation.4 = 1.0},
+//                 _ => {},
+//             }
+//             temp.push(mutation);
+//         }
+//         temp.sort_by(|a, b| a.0.cmp(&b.0));
+//         temp.dedup_by(|a, b| a.0.eq(&b.0));
+//         output.push(temp);
+//     }
+
+//     for _ in 0..(n_leaves + 1) {
+//         output.push(Vec::new());
+//     }
+
+//     output
+// }
+
+// Combines two vectors of Mutations into a single vector
+// pub fn combine_lists(
+//     seq1: Option<&Vec<Mutation>>,
+//     seq2: Option<&Vec<Mutation>>,
+//     branchlengths: (f64, f64),
+//     rate_matrix: &na::Matrix4<f64>,
+// ) -> Vec<Mutation> {
+//     let mut out: Vec<Mutation> = Vec::new();
+
+//     // Probability matrices
+//     let p1 = na::Matrix::exp(&(rate_matrix * branchlengths.0));
+//     let p2 = na::Matrix::exp(&(rate_matrix * branchlengths.1));
+
+//     let mut s1 = seq1.unwrap().iter();
+//     let mut s2 = seq2.unwrap().iter();
+
+//     let mut mut1 = s1.next();
+//     let mut mut2 = s2.next();
+
+//     while mut1.is_some() | mut2.is_some() {
+//         if mut1.is_none() {
+//             // First iterator empty, push second
+//             out.push(mut2.unwrap().child_log_likelihood(&p2));
+//             mut2 = s2.next();
+//         } else if mut2.is_none() {
+//             // Second iterator empty, push first
+//             out.push(mut1.unwrap().child_log_likelihood(&p1));
+//             mut1 = s1.next();
+//         } else {
+//             // println!("mut1 = {:?} mut2 = {:?}", mut1.unwrap(), mut2.unwrap());
+//             // Neither iterator empty, compare indices of mutations and push highest
+//             // or combine likelihood if mutations at same location
+//             match mut1.unwrap().0.cmp(&mut2.unwrap().0) {
+//                 Ordering::Equal => {
+//                     // println!("mut1 == mut2 so pushing {:?}", mut1.unwrap());
+//                     out.push(
+//                         mut1.unwrap()
+//                             .child_log_likelihood(&p1)
+//                             .sum(mut2.unwrap().child_log_likelihood(&p2)),
+//                     );
+//                     mut1 = s1.next();
+//                     mut2 = s2.next();
+//                 }
+//                 Ordering::Greater => {
+//                     // println!("mut1 > mut2 so pushing {:?}", mut2.unwrap());
+//                     out.push(mut2.unwrap().child_log_likelihood(&p2));
+//                     mut2 = s2.next();
+//                 }
+//                 Ordering::Less => {
+//                     // println!("mut2 > mut1 so pushing {:?}", mut1.unwrap());
+//                     out.push(mut1.unwrap().child_log_likelihood(&p1));
+//                     mut1 = s1.next();
+//                 }
+//             }
+//         }
+//     }
+//     out
+// }

src/hillclimb.rs

@@ -45,15 +45,15 @@ impl Tree {
 
             println!("Optimisation step {} out of {}", k, iterations);
             // println!("Old vector {:?}", self.tree_vec);
-            println!("Tree log likelihood: {}", self.get_tree_likelihood());
+            println!("Current best likelihood: {}", best_likelihood);
 
             candidate_vec = hill_peturb(&self.tree_vec, self.tree_vec.len());
             working_tree.update(&candidate_vec);
             working_tree.update_likelihood(q);
             let new_likelihood = working_tree.get_tree_likelihood();
 
             // println!("New vector {:?}", candidate_vec);
-            println!("New likelihood {}", new_likelihood);
+            println!("Candidate likelihood {}", new_likelihood);
 
             if new_likelihood > best_likelihood {
                 println!("Climbing hill!");

src/lib.rs

@@ -1,4 +1,4 @@
-mod gen_list;
+mod mutation;
 mod import;
 mod likelihoods;
 mod node;
@@ -9,12 +9,9 @@ mod hillclimb;
 mod tree_iterators;
 mod tree_to_newick;
 
-use crate::gen_list::*;
 use crate::build_tree::*;
 use crate::tree::Tree;
 extern crate nalgebra as na;
-
-
 pub mod cli;
 use crate::cli::*;
 
@@ -46,7 +43,7 @@ pub fn main() {
     println!("{:?}", tr.tree_vec);
 
     if !args.no_optimise {
-        // tr.hillclimb(&q, 100);
+        tr.hillclimb(&q, 25);
     }
 
     // let end = Instant::now();