tests/test-action_spawn.R

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

library(magrittr)
library(unittest)

library(gadget3)

stock_imm1 <- g3_stock('stock_imm1', seq(10, 40, 10)) %>% g3s_age(3, 7)
stock_imm2 <- g3_stock('stock_imm2', seq(10, 40, 10)) %>% g3s_age(4, 7)
stock_mat <- g3_stock('stock_mat', seq(30, 40, 10)) %>% g3s_age(5, 7)

ok_group('g3a_spawn_recruitment_fecundity', {
    ok(ut_cmp_identical(
        rlang::f_rhs( g3a_spawn_recruitment_fecundity(90, 91, 92, 93, 94)$s ),
        quote( sum(stock__midlen^91 * age^92 * stock_ss(stock__spawningnum)^93 * stock_ss(stock__wgt)^94) )), "g3a_spawn_recruitment_fecundity$s")
    ok(ut_cmp_identical(
        rlang::f_rhs( g3a_spawn_recruitment_fecundity(90, 91, 92, 93, 94)$r ),
        quote( 90 * s  )), "g3a_spawn_recruitment_fecundity$r")
})

ok_group('g3a_spawn_recruitment_simplessb', {
    ok(ut_cmp_identical(
        rlang::f_rhs( g3a_spawn_recruitment_simplessb(91)$s ),
        quote( sum(stock_ss(stock__wgt) * stock_ss(stock__spawningnum)) )), "g3a_spawn_recruitment_simplessb$s")
    ok(ut_cmp_identical(
        rlang::f_rhs( g3a_spawn_recruitment_simplessb(91)$r ),
        quote( 91 * s )), "g3a_spawn_recruitment_simplessb$r")
})

ok(ut_cmp_error(
    g3a_spawn(stock_mat, output_stocks = list(stock_imm1, stock_imm2), output_ratios = c(9,9,9), recruitment_f = list(s = 1, r = 1)),
    "output_ratios"), "Length of output_ratios must match")
ok(ut_cmp_error(
    g3a_spawn(stock_mat, output_stocks = list(stock_imm1, stock_imm2), output_ratios = c(9,9), recruitment_f = list(s = 1, r = 1)),
    "output_ratios"), "output_ratios must sum to 1")

