rework on adjust bl for simulating trees, fix #25

[shajoezhu]

rm(list=ls())

dac <- function(n = 1000, m = 1000, x = 1000, y = 1000, z = 1000, w = 2000, l = 1:20000){
    return(1/(2*m) * exp(-(l-y-z-w)/(2*m)))
}

dab <- function(n = 1000, m = 1000, x = 1000, y = 1000, z = 1000, w = 2000, l = 1:20000){
    ret = rep(0, length(l))
    idx1 = which(l > (x+y) & l < (x+y+2*z))
    ret[idx1] = 1/(2*n) * exp(-(l[idx1]-x-y)/(2*n))
    idx2 = which(l > (x+y+2*z))
    ret[idx2] = exp(-z/n) * 1/(2*m) * exp(-(l[idx2]-y-z-2*z)/(2*m))
    return(ret)
}

library(ape)

# case1
#./hybrid-Lambda -spng "((t1:1000,t2:1000)int:1000,t3:2000);" -pop "((t1:1000,t2:1000)int:1000,t3:1000)r:1000;" -sim_num_gener -num 10000 -o sim -f
#joezhu@joezhu-M3800:~/hybrid-Lambda$ cat sim_frequencies
#1 (t2_1,(t3_1,t1_1));  1212
#2 (t1_1,(t2_1,t3_1));  1251
#3 ((t1_1,t2_1),t3_1);  7537
#joezhu@joezhu-M3800:~/hybrid-Lambda$ hybrid-Lambda -gt sim_num_gener -f
#Frequency file is saved at: OUT_frequencies
#joezhu@joezhu-M3800:~/hybrid-Lambda$ cat OUT_frequencies
#1 (t2_1,(t3_1,t1_1));  1212
#2 (t1_1,(t2_1,t3_1));  1251
#3 ((t1_1,t2_1),t3_1);  7537

gt <- read.tree("sim_num_gener")
png("case1.png", width = 1500, height =500)
par(mfrow=c(1,3))
hist(unlist(lapply(gt, function(x){sum(x$edge.length[which.edge(x,c("t1_1", "t2_1"))])})), probability=T)
lines(dab(),col="red")

hist(unlist(lapply(gt, function(x){sum(x$edge.length[which.edge(x,c("t1_1", "t3_1"))])})), probability=T)
lines(dac(),col="red")

hist(unlist(lapply(gt, function(x){sum(x$edge.length[which.edge(x,c("t2_1", "t3_1"))])})), probability=T)
lines(dac(),col="red")

dev.off()

#case2
#./hybrid-Lambda -spng "((t1:1000,t2:1000)int:1000,t3:2000);" -pop "((t1:1000,t2:1000)int:1000,t3:1000)r:2000;" -sim_num_gener -num 10000 -o sim -f
#joezhu@joezhu-M3800:~/hybrid-Lambda$ cat sim_frequencies
#1 ((t2_1,t1_1),t3_1);  7557
#2 (t2_1,(t1_1,t3_1));  1242
#3 ((t3_1,t2_1),t1_1);  1201
#joezhu@joezhu-M3800:~/hybrid-Lambda$ hybrid-Lambda -gt sim_num_gener -f
#Frequency file is saved at: OUT_frequencies
#joezhu@joezhu-M3800:~/hybrid-Lambda$ cat OUT_frequencies
#1 ((t2_1,t1_1),t3_1);  7557
#2 (t2_1,(t1_1,t3_1));  1242
#3 ((t3_1,t2_1),t1_1);  1201

gt <- read.tree("sim_num_gener")
png("case2.png", width = 1500, height =500)
par(mfrow=c(1,3))
hist(unlist(lapply(gt, function(x){sum(x$edge.length[which.edge(x,c("t1_1", "t2_1"))])})), probability=T)
lines(dab(m=2000),col="red")

hist(unlist(lapply(gt, function(x){sum(x$edge.length[which.edge(x,c("t1_1", "t3_1"))])})), probability=T)
lines(dac(m=2000),col="red")

hist(unlist(lapply(gt, function(x){sum(x$edge.length[which.edge(x,c("t2_1", "t3_1"))])})), probability=T)
lines(dac(m=2000),col="red")

dev.off()

#case3
#./hybrid-Lambda -spng "((t1:1000,t2:1000)int:1000,t3:2000);" -pop "((t1:1000,t2:1000)int:2000,t3:1000)r:1000;" -sim_num_gener -num 10000 -o sim -f
#hybrid-Lambda -gt sim_num_gener -f; cat sim_frequencies ; cat OUT_frequencies
#Frequency file is saved at: OUT_frequencies
#1 (t3_1,(t2_1,t1_1));  5981
#2 ((t3_1,t2_1),t1_1);  1990
#3 ((t1_1,t3_1),t2_1);  2029
#1 (t3_1,(t2_1,t1_1));  5981
#2 ((t3_1,t2_1),t1_1);  1990
#3 ((t1_1,t3_1),t2_1);  2029

gt <- read.tree("sim_num_gener")
png("case3.png", width = 1500, height =500)
par(mfrow=c(1,3))
hist(unlist(lapply(gt, function(x){sum(x$edge.length[which.edge(x,c("t1_1", "t2_1"))])})), probability=T)
lines(dab(n=2000),col="red")

hist(unlist(lapply(gt, function(x){sum(x$edge.length[which.edge(x,c("t1_1", "t3_1"))])})), probability=T)
lines(dac(n=2000),col="red")

hist(unlist(lapply(gt, function(x){sum(x$edge.length[which.edge(x,c("t2_1", "t3_1"))])})), probability=T)
lines(dac(n=2000),col="red")

dev.off()

#case4
#./hybrid-Lambda -spng "((t1:1000,t2:1000)int:1000,t3:2000);" -pop "((t1:1000,t2:1000)int:2000,t3:1000)r:2000;" -sim_num_gener -num 10000 -o sim -f
#hybrid-Lambda -gt sim_num_gener -f; cat sim_frequencies ; cat OUT_frequencies Frequency file is saved at: OUT_frequencies
#1 ((t2_1,t1_1),t3_1);  5991
#2 (t2_1,(t3_1,t1_1));  1969
#3 (t1_1,(t3_1,t2_1));  2040
#1 ((t2_1,t1_1),t3_1);  5991
#2 (t2_1,(t3_1,t1_1));  1969
#3 (t1_1,(t3_1,t2_1));  2040

gt <- read.tree("sim_num_gener")
png("case4.png", width = 1500, height =500)
par(mfrow=c(1,3))
hist(unlist(lapply(gt, function(x){sum(x$edge.length[which.edge(x,c("t1_1", "t2_1"))])})), breaks = 30, probability=T)
lines(dab(m = 2000, n=2000),col="red")

hist(unlist(lapply(gt, function(x){sum(x$edge.length[which.edge(x,c("t1_1", "t3_1"))])})), probability=T)
lines(dac(m = 2000, n=2000),col="red")

hist(unlist(lapply(gt, function(x){sum(x$edge.length[which.edge(x,c("t2_1", "t3_1"))])})), probability=T)
lines(dac(m = 2000, n=2000),col="red")

dev.off()