rm(list=ls())
setwd(dirname(rstudioapi::getActiveDocumentContext()$path)) ##Set path to current dir
library(tidyverse) #Lib for data analysis


# List all .RDS files in your Data folder
rds_files <- list.files("Data", pattern = "\\.RDS$", full.names = TRUE)

# Read each file and assign it to a variable named after the file (without extension)
for (file in rds_files) {
  var_name <- tools::file_path_sans_ext(basename(file))  # Remove .RDS
  assign(var_name, readRDS(file))
}


# List of all your dataset names as strings
# Get all variable names in the environment that start with "Scenario1_"
Scenario1_names=ls(pattern = "^Fig2_Scenario1_")
Scenario2_names=ls(pattern = "^Fig2_Scenario2_")


Fig3_names <- ls(pattern = "^Fig3_")

Fig4_names <- ls(pattern = "^Fig4_")

#Fig5_names <- ls(pattern = "^Fig5_") #dispersal analysis, changing parameter dist_glob

Fig6_names <- ls(pattern = "^Fig6_")
# Years you want to subset and analyze
years <- 2



data_names <- c(Scenario1_names, Scenario2_names, Fig3_names, Fig4_names, Fig6_names)


# Function to process each dataset
process_dataset <- function(data_names) {
  data <- get(data_names)  # Get the actual data frame from its name
  results <- list()
  
  for (year in years) {
    subset_data <- subset(data, Year == year)
    n2_mean <- mean(data$Aphids_N2, na.rm = TRUE)
    n3_mean <- mean(data$Aphids_N3, na.rm = TRUE)
    p_mean <- mean(data$ill+data$well, na.rm=TRUE)
    
    # Store result in a named list
    results[[paste0(data_names, "_year", year)]] <- list(
      N2 = n2_mean,
      N3 = n3_mean,
      P = p_mean
    )
  }
  
  return(results)
}



# Apply the function to all datasets
all_results <- lapply(data_names, process_dataset)

# Flatten the list if needed
flat_results <- do.call(c, all_results)

Scenario1_N2 <- c()
Scenario1_N2 <- numeric(length(Scenario1_names))
for (i in seq_along(Scenario1_names)) {
  key <- paste0(Scenario1_names[i], "_year2")
  Scenario1_N2[i] <- flat_results[[key]][[1]]
}


Scenario1_N3 <- c()
Scenario1_N3 <- numeric(length(Scenario1_names))
for (i in seq_along(Scenario1_names)) {
  key <- paste0(Scenario1_names[i], "_year2")
  Scenario1_N3[i] <- flat_results[[key]][[1]]
}


Scenario2_N2 <- c()
Scenario2_N2 <- numeric(length(Scenario2_names))
for (i in seq_along(Scenario2_names)) {
  key <- paste0(Scenario2_names[i], "_year2")
  Scenario2_N2[i] <- flat_results[[key]][[1]]
}


Scenario2_N3 <- c()
Scenario2_N3 <- numeric(length(Scenario2_names))
for (i in seq_along(Scenario2_names)) {
  key <- paste0(Scenario2_names[i], "_year2")
  Scenario2_N3[i] <- flat_results[[key]][[1]]
}


Fig3_N2 <- c()
Fig3_N2 <- numeric(length(Fig3_names))
for (i in seq_along(Fig3_names)) {
  key <- paste0(Fig3_names[i], "_year2")
  Fig3_N2[i] <- flat_results[[key]][[1]]
}

Fig3_N3 <- c()
Fig3_N3 <- numeric(length(Fig3_names))
for (i in seq_along(Fig3_names)) {
  key <- paste0(Fig3_names[i], "_year2")
  Fig3_N3[i] <- flat_results[[key]][[1]]
}


Fig4_N2 <- c()
Fig4_N2 <- numeric(length(Fig4_names))
for (i in seq_along(Fig4_names)) {
  key <- paste0(Fig4_names[i], "_year2")
  Fig4_N2[i] <- flat_results[[key]][[1]]
}

Fig4_N3 <- c()
Fig4_N3 <- numeric(length(Fig4_names))
for (i in seq_along(Fig4_names)) {
  key <- paste0(Fig4_names[i], "_year2")
  Fig4_N3[i] <- flat_results[[key]][[1]]
}


Fig6_N2 <- c()
Fig6_N2 <- numeric(length(Fig6_names))
for (i in seq_along(Fig6_names)) {
  key <- paste0(Fig6_names[i], "_year2")
  Fig6_N2[i] <- flat_results[[key]][[1]]
}

Fig6_N3 <- c()
Fig6_N3 <- numeric(length(Fig6_names))
for (i in seq_along(Fig6_names)) {
  key <- paste0(Fig6_names[i], "_year2")
  Fig6_N3[i] <- flat_results[[key]][[1]]
}



