mi = 0

# read the data 
contData <- readContinuousCharacterData("data/primates_lht.nex")

contData.excludeAll()
contData.includeCharacter(1) 
#contData.excludeCharacter(2:11)

# work under fixed topology
treeArray <- readTrees("data/Primates_tree.nex")
psi <- treeArray[1]

# sigma: variance per unit of time of the Brownian motion
logSigma ~ dnUniform(-5,5)
sigma := 10^logSigma

moves[++mi] = mvSlide(logSigma, delta=1.0, tune=true, weight=2.0)

# univariate Brownian process along the tree
# parameterized by sigma
logmass ~ dnPhyloBrownianREML(psi, branchRates=1.0, siteRates=sigma, nSites=1)

# let us have a look at the stochastic variable logmass
#logmass
# you can see the actual data by
#logmass.show()

logmass.clamp( contData )

# create the model
mymodel = model(sigma)

# on screen, we will monitor only the correlation coefficient and the mean value of each trait
monitors[1] = mnScreen(printgen=10000, sigma)

# a model monitor
monitors[2] = mnFile(filename="output/primates_mass_REML.log", printgen=10, separator = TAB, sigma)
 
mymcmc = mcmc(mymodel, monitors, moves)

mymcmc.burnin(generations=10000, tuningInterval=250)
mymcmc.run(300000)

#q()