sizespectrum / mizer

Compare fd17403 ... +4 ... 2a2e348


@@ -729,7 +729,7 @@
Loading
729 729
730 730
#' Get density independent rate of egg production
731 731
#' 
732 -
#' Calculates the density-independent rate of total egg production \eqn{R_{p.i}}
732 +
#' Calculates the density-independent rate of total egg production \eqn{R_{di}}{R_di}
733 733
#' (units 1/year) before density dependence, by species. 
734 734
#'
735 735
#' @inherit mizerRDI

@@ -88,13 +88,13 @@
Loading
88 88
        f_mort = r$f_mort, pred_mort = r$pred_mort, t = t, ...)
89 89
    
90 90
    ## Reproduction ----
91 -
    # R_{p,i}
91 +
    # R_di
92 92
    r$rdi <- rates_fns$RDI(
93 93
        params, n = n, n_pp = n_pp, n_other = n_other,
94 94
        e_growth = r$e_growth,
95 95
        mort = r$mort,
96 96
        e_repro = r$e_repro, t = t, ...)
97 -
    # R_i
97 +
    # R_dd
98 98
    r$rdd <- rates_fns$RDD(
99 99
        rdi = r$rdi, species_params = params@species_params, ...)
100 100
    
@@ -643,17 +643,19 @@
Loading
643 643
#' Get density-independent rate of reproduction needed to project standard
644 644
#' mizer model
645 645
#'
646 -
#' Calculates the density-independent rate of total egg production \eqn{R_{p.i}}
647 -
#' (units 1/year) before density dependence, by species. 
646 +
#' Calculates the density-independent rate of total egg production 
647 +
#' \eqn{R_{di}}{R_di} (units 1/year) before density dependence, by species. 
648 648
#' You would not usually call this
649 649
#' function directly but instead use [getRDI()], which then calls this
650 650
#' function unless an alternative function has been registered, see below.
651 651
#' 
652 -
#' This rate is obtained by
653 -
#' taking the total rate at which energy is invested in reproduction, as
654 -
#' calculated by [getERepro()],
655 -
#' multiplying by the reproductive efficiency `erepro` and dividing by the egg
656 -
#' size `w_min`, and by a factor of two to account for the two sexes.
652 +
#' This rate is obtained by taking the per capita rate \eqn{E_r(w)\psi(w)} at
653 +
#' which energy is invested in reproduction, as calculated by [getERepro()],
654 +
#' multiplying it by the number of individuals\eqn{N(w)} and integrating over all sizes
655 +
#' \eqn{w} and then multiplying by the reproductive efficiency \eqn{\epsilon}
656 +
#' and dividing by the egg size `w_min`, and by a factor of two to account for
657 +
#' the two sexes:
658 +
#' \deqn{R_{di} = \frac{\epsilon}{2 w_{min}} \int N(w)  E_r(w) \psi(w) \, dw}{R_di = (\epsilon/(2 w_min)) \int N(w)  E_r(w) \psi(w) dw}
657 659
#' 
658 660
#' Used by [getRDD()] to calculate the actual, density dependent rate.
659 661
#' See [setReproduction()] for more details.
@@ -693,22 +695,26 @@
Loading
693 695
694 696
#' Beverton Holt function to calculate density-dependent reproduction rate
695 697
#'
696 -
#' Takes the density-independent rates \eqn{R_p} of egg production and returns
697 -
#' reduced, density-dependent reproduction rates \eqn{R} given as \deqn{R = R_p
698 -
#' \frac{R_{max}}{R_p + R_{max}}}{R = R_p R_{max}/(R_p + R_{max})} where
699 -
#' \eqn{R_{max}} are the maximum possible reproduction rates that must be
698 +
#' Takes the density-independent rates \eqn{R_{di}}{R_di} of egg production (as
699 +
#' calculated by [getRDI()]) and returns
700 +
#' reduced, density-dependent reproduction rates \eqn{R_{dd}}{R_dd} given as 
701 +
#' \deqn{R_{dd} = R_{di}
702 +
#' \frac{R_{max}}{R_{di} + R_{max}}}{R_dd = R_di R_max/(R_di + R_max)} where
703 +
#' \eqn{R_{max}}{R_max} are the maximum possible reproduction rates that must be
700 704
#' specified in a column in the species parameter dataframe.
705 +
#' (All quantities in the above equation are species-specific but we dropped
706 +
#' the species index for simplicity.)
701 707
#'
702 708
#' This is only one example of a density-dependence. You can write your own
703 709
#' function based on this example, returning different density-dependent
704 710
#' reproduction rates. Two other examples provided are [RickerRDD()]
705 711
#' and [SheperdRDD()]. For more explanation see
706 712
#' [setReproduction()].
707 713
#'
708 -
#' @param rdi Vector of density-independent reproduction rates \eqn{R_p} for all
709 -
#'   species.
714 +
#' @param rdi Vector of density-independent reproduction rates 
715 +
#'   \eqn{R_{di}}{R_di} for all species.
710 716
#' @param species_params A species parameter dataframe. Must contain a column
711 -
#'   R_max holding the maximum reproduction rate \eqn{R_{max}} for each species.
717 +
#'   R_max holding the maximum reproduction rate \eqn{R_{max}}{R_max} for each species.
712 718
#' @param ... Unused
713 719
#'
714 720
#' @return Vector of density-dependent reproduction rates.
@@ -723,12 +729,12 @@
Loading
723 729
724 730
#' Ricker function to calculate density-dependent reproduction rate
725 731
#' 
726 -
#' Takes the density-independent rates \eqn{R_p} of egg production and returns
727 -
#' reduced, density-dependent rates \eqn{R} given as
728 -
#' \deqn{R = R_p \exp{- b R_p}}
732 +
#' Takes the density-independent rates \eqn{R_{di}}{R_di} of egg production and returns
733 +
#' reduced, density-dependent rates \eqn{R_{dd}}{R_dd} given as
734 +
#' \deqn{R_{dd} = R_{di} \exp{- b R_{di}}}{R_dd = R_di \exp{- b R_di}}
729 735
#' 
730 -
#' @param rdi Vector of density-independent reproduction rates \eqn{R_p} for
731 -
#'   all species.
736 +
#' @param rdi Vector of density-independent reproduction rates 
737 +
#'   \eqn{R_{di}}{R_di} for all species.
732 738
#' @param species_params A species parameter dataframe. Must contain a column
733 739
#'   `ricker_b` holding the coefficient b.
734 740
#' @param ... Unused
@@ -745,12 +751,12 @@
Loading
745 751
746 752
#' Sheperd function to calculate density-dependent reproduction rate
747 753
#' 
748 -
#' Takes the density-independent rates \eqn{R_p} of egg production and returns
749 -
#' reduced, density-dependent rates \eqn{R} as
750 -
#' \deqn{R = \frac{R_p}{1+(b\ R_p)^c}}{R = R_p / (1 + (b R_p)^c)}
754 +
#' Takes the density-independent rates \eqn{R_{di}}{R_di} of egg production and returns
755 +
#' reduced, density-dependent rates \eqn{R_{dd}}{R_dd} given as
756 +
#' \deqn{R_{dd} = \frac{R_{di}}{1+(b\ R_{di})^c}}{R_dd = R_di / (1 + (b R_di)^c)}
751 757
#' 
752 -
#' @param rdi Vector of density-independent reproduction rates \eqn{R_p} for
753 -
#'   all species.
758 +
#' @param rdi Vector of density-independent reproduction rates 
759 +
#'   \eqn{R_{di}}{R_di} for all species.
754 760
#' @param species_params A species parameter dataframe. Must contain columns
755 761
#'   `sheperd_b` and `sheperd_c` with the parameters b and c.
756 762
#' @param ... Unused
@@ -769,8 +775,8 @@
Loading
769 775
#' 
770 776
#' Simply returns its `rdi` argument.
771 777
#' 
772 -
#' @param rdi Vector of density-independent reproduction rates \eqn{R_p} for
773 -
#'   all species.
778 +
#' @param rdi Vector of density-independent reproduction rates 
779 +
#'   \eqn{R_{di}}{R_di} for all species.
774 780
#' @param ... Not used.
775 781
#' 
776 782
#' @return Vector of density-dependent reproduction rates.
@@ -784,8 +790,8 @@
Loading
784 790
#'
785 791
#' Simply returns the value from `species_params$constant_reproduction`.
786 792
#'
787 -
#' @param rdi Vector of density-independent reproduction rates \eqn{R_p} for all
788 -
#'   species.
793 +
#' @param rdi Vector of density-independent reproduction rates 
794 +
#'   \eqn{R_{di}}{R_di} for all species.
789 795
#' @param species_params A species parameter dataframe. Must contain a column
790 796
#'   `constant_reproduction`.
791 797
#' @param ... Unused

