###################################################################################################################
# Sterre Witte
# June 2023

Foodwebmetrics=function(year, sample, community, mat, network){

  ## Required packages
  
  
##########################################################################################
#calculate metrics
	
# TC=matrix(data=NA,sp,sp)
	
#species number
sp=length(mat)

#total number of links in the food web
totallinks=sum(mat) 

# #mean links per species
predators <- rowSums(mat)
sources <-  colSums(mat)
meanlinks <- mean(sources+predators)

#linkdensity: number of links per species
linkdens=totallinks/sp 

#connectance: saturation of the network
connectance=totallinks/(sp*sp) 
  
#vulnerability:
vulnerability=mean(predators)

#generality: number of consumed species for all species
generality=mean(sources)

#topspecies
tops <- which(predators==0)
top =length(tops) # number of top species in matrix (species without consumers)
topspecies=top/sp

#basal species
basal <- which(sources==0)#define basal species
base=length(basal) # number of basal species  (species without prey) --> !!hoe om te gaan met POM etc?
basalspecies=base/sp # share of basal species

#intermediate species
int=sp-top-base# number of intermediate species (species with both consumers and prey)
intermediatespecies=int/sp

#Chain Length: number of links from any species to a basal species. 
chain=data.frame(species=colnames(mat), level=0)
lp=1 # start loop
lower=as.data.frame(as.list(basal))#define the starting level
repeat{
  consumed<- mat%>%
    filter(rownames(mat) %in% colnames(lower))%>%
    select_if(colSums(.) != 0) 
  sp_consumed <- colnames(consumed)
  if (identical(consumed,lower)) {
    break
  }
  chain[match(sp_consumed, chain$species),]$level=ifelse(chain[match(sp_consumed, chain$species),]$level==0 , lp , chain[match(sp_consumed, chain$species),]$level) #fill in chainlength for any consumer feeding on the lower level
  lower=consumed # for level+1, level becomes lower
  lp=lp+1
}

chainlength=mean(chain$level) #!! Why mean, and not max?

#herbivores <- & detritivores (primary consumers)
eatbasal <- mat %>%
  filter(rownames(mat) %in% names(basal)) %>%
  select_if(colSums(.) != 0)

herbivore <- mat %>%
  filter(rownames(mat) %in% names(eatbasal)) %>%
  select_if(colSums(.) == 0 & !colnames(mat) %in% names(basal))

#omnivores
omnivore <- mat %>%
  filter(rownames(mat) %in% names(eatbasal)) %>%
  select_if(colSums(.) != 0 & !colnames(mat) %in% names(basal))

#carnivores
carnivore <- mat %>%
  filter(!rownames(mat) %in% names(basal)) %>%
  # select_if(!colnames(mat) %in% names(basal) & !colnames(mat) %in% colnames(omnivore) & !colnames(mat) %in% colnames(herbivore))
  select_if(colSums(.) != 0)

#cannibals
cannibal <- colnames(mat[,which(diag(as.matrix(mat))==1)]) # Which species have cannibalistic interaction
# mat %>%
#   select(which(row(mat)==col(mat)))
# 
# cannibal1 <- which(col(mat) == row(mat) & mat==1)

#   H=NULL
#   O=NULL
#   C=NULL
#   cann=0
#   herb=0
#   omn=0
#   carn=0
#   
# for (i in 1:sp){
#   anybase=any(mat[chn1,i]==1) # find species that eat basal species
#   noother=all(mat[-chn1,i]==0)		#find which eat nothing but basal species
#   H[i]=ifelse(anybase==T&noother==T,1,0)  #Which are herbivores: species only eating basal species !! this also includes detritivores
#   herb=ifelse(anybase==T&noother==T,herb+1,herb) # herbivores should eat basal species, nothing else
#   O[i]=ifelse(anybase==T&noother==F,1,0) #which are omnivores
#   omn=ifelse(anybase==T&noother==F,omn+1,omn)  # omnivores should eat a basal and a higher level food
#   C[i]=ifelse(H[i]==0&O[i]==0&con[i]>0,1,0) #which are carnivores
#   carn=sp-herb-omn-length(which(con==0)) #carnivores , everything besides basal, herbivores and omnivores. !! Why? Shouldn't this be the species that eat only other animals?
#   cann=ifelse(mat[i,i]==1,cann+1,cann) #number of cannibals of matrix
# }

herbivores=length(herbivore)/sp
omnivores=length(omnivore)/sp
carnivores=length(carnivore)/sp
cannibals=length(cannibal)/sp

#Trophic level
TL=TrophInd(mat)$TL #trophic level (with cannibals included)
trophiclevel=mean(TL)	#!! why mean?

# #Similarity: shared predators + shared prey/ total predators + total prey of species i and j 
# similarity=matrix(NA,nrow=dim(mat)[1],ncol=dim(mat)[1]) #define matrix to store similarity scores
#   for (i in 1:dim(mat)[1]){
#   for(j in 1:dim(mat)[2]){
#     if (i==j){ # similarity of species i with itself is 0
#       similarity[i,j]=0
#     }
#     else{		
#       conni=which(mat[i,]==1) #consumers of species i
#       connj=which(mat[j,]==1) #consumers of species j
#       comcon=length(which(conni%in%connj)) #shared consumers
#       srci=which(mat[,i]==1) #sources of species i
#       srcj=which(mat[,j]==1) #sources of species j
#       comsrc=length(which(srci%in%srcj)) #shared sources
#       tot=length(unique(c(srci,srcj)))+length(unique(c(conni,connj))) #total unique sources and consumers of i and j
#       similarity[i,j]=(comsrc+comcon)/tot #similarity index of species i and species j
#     }
#   }
# }
# 
# maxsimilarity=mean(apply(similarity, 2, max)) #maximum similarity of each species
# meansimilarity=mean(apply(similarity, 2, mean)) #mean similarity of each species
# 
# #Cluster Coefficient: how much neighboring species are linked with each other compared with maximum possible links
# 
# CC <- numeric(sp)
# Madj <- replace(mat, col(mat) == row(mat) & mat==1, 0)
# 
# # mat%>%
# #   rowwise() %>% 
# #   filter(any(c_across() == 1))
# 
# 
# ClusCoeff=function(mat){
#   CC=vector(length=dim(mat)[1]) #define vector to store data
#   Madj=mat*(matrix (1,dim(mat)[1],dim(mat)[1])-diag(dim(mat)[1])) #exclude cannibalism
#   for (i in 1:dim(Madj)[1]){ #run by every species in the matrix
#     #define neighbours
#     neighsrc=which(Madj[i,]==1) # which are the sources of species i
#     neighcon=which(Madj[,i]==1) # which are consumers of species i
#     neigh=unique(sort(c(neighsrc,neighcon))) # define unique neighbours
#     M=Madj[neigh,neigh] # make a submatrix including only the direct neighbours of species i
#     curv=sum(M)	#total actual links between neighbours
#     if(curv==0){Clus=0}
#     else{
#       for( j in 1:dim(M)[1]){
#         for (k in 1 :dim(M)[1]){
#           curv=ifelse(M[j,k]==1&M[k,j]==1,curv-0.5,curv)#curves are links between species, doesn't matter which way they go.(wattz& strogatz,nature 1998) So links going both ways have to be halved. 
#         }
#       }
#       
#       Clus=curv/((dim(M)[1]^2-dim(M)[1])/2) #clustering coefficient: divide actual links between neighbours by amount of possible links.
#     }
#     CC[i]=Clus #store clustering coefficient of species i
#   }
#   ClustCoeff=mean(CC) #Clustering coefficient of whole matrix is average of clustering coefficient in 
#   return (ClustCoeff)
#   
# }
# CC=ClusCoeff(mat)
# clustering=mean(CC) 

#   # CC=vector(length=dim(mat)[1]) #define vector to store data
#   Madj=mat*(matrix (1,dim(mat)[1],dim(mat)[1])-diag(dim(mat)[1])) #exclude cannibalism !! Why do you do that here?
#   for (i in 1:dim(Madj)[1]){ #run by every species in the matrix
#     #define neighbours
#     neighsrc=which(Madj[i,]==1) # which are the sources of species i
#     neighcon=which(Madj[,i]==1) # which are consumers of species i
#     neigh=unique(sort(c(neighsrc,neighcon))) # define unique neighbors
#     M=Madj[neigh,neigh] # make a submatrix including only the direct neighbors of species i
#     curv=sum(M)	#total actual links between neighbors
#     if(curv==0){Clus=0}
#     else{
#       for( j in 1:dim(M)[1]){
#         for (k in 1 :dim(M)[1]){
#           curv=ifelse(M[j,k]==1&M[k,j]==1,curv-0.5,curv)#curves are links between species, doesn't matter which way they go.(wattz& strogatz,nature 1998) So links going both ways have to be halved. 
#         }
#       }
#       
#       Clus=curv/((dim(M)[1]^2-dim(M)[1])/2) #clustering coefficient: divide actual links between neighbours by amount of possible links.
#     }
#     CC[i]=Clus #store clustering coefficient of species i
#   }
# clustering=mean(CC) #Clustering coefficient of whole matrix is average of clustering coefficient in 

#Path length: Average path length between all species.
disthist <- path.length.hist(network, directed=FALSE)$res
diameter <- length(disthist)
pathlength <- weighted.mean(1:diameter, disthist) #Average path length between nodes

#Compartmentalisation
multilevel=cluster_louvain(network)
compartmentalisation=modularity(network,membership(multilevel))

#Save all metrics in a vector
metrics=c('Year'=year, 'Field.ID'=sample, 'Community'=community, 'Species.number'=sp,'Totallinks'=totallinks, 'Meanlinks'=meanlinks, 'Linkdensity'=linkdens, 
                   'Connectance'=connectance, "Vulnerability"=vulnerability,"Generality"=generality,
                   'Topspecies'=topspecies,'Intermediatespecies'=intermediatespecies,'Basalspecies'=basalspecies,
                   'Herbivores'=herbivores,'Omnivores'=omnivores,'Carnivores'=carnivores,"Cannibals"=cannibals,
                   'Chain Length'=chainlength,"Trophic Level"=trophiclevel,
                   #'Max Similarity'=maxsimilarity,'Mean Similarity'=meansimilarity,
                   #'Clustering'=clustering,
                   'Path Length'= pathlength,'Compartmentalisation'=compartmentalisation)

return(metrics)
}