ok_group('g3a_spawn', {
    year_range <- 1982:1990

    ling_imm <- g3_stock('ling_imm', seq(20, 156, 4)) %>%
        g3s_livesonareas(c(1)) %>%
        g3s_age(3, 10)

    ling_mat <- g3_stock('ling_mat', seq(20, 156, 4)) %>%
        g3s_livesonareas(c(1)) %>%
        g3s_age(5, 15)

    igfs <- g3_fleet('igfs') %>% g3s_livesonareas(c(1))

    imm_report <- g3s_clone(ling_imm, 'imm_report') %>% g3s_time(year = local(year_range), step = 1:4)
    mat_report <- g3s_clone(ling_mat, 'mat_report') %>% g3s_time(year = local(year_range), step = 1:4)

    ling_imm_actions <- list(
        g3a_initialconditions_normalparam(ling_imm,
            factor_f = ~age * g3_param("lingimm.init") * g3_param("lingimm.init.scalar"),
            mean_f = ~g3_param("ling.Linf"),
            stddev_f = ~10,
            alpha_f = ~g3_param("lingimm.walpha"),
            beta_f = ~g3_param("lingimm.wbeta")),
        g3a_naturalmortality(ling_imm, g3a_naturalmortality_exp(~g3_param("lingimm.M"))),
        g3a_age(ling_imm),
        list())

    ling_mat_actions <- list(
        g3a_initialconditions_normalparam(ling_mat,
            factor_f = ~age * g3_param("lingmat.init") * g3_param("lingmat.init.scalar"),
            mean_f = ~g3_param("ling.Linf"),
            stddev_f = ~10,
            alpha_f = ~g3_param("lingmat.walpha"),
            beta_f = ~g3_param("lingmat.wbeta")),
        g3a_naturalmortality(ling_mat, g3a_naturalmortality_exp(~g3_param("lingmat.M"))),
        g3a_age(ling_mat),
        g3a_spawn(
            ling_mat,
            recruitment_f = g3a_spawn_recruitment_ricker(
                ~g3_param("ricker.mu"),
                ~g3_param("ricker.lambda")),
            proportion_f = g3_suitability_exponentiall50(alpha = ~-g3_param("spawn.prop.alpha"), l50 =  ~g3_param("spawn.prop.l50")),
            mortality_f = g3_suitability_straightline(alpha = ~g3_param("spawn.mort.alpha"), beta =  ~g3_param("spawn.mort.beta")),
            weightloss_f = ~g3_param("spawn.weightloss"),
            output_stocks = list(ling_imm),
            mean_f = 50,
            stddev_f = 0.9,
            alpha_f = 1,
            beta_f = 1,
            run_f = ~cur_step==1),
        list())

    report_actions <- list(
           g3a_report_stock(imm_report,ling_imm, ~stock_ss(ling_imm__num)),
           g3a_report_stock(imm_report,ling_imm, ~stock_ss(ling_imm__wgt)),
           g3a_report_stock(mat_report,ling_mat, ~stock_ss(ling_mat__num)),
           g3a_report_stock(mat_report,ling_mat, ~stock_ss(ling_mat__wgt)),
           list())

    time_actions <- list(
        g3a_time(min(year_range), max(year_range), c(3,3,3,3), project_years = 0),
        list())

    # Add steps to exercise rest of recruitment functions, and check they produce identical TMB results
    recruitment_test_step <- function (recruitment_f) {
        action_name <- gadget3:::unique_action_name()

        # Re-implement enough of spawning to test recruitment
        stock <- ling_mat
        out_var_name <- paste0('stock__rf_', sys.call()[[2]][[1]])
        assign(out_var_name, g3_stock_instance(stock))
        out <- list()
        out[[gadget3:::step_id(999, action_name)]] <- gadget3:::g3_step(gadget3:::f_substitute(
            ~g3_with(s := 0 * nll, {  # TODO: Ugly mess to get type right
                stock_iterate(stock, if (run_f) {
                    s <- s + recruitment_s_f
                    stock_ss(stock__outvar) <- 1
                } else {
                    stock_ss(stock__outvar) <- 0
                })
                g3_with(r := recruitment_r_f,
                    stock_with(stock, stock__outvar <- r * stock__outvar / avoid_zero(sum(stock__outvar))))
            }),
            list(
                recruitment_r_f = recruitment_f$r,
                recruitment_s_f = recruitment_f$s,
                stock__outvar = as.symbol(out_var_name))))
        return(out)
    }
    recruitment_f_actions <- list(
        recruitment_test_step(g3a_spawn_recruitment_fecundity(
            p0 = runif(1, min=0.1, max=0.9),
            p1 = runif(1, min=0.1, max=0.9),
            p2 = runif(1, min=0.1, max=0.9),
            p3 = runif(1, min=0.1, max=0.9),
            p4 = runif(1, min=0.1, max=0.9))),
        recruitment_test_step(g3a_spawn_recruitment_simplessb(runif(1, min=0.1, max=0.9))),
        recruitment_test_step(g3a_spawn_recruitment_ricker(runif(1, min=0.1, max=0.9), runif(1, min=0.1, max=0.9))),
        recruitment_test_step(g3a_spawn_recruitment_bevertonholt(runif(1, min=0.1, max=0.9), runif(1, min=0.1, max=0.9))),
        recruitment_test_step(g3a_spawn_recruitment_hockeystick(runif(1, min=0.1, max=0.9), runif(1, min=0.1, max=0.9))),
        list())

    actions <- c(
      ling_imm_actions,
      ling_mat_actions,
      report_actions,
      recruitment_f_actions,
      time_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")
    }

    params <- attr(model_fn, 'parameter_template')
    params$lingimm.init <- 0
    params$lingimm.init.scalar <- 0
    params$lingimm.rec.scalar <- 100
    params$lingimm.M <- 0
    params$lingimm.walpha <- 1e-1
    params$lingimm.wbeta <- 2
    params$lingmat.init <- 1
    params$lingmat.init.scalar <- 100
    params$lingmat.rec.scalar <- 100
    params$lingmat.M <- 0.95
    params$lingmat.walpha <- 1e-6
    params$lingmat.wbeta <- 2
    params$ling.init.F <- 0.4
    params$ling.mat.alpha <- 0.01
    params$ling.mat.l50 <- 75
    params$ling.mat.beta <- 0.01
    params$ling.mat.a50 <- 7
    params$ling.Linf <- 160
    params$ling.bbin <- 6
    params$ling.k <- 10
    params$ricker.mu <- 1
    params$ricker.lambda <- 1e-6
    params$spawn.prop.alpha <- 0.5
    params$spawn.prop.l50 <- 120
    params$spawn.mort.alpha <- 0
    params$spawn.mort.beta <- 0
    params$spawn.weightloss <- 0.1

    # Make sure the inttest model produces identical output in TMB and R
    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)
    }
})