#Fig2

x_scenarios=c(1,2,3,4,5,6,7,8, 9, 10 )
label_scenarios=c(0, 2.5, 5, 7.5, 10, 0, 2.5, 5, 7.5, 10)



#N2
pdf("Figures/Scenarios_N2.pdf", width=12, height=8.5)

#Plotting
par(mgp=c(4,2.3, 0), mar=c(10, 10, 2, 2))
plot(x_scenarios, Scenario1_N2, xlab="", ylab="", col='#a1d99b',  cex.axis=3, cex=4, pch=16, xaxt="n", ylim=range(0, 500))
axis(1, at=x_scenarios, labels=label_scenarios, cex.axis=2.8)
points(x_scenarios, Scenario2_N2, cex=4, col="#a1d99b", pch=17)

#x-axis lines and cross
abline(v=5.5, lwd=4)


#x-axis texts
mtext(c("Flower strip", "Woody habitat"), side=1, at=c(3, 8), cex=3.5, line=5)
mtext("Percentage in landscape cover", side=1, cex=3.5, line=7.8, at=5)
mtext(expression(paste("Aphid density " (m^{-2}))), side=2, cex=3.5, line=4.5, at=220)



dev.off()



#N3
pdf("Figures/Scenarios_N3.pdf", width=12, height=8.5)

#Plotting
par(new=TRUE)

layout(matrix(c(1,2,2,2), nrow=4, ncol=1, byrow=TRUE))

par(mgp=c(4,2.3, 0), mar = c(1*1.50, 8*1.50, 2*1.50, 2*1.50))
plot(x_scenarios, Scenario1_N3, xlab = "", ylab = "", col = '#238b45', yaxt="n",
     xaxt = "n", cex.axis = 3,  pch = 16, cex=4*1.5,
     ylim = c(1550, 1560))
axis(2, at = c(1550, 1560), label_scenarioss=FALSE, cex.axis=2.8*1.5)
mtext(c(1550,1560), side=2, at=c(1550,1558), cex=2.8, line=3)
points(x_scenarios, Scenario2_N3, cex=4, col="#238b45", pch=17)

#x-axis lines and cross
abline(v=5.5, lwd=4)

legend("topright", legend=c('Scenario 1', "Scenario 2"), cex=4, pch=c(16, 17),
       col=c('black'))


par(mar=c(10*1.50,8*1.50,0,2*1.50))
plot(x_scenarios, Scenario1_N3, xlab = "", ylab = "", col = '#238b45', yaxt="n", xaxt = "n",
     pch = 16, cex=4*1.5, ylim = c(0, 50))
axis(1, at = c(1,2,3,5,6,7), label_scenarioss = FALSE, cex.axis=3*1.50)  # Hide x-axis label_scenarioss
axis(2, at=c(10,20,30,40,50), label_scenarioss=FALSE, cex.axis=3*1.50)
points(x_scenarios, Scenario2_N3, cex=4, col="#238b45", pch=17)

axis(1, at=x_scenarios, label_scenarioss=label_scenarios, cex.axis=2.8*1.5)
mtext(c(10, 20,30,40), side=2, at=c(10,20,30,40), cex=3, line=3)

#x-axis lines and cross
abline(v=5.5, lwd=4)

#x-axis texts
mtext(c("Flower strip", "Woody habitat"), side=1, at=c(3, 8), cex=3.5, line=7)
mtext("Percentage in landscape cover", side=1, cex=3.5, line=11.5, at=5.5)
#mtext(expression(paste("Aphid density " (m^{-2}))), side=2, cex=3.5, line=4.5, at=3)


dev.off()







#Fig3

x_Fig3=c(1,2,3,4,5,6,7,8)
label_Fig3=c("A","B","C","D", "E", "F", "G", "H")


#N2
pdf("Figures/Fig3_N2.pdf", width=12, height=8.5)

#Plotting
par(mgp=c(4,2.3, 0), mar=c(10, 10, 2, 2))
plot(x_Fig3, Fig3_N2, xlab="", ylab="", col='#a1d99b',  cex.axis=3, cex=4, pch=16, xaxt="n", ylim=range(250, 500))
axis(1, at=x_Fig3, labels=label_Fig3, cex.axis=2.8)

abline(h= Fig3_N2[[3]], lwd=4, col="red")


mtext("Case", side=1, cex=3.5, line=6, at=4.5)
mtext(expression(paste("Aphid density " (m^{-2}))), side=2, cex=3.5, line=4.5, at=350)



dev.off()




#N3
pdf("Figures/Fig3_N3.pdf", width=12, height=8.5)

