R/extras.R

In imarpe/path: Pelagic Assessment Tools for Hydroacoustics

brody = function(l, Linf, k) {
  # assumes k is dt^-1
  out = Linf - (Linf-l)*exp(-k)
  return(out)
}

.lengthProj = function(l, marcas) {
  
  x = sort(c(l, marcas))
  dif = diff(x)
  pos = findInterval(l, marcas) + 1
  pos = seq(from=pos[1], to=pos[2])
  props = dif[pos]
  props = props/sum(props)
  out = numeric(length(marcas)-1)
  out[pos-1] = props
  return(out)
}

lengthProjMatrix = function(sp, scenario, freq) {
  
  if(!(scenario %in% c("neutro", "favorable", "desfavorable")))
    stop("Wrong scenario.")
  
  dt = 1/freq
  species = getSpeciesInfo(sp)
  bin = species$bin
  
  marcas = .createMarks(species)
  marcas_phi = .createMarks(species, phi=TRUE)
  marcas_inf = marcas - 0.5*bin
  marcas_sup = marcas + 0.5*bin
  
  k    = LHT[[sp]]$G["k", scenario]*dt
  Linf = LHT[[sp]]$G["Linf", scenario]
  
  l_inf = brody(marcas_inf, Linf=Linf, k=k)
  l_sup = brody(marcas_sup, Linf=Linf, k=k)
  
  newMarcas = cbind(l_inf, l_sup)
  A = t(apply(newMarcas, 1, .lengthProj, marcas=marcas_phi))
  
  return(A)
  
}


naturalMortality = function(sp, scenario, freq) {
  
  if(!(scenario %in% c("neutro", "favorable", "desfavorable")))
    stop("Wrong scenario.")  
  dt = 1/freq
  species = getSpeciesInfo(sp)
  bin = species$bin
  
  marcas = .createMarks(species)
  
  M_table    = LHT[[sp]]$M
  
  mPos = findInterval(marcas, M_table$size)
  
  M = M_table[mPos, scenario]*dt
  names(M) = marcas
  
  return(M)
  
}


fishingMortality = function(F, sp) {
  
  species = getSpeciesInfo(sp)
  
  L50 = species$L50
  L75 = species$L75
  
  marcas = .createMarks(species)
  
  selec = .selectivity.log(marcas, L50=L50, L75=L75)
  F = outer(F, selec)

  colnames(F) = marcas
  
  return(F)
  
}


.selectivity.log = function(x, L50, L75, tiny=1e-6) {
  
  s1 = (L50*log(3))/(L75-L50)
  s2 = s1/L50
  selec = 1/(1+exp(s1-(s2*x)))
  selec[selec<tiny] = 0
  names(selec) = x
  return(selec)
  
}

.projectionTimeWeeks = function(start, end){
  # TO_DO: 
  weeks = as.numeric(as.character(difftime(end, start, units = "weeks")))
  weeks = ceiling(weeks)
  return(weeks)
}

.projectPopulation = function (N0, F, scenario, freq, sp, T) {

  if(length(F)==1) F = rep(F, T)/freq
  if(length(F)!=T) stop("Incorrect temporal explotation pattern.")
    
  A = lengthProjMatrix(sp=sp, scenario=scenario, freq=freq)
  M = naturalMortality(sp=sp, scenario=scenario, freq=freq)
  F = fishingMortality(F=F, sp=sp)
  Z = colSums(t(M + t(F)))
  E = max((colSums(F)/Z)*(1-exp(-Z)))
  
  N = matrix(ncol=length(M), nrow=T+1)
  C = matrix(ncol=length(M), nrow=T)
  
  N[1, ] = N0
  
  for(t in seq_len(T)) {
    
    Ft = F[t, ]
    nNew = N[t,]*exp(-(M+Ft))
    C[t, ] = (Ft/(Ft+M))*(N[t,] - nNew)
    N[t+1, ] = as.numeric(nNew %*% A)
    
  }

  species = getSpeciesInfo(sp)
  marcas = .createMarks(species)
  
  #Use "a" and "b" values of the current survey.
  #Change the values in the species.csv (auxiliar folder, 
  #then run the script species.R and build & reload the package.)
  weights = species$a*marcas^species$b
  
  # maturity TO_DO: get values from survey (baseBiologia)
  maturity = .maturity.ojive(sp)
  
  B   = N %*% weights
  Y   = C %*% weights
  BD  = N %*% (maturity*weights)
  Q   = sum(Y)
  BDR = tail(BD, 1)
  
  return(list(N=N, C=C, B=B, Y=Y, BD=BD, Q=Q, BDR=BDR, E=E))
  
}

