library(dplyr)
library(ggplot2)
library(ggnewscale)
library(readr)



###### Read in data and build results DF #########


#Depending on preference you can run either the rds lines below or the csv lines. Both produce the same results in R, but
#we included the csv files as well in case anyone wants to look at the data using non-R tools. RDS files are R data 
#structures and so don't transport as easily between languages.


#Be sure the working directory is set to where the scripts are, or change the working directory below so that it reads
#into the data directory.

setwd("Data/")


file_list <- list.files(pattern = '*no_age.rds')


frame_list <- lapply(file_list, readRDS)


all_results_df <- bind_rows(frame_list, .id='column names')


file_list <- list.files(pattern = '*.csv')


frame_list <- lapply(file_list, read_csv)


all_results_df <- bind_rows(frame_list, .id='column names')


names(frame_list[[1]])
######### Summarize data post T2000, grouped by mortality and iteration $$$$$$$






gradient_means_post_2000 <- all_results_df%>%
  group_by(Mort, iteration)%>%
  filter(plants.generation >2000)%>%
  summarise(plant_mean_diff = mean(plants.sel_differential),
            plant_var_diff = var(plants.sel_differential),
            poll_mean_diff = mean(polls.sel_differential),
            poll_var_diff = var(polls.sel_differential),
            plant_slope = mean(plant_slope),
            poll_slope = mean(poll_slope),
            plant_int_means = mean(plants.int_means),
            poll_int_means = mean(polls.int_means),
            mean_diff = mean(plants.trait_means - polls.trait_means),
            mean_demo_change = mean(demo_change),
            var_demo_change = var(demo_change),
            scaled_demo_change = mean(demo_change)/3000,
            scaled_var_demo_change = var(demo_change)/3000)


names(gradient_means_post_2000)


###### Plot summaries across Mortality rates #######



grad_means_across_M <- ggplot(data=gradient_means_post_2000, aes(x=Mort))+
  geom_point(aes(y=plant_mean_diff, colour="Species L"))+geom_point(aes(y=poll_mean_diff, colour="Species S"))+ 
  theme_bw() + theme(panel.border = element_blank(), 
                     panel.grid.major = element_blank(),
                     panel.grid.minor = element_blank(), 
                     axis.line = element_line(colour = "black"))+
  ylab("Mean Selection Differential")+xlab(expression(paste("Annual Mortality of Species ", italic("L"))))+
  labs(colour="Species")+
  scale_colour_discrete(type=c('orange', 'purple'))+
  ggtitle("Mean Selection Differential", subtitle = "64 Simulations across M  @ T 2000 - 6000")
grad_means_across_M

ggsave("Selection_differential_means_across_M.pdf")



log_grad_var_across_A <- ggplot(data=gradient_means_post_2000, aes(x=Mort))+
  geom_point(aes(y=log(plant_var_diff), colour="Species L"))+geom_point(aes(y=log(poll_var_diff), colour="Species S"))+ 
  theme_bw() + theme(panel.border = element_blank(), 
                     panel.grid.major = element_blank(),
                     panel.grid.minor = element_blank(), 
                     axis.line = element_line(colour = "black"))+
  ylab("Log of Variance of Selection Differential")+xlab(expression(paste("Annual Mortality of Species", italic("L"))))+
  labs(colour="Species")+
  scale_colour_discrete(type=c('orange', 'purple'))+
  ggtitle("Log Variance of Selection Differntial, within simulations", subtitle = "64 Simulations across M @ T 2000 - 6000")
log_grad_var_across_A


ggsave("Log_selection_differential_variance_across_M.pdf")



raw_slopes_across_A <- ggplot(data=gradient_means_post_2000, aes(x=Mort))+
  geom_point(aes(y=plant_slope, colour="Species L"))+geom_point(aes(y=poll_slope, colour="Species S"))+ 
  theme_bw() + theme(panel.border = element_blank(), 
                     panel.grid.major = element_blank(),
                     panel.grid.minor = element_blank(), 
                     axis.line = element_line(colour = "black"))+
  ylab("Trait Slope")+xlab(expression(paste("Annual Mortality of Species ", italic("L"))))+
  labs(colour="Species")+
  scale_colour_discrete(type=c('orange', 'purple'))+
  ggtitle("Trait Slopes", subtitle = "10 Simulations across M @ T 2000 - 6000")

raw_slopes_across_A

ggsave("Raw_slopes_M1_to_M10.pdf")

fit_means_across_A <- ggplot(data=gradient_means_post_2000, aes(x=Mort))+
  geom_point(aes(y=plant_int_means, colour="Species L"))+geom_point(aes(y=poll_int_means, colour="Species S"))+ 
  theme_bw() + theme(panel.border = element_blank(), 
                     panel.grid.major = element_blank(),
                     panel.grid.minor = element_blank(), 
                     axis.line = element_line(colour = "black"))+
  ylab("Mean Fitness Values")+xlab(expression(paste("Annual Mortality of Species ", italic("L"))))+
  labs(colour="Species")+
  scale_colour_discrete(type=c('orange', 'purple'))+
  ggtitle("Mean Fitness values", subtitle = "10 Simulations across M @ T 2000 - 6000")

