tests/test-action_renewal-otherfood.R

In gadget-framework/gadget3: Globally-Applicable Area Disaggregated General Ecosystem Toolbox V3

if (!interactive()) options(warn=2, error = function() { sink(stderr()) ; traceback(3) ; q(status = 1) })
library(unittest)

library(gadget3)

other_wgt <- g3_stock('other_wgt', 0)
other_np <- g3_stock('other_np', seq(10, 50, by = 10)) |> g3s_age(3,3)
other_cv <- g3_stock('other_cv', seq(50, 100, by = 10)) |> g3s_age(5,10)
prey_stocks <- list(other_wgt, other_np, other_cv)
pred_a <- g3_stock('pred_a', seq(50, 80, by = 10)) |> g3s_age(0, 10)

pred_a_catch_obs <- expand.grid(
    year = 2000:2005,
    length = seq(0, 100, 10),
    stock = c('other_wgt', 'other_np', 'other_cv'),
    number = 0 )

other_abund_obs <- expand.grid(
    year = 2000:2005,
    length = seq(0, 100, 10),
    stock = c('other_wgt', 'other_np', 'other_cv'),
    number = 0 )

actions <- list(
    g3a_time(2000, 2005, step_lengths = c(6,6), project_years = 0),
    # Parameterized abundance, mean weight

    g3a_otherfood(
       other_wgt,
       g3_parameterized('of_abund', by_year = TRUE, by_stock = TRUE),
       g3_parameterized('of_meanwgt', by_stock = TRUE)),
    g3a_otherfood_normalparam(other_np),
    g3a_otherfood_normalcv(other_cv),
    g3a_age(pred_a),
    g3a_initialconditions(pred_a, ~1e5 + 0 * pred_a__midlen, ~1000),

    g3a_predate(
        pred_a,
        prey_stocks,
        suitabilities = list(
            other_wgt = g3_parameterized('suit', by_stock = TRUE, value = 1),
            other_np = g3_suitability_exponentiall50(),
            other_cv = g3_suitability_exponentiall50() ),
        catchability_f = g3a_predate_catchability_predator(
            temperature = g3_parameterized('temp', value = 0, by_year = TRUE, optimise = FALSE)) ),

    g3l_understocking(prey_stocks),
    g3l_catchdistribution(
        'pred_a_catch',
        pred_a_catch_obs,
        fleets = list(pred_a),
        stocks = prey_stocks,
        g3l_distribution_sumofsquares(),
        nll_breakdown = TRUE,
        report = TRUE ),

    g3l_abundancedistribution(
        'other_abund',
        other_abund_obs,
        fleets = list(),
        stocks = prey_stocks,
        g3l_distribution_sumofsquares(),
        nll_breakdown = TRUE,
        report = TRUE ),

    # NB: Dummy parameter so model will compile in TMB
    ~{nll <- nll + g3_param("x", value = 0)} )
actions <- c(actions, list(
    g3a_report_detail(actions) ))
model_fn <- g3_to_r(actions)
model_cpp <- g3_to_tmb(actions)

ok_group("Default params") ########
attr(model_fn, 'parameter_template') |>
    g3_init_val("other_np.Linf", max(g3_stock_def(other_np, "midlen"))) |>
    g3_init_val("other_cv.Linf", max(g3_stock_def(other_cv, "midlen"))) |>
    g3_init_val("other_np.t0", -1) |>
    g3_init_val("other_cv.t0", -1) |>
    g3_init_val("*.walpha", 0.01) |>
    g3_init_val("*.wbeta", 3) |>
    g3_init_val("other_wgt.of_abund.#", 2000:2005 * 1e3) |>
    g3_init_val("other_np.of_abund.#", 1) |>
    g3_init_val("other_cv.of_abund.#", 1) |>
    g3_init_val("other_wgt.of_meanwgt", 100) |>

    g3_init_val("*.*.alpha", 0.07, lower = 0.01, upper = 0.2) |>
    g3_init_val("other_np.*.l50", 45) |>
    g3_init_val("other_cv.*.l50", 55) |>

    identity() -> params
result <- model_fn(params)

r <- attributes(result)

