Running ‘test-aab_env.R’ [1s/1s]
  Running ‘test-action_age.R’ [2s/3s]
  Running ‘test-action_grow.R’ [1s/3s]
  Running ‘test-action_mature.R’ [1s/2s]
  Running ‘test-action_migrate.R’ [2s/3s]
  Running ‘test-action_naturalmortality.R’ [1s/2s]
  Running ‘test-action_predate-catchability.R’ [4s/7s]
  Running ‘test-action_predate-numberfleet.R’ [2s/3s]
  Running ‘test-action_predate.R’ [3s/6s]
  Running ‘test-action_renewal.R’ [1s/3s]
  Running ‘test-action_report.R’ [1s/2s]
  Running ‘test-action_spawn.R’ [3s/5s]
  Running ‘test-action_tagging.R’ [2s/3s]
  Running ‘test-action_time.R’ [1s/2s]
  Running ‘test-formula_utils.R’ [1s/1s]
  Running ‘test-likelihood_data.R’ [1s/2s]
  Running ‘test-likelihood_distribution.R’ [10s/16s]
  Running ‘test-likelihood_random.R’ [1s/2s]
  Running ‘test-likelihood_tagging_ckmr.R’ [3s/4s]
  Running ‘test-likelihood_understocking.R’ [0s/1s]
  Running ‘test-params.R’ [1s/1s]
  Running ‘test-run.R’ [1s/1s]
  Running ‘test-run_r.R’ [1s/1s]
  Running ‘test-run_tmb.R’ [2s/3s]
  Running ‘test-step.R’ [1s/3s]
  Running ‘test-stock.R’ [1s/1s]
  Running ‘test-stock_age.R’ [1s/2s]
  Running ‘test-stock_areas.R’ [1s/3s]
  Running ‘test-stock_tag.R’ [0s/1s]
  Running ‘test-stock_time.R’ [1s/1s]
  Running ‘test-suitability.R’ [2s/3s]
  Running ‘test-timedata.R’ [1s/2s]
  Running ‘test-timevariable.R’ [1s/2s]
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_logspace_add_vec <- g3_env$logspace_add_vec
  > g3_lgamma_vec <- lgamma
  > 
  > 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) ),
  +         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: 3.131472e-08
  # --- as.numeric(r)
  # +++ expected_nll
  #  num [-2.17527312-]{+2.175273188+}
  > 
  > 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))[modelstock__preage_idx, modelstock__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')
  +         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(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 - age1: Identity matrix
  ok - age2: Identity matrix
  ok - age3: Identity matrix
  ok - age4: Identity matrix
  ok - age5: Identity matrix
  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
  > 
  > # Can't make g3l_distribution_surveyindices without report
  > ok(ut_cmp_error(g3l_abundancedistribution(
  +     'surveyindices',
  +     surveyindices_data,
  +     fleets = list(),
  +     stocks = list(prey_b),
  +     area_group = areas,
  +     g3l_distribution_surveyindices_linear(alpha = ~g3_param("si_alpha", value = 0), beta = ~g3_param("si_beta", value = 0)),
  +     report = FALSE), "report"), "Can't have report = FALSE & g3l_distribution_surveyindices_linear")
  ok - Can't have report = FALSE & g3l_distribution_surveyindices_linear
  > ok(ut_cmp_error(g3l_abundancedistribution(
  +     'surveyindices',
  +     surveyindices_data,
  +     fleets = list(),
  +     stocks = list(prey_b),
  +     area_group = areas,
  +     g3l_distribution_surveyindices_log(alpha = ~g3_param("si_alpha", value = 0), beta = ~g3_param("si_beta", value = 0)),
  +     report = FALSE), "report"), "Can't have report = FALSE & g3l_distribution_surveyindices_log")
  ok - Can't have report = FALSE & g3l_distribution_surveyindices_log
  > 
  > # 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_substitute(~if (cur_time == i) {
  +             comment(report_comment)
  +             step_var[] <- var
  +             REPORT(step_var)
  +         } else rest, list(
  +             report_comment = paste0("Reporting on ", var_name, " at step ", i),
  +             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 
    9.702   0.160  15.699 
  # Looks like you failed 1 of 65 tests.