From 29c1db8c319ca3dde9889f0fc691da8a1ebd37b1 Mon Sep 17 00:00:00 2001 From: Joel Hellewell Date: Thu, 21 Mar 2024 13:57:37 +0000 Subject: [PATCH] Added very simple hill climbing optimiser --- src/dspsa.rs | 70 ++++++++++++++++++++++++++++++++++++++++-------- src/lib.rs | 3 ++- src/phylo2vec.rs | 6 +---- 3 files changed, 62 insertions(+), 17 deletions(-) diff --git a/src/dspsa.rs b/src/dspsa.rs index 04336bb..86ba93e 100644 --- a/src/dspsa.rs +++ b/src/dspsa.rs @@ -51,12 +51,56 @@ pub fn theta_change(pivec: &Vec, delta: &Vec, plus: bool) -> Vec, iterations: usize) { + + let mut working_tree: Tree = Tree { + tree_vec: self.tree_vec.clone(), + nodes: self.nodes.clone(), + max_depth: self.max_depth, + leaf_permutations: self.leaf_permutations.clone(), + changes: self.changes.clone(), + mutation_lists: self.mutation_lists.clone() + }; + + let mut theta: Vec = self.tree_vec.iter().map(|x| *x as f64).collect(); + let n: usize = theta.len(); + let mut delta: Vec = Vec::with_capacity(n); + let mut pivec: Vec = Vec::with_capacity(n); + let mut thetaplus: Vec = Vec::with_capacity(n); + + for k in 0..=iterations { + println!("Optimisation step {} out of {}", k, iterations); + println!("Tree log likelihood: {}", self.get_tree_likelihood()); + // Generate peturbation vector + delta = peturbation_vec(n); + // println!("Peturbation vector: {:?}", delta); + + // Generate pi vector + pivec = piv(&theta); + // println!("Pi vector: {:?}", pivec); + + // Calculate theta+ and theta-, + // New tree vectors based on peturbation + thetaplus = theta_change(&pivec, &delta, true); + + working_tree.update_quad(thetaplus); + working_tree.update_likelihood(&q); + + if(working_tree.get_tree_likelihood() > self.get_tree_likelihood()) { + println!("Climbing hill!"); + self.update_quad(working_tree.tree_vec.clone()); + self.update_likelihood(&q); + theta = self.tree_vec.iter().map(|x| *x as f64).collect(); + } + } + } + pub fn optimise(&mut self, q: &na::Matrix4, iterations: usize) { // Update likelihood if not done already - if self.get_tree_likelihood().eq(&0.0) { - self.update_likelihood(&q); - } + // if self.get_tree_likelihood().eq(&0.0) { + // self.update_likelihood(&q); + // } // Convert tree vector to Vec let mut theta: Vec = self.tree_vec.iter().map(|x| *x as f64).collect(); @@ -67,7 +111,7 @@ impl Tree { // Tuning parameters for optimisation, will // eventually have defaults or be passed in let a: f64 = 1.5; - let cap_a: f64 = 10000.0; + let cap_a: f64 = 1000.0; let alpha: f64 = 0.51; // Pre-allocate vectors @@ -76,11 +120,12 @@ impl Tree { let mut thetaplus: Vec = Vec::with_capacity(n); let mut thetaminus: Vec = Vec::with_capacity(n); let mut ghat: Vec = Vec::with_capacity(n); + let mut new_tree_vec: Vec = Vec::with_capacity(n); // Optimisation loop for k in 0..=iterations { println!("Optimisation step {} out of {}", k, iterations); - println!("Tree likelihood: {}", -self.get_tree_likelihood()); + println!("Negative tree log likelihood: {}", -self.get_tree_likelihood()); // Generate peturbation vector delta = peturbation_vec(n); // println!("Peturbation vector: {:?}", delta); @@ -121,16 +166,19 @@ impl Tree { // Set new theta theta = theta.iter().zip(ghat.iter()) - .map(|(theta, g)| *theta - ak * g).collect(); + .map(|(theta, g)| *theta - (ak * g)).collect(); + new_tree_vec = phi(&theta).iter().map(|x| x.round() as usize).collect(); + self.update_quad(new_tree_vec); + self.update_likelihood(&q) // println!("New theta is: {:?}", theta); } // Update final tree after finishing optimisation - let new_tree_vec: Vec = phi(&theta).iter().map(|x| x.round() as usize).collect(); - println!("New tree vector is: {:?}", new_tree_vec); - self.update_quad(new_tree_vec); - self.update_likelihood(&q); - println!("New tree likelihood is {}", -self.get_tree_likelihood()); + // let new_tree_vec: Vec = phi(&theta).iter().map(|x| x.round() as usize).collect(); + // println!("New tree vector is: {:?}", new_tree_vec); + // self.update_quad(new_tree_vec); + // self.update_likelihood(&q); + println!("New tree likelihood is {}", self.get_tree_likelihood()); } } \ No newline at end of file diff --git a/src/lib.rs b/src/lib.rs index ec9a1f3..cd7b353 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -52,7 +52,8 @@ pub fn main() { println!("{:?}", tr.tree_vec); if !args.no_optimise { - tr.optimise(&q, 25); + // tr.optimise(&q, 10); + tr.hillclimb(&q, 20); } // let end = Instant::now(); diff --git a/src/phylo2vec.rs b/src/phylo2vec.rs index 7cb6881..32de3b4 100644 --- a/src/phylo2vec.rs +++ b/src/phylo2vec.rs @@ -109,7 +109,7 @@ pub fn random_tree(k: usize) -> Vec { vec![0; k] .iter() .enumerate() - .map(|(i, _el)| if i > 0 { rng.gen_range(0..i) } else { 0 }) + .map(|(i, _el)| if i > 0 { rng.gen_range(0..((2 * i) - 1)) } else { 0 }) .collect() } @@ -188,10 +188,6 @@ impl Tree { pub fn update_quad(&mut self, new_vec: Vec) { - // if !self.changes.is_empty() { - // panic!("There are already changes that need updating"); - // } - let new_tree: Tree = phylo2vec_quad(new_vec); let k: usize = new_tree.nodes.len(); let mut old_parent: Option;