@@ -48,6 +48,12 @@
Loading
48 48
        msg <- "The later entries of dw_full should be equal to those of dw."
49 49
        errors <- c(errors, msg)
50 50
    }
51 +
    # Check w_min_idx
52 +
    if (any(object@species_params$w_min < object@w[object@w_min_idx]) ||
53 +
        any(object@species_params$w_min > object@w[object@w_min_idx + 1])) {
54 +
        msg <- "The `w_min_idx` should point to the start of the size bin containing the egg size `w_min`."
55 +
        errors <- c(errors, msg)
56 +
    }
51 57
52 58
    # Check the array dimensions are good ----
53 59
    # Bit tricky this one as I don't know of a way to compare lots of vectors 
@@ -426,14 +432,9 @@
Loading
426 432
#' 
427 433
#' @section Changes to species params:
428 434
#' The `species_params` slot of the returned MizerParams object may differ
429 -
#' slightly from the data frame supplied as argument to this function in the
430 -
#' following ways:
431 -
#' \itemize{
432 -
#'   \item Default values are set for \code{w_min, w_inf, alpha, gear, interaction_resource}.
433 -
#'   \item The egg sizes in `w_min` are rounded down to lie on a grid point.
434 -
#' }
435 -
#' Note that the other characteristic sizes of the species, like `w_mat` and
436 -
#' `w_inf`, are not modified to lie on grid points.
435 +
#' slightly from the data frame supplied as argument to this function because
436 +
#' default values are set for 
437 +
#' `w_min, w_inf, alpha, gear, interaction_resource`.
437 438
#' 
438 439
#' @param species_params A data frame of species-specific parameter values.
439 440
#' @param gear_params A data frame with gear-specific parameter values.
@@ -560,15 +561,7 @@
Loading
560 561
    vec1 <- as.numeric(rep(NA, no_w_full))
