From 99a2329f40f0f558151c683cd7f07b5ecf5c4e40 Mon Sep 17 00:00:00 2001 From: Gustav Delius Date: Wed, 27 Jan 2021 21:54:37 +0100 Subject: [PATCH] Some improvement to calculation of yield curve, see #6 --- .Rbuildignore | 1 + NAMESPACE | 3 +- R/getYieldVsFcurve.R | 107 ++++++++++++++++++++++++---------------- R/projectToSteady.R | 9 ++++ man/displayFrames.Rd | 3 +- man/getYieldVsF.Rd | 46 +++++++++++++++++ man/getYieldVsFcurve.Rd | 46 ----------------- man/plotYieldVsF.Rd | 50 +++++++++++++++++++ man/projectToSteady.Rd | 4 ++ man/steady.Rd | 53 ++++++++++++++++++++ 10 files changed, 232 insertions(+), 90 deletions(-) create mode 100644 man/getYieldVsF.Rd delete mode 100644 man/getYieldVsFcurve.Rd create mode 100644 man/plotYieldVsF.Rd create mode 100644 man/steady.Rd diff --git a/.Rbuildignore b/.Rbuildignore index 517ca020..138ed169 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -4,3 +4,4 @@ ^_pkgdown\.yml$ ^docs$ ^pkgdown$ +^julia$ diff --git a/NAMESPACE b/NAMESPACE index bada86f5..03d1494c 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -10,9 +10,10 @@ export(distanceMaxRelRDI) export(distanceSSLogN) export(getBiomassFrame) export(getSSBFrame) -export(getYieldVsFcurve) +export(getYieldVsF) export(markBackground) export(newSheldonParams) +export(plotYieldVsF) export(projectToSteady) export(pruneSpecies) export(removeSpecies) diff --git a/R/getYieldVsFcurve.R b/R/getYieldVsFcurve.R index f49a7f61..5b8461b7 100644 --- a/R/getYieldVsFcurve.R +++ b/R/getYieldVsFcurve.R @@ -1,19 +1,18 @@ -#' Calculate a yield-vs-F curve for a species. +#' Calculate a yield-vs-F curve #' -#' Calculates the yield of a species for a range of fishing mortalities -#' for species across all gears at each -#' simulation time step. **Note:** The function just returns the yield -#' at the last time step of the simulation, which might not have converged, -#' or might oscillate. +#' Calculates the yield of a species for a range of fishing mortalities for +#' that species while the fishing mortalities for the other species are held +#' fixed. #' #' @param params An object of class `MizerParams`. -#' @param ixSpecies The species number to make the calculation over -#' @param nSteps The number of steps in F to calculate -#' @param Fmax The maximum fishing mortality -#' @param Frange The range of fishing mortalities to run over -#' @param bPlot Boolean that indicates whether a plot is to be made +#' @param species Name of the target species +#' @param F_range A sequence of fishing mortalities at which to evaluate the +#' yield. If missing, it is set to +#' `seq(from = 0, to = F_max, length.out = no_steps)`. +#' @param F_max The maximum fishing mortality. Used only if `F_range` is missing. +#' @param no_steps The number of steps to use. Only used if `F_range` is missing. #' -#' @return An list with yields and fishing mortalities +#' @return A data frame with columns `F` and `yield`. #' @export #' @family summary functions #' @concept summary_function @@ -21,23 +20,22 @@ #' @examples #' \dontrun{ #' params <- newMultispeciesParams(NS_species_params_gears, inter) -#' y <- getYieldVsFcurve(params,11,bPlot=TRUE) +#' y <- getYieldVsF(params, "Cod", F_max = 1, no_steps = 5) #' } -getYieldVsFcurve <- function(params, - ixSpecies, - nSteps = 10, - Fmax=3, - Frange = seq(0, Fmax, length.out = nSteps), - bPlot = FALSE) { +getYieldVsF <- function(params, + species, + no_steps = 10, + F_max = 3, + F_range = seq(0, F_max, length.out = no_steps)) { # Check parameters params <- validParams(params) - ixSpecies <- as.integer(ixSpecies) - assert_that( - ixSpecies > 0, - ixSpecies <= nrow(params@species_params), - is.numeric(Frange)) + idx_species <- which(params@species_params$species == species) + if (length(idx_species) != 1) { + stop("Invalid species name") + } + assert_that(is.numeric(F_range)) # First make a new gear for that specific species - sp_name <- params@species_params$species[[ixSpecies]] + sp_name <- params@species_params$species[[idx_species]] gp <- gear_params(params) gp$gear <- as.character(gp$gear) gps <- gp$species == sp_name @@ -50,24 +48,49 @@ getYieldVsFcurve <- function(params, gp_extra$catchability <- 1 gp$catchability[gps] <- 0 gear_params(params) <- rbind(gp, gp_extra) - # First run with zero fishing mortality for 150 years - s <- project(params, t_max = 150, c(getInitialEffort(params), tmp = 0), - progress_bar = FALSE) - yield <- 0 * Frange + effort <- getInitialEffort(params) + yield <- 0 * F_range # Loop over fishing mortalities: - for (i in 2:length(Frange)) { - efforts <- c(getInitialEffort(params), tmp = Frange[i]) - s <- project(params, t_max = 40, effort = efforts, - progress_bar = FALSE) - y = getYield(s) - yield[i] = y[dim(y)[[1]], ixSpecies] - } - # Make a plot if requested - if (bPlot) { - plot(Frange, yield, type = "b", lwd = 3, - xlab = "Fishing mortality (1/yr)", ylab = "Yield", - main = species_params(params)$species[[ixSpecies]]) + for (i in seq_along(F_range)) { + params@initial_effort["tmp"] <- F_range[i] + sim <- projectToSteady(params, t_max = 100, + return_sim = TRUE, progress_bar = FALSE) + y <- getYield(sim) + yield[i] <- y[dim(y)[[1]], idx_species] + params <- setInitialValues(params, sim) } - return(list(F = Frange, yield = yield)) + return(data.frame(F = F_range, yield = yield)) +} + +#' Plot a yield-vs-F curve +#' +#' @inherit getYieldVsF +#' +#' @inheritParams getYieldVsF +#' +#' @return A ggplot object +#' @export +#' @family plotting functions +#' @md +#' @examples +#' \dontrun{ +#' params <- newMultispeciesParams(NS_species_params_gears, inter) +#' plotYieldVsFcurve(params, "Cod") +#' } +plotYieldVsF <- function(params, + species, + no_steps = 10, + F_max = 3, + F_range = seq(0, F_max, length.out = no_steps)) { + + curve <- getYieldVsF(params, + species = species, + F_range = F_range) + + ggplot(curve, aes(x = F, y = yield)) + + geom_line() + + xlab("Fishing mortality (1/yr)") + + ylab("Yield") + + ggtitle(species) } diff --git a/R/projectToSteady.R b/R/projectToSteady.R index 5128f6da..7c759f9f 100644 --- a/R/projectToSteady.R +++ b/R/projectToSteady.R @@ -47,6 +47,8 @@ distanceSSLogN <- function(params, current, previous) { #' calculated. #' #' @inheritParams steady +#' @param effort The fishing effort to be used throughout the simulation. +#' This must be a vector or list with one named entry per fishing gear. #' @param distance_func A function that will be called after every `t_per` years #' with both the previous and the new state and that should return a number #' that in some sense measures the distance between the states. By default @@ -56,6 +58,7 @@ distanceSSLogN <- function(params, current, previous) { #' @seealso [distanceSSLogN()], [distanceMaxRelRDI()] #' @export projectToSteady <- function(params, + effort = params@initial_effort, distance_func = distanceSSLogN, t_per = 1.5, t_max = 100, @@ -64,6 +67,7 @@ projectToSteady <- function(params, return_sim = FALSE, progress_bar = TRUE, ...) { params <- validParams(params) + effort <- validEffortVector(effort, params = params) assert_that(t_max >= t_per, tol > 0) if ((t_per < dt) || !isTRUE(all.equal((t_per - round(t_per / dt) * dt), 0))) { @@ -152,6 +156,11 @@ projectToSteady <- function(params, if (return_sim) { sim@params <- params + sel <- 1:i + sim@n <- sim@n[sel, , , drop = FALSE] + sim@n_pp <- sim@n_pp[sel, , drop = FALSE] + sim@n_other <- sim@n_other[sel, , drop = FALSE] + sim@effort <- sim@effort[sel, , drop = FALSE] return(sim) } else { return(params) diff --git a/man/displayFrames.Rd b/man/displayFrames.Rd index 7b05cd4a..667c59f7 100644 --- a/man/displayFrames.Rd +++ b/man/displayFrames.Rd @@ -56,7 +56,8 @@ Other frame functions: \code{\link{getSSBFrame}()} Other plotting functions: -\code{\link{animateSpectra}()} +\code{\link{animateSpectra}()}, +\code{\link{plotYieldVsF}()} } \concept{frame functions} \concept{plotting functions} diff --git a/man/getYieldVsF.Rd b/man/getYieldVsF.Rd new file mode 100644 index 00000000..72aefcaf --- /dev/null +++ b/man/getYieldVsF.Rd @@ -0,0 +1,46 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/getYieldVsFcurve.R +\name{getYieldVsF} +\alias{getYieldVsF} +\title{Calculate a yield-vs-F curve} +\usage{ +getYieldVsF( + params, + species, + no_steps = 10, + F_max = 3, + F_range = seq(0, F_max, length.out = no_steps) +) +} +\arguments{ +\item{params}{An object of class \code{MizerParams}.} + +\item{species}{Name of the target species} + +\item{no_steps}{The number of steps to use. Only used if \code{F_range} is missing.} + +\item{F_max}{The maximum fishing mortality. Used only if \code{F_range} is missing.} + +\item{F_range}{A sequence of fishing mortalities at which to evaluate the +yield. If missing, it is set to +\code{seq(from = 0, to = F_max, length.out = no_steps)}.} +} +\value{ +A data frame with columns \code{F} and \code{yield}. +} +\description{ +Calculates the yield of a species for a range of fishing mortalities for +that species while the fishing mortalities for the other species are held +fixed. +} +\details{ +If +} +\examples{ +\dontrun{ +params <- newMultispeciesParams(NS_species_params_gears, inter) +y <- getYieldVsFcurve(params, "Cod") +} +} +\concept{summary functions} +\concept{summary_function} diff --git a/man/getYieldVsFcurve.Rd b/man/getYieldVsFcurve.Rd deleted file mode 100644 index 070f7779..00000000 --- a/man/getYieldVsFcurve.Rd +++ /dev/null @@ -1,46 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/getYieldVsFcurve.R -\name{getYieldVsFcurve} -\alias{getYieldVsFcurve} -\title{Calculate a yield-vs-F curve for a species.} -\usage{ -getYieldVsFcurve( - params, - ixSpecies, - nSteps = 10, - Fmax = 3, - Frange = seq(0, Fmax, length.out = nSteps), - bPlot = FALSE -) -} -\arguments{ -\item{params}{An object of class \code{MizerParams}.} - -\item{ixSpecies}{The species number to make the calculation over} - -\item{nSteps}{The number of steps in F to calculate} - -\item{Fmax}{The maximum fishing mortality} - -\item{Frange}{The range of fishing mortalities to run over} - -\item{bPlot}{Boolean that indicates whether a plot is to be made} -} -\value{ -An list with yields and fishing mortalities -} -\description{ -Calculates the yield of a species for a range of fishing mortalities -for species across all gears at each -simulation time step. \strong{Note:} The function just returns the yield -at the last time step of the simulation, which might not have converged, -or might oscillate. -} -\examples{ -\dontrun{ -params <- newMultispeciesParams(NS_species_params_gears, inter) -y <- getYieldVsFcurve(params,11,bPlot=TRUE) -} -} -\concept{summary functions} -\concept{summary_function} diff --git a/man/plotYieldVsF.Rd b/man/plotYieldVsF.Rd new file mode 100644 index 00000000..6f7057e8 --- /dev/null +++ b/man/plotYieldVsF.Rd @@ -0,0 +1,50 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/getYieldVsFcurve.R +\name{plotYieldVsF} +\alias{plotYieldVsF} +\title{Plot a yield-vs-F curve} +\usage{ +plotYieldVsF( + params, + species, + no_steps = 10, + F_max = 3, + F_range = seq(0, F_max, length.out = no_steps) +) +} +\arguments{ +\item{params}{An object of class \code{MizerParams}.} + +\item{species}{Name of the target species} + +\item{no_steps}{The number of steps to use. Only used if \code{F_range} is missing.} + +\item{F_max}{The maximum fishing mortality. Used only if \code{F_range} is missing.} + +\item{F_range}{A sequence of fishing mortalities at which to evaluate the +yield. If missing, it is set to +\code{seq(from = 0, to = F_max, length.out = no_steps)}.} +} +\value{ +A ggplot object +} +\description{ +Calculates the yield of a species for a range of fishing mortalities for +that species while the fishing mortalities for the other species are held +fixed. +} +\details{ +If +} +\examples{ +\dontrun{ +params <- newMultispeciesParams(NS_species_params_gears, inter) +y <- plotYieldVsFcurve(params, "Cod") +} +} +\seealso{ +Other plotting functions: +\code{\link{animateSpectra}()}, +\code{\link{displayFrames}()} +} +\concept{plotting functions} diff --git a/man/projectToSteady.Rd b/man/projectToSteady.Rd index 03a42719..2243513b 100644 --- a/man/projectToSteady.Rd +++ b/man/projectToSteady.Rd @@ -6,6 +6,7 @@ \usage{ projectToSteady( params, + effort = params@initial_effort, distance_func = distanceSSLogN, t_per = 1.5, t_max = 100, @@ -19,6 +20,9 @@ projectToSteady( \arguments{ \item{params}{A \linkS4class{MizerParams} object} +\item{effort}{The fishing effort to be used throughout the simulation. +This must be a vector or list with one named entry per fishing gear.} + \item{distance_func}{A function that will be called after every \code{t_per} years with both the previous and the new state and that should return a number that in some sense measures the distance between the states. By default diff --git a/man/steady.Rd b/man/steady.Rd new file mode 100644 index 00000000..74bd4da5 --- /dev/null +++ b/man/steady.Rd @@ -0,0 +1,53 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/projectToSteady.R +\name{steady} +\alias{steady} +\title{Set initial values to a steady state for the model} +\usage{ +steady( + params, + t_max = 100, + t_per = 1.5, + dt = 0.1, + tol = 0.1 * dt, + return_sim = FALSE, + progress_bar = TRUE +) +} +\arguments{ +\item{params}{A \linkS4class{MizerParams} object} + +\item{t_max}{The maximum number of years to run the simulation. Default is 100.} + +\item{t_per}{The simulation is broken up into shorter runs of \code{t_per} years, +after each of which we check for convergence. Default value is 1.5. This +should be chosen as an odd multiple of the timestep \code{dt} in order to be +able to detect period 2 cycles.} + +\item{dt}{The time step to use in \code{project()}.} + +\item{tol}{The simulation stops when the relative change in the egg +production RDI over \code{t_per} years is less than \code{tol} for every species.} + +\item{return_sim}{If TRUE, the function returns the MizerSim object holding +the result of the simulation run. If FALSE (default) the function returns +a MizerParams object with the "initial" slots set to the steady state.} + +\item{progress_bar}{A shiny progress object to implement a progress bar in a +shiny app. Default FALSE.} +} +\description{ +The steady state is found by running the dynamics while keeping reproduction +and other components constant until the size spectra no longer change (or +until time \code{t_max} is reached, if earlier). Then the reproductive efficiencies +are set to the values that give the level of reproduction observed in that +steady state. +} +\examples{ +\dontrun{ +params <- newTraitParams() +species_params(params)$gamma[5] <- 3000 +params <- steady(params) +plotSpectra(params) +} +}