diff --git a/nimbleModel/.Rbuildignore b/nimbleModel/.Rbuildignore new file mode 100644 index 0000000..b2ba255 --- /dev/null +++ b/nimbleModel/.Rbuildignore @@ -0,0 +1 @@ +^include$ diff --git a/nimbleModel/DESCRIPTION b/nimbleModel/DESCRIPTION index 4c82d3d..12acd75 100644 --- a/nimbleModel/DESCRIPTION +++ b/nimbleModel/DESCRIPTION @@ -54,6 +54,7 @@ Collate: declFunBaseClass.R modelBaseClass.R instructions.R + instrRunners.R nimbleModel.R zzz.R diff --git a/nimbleModel/NAMESPACE b/nimbleModel/NAMESPACE index bc086ab..0e99885 100644 --- a/nimbleModel/NAMESPACE +++ b/nimbleModel/NAMESPACE @@ -90,6 +90,8 @@ export(makeMultiCopier) export(values) export("values<-") +export("values_indexed<-") +export(values_indexed) export(modelValues) export(modelValuesConf) @@ -110,4 +112,5 @@ export(setdiff_nodes) export(setNimbleModelOption) export(getNimbleModelOption) +export(withNimbleModelOptions) diff --git a/nimbleModel/R/copier.R b/nimbleModel/R/copier.R index 30d39b4..ac9f39c 100644 --- a/nimbleModel/R/copier.R +++ b/nimbleModel/R/copier.R @@ -65,12 +65,12 @@ multiCopier_nClass <- nCompiler::nClass( returnType = "numericVector", compileInfo = list( C_fun = function() { - cppLiteral("return flatViewGroup.copyIntoVector()") + cppLiteral("return getValues_()") # refreshes data pointers; see multiCopier_model_nC_base } ) ), # Note that this setValues works but is not actually used by nimble2 keyword processing - # which instead uses copiers->flatViewGroup.setValues_() which provides an operator=. + # which instead uses copiers->setValues_() which provides an operator=. # It is possible that this setValues() will be useful in other contexts, or if not that we # can eventually remove it to avoid confusion. setValues = nFunction( @@ -80,7 +80,33 @@ multiCopier_nClass <- nCompiler::nClass( }, compileInfo = list( C_fun = function(v = "numericVector") { - cppLiteral("flatViewGroup.copyFromVector(v)") + cppLiteral("setValues_() = v") # refreshes data pointers; see multiCopier_model_nC_base + } + ) + ), + # Get or set one entry (mode = "raw") or one flattened element (mode = "scalar"); i is 1-based. + # These choose the mode at run time, for use from R. nimble2 keyword processing instead + # calls getValues_raw_(i) etc. in multiCopier_model_nC_base directly, since its mode is fixed. + getValues_indexed = nFunction( + name = "getValues_indexed", + function(i = "integerScalar", mode = "string") { + cat("need uncompiled multiCopier getValues_indexed() implementation.") + }, + returnType = "numericVector", + compileInfo = list( + C_fun = function(i = "integerScalar", mode = "string") { + cppLiteral("return getValues_indexed_(i, mode)") # see multiCopier_model_nC_base + } + ) + ), + setValues_indexed = nFunction( + name = "setValues_indexed", + function(i = "integerScalar", v = "numericVector", mode = "string") { + cat("need uncompiled multiCopier setValues_indexed() implementation.") + }, + compileInfo = list( + C_fun = function(i = "integerScalar", v = "numericVector", mode = "string") { + cppLiteral("setValues_indexed_(i, v, mode)") # see multiCopier_model_nC_base } ) ) @@ -203,11 +229,15 @@ makeMultiCopier <- function(source, nodes, model = multiCopier_nClass$new(), modelValues = multiCopierMV_nClass$new() ) + # A single index i is the block i:i, which flattens identically. getRange <- function(indexRange) { - if (!inherits(indexRange, "indexRangeSequenceClass")) { - stop("In a copy operation, only contiguous index blocks are supported") + if (inherits(indexRange, "indexRangeSequenceClass")) { + return(c(indexRange$start, indexRange$end)) + } + if (inherits(indexRange, "indexRangeScalarClass")) { + return(c(indexRange$value, indexRange$value)) } - c(indexRange$start, indexRange$end) + stop("In a copy operation, only contiguous index blocks are supported") } for (i in seq_along(nodes)) { thisNode <- nodes[[i]] @@ -225,8 +255,23 @@ makeMultiCopier <- function(source, nodes, # Uncompiled implementations of values(model, nodes) and `values<-`(model, nodes, value). +# With mode = "scalar" or "node" and indexed `nodes`, as in values(model, nodes[i], mode = "scalar"), +# `nodes` and `i` are evaluated separately in the caller and passed to values_indexed(). +# With mode = "raw" (the default, as in nimble), nodes[i] is evaluated as usual, giving entry i of `nodes`. +# Without an index, "scalar" is the same as "raw", and "node" gives all the nodes' values in sorted order. + #' @export -values <- function(model, nodes) { +values <- function(model, nodes, mode = "raw") { + # Unindexed scalar mode is the same as raw mode: all the values, flattened. + if (identical(mode, "scalar") && !is_singly_indexed(substitute(nodes))) mode <- "raw" + # Unindexed node mode: all the nodes, in sorted order. + if (identical(mode, "node") && !is_singly_indexed(substitute(nodes))) { + return(model$getValues_nodes(nodes)) + } + if (!identical(mode, "raw")) { + base_and_index <- split_indexed_nodes(substitute(nodes), parent.frame(), mode) + return(values_indexed(model, base_and_index$nodes, base_and_index$i, mode)) + } if (is.character(nodes)) { nodes <- lapply(nodes, \(x) varRangeClass$new(x)) } @@ -241,17 +286,32 @@ values <- function(model, nodes) { } #' @export -`values<-` <- function(model, nodes, value) { +`values<-` <- function(model, nodes, mode = "raw", value) { + if (identical(mode, "scalar") && !is_singly_indexed(substitute(nodes))) mode <- "raw" + if (identical(mode, "node") && !is_singly_indexed(substitute(nodes))) { + model$setValues_nodes(nodes, value) + return(model) + } + if (!identical(mode, "raw")) { + base_and_index <- split_indexed_nodes(substitute(nodes), parent.frame(), mode) + values_indexed(model, base_and_index$nodes, base_and_index$i, mode) <- value + return(model) + } if (is.character(nodes)) { nodes <- lapply(nodes, \(x) varRangeClass$new(x)) } - iV_start <- 1 - replExprTemplate <- quote(model_stuff <- value[iV_start:iV_end]) + # First find each target expression and length, so a wrong-length value + # is an error before anything is assigned. + targetExprs <- vector("list", length(nodes)) + targetLengths <- numeric(length(nodes)) for (i in seq_along(nodes)) { thisNode <- nodes[[i]] expr <- thisNode$toExpr() if (is.name(expr)) { - thisLength <- length(model[[expr]]) + thisLength <- length(model[[as.character(expr)]]) + model_stuff <- quote(model$VAR) + model_stuff[[3]] <- expr + expr <- model_stuff } else { mm <- thisNode$getMinMax() sizes <- mm[, 2] - mm[, 1] + 1 @@ -271,11 +331,136 @@ values <- function(model, nodes) { model_stuff[[3]] <- varExpr expr[[2]] <- model_stuff } - iV_end <- iV_start + thisLength - 1 + targetExprs[[i]] <- expr + targetLengths[i] <- thisLength + } + if (length(value) != sum(targetLengths)) { + stop("values(model, nodes) <- value: value has length ", length(value), + " but nodes have total length ", sum(targetLengths), ".") + } + iV_start <- 1 + replExprTemplate <- quote(model_stuff <- value[iV_start:iV_end]) + for (i in seq_along(nodes)) { + iV_end <- iV_start + targetLengths[i] - 1 replExpr <- replExprTemplate - replExpr[[2]] <- expr + replExpr[[2]] <- targetExprs[[i]] eval(replExpr) iV_start <- iV_end + 1 } model } + +# From the unevaluated `nodes` argument of values(model, nodes[i], mode = ...), +# evaluate `nodes` and `i` separately in `env`. +is_singly_indexed <- function(nodesExpr) { + is.call(nodesExpr) && identical(nodesExpr[[1]], as.name("[")) && length(nodesExpr) == 3 +} + +split_indexed_nodes <- function(nodesExpr, env, mode) { + if (!is_singly_indexed(nodesExpr)) { + stop("values: mode = \"", mode, "\" requires a singly indexed `nodes` argument, as in nodes[i].") + } + list(nodes = eval(nodesExpr[[2]], env), i = eval(nodesExpr[[3]], env)) +} + +as_varRange_list <- function(nodes) { + if (inherits(nodes, "varRangeClass")) { + return(list(nodes)) + } + if (is.character(nodes)) { + return(lapply(nodes, \(x) varRangeClass$new(x))) + } + nodes +} + +check_values_index <- function(i, n, what) { + if (!(is.numeric(i) && length(i) == 1 && !is.na(i) && i >= 1 && i <= n && i == round(i))) { + stop("values: index ", deparse(i), " is out of range for ", n, " ", what, ".") + } +} + +# In node mode, `nodes` is anything makeInstrList accepts, and node i is in sorted order, +# as for model$calculate_indexed(nodes, i). Elements of `nodes` that are constants are not nodes and +# are omitted. RHS-only elements and logProb variables are errors: they have no instruction to +# align with calculate(nodes[i]) (and logProbs are available from getLogProb(nodes[i])). +# The check uses getNodes, so for large node sets it can be skipped with +# nimbleModelOptions(check_values_node_inputs = FALSE). +check_values_node_input <- function(model, nodes) { + if (!isTRUE(getNimbleModelOption("check_values_node_inputs"))) { + return(invisible(TRUE)) + } + if (!(is.character(nodes) || inherits(nodes, "varRangeClass") || + (is.list(nodes) && length(nodes) && all(vapply(nodes, \(x) inherits(x, "varRangeClass"), logical(1)))))) { + return(invisible(TRUE)) # instructions + } + vr <- as_varRange_list(nodes) + varNames <- vapply(vr, \(x) x$varName, character(1)) + if (any(startsWith(varNames, "logProb_"))) { + stop("values: mode = \"node\" does not support logProb variables (", + paste(unique(varNames[startsWith(varNames, "logProb_")]), collapse = ", "), ").") + } + nNodes <- length(model$getNodes(nodes, returnScalarComponents = TRUE, nodesAsChars = TRUE)) + nWithRHSonly <- length(model$getNodes(nodes, includeRHSonly = TRUE, returnScalarComponents = TRUE, nodesAsChars = TRUE)) + if (nWithRHSonly > nNodes) { + stop("values: mode = \"node\" does not support RHS-only elements of `nodes`.") + } + invisible(TRUE) +} + +# Uncompiled values(model, nodes[i], mode = mode) and its assignment form. +# mode = "raw": entry i of `nodes`. +# mode = "scalar": element i of values(model, nodes), as a length-1 vector. +# mode = "node": node i (see check_values_node_input), flattened. + +#' @export +values_indexed <- function(model, nodes, i, mode = c("raw", "scalar", "node")) { + mode <- match.arg(mode) + if (mode == "node") { + return(model$getValues_indexed(nodes, i)) + } + nodes <- as_varRange_list(nodes) + switch(mode, + raw = { + check_values_index(i, length(nodes), "entries of `nodes`") + values(model, nodes[i]) + }, + scalar = { + vals <- values(model, nodes) + check_values_index(i, length(vals), "elements of values(model, nodes)") + vals[i] + } + ) +} + +#' @export +`values_indexed<-` <- function(model, nodes, i, mode = c("raw", "scalar", "node"), value) { + mode <- match.arg(mode) + if (mode == "node") { + model$setValues_indexed(nodes, i, value) + return(model) + } + nodes <- as_varRange_list(nodes) + switch(mode, + raw = { + check_values_index(i, length(nodes), "entries of `nodes`") + n <- length(values(model, nodes[i])) + if (length(value) != n) { + stop("values: entry ", i, " of `nodes` has ", n, " elements, but ", length(value), " values were given.") + } + values(model, nodes[i]) <- value + }, + scalar = { + if (length(value) != 1) { + stop("values: mode = \"scalar\" requires a single value, but ", length(value), " were given.") + } + lens <- vapply(nodes, \(x) length(values(model, list(x))), numeric(1)) + check_values_index(i, sum(lens), "elements of values(model, nodes)") + k <- findInterval(i - 1, cumsum(c(0, lens)), left.open = FALSE, rightmost.closed = FALSE) + local_i <- i - sum(lens[seq_len(k - 1)]) + vals <- values(model, nodes[k]) + vals[local_i] <- value + values(model, nodes[k]) <- vals + } + ) + model +} diff --git a/nimbleModel/R/declFunBaseClass.R b/nimbleModel/R/declFunBaseClass.R index 0a7ec3a..b729403 100644 --- a/nimbleModel/R/declFunBaseClass.R +++ b/nimbleModel/R/declFunBaseClass.R @@ -3,765 +3,77 @@ declFunBase_nClass <- nClass( classname = "declFunBase_nClass", Rpublic = list( calculate = function(instr) { - calc_op(instr, "calc_one") + run_instr(instr, self$calc_one) }, calculateDiff = function(instr) { - calc_op(instr, "calcDiff_one") + run_instr(instr, self$calcDiff_one) }, getLogProb = function(instr) { - calc_op(instr, "getLogProb_one") + run_instr(instr, self$getLogProb_one) }, getParamID = function(param) { getParamID_one(param) }, - calc_op = function(instr, fn) { - if (instr$instr_type == 0) { - return(calc_0(instr, fn)) - } - if (instr$instr_type == 1) { - return(calc_1_seq(instr, fn)) - } - if (instr$instr_type == 2) { - return(calc_1_mat(instr, fn)) - } - if (instr$instr_type == 3) { - return(calc_1_matp(instr, fn)) - } - if (instr$instr_type == 4) { - return(calc_2_seq_seq(instr, fn)) - } - if (instr$instr_type == 5) { - return(calc_2_seq_mat(instr, fn)) - } - if (instr$instr_type == 6) { - return(calc_2_mat_seq(instr, fn)) - } - if (instr$instr_type == 7) { - return(calc_2_mat_mat(instr, fn)) - } - if (instr$instr_type == 8) { - return(calc_2_seq_matp(instr, fn)) - } - if (instr$instr_type == 9) { - return(calc_2_matp_seq(instr, fn)) - } - if (instr$instr_type == 10) { - return(calc_2_matp_matp(instr, fn)) - } - if (instr$instr_type == 11) { - return(calc_2_x_y_ord(instr, fn)) - } - if (instr$instr_type == 12) { - return(calc_3_allseq(instr, fn)) - } - if (instr$instr_type == 13) { - return(calc_3_generic(instr, fn)) - } - if (instr$instr_type == 14) { - return(calc_4_allseq(instr, fn)) - } - if (instr$instr_type == 15) { - return(calc_4_generic(instr, fn)) - } - if (instr$instr_type == 16) { - return(calc_5_allseq(instr, fn)) - } - if (instr$instr_type == 17) { - return(calc_5_generic(instr, fn)) - } - if (instr$instr_type == 18) { - return(calc_1_matp_ord(instr, fn)) - } - stop("declaration for type ", instr$instr_type, " no implemented") + # sim_one returns a dummy 0, so simulate shares the calc runners and ignores their sum. + simulate = function(instr) { + run_instr(instr, self$sim_one) + invisible(NULL) }, - calc_0 = function(instr, fn) { - return(self[[fn]](0)) + # One step (node) of an instruction, with no loop. `step` is 1-based within `instr`. + calculate_singleton = function(instr, step) { + calc_one(set_instr_idx(instr, step)) }, - calc_1_seq = - function(instr, fn) { - logProb <- 0 - iStart <- instr$values[[1]][1] - for (i in iStart:(iStart + instr$lens[1] - 1)) { - logProb <- logProb + self[[fn]](i) - } - return(logProb) - }, - calc_1_mat = - function(instr, fn) { - logProb <- 0 - for (i in 1:instr$lens[1]) { - logProb <- logProb + self[[fn]](instr$values[[1]][i]) - } - return(logProb) - }, - calc_1_matp = - function(instr, fn) { - logProb <- 0 - for (i in 1:instr$lens[1]) { - logProb <- logProb + self[[fn]](instr$values[[1]][(instr$dims[1] * (i - 1) + 1):(instr$dims[1] * i)]) - } - return(logProb) - }, - calc_1_matp_ord = - function(instr, fn) { - logProb <- 0 - idx <- rep(0, instr$dims[1]) - for (i in 1:instr$lens[1]) { - idx[instr$slots] <- instr$values[[1]][(instr$dims[1] * (i - 1) + 1):(instr$dims[1] * i)] - logProb <- logProb + self[[fn]](idx) - } - return(logProb) - }, - calc_2_seq_seq = - function(instr, fn) { - logProb <- 0 - idx <- rep(0, 2) - iStart1 <- instr$values[[1]][1] - iStart2 <- instr$values[[2]][1] - for (i in iStart1:(iStart1 + instr$lens[1] - 1)) { - idx[1] <- i - for (j in iStart2:(iStart2 + instr$lens[2] - 1)) { - idx[2] <- j - logProb <- logProb + self[[fn]](idx) - } - } - return(logProb) - }, - calc_2_seq_mat = - function(instr, fn) { - logProb <- 0 - idx <- rep(0, 2) - iStart1 <- instr$values[[1]][1] - for (i in iStart1:(iStart1 + instr$lens[1] - 1)) { - idx[1] <- i - for (j in 1:instr$lens[2]) { - idx[2] <- instr$values[[2]][j] - logProb <- logProb + self[[fn]](idx) - } - } - return(logProb) - }, - calc_2_mat_seq = - function(instr, fn) { - logProb <- 0 - idx <- rep(0, 2) - iStart2 <- instr$values[[2]][1] - for (i in 1:instr$lens[1]) { - idx[1] <- instr$values[[1]][i] - for (j in iStart2:(iStart2 + instr$lens[2] - 1)) { - idx[2] <- j - logProb <- logProb + self[[fn]](idx) - } - } - return(logProb) - }, - calc_2_mat_mat = - function(instr, fn) { - logProb <- 0 - idx <- rep(0, 2) - for (i in 1:instr$lens[1]) { - idx[1] <- instr$values[[1]][i] - for (j in 1:instr$lens[2]) { - idx[2] <- instr$values[[2]][j] - logProb <- logProb + self[[fn]](idx) - } - } - return(logProb) - }, - calc_2_seq_matp = - function(instr, fn) { - logProb <- 0 - idx <- rep(0, instr$nDim) - iStart1 <- instr$values[[1]][1] - for (i in iStart1:(iStart1 + instr$lens[1] - 1)) { - idx[instr$slots[1]] <- i - for (j in 1:instr$lens[2]) { - idx[instr$slots[2:instr$nDim]] <- instr$values[[2]][(instr$dims[2] * (j - 1) + 1):(instr$dims[2] * j)] - logProb <- logProb + self[[fn]](idx) - } - } - return(logProb) - }, - calc_2_matp_seq = - function(instr, fn) { - logProb <- 0 - idx <- rep(0, instr$nDim) - iStart2 <- instr$values[[2]][1] - for (i in 1:instr$lens[1]) { - idx[instr$slots[1:instr$dims[1]]] <- instr$values[[1]][(instr$dims[1] * (i - 1) + 1):(instr$dims[1] * i)] - for (j in iStart2:(iStart2 + instr$lens[2] - 1)) { - idx[instr$slots[instr$nDim]] <- j - logProb <- logProb + self[[fn]](idx) - } - } - return(logProb) - }, - calc_2_matp_matp = - function(instr, fn) { - logProb <- 0 - idx <- rep(0, instr$nDim) - for (i in 1:instr$lens[1]) { - idx[instr$slots[1:instr$dims[1]]] <- instr$values[[1]][(instr$dims[1] * (i - 1) + 1):(instr$dims[1] * i)] - for (j in 1:instr$lens[2]) { - idx[instr$slots[(instr$dims[1] + 1):instr$nDim]] <- instr$values[[2]][(instr$dims[2] * (j - 1) + 1):(instr$dims[2] * j)] - logProb <- logProb + self[[fn]](idx) - } - } - return(logProb) - }, - calc_2_x_y_ord = - function(instr, fn) { - logProb <- 0 - idx <- rep(0, 2) - iStart1 <- instr$values[[1]][1] - 1 - iStart2 <- instr$values[[2]][1] - 1 - index_types1 <- instr$index_types[1] - index_types2 <- instr$index_types[2] - for (i in 1:instr$lens[1]) { - if (index_types1 == 1) idx[2] <- iStart1 + i else idx[2] <- instr$values[[1]][i] - for (j in 1:instr$lens[2]) { - if (index_types2 == 1) idx[1] <- iStart2 + j else idx[1] <- instr$values[[2]][j] - logProb <- logProb + self[[fn]](idx) - } - } - return(logProb) - }, - calc_3_allseq = - function(instr, fn) { - logProb <- 0 - idx <- rep(0, 3) - iStart1 <- instr$values[[1]][1] - iStart2 <- instr$values[[2]][1] - iStart3 <- instr$values[[3]][1] - for (i in iStart1:(iStart1 + instr$lens[1] - 1)) { - idx[1] <- i - for (j in iStart2:(iStart2 + instr$lens[2] - 1)) { - idx[2] <- j - for (k in iStart3:(iStart3 + instr$lens[3] - 1)) { - idx[3] <- k - logProb <- logProb + self[[fn]](idx) - } - } - } - return(logProb) - }, - calc_3_generic = - function(instr, fn) { - logProb <- 0 - idx <- rep(0, instr$nDim) - iStart1 <- instr$values[[1]][1] - 1 - iStart2 <- instr$values[[2]][1] - 1 - iStart3 <- instr$values[[3]][1] - 1 - index_types1 <- instr$index_types[1] - index_types2 <- instr$index_types[2] - index_types3 <- instr$index_types[3] - cumdim2 <- instr$dims[1] + instr$dims[2] - for (i in 1:instr$lens[1]) { - if (index_types1 == 1) { - idx[instr$slots[1]] <- iStart1 + i - } else { - idx[instr$slots[1:instr$dims[1]]] <- instr$values[[1]][(instr$dims[1] * (i - 1) + 1):(instr$dims[1] * i)] - } - for (j in 1:instr$lens[2]) { - if (index_types2 == 1) { - idx[instr$slots[instr$dims[1] + 1]] <- iStart2 + j - } else { - idx[instr$slots[(instr$dims[1] + 1):cumdim2]] <- instr$values[[2]][(instr$dims[2] * (j - 1) + 1):(instr$dims[2] * j)] - } - for (k in 1:instr$lens[3]) { - if (index_types3 == 1) { - idx[instr$slots[instr$nDim]] <- iStart3 + k - } else { - idx[instr$slots[(cumdim2 + 1):instr$nDim]] <- instr$values[[3]][(instr$dims[3] * (k - 1) + 1):(instr$dims[3] * k)] - } - logProb <- logProb + self[[fn]](idx) - } - } - } - return(logProb) - }, - calc_4_allseq = - function(instr, fn) { - logProb <- 0 - idx <- rep(0, 4) - iStart <- sapply(instr$values, `[`, 1) - for (i1 in iStart[1]:(iStart[1] + instr$lens[1] - 1)) { - idx[1] <- i1 - for (i2 in iStart[2]:(iStart[2] + instr$lens[2] - 1)) { - idx[2] <- i2 - for (i3 in iStart[3]:(iStart[3] + instr$lens[3] - 1)) { - idx[3] <- i3 - for (i4 in iStart[4]:(iStart[4] + instr$lens[4] - 1)) { - idx[4] <- i4 - logProb <- logProb + self[[fn]](idx) - } - } - } - } - return(logProb) - }, - calc_4_generic = - function(instr, fn) { - logProb <- 0 - idx <- rep(0, instr$nDim) - iStart <- sapply(instr$values, function(x) x[1] - 1) - cumdim2 <- instr$dims[1] + instr$dims[2] - cumdim3 <- instr$dims[1] + instr$dims[2] + instr$dims[3] - for (i1 in 1:instr$lens[1]) { - if (instr$index_types[1] == 1) { - idx[instr$slots[1]] <- iStart[1] + i1 - } else { - idx[instr$slots[1:instr$dims[1]]] <- instr$values[[1]][(instr$dims[1] * (i1 - 1) + 1):(instr$dims[1] * i1)] - } - for (i2 in 1:instr$lens[2]) { - if (instr$index_types[2] == 1) { - idx[instr$slots[instr$dims[1] + 1]] <- iStart[2] + i2 - } else { - idx[instr$slots[(instr$dims[1] + 1):cumdim2]] <- instr$values[[2]][(instr$dims[2] * (i2 - 1) + 1):(instr$dims[2] * i2)] - } - for (i3 in 1:instr$lens[3]) { - if (instr$index_types[3] == 1) { - idx[instr$slots[cumdim2 + 1]] <- iStart[3] + i3 - } else { - idx[instr$slots[(cumdim2 + 1):cumdim3]] <- instr$values[[3]][(instr$dims[3] * (i3 - 1) + 1):(instr$dims[3] * i3)] - } - for (i4 in 1:instr$lens[4]) { - if (instr$index_types[4] == 1) { - idx[instr$slots[cumdim3 + 1]] <- iStart[4] + i4 - } else { - idx[instr$slots[(cumdim3 + 1):instr$nDim]] <- instr$values[[4]][(instr$dims[4] * (i4 - 1) + 1):(instr$dims[4] * i4)] - } - logProb <- logProb + self[[fn]](idx) - } - } - } - } - return(logProb) - }, - calc_5_allseq = - function(instr, fn) { - logProb <- 0 - idx <- rep(0, 4) - iStart <- sapply(instr$values, `[`, 1) - for (i1 in iStart[1]:(iStart[1] + instr$lens[1] - 1)) { - idx[1] <- i1 - for (i2 in iStart[2]:(iStart[2] + instr$lens[2] - 1)) { - idx[2] <- i2 - for (i3 in iStart[3]:(iStart[3] + instr$lens[3] - 1)) { - idx[3] <- i3 - for (i4 in iStart[4]:(iStart[4] + instr$lens[4] - 1)) { - idx[4] <- i4 - for (i5 in iStart[5]:(iStart[5] + instr$lens[5] - 1)) { - idx[5] <- i5 - logProb <- logProb + self[[fn]](idx) - } - } - } - } - } - return(logProb) - }, - calc_5_generic = - function(instr, fn) { - logProb <- 0 - idx <- rep(0, instr$nDim) - iStart <- sapply(instr$values, function(x) x[1] - 1) - cumdim2 <- instr$dims[1] + instr$dims[2] - cumdim3 <- instr$dims[1] + instr$dims[2] + instr$dims[3] - cumdim4 <- instr$dims[1] + instr$dims[2] + instr$dims[3] + instr$dims[4] - for (i1 in 1:instr$lens[1]) { - if (instr$index_types[1] == 1) { - idx[instr$slots[1]] <- iStart[1] + i1 - } else { - idx[instr$slots[1:instr$dims[1]]] <- instr$values[[1]][(instr$dims[1] * (i1 - 1) + 1):(instr$dims[1] * i1)] - } - for (i2 in 1:instr$lens[2]) { - if (instr$index_types[2] == 1) { - idx[instr$slots[instr$dims[1] + 1]] <- iStart[2] + i2 - } else { - idx[instr$slots[(instr$dims[1] + 1):cumdim2]] <- instr$values[[2]][(instr$dims[2] * (i2 - 1) + 1):(instr$dims[2] * i2)] - } - for (i3 in 1:instr$lens[3]) { - if (instr$index_types[3] == 1) { - idx[instr$slots[cumdim2 + 1]] <- iStart[3] + i3 - } else { - idx[instr$slots[(cumdim2 + 1):cumdim3]] <- instr$values[[3]][(instr$dims[3] * (i3 - 1) + 1):(instr$dims[3] * i3)] - } - for (i4 in 1:instr$lens[4]) { - if (instr$index_types[4] == 1) { - idx[instr$slots[cumdim3 + 1]] <- iStart[4] + i4 - } else { - idx[instr$slots[(cumdim3 + 1):instr$nDim]] <- instr$values[[4]][(instr$dims[4] * (i4 - 1) + 1):(instr$dims[4] * i4)] - } - for (i5 in 1:instr$lens[5]) { - if (instr$index_types[5] == 1) { - idx[instr$slots[cumdim4 + 1]] <- iStart[5] + i5 - } else { - idx[instr$slots[(cumdim4 + 1):instr$nDim]] <- instr$values[[5]][(instr$dims[5] * (i5 - 1) + 1):(instr$dims[5] * i5)] - } - logProb <- logProb + self[[fn]](idx) - } - } - } - } - } - return(logProb) - }, - simulate = function(instr) { - if (instr$instr_type == 0) { - return(sim_0(instr)) - } - if (instr$instr_type == 1) { - return(sim_1_seq(instr)) - } - if (instr$instr_type == 2) { - return(sim_1_mat(instr)) - } - if (instr$instr_type == 3) { - return(sim_1_matp(instr)) - } - if (instr$instr_type == 4) { - return(sim_2_seq_seq(instr)) - } - if (instr$instr_type == 5) { - return(sim_2_seq_mat(instr)) - } - if (instr$instr_type == 6) { - return(sim_2_mat_seq(instr)) - } - if (instr$instr_type == 7) { - return(sim_2_mat_mat(instr)) - } - if (instr$instr_type == 8) { - return(sim_2_seq_matp(instr)) - } - if (instr$instr_type == 9) { - return(sim_2_matp_seq(instr)) - } - if (instr$instr_type == 10) { - return(sim_2_matp_matp(instr)) - } - if (instr$instr_type == 11) { - return(sim_2_x_y_ord(instr)) - } - if (instr$instr_type == 12) { - return(sim_3_allseq(instr)) - } - if (instr$instr_type == 13) { - return(sim_3_generic(instr)) - } - if (instr$instr_type == 14) { - return(sim_4_allseq(instr)) - } - if (instr$instr_type == 15) { - return(sim_4_generic(instr)) - } - if (instr$instr_type == 16) { - return(sim_5_allseq(instr)) - } - if (instr$instr_type == 17) { - return(sim_5_generic(instr)) - } - if (instr$instr_type == 18) { - return(sim_1_matp_ord(instr)) - } - stop("declaration for type ", instr$instr_type, " no implemented") + calculateDiff_singleton = function(instr, step) { + calcDiff_one(set_instr_idx(instr, step)) }, - sim_0 = function(instr) { - sim_one(0) + getLogProb_singleton = function(instr, step) { + getLogProb_one(set_instr_idx(instr, step)) }, - sim_1_seq = function(instr) { - iStart <- instr$values[[1]][1] - for (i in iStart:(iStart + instr$lens[1] - 1)) { - sim_one(i) - } + simulate_singleton = function(instr, step) { + sim_one(set_instr_idx(instr, step)) + invisible(NULL) }, - sim_1_mat = function(instr) { - for (i in 1:instr$lens[1]) { - sim_one(instr$values[[1]][i]) - } + getParam_singleton = function(instr, step, paramID) { + getParam_one(set_instr_idx(instr, step), paramID) }, - sim_1_matp = function(instr) { - for (i in 1:instr$lens[1]) { - sim_one(instr$values[[1]][(instr$dims[1] * (i - 1) + 1):(instr$dims[1] * i)]) - } + getBound_singleton = function(instr, step, boundID) { + getBound_one(set_instr_idx(instr, step), boundID) + }, + # The values of one node, flattened, for values(model, nodes[i], mode = "node"). + getValues_singleton = function(instr, step) { + idx <- set_instr_idx(instr, step) + out <- numeric(valuesLength_one(idx)) + getValues_into_one(idx, out, 0) + }, + setValues_singleton = function(instr, step, v) { + idx <- set_instr_idx(instr, step) + if (length(v) != valuesLength_one(idx)) stop("setValues: value has the wrong length for this node.") + setValues_from_one(idx, v, 0) + invisible(NULL) + }, + # All the nodes of an instr, for values(model, nodes, mode = "node"). + valuesLength = function(instr) { + total <- 0 + for (k in seq_len(instr$nNodes)) { + idx <- set_instr_idx(instr, k) # a name, since idx is a refArg + total <- total + valuesLength_one(idx) + } + total + }, + getValues_instr = function(instr) { + unlist(lapply(seq_len(instr$nNodes), \(k) getValues_singleton(instr, k)), use.names = FALSE) + }, + # Set from v starting after offset (0-based), and return the offset after the last value used. + setValues_instr = function(instr, v, offset) { + for (k in seq_len(instr$nNodes)) { + idx <- set_instr_idx(instr, k) + n <- valuesLength_one(idx) + if (offset + n > length(v)) stop("setValues: value is too short for these nodes.") + setValues_from_one(idx, v, offset) + offset <- offset + n + } + offset }, - sim_1_matp_ord = - function(instr) { - logProb <- 0 - idx <- rep(0, instr$dims[1]) - for (i in 1:instr$lens[1]) { - idx[instr$slots] <- instr$values[[1]][(instr$dims[1] * (i - 1) + 1):(instr$dims[1] * i)] - sim_one(idx) - } - }, - sim_2_seq_seq = - function(instr) { - idx <- rep(0, 2) - iStart1 <- instr$values[[1]][1] - 1 - iStart2 <- instr$values[[2]][1] - 1 - for (i in 1:instr$lens[1]) { - idx[1] <- iStart1 + i - for (j in 1:instr$lens[2]) { - idx[2] <- iStart2 + j - sim_one(idx) - } - } - }, - sim_2_seq_mat = - function(instr, fn) { - idx <- rep(0, 2) - iStart1 <- instr$values[[1]][1] - for (i in iStart1:(iStart1 + instr$lens[1] - 1)) { - idx[1] <- i - for (j in 1:instr$lens[2]) { - idx[2] <- instr$values[[2]][j] - sim_one(idx) - } - } - }, - sim_2_mat_seq = - function(instr, fn) { - idx <- rep(0, 2) - iStart2 <- instr$values[[2]][1] - for (i in 1:instr$lens[1]) { - idx[1] <- instr$values[[1]][i] - for (j in iStart2:(iStart2 + instr$lens[2] - 1)) { - idx[2] <- j - sim_one(idx) - } - } - }, - sim_2_mat_mat = - function(instr, fn) { - idx <- rep(0, 2) - for (i in 1:instr$lens[1]) { - idx[1] <- instr$values[[1]][i] - for (j in 1:instr$lens[2]) { - idx[2] <- instr$values[[2]][j] - sim_one(idx) - } - } - }, - sim_2_seq_matp = - function(instr, fn) { - idx <- rep(0, instr$nDim) - iStart1 <- instr$values[[1]][1] - for (i in iStart1:(iStart1 + instr$lens[1] - 1)) { - idx[instr$slots[1]] <- i - for (j in 1:instr$lens[2]) { - idx[instr$slots[2:instr$nDim]] <- instr$values[[2]][(instr$dims[2] * (j - 1) + 1):(instr$dims[2] * j)] - sim_one(idx) - } - } - }, - sim_2_matp_seq = - function(instr, fn) { - idx <- rep(0, instr$nDim) - iStart2 <- instr$values[[2]][1] - for (i in 1:instr$lens[1]) { - idx[instr$slots[1:instr$dims[1]]] <- instr$values[[1]][(instr$dims[1] * (i - 1) + 1):(instr$dims[1] * i)] - for (j in iStart2:(iStart2 + instr$lens[2] - 1)) { - idx[instr$slots[instr$nDim]] <- j - sim_one(idx) - } - } - }, - sim_2_matp_matp = - function(instr, fn) { - idx <- rep(0, instr$nDim) - for (i in 1:instr$lens[1]) { - idx[instr$slots[1:instr$dims[1]]] <- instr$values[[1]][(instr$dims[1] * (i - 1) + 1):(instr$dims[1] * i)] - for (j in 1:instr$lens[2]) { - idx[instr$slots[(instr$dims[1] + 1):instr$nDim]] <- instr$values[[2]][(instr$dims[2] * (j - 1) + 1):(instr$dims[2] * j)] - sim_one(idx) - } - } - }, - sim_2_x_y_ord = - function(instr, fn) { - idx <- rep(0, 2) - iStart1 <- instr$values[[1]][1] - 1 - iStart2 <- instr$values[[2]][1] - 1 - index_types1 <- instr$index_types[1] - index_types2 <- instr$index_types[2] - for (i in 1:instr$lens[1]) { - if (index_types1 == 1) idx[2] <- iStart1 + i else idx[2] <- instr$values[[1]][i] - for (j in 1:instr$lens[2]) { - if (index_types2 == 1) idx[1] <- iStart2 + j else idx[1] <- instr$values[[2]][j] - sim_one(idx) - } - } - }, - sim_3_allseq = - function(instr, fn) { - idx <- rep(0, 3) - iStart1 <- instr$values[[1]][1] - iStart2 <- instr$values[[2]][1] - iStart3 <- instr$values[[3]][1] - for (i in iStart1:(iStart1 + instr$lens[1] - 1)) { - idx[1] <- i - for (j in iStart2:(iStart2 + instr$lens[2] - 1)) { - idx[2] <- j - for (k in iStart3:(iStart3 + instr$lens[3] - 1)) { - idx[3] <- k - sim_one(idx) - } - } - } - }, - sim_3_generic = - function(instr, fn) { - idx <- rep(0, instr$nDim) - iStart1 <- instr$values[[1]][1] - 1 - iStart2 <- instr$values[[2]][1] - 1 - iStart3 <- instr$values[[3]][1] - 1 - index_types1 <- instr$index_types[1] - index_types2 <- instr$index_types[2] - index_types3 <- instr$index_types[3] - cumdim2 <- instr$dims[1] + instr$dims[2] - for (i in 1:instr$lens[1]) { - if (index_types1 == 1) { - idx[instr$slots[1]] <- iStart1 + i - } else { - idx[instr$slots[1:instr$dims[1]]] <- instr$values[[1]][(instr$dims[1] * (i - 1) + 1):(instr$dims[1] * i)] - } - for (j in 1:instr$lens[2]) { - if (index_types2 == 1) { - idx[instr$slots[instr$dims[1] + 1]] <- iStart2 + j - } else { - idx[instr$slots[(instr$dims[1] + 1):cumdim2]] <- instr$values[[2]][(instr$dims[2] * (j - 1) + 1):(instr$dims[2] * j)] - } - for (k in 1:instr$lens[3]) { - if (index_types3 == 1) { - idx[instr$slots[instr$nDim]] <- iStart3 + k - } else { - idx[instr$slots[(cumdim2 + 1):instr$nDim]] <- instr$values[[3]][(instr$dims[3] * (k - 1) + 1):(instr$dims[3] * k)] - } - sim_one(idx) - } - } - } - }, - sim_4_allseq = - function(instr, fn) { - idx <- rep(0, 4) - iStart <- sapply(instr$values, `[`, 1) - for (i1 in iStart[1]:(iStart[1] + instr$lens[1] - 1)) { - idx[1] <- i1 - for (i2 in iStart[2]:(iStart[2] + instr$lens[2] - 1)) { - idx[2] <- i2 - for (i3 in iStart[3]:(iStart[3] + instr$lens[3] - 1)) { - idx[3] <- i3 - for (i4 in iStart[4]:(iStart[4] + instr$lens[4] - 1)) { - idx[4] <- i4 - sim_one(idx) - } - } - } - } - }, - sim_4_generic = - function(instr, fn) { - idx <- rep(0, instr$nDim) - iStart <- sapply(instr$values, function(x) x[1] - 1) - cumdim2 <- instr$dims[1] + instr$dims[2] - cumdim3 <- instr$dims[1] + instr$dims[2] + instr$dims[3] - for (i1 in 1:instr$lens[1]) { - if (instr$index_types[1] == 1) { - idx[instr$slots[1]] <- iStart[1] + i1 - } else { - idx[instr$slots[1:instr$dims[1]]] <- instr$values[[1]][(instr$dims[1] * (i1 - 1) + 1):(instr$dims[1] * i1)] - } - for (i2 in 1:instr$lens[2]) { - if (instr$index_types[2] == 1) { - idx[instr$slots[instr$dims[1] + 1]] <- iStart[2] + i2 - } else { - idx[instr$slots[(instr$dims[1] + 1):cumdim2]] <- instr$values[[2]][(instr$dims[2] * (i2 - 1) + 1):(instr$dims[2] * i2)] - } - for (i3 in 1:instr$lens[3]) { - if (instr$index_types[3] == 1) { - idx[instr$slots[cumdim2 + 1]] <- iStart[3] + i3 - } else { - idx[instr$slots[(cumdim2 + 1):cumdim3]] <- instr$values[[3]][(instr$dims[3] * (i3 - 1) + 1):(instr$dims[3] * i3)] - } - for (i4 in 1:instr$lens[4]) { - if (instr$index_types[4] == 1) { - idx[instr$slots[cumdim3 + 1]] <- iStart[4] + i4 - } else { - idx[instr$slots[(cumdim3 + 1):instr$nDim]] <- instr$values[[4]][(instr$dims[4] * (i4 - 1) + 1):(instr$dims[4] * i4)] - } - sim_one(idx) - } - } - } - } - }, - sim_5_allseq = - function(instr, fn) { - idx <- rep(0, 4) - iStart <- sapply(instr$values, `[`, 1) - for (i1 in iStart[1]:(iStart[1] + instr$lens[1] - 1)) { - idx[1] <- i1 - for (i2 in iStart[2]:(iStart[2] + instr$lens[2] - 1)) { - idx[2] <- i2 - for (i3 in iStart[3]:(iStart[3] + instr$lens[3] - 1)) { - idx[3] <- i3 - for (i4 in iStart[4]:(iStart[4] + instr$lens[4] - 1)) { - idx[4] <- i4 - for (i5 in iStart[5]:(iStart[5] + instr$lens[5] - 1)) { - idx[5] <- i5 - sim_one(idx) - } - } - } - } - } - }, - sim_5_generic = - function(instr, fn) { - idx <- rep(0, instr$nDim) - iStart <- sapply(instr$values, function(x) x[1] - 1) - cumdim2 <- instr$dims[1] + instr$dims[2] - cumdim3 <- instr$dims[1] + instr$dims[2] + instr$dims[3] - cumdim4 <- instr$dims[1] + instr$dims[2] + instr$dims[3] + instr$dims[4] - for (i1 in 1:instr$lens[1]) { - if (instr$index_types[1] == 1) { - idx[instr$slots[1]] <- iStart[1] + i1 - } else { - idx[instr$slots[1:instr$dims[1]]] <- instr$values[[1]][(instr$dims[1] * (i1 - 1) + 1):(instr$dims[1] * i1)] - } - for (i2 in 1:instr$lens[2]) { - if (instr$index_types[2] == 1) { - idx[instr$slots[instr$dims[1] + 1]] <- iStart[2] + i2 - } else { - idx[instr$slots[(instr$dims[1] + 1):cumdim2]] <- instr$values[[2]][(instr$dims[2] * (i2 - 1) + 1):(instr$dims[2] * i2)] - } - for (i3 in 1:instr$lens[3]) { - if (instr$index_types[3] == 1) { - idx[instr$slots[cumdim2 + 1]] <- iStart[3] + i3 - } else { - idx[instr$slots[(cumdim2 + 1):cumdim3]] <- instr$values[[3]][(instr$dims[3] * (i3 - 1) + 1):(instr$dims[3] * i3)] - } - for (i4 in 1:instr$lens[4]) { - if (instr$index_types[4] == 1) { - idx[instr$slots[cumdim3 + 1]] <- iStart[4] + i4 - } else { - idx[instr$slots[(cumdim3 + 1):instr$nDim]] <- instr$values[[4]][(instr$dims[4] * (i4 - 1) + 1):(instr$dims[4] * i4)] - } - for (i5 in 1:instr$lens[5]) { - if (instr$index_types[5] == 1) { - idx[instr$slots[cumdim4 + 1]] <- iStart[5] + i5 - } else { - idx[instr$slots[(cumdim4 + 1):instr$nDim]] <- instr$values[[5]][(instr$dims[5] * (i5 - 1) + 1):(instr$dims[5] * i5)] - } - sim_one(idx) - } - } - } - } - } - }, # getParam: use only the first entry in the instr. getParam = function(instr, paramID) { fi <- first_idx(instr) @@ -987,6 +299,162 @@ declFunBase_nClass <- nClass( return(0) } ) + ), + # One step (node) of an instruction, with no loop; see the uncompiled `*_singleton` methods. + calculate_singleton_cpp = nFunction( + name = "calculate_singleton_cpp", + function(instr, step) { + stop("Uncompiled version of calculate_singleton_cpp should not be called.") + }, + returnType = "numericScalar", + compileInfo = list( + virtual = TRUE, + C_fun = function(instr = "instr_nClass", step = "integerScalar") { + cppLiteral('Rprintf("declFunBase_nClass virtual base calculate_singleton_cpp should never be called (something is wrong)\\n");') + return(0) + } + ) + ), + calculateDiff_singleton_cpp = nFunction( + name = "calculateDiff_singleton_cpp", + function(instr, step) { + stop("Uncompiled version of calculateDiff_singleton_cpp should not be called.") + }, + returnType = "numericScalar", + compileInfo = list( + virtual = TRUE, + C_fun = function(instr = "instr_nClass", step = "integerScalar") { + cppLiteral('Rprintf("declFunBase_nClass virtual base calculateDiff_singleton_cpp should never be called (something is wrong)\\n");') + return(0) + } + ) + ), + getLogProb_singleton_cpp = nFunction( + name = "getLogProb_singleton_cpp", + function(instr, step) { + stop("Uncompiled version of getLogProb_singleton_cpp should not be called.") + }, + returnType = "numericScalar", + compileInfo = list( + virtual = TRUE, + C_fun = function(instr = "instr_nClass", step = "integerScalar") { + cppLiteral('Rprintf("declFunBase_nClass virtual base getLogProb_singleton_cpp should never be called (something is wrong)\\n");') + return(0) + } + ) + ), + simulate_singleton_cpp = nFunction( + name = "simulate_singleton_cpp", + function(instr, step) { + stop("Uncompiled version of simulate_singleton_cpp should not be called.") + }, + returnType = "void", + compileInfo = list( + virtual = TRUE, + C_fun = function(instr = "instr_nClass", step = "integerScalar") { + cppLiteral('Rprintf("declFunBase_nClass virtual base simulate_singleton_cpp should never be called (something is wrong)\\n");') + } + ) + ), + getParam_singleton_cpp = nFunction( + name = "getParam_singleton_cpp", + function(instr, step, paramID) { + stop("Uncompiled version of getParam_singleton_cpp should not be called.") + }, + returnType = "ETaccessor()", + compileInfo = list( + virtual = TRUE, + C_fun = function(instr = "instr_nClass", step = "integerScalar", paramID = "integerScalar") { + cppLiteral('Rprintf("declFunBase_nClass virtual base getParam_singleton_cpp should never be called (something is wrong)\\n");') + return(ETaccess(0, copy = TRUE)) + } + ) + ), + getBound_singleton_cpp = nFunction( + name = "getBound_singleton_cpp", + function(instr, step, boundID) { + stop("Uncompiled version of getBound_singleton_cpp should not be called.") + }, + returnType = "numericScalar", + compileInfo = list( + virtual = TRUE, + C_fun = function(instr = "instr_nClass", step = "integerScalar", boundID = "integerScalar") { + cppLiteral('Rprintf("declFunBase_nClass virtual base getBound_singleton_cpp should never be called (something is wrong)\\n");') + return(0) + } + ) + ), + getValues_singleton_cpp = nFunction( + name = "getValues_singleton_cpp", + function(instr, step) { + stop("Uncompiled version of getValues_singleton_cpp should not be called.") + }, + returnType = "numericVector", + compileInfo = list( + virtual = TRUE, + C_fun = function(instr = "instr_nClass", step = "integerScalar") { + cppLiteral('Rprintf("declFunBase_nClass virtual base getValues_singleton_cpp should never be called (something is wrong)\\n");') + return(nNumeric(length = 0)) + } + ) + ), + setValues_singleton_cpp = nFunction( + name = "setValues_singleton_cpp", + function(instr, step, v) { + stop("Uncompiled version of setValues_singleton_cpp should not be called.") + }, + returnType = "void", + compileInfo = list( + virtual = TRUE, + C_fun = function(instr = "instr_nClass", step = "integerScalar", v = "numericVector") { + cppLiteral('Rprintf("declFunBase_nClass virtual base setValues_singleton_cpp should never be called (something is wrong)\\n");') + } + ) + ), + # All the nodes of an instr, for values(model, nodes, mode = "node"). See declFunClass_.h. + valuesLength_cpp = nFunction( + name = "valuesLength_cpp", + function(instr) { + stop("Uncompiled version of valuesLength_cpp should not be called.") + }, + returnType = "integerScalar", + compileInfo = list( + virtual = TRUE, + C_fun = function(instr = "instr_nClass") { + cppLiteral('Rprintf("declFunBase_nClass virtual base valuesLength_cpp should never be called (something is wrong)\\n");') + return(0L) + } + ) + ), + getValues_into_cpp = nFunction( + name = "getValues_into_cpp", + function(instr, out, offset) { + stop("Uncompiled version of getValues_into_cpp should not be called.") + }, + returnType = "integerScalar", + refArgs = "out", + compileInfo = list( + virtual = TRUE, + C_fun = function(instr = "instr_nClass", out = "numericVector", offset = "integerScalar") { + cppLiteral('Rprintf("declFunBase_nClass virtual base getValues_into_cpp should never be called (something is wrong)\\n");') + return(0L) + } + ) + ), + setValues_from_cpp = nFunction( + name = "setValues_from_cpp", + function(instr, v, offset) { + stop("Uncompiled version of setValues_from_cpp should not be called.") + }, + returnType = "integerScalar", + refArgs = "v", + compileInfo = list( + virtual = TRUE, + C_fun = function(instr = "instr_nClass", v = "numericVector", offset = "integerScalar") { + cppLiteral('Rprintf("declFunBase_nClass virtual base setValues_from_cpp should never be called (something is wrong)\\n");') + return(0L) + } + ) ) ), # We haven't dealt with ensuring a virtual destructor when any method is virtual @@ -997,7 +465,7 @@ declFunBase_nClass <- nClass( interface = "full", createFromR = FALSE, exportName = "declFunBase_nClass_new", - interfaceExclude = "getParam_cpp", + interfaceExclude = c("getParam_cpp", "getParam_singleton_cpp"), needed_units = list("instr_nClass"), packageNames = c(uncompiled = "declFunBase_nClass_R", compiled = "declFunBase_nClass"), depends = "nimbleModel", diff --git a/nimbleModel/R/instrRunners.R b/nimbleModel/R/instrRunners.R new file mode 100644 index 0000000..c76c5eb --- /dev/null +++ b/nimbleModel/R/instrRunners.R @@ -0,0 +1,342 @@ +# Runners: each runs the iterations for one instruction. +# At each iteration, the "f" called is calc_one, sim_one, etc. +run_instr <- function(instr, f) { + switch(as.character(instr$instr_type), + "0" = calc_0(instr, f), + "1" = calc_1_seq(instr, f), + "2" = calc_1_mat(instr, f), + "3" = calc_1_matp(instr, f), + "4" = calc_2_seq_seq(instr, f), + "5" = calc_2_seq_mat(instr, f), + "6" = calc_2_mat_seq(instr, f), + "7" = calc_2_mat_mat(instr, f), + "8" = calc_2_seq_matp(instr, f), + "9" = calc_2_matp_seq(instr, f), + "10" = calc_2_matp_matp(instr, f), + "11" = calc_2_x_y_ord(instr, f), + "12" = calc_3_allseq(instr, f), + "13" = calc_3_generic(instr, f), + "14" = calc_4_allseq(instr, f), + "15" = calc_4_generic(instr, f), + "16" = calc_5_allseq(instr, f), + "17" = calc_5_generic(instr, f), + "18" = calc_1_matp_ord(instr, f), + stop("declaration for type ", instr$instr_type, " no implemented") + ) +} + +calc_0 <- function(instr, f) { + return(f(0)) +} +calc_1_seq <- function(instr, f) { + logProb <- 0 + iStart <- instr$values[[1]][1] + for (i in iStart:(iStart + instr$lens[1] - 1)) { + logProb <- logProb + f(i) + } + return(logProb) +} +calc_1_mat <- function(instr, f) { + logProb <- 0 + for (i in 1:instr$lens[1]) { + logProb <- logProb + f(instr$values[[1]][i]) + } + return(logProb) +} +calc_1_matp <- function(instr, f) { + logProb <- 0 + for (i in 1:instr$lens[1]) { + logProb <- logProb + f(instr$values[[1]][(instr$dims[1] * (i - 1) + 1):(instr$dims[1] * i)]) + } + return(logProb) +} +calc_1_matp_ord <- function(instr, f) { + logProb <- 0 + idx <- rep(0, instr$dims[1]) + for (i in 1:instr$lens[1]) { + idx[instr$slots] <- instr$values[[1]][(instr$dims[1] * (i - 1) + 1):(instr$dims[1] * i)] + logProb <- logProb + f(idx) + } + return(logProb) +} +calc_2_seq_seq <- function(instr, f) { + logProb <- 0 + idx <- rep(0, 2) + iStart1 <- instr$values[[1]][1] + iStart2 <- instr$values[[2]][1] + for (i in iStart1:(iStart1 + instr$lens[1] - 1)) { + idx[1] <- i + for (j in iStart2:(iStart2 + instr$lens[2] - 1)) { + idx[2] <- j + logProb <- logProb + f(idx) + } + } + return(logProb) +} +calc_2_seq_mat <- function(instr, f) { + logProb <- 0 + idx <- rep(0, 2) + iStart1 <- instr$values[[1]][1] + for (i in iStart1:(iStart1 + instr$lens[1] - 1)) { + idx[1] <- i + for (j in 1:instr$lens[2]) { + idx[2] <- instr$values[[2]][j] + logProb <- logProb + f(idx) + } + } + return(logProb) +} +calc_2_mat_seq <- function(instr, f) { + logProb <- 0 + idx <- rep(0, 2) + iStart2 <- instr$values[[2]][1] + for (i in 1:instr$lens[1]) { + idx[1] <- instr$values[[1]][i] + for (j in iStart2:(iStart2 + instr$lens[2] - 1)) { + idx[2] <- j + logProb <- logProb + f(idx) + } + } + return(logProb) +} +calc_2_mat_mat <- function(instr, f) { + logProb <- 0 + idx <- rep(0, 2) + for (i in 1:instr$lens[1]) { + idx[1] <- instr$values[[1]][i] + for (j in 1:instr$lens[2]) { + idx[2] <- instr$values[[2]][j] + logProb <- logProb + f(idx) + } + } + return(logProb) +} +calc_2_seq_matp <- function(instr, f) { + logProb <- 0 + idx <- rep(0, instr$nDim) + iStart1 <- instr$values[[1]][1] + for (i in iStart1:(iStart1 + instr$lens[1] - 1)) { + idx[instr$slots[1]] <- i + for (j in 1:instr$lens[2]) { + idx[instr$slots[2:instr$nDim]] <- instr$values[[2]][(instr$dims[2] * (j - 1) + 1):(instr$dims[2] * j)] + logProb <- logProb + f(idx) + } + } + return(logProb) +} +calc_2_matp_seq <- function(instr, f) { + logProb <- 0 + idx <- rep(0, instr$nDim) + iStart2 <- instr$values[[2]][1] + for (i in 1:instr$lens[1]) { + idx[instr$slots[1:instr$dims[1]]] <- instr$values[[1]][(instr$dims[1] * (i - 1) + 1):(instr$dims[1] * i)] + for (j in iStart2:(iStart2 + instr$lens[2] - 1)) { + idx[instr$slots[instr$nDim]] <- j + logProb <- logProb + f(idx) + } + } + return(logProb) +} +calc_2_matp_matp <- function(instr, f) { + logProb <- 0 + idx <- rep(0, instr$nDim) + for (i in 1:instr$lens[1]) { + idx[instr$slots[1:instr$dims[1]]] <- instr$values[[1]][(instr$dims[1] * (i - 1) + 1):(instr$dims[1] * i)] + for (j in 1:instr$lens[2]) { + idx[instr$slots[(instr$dims[1] + 1):instr$nDim]] <- instr$values[[2]][(instr$dims[2] * (j - 1) + 1):(instr$dims[2] * j)] + logProb <- logProb + f(idx) + } + } + return(logProb) +} +calc_2_x_y_ord <- function(instr, f) { + logProb <- 0 + idx <- rep(0, 2) + iStart1 <- instr$values[[1]][1] - 1 + iStart2 <- instr$values[[2]][1] - 1 + index_types1 <- instr$index_types[1] + index_types2 <- instr$index_types[2] + for (i in 1:instr$lens[1]) { + if (index_types1 == 1) idx[2] <- iStart1 + i else idx[2] <- instr$values[[1]][i] + for (j in 1:instr$lens[2]) { + if (index_types2 == 1) idx[1] <- iStart2 + j else idx[1] <- instr$values[[2]][j] + logProb <- logProb + f(idx) + } + } + return(logProb) +} +calc_3_allseq <- function(instr, f) { + logProb <- 0 + idx <- rep(0, 3) + iStart1 <- instr$values[[1]][1] + iStart2 <- instr$values[[2]][1] + iStart3 <- instr$values[[3]][1] + for (i in iStart1:(iStart1 + instr$lens[1] - 1)) { + idx[1] <- i + for (j in iStart2:(iStart2 + instr$lens[2] - 1)) { + idx[2] <- j + for (k in iStart3:(iStart3 + instr$lens[3] - 1)) { + idx[3] <- k + logProb <- logProb + f(idx) + } + } + } + return(logProb) +} +calc_3_generic <- function(instr, f) { + logProb <- 0 + idx <- rep(0, instr$nDim) + iStart1 <- instr$values[[1]][1] - 1 + iStart2 <- instr$values[[2]][1] - 1 + iStart3 <- instr$values[[3]][1] - 1 + index_types1 <- instr$index_types[1] + index_types2 <- instr$index_types[2] + index_types3 <- instr$index_types[3] + cumdim2 <- instr$dims[1] + instr$dims[2] + for (i in 1:instr$lens[1]) { + if (index_types1 == 1) { + idx[instr$slots[1]] <- iStart1 + i + } else { + idx[instr$slots[1:instr$dims[1]]] <- instr$values[[1]][(instr$dims[1] * (i - 1) + 1):(instr$dims[1] * i)] + } + for (j in 1:instr$lens[2]) { + if (index_types2 == 1) { + idx[instr$slots[instr$dims[1] + 1]] <- iStart2 + j + } else { + idx[instr$slots[(instr$dims[1] + 1):cumdim2]] <- instr$values[[2]][(instr$dims[2] * (j - 1) + 1):(instr$dims[2] * j)] + } + for (k in 1:instr$lens[3]) { + if (index_types3 == 1) { + idx[instr$slots[instr$nDim]] <- iStart3 + k + } else { + idx[instr$slots[(cumdim2 + 1):instr$nDim]] <- instr$values[[3]][(instr$dims[3] * (k - 1) + 1):(instr$dims[3] * k)] + } + logProb <- logProb + f(idx) + } + } + } + return(logProb) +} +calc_4_allseq <- function(instr, f) { + logProb <- 0 + idx <- rep(0, 4) + iStart <- sapply(instr$values, `[`, 1) + for (i1 in iStart[1]:(iStart[1] + instr$lens[1] - 1)) { + idx[1] <- i1 + for (i2 in iStart[2]:(iStart[2] + instr$lens[2] - 1)) { + idx[2] <- i2 + for (i3 in iStart[3]:(iStart[3] + instr$lens[3] - 1)) { + idx[3] <- i3 + for (i4 in iStart[4]:(iStart[4] + instr$lens[4] - 1)) { + idx[4] <- i4 + logProb <- logProb + f(idx) + } + } + } + } + return(logProb) +} +calc_4_generic <- function(instr, f) { + logProb <- 0 + idx <- rep(0, instr$nDim) + iStart <- sapply(instr$values, function(x) x[1] - 1) + cumdim2 <- instr$dims[1] + instr$dims[2] + cumdim3 <- instr$dims[1] + instr$dims[2] + instr$dims[3] + for (i1 in 1:instr$lens[1]) { + if (instr$index_types[1] == 1) { + idx[instr$slots[1]] <- iStart[1] + i1 + } else { + idx[instr$slots[1:instr$dims[1]]] <- instr$values[[1]][(instr$dims[1] * (i1 - 1) + 1):(instr$dims[1] * i1)] + } + for (i2 in 1:instr$lens[2]) { + if (instr$index_types[2] == 1) { + idx[instr$slots[instr$dims[1] + 1]] <- iStart[2] + i2 + } else { + idx[instr$slots[(instr$dims[1] + 1):cumdim2]] <- instr$values[[2]][(instr$dims[2] * (i2 - 1) + 1):(instr$dims[2] * i2)] + } + for (i3 in 1:instr$lens[3]) { + if (instr$index_types[3] == 1) { + idx[instr$slots[cumdim2 + 1]] <- iStart[3] + i3 + } else { + idx[instr$slots[(cumdim2 + 1):cumdim3]] <- instr$values[[3]][(instr$dims[3] * (i3 - 1) + 1):(instr$dims[3] * i3)] + } + for (i4 in 1:instr$lens[4]) { + if (instr$index_types[4] == 1) { + idx[instr$slots[cumdim3 + 1]] <- iStart[4] + i4 + } else { + idx[instr$slots[(cumdim3 + 1):instr$nDim]] <- instr$values[[4]][(instr$dims[4] * (i4 - 1) + 1):(instr$dims[4] * i4)] + } + logProb <- logProb + f(idx) + } + } + } + } + return(logProb) +} +calc_5_allseq <- function(instr, f) { + logProb <- 0 + idx <- rep(0, 4) + iStart <- sapply(instr$values, `[`, 1) + for (i1 in iStart[1]:(iStart[1] + instr$lens[1] - 1)) { + idx[1] <- i1 + for (i2 in iStart[2]:(iStart[2] + instr$lens[2] - 1)) { + idx[2] <- i2 + for (i3 in iStart[3]:(iStart[3] + instr$lens[3] - 1)) { + idx[3] <- i3 + for (i4 in iStart[4]:(iStart[4] + instr$lens[4] - 1)) { + idx[4] <- i4 + for (i5 in iStart[5]:(iStart[5] + instr$lens[5] - 1)) { + idx[5] <- i5 + logProb <- logProb + f(idx) + } + } + } + } + } + return(logProb) +} +calc_5_generic <- function(instr, f) { + logProb <- 0 + idx <- rep(0, instr$nDim) + iStart <- sapply(instr$values, function(x) x[1] - 1) + cumdim2 <- instr$dims[1] + instr$dims[2] + cumdim3 <- instr$dims[1] + instr$dims[2] + instr$dims[3] + cumdim4 <- instr$dims[1] + instr$dims[2] + instr$dims[3] + instr$dims[4] + for (i1 in 1:instr$lens[1]) { + if (instr$index_types[1] == 1) { + idx[instr$slots[1]] <- iStart[1] + i1 + } else { + idx[instr$slots[1:instr$dims[1]]] <- instr$values[[1]][(instr$dims[1] * (i1 - 1) + 1):(instr$dims[1] * i1)] + } + for (i2 in 1:instr$lens[2]) { + if (instr$index_types[2] == 1) { + idx[instr$slots[instr$dims[1] + 1]] <- iStart[2] + i2 + } else { + idx[instr$slots[(instr$dims[1] + 1):cumdim2]] <- instr$values[[2]][(instr$dims[2] * (i2 - 1) + 1):(instr$dims[2] * i2)] + } + for (i3 in 1:instr$lens[3]) { + if (instr$index_types[3] == 1) { + idx[instr$slots[cumdim2 + 1]] <- iStart[3] + i3 + } else { + idx[instr$slots[(cumdim2 + 1):cumdim3]] <- instr$values[[3]][(instr$dims[3] * (i3 - 1) + 1):(instr$dims[3] * i3)] + } + for (i4 in 1:instr$lens[4]) { + if (instr$index_types[4] == 1) { + idx[instr$slots[cumdim3 + 1]] <- iStart[4] + i4 + } else { + idx[instr$slots[(cumdim3 + 1):cumdim4]] <- instr$values[[4]][(instr$dims[4] * (i4 - 1) + 1):(instr$dims[4] * i4)] + } + for (i5 in 1:instr$lens[5]) { + if (instr$index_types[5] == 1) { + idx[instr$slots[cumdim4 + 1]] <- iStart[5] + i5 + } else { + idx[instr$slots[(cumdim4 + 1):instr$nDim]] <- instr$values[[5]][(instr$dims[5] * (i5 - 1) + 1):(instr$dims[5] * i5)] + } + logProb <- logProb + f(idx) + } + } + } + } + } + return(logProb) +} diff --git a/nimbleModel/R/instructions.R b/nimbleModel/R/instructions.R index d5850ad..4b56863 100644 --- a/nimbleModel/R/instructions.R +++ b/nimbleModel/R/instructions.R @@ -52,13 +52,33 @@ range2instr <- function(range) { ) }) # in calcRange, column major; need row major here for simpler/more efficient determination of indices } - instr$instr_type <- determineInstrType(instr) + instr <- add_instr_precalcs(instr) instr$sortID <- range$sortID instr$declID <- range$declID class(instr) <- "Rinstr" return(instr) } +# Add precalculated (once) fields to an instruction: +# nNodes: number of nodes in the instruction. +# strides: these are strides for the *indexing* (nested as written in the *model code*), not the *data*. E.g. in for(i in 1:3) for(j in 1:6) {...}: j goes fastest so has stride=1; i has stride=6. +# idxSlotStarts: number of idx slots filled by the ranges before each range, i.e. where its slots start. +# instr_type: the iteration plan, from `determineInstrType`. +# Every place that builds instructions calls this, so they can't drift apart. +add_instr_precalcs <- function(instr) { + lens <- as.numeric(instr$lens) + nNodes <- prod(lens) # computed as double so the overflow check below can't itself overflow + if (nNodes > .Machine$integer.max) { + stop("an instruction has ", nNodes, " nodes, more than the maximum of ", .Machine$integer.max) + } + instr$nNodes <- nNodes + nRanges <- length(lens) + instr$strides <- rev(cumprod(c(1, rev(lens)[-nRanges]))) + instr$idxSlotStarts <- if (length(instr$dims)) c(0, cumsum(instr$dims)[-length(instr$dims)]) else numeric(0) + instr$instr_type <- determineInstrType(instr) + return(instr) +} + # Eventually think about reordering order of looping for efficiency (and take parallelization into account). # For the moment, we determine mat vs. seq here and then in declFunClass calculate we will determine whether to # vectorize based on whether possible based on the declaration. @@ -130,29 +150,30 @@ makeInstrList <- function(model, input, includeData = TRUE, use_vec = FALSE) { # (1) a char vector of "nodes" # (2) a list of (or single) varRanges # (3) an nList of (or single) instr_nClass objects (assumed to be in sort order) - # (4) an R list of Rinstr objects (created by `range2instr`) (not assumed to be in sort order) + # (4) an R list of (or single) Rinstr objects (created by `range2instr`) (not assumed to be in sort order) + # (5) an `RinstrList` or `instrList_nClass` (returned as is) + # It returns an `RinstrList` (see `make_RinstrList`), or NULL if there are no nodes. + # Already a list of instructions with offsets. + if (inherits(input, "RinstrList") || inherits(input, "instrList_nClass")) { + return(input) + } # A single instruction. - if (inherits(input, "instr_nClass")) { - return(list(input)) + if (inherits(input, "instr_nClass") || inherits(input, "Rinstr")) { + return(make_RinstrList(list(input))) } # An nList of instructions. if (inherits(input, "nList")) { if (!inherits(input[[1]], "instr_nClass")) { stop("nList input to `makeInstrList` should contain `instr_nClass` objects") } else { - return(input) + return(make_RinstrList(input)) } } - # An R list of instr_nClass-like R lists - if (inherits(input, "Rlist_Rinstr")) { - return(input) - } - # An R list of `Rinstr` elements, not assumed to be in sort order. - if (inherits(input, "Rlist_Rinstr") || (is.list(input) && all(sapply(input, \(x) inherits(x, "Rinstr"))))) { - # Create sort-ordered `Rlist_Rinstr`. + if (is.list(input) && all(sapply(input, \(x) inherits(x, "Rinstr")))) { + # Sort by sortID. sortIDs <- lapply(input, \(x) x$sortID) sortIDranges <- sapply(sortIDs, \(x) range(x, na.rm = TRUE)) multiSortID <- which(sortIDranges[1, ] != sortIDranges[2, ]) @@ -167,8 +188,7 @@ makeInstrList <- function(model, input, includeData = TRUE, use_vec = FALSE) { } ord <- order(sortIDranges[1, ]) input <- input[ord] - class(input) <- "Rlist_Rinstr" - return(input) + return(make_RinstrList(input)) } # Finally handle character vectors or varRanges. @@ -202,6 +222,8 @@ makeInstrList <- function(model, input, includeData = TRUE, use_vec = FALSE) { if (!all(diff(sortIDs[[i]]) == 1, na.rm = TRUE)) { # This quickly creates a list of R lists, where the elements mimic instr_nClass objects, # from a calcRange with multiple sortID values. Creating many calcRanges or instr_nClass objects is slow. + # Note: `makeScalarInstrInfoLists` builds instructions with nDim = 1 (one idx slot), so it looks like it + # will only work for declarations with a single loop index (checked by `check_instr_nDims` below). newInstrs <- c(newInstrs, ranges[[i]]$makeScalarInstrInfoLists()) rangesToRemove <- c(rangesToRemove, i) } else { # Check for any overlapping sortID values for the sequential backward dependence calcRange. @@ -219,8 +241,88 @@ makeInstrList <- function(model, input, includeData = TRUE, use_vec = FALSE) { Rlist <- c(newInstrs, lapply(ranges, \(x) range2instr(x))) sortIDs <- sapply(Rlist, \(x) min(x$sortID, na.rm = TRUE)) # `min` still needed for case of ascending sortIDs (e.g., dependence on the past), which are not split. Rlist <- Rlist[order(sortIDs)] - class(Rlist) <- "Rlist_Rinstr" # For checking idempotency. - return(Rlist) + check_instr_nDims(model, Rlist) + return(make_RinstrList(Rlist)) +} + +# Check that each instruction fills as many idx slots as its declaration has loop indices. +# That is the length of `idx` the declaration's per-node methods expect, and in compiled code +# it determines which runners exist for the declaration (see declFunClass_.h). +# The loop indices come from `loopIndexingRule`, not `context`: a lifted node is given the same +# for-loop context as its original node, even if it doesn't use all of those indices +# (see the test "lifted node that uses only some of its loop indices" in test-nimbleModel.R). +check_instr_nDims <- function(model, instrs) { + declInfo <- model$modelDef$declInfo + for (instr in instrs) { + nIdx <- length(declInfo[[instr$declID]]$declRule$loopIndexingRule$indexVarNames) + if (instr$nDim != nIdx) { + stop( + "an instruction for declaration ", instr$declID, " fills ", instr$nDim, + " index slots, but the declaration has ", nIdx, " loop indices" + ) + } + } +} + +# Wrap a sort-ordered list of instructions (Rinstr lists or instr_nClass objects, or an nList of them) +# in the R-list form of `instrList_nClass`. Uncompiled code uses this directly; compiled code copies it +# into an `instrList_nClass` (see `makeCompiledInstrList`). The steps of the whole list are numbered in +# list order, so instruction j covers list steps offsets[j] + 1 through offsets[j] + nNodes of instruction j. +make_RinstrList <- function(instrs) { + n <- length(instrs) + nNodes <- vapply(seq_len(n), \(j) as.numeric(instrs[[j]]$nNodes), numeric(1)) + nNodes_total <- sum(nNodes) # double, so the overflow check below can't itself overflow + if (nNodes_total > .Machine$integer.max) { + stop("a list of instructions has ", nNodes_total, " nodes, more than the maximum of ", .Machine$integer.max) + } + result <- list( + instr_nList = instrs, + offsets = if (n) c(0, cumsum(nNodes)[-n]) else numeric(0), + nNodes_total = nNodes_total + ) + class(result) <- "RinstrList" + return(result) +} + +# Fill `idx` for one step (node) of an instruction, with no loop. `step` is 1-based. +# The counter of range r is ((step - 1) %/% strides[r]) %% lens[r] + 1, and each counter +# is turned into idx values as the runners do: start + counter - 1 for a seq range, or the +# counter's dims[r] values for a mat/matp range, written through that range's slots. +set_instr_idx <- function(instr, step) { + if (!instr$nDim) { + return(0) # as passed by calc_0 + } + idx <- numeric(instr$nDim) + k <- step - 1 + for (r in seq_along(instr$lens)) { + counter <- (k %/% instr$strides[r]) %% instr$lens[r] + 1 + d <- instr$dims[r] + slots <- instr$slots[instr$idxSlotStarts[r] + seq_len(d)] + if (instr$index_types[r] == 1) { + idx[slots] <- instr$values[[r]][1] + counter - 1 + } else { + idx[slots] <- instr$values[[r]][(counter - 1) * d + seq_len(d)] + } + } + idx +} + +# Find which instruction of an instruction list holds list step `step` (1-based), and the step +# within that instruction. Instruction j holds list steps offsets[j] + 1 through offsets[j] + nNodes, +# so j is the last instruction with offsets[j] < step (which skips any instructions with no nodes). +instrList_find_singleton <- function(instrList, step) { + if (step < 1 || step > instrList$nNodes_total) { + stop("step ", step, " is outside the instruction list's ", instrList$nNodes_total, " steps") + } + offsets <- instrList$offsets + j <- findInterval(step - 1, offsets) + list(j = j, k = step - offsets[j]) +} + +# Fill `idx` for list step `step` (1-based) of an instruction list. +set_instr_idx_from_instrList <- function(instrList, step) { + loc <- instrList_find_singleton(instrList, step) + set_instr_idx(instrList$instr_nList[[loc$j]], loc$k) } @@ -248,6 +350,9 @@ instr_nClass <- nClass( self$instr_type <- instr$instr_type %||% 0L # Use integer for compilation (would char be ok?). self$sortID <- instr$sortID %||% integer() self$declID <- instr$declID %||% 0L + self$nNodes <- instr$nNodes %||% 0L + self$strides <- instr$strides %||% integer() + self$idxSlotStarts <- instr$idxSlotStarts %||% integer() } } ), @@ -261,6 +366,9 @@ instr_nClass <- nClass( instr_type = "integerScalar", sortID = "integerVector", declID = "integerScalar", + nNodes = "integerScalar", + strides = "integerVector", + idxSlotStarts = "integerVector", instr_nClass = nFunction( function() { values <- nList(integerVector())$new() @@ -277,3 +385,28 @@ instr_nClass <- nClass( packageNames = c(uncompiled = "instr_nClass_R", compiled = "instr_nClass") ) ) + +# An ordered list of instructions, with the offset of each instruction's first step +# within the whole list. Filled from an `RinstrList` (see `make_RinstrList`). +instrList_nClass <- nClass( + classname = "instrList_nClass", + Cpublic = list( + instr_nList = "nList(instr_nClass)", + offsets = "integerVector", + nNodes_total = "integerScalar", + instrList_nClass = nFunction( + function() { + instr_nList <- nList(instr_nClass)$new() + }, + compileInfo = list(constructor = TRUE) + ) + ), + predefined = quote(system.file(file.path("include", "nimbleModel", "predef"), package = "nimbleModel") |> file.path("instrList_nClass")), + compileInfo = list( + interface = "full", + createFromR = TRUE, + exportName = "instrList_nClass_new", + needed_units = list("instr_nClass", "nList(instr_nClass)"), + packageNames = c(uncompiled = "instrList_nClass_R", compiled = "instrList_nClass") + ) +) diff --git a/nimbleModel/R/loopIndexingRules.R b/nimbleModel/R/loopIndexingRules.R index 6720828..f643d2f 100644 --- a/nimbleModel/R/loopIndexingRules.R +++ b/nimbleModel/R/loopIndexingRules.R @@ -14,6 +14,11 @@ loopIndexingRuleClass <- R6Class( internalRule = NULL, decl = NULL, varName = character(), + # Loop indices used in the LHS, in context order. These are the idx slots that instructions for this + # declaration fill, so they determine the length of `idx` and which `idx[k]` each index becomes + # (see the test "lifted node that uses only some of its loop indices" in test-nimbleModel.R). + # (This is a regression test of a bug that was fixed 2026-10-29.) + indexVarNames = character(), initialize = function(LHS, context, constants = list(), @@ -24,6 +29,7 @@ loopIndexingRuleClass <- R6Class( # Exclude indices not used in lifted expression, e.g., `i` in `y[i,j] ~ dnorm(mu[i], var = sigma2[j])` indexVarNames <- context$indexVarNames indexVarNames <- indexVarNames[indexVarNames %in% all.vars(LHS)] + indexVarNames <<- indexVarNames indexing <- if (length(indexVarNames)) { paste0("[", paste(indexVarNames, collapse = ","), "]") } else { diff --git a/nimbleModel/R/modelBaseClass.R b/nimbleModel/R/modelBaseClass.R index acc5ac7..77b9d29 100644 --- a/nimbleModel/R/modelBaseClass.R +++ b/nimbleModel/R/modelBaseClass.R @@ -1,3 +1,63 @@ +update_model_calc_methods_LAT <- function(code, symTab, auxEnv, info) { + # On input, model$calculate(instrList) + # has become MethodCall(model, calculate(instrList)) + # and here the input will be code = calculate(instrList) + # as exprClasses + method <- switch(code$name, + calculate = "calculate_impl", + calculate_indexed = "calculate_indexed_impl", + calculateDiff = "calculateDiff_impl", + calculateDiff_indexed = "calculateDiff_indexed_impl", + simulate = "simulate_impl", + simulate_indexed = "simulate_indexed_impl", + getLogProb = "getLogProb_impl", + getLogProb_indexed = "getLogProb_indexed_impl", + getBound = "getBound_impl", + getBound_indexed = "getBound_indexed_impl", + getParam = "getParam_impl", + getParam_indexed = "getParam_indexed_impl", + getParamID = "getParamID_impl", + getValues_indexed = "getValues_indexed_impl", + setValues_indexed = "setValues_indexed_impl", + getValues_nodes = "getValues_nodes_impl", + setValues_nodes = "setValues_nodes_impl" + ) + returnType <- switch(code$name, + calculate = "numericScalar", + calculate_indexed = "numericScalar", + calculateDiff = "numericScalar", + calculateDiff_indexed = "numericScalar", + simulate = "void", + simulate_indexed = "void", + getLogProb = "numericScalar", + getLogProb_indexed = "numericScalar", + getBound = "numericScalar", + getBound_indexed = "numericScalar", + getParam = "ETaccessor()", + getParam_indexed = "ETaccessor()", + getParamID = "integerScalar", + getValues_indexed = "numericVector", + setValues_indexed = "void", + getValues_nodes = "numericVector", + setValues_nodes = "void" + ) + code$name <- method + code$type <- nCompiler:::type2symbol(!!returnType) + invisible(NULL) +} + +# `boundID` is a numeric boundID, 1 for lower and 2 for upper (as in nimble), which is what compiled code +# must pass, or (uncompiled only) "lower" or "upper". +bound_to_boundID <- function(boundID, caller) { + if (is.numeric(boundID) && length(boundID) == 1 && boundID %in% c(1, 2)) { + return(boundID) + } + if (!(is.character(boundID) && length(boundID) == 1 && boundID %in% c("lower", "upper"))) { + stop(caller, ": 'boundID' must be 1 (lower) or 2 (upper), or \"lower\" or \"upper\"") + } + if (boundID == "lower") 1 else 2 +} + #' @export modelBase_nClass <- nClass( classname = "modelBase_nClass", @@ -495,45 +555,47 @@ modelBase_nClass <- nClass( } instrList <- makeInstrList(self, instr) if (isCompiled()) { - if (inherits(instrList, "Rlist_Rinstr") || !nCompiler:::isCompiled(instrList)) { + if (inherits(instrList, "RinstrList") || !nCompiler:::isCompiled(instrList)) { instrList <- makeCompiledInstrList(instrList) } return(self[[fn_cpp]](instrList)) } + instrs <- instrList$instr_nList logProb <- 0 - for (i in 1:length(instrList)) { - logProb <- logProb + declFunList[[instrList[[i]]$declID]][[fn]](instrList[[i]]) + for (i in 1:length(instrs)) { + logProb <- logProb + declFunList[[instrs[[i]]$declID]][[fn]](instrs[[i]]) } return(logProb) }, - calculate = function(instr) { - logProb <- calc_op(instr, "calculate", "calculate_impl") + calculate = function(instrList) { + logProb <- calc_op(instrList, "calculate", "calculate_impl") return(logProb) }, - calculateDiff = function(instr) { - logProb <- calc_op(instr, "calculateDiff", "calculateDiff_impl") + calculateDiff = function(instrList) { + logProb <- calc_op(instrList, "calculateDiff", "calculateDiff_impl") return(logProb) }, - getLogProb = function(instr) { - logProb <- calc_op(instr, "getLogProb", "getLogProb_impl") + getLogProb = function(instrList) { + logProb <- calc_op(instrList, "getLogProb", "getLogProb_impl") return(logProb) }, - simulate = function(instr, includeData = FALSE) { - if (missing(instr)) { - instr <- getVarNames() + simulate = function(instrList, includeData = FALSE) { + if (missing(instrList)) { + instrList <- getVarNames() } - instrList <- makeInstrList(self, instr, includeData = includeData) + instrList <- makeInstrList(self, instrList, includeData = includeData) if (is.null(instrList)) { return(invisible(NULL)) } if (isCompiled()) { - if (inherits(instrList, "Rlist_Rinstr") || !nCompiler:::isCompiled(instrList)) { + if (inherits(instrList, "RinstrList") || !nCompiler:::isCompiled(instrList)) { instrList <- makeCompiledInstrList(instrList) } self$simulate_impl(instrList) } else { - for (i in 1:length(instrList)) { - declFunList[[instrList[[i]]$declID]]$simulate(instrList[[i]]) + instrs <- instrList$instr_nList + for (i in 1:length(instrs)) { + declFunList[[instrs[[i]]$declID]]$simulate(instrs[[i]]) } } return(invisible(NULL)) @@ -541,62 +603,174 @@ modelBase_nClass <- nClass( # `one_instr` is a single instr_nClass (if an instList is provided, the first one is used). # Within it, one node will be used (see declFunBase_nClass$first_idx). # `param` is an integer paramID or a character param name. - getParamID = function(one_instr, param) { - instrList <- makeInstrList(self, one_instr) + getParamID = function(instr, param) { + instrList <- makeInstrList(self, instr) + first_instr <- instrList$instr_nList[[1]] if (is.numeric(param)) { paramID <- param } else if (is.character(param)) { if (isCompiled()) { - if (!nCompiler:::isCompiled(instrList)) instr <- makeCompiledInstr(instrList[[1]]) - paramID <- self$getParamID_impl(instr, param) + one_instr <- if (nCompiler:::isCompiled(first_instr)) first_instr else makeCompiledInstr(first_instr) + paramID <- self$getParamID_impl(one_instr, param) } else { - paramID <- declFunList[[instrList[[1]]$declID]]$getParamID(param) + paramID <- declFunList[[first_instr$declID]]$getParamID(param) } } else { stop("getParamID: argument `param` must be a character param name or (redundantly) a numeric paramID.") } return(paramID) }, - getParam = function(one_instr, param) { - instrList <- makeInstrList(self, one_instr) + getParam = function(instr, paramID) { + input_paramID <- paramID # an integer paramID or a character param name + instrList <- makeInstrList(self, instr) + first_instr <- instrList$instr_nList[[1]] if (isCompiled()) { - if (!nCompiler:::isCompiled(instrList)) instr <- makeCompiledInstr(instrList[[1]]) + one_instr <- if (nCompiler:::isCompiled(first_instr)) first_instr else makeCompiledInstr(first_instr) } - if (is.numeric(param)) { - paramID <- param - } else if (is.character(param)) { - paramID <- self$getParamID(instrList, param) + if (is.numeric(input_paramID)) { + paramID <- input_paramID + } else if (is.character(input_paramID)) { + paramID <- self$getParamID(instrList, input_paramID) # getParamID( # self$modelDef$declInfo[[one_instr$declID]]$distributionName, # param # ) # FIXME: Upon switching from the general function getParamID to the constructed method, # I am not sure that we return NA as a failure value. - if (is.na(paramID)) stop("getParam: argument `param` is not a valid parameter name for this node.") + if (is.na(paramID)) stop("getParam: argument `paramID` is not a valid parameter name for this node.") } else { - stop("getParam: argument `param` must be an integer paramID or a character param name.") + stop("getParam: argument `paramID` must be an integer paramID or a character param name.") } if (isCompiled()) { - # if (!nCompiler:::isCompiled(instrList)) instr <- makeCompiledInstr(instrList[[1]]) - return(self$getParam_impl(instr, paramID)) + return(self$getParam_impl(one_instr, paramID)) } - declFunList[[instrList[[1]]$declID]]$getParam(instrList[[1]], paramID) + declFunList[[first_instr$declID]]$getParam(first_instr, paramID) }, - getBound = function(one_instr, bound) { - instrList <- makeInstrList(self, one_instr) - if (!bound %in% c("lower", "upper")) { - stop("getBound: 'bound' must be either 'lower' or 'upper'") + getBound = function(instr, boundID) { + instrList <- makeInstrList(self, instr) + first_instr <- instrList$instr_nList[[1]] + boundID <- bound_to_boundID(boundID, "getBound") + if (isCompiled()) { + one_instr <- if (nCompiler:::isCompiled(first_instr)) first_instr else makeCompiledInstr(first_instr) + return(self$getBound_impl(one_instr, boundID)) } - boundID <- if (bound == "lower") 0 else 1 + declFunList[[first_instr$declID]]$getBound(first_instr, boundID) + }, + # One node of an instruction list, with no loop: `step` (1-based) numbers the nodes of the whole list + # in calculation order, i.e. this is `nodes[step]`. `instr` is anything `makeInstrList` accepts. + calc_indexed_op = function(instr, step, fn, fn_cpp, includeData = TRUE) { + instrList <- makeInstrList(self, instr, includeData = includeData) if (isCompiled()) { - if (!nCompiler:::isCompiled(instrList)) { - instr <- makeCompiledInstr(instrList[[1]]) - } else { - instr <- instrList[[1]] + if (inherits(instrList, "RinstrList") || !nCompiler:::isCompiled(instrList)) { + instrList <- makeCompiledInstrList(instrList) + } + return(self[[fn_cpp]](instrList, step)) + } + loc <- instrList_find_singleton(instrList, step) + instr_j <- instrList$instr_nList[[loc$j]] + declFunList[[instr_j$declID]][[fn]](instr_j, loc$k) + }, + calculate_indexed = function(instrList, i) { + calc_indexed_op(instrList, i, "calculate_singleton", "calculate_indexed_impl") + }, + calculateDiff_indexed = function(instrList, i) { + calc_indexed_op(instrList, i, "calculateDiff_singleton", "calculateDiff_indexed_impl") + }, + getLogProb_indexed = function(instrList, i) { + calc_indexed_op(instrList, i, "getLogProb_singleton", "getLogProb_indexed_impl") + }, + # As with `simulate`, data nodes are omitted (and so not numbered) unless `includeData = TRUE`. + simulate_indexed = function(instrList, i, includeData = FALSE) { + calc_indexed_op(instrList, i, "simulate_singleton", "simulate_indexed_impl", includeData = includeData) + return(invisible(NULL)) + }, + getParam_indexed = function(instrList, i, paramID) { + input_paramID <- paramID # an integer paramID or a character param name + instrList <- makeInstrList(self, instrList) + loc <- instrList_find_singleton(instrList, i) + instr_j <- instrList$instr_nList[[loc$j]] + if (is.numeric(input_paramID)) { + paramID <- input_paramID + } else if (is.character(input_paramID)) { + paramID <- self$getParamID(instr_j, input_paramID) + if (is.na(paramID)) stop("getParam_indexed: argument `paramID` is not a valid parameter name for this node.") + } else { + stop("getParam_indexed: argument `paramID` must be an integer paramID or a character param name.") + } + if (isCompiled()) { + if (inherits(instrList, "RinstrList") || !nCompiler:::isCompiled(instrList)) { + instrList <- makeCompiledInstrList(instrList) + } + return(self$getParam_indexed_impl(instrList, i, paramID)) + } + declFunList[[instr_j$declID]]$getParam_singleton(instr_j, loc$k, paramID) + }, + getBound_indexed = function(instrList, i, boundID) { + boundID <- bound_to_boundID(boundID, "getBound_indexed") + instrList <- makeInstrList(self, instrList) + if (isCompiled()) { + if (inherits(instrList, "RinstrList") || !nCompiler:::isCompiled(instrList)) { + instrList <- makeCompiledInstrList(instrList) + } + return(self$getBound_indexed_impl(instrList, i, boundID)) + } + loc <- instrList_find_singleton(instrList, i) + instr_j <- instrList$instr_nList[[loc$j]] + declFunList[[instr_j$declID]]$getBound_singleton(instr_j, loc$k, boundID) + }, + # The flattened values of node i from an instrList. + # This is used by values(model, nodes[i], mode = "node"). + getValues_indexed = function(instrList, i) { + check_values_node_input(self, instrList) + calc_indexed_op(instrList, i, "getValues_singleton", "getValues_indexed_impl") + }, + setValues_indexed = function(instrList, i, value) { + check_values_node_input(self, instrList) + instrList <- makeInstrList(self, instrList) + if (isCompiled()) { + if (inherits(instrList, "RinstrList") || !nCompiler:::isCompiled(instrList)) { + instrList <- makeCompiledInstrList(instrList) + } + self$setValues_indexed_impl(instrList, i, as.numeric(value)) + return(invisible(NULL)) + } + loc <- instrList_find_singleton(instrList, i) + instr_j <- instrList$instr_nList[[loc$j]] + declFunList[[instr_j$declID]]$setValues_singleton(instr_j, loc$k, as.numeric(value)) + invisible(NULL) + }, + # The flattened values of all the nodes of an instrList, in its (sorted) order. + # This is used by values(model, nodes, mode = "node"). + getValues_nodes = function(instrList) { + check_values_node_input(self, instrList) + instrList <- makeInstrList(self, instrList) + if (isCompiled()) { + if (inherits(instrList, "RinstrList") || !nCompiler:::isCompiled(instrList)) { + instrList <- makeCompiledInstrList(instrList) } - return(self$getBound_impl(instr, boundID)) + return(self$getValues_nodes_impl(instrList)) + } + unlist(lapply(instrList$instr_nList, \(instr) declFunList[[instr$declID]]$getValues_instr(instr)), use.names = FALSE) + }, + setValues_nodes = function(instrList, value) { + check_values_node_input(self, instrList) + instrList <- makeInstrList(self, instrList) + value <- as.numeric(value) + if (isCompiled()) { + if (inherits(instrList, "RinstrList") || !nCompiler:::isCompiled(instrList)) { + instrList <- makeCompiledInstrList(instrList) + } + self$setValues_nodes_impl(instrList, value) + return(invisible(NULL)) + } + instrs <- instrList$instr_nList + total <- sum(vapply(instrs, \(instr) declFunList[[instr$declID]]$valuesLength(instr), numeric(1))) + if (length(value) != total) { + stop("setValues: value has length ", length(value), " but the nodes have ", total, " values.") } - declFunList[[instrList[[1]]$declID]]$getBound(instrList[[1]], boundID) + offset <- 0 + for (instr in instrs) offset <- declFunList[[instr$declID]]$setValues_instr(instr, value, offset) + invisible(NULL) } ), Cpublic = list( @@ -611,11 +785,11 @@ modelBase_nClass <- nClass( makeCompiledInstrList = nFunction( name = "makeCompiledInstrList", function(input = "SEXP") { - ans <- nList(instr_nClass)$new() + ans <- instrList_nClass$new() cppLiteral("ans->set_all_values(input);") return(ans) }, - returnType = "nList(instr_nClass)" + returnType = "instrList_nClass" ), # Like `makeCompiledInstrList` but for a single instr_nClass (used by # `getParam`, which takes one node rather than a list of them). @@ -636,7 +810,7 @@ modelBase_nClass <- nClass( }, returnType = "numericScalar", compileInfo = list( - C_fun = function(instrList = "nList(instr_nClass)") { + C_fun = function(instrList = "instrList_nClass") { # NOTE: instrList input will be ordered. cppLiteral('Rprintf("modelBase_nClass calculate_impl (should not see this)\\n");') return(0) @@ -652,7 +826,7 @@ modelBase_nClass <- nClass( }, returnType = "numericScalar", compileInfo = list( - C_fun = function(instrList = "nList(instr_nClass)") { + C_fun = function(instrList = "instrList_nClass") { # NOTE: instrList input will be ordered. cppLiteral('Rprintf("modelBase_nClass calculateDiff_impl (should not see this)\\n");') return(0) @@ -668,7 +842,7 @@ modelBase_nClass <- nClass( }, returnType = "numericScalar", compileInfo = list( - C_fun = function(instrList = "nList(instr_nClass)") { + C_fun = function(instrList = "instrList_nClass") { # NOTE: instrList input will be ordered. cppLiteral('Rprintf("modelBase_nClass getLogProb_impl (should not see this)\\n");') return(0) @@ -684,7 +858,7 @@ modelBase_nClass <- nClass( }, returnType = "void", compileInfo = list( - C_fun = function(instrList = "nList(instr_nClass)") { + C_fun = function(instrList = "instrList_nClass") { # NOTE: instrList input will be ordered. cppLiteral('Rprintf("modelBase_nClass simulate_impl (should not see this)\\n");') }, @@ -736,25 +910,206 @@ modelBase_nClass <- nClass( }, virtual = TRUE ) + ), + # One node (list step) of an instruction list; see `calc_indexed_op`. + calculate_indexed_impl = nFunction( + name = "calculate_indexed_impl", + function(instrList, step) { + cat("Uncompiled `calculate_indexed_impl` should never be called.\n") + return(0) + }, + returnType = "numericScalar", + compileInfo = list( + C_fun = function(instrList = "instrList_nClass", step = "integerScalar") { + cppLiteral('Rprintf("modelBase_nClass calculate_indexed_impl (should not see this)\\n");') + return(0) + }, + virtual = TRUE + ) + ), + calculateDiff_indexed_impl = nFunction( + name = "calculateDiff_indexed_impl", + function(instrList, step) { + cat("Uncompiled `calculateDiff_indexed_impl` should never be called.\n") + return(0) + }, + returnType = "numericScalar", + compileInfo = list( + C_fun = function(instrList = "instrList_nClass", step = "integerScalar") { + cppLiteral('Rprintf("modelBase_nClass calculateDiff_indexed_impl (should not see this)\\n");') + return(0) + }, + virtual = TRUE + ) + ), + getLogProb_indexed_impl = nFunction( + name = "getLogProb_indexed_impl", + function(instrList, step) { + cat("Uncompiled `getLogProb_indexed_impl` should never be called.\n") + return(0) + }, + returnType = "numericScalar", + compileInfo = list( + C_fun = function(instrList = "instrList_nClass", step = "integerScalar") { + cppLiteral('Rprintf("modelBase_nClass getLogProb_indexed_impl (should not see this)\\n");') + return(0) + }, + virtual = TRUE + ) + ), + simulate_indexed_impl = nFunction( + name = "simulate_indexed_impl", + function(instrList, step) { + cat("Uncompiled `simulate_indexed_impl` should never be called.\n") + return(invisible(NULL)) + }, + returnType = "void", + compileInfo = list( + C_fun = function(instrList = "instrList_nClass", step = "integerScalar") { + cppLiteral('Rprintf("modelBase_nClass simulate_indexed_impl (should not see this)\\n");') + }, + virtual = TRUE + ) + ), + getParam_indexed_impl = nFunction( + name = "getParam_indexed_impl", + function(instrList, step, paramID) { + cat("Uncompiled `getParam_indexed_impl` should never be called.\n") + return(NULL) + }, + # We must return a general object because it may have different nDim + returnType = "ETaccessor()", + compileInfo = list( + C_fun = function(instrList = "instrList_nClass", step = "integerScalar", paramID = "integerScalar") { + cppLiteral('Rprintf("modelBase_nClass getParam_indexed_impl (should not see this)\\n"); return nullptr;') + }, + virtual = TRUE + ) + ), + getBound_indexed_impl = nFunction( + name = "getBound_indexed_impl", + function(instrList, step, boundID) { + cat("Uncompiled `getBound_indexed_impl` should never be called.\n") + return(0) + }, + returnType = "numericScalar", + compileInfo = list( + C_fun = function(instrList = "instrList_nClass", step = "integerScalar", boundID = "integerScalar") { + cppLiteral('Rprintf("modelBase_nClass getBound_indexed_impl (should not see this)\\n");') + return(0) + }, + virtual = TRUE + ) + ), + getValues_indexed_impl = nFunction( + name = "getValues_indexed_impl", + function(instrList, step) { + cat("Uncompiled `getValues_indexed_impl` should never be called.\n") + return(0) + }, + returnType = "numericVector", + compileInfo = list( + C_fun = function(instrList = "instrList_nClass", step = "integerScalar") { + cppLiteral('Rprintf("modelBase_nClass getValues_indexed_impl (should not see this)\\n");') + return(nNumeric(length = 0)) + }, + virtual = TRUE + ) + ), + setValues_indexed_impl = nFunction( + name = "setValues_indexed_impl", + function(instrList, step, value) { + cat("Uncompiled `setValues_indexed_impl` should never be called.\n") + }, + returnType = "void", + compileInfo = list( + C_fun = function(instrList = "instrList_nClass", step = "integerScalar", value = "numericVector") { + cppLiteral('Rprintf("modelBase_nClass setValues_indexed_impl (should not see this)\\n");') + }, + virtual = TRUE + ) + ), + getValues_nodes_impl = nFunction( + name = "getValues_nodes_impl", + function(instrList) { + cat("Uncompiled `getValues_nodes_impl` should never be called.\n") + return(0) + }, + returnType = "numericVector", + compileInfo = list( + C_fun = function(instrList = "instrList_nClass") { + cppLiteral('Rprintf("modelBase_nClass getValues_nodes_impl (should not see this)\\n");') + return(nNumeric(length = 0)) + }, + virtual = TRUE + ) + ), + setValues_nodes_impl = nFunction( + name = "setValues_nodes_impl", + function(instrList, value) { + cat("Uncompiled `setValues_nodes_impl` should never be called.\n") + }, + returnType = "void", + compileInfo = list( + C_fun = function(instrList = "instrList_nClass", value = "numericVector") { + cppLiteral('Rprintf("modelBase_nClass setValues_nodes_impl (should not see this)\\n");') + }, + virtual = TRUE + ) ) ), # See comment above about needing to ensure a virtual destructor predefined = quote(system.file(file.path("include", "nimbleModel", "predef"), package = "nimbleModel") |> file.path("modelBase_nC")), compileInfo = list( + # Each of these model methods, called from compiled code as model$method(...), becomes the + # corresponding *_impl method (see update_model_calc_methods_LAT). + opDefs = sapply( + c( + "calculate", "calculate_indexed", "calculateDiff", "calculateDiff_indexed", + "simulate", "simulate_indexed", "getLogProb", "getLogProb_indexed", + "getBound", "getBound_indexed", "getParam", "getParam_indexed", "getParamID", + "getValues_indexed", "setValues_indexed", "getValues_nodes", "setValues_nodes" + ), + function(method) { + matchDef <- if(method %in% c("calculate", "calculateDiff", "simulate", "getLogProb")) { + function(instrList) {} + } else if (method %in% c("calculate_indexed", "calculateDiff_indexed", "simulate_indexed", "getLogProb_indexed")) { + function(instrList, i) {} + } else if (method == "getParam") { + function(instr, paramID) {} + } else if (method == "getBound") { + function(instr, boundID) {} + } else if (method == "getParam_indexed") { + function(instrList, i, paramID) {} + } else if (method == "getBound_indexed") { + function(instrList, i, boundID) {} + } else if (method == "getParamID") { + function(instr, param) {} + } else if (method == "getValues_indexed") { + function(instrList, i) {} + } else if (method == "setValues_indexed") { + function(instrList, i, value) {} + } else if (method == "getValues_nodes") { + function(instrList) {} + } else if (method == "setValues_nodes") { + function(instrList, value) {} + } + # matchDef belongs at the top level of an opDef (see nCompiler:::normalizeCallArgs). For model$method(...), + # labelAbstractTypesEnv$MethodCall uses it to put the arguments in order before calling the handler. + list( + matchDef = matchDef, + labelAbstractTypes = list(handler = update_model_calc_methods_LAT) + ) + }, + simplify = FALSE + ), interface = "full", createFromR = FALSE, - Hincludes = c('"declFunBase_nClass_c_.h"', '"instr_nClass_c_.h"'), - needed_units = list("declFunBase_nClass", "instr_nClass", "nList(instr_nClass)"), + Hincludes = c('"declFunBase_nClass_c_.h"', '"instr_nClass_c_.h"', '"instrList_nClass_c_.h"'), + needed_units = list("declFunBase_nClass", "instr_nClass", "nList(instr_nClass)", "instrList_nClass"), exportName = "modelBase_nClass_new", packageNames = c(uncompiled = "modelBase_nClass_R", compiled = "modelBase_nClass"), depends = "nimbleModel", HpostCode = "#include " ) ) - -# Manually add -# # "#include " to that file, -# after the header content. - - -# nCompile(modelBase_nClass, control=list(generate_predefined=TRUE)) diff --git a/nimbleModel/R/nimbleModel.R b/nimbleModel/R/nimbleModel.R index 90347f0..c541bb6 100644 --- a/nimbleModel/R/nimbleModel.R +++ b/nimbleModel/R/nimbleModel.R @@ -89,6 +89,9 @@ make_modelClass_from_nimbleModel <- function(modelDef, data, inits, name = NULL) unlist() |> unique() |> setdiff(c("idx", "paramID_", "boundID_", "LocalNewLogProb_", "LocalAns_", "model", "param_name_")) %||% character() + # Other method locals (e.g. in getValues_into_one) are named Local*_, and + # nFunction's handling of refArgs adds *_Ref__ names to the uncompiled bodies. + declVars <- declVars[!grepl("^Local[[:alnum:]]*_$", declVars) & !grepl("_Ref__$", declVars)] declVarInfo <- modelVarInfo$vars[declVars] declID <- as.numeric(declInfo$declRule$ID) # Formerly `sourceLineNumber`, which may not be unique. declFun_membername <- declFunNames[i] @@ -96,7 +99,10 @@ make_modelClass_from_nimbleModel <- function(modelDef, data, inits, name = NULL) declFun_RvarName <- sub("declFun", "declFunClassGen", declFun_membername) # name of R var holding the nClass generator # Currently, we can't just make a list of these but need them as named objects in the environment, # which is passed into the nClass() call so that `initialize()` can use them via R's scoping. - assign(declFun_RvarName, make_declFun_nClass(declVarInfo, decl_methods, declFun_classname, declID)) + # Length of `idx` for this declaration. Not `context$singleContexts`, which for a lifted node + # can include loop indices the node doesn't use (see the test "lifted node that uses only some of its loop indices" in test-nimbleModel.R). + nIdx <- length(declInfo$declRule$loopIndexingRule$indexVarNames) + assign(declFun_RvarName, make_declFun_nClass(declVarInfo, decl_methods, declFun_classname, declID, nIdx)) declInfoList[[i]] <- make_decl_info_for_model_nClass(declFun_membername, declFun_RvarName, declFun_classname, declVarInfo) } # We have a canonical ordering of decls, but it does arise from a couple of places that should match. @@ -150,7 +156,8 @@ nm_addModelDollarSign <- function(expr, exceptionNames = character(0)) { make_declFun_nClass <- function(varInfo = list(), methods = list(), classname, - declID) { + declID, + nIdx) { # varInfo will be a list (names not used) of name, nDim, sizes. # These are the model member variables to be used by the declFxn. # They will be used in a constructor to set up C++ references to model variables. @@ -190,7 +197,8 @@ make_declFun_nClass <- function(varInfo = list(), classname <- declLabelCreator() } - baseclass <- paste0("declFunClass_<", classname, ">") + # nIdx lets the C++ base class compile only the runners this declaration's instructions can use. + baseclass <- paste0("declFunClass_<", classname, ", ", nIdx, ">") # Rpublic method to set the model pointer/reference. setModel <- function(model) { @@ -528,19 +536,19 @@ make_determ_calc_line <- function(LHSrep, RHSrep) { calc_line } -make_nFxn_from_Cfun <- function(Cfun) { - Rfun <- Cfun - body(calc1Rfun) <- nm_addModelDollarSign(body(Cfun), exceptionNames = c("idx")) - nFxn <- nFunction( - name = "calc_one", - fun = Rfun, - compileInfo = list(C_fun = Cfun), - argTypes = list(idx = "integerVector"), - returnType = "numericScalar" - ) - # declVars <- all.vars(body(calc1Cfun)) |> setdiff("idx") - nFxn -} +# make_nFxn_from_Cfun <- function(Cfun) { +# Rfun <- Cfun +# body(calc1Rfun) <- nm_addModelDollarSign(body(Cfun), exceptionNames = c("idx")) +# nFxn <- nFunction( +# name = "calc_one", +# fun = Rfun, +# compileInfo = list(C_fun = Cfun), +# argTypes = list(idx = "integerVector"), +# returnType = "numericScalar" +# ) +# # declVars <- all.vars(body(calc1Cfun)) |> setdiff("idx") +# nFxn +# } make_decl_method_nFxn <- function(f, name, returnType = "numericScalar") { Cfun <- f @@ -567,22 +575,28 @@ make_decl_methods_from_declInfo <- function(declInfo) { LHS <- modelCode[[2]] RHS <- modelCode[[3]] type <- if (modelCode[[1]] == "~") "stoch" else "determ" # or use declInfo$stoch (logical) - context <- declInfo$declRule$context + # Use the loop indices that instructions fill (see loopIndexingRuleClass), not all of `context`, + # which for a lifted node can include loop indices the node doesn't use + # (see the test "lifted node that uses only some of its loop indices" in test-nimbleModel.R). + # (This is a regression test of a bug fixed on 2026-10-29.) + indexVarNames <- declInfo$declRule$loopIndexingRule$indexVarNames replacements <- sapply( - seq_along(context$singleContexts), + seq_along(indexVarNames), function(i) parse(text = paste0("idx[", i, "]"))[[1]] ) - names(replacements) <- context$indexVarNames + names(replacements) <- indexVarNames LHSrep <- eval(substitute(substitute(e, replacements), list(e = LHS))) RHSrep <- eval(substitute(substitute(e, replacements), list(e = RHS))) if (type == "determ") { methodList <- eval(substitute( list( + # sim_one returns a dummy 0 so that simulate can share the calc runners. sim_one = (function(idx) { calc_one(idx) + return(0) }) |> - make_decl_method_nFxn("sim_one", NULL), + make_decl_method_nFxn("sim_one"), calc_one = (function(idx) { DETERMCALC return(0) @@ -612,8 +626,9 @@ make_decl_methods_from_declInfo <- function(declInfo) { list( sim_one = (function(idx) { STOCHSIM + return(0) }) |> - make_decl_method_nFxn("sim_one", NULL), + make_decl_method_nFxn("sim_one"), calc_one = (function(idx) { STOCHCALC return(LOGPROB) @@ -660,7 +675,147 @@ make_decl_methods_from_declInfo <- function(declInfo) { getBoundMethods <- make_getBound_methods(declInfo$distributionName, LHSrep, RHSrep, boundExprsRep) getParamIDMethod <- make_getParamID_method(declInfo$distributionName) } - c(methodList, getParamMethods, getBoundMethods, getParamIDMethod) + valuesMethods <- make_values_one_methods(LHSrep) + c(methodList, getParamMethods, getBoundMethods, getParamIDMethod, valuesMethods) +} + +# Generated per-node methods for values(model, nodes[i], mode = "node") and +# values(model, nodes, mode = "node"). For one node (its LHS, as LHSrep with idx filled in): +# valuesLength_one(idx): the number of values (as a double, since that is the type of +# index arithmetic in the DSL). +# getValues_into_one(idx, LocalOut_, LocalOff_): copy the values, flattened in column-major +# order, into LocalOut_ starting after LocalOff_ (a 0-based offset). +# setValues_from_one(idx, LocalV_, LocalOff_): the reverse. +# These do no checking. The loops over nodes, the offsets and the length checks are hand-coded +# in declFunClass_.h, which calls these without virtual dispatch. +# The vector arguments (and idx) are passed by reference, so a loop over nodes doesn't copy them. +# +# The compiled code depends on D, the number of ranges in the LHS: +# D = 0 (scalar) and D = 1 use the LHS directly in DSL code. +# D >= 2 needs a flattened view. In that case, we put +# the index expressions into local variables and call +# nm_getBlock_into_ or nm_setBlock_from_ +# (in declFunClass_.h), which use a RuntimeFlatView. +# All locals are named Local*_, which is used above as a pattern +# to exclude from the model$ prefixing in nm_addModelDollarSign. +make_values_one_methods <- function(LHSrep) { + if (is.name(LHSrep)) { + # A scalar variable. (A non-scalar variable on the LHS is always indexed.) + varName <- as.character(LHSrep) + indexExprs <- list() + } else { + varName <- as.character(LHSrep[[2]]) + indexExprs <- as.list(LHSrep)[-(1:2)] + } + isRange <- vapply(indexExprs, \(e) is.call(e) && identical(e[[1]], as.name(":")), logical(1)) + D <- sum(isRange) + if (D == 0) { + C_len <- function(idx) { + return(1) + } + C_get <- substitute(function(idx, LocalOut_, LocalOff_) { + LocalOut_[LocalOff_ + 1] <- LHS + }, list(LHS = LHSrep)) |> eval() + C_set <- substitute(function(idx, LocalV_, LocalOff_) { + LHS <- LocalV_[LocalOff_ + 1] + }, list(LHS = LHSrep)) |> eval() + } else if (D == 1) { + e <- indexExprs[[which(isRange)]] + LEN <- substitute((B - A + 1), list(A = e[[2]], B = e[[3]])) + C_len <- substitute(function(idx) { + return(LEN) + }, list(LEN = LEN)) |> eval() + C_get <- substitute(function(idx, LocalOut_, LocalOff_) { + LocalOut_[(LocalOff_ + 1):(LocalOff_ + LEN)] <- LHS + }, list(LHS = LHSrep, LEN = LEN)) |> eval() + C_set <- substitute(function(idx, LocalV_, LocalOff_) { + LHS <- LocalV_[(LocalOff_ + 1):(LocalOff_ + LEN)] + }, list(LHS = LHSrep, LEN = LEN)) |> eval() + } else { + # D >= 2. + # Create local variables to calculate index starts and ends (which could be expressions like i + 1). + # Then call nm_getBlock_into_/nm_setBlock_from_ (in declFunClass_.h) in C++. + boundLines <- list() + blocks <- character() + lenExprs <- list() + for (k in seq_along(indexExprs)) { + e <- indexExprs[[k]] + if (isRange[k]) { + S <- paste0("LocalS", k, "_") + E <- paste0("LocalE", k, "_") + boundLines <- c(boundLines, + substitute(S <- A, list(S = as.name(S), A = e[[2]])), + substitute(E <- B, list(E = as.name(E), B = e[[3]]))) + blocks <- c(blocks, paste0("nm_r_(", S, ", ", E, ")")) + lenExprs <- c(lenExprs, substitute((B - A + 1), list(A = e[[2]], B = e[[3]]))) + } else { + I <- paste0("LocalI", k, "_") + boundLines <- c(boundLines, substitute(I <- A, list(I = as.name(I), A = e))) + blocks <- c(blocks, paste0("nm_i_(", I, ")")) + } + } + blocksCode <- paste0("{", paste(blocks, collapse = ", "), "}") + lenExpr <- Reduce(\(a, b) substitute(A * B, list(A = a, B = b)), lenExprs) + C_len <- substitute(function(idx) { + return(LEN) + }, list(LEN = lenExpr)) |> eval() + C_get <- as.call(c( + as.name("{"), boundLines, + substitute(nCpp(CODE), list(CODE = paste0( + "nm_getBlock_into_(", varName, ", ", blocksCode, ", LocalOut_.data() + LocalOff_)" + ))) + )) + C_get <- eval(substitute(function(idx, LocalOut_, LocalOff_) BODY, list(BODY = C_get))) + C_set <- as.call(c( + as.name("{"), boundLines, + substitute(nCpp(CODE), list(CODE = paste0( + "nm_setBlock_from_(", varName, ", ", blocksCode, ", LocalV_.data() + LocalOff_)" + ))) + )) + C_set <- eval(substitute(function(idx, LocalV_, LocalOff_) BODY, list(BODY = C_set))) + } + + # Uncompiled versions. R has no pass by reference here, so getValues_into_one returns LocalOut_. + R_len <- substitute(function(idx) { + return(length(LHS)) + }, list(LHS = LHSrep)) |> eval() + R_get <- substitute(function(idx, LocalOut_, LocalOff_) { + LocalOut_[LocalOff_ + seq_along(LHS)] <- as.numeric(LHS) + return(LocalOut_) + }, list(LHS = LHSrep)) |> eval() + R_set <- substitute(function(idx, LocalV_, LocalOff_) { + LHS <- LocalV_[LocalOff_ + seq_along(LHS)] + }, list(LHS = LHSrep)) |> eval() + exceptionNames <- c("idx", "LocalOut_", "LocalV_", "LocalOff_") + body(R_len) <- nm_addModelDollarSign(body(R_len), exceptionNames = exceptionNames) + body(R_get) <- nm_addModelDollarSign(body(R_get), exceptionNames = exceptionNames) + body(R_set) <- nm_addModelDollarSign(body(R_set), exceptionNames = exceptionNames) + list( + valuesLength_one = nFunction( + name = "valuesLength_one", + R_len, + compileInfo = list(C_fun = C_len), + returnType = "numericScalar", # index arithmetic in the DSL is double; declFunClass_ casts to int + argTypes = list(idx = "integerVector"), + refArgs = "idx" + ), + getValues_into_one = nFunction( + name = "getValues_into_one", + R_get, + compileInfo = list(C_fun = C_get), + returnType = "void", + argTypes = list(idx = "integerVector", LocalOut_ = "numericVector", LocalOff_ = "integerScalar"), + refArgs = c("idx", "LocalOut_") + ), + setValues_from_one = nFunction( + name = "setValues_from_one", + R_set, + compileInfo = list(C_fun = C_set), + returnType = "void", + argTypes = list(idx = "integerVector", LocalV_ = "numericVector", LocalOff_ = "integerScalar"), + refArgs = c("idx", "LocalV_") + ) + ) } make_getParamID_method <- function(distName) { @@ -816,7 +971,7 @@ make_getBound_one_nFxn <- function(boundExprs = list()) { ) - switch_line <- bquote(nSwitch(boundID_, 0:1)) + switch_line <- bquote(nSwitch(boundID_, 1:2)) # boundID 1 is the lower bound, 2 the upper (as in nimble) switch_line[3 + (1:(length(case_lines)))] <- case_lines Rfun <- bquote( diff --git a/nimbleModel/R/nodeRules.R b/nimbleModel/R/nodeRules.R index f31ad3b..2fb4663 100644 --- a/nimbleModel/R/nodeRules.R +++ b/nimbleModel/R/nodeRules.R @@ -541,15 +541,17 @@ calcRangeClass <- R6Class( stop("mismatch between indexing values and node-based sortIDs") } + # All of these instructions share their derived fields, so compute them once. + template <- add_instr_precalcs(list( + dims = 1, index_types = 2, lens = 1, nDim = 1, slots = 1, declID = declID + )) + class(template) <- "Rinstr" results <- lapply(seq_along(indices), \(i) { - tmp <- list( - dims = 1, index_types = 2, lens = 1, nDim = 1, slots = 1, sortID = sortID[i], - declID = declID, instr_type = 2, values = list(indices[i]) - ) - class(tmp) <- "Rinstr" + tmp <- template + tmp$sortID <- sortID[i] + tmp$values <- list(indices[i]) return(tmp) }) - class(results) <- "Rlist_Rinstr" # For ease of determining the type of the object in `makeInstrList`. return(results) } ) diff --git a/nimbleModel/R/options.R b/nimbleModel/R/options.R index db35b5f..4999db0 100644 --- a/nimbleModel/R/options.R +++ b/nimbleModel/R/options.R @@ -12,7 +12,8 @@ nodesAsChars = FALSE, verbose = TRUE, getDependenciesIncludesPredictiveNodes = TRUE, ## may be toggled off during MCMC sampler building - includeUnneededLatents = FALSE + includeUnneededLatents = FALSE, + check_values_node_inputs = TRUE # FALSE skips check_values_node_input (RHS-only and logProb checks for values node mode) ) ) diff --git a/nimbleModel/inst/generate_predefined.R b/nimbleModel/inst/generate_predefined.R index f321d0b..46cc28a 100644 --- a/nimbleModel/inst/generate_predefined.R +++ b/nimbleModel/inst/generate_predefined.R @@ -4,6 +4,10 @@ nCompiler::nCompile(instr_nClass = nimbleModel:::instr_nClass, control = list(generate_predefined = TRUE)) test <- nCompiler::nCompile(instr_nClass = nimbleModel:::instr_nClass) ## # +## # generate new predef/instrList_nClass (after instr_nClass is in place). Move to package. +nCompiler::nCompile(instrList_nClass = nimbleModel:::instrList_nClass, control = list(generate_predefined = TRUE)) +test <- nCompiler::nCompile(instrList_nClass = nimbleModel:::instrList_nClass) +## # ## # generate new predef/declFunBase_nC. Move to package. ## # You should see "#include " in the hContent ## # after declaration of declFunBase_nClass diff --git a/nimbleModel/inst/include/nimbleModel/predef/declFunBase_nC/declFunBase_nClass_cppContent.cpp b/nimbleModel/inst/include/nimbleModel/predef/declFunBase_nC/declFunBase_nClass_cppContent.cpp index 3b41185..df23380 100644 --- a/nimbleModel/inst/include/nimbleModel/predef/declFunBase_nC/declFunBase_nClass_cppContent.cpp +++ b/nimbleModel/inst/include/nimbleModel/predef/declFunBase_nC/declFunBase_nClass_cppContent.cpp @@ -51,6 +51,59 @@ return(ETaccessPtr(0.0)); RESET_EIGEN_ERRORS Rprintf("declFunBase_nClass virtual base getBound_cpp should never be called (something is wrong)\n");; return(0.0); +} + double declFunBase_nClass::calculate_singleton_cpp ( std::shared_ptr instr, int step ) { +RESET_EIGEN_ERRORS +Rprintf("declFunBase_nClass virtual base calculate_singleton_cpp should never be called (something is wrong)\n");; +return(0.0); +} + double declFunBase_nClass::calculateDiff_singleton_cpp ( std::shared_ptr instr, int step ) { +RESET_EIGEN_ERRORS +Rprintf("declFunBase_nClass virtual base calculateDiff_singleton_cpp should never be called (something is wrong)\n");; +return(0.0); +} + double declFunBase_nClass::getLogProb_singleton_cpp ( std::shared_ptr instr, int step ) { +RESET_EIGEN_ERRORS +Rprintf("declFunBase_nClass virtual base getLogProb_singleton_cpp should never be called (something is wrong)\n");; +return(0.0); +} + void declFunBase_nClass::simulate_singleton_cpp ( std::shared_ptr instr, int step ) { +RESET_EIGEN_ERRORS +Rprintf("declFunBase_nClass virtual base simulate_singleton_cpp should never be called (something is wrong)\n");; +} + std::unique_ptr declFunBase_nClass::getParam_singleton_cpp ( std::shared_ptr instr, int step, int paramID ) { +RESET_EIGEN_ERRORS +Rprintf("declFunBase_nClass virtual base getParam_singleton_cpp should never be called (something is wrong)\n");; +return(ETaccessPtr(0.0)); +} + double declFunBase_nClass::getBound_singleton_cpp ( std::shared_ptr instr, int step, int boundID ) { +RESET_EIGEN_ERRORS +Rprintf("declFunBase_nClass virtual base getBound_singleton_cpp should never be called (something is wrong)\n");; +return(0.0); +} + Eigen::Tensor declFunBase_nClass::getValues_singleton_cpp ( std::shared_ptr instr, int step ) { +RESET_EIGEN_ERRORS +Rprintf("declFunBase_nClass virtual base getValues_singleton_cpp should never be called (something is wrong)\n");; +return(createTensor(0.0, 0.0)); +} + void declFunBase_nClass::setValues_singleton_cpp ( std::shared_ptr instr, int step, Eigen::Tensor v ) { +RESET_EIGEN_ERRORS +Rprintf("declFunBase_nClass virtual base setValues_singleton_cpp should never be called (something is wrong)\n");; +} + int declFunBase_nClass::valuesLength_cpp ( std::shared_ptr instr ) { +RESET_EIGEN_ERRORS +Rprintf("declFunBase_nClass virtual base valuesLength_cpp should never be called (something is wrong)\n");; +return(0); +} + int declFunBase_nClass::getValues_into_cpp ( std::shared_ptr instr, Eigen::Tensor & out, int offset ) { +RESET_EIGEN_ERRORS +Rprintf("declFunBase_nClass virtual base getValues_into_cpp should never be called (something is wrong)\n");; +return(0); +} + int declFunBase_nClass::setValues_from_cpp ( std::shared_ptr instr, Eigen::Tensor & v, int offset ) { +RESET_EIGEN_ERRORS +Rprintf("declFunBase_nClass virtual base setValues_from_cpp should never be called (something is wrong)\n");; +return(0); } declFunBase_nClass::declFunBase_nClass ( ) { RESET_EIGEN_ERRORS @@ -78,7 +131,17 @@ method("calculateDiff_cpp", &declFunBase_nClass::calculateDiff_cpp, args({{arg(" method("getLogProb_cpp", &declFunBase_nClass::getLogProb_cpp, args({{arg("instr",copy)}})), method("simulate_cpp", &declFunBase_nClass::simulate_cpp, args({{arg("instr",copy)}})), method("getParamID_cpp", &declFunBase_nClass::getParamID_cpp, args({{arg("param",copy)}})), -method("getBound_cpp", &declFunBase_nClass::getBound_cpp, args({{arg("instr",copy)},{arg("boundID",copy)}})) +method("getBound_cpp", &declFunBase_nClass::getBound_cpp, args({{arg("instr",copy)},{arg("boundID",copy)}})), +method("calculate_singleton_cpp", &declFunBase_nClass::calculate_singleton_cpp, args({{arg("instr",copy)},{arg("step",copy)}})), +method("calculateDiff_singleton_cpp", &declFunBase_nClass::calculateDiff_singleton_cpp, args({{arg("instr",copy)},{arg("step",copy)}})), +method("getLogProb_singleton_cpp", &declFunBase_nClass::getLogProb_singleton_cpp, args({{arg("instr",copy)},{arg("step",copy)}})), +method("simulate_singleton_cpp", &declFunBase_nClass::simulate_singleton_cpp, args({{arg("instr",copy)},{arg("step",copy)}})), +method("getBound_singleton_cpp", &declFunBase_nClass::getBound_singleton_cpp, args({{arg("instr",copy)},{arg("step",copy)},{arg("boundID",copy)}})), +method("getValues_singleton_cpp", &declFunBase_nClass::getValues_singleton_cpp, args({{arg("instr",copy)},{arg("step",copy)}})), +method("setValues_singleton_cpp", &declFunBase_nClass::setValues_singleton_cpp, args({{arg("instr",copy)},{arg("step",copy)},{arg("v",copy)}})), +method("valuesLength_cpp", &declFunBase_nClass::valuesLength_cpp, args({{arg("instr",copy)}})), +method("getValues_into_cpp", &declFunBase_nClass::getValues_into_cpp, args({{arg("instr",copy)},{arg("out",ref)},{arg("offset",copy)}})), +method("setValues_from_cpp", &declFunBase_nClass::setValues_from_cpp, args({{arg("instr",copy)},{arg("v",ref)},{arg("offset",copy)}})) ) ) #endif diff --git a/nimbleModel/inst/include/nimbleModel/predef/declFunBase_nC/declFunBase_nClass_hContent.h b/nimbleModel/inst/include/nimbleModel/predef/declFunBase_nC/declFunBase_nClass_hContent.h index 4f80a92..bd27ced 100644 --- a/nimbleModel/inst/include/nimbleModel/predef/declFunBase_nC/declFunBase_nClass_hContent.h +++ b/nimbleModel/inst/include/nimbleModel/predef/declFunBase_nC/declFunBase_nClass_hContent.h @@ -20,6 +20,17 @@ class declFunBase_nClass : public interface_resolver< genericInterfaceC getParam_cpp ( std::shared_ptr instr, int paramID ) ; virtual double getBound_cpp ( std::shared_ptr instr, int boundID ) ; + virtual double calculate_singleton_cpp ( std::shared_ptr instr, int step ) ; + virtual double calculateDiff_singleton_cpp ( std::shared_ptr instr, int step ) ; + virtual double getLogProb_singleton_cpp ( std::shared_ptr instr, int step ) ; + virtual void simulate_singleton_cpp ( std::shared_ptr instr, int step ) ; + virtual std::unique_ptr getParam_singleton_cpp ( std::shared_ptr instr, int step, int paramID ) ; + virtual double getBound_singleton_cpp ( std::shared_ptr instr, int step, int boundID ) ; + virtual Eigen::Tensor getValues_singleton_cpp ( std::shared_ptr instr, int step ) ; + virtual void setValues_singleton_cpp ( std::shared_ptr instr, int step, Eigen::Tensor v ) ; + virtual int valuesLength_cpp ( std::shared_ptr instr ) ; + virtual int getValues_into_cpp ( std::shared_ptr instr, Eigen::Tensor & out, int offset ) ; + virtual int setValues_from_cpp ( std::shared_ptr instr, Eigen::Tensor & v, int offset ) ; declFunBase_nClass ( ) ; }; diff --git a/nimbleModel/inst/include/nimbleModel/predef/declFunBase_nC/declFunBase_nClass_manifest.txt b/nimbleModel/inst/include/nimbleModel/predef/declFunBase_nC/declFunBase_nClass_manifest.txt index 89e7014..da3f2c3 100644 --- a/nimbleModel/inst/include/nimbleModel/predef/declFunBase_nC/declFunBase_nClass_manifest.txt +++ b/nimbleModel/inst/include/nimbleModel/predef/declFunBase_nC/declFunBase_nClass_manifest.txt @@ -1,4 +1,4 @@ -list(saved_at = structure(1793056375.6661, class = c("POSIXct", +list(saved_at = structure(1793832597.19755, class = c("POSIXct", "POSIXt")), packet_name = "declFunBase_nClass", elements = c("preamble", "cppContent", "hContent", "filebase", "post_cpp_compiler", "copyFiles", "depends"), files = list(preamble = "declFunBase_nClass_preamble.cpp", diff --git a/nimbleModel/inst/include/nimbleModel/predef/declFunClass_/declFunClass_.h b/nimbleModel/inst/include/nimbleModel/predef/declFunClass_/declFunClass_.h index 93cad35..eca1920 100644 --- a/nimbleModel/inst/include/nimbleModel/predef/declFunClass_/declFunClass_.h +++ b/nimbleModel/inst/include/nimbleModel/predef/declFunClass_/declFunClass_.h @@ -1,12 +1,183 @@ // to be included from the predefined nodeFxnBase_nClass. -// Add "#include " to that file, -// after the declaration of nodeFxnBase_nClass. +// Add "#include " the compileInfo$HpostCode +// so it appears after the declaration of nodeFxnBase_nClass. -template +// Helpers for getValues_into_one/setValues_from_one of a node whose LHS has two or more ranges +// (see make_values_one_methods). Indices and range bounds are 1-based, as in the DSL. +// The generated code looks like +// nm_getBlock_into_(W, {nm_i_(LocalI1_), nm_r_(LocalS2_, LocalE2_), nm_r_(LocalS3_, LocalE3_)}, LocalOut_.data() + LocalOff_); +// There is no checking here: the caller (in declFunClass_) checks lengths. +inline b__ nm_i_(int i) { return b__(static_cast(i) - 1); } +inline b__ nm_r_(int s, int e) { return b__(static_cast(s) - 1, static_cast(e) - 1); } + +template +RuntimeFlatView nm_blockView_(T &var, const std::vector &blocks) { + return RuntimeFlatView(var.data(), + std::vector(var.dimensions().begin(), var.dimensions().end()), + blocks); +} +template +void nm_getBlock_into_(T &var, const std::vector &blocks, double *dest) { + nm_blockView_(var, blocks).copyIntoVector(dest); +} +template +void nm_setBlock_from_(T &var, const std::vector &blocks, const double *src) { + RuntimeFlatView view = nm_blockView_(var, blocks); + view.copyFromVector(src); +} + +// NIDX is the number of loop indices of the declaration, i.e. the length of `idx`. +template class declFunClass_ : public declFunBase_nClass { public: double v; - declFunClass_() {}; + // Reusable idx tensor for singleton runners. + Eigen::Tensor singleton_idx_; + declFunClass_() : singleton_idx_(NIDX) {}; + + // Fill `idx` for one step (node) of `instr`, with no loop. `step` is 1-based; it becomes + // 0-based here, where the counters are used as C++ indices into each range's values. + // Slots stay 1-based on the instruction and are converted at the moment of use. + // This is important for decls like z[i] <- 10*i, where i on the RHS must be 1-based. + void set_instr_idx_(const instr_nClass& instr, int step, Eigen::Tensor& idx) { + const int k = step - 1; + const int nRanges = instr.lens.size(); + for(int r = 0; r < nRanges; ++r) { + const int counter = (k / instr.strides[r]) % instr.lens[r]; + const auto& vals = instr.values->operator[](r); + const int start = instr.idxSlotStarts[r]; + if(instr.index_types[r] == 1) { + idx[instr.slots[start] - 1] = vals[0] + counter; + } else { + const int d = instr.dims[r]; + for(int p = 0; p < d; ++p) + idx[instr.slots[start + p] - 1] = vals[counter * d + p]; + } + } + } + // Run one per-node method for one step. A declaration with no loop indices + // has one node and gets a dummy idx, as in calc_0_. + template + double singleton_op_(std::shared_ptr instr, int step) { + RESET_EIGEN_ERRORS; + if constexpr (NIDX == 0) { + return (static_cast(this)->*Method)(instr->lens); + } else { + set_instr_idx_(*instr, step, singleton_idx_); + return (static_cast(this)->*Method)(singleton_idx_); + } + } + double calculate_singleton_cpp(std::shared_ptr instr, int step) override { + return singleton_op_< &Derived::calc_one >(instr, step); + } + double calculateDiff_singleton_cpp(std::shared_ptr instr, int step) override { + return singleton_op_< &Derived::calcDiff_one >(instr, step); + } + double getLogProb_singleton_cpp(std::shared_ptr instr, int step) override { + return singleton_op_< &Derived::getLogProb_one >(instr, step); + } + void simulate_singleton_cpp(std::shared_ptr instr, int step) override { + singleton_op_< &Derived::sim_one >(instr, step); + } + std::unique_ptr getParam_singleton_cpp(std::shared_ptr instr, int step, int paramID) override { + RESET_EIGEN_ERRORS; + if constexpr (NIDX == 0) { + return static_cast(this)->getParam_one(instr->lens, paramID); + } else { + set_instr_idx_(*instr, step, singleton_idx_); + return static_cast(this)->getParam_one(singleton_idx_, paramID); + } + } + double getBound_singleton_cpp(std::shared_ptr instr, int step, int boundID) override { + RESET_EIGEN_ERRORS; + if constexpr (NIDX == 0) { + return static_cast(this)->getBound_one(instr->lens, boundID); + } else { + set_instr_idx_(*instr, step, singleton_idx_); + return static_cast(this)->getBound_one(singleton_idx_, boundID); + } + } + + // values(model, nodes[i], mode = "node") and values(model, nodes, mode = "node"). + // The generated valuesLength_one, getValues_into_one and setValues_from_one handle one node + // at an offset into a vector, with no checking. The loops over nodes, the offsets and the + // length checks are here. Calls to the generated methods are not virtual. + // + // Call f(idx) for each node of `instr`, in order. A declaration with no loop indices + // has one node and gets a dummy idx, as in calc_0_. + template + void for_each_node_(const std::shared_ptr &instr, F &&f) { + if constexpr (NIDX == 0) { + f(instr->lens); + } else { + const int n = instr->nNodes; + for(int k = 1; k <= n; ++k) { + set_instr_idx_(*instr, k, singleton_idx_); + f(singleton_idx_); + } + } + } + // The total number of values in the nodes of `instr`. + int valuesLength_cpp(std::shared_ptr instr) override { + RESET_EIGEN_ERRORS; + Derived *self = static_cast(this); + int total = 0; + for_each_node_(instr, [&](Eigen::Tensor &idx) { total += static_cast(self->valuesLength_one(idx)); }); + return total; + } + // Copy the values of the nodes of `instr` into `out`, starting after `offset` (0-based). + // Returns the offset after the last value copied. + int getValues_into_cpp(std::shared_ptr instr, Eigen::Tensor &out, int offset) override { + RESET_EIGEN_ERRORS; + Derived *self = static_cast(this); + for_each_node_(instr, [&](Eigen::Tensor &idx) { + const int n = static_cast(self->valuesLength_one(idx)); + if(offset + n > out.size()) + Rcpp::stop("getValues: the output vector is too short for these nodes."); + self->getValues_into_one(idx, out, offset); + offset += n; + }); + return offset; + } + // Set the values of the nodes of `instr` from `v`, starting after `offset` (0-based). + // Returns the offset after the last value used. + int setValues_from_cpp(std::shared_ptr instr, Eigen::Tensor &v, int offset) override { + RESET_EIGEN_ERRORS; + Derived *self = static_cast(this); + for_each_node_(instr, [&](Eigen::Tensor &idx) { + const int n = static_cast(self->valuesLength_one(idx)); + if(offset + n > v.size()) + Rcpp::stop("setValues: value is too short for these nodes."); + self->setValues_from_one(idx, v, offset); + offset += n; + }); + return offset; + } + // One node (step) of `instr`. + Eigen::Tensor &singleton_idx_for_(const std::shared_ptr &instr, int step) { + if constexpr (NIDX == 0) { + return instr->lens; + } else { + set_instr_idx_(*instr, step, singleton_idx_); + return singleton_idx_; + } + } + Eigen::Tensor getValues_singleton_cpp(std::shared_ptr instr, int step) override { + RESET_EIGEN_ERRORS; + Derived *self = static_cast(this); + Eigen::Tensor &idx = singleton_idx_for_(instr, step); + Eigen::Tensor out(static_cast(self->valuesLength_one(idx))); + self->getValues_into_one(idx, out, 0); + return out; + } + void setValues_singleton_cpp(std::shared_ptr instr, int step, Eigen::Tensor v) override { + RESET_EIGEN_ERRORS; + Derived *self = static_cast(this); + Eigen::Tensor &idx = singleton_idx_for_(instr, step); + if(v.size() != static_cast(self->valuesLength_one(idx))) + Rcpp::stop("setValues: value has the wrong length for this node."); + self->setValues_from_one(idx, v, 0); + } double calculate_cpp( std::shared_ptr instr) override { return calc_op_< &Derived::calc_one >(instr); @@ -17,31 +188,40 @@ class declFunClass_ : public declFunBase_nClass { double getLogProb_cpp( std::shared_ptr instr) override { return calc_op_< &Derived::getLogProb_one >(instr); } + // sim_one returns a dummy 0, so simulate shares the calc runners and ignores their sum. + void simulate_cpp( std::shared_ptr instr) override { + calc_op_< &Derived::sim_one >(instr); + } + // Each case is compiled only if its instruction type can occur for this declaration + // (see NIDX), so a declaration doesn't instantiate runners it can never use. + // Types with only single-slot ranges fill exactly as many slots as they have ranges; + // types with a matp range fill more. template double calc_op_ ( std::shared_ptr instr ) { RESET_EIGEN_ERRORS; - int instr_type = instr->instr_type; - if(instr_type == 0) return calc_0_< Method >(instr); - if(instr_type == 1) return calc_1_seq_< Method >(instr); - if(instr_type == 2) return calc_1_mat_< Method >(instr); - if(instr_type == 3) return calc_1_matp_< Method >(instr); - if(instr_type == 4) return calc_2_seq_seq_< Method >(instr); - if(instr_type == 5) return calc_2_seq_mat_< Method >(instr); - if(instr_type == 6) return calc_2_mat_seq_< Method >(instr); - if(instr_type == 7) return calc_2_mat_mat_< Method >(instr); - if(instr_type == 8) return calc_2_seq_matp_< Method >(instr); - if(instr_type == 9) return calc_2_matp_seq_< Method >(instr); - if(instr_type == 10) return calc_2_matp_matp_< Method >(instr); - if(instr_type == 11) return calc_2_x_y_ord_< Method >(instr); - if(instr_type == 12) return calc_3_allseq_< Method >(instr); - if(instr_type == 13) return calc_3_generic_< Method >(instr); - if(instr_type == 14) return calc_4_allseq_< Method >(instr); - if(instr_type == 15) return calc_4_generic_< Method >(instr); - if(instr_type == 16) return calc_5_allseq_< Method >(instr); - if(instr_type == 17) return calc_5_generic_< Method >(instr); - if(instr_type == 18) return calc_1_matp_ord_< Method >(instr); - // TODO: we should probably error out if no method found. - return(0); + switch(instr->instr_type) { + case 0: if constexpr (NIDX == 0) return calc_0_< Method >(instr); break; + case 1: if constexpr (NIDX == 1) return calc_1_seq_< Method >(instr); break; + case 2: if constexpr (NIDX == 1) return calc_1_mat_< Method >(instr); break; + case 3: if constexpr (NIDX >= 2) return calc_1_matp_< Method >(instr); break; + case 4: if constexpr (NIDX == 2) return calc_2_seq_seq_< Method >(instr); break; + case 5: if constexpr (NIDX == 2) return calc_2_seq_mat_< Method >(instr); break; + case 6: if constexpr (NIDX == 2) return calc_2_mat_seq_< Method >(instr); break; + case 7: if constexpr (NIDX == 2) return calc_2_mat_mat_< Method >(instr); break; + case 8: if constexpr (NIDX >= 3) return calc_2_seq_matp_< Method >(instr); break; + case 9: if constexpr (NIDX >= 3) return calc_2_matp_seq_< Method >(instr); break; + case 10: if constexpr (NIDX >= 3) return calc_2_matp_matp_< Method >(instr); break; + case 11: if constexpr (NIDX == 2) return calc_2_x_y_ord_< Method >(instr); break; + case 12: if constexpr (NIDX == 3) return calc_3_allseq_< Method >(instr); break; + case 13: if constexpr (NIDX >= 3) return calc_3_generic_< Method >(instr); break; + case 14: if constexpr (NIDX == 4) return calc_4_allseq_< Method >(instr); break; + case 15: if constexpr (NIDX >= 4) return calc_4_generic_< Method >(instr); break; + case 16: if constexpr (NIDX == 5) return calc_5_allseq_< Method >(instr); break; + case 17: if constexpr (NIDX >= 5) return calc_5_generic_< Method >(instr); break; + case 18: if constexpr (NIDX >= 2) return calc_1_matp_ord_< Method >(instr); break; + } + Rcpp::stop("instruction type %d cannot occur for a declaration with %d loop indices", + instr->instr_type, NIDX); } template double calc_0_ (std::shared_ptr instr) { @@ -767,528 +947,5 @@ class declFunClass_ : public declFunBase_nClass { return idx; } - // simulate - void simulate_cpp ( std::shared_ptr instr ) override { - RESET_EIGEN_ERRORS; - int instr_type = instr->instr_type; - if(instr_type == 0) return sim_0_(instr); - if(instr_type == 1) return sim_1_seq_(instr); - if(instr_type == 2) return sim_1_mat_(instr); - if(instr_type == 3) return sim_1_matp_(instr); - if(instr_type == 4) return sim_2_seq_seq_(instr); - if(instr_type == 5) return sim_2_seq_mat_(instr); - if(instr_type == 6) return sim_2_mat_seq_(instr); - if(instr_type == 7) return sim_2_mat_mat_(instr); - if(instr_type == 8) return sim_2_seq_matp_(instr); - if(instr_type == 9) return sim_2_matp_seq_(instr); - if(instr_type == 10) return sim_2_matp_matp_(instr); - if(instr_type == 11) return sim_2_x_y_ord_(instr); - if(instr_type == 12) return sim_3_allseq_(instr); - if(instr_type == 13) return sim_3_generic_(instr); - if(instr_type == 14) return sim_4_allseq_(instr); - if(instr_type == 15) return sim_4_generic_(instr); - if(instr_type == 16) return sim_5_allseq_(instr); - if(instr_type == 17) return sim_5_generic_(instr); - if(instr_type == 18) return sim_1_matp_ord_(instr); - // TODO: we should probably error out if no method found. - } - void sim_0_ (std::shared_ptr instr) { - static_cast(this)->sim_one(instr->lens); // lens serves as a dummy here, to have the right type to pass - } - void sim_1_seq_(std::shared_ptr instr) { - int len = instr->lens[0]; - if(len < 1) return; - int iStart = instr->values->operator[](0)[0]; - int iEnd = iStart + len; - Eigen::Tensor idx(1); - for(int i = iStart; i < iEnd; ++i) { - idx[0] = i; - static_cast(this)->sim_one(idx); - } - } - void sim_1_mat_(std::shared_ptr instr) { - int len = instr->lens[0]; - const auto& vals = instr->values->operator[](0); - if(len != vals.size()) std::cout<<"len != vals.size() in sim_1_mat_"< idx(1); - for(int i = 0; i < len; ++i) { - idx[0] = vals[i]; - static_cast(this)->sim_one(idx); - } - } - void sim_1_matp_(std::shared_ptr instr) { - int len = instr->lens[0]; - int nDim = instr->nDim; - const auto& vals = instr->values->operator[](0); - if(len*nDim != vals.size()) std::cout<<"len*nDim != vals.size() in sim_1_matp_"< idx(nDim); - for(int i = 0; i < len; ++i) { - for(int p = 0; p < nDim; ++p) - idx[p] = vals[i*nDim+p]; - static_cast(this)->sim_one(idx); - } - } - void sim_1_matp_ord_(std::shared_ptr instr) { - int len = instr->lens[0]; - int nDim = instr->nDim; - const auto& vals = instr->values->operator[](0); - if(len*nDim != vals.size()) std::cout<<"len*nDim != vals.size() in sim_1_matp_"< idx(nDim); - Eigen::Tensor slots(nDim); - for(int p = 0; p < nDim; ++p) - slots[p] = instr->slots[p]-1; - for(int i = 0; i < len; ++i) { - for(int p = 0; p < nDim; ++p) - idx[slots[p]] = vals[i*nDim+p]; - static_cast(this)->sim_one(idx); - } - } - void sim_2_seq_seq_(std::shared_ptr instr) { - int len1 = instr->lens[0]; - int len2 = instr->lens[1]; - if(len1 < 1) return; - if(len2 < 1) return; - int iStart1 = instr->values->operator[](0)[0]; - int iEnd1 = iStart1 + len1; - int iStart2 = instr->values->operator[](1)[0]; - int iEnd2 = iStart2 + len2; - Eigen::Tensor idx(2); - for(int i = iStart1; i < iEnd1; ++i) { - idx[0] = i; - for(int j = iStart2; j < iEnd2; ++j) { - idx[1] = j; - static_cast(this)->sim_one(idx); - } - } - } - void sim_2_seq_mat_(std::shared_ptr instr) { - int len1 = instr->lens[0]; - int len2 = instr->lens[1]; - const auto& vals2 = instr->values->operator[](1); - if(len2 != vals2.size()) std::cout<<"len2 != vals2.size() in sim_2_seq_mat_"<values->operator[](0)[0]; - int iEnd1 = iStart1 + len1; - Eigen::Tensor idx(2); - for(int i = iStart1; i < iEnd1; ++i) { - idx[0] = i; - for(int j = 0; j < len2; ++j) { - idx[1] = vals2[j]; - static_cast(this)->sim_one(idx); - } - } - } - void sim_2_mat_seq_(std::shared_ptr instr) { - int len1 = instr->lens[0]; - int len2 = instr->lens[1]; - const auto& vals1 = instr->values->operator[](0); - if(len1 != vals1.size()) std::cout<<"len1 != vals1.size() in sim_2_mat_seq_"<values->operator[](1)[0]; - int iEnd2 = iStart2 + len2; - Eigen::Tensor idx(2); - for(int i = 0; i < len1; ++i) { - idx[0] = vals1[i]; - for(int j = iStart2; j < iEnd2; ++j) { - idx[1] = j; - static_cast(this)->sim_one(idx); - } - } - } - void sim_2_mat_mat_(std::shared_ptr instr) { - int len1 = instr->lens[0]; - int len2 = instr->lens[1]; - const auto& vals1 = instr->values->operator[](0); - const auto& vals2 = instr->values->operator[](1); - if(len1 != vals1.size()) std::cout<<"len1 != vals1.size() in sim_2_mat_mat_"< idx(2); - for(int i = 0; i < len1; ++i) { - idx[0] = vals1[i]; - for(int j = 0; j < len2; ++j) { - idx[1] = vals2[j]; - static_cast(this)->sim_one(idx); - } - } - } - void sim_2_seq_matp_(std::shared_ptr instr) { - int nDim = instr->nDim; - int dim2 = instr->dims[1]; - int len1 = instr->lens[0]; - int len2 = instr->lens[1]; - const auto& vals2 = instr->values->operator[](1); - if(len2*dim2 != vals2.size()) std::cout<<"len2*dim2 != vals2.size() in sim_2_seq_matp_"<values->operator[](0)[0]; - int iEnd1 = iStart1 + len1; - Eigen::Tensor idx(nDim); - Eigen::Tensor slots(nDim); - for(int p = 0; p < nDim; ++p) - slots[p] = instr->slots[p]-1; - for(int i = iStart1; i < iEnd1; ++i) { - // TODO: because of reordering in nodeRange determination, we should be able - // to assume that slots[0] is 1 and slots[1:(nDim-1)] is 2:nDim. - idx[slots[0]] = i; - for(int j = 0; j < len2; ++j) { - for(int p = 0; p < dim2; ++p) - idx[slots[p+1]] = vals2[j*dim2+p]; - static_cast(this)->sim_one(idx); - } - } - } - void sim_2_matp_seq_(std::shared_ptr instr) { - int nDim = instr->nDim; - int dim1 = instr->dims[0]; - int len1 = instr->lens[0]; - int len2 = instr->lens[1]; - const auto& vals1 = instr->values->operator[](0); - if(len1*dim1 != vals1.size()) std::cout<<"len1*dim1 != vals1.size() in sim_2_matp_seq_"<values->operator[](1)[0]; - int iEnd2 = iStart2 + len2; - Eigen::Tensor idx(nDim); - Eigen::Tensor slots(nDim); - for(int p = 0; p < nDim; ++p) - slots[p] = instr->slots[p]-1; - for(int i = 0; i < len1; ++i) { - for(int p = 0; p < dim1; ++p) - idx[slots[p]] = vals1[i*dim1+p]; - for(int j = iStart2; j < iEnd2; ++j) { - idx[slots[dim1]] = j; - static_cast(this)->sim_one(idx); - } - } - } - void sim_2_matp_matp_(std::shared_ptr instr) { - int nDim = instr->nDim; - int dim1 = instr->dims[0]; - int dim2 = instr->dims[1]; - int len1 = instr->lens[0]; - int len2 = instr->lens[1]; - const auto& vals1 = instr->values->operator[](0); - const auto& vals2 = instr->values->operator[](1); - if(len1*dim1 != vals1.size()) std::cout<<"len1*dim1 != vals1.size() in sim_2_matp_matp_"< idx(nDim); - Eigen::Tensor slots(nDim); - for(int p = 0; p < nDim; ++p) - slots[p] = instr->slots[p]-1; - for(int i = 0; i < len1; ++i) { - for(int p = 0; p < dim1; ++p) - idx[slots[p]] = vals1[i*dim1+p]; - for(int j = 0; j < len2; ++j) { - for(int p = 0; p < dim2; ++p) - idx[slots[p+dim1]] = vals2[j*dim2+p]; - static_cast(this)->sim_one(idx); - } - } - } - // 2 indexRanges, single slot; assumes need to reorder. - void sim_2_x_y_ord_(std::shared_ptr instr) { - if(instr->slots[0] != 2 || instr->slots[1] != 1) - std::cout<<"slots not equal to 2,1 in sim_2_seq_seq_ord_"<dims[0]; - int dim2 = instr->dims[1]; - int index_type1 = instr->index_types[0]; - int index_type2 = instr->index_types[1]; - int len1 = instr->lens[0]; - int len2 = instr->lens[1]; - const auto& vals1 = instr->values->operator[](0); - const auto& vals2 = instr->values->operator[](1); - if(len1 != vals1.size()) std::cout<<"len1 != vals1.size() in sim_3_generic_"<values->operator[](0)[0]; - int iStart2 = instr->values->operator[](1)[0]; - Eigen::Tensor idx(2); - // We could manually write out the four cases (seq-mat, seq-seq, mat-seq, mat-mat) to avoid - // conditionals inside of loops. - for(int i = 0; i < len1; ++i) { - if(index_type1 == 1) { - idx[1] = iStart1 + i; - } else idx[1] = vals1[i]; - for(int j = 0; j < len2; ++j) { - if(index_type2 == 1) { - idx[0] = iStart2 + j; - } else idx[0] = vals2[j]; - static_cast(this)->sim_one(idx); - } - } - } - void sim_3_allseq_(std::shared_ptr instr) { - int len1 = instr->lens[0]; - int len2 = instr->lens[1]; - int len3 = instr->lens[2]; - if(len1 < 1) return; - if(len2 < 1) return; - if(len3 < 1) return; - int iStart1 = instr->values->operator[](0)[0]; - int iEnd1 = iStart1 + len1; - int iStart2 = instr->values->operator[](1)[0]; - int iEnd2 = iStart2 + len2; - int iStart3 = instr->values->operator[](2)[0]; - int iEnd3 = iStart3 + len3; - Eigen::Tensor idx(3); - for(int i = iStart1; i < iEnd1; ++i) { - idx[0] = i; - for(int j = iStart2; j < iEnd2; ++j) { - idx[1] = j; - for(int k = iStart3; k < iEnd3; ++k) { - idx[2] = k; - static_cast(this)->sim_one(idx); - } - } - } - } - void sim_3_generic_(std::shared_ptr instr) { - int nDim = instr->nDim; - int dim1 = instr->dims[0]; - int dim2 = instr->dims[1]; - int dim3 = instr->dims[2]; - int cumdim2 = dim1+dim2; - int index_type1 = instr->index_types[0]; - int index_type2 = instr->index_types[1]; - int index_type3 = instr->index_types[2]; - int len1 = instr->lens[0]; - int len2 = instr->lens[1]; - int len3 = instr->lens[2]; - const auto& vals1 = instr->values->operator[](0); - const auto& vals2 = instr->values->operator[](1); - const auto& vals3 = instr->values->operator[](2); - if(index_type1 == 2 && len1*dim1 != vals1.size()) std::cout<<"len1*dim1 != vals1.size() in sim_3_generic_"<values->operator[](0)[0]; - int iStart2 = instr->values->operator[](1)[0]; - int iStart3 = instr->values->operator[](2)[0]; - Eigen::Tensor idx(nDim); - Eigen::Tensor slots(nDim); - for(int p = 0; p < nDim; ++p) - slots[p] = instr->slots[p]-1; - for(int i = 0; i < len1; ++i) { - if(index_type1 == 1) { - idx[slots[0]] = iStart1 + i; - } else - for(int p = 0; p < dim1; ++p) - idx[slots[p]] = vals1[i*dim1+p]; - for(int j = 0; j < len2; ++j) { - if(index_type2 == 1) { - idx[slots[dim1]] = iStart2 + j; - } else - for(int p = 0; p < dim2; ++p) - idx[slots[p+dim1]] = vals2[j*dim2+p]; - for(int k = 0; k < len3; ++k) { - if(index_type3 == 1) { - idx[slots[cumdim2]] = iStart3 + k; - } else - for(int p = 0; p < dim3; ++p) - idx[slots[p+cumdim2]] = vals3[k*dim3+p]; - static_cast(this)->sim_one(idx); - } - } - } - } - void sim_4_allseq_(std::shared_ptr instr) { - int nDim = instr->nDim; - Eigen::Tensor lens(4); - Eigen::Tensor iStart(4); - Eigen::Tensor iEnd(4); - for(int p = 0; p < 4; ++p) { - lens[p] = instr->lens[p]; - iStart[p] = instr->values->operator[](p)[0]; - iEnd[p] = iStart[p] + lens[p]; - if(lens[p] < 1) return; - } - Eigen::Tensor idx(4); - for(int i0 = iStart[0]; i0 < iEnd[0]; ++i0) { - idx[0] = i0; - for(int i1 = iStart[1]; i1 < iEnd[1]; ++i1) { - idx[1] = i1; - for(int i2 = iStart[2]; i2 < iEnd[2]; ++i2) { - idx[2] = i2; - for(int i3 = iStart[3]; i3 < iEnd[3]; ++i3) { - idx[3] = i3; - static_cast(this)->sim_one(idx); - } - } - } - } - } - void sim_4_generic_(std::shared_ptr instr) { - int nDim = instr->nDim; - Eigen::Tensor lens(4); - // Some of the indexRanges might be seq ranges. - Eigen::Tensor iStart(4); - Eigen::Tensor index_types(4); - Eigen::Tensor dims(4); - for(int p = 0; p < 4; ++p) { - lens[p] = instr->lens[p]; - iStart[p] = instr->values->operator[](p)[0]; - index_types[p] = instr->index_types[p]; - dims[p] = instr->dims[p]; - if(lens[p] < 1) return; - } - Eigen::Tensor slots(nDim); - for(int p = 0; p < nDim; ++p) { - slots[p] = instr->slots[p]-1; - } - int cumdim1 = dims[0]+dims[1]; - int cumdim2 = dims[0]+dims[1]+dims[2]; - const auto& vals0 = instr->values->operator[](0); - const auto& vals1 = instr->values->operator[](1); - const auto& vals2 = instr->values->operator[](2); - const auto& vals3 = instr->values->operator[](3); - if(index_types[0] == 2 && lens[0]*dims[0] != vals0.size()) std::cout<<"lens[0]*dims[0] != vals0.size() in sim_4_generic_"< idx(nDim); - for(int i0 = 0; i0 < lens[0]; ++i0) { - if(index_types[0] == 1) { - idx[slots[0]] = iStart[0] + i0; - } else - for(int p = 0; p < dims[0]; ++p) - idx[slots[p]] = vals0[i0*dims[0]+p]; - for(int i1 = 0; i1 < lens[1]; ++i1) { - if(index_types[1] == 1) { - idx[slots[dims[0]]] = iStart[1] + i1; - } else - for(int p = 0; p < dims[1]; ++p) - idx[slots[p+dims[0]]] = vals1[i1*dims[1]+p]; - for(int i2 = 0; i2 < lens[2]; ++i2) { - if(index_types[2] == 1) { - idx[slots[cumdim1]] = iStart[2] + i2; - } else - for(int p = 0; p < dims[2]; ++p) - idx[slots[p+cumdim1]] = vals2[i2*dims[2]+p]; - for(int i3 = 0; i3 < lens[3]; ++i3) { - if(index_types[3] == 1) { - idx[slots[cumdim2]] = iStart[3] + i3; - } else - for(int p = 0; p < dims[3]; ++p) - idx[slots[p+cumdim2]] = vals3[i3*dims[3]+p]; - static_cast(this)->sim_one(idx); - } - } - } - } - } - void sim_5_allseq_(std::shared_ptr instr) { - int nDim = instr->nDim; - Eigen::Tensor lens(5); - Eigen::Tensor iStart(5); - Eigen::Tensor iEnd(5); - for(int p = 0; p < 5; ++p) { - lens[p] = instr->lens[p]; - iStart[p] = instr->values->operator[](p)[0]; - iEnd[p] = iStart[p] + lens[p]; - if(lens[p] < 1) return; - } - Eigen::Tensor idx(5); - for(int i0 = iStart[0]; i0 < iEnd[0]; ++i0) { - idx[0] = i0; - for(int i1 = iStart[1]; i1 < iEnd[1]; ++i1) { - idx[1] = i1; - for(int i2 = iStart[2]; i2 < iEnd[2]; ++i2) { - idx[2] = i2; - for(int i3 = iStart[3]; i3 < iEnd[3]; ++i3) { - idx[3] = i3; - for(int i4 = iStart[4]; i4 < iEnd[4]; ++i4) { - idx[4] = i4; - static_cast(this)->sim_one(idx); - } - } - } - } - } - } - void sim_5_generic_(std::shared_ptr instr) { - int nDim = instr->nDim; - Eigen::Tensor lens(5); - // Some of the indexRanges might be seq ranges. - Eigen::Tensor iStart(5); - Eigen::Tensor index_types(5); - Eigen::Tensor dims(5); - for(int p = 0; p < 5; ++p) { - lens[p] = instr->lens[p]; - iStart[p] = instr->values->operator[](p)[0]; - index_types[p] = instr->index_types[p]; - dims[p] = instr->dims[p]; - if(lens[p] < 1) return; - } - Eigen::Tensor slots(nDim); - for(int p = 0; p < nDim; ++p) { - slots[p] = instr->slots[p]-1; - } - int cumdim1 = dims[0]+dims[1]; - int cumdim2 = dims[0]+dims[1]+dims[2]; - int cumdim3 = dims[0]+dims[1]+dims[2]+dims[3]; - const auto& vals0 = instr->values->operator[](0); - const auto& vals1 = instr->values->operator[](1); - const auto& vals2 = instr->values->operator[](2); - const auto& vals3 = instr->values->operator[](3); - const auto& vals4 = instr->values->operator[](4); - if(index_types[0] == 2 && lens[0]*dims[0] != vals0.size()) std::cout<<"lens[0]*dims[0] != vals0.size() in sim_5_generic_"< idx(nDim); - for(int i0 = 0; i0 < lens[0]; ++i0) { - if(index_types[0] == 1) { - idx[slots[0]] = iStart[0] + i0; - } else - for(int p = 0; p < dims[0]; ++p) - idx[slots[p]] = vals0[i0*dims[0]+p]; - for(int i1 = 0; i1 < lens[1]; ++i1) { - if(index_types[1] == 1) { - idx[slots[dims[0]]] = iStart[1] + i1; - } else - for(int p = 0; p < dims[1]; ++p) - idx[slots[p+dims[0]]] = vals1[i1*dims[1]+p]; - for(int i2 = 0; i2 < lens[2]; ++i2) { - if(index_types[2] == 1) { - idx[slots[cumdim1]] = iStart[2] + i2; - } else - for(int p = 0; p < dims[2]; ++p) - idx[slots[p+cumdim1]] = vals2[i2*dims[2]+p]; - for(int i3 = 0; i3 < lens[3]; ++i3) { - if(index_types[3] == 1) { - idx[slots[cumdim2]] = iStart[3] + i3; - } else - for(int p = 0; p < dims[3]; ++p) - idx[slots[p+cumdim2]] = vals3[i3*dims[3]+p]; - for(int i4 = 0; i4 < lens[4]; ++i4) { - if(index_types[4] == 1) { - idx[slots[cumdim3]] = iStart[4] + i4; - } else - for(int p = 0; p < dims[4]; ++p) - idx[slots[p+cumdim3]] = vals4[i4*dims[4]+p]; - static_cast(this)->sim_one(idx); - } - } - } - } - } - } virtual ~declFunClass_() {}; }; diff --git a/nimbleModel/inst/include/nimbleModel/predef/instrList_nClass/instrList_nClass_copyFiles.txt b/nimbleModel/inst/include/nimbleModel/predef/instrList_nClass/instrList_nClass_copyFiles.txt new file mode 100644 index 0000000..e69de29 diff --git a/nimbleModel/inst/include/nimbleModel/predef/instrList_nClass/instrList_nClass_cppContent.cpp b/nimbleModel/inst/include/nimbleModel/predef/instrList_nClass/instrList_nClass_cppContent.cpp new file mode 100644 index 0000000..1843556 --- /dev/null +++ b/nimbleModel/inst/include/nimbleModel/predef/instrList_nClass/instrList_nClass_cppContent.cpp @@ -0,0 +1,48 @@ +/* OPENER (Do not edit this comment) */ +#ifndef __instrList_nClass_CPP +#define __instrList_nClass_CPP +/* BODY (Do not edit this comment) */ +#ifndef R_NO_REMAP +#define R_NO_REMAP +#endif +#include +#include "instrList_nClass_c_.h" +using namespace Rcpp; +// [[Rcpp::plugins(nCompiler_Eigen_plugin)]] +// [[Rcpp::depends(RcppParallel)]] +// [[Rcpp::depends(nCompiler)]] +// [[Rcpp::depends(Rcereal)]] + + instrList_nClass::instrList_nClass ( ) { +RESET_EIGEN_ERRORS +instr_nList = nClass_builder()(); +} + +// [[Rcpp::export(name = "instrList_nClass_new")]] + SEXP new_instrList_nClass ( ) { +RESET_EIGEN_ERRORS +return CREATE_NEW_NCOMP_OBJECT(instrList_nClass);; +} + +// [[Rcpp::export(name = "set_CnClass_env_instrList_nClass_new")]] + void set_CnClass_env_instrList_nClass ( SEXP env ) { +RESET_EIGEN_ERRORS +SET_CNCLASS_ENV(instrList_nClass, env);; +} + +// [[Rcpp::export(name = "get_CnClass_env_instrList_nClass_new")]] + Rcpp::Environment get_CnClass_env_instrList_nClass ( ) { +RESET_EIGEN_ERRORS +return GET_CNCLASS_ENV(instrList_nClass);; +} + +NCOMPILER_INTERFACE( +instrList_nClass, +NCOMPILER_FIELDS( +field("instr_nList", &instrList_nClass::instr_nList), +field("offsets", &instrList_nClass::offsets), +field("nNodes_total", &instrList_nClass::nNodes_total) +), +NCOMPILER_METHODS() +) +#endif diff --git a/nimbleModel/inst/include/nimbleModel/predef/instrList_nClass/instrList_nClass_depends.txt b/nimbleModel/inst/include/nimbleModel/predef/instrList_nClass/instrList_nClass_depends.txt new file mode 100644 index 0000000..e69de29 diff --git a/nimbleModel/inst/include/nimbleModel/predef/instrList_nClass/instrList_nClass_filebase.txt b/nimbleModel/inst/include/nimbleModel/predef/instrList_nClass/instrList_nClass_filebase.txt new file mode 100644 index 0000000..9f932a3 --- /dev/null +++ b/nimbleModel/inst/include/nimbleModel/predef/instrList_nClass/instrList_nClass_filebase.txt @@ -0,0 +1 @@ +instrList_nClass_c_ diff --git a/nimbleModel/inst/include/nimbleModel/predef/instrList_nClass/instrList_nClass_hContent.h b/nimbleModel/inst/include/nimbleModel/predef/instrList_nClass/instrList_nClass_hContent.h new file mode 100644 index 0000000..415134f --- /dev/null +++ b/nimbleModel/inst/include/nimbleModel/predef/instrList_nClass/instrList_nClass_hContent.h @@ -0,0 +1,29 @@ +/* OPENER (Do not edit this comment) */ +#ifndef __instrList_nClass_H +#define __instrList_nClass_H +/* BODY (Do not edit this comment) */ +#ifndef R_NO_REMAP +#define R_NO_REMAP +#endif +#include +#include "nList_instr_nClass_c_.h" +#include "instrList_nClass_c_.h" + +class instrList_nClass : public interface_resolver< genericInterfaceC >, public loadedObjectHookC { +public: + static const char* nc_class_key() {return "instrList_nClass";} + instrList_nClass ( ) ; + std::shared_ptr instr_nList; + Eigen::Tensor offsets; + int nNodes_total; + +}; + + SEXP new_instrList_nClass ( ) ; + + void set_CnClass_env_instrList_nClass ( SEXP env ) ; + + Rcpp::Environment get_CnClass_env_instrList_nClass ( ) ; + + +#endif diff --git a/nimbleModel/inst/include/nimbleModel/predef/instrList_nClass/instrList_nClass_manifest.txt b/nimbleModel/inst/include/nimbleModel/predef/instrList_nClass/instrList_nClass_manifest.txt new file mode 100644 index 0000000..1c83b72 --- /dev/null +++ b/nimbleModel/inst/include/nimbleModel/predef/instrList_nClass/instrList_nClass_manifest.txt @@ -0,0 +1,7 @@ +list(saved_at = structure(1793298692.68275, class = c("POSIXct", +"POSIXt")), packet_name = "instrList_nClass", elements = c("preamble", +"cppContent", "hContent", "filebase", "post_cpp_compiler", "copyFiles", +"depends"), files = list(preamble = "instrList_nClass_preamble.cpp", + cppContent = "instrList_nClass_cppContent.cpp", hContent = "instrList_nClass_hContent.h", + filebase = "instrList_nClass_filebase.txt", post_cpp_compiler = "instrList_nClass_post_cpp_compiler.txt", + copyFiles = "instrList_nClass_copyFiles.txt", depends = "instrList_nClass_depends.txt")) diff --git a/nimbleModel/inst/include/nimbleModel/predef/instrList_nClass/instrList_nClass_post_cpp_compiler.txt b/nimbleModel/inst/include/nimbleModel/predef/instrList_nClass/instrList_nClass_post_cpp_compiler.txt new file mode 100644 index 0000000..e69de29 diff --git a/nimbleModel/inst/include/nimbleModel/predef/instrList_nClass/instrList_nClass_preamble.cpp b/nimbleModel/inst/include/nimbleModel/predef/instrList_nClass/instrList_nClass_preamble.cpp new file mode 100644 index 0000000..0c6b4c1 --- /dev/null +++ b/nimbleModel/inst/include/nimbleModel/predef/instrList_nClass/instrList_nClass_preamble.cpp @@ -0,0 +1,6 @@ +#define NCOMPILER_HANDLE_EIGEN_ERRORS +#define NCOMPILER_USES_EIGEN +#define NCOMPILER_USES_NCPPVEC +#define NCOMPILER_USES_MVDISTS +#define USES_NCOMPILER +#define NCOMPILER_USES_NCLASS_INTERFACE diff --git a/nimbleModel/inst/include/nimbleModel/predef/instr_nClass/instr_nClass_cppContent.cpp b/nimbleModel/inst/include/nimbleModel/predef/instr_nClass/instr_nClass_cppContent.cpp index eec485f..b62c359 100644 --- a/nimbleModel/inst/include/nimbleModel/predef/instr_nClass/instr_nClass_cppContent.cpp +++ b/nimbleModel/inst/include/nimbleModel/predef/instr_nClass/instr_nClass_cppContent.cpp @@ -47,7 +47,10 @@ field("slots", &instr_nClass::slots), field("values", &instr_nClass::values), field("instr_type", &instr_nClass::instr_type), field("sortID", &instr_nClass::sortID), -field("declID", &instr_nClass::declID) +field("declID", &instr_nClass::declID), +field("nNodes", &instr_nClass::nNodes), +field("strides", &instr_nClass::strides), +field("idxSlotStarts", &instr_nClass::idxSlotStarts) ), NCOMPILER_METHODS() ) diff --git a/nimbleModel/inst/include/nimbleModel/predef/instr_nClass/instr_nClass_hContent.h b/nimbleModel/inst/include/nimbleModel/predef/instr_nClass/instr_nClass_hContent.h index 489c5ce..80ea76d 100644 --- a/nimbleModel/inst/include/nimbleModel/predef/instr_nClass/instr_nClass_hContent.h +++ b/nimbleModel/inst/include/nimbleModel/predef/instr_nClass/instr_nClass_hContent.h @@ -22,6 +22,9 @@ class instr_nClass : public interface_resolver< genericInterfaceC int instr_type; Eigen::Tensor sortID; int declID; + int nNodes; + Eigen::Tensor strides; + Eigen::Tensor idxSlotStarts; }; diff --git a/nimbleModel/inst/include/nimbleModel/predef/mCopy_nC/mCopy_nC_cppContent.cpp b/nimbleModel/inst/include/nimbleModel/predef/mCopy_nC/mCopy_nC_cppContent.cpp index 2c4ed20..db80f42 100644 --- a/nimbleModel/inst/include/nimbleModel/predef/mCopy_nC/mCopy_nC_cppContent.cpp +++ b/nimbleModel/inst/include/nimbleModel/predef/mCopy_nC/mCopy_nC_cppContent.cpp @@ -20,11 +20,19 @@ multiCopier_model_nC_base::init(this->copiers, model); } Eigen::Tensor mCopy_nC::getValues ( ) { RESET_EIGEN_ERRORS -return flatViewGroup.copyIntoVector(); +return getValues_(); } void mCopy_nC::setValues ( Eigen::Tensor v ) { RESET_EIGEN_ERRORS -flatViewGroup.copyFromVector(v); +setValues_() = v; +} + Eigen::Tensor mCopy_nC::getValues_indexed ( int i, std::string mode ) { +RESET_EIGEN_ERRORS +return getValues_indexed_(i, mode); +} + void mCopy_nC::setValues_indexed ( int i, Eigen::Tensor v, std::string mode ) { +RESET_EIGEN_ERRORS +setValues_indexed_(i, v, mode); } mCopy_nC::mCopy_nC ( ) { RESET_EIGEN_ERRORS @@ -56,7 +64,9 @@ field("copiers", &mCopy_nC::copiers) NCOMPILER_METHODS( method("init", &mCopy_nC::init, args({{arg("model",copy)}})), method("getValues", &mCopy_nC::getValues, args({{}})), -method("setValues", &mCopy_nC::setValues, args({{arg("v",copy)}})) +method("setValues", &mCopy_nC::setValues, args({{arg("v",copy)}})), +method("getValues_indexed", &mCopy_nC::getValues_indexed, args({{arg("i",copy)},{arg("mode",copy)}})), +method("setValues_indexed", &mCopy_nC::setValues_indexed, args({{arg("i",copy)},{arg("v",copy)},{arg("mode",copy)}})) ) ) #endif diff --git a/nimbleModel/inst/include/nimbleModel/predef/mCopy_nC/mCopy_nC_hContent.h b/nimbleModel/inst/include/nimbleModel/predef/mCopy_nC/mCopy_nC_hContent.h index 4b5c95d..b33a1fc 100644 --- a/nimbleModel/inst/include/nimbleModel/predef/mCopy_nC/mCopy_nC_hContent.h +++ b/nimbleModel/inst/include/nimbleModel/predef/mCopy_nC/mCopy_nC_hContent.h @@ -17,6 +17,8 @@ class mCopy_nC : public interface_resolver< genericInterfaceC, multiCo void init ( std::shared_ptr model ) ; Eigen::Tensor getValues ( ) ; void setValues ( Eigen::Tensor v ) ; + Eigen::Tensor getValues_indexed ( int i, std::string mode ) ; + void setValues_indexed ( int i, Eigen::Tensor v, std::string mode ) ; mCopy_nC ( ) ; std::shared_ptr copiers; diff --git a/nimbleModel/inst/include/nimbleModel/predef/mCopy_nC/mCopy_nC_manifest.txt b/nimbleModel/inst/include/nimbleModel/predef/mCopy_nC/mCopy_nC_manifest.txt index d6a3595..c2676e7 100644 --- a/nimbleModel/inst/include/nimbleModel/predef/mCopy_nC/mCopy_nC_manifest.txt +++ b/nimbleModel/inst/include/nimbleModel/predef/mCopy_nC/mCopy_nC_manifest.txt @@ -1,4 +1,4 @@ -list(saved_at = structure(1793056548.82963, class = c("POSIXct", +list(saved_at = structure(1793756732.99272, class = c("POSIXct", "POSIXt")), packet_name = "mCopy_nC", elements = c("preamble", "cppContent", "hContent", "filebase", "post_cpp_compiler", "copyFiles", "depends"), files = list(preamble = "mCopy_nC_preamble.cpp", diff --git a/nimbleModel/inst/include/nimbleModel/predef/mCopy_nC_base/mCopy_model_nC_base.h b/nimbleModel/inst/include/nimbleModel/predef/mCopy_nC_base/mCopy_model_nC_base.h index ba22213..76561a0 100644 --- a/nimbleModel/inst/include/nimbleModel/predef/mCopy_nC_base/mCopy_model_nC_base.h +++ b/nimbleModel/inst/include/nimbleModel/predef/mCopy_nC_base/mCopy_model_nC_base.h @@ -45,5 +45,54 @@ class multiCopier_model_nC_base : public multiCopier_nC_base { for (auto& acc : accessors) flatViewGroup.setNextData(acc->template S().data()); // update each one's data pointer. } + + // Get and set all values, refreshing the data pointers first in case a + // model variable has been reallocated since init(). setValues_() returns + // a proxy so that `setValues_() = v;` does the copy. + Eigen::Tensor getValues_() { + setActiveRow(0); + return flatViewGroup.copyIntoVector(); + } + RuntimeFlatViewGroup::SetValuesProxy setValues_() { + setActiveRow(0); + return flatViewGroup.setValues_(); + } + + // Get and set one entry, for values(model, nodes[i], mode = ...). i is 1-based. + // "raw" mode: entry i of nodes (one view each). + // "scalar" mode: element i of the concatenated values, returned as a length-1 vector. + // The setters return proxies, as setValues_() does. + Eigen::Tensor getValues_raw_(int i) { + setActiveRow(0); + return flatViewGroup.copyIntoVector_one(i - 1); + } + RuntimeFlatViewGroup::SetValuesOneProxy setValues_raw_(int i) { + setActiveRow(0); + return flatViewGroup.setValues_one_(i - 1); + } + Eigen::Tensor getValues_scalar_(int i) { + setActiveRow(0); + Eigen::Tensor ans(1); + ans(0) = flatViewGroup.at(i - 1); + return ans; + } + RuntimeFlatViewGroup::SetValuesAtProxy setValues_scalar_(int i) { + setActiveRow(0); + return flatViewGroup.setValues_at_(i - 1); + } + + // For use from R, via + // multiCopier_nClass$getValues_indexed and $setValues_indexed, + // allowing run-time choice of mode. + Eigen::Tensor getValues_indexed_(int i, const std::string &mode) { + if (mode == "raw") return getValues_raw_(i); + if (mode == "scalar") return getValues_scalar_(i); + Rcpp::stop("getValues_indexed: mode must be \"raw\" or \"scalar\"."); + } + void setValues_indexed_(int i, const Eigen::Tensor &v, const std::string &mode) { + if (mode == "raw") {setValues_raw_(i) = v; return;} + if (mode == "scalar") {setValues_scalar_(i) = v; return;} + Rcpp::stop("setValues_indexed: mode must be \"raw\" or \"scalar\"."); + } }; #endif diff --git a/nimbleModel/inst/include/nimbleModel/predef/modelBase_nC/modelBase_nClass_cppContent.cpp b/nimbleModel/inst/include/nimbleModel/predef/modelBase_nC/modelBase_nClass_cppContent.cpp index 1078d89..925ea25 100644 --- a/nimbleModel/inst/include/nimbleModel/predef/modelBase_nC/modelBase_nClass_cppContent.cpp +++ b/nimbleModel/inst/include/nimbleModel/predef/modelBase_nC/modelBase_nClass_cppContent.cpp @@ -18,10 +18,10 @@ using namespace Rcpp; RESET_EIGEN_ERRORS return(true); } - std::shared_ptr modelBase_nClass::makeCompiledInstrList ( SEXP input ) { + std::shared_ptr modelBase_nClass::makeCompiledInstrList ( SEXP input ) { RESET_EIGEN_ERRORS -std::shared_ptr ans; -ans = nClass_builder()(); +std::shared_ptr ans; +ans = nClass_builder()(); ans->set_all_values(input);; return(ans); } @@ -32,22 +32,22 @@ ans = nClass_builder()(); ans->set_all_values(input);; return(ans); } - double modelBase_nClass::calculate_impl ( std::shared_ptr instrList ) { + double modelBase_nClass::calculate_impl ( std::shared_ptr instrList ) { RESET_EIGEN_ERRORS Rprintf("modelBase_nClass calculate_impl (should not see this)\n");; return(0.0); } - double modelBase_nClass::calculateDiff_impl ( std::shared_ptr instrList ) { + double modelBase_nClass::calculateDiff_impl ( std::shared_ptr instrList ) { RESET_EIGEN_ERRORS Rprintf("modelBase_nClass calculateDiff_impl (should not see this)\n");; return(0.0); } - double modelBase_nClass::getLogProb_impl ( std::shared_ptr instrList ) { + double modelBase_nClass::getLogProb_impl ( std::shared_ptr instrList ) { RESET_EIGEN_ERRORS Rprintf("modelBase_nClass getLogProb_impl (should not see this)\n");; return(0.0); } - void modelBase_nClass::simulate_impl ( std::shared_ptr instrList ) { + void modelBase_nClass::simulate_impl ( std::shared_ptr instrList ) { RESET_EIGEN_ERRORS Rprintf("modelBase_nClass simulate_impl (should not see this)\n");; } @@ -63,6 +63,52 @@ Rprintf("modelBase_nClass getParam_impl (should not see this)\n"); return nullpt RESET_EIGEN_ERRORS Rprintf("modelBase_nClass getBound_impl (should not see this)\n");; return(0.0); +} + double modelBase_nClass::calculate_indexed_impl ( std::shared_ptr instrList, int step ) { +RESET_EIGEN_ERRORS +Rprintf("modelBase_nClass calculate_indexed_impl (should not see this)\n");; +return(0.0); +} + double modelBase_nClass::calculateDiff_indexed_impl ( std::shared_ptr instrList, int step ) { +RESET_EIGEN_ERRORS +Rprintf("modelBase_nClass calculateDiff_indexed_impl (should not see this)\n");; +return(0.0); +} + double modelBase_nClass::getLogProb_indexed_impl ( std::shared_ptr instrList, int step ) { +RESET_EIGEN_ERRORS +Rprintf("modelBase_nClass getLogProb_indexed_impl (should not see this)\n");; +return(0.0); +} + void modelBase_nClass::simulate_indexed_impl ( std::shared_ptr instrList, int step ) { +RESET_EIGEN_ERRORS +Rprintf("modelBase_nClass simulate_indexed_impl (should not see this)\n");; +} + std::unique_ptr modelBase_nClass::getParam_indexed_impl ( std::shared_ptr instrList, int step, int paramID ) { +RESET_EIGEN_ERRORS +Rprintf("modelBase_nClass getParam_indexed_impl (should not see this)\n"); return nullptr;; +} + double modelBase_nClass::getBound_indexed_impl ( std::shared_ptr instrList, int step, int boundID ) { +RESET_EIGEN_ERRORS +Rprintf("modelBase_nClass getBound_indexed_impl (should not see this)\n");; +return(0.0); +} + Eigen::Tensor modelBase_nClass::getValues_indexed_impl ( std::shared_ptr instrList, int step ) { +RESET_EIGEN_ERRORS +Rprintf("modelBase_nClass getValues_indexed_impl (should not see this)\n");; +return(createTensor(0.0, 0.0)); +} + void modelBase_nClass::setValues_indexed_impl ( std::shared_ptr instrList, int step, Eigen::Tensor value ) { +RESET_EIGEN_ERRORS +Rprintf("modelBase_nClass setValues_indexed_impl (should not see this)\n");; +} + Eigen::Tensor modelBase_nClass::getValues_nodes_impl ( std::shared_ptr instrList ) { +RESET_EIGEN_ERRORS +Rprintf("modelBase_nClass getValues_nodes_impl (should not see this)\n");; +return(createTensor(0.0, 0.0)); +} + void modelBase_nClass::setValues_nodes_impl ( std::shared_ptr instrList, Eigen::Tensor value ) { +RESET_EIGEN_ERRORS +Rprintf("modelBase_nClass setValues_nodes_impl (should not see this)\n");; } modelBase_nClass::modelBase_nClass ( ) { RESET_EIGEN_ERRORS @@ -93,7 +139,17 @@ method("getLogProb_impl", &modelBase_nClass::getLogProb_impl, args({{arg("instrL method("simulate_impl", &modelBase_nClass::simulate_impl, args({{arg("instrList",copy)}})), method("getParamID_impl", &modelBase_nClass::getParamID_impl, args({{arg("node",copy)},{arg("param",copy)}})), method("getParam_impl", &modelBase_nClass::getParam_impl, args({{arg("node",copy)},{arg("paramID",copy)}})), -method("getBound_impl", &modelBase_nClass::getBound_impl, args({{arg("node",copy)},{arg("boundID",copy)}})) +method("getBound_impl", &modelBase_nClass::getBound_impl, args({{arg("node",copy)},{arg("boundID",copy)}})), +method("calculate_indexed_impl", &modelBase_nClass::calculate_indexed_impl, args({{arg("instrList",copy)},{arg("step",copy)}})), +method("calculateDiff_indexed_impl", &modelBase_nClass::calculateDiff_indexed_impl, args({{arg("instrList",copy)},{arg("step",copy)}})), +method("getLogProb_indexed_impl", &modelBase_nClass::getLogProb_indexed_impl, args({{arg("instrList",copy)},{arg("step",copy)}})), +method("simulate_indexed_impl", &modelBase_nClass::simulate_indexed_impl, args({{arg("instrList",copy)},{arg("step",copy)}})), +method("getParam_indexed_impl", &modelBase_nClass::getParam_indexed_impl, args({{arg("instrList",copy)},{arg("step",copy)},{arg("paramID",copy)}})), +method("getBound_indexed_impl", &modelBase_nClass::getBound_indexed_impl, args({{arg("instrList",copy)},{arg("step",copy)},{arg("boundID",copy)}})), +method("getValues_indexed_impl", &modelBase_nClass::getValues_indexed_impl, args({{arg("instrList",copy)},{arg("step",copy)}})), +method("setValues_indexed_impl", &modelBase_nClass::setValues_indexed_impl, args({{arg("instrList",copy)},{arg("step",copy)},{arg("value",copy)}})), +method("getValues_nodes_impl", &modelBase_nClass::getValues_nodes_impl, args({{arg("instrList",copy)}})), +method("setValues_nodes_impl", &modelBase_nClass::setValues_nodes_impl, args({{arg("instrList",copy)},{arg("value",copy)}})) ) ) #endif diff --git a/nimbleModel/inst/include/nimbleModel/predef/modelBase_nC/modelBase_nClass_hContent.h b/nimbleModel/inst/include/nimbleModel/predef/modelBase_nC/modelBase_nClass_hContent.h index d1cefef..ac39cb1 100644 --- a/nimbleModel/inst/include/nimbleModel/predef/modelBase_nC/modelBase_nClass_hContent.h +++ b/nimbleModel/inst/include/nimbleModel/predef/modelBase_nC/modelBase_nClass_hContent.h @@ -7,7 +7,7 @@ #endif #include #include "modelBase_nClass_c_.h" -#include "nList_instr_nClass_c_.h" +#include "instrList_nClass_c_.h" #include "instr_nClass_c_.h" #include "declFunBase_nClass_c_.h" @@ -15,15 +15,25 @@ class modelBase_nClass : public interface_resolver< genericInterfaceC makeCompiledInstrList ( SEXP input ) ; + std::shared_ptr makeCompiledInstrList ( SEXP input ) ; std::shared_ptr makeCompiledInstr ( SEXP input ) ; - virtual double calculate_impl ( std::shared_ptr instrList ) ; - virtual double calculateDiff_impl ( std::shared_ptr instrList ) ; - virtual double getLogProb_impl ( std::shared_ptr instrList ) ; - virtual void simulate_impl ( std::shared_ptr instrList ) ; + virtual double calculate_impl ( std::shared_ptr instrList ) ; + virtual double calculateDiff_impl ( std::shared_ptr instrList ) ; + virtual double getLogProb_impl ( std::shared_ptr instrList ) ; + virtual void simulate_impl ( std::shared_ptr instrList ) ; virtual int getParamID_impl ( std::shared_ptr node, std::string param ) ; virtual std::unique_ptr getParam_impl ( std::shared_ptr node, int paramID ) ; virtual double getBound_impl ( std::shared_ptr node, int boundID ) ; + virtual double calculate_indexed_impl ( std::shared_ptr instrList, int step ) ; + virtual double calculateDiff_indexed_impl ( std::shared_ptr instrList, int step ) ; + virtual double getLogProb_indexed_impl ( std::shared_ptr instrList, int step ) ; + virtual void simulate_indexed_impl ( std::shared_ptr instrList, int step ) ; + virtual std::unique_ptr getParam_indexed_impl ( std::shared_ptr instrList, int step, int paramID ) ; + virtual double getBound_indexed_impl ( std::shared_ptr instrList, int step, int boundID ) ; + virtual Eigen::Tensor getValues_indexed_impl ( std::shared_ptr instrList, int step ) ; + virtual void setValues_indexed_impl ( std::shared_ptr instrList, int step, Eigen::Tensor value ) ; + virtual Eigen::Tensor getValues_nodes_impl ( std::shared_ptr instrList ) ; + virtual void setValues_nodes_impl ( std::shared_ptr instrList, Eigen::Tensor value ) ; modelBase_nClass ( ) ; }; diff --git a/nimbleModel/inst/include/nimbleModel/predef/modelBase_nC/modelBase_nClass_manifest.txt b/nimbleModel/inst/include/nimbleModel/predef/modelBase_nC/modelBase_nClass_manifest.txt index cc6d6fb..8225cd6 100644 --- a/nimbleModel/inst/include/nimbleModel/predef/modelBase_nC/modelBase_nClass_manifest.txt +++ b/nimbleModel/inst/include/nimbleModel/predef/modelBase_nC/modelBase_nClass_manifest.txt @@ -1,4 +1,4 @@ -list(saved_at = structure(1793056441.24624, class = c("POSIXct", +list(saved_at = structure(1793832608.695, class = c("POSIXct", "POSIXt")), packet_name = "modelBase_nClass", elements = c("preamble", "cppContent", "hContent", "filebase", "post_cpp_compiler", "copyFiles", "depends"), files = list(preamble = "modelBase_nClass_preamble.cpp", diff --git a/nimbleModel/inst/include/nimbleModel/predef/modelClass_/modelClass_.h b/nimbleModel/inst/include/nimbleModel/predef/modelClass_/modelClass_.h index add2a38..79b6446 100644 --- a/nimbleModel/inst/include/nimbleModel/predef/modelClass_/modelClass_.h +++ b/nimbleModel/inst/include/nimbleModel/predef/modelClass_/modelClass_.h @@ -9,26 +9,26 @@ class modelClass_ : public modelBase_nClass { modelClass_() {}; std::vector< std::shared_ptr > declFunPtrs; std::map name2index_map; - double calculate_impl(std::shared_ptr instrList) override { + double calculate_impl(std::shared_ptr instrList) override { // double logProb(0.0); - // const auto& instrVec = instrList->contents(); + // const auto& instrVec = instrList->instr_nList->contents(); // for (const auto& instr : instrVec) { // logProb += declFunPtrs[instr->declID - 1 ]->calculate_cpp(instr); // } // return(logProb); return calc_op_impl< &declFunBase_nClass::calculate_cpp >(instrList); } - double calculateDiff_impl(std::shared_ptr instrList) override { + double calculateDiff_impl(std::shared_ptr instrList) override { return calc_op_impl< &declFunBase_nClass::calculateDiff_cpp >(instrList); } - double getLogProb_impl(std::shared_ptr instrList) override { + double getLogProb_impl(std::shared_ptr instrList) override { return calc_op_impl< &declFunBase_nClass::getLogProb_cpp >(instrList); } template - double calc_op_impl(std::shared_ptr instrList) { + double calc_op_impl(std::shared_ptr instrList) { double logProb(0.0); - const auto& instrVec = instrList->contents(); + const auto& instrVec = instrList->instr_nList->contents(); for (const auto& instr : instrVec) { logProb += ((*(declFunPtrs[instr->declID - 1 ])).*Method)(instr); } @@ -49,13 +49,96 @@ class modelClass_ : public modelBase_nClass { return declFunPtrs[node->declID - 1]->getBound_cpp(node, boundID); } - void simulate_impl(std::shared_ptr instrList) { - const auto& instrVec = instrList->contents(); + void simulate_impl(std::shared_ptr instrList) { + const auto& instrVec = instrList->instr_nList->contents(); for (const auto& instr : instrVec) { declFunPtrs[instr->declID - 1 ]->simulate_cpp(instr); } } + // One node (list step) of an instruction list, with no loop. `step` (1-based) numbers the + // nodes of the whole list in calculation order. Instruction j holds list steps + // offsets[j] + 1 through offsets[j] + nNodes, so a binary search finds the last j with + // offsets[j] < step (skipping any instructions with no nodes). Returns that instruction + // and sets `k` to the step within it (1-based). + std::shared_ptr find_instr_for_step_(const instrList_nClass& instrList, int step, int& k) { + if(step < 1 || step > instrList.nNodes_total) + Rcpp::stop("step %d is outside the instruction list's %d steps", step, instrList.nNodes_total); + const auto& offsets = instrList.offsets; + int lo = 0, hi = offsets.size() - 1; + while(lo < hi) { + int mid = (lo + hi + 1) / 2; + if(offsets[mid] < step) lo = mid; else hi = mid - 1; + } + k = step - offsets[lo]; + return instrList.instr_nList->contents()[lo]; + } + template + double calc_indexed_op_impl(std::shared_ptr instrList, int step) { + int k; + auto instr = find_instr_for_step_(*instrList, step, k); + return ((*(declFunPtrs[instr->declID - 1])).*Method)(instr, k); + } + double calculate_indexed_impl(std::shared_ptr instrList, int step) override { + return calc_indexed_op_impl< &declFunBase_nClass::calculate_singleton_cpp >(instrList, step); + } + double calculateDiff_indexed_impl(std::shared_ptr instrList, int step) override { + return calc_indexed_op_impl< &declFunBase_nClass::calculateDiff_singleton_cpp >(instrList, step); + } + double getLogProb_indexed_impl(std::shared_ptr instrList, int step) override { + return calc_indexed_op_impl< &declFunBase_nClass::getLogProb_singleton_cpp >(instrList, step); + } + void simulate_indexed_impl(std::shared_ptr instrList, int step) override { + int k; + auto instr = find_instr_for_step_(*instrList, step, k); + declFunPtrs[instr->declID - 1]->simulate_singleton_cpp(instr, k); + } + std::unique_ptr getParam_indexed_impl(std::shared_ptr instrList, int step, int paramID) override { + int k; + auto instr = find_instr_for_step_(*instrList, step, k); + return declFunPtrs[instr->declID - 1]->getParam_singleton_cpp(instr, k, paramID); + } + double getBound_indexed_impl(std::shared_ptr instrList, int step, int boundID) override { + int k; + auto instr = find_instr_for_step_(*instrList, step, k); + return declFunPtrs[instr->declID - 1]->getBound_singleton_cpp(instr, k, boundID); + } + Eigen::Tensor getValues_indexed_impl(std::shared_ptr instrList, int step) override { + int k; + auto instr = find_instr_for_step_(*instrList, step, k); + return declFunPtrs[instr->declID - 1]->getValues_singleton_cpp(instr, k); + } + void setValues_indexed_impl(std::shared_ptr instrList, int step, Eigen::Tensor value) override { + int k; + auto instr = find_instr_for_step_(*instrList, step, k); + declFunPtrs[instr->declID - 1]->setValues_singleton_cpp(instr, k, value); + } + // All the nodes of an instruction list, for values(model, nodes, mode = "node"). + // The total length is found first (cheaply: no values are copied), so the result is + // allocated once and a wrong-length value is an error before anything is set. + // There is one virtual call per instruction (not per node) for each pass. + int valuesLength_impl_(const instrList_nClass &instrList) { + int total = 0; + for (const auto& instr : instrList.instr_nList->contents()) + total += declFunPtrs[instr->declID - 1]->valuesLength_cpp(instr); + return total; + } + Eigen::Tensor getValues_nodes_impl(std::shared_ptr instrList) override { + Eigen::Tensor out(valuesLength_impl_(*instrList)); + int offset = 0; + for (const auto& instr : instrList->instr_nList->contents()) + offset = declFunPtrs[instr->declID - 1]->getValues_into_cpp(instr, out, offset); + return out; + } + void setValues_nodes_impl(std::shared_ptr instrList, Eigen::Tensor value) override { + const int total = valuesLength_impl_(*instrList); + if (value.size() != total) + Rcpp::stop("setValues: value has length %d but the nodes have %d values.", (int)value.size(), total); + int offset = 0; + for (const auto& instr : instrList->instr_nList->contents()) + offset = declFunPtrs[instr->declID - 1]->setValues_from_cpp(instr, value, offset); + } + // This version takes a character vector of names from R so that // the ordering of nodeFxns matches that in R, which is important for // the calculation instructions. @@ -186,3 +269,48 @@ class modelClass_ : public modelBase_nClass { } } }; + +// Support for nimble2's values(model, nodes[i], mode = "node") <- y, which becomes +// setValues_indexed_(model, instrList, i) = y. When generating code, +// we know we are on the LHS, and it is easiest to modify the LHS alone. +// Using a proxy object with an operator=() makes that simple. +struct SetValuesIndexedProxy_ { + std::shared_ptr model; + std::shared_ptr instrList; + int step; + void operator=(const Eigen::Tensor &v) { + model->setValues_indexed_impl(instrList, step, v); + } + void operator=(double x) { + Eigen::Tensor v(1); + v(0) = x; + model->setValues_indexed_impl(instrList, step, v); + } +}; + +template +SetValuesIndexedProxy_ setValues_indexed_(const std::shared_ptr &model, + std::shared_ptr instrList, int step) { + return SetValuesIndexedProxy_{model, instrList, step}; +} + +// The same for values(model, nodes, mode = "node") <- y, which becomes +// setValues_nodes_(model, instrList) = y. +struct SetValuesNodesProxy_ { + std::shared_ptr model; + std::shared_ptr instrList; + void operator=(const Eigen::Tensor &v) { + model->setValues_nodes_impl(instrList, v); + } + void operator=(double x) { + Eigen::Tensor v(1); + v(0) = x; + model->setValues_nodes_impl(instrList, v); + } +}; + +template +SetValuesNodesProxy_ setValues_nodes_(const std::shared_ptr &model, + std::shared_ptr instrList) { + return SetValuesNodesProxy_{model, instrList}; +} diff --git a/nimbleModel/tests/testthat/test-copier.R b/nimbleModel/tests/testthat/test-copier.R index e45f34f..75cf5a2 100644 --- a/nimbleModel/tests/testthat/test-copier.R +++ b/nimbleModel/tests/testthat/test-copier.R @@ -47,7 +47,7 @@ test_that("multiCopier basics works", { function(v = "numericVector") { # Work is needed on opDef dispatch to support # some assignment using multCopy in the syntax via an opDef. - cppLiteral("multCopy->flatViewGroup.setValues_() = v;") + cppLiteral("multCopy->setValues_() = v;") } ) ) @@ -73,3 +73,141 @@ test_that("multiCopier basics works", { check2 <- c(m$y[2:4], as.numeric(m$z[2:4, 3:5])) expect_equal(val2, check2) }) + +# values(model, nodes[i], mode = "raw" or "scalar"): one entry of nodes, or one flattened element. +test_that("multiCopier getValues_indexed and setValues_indexed work, compiled", { + code <- quote({ + tau ~ dunif(0, 100) + mu ~ dnorm(0, 1) + for (i in 1:5) { + y[i] ~ dnorm(mu, var = tau) + } + for (i in 1:5) { + for (j in 1:5) { + z[i, j] ~ dnorm(mu, var = tau) + } + } + }) + inits <- list(tau = 25, mu = 0, z = matrix(rnorm(25), nrow = 5)) + data <- list(y = rnorm(5)) + mclass <- nimbleModel(code, inits = inits, data = data, returnClass = TRUE) + + nf <- nClass( + Cpublic = list( + m = "nimbleModel:::modelBase_nClass()", + multCopy = "nimbleModel:::multiCopier_nClass()", + init = nFunction(function() { + multCopy$init(m) + }), + get = nFunction(function() { + return(multCopy$getValues()) + }, returnType = "numericVector"), + get_i = nFunction(function(i = "integerScalar", mode = "string") { + return(multCopy$getValues_indexed(i, mode)) + }, returnType = "numericVector"), + set_i = nFunction(function(i = "integerScalar", v = "numericVector", mode = "string") { + multCopy$setValues_indexed(i, v, mode) + }), + set_raw_scalar = nFunction(function(i = "integerScalar", x = "numericScalar") { + # A scalar into a length-1 entry, as values(model, nodes[i]) <- x would do. + cppLiteral("multCopy->setValues_raw_(i) = x;") + }) + ) + ) + comp <- nCompile(mclass, nf) + m <- comp$mclass$new() + set.seed(1) + m$simulate() + m$calculate() + obj <- comp$nf$new() + obj$m <- m + # Includes a whole scalar variable, a duplicated (overlapping) entry, and a logProb variable. + nodes <- c("y[2:4]", "z[2:4, 3:5]", "mu", "y[3]", "logProb_y[1:2]") + obj$multCopy <- makeMultiCopier(m, nodes) + obj$init() + + entries <- function() { + list(m$y[2:4], as.numeric(m$z[2:4, 3:5]), m$mu, m$y[3], m$logProb_y[1:2]) + } + vals <- entries() + for (k in seq_along(nodes)) { + expect_equal(obj$get_i(k, "raw"), vals[[k]], info = nodes[k]) + } + flat <- unlist(vals) + expect_equal(obj$get(), flat) + for (k in seq_along(flat)) { + expect_equal(obj$get_i(k, "scalar"), flat[k], info = k) + } + + # raw set, including into the duplicated entry + new_z <- rnorm(9) + obj$set_i(2L, new_z, "raw") + expect_equal(as.numeric(m$z[2:4, 3:5]), new_z) + obj$set_i(4L, 10.5, "raw") + expect_equal(m$y[3], 10.5) + obj$set_raw_scalar(3L, -2.5) + expect_equal(m$mu, -2.5) + + # scalar set: element 5 is z[3, 3]; element 13 is mu + obj$set_i(5L, 7.25, "scalar") + expect_equal(m$z[3, 3], 7.25) + obj$set_i(13L, 3.5, "scalar") + expect_equal(m$mu, 3.5) + + # errors for bad indices, sizes, and modes + expect_error(obj$get_i(0L, "raw")) + expect_error(obj$get_i(6L, "raw")) + expect_error(obj$get_i(length(flat) + 1L, "scalar")) + expect_error(obj$set_i(1L, c(1, 2), "raw")) + expect_error(obj$set_i(1L, c(1, 2), "scalar")) + expect_error(obj$set_raw_scalar(1L, 1)) + expect_error(obj$get_i(1L, "node")) +}) + +test_that("uncompiled values(model, nodes[i], mode = ...) and values_indexed work", { + code <- quote({ + mu ~ dnorm(0, 1) + for (i in 1:5) { + y[i] ~ dnorm(mu, 1) + } + for (i in 1:5) { + for (j in 1:5) { + z[i, j] ~ dnorm(mu, 1) + } + } + }) + m <- nimbleModel(code, inits = list(mu = 1, y = rnorm(5), z = matrix(rnorm(25), nrow = 5))) + nodes <- c("y[2:4]", "z[2:4, 3:5]", "mu") + flat <- c(m$y[2:4], as.numeric(m$z[2:4, 3:5]), m$mu) + + # raw is the default and is ordinary R indexing of nodes + expect_equal(values(m, nodes[2]), as.numeric(m$z[2:4, 3:5])) + expect_equal(values_indexed(m, nodes, 2, "raw"), as.numeric(m$z[2:4, 3:5])) + # nodes as a list of varRanges + vr <- lapply(nodes, \(x) varRangeClass$new(x)) + expect_equal(values(m, vr[3]), m$mu) + for (k in seq_along(flat)) { + expect_equal(values(m, nodes[k], mode = "scalar"), flat[k], info = k) + } + i <- 4 + expect_equal(values(m, nodes[i + 1], mode = "scalar"), flat[5]) + + values(m, nodes[3]) <- 2.5 # whole-variable entry + expect_equal(m$mu, 2.5) + values(m, nodes[1], mode = "raw") <- c(1, 2, 3) + expect_equal(m$y[2:4], c(1, 2, 3)) + values(m, nodes[5], mode = "scalar") <- 7.25 # z[3, 3] + expect_equal(m$z[3, 3], 7.25) + values_indexed(m, vr, 13, "scalar") <- -1 + expect_equal(m$mu, -1) + + # unindexed scalar mode is raw mode + expect_equal(values(m, nodes, mode = "scalar"), values(m, nodes)) + allVals <- values(m, nodes) + values(m, nodes, mode = "scalar") <- allVals + 1 + expect_equal(values(m, nodes), allVals + 1) + values(m, nodes) <- allVals + expect_error(values(m, nodes[length(flat) + 1], mode = "scalar")) + expect_error(values_indexed(m, nodes, 1, "raw") <- 1) + expect_equal(values_indexed(m, nodes, 1, "node"), m$mu) # node mode uses sorted order +}) diff --git a/nimbleModel/tests/testthat/test-getParam_getBound.R b/nimbleModel/tests/testthat/test-getParam_getBound.R index 686cf37..73a22db 100644 --- a/nimbleModel/tests/testthat/test-getParam_getBound.R +++ b/nimbleModel/tests/testthat/test-getParam_getBound.R @@ -36,10 +36,10 @@ test_that("basic getParam and getBound works from full interface", { # Now that is supported (and tested) inside the functions # vr <- varRangeClass$new('y[3]') # instrList <- makeInstrList(m, vr) - # y3 <- instrList[[1]] + # y3 <- instrList$instr_nList[[1]] # vr <- varRangeClass$new("v") # instrList <- makeInstrList(m, vr) - # v <- instrList[[1]] + # v <- instrList$instr_nList[[1]] y3 <- "y[3]" v <- "v" for (case in c("comp", "uncomp")) { diff --git a/nimbleModel/tests/testthat/test-instrSteps.R b/nimbleModel/tests/testthat/test-instrSteps.R new file mode 100644 index 0000000..c03c229 --- /dev/null +++ b/nimbleModel/tests/testthat/test-instrSteps.R @@ -0,0 +1,251 @@ +library(nimbleModel) +library(testthat) +library(nCompiler) + + +# Tests of numbering the nodes (steps) of instructions and instruction lists: +# `set_instr_idx` must give, for every step k, the k-th idx visited by the full runners, +# and the one-step (`*_indexed`, and through them `*_singleton`) entry points must act on exactly that node. + +# Build an instruction (in its R-list form) directly from its indexing. +# values: start for a seq range; for a mat/matp range, the counter values, each counter's dims[r] values together. +make_test_instr <- function(index_types, dims, lens, values, slots = seq_len(sum(dims))) { + instr <- list( + lens = lens, index_types = index_types, nDim = sum(dims), dims = dims, + slots = as.numeric(slots), values = values, declID = 1, sortID = 1 + ) + instr <- nimbleModel:::add_instr_precalcs(instr) + class(instr) <- "Rinstr" + instr +} + +# Run the full runner for `instr` with a per-node function that records each idx it receives. +record_runner_idx <- function(instr) { + seen <- list() + nimbleModel:::run_instr(instr, function(idx) { + seen[[length(seen) + 1]] <<- as.numeric(idx) + 0 + }) + seen +} + +test_instrs <- list( + "0" = list(lens = 1, index_types = 0, nDim = 0) |> nimbleModel:::add_instr_precalcs() |> structure(class = "Rinstr"), + "1" = make_test_instr(1, 1, 4, list(3)), + "2" = make_test_instr(2, 1, 3, list(c(5, 2, 7))), + "3" = make_test_instr(2, 2, 3, list(c(1, 4, 2, 5, 3, 6))), + "18" = make_test_instr(2, 2, 3, list(c(1, 4, 2, 5, 3, 6)), slots = c(2, 1)), + "4" = make_test_instr(c(1, 1), c(1, 1), c(3, 4), list(2, 5)), + "5" = make_test_instr(c(1, 2), c(1, 1), c(3, 2), list(2, c(9, 4))), + "6" = make_test_instr(c(2, 1), c(1, 1), c(2, 3), list(c(9, 4), 2)), + "7" = make_test_instr(c(2, 2), c(1, 1), c(2, 3), list(c(9, 4), c(1, 8, 3))), + "8" = make_test_instr(c(1, 2), c(1, 2), c(3, 2), list(2, c(1, 4, 2, 5))), + "9" = make_test_instr(c(2, 1), c(2, 1), c(2, 3), list(c(1, 4, 2, 5), 2)), + "10" = make_test_instr(c(2, 2), c(2, 2), c(2, 2), list(c(1, 4, 2, 5), c(7, 8, 9, 6)), slots = c(3, 1, 4, 2)), + "11" = make_test_instr(c(2, 1), c(1, 1), c(2, 3), list(c(9, 4), 2), slots = c(2, 1)), + "12" = make_test_instr(c(1, 1, 1), c(1, 1, 1), c(2, 3, 2), list(1, 4, 2)), + "13" = make_test_instr(c(1, 2, 2), c(1, 2, 1), c(2, 2, 3), list(3, c(1, 4, 2, 5), c(7, 9, 8)), slots = c(3, 1, 4, 2)), + "14" = make_test_instr(rep(1, 4), rep(1, 4), c(2, 2, 3, 2), list(1, 2, 3, 4)), + "15" = make_test_instr(c(2, 1, 2, 1), c(1, 1, 2, 1), c(2, 2, 2, 3), list(c(5, 3), 2, c(1, 4, 2, 5), 6), slots = c(5, 2, 1, 3, 4)), + "16" = make_test_instr(rep(1, 5), rep(1, 5), c(2, 2, 2, 2, 2), list(1, 2, 3, 4, 5)), + "17" = make_test_instr(c(1, 2, 1, 2, 1), c(1, 1, 1, 2, 1), c(2, 2, 2, 2, 3), list(1, c(6, 3), 2, c(1, 4, 2, 5), 7), slots = c(6, 4, 5, 1, 2, 3)) +) + +test_that("test instructions cover every instruction type", { + types <- sapply(test_instrs, \(x) x$instr_type) + expect_identical(as.character(types), names(test_instrs)) + expect_setequal(types, unlist(nimbleModel:::type2itype) |> unname()) +}) + +test_that("set_instr_idx gives the k-th idx of the full runner for every instruction type", { + for (type in names(test_instrs)) { + instr <- test_instrs[[type]] + full <- record_runner_idx(instr) + expect_identical(length(full), as.integer(instr$nNodes), info = paste("type", type)) + for (k in seq_along(full)) { + expect_identical(nimbleModel:::set_instr_idx(instr, k), full[[k]], info = paste("type", type, "step", k)) + } + } +}) + +test_that("list steps map to the right instruction and step, including empty instructions", { + empty <- make_test_instr(1, 1, 0, list(1)) + instrs <- list(test_instrs[["4"]], empty, test_instrs[["1"]], empty, empty, test_instrs[["2"]]) + instrList <- nimbleModel:::make_RinstrList(instrs) + expect_equal(instrList$offsets, c(0, 12, 12, 16, 16, 16)) + expect_equal(instrList$nNodes_total, 19) + expected <- list() + for (j in seq_along(instrs)) { + for (k in seq_len(instrs[[j]]$nNodes)) expected[[length(expected) + 1]] <- list(j = j, k = k) + } + for (step in 1:19) { + expect_equal(nimbleModel:::instrList_find_singleton(instrList, step), expected[[step]], info = paste("step", step)) + loc <- expected[[step]] + expect_identical( + nimbleModel:::set_instr_idx_from_instrList(instrList, step), + nimbleModel:::set_instr_idx(instrs[[loc$j]], loc$k) + ) + } + expect_error(nimbleModel:::instrList_find_singleton(instrList, 0), "outside") + expect_error(nimbleModel:::instrList_find_singleton(instrList, 20), "outside") +}) + +# For each step of an instruction list, check that the one-step methods act on exactly the +# node of `z` that `set_instr_idx` names, in model `m` (compiled or not). +check_indexed_steps <- function(m, nodes) { + instrList <- makeInstrList(m, nodes) + n <- instrList$nNodes_total + for (step in seq_len(n)) { + idx <- nimbleModel:::set_instr_idx_from_instrList(instrList, step) + m$z <- array(0, dim(m$z)) + m$calculate_indexed(instrList, step) + changed <- which(m$z != 0, arr.ind = TRUE) + expect_identical(nrow(changed), 1L, info = paste("step", step)) + expect_equal(as.numeric(changed), idx, info = paste("step", step)) + # idx holds 1-based (R) index values, used both to index z and in z's expression. + expect_equal(m$z[changed], idx[1] * 10 + idx[2], info = paste("step", step)) + } + # Running every step in turn gives the same result as a full calculate. + m$z <- array(0, dim(m$z)) + m$calculate(instrList) + full <- m$z + m$z <- array(0, dim(m$z)) + for (step in seq_len(n)) m$calculate_indexed(instrList, step) + expect_equal(m$z, full) +} + +test_that("calculate_indexed acts on the node set_instr_idx names, uncompiled and compiled", { + code <- nimbleCode({ + for (i in 1:4) { + for (j in 1:3) { + z[i, j] <- i * 10 + j # ensure that i and j are used in R-based indexing on RHS, even when compiled. + } + } + }) + m <- nimbleModel(code) + cm <- nimbleModel(code, compile = TRUE) + for (model in list(m, cm)) { + check_indexed_steps(model, "z") # 2_seq_seq + check_indexed_steps(model, "z[2:3, 2]") + check_indexed_steps(model, "z[c(4, 1), 2:3]") # mat outer range + check_indexed_steps(model, "z[2:4, c(3, 1)]") # mat inner range + } + expect_error(cm$calculate_indexed("z", 13), "outside") +}) + +test_that("getLogProb_indexed, getParam_indexed, getBound_indexed and simulate_indexed act on the k-th node", { + code <- nimbleCode({ + for (i in 1:3) { + for (j in 1:2) { + y[i, j] ~ dunif(lo[i, j], lo[i, j] + w[i, j]) + } + } + }) + lo <- matrix(1:6, 3, 2) + w <- matrix(c(2, 4, 5, 8, 10, 20), 3, 2) + y <- lo + 1 + m <- nimbleModel(code, data = list(lo = lo, w = w), inits = list(y = y)) + cm <- nimbleModel(code, data = list(lo = lo, w = w), inits = list(y = y), compile = TRUE) + for (model in list(m, cm)) { + model$calculate() # also calculates the lifted node for lo + w + instrList <- makeInstrList(model, "y") + for (step in seq_len(instrList$nNodes_total)) { + idx <- nimbleModel:::set_instr_idx_from_instrList(instrList, step) + i <- idx[1] + j <- idx[2] + expect_equal(model$getParam_indexed(instrList, step, "min"), lo[i, j]) + expect_equal(model$getBound_indexed(instrList, step, "upper"), lo[i, j] + w[i, j]) + expect_equal(model$getLogProb_indexed(instrList, step), -log(w[i, j])) + before <- model$y + model$simulate_indexed(instrList, step) + changed <- which(model$y != before, arr.ind = TRUE) + expect_identical(nrow(changed), 1L) + expect_equal(as.numeric(changed), c(i, j)) + } + } +}) + +# Until compile-time handling maps `model$calculate_indexed(nodes, i)` and its siblings to the +# `*_indexed_impl` methods, compiled code calls those methods directly. This checks that the new +# virtual methods and `instrList_nClass` work from generated C++. (The uncompiled `*_impl` +# methods are stubs, so this is compiled only.) +test_that("compiled code can call the *_indexed_impl methods with an instrList_nClass", { + code <- nimbleCode({ + for (i in 1:3) { + for (j in 1:2) { + y[i, j] ~ dunif(lo[i, j], lo[i, j] + w[i, j]) + } + } + }) + lo <- matrix(1:6, 3, 2) + w <- matrix(c(2, 4, 5, 8, 10, 20), 3, 2) + y <- lo + 1 + mclass <- nimbleModel(code, data = list(lo = lo, w = w), inits = list(y = y), returnClass = TRUE) + indexedCaller <- nClass( + classname = "indexedCaller", + Cpublic = list( + model = "nimbleModel:::modelBase_nClass()", + nodes = "nimbleModel:::instrList_nClass()", + calc = nFunction( + function(i = "integerScalar") { + return(model$calculate_indexed_impl(nodes, i)) + }, + returnType = "numericScalar" + ), + calcDiff = nFunction( + function(i = "integerScalar") { + return(model$calculateDiff_indexed_impl(nodes, i)) + }, + returnType = "numericScalar" + ), + logProb = nFunction( + function(i = "integerScalar") { + return(model$getLogProb_indexed_impl(nodes, i)) + }, + returnType = "numericScalar" + ), + sim = nFunction( + function(i = "integerScalar") { + model$simulate_indexed_impl(nodes, i) + }, + returnType = "void" + ), + param = nFunction( + function(i = "integerScalar", paramID = "integerScalar") { + return(model$getParam_indexed_impl(nodes, i, paramID)) + }, + returnType = "ETaccessor()" + ), + bound = nFunction( + function(i = "integerScalar", boundID = "integerScalar") { + return(model$getBound_indexed_impl(nodes, i, boundID)) + }, + returnType = "numericScalar" + ) + ) + ) + comp <- nCompile(mclass, indexedCaller) + cm <- comp$mclass$new() + caller <- comp$indexedCaller$new() + caller$model <- cm + instrList <- makeInstrList(cm, "y") + caller$nodes <- cm$makeCompiledInstrList(instrList) + minID <- cm$getParamID(instrList, "min") + cm$calculate() + for (step in seq_len(instrList$nNodes_total)) { + idx <- nimbleModel:::set_instr_idx_from_instrList(instrList, step) + i <- idx[1] + j <- idx[2] + expect_equal(caller$calc(step), -log(w[i, j])) + expect_equal(caller$logProb(step), -log(w[i, j])) + expect_equal(caller$calcDiff(step), 0) + expect_equal(caller$param(step, minID), lo[i, j]) + expect_equal(caller$bound(step, 2L), lo[i, j] + w[i, j]) # boundID 2 is the upper bound + before <- cm$y + caller$sim(step) + changed <- which(cm$y != before, arr.ind = TRUE) + expect_identical(nrow(changed), 1L) + expect_equal(as.numeric(changed), c(i, j)) + } + expect_error(caller$calc(7L), "outside") +}) diff --git a/nimbleModel/tests/testthat/test-nimbleModel.R b/nimbleModel/tests/testthat/test-nimbleModel.R index a349030..8c77dc0 100644 --- a/nimbleModel/tests/testthat/test-nimbleModel.R +++ b/nimbleModel/tests/testthat/test-nimbleModel.R @@ -117,7 +117,7 @@ test_that("basic testing of models, compiled and uncompiled", { ## Check that instrList is in correct order. instrList <- makeInstrList(m, c("y", "lifted_sqrt_oPtau_cP")) - expect_equal(instrList[[1]]$lens, 1) # lifted node first + expect_equal(instrList$instr_nList[[1]]$lens, 1) # lifted node first lp <- m$calculate(instrList) expect_equal(m$lifted_sqrt_oPtau_cP, 5) expect_equal(lp, lp_y) @@ -221,6 +221,163 @@ test_that("basic testing of models, compiled and uncompiled", { gc() }) +test_that("model methods called from nClass methods, uncompiled and compiled", { + # Each model method is called from its own method of one nClass, so that there is one compilation. + # In compiled code, model$calculate(nodes) and so on become model$calculate_impl(nodes) and so on + # (see update_model_calc_methods_LAT in modelBaseClass.R). + code <- nimbleCode({ + mu ~ dnorm(0, sd = 2) + for (i in 1:4) { + y[i] ~ dnorm(mu, sd = 1.5) + } + for (i in 1:3) { + u[i] ~ dunif(lo[i], lo[i] + w[i]) + } + }) + y <- c(0.3, -1.2, 2.1, 0.8) + lo <- c(-1, 0, 2) + w <- c(1, 2, 3) + mclass <- nimbleModel(code, + inits = list(mu = 0.5, lo = lo, w = w, u = lo + w / 2), + data = list(y = y), returnClass = TRUE + ) + modelCaller <- nClass( + classname = "modelCaller", + Cpublic = list( + model = "nimbleModel:::modelBase_nClass()", + ynodes = "nimbleModel:::instrList_nClass()", # y[1:4], data + unodes = "nimbleModel:::instrList_nClass()", # u[1:3] + unode = "nimbleModel:::instr_nClass()", # the instruction for u[1:3]; getParam and getBound use u[1] + calc = nFunction(function() { + return(model$calculate(ynodes)) + }, returnType = "numericScalar"), + calcDiff = nFunction(function() { + return(model$calculateDiff(ynodes)) + }, returnType = "numericScalar"), + logProb = nFunction(function() { + return(model$getLogProb(ynodes)) + }, returnType = "numericScalar"), + sim = nFunction(function() { + model$simulate(unodes) + }, returnType = "void"), + paramID = nFunction(function(param = "character()") { + return(model$getParamID(unode, param)) + }, returnType = "integerScalar"), + param = nFunction(function(paramID = "integerScalar") { + return(model$getParam(unode, paramID)) + }, returnType = "ETaccessor()"), + bound = nFunction(function(boundID = "integerScalar") { + return(model$getBound(unode, boundID)) + }, returnType = "numericScalar"), + calc_indexed = nFunction(function(i = "integerScalar") { + return(model$calculate_indexed(ynodes, i)) + }, returnType = "numericScalar"), + calcDiff_indexed = nFunction(function(i = "integerScalar") { + return(model$calculateDiff_indexed(ynodes, i)) + }, returnType = "numericScalar"), + logProb_indexed = nFunction(function(i = "integerScalar") { + return(model$getLogProb_indexed(ynodes, i)) + }, returnType = "numericScalar"), + sim_indexed = nFunction(function(i = "integerScalar") { + model$simulate_indexed(unodes, i) + }, returnType = "void"), + param_indexed = nFunction(function(i = "integerScalar", paramID = "integerScalar") { + return(model$getParam_indexed(unodes, i, paramID)) + }, returnType = "ETaccessor()"), + # Named arguments out of order, which the matchDef in the opDef must put in order. + param_indexed_named = nFunction(function(i = "integerScalar", paramID = "integerScalar") { + return(model$getParam_indexed(paramID = paramID, i = i, instrList = unodes)) + }, returnType = "ETaccessor()"), + bound_indexed = nFunction(function(i = "integerScalar", boundID = "integerScalar") { + return(model$getBound_indexed(unodes, i, boundID)) + }, returnType = "numericScalar") + ) + ) + comp <- nCompile(mclass, modelCaller) + + for (case in c("uncompiled", "compiled")) { + if (case == "uncompiled") { + m <- mclass$new() + caller <- modelCaller$new() + } else { + m <- comp$mclass$new() + caller <- comp$modelCaller$new() + } + caller$model <- m + ynodes <- makeInstrList(m, "y") + unodes <- makeInstrList(m, "u") + if (case == "uncompiled") { + caller$ynodes <- ynodes + caller$unodes <- unodes + caller$unode <- unodes$instr_nList[[1]] + } else { + caller$ynodes <- m$makeCompiledInstrList(ynodes) + caller$unodes <- m$makeCompiledInstrList(unodes) + caller$unode <- m$makeCompiledInstr(unodes$instr_nList[[1]]) + } + m$calculate() # also calculates the lifted nodes for lo + w + + # calculate, getLogProb and calculateDiff + lp_y <- dnorm(y, 0.5, 1.5, log = TRUE) + expect_equal(caller$calc(), sum(lp_y), info = case) + expect_equal(caller$logProb(), sum(lp_y), info = case) + expect_equal(caller$calcDiff(), 0, info = case) + m$mu <- 1 + new_lp_y <- dnorm(y, 1, 1.5, log = TRUE) + expect_equal(caller$calcDiff(), sum(new_lp_y) - sum(lp_y), info = case) + expect_equal(caller$logProb(), sum(new_lp_y), info = case) + for (k in 1:4) { + expect_equal(caller$logProb_indexed(k), new_lp_y[k], info = case) + } + m$mu <- 0.5 + for (k in 1:4) { + expect_equal(caller$calcDiff_indexed(k), lp_y[k] - new_lp_y[k], info = case) + expect_equal(caller$calc_indexed(k), lp_y[k], info = case) + } + expect_equal(caller$logProb(), sum(lp_y), info = case) + + # getParamID, getParam and getBound: the non-indexed versions use the first node, u[1] + minID <- caller$paramID("min") + expect_equal(minID, m$getParamID("u", "min"), info = case) + expect_warning(expect_equal(caller$param(minID), lo[1], info = case)) + expect_warning(expect_equal(caller$bound(1L), lo[1], info = case)) # boundID 1 is the lower bound + expect_warning(expect_equal(caller$bound(2L), lo[1] + w[1], info = case)) # boundID 2 is the upper bound + for (k in 1:3) { + expect_equal(caller$param_indexed(k, minID), lo[k], info = case) + expect_equal(caller$param_indexed_named(k, minID), lo[k], info = case) + expect_equal(caller$bound_indexed(k, 1L), lo[k], info = case) + expect_equal(caller$bound_indexed(k, 2L), lo[k] + w[k], info = case) + } + + # simulate_indexed changes only the k-th node, within its bounds + set.seed(1) + for (k in 1:3) { + before <- m$u + caller$sim_indexed(k) + expect_identical(which(m$u != before), k, info = case) + expect_true(m$u[k] >= lo[k] && m$u[k] <= lo[k] + w[k], info = case) + } + if(case == "uncompiled") { + sim_u <- m$u # record uncompiled + } else { + expect_equal(sim_u, m$u) # then check for compiled match + } + # simulate changes all of them, within their bounds + set.seed(2) + before <- m$u + caller$sim() + expect_true(all(m$u != before), info = case) + expect_true(all(m$u >= lo & m$u <= lo + w), info = case) + if(case == "uncompiled") { + sim_u2 <- m$u # record uncompiled + } else { + expect_equal(sim_u2, m$u) # then check for compiled match + } + # and the data nodes were not touched + expect_identical(m$y, y, info = case) + } +}) + test_that("no indices", { code <- quote({ mu ~ dnorm(0, 1) @@ -917,7 +1074,7 @@ test_that("calculate works correctly for time series/SSM recursion", { m <- mclass$new() instrList <- makeInstrList(m, c("y", "mu")) - sortIDs <- lapply(instrList, \(x) x$sortID) + sortIDs <- lapply(instrList$instr_nList, \(x) x$sortID) expect_true(all(diff(sapply(sortIDs, \(x) min(x, na.rm = TRUE))) >= 0), "incorrect order") expect_true( all(sapply(sortIDs, \(x) length(x) == 1 || all(diff(x) >= 1, na.rm = TRUE))), @@ -925,8 +1082,8 @@ test_that("calculate works correctly for time series/SSM recursion", { ) # Check proper handling when providing list of instructions. - instrList <- makeInstrList(m, instrList[c(5, 3, 1, 2, 4)]) - sortIDs <- lapply(instrList, \(x) x$sortID) + instrList <- makeInstrList(m, instrList$instr_nList[c(5, 3, 1, 2, 4)]) + sortIDs <- lapply(instrList$instr_nList, \(x) x$sortID) expect_true(all(diff(sapply(sortIDs, \(x) min(x, na.rm = TRUE))) >= 0)) expect_true(all(sapply(sortIDs, \(x) length(x) == 1 || all(diff(x) >= 1, na.rm = TRUE)))) @@ -958,7 +1115,7 @@ test_that("calculate works correctly for time series/SSM recursion", { m <- mclass$new() instrList <- makeInstrList(m, c("y", "mu")) - sortIDs <- lapply(instrList, \(x) x$sortID) + sortIDs <- lapply(instrList$instr_nList, \(x) x$sortID) expect_true(all(diff(sapply(sortIDs, \(x) min(x, na.rm = TRUE))) >= 0), "incorrect order") expect_true( all(sapply(sortIDs, \(x) length(x) == 1 || all(diff(x) >= 1, na.rm = TRUE))), @@ -967,7 +1124,7 @@ test_that("calculate works correctly for time series/SSM recursion", { # Check idempotency instrList <- makeInstrList(m, instrList) - sortIDs <- lapply(instrList, \(x) x$sortID) + sortIDs <- lapply(instrList$instr_nList, \(x) x$sortID) expect_true(all(diff(sapply(sortIDs, \(x) min(x, na.rm = TRUE))) >= 0), "incorrect order") expect_true( all(sapply(sortIDs, \(x) length(x) == 1 || all(diff(x) >= 1, na.rm = TRUE))), @@ -976,8 +1133,8 @@ test_that("calculate works correctly for time series/SSM recursion", { # Check proper handling when providing list of instructions. # This would need to be created as an R list of instrClass objects, - # instrList <- makeInstrList(m, lapply(c(5,3,1,2,4), \(i) instrList[[i]])) - # sortIDs <- lapply(instrList, \(x) x$sortID) + # instrList <- makeInstrList(m, lapply(c(5,3,1,2,4), \(i) instrList$instr_nList[[i]])) + # sortIDs <- lapply(instrList$instr_nList, \(x) x$sortID) # expect_true(all(diff(sapply(sortIDs, \(x) min(x,na.rm=TRUE))) >= 0)) # expect_true(all(sapply(sortIDs, \(x) length(x) == 1 || all(diff(x) >= 1, na.rm=TRUE)))) @@ -1002,15 +1159,15 @@ test_that("calculate works correctly for time series/SSM recursion", { m <- mclass$new() instrList <- makeInstrList(m, c("y", "mu")) - sortIDs <- lapply(instrList, \(x) x$sortID) + sortIDs <- lapply(instrList$instr_nList, \(x) x$sortID) expect_true(all(diff(sapply(sortIDs, \(x) min(x, na.rm = TRUE))) >= 0), "incorrect order") expect_true( all(sapply(sortIDs, \(x) length(x) == 1 || all(diff(x) >= 1, na.rm = TRUE))), "something other than sequential dependence on the past" ) - instrList <- makeInstrList(m, instrList[c(5, 3, 1, 2, 4, 6, 9, 8, 7)]) - sortIDs <- lapply(instrList, \(x) x$sortID) + instrList <- makeInstrList(m, instrList$instr_nList[c(5, 3, 1, 2, 4, 6, 9, 8, 7)]) + sortIDs <- lapply(instrList$instr_nList, \(x) x$sortID) expect_true(all(diff(sapply(sortIDs, \(x) min(x, na.rm = TRUE))) >= 0)) expect_true(all(sapply(sortIDs, \(x) length(x) == 1 || all(diff(x) >= 1, na.rm = TRUE)))) @@ -1091,7 +1248,7 @@ test_that("calculate works correctly for time series/SSM recursion", { m <- mclass$new() instrList <- makeInstrList(m, c("y", "mu")) - sortIDs <- lapply(instrList, \(x) x$sortID) + sortIDs <- lapply(instrList$instr_nList, \(x) x$sortID) expect_true(all(diff(sapply(sortIDs, \(x) min(x, na.rm = TRUE))) >= 0), "incorrect order") expect_true( all(sapply(sortIDs, \(x) length(x) == 1 || all(diff(x) >= 1, na.rm = TRUE))), @@ -1099,8 +1256,8 @@ test_that("calculate works correctly for time series/SSM recursion", { ) # Check proper handling when providing list of instructions. - instrList <- makeInstrList(m, lapply(c(5, 3, 1, 2, 4), \(i) instrList[[i]])) - sortIDs <- lapply(instrList, \(x) x$sortID) + instrList <- makeInstrList(m, lapply(c(5, 3, 1, 2, 4), \(i) instrList$instr_nList[[i]])) + sortIDs <- lapply(instrList$instr_nList, \(x) x$sortID) expect_true(all(diff(sapply(sortIDs, \(x) min(x, na.rm = TRUE))) >= 0)) expect_true(all(sapply(sortIDs, \(x) length(x) == 1 || all(diff(x) >= 1, na.rm = TRUE)))) @@ -1204,25 +1361,25 @@ test_that("basic creation of list of instr_nClass objects", { m <- nimbleModel(code, data = data, returnClass = FALSE) - instr0 <- makeInstrList(m, "mu")[[1]] + instr0 <- makeInstrList(m, "mu")$instr_nList[[1]] expect_identical(instr0$lens, 1) expect_identical(length(instr0$values), 0L) expect_identical(instr0$index_types, 0) expect_identical(instr0$instr_type, 0) - instr1 <- makeInstrList(m, "y[3:4]")[[1]] + instr1 <- makeInstrList(m, "y[3:4]")$instr_nList[[1]] expect_identical(instr1$lens, 2) expect_identical(instr1$values[[1]], 3) # offset expect_identical(instr1$index_types, 1) expect_identical(instr1$instr_type, 1) - instr2 <- makeInstrList(m, c("y[c(2,5)]"))[[1]] + instr2 <- makeInstrList(m, c("y[c(2,5)]"))$instr_nList[[1]] expect_identical(instr2$lens, 2) expect_identical(instr2$values[[1]], c(2, 5)) expect_identical(instr2$index_types, 2) expect_identical(instr2$instr_type, 2) - instr2 <- makeInstrList(m, varRangeClass$new(list(newIndexRange(matrix(c(2, 5), ncol = 1))), varName = "y"))[[1]] + instr2 <- makeInstrList(m, varRangeClass$new(list(newIndexRange(matrix(c(2, 5), ncol = 1))), varName = "y"))$instr_nList[[1]] expect_identical(instr2$lens, 2) expect_identical(instr2$values[[1]], c(2, 5)) expect_identical(instr2$index_types, 2) @@ -1240,22 +1397,22 @@ test_that("basic creation of list of instr_nClass objects", { data <- list(y = matrix(rnorm(20), 5)) m <- nimbleModel(code, data = data) - instr <- makeInstrList(m, varRangeClass$new(list(newIndexRange(quote(2:5)), newIndexRange(quote(1:3))), varName = "y"))[[1]] + instr <- makeInstrList(m, varRangeClass$new(list(newIndexRange(quote(2:5)), newIndexRange(quote(1:3))), varName = "y"))$instr_nList[[1]] expect_identical(instr$instr_type, 4) - instr <- makeInstrList(m, varRangeClass$new(list(newIndexRange(quote(2:5)), newIndexRange(matrix(c(1, 4), ncol = 1))), varName = "y"))[[1]] + instr <- makeInstrList(m, varRangeClass$new(list(newIndexRange(quote(2:5)), newIndexRange(matrix(c(1, 4), ncol = 1))), varName = "y"))$instr_nList[[1]] expect_identical(instr$instr_type, 5) - instr <- makeInstrList(m, varRangeClass$new(list(newIndexRange(matrix(c(1, 4), ncol = 1)), newIndexRange(quote(2:5))), varName = "y"))[[1]] + instr <- makeInstrList(m, varRangeClass$new(list(newIndexRange(matrix(c(1, 4), ncol = 1)), newIndexRange(quote(2:5))), varName = "y"))$instr_nList[[1]] expect_identical(instr$instr_type, 6) - instr <- makeInstrList(m, varRangeClass$new(list(newIndexRange(matrix(c(1, 4), ncol = 1)), newIndexRange(matrix(c(2, 4), ncol = 1))), varName = "y"))[[1]] + instr <- makeInstrList(m, varRangeClass$new(list(newIndexRange(matrix(c(1, 4), ncol = 1)), newIndexRange(matrix(c(2, 4), ncol = 1))), varName = "y"))$instr_nList[[1]] expect_identical(instr$instr_type, 7) ## order is shuffled to put first index slot in first range instr <- makeInstrList(m, varRangeClass$new(list(newIndexRange(quote(2:5)), newIndexRange(matrix(c(1, 4), ncol = 1))), rangeToIndexSlot = list(2, 1), varName = "y" - ))[[1]] + ))$instr_nList[[1]] expect_identical(instr$instr_type, 6) expect_identical(instr$slots, c(1, 2)) expect_identical(instr$index_types, c(2, 1)) @@ -1263,7 +1420,7 @@ test_that("basic creation of list of instr_nClass objects", { ## order is shuffled to put first index slot in first range instr <- makeInstrList(m, varRangeClass$new(list(newIndexRange(matrix(c(1, 4), ncol = 1)), newIndexRange(quote(2:5))), rangeToIndexSlot = list(2, 1), varName = "y" - ))[[1]] + ))$instr_nList[[1]] expect_identical(instr$instr_type, 5) expect_identical(instr$slots, c(1, 2)) expect_identical(instr$index_types, c(1, 2)) @@ -1281,22 +1438,22 @@ test_that("basic creation of list of instr_nClass objects", { data <- list(y = matrix(rnorm(20), 4)) m <- nimbleModel(code, data = data) - instr <- makeInstrList(m, varRangeClass$new(list(newIndexRange(quote(1:3)), newIndexRange(quote(2:5))), varName = "y"))[[1]] + instr <- makeInstrList(m, varRangeClass$new(list(newIndexRange(quote(1:3)), newIndexRange(quote(2:5))), varName = "y"))$instr_nList[[1]] expect_identical(instr$instr_type, 4) - instr <- makeInstrList(m, varRangeClass$new(list(newIndexRange(quote(1:3)), newIndexRange(matrix(c(2, 5), ncol = 1))), varName = "y"))[[1]] + instr <- makeInstrList(m, varRangeClass$new(list(newIndexRange(quote(1:3)), newIndexRange(matrix(c(2, 5), ncol = 1))), varName = "y"))$instr_nList[[1]] expect_identical(instr$instr_type, 6) # shuffled - instr <- makeInstrList(m, varRangeClass$new(list(newIndexRange(matrix(c(1, 4), ncol = 1)), newIndexRange(quote(2:5))), varName = "y"))[[1]] + instr <- makeInstrList(m, varRangeClass$new(list(newIndexRange(matrix(c(1, 4), ncol = 1)), newIndexRange(quote(2:5))), varName = "y"))$instr_nList[[1]] expect_identical(instr$instr_type, 5) # shuffled - instr <- makeInstrList(m, varRangeClass$new(list(newIndexRange(matrix(c(1, 4), ncol = 1)), newIndexRange(matrix(c(2, 5), ncol = 1))), varName = "y"))[[1]] + instr <- makeInstrList(m, varRangeClass$new(list(newIndexRange(matrix(c(1, 4), ncol = 1)), newIndexRange(matrix(c(2, 5), ncol = 1))), varName = "y"))$instr_nList[[1]] expect_identical(instr$instr_type, 7) ## order is shuffled to put first index slot in first range instr <- makeInstrList(m, varRangeClass$new(list(newIndexRange(quote(2:5)), newIndexRange(matrix(c(1, 4), ncol = 1))), rangeToIndexSlot = list(2, 1), varName = "y" - ))[[1]] + ))$instr_nList[[1]] expect_identical(instr$instr_type, 5) expect_identical(instr$slots, c(1, 2)) expect_identical(instr$index_types, c(1, 2)) @@ -1304,7 +1461,7 @@ test_that("basic creation of list of instr_nClass objects", { ## order is shuffled to put first index slot in first range instr <- makeInstrList(m, varRangeClass$new(list(newIndexRange(matrix(c(1, 5), ncol = 1)), newIndexRange(quote(2:4))), rangeToIndexSlot = list(2, 1), varName = "y" - ))[[1]] + ))$instr_nList[[1]] expect_identical(instr$instr_type, 6) expect_identical(instr$slots, c(1, 2)) expect_identical(instr$index_types, c(2, 1)) @@ -1325,7 +1482,7 @@ test_that("basic creation of list of instr_nClass objects", { m <- nimbleModel(code, data = data) instr <- makeInstrList(m, varRangeClass$new(list(newIndexRange(quote(2:5)), newIndexRange(matrix(c(1, 2, 2, 3), ncol = 2, byrow = TRUE))), varName = "y" - ))[[1]] + ))$instr_nList[[1]] expect_identical(instr$slots, c(1, 2, 3)) expect_identical(instr$index_types, c(1, 2)) expect_identical(instr$values[[2]], c(1, 2, 2, 3)) @@ -1339,14 +1496,14 @@ test_that("basic creation of list of instr_nClass objects", { newIndexRange(quote(1:2)) ), varName = "y" - ))[[1]] + ))$instr_nList[[1]] expect_identical(instr$instr_type, 12) # allseq ## indexRange order is shuffled to put first index slot in first indexRange. instr <- makeInstrList(m, varRangeClass$new(list(newIndexRange(quote(2:5)), newIndexRange(matrix(c(1, 2, 2, 3), ncol = 2, byrow = TRUE))), varName = "y", rangeToIndexSlot = list(2, c(1, 3)) - ))[[1]] + ))$instr_nList[[1]] expect_identical(instr$slots, c(1, 3, 2)) # shuffled expect_identical(instr$index_types, c(2, 1)) expect_identical(instr$values[[1]], c(1, 2, 2, 3)) @@ -1355,7 +1512,7 @@ test_that("basic creation of list of instr_nClass objects", { ## order within the matrix indexRange is shuffled to be ascending. instr <- makeInstrList(m, varRangeClass$new(list(newIndexRange(quote(2:5)), newIndexRange(matrix(c(1, 2, 2, 3), ncol = 2, byrow = TRUE))), varName = "y", rangeToIndexSlot = list(2, c(3, 1)) - ))[[1]] + ))$instr_nList[[1]] expect_identical(instr$slots, c(1, 3, 2)) # shuffled expect_identical(instr$index_types, c(2, 1)) expect_identical(instr$values[[1]], c(2, 1, 3, 2)) # shuffled @@ -1363,7 +1520,7 @@ test_that("basic creation of list of instr_nClass objects", { instr <- makeInstrList(m, varRangeClass$new(list(newIndexRange(matrix(c(1, 3), ncol = 1)), newIndexRange(matrix(c(1, 2, 2, 3), ncol = 2, byrow = TRUE))), varName = "y" - ))[[1]] + ))$instr_nList[[1]] expect_identical(instr$instr_type, 10) ## Check technique of building and copying nList(instr_nClass) as a method. @@ -1381,14 +1538,46 @@ test_that("basic creation of list of instr_nClass objects", { instrList <- makeInstrList(cm, "y[2:5]") cinstrList <- cm$makeCompiledInstrList(instrList) expect_true(cinstrList$isCompiled()) - expect_true(inherits(cinstrList, "nList")) - expect_identical(length(cinstrList), 1L) - cinstr <- cinstrList[[1]] + expect_true(inherits(cinstrList, "instrList_nClass")) + expect_identical(cinstrList$offsets, 0L) + expect_identical(cinstrList$nNodes_total, 4L) + expect_true(inherits(cinstrList$instr_nList, "nList")) + expect_identical(length(cinstrList$instr_nList), 1L) + cinstr <- cinstrList$instr_nList[[1]] expect_true(cinstr$isCompiled()) expect_identical(cinstr$nDim, 1L) expect_identical(cinstr$dims, 1L) expect_identical(cinstr$lens, 4L) expect_identical(cinstr$values[[1]], 2L) + expect_identical(cinstr$nNodes, 4L) + expect_identical(cinstr$strides, 1L) + expect_identical(cinstr$idxSlotStarts, 0L) + # The compiled list can be used directly. + expect_equal(cm$calculate(cinstrList), cm$calculate("y[2:5]")) +}) + +test_that("instruction lists record each instruction's node count and offset", { + code <- quote({ + mu ~ dnorm(0, 1) + for (i in 1:5) { + for (j in 1:3) { + y[i, j] ~ dnorm(mu, 1) + } + } + }) + m <- nimbleModel(code, returnClass = FALSE) + instrList <- makeInstrList(m, c("mu", "y[2:4, 1:3]")) + expect_true(inherits(instrList, "RinstrList")) + instrs <- instrList$instr_nList + expect_identical(length(instrs), 2L) + expect_identical(instrs[[1]]$nNodes, 1) # mu + expect_identical(instrs[[2]]$nNodes, 9) # y[2:4, 1:3] + expect_identical(instrs[[2]]$strides, c(3, 1)) # last range fastest + expect_identical(instrs[[2]]$idxSlotStarts, c(0, 1)) + expect_identical(instrList$offsets, c(0, 1)) + expect_identical(instrList$nNodes_total, 10) + # Idempotent. + expect_identical(makeInstrList(m, instrList), instrList) }) # Note that whether `k` is considered as a RHSonly node might change; see Issue 28. @@ -2052,3 +2241,35 @@ test_that("duplication cases", { expect_identical(expected, m$y[2:8]) expect_identical(expected, cm$y[2:8]) }) + +test_that("lifted node that uses only some of its loop indices", { + # Regression test: `sqrt(s2[j])` is lifted from `y[i, j]` and keeps the for-loop context `i, j`, + # but its instructions fill only one idx slot (for `j`). Its generated code must use `idx[1]` for `j` + # (it used `idx[2]`, leaving the lifted node NA), and the compiled declaration must be built for + # one loop index (otherwise compiled calculation errors with "instruction type ... cannot occur"). + code <- nimbleCode({ + for (j in 1:3) { + s2[j] ~ dunif(0, 10) + } + for (i in 1:4) { + for (j in 1:3) { + y[i, j] ~ dnorm(0, var = s2[j]) + } + } + }) + s2 <- c(0.25, 1, 4) + y <- matrix(rnorm(12), 4, 3) + mclass <- nimbleModel(code, inits = list(s2 = s2), data = list(y = y), returnClass = TRUE) + cmclass <- nCompile(mclass) + m <- mclass$new() + cm <- cmclass$new() + expected <- dnorm(y, 0, rep(sqrt(s2), each = 4), log = TRUE) + for (obj in list(m, cm)) { + obj$calculate() + expect_equal(obj$calculate("y"), sum(expected)) + for (k in 1:12) { + node <- paste0("y[", (k - 1) %% 4 + 1, ", ", (k - 1) %/% 4 + 1, "]") + expect_equal(obj$getLogProb(node), expected[k]) + } + } +}) diff --git a/nimbleModel/tests/testthat/test-valuesNode.R b/nimbleModel/tests/testthat/test-valuesNode.R new file mode 100644 index 0000000..7d89046 --- /dev/null +++ b/nimbleModel/tests/testthat/test-valuesNode.R @@ -0,0 +1,237 @@ +library(nimbleModel) +library(testthat) +library(nCompiler) + +# values(model, nodes[i], mode = "node"): the values of node i, in sorted order +# (as for calculate_indexed), flattened. This uses getValues_one/setValues_one in each declFun. + +values_node_code <- nimbleCode({ + mu ~ dnorm(0, 1) # scalar, no loop index + tau ~ dunif(0, 10) + for (i in 1:3) { + y[i] ~ dnorm(mu + x[i], tau = tau) # x is RHS-only; tau gives a lifted node + z[i] <- 2 * y[i] + } + for (i in 1:2) { + v[i, 1:3] <- m0[1:3] * i # one range and a fixed index (D = 1) + W[i, 1:2, 1:2] <- B[1:2, 1:2] * i # two ranges and a fixed index (D = 2) + } +}) + +make_values_node_inputs <- function() { + list( + data = list(x = c(0.1, 0.2, 0.3), m0 = c(1, 2, 3), B = matrix(1:4, 2)), + inits = list(mu = 0.5, tau = 2, y = c(1, 2, 3)) + ) +} + +check_values_node <- function(m) { + m$calculate() + liftedNames <- grep("^lifted", m$getVarNames(), value = TRUE) + expect_true(length(liftedNames) > 0) + + # get + expect_equal(values_indexed(m, "mu", 1, "node"), m$mu) + expect_equal(values(m, c("y[2:3]")[1], mode = "node"), m$y[2]) + nodes <- "y[2:3]" + expect_equal(values(m, nodes[2], mode = "node"), m$y[3]) + expect_equal(values_indexed(m, "v", 2, "node"), m$v[2, ]) + expect_equal(values_indexed(m, "W", 1, "node"), as.numeric(m$W[1, , ])) + expect_equal(values_indexed(m, "W", 2, "node"), as.numeric(m$W[2, , ])) + for (ln in liftedNames) { + expect_equal(values_indexed(m, ln, 1, "node"), m[[ln]][1], info = ln) + } + # a partial overlap gives the whole node + expect_equal(values_indexed(m, "v[1, 2]", 1, "node"), m$v[1, ]) + # sorted order: mu comes before z + expect_equal(values_indexed(m, c("z[1:2]", "mu"), 1, "node"), m$mu) + expect_equal(values_indexed(m, c("z[1:2]", "mu"), 3, "node"), m$z[2]) + # varRange and instrList inputs + expect_equal(values_indexed(m, varRangeClass$new("v[1:2, 1:3]"), 1, "node"), m$v[1, ]) + il <- makeInstrList(m, c("z[1:2]", "mu")) + expect_equal(values_indexed(m, il, 2, "node"), m$z[1]) + + # set + values(m, nodes[1], mode = "node") <- 4.5 + expect_equal(m$y[2], 4.5) + values_indexed(m, "mu", 1, "node") <- -1 + expect_equal(m$mu, -1) + values_indexed(m, "v", 2, "node") <- c(7, 8, 9) + expect_equal(m$v[2, ], c(7, 8, 9)) + expect_equal(m$v[1, ], c(1, 2, 3)) + newW <- c(11, 12, 13, 14) + values_indexed(m, "W", 2, "node") <- newW + expect_equal(as.numeric(m$W[2, , ]), newW) + expect_equal(as.numeric(m$W[1, , ]), as.numeric(matrix(1:4, 2))) + + # errors + expect_error(values_indexed(m, "x[1:2]", 1, "node"), "RHS-only") + expect_error(values_indexed(m, "logProb_y[1:2]", 1, "node"), "logProb") + expect_error(values_indexed(m, "W", 3, "node")) + expect_error(values_indexed(m, "W", 1, "node") <- 1:3, "wrong length") + expect_error(values_indexed(m, "mu", 1, "node") <- c(1, 2), "wrong length") + expect_error(values_indexed(m, "v", 1, "node") <- c(1, 2), "wrong length") +} + +test_that("values(model, nodes[i], mode = \"node\") works, uncompiled", { + inputs <- make_values_node_inputs() + m <- nimbleModel(values_node_code, data = inputs$data, inits = inputs$inits) + check_values_node(m) +}) + +test_that("values(model, nodes[i], mode = \"node\") works, compiled", { + inputs <- make_values_node_inputs() + cm <- nimbleModel(values_node_code, data = inputs$data, inits = inputs$inits, compile = TRUE) + check_values_node(cm) +}) + +test_that("compiled code can call model$getValues_indexed and model$setValues_indexed", { + inputs <- make_values_node_inputs() + mclass <- nimbleModel(values_node_code, data = inputs$data, inits = inputs$inits, returnClass = TRUE) + valuesCaller <- nClass( + classname = "valuesCaller", + Cpublic = list( + model = "nimbleModel:::modelBase_nClass()", + nodes = "nimbleModel:::instrList_nClass()", + get = nFunction( + function(i = "integerScalar") { + return(model$getValues_indexed(nodes, i)) + }, + returnType = "numericVector" + ), + set = nFunction( + function(i = "integerScalar", v = "numericVector") { + model$setValues_indexed(nodes, i, v) + }, + returnType = "void" + ), + getAll = nFunction( + function() { + return(model$getValues_nodes(nodes)) + }, + returnType = "numericVector" + ), + setAll = nFunction( + function(v = "numericVector") { + model$setValues_nodes(nodes, v) + }, + returnType = "void" + ) + ) + ) + comp <- nCompile(mclass, valuesCaller) + cm <- comp$mclass$new() + caller <- comp$valuesCaller$new() + caller$model <- cm + instrList <- makeInstrList(cm, c("W", "mu")) + caller$nodes <- cm$makeCompiledInstrList(instrList) + cm$calculate() + # sorted order: mu, W[1, , ], W[2, , ] + expect_equal(caller$get(1L), cm$mu) + expect_equal(caller$get(2L), as.numeric(cm$W[1, , ])) + caller$set(3L, c(4, 3, 2, 1)) + expect_equal(as.numeric(cm$W[2, , ]), c(4, 3, 2, 1)) + caller$set(1L, 2.5) + expect_equal(cm$mu, 2.5) + expect_error(caller$set(1L, c(1, 2)), "wrong length") + expect_error(caller$set(2L, c(1, 2)), "wrong length") + expect_error(caller$get(4L), "outside") + # all the nodes: mu, W[1, , ], W[2, , ] + expect_equal(caller$getAll(), c(cm$mu, as.numeric(cm$W[1, , ]), as.numeric(cm$W[2, , ]))) + caller$setAll(1:9) + expect_equal(cm$mu, 1) + expect_equal(as.numeric(cm$W[1, , ]), 2:5) + expect_equal(as.numeric(cm$W[2, , ]), 6:9) + expect_error(caller$setAll(1:8), "length") +}) + +# Alignment: node i of values(model, nodes[i], mode = "node") is the node that calculate(nodes[i]) +# and getLogProb(nodes[i]) act on. The reference is the sorted, expanded node names. +check_values_node_alignment <- function(m, nodes) { + nodeNames <- m$topologicallySortNodes(nodes) + instrList <- makeInstrList(m, nodes, includeData = TRUE) + expect_equal(instrList$nNodes_total, length(nodeNames)) + for (k in seq_along(nodeNames)) { + info <- paste("node", k, nodeNames[k]) + expect_equal(values_indexed(m, nodes, k, "node"), values(m, nodeNames[k]), info = info) + expect_equal(m$getLogProb_indexed(nodes, k), m$getLogProb(nodeNames[k]), info = info) + # Setting node k changes exactly that node. + old <- lapply(nodeNames, \(nn) values(m, nn)) + new_k <- old[[k]] + 100 + values_indexed(m, nodes, k, "node") <- new_k + expected <- old + expected[[k]] <- new_k + expect_equal(lapply(nodeNames, \(nn) values(m, nn)), expected, info = info) + # calculate_indexed acts on the same node: its logProb now matches a calculate of that node. + lp <- m$calculate_indexed(nodes, k) + expect_equal(lp, m$getLogProb(nodeNames[k]), info = info) + values_indexed(m, nodes, k, "node") <- old[[k]] + m$calculate(nodeNames[k]) + } + # Unindexed node mode: all the nodes, concatenated in the same order. + allOld <- lapply(nodeNames, \(nn) values(m, nn)) + expect_equal(values(m, nodes, mode = "node"), unlist(allOld, use.names = FALSE)) + newAll <- unlist(allOld, use.names = FALSE) + 1000 + values(m, nodes, mode = "node") <- newAll + expect_equal(unlist(lapply(nodeNames, \(nn) values(m, nn)), use.names = FALSE), newAll) + # A wrong length is an error before anything is set. + expect_error(values(m, nodes, mode = "node") <- newAll[-1], "length") + expect_error(values(m, nodes, mode = "node") <- c(newAll, 1), "length") + expect_equal(values(m, nodes, mode = "node"), newAll) + values(m, nodes, mode = "node") <- unlist(allOld, use.names = FALSE) + expect_equal(lapply(nodeNames, \(nn) values(m, nn)), allOld) + m$calculate(nodeNames) +} + +check_values_edge_cases <- function(m) { + y <- m$y + # duplicates and overlaps: raw and scalar modes keep them, node mode has each node once + expect_equal(values(m, c("mu", "mu")), c(m$mu, m$mu)) + expect_equal(values_indexed(m, c("y[1:2]", "y[2:3]"), 2, "raw"), y[2:3]) + expect_equal(values_indexed(m, c("y[1:2]", "y[2:3]"), 3, "scalar"), y[2]) + dupNodes <- c("mu", "mu", "y[1:2]", "y[2:3]") + expect_equal(makeInstrList(m, dupNodes, includeData = TRUE)$nNodes_total, 4) + expect_equal(values_indexed(m, dupNodes, 3, "node"), y[2]) + expect_error(values_indexed(m, dupNodes, 5, "node")) + # raw assignment with an overlap: later entries win + values(m, c("y[1:2]", "y[2:3]")) <- c(11, 12, 13, 14) + expect_equal(m$y, c(11, 13, 14)) + m$y <- y + + # i given as an expression, in each mode + nodes <- c("y[1:2]", "mu", "y[3]") + j <- 1 + expect_equal(values(m, nodes[j + 1]), m$mu) + expect_equal(values(m, nodes[j + 1], mode = "scalar"), y[2]) + expect_equal(values(m, nodes[2 * j + 1], mode = "node"), y[2]) # sorted: mu, y[1], y[2], y[3] + values(m, nodes[j + 2], mode = "raw") <- 30 + expect_equal(m$y[3], 30) + values(m, nodes[j + 1], mode = "scalar") <- 20 + expect_equal(m$y[2], 20) + values(m, nodes[j], mode = "node") <- -2 + expect_equal(m$mu, -2) + m$y <- y +} + +test_that("values(nodes[i], mode = \"node\") aligns with calculate and getLogProb, and edge cases work", { + data <- list(x = c(0.1, 0.2, 0.3), m0 = c(1, 2, 3), B = matrix(1:4, 2), y = c(1, 2, 3)) # y is data + inits <- list(mu = 0.5, tau = 2) + for (compile in c(FALSE, TRUE)) { + m <- nimbleModel(values_node_code, data = data, inits = inits, compile = compile) + m$calculate() + check_values_node_alignment(m, m$getNodeNames()) # every node, including lifted and data nodes + check_values_node_alignment(m, c("z", "W[2, 1:2, 1:2]", "y[2:3]", "mu")) # input order is not sorted + check_values_edge_cases(m) + } +}) + +test_that("nimbleModelOptions(check_values_node_inputs = FALSE) skips check_values_node_input", { + inputs <- make_values_node_inputs() + m <- nimbleModel(values_node_code, data = inputs$data, inits = inputs$inits) + expect_error(nimbleModel:::check_values_node_input(m, "x[1:2]"), "RHS-only") + expect_true(nimbleModel:::check_values_node_input(m, "mu")) + setNimbleModelOption("check_values_node_inputs", FALSE) + on.exit(setNimbleModelOption("check_values_node_inputs", TRUE)) + expect_true(nimbleModel:::check_values_node_input(m, "x[1:2]")) + expect_true(nimbleModel:::check_values_node_input(m, "logProb_y[1:2]")) +})