## other_wgt catch
ok(ut_cmp_equal(
    sum(r$cdist_sumofsquares_pred_a_catch_model__num[-1, stock="other_wgt",]),
    0), "cdist_sumofsquares_pred_a_catch_model__num: No other_wgt caught outside it's single lengthgroup")
ok(ut_cmp_equal(
    r$cdist_sumofsquares_pred_a_catch_model__num[1,stock="other_wgt",],
    c("2000" = 11286.74, "2001" = 11290.94, "2002" = 11295.13, "2003" = 11299.33, "2004" = 11303.52, "2005" = 11307.71),
    tolerance = 1e-6 ), "cdist_sumofsquares_pred_a_catch_model__num: other_wgt catch steadily increasing")

## other_np catch
ok(ut_cmp_equal(
    r$cdist_sumofsquares_pred_a_catch_model__num[,stock="other_np",time="2000"],
    c("0:10" = 0, "10:20" = 0.0016687, "20:30" = 0.090596, "30:40" = 1.6739, "40:50" = 10.221,
    "50:60" = 20.362, "60:70" = 0, "70:80" = 0, "80:90" = 0, "90:100" = 0, "100:Inf" = 0),
    tolerance = 1e-3 ), "cdist_sumofsquares_pred_a_catch_model__num: other_np catch has selectivity curve")
for (yr in 2000:2005) ok(ut_cmp_equal(
    r$cdist_sumofsquares_pred_a_catch_model__num[,stock="other_np",time=as.character(yr)],
    r$cdist_sumofsquares_pred_a_catch_model__num[,stock="other_np",time="2000"],
    tolerance = 1e-3), paste0("cdist_sumofsquares_pred_a_catch_model__num: other_np catch ~stable across year ", yr))

## other_cv catch
ok(ut_cmp_equal(
    r$cdist_sumofsquares_pred_a_catch_model__num[,stock="other_cv",time="2000"],
    c("0:10" = 0, "10:20" = 0, "20:30" = 0, "30:40" = 0, "40:50" = 0, "50:60" = 0.0011257,
    "60:70" = 0.08904, "70:80" = 2.5512, "80:90" = 27.268, "90:100" = 111.97, "100:Inf" = 180.45),
    tolerance = 1e-3 ), "cdist_sumofsquares_pred_a_catch_model__num: other_cv catch has selectivity curve")
for (yr in 2000:2005) ok(ut_cmp_equal(
    r$cdist_sumofsquares_pred_a_catch_model__num[,stock="other_cv",time=as.character(yr)],
    r$cdist_sumofsquares_pred_a_catch_model__num[,stock="other_cv",time="2000"],
    tolerance = 1e-3), paste0("cdist_sumofsquares_pred_a_catch_model__num: other_cv catch ~stable across year ", yr))

