tests/test-action_age.R

In gadget-framework/gadget3: Globally-Applicable Area Disaggregated General Ecosystem Toolbox V3

library(magrittr)
library(unittest)

library(gadget3)

ok_group("g3a_age:no_age")
ok(ut_cmp_error({
    g3a_age(g3_stock('prey_a', seq(20, 40, 4)))
}, "age.*prey_a"), "Cannot use g3a_age() without an age dimension")
#### g3a_age:no_age

ok_group("g3a_age:single_age", {
    # prey_a, with single age, outputs into b & c
    prey_a <- g3_stock('prey_a', seq(20, 40, 4)) %>% g3s_age(11, 11)
    prey_b <- g3_stock('prey_b', seq(20, 40, 4)) %>% g3s_age(11, 15)
    prey_c <- g3_stock('prey_c', seq(20, 40, 4)) %>% g3s_age(11, 15)

    model_fn <- g3_to_r(c(
        list("99999" = g3_formula({
            # Keep TMB happy
            nll <- nll + g3_param("dummy", value = 0)
            REPORT(prey_a__num)
            REPORT(prey_b__num)
            REPORT(prey_c__num)
        })),
        g3a_initialconditions(prey_a, ~10 * (age-10) + prey_a__midlen * 0, ~100 * (age-10) + prey_a__midlen * 0),
        g3a_initialconditions(prey_b, ~10 * (age-10) + prey_b__midlen * 0, ~100 * (age-10) + prey_b__midlen * 0),
        g3a_initialconditions(prey_c, ~prey_c__midlen * 0, ~prey_c__midlen * 0),
        g3a_age(prey_a, output_stocks = list(prey_b, prey_c), output_ratios = c(0.75, 0.25)),
        g3a_time(2000, 2002)))
    model_cpp <- g3_to_tmb(attr(model_fn, 'actions'), trace = FALSE)
    if (Sys.getenv('G3_TEST_TMB') == "2") {
        #model_cpp <- edit(model_cpp)
        #writeLines(TMB::gdbsource(g3_tmb_adfun(model_cpp, compile_flags = "-g", output_script = TRUE)))
        model_tmb <- g3_tmb_adfun(model_cpp, trace = TRUE, compile_flags = c("-O0", "-g"))
    }

    params <- attr(model_fn, 'parameter_template')
    r <- model_fn(params)
    for (len_idx in 1:6) {
        ok(ut_cmp_equal(
            attr(r, 'prey_a__num')[length=len_idx,],
            c(0)), paste0("prey_a__num length=", len_idx))
        # prey_b got 75% of prey_a
        ok(ut_cmp_equal(
            attr(r, 'prey_b__num')[length=len_idx,],
            c(age11=10, age12=20 + 7.5, age13=30, age14=40, age15=50)), paste0("prey_b__num length=", len_idx))
        # prey_c got 25% of prey_a
        ok(ut_cmp_equal(
            attr(r, 'prey_c__num')[length=len_idx,],
            c(age11=0, age12=2.5, age13=0, age14=0, age15=0)), paste0("prey_c__num length=", len_idx))
    }
    if (Sys.getenv('G3_TEST_TMB') == "2") gadget3:::ut_tmb_r_compare(model_fn, model_tmb, params, model_cpp = model_cpp)
})

# NB: We start at 11 to make sure we age into the right bracket
prey_a <- g3_stock('prey_a', c(1)) %>% g3s_age(11, 15)
prey_b <- g3_stock('prey_b', c(1)) %>% g3s_age(11, 13)
prey_c <- g3_stock('prey_c', c(1)) %>% g3s_age(11, 17)

# Store stock state in temporary variables labelled stock 0..n
report_action <- list()
for (step in 0:3) for (s in list(prey_a, prey_b, prey_c)) {
    assign(paste0("step", step, '_', s$name, '__num'), g3_stock_instance(s))
    assign(paste0("step", step, '_', s$name, '__wgt'), g3_stock_instance(s))
    report_action[[paste0("999:step", step, ':', s$name)]] <- gadget3:::f_substitute(~if (cur_time == step) {
        stepnum[] <- curnum
        stepwgt[] <- curwgt
        REPORT(stepnum)
        REPORT(stepwgt)
    }, list(
        stepnum = as.symbol(paste0("step", step, '_', s$name, '__num')),
        stepwgt = as.symbol(paste0("step", step, '_', s$name, '__wgt')),
        curnum = as.symbol(paste0(s$name, '__num')),
        curwgt = as.symbol(paste0(s$name, '__wgt')),
        step = step))
}

actions <- list(
    g3a_time(2000, 2002, step_lengths = c(6, 6), project_years = 0),
    g3a_initialconditions(prey_a, ~10 * (age-10) + prey_a__midlen * 0, ~100 * (age-10) + prey_a__midlen * 0),
    g3a_initialconditions(prey_b, ~10 * (age-10) + prey_b__midlen * 0, ~100 * (age-10) + prey_b__midlen * 0),
    g3a_initialconditions(prey_c, ~prey_b__midlen * 0, ~prey_b__midlen * 0),
    g3a_age(prey_a),
    g3a_age(prey_b, output_stocks = list(prey_c)),
    g3a_age(prey_c),
    report_action,
    list('999' = ~{ nll <- nll + g3_param('x', value = 1.0) }))

