# rstudioapi::addTheme("https://raw.githubusercontent.com/roshandarji/synthwaveBLACK/master/SynthwaveBlack.rstheme", TRUE, TRUE, FALSE)
# para usar directorio donde se encuentra script
library(rstudioapi)
current_path = rstudioapi::getActiveDocumentContext()$path 
setwd(dirname(current_path ))
getwd()

library(tidyverse) # Easily Install and Load the 'Tidyverse'
library(phytools)  # Phylogenetic Tools for Comparative Biology (and Other Things)

# cargamos el arbol 

 eltree <- read.tree("./TREES/Accipitridae.tree") 
# eltree <- read.tree("./TREES/Fringillidae.tree") 
# eltree <- read.tree("./TREES/Muscicapidae.tree") 
# eltree <- read.tree("./TREES/Trochilidae.tree")  
# eltree <- read.tree("./TREES/Columbidae.tree")   
# eltree <- read.tree("./TREES/Furnariidae.tree")  
# eltree <- read.tree("./TREES/Strigidae.tree")    
# eltree <- read.tree("./TREES/Tyrannidae.tree")   

# hacemos un plot de lineajes en el tiempo
ltt.phylo(tree = eltree)

# ahora hacemos el mccr test [nos fijamos el numero de spp en nuestro arbol
# y el total de spp en la familia (ya lo tenemos anotado)]
# este total lo sustituimos abajo donde dice total
# de esta forma corregimos por esa falta de representatividad
total=300

mccr.result<-mccr(ltt(eltree),rho=length(eltree$tip.label)/total,nsim=100)
mccr.result 

# le agregamos una linea roja con lo esperado segun 
ln_ltt_plot <- ltt(eltree,log=TRUE)
lines(c(0, max(nodeHeights(eltree))), c(log(2), log(length(eltree$tip.label))),lty = "dashed", lwd = 2, col = "red")

# PODEMOS COMPARAR CON OTROS ARBOLES
# Acá podemos ver todos los LTTs
read.tree("./TREES/Accipitridae.tree") %>% ltt.phylo(main="LTT para Accipitriddae" )
read.tree("./TREES/Fringillidae.tree") %>% ltt.phylo(main="LTT para Fringilliddae" )
read.tree("./TREES/Muscicapidae.tree") %>% ltt.phylo(main="LTT para Muscicapiddae" )
read.tree("./TREES/Trochilidae.tree")  %>% ltt.phylo(main="LTT para Trochiliddae" )
read.tree("./TREES/Columbidae.tree")   %>% ltt.phylo(main="LTT para Columbiddae" )
read.tree("./TREES/Furnariidae.tree")  %>% ltt.phylo(main="LTT para Furnariiddae" ) 
read.tree("./TREES/Strigidae.tree")    %>% ltt.phylo(main="LTT para Strigiddae" )  
read.tree("./TREES/Tyrannidae.tree")   %>% ltt.phylo(main="LTT para Tyranniddae" ) 

# hacemos el plot de árbol más LTT
obj<-ltt(eltree,plot=FALSE)
plot(obj,show.tree=TRUE, log.lineages=FALSE,
     main="LTT para UN ARBOL",lwd=2)
#esto es necesario para lo que sigue...
h<-max(nodeHeights(eltree))
x<-seq(0,h,by=h/100)
b<-(log(Ntip(eltree))-log(2))/h

# simulamos 20 (VEINTE) arboles de tamaño similar a nuestro arbol 
# e.g. igual tamaño arbol y tiempo o height
# ponemos nsim=2 o 3 si se complica
trees <- pbtree(b=b,n=Ntip(eltree), t=max(nodeHeights(eltree)), nsim=20,method="direct",quiet=TRUE)

multobj<-ltt(trees,plot=FALSE) # aca hacemos multiples LTTs
# plot de 20 LTTs, 20 historias
plot(multobj,col="grey", main="LTT obtenido de ARBOL comparado a los 20 LTTs simulados")
## agregamos "nuestro" LTT
ltt(eltree,add=TRUE,lwd=2,col="#ef476f")

# AJUSTAMOS DOS MODELOS Puro nacimiento y nacimiento y muerte...
# Puros Nacimientos
yule.model<-fit.yule(eltree)
yule.model
# Nacimiento y muerte de linajes
bd.model<-fit.bd(eltree)
bd.model
# vemos si el modelo más complejo esta mas ajustado
#install.packages("lmtest")
#library(lmtest)
lrtest(yule.model,bd.model)