561 562
    names(vec1) <- signif(w_full, 3)
562 563
    
563 -
    # Round down w_min to lie on grid points and store the indices of these
564 -
    # grid points in w_min_idx
565 -
    w_min_idx <- as.vector(suppressWarnings(
566 -
        tapply(species_params$w_min, 1:no_sp,
567 -
               function(w_min, wx) max(which(wx <= w_min)), wx = w)))
568 -
    # Due to rounding errors this might happen:
569 -
    w_min_idx[w_min_idx == -Inf] <- 1
570 -
    names(w_min_idx) <- species_names
571 -
    species_params$w_min <- w[w_min_idx]
564 +
    w_min_idx <- get_w_min_idx(species_params, w)
572 565
    
573 566
    # Colour and linetype scales ----
574 567
    # for use in plots
@@ -814,6 +807,20 @@
Loading
814 807
        params <- upgradeParams(params)
815 808
        warning("You need to upgrade your MizerParams object with `upgradeParams()`.")
816 809
    }
810 +
    params@w_min_idx <- get_w_min_idx(params@species_params, params@w)
817 811
    validObject(params)
818 812
    params
813 +
}
814 +
815 +
# helper function to calculate w_min_idx slot
816 +
get_w_min_idx <- function(species_params, w) {
817 +
    # Round down w_min to lie on grid points and store the indices of these
818 +
    # grid points in w_min_idx
819 +
    w_min_idx <- as.vector(suppressWarnings(
820 +
        tapply(species_params$w_min, seq_len(nrow(species_params)),
821 +
               function(w_min, wx) max(which(wx <= w_min)), wx = w)))
822 +
    # Due to rounding errors this might happen:
823 +
    w_min_idx[w_min_idx == -Inf] <- 1
824 +
    names(w_min_idx) <- as.character(species_params$species)
825 +
    w_min_idx
819 826
}

