inttest/codegeneration-defaults/model.R

structure(function (param = parameter_template) 
{
    if (is.data.frame(param)) {
        param_lower <- structure(param$lower, names = param$switch)
        param_upper <- structure(param$upper, names = param$switch)
        param <- structure(param$value, names = param$switch)
    }
    else {
        param_lower <- lapply(param, function(x) NA)
        param_upper <- lapply(param, function(x) NA)
    }
    stopifnot("retro_years" %in% names(param))
    stopifnot("fish.Linf" %in% names(param))
    stopifnot("fish.K" %in% names(param))
    stopifnot("fish.t0" %in% names(param))
    stopifnot("fish.lencv" %in% names(param))
    stopifnot("fish.init.scalar" %in% names(param))
    stopifnot("fish.init.1" %in% names(param))
    stopifnot("fish.init.2" %in% names(param))
    stopifnot("fish.init.3" %in% names(param))
    stopifnot("fish.init.4" %in% names(param))
    stopifnot("fish.init.5" %in% names(param))
    stopifnot("fish.init.6" %in% names(param))
    stopifnot("fish.init.7" %in% names(param))
    stopifnot("fish.init.8" %in% names(param))
    stopifnot("fish.init.9" %in% names(param))
    stopifnot("fish.init.10" %in% names(param))
    pt.fish.init <- list(`1` = param[["fish.init.1"]], `2` = param[["fish.init.2"]], `3` = param[["fish.init.3"]], `4` = param[["fish.init.4"]], `5` = param[["fish.init.5"]], `6` = param[["fish.init.6"]], `7` = param[["fish.init.7"]], `8` = param[["fish.init.8"]], `9` = param[["fish.init.9"]], `10` = param[["fish.init.10"]])
    stopifnot("fish.M.1" %in% names(param))
    stopifnot("fish.M.2" %in% names(param))
    stopifnot("fish.M.3" %in% names(param))
    stopifnot("fish.M.4" %in% names(param))
    stopifnot("fish.M.5" %in% names(param))
    stopifnot("fish.M.6" %in% names(param))
    stopifnot("fish.M.7" %in% names(param))
    stopifnot("fish.M.8" %in% names(param))
    stopifnot("fish.M.9" %in% names(param))
    stopifnot("fish.M.10" %in% names(param))
    pt.fish.M <- list(`1` = param[["fish.M.1"]], `2` = param[["fish.M.2"]], `3` = param[["fish.M.3"]], `4` = param[["fish.M.4"]], `5` = param[["fish.M.5"]], `6` = param[["fish.M.6"]], `7` = param[["fish.M.7"]], `8` = param[["fish.M.8"]], `9` = param[["fish.M.9"]], `10` = param[["fish.M.10"]])
    stopifnot("init.F" %in% names(param))
    stopifnot("recage" %in% names(param))
    stopifnot("fish.walpha" %in% names(param))
    stopifnot("fish.wbeta" %in% names(param))
    stopifnot("report_detail" %in% names(param))
    stopifnot("fish.comm.alpha" %in% names(param))
    stopifnot("fish.comm.l50" %in% names(param))
    stopifnot("fish.bbin" %in% names(param))
    stopifnot("fish.rec.proj" %in% names(param))
    stopifnot("fish.rec.1990" %in% names(param))
    stopifnot("fish.rec.1991" %in% names(param))
    stopifnot("fish.rec.1992" %in% names(param))
    stopifnot("fish.rec.1993" %in% names(param))
    stopifnot("fish.rec.1994" %in% names(param))
    stopifnot("fish.rec.1995" %in% names(param))
    stopifnot("fish.rec.1996" %in% names(param))
    stopifnot("fish.rec.1997" %in% names(param))
    stopifnot("fish.rec.1998" %in% names(param))
    stopifnot("fish.rec.1999" %in% names(param))
    stopifnot("fish.rec.2000" %in% names(param))
    pt.fish.rec <- list(`1990` = param[["fish.rec.1990"]], `1991` = param[["fish.rec.1991"]], `1992` = param[["fish.rec.1992"]], `1993` = param[["fish.rec.1993"]], `1994` = param[["fish.rec.1994"]], `1995` = param[["fish.rec.1995"]], `1996` = param[["fish.rec.1996"]], `1997` = param[["fish.rec.1997"]], `1998` = param[["fish.rec.1998"]], `1999` = param[["fish.rec.1999"]], `2000` = param[["fish.rec.2000"]])
    stopifnot("fish.rec.scalar" %in% names(param))
    stopifnot("adist_surveyindices_log_acoustic_dist_weight" %in% names(param))
    stopifnot("cdist_sumofsquares_comm_ldist_weight" %in% names(param))
    assert_msg <- function(expr, message) {
        if (isFALSE(expr)) {
            warning(message)
            return(TRUE)
        }
        return(FALSE)
    }
    dif_pmax <- function(a, b, scale) {
        logspace_add <- function(a, b) pmax(a, b) + log1p(exp(pmin(a, b) - pmax(a, b)))
        b <- as.vector(b)
        logspace_add(a * scale, b * scale)/scale
    }
    avoid_zero <- function(a) {
        dif_pmax(a, 0, 1000)
    }
    nvl <- function(...) {
        for (i in seq_len(...length())) if (!is.null(...elt(i))) 
            return(...elt(i))
        return(NULL)
    }
    normalize_vec <- function(a) {
        a/avoid_zero(sum(a))
    }
    as_numeric_arr <- function(x) x
    REPORT <- function(var) {
        var_name <- as.character(sys.call()[[2]])
        attr(nll, var_name) <<- if (var_name == "nll") 
            as.vector(var)
        else var
    }
    g3_cast_vector <- function(x) x
    intlookup_getdefault <- function(lookup, key, def) {
        if (is.environment(lookup)) {
            out <- lookup[[as.character(key)]]
        }
        else {
            out <- if (key < 1) 
                NULL
            else lookup[key][[1]]
        }
        return(if (is.null(out)) def else out)
    }
    nonconform_add <- function(base_ar, extra_ar) {
        base_ar + as.vector(extra_ar)
    }
    dif_pmin <- function(a, b, scale) {
        dif_pmax(a, b, -scale)
    }
    nonconform_mult <- function(base_ar, extra_ar) {
        base_ar * as.vector(extra_ar)
    }
    growth_bbinom <- function(delt_l, binn, beta) {
        alpha <- (beta * delt_l)/(binn - delt_l)
        x <- 0:binn
        na <- length(alpha)
        n <- length(x) - 1
        alpha <- rep(alpha, n + 1)
        x <- rep(x, each = na)
        val <- exp(lgamma(n + 1) + lgamma(alpha + beta) + lgamma(n - x + beta) + lgamma(x + alpha) - lgamma(n - x + 1) - lgamma(x + 1) - lgamma(n + alpha + beta) - lgamma(beta) - lgamma(alpha))
        dim(val) <- c(na, n + 1)
        return(val)
    }
    g3a_grow_vec_rotate <- function(vec, a) {
        out <- vapply(seq_len(a), function(i) vec[i:(i + length(vec) - 1)], numeric(length(vec)))
        out[is.na(out)] <- vec[length(vec)]
        out
    }
    g3a_grow_vec_extrude <- function(vec, a) {
        array(vec, dim = c(length(vec), a))
    }
    g3a_grow_matrix_wgt <- function(delta_w) {
        na <- dim(delta_w)[[1]]
        n <- dim(delta_w)[[2]] - 1
        wgt.matrix <- array(0, c(na, na))
        for (lg in 1:na) {
            if (lg == na) {
                wgt.matrix[lg, lg:na] <- delta_w[lg, 1:(na - lg + 1)]
            }
            else if (lg + n > na) {
                wgt.matrix[lg, lg:na] <- delta_w[lg, 1:(na - lg + 1)]
            }
            else {
                wgt.matrix[lg, lg:(n + lg)] <- delta_w[lg, ]
            }
        }
        return(wgt.matrix)
    }
    g3a_grow_matrix_len <- function(delta_l) {
        na <- dim(delta_l)[[1]]
        n <- dim(delta_l)[[2]] - 1
        growth.matrix <- array(0, c(na, na))
        for (lg in 1:na) {
            if (lg == na) {
                growth.matrix[na, na] <- sum(delta_l[lg, ])
            }
            else if (lg + n > na) {
                growth.matrix[lg, lg:(na - 1)] <- delta_l[lg, 1:(na - lg)]
                growth.matrix[lg, na] <- sum(delta_l[lg, (na - lg + 1):(n + 1)])
            }
            else {
                growth.matrix[lg, lg:(n + lg)] <- delta_l[lg, ]
            }
        }
        return(growth.matrix)
    }
    g3a_grow_apply <- function(growth.matrix, wgt.matrix, input_num, input_wgt) {
        na <- dim(growth.matrix)[[1]]
        growth.matrix <- growth.matrix * as.vector(input_num)
        wgt.matrix <- growth.matrix * (wgt.matrix + as.vector(input_wgt))
        growth.matrix.sum <- colSums(growth.matrix)
        return(array(c(growth.matrix.sum, colSums(wgt.matrix)/avoid_zero(growth.matrix.sum)), dim = c(na, 2)))
    }
    ratio_add_vec <- function(orig_vec, orig_amount, new_vec, new_amount) {
        (orig_vec * orig_amount + new_vec * new_amount)/avoid_zero(orig_amount + new_amount)
    }
    surveyindices_linreg <- function(N, I, fixed_alpha, fixed_beta) {
        meanI <- mean(I)
        meanN <- mean(N)
        beta <- if (is.nan(fixed_beta)) 
            sum((I - meanI) * (N - meanN))/avoid_zero(sum((N - meanN)^2))
        else fixed_beta
        alpha <- if (is.nan(fixed_alpha)) 
            meanI - beta * meanN
        else fixed_alpha
        return(c(alpha, beta))
    }
    cur_time <- -1L
    stopifnot("project_years" %in% names(param))
    project_years <- param[["project_years"]]
    cur_year <- 0L
    start_year <- 1990L
    step_lengths <- 12L
    step_count <- length(step_lengths)
    cur_year_projection <- FALSE
    end_year <- 2000L
    cur_step <- 0L
    cur_step_size <- step_lengths[[1]]/12
    cur_step_final <- FALSE
    fish__minage <- 1L
    fish__maxage <- 10L
    fish__area <- 1L
    fish__num <- array(0, dim = c(length = 6L, area = 1L, age = 10L), dimnames = list(length = c("50:60", "60:70", "70:80", "80:90", "90:100", "100:Inf"), area = "all", age = c("age1", "age2", "age3", "age4", "age5", "age6", "age7", "age8", "age9", "age10")))
    fish__wgt <- array(1, dim = c(length = 6L, area = 1L, age = 10L), dimnames = list(length = c("50:60", "60:70", "70:80", "80:90", "90:100", "100:Inf"), area = "all", age = c("age1", "age2", "age3", "age4", "age5", "age6", "age7", "age8", "age9", "age10")))
    as_integer <- as.integer
    retro_years <- param[["retro_years"]]
    total_steps <- length(step_lengths) * (end_year - as_integer(retro_years) - start_year + as_integer(project_years)) + length(step_lengths) - 1L
    dstart_fish__num <- array(0, dim = c(length = 6L, area = 1L, age = 10L, time = as_integer(total_steps + 1)), dimnames = list(length = c("50:60", "60:70", "70:80", "80:90", "90:100", "100:Inf"), area = "all", age = c("age1", "age2", "age3", "age4", "age5", "age6", "age7", "age8", "age9", "age10"), time = attributes(gen_dimnames(param))[["time"]]))
    dstart_fish__wgt <- array(1, dim = c(length = 6L, area = 1L, age = 10L, time = as_integer(total_steps + 1)), dimnames = list(length = c("50:60", "60:70", "70:80", "80:90", "90:100", "100:Inf"), area = "all", age = c("age1", "age2", "age3", "age4", "age5", "age6", "age7", "age8", "age9", "age10"), time = attributes(gen_dimnames(param))[["time"]]))
    suit_fish_comm__report <- array(NA, dim = c(length = 6L), dimnames = list(length = c("50:60", "60:70", "70:80", "80:90", "90:100", "100:Inf")))
    adist_surveyindices_log_acoustic_dist_model__wgt <- array(0, dim = c(length = 1L, time = 11L, area = 1L), dimnames = list(length = "0:Inf", time = c("1990", "1991", "1992", "1993", "1994", "1995", "1996", "1997", "1998", "1999", "2000"), area = "all"))
    cdist_sumofsquares_comm_ldist_model__wgt <- array(0, dim = c(length = 5L, time = 11L, area = 1L), dimnames = list(length = c("50:60", "60:70", "70:80", "80:90", "90:Inf"), time = c("1990", "1991", "1992", "1993", "1994", "1995", "1996", "1997", "1998", "1999", "2000"), area = "all"))
    detail_fish__predby_comm <- array(NA, dim = c(length = 6L, area = 1L, age = 10L, time = as_integer(total_steps + 1)), dimnames = list(length = c("50:60", "60:70", "70:80", "80:90", "90:100", "100:Inf"), area = "all", age = c("age1", "age2", "age3", "age4", "age5", "age6", "age7", "age8", "age9", "age10"), time = attributes(gen_dimnames(param))[["time"]]))
    detail_fish__renewalnum <- array(0, dim = c(length = 6L, area = 1L, age = 10L, time = as_integer(total_steps + 1)), dimnames = list(length = c("50:60", "60:70", "70:80", "80:90", "90:100", "100:Inf"), area = "all", age = c("age1", "age2", "age3", "age4", "age5", "age6", "age7", "age8", "age9", "age10"), time = attributes(gen_dimnames(param))[["time"]]))
    detail_fish_comm__cons <- array(NA, dim = c(length = 6L, area = 1L, age = 10L, time = as_integer(total_steps + 1)), dimnames = list(length = c("50:60", "60:70", "70:80", "80:90", "90:100", "100:Inf"), area = "all", age = c("age1", "age2", "age3", "age4", "age5", "age6", "age7", "age8", "age9", "age10"), time = attributes(gen_dimnames(param))[["time"]]))
    detail_fish_comm__suit <- array(NA, dim = c(length = 6L, area = 1L, age = 10L, time = as_integer(total_steps + 1)), dimnames = list(length = c("50:60", "60:70", "70:80", "80:90", "90:100", "100:Inf"), area = "all", age = c("age1", "age2", "age3", "age4", "age5", "age6", "age7", "age8", "age9", "age10"), time = attributes(gen_dimnames(param))[["time"]]))
    nll <- 0
    nll_adist_surveyindices_log_acoustic_dist__weight <- array(0, dim = c(time = as_integer(total_steps + 1L)), dimnames = list(time = attributes(gen_dimnames(param))[["time"]]))
    nll_adist_surveyindices_log_acoustic_dist__wgt <- array(0, dim = c(time = as_integer(total_steps + 1L)), dimnames = list(time = attributes(gen_dimnames(param))[["time"]]))
    nll_cdist_sumofsquares_comm_ldist__weight <- array(0, dim = c(time = as_integer(total_steps + 1L)), dimnames = list(time = attributes(gen_dimnames(param))[["time"]]))
    nll_cdist_sumofsquares_comm_ldist__wgt <- array(0, dim = c(time = as_integer(total_steps + 1L)), dimnames = list(time = attributes(gen_dimnames(param))[["time"]]))
    nll_understocking__wgt <- array(0, dim = c(time = as_integer(total_steps + 1L)), dimnames = list(time = attributes(gen_dimnames(param))[["time"]]))
    fish__totalpredate <- array(NA, dim = c(length = 6L, area = 1L, age = 10L), dimnames = list(length = c("50:60", "60:70", "70:80", "80:90", "90:100", "100:Inf"), area = "all", age = c("age1", "age2", "age3", "age4", "age5", "age6", "age7", "age8", "age9", "age10")))
    comm__totalsuit <- array(NA, dim = c(area = 1L), dimnames = list(area = "all"))
    fish_comm__suit <- array(NA, dim = c(length = 6L, area = 1L, age = 10L), dimnames = list(length = c("50:60", "60:70", "70:80", "80:90", "90:100", "100:Inf"), area = "all", age = c("age1", "age2", "age3", "age4", "age5", "age6", "age7", "age8", "age9", "age10")))
    comm__area <- 1L
    fish_comm__cons <- array(NA, dim = c(length = 6L, area = 1L, age = 10L), dimnames = list(length = c("50:60", "60:70", "70:80", "80:90", "90:100", "100:Inf"), area = "all", age = c("age1", "age2", "age3", "age4", "age5", "age6", "age7", "age8", "age9", "age10")))
    intlookup_zip <- function(keys, values) {
        if (min(keys) > 0 && max(keys) < 1e+05) {
            out <- list()
            out[as.integer(keys)] <- as.list(values)
        }
        else {
            out <- as.list(values)
            names(out) <- keys
            out <- as.environment(out)
        }
        attr(out, "key_var") <- deparse(sys.call()[[2]])
        attr(out, "value_var") <- deparse(sys.call()[[3]])
        return(out)
    }
    comm_landings <- intlookup_zip(comm_landings_keys, comm_landings_values)
    fish__consratio <- array(NA, dim = c(length = 6L, area = 1L, age = 10L), dimnames = list(length = c("50:60", "60:70", "70:80", "80:90", "90:100", "100:Inf"), area = "all", age = c("age1", "age2", "age3", "age4", "age5", "age6", "age7", "age8", "age9", "age10")))
    fish__overconsumption <- structure(0, desc = "Total overconsumption of fish")
    fish__consconv <- array(NA, dim = c(length = 6L, area = 1L, age = 10L), dimnames = list(length = c("50:60", "60:70", "70:80", "80:90", "90:100", "100:Inf"), area = "all", age = c("age1", "age2", "age3", "age4", "age5", "age6", "age7", "age8", "age9", "age10")))
    fish__predby_comm <- array(NA, dim = c(length = 6L, area = 1L, age = 10L), dimnames = list(length = c("50:60", "60:70", "70:80", "80:90", "90:100", "100:Inf"), area = "all", age = c("age1", "age2", "age3", "age4", "age5", "age6", "age7", "age8", "age9", "age10")))
    fish__growth_lastcalc <- -1L
    fish__growth_l <- array(NA, dim = c(length = 6L, delta = 6L), dimnames = list(length = c("50:60", "60:70", "70:80", "80:90", "90:100", "100:Inf"), delta = c("0", "1", "2", "3", "4", "5")))
    fish__plusdl <- 10
    fish__growth_w <- array(NA, dim = c(length = 6L, delta = 6L), dimnames = list(length = c("50:60", "60:70", "70:80", "80:90", "90:100", "100:Inf"), delta = c("0", "1", "2", "3", "4", "5")))
    fish__prevtotal <- 0
    fish__renewalnum <- array(0, dim = c(length = 6L, area = 1L, age = 10L), dimnames = list(length = c("50:60", "60:70", "70:80", "80:90", "90:100", "100:Inf"), area = "all", age = c("age1", "age2", "age3", "age4", "age5", "age6", "age7", "age8", "age9", "age10")))
    fish__renewalwgt <- array(0, dim = c(length = 6L, area = 1L, age = 10L), dimnames = list(length = c("50:60", "60:70", "70:80", "80:90", "90:100", "100:Inf"), area = "all", age = c("age1", "age2", "age3", "age4", "age5", "age6", "age7", "age8", "age9", "age10")))
    adist_surveyindices_log_acoustic_dist_model__area <- 1L
    adist_surveyindices_log_acoustic_dist_model__times <- intlookup_zip(adist_surveyindices_log_acoustic_dist_model__times_keys, adist_surveyindices_log_acoustic_dist_model__times_values)
    fish_adist_surveyindices_log_acoustic_dist_model_lgmatrix <- array(1, dim = c(1L, 6L), dimnames = NULL)
    adist_surveyindices_log_acoustic_dist_obs__area <- 1L
    cdist_sumofsquares_comm_ldist_model__area <- 1L
    cdist_sumofsquares_comm_ldist_model__times <- intlookup_zip(cdist_sumofsquares_comm_ldist_model__times_keys, cdist_sumofsquares_comm_ldist_model__times_values)
    cdist_sumofsquares_comm_ldist_obs__area <- 1L
    cdist_sumofsquares_comm_ldist_obs__times <- intlookup_zip(cdist_sumofsquares_comm_ldist_obs__times_keys, cdist_sumofsquares_comm_ldist_obs__times_values)
    g3l_understocking_total <- 0
    nll_understocking__weight <- array(0, dim = c(time = as_integer(total_steps + 1L)), dimnames = list(time = attributes(gen_dimnames(param))[["time"]]))
    while (TRUE) {
        {
            comment("g3a_time: Start of time period")
            cur_time <- cur_time + 1L
            if (cur_time == 0 && assert_msg(param[["retro_years"]] >= 0, "retro_years must be >= 0")) 
                return(NaN)
            if (cur_time == 0 && assert_msg(project_years >= 0, "project_years must be >= 0")) 
                return(NaN)
            cur_year <- start_year + (cur_time%/%step_count)
            cur_year_projection <- cur_year > end_year - param[["retro_years"]]
            cur_step <- (cur_time%%step_count) + 1L
            cur_step_size <- step_lengths[[cur_step]]/12
            cur_step_final <- cur_step == step_count
        }
        {
            comment("g3a_initialconditions for fish")
            for (age in seq(fish__minage, fish__maxage, by = 1)) if (cur_time == 0L) {
                fish__age_idx <- age - fish__minage + 1L
                area <- fish__area
                fish__area_idx <- (1L)
                ren_dnorm <- dnorm(fish__midlen, (param[["fish.Linf"]] * (1 - exp(-1 * param[["fish.K"]] * ((age - cur_step_size) - param[["fish.t0"]])))), avoid_zero(((param[["fish.Linf"]] * (1 - exp(-1 * param[["fish.K"]] * ((age - cur_step_size) - param[["fish.t0"]])))) * param[["fish.lencv"]])))
                factor <- (param[["fish.init.scalar"]] * nvl(pt.fish.init[[paste(age, sep = ".")]], {
                  warning("No value found in g3_param_table fish.init, ifmissing not specified")
                  NaN
                }) * exp(-1 * (nvl(pt.fish.M[[paste(age, sep = ".")]], {
                  warning("No value found in g3_param_table fish.M, ifmissing not specified")
                  NaN
                }) + param[["init.F"]]) * (age - param[["recage"]])))
                {
                  fish__num[, fish__area_idx, fish__age_idx] <- normalize_vec(ren_dnorm) * 10000 * factor
                  fish__wgt[, fish__area_idx, fish__age_idx] <- param[["fish.walpha"]] * fish__midlen^param[["fish.wbeta"]]
                }
            }
        }
        if ((cur_time <= total_steps && param[["report_detail"]] == 1)) 
            dstart_fish__num[, , , cur_time + 1] <- as_numeric_arr(fish__num)
        if ((cur_time <= total_steps && param[["report_detail"]] == 1)) 
            dstart_fish__wgt[, , , cur_time + 1] <- as_numeric_arr(fish__wgt)
        if (cur_time == 0L) {
            suit_fish_comm__report[] <- 1/(1 + exp(-param[["fish.comm.alpha"]] * (fish__midlen - param[["fish.comm.l50"]])))
            REPORT(suit_fish_comm__report)
        }
        if (reporting_enabled > 0L && cur_time > total_steps) 
            REPORT(adist_surveyindices_log_acoustic_dist_model__params)
        if (reporting_enabled > 0L && cur_time > total_steps) 
            REPORT(adist_surveyindices_log_acoustic_dist_model__wgt)
        if (reporting_enabled > 0L && cur_time > total_steps) 
            REPORT(adist_surveyindices_log_acoustic_dist_obs__wgt)
        if (reporting_enabled > 0L && cur_time > total_steps) 
            REPORT(cdist_sumofsquares_comm_ldist_model__wgt)
        if (reporting_enabled > 0L && cur_time > total_steps) 
            REPORT(cdist_sumofsquares_comm_ldist_obs__wgt)
        if (reporting_enabled > 0L && cur_time > total_steps) 
            REPORT(detail_fish__predby_comm)
        if (reporting_enabled > 0L && cur_time > total_steps) 
            REPORT(detail_fish__renewalnum)
        if (reporting_enabled > 0L && cur_time > total_steps) 
            REPORT(detail_fish_comm__cons)
        if (reporting_enabled > 0L && cur_time > total_steps) 
            REPORT(detail_fish_comm__suit)
        if (reporting_enabled > 0L && cur_time > total_steps) 
            REPORT(dstart_fish__num)
        if (reporting_enabled > 0L && cur_time > total_steps) 
            REPORT(dstart_fish__wgt)
        if (reporting_enabled > 0L && cur_time > total_steps) 
            REPORT(nll)
        if (reporting_enabled > 0L && cur_time > total_steps) 
            REPORT(nll_adist_surveyindices_log_acoustic_dist__weight)
        if (reporting_enabled > 0L && cur_time > total_steps) 
            REPORT(nll_adist_surveyindices_log_acoustic_dist__wgt)
        if (reporting_enabled > 0L && cur_time > total_steps) 
            REPORT(nll_cdist_sumofsquares_comm_ldist__weight)
        if (reporting_enabled > 0L && cur_time > total_steps) 
            REPORT(nll_cdist_sumofsquares_comm_ldist__wgt)
        if (reporting_enabled > 0L && cur_time > total_steps) 
            REPORT(nll_understocking__wgt)
        if (reporting_enabled > 0L && cur_time > total_steps) 
            REPORT(step_lengths)
        {
            if (cur_time > total_steps) 
                return(nll)
        }
        fish__totalpredate[] <- 0
        comm__totalsuit[] <- 0
        {
            suitability <- g3_cast_vector(suit_fish_comm__report[])
            {
                comment("g3a_predate for comm predating fish")
                fish_comm__suit[] <- 0
                for (age in seq(fish__minage, fish__maxage, by = 1)) {
                  fish__age_idx <- age - fish__minage + 1L
                  area <- fish__area
                  fish__area_idx <- (1L)
                  if (area == comm__area) {
                    catchability <- (suitability * fish__num[, fish__area_idx, fish__age_idx])
                    comm__area_idx <- (1L)
                    predator_area <- area
                    {
                      comment("Collect all suitable fish biomass for comm")
                      fish_comm__suit[, fish__area_idx, fish__age_idx] <- catchability
                      comm__totalsuit[comm__area_idx] <- comm__totalsuit[comm__area_idx] + sum(fish_comm__suit[, fish__area_idx, fish__age_idx])
                    }
                  }
                }
            }
        }
        {
            comment("Scale comm catch of fish by total expected catch")
            fish_comm__cons[] <- 0
            for (age in seq(fish__minage, fish__maxage, by = 1)) {
                fish__age_idx <- age - fish__minage + 1L
                area <- fish__area
                fish__area_idx <- (1L)
                if (area == comm__area) {
                  comm__area_idx <- (1L)
                  predator_area <- area
                  total_predsuit <- comm__totalsuit[comm__area_idx]
                  fish_comm__cons[, fish__area_idx, fish__age_idx] <- fish_comm__suit[, fish__area_idx, fish__age_idx] * ((if (area != 1L) 
                    0
                  else intlookup_getdefault(comm_landings, cur_year, 0))/total_predsuit) * fish__wgt[, fish__area_idx, fish__age_idx]
                }
            }
            fish__totalpredate[] <- nonconform_add(fish__totalpredate[], fish_comm__cons[, , ])
        }
        {
            comment("Calculate fish overconsumption coefficient")
            comment("Apply overconsumption to fish")
            fish__consratio <- fish__totalpredate/avoid_zero(fish__num * fish__wgt)
            fish__consratio <- dif_pmin(fish__consratio, 0.95, 1000)
            fish__overconsumption <- sum(fish__totalpredate)
            fish__consconv <- 1/avoid_zero(fish__totalpredate)
            fish__totalpredate <- (fish__num * fish__wgt) * fish__consratio
            fish__overconsumption <- fish__overconsumption - sum(fish__totalpredate)
            fish__consconv <- fish__consconv * fish__totalpredate
            fish__num <- fish__num * (1 - fish__consratio)
        }
        {
            comment("Apply overconsumption to fish_comm__cons")
            fish_comm__cons <- nonconform_mult(fish_comm__cons, fish__consconv)
        }
        {
            fish__predby_comm[] <- 0
            fish__predby_comm[] <- nonconform_add(fish__predby_comm, fish_comm__cons[, , ])
        }
        {
            comment("Natural mortality for fish")
            for (age in seq(fish__minage, fish__maxage, by = 1)) {
                fish__age_idx <- age - fish__minage + 1L
                area <- fish__area
                fish__area_idx <- (1L)
                fish__num[, fish__area_idx, fish__age_idx] <- fish__num[, fish__area_idx, fish__age_idx] * exp(-(nvl(pt.fish.M[[paste(age, sep = ".")]], {
                  warning("No value found in g3_param_table fish.M, ifmissing not specified")
                  NaN
                })) * cur_step_size)
            }
        }
        {
            growth_delta_l <- if (fish__growth_lastcalc == floor(cur_step_size * 12L)) 
                fish__growth_l
            else (fish__growth_l[] <- growth_bbinom(avoid_zero(avoid_zero((param[["fish.Linf"]] - fish__midlen) * (1 - exp(-(param[["fish.K"]]) * cur_step_size)))/fish__plusdl), 5L, avoid_zero(param[["fish.bbin"]])))
            growth_delta_w <- if (fish__growth_lastcalc == floor(cur_step_size * 12L)) 
                fish__growth_w
            else (fish__growth_w[] <- (g3a_grow_vec_rotate(fish__midlen^param[["fish.wbeta"]], 5L + 1) - g3a_grow_vec_extrude(fish__midlen^param[["fish.wbeta"]], 5L + 1)) * param[["fish.walpha"]])
            growthmat_w <- g3a_grow_matrix_wgt(growth_delta_w)
            growthmat_l <- g3a_grow_matrix_len(growth_delta_l)
            {
                comment("g3a_grow for fish")
                for (age in seq(fish__minage, fish__maxage, by = 1)) {
                  fish__age_idx <- age - fish__minage + 1L
                  area <- fish__area
                  fish__area_idx <- (1L)
                  growthresult <- g3a_grow_apply(growthmat_l, growthmat_w, fish__num[, fish__area_idx, fish__age_idx], fish__wgt[, fish__area_idx, fish__age_idx])
                  {
                    if (FALSE) 
                      fish__prevtotal <- sum(fish__num[, fish__area_idx, fish__age_idx])
                    comment("Update fish using delta matrices")
                    fish__num[, fish__area_idx, fish__age_idx] <- growthresult[, (1)]
                    fish__wgt[, fish__area_idx, fish__age_idx] <- growthresult[, (2)]
                    if (FALSE) 
                      assert_msg(~abs(fish__prevtotal - sum(fish__num[, fish__area_idx, fish__age_idx])) < 1e-04, "g3a_growmature: fish__num totals are not the same before and after growth")
                  }
                }
                fish__growth_lastcalc <- floor(cur_step_size * 12L)
            }
        }
        {
            factor <- (nvl(pt.fish.rec[[paste(cur_year, sep = ".")]], param[["fish.rec.proj"]]) * param[["fish.rec.scalar"]])
            {
                comment("g3a_renewal for fish")
                for (age in seq(fish__minage, fish__maxage, by = 1)) if (age == fish__minage && cur_step == 1) {
                  fish__age_idx <- age - fish__minage + 1L
                  area <- fish__area
                  fish__area_idx <- (1L)
                  ren_dnorm <- dnorm(fish__midlen, (param[["fish.Linf"]] * (1 - exp(-1 * param[["fish.K"]] * (age - param[["fish.t0"]])))), avoid_zero(((param[["fish.Linf"]] * (1 - exp(-1 * param[["fish.K"]] * (age - param[["fish.t0"]])))) * param[["fish.lencv"]])))
                  {
                    fish__renewalnum[, fish__area_idx, fish__age_idx] <- normalize_vec(ren_dnorm) * 10000 * factor
                    fish__renewalwgt[, fish__area_idx, fish__age_idx] <- param[["fish.walpha"]] * fish__midlen^param[["fish.wbeta"]]
                    comment("Add result to fish")
                    fish__wgt[, fish__area_idx, fish__age_idx] <- ratio_add_vec(fish__wgt[, fish__area_idx, fish__age_idx], fish__num[, fish__area_idx, fish__age_idx], fish__renewalwgt[, fish__area_idx, fish__age_idx], fish__renewalnum[, fish__area_idx, fish__age_idx])
                    fish__num[, fish__area_idx, fish__age_idx] <- fish__num[, fish__area_idx, fish__age_idx] + fish__renewalnum[, fish__area_idx, fish__age_idx]
                  }
                }
            }
        }
        if (param[["adist_surveyindices_log_acoustic_dist_weight"]] > 0) {
            comment("g3l_abundancedistribution_surveyindices_log: Collect abundance from fish for adist_surveyindices_log_acoustic_dist")
            for (age in seq(fish__minage, fish__maxage, by = 1)) {
                fish__age_idx <- age - fish__minage + 1L
                area <- fish__area
                fish__area_idx <- (1L)
                if (area == adist_surveyindices_log_acoustic_dist_model__area) {
                  adist_surveyindices_log_acoustic_dist_model__area_idx <- (1L)
                  adist_surveyindices_log_acoustic_dist_model__time_idx <- intlookup_getdefault(adist_surveyindices_log_acoustic_dist_model__times, (cur_year * 100L + cur_step * 0L), -1L)
                  if (adist_surveyindices_log_acoustic_dist_model__time_idx >= (1L)) {
                    comment("Convert fish to wgt")
                    adist_surveyindices_log_acoustic_dist_model__wgt[, adist_surveyindices_log_acoustic_dist_model__time_idx, adist_surveyindices_log_acoustic_dist_model__area_idx] <- adist_surveyindices_log_acoustic_dist_model__wgt[, adist_surveyindices_log_acoustic_dist_model__time_idx, adist_surveyindices_log_acoustic_dist_model__area_idx] + as.vector(fish_adist_surveyindices_log_acoustic_dist_model_lgmatrix %*% (fish__num[, fish__area_idx, fish__age_idx] * fish__wgt[, fish__area_idx, fish__age_idx]))
                  }
                }
            }
        }
        {
            adist_surveyindices_log_acoustic_dist_model__max_time_idx <- (11L)
            {
                comment("g3l_abundancedistribution_surveyindices_log: Compare adist_surveyindices_log_acoustic_dist_model to adist_surveyindices_log_acoustic_dist_obs")
                if (param[["adist_surveyindices_log_acoustic_dist_weight"]] > 0 && cur_step_final) {
                  area <- adist_surveyindices_log_acoustic_dist_model__area
                  adist_surveyindices_log_acoustic_dist_model__area_idx <- (1L)
                  adist_surveyindices_log_acoustic_dist_model__time_idx <- intlookup_getdefault(adist_surveyindices_log_acoustic_dist_model__times, (cur_year * 100L + cur_step * 0L), -1L)
                  if (adist_surveyindices_log_acoustic_dist_model__time_idx >= (1L)) 
                    if (area == adist_surveyindices_log_acoustic_dist_obs__area) {
                      adist_surveyindices_log_acoustic_dist_obs__area_idx <- (1L)
                      {
                        adist_surveyindices_log_acoustic_dist_model__params <- if (adist_surveyindices_log_acoustic_dist_model__time_idx != adist_surveyindices_log_acoustic_dist_model__max_time_idx) 
                          adist_surveyindices_log_acoustic_dist_model__params
                        else surveyindices_linreg(log(avoid_zero(adist_surveyindices_log_acoustic_dist_model__wgt[, , adist_surveyindices_log_acoustic_dist_model__area_idx])), log(avoid_zero(adist_surveyindices_log_acoustic_dist_obs__wgt[, , adist_surveyindices_log_acoustic_dist_obs__area_idx])), NaN, 1)
                        {
                          cur_cdist_nll <- if (adist_surveyindices_log_acoustic_dist_model__time_idx != adist_surveyindices_log_acoustic_dist_model__max_time_idx) 
                            0
                          else sum((adist_surveyindices_log_acoustic_dist_model__params[[1]] + adist_surveyindices_log_acoustic_dist_model__params[[2]] * log(avoid_zero(adist_surveyindices_log_acoustic_dist_model__wgt[, , adist_surveyindices_log_acoustic_dist_model__area_idx])) - log(avoid_zero(adist_surveyindices_log_acoustic_dist_obs__wgt[, , adist_surveyindices_log_acoustic_dist_obs__area_idx])))^2)
                          {
                            nll <- nll + param[["adist_surveyindices_log_acoustic_dist_weight"]] * cur_cdist_nll
                            nll_adist_surveyindices_log_acoustic_dist__wgt[cur_time + 1L] <- nll_adist_surveyindices_log_acoustic_dist__wgt[cur_time + 1L] + cur_cdist_nll
                            nll_adist_surveyindices_log_acoustic_dist__weight[cur_time + 1L] <- param[["adist_surveyindices_log_acoustic_dist_weight"]]
                          }
                        }
                      }
                    }
                }
            }
        }
        if (param[["cdist_sumofsquares_comm_ldist_weight"]] > 0) {
            comment("g3l_catchdistribution_sumofsquares: Collect catch from comm/fish for cdist_sumofsquares_comm_ldist")
            for (age in seq(fish__minage, fish__maxage, by = 1)) {
                fish__age_idx <- age - fish__minage + 1L
                area <- fish__area
                fish__area_idx <- (1L)
                if (area == cdist_sumofsquares_comm_ldist_model__area) {
                  cdist_sumofsquares_comm_ldist_model__area_idx <- (1L)
                  cdist_sumofsquares_comm_ldist_model__time_idx <- intlookup_getdefault(cdist_sumofsquares_comm_ldist_model__times, (cur_year * 100L + cur_step * 0L), -1L)
                  if (cdist_sumofsquares_comm_ldist_model__time_idx >= (1L)) {
                    comment("Convert fish_comm to wgt")
                    cdist_sumofsquares_comm_ldist_model__wgt[, cdist_sumofsquares_comm_ldist_model__time_idx, cdist_sumofsquares_comm_ldist_model__area_idx] <- cdist_sumofsquares_comm_ldist_model__wgt[, cdist_sumofsquares_comm_ldist_model__time_idx, cdist_sumofsquares_comm_ldist_model__area_idx] + as.vector(fish_comm_cdist_sumofsquares_comm_ldist_model_lgmatrix %*% fish_comm__cons[, fish__area_idx, fish__age_idx])
                  }
                }
            }
        }
        {
            comment("g3l_catchdistribution_sumofsquares: Compare cdist_sumofsquares_comm_ldist_model to cdist_sumofsquares_comm_ldist_obs")
            if (param[["cdist_sumofsquares_comm_ldist_weight"]] > 0 && cur_step_final) {
                area <- cdist_sumofsquares_comm_ldist_model__area
                cdist_sumofsquares_comm_ldist_model__area_idx <- (1L)
                cdist_sumofsquares_comm_ldist_model__time_idx <- intlookup_getdefault(cdist_sumofsquares_comm_ldist_model__times, (cur_year * 100L + cur_step * 0L), -1L)
                if (cdist_sumofsquares_comm_ldist_model__time_idx >= (1L)) 
                  if (area == cdist_sumofsquares_comm_ldist_obs__area) {
                    cdist_sumofsquares_comm_ldist_model__sstotal <- avoid_zero(sum(cdist_sumofsquares_comm_ldist_model__wgt[, cdist_sumofsquares_comm_ldist_model__time_idx, cdist_sumofsquares_comm_ldist_model__area_idx]))
                    {
                      cdist_sumofsquares_comm_ldist_obs__area_idx <- (1L)
                      cdist_sumofsquares_comm_ldist_obs__time_idx <- intlookup_getdefault(cdist_sumofsquares_comm_ldist_obs__times, (cur_year * 100L + cur_step * 0L), -1L)
                      if (cdist_sumofsquares_comm_ldist_obs__time_idx >= (1L)) {
                        cdist_sumofsquares_comm_ldist_obs__sstotal <- avoid_zero(sum(cdist_sumofsquares_comm_ldist_obs__wgt[, cdist_sumofsquares_comm_ldist_obs__time_idx, cdist_sumofsquares_comm_ldist_obs__area_idx]))
                        cur_cdist_nll <- sum((((cdist_sumofsquares_comm_ldist_model__wgt[, cdist_sumofsquares_comm_ldist_model__time_idx, cdist_sumofsquares_comm_ldist_model__area_idx]/cdist_sumofsquares_comm_ldist_model__sstotal) - (cdist_sumofsquares_comm_ldist_obs__wgt[, cdist_sumofsquares_comm_ldist_obs__time_idx, cdist_sumofsquares_comm_ldist_obs__area_idx]/cdist_sumofsquares_comm_ldist_obs__sstotal))^2))
                        {
                          nll <- nll + param[["cdist_sumofsquares_comm_ldist_weight"]] * cur_cdist_nll
                          nll_cdist_sumofsquares_comm_ldist__wgt[cur_time + 1L] <- nll_cdist_sumofsquares_comm_ldist__wgt[cur_time + 1L] + cur_cdist_nll
                          nll_cdist_sumofsquares_comm_ldist__weight[cur_time + 1L] <- param[["cdist_sumofsquares_comm_ldist_weight"]]
                        }
                      }
                    }
                  }
            }
        }
        {
            comment("Reset understocking total")
            g3l_understocking_total <- 0
        }
        {
            comment("g3l_understocking for fish")
            comment("Add understocking from fish as biomass to nll")
            g3l_understocking_total <- g3l_understocking_total + fish__overconsumption
        }
        {
            comment("g3l_understocking: Combine and add to nll")
            g3l_understocking_total <- g3l_understocking_total^2
            nll <- nll + 1e+08 * g3l_understocking_total
            nll_understocking__wgt[cur_time + 1L] <- nll_understocking__wgt[cur_time + 1L] + g3l_understocking_total
            nll_understocking__weight[cur_time + 1L] <- 1e+08
        }
        if (param[["report_detail"]] == 1) 
            detail_fish__predby_comm[, , , cur_time + 1] <- as_numeric_arr(fish__predby_comm)
        if (param[["report_detail"]] == 1) 
            detail_fish__renewalnum[, , , cur_time + 1] <- as_numeric_arr(fish__renewalnum)
        if (param[["report_detail"]] == 1) 
            detail_fish_comm__cons[, , , cur_time + 1] <- as_numeric_arr(fish_comm__cons)
        if (param[["report_detail"]] == 1) 
            detail_fish_comm__suit[, , , cur_time + 1] <- as_numeric_arr(fish_comm__suit)
        if (cur_step_final) {
            comment("g3a_age for fish")
            for (age in seq(fish__maxage, fish__minage, by = -1)) {
                fish__age_idx <- age - fish__minage + 1L
                {
                  comment("Check stock has remained finite for this step")
                  if (age == fish__maxage) {
                    comment("Oldest fish is a plus-group, combine with younger individuals")
                    fish__wgt[, , fish__age_idx] <- ratio_add_vec(fish__wgt[, , fish__age_idx], fish__num[, , fish__age_idx], fish__wgt[, , fish__age_idx - 1], fish__num[, , fish__age_idx - 1])
                    fish__num[, , fish__age_idx] <- fish__num[, , fish__age_idx] + fish__num[, , fish__age_idx - 1]
                  }
                  else if (age == fish__minage) {
                    comment("Empty youngest fish age-group")
                    fish__num[, , fish__age_idx] <- 0
                  }
                  else {
                    comment("Move fish age-group to next one up")
                    fish__num[, , fish__age_idx] <- fish__num[, , fish__age_idx - 1]
                    fish__wgt[, , fish__age_idx] <- fish__wgt[, , fish__age_idx - 1]
                  }
                }
            }
        }
    }
}, class = c("g3_r", "function"), parameter_template = list(retro_years = 0, fish.Linf = 1, fish.K = 1, fish.t0 = 0, fish.lencv = 0.1, fish.init.scalar = 1, fish.init.1 = 1, fish.init.2 = 1, fish.init.3 = 1, fish.init.4 = 1, fish.init.5 = 1, fish.init.6 = 1, fish.init.7 = 1, fish.init.8 = 1, fish.init.9 = 1, fish.init.10 = 1, fish.M.1 = 0, fish.M.2 = 0, fish.M.3 = 0, fish.M.4 = 0, fish.M.5 = 0, fish.M.6 = 0, fish.M.7 = 0, fish.M.8 = 0, fish.M.9 = 0, fish.M.10 = 0, init.F = 0, recage = 0, fish.walpha = 0, 
    fish.wbeta = 0, report_detail = 1L, fish.comm.alpha = 0, fish.comm.l50 = 0, fish.bbin = 0, fish.rec.proj = 0, fish.rec.1990 = 0, fish.rec.1991 = 0, fish.rec.1992 = 0, fish.rec.1993 = 0, fish.rec.1994 = 0, fish.rec.1995 = 0, fish.rec.1996 = 0, fish.rec.1997 = 0, fish.rec.1998 = 0, fish.rec.1999 = 0, fish.rec.2000 = 0, fish.rec.scalar = 0, adist_surveyindices_log_acoustic_dist_weight = 1, cdist_sumofsquares_comm_ldist_weight = 1, project_years = 0))
gadget-framework/gadget3 documentation built on June 13, 2025, 5:06 a.m.