R/F_pseudo_rectif.R
In polehalieutique/demerstem: This package aims to falicitate Fish stock assessement for west african partners involved in the EU DEMERSTEM project
Defines functions
F_pseudo_rectif
Documented in
F_pseudo_rectif
#' Rectified pseudo cohort
#'
#' \code{F_pseudo_rectif} From terminal F, returns N and F estimate and plots. Trimestrial only works when group = 3 or 4.
#'
#' @param FT Terminal F
#' @param Rinit Initial value for pseudo cohort analysis, equivalent to the last year recruitment
#' @param Mat_C data frame of catch, each column is a different gear
#' @param Mat_R vector of recruitment, if unknown, considered constant
#' @param Mat_E data frame of effort, each column is a different gear
#' @param Mat_M vector of mortality biological parameter
#' @param step_time Time step for the analysis
#' @param Prop_R Proportion of recruitment if quarter analysis (R can be continuons along the year)
#' @param quarter Boolean indicating if quarterly data is used
#' @param ngroup Number of groups
#'
#' @example
#' espece <- "Pseudotolithus senegalensis"
#' Linf <- 44.4
#' K <- 0.35
#' t0 <- -0.24
#' a <- 0.006
#' b <- 3.08
#' M_inf <- 0.2
#' list_age <- seq(0.125,2.875, by = 0.25)
#' age <- length(list_age)
#' Mat_C <- data.frame(indutrial = c(1,1,1,50000,100000,1e6,4e6,5e6,10e6,5e6,3e6,1e6), artisanal = c(1,1,1,30000,100000,2e6,6e6,10e6,4e6,2e6,3e6,1e6))
#' Mat_E <- data.frame(indutrial = c(200,400,200,180,488,174,246,237,342,108,206,150), artisanal = c(439,230,279,279,306,318,524,270,270,398,248,65))
#' Mat_M <- rep(0.2,12)
#' Mat_R <- rep(c(1),4) # Supposed constant here
#' Prop_R <- c(0.3,0.3,0.2,0.2, rep(0,8))
#'step_time <- 4
#'quarter = T
#'FT <- 0.343
#'res <- F_pseudo_rectif(FT = 0.343, age, Mat_C, Mat_R, Mat_E, Mat_M, step_time = 4, Prop_R, quarter = quarter, ngroup = F)
#'res
#' @export
F_pseudo_rectif <- function(FT, age, Mat_C, Mat_R, Mat_E, Mat_M, step_time, Prop_R, quarter, ngroup= F) {
plot <- list()
FT_init <- FT
Mat_q <<- Mat_C #initialisation des capturabilites par les captures
Mat_N<<-rep(NA,age) # effectifs methode de Pope
Mat_N2<<-rep(NA,age) # effectifs methode rectifiee
Mat_Cage <<- apply(Mat_C,1,sum)
# Approximate F from NT and catches
FT_min<- function (x)
{
res <- (Mat_Cage[indice] - (NT * (x/(x + Mat_M[indice]) *
(1 - exp(-x - Mat_M[indice])))))^2
res
}
# returns F and N
R_pseudo_rectif_bis <-function(init) # init <- Rinit
{
Mat_F2 <<- (Mat_C/Mat_Cage)#/step_time # initialisation of F2 as metier ratio (= 1 if only one metier/fleet)
nb_flot <<- ncol(Mat_C) # fleet size
if (quarter == T) {
print(init)
Mat_N2[1] <-round(init * (Prop_R[1])) # N1
assign('Mat_N2', Mat_N2, envir=globalenv())
indice <<-1
# on récupère q1
Mat_F2[indice,] <- Mat_F2[indice,]*optimize(F_a_min,interval=c(0,2))$minimum
assign('Mat_F2', Mat_F2, envir=globalenv())
Mat_q[indice,] <- Mat_F2[indice,]/Mat_E[4,]
assign('Mat_q', Mat_q, envir=globalenv())
#boucle pour apppliquer q1 à q4
if (nb_flot < 2) {
#for (i in 1:4)
#----
Mat_N2[2] <- Mat_N2[1] * exp(-(Mat_q[1, ] * Mat_E[4, ]) - Mat_M[1]) #+ init*Prop_R[1]
assign('Mat_N2', Mat_N2, envir=globalenv())
indice <<- 2
Mat_F2[indice,] <- Mat_F2[indice,]*optimize(F_a_min,interval=c(0,2))$minimum
assign('Mat_F2', Mat_F2, envir=globalenv())
Mat_q[indice,] <- Mat_F2[indice,]/Mat_E[4,]
assign('Mat_q', Mat_q, envir=globalenv())
#---
#---
Mat_N2[3] <- Mat_N2[2] * exp(-(Mat_q[2, ] * Mat_E[3, ]) - Mat_M[2]) + init*Prop_R[2]
assign('Mat_N2', Mat_N2, envir=globalenv())
indice <<- 3
Mat_F2[indice,] <- Mat_F2[indice,]*optimize(F_a_min,interval=c(0,2))$minimum
assign('Mat_F2', Mat_F2, envir=globalenv())
Mat_q[indice,] <- Mat_F2[indice,]/Mat_E[4,]
assign('Mat_q', Mat_q, envir=globalenv())
#---
#---
Mat_N2[4] <- Mat_N2[3] * exp(-(Mat_q[3, ] * Mat_E[2, ]) - Mat_M[3]) + init*Prop_R[3]
assign('Mat_N2', Mat_N2, envir=globalenv())
indice <<- 4
Mat_F2[indice,] <- Mat_F2[indice,]*optimize(F_a_min,interval=c(0,2))$minimum
assign('Mat_F2', Mat_F2, envir=globalenv())
Mat_q[indice,] <- Mat_F2[indice,]/Mat_E[4,]
assign('Mat_q', Mat_q, envir=globalenv())
#---
#---
Mat_N2[5] <- Mat_N2[4] * exp(-(Mat_q[4, ] * Mat_E[1, ]) - Mat_M[4]) + init*Prop_R[4]
assign('Mat_N2', Mat_N2, envir=globalenv())
### Fin 2020
NT2 <- Mat_N2[1] * Mat_R[2]/Mat_R[1] * exp(-(Mat_q[1, ] * Mat_E[8, ]) - Mat_M[1]) + init*Prop_R[1]
NT3 <- NT2 * exp(-(Mat_q[2, ] * Mat_E[7, ]) - Mat_M[2]) + init*Prop_R[2]
NT4 <- NT3 * exp(-(Mat_q[3, ] * Mat_E[6, ]) - Mat_M[3]) + init*Prop_R[3]
NT <- NT4 * exp(-(Mat_q[4, ] * Mat_E[5, ]) - Mat_M[4]) + init*Prop_R[4]
assign('NT', NT, envir=globalenv())
indice <<- 5
Mat_F2[indice,] <- Mat_F2[indice,]*optimize(FT_min,interval=c(0,2))$minimum
assign('Mat_F2', Mat_F2, envir=globalenv())
Mat_q[indice,] <- Mat_F2[indice,]/Mat_E[4,]
assign('Mat_q', Mat_q, envir=globalenv())
#----
#----
Mat_N2[6] <- NT * exp(-(Mat_q[5, ] * Mat_E[4, ]) - Mat_M[5])
assign('Mat_N2', Mat_N2, envir=globalenv())
indice <<- 6
Mat_F2[indice,] <- Mat_F2[indice,]*optimize(F_a_min,interval=c(0,2))$minimum
assign('Mat_F2', Mat_F2, envir=globalenv())
Mat_q[indice,] <- Mat_F2[indice,]/Mat_E[1,]
assign('Mat_q', Mat_q, envir=globalenv())
#----
#----
Mat_N2[7] <- Mat_N2[6] * exp(-(Mat_q[6, ] * Mat_E[3, ]) - Mat_M[6])
assign('Mat_N2', Mat_N2, envir=globalenv())
indice <<- 7
Mat_F2[indice,] <- Mat_F2[indice,]*optimize(F_a_min,interval=c(0,2))$minimum
assign('Mat_F2', Mat_F2, envir=globalenv())
Mat_q[indice,] <- Mat_F2[indice,]/Mat_E[1,]
assign('Mat_q', Mat_q, envir=globalenv())
#----
#----
Mat_N2[8] <- Mat_N2[7] * exp(-(Mat_q[7, ] * Mat_E[2, ]) - Mat_M[7])
assign('Mat_N2', Mat_N2, envir=globalenv())
indice <<- 8
Mat_F2[indice,] <- Mat_F2[indice,]*optimize(F_a_min,interval=c(0,2))$minimum
assign('Mat_F2', Mat_F2, envir=globalenv())
Mat_q[indice,] <- Mat_F2[indice,]/Mat_E[1,]
assign('Mat_q', Mat_q, envir=globalenv())
#----
#----
Mat_N2[9] <- Mat_N2[8] * exp(-(Mat_q[8, ] * Mat_E[1, ]) - Mat_M[8])
assign('Mat_N2', Mat_N2, envir=globalenv())
# fin 2020-1
NT2 <- Mat_N2[1] * Mat_R[3]/Mat_R[1] * exp(-(Mat_q[1, ] * Mat_E[12, ]) - Mat_M[1]) #+ init*Prop_R[2]
NT3 <- NT2 * exp(-(Mat_q[2, ] * Mat_E[11, ]) - Mat_M[2]) + init*Prop_R[2]
NT4 <- NT3 * exp(-(Mat_q[3, ] * Mat_E[10, ]) - Mat_M[3]) + init*Prop_R[3]
NT5 <- NT3 * exp(-(Mat_q[4, ] * Mat_E[9, ]) - Mat_M[4]) + init*Prop_R[4]
NT6 <- NT4 * exp(-(Mat_q[5, ] * Mat_E[9, ]) - Mat_M[4])
NT7 <- Mat_N2[5] * exp(-(Mat_q[6, ] * Mat_E[8, ]) - Mat_M[5])
NT8 <- Mat_N2[6] * exp(-(Mat_q[7, ] * Mat_E[7, ]) - Mat_M[6])
NT <- Mat_N2[7] * exp(-(Mat_q[8, ] * Mat_E[6, ]) - Mat_M[7])
assign('NT', NT, envir=globalenv())
indice <<- 9
Mat_F2[indice,] <- Mat_F2[indice,]*optimize(FT_min,interval=c(0,2))$minimum
assign('Mat_F2', Mat_F2, envir=globalenv())
Mat_q[indice,] <- Mat_F2[indice,]/Mat_E[4,]
assign('Mat_q', Mat_q, envir=globalenv())
Mat_N2[10] <- NT9 * exp(-(Mat_q[8, ] * Mat_E[5, ]) - Mat_M[8])
assign('Mat_N2', Mat_N2, envir=globalenv())
indice <<- 9
Mat_F2[indice,] <- Mat_F2[indice,]*optimize(F_a_min,interval=c(0,2))$minimum
assign('Mat_F2', Mat_F2, envir=globalenv())
Mat_q[indice,] <- Mat_F2[indice,]/Mat_E[1,]
assign('Mat_q', Mat_q, envir=globalenv())
Mat_N2[10] <- Mat_N2[9] * exp(-(Mat_q[9, ] * Mat_E[4, ]) - Mat_M[9])
assign('Mat_N2', Mat_N2, envir=globalenv())
indice <<- 10
Mat_F2[indice,] <- Mat_F2[indice,]*optimize(F_a_min,interval=c(0,2))$minimum
assign('Mat_F2', Mat_F2, envir=globalenv())
Mat_q[indice,] <- Mat_F2[indice,]/Mat_E[1,]
assign('Mat_q', Mat_q, envir=globalenv())
Mat_N2[11] <- Mat_N2[10] * exp(-(Mat_q[10, ] * Mat_E[3, ]) - Mat_M[10])
assign('Mat_N2', Mat_N2, envir=globalenv())
indice <<- 11
Mat_F2[indice,] <- Mat_F2[indice,]*optimize(F_a_min,interval=c(0,2))$minimum
assign('Mat_F2', Mat_F2, envir=globalenv())
Mat_q[indice,] <- Mat_F2[indice,]/Mat_E[1,]
assign('Mat_q', Mat_q, envir=globalenv())
Mat_N2[12] <- Mat_N2[11] * exp(-(Mat_q[11, ] * Mat_E[2, ]) - Mat_M[11])
assign('Mat_N2', Mat_N2, envir=globalenv())
indice <<- 12
Mat_F2[indice,] <- Mat_F2[indice,]*optimize(F_a_min,interval=c(0,2))$minimum
assign('Mat_F2', Mat_F2, envir=globalenv())
Mat_q[indice,] <- Mat_F2[indice,]/Mat_E[1,]
assign('Mat_q', Mat_q, envir=globalenv())
assign('Mat_N2', Mat_N2, envir=globalenv())
}
else {
Mat_N2[2] <- Mat_N2[1] * exp(-apply((Mat_q[1, ] * Mat_E[4, ]), 1, sum) - Mat_M[1]) + init*Prop_R[1]
assign('Mat_N2', Mat_N2, envir=globalenv())
indice <<- 2
Mat_F2[indice,] <- Mat_F2[indice,]*optimize(F_a_min,interval=c(0,2))$minimum
assign('Mat_F2', Mat_F2, envir=globalenv())
Mat_q[indice,] <- Mat_F2[indice,]/Mat_E[4,]
assign('Mat_q', Mat_q, envir=globalenv())
#---
#---
Mat_N2[3] <- Mat_N2[2] * exp(-apply((Mat_q[2, ] * Mat_E[3, ]), 1, sum) - Mat_M[2]) + init*Prop_R[2]
assign('Mat_N2', Mat_N2, envir=globalenv())
indice <<- 3
Mat_F2[indice,] <- Mat_F2[indice,]*optimize(F_a_min,interval=c(0,2))$minimum
assign('Mat_F2', Mat_F2, envir=globalenv())
Mat_q[indice,] <- Mat_F2[indice,]/Mat_E[4,]
assign('Mat_q', Mat_q, envir=globalenv())
#---
#---
Mat_N2[4] <- Mat_N2[3] * exp(-apply((Mat_q[3, ] * Mat_E[2, ]), 1, sum) - Mat_M[3]) + init*Prop_R[3]
assign('Mat_N2', Mat_N2, envir=globalenv())
indice <<- 4
Mat_F2[indice,] <- Mat_F2[indice,]*optimize(F_a_min,interval=c(0,2))$minimum
assign('Mat_F2', Mat_F2, envir=globalenv())
Mat_q[indice,] <- Mat_F2[indice,]/Mat_E[4,]
assign('Mat_q', Mat_q, envir=globalenv())
#---
#---
Mat_N2[5] <- Mat_N2[4] * exp(-apply((Mat_q[4, ] * Mat_E[1, ]), 1, sum) - Mat_M[4]) + init*Prop_R[4]
assign('Mat_N2', Mat_N2, envir=globalenv())
### Fin 2020
NT2 <- Mat_N2[1] * Mat_R[2]/Mat_R[1] * exp(-apply((Mat_q[1, ] * Mat_E[8, ]), 1, sum) - Mat_M[1]) + init*Prop_R[1]
NT3 <- NT2 * exp(-apply((Mat_q[2, ] * Mat_E[7, ]), 1, sum) - Mat_M[2]) + init*Prop_R[2]
NT4 <- NT3 * exp(-apply((Mat_q[3, ] * Mat_E[6, ]), 1, sum) - Mat_M[3]) + init*Prop_R[3]
NT <- NT4 * exp(-apply((Mat_q[4, ] * Mat_E[5, ]), 1, sum) - Mat_M[4]) + init*Prop_R[4]
assign('NT', NT, envir=globalenv())
indice <<- 5
Mat_F2[indice,] <- Mat_F2[indice,]*optimize(FT_min,interval=c(0,2))$minimum
assign('Mat_F2', Mat_F2, envir=globalenv())
Mat_q[indice,] <- Mat_F2[indice,]/Mat_E[4,]
assign('Mat_q', Mat_q, envir=globalenv())
#---
#---
Mat_N2[6] <- NT * exp(-apply((Mat_q[5, ] * Mat_E[4, ]), 1, sum) - Mat_M[5])
assign('Mat_N2', Mat_N2, envir=globalenv())
indice <<- 6
Mat_F2[indice,] <- Mat_F2[indice,]*optimize(F_a_min,interval=c(0,2))$minimum
assign('Mat_F2', Mat_F2, envir=globalenv())
Mat_q[indice,] <- Mat_F2[indice,]/Mat_E[4,]
assign('Mat_q', Mat_q, envir=globalenv())
#---
#---
Mat_N2[7] <- Mat_N2[6] * exp(-apply((Mat_q[6, ] * Mat_E[3, ]), 1, sum) - Mat_M[6])
assign('Mat_N2', Mat_N2, envir=globalenv())
indice <<- 7
Mat_F2[indice,] <- Mat_F2[indice,]*optimize(F_a_min,interval=c(0,2))$minimum
assign('Mat_F2', Mat_F2, envir=globalenv())
Mat_q[indice,] <- Mat_F2[indice,]/Mat_E[4,]
assign('Mat_q', Mat_q, envir=globalenv())
#---
#---
Mat_N2[8] <- Mat_N2[7] * exp(-apply((Mat_q[7, ] * Mat_E[2, ]), 1, sum) - Mat_M[7])
assign('Mat_N2', Mat_N2, envir=globalenv())
indice <<- 8
Mat_F2[indice,] <- Mat_F2[indice,]*optimize(F_a_min,interval=c(0,2))$minimum
assign('Mat_F2', Mat_F2, envir=globalenv())
Mat_q[indice,] <- Mat_F2[indice,]/Mat_E[4,]
assign('Mat_q', Mat_q, envir=globalenv())
#---
#---
Mat_N2[9] <- Mat_N2[8] * exp(-apply((Mat_q[8, ] * Mat_E[1, ]), 1, sum) - Mat_M[8])
assign('Mat_N2', Mat_N2, envir=globalenv())
#---
# fin 2019
#---
NT2 <- Mat_N2[1] * Mat_R[3]/Mat_R[1] * exp(-apply((Mat_q[1, ] * Mat_E[12, ]), 1, sum) - Mat_M[1]) + init*Prop_R[1]
NT3 <- NT2 * exp(-apply((Mat_q[2, ] * Mat_E[11, ]), 1, sum) - Mat_M[2]) + init*Prop_R[2]
NT4 <- NT3 * exp(-apply((Mat_q[3, ] * Mat_E[10, ]), 1, sum) - Mat_M[3]) + init*Prop_R[3]
NT5 <- NT4 * exp(-apply((Mat_q[4, ] * Mat_E[9, ]), 1, sum) - Mat_M[4]) + init*Prop_R[4]
NT6 <- NT5 * exp(-apply((Mat_q[5, ] * Mat_E[8, ]), 1, sum) - Mat_M[5])
NT7 <- NT6 * exp(-apply((Mat_q[6, ] * Mat_E[7, ]), 1, sum) - Mat_M[6])
NT8 <- NT7 * exp(-apply((Mat_q[7, ] * Mat_E[6, ]), 1, sum) - Mat_M[7])
NT <- NT8 * exp(-apply((Mat_q[8, ] * Mat_E[5, ]), 1, sum) - Mat_M[8])
assign('NT', NT, envir=globalenv())
indice <<- 9
Mat_F2[indice,] <- Mat_F2[indice,]*optimize(FT_min,interval=c(0,2))$minimum
assign('Mat_F2', Mat_F2, envir=globalenv())
Mat_q[indice,] <- Mat_F2[indice,]/Mat_E[4,]
assign('Mat_q', Mat_q, envir=globalenv())
#---
#---
Mat_N2[10] <- NT * exp(-apply((Mat_q[9, ] * Mat_E[4, ]), 1, sum) - Mat_M[9])
assign('Mat_N2', Mat_N2, envir=globalenv())
indice <<- 10
Mat_F2[indice,] <- Mat_F2[indice,]*optimize(F_a_min,interval=c(0,2))$minimum
assign('Mat_F2', Mat_F2, envir=globalenv())
Mat_q[indice,] <- Mat_F2[indice,]/Mat_E[4,]
assign('Mat_q', Mat_q, envir=globalenv())
#---
#---
Mat_N2[11] <- Mat_N2[10] * exp(-apply((Mat_q[10, ] * Mat_E[3, ]), 1, sum) - Mat_M[10])
assign('Mat_N2', Mat_N2, envir=globalenv())
indice <<- 11
Mat_F2[indice,] <- Mat_F2[indice,]*optimize(F_a_min,interval=c(0,2))$minimum
assign('Mat_F2', Mat_F2, envir=globalenv())
Mat_q[indice,] <- Mat_F2[indice,]/Mat_E[4,]
assign('Mat_q', Mat_q, envir=globalenv())
#---
#---
Mat_N2[12] <- Mat_N2[11] * exp(-apply((Mat_q[11, ] * Mat_E[2, ]), 1, sum) - Mat_M[11])
assign('Mat_N2', Mat_N2, envir=globalenv())
indice <<- 12
Mat_F2[indice,] <- Mat_F2[indice,]*optimize(F_a_min,interval=c(0,2))$minimum
assign('Mat_F2', Mat_F2, envir=globalenv())
Mat_q[indice,] <- Mat_F2[indice,]/Mat_E[4,]
assign('Mat_q', Mat_q, envir=globalenv())
#---
#---
if(ngroup == 4){
Mat_N2[13] <- Mat_N2[12] * exp(-apply((Mat_q[12, ] * Mat_E[1, ]), 1, sum) - Mat_M[12])
assign('Mat_N2', Mat_N2, envir=globalenv())
NT2 <- Mat_N2[1] * Mat_R[4]/Mat_R[1] * exp(-apply((Mat_q[1, ] * Mat_E[16, ]), 1, sum) - Mat_M[1]) + init*Prop_R[1]
NT3 <- NT2 * exp(-apply((Mat_q[2, ] * Mat_E[15, ]), 1, sum) - Mat_M[2]) + init*Prop_R[2]
NT4 <- NT3 * exp(-apply((Mat_q[3, ] * Mat_E[14, ]), 1, sum) - Mat_M[3]) + init*Prop_R[3]
NT5 <- NT4 * exp(-apply((Mat_q[4, ] * Mat_E[13, ]), 1, sum) - Mat_M[4]) + init*Prop_R[4]
NT6 <- NT5 * exp(-apply((Mat_q[5, ] * Mat_E[12, ]), 1, sum) - Mat_M[5])
NT7 <- NT6 * exp(-apply((Mat_q[6, ] * Mat_E[11, ]), 1, sum) - Mat_M[6])
NT8 <- NT7 * exp(-apply((Mat_q[7, ] * Mat_E[10, ]), 1, sum) - Mat_M[7])
NT9 <- NT8 * exp(-apply((Mat_q[8, ] * Mat_E[9, ]), 1, sum) - Mat_M[8])
NT9 <- NT8 * exp(-apply((Mat_q[9, ] * Mat_E[8, ]), 1, sum) - Mat_M[9])
NT10 <- NT9 * exp(-apply((Mat_q[10, ] * Mat_E[7, ]), 1, sum) - Mat_M[10])
NT11 <- NT10 * exp(-apply((Mat_q[11, ] * Mat_E[6, ]), 1, sum) - Mat_M[11])
NT <- NT11 * exp(-apply((Mat_q[12, ] * Mat_E[5, ]), 1, sum) - Mat_M[12])
assign('NT', NT, envir=globalenv())
indice <<- 13
Mat_F2[indice,] <- Mat_F2[indice,]*optimize(FT_min,interval=c(0,2))$minimum
assign('Mat_F2', Mat_F2, envir=globalenv())
Mat_q[indice,] <- Mat_F2[indice,]/Mat_E[4,]
assign('Mat_q', Mat_q, envir=globalenv())
#---
#---
Mat_N2[14] <- NT * exp(-apply((Mat_q[13, ] * Mat_E[4, ]), 1, sum) - Mat_M[13])
assign('Mat_N2', Mat_N2, envir=globalenv())
indice <<- 14
Mat_F2[indice,] <- Mat_F2[indice,]*optimize(F_a_min,interval=c(0,2))$minimum
assign('Mat_F2', Mat_F2, envir=globalenv())
Mat_q[indice,] <- Mat_F2[indice,]/Mat_E[4,]
assign('Mat_q', Mat_q, envir=globalenv())
#---
#---
Mat_N2[15] <- Mat_N2[14] * exp(-apply((Mat_q[14, ] * Mat_E[3, ]), 1, sum) - Mat_M[14])
assign('Mat_N2', Mat_N2, envir=globalenv())
indice <<- 14
Mat_F2[indice,] <- Mat_F2[indice,]*optimize(F_a_min,interval=c(0,2))$minimum
assign('Mat_F2', Mat_F2, envir=globalenv())
Mat_q[indice,] <- Mat_F2[indice,]/Mat_E[1,]
assign('Mat_q', Mat_q, envir=globalenv())
#---
#---
Mat_N2[15] <- Mat_N2[14] * exp(-apply((Mat_q[14, ] * Mat_E[3, ]), 1, sum) - Mat_M[14])
assign('Mat_N2', Mat_N2, envir=globalenv())
indice <<- 15
Mat_F2[indice,] <- Mat_F2[indice,]*optimize(F_a_min,interval=c(0,2))$minimum
assign('Mat_F2', Mat_F2, envir=globalenv())
Mat_q[indice,] <- Mat_F2[indice,]/Mat_E[1,]
assign('Mat_q', Mat_q, envir=globalenv())
#---
#---
Mat_N2[16] <- Mat_N2[15] * exp(-apply((Mat_q[15, ] * Mat_E[2, ]), 1, sum) - Mat_M[15])
assign('Mat_N2', Mat_N2, envir=globalenv())
indice <<- 16
Mat_F2[indice,] <- Mat_F2[indice,]*optimize(F_a_min,interval=c(0,2))$minimum
assign('Mat_F2', Mat_F2, envir=globalenv())
Mat_q[indice,] <- Mat_F2[indice,]/Mat_E[1,]
assign('Mat_q', Mat_q, envir=globalenv())
assign('Mat_N2', Mat_N2, envir=globalenv())
}
}
}
else{
Mat_F2 <<- (Mat_C/Mat_Cage)#/step_time # initialisation of F2 as metier ratio (= 1 if only one metier/fleet)
indice <<-1
nb_flot <<- ncol(Mat_C) # fleet size
Mat_N2[indice] <- init
assign('Mat_N2', Mat_N2, envir=globalenv())
Mat_F2[indice,] <-Mat_F2[indice,]*optimize(F_a_min,interval=c(0,2))$minimum # Equation de ventilation : optimize (C_th - C_obs) and returns an approximation of F, used to multiply F2
assign('Mat_F2', Mat_F2, envir=globalenv())
# reallocating estimated F by fishery/gear
Mat_q[indice,] <-Mat_F2[indice,]/Mat_E[1,]
assign('Mat_q', Mat_q, envir=globalenv())
Mat_N2[indice+1] <-F_abundance(indice)
assign('Mat_N2', Mat_N2, envir=globalenv())
for(compteur in 2:(age))
{
indice <<- compteur
Mat_F2[indice,] <- Mat_F2[indice,]*optimize(F_a_min,interval=c(0,2))$minimum
assign('Mat_F2', Mat_F2, envir=globalenv())
Mat_q[indice,] <- Mat_F2[indice,]/Mat_E[1,]
assign('Mat_q', Mat_q, envir=globalenv())
if (indice<age) {
Mat_N2[indice+1] <- F_abundance(indice)
assign('Mat_N2', Mat_N2, envir=globalenv())
}
}
}
#on renvoie la somme
if (nb_flot < 2)
{
Mat_F2[age,] # si 1 seul engin
}
else
{
apply(Mat_F2[age,],1,sum) # si multi-engins
}
}
F_pseudo_rectif_global<-function(init) # Return difference between FT and R_pseudo FT and by extracting the minimum using an optimize function, it will return results from an initialization with a FT
{
res<-(FT-R_pseudo_rectif_bis(init))^2
res
}
####
# Repeat for different values of F (0 -> 1)
####
varFT<-seq(0.2/step_time,1/step_time,by=0.2/step_time)
matF3<-as.data.frame(matrix(NA, nrow = age, ncol = length(varFT)))
colnames(matF3) <- varFT
Rinit<-rep(NA,length(varFT))
for (i in 1:length(varFT)) {
VPA_Pope(varFT[i], age, Mat_Cage, Mat_M)
Rinit<-Mat_N[1]
FT <- varFT[i]
optimize(F_pseudo_rectif_global,interval=c(Rinit/10,Rinit*10))$minimum
for (j in 1:age)
if (nb_flot < 2)
{
matF3[j,i]<-Mat_F2[j,]
}
else
{
matF3[j,i]<-apply(Mat_F2[j,],1,sum)
}
}
ylim_max <- max(matF3)+0.2
matF3$Age <- list_age
matF3 <- matF3 %>% pivot_longer(cols = c(1:5), values_to = "Fishing_mortality", names_to = "FT")
matF3$FT <- as.factor(matF3$FT)
matF3$title <- "Variation of F at age for different FT"
####
# TRUE beginning of function, Mat_F is updated with final calculated # this is a cheap trick because I have not enough time
####
# get Rinit estimation to restrict iterations to an intervall [Rinit; 10*Rinit]
FT<- FT_init
Mat_Cage <<- apply(Mat_C,1,sum) # somme des captures de chaque flottille
VPA_Pope(FT, age, Mat_Cage, Mat_M) # Returns Mat_F
Rinit <- Mat_N[1]
# Mat_N<<-rep(NA,age) # effectifs methode de Pope
# Mat_N2<<-rep(NA,age) # effectifs methode rectifiee
# Mat_C <<- Mat_C
# Mat_E <<- Mat_E
# Mat_R <<- Mat_R
# Mat_M <<- Mat_M
indice <<- 1 # age de recrutement
optimize(F_pseudo_rectif_global,interval=c(Rinit/10,Rinit*10))$minimum # Examine for interval on R_init which one optimize FT
Mat_Cage_plot <- as.data.frame(Mat_Cage)[1]
Mat_Cage_plot$Age <- list_age
Mat_Cage_plot$title <- paste0("Catch at age - ", espece)
plot[[1]] <- ggplot(data = Mat_Cage_plot, aes(x = Age, y = Mat_Cage)) + labs(x = "Age", y = "Catch") + geom_line(size = 1.01)+ geom_point() + theme_nice() + facet_grid(~title)
print(plot[[1]])
print(paste0("Matrice des mortalités par pêche F pour FT = ",FT))
Mat_F <- apply(Mat_F2,1,sum)
assign('Mat_F', Mat_F, envir=globalenv())
print(Mat_F)
Mat_F_plot <- as.data.frame(Mat_F)
Mat_F_plot$Age <- list_age
Mat_F_plot$title <- paste0("F at age - ", espece)
Mat_F_plot$F_mean[1] <- mean(c(Mat_F_plot$Mat_F[1], Mat_F_plot$Mat_F[1], Mat_F_plot$Mat_F[2]), na.rm=T)
z <- nrow(Mat_F_plot)
Mat_F_plot$F_mean[z] <- mean(c(Mat_F_plot$Mat_F[z], Mat_F_plot$Mat_F[z], Mat_F_plot$Mat_F[z-1]), na.rm=T)
for (i in 2:(nrow(Mat_F_plot)-1)){
Mat_F_plot$F_mean[i] <- mean(c(Mat_F_plot$Mat_F[i-1], Mat_F_plot$Mat_F[i], Mat_F_plot$Mat_F[i+1]), na.rm=T)
}
plot[[2]] <- ggplot(data = Mat_F_plot) +
labs(x = "Age", y = "Fishing mortality") + geom_line(aes(x = Age, y = Mat_F)) + geom_line(size = 1, aes(x = Age, y = F_mean)) +
geom_point(aes(x = Age, y = F_mean)) + theme_nice() + facet_grid(~title)
# if(length(Mat_C) > 1) { # add Mat_F2 to plot[[2]] with F for each fishery
# Mat_F2_plot <- as.data.frame(Mat_F2)
# Mat_F2_plot$Age <- list_age
# Mat_F2_plot <- Mat_F2_plot %>% pivot_longer(cols=c(1,2),names_to = "Fleet", values_to = "Mat_F")
# plot[[2]] <- plot[[2]] + geom_line(data = Mat_F2_plot, aes(x = Age, y = Mat_F, col = Fleet, linetype = Fleet), size = 1.05)
# }
print(plot[[2]])
print("Mat_N2")
print(Mat_N2)
Mat_N2_plot <- as.data.frame(Mat_N2)
Mat_N2_plot$Age <- list_age
Mat_N2_plot$title <- paste0("Abdundance at age - ", espece)
plot[[3]] <- ggplot(data = Mat_N2_plot, aes(x = Age, y = Mat_N2)) + ylim(c(0, max(Mat_N2))) + labs(x = "Age", y = "Abundance") + geom_line(size = 1) + geom_point() + theme_nice() + facet_grid(~title)
print(plot[[3]])
plot[[4]] <- ggplot(data = matF3, aes(x = Age, y = Fishing_mortality, color = FT)) + geom_line(size = 1.01) + theme_nice() + facet_grid(~title)
print(plot[[4]])
return(plot)
}
polehalieutique/demerstem documentation built on Aug. 4, 2024, 5:12 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions F_pseudo_rectif
Documented in F_pseudo_rectif
#' Rectified pseudo cohort
#'
#' \code{F_pseudo_rectif} From terminal F, returns N and F estimate and plots. Trimestrial only works when group = 3 or 4.
#'
#' @param FT Terminal F
#' @param Rinit Initial value for pseudo cohort analysis, equivalent to the last year recruitment
#' @param Mat_C data frame of catch, each column is a different gear
#' @param Mat_R vector of recruitment, if unknown, considered constant
#' @param Mat_E data frame of effort, each column is a different gear
#' @param Mat_M vector of mortality biological parameter
#' @param step_time Time step for the analysis
#' @param Prop_R Proportion of recruitment if quarter analysis (R can be continuons along the year)
#' @param quarter Boolean indicating if quarterly data is used
#' @param ngroup Number of groups
#'
#' @example
#' espece <- "Pseudotolithus senegalensis"
#' Linf <- 44.4
#' K <- 0.35
#' t0 <- -0.24
#' a <- 0.006
#' b <- 3.08
#' M_inf <- 0.2
#' list_age <- seq(0.125,2.875, by = 0.25)
#' age <- length(list_age)
#' Mat_C <- data.frame(indutrial = c(1,1,1,50000,100000,1e6,4e6,5e6,10e6,5e6,3e6,1e6), artisanal = c(1,1,1,30000,100000,2e6,6e6,10e6,4e6,2e6,3e6,1e6))
#' Mat_E <- data.frame(indutrial = c(200,400,200,180,488,174,246,237,342,108,206,150), artisanal = c(439,230,279,279,306,318,524,270,270,398,248,65))
#' Mat_M <- rep(0.2,12)
#' Mat_R <- rep(c(1),4) # Supposed constant here
#' Prop_R <- c(0.3,0.3,0.2,0.2, rep(0,8))
#'step_time <- 4
#'quarter = T
#'FT <- 0.343
#'res <- F_pseudo_rectif(FT = 0.343, age, Mat_C, Mat_R, Mat_E, Mat_M, step_time = 4, Prop_R, quarter = quarter, ngroup = F)
#'res
#' @export
F_pseudo_rectif <- function(FT, age, Mat_C, Mat_R, Mat_E, Mat_M, step_time, Prop_R, quarter, ngroup= F) {
plot <- list()
FT_init <- FT
Mat_q <<- Mat_C #initialisation des capturabilites par les captures
Mat_N<<-rep(NA,age) # effectifs methode de Pope
Mat_N2<<-rep(NA,age) # effectifs methode rectifiee
Mat_Cage <<- apply(Mat_C,1,sum)
# Approximate F from NT and catches
FT_min<- function (x)
{
res <- (Mat_Cage[indice] - (NT * (x/(x + Mat_M[indice]) *
(1 - exp(-x - Mat_M[indice])))))^2
res
}
# returns F and N
R_pseudo_rectif_bis <-function(init) # init <- Rinit
{
Mat_F2 <<- (Mat_C/Mat_Cage)#/step_time # initialisation of F2 as metier ratio (= 1 if only one metier/fleet)
nb_flot <<- ncol(Mat_C) # fleet size
if (quarter == T) {
print(init)
Mat_N2[1] <-round(init * (Prop_R[1])) # N1
assign('Mat_N2', Mat_N2, envir=globalenv())
indice <<-1
# on récupère q1
Mat_F2[indice,] <- Mat_F2[indice,]*optimize(F_a_min,interval=c(0,2))$minimum
assign('Mat_F2', Mat_F2, envir=globalenv())
Mat_q[indice,] <- Mat_F2[indice,]/Mat_E[4,]
assign('Mat_q', Mat_q, envir=globalenv())
#boucle pour apppliquer q1 à q4
if (nb_flot < 2) {
#for (i in 1:4)
#----
Mat_N2[2] <- Mat_N2[1] * exp(-(Mat_q[1, ] * Mat_E[4, ]) - Mat_M[1]) #+ init*Prop_R[1]
assign('Mat_N2', Mat_N2, envir=globalenv())
indice <<- 2
Mat_F2[indice,] <- Mat_F2[indice,]*optimize(F_a_min,interval=c(0,2))$minimum
assign('Mat_F2', Mat_F2, envir=globalenv())
Mat_q[indice,] <- Mat_F2[indice,]/Mat_E[4,]
assign('Mat_q', Mat_q, envir=globalenv())
#---
#---
Mat_N2[3] <- Mat_N2[2] * exp(-(Mat_q[2, ] * Mat_E[3, ]) - Mat_M[2]) + init*Prop_R[2]
assign('Mat_N2', Mat_N2, envir=globalenv())
indice <<- 3
Mat_F2[indice,] <- Mat_F2[indice,]*optimize(F_a_min,interval=c(0,2))$minimum
assign('Mat_F2', Mat_F2, envir=globalenv())
Mat_q[indice,] <- Mat_F2[indice,]/Mat_E[4,]
assign('Mat_q', Mat_q, envir=globalenv())
#---
#---
Mat_N2[4] <- Mat_N2[3] * exp(-(Mat_q[3, ] * Mat_E[2, ]) - Mat_M[3]) + init*Prop_R[3]
assign('Mat_N2', Mat_N2, envir=globalenv())
indice <<- 4
Mat_F2[indice,] <- Mat_F2[indice,]*optimize(F_a_min,interval=c(0,2))$minimum
assign('Mat_F2', Mat_F2, envir=globalenv())
Mat_q[indice,] <- Mat_F2[indice,]/Mat_E[4,]
assign('Mat_q', Mat_q, envir=globalenv())
#---
#---
Mat_N2[5] <- Mat_N2[4] * exp(-(Mat_q[4, ] * Mat_E[1, ]) - Mat_M[4]) + init*Prop_R[4]
assign('Mat_N2', Mat_N2, envir=globalenv())
### Fin 2020
NT2 <- Mat_N2[1] * Mat_R[2]/Mat_R[1] * exp(-(Mat_q[1, ] * Mat_E[8, ]) - Mat_M[1]) + init*Prop_R[1]
NT3 <- NT2 * exp(-(Mat_q[2, ] * Mat_E[7, ]) - Mat_M[2]) + init*Prop_R[2]
NT4 <- NT3 * exp(-(Mat_q[3, ] * Mat_E[6, ]) - Mat_M[3]) + init*Prop_R[3]
NT <- NT4 * exp(-(Mat_q[4, ] * Mat_E[5, ]) - Mat_M[4]) + init*Prop_R[4]
assign('NT', NT, envir=globalenv())
indice <<- 5
Mat_F2[indice,] <- Mat_F2[indice,]*optimize(FT_min,interval=c(0,2))$minimum
assign('Mat_F2', Mat_F2, envir=globalenv())
Mat_q[indice,] <- Mat_F2[indice,]/Mat_E[4,]
assign('Mat_q', Mat_q, envir=globalenv())
#----
#----
Mat_N2[6] <- NT * exp(-(Mat_q[5, ] * Mat_E[4, ]) - Mat_M[5])
assign('Mat_N2', Mat_N2, envir=globalenv())
indice <<- 6
Mat_F2[indice,] <- Mat_F2[indice,]*optimize(F_a_min,interval=c(0,2))$minimum
assign('Mat_F2', Mat_F2, envir=globalenv())
Mat_q[indice,] <- Mat_F2[indice,]/Mat_E[1,]
assign('Mat_q', Mat_q, envir=globalenv())
#----
#----
Mat_N2[7] <- Mat_N2[6] * exp(-(Mat_q[6, ] * Mat_E[3, ]) - Mat_M[6])
assign('Mat_N2', Mat_N2, envir=globalenv())
indice <<- 7
Mat_F2[indice,] <- Mat_F2[indice,]*optimize(F_a_min,interval=c(0,2))$minimum
assign('Mat_F2', Mat_F2, envir=globalenv())
Mat_q[indice,] <- Mat_F2[indice,]/Mat_E[1,]
assign('Mat_q', Mat_q, envir=globalenv())
#----
#----
Mat_N2[8] <- Mat_N2[7] * exp(-(Mat_q[7, ] * Mat_E[2, ]) - Mat_M[7])
assign('Mat_N2', Mat_N2, envir=globalenv())
indice <<- 8
Mat_F2[indice,] <- Mat_F2[indice,]*optimize(F_a_min,interval=c(0,2))$minimum
assign('Mat_F2', Mat_F2, envir=globalenv())
Mat_q[indice,] <- Mat_F2[indice,]/Mat_E[1,]
assign('Mat_q', Mat_q, envir=globalenv())
#----
#----
Mat_N2[9] <- Mat_N2[8] * exp(-(Mat_q[8, ] * Mat_E[1, ]) - Mat_M[8])
assign('Mat_N2', Mat_N2, envir=globalenv())
# fin 2020-1
NT2 <- Mat_N2[1] * Mat_R[3]/Mat_R[1] * exp(-(Mat_q[1, ] * Mat_E[12, ]) - Mat_M[1]) #+ init*Prop_R[2]
NT3 <- NT2 * exp(-(Mat_q[2, ] * Mat_E[11, ]) - Mat_M[2]) + init*Prop_R[2]
NT4 <- NT3 * exp(-(Mat_q[3, ] * Mat_E[10, ]) - Mat_M[3]) + init*Prop_R[3]
NT5 <- NT3 * exp(-(Mat_q[4, ] * Mat_E[9, ]) - Mat_M[4]) + init*Prop_R[4]
NT6 <- NT4 * exp(-(Mat_q[5, ] * Mat_E[9, ]) - Mat_M[4])
NT7 <- Mat_N2[5] * exp(-(Mat_q[6, ] * Mat_E[8, ]) - Mat_M[5])
NT8 <- Mat_N2[6] * exp(-(Mat_q[7, ] * Mat_E[7, ]) - Mat_M[6])
NT <- Mat_N2[7] * exp(-(Mat_q[8, ] * Mat_E[6, ]) - Mat_M[7])
assign('NT', NT, envir=globalenv())
indice <<- 9
Mat_F2[indice,] <- Mat_F2[indice,]*optimize(FT_min,interval=c(0,2))$minimum
assign('Mat_F2', Mat_F2, envir=globalenv())
Mat_q[indice,] <- Mat_F2[indice,]/Mat_E[4,]
assign('Mat_q', Mat_q, envir=globalenv())
Mat_N2[10] <- NT9 * exp(-(Mat_q[8, ] * Mat_E[5, ]) - Mat_M[8])
assign('Mat_N2', Mat_N2, envir=globalenv())
indice <<- 9
Mat_F2[indice,] <- Mat_F2[indice,]*optimize(F_a_min,interval=c(0,2))$minimum
assign('Mat_F2', Mat_F2, envir=globalenv())
Mat_q[indice,] <- Mat_F2[indice,]/Mat_E[1,]
assign('Mat_q', Mat_q, envir=globalenv())
Mat_N2[10] <- Mat_N2[9] * exp(-(Mat_q[9, ] * Mat_E[4, ]) - Mat_M[9])
assign('Mat_N2', Mat_N2, envir=globalenv())
indice <<- 10
Mat_F2[indice,] <- Mat_F2[indice,]*optimize(F_a_min,interval=c(0,2))$minimum
assign('Mat_F2', Mat_F2, envir=globalenv())
Mat_q[indice,] <- Mat_F2[indice,]/Mat_E[1,]
assign('Mat_q', Mat_q, envir=globalenv())
Mat_N2[11] <- Mat_N2[10] * exp(-(Mat_q[10, ] * Mat_E[3, ]) - Mat_M[10])
assign('Mat_N2', Mat_N2, envir=globalenv())
indice <<- 11
Mat_F2[indice,] <- Mat_F2[indice,]*optimize(F_a_min,interval=c(0,2))$minimum
assign('Mat_F2', Mat_F2, envir=globalenv())
Mat_q[indice,] <- Mat_F2[indice,]/Mat_E[1,]
assign('Mat_q', Mat_q, envir=globalenv())
Mat_N2[12] <- Mat_N2[11] * exp(-(Mat_q[11, ] * Mat_E[2, ]) - Mat_M[11])
assign('Mat_N2', Mat_N2, envir=globalenv())
indice <<- 12
Mat_F2[indice,] <- Mat_F2[indice,]*optimize(F_a_min,interval=c(0,2))$minimum
assign('Mat_F2', Mat_F2, envir=globalenv())
Mat_q[indice,] <- Mat_F2[indice,]/Mat_E[1,]
assign('Mat_q', Mat_q, envir=globalenv())
assign('Mat_N2', Mat_N2, envir=globalenv())
}
else {
Mat_N2[2] <- Mat_N2[1] * exp(-apply((Mat_q[1, ] * Mat_E[4, ]), 1, sum) - Mat_M[1]) + init*Prop_R[1]
assign('Mat_N2', Mat_N2, envir=globalenv())
indice <<- 2
Mat_F2[indice,] <- Mat_F2[indice,]*optimize(F_a_min,interval=c(0,2))$minimum
assign('Mat_F2', Mat_F2, envir=globalenv())
Mat_q[indice,] <- Mat_F2[indice,]/Mat_E[4,]
assign('Mat_q', Mat_q, envir=globalenv())
#---
#---
Mat_N2[3] <- Mat_N2[2] * exp(-apply((Mat_q[2, ] * Mat_E[3, ]), 1, sum) - Mat_M[2]) + init*Prop_R[2]
assign('Mat_N2', Mat_N2, envir=globalenv())
indice <<- 3
Mat_F2[indice,] <- Mat_F2[indice,]*optimize(F_a_min,interval=c(0,2))$minimum
assign('Mat_F2', Mat_F2, envir=globalenv())
Mat_q[indice,] <- Mat_F2[indice,]/Mat_E[4,]
assign('Mat_q', Mat_q, envir=globalenv())
#---
#---
Mat_N2[4] <- Mat_N2[3] * exp(-apply((Mat_q[3, ] * Mat_E[2, ]), 1, sum) - Mat_M[3]) + init*Prop_R[3]
assign('Mat_N2', Mat_N2, envir=globalenv())
indice <<- 4
Mat_F2[indice,] <- Mat_F2[indice,]*optimize(F_a_min,interval=c(0,2))$minimum
assign('Mat_F2', Mat_F2, envir=globalenv())
Mat_q[indice,] <- Mat_F2[indice,]/Mat_E[4,]
assign('Mat_q', Mat_q, envir=globalenv())
#---
#---
Mat_N2[5] <- Mat_N2[4] * exp(-apply((Mat_q[4, ] * Mat_E[1, ]), 1, sum) - Mat_M[4]) + init*Prop_R[4]
assign('Mat_N2', Mat_N2, envir=globalenv())
### Fin 2020
NT2 <- Mat_N2[1] * Mat_R[2]/Mat_R[1] * exp(-apply((Mat_q[1, ] * Mat_E[8, ]), 1, sum) - Mat_M[1]) + init*Prop_R[1]
NT3 <- NT2 * exp(-apply((Mat_q[2, ] * Mat_E[7, ]), 1, sum) - Mat_M[2]) + init*Prop_R[2]
NT4 <- NT3 * exp(-apply((Mat_q[3, ] * Mat_E[6, ]), 1, sum) - Mat_M[3]) + init*Prop_R[3]
NT <- NT4 * exp(-apply((Mat_q[4, ] * Mat_E[5, ]), 1, sum) - Mat_M[4]) + init*Prop_R[4]
assign('NT', NT, envir=globalenv())
indice <<- 5
Mat_F2[indice,] <- Mat_F2[indice,]*optimize(FT_min,interval=c(0,2))$minimum
assign('Mat_F2', Mat_F2, envir=globalenv())
Mat_q[indice,] <- Mat_F2[indice,]/Mat_E[4,]
assign('Mat_q', Mat_q, envir=globalenv())
#---
#---
Mat_N2[6] <- NT * exp(-apply((Mat_q[5, ] * Mat_E[4, ]), 1, sum) - Mat_M[5])
assign('Mat_N2', Mat_N2, envir=globalenv())
indice <<- 6
Mat_F2[indice,] <- Mat_F2[indice,]*optimize(F_a_min,interval=c(0,2))$minimum
assign('Mat_F2', Mat_F2, envir=globalenv())
Mat_q[indice,] <- Mat_F2[indice,]/Mat_E[4,]
assign('Mat_q', Mat_q, envir=globalenv())
#---
#---
Mat_N2[7] <- Mat_N2[6] * exp(-apply((Mat_q[6, ] * Mat_E[3, ]), 1, sum) - Mat_M[6])
assign('Mat_N2', Mat_N2, envir=globalenv())
indice <<- 7
Mat_F2[indice,] <- Mat_F2[indice,]*optimize(F_a_min,interval=c(0,2))$minimum
assign('Mat_F2', Mat_F2, envir=globalenv())
Mat_q[indice,] <- Mat_F2[indice,]/Mat_E[4,]
assign('Mat_q', Mat_q, envir=globalenv())
#---
#---
Mat_N2[8] <- Mat_N2[7] * exp(-apply((Mat_q[7, ] * Mat_E[2, ]), 1, sum) - Mat_M[7])
assign('Mat_N2', Mat_N2, envir=globalenv())
indice <<- 8
Mat_F2[indice,] <- Mat_F2[indice,]*optimize(F_a_min,interval=c(0,2))$minimum
assign('Mat_F2', Mat_F2, envir=globalenv())
Mat_q[indice,] <- Mat_F2[indice,]/Mat_E[4,]
assign('Mat_q', Mat_q, envir=globalenv())
#---
#---
Mat_N2[9] <- Mat_N2[8] * exp(-apply((Mat_q[8, ] * Mat_E[1, ]), 1, sum) - Mat_M[8])
assign('Mat_N2', Mat_N2, envir=globalenv())
#---
# fin 2019
#---
NT2 <- Mat_N2[1] * Mat_R[3]/Mat_R[1] * exp(-apply((Mat_q[1, ] * Mat_E[12, ]), 1, sum) - Mat_M[1]) + init*Prop_R[1]
NT3 <- NT2 * exp(-apply((Mat_q[2, ] * Mat_E[11, ]), 1, sum) - Mat_M[2]) + init*Prop_R[2]
NT4 <- NT3 * exp(-apply((Mat_q[3, ] * Mat_E[10, ]), 1, sum) - Mat_M[3]) + init*Prop_R[3]
NT5 <- NT4 * exp(-apply((Mat_q[4, ] * Mat_E[9, ]), 1, sum) - Mat_M[4]) + init*Prop_R[4]
NT6 <- NT5 * exp(-apply((Mat_q[5, ] * Mat_E[8, ]), 1, sum) - Mat_M[5])
NT7 <- NT6 * exp(-apply((Mat_q[6, ] * Mat_E[7, ]), 1, sum) - Mat_M[6])
NT8 <- NT7 * exp(-apply((Mat_q[7, ] * Mat_E[6, ]), 1, sum) - Mat_M[7])
NT <- NT8 * exp(-apply((Mat_q[8, ] * Mat_E[5, ]), 1, sum) - Mat_M[8])
assign('NT', NT, envir=globalenv())
indice <<- 9
Mat_F2[indice,] <- Mat_F2[indice,]*optimize(FT_min,interval=c(0,2))$minimum
assign('Mat_F2', Mat_F2, envir=globalenv())
Mat_q[indice,] <- Mat_F2[indice,]/Mat_E[4,]
assign('Mat_q', Mat_q, envir=globalenv())
#---
#---
Mat_N2[10] <- NT * exp(-apply((Mat_q[9, ] * Mat_E[4, ]), 1, sum) - Mat_M[9])
assign('Mat_N2', Mat_N2, envir=globalenv())
indice <<- 10
Mat_F2[indice,] <- Mat_F2[indice,]*optimize(F_a_min,interval=c(0,2))$minimum
assign('Mat_F2', Mat_F2, envir=globalenv())
Mat_q[indice,] <- Mat_F2[indice,]/Mat_E[4,]
assign('Mat_q', Mat_q, envir=globalenv())
#---
#---
Mat_N2[11] <- Mat_N2[10] * exp(-apply((Mat_q[10, ] * Mat_E[3, ]), 1, sum) - Mat_M[10])
assign('Mat_N2', Mat_N2, envir=globalenv())
indice <<- 11
Mat_F2[indice,] <- Mat_F2[indice,]*optimize(F_a_min,interval=c(0,2))$minimum
assign('Mat_F2', Mat_F2, envir=globalenv())
Mat_q[indice,] <- Mat_F2[indice,]/Mat_E[4,]
assign('Mat_q', Mat_q, envir=globalenv())
#---
#---
Mat_N2[12] <- Mat_N2[11] * exp(-apply((Mat_q[11, ] * Mat_E[2, ]), 1, sum) - Mat_M[11])
assign('Mat_N2', Mat_N2, envir=globalenv())
indice <<- 12
Mat_F2[indice,] <- Mat_F2[indice,]*optimize(F_a_min,interval=c(0,2))$minimum
assign('Mat_F2', Mat_F2, envir=globalenv())
Mat_q[indice,] <- Mat_F2[indice,]/Mat_E[4,]
assign('Mat_q', Mat_q, envir=globalenv())
#---
#---
if(ngroup == 4){
Mat_N2[13] <- Mat_N2[12] * exp(-apply((Mat_q[12, ] * Mat_E[1, ]), 1, sum) - Mat_M[12])
assign('Mat_N2', Mat_N2, envir=globalenv())
NT2 <- Mat_N2[1] * Mat_R[4]/Mat_R[1] * exp(-apply((Mat_q[1, ] * Mat_E[16, ]), 1, sum) - Mat_M[1]) + init*Prop_R[1]
NT3 <- NT2 * exp(-apply((Mat_q[2, ] * Mat_E[15, ]), 1, sum) - Mat_M[2]) + init*Prop_R[2]
NT4 <- NT3 * exp(-apply((Mat_q[3, ] * Mat_E[14, ]), 1, sum) - Mat_M[3]) + init*Prop_R[3]
NT5 <- NT4 * exp(-apply((Mat_q[4, ] * Mat_E[13, ]), 1, sum) - Mat_M[4]) + init*Prop_R[4]
NT6 <- NT5 * exp(-apply((Mat_q[5, ] * Mat_E[12, ]), 1, sum) - Mat_M[5])
NT7 <- NT6 * exp(-apply((Mat_q[6, ] * Mat_E[11, ]), 1, sum) - Mat_M[6])
NT8 <- NT7 * exp(-apply((Mat_q[7, ] * Mat_E[10, ]), 1, sum) - Mat_M[7])
NT9 <- NT8 * exp(-apply((Mat_q[8, ] * Mat_E[9, ]), 1, sum) - Mat_M[8])
NT9 <- NT8 * exp(-apply((Mat_q[9, ] * Mat_E[8, ]), 1, sum) - Mat_M[9])
NT10 <- NT9 * exp(-apply((Mat_q[10, ] * Mat_E[7, ]), 1, sum) - Mat_M[10])
NT11 <- NT10 * exp(-apply((Mat_q[11, ] * Mat_E[6, ]), 1, sum) - Mat_M[11])
NT <- NT11 * exp(-apply((Mat_q[12, ] * Mat_E[5, ]), 1, sum) - Mat_M[12])
assign('NT', NT, envir=globalenv())
indice <<- 13
Mat_F2[indice,] <- Mat_F2[indice,]*optimize(FT_min,interval=c(0,2))$minimum
assign('Mat_F2', Mat_F2, envir=globalenv())
Mat_q[indice,] <- Mat_F2[indice,]/Mat_E[4,]
assign('Mat_q', Mat_q, envir=globalenv())
#---
#---
Mat_N2[14] <- NT * exp(-apply((Mat_q[13, ] * Mat_E[4, ]), 1, sum) - Mat_M[13])
assign('Mat_N2', Mat_N2, envir=globalenv())
indice <<- 14
Mat_F2[indice,] <- Mat_F2[indice,]*optimize(F_a_min,interval=c(0,2))$minimum
assign('Mat_F2', Mat_F2, envir=globalenv())
Mat_q[indice,] <- Mat_F2[indice,]/Mat_E[4,]
assign('Mat_q', Mat_q, envir=globalenv())
#---
#---
Mat_N2[15] <- Mat_N2[14] * exp(-apply((Mat_q[14, ] * Mat_E[3, ]), 1, sum) - Mat_M[14])
assign('Mat_N2', Mat_N2, envir=globalenv())
indice <<- 14
Mat_F2[indice,] <- Mat_F2[indice,]*optimize(F_a_min,interval=c(0,2))$minimum
assign('Mat_F2', Mat_F2, envir=globalenv())
Mat_q[indice,] <- Mat_F2[indice,]/Mat_E[1,]
assign('Mat_q', Mat_q, envir=globalenv())
#---
#---
Mat_N2[15] <- Mat_N2[14] * exp(-apply((Mat_q[14, ] * Mat_E[3, ]), 1, sum) - Mat_M[14])
assign('Mat_N2', Mat_N2, envir=globalenv())
indice <<- 15
Mat_F2[indice,] <- Mat_F2[indice,]*optimize(F_a_min,interval=c(0,2))$minimum
assign('Mat_F2', Mat_F2, envir=globalenv())
Mat_q[indice,] <- Mat_F2[indice,]/Mat_E[1,]
assign('Mat_q', Mat_q, envir=globalenv())
#---
#---
Mat_N2[16] <- Mat_N2[15] * exp(-apply((Mat_q[15, ] * Mat_E[2, ]), 1, sum) - Mat_M[15])
assign('Mat_N2', Mat_N2, envir=globalenv())
indice <<- 16
Mat_F2[indice,] <- Mat_F2[indice,]*optimize(F_a_min,interval=c(0,2))$minimum
assign('Mat_F2', Mat_F2, envir=globalenv())
Mat_q[indice,] <- Mat_F2[indice,]/Mat_E[1,]
assign('Mat_q', Mat_q, envir=globalenv())
assign('Mat_N2', Mat_N2, envir=globalenv())
}
}
}
else{
Mat_F2 <<- (Mat_C/Mat_Cage)#/step_time # initialisation of F2 as metier ratio (= 1 if only one metier/fleet)
indice <<-1
nb_flot <<- ncol(Mat_C) # fleet size
Mat_N2[indice] <- init
assign('Mat_N2', Mat_N2, envir=globalenv())
Mat_F2[indice,] <-Mat_F2[indice,]*optimize(F_a_min,interval=c(0,2))$minimum # Equation de ventilation : optimize (C_th - C_obs) and returns an approximation of F, used to multiply F2
assign('Mat_F2', Mat_F2, envir=globalenv())
# reallocating estimated F by fishery/gear
Mat_q[indice,] <-Mat_F2[indice,]/Mat_E[1,]
assign('Mat_q', Mat_q, envir=globalenv())
Mat_N2[indice+1] <-F_abundance(indice)
assign('Mat_N2', Mat_N2, envir=globalenv())
for(compteur in 2:(age))
{
indice <<- compteur
Mat_F2[indice,] <- Mat_F2[indice,]*optimize(F_a_min,interval=c(0,2))$minimum
assign('Mat_F2', Mat_F2, envir=globalenv())
Mat_q[indice,] <- Mat_F2[indice,]/Mat_E[1,]
assign('Mat_q', Mat_q, envir=globalenv())
if (indice<age) {
Mat_N2[indice+1] <- F_abundance(indice)
assign('Mat_N2', Mat_N2, envir=globalenv())
}
}
}
#on renvoie la somme
if (nb_flot < 2)
{
Mat_F2[age,] # si 1 seul engin
}
else
{
apply(Mat_F2[age,],1,sum) # si multi-engins
}
}
F_pseudo_rectif_global<-function(init) # Return difference between FT and R_pseudo FT and by extracting the minimum using an optimize function, it will return results from an initialization with a FT
{
res<-(FT-R_pseudo_rectif_bis(init))^2
res
}
####
# Repeat for different values of F (0 -> 1)
####
varFT<-seq(0.2/step_time,1/step_time,by=0.2/step_time)
matF3<-as.data.frame(matrix(NA, nrow = age, ncol = length(varFT)))
colnames(matF3) <- varFT
Rinit<-rep(NA,length(varFT))
for (i in 1:length(varFT)) {
VPA_Pope(varFT[i], age, Mat_Cage, Mat_M)
Rinit<-Mat_N[1]
FT <- varFT[i]
optimize(F_pseudo_rectif_global,interval=c(Rinit/10,Rinit*10))$minimum
for (j in 1:age)
if (nb_flot < 2)
{
matF3[j,i]<-Mat_F2[j,]
}
else
{
matF3[j,i]<-apply(Mat_F2[j,],1,sum)
}
}
ylim_max <- max(matF3)+0.2
matF3$Age <- list_age
matF3 <- matF3 %>% pivot_longer(cols = c(1:5), values_to = "Fishing_mortality", names_to = "FT")
matF3$FT <- as.factor(matF3$FT)
matF3$title <- "Variation of F at age for different FT"
####
# TRUE beginning of function, Mat_F is updated with final calculated # this is a cheap trick because I have not enough time
####
# get Rinit estimation to restrict iterations to an intervall [Rinit; 10*Rinit]
FT<- FT_init
Mat_Cage <<- apply(Mat_C,1,sum) # somme des captures de chaque flottille
VPA_Pope(FT, age, Mat_Cage, Mat_M) # Returns Mat_F
Rinit <- Mat_N[1]
# Mat_N<<-rep(NA,age) # effectifs methode de Pope
# Mat_N2<<-rep(NA,age) # effectifs methode rectifiee
# Mat_C <<- Mat_C
# Mat_E <<- Mat_E
# Mat_R <<- Mat_R
# Mat_M <<- Mat_M
indice <<- 1 # age de recrutement
optimize(F_pseudo_rectif_global,interval=c(Rinit/10,Rinit*10))$minimum # Examine for interval on R_init which one optimize FT
Mat_Cage_plot <- as.data.frame(Mat_Cage)[1]
Mat_Cage_plot$Age <- list_age
Mat_Cage_plot$title <- paste0("Catch at age - ", espece)
plot[[1]] <- ggplot(data = Mat_Cage_plot, aes(x = Age, y = Mat_Cage)) + labs(x = "Age", y = "Catch") + geom_line(size = 1.01)+ geom_point() + theme_nice() + facet_grid(~title)
print(plot[[1]])
print(paste0("Matrice des mortalités par pêche F pour FT = ",FT))
Mat_F <- apply(Mat_F2,1,sum)
assign('Mat_F', Mat_F, envir=globalenv())
print(Mat_F)
Mat_F_plot <- as.data.frame(Mat_F)
Mat_F_plot$Age <- list_age
Mat_F_plot$title <- paste0("F at age - ", espece)
Mat_F_plot$F_mean[1] <- mean(c(Mat_F_plot$Mat_F[1], Mat_F_plot$Mat_F[1], Mat_F_plot$Mat_F[2]), na.rm=T)
z <- nrow(Mat_F_plot)
Mat_F_plot$F_mean[z] <- mean(c(Mat_F_plot$Mat_F[z], Mat_F_plot$Mat_F[z], Mat_F_plot$Mat_F[z-1]), na.rm=T)
for (i in 2:(nrow(Mat_F_plot)-1)){
Mat_F_plot$F_mean[i] <- mean(c(Mat_F_plot$Mat_F[i-1], Mat_F_plot$Mat_F[i], Mat_F_plot$Mat_F[i+1]), na.rm=T)
}
plot[[2]] <- ggplot(data = Mat_F_plot) +
labs(x = "Age", y = "Fishing mortality") + geom_line(aes(x = Age, y = Mat_F)) + geom_line(size = 1, aes(x = Age, y = F_mean)) +
geom_point(aes(x = Age, y = F_mean)) + theme_nice() + facet_grid(~title)
# if(length(Mat_C) > 1) { # add Mat_F2 to plot[[2]] with F for each fishery
# Mat_F2_plot <- as.data.frame(Mat_F2)
# Mat_F2_plot$Age <- list_age
# Mat_F2_plot <- Mat_F2_plot %>% pivot_longer(cols=c(1,2),names_to = "Fleet", values_to = "Mat_F")
# plot[[2]] <- plot[[2]] + geom_line(data = Mat_F2_plot, aes(x = Age, y = Mat_F, col = Fleet, linetype = Fleet), size = 1.05)
# }
print(plot[[2]])
print("Mat_N2")
print(Mat_N2)
Mat_N2_plot <- as.data.frame(Mat_N2)
Mat_N2_plot$Age <- list_age
Mat_N2_plot$title <- paste0("Abdundance at age - ", espece)
plot[[3]] <- ggplot(data = Mat_N2_plot, aes(x = Age, y = Mat_N2)) + ylim(c(0, max(Mat_N2))) + labs(x = "Age", y = "Abundance") + geom_line(size = 1) + geom_point() + theme_nice() + facet_grid(~title)
print(plot[[3]])
plot[[4]] <- ggplot(data = matF3, aes(x = Age, y = Fishing_mortality, color = FT)) + geom_line(size = 1.01) + theme_nice() + facet_grid(~title)
print(plot[[4]])
return(plot)
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.