## Creation of non-trophic interaction matrix
## Matrix is formated as [from, to] or [source, utilizer], e.g ["Astrophorina", "polychates"]
fillNTImatrix <- function(NTI  = NULL,    # empty non-trophic interaction matrix in correct dimensions
                          Data = NULL){   # excel sheet with data
  for (i in 1:nrow(Data)){
    #fill in Astrophorina non-trophic dependence
    if (Data[i, "NTI_Astrophorina"] == "obligatory"){ 
      NTI <- nti(NTI, "Astrophorina", Data[i, "modellabel"], ity = 1)   
    }else if (Data[i, "NTI_Astrophorina"] == "facultative"){
      NTI <- nti(NTI, "Astrophorina", Data[i, "modellabel"], ity = 2)
    }
    
    #fill in organisms non-trophic dependence
    if (Data[i, "NTI_organisms"] == "obligatory"){ 
      NTI <- nti(NTI, Data[i, "NTI_organismspecification"], Data[i, "modellabel"], ity = 1)    
    }else if (Data[i, "NTI_organisms"] == "facultative"){
      NTI <- nti(NTI, Data[i, "NTI_organismspecification"], Data[i, "modellabel"], ity = 2)
    }
  }
  return(NTI)
}

## Function that adds a single nti to the NTI matrix
nti <- function(NTI  = NULL,
                from = NULL,
                to   = NULL,
                ity  = NULL){ #interaction type (1:3)
  if(from %in% rownames(NTI) & to %in% colnames(NTI)){
    NTI[which(rownames(NTI) == from), which(colnames(NTI) == to)] <- ity
  }else{
    message("Wrong link from, ", from, " to ", to)
  }
  return(NTI)
}

## Creation of the trophic interaction matrix 
## Matrix is formatted as [from, to] or [prey, predator]
fillTImatrix <- function(TI   = NULL,   # empty matrix in correct dimensions
                         Data  = NULL){  # excel sheet with data
  for (i in 1:nrow(Data)){
    from <- Data[i, "feeding"]
    to   <- Data[i, "modellabel"]
    if(from %in% rownames(TI) & to %in% colnames(TI)){
      TI[which(rownames(TI) == from), which(colnames(TI) == to)] <- 1
    }else{
      message("Wrong species names at i = ", i)}
  }
  return(TI)
}


  findNameIndex <- function(names, name){
  	return (match(name, names))
  }

  initiateRemoved <- function(ti, nti, toremove){
  	i_len <- length(toremove)
  	ii <- sort(findNameIndex(rownames(nti), toremove))
  	removed <- data.frame(index = ii,
  		                level = rep(0, i_len),
  		                NTI = rep(0, i_len),
  		                TI = rep(0, i_len),
  		                row.names = rownames(nti)[ii])
  	return(removed)
    }
 	 
  gatherRemovals <- function(removed, matrix, type){
  	i_old <- removed$index
  	ii <- c()
  	for(i in 1:length(i_old)){
  		if(any(matrix[i_old[i], ] == type)) ii <- c(ii, which(matrix[i_old[i], ]  ==  type))
  	}
  	return(ii)
  }
  
  updateObligateNTI <- function(Removed, TI, NTI, level, obl){
  	for(i in 1:length(obl)){
  	  ii <- obl[i]
  		if(!(ii %in% Removed$index)){
  			Removed <- rbind(Removed, c(ii, level, 1, 0))
  		 	rownames(Removed)[nrow(Removed)] <- rownames(NTI)[ii]
  		}
  	}
  	return(Removed)
  }
  
  updateFaculatativeNTI <- function(Removed, TI, NTI, level, fac){
  		for(i in 1:length(fac)){
  			ii <- fac[i]
  			if(!(ii %in% Removed$index)){
  			  jj <- which(NTI[, ii] > 0)  
  			  if(all(ii == jj) | all(jj %in% c(ii, Removed$index[which(Removed$level < level)]))){
  			  	Removed <- rbind(Removed, c(ii, level, 2, 0))
  			    rownames(Removed)[nrow(Removed)] <- rownames(NTI)[ii]
  		   	}
  			}
  		}	
  	  return(Removed)
  	  }
  
  updateTI <- function(Removed, TI, NTI, level, tro){
 		for(i in 1:length(tro)){
  			ii <- tro[i]
  			jj <- which(TI[, ii] > 0)
  			if(ii %in% Removed$index & all(jj %in% c(ii, Removed$index[which(Removed$level < level)]))){
  				  j <- which(ii == Removed$index)
  			    if(level == Removed$level[j]) Removed[j, "TI"] <- 1
  			 }else if (all(jj %in% c(ii, Removed$index[which(Removed$level < level)]))){
  			  	Removed <- rbind(Removed, c(ii, level, 0, 1))
  			  	rownames(Removed)[nrow(Removed)] <- rownames(NTI)[ii]
  		  }
  	}
  	return(Removed)
  	}
  
  updateRemoved <- function(Removed, TI, NTI, level, obl, fac, tro){
  	if(!is.null(obl)){
  		Removed <- updateObligateNTI(Removed, TI, NTI, level, obl)
  	}	
  	if(!is.null(fac)){
  		Removed <- updateFaculatativeNTI(Removed, TI, NTI, level, fac)
  	}
  	if(!is.null(tro)){
  		Removed <- updateTI(Removed, TI, NTI, level, tro)
  	}
    return(Removed)
  }
  
  Remove <- function(TI       = NULL,
	                   NTI      = NULL,
	                   ToRemove = "Astrophorina"){
	  Removed <- initiateRemoved(TI, NTI, ToRemove)
	  checkRemovals <- TRUE
	  level <- 0
	  while(checkRemovals){
      level <- level + 1
      obl <- gatherRemovals(Removed, NTI, 1)
	  	fac <- gatherRemovals(Removed, NTI, 2)
	    tro <- gatherRemovals(Removed, TI, 1)
	    Removed_len <- nrow(Removed)
	    Removed <- updateRemoved(Removed, TI, NTI, level, obl, fac, tro)
	    if (Removed_len == nrow(Removed)) checkRemovals <- FALSE
	  }
	  return(Removed)
  }
 
  ## Calculate food web properties
  FoodWebProperties <- function(fw = NULL){
    S <- nrow(fw)
    L <- sum(fw)
    D <- L/S
    C <- L/(S^2)
    return(data.frame(S = S           ,
                      L = L           ,
                      D = signif(D, 2),
                      C = signif(C, 2))) 	
  }
  
  
  ## Calculate S
  Cal_S <- function(fw = NULL){
    S <- nrow(fw)
  }
  
  Cal_L <- function(fw = NULL){
    L <- sum(fw)
  }
  
  Cal_D <- function(fw = NULL){
    S <- nrow(fw)
    L <- sum(fw)
    D <- L/S
  }
  
  Cal_C <- function(fw = NULL){
    S <- nrow(fw)
    L <- sum(fw)
    C <- L/(S^2)
  }