.maturity.ojive = function(sp) {
  
  species = getSpeciesInfo(sp)
  marcas = .createMarks(species)
  
  #Use "mat1" and "mat2" values of the current survey.
  #Change the values in the species.csv (auxiliar folder, 
  #then run the script species.R and build & reload the package.)
  out = 1/(1+exp(species$mat1+species$mat2*marcas))
  return(out)
  
}

.getScenarios = function(scenario, n) {
  
  scenarios = c("neutro", "desfavorable", "favorable")
  
  if(is.character(scenario)) {
    if(length(scenario)!=1) stop("You can only provide one scenario")
    return(rep(scenario, n))
  }
  
  if(length(scenario)!=length(scenarios)) 
    stop("You must indicate one probability per scenario.")
  
  scenario = sample(x=scenarios, size=n, prob = scenario, replace=TRUE)
  
  return(scenario)
  
}

projectPopulation = function(N, F, scenario, freq, sp, T) {
  
  matrixN    = array(dim=c(T+1, dim(N)))
  matrixC    = array(dim=c(T, dim(N)))
  matrixB    = matrix(nrow=T+1, ncol=ncol(N))
  matrixBD   = matrix(nrow=T+1, ncol=ncol(N))
  matrixY    = matrix(nrow=T, ncol=ncol(N))
  matrixQ    = numeric(ncol(N))
  matrixBDR  = numeric(ncol(N))
  matrixE    = numeric(ncol(N))
  
  scenario = .getScenarios(scenario, ncol(N))
  
  for(i in seq_len(ncol(N))) {
    
    N0 = N[, i]
    sim = .projectPopulation(N=N0, F=F, scenario=scenario[i], freq=freq, sp=sp, T=T)
    
    matrixN[, ,i] = sim$N 
    matrixC[, ,i] = sim$C 
    matrixB[, i]  = sim$B
    matrixBD[, i] = sim$BD
    matrixY[, i]  = sim$Y
    matrixQ[i]    = sim$Q
    matrixBDR[i]  = sim$BDR
    matrixE[i]    = sim$E
    
  }
  
  
  N    = apply(matrixN, c(1,2), median)
  C    = apply(matrixC, c(1,2), median)
  B    = apply(matrixB, 1, median)
  BD   = apply(matrixBD, 1, median)
  Y    = apply(matrixY, 1, median)
  Q    = median(matrixQ)
  BDR  = median(matrixBDR)
  E    = mean(matrixE)
  
  rawData = list(N=matrixN, C=matrixC, B=matrixB, BD=matrixBD, Y=matrixY, 
         Q=matrixQ, BDR=matrixBDR, E=matrixE)
  
  output = list(N=N, C=C, B=B, Y=Y, BD=BD, Q=Q, BDR=BDR, E=E, F=sum(F), 
                raw=rawData)
  
  attr(output, which="sp") = sp
  attr(output, which="freq") = freq
  attr(output, which="T") = T
  
  class(output) = "surveyProj"
  
  return(output)
  
}

summary.surveyProj = function(object, factor=1e-6, BDlim=NULL, ...) {
  
  sp = attr(object, which="sp")
  species = getSpeciesInfo(sp)
  
  if(is.null(BDlim)) BDlim = species$BDlim
  
  matrixBDR = object$raw$BDR
  
  risks = .calculateRisk(x=matrixBDR, ref=BDlim)
  
  out = c(object$F, object$E, factor*object$Q, factor*object$BDR, 
          as.numeric(risks))
  
  riskNames = c("risk", "CTE")
  riskNames = if(length(BDlim)==1) riskNames else 
    paste(riskNames, rep(seq_len(length(BDlim)), each=2), sep=".")
  
  names(out) = c("F", "E", "Q", "BDR", riskNames)
  
  class(out) = "summary.surveyProj"
  return(out)
}

.calculateRisk = function(x, ref) {
  
  output = matrix(nrow=2, ncol=length(ref))
  
  for(iref in seq_along(ref)) {
    output[1, iref] = round(sum(x<ref[iref])/length(x),3)
    output[2, iref] = mean(x[x<ref[iref]])
  }
  
  rownames(output) = c("risk", "CTE")
  colnames(output) = paste("ref", seq_along(ref))
  
  return(output)
}

print.summary.surveyProj = function(x, factor=1e-6, ...) {

  print.default(c(x), ...)
  
}

