source("https://tonyjhwueng.info/phyrates/ComparingBM.r")
library(ape)
library(geiger)

num.spe.array<-c(16,32,64,128)
sigma<-1
sigma1.array<-rep(sigma,5)
sigma2.array<-sigma*c(1.25,1.5,2.0,3.0,4.0)
simstrait<-1000

for(taxaIndex in 1: length(num.spe.array)){
  #taxaIndex<-1
  taxa<-num.spe.array[taxaIndex]
  print(taxa)
  #taxa<-compute.brlen(stree(taxa,"balanced"))#balanced tree
  tree<-rcoal(taxa)
  #plot(tree)
  for(s2Index in 1:length(sigma2.array)){
   # s2Index<-1
    sigma2<-sigma2.array[s2Index]
    
  for(traitIndex in 1:simstrait){
    #traitIndex<-1
#?sim.char 
  x1<-sim.char(phy=tree,par=sigma, model="BM")
  x2<-sim.char(phy=tree,par=sigma2, model="BM")
  x<-cbind(x1,x2)
  rownames(x)<-tree$tip.label
  assign(paste("taxa",taxa,"sgmratio",s2Index,"sim",traitIndex,sep="") ,CompareRates.multTrait(phy=tree,x=x,TraitCov=T,ms.err=NULL,ms.cov=NULL))
  }
 }
}

TypeIIerr<-array(0,c(length(num.spe.array),length(sigma2.array)))
rownames(TypeIIerr)<-num.spe.array
colnames(TypeIIerr)<-sigma2.array
TypeIIerr


for(taxaIndex in 1: length(num.spe.array)){
  #taxaIndex<-1
  taxa<-num.spe.array[taxaIndex]
#  print(taxa)
  #compute.brnlen(rtree(taxa,))
#  tree<-rcoal(taxa)
  #plot(tree)
  for(s2Index in 1:length(sigma2.array)){
    # s2Index<-1
    sigma2<-sigma2.array[s2Index]
    for(traitIndex in 1:simstrait){
      if(get(paste("taxa",taxa,"sgmratio",s2Index,"sim",traitIndex,sep=""))$Prob>0.05){
        TypeIIerr[taxaIndex,s2Index]<-TypeIIerr[taxaIndex,s2Index]+1
      }
    }
  }
}   
TypeIIerr/simstrait
save.image("simsBM.rda")

#nohup R CMD BATCH simsBM.R >/dev/null &

load(url("https://tonyjhwueng.info/phyrates/simsBM.rda"))

power.seq<-1-TypeIIerr/simstrait
power.seq

##      1.25   1.5     2     3     4
## 16  0.222 0.244 0.372 0.659 0.812
## 32  0.196 0.312 0.565 0.862 0.970
## 64  0.224 0.435 0.809 0.993 0.999
## 128 0.288 0.691 0.983 1.000 1.000

plot(power.seq[1,],type="b",ylim=c(0,1),xaxt="n",xlab="ratio",ylab="Power",main="BM")
axis(1, at=1:length(sigma2.array),labels=sigma2.array, las=2)
points(1:length(sigma2.array),power.seq[2,],type="b",col="blue",pch=3)
points(1:length(sigma2.array),power.seq[3,],type="b",col="red",pch=4)
points(1:length(sigma2.array),power.seq[4,],type="b",col="purple",pch=5)