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)
> 
> # 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_logspace_add_vec <- g3_env$logspace_add_vec
> 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")
+ 
+ })
# g3_distribution_preview
ok - Returned number array
ok - Returned weight array
ok - Number array wins if both present
> 
> 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()
ok - g3l_distribution_sumofsquares statified over length
> 
> 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("wtperstock",
+             obsdata,
+             function_f = g3l_distribution_sumofsquares(),
+             stocks = list(prey_a),
+             transform_fs = list(age = g3_formula( stock_switch(stock,
+                 prey_a = g3_param_array('reader1matrix', value = diag(prey_a__maxage - prey_a__minage + 1)),
+                 x = 0)[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)
+ 
+     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 - 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 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 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"),
+             prey_b = ~g3_param_vector("fleet_abc_b"),
+             prey_c = ~g3_param_vector("fleet_abc_c")),
+         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__predby_fleet_abc'), report_step('prey_a__predby_fleet_abc', 4, g3_stock_instance(prey_a)),
+         gadget3:::step_id(990, 'prey_b__predby_fleet_abc'), report_step('prey_b__predby_fleet_abc', 4, g3_stock_instance(prey_b)),
+         gadget3:::step_id(990, 'prey_c__predby_fleet_abc'), report_step('prey_c__predby_fleet_abc', 4, g3_stock_instance(prey_c)),
+         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))
> 
> # 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$step0_prey_a__predby_fleet_abc[1:5,] / r$prey_a__wgt[1:5,]),
+         sum(r$step0_prey_a__predby_fleet_abc[6:10,] / 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$step0_prey_c__predby_fleet_abc[1:5,] / r$prey_c__wgt[1:5,]),
+         sum(r$step0_prey_c__predby_fleet_abc[6:10,] / 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$step1_prey_a__predby_fleet_abc[1:5,] / r$prey_a__wgt[1:5,]),
+         sum(r$step1_prey_a__predby_fleet_abc[6:10,] / 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$step1_prey_c__predby_fleet_abc[1:5,] / r$prey_c__wgt[1:5,]),
+         sum(r$step1_prey_c__predby_fleet_abc[6:10,] / 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$step2_prey_a__predby_fleet_abc[1:5,] / r$prey_a__wgt[1:5,]),
+         sum(r$step2_prey_a__predby_fleet_abc[6:10,] / 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$step2_prey_c__predby_fleet_abc[1:5,] / r$prey_c__wgt[1:5,]),
+         sum(r$step2_prey_c__predby_fleet_abc[6:10,] / 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$step3_prey_a__predby_fleet_abc[1:5,] / r$prey_a__wgt[1:5,]),
+         sum(r$step3_prey_a__predby_fleet_abc[6:10,] / 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$step3_prey_c__predby_fleet_abc[1:5,] / r$prey_c__wgt[1:5,]),
+         sum(r$step3_prey_c__predby_fleet_abc[6:10,] / 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$step0_prey_a__predby_fleet_abc[1:5,1] / r$prey_a__wgt[1:5,1]),
+              sum(r$step0_prey_b__predby_fleet_abc[1:5,1] / r$prey_b__wgt[1:5,1])),
+          sum(
+              sum(r$step0_prey_a__predby_fleet_abc[6:10,1] / r$prey_a__wgt[6:10,1]),
+              sum(r$step0_prey_b__predby_fleet_abc[6:10,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$step0_prey_c__predby_fleet_abc[1:5,1] / r$prey_c__wgt[1:5,1])),
+          sum(
+              sum(r$step0_prey_c__predby_fleet_abc[6:10,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$step1_prey_a__predby_fleet_abc[1:5,1] / r$prey_a__wgt[1:5,1]),
+              sum(r$step1_prey_b__predby_fleet_abc[1:5,1] / r$prey_b__wgt[1:5,1])),
+          sum(
+              sum(r$step1_prey_a__predby_fleet_abc[6:10,1] / r$prey_a__wgt[6:10,1]),
+              sum(r$step1_prey_b__predby_fleet_abc[6:10,1] / 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$step1_prey_c__predby_fleet_abc[1:5,1] / r$prey_c__wgt[1:5,1])),
+          sum(
+              sum(r$step1_prey_c__predby_fleet_abc[6:10,1] / 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$step2_prey_a__predby_fleet_abc[1:5,1] / r$prey_a__wgt[1:5,1]),
+              sum(r$step2_prey_b__predby_fleet_abc[1:5,1] / r$prey_b__wgt[1:5,1])),
+          sum(
+              sum(r$step2_prey_a__predby_fleet_abc[6:10,1] / r$prey_a__wgt[6:10,1]),
+              sum(r$step2_prey_b__predby_fleet_abc[6:10,1] / 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$step2_prey_c__predby_fleet_abc[1:5,1] / r$prey_c__wgt[1:5,1])),
+          sum(
+              sum(r$step2_prey_c__predby_fleet_abc[6:10,1] / 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$step3_prey_a__predby_fleet_abc[1:5,1] / r$prey_a__wgt[1:5,1]),
+              sum(r$step3_prey_b__predby_fleet_abc[1:5,1] / r$prey_b__wgt[1:5,1])),
+          sum(
+              sum(r$step3_prey_a__predby_fleet_abc[6:10,1] / r$prey_a__wgt[6:10,1]),
+              sum(r$step3_prey_b__predby_fleet_abc[6:10,1] / 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$step3_prey_c__predby_fleet_abc[1:5,1] / r$prey_c__wgt[1:5,1])),
+          sum(
+              sum(r$step3_prey_c__predby_fleet_abc[6:10,1] / 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_vec(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_vec(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_vec(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_vec(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)))))
+ 
+     # 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$step0_prey_a__predby_fleet_abc[1:5,] / r$prey_a__wgt[1:5,]),
+             (r$step1_prey_a__predby_fleet_abc[1:5,] / r$prey_a__wgt[1:5,])),
+         sum(
+             sum(r$step0_prey_a__predby_fleet_abc[6:10,] / r$prey_a__wgt[6:10,]),
+             sum(r$step1_prey_a__predby_fleet_abc[6:10,] / 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$step0_prey_c__predby_fleet_abc[1:5,] / r$prey_c__wgt[1:5,]),
+             (r$step1_prey_c__predby_fleet_abc[1:5,] / r$prey_c__wgt[1:5,])),
+         sum(
+             sum(r$step0_prey_c__predby_fleet_abc[6:10,] / r$prey_c__wgt[6:10,]),
+             sum(r$step1_prey_c__predby_fleet_abc[6:10,] / 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$step0_prey_a__predby_fleet_abc[1:5,1] / r$prey_a__wgt[1:5,1]),
+              sum(r$step0_prey_b__predby_fleet_abc[1:5,1] / r$prey_b__wgt[1:5,1]),
+              sum(r$step1_prey_a__predby_fleet_abc[1:5,1] / r$prey_a__wgt[1:5,1]),
+              sum(r$step1_prey_b__predby_fleet_abc[1:5,1] / r$prey_b__wgt[1:5,1])),
+          sum(
+              sum(r$step0_prey_a__predby_fleet_abc[6:10,1] / r$prey_a__wgt[6:10,1]),
+              sum(r$step0_prey_b__predby_fleet_abc[6:10,1] / r$prey_b__wgt[6:10,1]),
+              sum(r$step1_prey_a__predby_fleet_abc[6:10,1] / r$prey_a__wgt[6:10,1]),
+              sum(r$step1_prey_b__predby_fleet_abc[6:10,1] / 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$step0_prey_c__predby_fleet_abc[1:5,1] / r$prey_c__wgt[1:5,1]),
+              sum(r$step1_prey_c__predby_fleet_abc[1:5,1] / r$prey_c__wgt[1:5,1])),
+          sum(
+              sum(r$step0_prey_c__predby_fleet_abc[6:10,1] / r$prey_c__wgt[6:10,1]),
+              sum(r$step1_prey_c__predby_fleet_abc[6:10,1] / 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$step0_prey_a__predby_fleet_abc[1:5,1]),
+              sum(r$step0_prey_b__predby_fleet_abc[1:5,1]),
+              sum(r$step1_prey_a__predby_fleet_abc[1:5,1]),
+              sum(r$step1_prey_b__predby_fleet_abc[1:5,1])),
+          sum(
+              sum(r$step0_prey_a__predby_fleet_abc[6:10,1]),
+              sum(r$step0_prey_b__predby_fleet_abc[6:10,1]),
+              sum(r$step1_prey_a__predby_fleet_abc[6:10,1]),
+              sum(r$step1_prey_b__predby_fleet_abc[6:10,1])),
+          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$step0_prey_c__predby_fleet_abc[1:5,1]),
+              sum(r$step1_prey_c__predby_fleet_abc[1:5,1])),
+          sum(
+              sum(r$step0_prey_c__predby_fleet_abc[6:10,1]),
+              sum(r$step1_prey_c__predby_fleet_abc[6:10,1])),
+          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_vec(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_vec(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 
   6.25    0.14    6.39 
# Looks like you passed all 68 tests.