decisionTable = function(N, f, scenario, freq, sp, T, Fmin=0, 
                         Fmax=0.5, dF=0.01, BDlim=NULL) {

  # TO_DO: save information (f, scenario, etc.)
  # TO_DO: print more information (f, scenario) 
  
  species = getSpeciesInfo(sp)
  if(is.null(BDlim)) BDlim = species$BDlim
  
  Fs = seq(from=Fmin, to=Fmax, by=dF)
  
  output = matrix(nrow=length(Fs), ncol=4 + 2*length(BDlim))
  
  for(i in seq_along(Fs)) {
    F0 = Fs[i]
    sceText = paste(scenario, collapse=", ")
    cat(sprintf("\tSimulation for F=%0.2f and scenario %s\n", F0, sceText))
    sim = projectPopulation(N=N, F=F0*f, scenario=scenario, freq=freq, 
                            sp=sp, T=T)
    output[i, ] = summary(sim, BDlim=BDlim)
  }
  
  output = as.data.frame(output)
  
  riskNames = c("risk", "CTE")
  riskNames = if(length(BDlim)==1) riskNames else 
    paste(riskNames, rep(seq_len(length(BDlim)), each=2), sep=".")
  
  names(output) = c("F", "E", "Q", "BDR", riskNames)
  
  attr(output, which="BDlim") = BDlim
  attr(output, which="Eref")  = species$Eref
  
  class(output) = c("decisionTable", class(output))
  
  return(output)
  
}


plot.decisionTable = function(x, BDlim=NULL, Eref=NULL, factor=1e-6, 
                              xlab = "Tasa de explotación (E)",
                              ylab = "MCTP y Biomasa Desovante (millones t)", 
                              colQ = "blue", colBDR = "red", colBDlim = "grey62", ...) {
  ylim = range(c(x$Q, x$BDR))*c(0.8, 1.2)
  plot(x$E, x$Q, type="l", col=colQ, lwd=2, ylim=ylim,
       xlab = xlab , ylab = ylab)
  
  lines(x$E, x$BDR, col = colBDR, lwd=2)
  if(is.null(BDlim)) BDlim = attributes(x)$BDlim*factor
  #if(is.null(Eref)) Eref = attributes(x)$Eref
  abline(h=BDlim, v=Eref, lty=3, col = colBDlim, lwd = 2)
  legend("topright", c("Biomasa Desovante", "MTCP", "BDlim"), 
         col =c(colBDR, colQ, colBDlim), lty=c(1,1,3), lwd =c(2,2,2), bty = "n")
  return(invisible())
}


plot.decisionTable2 = function(x, BDlim = 5, xlab = "Cuota (mt)", ylab = "Biomasa desovante (mt)", 
                               colBD = "red", colRisk = "blue", colBDlim = "grey62", ...) {
  par(mar = c(4,4,4,4))
  plot(x$Q, x$BDR, type = "l", col = colBD, 
       lwd = 2, ylab = ylab , xlab = xlab,
       ylim = c(0, 1.1*max(x$BDR)))
  par(new = TRUE) 
  plot(x$Q, x$risk*100, type = "l", col = colRisk,
       xaxt = "n", yaxt = "n", xlab = "", ylab = "", lwd = 2) 
  axis(4) 
  par(new=TRUE)
  plot(x$Q, rep(BDlim, nrow(x)), type = "l", col = colBDlim, lwd = 2, lty = 3, 
       ylab = "", xlab = "", ylim = c(0, 1.1*max(x$BDR)))
  mtext("Riesgo (%)", side = 4, line = 3) 
  legend("top", c("Biomasa desovante", "Riesgo", "BDlim"), 
         col = c(colBD, colRisk, colBDlim), lty = c(1,1,3), lwd = c(2,2,2), bty = "n")
}



SEP = function(Tf, T) {
  ft = c(rep(1, Tf), rep(0, T-Tf)) # TO_DO: functions for temporal explotation pattern
  ft = ft/sum(ft)
  return(ft)
}


loadInput = function(file) {
  cat("Loading", file, "...\n")
  load(file)
  objs = ls()
  ind = sapply(mget(objs), inherits, what = "imarpeSurvey.output")
  n = length(which(ind))
  msg1 = "The file %s does not contain valid survey data."
  msg2 = "The file %s contains multiple survey data, only one allowed."
  if(n==0) stop(sprintf(msg1, file))
  if(n>1)  stop(sprintf(msg2, file))
  output = get(x = objs[which(ind)])
  return(output)
}


getAbundanceByLength = function(x){
  n = x$raw$boot$bootMatrixAreaN
  N = apply(n, c(2,3), sum, na.rm = TRUE)
  return(N)
}

assignNames = function(x, labels) {
  names(x) = labels
  return(x)
}


