#' @author Dr. Martin C. Arostegui \email{martin.arostegui@@whoi.edu}
# Code for assessing overlap with scattering layers from light attenuation profiles

#### Create dataframe of all K profiles, depth matched ####
profiles <- list.files("~/Dropbox/WHOI/Albacore/R/Profiles_Kdaytime/")
potential.max.depth <- 500
for (i in 1:length(profiles)){
  profile <- readRDS(paste("~/Dropbox/WHOI/Albacore/R/Profiles_Kdaytime/",profiles[i],sep=""))
  latitude <- profile$meta$latitude
  longitude <- profile$meta$longitude
  date <- profile$meta$date
  instrument_name <- profile$meta$instrument_name
  tag.manufacturer <- as.character(profile$meta$tag.manufacturer)
  min.depth <- min(profile$attenuation$depth)
  max.depth <- max(profile$attenuation$depth)
  mean.depth <- (min.depth + max.depth) / 2
  mae <- profile$meta$mae
  kdz_extract <- profile$attenuation$kdz
  kdz <- rep(NA,potential.max.depth/10)
  kdz[which(sprintf("%sm",seq(10,potential.max.depth,10))==paste(min.depth,"m",sep="")):which(sprintf("%sm",seq(10,potential.max.depth,10))==paste(max.depth,"m",sep=""))] <- 
    kdz_extract
  
  if (i == 1){
    profile_info <- c(instrument_name,tag.manufacturer,date,latitude,longitude,min.depth,max.depth,mean.depth,mae,
                      kdz)
  } else {
    add_info <- c(instrument_name,tag.manufacturer,date,latitude,longitude,min.depth,max.depth,mean.depth,mae,
                  kdz)
    profile_info <- rbind(profile_info,add_info)
  }
}
profile_info <- as.data.frame(profile_info,stringsAsFactors = FALSE)
colnames(profile_info) <- c("instrument_name","tag.manufacturer","date","lat","lon","min","max","mean","mae",
                            sprintf("%sm",seq(10,potential.max.depth,10)))

# Assign date format
profile_info$date <- as.Date(as.numeric(profile_info$date),origin='1970-01-01')

# Assign numeric format to coordinates and profile metrics/kdz
profile_info[,4:ncol(profile_info)] <- sapply(profile_info[,4:ncol(profile_info)],as.numeric)

#Filter profiles with high MAE (> .1 * units/decade)
I <- profile_info$mae > 2 & profile_info$tag.manufacturer=="Wildlife Computers"|profile_info$mae > 3.2 & profile_info$tag.manufacturer=="Lotek"
profile_info <- profile_info[!I,]

#### Align SL Overlap ####
# Long reshape of dataframe to have single column with all kdz values
profile_long <- reshape(profile_info,varying=c(sprintf("%sm",seq(10,potential.max.depth,10))),direction="long",v.names="kdz",idvar=c("instrument_name","date"),timevar="depth")
# Sort by date within fish
profile_long <- dplyr::arrange(profile_long,instrument_name,date)
# Set corresponding depths to kdz values
profile_long$depth <- seq(-10,-potential.max.depth,-10)
# Calculate difference from bin to bin in kdz
profile_long$delta.kdz <- c(diff(profile_long$kdz),NA)

# Define DCM metadata (Cornec et al. 2021 GBC 35:e2020GB006759 - Table S5)
DCM <- data.frame(lat.band = seq(10,50,10), lat.min = seq(5,45,10), lat.max = seq(15,55,10),
                  DCM.depth.mean = c(59,84,94,65,41), DCM.depth.sd = c(16,38,24,24,20))
# 99.7% DCM depth cutoff (potential SL identification must occur below the cutoff)
DCM$depth.cutoff <- DCM$DCM.depth.mean + 3*DCM$DCM.depth.sd

# Create dataframe of all fish profiles and determine scattering layer overlap
sl_prop <- vector(length=25)
for (j in 1:25){
  alb <- subset(profile_long,instrument_name==unique(profile_info$instrument_name)[j])
  # Determine if Scattering Layer was Overlapped
  alb$sl <- NA
  for (i in 1:length(unique(alb$date))){
    ind.date <- which(alb$date==unique(alb$date)[i])
    ind.lat <- which.min(abs(DCM$lat.band-alb$lat[ind.date][1]))
    kdz <- alb$kdz[ind.date]
    depth <- alb$depth[ind.date]
    # Take the difference of the raw data minus the rolling mean (right-aligned) over a 3 depth bin window
    diff.rm <- kdz-zoo::rollmean(kdz,align="right",k=3,fill=NA)
    # Assign NAs at start of df as difference of zero
    diff.rm[is.na(diff.rm)] <- 0
    # If no positive peak differences
    if (all(diff.rm <= 0)) {
      alb$sl[ind.date] <- 0
      # If any positive peak differences
    } else {
      # Find peaks in the raw minus rolling mean
      peaks.rm <- data.frame(pracma::findpeaks(diff.rm)) ; colnames(peaks.rm) <- c("peak.value","peak","start","end")
      # Eliminate any peak differences that are negative
      ind.positive <- which(peaks.rm[,1]>0)
      peaks.rm <- peaks.rm[ind.positive,]
      if (nrow(peaks.rm)==0) {
        alb$sl[ind.date] <- 0
      } else {
        # If profile extends below DCM cutoff
        if (any(depth[peaks.rm[,2]] < -DCM$depth.cutoff[ind.lat])) {
          # Find min kdz below DCM cutoff and its index within the profile
          kdz.minbelowDCM <- min(kdz[which(depth < -DCM$depth.cutoff[ind.lat])],na.rm=T)
          ind.kdz.minbelowDCM <- which(kdz==kdz.minbelowDCM)
          # If min kdz below DCM cutoff is shallower than any kdz below cutoff and their difference is greater than a threshold
          if (any(kdz[ind.kdz.minbelowDCM:length(kdz)] - kdz.minbelowDCM >= 0.005, na.rm=T)) {
            alb$sl[ind.date] <- 1
          } else {
            # Index of peak differences occuring at depths below DCM cutoff
            ind.peaksbelowDCM <- which(depth[peaks.rm[,2]] < -DCM$depth.cutoff[ind.lat])
            # Filter out peaks shallower than DCM cutoff
            peaks.rm <- peaks.rm[ind.peaksbelowDCM,]
            # Index of max peak difference occuring below DCM cutoff
            ind.profile <- peaks.rm[which.max(peaks.rm$peak.value),2]
            # Maximum kdz at/after peak difference minus minimum kdz before [deeper than 100 m]/at peak difference
            diff.sl <- max(kdz[ind.profile:length(kdz)],na.rm=T) - min(kdz[11:ind.profile],na.rm=T)
            # Difference threshold for scattering layer overlap assignment
            if (diff.sl >= 0.005) {
              alb$sl[ind.date] <- 1
            } else {
              alb$sl[ind.date] <- 0
            }
          }
          # If profile does not extend below DCM cutoff
        } else {
          alb$sl[ind.date] <- 0
        }
      }
    }
  }
  if (j == 1){
    profile_SL <- alb
  } else {
    add_info <- alb
    profile_SL <- rbind(profile_SL,add_info)
  }
  sl_prop[j] <- sum(alb$sl)/nrow(alb)
}

# Wide reshape of dataframe to have single row per fish/date with all depth matched kdz values
profile_info <- reshape(profile_SL,direction="wide",v.names="kdz",idvar=c("instrument_name","date"),timevar="depth")
# Remove delta.kdz column, which no longer matches and is not needed
profile_info <- subset(profile_info, select= -delta.kdz)