####################################################################################################################
# Function to compute portfolio effect components
# input SADdat: community matrix (make sure the time index is matching up the columns)
####################################################################################################################
portfolio_fn <- function(SADmat) {
  mat <- SADmat[rowSums(SADmat)!=0,colSums(SADmat)!=0]
  abund.mean <- rowMeans(mat)
  abund.relative <- rowMeans(mat)/sum(rowMeans(mat))
  abund.vcov <- cov(t(mat))
  abund.var <- diag(abund.vcov)
  phi <-  sum(abund.vcov)/(sum(sqrt(abund.var))^2)
  cvpop <- sum(abund.relative*(sqrt(abund.var)/abund.mean)) 
  cvcomm <- sd(colSums(mat))/mean(colSums(mat))
  #Alternative: cvcomm <- sqrt(sum(abund.vcov))/sum(rowMeans(mat)) 
  return(list(cvcomm=cvcomm,cvpop=cvpop,phi=phi))
}
####################################################################################################################
#library(tidyverse)
#load("dsrltmp2024_meta.RData")
#str(gbrfishes) 
# See README FILE for details of data structure 
#ltmp_transectlarge_sprowindex 
# See README FILE: an index to indicate large-transect species in a community matrix

####################################################################################################################
port39reefs_list_fn <- function(gbrfishes, ltmp_transectlarge_sprowindex, ltmp_reefnames, portfolio_fn, simulations = 100){
  # Validate Inputs
  if(length(gbrfishes$ltmp_fishcomm_list) < 39){
    stop("gbrfishes$ltmp_fishcomm_list must contain at least 39 reefs")
  }
  
  if(length(ltmp_reefnames) != 39){
    stop("ltmp_reefnames must have 39 names corresponding to the reefs")
  }
  
  port39reefs.list <- vector("list", simulations) # Re-sampling for specified number of simulations
  
  for(sim in 1:simulations){
    reef_results <- vector("list", 39)
    
    for(i in 1:39){
      mat <- gbrfishes$ltmp_fishcomm_list[[i]] # 39 reefs
      
      # Validate mat structure
      if(!is.matrix(mat)){
        stop(paste("ltmp_fishcomm_list[[", i, "]] is not a matrix", sep=""))
      }
      
      for(j in 1:ncol(mat)){
        # Ensure indices are valid
        if(any(ltmp_transectlarge_sprowindex > nrow(mat) | ltmp_transectlarge_sprowindex < 1)){
          stop("Invalid indices in ltmp_transectlarge_sprowindex")
        }
        
        lamda <- mat[ltmp_transectlarge_sprowindex, j] * 0.2
        
        # Validate lamda
        if(any(is.na(lamda) | lamda < 0)){
          stop("Invalid lambda values for Poisson sampling")
        }
        
        # Perform Poisson sampling
        mat[ltmp_transectlarge_sprowindex, j] <- rpois(length(lamda), lamda)
      }
      
      # Store the resampled matrix locally
      local_fishcomm_list_rpois <- vector("list", 39)
      local_fishcomm_list_rpois[[i]] <- mat
      
      # Compute portfolio effect components
      out <- portfolio_fn(SADmat = local_fishcomm_list_rpois[[i]])
      
      # Ensure 'out' is in the expected format
      if(!is.list(out) && !is.vector(out)){
        stop("portfolio_fn must return a list or a vector")
      }
      
      reef_results[[i]] <- unlist(out)
    }
    
    # Combine reef results into a data frame
    port39reefs <- do.call(rbind, reef_results)
    
    # Validate dimensions before creating the data frame
    if(nrow(port39reefs) != length(ltmp_reefnames)){
      stop("Mismatch between port39reefs rows and ltmp_reefnames length")
    }
    
    port39reefs.list[[sim]] <- data.frame(port39reefs, ReefName = ltmp_reefnames)
  }
  
  return(port39reefs.list)
}