R Under development (unstable) (2024-07-01 r86857 ucrt) -- "Unsuffered Consequences" Copyright (C) 2024 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > library(magrittr) > library(unittest) > > library(gadget3) > > capture_warnings <- function(x, full_object = FALSE) { + all_warnings <- list() + rv <- withCallingHandlers(x, warning = function (w) { + all_warnings <<- c(all_warnings, list(w)) + invokeRestart("muffleWarning") + }) + if (!full_object) all_warnings <- vapply(all_warnings, function (w) w$message, character(1)) + return(list(rv = rv, warnings = all_warnings)) + } > > ok(ut_cmp_error({ + invalid_subset <- array(dim = c(2,2,2)) + g3_to_tmb(list(~{invalid_subset[,g3_idx(1),]})) + }, "invalid_subset"), "Complained when trying to subset in middle") ok - Complained when trying to subset in middle > ok(ut_cmp_error({ + invalid_subset <- array(dim = c(2,2,2)) + g3_to_tmb(list(~{invalid_subset[g3_idx(1),,] <- 0})) + }, "invalid_subset"), "Complained when trying to subset by row in an lvalue") ok - Complained when trying to subset by row in an lvalue > > ok(ut_cmp_error({ + g3_to_tmb(list(~{ 2 + NA })) + }, "NA"), "No such thing as NA in C++") ok - No such thing as NA in C++ > > ok(grepl( + "unknown_func(.*2.*)", + paste(g3_to_tmb(list(~{ + unknown_func(2) + })), collapse = "\n"), + perl = TRUE), "Unknown functions that are valid C++ get included") ok - Unknown functions that are valid C++ get included > > ok(ut_cmp_error({ + g3_to_tmb(list(~`0unknown0`(2))) + }, "0unknown0"), "An unknown function has to at least be a valid C++ function") ok - An unknown function has to at least be a valid C++ function > > ok(ut_cmp_error({ + g3_to_tmb(list(~not.a.function(2))) + }, "not\\.a\\.function"), "An unknown function has to at least be a valid C++ function") ok - An unknown function has to at least be a valid C++ function > > ok(ut_cmp_error({ + g3_tmb_adfun(g3_to_tmb(list(~{ + g3_param('randomandoptimise', random = TRUE, optimise = TRUE) + g3_param('ro2', random = TRUE, optimise = TRUE) + }))) + }, "randomandoptimise,ro2"), "Specifying random and optimise isn't allowed") ok - Specifying random and optimise isn't allowed > > ok(ut_cmp_error({ + g3_tmb_adfun(g3_to_tmb(list(~{ g3_param_table('randomandoptimise', expand.grid(cur_step = 2:3), random = TRUE, optimise = TRUE) }))) + }, "randomandoptimise"), "Specifying random and optimise isn't allowed") ok - Specifying random and optimise isn't allowed > > ok(ut_cmp_error({ + g3_to_tmb(list(~g3_param("camel", optimize = FALSE))) + }, "optimise"), "Optimise is spelt with an s in g3_param()") ok - Optimise is spelt with an s in g3_param() > > ok_group("Exponentiate params") # Exponentiate params > params.in <- attr(g3_to_tmb(list( g3a_time(1990, 2000), g3_formula( + quote(d), + d = g3_parameterized('par.years', value = 0, by_year = TRUE, exponentiate = TRUE), + x = NA) )), 'parameter_template') > ok(ut_cmp_identical(params.in[grep('^par', params.in$switch), 'switch'], c( + paste0('par.years.', 1990:2000, '_exp'), + NULL)), "exponentiate prefix ends up at the end of parameters") ok - exponentiate prefix ends up at the end of parameters > > ok_group('g3_tmb_par', { + param <- attr(g3_to_tmb(list(~{ + g3_param('param.b') + g3_param_vector('param_vec') + g3_param('aaparam') + g3_param('randomparam', random = TRUE) + })), 'parameter_template') + param$value <- I(list( + aaparam = 55, + param.b = 66, + randomparam = 2, + param_vec = 6:10)[rownames(param)]) + + ok(ut_cmp_identical(g3_tmb_par(param, include_random = FALSE), c( + param__b = 66, + param_vec1 = 6, param_vec2 = 7, param_vec3 = 8, param_vec4 = 9, param_vec5 = 10, + aaparam = 55)), "g3_tmb_par: Flattened parameters in right order") + + param['param_vec', 'optimise'] <- FALSE + ok(ut_cmp_identical(g3_tmb_par(param, include_random = FALSE), c( + param__b = 66, + aaparam = 55)), "g3_tmb_par: Turning off optimise removed values") + + ok(ut_cmp_identical(g3_tmb_par(param, include_random = TRUE), c( + param__b = 66, + aaparam = 55, + randomparam = 2)), "g3_tmb_par: randomparam visible if include_random on") + }) # g3_tmb_par ok - g3_tmb_par: Flattened parameters in right order ok - g3_tmb_par: Turning off optimise removed values ok - g3_tmb_par: randomparam visible if include_random on > > ok_group('g3_tmb_lower', { + param <- attr(g3_to_tmb(list(~{ + g3_param('param.b') + g3_param_vector('param_vec') + g3_param('aaparam') + # NB: Never visible + g3_param('randomparam', random = TRUE) + })), 'parameter_template') + param$value <- I(list( + aaparam = 55, + param.b = 66, + randomparam = 2, + param_vec = 6:10)[rownames(param)]) + param$lower <- c( + aaparam = 500, + param.b = 600, + param_vec = 100)[rownames(param)] + + ok(ut_cmp_identical(g3_tmb_lower(param), c( + param__b = 600, + param_vec1 = 100, param_vec2 = 100, param_vec3 = 100, param_vec4 = 100, param_vec5 = 100, + aaparam = 500)), "g3_tmb_lower: All lower bounds in right order") + + param['param_vec', 'lower'] <- 3 + ok(ut_cmp_identical(g3_tmb_lower(param), c( + param__b = 600, + param_vec1 = 3, param_vec2 = 3, param_vec3 = 3, param_vec4 = 3, param_vec5 = 3, + aaparam = 500)), "g3_tmb_lower: Set all lower values of param_vec in one go") + ok(ut_cmp_identical( + names(g3_tmb_par(param, include_random = FALSE)), + names(g3_tmb_lower(param))), "g3_tmb_lower: Structure matches par after setting param_vec") + + param['param.b', 'optimise'] <- FALSE + ok(ut_cmp_identical(g3_tmb_lower(param), c( + param_vec1 = 3, param_vec2 = 3, param_vec3 = 3, param_vec4 = 3, param_vec5 = 3, + aaparam = 500)), "g3_tmb_lower: Cleared param.b by setting optimise = F") + ok(ut_cmp_identical( + names(g3_tmb_par(param, include_random = FALSE)), + names(g3_tmb_lower(param))), "g3_tmb_lower: Structure matches par after setting param.b") + }) # g3_tmb_lower ok - g3_tmb_lower: All lower bounds in right order ok - g3_tmb_lower: Set all lower values of param_vec in one go ok - g3_tmb_lower: Structure matches par after setting param_vec ok - g3_tmb_lower: Cleared param.b by setting optimise = F ok - g3_tmb_lower: Structure matches par after setting param.b > > ok_group('g3_tmb_upper', { + param <- attr(g3_to_tmb(list(~{ + g3_param('param.b') + g3_param_vector('param_vec') + g3_param('aaparam') + # NB: Never visible + g3_param('randomparam', random = TRUE) + })), 'parameter_template') + param$value <- I(list( + aaparam = 55, + param.b = 66, + randomparam = 2, + param_vec = 6:10)[rownames(param)]) + param$upper <- c( + aaparam = 500, + param.b = 600, + param_vec = 100)[rownames(param)] + + ok(ut_cmp_identical(g3_tmb_upper(param), c( + param__b = 600, + param_vec1 = 100, param_vec2 = 100, param_vec3 = 100, param_vec4 = 100, param_vec5 = 100, + aaparam = 500)), "g3_tmb_upper: All upper bounds in right order") + + param['param_vec', 'upper'] <- 3 + ok(ut_cmp_identical(g3_tmb_upper(param), c( + param__b = 600, + param_vec1 = 3, param_vec2 = 3, param_vec3 = 3, param_vec4 = 3, param_vec5 = 3, + aaparam = 500)), "g3_tmb_upper: Set all lower values of param_vec in one go") + ok(ut_cmp_identical( + names(g3_tmb_par(param, include_random = FALSE)), + names(g3_tmb_lower(param))), "g3_tmb_upper: Structure matches par after setting param_vec") + + param['param.b', 'optimise'] <- FALSE + ok(ut_cmp_identical(g3_tmb_upper(param), c( + param_vec1 = 3, param_vec2 = 3, param_vec3 = 3, param_vec4 = 3, param_vec5 = 3, + aaparam = 500)), "g3_tmb_upper: Cleared param.b by setting optimise = F") + ok(ut_cmp_identical( + names(g3_tmb_par(param, include_random = FALSE)), + names(g3_tmb_lower(param))), "g3_tmb_lower: Structure matches par after setting param.b") + }) # g3_tmb_upper ok - g3_tmb_upper: All upper bounds in right order ok - g3_tmb_upper: Set all lower values of param_vec in one go ok - g3_tmb_upper: Structure matches par after setting param_vec ok - g3_tmb_upper: Cleared param.b by setting optimise = F ok - g3_tmb_lower: Structure matches par after setting param.b > > ok_group('g3_tmb_parscale', { + param <- attr(g3_to_tmb(list(~{ + g3_param('param.lu', lower = 4, upper = 8) + g3_param('param.ps', parscale = 22) + g3_param('param.lups', lower = 4, upper = 8, parscale = 44) + })), 'parameter_template') + ok(ut_cmp_identical(g3_tmb_parscale(param), c( + param__lu = NA, + param__ps = 22, + param__lups = 44)), "We populate parscale") + }) # g3_tmb_parscale ok - We populate parscale > > ok_group('g3_tmb_relist', { + param <- attr(g3_to_tmb(list(~{ + g3_param('param.b') + g3_param_vector('param_vec') + g3_param('unopt_param', optimise = FALSE) + g3_param('randomparam', random = TRUE) + g3_param('aaparam') + })), 'parameter_template') + param$value <- I(list( + aaparam = 55, + param.b = 66, + unopt_param = 95, + randomparam = 2, + param_vec = 6:10)[rownames(param)]) + + ok(ut_cmp_identical( + g3_tmb_relist(param, c( + param__b = 660, + param_vec1 = 60, param_vec2 = 70, param_vec3 = 80, param_vec4 = 90, param_vec5 = 100, + aaparam = 550)), + list( + "param.b" = 660, + "param_vec" = c(60, 70, 80, 90, 100), + "unopt_param" = 95, + "randomparam" = 2, + "aaparam" = 550)), "g3_tmb_relist: Put parameters back in right slots, used old unopt_param value") + ok(ut_cmp_identical( + g3_tmb_relist(param, c( + param__b = 660, + param_vec1 = 60, param_vec2 = 70, param_vec3 = 80, param_vec4 = 90, param_vec5 = 100, + randomparam = 5, # NB: randomparam included, so we update it + aaparam = 550)), + list( + "param.b" = 660, + "param_vec" = c(60, 70, 80, 90, 100), + "unopt_param" = 95, + "randomparam" = 5, + "aaparam" = 550)), "g3_tmb_relist: Put parameters back in right slots, including random") + + param['param.b', 'optimise'] <- FALSE + ok(ut_cmp_error( + g3_tmb_relist(param, c( + param__b = 660, + param_vec1 = 60, param_vec2 = 70, param_vec3 = 80, param_vec4 = 90, param_vec5 = 100, + randomparam = 6, # NB: randomparam included, so we update it + aaparam = 550)), + "par"), "g3_tmb_relist: Still included param__b in par, now an error") + ok(ut_cmp_identical( + g3_tmb_relist(param, c( + param_vec1 = 60, param_vec2 = 70, param_vec3 = 80, param_vec4 = 90, param_vec5 = 100, + aaparam = 550)), + list( + "param.b" = 66, + "param_vec" = c(60, 70, 80, 90, 100), + "unopt_param" = 95, + "randomparam" = 2, + "aaparam" = 550)), "g3_tmb_relist: Works without param.b set, use initial table value") + }) # g3_tmb_relist ok - g3_tmb_relist: Put parameters back in right slots, used old unopt_param value ok - g3_tmb_relist: Put parameters back in right slots, including random ok - g3_tmb_relist: Still included param__b in par, now an error ok - g3_tmb_relist: Works without param.b set, use initial table value > > ok_group('g3_param', { + param <- attr(g3_to_tmb(list(g3a_time(2000, 2004), ~{ + g3_param('a') + g3_param('b', value = 4, optimise = FALSE, random = TRUE, lower = 5, upper = 10) + })), 'parameter_template') + ok(ut_cmp_identical( + param[c('a', 'b'),], + data.frame( + row.names = c('a', 'b'), + switch = c('a', 'b'), + type = c("", ""), + value = I(list(a = 0, b = 4)), + optimise = c(TRUE, FALSE), + random = c(FALSE, TRUE), + lower = c(NA, 5), + upper = c(NA, 10), + parscale = as.numeric(c(NA, NA)), # NB: Not derived yet, only when calling g3_tmb_parscale() + stringsAsFactors = FALSE)), "Param table included custom values") + }) # g3_param ok - Param table included custom values > > ok_group('g3_param_table', { + param <- attr(g3_to_tmb(list(g3a_time(2000, 2004, step_lengths = c(3,3,3,3)), ~{ + g3_param_table('pt', expand.grid( # NB: We can use base R + cur_year = seq(start_year, end_year), # NB: We can use g3a_time's vars + cur_step = 2:3)) + g3_param_table('pg', expand.grid( + cur_year = start_year, + cur_step = 1:2), value = 4, optimise = FALSE, random = TRUE, lower = 5, upper = 10) + })), 'parameter_template') + ok(ut_cmp_identical( + param[c(paste('pt', 2000:2004, 2, sep = '.'), paste('pt', 2000:2004, 3, sep = '.'), 'pg.2000.1', 'pg.2000.2'),], + data.frame( + row.names = c(paste('pt', 2000:2004, 2, sep = '.'), paste('pt', 2000:2004, 3, sep = '.'), 'pg.2000.1', 'pg.2000.2'), + switch = c(paste('pt', 2000:2004, 2, sep = '.'), paste('pt', 2000:2004, 3, sep = '.'), 'pg.2000.1', 'pg.2000.2'), + type = c("", "", "", "", "", "", "", "", "", "", "", ""), + value = I(list( + "pt.2000.2" = 0, + "pt.2001.2" = 0, + "pt.2002.2" = 0, + "pt.2003.2" = 0, + "pt.2004.2" = 0, + "pt.2000.3" = 0, + "pt.2001.3" = 0, + "pt.2002.3" = 0, + "pt.2003.3" = 0, + "pt.2004.3" = 0, + "pg.2000.1" = 4, + "pg.2000.2" = 4)), + optimise = c(TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, FALSE, FALSE), + random = c(FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, TRUE, TRUE), + lower = c(NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 5, 5), + upper = c(NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 10, 10), + # NB: Not derived yet, only when calling g3_tmb_parscale() + parscale = as.numeric(c(NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA)), + stringsAsFactors = FALSE)), "Param table included custom values") + ok(ut_cmp_identical( + attr(g3_to_tmb(list(~{ + g3_param_table('moo', expand.grid(cur_year=1990:1994)) + g3_param('moo.1990') + g3_param('oink.2') + })), 'parameter_template')$switch, + c("moo.1990", "moo.1991", "moo.1992", "moo.1993", "moo.1994", "oink.2")), "Refering to individual parameters doesn't generate duplicate table entries") + }) # g3_param_table ok - Param table included custom values ok - Refering to individual parameters doesn't generate duplicate table entries > > ok_group("g3_to_tmb: attr.actions", { + actions <- list( + list("001" = ~{ 1 + 1 }, "002" = ~{2 + 2}), + "003" = ~{3 + 3}) + model_fn <- g3_to_tmb(actions) + ok(ut_cmp_identical(attr(model_fn, 'actions'), actions), "actions returned as attribute uncollated") + }) # g3_to_tmb: attr.actions ok - actions returned as attribute uncollated > > ok_group("g3_to_tmb: Can use random parameters without resorting to include_random", local({ + # ./TMB/inst/examples/randomregression.R converted into gadget3 + actions <- local({ + set.seed(123) + n <- 1 + ng <- 1 + f <- gl(ng, n/ng) + t <- rnorm(n, mean=2, sd=5) + a <- rnorm(ng, mean=3, sd=4) + b <- rnorm(ng, mean=8, sd=7) + x <- a[f] * t + b[f] + rnorm(n, mean=0, sd=1) + set.seed(NULL) + + list( + g3a_time(1990, 1991), + g3l_random_dnorm("a", + ~g3_param('a', value=1, random=TRUE), + ~g3_param('mu.a', value=1), + ~g3_param('sigma.a', value=1), period="single"), + g3l_random_dnorm("b", + ~g3_param('b', value=1, random=TRUE), + ~g3_param('mu.b', value=1), + ~g3_param('sigma.b', value=1), period="single"), + g3l_random_dnorm("x", + ~x, + ~g3_param('a', value=1, random=TRUE) * t + g3_param('b', value=1, random=TRUE), + ~g3_param('sigma0', value=1), period="single" )) + }) + model_fn <- g3_to_r(actions) + + if (nzchar(Sys.getenv('G3_TEST_TMB'))) { + model_cpp <- g3_to_tmb(actions) + model_tmb <- g3_tmb_adfun(model_cpp, compile_flags = c("-O0", "-g"), inner.control = list(fnscale = -1)) + param_tbl <- attr(model_cpp, 'parameter_template') + param_tbl[,'lower'] <- -2 + param_tbl[,'upper'] <- 3 + param_tbl[,'parscale'] <- 1 + + ok(ut_cmp_identical(names(model_tmb$env$last.par)[model_tmb$env$random], c( + "a", + "b", + NULL)), "env$random: TMB got the random parameters we expected") + + ok(ut_cmp_error({ + nlminb(g3_tmb_par(param_tbl), model_tmb$fn, model_tmb$gr) + }, "\\$par"), "g3_tmb_par: Not allowed in nlminb() call") + ok(ut_cmp_error({ + optim(g3_tmb_par(param_tbl), model_tmb$fn, model_tmb$gr) + }, "\\$par"), "g3_tmb_par: Not allowed in optim() call") + + res <- suppressWarnings({ + nlminb(model_tmb$par, model_tmb$fn, model_tmb$gr, + upper = g3_tmb_upper(param_tbl), + lower = g3_tmb_lower(param_tbl), + control = list( + fnscale = -1, + parscale = g3_tmb_parscale(param_tbl))) + }) + ok(res$convergence == 0, "Model ran successfully and converged") + + # The first won't have random parameters, the latter will. Both should work + value_no_random <- g3_tmb_relist(param_tbl, res$par) + value_inc_random <- g3_tmb_relist(param_tbl, model_tmb$env$last.par) + ok(ut_cmp_equal( + value_no_random[param_tbl[!param_tbl$random, 'switch']], + value_inc_random[param_tbl[!param_tbl$random, 'switch']]), "Non-random parameters match") + ok(ut_cmp_equal( + I(value_no_random[param_tbl[param_tbl$random, 'switch']]), + param_tbl$value[param_tbl[param_tbl$random, 'switch']]), "value_no_random: Random parameters as default") + ok(!isTRUE(all.equal( + I(value_inc_random[param_tbl[param_tbl$random, 'switch']]), + param_tbl$value[param_tbl[param_tbl$random, 'switch']])), "value_random: Random parameters have been updated") + param_tbl$value <- value_inc_random + + # NB: ut_tmb_r_compare will be using g3_tmb_par() + gadget3:::ut_tmb_r_compare(model_fn, model_tmb, param_tbl) + } else { + writeLines("# skip: not compiling TMB model") + } + })) # g3_to_tmb: Can use random parameters without resorting to include_random # skip: not compiling TMB model > > ok_group("cpp_code", { + ns <- function (s) gsub("\\s+", "", s) + ok(ut_cmp_identical( + ns(gadget3:::cpp_code(quote( 4 + 4 ))), + ns("(double)(4) + (double)(4)")), + "4 + 4: Defaulted to double") + ok(ut_cmp_identical( + ns(gadget3:::cpp_code(quote( (4 + 4) ))), + ns("((double)(4) + (double)(4))")), + "(4 + 4): Defaulted to double") + ok(ut_cmp_identical( + ns(gadget3:::cpp_code(quote( x[[4]] + 4 ))), + ns("x(3) + (double)(4)")), + "x[[4]] + 4: Subtracted one, assumed int") + # NB: We rely on this in g3_param_table(ifmissing_param_table) below + ok(ut_cmp_identical( + ns(gadget3:::cpp_code(quote( x[[(-4 + 1 - 1) / (2*8)]] + 4 ))), + ns("x( (-4 + 1 - 1) / (2*8) - 1 ) + (double)(4)")), + "x[[(-4 + 1 - 1) / (2*8)]] + 4: Assumed int passses through +, -, *, /, ( ") + }) # cpp_code ok - 4 + 4: Defaulted to double ok - (4 + 4): Defaulted to double ok - x[[4]] + 4: Subtracted one, assumed int ok - x[[(-4 + 1 - 1) / (2*8)]] + 4: Assumed int passses through +, -, *, /, ( > > > ############################################################################### > actions <- list() > expecteds <- new.env(parent = emptyenv()) > expected_warnings_r <- c() > expected_warnings_tmb <- c() > params <- list(rv=0) > > # Check constants can pass through cleanly > constants_integer <- 999 > constants_nan <- 999 > actions <- c(actions, ~{ + comment('constants') + constants_integer <- 5L ; REPORT(constants_integer) + constants_nan <- NaN ; REPORT(constants_nan) + }) > expecteds$constants_integer <- 5L > expecteds$constants_nan <- NaN > > # Can assign a single value to 1x1 array > assign_to_1_1 <- array(dim = c(1,1)) > actions <- c(actions, ~{ + comment('assign_to_1_1') + assign_to_1_1[g3_idx(1)] <- 100.1 + REPORT(assign_to_1_1) + }) > expecteds$assign_to_1_1 <- array(c(100.1), dim = c(1,1)) > > assign_to_2_1 <- array(dim = c(2,1)) > data_to_2_1 <- c(100.1, 200.2) > actions <- c(actions, ~{ + comment('assign_to_2_1') + assign_to_2_1[,g3_idx(1)] <- data_to_2_1 + REPORT(assign_to_2_1) + }) > expecteds$assign_to_2_1 <- array(c(100.1, 200.2), dim = c(2,1)) > > assign_to_2_2 <- array(dim = c(2,2)) > data_to_2_2 <- c(110.1, 220.2) > actions <- c(actions, ~{ + comment('assign_to_2_2') + assign_to_2_2[,g3_idx(1)] <- data_to_2_1 + assign_to_2_2[,g3_idx(2)] <- data_to_2_2 + REPORT(assign_to_2_2) + }) > expecteds$assign_to_2_2 <- array(c(100.1, 200.2, 110.1, 220.2), dim = c(2,2)) > > # Assign single value, horizontally > assign_to_2_2a <- array(dim = c(2,2)) > actions <- c(actions, ~{ + comment('assign_to_2_2a') + assign_to_2_2a[[g3_idx(1),g3_idx(2)]] <- 99 + assign_to_2_2a[[g3_idx(2),g3_idx(1)]] <- 88 + assign_to_2_2a[[g3_idx(2),g3_idx(2)]] <- 0 + assign_to_2_2a[[g3_idx(1),g3_idx(1)]] <- 0 + REPORT(assign_to_2_2a) + }) > expecteds$assign_to_2_2a <- array(c(0, 88, 99, 0), dim = c(2,2)) > > # Assign zero and other scalars to column, arrays > assign_scalar <- array(dim = c(2,4)) > params$assign_scalar_const <- runif(1, 100, 200) > actions <- c(actions, ~{ + comment('assign_scalar') + assign_scalar[] <- 0 # TODO: TMB auto-setZeros, R doesn't + assign_scalar[g3_idx(1),g3_idx(1)] <- 99 # NB: Overwritten + assign_scalar[,g3_idx(1)] <- 0 + assign_scalar[,g3_idx(2)] <- 88 + assign_scalar[g3_idx(2),g3_idx(3)] <- 27 + assign_scalar[,g3_idx(4)] <- g3_param("assign_scalar_const") + REPORT(assign_scalar) + }) > expecteds$assign_scalar <- array(c( + 0, 0, + 88, 88, + 0, 27, + params$assign_scalar_const, params$assign_scalar_const, + NULL), dim = c(2,4)) > > # Can subset from left and right, using transpose > assign_subset <- array( + rep(seq(100, 200, 100), each = 3 * 5) + + rep(seq(10, 30, 10), each = 5, times = 2) + + rep(seq(1, 5, 1), each = 1, times = 2), + dim = c(5, 3, 2)) > assign_right <- assign_subset[,2,2] ; assign_right[] <- 0 > assign_leftright <- assign_subset[4,,2] ; assign_leftright[] <- 0 > actions <- c(actions, ~{ + comment('assign_subset') + assign_right <- assign_subset[,g3_idx(2),g3_idx(2)] + REPORT(assign_right) + assign_leftright <- assign_subset[g3_idx(4),,g3_idx(2)] + REPORT(assign_leftright) + }) > expecteds$assign_right <- assign_subset[,2,2] > expecteds$assign_leftright <- assign_subset[4,,2] > > # Arrays with dynamic dimensions > dynamic_dim_array <- array(0, dim = c(2,1)) > dynamic_dim_array_na <- array(NA, dim = c(2,1)) > attr(dynamic_dim_array, 'dynamic_dim') <- list(2, quote(dynamic_dim_array_dim_2)) > attr(dynamic_dim_array_na, 'dynamic_dim') <- list(2, quote(dynamic_dim_array_dim_2)) > dynamic_dim_array_dim_2 <- 4L > dynamic_dim_array_dim_na_2 <- 4L > actions <- c(actions, ~{ + comment('assign_scalar') + # NB: Under TMB this won't be NA, so set it to hide differences. + dynamic_dim_array_na[] <- 5 + REPORT(dynamic_dim_array) + REPORT(dynamic_dim_array_na) + }) > expecteds$dynamic_dim_array <- array(0, dim = c(2, 4)) > expecteds$dynamic_dim_array_na <- array(5, dim = c(2, 4)) > > # Data variable names are escaped > escaped.data.scalar <- 44 > escaped.data.array <- c(1,2,3,4,5) > escaped_data_output <- 0.0 > actions <- c(actions, ~{ + comment('escaped.data.array') + escaped_data_output <- escaped.data.scalar + sum(escaped.data.array) + REPORT(escaped_data_output) + }) > # NB: In theory we could also report the escaped name, but we can't rename > # items in the produced report, so reports will differ TMB/R > expecteds$escaped_data_output <- escaped.data.scalar + sum(escaped.data.array) > > # is.nan / is.finite > is_nan_nan_scalar_input <- NaN > is_nan_finite_scalar_input <- 4 > is_nan_finite_array_input <- as.array(1:4) > is_nan_nan_array_input <- as.array(rep(NaN, 5)) > is_nan_finite_vector_input <- 1:4 > is_nan_nan_vector_input <- rep(NaN, 5) > is_nan_finite_single_array_input <- as.array(1) > is_nan_output <- 0L > is_finite_output <- 0L > actions <- c(actions, ~{ + comment('is.nan / is.finite') + is_nan_output <- is_nan_output + (if (is.nan(is_nan_nan_scalar_input)) 1 else 0) + is_nan_output <- is_nan_output + (if (is.nan(is_nan_finite_scalar_input)) 2 else 0) + is_nan_output <- is_nan_output + (if (any(is.nan(is_nan_finite_array_input))) 4 else 0) + is_nan_output <- is_nan_output + (if (any(is.nan(is_nan_nan_array_input))) 8 else 0) + is_nan_output <- is_nan_output + (if (any(is.nan(is_nan_nan_vector_input))) 16 else 0) + is_nan_output <- is_nan_output + (if (any(is.nan(is_nan_finite_single_array_input))) 32 else 0) + REPORT(is_nan_output) + is_finite_output <- is_finite_output + (if (is.finite(is_nan_nan_scalar_input)) 1 else 0) + is_finite_output <- is_finite_output + (if (is.finite(is_nan_finite_scalar_input)) 2 else 0) + is_finite_output <- is_finite_output + (if (any(is.finite(is_nan_finite_array_input))) 4 else 0) + is_finite_output <- is_finite_output + (if (any(is.finite(is_nan_nan_array_input))) 8 else 0) + is_finite_output <- is_finite_output + (if (any(is.finite(is_nan_nan_vector_input))) 16 else 0) + is_finite_output <- is_finite_output + (if (any(is.finite(is_nan_finite_single_array_input))) 32 else 0) + REPORT(is_finite_output) + }) > expecteds$is_nan_output <- 1 + 8 + 16 > expecteds$is_finite_output <- 2 + 4 + 32 > > # as.vector() --> .vec() > as_vector_array <- array(runif(20), dim=c(10, 2)) > as_vector_result1 <- rep(NaN, 10) > as_vector_result2 <- rep(NaN, 10) > as_vector_mean <- 0.5 > as_vector_sigma <- 1 > actions <- c(actions, ~{ + comment('as_vector') + as_vector_result1 <- pnorm(as.vector(as_vector_array[,g3_idx(1)]), as_vector_mean, as_vector_sigma) + as_vector_result2 <- pnorm(as.vector(as_vector_array[,g3_idx(2)]), as_vector_mean, as_vector_sigma) + REPORT(as_vector_result1) + REPORT(as_vector_result2) + }) > expecteds$as_vector_result1 <- pnorm(as_vector_array[,1], as_vector_mean, as_vector_sigma) > expecteds$as_vector_result2 <- pnorm(as_vector_array[,2], as_vector_mean, as_vector_sigma) > > # mean() --> .mean() > mean_vector <- array(c(1, 2, 88, 99)) > mean_vector_result <- 0 > actions <- c(actions, ~{ + comment('mean_vector') + mean_vector_result <- mean(mean_vector) + REPORT(mean_vector_result) + }) > expecteds$mean_vector_result <- mean(mean_vector) > > # colsums > colsums_in <- array(1:6, dim = c(3,2)) > colsums_result <- c(0L, 0L) > actions <- c(actions, ~{ + comment('colsums') + colsums_result <- colSums(colsums_in) + REPORT(colsums_result) + }) > expecteds$colsums_result <- colSums(colsums_in) > > # rowsums > rowsums_in <- array(c(1,2,3,4,5,6), dim = c(3,2)) > rowsums_result <- c(0, 0, 0) > actions <- c(actions, ~{ + comment('rowsums') + rowsums_result <- rowSums(rowsums_in) + REPORT(rowsums_result) + }) > expecteds$rowsums_result <- rowSums(rowsums_in) > > # if statement without braces > if_no_brace_result <- 0.0 > actions <- c(actions, ~{ + comment('if_without_braces') + if (FALSE) if_no_brace_result <- 1 else if_no_brace_result <- 0.2 + REPORT(if_no_brace_result) + }) > expecteds$if_no_brace_result <- 0.2 > > # if expression (not statement) turns into tertiary > tertiary_result_0 <- 3 > tertiary_result_1 <- 6 > tertiary_result_2 <- 6 > ok(ut_cmp_error({ + g3_to_tmb(list(~{ tertiary_result_0 <- if (tertiary_result_0 == 3) 9 })) + }, "tertiary_result_0"), "Complained about missing else for tertiary") ok - Complained about missing else for tertiary > actions <- c(actions, ~{ + comment('tertiary') + tertiary_result_0 <- if (tertiary_result_0 == 3) 9 else 4 + tertiary_result_1 <- if (tertiary_result_1 == 3) 9 else 4 + tertiary_result_2 <- 100 - if (tertiary_result_0 == -30) 20 else 40 + REPORT(tertiary_result_0) + REPORT(tertiary_result_1) + REPORT(tertiary_result_2) + }) > expecteds$tertiary_result_0 <- 9 > expecteds$tertiary_result_1 <- 4 > expecteds$tertiary_result_2 <- 100 - 40 # NB: Operator precedence preserved > > # g3_with() > g3_with_result <- 0L > # NB: We don't define g3_with_iterator, it's defined within the block > actions <- c(actions, ~{ + comment('g3_with') + g3_with( + g3_with_iterator := g3_idx(2L), # NB: Tests we can use g3 functions in definition + g3_with_other_exp := 1L + 2L + 3L, + { + g3_with_result <- g3_with_iterator - g3_idx(1) # NB: Reverse g3_idx from definition + REPORT(g3_with_result) + REPORT(g3_with_other_exp) + }) + }) > expecteds$g3_with_result <- 1L # i.e. 2 - 1 in R or 1 - 0 in TMB > expecteds$g3_with_other_exp <- 1L + 2L + 3L > > # min() & max() > min_result <- 0.0 > max_result <- 0.0 > actions <- c(actions, ~{ + comment('min/max') + min_result <- min(4, 9) + max_result <- max(sum(mean_vector), 2) # NB: sum gets cast to Type + REPORT(min_result) + REPORT(max_result) + }) > expecteds$min_result <- 4 > expecteds$max_result <- sum(mean_vector) > > # negate single value > negate_x <- 10 > actions <- c(actions, ~{ + comment('negate') + negate_x <- -negate_x + REPORT(negate_x) + }) > expecteds$negate_x <- -10 > > # Modulo assumes integer > modulo_x <- 10 > modulo_y <- 11 > actions <- c(actions, ~{ + comment('modulo') + modulo_x <- as_integer(modulo_x) %% 2 + modulo_y <- as_integer(modulo_y) %% 2 + REPORT(modulo_x) + REPORT(modulo_y) + }) > expecteds$modulo_x <- 0 > expecteds$modulo_y <- 1 > > # sum() & prod() > sumprod_input <- runif(10) > sumprod_sum <- 0 > sumprod_prod <- 0 > actions <- c(actions, ~{ + comment('sum/prod') + sumprod_sum <- sum(sumprod_input) + sumprod_prod <- prod(sumprod_input) + REPORT(sumprod_sum) + REPORT(sumprod_prod) + }) > expecteds$sumprod_sum <- sum(sumprod_input) > expecteds$sumprod_prod <- prod(sumprod_input) > > # g3_param_table() > param_table_out <- array(rep(0, 6)) > actions <- c(actions, ~{ + for (pt_a in seq(1, 3, by = 1)) for (pt_b in seq(7, 8, by = 1)) { + param_table_out[[((pt_a - 1L) * 2L) + (pt_b - 7L) + 1L]] <- g3_param_table('param_table', expand.grid(pt_a = 1:2, pt_b = c(8, 7))) + } + REPORT(param_table_out) + }) > params[['param_table.1.8']] <- 18 > params[['param_table.1.7']] <- 17 > params[['param_table.2.8']] <- 28 > params[['param_table.2.7']] <- 27 > expecteds$param_table_out <- array(c( + 17, + 18, + 27, + 28, + NaN, + NaN, + NULL)) > expected_warnings_r <- c(expected_warnings_r, + "No value found in g3_param_table param_table, ifmissing not specified", + "No value found in g3_param_table param_table, ifmissing not specified", + NULL) > expected_warnings_tmb <- c(expected_warnings_tmb, + "No value found in g3_param_table param_table, ifmissing not specified", + "No value found in g3_param_table param_table, ifmissing not specified", + NULL) > > # g3_param_table(ifmissing) > param_table_ifmissing_out <- array(c(1,2,3,4,5,6)) > actions <- c(actions, ~{ + for (ifmissing in seq(1, 6, by = 1)) { + param_table_ifmissing_out[[ifmissing]] <- g3_param_table( + 'param_table_ifmissing', + expand.grid(ifmissing = 2:4), ifmissing = -1) + } + REPORT(param_table_ifmissing_out) + }) > params[['param_table_ifmissing.2']] <- 27 > params[['param_table_ifmissing.3']] <- 47 > params[['param_table_ifmissing.4']] <- 22 > expecteds$param_table_ifmissing_out <- array(c(-1, 27, 47, 22, -1, -1)) > > # g3_param_table(ifmissing_param) > param_table_ifmparam_out <- array(c(1,2,3,4,5,6)) > actions <- c(actions, ~{ + for (ifmparam in seq(1, 6, by = 1)) { + param_table_ifmparam_out[[ifmparam]] <- g3_param_table( + 'param_table_ifmparam', + expand.grid(ifmparam = 2:4), ifmissing = g3_param("param_table_ifmparam.missing")) + } + REPORT(param_table_ifmparam_out) + }) > params[['param_table_ifmparam.2']] <- floor(runif(1, 100, 200)) > params[['param_table_ifmparam.3']] <- floor(runif(1, 100, 200)) > params[['param_table_ifmparam.4']] <- floor(runif(1, 100, 200)) > params[['param_table_ifmparam.missing']] <- floor(runif(1, 100, 200)) > expecteds$param_table_ifmparam_out <- array(c( + params[['param_table_ifmparam.missing']], + params[['param_table_ifmparam.2']], + params[['param_table_ifmparam.3']], + params[['param_table_ifmparam.4']], + params[['param_table_ifmparam.missing']], + params[['param_table_ifmparam.missing']])) > > # g3_param_table(ifmissing_param_table) > param_table_ifmpartab_out <- array(c(1,2,3,4,5,6,7,8,9)) > actions <- c(actions, ~{ + for (ifmpartab in seq(1, 3, by = 1)) for (ifmpartab_m in seq(1, 3, by = 1)) { + param_table_ifmpartab_out[[((ifmpartab - 1) * 3) + (ifmpartab_m - 1) + 1]] <- g3_param_table( + 'param_table_ifmpartab', + expand.grid(ifmpartab = 2:3), ifmissing = g3_param_table( + "param_table_ifmpartab_m", + expand.grid(ifmpartab_m = 1:3))) + } + REPORT(param_table_ifmpartab_out) + }) > params[['param_table_ifmpartab.2']] <- floor(runif(1, 100, 200)) > params[['param_table_ifmpartab.3']] <- floor(runif(1, 100, 200)) > params[['param_table_ifmpartab_m.1']] <- floor(runif(1, 100, 200)) > params[['param_table_ifmpartab_m.2']] <- floor(runif(1, 100, 200)) > params[['param_table_ifmpartab_m.3']] <- floor(runif(1, 100, 200)) > expecteds$param_table_ifmpartab_out <- array(c( + params[['param_table_ifmpartab_m.1']], + params[['param_table_ifmpartab_m.2']], + params[['param_table_ifmpartab_m.3']], + + params[['param_table_ifmpartab.2']], + params[['param_table_ifmpartab.2']], + params[['param_table_ifmpartab.2']], + + params[['param_table_ifmpartab.3']], + params[['param_table_ifmpartab.3']], + params[['param_table_ifmpartab.3']], + NULL)) > > ############################################################################### > > nll <- 0.0 > actions <- c(actions, ~{ + comment('done') + nll <- nll + g3_param('rv') + return(nll) + }) > > # 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") + } # skip: not compiling TMB model > > # Compare everything we've been told to compare > result <- capture_warnings(model_fn(params)) > # str(attributes(result), vec.len = 10000) > for (n in ls(expecteds)) { + ok(ut_cmp_equal( + attr(result$rv, n), + expecteds[[n]], tolerance = 1e-6), n) + } ok - as_vector_result1 ok - as_vector_result2 ok - assign_leftright ok - assign_right ok - assign_scalar ok - assign_to_1_1 ok - assign_to_2_1 ok - assign_to_2_2 ok - assign_to_2_2a ok - colsums_result ok - constants_integer ok - constants_nan ok - dynamic_dim_array ok - dynamic_dim_array_na ok - escaped_data_output ok - g3_with_other_exp ok - g3_with_result ok - if_no_brace_result ok - is_finite_output ok - is_nan_output ok - max_result ok - mean_vector_result ok - min_result ok - modulo_x ok - modulo_y ok - negate_x ok - param_table_ifmissing_out ok - param_table_ifmparam_out ok - param_table_ifmpartab_out ok - param_table_out ok - rowsums_result ok - sumprod_prod ok - sumprod_sum ok - tertiary_result_0 ok - tertiary_result_1 ok - tertiary_result_2 > ok(ut_cmp_identical(result$warnings, expected_warnings_r), "Warnings from R as expected") ok - Warnings from R as expected > if (nzchar(Sys.getenv('G3_TEST_TMB'))) { + param_template <- attr(model_cpp, "parameter_template") + param_template$value <- params[param_template$switch] + generated_warnings <- capture_warnings(gadget3:::ut_tmb_r_compare(model_fn, model_tmb, param_template))$warnings + ok(ut_cmp_identical(generated_warnings, c(expected_warnings_r, expected_warnings_tmb)), "Warnings from R+TMB as expected") + } > > proc.time() user system elapsed 1.21 0.04 1.25 # Looks like you passed all 74 tests.