Last updated on 2025-07-09 07:49:51 CEST.
| Package | ERROR | NOTE | OK |
|---|---|---|---|
| gadget3 | 1 | 12 | |
| mfdb | 4 | 9 | |
| unittest | 13 |
Current CRAN status: ERROR: 1, OK: 12
Version: 0.13-0
Check: tests
Result: ERROR
Running 'test-aab_env.R' [1s]
Running 'test-action_age.R' [2s]
Running 'test-action_grow-methods.R' [4s]
Running 'test-action_grow.R' [2s]
Running 'test-action_mature.R' [2s]
Running 'test-action_migrate.R' [2s]
Running 'test-action_naturalmortality.R' [1s]
Running 'test-action_predate-catchability.R' [5s]
Running 'test-action_predate-numberfleet.R' [2s]
Running 'test-action_predate-predator.R' [11s]
Running 'test-action_predate-timebasedsuitability.R' [2s]
Running 'test-action_predate.R' [5s]
Running 'test-action_renewal-otherfood.R' [7s]
Running 'test-action_renewal.R' [4s]
Running 'test-action_report.R' [1s]
Running 'test-action_spawn-multipleoutputs.R' [4s]
Running 'test-action_spawn.R' [4s]
Running 'test-action_spmodel.R' [2s]
Running 'test-action_tagging.R' [2s]
Running 'test-action_time.R' [1s]
Running 'test-action_weightloss.R' [2s]
Running 'test-array_utils.R' [4s]
Running 'test-env_dif.R' [1s]
Running 'test-eval.R' [1s]
Running 'test-formula_utils.R' [1s]
Running 'test-init_val.R' [1s]
Running 'test-likelihood_bounds.R' [3s]
Running 'test-likelihood_data.R' [3s]
Running 'test-likelihood_distribution-surveyindices.R' [2s]
Running 'test-likelihood_distribution.R' [13s]
Running 'test-likelihood_random.R' [1s]
Running 'test-likelihood_sparsesample.R' [3s]
Running 'test-likelihood_tagging_ckmr.R' [3s]
Running 'test-likelihood_understocking.R' [0s]
Running 'test-param_project-ar1.R' [3s]
Running 'test-param_project-logar1.R' [3s]
Running 'test-param_project.R' [2s]
Running 'test-params.R' [2s]
Running 'test-quota-assess.R' [5s]
Running 'test-quota-hockeyfleet.R' [3s]
Running 'test-quota.R' [1s]
Running 'test-regression.R' [1s]
Running 'test-run.R' [1s]
Running 'test-run_r.R' [2s]
Running 'test-run_tmb-reporting_enabled.R' [1s]
Running 'test-run_tmb.R' [3s]
Running 'test-step.R' [2s]
Running 'test-stock.R' [1s]
Running 'test-stock_age.R' [2s]
Running 'test-stock_areas.R' [2s]
Running 'test-stock_product.R' [0s]
Running 'test-stock_tag.R' [1s]
Running 'test-stock_time-fishingyear.R' [2s]
Running 'test-stock_time.R' [1s]
Running 'test-suitability-report.R' [6s]
Running 'test-suitability.R' [2s]
Running 'test-timedata.R' [1s]
Running 'test-timevariable.R' [1s]
Running the tests in 'tests/test-likelihood_distribution.R' failed.
Complete output:
> library(magrittr)
> library(unittest)
>
> library(gadget3)
>
> # Zip name/value arguments together into a list
> named_list <- function(...) {
+ x <- list(...)
+ structure(
+ x[seq_along(x) %% 2 == 0],
+ names = as.character(x[seq_along(x) %% 2 == 1]))
+ }
>
> cmp_grep <- function (a, ...) {
+ re <- paste0('\\Q', c(...), '\\E', collapse = ".*")
+ if (grepl(re, a, perl = TRUE)) return(TRUE)
+ c(re, "Not found in:", a)
+ }
>
> # NB: Name has to be different, or it gets sucked into the model
> g3_avoid_zero <- g3_env$avoid_zero
> g3_logspace_add <- g3_env$logspace_add
> g3_lgamma_vec <- lgamma
>
> ok_group("g3_distribution_preview", {
+ dat <- expand.grid(year = 1990:1994, step = 2, area = 'IXa')
+ dat$number <- seq_len(nrow(dat))
+ out <- g3_distribution_preview(dat, area_group = c(IXa = 1))
+ ok(ut_cmp_equal(out, structure(
+ 1:5,
+ dim = c(length = 1L, time = 5L, area = 1L),
+ dimnames = list(
+ length = "0:Inf",
+ time = c("1990-02", "1991-02", "1992-02", "1993-02", "1994-02"),
+ area = "IXa") )), "Returned number array")
+
+ dat <- expand.grid(year = 1990:1994, step = 2, area = 'IXa')
+ dat$weight <- seq_len(nrow(dat)) * 40
+ out <- g3_distribution_preview(dat, area_group = c(IXa = 1))
+ ok(ut_cmp_equal(out, structure(
+ 1:5 * 40,
+ dim = c(length = 1L, time = 5L, area = 1L),
+ dimnames = list(
+ length = "0:Inf",
+ time = c("1990-02", "1991-02", "1992-02", "1993-02", "1994-02"),
+ area = "IXa") )), "Returned weight array")
+
+ dat <- expand.grid(year = 1990:1994, step = 2, area = 'IXa')
+ dat$number <- seq_len(nrow(dat)) * 9
+ dat$weight <- seq_len(nrow(dat)) * 40
+ out <- g3_distribution_preview(dat, area_group = c(IXa = 1))
+ ok(ut_cmp_equal(out, structure(
+ 1:5 * 9,
+ dim = c(length = 1L, time = 5L, area = 1L),
+ dimnames = list(
+ length = "0:Inf",
+ time = c("1990-02", "1991-02", "1992-02", "1993-02", "1994-02"),
+ area = "IXa") )), "Number array wins if both present")
+
+ stocks <- list(
+ g3_stock(c(species = 'fish', 'imm', 'f'), 1:10),
+ g3_stock(c(species = 'fish', 'imm', 'm'), 1:10),
+ g3_stock(c(species = 'fish', 'mat', 'f'), 1:10),
+ g3_stock(c(species = 'fish', 'mat', 'm'), 1:10) )
+ dat <- expand.grid(year = 1990:1994, stock = c("imm", "mat"), number = 0)
+ out <- drop(g3_distribution_preview(dat, stocks = stocks, area_group = c(IXa = 1)))
+ ok(ut_cmp_identical(out, structure(
+ array(0, dim = c(stock = 2L, time = 5L), dimnames = list(stock = c("imm", "mat"), time = c("1990", "1991", "1992", "1993", "1994"))),
+ stock_map = c(fish_imm_f = "imm", fish_imm_m = "imm", fish_mat_f = "mat", fish_mat_m = "mat") )), "Included stock map for stock column")
+
+ fleets <- list(
+ g3_fleet(c('comm', 'se')),
+ g3_fleet(c('comm', 'no')),
+ g3_fleet(c('surv')) )
+ dat <- expand.grid(year = 1990:1994, fleet = c("comm", "surv"), number = 0)
+ out <- drop(g3_distribution_preview(dat, fleets = fleets, area_group = c(IXa = 1)))
+ ok(ut_cmp_identical(out, structure(
+ array(0, dim = c(fleet = 2L, time = 5L), dimnames = list(fleet = c("comm", "surv"), time = c("1990", "1991", "1992", "1993", "1994"))),
+ fleet_map = c(comm_se = "comm", comm_no = "comm", surv = "surv") )), "Included fleet map for fleet column")
+ })
# g3_distribution_preview
ok - Returned number array
ok - Returned weight array
ok - Number array wins if both present
ok - Included stock map for stock column
ok - Included fleet map for fleet column
>
> ok_group("g3l_distribution_sumofsquares", {
+ ok(cmp_grep(
+ deparse1(environment(g3l_distribution_sumofsquares(c('area', 'age')))$modelstock__sstotal),
+ 'sum(stock_ssinv(modelstock__x, "time", "area", "age"))',
+ NULL), "Added custom totals to sumofsquares modelstock__x")
+ ok(cmp_grep(
+ deparse1(environment(g3l_distribution_sumofsquares(c('area', 'age')))$obsstock__sstotal),
+ 'sum(stock_ssinv(obsstock__x, "time", "area", "age"))',
+ NULL), "Added custom totals to sumofsquares modelstock__x")
+ ok(cmp_grep(
+ deparse1(g3l_distribution_sumofsquares(c('area', 'length'))),
+ 'stock_ss(modelstock__x, length = default)',
+ 'modelstock__sstotal[[modelstock__length_idx]]',
+ 'stock_ss(obsstock__x, length = default)',
+ 'obsstock__sstotal[[obsstock__length_idx]]',
+ NULL), "Adding length also adds to the stock_ss()")
+
+ # Stratified sumofsquares
+ prey_a <- g3_stock('prey_a', seq(1, 5, by = 1)) %>% g3s_age(1,5)
+ prey_a__init <- g3_stock_instance(prey_a)
+ prey_a__init[] <- runif(length(prey_a__init))
+ obsdata <- expand.grid(
+ year = 2000,
+ length = seq(1, 5, by = 1),
+ age = 3:4) # NB: Only report age 3,4
+ obsdata$number <- runif(nrow(obsdata))
+ model_fn <- g3_to_r(list(
+ # Keep TMB happy
+ g3_formula({
+ nll <- nll + g3_param("dummy", value = 0)
+ REPORT(prey_a__num)
+ }),
+ g3a_time(2000, 2001),
+ g3a_initialconditions(prey_a,
+ num_f = g3_formula(stock_ss(prey_a__init), prey_a__init = prey_a__init),
+ wgt_f = 10),
+ g3l_abundancedistribution("adist",
+ obsdata,
+ function_f = g3l_distribution_sumofsquares(over = c("area", "length")),
+ stocks = list(prey_a),
+ report = TRUE)))
+ model_cpp <- g3_to_tmb(attr(model_fn, 'actions'), trace = FALSE)
+ if (Sys.getenv('G3_TEST_TMB') == "2") {
+ #model_cpp <- edit(model_cpp)
+ #writeLines(TMB::gdbsource(g3_tmb_adfun(model_cpp, compile_flags = "-g", output_script = TRUE)))
+ model_tmb <- g3_tmb_adfun(model_cpp, trace = TRUE, compile_flags = c("-O0", "-g"))
+ }
+
+ params <- attr(model_fn, 'parameter_template')
+ r <- model_fn(params)
+ modeldata <- attr(r, 'prey_a__num')
+ expected_nll <- 0
+ for (age in 3:4) for (length in 1:5) {
+ # Proportion compared to others in age group, for that length
+ expected_nll <- expected_nll + (
+ # Proportion compared to others in this length group
+ modeldata[length, age] / sum(modeldata[length, c(3,4)]) -
+ obsdata[obsdata$age == age & obsdata$length == length, 'number'] / sum(obsdata[obsdata$length == length, 'number'])
+ ) ** 2
+ }
+ ok(ut_cmp_equal(as.numeric(r), expected_nll), "g3l_distribution_sumofsquares statified over length")
+ if (Sys.getenv('G3_TEST_TMB') == "2") gadget3:::ut_tmb_r_compare(model_fn, model_tmb, params, model_cpp = model_cpp)
+ })
# g3l_distribution_sumofsquares
ok - Added custom totals to sumofsquares modelstock__x
ok - Added custom totals to sumofsquares modelstock__x
ok - Adding length also adds to the stock_ss()
not ok - g3l_distribution_sumofsquares statified over length
# Test returned non-TRUE value:
# Mean relative difference: 0.001187689
# --- as.numeric(r)
# +++ expected_nll
# [1] [-1.369097-]{+1.370723+}
>
> ok_group("g3l_distribution:transform_fs", {
+ prey_a <- g3_stock('prey_a', seq(1, 5, by = 1)) %>% g3s_age(1,5) %>% g3s_livesonareas(1:2)
+ prey_a__init <- g3_stock_instance(prey_a)
+ prey_a__init[] <- runif(length(prey_a__init))
+ obsdata <- expand.grid(
+ year = 2000,
+ length = seq(1, 5, by = 1),
+ age = 1:5)
+ obsdata$number <- runif(nrow(obsdata))
+ model_fn <- g3_to_r(list(
+ g3a_time(2000, 2001),
+ g3a_initialconditions(prey_a,
+ num_f = g3_formula(stock_ss(prey_a__init), prey_a__init = prey_a__init),
+ wgt_f = 10),
+ g3l_abundancedistribution("wt",
+ obsdata,
+ function_f = g3l_distribution_sumofsquares(),
+ stocks = list(prey_a),
+ transform_fs = list(
+ age = g3_formula( g3_param_array('reader1matrix', value = diag(5))[g3_idx(preage), g3_idx(age)] )),
+ report = TRUE),
+ g3l_abundancedistribution("len",
+ obsdata,
+ function_f = g3l_distribution_sumofsquares(),
+ stocks = list(prey_a),
+ transform_fs = list(
+ length = quote(diag(10 * prey_a__midlen)) ),
+ report = TRUE),
+ g3l_abundancedistribution("wtperstock",
+ obsdata,
+ function_f = g3l_distribution_sumofsquares(),
+ stocks = list(prey_a),
+ transform_fs = list(age = list(
+ prey_a = g3_formula(g3_param_array('reader1matrix', value = diag(prey_a__maxage - prey_a__minage + 1))[prey_a__preage_idx, prey_a__age_idx] ))),
+ report = TRUE),
+ g3l_abundancedistribution("nt",
+ obsdata,
+ function_f = g3l_distribution_sumofsquares(),
+ stocks = list(prey_a),
+ report = TRUE),
+ # Keep TMB happy
+ g3_formula( nll <- nll + g3_param("dummy", value = 0) )))
+ model_cpp <- g3_to_tmb(attr(model_fn, 'actions'), trace = FALSE)
+ if (Sys.getenv('G3_TEST_TMB') == "2") {
+ #model_cpp <- edit(model_cpp)
+ #writeLines(TMB::gdbsource(g3_tmb_adfun(model_cpp, compile_flags = "-g", output_script = TRUE)))
+ model_tmb <- g3_tmb_adfun(model_cpp, trace = TRUE, compile_flags = c("-O0", "-g"))
+ }
+
+ # Given results / params, apply matrix manually and make sure the results match
+ do_test <- function (r, params, msg) {
+ nt <- attr(r, 'adist_sumofsquares_nt_model__num')
+ wt <- attr(r, 'adist_sumofsquares_wt_model__num')
+ wtperstock <- attr(r, 'adist_sumofsquares_wtperstock_model__num')
+ m <- params$reader1matrix
+ apply_matrix <- function (destage) {
+ nt[,1,] * m[1,destage] +
+ nt[,2,] * m[2,destage] +
+ nt[,3,] * m[3,destage] +
+ nt[,4,] * m[4,destage] +
+ nt[,5,] * m[5,destage] +
+ 0
+ }
+ ok(ut_cmp_equal(wt, wtperstock), paste0("wt / wtperstock: ", msg))
+ ok(ut_cmp_equal(apply_matrix(1), wt[,1,]), paste0("age1: ", msg))
+ ok(ut_cmp_equal(apply_matrix(2), wt[,2,]), paste0("age2: ", msg))
+ ok(ut_cmp_equal(apply_matrix(3), wt[,3,]), paste0("age3: ", msg))
+ ok(ut_cmp_equal(apply_matrix(4), wt[,4,]), paste0("age4: ", msg))
+ ok(ut_cmp_equal(apply_matrix(5), wt[,5,]), paste0("age5: ", msg))
+ }
+
+ params <- attr(model_fn, 'parameter_template')
+ r <- model_fn(params)
+ do_test(r, params, "Identity matrix")
+ if (Sys.getenv('G3_TEST_TMB') == "2") gadget3:::ut_tmb_r_compare(model_fn, model_tmb, params, model_cpp = model_cpp)
+
+ # Length vector applied for len
+ ok(ut_cmp_equal(
+ attr(r, "adist_sumofsquares_len_model__num"),
+ attr(r, "adist_sumofsquares_nt_model__num") * (10 * g3_stock_def(prey_a, 'midlen')),
+ unused = NULL), "r$adist_sumofsquares_len_model__num: Length vector applied")
+
+ params <- attr(model_fn, 'parameter_template')
+ # Age 1 smeared over bottom 3 groups
+ params$reader1matrix[1,] <- c(0.5, 0.25, 0.25, 0, 0)
+ # Age 2 smeared over 2 groups
+ params$reader1matrix[2,] <- c(0, 0.75, 0.25, 0, 0)
+ r <- model_fn(params)
+ do_test(r, params, "age1, age2 smeared")
+ if (Sys.getenv('G3_TEST_TMB') == "2") gadget3:::ut_tmb_r_compare(model_fn, model_tmb, params, model_cpp = model_cpp)
+ })
# g3l_distribution:transform_fs
ok - wt / wtperstock: Identity matrix
ok - age1: Identity matrix
ok - age2: Identity matrix
ok - age3: Identity matrix
ok - age4: Identity matrix
ok - age5: Identity matrix
ok - r$adist_sumofsquares_len_model__num: Length vector applied
ok - wt / wtperstock: age1, age2 smeared
ok - age1: age1, age2 smeared
ok - age2: age1, age2 smeared
ok - age3: age1, age2 smeared
ok - age4: age1, age2 smeared
ok - age5: age1, age2 smeared
>
> # g3l_distribution_sumofsquaredlogratios
> ok(ut_cmp_identical(
+ deparse1(g3l_distribution_sumofsquaredlogratios()),
+ "~sum((log(stock_ss(obsstock__x) + 10) - log(stock_ss(modelstock__x) + 10))^2)"), "g3l_distribution_sumofsquaredlogratios: Formula code as expected")
ok - g3l_distribution_sumofsquaredlogratios: Formula code as expected
> ok(ut_cmp_identical(
+ length(environment(g3l_distribution_sumofsquaredlogratios())), 0L), "g3l_distribution_sumofsquaredlogratios: Environment empty")
ok - g3l_distribution_sumofsquaredlogratios: Environment empty
>
> # NB: Also added some aggregate areas for fleet data
> areas <- list(a = 1, b = 2, c = 3, x = 1:2, y = 3)
> prey_a <- g3_stock('prey_a', seq(1, 10)) %>% g3s_livesonareas(areas[c('a')])
> prey_b <- g3_stock('prey_b', seq(1, 10)) %>% g3s_livesonareas(areas[c('b')])
> prey_c <- g3_stock('prey_c', seq(1, 10)) %>% g3s_livesonareas(areas[c('c')])
> fleet_abc <- g3_fleet('fleet_abc') %>% g3s_livesonareas(areas[c('a', 'b', 'c')])
>
> # Generate observation data for stock distribution
> # NB: No prey_b, only compare prey_a and prey_c
> sd_data <- expand.grid(year = 1999:2000, step = c(1, 2), area = c('x', 'y'), stock = c("prey_a", "prey_c"), length = c(1,6))
> sd_data$number <- floor(runif(length(sd_data$year), min=100, max=999))
>
> # Generate observation data for catch distribution
> cd_data <- expand.grid(year = 1999:2000, step = c(1, 2), area = c('x', 'y'), length = c(1,6))
> cd_data$number <- floor(runif(length(cd_data$year), min=100, max=999))
>
> # Generate observation data for catch distribution by biomass
> cd_weight_data <- expand.grid(year = 1999:2000, step = c(1, 2), area = c('x', 'y'), length = c(1,6))
> cd_weight_data$weight <- floor(runif(length(cd_data$year), min=1000, max=9999))
>
> # Generate observation data for catch distribution (multinomial)
> multinomial_data <- expand.grid(year = 1999:2000, step = c(1, 2), length = c(1,6))
> multinomial_data$number <- floor(runif(length(multinomial_data$year), min=100, max=999))
>
> surveyindices_data <- expand.grid(year = 1999:2000, step = c(1, 2))
> surveyindices_data$number <- floor(runif(length(surveyindices_data$year), min=100, max=999))
>
> # Can't make catchdistribution without fleets
> ok(ut_cmp_error(g3l_catchdistribution(
+ 'utcd',
+ cd_data,
+ fleets = list(),
+ stocks = list(prey_a, prey_b, prey_c),
+ area_group = areas,
+ g3l_distribution_sumofsquares()), "fleets/predators must be supplied"), "g3l_catchdistribution: Invalid without fleets")
ok - g3l_catchdistribution: Invalid without fleets
> ok(ut_cmp_error(g3l_abundancedistribution(
+ 'utcd',
+ cd_data,
+ fleets = list(fleet_abc),
+ stocks = list(prey_a, prey_b, prey_c),
+ area_group = areas,
+ g3l_distribution_sumofsquares()), "fleets/predators must not be supplied"), "g3l_abundancedistribution: Invalid with fleets")
ok - g3l_abundancedistribution: Invalid with fleets
>
> # Generate a step that reports the value of (var_name) into separate variable (steps) times
> # (initial_val) provides a definition to use to set variable type
> report_step <- function (var_name, steps, initial_val) {
+ out <- ~{}
+ for (i in seq(steps - 1, 0, by = -1)) {
+ step_var_name <- paste0("step", i, "_", var_name)
+ assign(step_var_name, initial_val)
+ out <- gadget3:::f_optimize(gadget3:::f_substitute(~if (cur_time == i) {
+ comment(report_comment)
+ if (is_nll) {
+ # nll is a scalar, and should have attributes stripped
+ step_var <- as.numeric(nll)
+ } else {
+ step_var[] <- var
+ }
+ REPORT(step_var)
+ } else rest, list(
+ report_comment = paste0("Reporting on ", var_name, " at step ", i),
+ is_nll = var_name == 'nll',
+ step_var = as.symbol(step_var_name),
+ var = as.symbol(var_name),
+ i = i,
+ rest = out )))
+ }
+ return(out)
+ }
>
> base_actions <- list(
+ g3a_time(1999,2000, step_lengths = c(6, 6), project_years = 0),
+ g3a_initialconditions(prey_a, ~10 * prey_a__midlen, ~100 * prey_a__midlen),
+ g3a_initialconditions(prey_b, ~10 * prey_b__midlen, ~100 * prey_b__midlen),
+ g3a_initialconditions(prey_c, ~10 * prey_c__midlen, ~100 * prey_c__midlen),
+ g3a_predate_totalfleet(
+ fleet_abc,
+ list(prey_a, prey_b, prey_c),
+ suitabilities = list(
+ prey_a = ~g3_param_vector("fleet_abc_a") + 0 * stock__midlen,
+ prey_b = ~g3_param_vector("fleet_abc_b") + 0 * stock__midlen,
+ prey_c = ~g3_param_vector("fleet_abc_c") + 0 * stock__midlen),
+ amount_f = ~g3_param('amount_ab', value = 100) * area),
+ named_list(
+ # Capture data just before final step erases it
+ gadget3:::step_id(10, 'g3l_distribution', 'cdist_sumofsquares_utcd', 1, 'zzzz'), report_step('cdist_sumofsquares_utcd_model__num', 4, array(
+ # NB: Lift definition from deparse(r$cdist_sumofsquares_utcd_model__num)
+ dim = c(2L, 2L),
+ dimnames = list(
+ c("1:6", "6:10"),
+ c("x", "y")))),
+ gadget3:::step_id(10, 'g3l_distribution', 'cdist_sumofsquares_utsd', 1, 'zzzz'), report_step('cdist_sumofsquares_utsd_model__num', 4, array(
+ # NB: Lift definition from deparse(r$cdist_sumofsquares_utsd_model__num)
+ dim = c(2L, 2L, 2L),
+ dimnames = list(
+ c("1:6", "6:10"),
+ c("prey_a", "prey_c"),
+ c("x", "y")))),
+ gadget3:::step_id(10, 'g3l_distribution', 'cdist_sumofsquares_utcd_weight', 1, 'zzzz'), report_step('cdist_sumofsquares_utcd_weight_model__wgt', 4, array(
+ # NB: Lift definition from deparse(r$cdist_sumofsquares_utcd_weight_model__wgt)
+ dim = c(2L, 2L),
+ dimnames = list(
+ c("1:6", "6:10"),
+ c("x", "y")))),
+ gadget3:::step_id(10, 'g3l_distribution', 'cdist_multinomial_multinom', 1, 'zzzz'), report_step('cdist_multinomial_multinom_model__num', 4, array(
+ # NB: Lift definition from deparse(r$cdist_multinomial_multinom_model__num)
+ dim = c(2L),
+ dimnames = list(
+ c("1:6", "6:10")))),
+ gadget3:::step_id(990, 'prey_a__num'), report_step('prey_a__num', 4, g3_stock_instance(prey_a)),
+ gadget3:::step_id(990, 'prey_b__num'), report_step('prey_b__num', 4, g3_stock_instance(prey_b)),
+ gadget3:::step_id(990, 'prey_c__num'), report_step('prey_c__num', 4, g3_stock_instance(prey_c)),
+ gadget3:::step_id(990, 'nll'), report_step('nll', 4, 0.0),
+ gadget3:::step_id(999), ~{
+ REPORT(cdist_sumofsquares_utcd_model__num)
+ REPORT(cdist_sumofsquares_utcd_obs__num)
+ REPORT(cdist_sumofsquares_utcd_weight_model__wgt)
+ REPORT(cdist_sumofsquares_utcd_weight_obs__wgt)
+ REPORT(cdist_sumofsquares_utsd_model__num)
+ REPORT(cdist_sumofsquares_utsd_obs__num)
+ REPORT(cdist_multinomial_multinom_model__num)
+ REPORT(cdist_multinomial_multinom_obs__num)
+ REPORT(prey_a__wgt)
+ REPORT(prey_b__wgt)
+ REPORT(prey_c__wgt)
+ REPORT(prey_a__num)
+ REPORT(prey_b__num)
+ REPORT(prey_c__num)
+
+ # NB: In theory we could inspect the return value, but TMB doesn't give an easy public method for that
+ REPORT(nll)
+ }))
> actions <- c(base_actions, list(
+ g3l_catchdistribution(
+ 'utcd',
+ cd_data,
+ list(fleet_abc),
+ list(prey_a, prey_b, prey_c),
+ area_group = areas,
+ g3l_distribution_sumofsquares()),
+ g3l_catchdistribution(
+ 'utcd_weight',
+ cd_weight_data,
+ list(fleet_abc),
+ list(prey_a, prey_b, prey_c),
+ area_group = areas,
+ g3l_distribution_sumofsquares()),
+ g3l_catchdistribution(
+ 'utsd',
+ sd_data,
+ list(fleet_abc),
+ list(prey_a, prey_b, prey_c),
+ area_group = areas,
+ g3l_distribution_sumofsquares()),
+ g3l_catchdistribution(
+ 'multinom',
+ multinomial_data,
+ list(fleet_abc),
+ list(prey_a),
+ area_group = areas,
+ g3l_distribution_multinomial()),
+ g3l_abundancedistribution(
+ 'surveyindices',
+ surveyindices_data,
+ fleets = list(),
+ stocks = list(prey_b),
+ area_group = areas,
+ report = TRUE, # NB: Using built-in reporting vs. version hacked in tests
+ g3l_distribution_surveyindices_log(alpha = ~g3_param("si_alpha", value = 0), beta = ~g3_param("si_beta", value = 0))),
+ NULL))
> actions <- c(actions, list(
+ # NB: Don't use _detail, since it doesn't play well with report_step()
+ # Ideally we'd remove report_step, but we have to time it exactly between data-collect and reset
+ g3a_report_history(actions, '__cons$', out_prefix = "detail_") ))
>
> # Compile model
> model_fn <- g3_to_r(actions, trace = FALSE)
> # model_fn <- edit(model_fn)
> if (nzchar(Sys.getenv('G3_TEST_TMB'))) {
+ params <- attr(model_fn, 'parameter_template')
+ params$fleet_abc_a <- c(0, 0, 0, 0.1, 0.2, 0.1, 0, 0, 0, 0)
+ params$fleet_abc_b <- c(0, 0, 0, 0, 0, 0.1, 0.2, 0.1, 0, 0)
+ params$fleet_abc_c <- c(0, 0, 0, 0, 0, 0, 0.1, 0.2, 0.1, 0)
+
+ model_cpp <- g3_to_tmb(actions, trace = FALSE)
+ # model_cpp <- edit(model_cpp)
+ model_tmb <- g3_tmb_adfun(model_cpp, params, compile_flags = c("-O0", "-g"))
+ } else {
+ writeLines("# skip: not compiling TMB model")
+ model_cpp <- c()
+ }
# skip: not compiling TMB model
>
> ok_group("Likelihood per step", {
+ params <- attr(model_fn, 'parameter_template')
+ # Randomly catch, but get something everywhere
+ params$fleet_abc_a <- runif(10, min=0.1, max=0.9)
+ params$fleet_abc_b <- runif(10, min=0.1, max=0.9)
+ params$fleet_abc_c <- runif(10, min=0.1, max=0.9)
+ params$amount_ab <- 1000000
+ params$si_alpha <- 4
+ params$si_beta <- 2
+ params$cdist_sumofsquares_utcd_weight <- 1
+ params$cdist_sumofsquares_utcd_weight_weight <- 1
+ params$cdist_sumofsquares_utsd_weight <- 1
+ params$cdist_multinomial_multinom_weight <- 1
+ params$adist_surveyindices_log_surveyindices_weight <- 1
+
+ result <- model_fn(params)
+ r <- attributes(result)
+ # str(result)
+ # str(as.list(r), vec.len = 10000)
+
+ ok(ut_cmp_equal(
+ sort(as.vector(r$cdist_sumofsquares_utsd_obs__num)),
+ sort(sd_data$number)), "cdist_sumofsquares_utsd_obs__num: Imported from data.frame, order not necessarily the same")
+
+ ok(ut_cmp_equal(
+ r$adist_surveyindices_log_surveyindices_model__params,
+ c(params$si_alpha, params$si_beta)), "adist_surveyindices_log_surveyindices_model__params: Reported our hard-coded linear regression parameters")
+
+ ######## cdist_sumofsquares_utsd_model__num
+ ok(ut_cmp_equal(as.vector(r$step0_cdist_sumofsquares_utsd_model__num[,'prey_a', 1]), c(
+ sum(r$detail_prey_a_fleet_abc__cons[1:5,,1] / r$prey_a__wgt[1:5,]),
+ sum(r$detail_prey_a_fleet_abc__cons[6:10,,1] / r$prey_a__wgt[6:10,]))), "step0_cdist_sumofsquares_utsd_model__num[,'prey_a',1]: prey_a in area 1")
+ ok(ut_cmp_equal(as.vector(r$step0_cdist_sumofsquares_utsd_model__num[,'prey_c', 1]), c(
+ 0,
+ 0)), "step0_cdist_sumofsquares_utsd_model__num[,'prey_c',1]: No prey_c in area 1")
+ ok(ut_cmp_equal(as.vector(r$step0_cdist_sumofsquares_utsd_model__num[,'prey_a', 2]), c(
+ 0,
+ 0)), "step0_cdist_sumofsquares_utsd_model__num[,'prey_a',2]: No prey_a in area 2")
+ ok(ut_cmp_equal(as.vector(r$step0_cdist_sumofsquares_utsd_model__num[,'prey_c', 2]), c(
+ sum(r$detail_prey_c_fleet_abc__cons[1:5,,1] / r$prey_c__wgt[1:5,]),
+ sum(r$detail_prey_c_fleet_abc__cons[6:10,,1] / r$prey_c__wgt[6:10,]))), "step0_cdist_sumofsquares_utsd_model__num[,'prey_c',2]: prey_c in area 2")
+
+ ok(ut_cmp_equal(as.vector(r$step1_cdist_sumofsquares_utsd_model__num[,'prey_a', 1]), c(
+ sum(r$detail_prey_a_fleet_abc__cons[1:5,,2] / r$prey_a__wgt[1:5,]),
+ sum(r$detail_prey_a_fleet_abc__cons[6:10,,2] / r$prey_a__wgt[6:10,]))), "step1_cdist_sumofsquares_utsd_model__num[,'prey_a',1]: prey_a in area 1")
+ ok(ut_cmp_equal(as.vector(r$step1_cdist_sumofsquares_utsd_model__num[,'prey_c', 1]), c(
+ 0,
+ 0)), "step1_cdist_sumofsquares_utsd_model__num[,'prey_c',1]: No prey_c in area 1")
+ ok(ut_cmp_equal(as.vector(r$step1_cdist_sumofsquares_utsd_model__num[,'prey_a', 2]), c(
+ 0,
+ 0)), "step1_cdist_sumofsquares_utsd_model__num[,'prey_a',2]: No prey_a in area 2")
+ ok(ut_cmp_equal(as.vector(r$step1_cdist_sumofsquares_utsd_model__num[,'prey_c', 2]), c(
+ sum(r$detail_prey_c_fleet_abc__cons[1:5,,2] / r$prey_c__wgt[1:5,]),
+ sum(r$detail_prey_c_fleet_abc__cons[6:10,,2] / r$prey_c__wgt[6:10,]))), "step1_cdist_sumofsquares_utsd_model__num[,'prey_c',2]: prey_c in area 2")
+
+ ok(ut_cmp_equal(as.vector(r$step2_cdist_sumofsquares_utsd_model__num[,'prey_a', 1]), c(
+ sum(r$detail_prey_a_fleet_abc__cons[1:5,,3] / r$prey_a__wgt[1:5,]),
+ sum(r$detail_prey_a_fleet_abc__cons[6:10,,3] / r$prey_a__wgt[6:10,]))), "step2_cdist_sumofsquares_utsd_model__num[,'prey_a',1]: prey_a in area 1")
+ ok(ut_cmp_equal(as.vector(r$step2_cdist_sumofsquares_utsd_model__num[,'prey_c', 1]), c(
+ 0,
+ 0)), "step2_cdist_sumofsquares_utsd_model__num[,'prey_c',1]: No prey_c in area 1")
+ ok(ut_cmp_equal(as.vector(r$step2_cdist_sumofsquares_utsd_model__num[,'prey_a', 2]), c(
+ 0,
+ 0)), "step2_cdist_sumofsquares_utsd_model__num[,'prey_a',2]: No prey_a in area 2")
+ ok(ut_cmp_equal(as.vector(r$step2_cdist_sumofsquares_utsd_model__num[,'prey_c', 2]), c(
+ sum(r$detail_prey_c_fleet_abc__cons[1:5,,3] / r$prey_c__wgt[1:5,]),
+ sum(r$detail_prey_c_fleet_abc__cons[6:10,,3] / r$prey_c__wgt[6:10,]))), "step2_cdist_sumofsquares_utsd_model__num[,'prey_c',2]: prey_c in area 2")
+
+ ok(ut_cmp_equal(as.vector(r$step3_cdist_sumofsquares_utsd_model__num[,'prey_a', 1]), c(
+ sum(r$detail_prey_a_fleet_abc__cons[1:5,,4] / r$prey_a__wgt[1:5,]),
+ sum(r$detail_prey_a_fleet_abc__cons[6:10,,4] / r$prey_a__wgt[6:10,]))), "step3_cdist_sumofsquares_utsd_model__num[,'prey_a',1]: prey_a in area 1")
+ ok(ut_cmp_equal(as.vector(r$step3_cdist_sumofsquares_utsd_model__num[,'prey_c', 1]), c(
+ 0,
+ 0)), "step3_cdist_sumofsquares_utsd_model__num[,'prey_c',1]: No prey_c in area 1")
+ ok(ut_cmp_equal(as.vector(r$step3_cdist_sumofsquares_utsd_model__num[,'prey_a', 2]), c(
+ 0,
+ 0)), "step3_cdist_sumofsquares_utsd_model__num[,'prey_a',2]: No prey_a in area 2")
+ ok(ut_cmp_equal(as.vector(r$step3_cdist_sumofsquares_utsd_model__num[,'prey_c', 2]), c(
+ sum(r$detail_prey_c_fleet_abc__cons[1:5,,4] / r$prey_c__wgt[1:5,]),
+ sum(r$detail_prey_c_fleet_abc__cons[6:10,,4] / r$prey_c__wgt[6:10,]))), "step3_cdist_sumofsquares_utsd_model__num[,'prey_c',2]: prey_c in area 2")
+ ########
+
+ ######## cdist_sumofsquares_utcd_model__num
+ ok(ut_cmp_equal(as.vector(r$step0_cdist_sumofsquares_utcd_model__num[,1]), c(
+ sum(
+ sum(r$detail_prey_a_fleet_abc__cons[1:5,1,1] / r$prey_a__wgt[1:5,1]),
+ sum(r$detail_prey_b_fleet_abc__cons[1:5,1,1] / r$prey_b__wgt[1:5,1])),
+ sum(
+ sum(r$detail_prey_a_fleet_abc__cons[6:10,1,1] / r$prey_a__wgt[6:10,1]),
+ sum(r$detail_prey_b_fleet_abc__cons[6:10,1,1] / r$prey_b__wgt[6:10,1])),
+ NULL)), "step0_cdist_sumofsquares_utsd_model__num[,1]: all prey in area a/b")
+ ok(ut_cmp_equal(as.vector(r$step0_cdist_sumofsquares_utcd_model__num[,2]), c(
+ sum(
+ sum(r$detail_prey_c_fleet_abc__cons[1:5,1,1] / r$prey_c__wgt[1:5,1])),
+ sum(
+ sum(r$detail_prey_c_fleet_abc__cons[6:10,1,1] / r$prey_c__wgt[6:10,1])),
+ NULL)), "step0_cdist_sumofsquares_utsd_model__num[,2]: all prey in area c")
+
+ ok(ut_cmp_equal(as.vector(r$step1_cdist_sumofsquares_utcd_model__num[,1]), c(
+ sum(
+ sum(r$detail_prey_a_fleet_abc__cons[1:5,1,2] / r$prey_a__wgt[1:5,1]),
+ sum(r$detail_prey_b_fleet_abc__cons[1:5,1,2] / r$prey_b__wgt[1:5,1])),
+ sum(
+ sum(r$detail_prey_a_fleet_abc__cons[6:10,1,2] / r$prey_a__wgt[6:10,1]),
+ sum(r$detail_prey_b_fleet_abc__cons[6:10,1,2] / r$prey_b__wgt[6:10,1])),
+ NULL)), "step1_cdist_sumofsquares_utsd_model__num[,1]: all prey in area a/b")
+ ok(ut_cmp_equal(as.vector(r$step1_cdist_sumofsquares_utcd_model__num[,2]), c(
+ sum(
+ sum(r$detail_prey_c_fleet_abc__cons[1:5,1,2] / r$prey_c__wgt[1:5,1])),
+ sum(
+ sum(r$detail_prey_c_fleet_abc__cons[6:10,1,2] / r$prey_c__wgt[6:10,1])),
+ NULL)), "step1_cdist_sumofsquares_utsd_model__num[,2]: all prey in area c")
+
+ ok(ut_cmp_equal(as.vector(r$step2_cdist_sumofsquares_utcd_model__num[,1]), c(
+ sum(
+ sum(r$detail_prey_a_fleet_abc__cons[1:5,1,3] / r$prey_a__wgt[1:5,1]),
+ sum(r$detail_prey_b_fleet_abc__cons[1:5,1,3] / r$prey_b__wgt[1:5,1])),
+ sum(
+ sum(r$detail_prey_a_fleet_abc__cons[6:10,1,3] / r$prey_a__wgt[6:10,1]),
+ sum(r$detail_prey_b_fleet_abc__cons[6:10,1,3] / r$prey_b__wgt[6:10,1])),
+ NULL)), "step2_cdist_sumofsquares_utsd_model__num[,1]: all prey in area a/b")
+ ok(ut_cmp_equal(as.vector(r$step2_cdist_sumofsquares_utcd_model__num[,2]), c(
+ sum(
+ sum(r$detail_prey_c_fleet_abc__cons[1:5,1,3] / r$prey_c__wgt[1:5,1])),
+ sum(
+ sum(r$detail_prey_c_fleet_abc__cons[6:10,1,3] / r$prey_c__wgt[6:10,1])),
+ NULL)), "step2_cdist_sumofsquares_utsd_model__num[,2]: all prey in area c")
+
+ ok(ut_cmp_equal(as.vector(r$step3_cdist_sumofsquares_utcd_model__num[,1]), c(
+ sum(
+ sum(r$detail_prey_a_fleet_abc__cons[1:5,1,4] / r$prey_a__wgt[1:5,1]),
+ sum(r$detail_prey_b_fleet_abc__cons[1:5,1,4] / r$prey_b__wgt[1:5,1])),
+ sum(
+ sum(r$detail_prey_a_fleet_abc__cons[6:10,1,4] / r$prey_a__wgt[6:10,1]),
+ sum(r$detail_prey_b_fleet_abc__cons[6:10,1,4] / r$prey_b__wgt[6:10,1])),
+ NULL)), "step3_cdist_sumofsquares_utsd_model__num[,1]: all prey in area a/b")
+ ok(ut_cmp_equal(as.vector(r$step3_cdist_sumofsquares_utcd_model__num[,2]), c(
+ sum(
+ sum(r$detail_prey_c_fleet_abc__cons[1:5,1,4] / r$prey_c__wgt[1:5,1])),
+ sum(
+ sum(r$detail_prey_c_fleet_abc__cons[6:10,1,4] / r$prey_c__wgt[6:10,1])),
+ NULL)), "step3_cdist_sumofsquares_utsd_model__num[,2]: all prey in area c")
+ ########
+
+ ######## adist_surveyindices_log_surveyindices_model__num
+ ok(ut_cmp_equal(
+ as.vector(r$adist_surveyindices_log_surveyindices_model__num),
+ c(
+ sum(r$step0_prey_b__num),
+ sum(r$step1_prey_b__num),
+ sum(r$step2_prey_b__num),
+ sum(r$step3_prey_b__num),
+ NULL)), "adist_surveyindices_log_surveyindices_model__num: Built-in reporting gave us step 0..3 abundance")
+ ########
+
+ ok(ut_cmp_equal(r$step0_nll, sum(
+ # utsd: stock 1 / area 1
+ (r$step0_cdist_sumofsquares_utsd_model__num[,1,1] / sum(r$step0_cdist_sumofsquares_utsd_model__num[,,1]) -
+ r$cdist_sumofsquares_utsd_obs__num[,1,1,1] / sum(r$cdist_sumofsquares_utsd_obs__num[,,1,1])) ** 2,
+ # utsd: stock 2 / area 1
+ (r$step0_cdist_sumofsquares_utsd_model__num[,2,1] / sum(r$step0_cdist_sumofsquares_utsd_model__num[,,1]) -
+ r$cdist_sumofsquares_utsd_obs__num[,2,1,1] / sum(r$cdist_sumofsquares_utsd_obs__num[,,1,1])) ** 2,
+ # utsd: stock 1 / area 2
+ (r$step0_cdist_sumofsquares_utsd_model__num[,1,2] / sum(r$step0_cdist_sumofsquares_utsd_model__num[,,2]) -
+ r$cdist_sumofsquares_utsd_obs__num[,1,1,2] / sum(r$cdist_sumofsquares_utsd_obs__num[,,1,2])) ** 2,
+ # utsd: stock 2 / area 2
+ (r$step0_cdist_sumofsquares_utsd_model__num[,2,2] / sum(r$step0_cdist_sumofsquares_utsd_model__num[,,2]) -
+ r$cdist_sumofsquares_utsd_obs__num[,2,1,2] / sum(r$cdist_sumofsquares_utsd_obs__num[,,1,2])) ** 2,
+ # utcd: area 1
+ (r$step0_cdist_sumofsquares_utcd_model__num[,1] / sum(r$step0_cdist_sumofsquares_utcd_model__num[,1]) -
+ r$cdist_sumofsquares_utcd_obs__num[,1,1] / sum(r$cdist_sumofsquares_utcd_obs__num[,1,1])) ** 2,
+ # utcd: area 2
+ (r$step0_cdist_sumofsquares_utcd_model__num[,2] / sum(r$step0_cdist_sumofsquares_utcd_model__num[,2]) -
+ r$cdist_sumofsquares_utcd_obs__num[,1,2] / sum(r$cdist_sumofsquares_utcd_obs__num[,1,2])) ** 2,
+ # utcd_weight: area 1
+ (r$step0_cdist_sumofsquares_utcd_weight_model__wgt[,1] / g3_avoid_zero(sum(r$step0_cdist_sumofsquares_utcd_weight_model__wgt[,1])) -
+ r$cdist_sumofsquares_utcd_weight_obs__wgt[,1,1] / g3_avoid_zero(sum(r$cdist_sumofsquares_utcd_weight_obs__wgt[,1,1]))) ** 2,
+ # utcd_weight: area 2
+ (r$step0_cdist_sumofsquares_utcd_weight_model__wgt[,2] / g3_avoid_zero(sum(r$step0_cdist_sumofsquares_utcd_weight_model__wgt[,2])) -
+ r$cdist_sumofsquares_utcd_weight_obs__wgt[,1,2] / g3_avoid_zero(sum(r$cdist_sumofsquares_utcd_weight_obs__wgt[,1,2]))) ** 2,
+ # multinom:
+ (2 * (-sum(r$cdist_multinomial_multinom_obs__num[,1] * log(g3_logspace_add(r$step0_cdist_multinomial_multinom_model__num/g3_avoid_zero(sum(r$step0_cdist_multinomial_multinom_model__num)) *
+ 10000, (1/(length(r$cdist_multinomial_multinom_obs__num[,1]) * 10)) * 10000)/10000)) +
+ (sum(g3_lgamma_vec(1 + r$cdist_multinomial_multinom_obs__num[,1])) - lgamma(1 +
+ sum(r$cdist_multinomial_multinom_obs__num[,1]))))),
+ # surveyindices:
+ 0, # NB: We don't calculate until end
+ 0)), "step0_nll: Sum of squares")
+
+ ok(ut_cmp_equal(r$step1_nll, sum(
+ # utsd: stock 1 / area 1
+ (r$step1_cdist_sumofsquares_utsd_model__num[,1,1] / g3_avoid_zero(sum(r$step1_cdist_sumofsquares_utsd_model__num[,,1])) -
+ r$cdist_sumofsquares_utsd_obs__num[,1,2,1] / g3_avoid_zero(sum(r$cdist_sumofsquares_utsd_obs__num[,,2,1]))) ** 2,
+ # utsd: stock 2 / area 1
+ (r$step1_cdist_sumofsquares_utsd_model__num[,2,1] / g3_avoid_zero(sum(r$step1_cdist_sumofsquares_utsd_model__num[,,1])) -
+ r$cdist_sumofsquares_utsd_obs__num[,2,2,1] / g3_avoid_zero(sum(r$cdist_sumofsquares_utsd_obs__num[,,2,1]))) ** 2,
+ # utsd: stock 1 / area 2
+ (r$step1_cdist_sumofsquares_utsd_model__num[,1,2] / g3_avoid_zero(sum(r$step1_cdist_sumofsquares_utsd_model__num[,,2])) -
+ r$cdist_sumofsquares_utsd_obs__num[,1,2,2] / g3_avoid_zero(sum(r$cdist_sumofsquares_utsd_obs__num[,,2,2]))) ** 2,
+ # utsd: stock 2 / area 2
+ (r$step1_cdist_sumofsquares_utsd_model__num[,2,2] / g3_avoid_zero(sum(r$step1_cdist_sumofsquares_utsd_model__num[,,2])) -
+ r$cdist_sumofsquares_utsd_obs__num[,2,2,2] / g3_avoid_zero(sum(r$cdist_sumofsquares_utsd_obs__num[,,2,2]))) ** 2,
+ # utcd: area 1
+ (r$step1_cdist_sumofsquares_utcd_model__num[,1] / g3_avoid_zero(sum(r$step1_cdist_sumofsquares_utcd_model__num[,1])) -
+ r$cdist_sumofsquares_utcd_obs__num[,2,1] / g3_avoid_zero(sum(r$cdist_sumofsquares_utcd_obs__num[,2,1]))) ** 2,
+ # utcd: area 2
+ (r$step1_cdist_sumofsquares_utcd_model__num[,2] / g3_avoid_zero(sum(r$step1_cdist_sumofsquares_utcd_model__num[,2])) -
+ r$cdist_sumofsquares_utcd_obs__num[,2,2] / g3_avoid_zero(sum(r$cdist_sumofsquares_utcd_obs__num[,2,2]))) ** 2,
+ # utcd_weight: area 1
+ (r$step1_cdist_sumofsquares_utcd_weight_model__wgt[,1] / g3_avoid_zero(sum(r$step1_cdist_sumofsquares_utcd_weight_model__wgt[,1])) -
+ r$cdist_sumofsquares_utcd_weight_obs__wgt[,2,1] / g3_avoid_zero(sum(r$cdist_sumofsquares_utcd_weight_obs__wgt[,2,1]))) ** 2,
+ # utcd_weight: area 2
+ (r$step1_cdist_sumofsquares_utcd_weight_model__wgt[,2] / g3_avoid_zero(sum(r$step1_cdist_sumofsquares_utcd_weight_model__wgt[,2])) -
+ r$cdist_sumofsquares_utcd_weight_obs__wgt[,2,2] / g3_avoid_zero(sum(r$cdist_sumofsquares_utcd_weight_obs__wgt[,2,2]))) ** 2,
+ # multinom:
+ (2 * (-sum(r$cdist_multinomial_multinom_obs__num[,2] * log(g3_logspace_add(r$step1_cdist_multinomial_multinom_model__num/g3_avoid_zero(sum(r$step1_cdist_multinomial_multinom_model__num)) *
+ 10000, (1/(length(r$cdist_multinomial_multinom_obs__num[,2]) * 10)) * 10000)/10000)) +
+ (sum(g3_lgamma_vec(1 + r$cdist_multinomial_multinom_obs__num[,2])) - lgamma(1 +
+ sum(r$cdist_multinomial_multinom_obs__num[,2]))))),
+ # surveyindices:
+ 0, # NB: We don't calculate until end
+ r$step0_nll)), "step1_nll: Sum of squares, including step0_nll")
+
+ ok(ut_cmp_equal(r$step2_nll, sum(
+ # utsd: stock 1 / area 1
+ (r$step2_cdist_sumofsquares_utsd_model__num[,1,1] / g3_avoid_zero(sum(r$step2_cdist_sumofsquares_utsd_model__num[,,1])) -
+ r$cdist_sumofsquares_utsd_obs__num[,1,3,1] / g3_avoid_zero(sum(r$cdist_sumofsquares_utsd_obs__num[,,3,1]))) ** 2,
+ # utsd: stock 2 / area 1
+ (r$step2_cdist_sumofsquares_utsd_model__num[,2,1] / g3_avoid_zero(sum(r$step2_cdist_sumofsquares_utsd_model__num[,,1])) -
+ r$cdist_sumofsquares_utsd_obs__num[,2,3,1] / g3_avoid_zero(sum(r$cdist_sumofsquares_utsd_obs__num[,,3,1]))) ** 2,
+ # utsd: stock 1 / area 2
+ (r$step2_cdist_sumofsquares_utsd_model__num[,1,2] / g3_avoid_zero(sum(r$step2_cdist_sumofsquares_utsd_model__num[,,2])) -
+ r$cdist_sumofsquares_utsd_obs__num[,1,3,2] / g3_avoid_zero(sum(r$cdist_sumofsquares_utsd_obs__num[,,3,2]))) ** 2,
+ # utsd: stock 2 / area 2
+ (r$step2_cdist_sumofsquares_utsd_model__num[,2,2] / g3_avoid_zero(sum(r$step2_cdist_sumofsquares_utsd_model__num[,,2])) -
+ r$cdist_sumofsquares_utsd_obs__num[,2,3,2] / g3_avoid_zero(sum(r$cdist_sumofsquares_utsd_obs__num[,,3,2]))) ** 2,
+ # utcd: area 1
+ (r$step2_cdist_sumofsquares_utcd_model__num[,1] / g3_avoid_zero(sum(r$step2_cdist_sumofsquares_utcd_model__num[,1])) -
+ r$cdist_sumofsquares_utcd_obs__num[,3,1] / g3_avoid_zero(sum(r$cdist_sumofsquares_utcd_obs__num[,3,1]))) ** 2,
+ # utcd: area 2
+ (r$step2_cdist_sumofsquares_utcd_model__num[,2] / g3_avoid_zero(sum(r$step2_cdist_sumofsquares_utcd_model__num[,2])) -
+ r$cdist_sumofsquares_utcd_obs__num[,3,2] / g3_avoid_zero(sum(r$cdist_sumofsquares_utcd_obs__num[,3,2]))) ** 2,
+ # utcd_weight: area 1
+ (r$step2_cdist_sumofsquares_utcd_weight_model__wgt[,1] / g3_avoid_zero(sum(r$step2_cdist_sumofsquares_utcd_weight_model__wgt[,1])) -
+ r$cdist_sumofsquares_utcd_weight_obs__wgt[,3,1] / g3_avoid_zero(sum(r$cdist_sumofsquares_utcd_weight_obs__wgt[,3,1]))) ** 2,
+ # utcd_weight: area 2
+ (r$step2_cdist_sumofsquares_utcd_weight_model__wgt[,2] / g3_avoid_zero(sum(r$step2_cdist_sumofsquares_utcd_weight_model__wgt[,2])) -
+ r$cdist_sumofsquares_utcd_weight_obs__wgt[,3,2] / g3_avoid_zero(sum(r$cdist_sumofsquares_utcd_weight_obs__wgt[,3,2]))) ** 2,
+ # multinom:
+ (2 * (-sum(r$cdist_multinomial_multinom_obs__num[,3] * log(g3_logspace_add(r$step2_cdist_multinomial_multinom_model__num/g3_avoid_zero(sum(r$step2_cdist_multinomial_multinom_model__num)) *
+ 10000, (1/(length(r$cdist_multinomial_multinom_obs__num[,3]) * 10)) * 10000)/10000)) +
+ (sum(g3_lgamma_vec(1 + r$cdist_multinomial_multinom_obs__num[,3])) - lgamma(1 +
+ sum(r$cdist_multinomial_multinom_obs__num[,3]))))),
+ # surveyindices:
+ 0, # NB: We don't calculate until end
+ r$step1_nll)), "step2_nll: Sum of squares, including step1_nll")
+
+ ok(ut_cmp_equal(r$step3_nll, sum(
+ # utsd: stock 1 / area 1
+ (r$step3_cdist_sumofsquares_utsd_model__num[,1,1] / g3_avoid_zero(sum(r$step3_cdist_sumofsquares_utsd_model__num[,,1])) -
+ r$cdist_sumofsquares_utsd_obs__num[,1,4,1] / g3_avoid_zero(sum(r$cdist_sumofsquares_utsd_obs__num[,,4,1]))) ** 2,
+ # utsd: stock 2 / area 1
+ (r$step3_cdist_sumofsquares_utsd_model__num[,2,1] / g3_avoid_zero(sum(r$step3_cdist_sumofsquares_utsd_model__num[,,1])) -
+ r$cdist_sumofsquares_utsd_obs__num[,2,4,1] / g3_avoid_zero(sum(r$cdist_sumofsquares_utsd_obs__num[,,4,1]))) ** 2,
+ # utsd: stock 1 / area 2
+ (r$step3_cdist_sumofsquares_utsd_model__num[,1,2] / g3_avoid_zero(sum(r$step3_cdist_sumofsquares_utsd_model__num[,,2])) -
+ r$cdist_sumofsquares_utsd_obs__num[,1,4,2] / g3_avoid_zero(sum(r$cdist_sumofsquares_utsd_obs__num[,,4,2]))) ** 2,
+ # utsd: stock 2 / area 2
+ (r$step3_cdist_sumofsquares_utsd_model__num[,2,2] / g3_avoid_zero(sum(r$step3_cdist_sumofsquares_utsd_model__num[,,2])) -
+ r$cdist_sumofsquares_utsd_obs__num[,2,4,2] / g3_avoid_zero(sum(r$cdist_sumofsquares_utsd_obs__num[,,4,2]))) ** 2,
+ # utcd: area 1
+ (r$step3_cdist_sumofsquares_utcd_model__num[,1] / g3_avoid_zero(sum(r$step3_cdist_sumofsquares_utcd_model__num[,1])) -
+ r$cdist_sumofsquares_utcd_obs__num[,4,1] / g3_avoid_zero(sum(r$cdist_sumofsquares_utcd_obs__num[,4,1]))) ** 2,
+ # utcd: area 2
+ (r$step3_cdist_sumofsquares_utcd_model__num[,2] / g3_avoid_zero(sum(r$step3_cdist_sumofsquares_utcd_model__num[,2])) -
+ r$cdist_sumofsquares_utcd_obs__num[,4,2] / g3_avoid_zero(sum(r$cdist_sumofsquares_utcd_obs__num[,4,2]))) ** 2,
+ # utcd_weight: area 1
+ (r$step3_cdist_sumofsquares_utcd_weight_model__wgt[,1] / g3_avoid_zero(sum(r$step3_cdist_sumofsquares_utcd_weight_model__wgt[,1])) -
+ r$cdist_sumofsquares_utcd_weight_obs__wgt[,4,1] / g3_avoid_zero(sum(r$cdist_sumofsquares_utcd_weight_obs__wgt[,4,1]))) ** 2,
+ # utcd_weight: area 2
+ (r$step3_cdist_sumofsquares_utcd_weight_model__wgt[,2] / g3_avoid_zero(sum(r$step3_cdist_sumofsquares_utcd_weight_model__wgt[,2])) -
+ r$cdist_sumofsquares_utcd_weight_obs__wgt[,4,2] / g3_avoid_zero(sum(r$cdist_sumofsquares_utcd_weight_obs__wgt[,4,2]))) ** 2,
+ # multinom:
+ (2 * (-sum(r$cdist_multinomial_multinom_obs__num[,4] * log(g3_logspace_add(r$step3_cdist_multinomial_multinom_model__num/g3_avoid_zero(sum(r$step3_cdist_multinomial_multinom_model__num)) *
+ 10000, (1/(length(r$cdist_multinomial_multinom_obs__num[,4]) * 10)) * 10000)/10000)) +
+ (sum(g3_lgamma_vec(1 + r$cdist_multinomial_multinom_obs__num[,4])) - lgamma(1 +
+ sum(r$cdist_multinomial_multinom_obs__num[,4]))))),
+ # surveyindices:
+ sum((params$si_alpha +
+ params$si_beta * log(g3_avoid_zero(r$adist_surveyindices_log_surveyindices_model__num[,])) -
+ log(g3_avoid_zero(r$adist_surveyindices_log_surveyindices_obs__num[,])))**2),
+ r$step2_nll)), "step3_nll: Sum of squares, including step2_nll")
+
+ 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)
+ }
+ })
# Likelihood per step
ok - cdist_sumofsquares_utsd_obs__num: Imported from data.frame, order not necessarily the same
ok - adist_surveyindices_log_surveyindices_model__params: Reported our hard-coded linear regression parameters
ok - step0_cdist_sumofsquares_utsd_model__num[,'prey_a',1]: prey_a in area 1
ok - step0_cdist_sumofsquares_utsd_model__num[,'prey_c',1]: No prey_c in area 1
ok - step0_cdist_sumofsquares_utsd_model__num[,'prey_a',2]: No prey_a in area 2
ok - step0_cdist_sumofsquares_utsd_model__num[,'prey_c',2]: prey_c in area 2
ok - step1_cdist_sumofsquares_utsd_model__num[,'prey_a',1]: prey_a in area 1
ok - step1_cdist_sumofsquares_utsd_model__num[,'prey_c',1]: No prey_c in area 1
ok - step1_cdist_sumofsquares_utsd_model__num[,'prey_a',2]: No prey_a in area 2
ok - step1_cdist_sumofsquares_utsd_model__num[,'prey_c',2]: prey_c in area 2
ok - step2_cdist_sumofsquares_utsd_model__num[,'prey_a',1]: prey_a in area 1
ok - step2_cdist_sumofsquares_utsd_model__num[,'prey_c',1]: No prey_c in area 1
ok - step2_cdist_sumofsquares_utsd_model__num[,'prey_a',2]: No prey_a in area 2
ok - step2_cdist_sumofsquares_utsd_model__num[,'prey_c',2]: prey_c in area 2
ok - step3_cdist_sumofsquares_utsd_model__num[,'prey_a',1]: prey_a in area 1
ok - step3_cdist_sumofsquares_utsd_model__num[,'prey_c',1]: No prey_c in area 1
ok - step3_cdist_sumofsquares_utsd_model__num[,'prey_a',2]: No prey_a in area 2
ok - step3_cdist_sumofsquares_utsd_model__num[,'prey_c',2]: prey_c in area 2
ok - step0_cdist_sumofsquares_utsd_model__num[,1]: all prey in area a/b
ok - step0_cdist_sumofsquares_utsd_model__num[,2]: all prey in area c
ok - step1_cdist_sumofsquares_utsd_model__num[,1]: all prey in area a/b
ok - step1_cdist_sumofsquares_utsd_model__num[,2]: all prey in area c
ok - step2_cdist_sumofsquares_utsd_model__num[,1]: all prey in area a/b
ok - step2_cdist_sumofsquares_utsd_model__num[,2]: all prey in area c
ok - step3_cdist_sumofsquares_utsd_model__num[,1]: all prey in area a/b
ok - step3_cdist_sumofsquares_utsd_model__num[,2]: all prey in area c
ok - adist_surveyindices_log_surveyindices_model__num: Built-in reporting gave us step 0..3 abundance
ok - step0_nll: Sum of squares
ok - step1_nll: Sum of squares, including step0_nll
ok - step2_nll: Sum of squares, including step1_nll
ok - step3_nll: Sum of squares, including step2_nll
>
> ok_group("Likelihood per year", {
+ # Change model to aggregate over a year
+ # NB: No prey_b, only compare prey_a and prey_c
+ # NB: No step column now
+ sd_data <- expand.grid(year = 1999:2000, area = c('x', 'y'), stock = c("prey_a", "prey_c"), length = c(1,6))
+ sd_data$number <- floor(runif(length(sd_data$year), min=100, max=999))
+ attr(sd_data, 'area') <- list(
+ x = areas[c('a', 'b')],
+ y = areas[c('c')])
+
+ # Change model to aggregate over a year
+ # NB: No step column now
+ cd_data <- expand.grid(year = 1999:2000, area = c('x', 'y'), length = c(1,6))
+ cd_data$number <- floor(runif(length(cd_data$year), min=100, max=999))
+ attr(cd_data, 'area') <- list(
+ x = areas[c('a', 'b')],
+ y = areas[c('c')])
+
+ # Change model to aggregate over a year
+ cd_weight_data <- expand.grid(year = 1999:2000, area = c('x', 'y'), length = c(1,6))
+ cd_weight_data$weight <- floor(runif(length(cd_data$year), min=1000, max=9999))
+
+ # Change model to aggregate over a year
+ # NB: No step column now
+ multinomial_data <- expand.grid(year = 1999:2000, length = c(1,6))
+ multinomial_data$number <- floor(runif(length(multinomial_data$year), min=100, max=999))
+
+ # NB: No step column now
+ surveyindices_data <- expand.grid(year = 1999:2000)
+ surveyindices_data$number <- floor(runif(length(surveyindices_data$year), min=100, max=999))
+
+ actions <- c(base_actions, list(
+ g3l_catchdistribution(
+ 'utcd',
+ cd_data,
+ list(fleet_abc),
+ list(prey_a, prey_b, prey_c),
+ area_group = areas,
+ g3l_distribution_sumofsquares()),
+ g3l_catchdistribution(
+ 'utcd_weight',
+ cd_weight_data,
+ list(fleet_abc),
+ list(prey_a, prey_b, prey_c),
+ area_group = areas,
+ g3l_distribution_sumofsquares()),
+ g3l_catchdistribution(
+ 'utsd',
+ sd_data,
+ list(fleet_abc),
+ list(prey_a, prey_b, prey_c),
+ area_group = areas,
+ g3l_distribution_sumofsquares()),
+ g3l_catchdistribution(
+ 'multinom',
+ multinomial_data,
+ list(fleet_abc),
+ list(prey_a),
+ area_group = areas,
+ g3l_distribution_multinomial()),
+ g3l_abundancedistribution(
+ 'surveyindices',
+ surveyindices_data,
+ fleets = list(),
+ stocks = list(prey_b),
+ area_group = areas,
+ report = TRUE,
+ g3l_distribution_surveyindices_log(alpha = ~g3_param("si_alpha", value = 0), beta = ~g3_param("si_beta", value = 0)))))
+ actions <- c(actions, list(
+ g3a_report_history(actions, '__cons$', out_prefix = "detail_") ))
+
+ # 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, params, compile_flags = c("-O0", "-g"))
+ } else {
+ writeLines("# skip: not compiling TMB model")
+ }
+
+ params <- attr(model_fn, 'parameter_template')
+ params$fleet_abc_a <- runif(10, min=0.1, max=0.9)
+ params$fleet_abc_b <- runif(10, min=0.1, max=0.9)
+ params$fleet_abc_c <- runif(10, min=0.1, max=0.9)
+ params$amount_ab <- 1000000
+ params$si_alpha <- 1.82
+ params$si_beta <- 3.74
+ params$cdist_sumofsquares_utcd_weight <- 1
+ params$cdist_sumofsquares_utcd_weight_weight <- 1
+ params$cdist_sumofsquares_utsd_weight <- 1
+ params$cdist_multinomial_multinom_weight <- 1
+ params$adist_surveyindices_log_surveyindices_weight <- 1
+
+ result <- model_fn(params)
+ r <- attributes(result)
+ # str(result)
+ # str(as.list(r), vec.len = 10000)
+
+ ok(ut_cmp_equal(
+ sort(as.vector(r$cdist_sumofsquares_utsd_obs__num)),
+ sort(sd_data$number)), "cdist_sumofsquares_utsd_obs__num: Imported from data.frame, order not necessarily the same")
+
+ ok(ut_cmp_equal(
+ r$adist_surveyindices_log_surveyindices_model__params,
+ c(params$si_alpha, params$si_beta)), "adist_surveyindices_log_surveyindices_model__params: Reported our hard-coded linear regression parameters")
+
+
+ ######## cdist_sumofsquares_utsd_model__num
+ ok(ut_cmp_equal(as.vector(r$step1_cdist_sumofsquares_utsd_model__num[,'prey_a', 1]), c(
+ sum(
+ (r$detail_prey_a_fleet_abc__cons[1:5,,1] / r$prey_a__wgt[1:5,]),
+ (r$detail_prey_a_fleet_abc__cons[1:5,,2] / r$prey_a__wgt[1:5,])),
+ sum(
+ sum(r$detail_prey_a_fleet_abc__cons[6:10,,1] / r$prey_a__wgt[6:10,]),
+ sum(r$detail_prey_a_fleet_abc__cons[6:10,,2] / r$prey_a__wgt[6:10,])),
+ NULL)), "step1_cdist_sumofsquares_utsd_model__num[,'prey_a',1]: prey_a summed over year")
+ ok(ut_cmp_equal(as.vector(r$step1_cdist_sumofsquares_utsd_model__num[,'prey_c', 1]), c(
+ 0,
+ 0,
+ NULL)), "step1_cdist_sumofsquares_utsd_model__num[,'prey_c',1]: No prey_c in first area")
+ ok(ut_cmp_equal(as.vector(r$step1_cdist_sumofsquares_utsd_model__num[,'prey_a', 2]), c(
+ 0,
+ 0,
+ NULL)), "step1_cdist_sumofsquares_utsd_model__num[,'prey_a',2]: No prey_a in second area")
+ ok(ut_cmp_equal(as.vector(r$step1_cdist_sumofsquares_utsd_model__num[,'prey_c', 2]), c(
+ sum(
+ (r$detail_prey_c_fleet_abc__cons[1:5,,1] / r$prey_c__wgt[1:5,]),
+ (r$detail_prey_c_fleet_abc__cons[1:5,,2] / r$prey_c__wgt[1:5,])),
+ sum(
+ sum(r$detail_prey_c_fleet_abc__cons[6:10,,1] / r$prey_c__wgt[6:10,]),
+ sum(r$detail_prey_c_fleet_abc__cons[6:10,,2] / r$prey_c__wgt[6:10,])),
+ NULL)), "step1_cdist_sumofsquares_utsd_model__num[,'prey_c',2]: prey_c summed over year")
+ ########
+
+ ######## cdist_sumofsquares_utcd_model__num
+ ok(ut_cmp_equal(as.vector(r$step1_cdist_sumofsquares_utcd_model__num[,1]), c(
+ sum(
+ sum(r$detail_prey_a_fleet_abc__cons[1:5,1,1] / r$prey_a__wgt[1:5,1]),
+ sum(r$detail_prey_b_fleet_abc__cons[1:5,1,1] / r$prey_b__wgt[1:5,1]),
+ sum(r$detail_prey_a_fleet_abc__cons[1:5,1,2] / r$prey_a__wgt[1:5,1]),
+ sum(r$detail_prey_b_fleet_abc__cons[1:5,1,2] / r$prey_b__wgt[1:5,1])),
+ sum(
+ sum(r$detail_prey_a_fleet_abc__cons[6:10,1,1] / r$prey_a__wgt[6:10,1]),
+ sum(r$detail_prey_b_fleet_abc__cons[6:10,1,1] / r$prey_b__wgt[6:10,1]),
+ sum(r$detail_prey_a_fleet_abc__cons[6:10,1,2] / r$prey_a__wgt[6:10,1]),
+ sum(r$detail_prey_b_fleet_abc__cons[6:10,1,2] / r$prey_b__wgt[6:10,1])),
+ NULL)), "step1_cdist_sumofsquares_utcd_model__num[,1]: all prey from a/b and steps 0 & 1")
+
+ ok(ut_cmp_equal(as.vector(r$step1_cdist_sumofsquares_utcd_model__num[,2]), c(
+ sum(
+ sum(r$detail_prey_c_fleet_abc__cons[1:5,1,1] / r$prey_c__wgt[1:5,1]),
+ sum(r$detail_prey_c_fleet_abc__cons[1:5,1,2] / r$prey_c__wgt[1:5,1])),
+ sum(
+ sum(r$detail_prey_c_fleet_abc__cons[6:10,1,1] / r$prey_c__wgt[6:10,1]),
+ sum(r$detail_prey_c_fleet_abc__cons[6:10,1,2] / r$prey_c__wgt[6:10,1])),
+ NULL)), "step1_cdist_sumofsquares_utcd_model__num[,2]: all prey from c and steps 0/1")
+ ########
+
+ ######## cdist_sumofsquares_utcd_weight_model__wgt
+ ok(ut_cmp_equal(as.vector(r$step1_cdist_sumofsquares_utcd_weight_model__wgt[,1]), c(
+ sum(
+ sum(r$detail_prey_a_fleet_abc__cons[1:5,1,1]),
+ sum(r$detail_prey_b_fleet_abc__cons[1:5,1,1]),
+ sum(r$detail_prey_a_fleet_abc__cons[1:5,1,2]),
+ sum(r$detail_prey_b_fleet_abc__cons[1:5,1,2])),
+ sum(
+ sum(r$detail_prey_a_fleet_abc__cons[6:10,1,1]),
+ sum(r$detail_prey_b_fleet_abc__cons[6:10,1,1]),
+ sum(r$detail_prey_a_fleet_abc__cons[6:10,1,2]),
+ sum(r$detail_prey_b_fleet_abc__cons[6:10,1,2])),
+ NULL)), "step1_cdist_sumofsquares_utcd_weight_model__wgt[,1]: total biomass of prey from a/b and steps 0 & 1")
+
+ ok(ut_cmp_equal(as.vector(r$step1_cdist_sumofsquares_utcd_weight_model__wgt[,2]), c(
+ sum(
+ sum(r$detail_prey_c_fleet_abc__cons[1:5,1,1]),
+ sum(r$detail_prey_c_fleet_abc__cons[1:5,1,2])),
+ sum(
+ sum(r$detail_prey_c_fleet_abc__cons[6:10,1,1]),
+ sum(r$detail_prey_c_fleet_abc__cons[6:10,1,2])),
+ NULL)), "step1_cdist_sumofsquares_utcd_weight_model__wgt[,2]: total biomass of prey from c and steps 0/1")
+ ########
+
+ ######## adist_surveyindices_surveyindices_model__num
+ ok(ut_cmp_equal(
+ as.vector(r$adist_surveyindices_log_surveyindices_model__num),
+ c(
+ sum(r$step0_prey_b__num) + sum(r$step1_prey_b__num),
+ sum(r$step2_prey_b__num) + sum(r$step3_prey_b__num),
+ NULL)), "adist_surveyindices_log_surveyindices_model__num: Built-in reporting gave us step 0..3 abundance")
+ ########
+
+ ok(ut_cmp_equal(r$step0_nll, 0), "step0_nll: nll not calculated yet")
+
+ ok(ut_cmp_equal(r$step1_nll, sum(
+ # utsd: stock 1 / area 1
+ (r$step1_cdist_sumofsquares_utsd_model__num[,1,1] / sum(r$step1_cdist_sumofsquares_utsd_model__num[,,1]) -
+ # NB: Still using first time data, unlike per-step example
+ r$cdist_sumofsquares_utsd_obs__num[,1,1,1] / sum(r$cdist_sumofsquares_utsd_obs__num[,,1,1])) ** 2,
+ # utsd: stock 2 / area 1
+ (r$step1_cdist_sumofsquares_utsd_model__num[,2,1] / sum(r$step1_cdist_sumofsquares_utsd_model__num[,,1]) -
+ r$cdist_sumofsquares_utsd_obs__num[,2,1,1] / sum(r$cdist_sumofsquares_utsd_obs__num[,,1,1])) ** 2,
+ # utsd: stock 1 / area 2
+ (r$step1_cdist_sumofsquares_utsd_model__num[,1,2] / sum(r$step1_cdist_sumofsquares_utsd_model__num[,,2]) -
+ r$cdist_sumofsquares_utsd_obs__num[,1,1,2] / sum(r$cdist_sumofsquares_utsd_obs__num[,,1,2])) ** 2,
+ # utsd: stock 2 / area 2
+ (r$step1_cdist_sumofsquares_utsd_model__num[,2,2] / sum(r$step1_cdist_sumofsquares_utsd_model__num[,,2]) -
+ r$cdist_sumofsquares_utsd_obs__num[,2,1,2] / sum(r$cdist_sumofsquares_utsd_obs__num[,,1,2])) ** 2,
+ # utcd: area 1
+ (r$step1_cdist_sumofsquares_utcd_model__num[,1] / sum(r$step1_cdist_sumofsquares_utcd_model__num[,1]) -
+ r$cdist_sumofsquares_utcd_obs__num[,1,1] / sum(r$cdist_sumofsquares_utcd_obs__num[,1,1])) ** 2,
+ # utcd: area 2
+ (r$step1_cdist_sumofsquares_utcd_model__num[,2] / sum(r$step1_cdist_sumofsquares_utcd_model__num[,2]) -
+ r$cdist_sumofsquares_utcd_obs__num[,1,2] / sum(r$cdist_sumofsquares_utcd_obs__num[,1,2])) ** 2,
+ # utcd_weight: area 1
+ (r$step1_cdist_sumofsquares_utcd_weight_model__wgt[,1] / sum(r$step1_cdist_sumofsquares_utcd_weight_model__wgt[,1]) -
+ r$cdist_sumofsquares_utcd_weight_obs__wgt[,1,1] / sum(r$cdist_sumofsquares_utcd_weight_obs__wgt[,1,1])) ** 2,
+ # utcd_weight: area 2
+ (r$step1_cdist_sumofsquares_utcd_weight_model__wgt[,2] / sum(r$step1_cdist_sumofsquares_utcd_weight_model__wgt[,2]) -
+ r$cdist_sumofsquares_utcd_weight_obs__wgt[,1,2] / sum(r$cdist_sumofsquares_utcd_weight_obs__wgt[,1,2])) ** 2,
+ # multinom:
+ (2 * (-sum(r$cdist_multinomial_multinom_obs__num[,1] * log(g3_logspace_add(r$step1_cdist_multinomial_multinom_model__num/g3_avoid_zero(sum(r$step1_cdist_multinomial_multinom_model__num)) *
+ 10000, (1/(length(r$cdist_multinomial_multinom_obs__num[,1]) * 10)) * 10000)/10000)) +
+ (sum(g3_lgamma_vec(1 + r$cdist_multinomial_multinom_obs__num[,1])) - lgamma(1 +
+ sum(r$cdist_multinomial_multinom_obs__num[,1]))))),
+ # surveyindices:
+ 0, # NB: We don't calculate until end
+ 0)), "step1_nll: Sum of squares")
+
+ ok(ut_cmp_equal(r$step3_nll, sum(
+ # utsd: stock 1 / area 1
+ (r$step3_cdist_sumofsquares_utsd_model__num[,1,1] / g3_avoid_zero(sum(r$step3_cdist_sumofsquares_utsd_model__num[,,1])) -
+ # NB: Second time data, not 4th as in per-step example
+ r$cdist_sumofsquares_utsd_obs__num[,1,2,1] / g3_avoid_zero(sum(r$cdist_sumofsquares_utsd_obs__num[,,2,1]))) ** 2,
+ # utsd: stock 2 / area 1
+ (r$step3_cdist_sumofsquares_utsd_model__num[,2,1] / g3_avoid_zero(sum(r$step3_cdist_sumofsquares_utsd_model__num[,,1])) -
+ r$cdist_sumofsquares_utsd_obs__num[,2,2,1] / g3_avoid_zero(sum(r$cdist_sumofsquares_utsd_obs__num[,,2,1]))) ** 2,
+ # utsd: stock 1 / area 2
+ (r$step3_cdist_sumofsquares_utsd_model__num[,1,2] / g3_avoid_zero(sum(r$step3_cdist_sumofsquares_utsd_model__num[,,2])) -
+ r$cdist_sumofsquares_utsd_obs__num[,1,2,2] / g3_avoid_zero(sum(r$cdist_sumofsquares_utsd_obs__num[,,2,2]))) ** 2,
+ # utsd: stock 2 / area 2
+ (r$step3_cdist_sumofsquares_utsd_model__num[,2,2] / g3_avoid_zero(sum(r$step3_cdist_sumofsquares_utsd_model__num[,,2])) -
+ r$cdist_sumofsquares_utsd_obs__num[,2,2,2] / g3_avoid_zero(sum(r$cdist_sumofsquares_utsd_obs__num[,,2,2]))) ** 2,
+ # utcd: area 1
+ (r$step3_cdist_sumofsquares_utcd_model__num[,1] / g3_avoid_zero(sum(r$step3_cdist_sumofsquares_utcd_model__num[,1])) -
+ r$cdist_sumofsquares_utcd_obs__num[,2,1] / g3_avoid_zero(sum(r$cdist_sumofsquares_utcd_obs__num[,2,1]))) ** 2,
+ # utcd: area 2
+ (r$step3_cdist_sumofsquares_utcd_model__num[,2] / g3_avoid_zero(sum(r$step3_cdist_sumofsquares_utcd_model__num[,2])) -
+ r$cdist_sumofsquares_utcd_obs__num[,2,2] / g3_avoid_zero(sum(r$cdist_sumofsquares_utcd_obs__num[,2,2]))) ** 2,
+ # utcd_weight: area 1
+ (r$step3_cdist_sumofsquares_utcd_weight_model__wgt[,1] / g3_avoid_zero(sum(r$step3_cdist_sumofsquares_utcd_weight_model__wgt[,1])) -
+ r$cdist_sumofsquares_utcd_weight_obs__wgt[,2,1] / g3_avoid_zero(sum(r$cdist_sumofsquares_utcd_weight_obs__wgt[,2,1]))) ** 2,
+ # utcd_weight: area 2
+ (r$step3_cdist_sumofsquares_utcd_weight_model__wgt[,2] / g3_avoid_zero(sum(r$step3_cdist_sumofsquares_utcd_weight_model__wgt[,2])) -
+ r$cdist_sumofsquares_utcd_weight_obs__wgt[,2,2] / g3_avoid_zero(sum(r$cdist_sumofsquares_utcd_weight_obs__wgt[,2,2]))) ** 2,
+ # multinom:
+ (2 * (-sum(r$cdist_multinomial_multinom_obs__num[,2] * log(g3_logspace_add(r$step3_cdist_multinomial_multinom_model__num/g3_avoid_zero(sum(r$step3_cdist_multinomial_multinom_model__num)) *
+ 10000, (1/(length(r$cdist_multinomial_multinom_obs__num[,2]) * 10)) * 10000)/10000)) +
+ (sum(g3_lgamma_vec(1 + r$cdist_multinomial_multinom_obs__num[,2])) - lgamma(1 +
+ sum(r$cdist_multinomial_multinom_obs__num[,2]))))),
+ # surveyindices:
+ sum((params$si_alpha +
+ params$si_beta * log(g3_avoid_zero(r$adist_surveyindices_log_surveyindices_model__num[,])) -
+ log(g3_avoid_zero(r$adist_surveyindices_log_surveyindices_obs__num[,])))**2),
+ r$step1_nll)), "step3_nll: Sum of squares, including step1_nll")
+
+ 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)
+ }
+ })
# Likelihood per year
# skip: not compiling TMB model
ok - cdist_sumofsquares_utsd_obs__num: Imported from data.frame, order not necessarily the same
ok - adist_surveyindices_log_surveyindices_model__params: Reported our hard-coded linear regression parameters
ok - step1_cdist_sumofsquares_utsd_model__num[,'prey_a',1]: prey_a summed over year
ok - step1_cdist_sumofsquares_utsd_model__num[,'prey_c',1]: No prey_c in first area
ok - step1_cdist_sumofsquares_utsd_model__num[,'prey_a',2]: No prey_a in second area
ok - step1_cdist_sumofsquares_utsd_model__num[,'prey_c',2]: prey_c summed over year
ok - step1_cdist_sumofsquares_utcd_model__num[,1]: all prey from a/b and steps 0 & 1
ok - step1_cdist_sumofsquares_utcd_model__num[,2]: all prey from c and steps 0/1
ok - step1_cdist_sumofsquares_utcd_weight_model__wgt[,1]: total biomass of prey from a/b and steps 0 & 1
ok - step1_cdist_sumofsquares_utcd_weight_model__wgt[,2]: total biomass of prey from c and steps 0/1
ok - adist_surveyindices_log_surveyindices_model__num: Built-in reporting gave us step 0..3 abundance
ok - step0_nll: nll not calculated yet
ok - step1_nll: Sum of squares
ok - step3_nll: Sum of squares, including step1_nll
>
> proc.time()
user system elapsed
12.20 0.29 12.75
1..71
# Looks like you failed 1 of 71 tests.
# 9: g3l_distribution_sumofsquares statified over length
Flavor: r-devel-windows-x86_64
Current CRAN status: NOTE: 4, OK: 9
Version: 7.3-1
Check: for unstated dependencies in ‘demo’
Result: NOTE
'library' or 'require' call not declared from: ‘tibble’
Flavors: r-devel-linux-x86_64-debian-clang, r-devel-linux-x86_64-debian-gcc, r-patched-linux-x86_64, r-release-linux-x86_64
Current CRAN status: OK: 13
These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.