Proposal: Generalise reproduction density dependence balancing for multiple RDD functions

[gustavdelius]

Background & Motivation

Currently in mizer, density dependence in the resource spectrum is handled generically via setResource() and resource_level<-. Different resource dynamics (e.g. resource_semichemostat, resource_logistic) implement their own steady-state balancing functions (balance_resource_semichemostat(), balance_resource_logistic()). This preserves the steady state while tuning resource competition.

In contrast, density dependence in reproduction is split:

• setReproduction() sets the maturity/allocation ogives, erepro, and selects the RDD function name (BevertonHoltRDD, RickerRDD, SheperdRDD, constantRDD, noRDD).
• setBevertonHolt() (and reproduction_level<-) is specific to Beverton-Holt reproduction and forces RDD = "BevertonHoltRDD".

This issue proposes generalizing reproduction density-dependence tuning so that any RDD function can define a steady-state balancing method (e.g. balance_<RDD>()), allowing setReproduction() and reproduction_level<- to work across different stock-recruitment relationships in a unified way.

---

Mathematical Formulation of RDD Balancing

At steady state, the required density-dependent egg production rate is $R_{dd}^* = \text{getRequiredRDD}(\text{params})$.
The density-independent egg production rate is $R_{di} = \frac{\epsilon_{\text{repro}}}{2 w_{\min}} E_R$.

For any density-dependence function $f_{\mathrm{RDD}}(R_{di}, \mathrm{species.params})$, preserving the steady state requires satisfying:
$$f_{\mathrm{RDD}}(R_{di}, \mathrm{species.params}) = R_{dd}^*$$

---

1. Beverton-Holt (BevertonHoltRDD)

$$R_{dd} = \frac{R_{di}}{1 + R_{di} / R_{\max}} = \frac{R_{di} R_{\max}}{R_{di} + R_{\max}}$$

• Parameters: $\epsilon_{\text{repro}}$ (sets $R_{di}$) and $R_{\max}$.
• Reproduction Level: $L_R = \frac{R_{dd}^}{R_{\max}} = 1 - \frac{R_{dd}^}{R_{di}} \in [0, 1)$.
• Balancing:
  • Given $L_R \implies R_{\max} = R_{dd}^* / L_R$ and $R_{di} = R_{dd}^* / (1 - L_R)$.
  • Given $\epsilon_{\text{repro}} \implies R_{\max} = \frac{R_{di} R_{dd}^}{R_{di} - R_{dd}^}$.
  • Given $R_{\max} \implies R_{di} = \frac{R_{\max} R_{dd}^}{R_{\max} - R_{dd}^}$.

---

2. Ricker (RickerRDD)

$$R_{dd} = R_{di} \exp(- b_{\text{ricker}} R_{di})$$

• Parameters: $\epsilon_{\text{repro}}$ (sets $R_{di}$) and $b_{\text{ricker}}$ (ricker_b).
• Maximum Possible Recruitment: $R_{\max} = \max_{R_{di}} R_{dd} = \frac{1}{e \cdot b_{\text{ricker}}}$ (attained at $R_{di} = 1/b_{\text{ricker}}$).
• Reproduction Level: $L_R = \frac{R_{dd}^}{R_{\max}} = e \cdot b_{\text{ricker}} \cdot R_{dd}^ \in (0, 1]$.
• Balancing:
  • Given $\epsilon_{\text{repro}}$ ($R_{di} > R_{dd}^*$):
    $$b_{\text{ricker}} = \frac{\ln(R_{di} / R_{dd}^*)}{R_{di}}$$
    (Unique closed-form solution).
  • Given $L_R$ or $b_{\text{ricker}}$:
    The equation $-b_{\text{ricker}} R_{di} e^{-b_{\text{ricker}} R_{di}} = -b_{\text{ricker}} R_{dd}^*$ has two branches due to overcompensation, solvable via the Lambert $W$ function:
    1. Ascending limb ($R_{di} < 1/b$): $R_{di} = -\frac{1}{b} W_0(-b R_{dd}^*)$ (monotonic / undercompensating regime).
    2. Descending limb ($R_{di} > 1/b$): $R_{di} = -\frac{1}{b} W_{-1}(-b R_{dd}^*)$ (overcompensated / oscillatory regime).

---

3. Shepherd (SheperdRDD)

$$R_{dd} = \frac{R_{di}}{1 + (b_{\text{shep}} R_{di})^c}$$

• Parameters: $\epsilon_{\text{repro}}$ ($R_{di}$), $b_{\text{shep}}$ (sheperd_b), and shape exponent $c$ (sheperd_c).
• Balancing:
  • Given $\epsilon_{\text{repro}}$ and fixed $c$:
    $$b_{\text{shep}} = \frac{1}{R_{di}} \left( \frac{R_{di} - R_{dd}^}{R_{dd}^} \right)^{1/c}$$
    (Unique closed form).
  • When $c = 1$, this collapses identically to Beverton-Holt with $b_{\text{shep}} = 1/R_{\max}$.
  • When $c > 1$, Shepherd exhibits overcompensation (similar to Ricker).

---

4. Constant and No RDD (constantRDD, noRDD)

• noRDD ($R_{dd} = R_{di}$): Forced $R_{di} = R_{dd}^* \implies \epsilon_{\text{repro}} = \epsilon_{\min}$. No free density-dependence parameter; $L_R = 0$.
• constantRDD ($R_{dd} = \text{const}$): Forced $\mathrm{constant.reproduction} = R_{dd}^*$. Reproduction is decoupled from spawner biomass $R_{di}$.

---

Proposed Software Architecture

1. Modular Balancing Dispatch

Similar to balance_resource_<dynamics>(), define:

balance_BevertonHoltRDD(params, erepro = NULL, R_max = NULL, reproduction_level = NULL)
balance_RickerRDD(params, erepro = NULL, ricker_b = NULL, reproduction_level = NULL, limb = c("ascending", "descending"))
balance_SheperdRDD(params, erepro = NULL, sheperd_b = NULL, sheperd_c = NULL, reproduction_level = NULL)

2. Unified setReproduction() and reproduction_level<-

# Setting RDD dynamics and balancing simultaneously:
params <- setReproduction(params, RDD = "RickerRDD", reproduction_level = 0.5)

# Setting reproduction level for the active RDD:
reproduction_level(params) <- 0.5

3. Open Design Considerations

• Handling multi-valued solutions for overcompensating functions (Ricker, Shepherd $c > 1$) by defaulting to the ascending limb unless specified.
• Retaining setBevertonHolt() as a specialized wrapper for backwards compatibility.