setwd("C:/Users/Mike/Dropbox/To do/Chapter 1/R")

require(ape)
require(phytools)
require(geiger)
require (mvMORPH)
require(devtools)
require(nlme)
require(qpcR)
require(reshape)
require(ouch)
require(calibrate)
require(surface)


tree<- read.nexus("Chapter 1_trees.nexus")
dat<- read.csv("Chapter 1_data_all PCs.csv", row.names=1)
dd<- as.matrix(dat)

TreeOnly <- setdiff(tree$tip.label,rownames(dat)) #checks to make sure the tips of the phylogeny matches the morpho data

DataOnly <- setdiff(rownames(dat), tree$tip.label) #checks to make sure the tips of the phylogeny matches the morpho data

#########################################
#Prune Tree
########################################

phyloTime <- drop.tip(tree,TreeOnly)
phyloTimeLadderized <- (ladderize(phyloTime)) 
phyloTimeLadderized <- rescale(phyloTimeLadderized, "depth", 1)
#phyloTimeLadderized <- chronos(phyloTimeLadderized)

tip.col<-as.vector(dat$color)
names(tip.col)<- rownames(dat)
cols<-c(tip.col[phyloTime$tip.label],rep("black",phyloTime$Nnode))
names(cols)<-1:(length(phyloTime$tip)+phyloTime$Nnode)
phylomorphospace(phyloTime, dat, label="off", control= list(col.node=cols))


#phyloTime$root.time <-max(dist.nodes(phyloTime))
phyloTime$root.time <-110
 geoscalePhylo(tree=ladderize(phyloTime),label.offset=0)

dis<- dtt(phyloTime, dat[,c(1,2)], index=c("avg.sq", "avg.manhattan", "num.states"),
mdi.range=c(0,1), nsim=100, CI=0.95, plot=TRUE, calculateMDIp=F)

#########SURFACE

dd<- dat[,1:2]

tree3<-nameNodes(phyloTimeLadderized)
olist<-convertTreeData(tree3,dd[,1:2])
	otree<-olist[[1]] 
	odata<-olist[[2]]

fwd<-surfaceForward(otree, odata, aic_threshold = -4.5, exclude = 0,verbose = FALSE, plotaic = FALSE)
k<-length(fwd)

fsum<-surfaceSummary(fwd)
	fsum$aics

bwd<-surfaceBackward(otree, odata, starting_model = fwd[[k]], aic_threshold = -4.5, only_best = FALSE, verbose = FALSE, plotaic = FALSE)

	bsum<-surfaceSummary(bwd)
	kk<-length(bwd)

surfaceTreePlot(tree3, bwd[[kk]], labelshifts = T)

pdf("Surface_tree_claudio_2.pdf")
surfaceTreePlot(tree3, bwd[[kk]], labelshifts = T)
dev.off()

##################################
		#Plot Morphospace
##################################

pdf("Surface_trait_claudio_2.pdf")
surfaceTraitPlot(dat, bwd[[kk]], whattraits = c(1,2))
dev.off()

