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)
}
})
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.