"Parte 4"
# leemos los datos de las regiones/distribuciones pero solo para Trochilidae
eltree <- read.tree("./TREES/Trochilidae.tree"  )
south <-as.matrix(read.csv("South.tsv",row.names=1))[,1] # como una matriz primero
south<-as.factor(south)                         # ahora lo convertimos en factor

read_tsv("SUR.tsv", show_col_types = FALSE) # vemos los datos mejor cargandolos asi

# troch.maps<-make.simmap(eltree,south, model="ARD",nsim=100)
# write.simmap(troch.maps, "troci.tre", format="phylip",version=1, append = TRUE)
troch.maps <-  read.simmap("trochi.tre",format="phylip",version=1)

describe.simmap(troch.maps) # describo las transiciones entre estados
# CUAL ES EL ESTADO MAS FRECUENTE Y QUE CAMBIOS SON MAS FRECUENTES

# ahora el troch.maps "se despedaza" y nos sirve de input

# a continuacion reusamos y adaptamos algo de código tomado desde aqui del creador de phytools 
# http://blog.phytools.org/2022/07/creating-lineage-through-time-plot.html

# plot de la historia del grupo...
# metodo simmap, reconstruccion ancestral de caracteres

colors<-setNames(c("#F0D32D","#2a9d8f","#C34582","#BEC0BF","#1F386E"),c("CS","MA","NAm","NO","SA" ))
plot(summary(troch.maps), colors=colors, ftype="i",fsize=0.1)
add.simmap.legend(colors=colors,x=2,y=80,prompt=FALSE)

## transformamos un poco los arboles...
tt<-map.to.singleton(troch.maps[[1]])
## del primer arbol sacamos el largo de cada rama
H<-nodeHeights(tt)
## sacamos altura o largo de rama (tiempo) de estas ramas
h<-max(H)-branching.times(tt)
## obtenemos estados (en este caso que region ocupa)
ss<-setNames(as.factor(names(tt$edge.length)), tt$edge[,2])
## matriz de cada linaje
lineages<-matrix(0,length(h),length(levels(south)),dimnames=list(names(h),levels(south)))
## los contamos...
for(i in 1:length(h)){
  ii<-intersect(which(h[i]>H[,1]),which(h[i]<=H[,2]))
  lineages[i,]<-summary(ss[ii])
}
## ordenamos eventos...
ii<-order(h)
times<-h[ii]
lineages<-lineages[ii,]

## creamos area de ploteo
plot(NA,xlim=range(times),ylim=c(0,max(lineages)),
     xlab="time",ylab="lineages",bty="n",las=1)
# leyenda
legend("bottomleft",c("Cono Sur","MesoAm", "NAm", "Norte y MesoAm", "SAm"),
       pch=22,pt.bg=c("#F0D32D","#2a9d8f","#C34582","#BEC0BF","#1F386E"),
       pt.cex=2.8,cex=1.6,bty="n")

## agregamos linajes de cada region
lines(times,lineages[,1],type="s",col="#F0D32D",lwd=3.5)
lines(times,lineages[,2],type="s",col="#2a9d8f",lwd=3.5)
lines(times,lineages[,3],type="s",col="#C34582",lwd=3.5)
lines(times,lineages[,4],type="s",col="#BEC0BF",lwd=3.5)
lines(times,lineages[,5],type="s",col="#1F386E",lwd=5)

# esta es la paleta... "#a9cfdb","#2a9d8f","#D526C3","#f4e3b8","#516de8"
# la pueden ver acá https://coolors.co/a9cfdb-2a9d8f-d526c3-e57010-516de8
## ponemos el arbol encima
cols<-setNames(make.transparent(c("#F0D32D","#2a9d8f","#C34582","#BEC0BF","#1F386E"),0.72),levels(south))
plot(troch.maps[[1]],cols,ftype="off",add=TRUE,lwd=1.5, mar=c(5.1,4.1,4.1,2.1))
obj<-markChanges(troch.maps[[1]],plot=FALSE)

# le agregamos lineas
for(i in 1:nrow(obj)) lines(rep(obj[i,"x"],2),c(0,obj[i,"y"]),lty="dotted",col=make.transparent("#84a59d",0.8))
# leyenda final