Moving back to optimisation algorithm

bacpop · Feb 23, 2024 · 6378788 · 6378788
1 parent a7b38be
commit 6378788
Showing 1 changed file with 36 additions and 51 deletions.
diff --git a/src/main.rs b/src/main.rs
@@ -48,88 +48,73 @@ fn main() {
     println!("{:?}", tr.tree_vec);
     // println!("{:?}", tr.nodes);
 
-    // let mut out = it.newick;
+    let mut theta: Vec<f64> = tr.tree_vec.iter().map(|x| *x as f64).collect();
+    let n = theta.len();
 
-    // while it.next().is_some() {
-    //     it.next();
-
-    // }
-
-    // println!("{:?}", &it.newick);
-    // let prr: String = pr.chars().rev().collect();
-    // println!("{:?}", prr);
+    let a = 2.0;
+    let A = 2.0;
+    let alpha = 0.75;
+    // let k = 0;
 
-    // out.reverse();
-    // println!("{:?}", out.join(""));
-
-
-    // let mut theta: Vec<f64> = tr.tree_vec.iter().map(|x| *x as f64).collect();
-    // let n = theta.len();
+    let mut llvec: Vec<f64> = Vec::new();
 
-    // 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..=100 {
-    //     println!("k: {:?}", k);
-    //     // println!("theta: {:?}", theta);
+    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);
+        let delta = peturbation_vec(n);
+        // println!("delta: {:?}", delta);
 
     //     // Pi vector
-    //     let pivec: Vec<f64> = piv(&theta);
+        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();
+        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);
+        tr.update_tree(Some(thetaplus), false);
     //     // println!("tree changes: {:?}", tr.changes);
-    //     tr.update_likelihood(&q);
-    //     let x1 = tr.get_tree_likelihood();
+        tr.update_likelihood(&q);
+        let x1 = tr.get_tree_likelihood();
     //     // println!("thetaplus ll: {:?}", x1);
 
-    //     tr.update_tree(Some(thetaminus), false);
+        tr.update_tree(Some(thetaminus), false);
     //     // println!("tree changes: {:?}", tr.changes);
-    //     tr.update_likelihood(&q);
-    //     let x2 = tr.get_tree_likelihood();
+        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 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);
+        let ak = a / (1.0 + A + k as f64).powf(alpha);
 
-    //     theta = theta.iter().zip(ghat.iter()).map(|(theta, g)| *theta as f64 - ak * g).collect();
+        theta = theta.iter().zip(ghat.iter()).map(|(theta, g)| *theta as f64 - ak * g).collect();
 
-    //     llvec.push(x1);
+        llvec.push(x1);
 
     //     // println!("ghat: {:?}", ghat);
 
-    // }
+    }
 
-    // let out: Vec<f64> = phi(&theta).iter().map(|x| x.round()).collect();
-    // println!("final theta: {:?}", out);
-    // let end = Instant::now();
+    let out: Vec<f64> = phi(&theta).iter().map(|x| x.round()).collect();
+    println!("final theta: {:?}", out);
+    let end = Instant::now();
 
 
-    // // // println!("{:?}", &llvec[0..5]);
-    // // // println!("{:?}", &llvec[95..100]);
+    println!("{:?}", &llvec);
+    // println!("{:?}", &llvec[95..100]);
 
-    // eprintln!("Done in {}s", end.duration_since(start).as_secs());
-    // eprintln!("Done in {}ms", end.duration_since(start).as_millis());
-    // eprintln!("Done in {}ns", end.duration_since(start).as_nanos());
+    eprintln!("Done in {}s", end.duration_since(start).as_secs());
+    eprintln!("Done in {}ms", end.duration_since(start).as_millis());
+    eprintln!("Done in {}ns", end.duration_since(start).as_nanos());
 
 }