tests/test-action_spawn.R

if (!interactive()) options(warn=2, error = function() { sink(stderr()) ; traceback(3) ; q(status = 1) })
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

    # Add spawn step, spawning into a new stock with a name matching the recruitment function
    recruitment_test_step <- function (recruitment_f) {
        fn_name <- as.character(sys.call()[[2]][[1]])
        mat_stock <- g3_stock(c("mat", fn_name), seq(10, 40, 10)) %>% g3s_age(5, 5)
        imm_stock <- g3_stock(c("imm", fn_name), seq(10, 40, 10)) %>% g3s_age(0, 0)

        return(list(
            g3a_otherfood(  # Matures steadily grow over model run
                mat_stock,
                num_f = g3_formula(
                    cur_time * s1 + stock__midlen * 0,
                    s1 = g3_parameterized("abund_scalar", value = runif(1, min=90, max=110)) ),
                wgt_f = g3_formula(
                    cur_time * s2 + stock__midlen * 0,
                    s2 = g3_parameterized("mwgt_scalar", value = runif(1, min=9, max=11)) )),
            g3a_otherfood(  # Empty immature after each step
                imm_stock,
                num_f = quote( stock__midlen * 0 ),
                wgt_f = 1 ),
            g3a_spawn(
                mat_stock,
                recruitment_f,
                output_stocks = list(imm_stock) )))
    }
    actions <- list(
        g3a_time(min(year_range), max(year_range), c(3,3,3,3), project_years = 0),
        recruitment_test_step(g3a_spawn_recruitment_fecundity(
            g3_parameterized("fecundity_p0", value = runif(1, min=0.1, max=0.9)),
            g3_parameterized("fecundity_p1", value = runif(1, min=0.1, max=0.9)),
            g3_parameterized("fecundity_p2", value = runif(1, min=0.1, max=0.9)),
            g3_parameterized("fecundity_p3", value = runif(1, min=0.1, max=0.9)),
            g3_parameterized("fecundity_p4", value = runif(1, min=0.1, max=0.9)) )),
        recruitment_test_step(g3a_spawn_recruitment_simplessb(
            g3_parameterized("simplessb_mu", value = runif(1, min=0.1, max=0.9)) )),
        recruitment_test_step(g3a_spawn_recruitment_ricker(
            g3_parameterized("ricker_mu", value = runif(1, min=0.9, max=1.2)),
            g3_parameterized("ricker_lambda", value = runif(1, min=1e-6, max=1e-5)) )),
        recruitment_test_step(g3a_spawn_recruitment_bevertonholt(
            g3_parameterized("bevertonholt_mu", value = runif(1, min=0.4, max=0.6)),
            g3_parameterized("bevertonholt_lambda", value = runif(1, min=0.4, max=0.6)) )),
        recruitment_test_step(g3a_spawn_recruitment_bevertonholt_ss3()),
        recruitment_test_step(g3a_spawn_recruitment_hockeystick(
            g3_parameterized("hockeystick_r0", value = runif(1, min=0.4, max=0.6)),
            g3_parameterized("hockeystick_blim", value = runif(1, min=0.4, max=0.6)) )),
        list())

    all_actions <- c(actions, list(
        g3a_report_history(actions, '__num$|__wgt$|__offspringnum$'),
        g3a_report_history(actions, '__midlen$', out_prefix = NULL),
        NULL ))

    # Compile model
    model_fn <- g3_to_r(all_actions, trace = FALSE)
    # model_fn <- edit(model_fn)
    if (nzchar(Sys.getenv('G3_TEST_TMB'))) {
        model_cpp <- g3_to_tmb(all_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')

    # g3a_spawn_recruitment_bevertonholt_ss3
    params <- g3_init_val(params, 'mat_*.spawn.h', runif(1, 1.4, 1.5))
    params <- g3_init_val(params, 'mat_*.spawn.R0', runif(1, 1e4 - 10, 1e4 + 10))
    bholt_ss3_R <- runif(length(year_range), 1e2 - 10, 1e2 + 10)
    params <- g3_init_val(params, 'mat_*.spawn.R.#', bholt_ss3_R)
    params <- g3_init_val(params, 'mat_*.spawn.B0', runif(1, 1e3 - 10, 1e3 + 10))

    r <- lapply(attributes(model_fn(params)), drop)  # NB: Drop all redundant dimensions,
    age <- 5
    # Repeat midlen history for each time bin
    hist_midlen <- array(r$mat_g3a_spawn_recruitment_fecundity__midlen, dim = c(4, 36))
    hist_biomass <- colSums(r$hist_mat_g3a_spawn_recruitment_simplessb__num * r$hist_mat_g3a_spawn_recruitment_simplessb__wgt)
    ok(ut_cmp_equal(
        as.vector(colSums(r$hist_mat_g3a_spawn_recruitment_simplessb__num)),
        4 * 0:35 * params$abund_scalar,
        end = NULL ), "hist_mat_g3a_spawn_recruitment_simplessb__num: Test rig produced increasing abundance" )
    ok(ut_cmp_equal(
        as.vector(colSums(r$hist_mat_g3a_spawn_recruitment_simplessb__wgt)),
        4 * 0:35 * params$mwgt_scalar,
        end = NULL ), "hist_mat_g3a_spawn_recruitment_simplessb__wgt: Test rig produced increasing mean weight" )

    ok(ut_cmp_equal(
        as.vector(colSums(r$hist_mat_g3a_spawn_recruitment_fecundity__offspringnum) > 0),
        c(FALSE, rep(TRUE, 36 - 1)) ), "hist_mat_g3a_spawn_recruitment_fecundity__offspringnum: Some recruitment happened")
    ok(ut_cmp_equal(
        colSums(r$hist_mat_g3a_spawn_recruitment_fecundity__offspringnum),
        colSums(
            hist_midlen ^ params$fecundity_p1 *
            age ^ params$fecundity_p2 *
            r$hist_mat_g3a_spawn_recruitment_simplessb__num^params$fecundity_p3 *
            r$hist_mat_g3a_spawn_recruitment_simplessb__wgt^params$fecundity_p4 ) * params$fecundity_p0,
        end = NULL), "hist_mat_g3a_spawn_recruitment_fecundity__offspringnum: Matches formula")

    ok(ut_cmp_equal(
        as.vector(colSums(r$hist_mat_g3a_spawn_recruitment_simplessb__offspringnum) > 0),
        c(FALSE, rep(TRUE, 36 - 1)) ), "hist_mat_g3a_spawn_recruitment_simplessb__offspringnum: Some recruitment happened")
    ok(ut_cmp_equal(
        colSums(r$hist_mat_g3a_spawn_recruitment_simplessb__offspringnum),
        hist_biomass * params$simplessb_mu,
        end = NULL ), "hist_mat_g3a_spawn_recruitment_simplessb__offspringnum: Total biomass * mu")

    ok(ut_cmp_equal(
        as.vector(colSums(r$hist_mat_g3a_spawn_recruitment_ricker__offspringnum) > 0),
        c(FALSE, rep(TRUE, 36 - 1)) ), "hist_mat_g3a_spawn_recruitment_ricker__offspringnum: Some recruitment happened")
    ok(ut_cmp_equal(
        colSums(r$hist_mat_g3a_spawn_recruitment_ricker__offspringnum),
        params$ricker_mu * hist_biomass * exp(-params$ricker_lambda * hist_biomass),
        end = NULL ), "hist_mat_g3a_spawn_recruitment_ricker__offspringnum: Matches formula")

    ok(ut_cmp_equal(
        as.vector(colSums(r$hist_mat_g3a_spawn_recruitment_bevertonholt__offspringnum) > 0),
        c(FALSE, rep(TRUE, 36 - 1)) ), "hist_mat_g3a_spawn_recruitment_bevertonholt__offspringnum: Some recruitment happened")
    ok(ut_cmp_equal(
        colSums(r$hist_mat_g3a_spawn_recruitment_bevertonholt__offspringnum),
        params$bevertonholt_mu * hist_biomass/(params$bevertonholt_lambda + hist_biomass),
        end = NULL ), "hist_mat_g3a_spawn_recruitment_bevertonholt__offspringnum: Matches formula")

    ok(ut_cmp_equal(
        as.vector(colSums(r$hist_mat_g3a_spawn_recruitment_bevertonholt__offspringnum) > 0),
        c(FALSE, rep(TRUE, 36 - 1)) ), "hist_mat_g3a_spawn_recruitment_bevertonholt__offspringnum: Some recruitment happened")
    ok(ut_cmp_equal(
        colSums(r$hist_mat_g3a_spawn_recruitment_bevertonholt__offspringnum),
        params$bevertonholt_mu * hist_biomass/(params$bevertonholt_lambda + hist_biomass),
        end = NULL ), "hist_mat_g3a_spawn_recruitment_bevertonholt__offspringnum: Matches formula")

    ok(ut_cmp_equal(
        as.vector(colSums(r$hist_mat_g3a_spawn_recruitment_bevertonholt_ss3__offspringnum) > 0),
        c(FALSE, rep(TRUE, 36 - 1)) ), "hist_mat_g3a_spawn_recruitment_bevertonholt_ss3__offspringnum: Some recruitment happened")
    ok(ut_cmp_equal(
        colSums(r$hist_mat_g3a_spawn_recruitment_bevertonholt_ss3__offspringnum),
        4 * params$mat_g3a_spawn_recruitment_bevertonholt_ss3.spawn.h *
        params$mat_g3a_spawn_recruitment_bevertonholt_ss3.spawn.R0 *
        hist_biomass *
        rep(bholt_ss3_R, each = 4) / (
            params$mat_g3a_spawn_recruitment_bevertonholt_ss3.spawn.B0 *
            (1 - params$mat_g3a_spawn_recruitment_bevertonholt_ss3.spawn.h)
          +
            hist_biomass *
            (5 * params$mat_g3a_spawn_recruitment_bevertonholt_ss3.spawn.h - 1)),
        end = NULL ), "hist_mat_g3a_spawn_recruitment_bevertonholt_ss3__offspringnum: Matches formula")

    ok(ut_cmp_equal(
        as.vector(colSums(r$hist_mat_g3a_spawn_recruitment_hockeystick__offspringnum) > 0),
        c(FALSE, rep(TRUE, 36 - 1)) ), "hist_mat_g3a_spawn_recruitment_hockeystick__offspringnum: Some recruitment happened")
    ok(ut_cmp_equal(
        colSums(r$hist_mat_g3a_spawn_recruitment_hockeystick__offspringnum),
        params$hockeystick_r0 * pmin(hist_biomass/params$hockeystick_blim, 1),
        end = NULL ), "hist_mat_g3a_spawn_recruitment_hockeystick__offspringnum: Matches formula")

    # 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)
    }
})
gadget-framework/gadget3 documentation built on June 13, 2025, 5:06 a.m.