Moved optimisation fully to Tree impl, runs but needs checking

bacpop · Feb 26, 2024 · c2a8bea · c2a8bea
1 parent be110b7
commit c2a8bea
Show file tree

Hide file tree

Showing 2 changed files with 5 additions and 69 deletions.
diff --git a/src/dspsa.rs b/src/dspsa.rs
@@ -102,12 +102,10 @@ impl Tree {
             // Set new theta
             theta = theta.iter().zip(ghat.iter())
             .map(|(theta, g)| *theta - ak * g).collect();
-        println!("New tree vector is: {:?}", theta);
         }
 
         // Update final tree after finishing optimisation
-        println!("New tree vector is: {:?}", theta);
-        let new_tree_vec: Vec<usize> = theta.iter().map(|x| *x as usize).collect();
+        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_tree(Some(new_tree_vec), false);
         self.update_likelihood(&q);

diff --git a/src/lib.rs b/src/lib.rs
@@ -51,72 +51,10 @@ pub fn main() {
     println!("{:?}", tr.newick());
     println!("{:?}", tr.tree_vec);
 
-    tr.optimise(&q, 5);
-
-    // if !args.no_optimise {
-    //     let mut theta: Vec<f64> = tr.tree_vec.iter().map(|x| *x as f64).collect();
-    //     let n = theta.len();
-
-    //     let a = 2.0;
-    //     let A = 2.0;
-    //     let alpha = 0.75;
-    //     // let k = 0;
-
-    //     let mut llvec: Vec<f64> = Vec::new();
-
-    //     let start = Instant::now();
-    //     for k in 0..=200 {
-    //         println!("k: {:?}", k);
-    //         // println!("theta: {:?}", theta);
-
-    //     //     // Peturbation vector
-    //         let delta = peturbation_vec(n);
-    //         // println!("delta: {:?}", delta);
-
-    //     //     // Pi vector
-    //         let pivec: Vec<f64> = piv(&theta);
-    //     //     // println!("pivec: {:?}", pivec);
-
-    //     //     // theta+/-
-    //         let thetaplus: Vec<usize> = pivec.iter().zip(delta.iter()).map(|(x, y)| (x + (y / 2.0)).round() as usize).collect();
-    //         let thetaminus: Vec<usize> = pivec.iter().zip(delta.iter()).map(|(x, y)| (x - (y / 2.0)).round() as usize).collect();
-
-    //     //     // println!("thetaplus: {:?}", thetaplus);
-    //     //     // println!("thetaminus: {:?}", thetaminus);
-
-    //     //     // Calculate likelihood at theta trees
-    //         tr.update_tree(Some(thetaplus), false);
-    //     //     // println!("tree changes: {:?}", tr.changes);
-    //         tr.update_likelihood(&q);
-    //         let x1 = tr.get_tree_likelihood();
-    //     //     // println!("thetaplus ll: {:?}", x1);
-
-    //         tr.update_tree(Some(thetaminus), false);
-    //     //     // println!("tree changes: {:?}", tr.changes);
-    //         tr.update_likelihood(&q);
-    //         let x2 = tr.get_tree_likelihood();
-    //     //     // println!("thetaminus ll: {:?}", x2);
-
-    //     //     // Calculations to work out new theta
-    //         let ldiff = x1 - x2; 
-    //         let ghat: Vec<f64> = delta.iter().map(|el| if !el.eq(&0.0) {el * ldiff} else {0.0}).collect();
-
-    //         let ak = a / (1.0 + A + k as f64).powf(alpha);
-
-    //         theta = theta.iter().zip(ghat.iter()).map(|(theta, g)| *theta - ak * g).collect();
-
-    //         llvec.push(x1);
-
-    //     //     // println!("ghat: {:?}", ghat);
-
-    //     }
-
-    //     let out: Vec<f64> = phi(&theta).iter().map(|x| x.round()).collect();
-    //     println!("final theta: {:?}", out);
-
-    //     println!("{:?}", &llvec);
-    //     // println!("{:?}", &llvec[95..100]);
-    // }
+    if !args.no_optimise {
+        tr.optimise(&q, 5);
+    }
+
     // let end = Instant::now();
 
     // eprintln!("Done in {}s", end.duration_since(start).as_secs());