tests/test-action_predate-catchability.R

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

library(magrittr)
library(unittest)

library(gadget3)

prey_a <- g3_stock('prey_a', seq(1, 10))
prey_b <- g3_stock('prey_b', seq(1, 10))
fleet_totalfleet <- g3_fleet('fleet_totalfleet')
fleet_numberfleet <- g3_fleet('fleet_numberfleet')
fleet_linearfleet <- g3_fleet('fleet_linearfleet')
fleet_effortfleet <- g3_fleet('fleet_effortfleet')
fleet_quotafleet <- g3_fleet('fleet_quotafleet')

actions <- list(
    g3a_time(2000, 2005, step_lengths = c(6,6), project_years = 0),
    g3a_initialconditions(prey_a, ~1e10 * prey_a__midlen, ~100 * prey_a__midlen),
    g3a_initialconditions(prey_b, ~2e10 * prey_b__midlen, ~200 * prey_b__midlen),
    g3a_predate_fleet(
        fleet_totalfleet,
        list(prey_a, prey_b),
        suitabilities = list(prey_a = 0.1, prey_b = 0.1),
        catchability_f = g3a_predate_catchability_totalfleet(10)),
    g3a_predate_fleet(
        fleet_numberfleet,
        list(prey_a, prey_b),
        suitabilities = list(prey_a = 0.1, prey_b = 0.1),
        catchability_f = g3a_predate_catchability_numberfleet(10)),
    g3a_predate_fleet(
        fleet_linearfleet,
        list(prey_a, prey_b),
        suitabilities = list(prey_a = 0.1, prey_b = 0.1),
        catchability_f = g3a_predate_catchability_linearfleet(1e-5)),
    g3a_predate_fleet(
        fleet_effortfleet,
        list(prey_a, prey_b),
        suitabilities = list(prey_a = 0.1, prey_b = 0.1),
        catchability_f = g3a_predate_catchability_effortfleet(list(prey_a = 1e-3, prey_b= 1e-4), 10)),
    g3a_predate_fleet(
        fleet_quotafleet,
        list(prey_a, prey_b),
        suitabilities = list(prey_a = 0.1, prey_b = 0.1),
        catchability_f = g3a_predate_catchability_quotafleet(
            data.frame(
                biomass = c(100, Inf),
                quota = I(list(~g3_param('quota.low', value = 1e-9), ~g3_param('quota.high', value = 1e-5)))),
            E = 2,
            recalc_f = ~cur_step == 1)),
    list())
actions <- c(actions, list(
    g3a_report_detail(actions) ))

# Compile model
model_fn <- g3_to_r(actions, trace = FALSE)
# model_fn <- edit(model_fn)
if (nzchar(Sys.getenv('G3_TEST_TMB'))) {
    model_cpp <- g3_to_tmb(actions, trace = FALSE)
    # model_cpp <- edit(model_cpp)
    model_tmb <- g3_tmb_adfun(model_cpp, compile_flags = c("-O0", "-g"))
} else {
    writeLines("# skip: not compiling TMB model")
}