@@ -0,0 +1,41 @@
Loading
1 +
#' Set Beverton-Holt density dependence
2 +
#' 
3 +
#' Takes a MizerParams object (with arbitrary density dependence) and sets a
4 +
#' Beverton-Holt density-dependence with a maximum reproduction rate that is a
5 +
#' chosen factor `R_factor` higher than the initial-state reproduction rate. At
6 +
#' the same time it adjusts the reproductive efficiency `erepro` so that the
7 +
#' steady-state abundances do not change. Setting `R_factor = Inf` switches off
8 +
#' all density dependence.
9 +
#' 
10 +
#' @param params A MizerParams object
11 +
#' @param R_factor The factor by which the maximum reproduction rate should be
12 +
#'   higher than the initial-state reproduction rate
13 +
#' 
14 +
#' @return A MizerParams object
15 +
#' @export
16 +
setBevertonHolt <- function(params, R_factor) {
17 +
    assert_that(is(params, "MizerParams"),
18 +
                is.numeric(R_factor),
19 +
                length(R_factor) %in% c(1, nrow(params@species_params)),
20 +
                all(R_factor > 1))
21 +
    
22 +
    rdi <- getRDI(params)
23 +
    rdd <- getRDD(params)
24 +
    params@species_params$R_max <- R_factor * getRDD(params)
25 +
    
26 +
    # erepro needs to be changed to
27 +
    # compensate for using a Beverton Holt relationship
28 +
    # because RDD = (1-1/R_factor) RDI
29 +
    params@species_params$erepro <- 
30 +
        params@species_params$erepro / (1 - 1 / R_factor) *
31 +
        rdd / rdi
32 +
    
33 +
    return(setReproduction(params, RDD = "BevertonHoltRDD"))
34 +
}
35 +
36 +
#' Alias for setBevertonHolt
37 +
#' 
38 +
#' An alias provided for backward compatibility with mizer version <= 2.0.4
39 +
#' @inherit setBevertonHolt
40 +
#' @export
41 +
setRmax <- setBevertonHolt

@@ -58,12 +58,13 @@
Loading
58 58
#' }
59 59
#' 
60 60
#' \subsection{Reproductive efficiency}{
61 -
#' The reproductive efficiency, i.e., the proportion of energy allocated to
61 +
#' The reproductive efficiency \eqn{\epsilon}, i.e., the proportion of energy allocated to
62 62
#' reproduction that results in egg biomass, is set through the `erepro`
63 63
#' column in the species_params data frame. If that is not provided, the default
64 64
#' is set to 1 (which you will want to override). The offspring biomass divided
65 65
#' by the egg biomass gives the rate of egg production, returned by
66 -
#' [getRDI()].
66 +
#' [getRDI()]:
67 +
#' \deqn{R_{di} = \frac{\epsilon}{2 w_{min}} \int N(w)  E_r(w) \psi(w) \, dw}{R_di = (\epsilon/(2 w_min)) \int N(w)  E_r(w) \psi(w) dw}
67 68
#' }
68 69
#' 
69 70
#' \subsection{Density dependence}{

Click to load this diff.
Loading diff...

Learn more Showing 2 files with coverage changes found.

Changes in R/project_methods.R
-1
+1
Loading file...
New file R/setBevertonHolt.R
New
Loading file...
Files Coverage
R 0.01% 91.79%
src/inner_project_loop.cpp 100.00%
Project Totals (30 files) 91.81%
Loading