From fe5575756cbfe72b1cb81943abcee788e44f596b Mon Sep 17 00:00:00 2001 From: jhellewell14 Date: Thu, 31 Oct 2024 14:05:35 +0000 Subject: [PATCH] Added peturb vector MoveFn, put moves in their own file --- src/genetic_data.rs | 79 ++------------------------------------------- src/hillclimb.rs | 68 -------------------------------------- src/lib.rs | 2 ++ src/moves.rs | 76 +++++++++++++++++++++++++++++++++++++++++++ 4 files changed, 81 insertions(+), 144 deletions(-) delete mode 100644 src/hillclimb.rs create mode 100644 src/moves.rs diff --git a/src/genetic_data.rs b/src/genetic_data.rs index 5c7db49..e7d4428 100644 --- a/src/genetic_data.rs +++ b/src/genetic_data.rs @@ -4,7 +4,7 @@ use needletail::parse_fastx_file; use crate::topology::Topology; use ndarray::s; use logaddexp::LogAddExp; -use rand::Rng; +use crate::moves::MoveFn; const NEGINF: f64 = -f64::INFINITY; // (A, C, G, T) @@ -153,63 +153,9 @@ pub fn base_freq_logse(muta: ndarray::ArrayBase, ndarray .fold(0.0, |tot, (muta, bf)| tot + muta.exp() * bf) .ln() } - -// A move produces a vector of nodes that have a new parent node -// Now Likelihood update -// Iterate over nodes that need changing from Topology -// Calculate new likelihoods at internal nodes, add to HashMap -// Get whole Topology likelihood at root node from HashMap -// If better, can move HashMap likelihoods into GeneticData - - -// A TryTopology struct that contains old Topology, new tree_vec, changes, reference to gen_data?, HashMap? - pub struct CandidateTopology{ - new_topology: Topology, - changes: Option>, -} - -pub fn peturb_vector(topology: &Topology) -> CandidateTopology { - - let mut vout = topology.tree_vec.to_vec(); - let n = vout.len().div_ceil(5); - let mut rng = rand::thread_rng(); - let ind_rng = rand::thread_rng(); - let distr = rand::distributions::Bernoulli::new(0.5).unwrap(); - let ind_distr = rand::distributions::Uniform::new(0, vout.len()); - - let mut inds: Vec = ind_rng.sample_iter(ind_distr).take(n).collect(); - inds.sort(); - - for ind in inds { - if ind.eq(&0) { - continue; - } - - match rng.sample(distr) { - true => { - if vout[ind].lt(&(2 * (ind - 1))) { - vout[ind] += 1; - } - } - false => { - if vout[ind].gt(&0) { - vout[ind] -= 1; - } - } - }; - } - - let new_topology: Topology = Topology::from_vec(&vout); - let changes: Option> = topology.find_changes(&new_topology); - CandidateTopology{ - new_topology, - changes, - } -} - -pub fn hillclimb_accept(old_ll: &f64, new_ll: &f64) -> bool { - new_ll.gt(old_ll) + pub new_topology: Topology, + pub changes: Option>, } impl Topology { @@ -306,22 +252,3 @@ impl Topology { } } -pub trait MoveFn { - fn generate_move(&self, current_topology: &Topology) -> CandidateTopology; -} - -pub struct ExactMove { - pub target_vector: Vec, -} - -impl MoveFn for ExactMove { - fn generate_move(&self, current_topology: &Topology) -> CandidateTopology { - - let new_topology: Topology = Topology::from_vec(&self.target_vector); - let changes: Option> = current_topology.find_changes(&new_topology); - CandidateTopology{ - new_topology, - changes, - } - } -} \ No newline at end of file diff --git a/src/hillclimb.rs b/src/hillclimb.rs deleted file mode 100644 index cc52f58..0000000 --- a/src/hillclimb.rs +++ /dev/null @@ -1,68 +0,0 @@ -// use crate::Tree; -// use rand::Rng; -// use crate::rate_matrix::RateMatrix; - -// // Makes random moves +/- 1 moves on the integer vector (v) for a given number of elements (n) -// pub fn peturb_vector(v: &[usize], n: usize) -> Vec { -// let mut vout = v.to_vec(); -// let mut rng = rand::thread_rng(); -// let ind_rng = rand::thread_rng(); -// let distr = rand::distributions::Bernoulli::new(0.5).unwrap(); -// let ind_distr = rand::distributions::Uniform::new(0, v.len()); - -// let mut inds: Vec = ind_rng.sample_iter(ind_distr).take(n).collect(); -// inds.sort(); - -// for ind in inds { -// if ind.eq(&0) { -// continue; -// } - -// match rng.sample(distr) { -// true => { -// if vout[ind].lt(&(2 * (ind - 1))) { -// vout[ind] += 1; -// } -// } -// false => { -// if vout[ind].gt(&0) { -// vout[ind] -= 1; -// } -// } -// }; -// } - -// vout -// } - -// impl Tree { -// // Hill climbing optimisation algorithm -// pub fn hillclimb(&mut self, iterations: usize) { -// let mut candidate_vec: Vec = Vec::with_capacity(self.tree_vec.len()); -// let mut best_vec: Vec = self.tree_vec.clone(); -// let mut best_likelihood: f64 = self.get_tree_likelihood(); -// let mut new_likelihood: f64; - -// for k in 0..=iterations { -// println!("Optimisation step {} out of {}", k, iterations); -// candidate_vec = peturb_vector(&best_vec, self.tree_vec.len()); -// println!("new vec: {:?}", candidate_vec); -// self.update(&candidate_vec); -// self.update_likelihood(); -// new_likelihood = self.get_tree_likelihood(); -// println!( -// "Candidate likelihood: {} \n Current likelihood: {}", -// new_likelihood, best_likelihood -// ); - -// if new_likelihood > best_likelihood { -// println!("Climbing hill!"); -// best_vec = candidate_vec; -// best_likelihood = new_likelihood; -// } -// } - -// self.update(&best_vec); -// self.update_likelihood(); -// } -// } diff --git a/src/lib.rs b/src/lib.rs index 7ed2422..b4995ff 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -4,6 +4,7 @@ mod iterators; mod rate_matrix; mod topology; mod genetic_data; +mod moves; use rate_matrix::RateMatrix; use topology::Topology; @@ -15,6 +16,7 @@ use crate::cli::*; use std::env::args; use std::time::Instant; use crate::genetic_data::*; +use crate::moves::*; pub fn main() { let args = cli_args(); diff --git a/src/moves.rs b/src/moves.rs new file mode 100644 index 0000000..500e460 --- /dev/null +++ b/src/moves.rs @@ -0,0 +1,76 @@ +use crate::Topology; +use crate::CandidateTopology; +use rand::Rng; + +pub trait MoveFn { + fn generate_move(&self, current_topology: &Topology) -> CandidateTopology; +} + +pub struct ExactMove { + pub target_vector: Vec, +} + +impl MoveFn for ExactMove { + fn generate_move(&self, current_topology: &Topology) -> CandidateTopology { + + let new_topology: Topology = Topology::from_vec(&self.target_vector); + let changes: Option> = current_topology.find_changes(&new_topology); + CandidateTopology{ + new_topology, + changes, + } + } +} + +pub struct PeturbVec { + pub n: usize, +} + +impl MoveFn for PeturbVec { + fn generate_move(&self, current_topology: &Topology) -> CandidateTopology { + let mut vout = current_topology.tree_vec.to_vec(); + let mut rng = rand::thread_rng(); + let ind_rng = rand::thread_rng(); + let distr = rand::distributions::Bernoulli::new(0.5).unwrap(); + let ind_distr = rand::distributions::Uniform::new(0, vout.len()); + + let samp_n: usize = match self.n.gt(&vout.len()) { + true => {vout.len()}, + false => {self.n}, + }; + + let mut inds: Vec = ind_rng.sample_iter(ind_distr).take(samp_n).collect(); + inds.sort(); + + for ind in inds { + if ind.eq(&0) { + continue; + } + + match rng.sample(distr) { + true => { + if vout[ind].lt(&(2 * (ind - 1))) { + vout[ind] += 1; + } + } + false => { + if vout[ind].gt(&0) { + vout[ind] -= 1; + } + } + }; + }; + + let new_topology: Topology = Topology::from_vec(&vout); + let changes: Option> = current_topology.find_changes(&new_topology); + CandidateTopology{ + new_topology, + changes, + } + } +} + + +pub fn hillclimb_accept(old_ll: &f64, new_ll: &f64) -> bool { + new_ll.gt(old_ll) +}