## Sterre Witte September 2023

##---------------------------------
## Packages etc.
library(tidyr)
library(ggplot2)
library(igraph)
library(NetIndices)
library(viridis)
library(ggnetwork)
library(dplyr)
source("Foodwebmetrics functions.R")

##---------------------------------
## Species data from samples

data <- read.csv("... a file containing abundances of different species found in several samples and years.csv", sep=",", head=T, fileEncoding="UTF-8-BOM")

## Find the unique species in all the different samples
speciesobserved <-data %>%
  group_by(sampleyear, sample, species_name)%>%
  dplyr::summarise(abundance=sum(abundance))

##---------------------------------
## Food web data

## Interaction matrix
interactions <- read.csv("Interaction Matrix.csv", head=T, sep= ";",  check.names = FALSE) #the interaction matrix tab of the interaction matrix
rownames(interactions) <-interactions[,1]
interactions[,1] <- NULL     #matrix1 <- subset(matrix1, select=-c(X)), dit is het voor een matrix
interactions1=as.matrix(interactions)

##---------------------------------
## Function subset for substrates

fwmetrics <- NULL
plot_list = list()

pdf("~... folder directory fir plots")# a folder to save your plots to
for(s in unique(speciesobserved$sampleyear)){
  for (i in unique(speciesobserved$sample_id)){
      # subset the matrix per cage
      data <- speciesobserved%>%
        filter(sampleyear==s & sample_id==i)
      
      if(nrow(data)==0) {
        print(paste("skipped", s,i,sep=",")); next}
      
      ## Add algae, plankton and POM to every foodweb as base
      Macroalgae <- data.frame(sampleyear=s, fieldlabel= i, species_name= "Macroalgae", abundance = 1, mobile= 1, benthic= 1)
      Phytoplankton <- data.frame(sampleyear=s, fieldlabel= i, species_name= "Phytoplankton", abundance = 1, mobile= 1, benthic= 1)
      Zooplankton <- data.frame(sampleyear=s, fieldlabel= i, species_name= "Zooplankton", abundance = 1, mobile= 1, benthic= 1)
      wPOM <- data.frame(sampleyear=s, fieldlabel= i, species_name= "wPOM", abundance = 1, mobile= 1, benthic= 1)
      sPOM <- data.frame(sampleyear=s, fieldlabel= i, species_name= "sPOM", abundance = 1, mobile= 1, benthic= 1)
      
      data <- rbind(data, Macroalgae, Phytoplankton, Zooplankton, wPOM, sPOM)
      
      species=unique(data$species_name)
      
      interact<- interactions%>%
        filter(row.names(.) %in% species)%>%
        dplyr::select(all_of(species))%>%
        relocate(any_of(names(interactions)))
     
      matrix=as.matrix(interact)
      mode(matrix) = "numeric"
      
      network <- graph.adjacency(matrix, mode="undirected")

      #calculate metrics 
      metrics <- Foodwebmetrics(s,i, interact, network)
      #matrix requires a square interaction dataframe, network an adjacency object
      
      # Using rbind() to append the output of one iteration to the dataframe
      fwmetrics = rbind(fwmetrics,metrics)
      
      
      ## Setup the plot
      #calculate Trophic level
      TL <- TrophInd(matrix) #trophic level (with cannibals included)
      TL$round <- as.factor(findInterval(TL$TL, c(0, 0.5,1, 1.5, 2, 2.5,3, 3.5, 4), all.inside=T))
      V(network)$TL<- TL$TL[match(V(network)$name, row.names(TL))]
      V(network)$round <- TL$round[match(V(network)$name, row.names(TL))]
      
      coordsp <- cbind(runif(length(V(network))), TL$TL)
      options(ggrepel.max.overlaps = Inf)

      #Save plot as pdf
      # figurepath <- file.path("~... folder directory for plot pdfs",
      #                        paste("Foodweb", s, i, c, ".pdf", sep = ""))
      # 
      # pdf(file=figurepath)
      
      print(
          ggplot(ggnetwork(network, layout=coordsp, scale = FALSE, cell.jitter = 0.25), aes(x, y, xend = xend, yend = yend))+
          geom_edges(size = 0.3, color = "black", arrow = arrow(length = unit(6,"pt"), type = "closed"), curvature = 0.1) +
          geom_nodes(aes(color = TL), size=8)+
          # geom_nodelabel_repel(aes(label=name, color= sg), size=2, fontface = "bold", box.padding = unit(0.1, "lines"))+
          scale_colour_viridis()+
          theme_classic()+
            # ylim(c(1,4))+
          coord_cartesian(clip = "off") +
          theme(legend.position="none", text = element_text(size = 20),axis.line.x = element_blank(),
                axis.title.x=element_blank(),axis.text.x=element_blank(),axis.ticks.x=element_blank())+
            labs(y="Trophic level")
          # labs(title = paste("Foodweb", s, i, c, sep = " "))
          )
        # dev.off()

      print(paste(s,i,sep=","))
    }
  }

dev.off()

## Write away the metrics
write.csv(fwmetrics, "Foodweb metrics.csv", row.names=F)