R Under development (unstable) (2025-04-18 r88159 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2025 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.

> if (!interactive()) options(warn=2, error = function() { sink(stderr()) ; traceback(3) ; q(status = 1) })
> library(unittest)
> 
> library(gadget3)
> 
> st_imm <- g3_stock(c(species = "st", maturity = "imm"), c(10, 20, 30)) |> g3s_age(1, 5)
> st_mat <- g3_stock(c(species = "st", maturity = "mat"), c(10, 20, 30)) |> g3s_age(3, 15)
> stocks_st <- list(st_imm, st_mat)
> fl1 <- g3_fleet(c(type = "surv", 1))
> fl2 <- g3_fleet(c(type = "surv", 2))
> stocks_fl <- list(fl1, fl2)
> 
> assess_outputs <- list()
> 
> # Assessment function, gets pulled into model by g3_formula
> assess_fn <- function (
+         # The start model year, as defined by g3a_time,
+         start_year,
+         # The current model year, as defined by g3a_time,
+         cur_year,
+         # Nested list of pred -> prey -> detail_prey_pred__cons
+         cons,
+         # List of prey -> dstart_prey__num
+         abund,
+         # List of prey -> dstart_prey__wgt
+         meanwgt ) {
+     years <- seq(start_year, cur_year)
+ 
+     ## catch in numbers at age over all fleets
+     cn <- g3_array_combine(lapply(names(cons), function (pred_n) lapply(names(cons[[pred_n]]), function (prey_n) {
+         g3_array_agg(cons[[pred_n]][[prey_n]] / pmax(meanwgt[[prey_n]], 0.001), c("age", "year"), year = years)
+     })))
+ 
+     ## Abundance by age at step 1
+     smb <- g3_array_combine(lapply(names(abund), function (prey_n) {
+         g3_array_agg(abund[[prey_n]], c("age", "year"), year = years, step = 1)
+     }))
+ 
+     ## total abundance by maturity at step 1 by age
+     immtotal <- g3_array_agg(abund$st_imm, c("age", "year"), year = years, step = 1)
+     mattotal <- g3_array_agg(abund$st_mat, c("age", "year"), year = years, step = 1)
+ 
+     ## Log outputs in globalenv
+     assess_outputs[[as.character(cur_year)]] <<- list(
+         cn = cn,
+         smb = smb,
+         immtotal = immtotal,
+         mattotal = mattotal )
+ 
+     ## Perform the assessment
+     (sum(immtotal) + sum(mattotal)) / 1e10
+ }
> 
> fl_quota <- g3_quota(
+     g3_quota_assess(stocks_fl, stocks_st, g3_formula(
+         assess_fn(start_year, cur_year, cons, abund, meanwgt),
+         assess_fn = assess_fn )),
+     run_revstep = -3L,  # Run in spring
+     start_step = 4L,  # Skip first year, so we run 2000.2005
+     year_length = 5L )  # Run every years
> 
> actions <- list(
+     g3a_time(1990, 1994, c(3,3,3,3)),
+     g3a_otherfood_normalcv(
+         st_imm,
+         factor_f = g3_timeareadata('st_abund', data.frame(
+             year = 1990:2050,
+             abund = 1e6 - 1e4 * seq(0, 2050-1990)), "abund")),
+     g3a_otherfood_normalcv(
+         st_mat,
+         factor_f = g3_timeareadata('st_abund', data.frame(
+             year = 1990:2050,
+             abund = 1e6 - 1e4 * seq(0, 2050-1990)), "abund")),
+     g3a_predate(
+         fl1,
+         stocks_st,
+         catchability_f = g3a_predate_catchability_project(
+             fl_quota,
+             g3_parameterized("landings", value = 0, by_year = TRUE, by_predator = TRUE) ),
+         g3_suitability_exponentiall50(by_stock = 'species') ),
+     g3a_predate(
+         fl2,
+         stocks_st,
+         catchability_f = g3a_predate_catchability_project(
+             fl_quota,
+             g3_parameterized("landings", value = 0, by_year = TRUE, by_predator = TRUE) ),
+         g3_suitability_exponentiall50(by_stock = 'species') ),
+     # NB: Dummy parameter so model will compile in TMB
+     quote( nll <- nll + g3_param("x", value = 0) ) )
> full_actions <- c(actions, list(
+     g3a_report_detail(actions),
+     g3a_report_history(actions, 'proj_.*', out_prefix = NULL),
+     g3a_report_history(actions, 'quota_.*', out_prefix = NULL),
+     NULL))
> model_fn <- g3_to_r(full_actions)
> model_cpp <- g3_to_tmb(full_actions)
> 
> attr(model_cpp, 'parameter_template') |>
+     # Project for 30 years
+     g3_init_val("project_years", 30) |>
+ 
+     g3_init_val("*.K", 0.3, lower = 0.04, upper = 1.2) |>
+     g3_init_val("*.Linf", max(g3_stock_def(st_imm, "midlen")), spread = 0.2) |>
+     g3_init_val("*.t0", g3_stock_def(st_imm, "minage") - 0.8, spread = 2) |>
+     g3_init_val("*.lencv", 0.1, optimise = FALSE) |>
+     g3_init_val("*.walpha", 0.01, optimise = FALSE) |>
+     g3_init_val("*.wbeta", 3, optimise = FALSE) |>
+     g3_init_val("*.*.l50", g3_stock_def(st_imm, "midlen")[[length(g3_stock_def(st_imm, "midlen")) / 2]], spread = 0.25) |>
+ 
+     # surv_1 takes 80% of the quota
+     g3_init_val("surv_1.quota.prop", 0.8) |>
+     g3_init_val("surv_2.quota.prop", 0.2) |>
+     # surv_1 active at 1, surv_2 active at 2/3
+     g3_init_val("surv_1.cons.step.#", c(1, 0, 0, 0)) |>
+     g3_init_val("surv_2.cons.step.#", c(0, 0.6, 0.4, 0)) |>
+ 
+     identity() -> params.in
> nll <- model_fn(params.in) ; r <- attributes(nll) ; nll <- as.vector(nll)
> 
> ok(ut_cmp_identical(
+     names(assess_outputs),
+     c("1990", "1995", "2000", "2005", "2010", "2015", "2020") ), "assess_outputs: Ran assessment according to fishing year")
ok - assess_outputs: Ran assessment according to fishing year
> 
> ok(ut_cmp_identical(sapply(assess_outputs, function (x) dimnames(x$mattotal)$year), list(
+     "1990" = as.character(seq(1990, 1990)),
+     "1995" = as.character(seq(1990, 1995)),
+     "2000" = as.character(seq(1990, 2000)),
+     "2005" = as.character(seq(1990, 2005)),
+     "2010" = as.character(seq(1990, 2010)),
+     "2015" = as.character(seq(1990, 2015)),
+     "2020" = as.character(seq(1990, 2020)) )), "assess_outputs$mattotal: Correct years aggreagated")
ok - assess_outputs$mattotal: Correct years aggreagated
> 
> for (yr in seq(1995, 2020, by = 5)) ok_group(paste0("Year = ", yr), {
+     if (yr > 1995) {
+         x <- assess_outputs[[as.character(yr - 5)]]$cn[,as.character(seq(1990, yr - 6)), drop = FALSE]
+         y <- assess_outputs[[as.character(yr)]]$cn[,as.character(seq(1990, yr - 6)), drop = FALSE]
+         ok(ut_cmp_equal(x, y), paste0("assess_outputs$cn$", yr, ": Matches previous year, bar final year"))
+         ok(all(is.na(assess_outputs[[as.character(yr)]]$cn[,as.character(yr),drop = FALSE])), paste0("assess_outputs$cn$", yr, ": Final year NA"))
+     }
+     ok(ut_cmp_equal(
+         assess_outputs[[as.character(yr - 5)]]$smb,
+         assess_outputs[[as.character(yr)]]$smb[,as.character(seq(1990, yr - 5)), drop = FALSE] ), paste0("assess_outputs$smb$", yr, ": Matches previous year"))
+     ok(ut_cmp_equal(
+         assess_outputs[[as.character(yr - 5)]]$immtotal,
+         assess_outputs[[as.character(yr)]]$immtotal[,as.character(seq(1990, yr - 5)), drop = FALSE] ), paste0("assess_outputs$immtotal$", yr, ": Matches previous year"))
+     ok(ut_cmp_equal(
+         assess_outputs[[as.character(yr - 5)]]$mattotal,
+         assess_outputs[[as.character(yr)]]$mattotal[,as.character(seq(1990, yr - 5)), drop = FALSE] ), paste0("assess_outputs$mattotal$", yr, ": Matches previous year"))
+ })
# Year = 1995
ok - assess_outputs$smb$1995: Matches previous year
ok - assess_outputs$immtotal$1995: Matches previous year
ok - assess_outputs$mattotal$1995: Matches previous year
# Year = 2000
ok - assess_outputs$cn$2000: Matches previous year, bar final year
ok - assess_outputs$cn$2000: Final year NA
ok - assess_outputs$smb$2000: Matches previous year
ok - assess_outputs$immtotal$2000: Matches previous year
ok - assess_outputs$mattotal$2000: Matches previous year
# Year = 2005
ok - assess_outputs$cn$2005: Matches previous year, bar final year
ok - assess_outputs$cn$2005: Final year NA
ok - assess_outputs$smb$2005: Matches previous year
ok - assess_outputs$immtotal$2005: Matches previous year
ok - assess_outputs$mattotal$2005: Matches previous year
# Year = 2010
ok - assess_outputs$cn$2010: Matches previous year, bar final year
ok - assess_outputs$cn$2010: Final year NA
ok - assess_outputs$smb$2010: Matches previous year
ok - assess_outputs$immtotal$2010: Matches previous year
ok - assess_outputs$mattotal$2010: Matches previous year
# Year = 2015
ok - assess_outputs$cn$2015: Matches previous year, bar final year
ok - assess_outputs$cn$2015: Final year NA
ok - assess_outputs$smb$2015: Matches previous year
ok - assess_outputs$immtotal$2015: Matches previous year
ok - assess_outputs$mattotal$2015: Matches previous year
# Year = 2020
ok - assess_outputs$cn$2020: Matches previous year, bar final year
ok - assess_outputs$cn$2020: Final year NA
ok - assess_outputs$smb$2020: Matches previous year
ok - assess_outputs$immtotal$2020: Matches previous year
ok - assess_outputs$mattotal$2020: Matches previous year
> 
> ok(ut_cmp_equal(
+     g3_array_agg(r$detail_st_imm_surv_1__cons / r$dstart_st_imm__wgt, c("year"), age = 5, year = 1990:2019) +
+     g3_array_agg(r$detail_st_mat_surv_1__cons / r$dstart_st_mat__wgt, c("year"), age = 5, year = 1990:2019) +
+     g3_array_agg(r$detail_st_imm_surv_2__cons / r$dstart_st_imm__wgt, c("year"), age = 5, year = 1990:2019) +
+     g3_array_agg(r$detail_st_mat_surv_2__cons / r$dstart_st_mat__wgt, c("year"), age = 5, year = 1990:2019) +
+     0,
+     assess_outputs[["2020"]]$cn["age5", as.character(1990:2019)] ), "assess_outputs$2020$cn: Matches reporting output")
ok - assess_outputs$2020$cn: Matches reporting output
> 
> # NB: Not testing TMB, doesn't make sense to
> 
> proc.time()
   user  system elapsed 
   4.78    0.17    4.92 
1..31
# Looks like you passed all 31 tests.