## other_wgt abundance
ok(gadget3:::ut_cmp_array(r$dstart_other_wgt__num, "
   length    time Freq
1   0:Inf 2000-01 2000000
2   0:Inf 2000-02 2000000
3   0:Inf 2001-01 2001000
4   0:Inf 2001-02 2001000
5   0:Inf 2002-01 2002000
6   0:Inf 2002-02 2002000
7   0:Inf 2003-01 2003000
8   0:Inf 2003-02 2003000
9   0:Inf 2004-01 2004000
10  0:Inf 2004-02 2004000
11  0:Inf 2005-01 2005000
12  0:Inf 2005-02 2005000
"), "dstart_other_wgt__num: current year * 1000, not affected by catch")
ok(gadget3:::ut_cmp_array(r$dstart_other_wgt__wgt, "
   length    time Freq
1   0:Inf 2000-01  100
2   0:Inf 2000-02  100
3   0:Inf 2001-01  100
4   0:Inf 2001-02  100
5   0:Inf 2002-01  100
6   0:Inf 2002-02  100
7   0:Inf 2003-01  100
8   0:Inf 2003-02  100
9   0:Inf 2004-01  100
10  0:Inf 2004-02  100
11  0:Inf 2005-01  100
12  0:Inf 2005-02  100
"), "dstart_other_wgt__wgt: All 100, not affected by catch")

## other_np abundance
ok(gadget3:::ut_cmp_array(r$dstart_other_np__num[,,"2000-01", drop = FALSE], "
  length  age    time        Freq
1  10:20 age3 2000-01    2.710355
2  20:30 age3 2000-01   81.153587
3  30:40 age3 2000-01  893.912126
4  40:50 age3 2000-01 3622.325360
5 50:Inf age3 2000-01 5399.898572
"), "dstart_other_np__num: vonB applied")
for (t in dimnames(r$dstart_other_np__num)$time) ok(ut_cmp_equal(
    r$dstart_other_np__num[,,"2000-01"],
    r$dstart_other_np__num[,,t],
    ignore = NULL ), paste0("dstart_other_np__num: Stable across timestep ", t, ", not affected by catch"))
for (t in dimnames(r$dstart_other_np__wgt)$time) ok(ut_cmp_equal(
    r$dstart_other_np__wgt[,,"2000-01"],
    r$dstart_other_np__wgt[,,t],
    ignore = NULL ), paste0("dstart_other_np__wgt: Stable across timestep ", t, ", not affected by catch"))

## other_cv abundance
ok(gadget3:::ut_cmp_array(r$dstart_other_cv__num[,,"2000-01", drop = FALSE], "
  length  age    time        Freq
1    50:60  age5 2000-01 6.901043e-02
2    60:70  age5 2000-01 4.074388e+00
3    70:80  age5 2000-01 9.667898e+01
4    80:90  age5 2000-01 9.219846e+02
5   90:100  age5 2000-01 3.533761e+03
6  100:Inf  age5 2000-01 5.443432e+03
7    50:60  age6 2000-01 6.683220e-02
8    60:70  age6 2000-01 3.954349e+00
9    70:80  age6 2000-01 9.430333e+01
10   80:90  age6 2000-01 9.064455e+02
11  90:100  age6 2000-01 3.511713e+03
12 100:Inf  age6 2000-01 5.483517e+03
13   50:60  age7 2000-01 6.605127e-02
14   60:70  age7 2000-01 3.911228e+00
15   70:80  age7 2000-01 9.344646e+01
16   80:90  age7 2000-01 9.008037e+02
17  90:100  age7 2000-01 3.503606e+03
18 100:Inf  age7 2000-01 5.498167e+03
19   50:60  age8 2000-01 6.576669e-02
20   60:70  age8 2000-01 3.895503e+00
21   70:80  age8 2000-01 9.313351e+01
22   80:90  age8 2000-01 8.987383e+02
23  90:100  age8 2000-01 3.500623e+03
24 100:Inf  age8 2000-01 5.503543e+03
25   50:60  age9 2000-01 6.566236e-02
26   60:70  age9 2000-01 3.889736e+00
27   70:80  age9 2000-01 9.301869e+01
28   80:90  age9 2000-01 8.979798e+02
29  90:100  age9 2000-01 3.499526e+03
30 100:Inf  age9 2000-01 5.505520e+03
31   50:60 age10 2000-01 6.562402e-02
32   60:70 age10 2000-01 3.887617e+00
33   70:80 age10 2000-01 9.297649e+01
34   80:90 age10 2000-01 8.977010e+02
35  90:100 age10 2000-01 3.499123e+03
36 100:Inf age10 2000-01 5.506246e+03
", tolerance = 1e-7), "dstart_other_cv__num: vonB applied to all ages")
for (t in dimnames(r$dstart_other_cv__num)$time) ok(ut_cmp_equal(
    r$dstart_other_cv__num[,,"2000-01"],
    r$dstart_other_cv__num[,,t],
    ignore = NULL ), paste0("dstart_other_cv__num: Stable across timestep ", t, ", not affected by catch"))
for (t in dimnames(r$dstart_other_cv__wgt)$time) ok(ut_cmp_equal(
    r$dstart_other_cv__wgt[,,"2000-01"],
    r$dstart_other_cv__wgt[,,t],
    ignore = NULL ), paste0("dstart_other_cv__wgt: Stable across timestep ", t, ", not affected by catch"))

gadget3:::ut_tmb_r_compare2(model_fn, model_cpp, params)
######## Default params