Added very simple hill climbing optimiser

src/dspsa.rs

@@ -51,12 +51,56 @@ pub fn theta_change(pivec: &Vec<f64>, delta: &Vec<f64>, plus: bool) -> Vec<usize
 }
 
 impl Tree {
+    pub fn hillclimb(&mut self, q: &na::Matrix4<f64>, 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<f64> = self.tree_vec.iter().map(|x| *x as f64).collect();
+        let n: usize = theta.len();
+        let mut delta: Vec<f64> = Vec::with_capacity(n);
+        let mut pivec: Vec<f64> = Vec::with_capacity(n);
+        let mut thetaplus: Vec<usize> = 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<f64>, 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<f64>
         let mut theta: Vec<f64> = 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<usize> = Vec::with_capacity(n);
         let mut thetaminus: Vec<usize> = Vec::with_capacity(n);
         let mut ghat: Vec<f64> = Vec::with_capacity(n);
+        let mut new_tree_vec: Vec<usize> = 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<usize> = 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<usize> = 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());
     }
 }

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();

src/phylo2vec.rs

@@ -109,7 +109,7 @@ pub fn random_tree(k: usize) -> Vec<usize> {
     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<usize>) {
 
-        // 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<usize>;