R Under development (unstable) (2024-03-18 r86148 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2024 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> 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")
+ })
# g3a_spawn_recruitment_fecundity
ok - g3a_spawn_recruitment_fecundity$s
ok - 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")
+ })
# g3a_spawn_recruitment_simplessb
ok - g3a_spawn_recruitment_simplessb$s
ok - 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 - 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 - 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)
+     }
+ })
# g3a_spawn
# skip: not compiling TMB model
> 
> proc.time()
   user  system elapsed 
   1.43    0.06    1.48 
# Looks like you passed all 6 tests.