Surplus production models in place of otherfood

[lentinj]

Instead of ~fixed-abundance otherfood, add action to create a surplus production model:

• Based on total-biomass, i.e. no length/age, __num permanently 1.
• Total biomass varies based on provided model

Some examples for defaults to be found at: https://haddonm.github.io/URMQMF/simple-population-models.html

[bthe]

I think this could be fairly simple to achieve, something like this would go a long way:

st <- g3_stock('test',lengthgroups = 1) %>%
  g3s_livesonareas(1) 

test_spawn <- 
  function(r = g3_parameterized("spm_r", lower = 0.01, upper = 1, value = 0.5, by_stock = by_stock), 
           p = g3_parameterized("spm_p", lower = 0.01, upper = 10, value = 1, by_stock = by_stock), 
           K = g3_parameterized("spm_K", lower = 100, upper = 1e6, value = 1000, by_stock = by_stock))
           {
            list(s = ~sum(stock_ss(stock__wgt) * stock_ss(stock__spawningnum)), 
                r = ~(r*s*(1-(s/K)^p)))
        }

st_actions <- 
list(
  g3a_initialconditions_normalparam(st,
  alpha = 1,
beta = 1), 
  g3a_spawn(st,
    recruitment_f = test_spawn(),
    output_stocks = st,
    run_f = ~cur_step==1,
    mean_f = 1,
    stddev_f = 0.9,
    alpha_f = 1,
    beta_f = 1
  )                                                                                    #  ifmissing = g3_parameterized("t0_init", by_stock = TRUE)))),
  list()
)

where the assumption is that the mean weight of an individual is given by the constant length weight relationship and new individuals are spawned with the fixed l-w relationship.

[lentinj]

It could be even simpler with it's own action, but yes approximately along these lines. test_spawn as above is a bit fast and loose with abundance vs. total biomass, but they're equivalent in this case anyway.

We talked in Gothenburg about fixing abundance to 1 and varying mean(sic) weight as total biomass, but it feels like it be more useful to vary numbers using the discrete logistic model and fixing the mean weight. Then other actions (notably predation, but naturalmortality too, e.g.) would all do the right thing, removing numbers rather than fighting over the 1 giant fish.

The main action would be the same as naturalmortality, but adding to __num rather than scaling it.

You'd have to give a mean-weight value then, but worst-case it could just be 1.

[lentinj]

i.e. this sort of thing (untested, but produces some not-unreasonable R code)

g3a_surplusproductionmodel_logistic <- function(
        r = g3_parameterized("spm_r", lower = 0.01, upper = 1, value = 0.5, by_stock = by_stock), 
        p = g3_parameterized("spm_p", lower = 0.01, upper = 10, value = 1, by_stock = by_stock), 
        K = g3_parameterized("spm_K", lower = 100, upper = 1e6, value = 1000, by_stock = by_stock),
        by_stock = TRUE) {
    s <- ~sum(stock_ss(stock__num))
    ~r *  s * (1-(s/K)^p)
}

g3a_surplusproductionmodel <- function (
        stock,
        spm_num = g3a_surplusproductionmodel_logistic(),
        spm_num_init = g3_parameterized("spm_n0", by_stock = TRUE),
        spm_wgt = 1,
        run_f = TRUE,
        run_at = g3_action_order$initial) {
    out <- new.env(parent = emptyenv())
    action_name <- gadget3:::unique_action_name()

    stock__num <- g3_stock_instance(stock, 0)
    stock__wgt <- g3_stock_instance(stock, 1)

    out[[gadget3:::step_id(run_at, "g3a_surplusproductionmodel", stock, action_name)]] <- g3_step(gadget3:::f_substitute(~{
        debug_label("Surplus production model for ", stock)
        stock_iterate(stock, if (cur_step == 0) {
            stock_ss(stock__num) <- spm_num_init
            stock_ss(stock__wgt) <- spm_wgt
        } else {
            stock_ss(stock__num) <- stock_ss(stock__num) + spm_num
            stock_ss(stock__wgt) <- spm_wgt
        })
    }, list(
        spm_num = spm_num,
        spm_wgt = spm_wgt,
        run_f = run_f )))

    return(as.list(out))
}
g3_to_r(list(g3a_time(1990, 1994), g3a_surplusproductionmodel(g3_stock("moo", 1))))

[bthe]

Yes, that should just about cover it. It may be worthwhile to test a real example with data to see where the wheels fall off and that we get reasonable estimates for the parameters.

[lentinj]

I'm fairly confident putting this together with the other actions as they are won't be a problem. Adding age would require more thought, but I don't think you'd expect that to work anyway.

Presumably if we put this together with a fleet & some likelihood components, we've just replicated a simpler model that we can compare against.