# Compile model
model_fn <- g3_to_r(actions, trace = FALSE)
# model_fn <- edit(model_fn)
if (nzchar(Sys.getenv('G3_TEST_TMB'))) {
    model_cpp <- g3_to_tmb(actions, trace = FALSE)
    # model_cpp <- edit(model_cpp)
    model_tmb <- g3_tmb_adfun(model_cpp, compile_flags = c("-O0", "-g"))
} else {
    writeLines("# skip: not compiling TMB model")
}

ok_group("age", {
    params <- attr(model_fn, 'parameter_template')
    params$x <- 1.0

    result <- model_fn(params)
    r <- attributes(result)
    # str(result)
    # str(as.list(r), vec.len = 10000)

    # Step 0
    ok(ut_cmp_identical(as.vector(r$step0_prey_a__num), c(10,20,30,40,50)), "step0_prey_a__num: Populated by initialconditions, no ageing yet")
    ok(ut_cmp_identical(as.vector(r$step0_prey_a__wgt), c(100,200,300,400,500)), "step0_prey_a__wgt: Populated by initialconditions, no ageing yet")
    ok(ut_cmp_identical(as.vector(r$step0_prey_b__num), c(10,20,30)), "step0_prey_b__num: Populated by initialconditions, no ageing yet")
    ok(ut_cmp_identical(as.vector(r$step0_prey_b__wgt), c(100,200,300)), "step0_prey_b__wgt: Populated by initialconditions, no ageing yet")
    ok(ut_cmp_identical(as.vector(r$step0_prey_c__num), c(0, 0, 0, 0, 0, 0, 0)), "step0_prey_c__num: Starting off empty")
    ok(ut_cmp_identical(as.vector(r$step0_prey_c__wgt), c(0, 0, 0, 0, 0, 0, 0)), "step0_prey_c__wgt: Starting off empty")

    # Step 1
    ok(ut_cmp_identical(as.vector(r$step1_prey_a__num), c(0, 10,20,30,40 + 50)), "step1_prey_a__num: Numbers rotated by 1, final group a plus group")
    ok(ut_cmp_equal(as.vector(r$step1_prey_a__wgt), c(100, 100,200,300, (400 * 40 + 500 * 50) / (40 + 50))), "step1_prey_a__wgt: Numbers rotated by 1, final group averaged plus-group, first group kept __wgt value")
    ok(ut_cmp_identical(as.vector(r$step1_prey_b__num), c(0, 10,20)), "step1_prey_b__num: Numbers rotated by 1, final group transitioned")
    ok(ut_cmp_equal(as.vector(r$step1_prey_b__wgt), c(100, 100, 200)), "step1_prey_b__wgt: Numbers rotated by 1, final group a plus group")
    ok(ut_cmp_identical(as.vector(r$step1_prey_c__num), c(0, 0, 0, 30, 0, 0, 0)), "step1_prey_c__num: Final age of prey_b transferred")
    ok(ut_cmp_equal(as.vector(r$step1_prey_c__wgt), c(0, 0, 0, 300, 0, 0, 0)), "step1_prey_c__wgt: Final weight of prey_b transferred")

    # Step 2
    ok(ut_cmp_identical(as.vector(r$step2_prey_a__num), c(0, 10,20,30,40 + 50)), "step2_prey_a__num: Not final step, nothing changed")
    ok(ut_cmp_identical(as.vector(r$step2_prey_a__wgt), c(100, 100,200,300, (400 * 40 + 500 * 50) / (40 + 50))), "step2_prey_a__wgt: Not final step, nothing changed")
    ok(ut_cmp_identical(as.vector(r$step2_prey_b__num), c(0, 10, 20)), "step2_prey_b__num: Not final step, nothing changed")
    ok(ut_cmp_equal(as.vector(r$step2_prey_b__wgt), c(100, 100, 200)), "step2_prey_b__wgt: Not final step, nothing changed")
    ok(ut_cmp_identical(as.vector(r$step2_prey_c__num), c(0, 0, 0, 30, 0, 0, 0)), "step2_prey_c__num: Not final step, nothing changed")
    ok(ut_cmp_equal(as.vector(r$step2_prey_c__wgt), c(0, 0, 0, 300, 0, 0, 0)), "step2_prey_c__wgt: Not final step, nothing changed")

    # Step 3
    ok(ut_cmp_identical(as.vector(r$step3_prey_a__num), c(0, 0, 10,20,30 + 40 + 50)), "step3_prey_a__num: Numbers rotated by 1, final group a plus group")
    ok(ut_cmp_equal(as.vector(r$step3_prey_a__wgt), c(100, 100, 100,200, (300 * 30 + 400 * 40 + 500 * 50) / (30 + 40 + 50))), "step3_prey_a__wgt: Numbers rotated by 1, final group a averaged plus group")
    ok(ut_cmp_identical(as.vector(r$step3_prey_b__num), c(0, 0, 10)), "step3_prey_b__num: Numbers rotated by 1, final group transitioned")
    ok(ut_cmp_equal(as.vector(r$step3_prey_b__wgt), c(100, 100, 100)), "step3_prey_b__wgt: Numbers rotated by 1, final group a plus group")
    ok(ut_cmp_identical(as.vector(r$step3_prey_c__num), c(0, 0, 0, 20, 30, 0, 0)), "step3_prey_c__num: Final age of prey_b transferred, numbers rotated")
    ok(ut_cmp_equal(as.vector(r$step3_prey_c__wgt), c(0, 0, 0, 200, 300, 0, 0)), "step3_prey_c__wgt: Final age of prey_b transferred")

    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)
    }
})