readCatchAtLength = function(file, sp = "anchoveta"){
  if(is.null(file)) return(NULL)
  
  base =  read.csv(file, stringsAsFactors = FALSE)
  colnames(base) = tolower(colnames(base))
  specie = getSpeciesInfo(sp)
  marcas = path:::.createMarks(specie)
  base = as.matrix(base)
  y = base[ ,1]
  
  # Only for anchoveta
  
  if(min(y) != 2 & max(y) != 20){
    x = marcas
    base1 = rbind(matrix(0, nrow = length(seq(from = marcas[1], to = min(base[,1])-specie$bin, by = specie$bin)), 
                         ncol = ncol(base)), base)
    base2 = rbind(base1, matrix(0, nrow = length(seq(from = max(base1[,1])+specie$bin, to = rev(marcas)[1], by = specie$bin)),
                                ncol = ncol(base1)))
    base = base2[ , -1] 
    rownames(base) = NULL} else{
      x = base[, 1]
      base = base
      base = base[, -1]
      rownames(base) = NULL
    }
  
  return(base)
}


plotAbundaceLengthMonth = function(data, sp, monthLabel, xlim, ylim, xlab, ylab, stat,
                                   xtextpos, ytextpos, xjuvpos, yjuvpos, vline, col, ...){
  
  specie = getSpeciesInfo(sp)
  marcas = path:::.createMarks(specie)
  
  data = apply(data, c(1,2), stat)
  data = data.frame(data)
  colnames(data) = marcas
  data = data[nrow(data):1, ]
  data = data/rowSums(data)*100
  juv = round(rowSums(data[1:20])/rowSums(data)*100,1) #only anchoveta 
  data[data==0] = NA
  
  
  plot.new()
  par(mar = c(4,2,1,1))
  plot.window(xlim = xlim, ylim = ylim)
  for(i in seq_len(nrow(data))){
    abline(h = (10*(i-1)), col = "gray", lty = 2)
    lines(marcas, data[i, ] + (10*(i-1)), col = col, lwd = 2)
    text(x = xtextpos, y = ytextpos + 10*(i-1), labels = rev(monthLabel)[i])
    text(x = xjuvpos, y = yjuvpos + 10*(i-1), labels = paste0("juv ="," ", juv[i], "%"))
  }
  abline(v = vline, lty = 2, col = "red")
  axis(1, at = marcas, labels = marcas, cex.axis = 1)
  mtext(text = xlab, side = 1, line = 2.5)
  mtext(text = ylab, side = 2, line = 1)
  box()
}


AbundanceLengthMonth = function(data, stat, ...){
  output = t(apply(data, c(1,2), stat))
  output[output<0] = 0 #Some values are negative (?)
  return(output)
}

#You need "a" and "b" values. 
biomassProyection = function(data, a, b, sp, plot = FALSE, monthLabel, ...){

  if (any(is.null(a), is.null(b))) 
    stop("You need a and b values")
  
  specie = getSpeciesInfo(sp)
  marcas = path:::.createMarks(specie)
  
  biomassMonth = NULL
  biomassQuantil = NULL
  
  for(i in seq_len(dim(data)[1])){
    biomass  = round(median(apply((a*(marcas^b))*data[i, , ], 2, sum)), 2) 
    biomassQ = quantile(apply((a*(marcas^b))*data[i, , ], 2, sum), c(0.025, 0.975))
    biomassMonth   = c(biomassMonth, biomass)
    biomassQuantil = rbind(biomassQuantil, biomassQ)
  }
  
  if(plot == TRUE){
    par(mar = c(4,4,2,2))
    plot(biomassMonth/1e6, type = "l", lwd = 2, col = "black", axes = FALSE, 
         xlab = "Meses", ylab = "Biomasa (millones t)",
         ylim = c(0.9*min(biomassQuantil[ ,1]/1e6), 1.1*max(biomassQuantil[ ,2]/1e6)))
    lines(biomassQuantil[,1]/1e6, col = "red", lwd = 2, lty = 2)
    lines(biomassQuantil[,2]/1e6, col = "red", lwd = 2, lty = 2)
    axis(1, at = seq(1, by = 1, length.out = length(monthLabel)), labels = monthLabel, 
         las = 1, cex.axis = 0.8, line = 0)
    axis(2)
    box()
    legend("topright", c("Biomasa", "Quantiles (2.5 , 97.5)"), col = c("black", "red"),
           lwd = c(2,2), lty = c(1,2), bty = "n")
  }
  
  return(list(Month = biomassMonth, Quantile = biomassQuantil))
}