diff --git a/R/convert_output.R b/R/convert_output.R index 62fed9e0..f3a509c2 100644 --- a/R/convert_output.R +++ b/R/convert_output.R @@ -1981,7 +1981,10 @@ convert_output <- function( # can late take approach like BAM? # TODO: Do we want users to input the saved file or already loaded into the R environment? if (is.character(file)) { - if (!grepl(".rds", file)) { + # Anchored and case-insensitive: saveRDS() output is conventionally named + # either .rds or .RDS, and the real assessment files this targets use the + # latter. The unanchored, case-sensitive form rejected them. + if (!grepl("\\.rds$", file, ignore.case = TRUE)) { cli::cli_abort("File must be in an .rds format otherwise load data into your environment and run `convert_output()` using the loaded data.") } dat <- readRDS(file) @@ -1989,9 +1992,38 @@ convert_output <- function( dat <- file } + # The standardized output describes one stock and has no species dimension, + # but CEATTLE is a multispecies model. A multispecies fit's `sdrep` holds + # nspp * nyrs values per label, which fails the `n == length(styr:projyr)` + # test used to attach years below, so every biomass / spawning biomass / + # recruitment value would silently come back with `year = NA`. Refuse + # instead: a converted file that quietly merges two stocks is worse than no + # converted file. + # `[[` rather than `$`: `$` partial-matches, so a renamed element would + # return something else and skip the refusal. + if (!"data_list" %in% names(dat)) { + cli::cli_abort(c( + "Rceattle fit has no {.field data_list} element.", + i = "The catch, index and composition data are read from it." + )) + } + nspp <- dat[["data_list"]][["nspp"]] + if (is.null(nspp)) { + cli::cli_abort(c( + "Rceattle fit does not report {.field nspp}, so its species count is unknown.", + i = "The standardized output has no species dimension and must not merge stocks." + )) + } + if (nspp > 1) { + cli::cli_abort(c( + "Rceattle model has {nspp} species, and the standardized output has no species dimension.", + i = "Subset the fit to a single species before converting." + )) + } + # Extract or use fleet names if (is.null(fleet_names)) { - fleet_names <- names(dat$estimated_params$index_ln_q) + fleet_names <- names(dat$estimated_params$index_log_q) } # Output fleet names in console @@ -2005,9 +2037,25 @@ convert_output <- function( # units <- c("mt", "lbs", "eggs") ##### Loop #### - for (p in (2:length(dat))[-c(6, 9, 10)]) { - extract <- dat[p] - module_name <- names(extract) + # Select the modules by name rather than by position. An Rceattle fit's + # element list depends on how it was fitted -- `phase = TRUE` adds + # `.conv_phase` and `phase_params`, and `getsd = FALSE` drops `sdrep` -- so a + # fixed offset selects different modules for different fits of the same + # model. Under `phase = TRUE` the positional form dropped `data_list` and + # passed `obj` to the generic vector handler below, which aborts with + # "attempt to select less than one element in OneIndex". + # + # `sdrep` carries every derived time series, so without it the output holds + # only the input data. Say so rather than returning a table whose biomass + # and recruitment are simply absent. + if (!"sdrep" %in% names(dat)) { + cli::cli_warn(c( + "Rceattle fit has no {.field sdrep}, so the converted output has no biomass, spawning biomass or recruitment.", + i = "Refit with {.code fit_control(getsd = TRUE)} to standardize the derived time series." + )) + } + for (module_name in intersect(c("sdrep", "data_list", "quantities"), names(dat))) { + extract <- dat[module_name] cli::cli_alert_info("Processing {module_name}") if (module_name == "sdrep") { ##### sdrep #### @@ -2034,29 +2082,49 @@ convert_output <- function( ) # make year column year_col <- rep( - file[["data_list"]]$styr:file[["data_list"]]$projyr, + dat[["data_list"]]$styr:dat[["data_list"]]$projyr, length(unique( - dplyr::filter(values_count, n == length(file[["data_list"]]$styr:file[["data_list"]]$projyr)) |> + dplyr::filter(values_count, n == length(dat[["data_list"]]$styr:dat[["data_list"]]$projyr)) |> dplyr::pull(label) )) ) df2 <- values |> - dplyr::filter(n == length(file[["data_list"]]$styr:file[["data_list"]]$projyr)) |> + dplyr::filter(n == length(dat[["data_list"]]$styr:dat[["data_list"]]$projyr)) |> dplyr::mutate(year = year_col) df2 <- values |> dplyr::filter( - n != length(file[["data_list"]]$styr:file[["data_list"]]$projyr) + n != length(dat[["data_list"]]$styr:dat[["data_list"]]$projyr) ) |> dplyr::mutate(year = NA) |> rbind(df2) # Extract parameter values ts + # The standard errors of the fixed effects are the square root of the + # diagonal of cov.fixed, in the same order as par.fixed. + has_cov_fixed <- !is.null(extract[[1]]$cov.fixed) par_fixes <- data.frame( label = names(extract[[1]]$par.fixed), - estimate = extract[[1]]$par.fixed + estimate = extract[[1]]$par.fixed, + uncertainty = if (has_cov_fixed) sqrt(diag(extract[[1]]$cov.fixed)) else NA_real_, + uncertainty_label = if (has_cov_fixed) "sd" else NA_character_ ) + # Random effects, such as rec_dev under random_rec = TRUE, are in + # par.random rather than par.fixed. Their variances are in + # diag.cov.random. + if (length(extract[[1]]$par.random) > 0) { + has_cov_random <- !is.null(extract[[1]]$diag.cov.random) + par_fixes <- rbind( + par_fixes, + data.frame( + label = names(extract[[1]]$par.random), + estimate = extract[[1]]$par.random, + uncertainty = if (has_cov_random) sqrt(extract[[1]]$diag.cov.random) else NA_real_, + uncertainty_label = if (has_cov_random) "sd" else NA_character_ + ) + ) + } par_fixes_count <- par_fixes |> dplyr::group_by(label) |> dplyr::count() @@ -2067,26 +2135,22 @@ convert_output <- function( ) year_col_par_fix <- rep( - file[["data_list"]]$styr:file[["data_list"]]$endyr, + dat[["data_list"]]$styr:dat[["data_list"]]$endyr, length(unique( - dplyr::filter(par_fixes_count, n == length(file[["data_list"]]$styr:file[["data_list"]]$endyr)) |> + dplyr::filter(par_fixes_count, n == length(dat[["data_list"]]$styr:dat[["data_list"]]$endyr)) |> dplyr::pull(label) )) ) df3 <- par_fixes |> - dplyr::filter(n == length(file[["data_list"]]$styr:file[["data_list"]]$endyr)) |> + dplyr::filter(n == length(dat[["data_list"]]$styr:dat[["data_list"]]$endyr)) |> dplyr::mutate(year = year_col_par_fix) df3 <- par_fixes |> dplyr::filter( - n != length(file[["data_list"]]$styr:file[["data_list"]]$endyr) + n != length(dat[["data_list"]]$styr:dat[["data_list"]]$endyr) ) |> dplyr::mutate(year = NA) |> - rbind(df3) |> - dplyr::mutate( - uncertainty = NA, - uncertainty_label = NA - ) + rbind(df3) # not sure how pop_scalar is indexed # not sure how log_index_hat is indexes # Did not use r_sd for the error in rec bc used it from the other element in the list @@ -2121,7 +2185,10 @@ convert_output <- function( uncertainty_label = "log_sd", indices_predicted = dat$quantities$index_hat ) |> - dplyr::select(-c(fleet_code, q_block)) + # any_of(): `Q_block` is deprecated in Rceattle and no longer written + # to `index_data`, so naming it directly fails on a current fit while + # still needing to be dropped from an older one. + dplyr::select(-dplyr::any_of(c("fleet_code", "q_block"))) cli::cli_alert_info("'Selectivity_block' values located in 'block' columns.") # check if species > 1 @@ -2159,12 +2226,16 @@ convert_output <- function( # dplyr::filter(!is.na(Catch)) |> dplyr::rename_with(tolower) |> dplyr::mutate( + module_name = names(extract), + # Set the era here from the year, as the index and composition + # frames do, so the frame is complete where it is built. The mutate + # at the end of the branch recomputes it the same way. era = dplyr::if_else( year > dat$data_list$endyr, "fore", - NA_character_ + "time" ), - catch_h = dat$quantities$catch_h, + catch_h = dat$quantities$catch_hat, # TODO: follow up on this quantity # log_index_hat = dat$quantities$log_index_hat uncertainty_label = "log_sd" @@ -2203,8 +2274,15 @@ convert_output <- function( indexing_vars_cols <- colnames(df_comp_obs)[!grepl("comp", colnames(df_comp_obs))] + # `age_hat` holds one row per row of `comp_data`, in the same order, so + # it binds column-wise. cross_join() paired every observation with the + # whole prediction matrix instead, which is a Cartesian product: on the + # GOAatf fixture it returned the 2,478 correct predictions 59 times + # over, once per row of comp_data, for 146,202 in all. + age_hat <- as.data.frame(dat$quantities$age_hat) + names(age_hat) <- paste0("agehat_", seq_len(ncol(age_hat))) df_comp_pred <- dplyr::select(df_comp_obs, dplyr::all_of(indexing_vars_cols)) |> - dplyr::cross_join(as.data.frame(dat$quantities$age_hat)) |> + cbind(age_hat) |> tidyr::pivot_longer( cols = -dplyr::all_of(indexing_vars_cols), names_to = "age", @@ -2215,7 +2293,7 @@ convert_output <- function( dplyr::mutate( label = "composition_predicted", # NOTE: the below age mutate slows down code - age = stringr::str_replace(age, "V", "") + age = stringr::str_replace(age, "agehat_", "") ) # Finish adjusting comp_obs @@ -2234,15 +2312,77 @@ convert_output <- function( ) df_comp <- rbind(df_comp_obs, df_comp_pred) |> + # Match the index and catch frames: the standard names the fleet + # `fleet`, carries no `fleet_name` or `species`, and every row states + # which module it came from. + dplyr::rename(fleet = fleet_name) |> + dplyr::mutate(module_name = names(extract)) |> + dplyr::select(-dplyr::any_of("species")) |> # TODO: do I need to filter sample size by <3 for confidentiality or does this not apply? dplyr::select(-c(age0_length1, fleet_code, sample_size)) - df_catch[setdiff(tolower(names(out_new)), tolower(names(df_catch)))] <- NA - data_list_list[["comp_data"]] <- df_catch + df_comp[setdiff(tolower(names(out_new)), tolower(names(df_comp)))] <- NA + data_list_list[["comp_data"]] <- df_comp # final df for data_list module + # rbind() rather than bind_rows(): every frame has already been padded + # to the standard column set above, so rbind() succeeds and keeps the + # column contract -- it aborts on a frame carrying an extra column, + # where bind_rows() would silently widen the output. new_df <- Reduce(rbind, data_list_list) out_list[[names(extract)]] <- new_df + } else if (module_name == "quantities") { + ###### quantities #### + # Only selectivity and natural mortality at age for now. Both are + # [fleet or species, sex, age, year] arrays. The age dimension is the + # bin index, so the age itself is offset by minage. + ages <- dat[["data_list"]][["minage"]] - 1 + + seq_len(dat[["data_list"]][["nages"]]) + nsex <- dat[["data_list"]][["nsex"]] + # Rceattle stores sex 1 as females, or both sexes when nsex = 1. + sex_names <- if (nsex == 1) NA_character_ else c("female", "male") + quantities_list <- list() + + if ("sel_at_age" %in% names(extract[[1]])) { + sel <- extract[[1]][["sel_at_age"]][, seq_len(nsex), seq_along(ages), , drop = FALSE] + df_sel <- expand.grid( + fleet = dimnames(sel)[[1]], + sex = sex_names, + age = ages, + year = as.numeric(dimnames(sel)[[4]]), + stringsAsFactors = FALSE + ) |> + dplyr::mutate( + estimate = as.vector(sel), + label = "sel_at_age" + ) + df_sel[setdiff(tolower(names(out_new)), tolower(names(df_sel)))] <- NA + quantities_list[["sel_at_age"]] <- df_sel + } + + if ("M_at_age" %in% names(extract[[1]])) { + m_at_age <- extract[[1]][["M_at_age"]][1, seq_len(nsex), seq_along(ages), , drop = FALSE] + df_m <- expand.grid( + sex = sex_names, + age = ages, + year = as.numeric(dimnames(m_at_age)[[4]]), + stringsAsFactors = FALSE + ) |> + dplyr::mutate( + estimate = as.vector(m_at_age), + label = "m_at_age" + ) + df_m[setdiff(tolower(names(out_new)), tolower(names(df_m)))] <- NA + quantities_list[["m_at_age"]] <- df_m + } + + if (length(quantities_list) > 0) { + out_list[[names(extract)]] <- Reduce(rbind, quantities_list) |> + dplyr::mutate( + module_name = names(extract), + age = as.character(age) + ) + } } else if (is.list(extract[[1]])) { # indicates vector and list ##### remaining lmnts #### if (any(vapply(extract[[1]], is.matrix, FUN.VALUE = logical(1)))) { @@ -2301,6 +2441,8 @@ convert_output <- function( # } # close if statement } # close loop over objects listed in dat file # Finish out df + # rbind() as above, so a module that has drifted from the standard column + # set aborts here rather than widening the output. out_new <- Reduce(rbind, out_list) |> # Add era as factor into BAM conout dplyr::mutate( diff --git a/inst/resources/rceattle_var_names.csv b/inst/resources/rceattle_var_names.csv index bc31703a..a0467530 100644 --- a/inst/resources/rceattle_var_names.csv +++ b/inst/resources/rceattle_var_names.csv @@ -31,8 +31,8 @@ map,alpha,, map,beta,, map,beta_rec_pars,, map,catch_ln_sd,ln_catch_sd, -map,comp_weights,composition_weights,kg -map,diet_comp_weights,diet_composition_weights,kg +map,comp_weights,composition_weights, +map,diet_comp_weights,diet_composition_weights, map,dummy,, map,index_ln_q,, map,index_ln_sd,, @@ -72,8 +72,8 @@ phase_params,alpha,, phase_params,beta,, phase_params,beta_rec_pars,, phase_params,catch_ln_sd,ln_catch_sd, -phase_params,comp_weights,composition_weights,kg -phase_params,diet_comp_weights,diet_composition_weights,kg +phase_params,comp_weights,composition_weights, +phase_params,diet_comp_weights,diet_composition_weights, phase_params,dummy,, phase_params,index_ln_q,, phase_params,index_ln_sd,, @@ -109,36 +109,36 @@ phase_params,sel_curve_pen,, phase_params,sel_dev_ln_sd,, phase_params,sel_inf,, phase_params,sel_inf_dev,, -quantities,age_hat,catch,kg -quantities,age_obs_hat,catch_adjusted,kg -quantities,avail_food,,kg +quantities,age_hat,catch,thousands of fish +quantities,age_obs_hat,catch_adjusted,thousands of fish +quantities,avail_food,,mt quantities,avg_r,recruitment_mean, -quantities,avgn_at_age,abundance_average,numbers -quantities,b0,,kg -quantities,b_eaten,,kg -quantities,b_eaten_as_prey,,kg -quantities,biomass,,kg -quantities,biomass_at_age,biomass,kg -quantities,biomass_depletion,biomass_biomass_unfished,kg -quantities,catch_hat,landings,kg +quantities,avgn_at_age,abundance_average,thousands of fish +quantities,b0,,mt +quantities,b_eaten,,mt +quantities,b_eaten_as_prey,,mt +quantities,biomass,,mt +quantities,biomass_at_age,biomass,mt +quantities,biomass_depletion,biomass_biomass_unfished,proportion +quantities,catch_hat,landings, quantities,comp_hat,, quantities,comp_n,, -quantities,comp_obs,,kg -quantities,consumption_at_age,,kg +quantities,comp_obs,,proportion +quantities,consumption_at_age,,kg per fish per yr quantities,diet_hat,, -quantities,diet_obs,,kg -quantities,dynamicb0,,kg -quantities,dynamicsb0,,kg -quantities,dynamicsbf,,kg -quantities,exploitable_biomass,,kg +quantities,diet_obs,,proportion +quantities,dynamicb0,,mt +quantities,dynamicsb0,,mt +quantities,dynamicsbf,,mt +quantities,exploitable_biomass,,mt quantities,f_at_age,fishing_mortality, quantities,f_flt,fishing_mortality, quantities,f_flt_age,fishing_mortality, quantities,f_spp,fishing_mortality, -quantities,flimit,fishing_mortality_target, +quantities,flimit,fishing_mortality_limit, quantities,flt_sel_maxage,, quantities,ft,, -quantities,ftarget,fishing_mortality_limit, +quantities,ftarget,fishing_mortality_target, quantities,gam_a,, quantities,gam_b,, quantities,index_hat,, @@ -148,16 +148,16 @@ quantities,jnll_comp,, quantities,ln_catch_sd,, quantities,ln_index_sd,, quantities,log_index_hat,, -quantities,log_ssb,log_spawning_biomass,kg +quantities,log_ssb,log_spawning_biomass,log mt quantities,m1_at_age,m1, quantities,m2_at_age,m2, quantities,m2_prop,proportion_m2, quantities,m_at_age,natural_mortality, -quantities,max_catch_hat,,kg -quantities,n_at_age,abundance,numbers -quantities,nbyage0,abundance_rec_hat_f0,numbers +quantities,max_catch_hat,, +quantities,n_at_age,abundance,thousands of fish +quantities,nbyage0,abundance_rec_hat_f0,thousands of fish quantities,nbyagef,abundance_rec_hat_f_target, -quantities,nbyagespr,abundance_spr,numbers +quantities,nbyagespr,abundance_spr,fish per recruit quantities,pop_scalar,, quantities,proj_f,fishing_mortality, quantities,r,recruitment, @@ -166,15 +166,16 @@ quantities,r_hat,, quantities,r_init,recruitment_initial, quantities,r_sd,recruitment_sd, quantities,ration,, -quantities,sb0,spawning_biomass_0,kg -quantities,sbf,spawning_biomass_target,kg +quantities,sb0,spawning_biomass_0,mt +quantities,sbf,spawning_biomass_target,mt quantities,sel,selectivity, +quantities,sel_at_age,selectivity_age, quantities,spr0,, quantities,sprfinit,spr_initial, quantities,sprlimit,spr_limit, quantities,sprtarget,spr_target, -quantities,ssb,spawning_biomass,kg -quantities,ssb_depletion,spawning_biomass_depletion,kg +quantities,ssb,spawning_biomass,mt +quantities,ssb_depletion,spawning_biomass_depletion,proportion quantities,steepness,, quantities,stom_div_bio,, quantities,suit_other,, @@ -184,19 +185,19 @@ quantities,unweighted_jnll_comp,, quantities,vulnerability,, quantities,vulnerability_other,, quantities,z_at_age,total_mortality, -sdrep,biomass,,kg +sdrep,biomass,,mt sdrep,index_ln_q,, sdrep,init_dev,, sdrep,ln_f,, sdrep,ln_sel_slp,, sdrep,ln_sel_slp_dev,, sdrep,log_index_hat,, -sdrep,log_ssb,log_spawning_biomass,kg -sdrep,pop_scalar,,kg +sdrep,log_ssb,log_spawning_biomass,log mt +sdrep,pop_scalar,,multiplier sdrep,r,recruitment, sdrep,r_sd,recruitment_sd, sdrep,rec_dev,recruitment_deviations, sdrep,rec_pars,recruitment_parameters, sdrep,sel_inf,, sdrep,sel_inf_dev,, -sdrep,ssb,spawning_biomass,kg +sdrep,ssb,spawning_biomass,mt diff --git a/tests/testthat/fixtures/prepare_rceattle_output_file.R b/tests/testthat/fixtures/prepare_rceattle_output_file.R new file mode 100644 index 00000000..85ce614c --- /dev/null +++ b/tests/testthat/fixtures/prepare_rceattle_output_file.R @@ -0,0 +1,69 @@ +# Build the Rceattle fixture used by test-convert_output.R. +# +# Source data is `GOAatf`, the Gulf of Alaska arrowtooth flounder dataset that +# ships with Rceattle, so the fixture contains nothing that is not already +# public. It is a single-species, two-sex assessment that converges with a +# positive-definite Hessian, which is what makes it exercise the `sdrep` +# branch of convert_output(). `GOAcod` is the one-sex alternative. +# +# Run from the package root, so test_path() resolves. The estimates depend on +# the Rceattle version, so a rebuild reproduces the fixture's structure and +# the test assertions, not its exact values. + +utils::data("GOAatf", package = "Rceattle", envir = environment()) + +# estimateMode = 1 fits the hindcast only; getsd = TRUE is required because the +# converter reads its time series out of sdrep$value. +full <- Rceattle::fit_mod( + data_list = GOAatf, + estimateMode = 1, + msmMode = 0, + fit_control = Rceattle::fit_control( + getsd = TRUE, + verbose = 0, + phase = FALSE + ) +) + +stopifnot( + !is.null(full$sdrep), + isTRUE(full$sdrep$pdHess), + full$data_list$nspp == 1 +) + +# A saved Rceattle fit is ~50 MB, nearly all of it the TMB objective-function +# object and the Hessian. Keep only the elements convert_output() reads, which +# brings the fixture to ~40 KB. +data_list_keep <- c( + "styr", "endyr", "projyr", "nspp", "nsex", "minage", "nages", "spnames", + "index_data", "catch_data", "comp_data", "fleet_control" +) +quantities_keep <- c( + "index_hat", "catch_hat", "age_hat", "log_index_hat", "sel_at_age", + "M_at_age" +) + +fixture <- structure( + list( + data_list = full$data_list[intersect(data_list_keep, names(full$data_list))], + sdrep = list( + value = full$sdrep$value, + sd = full$sdrep$sd, + par.fixed = full$sdrep$par.fixed, + # convert_output() reads only the variances. Zeroing the covariances + # keeps the fixture small, since the full matrix is ~170 KB. + cov.fixed = diag(diag(full$sdrep$cov.fixed)) + ), + quantities = full$quantities[intersect(quantities_keep, names(full$quantities))], + estimated_params = full$estimated_params[ + intersect("index_log_q", names(full$estimated_params)) + ] + ), + class = "Rceattle" +) + +saveRDS( + fixture, + testthat::test_path("fixtures", "rceattle_goa_atf.rds"), + compress = "xz" +) diff --git a/tests/testthat/fixtures/rceattle_goa_atf.rds b/tests/testthat/fixtures/rceattle_goa_atf.rds new file mode 100644 index 00000000..ec103552 Binary files /dev/null and b/tests/testthat/fixtures/rceattle_goa_atf.rds differ diff --git a/tests/testthat/test-convert_output.R b/tests/testthat/test-convert_output.R index 1a546672..6e90889d 100644 --- a/tests/testthat/test-convert_output.R +++ b/tests/testthat/test-convert_output.R @@ -129,3 +129,172 @@ test_that("replace_empty_with_na_all replaces empty strings only", { expect_equal(out$b, c(NA, "y", NA)) expect_equal(out$c, 1:3) }) + +test_that("convert_output works for Rceattle", { + fit_path <- test_path("fixtures", "rceattle_goa_atf.rds") + + # Both documented input modes: a path, and an already-loaded fit. + expect_no_error(from_path <- convert_output( + file = fit_path, + model = "rceattle" + )) + expect_no_error(from_object <- convert_output( + file = readRDS(fit_path), + model = "rceattle" + )) + + expect_equal(dim(from_path)[2], 35) + expect_equal(from_path, from_object) + + # Every row says which module it came from. + expect_false(any(is.na(from_path[["module_name"]]))) + + # Predicted composition binds column-wise to the observations, one row of + # age_hat per row of comp_data, so reading age_hat row by row reproduces the + # series exactly. A cross join repeats it once per row of comp_data instead. + fit <- readRDS(fit_path) + age_hat <- as.vector(t(fit[["quantities"]][["age_hat"]])) + expect_equal( + as.numeric(from_path[from_path[["label"]] == "composition_predicted", ][["estimate"]]), + age_hat[age_hat != 0] + ) + + # The index series are the fit's own values, in order. + expect_equal( + as.numeric(from_path[from_path[["label"]] == "index_predicted", ][["estimate"]]), + unname(as.numeric(fit[["quantities"]][["index_hat"]])) + ) + + # A derived time series carries its years, rather than falling into the + # `year = NA` branch. + biomass <- from_path[from_path[["label"]] == "biomass", ] + expect_gt(nrow(biomass), 0) + expect_false(any(is.na(biomass[["year"]]))) + expect_true(all(is.finite(as.numeric(biomass[["estimate"]])))) +}) + +test_that("convert_output places Rceattle selectivity and M at age", { + fit <- readRDS(test_path("fixtures", "rceattle_goa_atf.rds")) + out <- convert_output(file = fit, model = "rceattle") + + # One row per fleet, sex, age and year of sel_at_age, and each cell lands + # on its own fleet, sex, age and year. + selectivity <- out[out[["label"]] == "selectivity_age", ] + expect_equal(nrow(selectivity), length(fit[["quantities"]][["sel_at_age"]])) + fishery_age_5 <- selectivity |> + dplyr::filter( + fleet == "GOA_atf_fishery", + sex == "male", + age == "5", + year == 2018 + ) + expect_equal( + fishery_age_5[["estimate"]], + unname(fit[["quantities"]][["sel_at_age"]]["GOA_atf_fishery", 2, 5, "2018"]) + ) + + # GOAatf fixes M at 0.2 for females and 0.35 for males at every age and year. + natural_mortality <- out[out[["label"]] == "natural_mortality", ] + expect_equal(nrow(natural_mortality), 2 * 21 * 90) + expect_equal( + unique(natural_mortality[natural_mortality[["sex"]] == "female", ][["estimate"]]), + 0.2 + ) + expect_equal( + unique(natural_mortality[natural_mortality[["sex"]] == "male", ][["estimate"]]), + 0.35 + ) + + # The age dimension is the bin index, so the ages start at minage. + fit[["data_list"]][["minage"]] <- 3 + shifted <- convert_output(file = fit, model = "rceattle") + expect_equal( + range(as.numeric(shifted[shifted[["label"]] == "natural_mortality", ][["age"]])), + c(3, 23) + ) + + expect_no_error(plot_selectivity(out, interactive = FALSE)) + expect_no_error(plot_natural_mortality(out, interactive = FALSE)) +}) + +test_that("convert_output gives Rceattle recruitment deviations a standard error", { + fit <- readRDS(test_path("fixtures", "rceattle_goa_atf.rds")) + out <- convert_output(file = fit, model = "rceattle") + + # The standard error of a fixed effect is the square root of its diagonal + # entry in cov.fixed. + rec_dev_sd <- sqrt(diag(fit[["sdrep"]][["cov.fixed"]]))[ + names(fit[["sdrep"]][["par.fixed"]]) == "rec_dev" + ] + rec_devs <- out[out[["label"]] == "recruitment_deviations", ] + expect_equal(rec_devs[["uncertainty"]], rec_dev_sd) + expect_true(all(rec_devs[["uncertainty_label"]] == "sd")) + + expect_no_error(plot_recruitment_deviations(out, interactive = FALSE)) + + # Under random_rec = TRUE the deviations are random effects, held in + # par.random with their variances in diag.cov.random. + is_rec_dev <- names(fit[["sdrep"]][["par.fixed"]]) == "rec_dev" + random_rec <- fit + random_rec[["sdrep"]][["par.random"]] <- fit[["sdrep"]][["par.fixed"]][is_rec_dev] + random_rec[["sdrep"]][["diag.cov.random"]] <- rec_dev_sd^2 + random_rec[["sdrep"]][["par.fixed"]] <- fit[["sdrep"]][["par.fixed"]][!is_rec_dev] + random_rec[["sdrep"]][["cov.fixed"]] <- + fit[["sdrep"]][["cov.fixed"]][!is_rec_dev, !is_rec_dev] + random_out <- convert_output(file = random_rec, model = "rceattle") + random_rec_devs <- random_out[random_out[["label"]] == "recruitment_deviations", ] + expect_equal(random_rec_devs[["estimate"]], rec_devs[["estimate"]]) + expect_equal(random_rec_devs[["uncertainty"]], rec_dev_sd) + expect_equal(random_rec_devs[["year"]], rec_devs[["year"]]) +}) + +test_that("convert_output accepts .RDS as well as .rds for Rceattle", { + # saveRDS() output is conventionally named either way, and real Rceattle + # assessment files use the upper-case form. + upper <- tempfile(fileext = ".RDS") + on.exit(unlink(upper), add = TRUE) + saveRDS(readRDS(test_path("fixtures", "rceattle_goa_atf.rds")), upper) + expect_no_error(convert_output(file = upper, model = "rceattle")) +}) + +test_that("convert_output says so when an Rceattle fit carries no sdrep", { + # getsd = FALSE leaves sdrep NULL, and every derived time series lives there, + # so the converted output holds only the input data. + no_sd <- readRDS(test_path("fixtures", "rceattle_goa_atf.rds")) + no_sd[["sdrep"]] <- NULL + + # Matched on a single word, since cli wraps a message and a phrase can be + # split across a line break. + expect_warning( + out <- convert_output(file = no_sd, model = "rceattle"), + "sdrep" + ) + expect_false("biomass" %in% out[["label"]]) + expect_true("catch" %in% out[["label"]]) +}) + +test_that("convert_output refuses an Rceattle fit of unknown species count", { + # The standardized output has no species dimension, so an unreadable species + # count must stop the conversion rather than risk merging two stocks. + no_nspp <- readRDS(test_path("fixtures", "rceattle_goa_atf.rds")) + no_nspp[["data_list"]][["nspp"]] <- NULL + + expect_error( + convert_output(file = no_nspp, model = "rceattle"), + "species count is unknown" + ) +}) + +test_that("convert_output refuses a multispecies Rceattle fit", { + # The standardized output describes one stock and has no species dimension. + # CEATTLE is a multispecies model, so this must fail rather than return rows + # in which two stocks share a year. + multi <- readRDS(test_path("fixtures", "rceattle_goa_atf.rds")) + multi[["data_list"]][["nspp"]] <- 3L + + # Matched on a phrase short enough that cli's wrapping cannot split it. + expect_error( + convert_output(file = multi, model = "rceattle"), + "Rceattle model has 3 species" + ) +})