tests/test-likelihood_sparsesample.R
In gadget-framework/gadget3: Globally-Applicable Area Disaggregated General Ecosystem Toolbox V3
if (!interactive()) options(warn=2, error = function() { sink(stderr()) ; traceback(3) ; q(status = 1) })
library(unittest)
library(gadget3)
ok_group("g3s_sparsedata:area", {
sd <- gadget3:::g3s_sparsedata("frank", data.frame(
year = 1994,
area = c("a", "a", "c"),
length = 2,
stringsAsFactors = FALSE), area_group = c(a=1,b=2,c=3))
model_fn <- g3_to_r(list(g3a_time(1990, 1999), g3_step(g3_formula(
stock_iterate(sd, {
print(c(year = cur_year, length = length, area = area, sd__i = stock_ss(sd__i, vec = single)))
}),
sd__i = gadget3:::g3_sparsedata_instance(sd, 50.1 * seq_along(g3_stock_def(sd, 'year'))),
sd = sd ))))
ok(ut_cmp_identical(capture.output(model_fn()), c(
" year length area sd__i ",
"1994.0 2.0 1.0 50.1 ",
" year length area sd__i ",
"1994.0 2.0 1.0 100.2 ",
" year length area sd__i ",
"1994.0 2.0 3.0 150.3 ",
"[1] 0")), "model_fn, iterated over numeric areas")
})
st <- g3_stock("stst", c(10, 20, 30)) |> g3s_age(3,5)
fl <- g3_fleet(c("fl", "surv"))
obs_linreg_df <- data.frame(
# NB: No 1993, just ignored
year = rep(c(1990, 1991, 1992, 1994), each = 2),
step = 1:2 )
obs_linreg_df$age <- floor(runif(nrow(obs_linreg_df), min = 3, max = 5.1))
obs_linreg_df$length <- floor(runif(nrow(obs_linreg_df), min = 10, max = 50))
obs_linreg_df$mean <- runif(nrow(obs_linreg_df), min = 10, max = 1000)
# Occasionally test with integer values
if (runif(1) < 0.5) obs_linreg_df$mean <- as.integer(obs_linreg_df$mean)
obs_ss_df <- data.frame(
# NB: No step, just year
year = c(1990, 1992, 1993) )
obs_ss_df$mean = runif(nrow(obs_ss_df), min = 10, max = 1000)
obs_flc_df <- data.frame(
year = rep(c(1990, 1991, 1992, 1994), each = 2) )
obs_flc_df$length <- floor(runif(nrow(obs_linreg_df), min = 10, max = 50))
obs_flc_df$mean <- runif(nrow(obs_flc_df), min = 10, max = 1000)
obs_flc_df$stddev <- runif(nrow(obs_flc_df), min = 0.01, max = 0.1)
obs_flc_df$number <- as.integer(runif(nrow(obs_flc_df), min = 10, max = 10))
actions <- list(
g3a_time(1990, 1994, c(6,6)),
g3a_otherfood(st,
quote( age * 100 + stock__minlen ),
quote( cur_year * 1e5 + cur_step * 1e4 + 0 * stock__minlen ) ),
g3a_predate(
fl,
list(st),
suitabilities = g3_suitability_exponentiall50(),
catchability_f = g3a_predate_catchability_numberfleet(g3_parameterized("predate_num", value = 0)) ),
g3l_sparsesample(
"bt",
obs_linreg_df,
list(st),
measurement_f = g3_formula(
wgt + length,
end = NULL ),
function_f = g3l_sparsesample_linreg(fit = "linear") ),
g3l_sparsesample(
"cs_model",
obs_ss_df,
list(st),
measurement_f = g3_formula(
wgt * age,
end = NULL ),
function_f = g3l_sparsesample_sumsquares(weighting = 'model_stddev') ),
g3l_sparsesample(
"flc",
obs_flc_df,
list(st),
predstocks = list(fl),
measurement_f = quote(wgt + length),
function_f = g3l_sparsesample_sumsquares(weighting = 'obs_stddev') ),
# NB: Dummy parameter so model will compile in TMB
quote( nll <- nll + g3_param("x", value = 0) ) )
full_actions <- c(actions, list(
g3a_report_detail(actions),
g3a_report_history(actions, var_re = "__num$|__wgt$|__cons$|^nll_.sparse_.+__(model|obs)_.+$"),
g3a_report_history(actions, var_re = "^nll_.sparse_.+__nll$", out_prefix = NULL),
NULL))
model_fn <- g3_to_r(full_actions)
model_cpp <- g3_to_tmb(full_actions)
attr(model_fn, 'parameter_template') |>
g3_init_val("asparse_linreg_bt_weight", runif(1, 0.5, 5)) |>
g3_init_val("asparse_sumsquares_cs_model_weight", runif(1, 0.5, 5)) |>
g3_init_val("csparse_sumsquares_flc_weight", 0) |>
identity() -> params
nll <- model_fn(params) ; r <- attributes(nll) ; nll <- as.vector(nll)
# Final result can be derived from stock abundance
expected_m <- r$hist_stst__wgt + c(15, 25, 35)
length_grp <- as.character(cut(obs_linreg_df$length, c(9, 20, 30, Inf), right = FALSE, labels = c("10:20", "20:30", "30:Inf")))
for (i in seq_len(nrow(obs_linreg_df))) {
m <- expected_m[
length = length_grp[[i]],
age = sprintf("age%d", obs_linreg_df[i, "age"]),
time = sprintf("%04d-%02d", obs_linreg_df[i, "year"], obs_linreg_df[i, "step"]) ]
count <- r$hist_stst__num[
length = length_grp[[i]],
age = sprintf("age%d", obs_linreg_df[i, "age"]),
time = sprintf("%04d-%02d", obs_linreg_df[i, "year"], obs_linreg_df[i, "step"]) ]
ok(ut_cmp_equal(
as.vector(r$hist_nll_asparse_linreg_bt__model_sum[i, '1994-02']),
sum(m * count) ), paste0("r$hist_nll_asparse_linreg_bt__model_sum[", i, ", '1994-02']"))
ok(ut_cmp_equal(
as.vector(r$hist_nll_asparse_linreg_bt__model_sqsum[i, '1994-02']),
sum(m**2 * count) ), paste0("r$hist_nll_asparse_linreg_bt__model_sqsum[", i, ", '1994-02']"))
ok(ut_cmp_equal(
as.vector(r$hist_nll_asparse_linreg_bt__model_n[i, '1994-02']),
sum(count) ), paste0("r$hist_nll_asparse_linreg_bt__model_n[", i, ", '1994-02']"))
}
ok(ut_cmp_equal(
r$hist_nll_asparse_linreg_bt__model_sum,
ifelse(r$hist_nll_asparse_linreg_bt__model_n > 0, r$hist_nll_asparse_linreg_bt__model_sum[,"1994-02"], 0),
end = NULL), "hist_nll_asparse_linreg_bt__model_sum: History is the same shape as __model_n")
ok(gadget3:::ut_cmp_df(as.data.frame(r$hist_nll_asparse_linreg_bt__model_n > 0), '
1990-01 1990-02 1991-01 1991-02 1992-01 1992-02 1993-01 1993-02 1994-01 1994-02
TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
FALSE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
FALSE FALSE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
FALSE FALSE FALSE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
FALSE FALSE FALSE FALSE TRUE TRUE TRUE TRUE TRUE TRUE
FALSE FALSE FALSE FALSE FALSE TRUE TRUE TRUE TRUE TRUE
FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE TRUE
FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE
'), "hist_nll_asparse_linreg_bt__model_n: Added value once per step, apart from 1993")
ok(ut_cmp_equal(
r$hist_nll_asparse_linreg_bt__obs_mean,
array(
obs_linreg_df$mean,
dim = c(vec = length(obs_linreg_df$mean), time = 10),
dimnames = list(vec = NULL, time = paste0(rep(1990:1994, each=2), c("-01", "-02"))) ),
end = NULL), "hist_nll_asparse_linreg_bt__obs_mean: Repeats the input observations")
nll_mean <- r$hist_nll_asparse_linreg_bt__model_sum[, '1994-02']/r$hist_nll_asparse_linreg_bt__model_n[, '1994-02']
ok(ut_cmp_equal(
r$nll_asparse_linreg_bt__nll,
gadget3:::regression_linear(nll_mean, obs_linreg_df$mean, NaN, 1) ), "r$nll_asparse_linreg_bt__nll: Matches derived version")
# cs_model
ok(ut_cmp_equal(colSums(r$hist_nll_asparse_sumsquares_cs_model__model_n) / colSums(r$hist_stst__num, dims = 2), c(
"1990-01" = 1,
"1990-02" = 2,
"1991-01" = 2,
"1991-02" = 2,
"1992-01" = 3,
"1992-02" = 4,
"1993-01" = 5,
"1993-02" = 6,
"1994-01" = 6,
"1994-02" = 6,
NULL )), "hist_nll_asparse_sumsquares_cs_model__model_n: Final data point has 6 step's worth of individuals")
expected_m <- r$hist_stst__wgt * rep(3:5, each = 3)
for (i in seq_len(nrow(obs_ss_df))) {
# Pick out all relevant measurements to this row
m <- expected_m[,,time = paste0(obs_ss_df[i, "year"], c("-01", "-02"))]
count <- r$hist_stst__num[,,time = paste0(obs_ss_df[i, "year"], c("-01", "-02"))]
ok(ut_cmp_equal(
as.vector(r$hist_nll_asparse_sumsquares_cs_model__model_sum[i, '1994-02']),
sum(m * count) ), paste0("r$hist_nll_asparse_sumsquares_cs_model__model_sum[", i, ", '1994-02']"))
ok(ut_cmp_equal(
as.vector(r$hist_nll_asparse_sumsquares_cs_model__model_sqsum[i, '1994-02']),
sum(m**2 * count) ), paste0("r$hist_nll_asparse_sumsquares_cs_model__model_sqsum[", i, ", '1994-02']"))
ok(ut_cmp_equal(
as.vector(r$hist_nll_asparse_sumsquares_cs_model__model_n[i, '1994-02']),
sum(count) ), paste0("r$hist_nll_asparse_cs__model_n[", i, ", '1994-02']"))
}
nll_var <- r$hist_nll_asparse_sumsquares_cs_model__model_sqsum[, '1994-02']/r$hist_nll_asparse_sumsquares_cs_model__model_n[, '1994-02'] -
(r$hist_nll_asparse_sumsquares_cs_model__model_sum[, '1994-02']/r$hist_nll_asparse_sumsquares_cs_model__model_n[, '1994-02'])**2
nll_mean <- r$hist_nll_asparse_sumsquares_cs_model__model_sum[, '1994-02']/r$hist_nll_asparse_sumsquares_cs_model__model_n[, '1994-02']
ok(ut_cmp_equal(
(r$nll_asparse_sumsquares_cs_model__nll),
array((nll_mean - obs_ss_df$mean)^2 * 1/nll_var, dim = c(row = 3L), dimnames = list(row = NULL)) ), "r$nll_asparse_sumsquares_cs_model__nll: Matches derived version")
# Overall nll
ok(ut_cmp_equal(nll,
params$asparse_linreg_bt_weight * r$nll_asparse_linreg_bt__nll[["nll"]] +
params$asparse_sumsquares_cs_model_weight * sum(r$nll_asparse_sumsquares_cs_model__nll) +
params$csparse_sumsquares_flc_weight * sum(r$nll_csparse_sumsquares_flc__nll) +
0 ), "nll: Overall value matches")
for (n in grep("^hist_nll_", names(r), value = TRUE)) {
ok(ut_cmp_equal(
as.vector(r[[n]][, "1994-02"]),
as.vector(r[[gsub("^hist_", "", n)]])), paste0(n, ": _detail reports match last column of our history") )
}
gadget3:::ut_tmb_r_compare2(model_fn, model_cpp, params)
ok_group("flc") ###############################################################
attr(model_fn, 'parameter_template') |>
g3_init_val("asparse_linreg_bt_weight", 0) |>
g3_init_val("asparse_sumsquares_cs_model_weight", 0) |>
g3_init_val("csparse_sumsquares_flc_weight", runif(1, 0.5, 5)) |>
g3_init_val("predate_num", 100) |>
g3_init_val("stst.fl_surv.alpha", 10) |>
g3_init_val("stst.fl_surv.l50", 25) |>
identity() -> params
nll <- model_fn(params) ; r <- attributes(nll) ; nll <- as.vector(nll)
ok(gadget3:::ut_cmp_df(as.data.frame((r$hist_stst_fl_surv__cons / r$hist_stst__wgt)[,,1]), '
age3 age4 age5
10:20 0.00000 0.00000 0.00000
20:30 8.33333 10.93750 13.54167
30:Inf 17.18750 22.39583 27.60417
', tolerance = 1e-5), "hist_stst_fl_surv__cons: Caught 100 fish, with a selectivity pattern")
expected_m <- r$hist_stst__wgt + c(15, 25, 35)
length_grp <- as.character(cut(obs_flc_df$length, c(9, 20, 30, Inf), right = FALSE, labels = c("10:20", "20:30", "30:Inf")))
for (i in seq_len(nrow(obs_ss_df))) {
m <- expected_m[
length = length_grp[[i]],
age = ,
time = sprintf("%04d-%02d", obs_flc_df[i, "year"], 1:2) ]
count <- (r$hist_stst_fl_surv__cons / r$hist_stst__wgt)[
length = length_grp[[i]],
age = ,
time = sprintf("%04d-%02d", obs_flc_df[i, "year"], 1:2) ]
ok(ut_cmp_equal(
as.vector(r$hist_nll_csparse_sumsquares_flc__model_n[i, '1994-02']),
sum(count) ), paste0("r$hist_nll_csparse_sumsquares_flc__model_n[", i, ", '1994-02']"))
ok(ut_cmp_equal(
as.vector(r$hist_nll_csparse_sumsquares_flc__model_sum[i, '1994-02']),
sum(m * count) ), paste0("r$hist_nll_csparse_sumsquares_flc__model_sum[", i, ", '1994-02']"))
ok(ut_cmp_equal(
as.vector(r$hist_nll_csparse_sumsquares_flc__model_sqsum[i, '1994-02']),
sum(m**2 * count) ), paste0("r$hist_nll_csparse_sumsquares_flc__model_sqsum[", i, ", '1994-02']"))
}
nll_var <- obs_flc_df$stddev^2 # Using weighting = 'obs_stddev'
nll_mean <- r$hist_nll_csparse_sumsquares_flc__model_sum[, '1994-02']/r$hist_nll_csparse_sumsquares_flc__model_n[, '1994-02']
ok(ut_cmp_equal(
r$nll_csparse_sumsquares_flc__nll,
array((nll_mean - obs_flc_df$mean)^2 * (1/nll_var) * obs_flc_df$number, dim = c(row = length(r$nll_csparse_sumsquares_flc__nll)), dimnames = list(row = NULL)) ), "r$nll_csparse_sumsquares_flc__nll: Matches derived version")
# Overall nll
ok(ut_cmp_equal(nll,
params$asparse_linreg_bt_weight * r$nll_asparse_linreg_bt__nll[["nll"]] +
params$asparse_sumsquares_cs_model_weight * sum(r$nll_asparse_sumsquares_cs_model__nll) +
params$csparse_sumsquares_flc_weight * sum(r$nll_csparse_sumsquares_flc__nll) +
0 ), "nll: Overall value matches")
for (n in grep("^hist_nll_", names(r), value = TRUE)) {
ok(ut_cmp_equal(
as.vector(r[[n]][, "1994-02"]),
as.vector(r[[gsub("^hist_", "", n)]])), paste0(n, ": _detail reports match last column of our history") )
}
gadget3:::ut_tmb_r_compare2(model_fn, model_cpp, params)
gadget-framework/gadget3 documentation built on June 13, 2025, 5:06 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
if (!interactive()) options(warn=2, error = function() { sink(stderr()) ; traceback(3) ; q(status = 1) })
library(unittest)
library(gadget3)
ok_group("g3s_sparsedata:area", {
sd <- gadget3:::g3s_sparsedata("frank", data.frame(
year = 1994,
area = c("a", "a", "c"),
length = 2,
stringsAsFactors = FALSE), area_group = c(a=1,b=2,c=3))
model_fn <- g3_to_r(list(g3a_time(1990, 1999), g3_step(g3_formula(
stock_iterate(sd, {
print(c(year = cur_year, length = length, area = area, sd__i = stock_ss(sd__i, vec = single)))
}),
sd__i = gadget3:::g3_sparsedata_instance(sd, 50.1 * seq_along(g3_stock_def(sd, 'year'))),
sd = sd ))))
ok(ut_cmp_identical(capture.output(model_fn()), c(
" year length area sd__i ",
"1994.0 2.0 1.0 50.1 ",
" year length area sd__i ",
"1994.0 2.0 1.0 100.2 ",
" year length area sd__i ",
"1994.0 2.0 3.0 150.3 ",
"[1] 0")), "model_fn, iterated over numeric areas")
})
st <- g3_stock("stst", c(10, 20, 30)) |> g3s_age(3,5)
fl <- g3_fleet(c("fl", "surv"))
obs_linreg_df <- data.frame(
# NB: No 1993, just ignored
year = rep(c(1990, 1991, 1992, 1994), each = 2),
step = 1:2 )
obs_linreg_df$age <- floor(runif(nrow(obs_linreg_df), min = 3, max = 5.1))
obs_linreg_df$length <- floor(runif(nrow(obs_linreg_df), min = 10, max = 50))
obs_linreg_df$mean <- runif(nrow(obs_linreg_df), min = 10, max = 1000)
# Occasionally test with integer values
if (runif(1) < 0.5) obs_linreg_df$mean <- as.integer(obs_linreg_df$mean)
obs_ss_df <- data.frame(
# NB: No step, just year
year = c(1990, 1992, 1993) )
obs_ss_df$mean = runif(nrow(obs_ss_df), min = 10, max = 1000)
obs_flc_df <- data.frame(
year = rep(c(1990, 1991, 1992, 1994), each = 2) )
obs_flc_df$length <- floor(runif(nrow(obs_linreg_df), min = 10, max = 50))
obs_flc_df$mean <- runif(nrow(obs_flc_df), min = 10, max = 1000)
obs_flc_df$stddev <- runif(nrow(obs_flc_df), min = 0.01, max = 0.1)
obs_flc_df$number <- as.integer(runif(nrow(obs_flc_df), min = 10, max = 10))
actions <- list(
g3a_time(1990, 1994, c(6,6)),
g3a_otherfood(st,
quote( age * 100 + stock__minlen ),
quote( cur_year * 1e5 + cur_step * 1e4 + 0 * stock__minlen ) ),
g3a_predate(
fl,
list(st),
suitabilities = g3_suitability_exponentiall50(),
catchability_f = g3a_predate_catchability_numberfleet(g3_parameterized("predate_num", value = 0)) ),
g3l_sparsesample(
"bt",
obs_linreg_df,
list(st),
measurement_f = g3_formula(
wgt + length,
end = NULL ),
function_f = g3l_sparsesample_linreg(fit = "linear") ),
g3l_sparsesample(
"cs_model",
obs_ss_df,
list(st),
measurement_f = g3_formula(
wgt * age,
end = NULL ),
function_f = g3l_sparsesample_sumsquares(weighting = 'model_stddev') ),
g3l_sparsesample(
"flc",
obs_flc_df,
list(st),
predstocks = list(fl),
measurement_f = quote(wgt + length),
function_f = g3l_sparsesample_sumsquares(weighting = 'obs_stddev') ),
# NB: Dummy parameter so model will compile in TMB
quote( nll <- nll + g3_param("x", value = 0) ) )
full_actions <- c(actions, list(
g3a_report_detail(actions),
g3a_report_history(actions, var_re = "__num$|__wgt$|__cons$|^nll_.sparse_.+__(model|obs)_.+$"),
g3a_report_history(actions, var_re = "^nll_.sparse_.+__nll$", out_prefix = NULL),
NULL))
model_fn <- g3_to_r(full_actions)
model_cpp <- g3_to_tmb(full_actions)
attr(model_fn, 'parameter_template') |>
g3_init_val("asparse_linreg_bt_weight", runif(1, 0.5, 5)) |>
g3_init_val("asparse_sumsquares_cs_model_weight", runif(1, 0.5, 5)) |>
g3_init_val("csparse_sumsquares_flc_weight", 0) |>
identity() -> params
nll <- model_fn(params) ; r <- attributes(nll) ; nll <- as.vector(nll)
# Final result can be derived from stock abundance
expected_m <- r$hist_stst__wgt + c(15, 25, 35)
length_grp <- as.character(cut(obs_linreg_df$length, c(9, 20, 30, Inf), right = FALSE, labels = c("10:20", "20:30", "30:Inf")))
for (i in seq_len(nrow(obs_linreg_df))) {
m <- expected_m[
length = length_grp[[i]],
age = sprintf("age%d", obs_linreg_df[i, "age"]),
time = sprintf("%04d-%02d", obs_linreg_df[i, "year"], obs_linreg_df[i, "step"]) ]
count <- r$hist_stst__num[
length = length_grp[[i]],
age = sprintf("age%d", obs_linreg_df[i, "age"]),
time = sprintf("%04d-%02d", obs_linreg_df[i, "year"], obs_linreg_df[i, "step"]) ]
ok(ut_cmp_equal(
as.vector(r$hist_nll_asparse_linreg_bt__model_sum[i, '1994-02']),
sum(m * count) ), paste0("r$hist_nll_asparse_linreg_bt__model_sum[", i, ", '1994-02']"))
ok(ut_cmp_equal(
as.vector(r$hist_nll_asparse_linreg_bt__model_sqsum[i, '1994-02']),
sum(m**2 * count) ), paste0("r$hist_nll_asparse_linreg_bt__model_sqsum[", i, ", '1994-02']"))
ok(ut_cmp_equal(
as.vector(r$hist_nll_asparse_linreg_bt__model_n[i, '1994-02']),
sum(count) ), paste0("r$hist_nll_asparse_linreg_bt__model_n[", i, ", '1994-02']"))
}
ok(ut_cmp_equal(
r$hist_nll_asparse_linreg_bt__model_sum,
ifelse(r$hist_nll_asparse_linreg_bt__model_n > 0, r$hist_nll_asparse_linreg_bt__model_sum[,"1994-02"], 0),
end = NULL), "hist_nll_asparse_linreg_bt__model_sum: History is the same shape as __model_n")
ok(gadget3:::ut_cmp_df(as.data.frame(r$hist_nll_asparse_linreg_bt__model_n > 0), '
1990-01 1990-02 1991-01 1991-02 1992-01 1992-02 1993-01 1993-02 1994-01 1994-02
TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
FALSE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
FALSE FALSE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
FALSE FALSE FALSE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
FALSE FALSE FALSE FALSE TRUE TRUE TRUE TRUE TRUE TRUE
FALSE FALSE FALSE FALSE FALSE TRUE TRUE TRUE TRUE TRUE
FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE TRUE
FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE
'), "hist_nll_asparse_linreg_bt__model_n: Added value once per step, apart from 1993")
ok(ut_cmp_equal(
r$hist_nll_asparse_linreg_bt__obs_mean,
array(
obs_linreg_df$mean,
dim = c(vec = length(obs_linreg_df$mean), time = 10),
dimnames = list(vec = NULL, time = paste0(rep(1990:1994, each=2), c("-01", "-02"))) ),
end = NULL), "hist_nll_asparse_linreg_bt__obs_mean: Repeats the input observations")
nll_mean <- r$hist_nll_asparse_linreg_bt__model_sum[, '1994-02']/r$hist_nll_asparse_linreg_bt__model_n[, '1994-02']
ok(ut_cmp_equal(
r$nll_asparse_linreg_bt__nll,
gadget3:::regression_linear(nll_mean, obs_linreg_df$mean, NaN, 1) ), "r$nll_asparse_linreg_bt__nll: Matches derived version")
# cs_model
ok(ut_cmp_equal(colSums(r$hist_nll_asparse_sumsquares_cs_model__model_n) / colSums(r$hist_stst__num, dims = 2), c(
"1990-01" = 1,
"1990-02" = 2,
"1991-01" = 2,
"1991-02" = 2,
"1992-01" = 3,
"1992-02" = 4,
"1993-01" = 5,
"1993-02" = 6,
"1994-01" = 6,
"1994-02" = 6,
NULL )), "hist_nll_asparse_sumsquares_cs_model__model_n: Final data point has 6 step's worth of individuals")
expected_m <- r$hist_stst__wgt * rep(3:5, each = 3)
for (i in seq_len(nrow(obs_ss_df))) {
# Pick out all relevant measurements to this row
m <- expected_m[,,time = paste0(obs_ss_df[i, "year"], c("-01", "-02"))]
count <- r$hist_stst__num[,,time = paste0(obs_ss_df[i, "year"], c("-01", "-02"))]
ok(ut_cmp_equal(
as.vector(r$hist_nll_asparse_sumsquares_cs_model__model_sum[i, '1994-02']),
sum(m * count) ), paste0("r$hist_nll_asparse_sumsquares_cs_model__model_sum[", i, ", '1994-02']"))
ok(ut_cmp_equal(
as.vector(r$hist_nll_asparse_sumsquares_cs_model__model_sqsum[i, '1994-02']),
sum(m**2 * count) ), paste0("r$hist_nll_asparse_sumsquares_cs_model__model_sqsum[", i, ", '1994-02']"))
ok(ut_cmp_equal(
as.vector(r$hist_nll_asparse_sumsquares_cs_model__model_n[i, '1994-02']),
sum(count) ), paste0("r$hist_nll_asparse_cs__model_n[", i, ", '1994-02']"))
}
nll_var <- r$hist_nll_asparse_sumsquares_cs_model__model_sqsum[, '1994-02']/r$hist_nll_asparse_sumsquares_cs_model__model_n[, '1994-02'] -
(r$hist_nll_asparse_sumsquares_cs_model__model_sum[, '1994-02']/r$hist_nll_asparse_sumsquares_cs_model__model_n[, '1994-02'])**2
nll_mean <- r$hist_nll_asparse_sumsquares_cs_model__model_sum[, '1994-02']/r$hist_nll_asparse_sumsquares_cs_model__model_n[, '1994-02']
ok(ut_cmp_equal(
(r$nll_asparse_sumsquares_cs_model__nll),
array((nll_mean - obs_ss_df$mean)^2 * 1/nll_var, dim = c(row = 3L), dimnames = list(row = NULL)) ), "r$nll_asparse_sumsquares_cs_model__nll: Matches derived version")
# Overall nll
ok(ut_cmp_equal(nll,
params$asparse_linreg_bt_weight * r$nll_asparse_linreg_bt__nll[["nll"]] +
params$asparse_sumsquares_cs_model_weight * sum(r$nll_asparse_sumsquares_cs_model__nll) +
params$csparse_sumsquares_flc_weight * sum(r$nll_csparse_sumsquares_flc__nll) +
0 ), "nll: Overall value matches")
for (n in grep("^hist_nll_", names(r), value = TRUE)) {
ok(ut_cmp_equal(
as.vector(r[[n]][, "1994-02"]),
as.vector(r[[gsub("^hist_", "", n)]])), paste0(n, ": _detail reports match last column of our history") )
}
gadget3:::ut_tmb_r_compare2(model_fn, model_cpp, params)
ok_group("flc") ###############################################################
attr(model_fn, 'parameter_template') |>
g3_init_val("asparse_linreg_bt_weight", 0) |>
g3_init_val("asparse_sumsquares_cs_model_weight", 0) |>
g3_init_val("csparse_sumsquares_flc_weight", runif(1, 0.5, 5)) |>
g3_init_val("predate_num", 100) |>
g3_init_val("stst.fl_surv.alpha", 10) |>
g3_init_val("stst.fl_surv.l50", 25) |>
identity() -> params
nll <- model_fn(params) ; r <- attributes(nll) ; nll <- as.vector(nll)
ok(gadget3:::ut_cmp_df(as.data.frame((r$hist_stst_fl_surv__cons / r$hist_stst__wgt)[,,1]), '
age3 age4 age5
10:20 0.00000 0.00000 0.00000
20:30 8.33333 10.93750 13.54167
30:Inf 17.18750 22.39583 27.60417
', tolerance = 1e-5), "hist_stst_fl_surv__cons: Caught 100 fish, with a selectivity pattern")
expected_m <- r$hist_stst__wgt + c(15, 25, 35)
length_grp <- as.character(cut(obs_flc_df$length, c(9, 20, 30, Inf), right = FALSE, labels = c("10:20", "20:30", "30:Inf")))
for (i in seq_len(nrow(obs_ss_df))) {
m <- expected_m[
length = length_grp[[i]],
age = ,
time = sprintf("%04d-%02d", obs_flc_df[i, "year"], 1:2) ]
count <- (r$hist_stst_fl_surv__cons / r$hist_stst__wgt)[
length = length_grp[[i]],
age = ,
time = sprintf("%04d-%02d", obs_flc_df[i, "year"], 1:2) ]
ok(ut_cmp_equal(
as.vector(r$hist_nll_csparse_sumsquares_flc__model_n[i, '1994-02']),
sum(count) ), paste0("r$hist_nll_csparse_sumsquares_flc__model_n[", i, ", '1994-02']"))
ok(ut_cmp_equal(
as.vector(r$hist_nll_csparse_sumsquares_flc__model_sum[i, '1994-02']),
sum(m * count) ), paste0("r$hist_nll_csparse_sumsquares_flc__model_sum[", i, ", '1994-02']"))
ok(ut_cmp_equal(
as.vector(r$hist_nll_csparse_sumsquares_flc__model_sqsum[i, '1994-02']),
sum(m**2 * count) ), paste0("r$hist_nll_csparse_sumsquares_flc__model_sqsum[", i, ", '1994-02']"))
}
nll_var <- obs_flc_df$stddev^2 # Using weighting = 'obs_stddev'
nll_mean <- r$hist_nll_csparse_sumsquares_flc__model_sum[, '1994-02']/r$hist_nll_csparse_sumsquares_flc__model_n[, '1994-02']
ok(ut_cmp_equal(
r$nll_csparse_sumsquares_flc__nll,
array((nll_mean - obs_flc_df$mean)^2 * (1/nll_var) * obs_flc_df$number, dim = c(row = length(r$nll_csparse_sumsquares_flc__nll)), dimnames = list(row = NULL)) ), "r$nll_csparse_sumsquares_flc__nll: Matches derived version")
# Overall nll
ok(ut_cmp_equal(nll,
params$asparse_linreg_bt_weight * r$nll_asparse_linreg_bt__nll[["nll"]] +
params$asparse_sumsquares_cs_model_weight * sum(r$nll_asparse_sumsquares_cs_model__nll) +
params$csparse_sumsquares_flc_weight * sum(r$nll_csparse_sumsquares_flc__nll) +
0 ), "nll: Overall value matches")
for (n in grep("^hist_nll_", names(r), value = TRUE)) {
ok(ut_cmp_equal(
as.vector(r[[n]][, "1994-02"]),
as.vector(r[[gsub("^hist_", "", n)]])), paste0(n, ": _detail reports match last column of our history") )
}
gadget3:::ut_tmb_r_compare2(model_fn, model_cpp, params)
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.