#Plotting
par(mgp=c(4,2.3, 0), mar=c(10, 10, 2, 2))
plot(x_Fig3, Fig3_N3, xlab="", ylab="", col='#238b45',  cex.axis=3, cex=4, pch=16, xaxt="n", ylim=range(500, 1600))
axis(1, at=x_Fig3, labels=label_Fig3, cex.axis=2.8)


abline(h= Fig3_N3[[8]], lwd=4, col="red")

mtext("Case", side=1, cex=3.5, line=6, at=4.5)
#mtext(expression(paste("Aphid density " (m^{-2}))), side=2, cex=3.5, line=4.5, at=200)



dev.off()





#Fig4


x_Fig4=c(1,2,3,4,5,6,7,8,9,10,11,12)
label_Fig4=c("Def","A","B","C", "D", "E", "F", "G", "H", "I", "J", "K")



#N2
pdf("Figures/Fig4_N2.pdf",width=14, height=8.5)

#Plotting
par(mgp=c(4,2.3, 0), mar=c(10, 10, 2, 2))
plot(x_Fig4, Fig4_N2, xlab="", ylab="", col='#a1d99b',  cex.axis=3, cex=4, pch=16, xaxt="n", ylim=range(300, 500))
axis(1, at=x_Fig4, labels=label_Fig4, cex.axis=2.8)

abline(h= Fig4_N2[[1]], lwd=4, col="red")

mtext("Case", side=1, cex=3.5, line=6, at=6)
mtext(expression(paste("Aphid density " (m^{-2}))), side=2, cex=3.5, line=6, at=400)



dev.off()



#N3
pdf("Figures/Fig4_N3.pdf", width=14, height=8.5)

#Plotting
par(mgp=c(4,2.3, 0), mar=c(10, 10, 2, 2))
plot(x_Fig4, Fig4_N3, xlab="", ylab="", col='#238b45',  cex.axis=3, cex=4, pch=16, xaxt="n", ylim=range(500, 1200))
axis(1, at=x_Fig4, labels=label_Fig4, cex.axis=2.8)


abline(h= Fig4_N3[[1]], lwd=4, col="red")

mtext("Case", side=1, cex=3.5, line=6, at=6)
#mtext(expression(paste("Aphid density " (m^{-2}))), side=2, cex=3.5, line=6.5, at=2)



dev.off()





#Fig6


x_Fig6=c(1,2,3,4)
label_Fig6=c("Default","A", "B", "C")



#N2
pdf("Figures/Fig6_N2.pdf", width=12, height=8.5)

#Plotting
par(mgp=c(4,2.3, 0), mar=c(10, 10, 2, 2))
plot(x_Fig6, Fig6_N2, xlab="", ylab="", col='#a1d99b',  cex.axis=3, cex=4, pch=16, xaxt="n", ylim=range(200, 350))
axis(1, at=x_Fig6, labels=label_Fig6, cex.axis=2.8)

abline(h= Fig6_N2[[1]], lwd=4, col="red")

mtext("Case", side=1, cex=3.5, line=6, at=2.5)
mtext(expression(paste("Aphid density " (m^{-2}))), side=2, cex=3.5, line=6, at=275)



dev.off()

#P
pdf("Figures/Fig6_P.pdf", width=12, height=8.5)

#Plotting
par(mgp=c(4,2.3, 0), mar=c(10, 10, 2, 2))
plot(x_Fig6, Fig6_P, xlab="", ylab="", col='#ef3b2c',  cex.axis=3, cex=4, pch=16, xaxt="n", ylim=range(0, 0.5))
axis(1, at=x_Fig6, labels=label_Fig6, cex.axis=2.8)

abline(h= Fig6_P[[1]], lwd=4, col="red")

mtext("Case", side=1, cex=3.5, line=6, at=2.5)
mtext(expression(paste("Hoverfly density " (m^{-2}))), side=2, cex=3.5, line=6, at=0.25)



dev.off()

#N3
pdf("Figures/Fig6_N3.pdf", width=12, height=8.5)

#Plotting
par(mgp=c(4,2.3, 0), mar=c(10, 10, 2, 2))
plot(x_Fig6, Fig6_N3, xlab="", ylab="", col='#238b45',  cex.axis=3, cex=4, pch=16, xaxt="n", ylim=range(0, 800))
axis(1, at=x_Fig6, labels=label_Fig6, cex.axis=2.8)

abline(h= Fig6_N3[[1]], lwd=4, col="red")

mtext("Case", side=1, cex=3.5, line=6, at=2.5)
# mtext(expression(paste("Aphid density " (m^{-2}))), side=2, cex=3.5, line=6.5, at=2)



dev.off()

