library(tidyverse)


#Figure 4 in Simulation 2


load("results/sim_2_mcrts_f.RData")
load("results/sim_2_mcrts_z.RData")
load("results/sim_2_mixed.RData")

plot_fig <- function(dat_m,dat_crt,dat_crt_f,to_do,index,weight,i){
  
  data_ <- data.frame(lag = c(0,1,2,3,4),
                      mixed = dat_m,
                      crtsf = dat_crt_f,
                      crts = dat_crt )
  if (index == 1){
      q <- "Quadratic"
  } else if (index==2){
      q <- "Exponential"
  }  else{
      q <- "Tanh"  
    }
  
  if (i==1){
    data_$title <- paste("No Interation",sep = "")
  }else{
    data_$title <- paste(q, " Interation",sep = "")
  }
  
  
  if(to_do == "coverage") {
    p <- ggplot(data_, aes(x=lag))+
      geom_hline(aes(yintercept = 0.90))+
      xlim(-0.1, 4.1)+  
      scale_y_continuous(breaks=c(0.00, 0.30, 0.60, 0.90),limits = c(0.0,1.0))+
      geom_line(aes(y = mixed,color = "Mixed"),size=2) +
      geom_line(aes(y = crtsf, color =  "MCRTs+F"),size=2) +
      geom_line(aes(y = crts, color =  "MCRTs+Z"),size=2) +
      geom_point(aes(y = mixed,color = "Mixed",shape = "Mixed"),size=5)+
      geom_point(aes(y = crtsf,color =  "MCRTs+F",shape = "MCRTs+F"),size=5)+
      geom_point(aes(y = crts,color =  "MCRTs+Z",shape = "MCRTs+Z"),size=5)+
      facet_grid(. ~ title)+
      scale_color_manual(name = NULL, values = c("Mixed" = "#999999", "MCRTs+F" = "#56B4E9", "MCRTs+Z" = "#E69F00"))+
      scale_shape_manual(name = NULL, values = c("Mixed" =15, "MCRTs+F" = 19, "MCRTs+Z" = 17))+
      xlab("Time lag") +
      ylab("Coverage rate")
    
      p <- p + theme_bw(base_line_size =0.5,base_rect_size = 1.0)
      p <-p + theme( legend.position="top",
      text = element_text(size=32, face="bold"),
      axis.title.x = element_text(color="black", size=32, face="bold",vjust = -2),
      axis.title.y = element_text(color="black", size=32, face="bold",vjust = 5),
      plot.margin=unit(c(1.0,1.0,1.0,1.0),"cm"))
  } else if (to_do == "length"){
    p <- ggplot(data_, aes(x=lag))+
      xlim(-0.1, 4.1)+  
      geom_line(aes(y = mixed,color = "Mixed"),size=2) +
      geom_line(aes(y = crtsf, color =  "MCRTs+F"),size=2) +
      geom_line(aes(y = crts, color =  "MCRTs+Z"),size=2) +
      geom_point(aes(y = mixed,color = "Mixed",shape = "Mixed"),size=5)+
      geom_point(aes(y = crtsf,color =  "MCRTs+F",shape = "MCRTs+F"),size=5)+
      geom_point(aes(y = crts,color =  "MCRTs+Z",shape = "MCRTs+Z"),size=5)+
      facet_grid(. ~ title)+
      scale_color_manual(name = NULL, values = c("Mixed" = "#999999", "MCRTs+F" = "#56B4E9", "MCRTs+Z" = "#E69F00"))+
      scale_shape_manual(name = NULL, values = c("Mixed" =15, "MCRTs+F" = 19, "MCRTs+Z" = 17))+
      xlab("Time lag") +
      ylab("Length of CIs")
      p <- p + theme_bw(base_line_size =0.5,base_rect_size = 1.0)
      p <-p + theme( legend.position="top",
                     text = element_text(size=32, face="bold"),
                     axis.title.x = element_text(color="black", size=32, face="bold",vjust = -2),
                     axis.title.y = element_text(color="black", size=32, face="bold",vjust = 5),
                     plot.margin=unit(c(1.0,1.0,1.0,1.0),"cm"))
  } else{
    stop("Error")
  }
  ggsave(p, width=6,height=6, dpi = 300, filename =paste("figures/sim_2/index_",index,"_weight_",weight,"_",to_do,".png"))
}


for (i in 1:length(list_mixed)){
  index<- list_mixed[[i]][1]
  weight<-list_mixed[[i]][2]
  dat_m <- list_mixed[[i]][3:7]
  dat_crt_f <-list_crts[[i]][3:7]
  dat_crt <-list_crts_f[[i]][3:7]
  plot_fig(dat_m,dat_crt,dat_crt_f,"coverage",index,weight,i)

  dat_m <- list_mixed[[i]][8:12]
  dat_crt_f <-list_crts[[i]][8:12]
  dat_crt <-list_crts_f[[i]][8:12]
  plot_fig(dat_m,dat_crt,dat_crt_f,"length",index,weight,i)
}