# # additional information about the food web #
# 
# #number of interactions per species
# totlinks=src+con[1:length(src)] 
# names(totlinks)=specnames
# hub=which(totlinks==max(totlinks)) #define species with highest number of interactions (the hub)
# hub=data.frame(paste(names(hub),"(",max(totlinks),")",sep=""))
# names(hub)[1]="Species with highest number of interactions (Number of links)"
# 
# top1=ifelse(src==0,1,0)
# basal1=ifelse(con==0,1,0)
# int1=ifelse(src==0|con==0,0,1)
# cann1=data.frame('Number of cannibales'=cannibal) 
# maxsi=apply(simil, 2, max)
# avsi=apply(simil, 2, mean)
# dat=data.frame('Foodweb_ID'=FW_ID,'Filename'=name,
# 'Top'=top1,'Intermediate'=int1,'Basal'=basal1,'Herbivore'=H,"Omnivore"=O,"Carnivore"=Carniv,'Cannibal'=diag(mat)==1,'Chain lengths'=ChnLg,"Trophic Level"=TL,'Number of sources-Generality'=src,
# 'Times eaten-Vulnerability'=con,'Total links'=totlinks,'Species with most links-Hub'=totlinks==max(totlinks),'Max.Similarity'=maxsi,'Average Similarity'=avsi)
# rownames(dat)=specnames
# space=data.frame('')
# 
# #save additional information as a csv
# write.table(mat, paste(filestem,"Matrix_info/",FW_ID,".csv",sep=''), col.names=NA,row.names=TRUE,sep=";")
# write.table(space, paste(filestem,"Matrix_info/",FW_ID,".csv",sep=''), col.names=F,row.names=F,sep=";",append=TRUE)
# write.table(dat, paste(filestem,"Matrix_info/",FW_ID,".csv",sep=''), col.names=NA,row.names=TRUE,sep=";",append=TRUE)
# write.table(space, paste(filestem,"Matrix_info/",FW_ID,".csv",sep=''), col.names=F,row.names=F,sep=";",append=TRUE)
# write.table(cann1, paste(filestem,"Matrix_info/",FW_ID,".csv",sep=''), col.names=TRUE,row.names=FALSE,sep=";",append=TRUE)
# write.table(space, paste(filestem,"Matrix_info/",FW_ID,".csv",sep=''), col.names=F,row.names=F,sep=";",append=TRUE)
# write.table(hub, paste(filestem,"Matrix_info/",FW_ID,".csv",sep=''), col.names=TRUE,row.names=FALSE,sep=";",append=TRUE)
#
##

# #define data.frame to store data
# met=data.frame()
# 
# #load workkbook
# wb=loadWorkbook(paste(filestem, name, ".xlsx", sep=''))
# sheets <- getSheets(wb) 
# 
# for(j in 1:length(sheets)){
# temp=IntMat(names(sheets[j]),name)
# met=rbind(met,temp)
# }
# 
# #write dataframe to csv
# write.table(met,paste(filestem,"Metrics_", name,".csv",sep=''),sep=";",col.names=NA,row.names=TRUE)