area_fn = "data/earth25.still.atlas.txt"
data_fn = "data/primates_bg_n25.nex"
tree_fn = "data/primates.tree"
out_fn = "output/bg_3"

psi <- readTrees( tree_fn )[1]

data = readDiscreteCharacterData( data_fn )

atlas = readAtlas( area_fn)
n_areas = atlas.nAreas()

mi = 0

ngen = 5000

write(psi,filename=tree_fn+".index.tre")

clock_bg ~ dnExponential(10)
moves[++mi] = mvScale(clock_bg, lambda=0.5, weight=5.0)

glr ~ dnDirichlet([1,1])
moves[++mi] = mvSimplexElementScale(glr, alpha=30.0, weight=5.0)

r_gain := glr[1]
r_loss := glr[2]

q_area := fnFreeBinary(glr)

dp ~ dnExponential(10.0)
moves[++mi] = mvScale(x=dp, lambda=0.5, tune=true, weight=5.0)

grm := fnBiogeoGRM(atlas=atlas, distancePower=dp, useDistance=true)

q_range := fnBiogeoDE(gainLossRates=q_area, branchRates=clock_bg, geoRateMod=grm,
numAreas=n_areas, forbidExtinction=true)

clado_prob ~ dnDirichlet( [1, 1, 1] )
widespread_sympatry := clado_prob[1]

subset_sympatry := clado_prob[2]
allopatry := clado_prob[3]
moves[++mi] = mvSimplexElementScale(clado_prob, alpha=20.0, weight=5.0)

m ~ dnPhyloDACTMC(tree=psi, Q=q_range, cladoProbs=clado_prob, type="Biogeo",
forbidExtinction=true, useCladogenesis=true)

m.clamp(data)
m.lnProbability()

moves[++mi] = mvCharacterHistory(ctmc=m, qmap=q_range, tree=psi, lambda=0.05, type="Biogeo", graph="node", proposal="rejection", weight=100.0)
moves[++mi] = mvCharacterHistory(ctmc=m, qmap=q_range, tree=psi, lambda=0.05, type="Biogeo", graph="branch", proposal="rejection", weight=100.0)

moves[++mi] = mvCharacterHistory(ctmc=m, qmap=q_range, tree=psi, lambda=1.0, type="Biogeo", graph="node", proposal="rejection", weight=40.0)
moves[++mi] = mvCharacterHistory(ctmc=m, qmap=q_range, tree=psi, lambda=1.0, type="Biogeo", graph="branch", proposal="rejection", weight=40.0)


my_model = model(m)

monitors[1] = mnScreen(clock_bg, r_gain, r_loss, dp, subset_sympatry, allopatry,
widespread_sympatry, printgen=100)
monitors[2] = mnFile(clock_bg, r_gain, r_loss, dp, subset_sympatry, allopatry,
widespread_sympatry, filename=out_fn+".params.txt", printgen=10)

monitors[3] = mnCharHistoryNewick(filename=out_fn+".events.txt", ctmc=m, tree=psi,
printgen=100, style="events")
monitors[4] = mnCharHistoryNewick(filename=out_fn+".counts.txt", ctmc=m, tree=psi,
printgen=100, style="counts")

monitors[5] = mnCharHistoryNhx(filename=out_fn+".phw.txt", ctmc=m, tree=psi, atlas=
atlas, samplegen=100, maxgen=ngen, burnin=0.5)

my_mcmc = mcmc(my_model, monitors, moves)

my_mcmc.run(generations=ngen)








q()