---
title: ''
output:
  html_document: default
  pdf_document: default
---

```{r}
rm(list=ls())
setwd("~/Dropbox/FCU/Teaching/Mentoring/2022Spring-Current_Professor/2022Fall/YourueiMin/rcode_dcj/simulation")
# Load required libraries
#install.packages("ggplot2")
library(ggplot2)
library(reshape2)  # for melting the data
library(xtable) 
#install.packages("gridExtra")
library(gridExtra)
# wirte a function to plot the power curve

T1PowerCurve<-function(power_matrix=power_matrix,T1err.all=T1err.all,alt.phi=alt.phi,taxa.all=taxa.all,treetype=treetype){
  power <- data.frame(
    alt.phi = rep(alt.phi, length(taxa.all)),
    power = c(power_matrix),
    taxa = rep(taxa.all, each = length(alt.phi))
  )
power
type_I_error <- T1err.all
type_I_error
# Plot
p1<-ggplot(power, aes(x = alt.phi, y = power, color = as.factor(taxa), group = taxa)) +  geom_line(size = 1) +  geom_hline(yintercept = type_I_error, linetype = "dashed", color = "red", size = 0.5) +  labs(title = "Power Across Alt.Phi for Different Taxa",       x = "Alt.Phi",       y = "Power",       color = "Taxa") +   ggtitle(paste("T1 Power Curve for ",treetype," Tree",sep="")) + scale_color_discrete(name = "Taxa") +  theme_minimal()
print(p1)
return(p1)
}


treetype.array<-c("coal","purebirth","balanced","birthdeath")
taxa.all<-c(16,32,64,128)
power_matrix_full<-NULL
T1err.df_full<-NULL
for (treetype in treetype.array){
  filetoload<-paste("onerateT1power_",treetype,"Tree.RData",sep="")
  print(filetoload)
  load(filetoload)
  print(T1err.all)
  print(power_matrix)
  power_matrix_full<-rbind(power_matrix_full,power_matrix)
#  p1<-T1PowerCurve(power_matrix=power_matrix,T1err.all=T1err.all, alt.phi=alt_phi_values,taxa.all=taxa.all,treetype=treetype)  
#  print(p1)

#---given power_matrix of rows names taxa.all and colname phi.alt 
# plot the power curve using alt.phi  as x axis and power as y axis, plot the line each taxa across alt.phi, each line represent a taxa for each tree type
# add T1err.all as horizontal line, 
# make latex table for each tree type
# Reshape data for ggplot

power_melted <- melt(power_matrix, id.vars = "taxa.all", variable.name = "alt.phi", value.name = "power")
power_melted
colnames(power_melted)<-c("taxa.all","alt.phi","power")
T1err.df <- data.frame(taxa.all = names(T1err.all), T1err = unname(T1err.all))

# Convert taxa.all to a factor
power_melted$taxa.all <- as.factor(power_melted$taxa.all)
T1err.df$taxa.all <- as.factor(T1err.df$taxa.all)

T1err.df_full<-rbind(T1err.df_full,T1err.df)

# Define colors (you can change these to any colors you want)
my_colors <- rainbow(length(unique(power_melted$taxa.all)))

# Plot
assign(paste("p.",treetype,sep = ""),  
ggplot(power_melted, aes(x = as.numeric(as.character(alt.phi)), y = power, group = taxa.all)) +   geom_line(aes(color = taxa.all), size = 1) + geom_hline(data = T1err.df, aes(yintercept = T1err, color = taxa.all), linetype = "dashed", size = 0.5) +   scale_color_manual(values = my_colors) +   labs(title = paste("Power curve: ", treetype," tree",sep="") , x = expression(phi[a]), y = "Power value", color = "Taxa") +   theme_minimal()
)
}

grid.arrange(p.coal,p.purebirth,p.balanced,p.birthdeath,ncol=2)


fulltable<-cbind(T1err.df_full,power_matrix_full)

rownames(fulltable)<-NULL

fulltable
# make latex table for each tree type
xtable(fulltable,digits=3)
```




