library(unittest) library(gadget3) read.matrix <- function (x, ...) unname(as.matrix(read.table(text=x, ...))) teststock <- g3_stock('teststock', seq(10, 35, 5)) ok_group("g3a_grow_lengthvbsimple", { # Run g3a_grow_lengthvbsimple for given linf and midlen lvs <- function (linf, kappa = 100) { as.numeric(g3_eval( g3a_grow_lengthvbsimple(linf, kappa), stock = g3_stock("s", seq(10, 100, 10) - 5), cur_step_size = 1)) } ok(ut_cmp_equal( lvs(linf = 150), c(140, 130, 120, 110, 100, 90, 80, 70, 60, 50)), "linf = 150") ok(ut_cmp_equal( lvs(linf = 55), c(45, 35, 25, 15, 5, 0, 0, 0, 0, 0), tolerance = 1e-7), "linf = 55 (linf doesn't fall negative)") }) ok_group("g3a_grow_impl_bbinom:zero-growth", { # Run g3a_grow_impl_bbinom for given linf and midlen gim <- function (linf, kappa = 20, beta = 0.5, maxlengthgroupgrowth = 10) { g3_eval( g3a_grow_impl_bbinom(g3a_grow_lengthvbsimple(linf, kappa), ~0, beta, maxlengthgroupgrowth)$len, stock = g3_stock("s", seq(10, 100, 10) - 5), cur_step_size = 1) } ok(ut_cmp_equal(round(gim(linf = 120, maxlengthgroupgrowth = 5), 2), read.matrix(' 8.29 8.44 0.40 0.17 0.12 0.15 9.00 10.00 0.00 0.00 0.00 0.00 11.43 14.29 1.02 0.36 0.22 0.26 26.56 39.35 7.50 2.00 1.11 1.18 26.18 50.91 21.82 2.18 0.91 0.82 21.00 70.00 80.00 32.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 1.00 0.02 0.05 0.08 0.13 0.21 0.51 0.14 0.12 0.12 0.13 0.17 0.32 0.37 0.14 0.10 0.10 0.11 0.19 ')),"linf=120, maxlengthgroupgrowth=5") # TODO: rowSums(gim(linf = 120, maxlengthgroupgrowth = 5)) should be 1? ok(ut_cmp_equal(round(gim(linf = 40, maxlengthgroupgrowth = 5), 2), read.matrix(' 0.14 0.12 0.12 0.13 0.17 0.32 0.37 0.14 0.10 0.10 0.11 0.19 0.66 0.09 0.06 0.05 0.05 0.09 1.00 0.00 0.00 0.00 0.00 0.00 1.00 0.00 0.00 0.00 0.00 0.00 1.00 0.00 0.00 0.00 0.00 0.00 1.00 0.00 0.00 0.00 0.00 0.00 1.00 0.00 0.00 0.00 0.00 0.00 1.00 0.00 0.00 0.00 0.00 0.00 1.00 0.00 0.00 0.00 0.00 0.00 ')),"linf=40, maxlengthgroupgrowth=5 - Linf lower than Lmax produces valid output") ok(ut_cmp_equal( rowSums(gim(linf = 40, maxlengthgroupgrowth = 5)), rep(1, 10), tolerance = 1e-2), "linf=40, maxlengthgroupgrowth=5 - Rows sum to 1") }) ok_group("g3a_grow_impl_bbinom", { actions <- list( # dmu, lengthgrouplen, binn, beta g3a_time(2000, 2001, project_years = 0), g3a_initialconditions(teststock, ~g3_param_vector("initial", value = c(10, 100, 1000, 1000, 10000, 100000)), ~0 * teststock__minlen), g3a_growmature(teststock, impl_f = g3a_grow_impl_bbinom( delta_len_f = g3a_grow_lengthvbsimple(~g3_param('linf', value = 160), ~g3_param('kappa', value = 90)), delta_wgt_f = g3a_grow_weightsimple(~g3_param('walpha', value = 3), ~g3_param('wbeta', value = 2)), beta = ~g3_param('beta', value = 30), maxlengthgroupgrowth = 4)), list( "999" = ~{ REPORT(teststock__growth_l) REPORT(teststock__growth_w) nll <- nll + g3_param('x', value = 1.0) return(nll) })) # Compile model # NB: Growth not valid, but not testing growth at this point model_fn <- g3_to_r(actions, trace = FALSE, strict = FALSE) # model_fn <- edit(model_fn) if (nzchar(Sys.getenv('G3_TEST_TMB'))) { model_cpp <- g3_to_tmb(actions, trace = FALSE, strict = FALSE) # model_cpp <- edit(model_cpp) model_tmb <- g3_tmb_adfun(model_cpp, attr(model_fn, 'parameter_template'), compile_flags = c("-O0", "-g")) } else { writeLines("# skip: not compiling TMB model") } params <- attr(model_fn, 'parameter_template') result <- model_fn(params) r <- attributes(result) # str(result) # str(as.list(r), vec.len = 10000) #writeLines(deparse(r$teststock__growth_l)) ok(ut_cmp_equal(r$teststock__growth_l, structure(c(12199.8976226175, 9352.23228375502, 7157.83523789257, 5463.49450549465, 4154.31593595054, 3143.22179218075, 51322.2074676947, 39560.4630926918, 30458.8733527337, 23399.2087912083, 17917.1342692155, 13660.121314133, 81087.2009530949, 62860.2639001552, 48695.7187344743, 37658.1016483506, 29043.7267350904, 22317.4314305221, 57032.8699935755, 44472.5676572519, 34669.6583623158, 26995.0549450546, 20974.8144063369, 16247.8860873661, 15068.9788855573, 11821.5345660349, 9275.97774268346, 7273.66758241751, 5694.90600451131, 4448.35417879727), .Dimnames = list(length = c("10:15", "15:20", "20:25", "25:30", "30:35", "35:Inf"), delta = c("0", "1", "2", "3", "4")), .Dim = c(length = 6, delta = 5)), tolerance = 1e-5), "Matches baseline") 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) } # Try comparing with a few different inputs for (x in 1:10) ok_group("g3a_grow_impl_bbinom random params", { params <- attr(model_fn, 'parameter_template') params$linf <- runif(1, 100, 200) params$kappa <- runif(1, 50, 100) params$walpha <- runif(1, 0.1, 3) params$wbeta <- runif(1, 0.1, 2) params$beta <- runif(1, 10, 30) params$initial <- runif(6, 1000, 9000) 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) } }) }) ok_group("g3a_growmature", { actions <- list( g3a_time(2000, 2001, project_years = 0), g3a_initialconditions(teststock, ~g3_param_vector("initial_num"), ~g3_param_vector("initial_wgt")), g3a_growmature(teststock, impl_f = list( delta_dim = as.character(0:6), len = ~g3_param_array('growth_matrix'), wgt = ~g3_param_array('weight_matrix'))), list( "999" = ~{ REPORT(teststock__num) REPORT(teststock__wgt) nll <- nll + g3_param('x') return(nll) })) params <- attr(model_fn, 'parameter_template') params$initial_num <- c(10, 100, 1000, 1000, 10000, 100000) params$initial_wgt <- c(100, 200, 300, 400, 500, 600) params$growth_matrix <- array(dim = c(6, 7)) params$weight_matrix <- array(0, dim = c(6, 7)) # Compile model # NB: Our tests don't preserve counts, so we can't enable strict here model_fn <- g3_to_r(actions, trace = FALSE, strict = FALSE) # model_fn <- edit(model_fn) if (nzchar(Sys.getenv('G3_TEST_TMB'))) { model_cpp <- g3_to_tmb(actions, trace = FALSE, strict = 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") } gm <- array(c( # 10 15 20 25 30 35 0, 0.25, 1, 0.5, 0.5, 1, # 0 0, 0, 0, 0, 0, 0, # +1 1, 0.75, 0, 0, 0, 0, # +2 0, 0, 0, 0, 0, 0, # +3 0, 0, 0, 0, 0.5, 0, # +4 0, 0, 0, 0, 0, 0, # +5 0, 0, 0, 0.5, 0, 0), dim = c(6,7)) params <- attr(model_fn, 'parameter_template') params$initial_num <- c(10, 100, 1000, 1000, 10000, 100000) params$initial_wgt <- c(100, 200, 300, 400, 500, 600) params$growth_matrix <- gm params$weight_matrix <- array(0, dim = c(6, 7)) result <- model_fn(params) r <- attributes(result) # str(result) # str(as.list(r), vec.len = 10000) # model_fn <- edit(model_fn) result <- model_fn(params) ok(ut_cmp_identical( as.vector(r$teststock__num), c(0, 25, 1010, 575, 5000, 105500)), "Stock individuals have been scaled by matrix, can't escape plus-group") ok(ut_cmp_equal(as.vector(r$teststock__wgt), c( 0 / g3_env$logspace_add(0, 0), (200 * 100*0.25) / 25, (300 * 1000 + 100 * 10) / 1010, (400 * 1000*0.5 + 200 * 100*0.75) / 575, (500 * 10000*0.5) / 5000, (600 * 100000 + 500 * 10000*0.5 + 400 * 1000*0.5) / 105500, NULL), tolerance = 1e-5), "Weight scaled, didn't let weight go to infinity when dividing by zero") 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) } })