Some improvement to calculation of yield curve, see #6

.Rbuildignore

@@ -4,3 +4,4 @@
 ^_pkgdown\.yml$
 ^docs$
 ^pkgdown$
+^julia$

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)

R/getYieldVsFcurve.R

@@ -1,43 +1,41 @@
-#' 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
 #' @md
 #' @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)
 }

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)