/working/checkVMAT.R

 1require(auteur)
 2bytens=function(x, num) sapply(0:num, function(y)x*(10^y))
 3vcv2 <- function(tree){
 4    res = cophenetic(tree) / 2
 5    max(res) - res
 6}
 7
 8n=bytens(20,2)
 9reps=5
10out=array(dim=c(length(n)*reps,5))
11for(i in n){
12	for(j in 1:reps){
13		cat(paste("working on rep", j, "of size", i, "\n",sep=" "))
14		phy=rcoal(i)
15		ii=which(n==i)
16		sA=system.time(vcv(phy)->vA)[3] 
17		sR=system.time(vmat(phy)->vR)[3] 
18		s2=system.time(vcv2(phy)->v2)[3]
19		match=all(vA==vR&vA==v2) 
20		out[(ii-1)*reps+j,]=c(match, sA,sR,s2,Ntip(phy))
21		}
22		
23		}
24out=data.frame(out)
25out=cbind(out, out[,2]/out[,3],out[,2]/out[,4])
26names(out)=c("match","ape","rjmcmc","klaus","treesize","speedupRJMCMC","speedupKLAUS")
27save(out,file="checkVMAT.rda")