ok_group("Catchability", {
    params <- attr(model_fn, 'parameter_template')

    result <- model_fn(params)
    r <- attributes(result)
    # str(as.list(r), vec.len = 10000)

    # Make sure everything is internally consistent
    ok(all(r$report_a__consratio > 0.9499), "report_a__consratio: No overconsumption")
    ok(all(r$report_b__consratio > 0.9499), "report_b__consratio: No overconsumption")
    for (time_idx in head(seq_along(dimnames(r$dstart_prey_b__num)$time), -1)) {
        time_lbl <- dimnames(r$dstart_prey_b__num)$time[[time_idx]]
        ok(ut_cmp_equal(
             r$dstart_prey_a__num[,time = time_idx] * r$dstart_prey_a__wgt[,time = time_idx] - 
                 r$detail_prey_a_fleet_totalfleet__cons[,time = time_idx] -
                 r$detail_prey_a_fleet_numberfleet__cons[,time = time_idx] -
                 r$detail_prey_a_fleet_linearfleet__cons[,time = time_idx] -
                 r$detail_prey_a_fleet_effortfleet__cons[,time = time_idx] -
                 r$detail_prey_a_fleet_quotafleet__cons[,time = time_idx],
             r$dstart_prey_a__num[,time = time_idx + 1] * r$dstart_prey_a__wgt[,time = time_idx + 1],
             tolerance = 1e-6), paste0(time_lbl, "/prey_a: Consumption and abundance consistent"))
        ok(ut_cmp_equal(
             r$dstart_prey_b__num[,time = time_idx] * r$dstart_prey_b__wgt[,time = time_idx] - 
                 r$detail_prey_b_fleet_totalfleet__cons[,time = time_idx] -
                 r$detail_prey_b_fleet_numberfleet__cons[,time = time_idx] -
                 r$detail_prey_b_fleet_linearfleet__cons[,time = time_idx] -
                 r$detail_prey_b_fleet_effortfleet__cons[,time = time_idx] -
                 r$detail_prey_b_fleet_quotafleet__cons[,time = time_idx],
             r$dstart_prey_b__num[,time = time_idx + 1] * r$dstart_prey_b__wgt[,time = time_idx + 1],
             tolerance = 1e-6), paste0(time_lbl, "/prey_b: Consumption and abundance consistent"))

        ok(ut_cmp_equal(
            colSums(r$detail_prey_a_fleet_totalfleet__cons + r$detail_prey_b_fleet_totalfleet__cons)[[time_idx]],
            10,
            tolerance = 1e6), paste0(time_lbl, "/totalfleet: Total consumption adds to 10"))
        ok(ut_cmp_equal(
            colSums(r$detail_prey_a_fleet_numberfleet__cons / r$dstart_prey_a__wgt + r$detail_prey_b_fleet_numberfleet__cons / r$dstart_prey_b__wgt)[[time_idx]],
            10,
            tolerance = 1e6), paste0(time_lbl, "/numberfleet: Total consumption in numbers adds to 10"))

        ok(ut_cmp_equal(
            r$detail_prey_a_fleet_linearfleet__cons[,time_idx],
            (r$dstart_prey_a__num * r$dstart_prey_a__wgt * 0.5 * 0.1 * 1e-5)[,time_idx],
            tolerance = 1e6), paste0(time_lbl, "/prey_a/linearfleet: Total consumption proportion of overall stock"))

        ok(ut_cmp_equal(
            r$detail_prey_b_fleet_linearfleet__cons[,time_idx],
            (r$dstart_prey_b__num * r$dstart_prey_b__wgt * 0.5 * 0.1 * 1e-5)[,time_idx],
            tolerance = 1e6), paste0(time_lbl, "/prey_b/linearfleet: Total consumption proportion of overall stock"))

        ok(ut_cmp_equal(
            r$detail_prey_a_fleet_effortfleet__cons[,time_idx],
            (r$dstart_prey_a__num * r$dstart_prey_a__wgt * 0.5 * 0.1 * 1e-3 * 10)[,time_idx],
            tolerance = 1e6), paste0(time_lbl, "/prey_a/effortfleet: Total consumption proportion of overall stock (using prey catchability)"))
        ok(ut_cmp_equal(
            r$detail_prey_b_fleet_effortfleet__cons[,time_idx],
            (r$dstart_prey_b__num * r$dstart_prey_b__wgt * 0.5 * 0.1 * 1e-4 * 10)[,time_idx],
            tolerance = 1e6), paste0(time_lbl, "/prey_b/effortfleet: Total consumption proportion of overall stock (using prey catchability)"))

        ok(ut_cmp_equal(
            r$detail_prey_a_fleet_quotafleet__cons[,time_idx],
            (r$dstart_prey_a__num * r$dstart_prey_a__wgt * 0.5 * 0.1 * 1e-5 * 2)[,time_idx],
            tolerance = 1e6), paste0(time_lbl, "/prey_a/quotafleet: Total consumption proportion of overall stock (high quota)"))
        ok(ut_cmp_equal(
            r$detail_prey_b_fleet_quotafleet__cons[,time_idx],
            (r$dstart_prey_b__num * r$dstart_prey_b__wgt * 0.5 * 0.1 * 1e-5 * 2)[,time_idx],
            tolerance = 1e6), paste0(time_lbl, "/prey_b/quotafleet: Total consumption proportion of overall stock (high quota)"))
    }
    if (nzchar(Sys.getenv('G3_TEST_TMB'))) {
        param_template <- attr(model_cpp, "parameter_template")
        param_template$value <- params[param_template$switch]
        gadget3:::ut_tmb_r_compare(model_fn, model_tmb, param_template)
    }
})