### Script to test using different number of arguments using switch

##TOOL 1: phy to chrono

#needs a tree and information about the length of sequence data to make the tree
#rtree(10)->phy

# conducts search for best lambda
library(geiger)
library(getopt)

args<-commandArgs(trailingOnly=TRUE)
#print(args)
options<-matrix(c(
		'arg', 'a', 1, "character",
		'phy', 'p', 1, "character", 
		'lambda', 'l', 2, "integer", 
		'seqLength', 's', 2, "integer", 
		'saveTree', 't', 2, "boolean"),
		ncol=4,byrow=TRUE)

ret.opts<-getopt(options,args)

if (is.null(ret.opts$seqLength)){
	ret.opts$seqLength=100
}
if (is.null(ret.opts$lambda)){
	ret.opts$lambda=NA
}

#cat("The arguments passed\n")
#print(args)
#cat("The list of option and their values\n")
#print(ret.opts)

#arg<-args[1]
#ntax<-as.numeric(args[2])

load.file<-function(arg){	
	file<-switch(ret.opts$arg,
		fL = findLambda(ret.opts$phy, ret.opts$seqLength),
		cCPL = calcChronoPL(ret.opts$phy, ret.opts$lambda, ret.opts$seqLength, ret.opts$saveTree),
		both = findLambdaANDcalcChronoPL(ret.opts$phy, ret.opts$seqLength, ret.opts$saveTree)
	)
return(file)
} 



findLambda<-function(phyFile, seqLength=ret.opts$seqLength) {
	read.tree(phyFile)->phy
	for(i in 1:length(phy$edge.length)) {
		if(phy$edge.length[i]<=0.00001) {
			phy$edge.length[i]=0.0001
		}
		lambda<-10^(-3:6) 
		loglik<-c()
		for(i in 1:length(lambda)) {
			loglik[i]=attr(chronopl(phy,lambda=lambda[i], S=seqLength, tol=0.00001),"ploglik")
		}
		bestLambda=lambda[which(loglik==max(loglik))]
	}	
return(bestLambda)
}

# transforms branchlengths with best lambda
calcChronoPL<-function(phyFile, lambda, seqLength=100, saveTree=TRUE) {
	read.tree(phyFile)->phy
	chrono=chronopl(phy,lambda,S=seqLength, tol=0.00001, CV=T)
		#if(saveTree) {
			write.tree(chrono, file="chrono.phy") 
		#}
return(chrono)
}

#does both find best lambda and then runs calcChronoPL
findLambdaANDcalcChronoPL<-function(phyFile, seqLength=100, saveTree=TRUE) {
	read.tree(phyFile)->phy
	for(i in 1:length(phy$edge.length)) {
		if(phy$edge.length[i]<=0.00001) {
			phy$edge.length[i]=0.0001
		}
		lambda<-10^(-3:6) 
		loglik<-c()
		for(i in 1:length(lambda)) {
			loglik[i]=attr(chronopl(phy,lambda=lambda[i], S=seqLength, tol=0.00001),"ploglik")
		}
		bestLambda=lambda[which(loglik==max(loglik))]
	}	
	chrono=chronopl(phy,bestLambda,S=seqLength, tol=0.00001, CV=T)
	write.tree(chrono, file="chrono.phy") 

}

load.file()