#setwd('~/Desktop/mimulus/')
library(phytools)
library(ape)

phy<-read.tree("mimtree4.tre")
d<-read.csv("mimulus_data.csv")
shape<-log(d$lxw)
names(shape)<-d$names
shape.sd<-1/50
library(diversitree)
p<-starting.point.quasse(phy,shape)
xr<-range(shape)+c(-1,1)*20*p["diffusion"]
linear.x<-make.linear.x(xr[1],xr[2])

make.shape<-function(lambda,mu)
make.quasse(phy,shape,shape.sd,lambda,mu)

nodrift<-function(f)
constrain(f,drift~0)


f.c<-make.shape(constant.x,constant.x)
f.l<-make.shape(linear.x,constant.x)
f.s<-make.shape(sigmoid.x,constant.x)
f.h<-make.shape(noroptimal.x,constant.x)
control<-list(parscale=0.1,reltol=0.001)
mle.c<-find.mle(nodrift(f.c),p,lower=0,control=control,verbose=0)
p.c <- mle.c$par 
p.l <- c(p.c[1], l.m=0, p.c[2:3]) 
p.s <- p.h <- c(p.c[1], p.c[1], mean(xr), 1, p.c[2:3]) 
names(p.s) <- argnames(nodrift(f.s)) 
names(p.h) <- argnames(nodrift(f.h))

mle.l<-find.mle(nodrift(f.l),p.l,control=control,verbose=0)
mle.s<-find.mle(nodrift(f.s),p.s,control=control,verbose=0)
mle.h<-find.mle(nodrift(f.h),p.h,control=control,verbose=0)
anova(mle.c,linear=mle.l,sigmoidal=mle.s,hump=mle.h)

coef(mle.c)
coef(mle.l)
coef(mle.s)
coef(mle.h)

mle.d.l<-find.mle(f.l,coef(mle.l,TRUE),control=control,verbose=0)
mle.d.s<-find.mle(f.s,coef(mle.s,TRUE),control=control,verbose=0)
mle.d.h<-find.mle(f.h,coef(mle.h,TRUE),control=control,verbose=0)
anova(mle.c, linear=mle.l, sigmoidal=mle.s, hump=mle.h, drift.linear=mle.d.l, drift.sigmoidal=mle.d.s, drift.hump=mle.d.h)