fit_means_across_A
ggsave("Mean_Fitness_Values_across_M.pdf")




Trait_differential_across_A <- ggplot(data=gradient_means_post_2000, aes(x=Mort))+
  geom_point(aes(y=mean_diff))+ theme_bw() + theme(panel.border = element_blank(), 
                                                   panel.grid.major = element_blank(),
                                                   panel.grid.minor = element_blank(), 
                                                   axis.line = element_line(colour = "black"))+
  ylab("Mean Trait Differential")+xlab(expression(paste("Annual Mortality of Species ", italic("L"))))+
  ggtitle("Mean Trait Differential", subtitle = "64 Simulations across M @ T 2000 - 6000")
Trait_differential_across_A

ggsave("Trait_differential_across_M.pdf")







########## Whole simulation plots - Post T2000 #######


P2_all_results_df <- all_results_df %>%
  filter(plants.generation>2000) %>%
  filter(Mort %in% c(0.1, 0.5, 1))%>%
  select(plants.generation, plants.trait_means, plants.int_means, plants.sel_differential,
         polls.trait_means, polls.int_means, polls.sel_differential, Mort, iteration)


names(P2_all_results_df)

p2000_pooled_mean_plots <- ggplot(data=P2_all_results_df, aes(x=plants.generation))+
  geom_line(aes(y=plants.trait_means, group=iteration, colour="Species L"), alpha=0.7)+
  geom_line(aes(y=polls.trait_means, group=iteration, colour="Species S"), alpha=0.7)+
  ylab("Mean Trait Value")+xlab("Time Step")+
  facet_wrap(~Mort, labeller = label_bquote(bold("m"[L] ~"="~ .(Mort))))+ theme_bw() + theme(panel.border = element_blank(), 
                                        panel.grid.major = element_blank(),
                                        panel.grid.minor = element_blank(), 
                                        axis.line = element_line(colour = "black"),
                                        axis.text.x = element_text(angle = 45, vjust=1,hjust=1),
                                        strip.text = element_text(size=13))+
  labs(colour="Species")+
  scale_colour_discrete(type=c('orange', 'purple'))
p2000_pooled_mean_plots
ggsave("P2000_Trait_mean_over_time_and_M.pdf")



p2000_pooled_diff_plots <- ggplot(P2_all_results_df, aes(x=plants.generation))+
  geom_line(aes(y=plants.trait_means-polls.trait_means, group=iteration), alpha=0.5)+
  ylab("Trait Differential")+xlab("Time Step")+
  facet_wrap(~Mort, labeller = label_bquote(bold("m"[L] ~"="~ .(Mort))))+ theme_bw() + theme(panel.border = element_blank(), 
                                        panel.grid.major = element_blank(),
                                        panel.grid.minor = element_blank(), 
                                        axis.line = element_line(colour = "black"),
                                        axis.text.x = element_text(angle = 45, vjust=1,hjust=1),
                                        strip.text = element_text(size=13))
  

p2000_pooled_diff_plots
ggsave("P2000_Trait_diff_over_time_and_M.pdf")





p2000_pooled_int_plots <- ggplot(data=P2_all_results_df, aes(x=plants.generation))+
  geom_line(aes(y=plants.int_means, group=iteration, colour="Species L"),  alpha=0.7)+
  geom_line(aes(y=polls.int_means, group=iteration, colour="Species S"),  alpha=0.7)+
  ylab("Mean Fitness Value")+xlab("Time Step")+
  facet_wrap(~Mort, labeller = label_bquote(bold("m"[L] ~"="~ .(Mort))))+ theme_bw() + theme(panel.border = element_blank(), 
                                        panel.grid.major = element_blank(),
                                        panel.grid.minor = element_blank(), 
                                        axis.line = element_line(colour = "black"),
                                        axis.text.x = element_text(angle = 45, vjust=1,hjust=1),
                                        strip.text = element_text(size=13))+
  labs(colour="Species")+
  scale_colour_discrete(type=c('orange', 'purple'))

p2000_pooled_int_plots
ggsave("P2000_Interaction_over_time_and_M.pdf")


p2000_pooled_sel_diff_plots <- ggplot(data=P2_all_results_df, aes(x=plants.generation))+
  geom_line(aes(y=plants.sel_differential, group=iteration, colour="Species L"), alpha=0.7)+
  geom_line(aes(y=polls.sel_differential, group=iteration, colour="Species S"), alpha=0.7)+
  ylab("Mean Selection Differential")+xlab("Time Step")+
  facet_wrap(~Mort, labeller = label_bquote(bold("m"[L] ~"="~ .(Mort))))+ theme_bw() + theme(panel.border = element_blank(), 
                                                                                             panel.grid.major = element_blank(),
                                                                                             panel.grid.minor = element_blank(), 
                                                                                             axis.line = element_line(colour = "black"),
                                                                                             axis.text.x = element_text(angle = 45, vjust=1,hjust=1),
                                                                                             strip.text = element_text(size=13))+
  labs(colour="Species")+
  scale_colour_discrete(type=c('orange', 'purple'))

p2000_pooled_sel_diff_plots

ggsave('P2000_sel_diff_over_time_and_M.pdf')





# 
# 
#