###Script analysis generated by RevScripter###

# setwd("E:\\Programas\\Analise filogenetica\\RevBayes\\Datasets\\Squamates(S18)") 

data = readDiscreteCharacterData("./revbayes_datafile.nex")

# Get some useful variables from the data. We need these later on.
num_taxa <- data.ntaxa()
num_branches <- 2 * num_taxa - 3
taxa <- data.taxa()

moves    = VectorMoves()
monitors = VectorMonitors()

num_char_states <- 4

######################
# Substitution Model #
######################

er_GTR ~ dnDirichlet(v(1, 1, 1, 1, 1, 1))
moves.append(mvBetaSimplex(er_GTR))
pi_GTR ~ dnDirichlet(v(1, 1, 1, 1))
moves.append(mvBetaSimplex(pi_GTR))
Q := fnGTR(er_GTR,pi_GTR)

site_rates_shape ~ dnUniform(0,1E7)
moves.append(mvScale(site_rates_shape))
num_rate_categories <- 4
site_rates := fnDiscretizeGamma(site_rates_shape, site_rates_shape, num_rate_categories)

##############
# Tree Model #
##############

topology ~ dnUniformTopology(taxa)
moves.append( mvNNI(topology, weight=num_taxa/2.0) )
moves.append( mvSPR(topology, weight=num_taxa/10.0) )

branch_hypershape ~ dnUniform(0,1E7)
moves.append(mvScale(branch_hypershape))

for(i in 1:num_branches){
   branch_lengths[i] ~ dnExponential(branch_hypershape)
    moves.append(mvScale(branch_lengths[i]))
}

psi := fnTreeAssembly(topology, branch_lengths)

phylo ~ dnPhyloCTMC(tree=psi, Q=Q, type="DNA", siteRates=site_rates)

phylo.clamp(data)

mymodel =  model(topology)

########
# MCMC #
########

monitors.append( mnModel(filename="SquamMol_Log", printgen=100) )
monitors.append( mnFile(filename="SquamMol_Trees", printgen=100, psi) )

monitors.append( mnScreen(printgen=100) )

mymcmc = mcmcmc(mymodel, monitors, moves, nruns=2, nchains=4, deltaHeat=0.05)
mymcmc.run(generations=1000)

# print the summary of the operators (now tuned)
mymcmc.operatorSummary()