diff --git a/src/lib.rs b/src/lib.rs index 16b30e1..4d630f7 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -40,6 +40,37 @@ pub fn main() { -1.0, ); + fn create_GTR_ratematrix(a: f64, b: f64, c: f64, d: f64, e: f64, f: f64, pv: Vec) -> na::Matrix4 { + // pv = pivec defined as (piA, piC, piG, piT) + let mut q = na::Matrix4::new( + -(a * pv[1] + b * pv[2] + c * pv[3]), + a * pv[1], + b * pv[2], + c * pv[3], + a * pv[0], + -(a * pv[0] + d * pv[2] + e * pv[3]), + d * pv[2], + e * pv[3], + b * pv[0], + d * pv[1], + -(b * pv[0] + d * pv[1] + f * pv[3]), + f * pv[3], + c * pv[0], + e * pv[1], + f * pv[2], + -(c * pv[0] + e * pv[1] * f * pv[2])); + + let mut diag = 0.0; + for i in 0..=3 { + diag -= q[(i, i)] * pv[i]; + } + + q / diag + } + + let q2 = create_GTR_ratematrix(1.0, 1.0, 1.0, 1.0, 1.0, 1.0, vec![0.25, 0.25, 0.25, 0.25]); + // println!("{:?}", q2); + // let mut tr = vector_to_tree(&random_vector(4)); // tr.add_genetic_data(&String::from("/Users/joel/Downloads/listeria0.aln")); let mut tr = vector_to_tree(&random_vector(27)); @@ -52,10 +83,10 @@ pub fn main() { if !args.no_optimise { let start = Instant::now(); - tr.hillclimb(&q, 50); + tr.hillclimb(&q, 0); let end = Instant::now(); eprintln!("Done in {}s", end.duration_since(start).as_secs()); - eprintln!("Done in {}s", end.duration_since(start).as_millis()); + eprintln!("Done in {}ms", end.duration_since(start).as_millis()); } }