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 f569193..b729403 100644 --- a/nimbleModel/R/declFunBaseClass.R +++ b/nimbleModel/R/declFunBaseClass.R @@ -39,6 +39,41 @@ declFunBase_nClass <- nClass( 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 + }, # getParam: use only the first entry in the instr. getParam = function(instr, paramID) { fi <- first_idx(instr) @@ -348,6 +383,78 @@ declFunBase_nClass <- nClass( 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 diff --git a/nimbleModel/R/modelBaseClass.R b/nimbleModel/R/modelBaseClass.R index 2b888e7..77b9d29 100644 --- a/nimbleModel/R/modelBaseClass.R +++ b/nimbleModel/R/modelBaseClass.R @@ -16,7 +16,11 @@ update_model_calc_methods_LAT <- function(code, symTab, auxEnv, info) { getBound_indexed = "getBound_indexed_impl", getParam = "getParam_impl", getParam_indexed = "getParam_indexed_impl", - getParamID = "getParamID_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", @@ -31,7 +35,11 @@ update_model_calc_methods_LAT <- function(code, symTab, auxEnv, info) { getBound_indexed = "numericScalar", getParam = "ETaccessor()", getParam_indexed = "ETaccessor()", - getParamID = "integerScalar" + getParamID = "integerScalar", + getValues_indexed = "numericVector", + setValues_indexed = "void", + getValues_nodes = "numericVector", + setValues_nodes = "void" ) code$name <- method code$type <- nCompiler:::type2symbol(!!returnType) @@ -709,6 +717,60 @@ modelBase_nClass <- nClass( 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$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.") + } + offset <- 0 + for (instr in instrs) offset <- declFunList[[instr$declID]]$setValues_instr(instr, value, offset) + invisible(NULL) } ), Cpublic = list( @@ -938,6 +1000,62 @@ modelBase_nClass <- nClass( }, 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 @@ -949,7 +1067,8 @@ modelBase_nClass <- nClass( c( "calculate", "calculate_indexed", "calculateDiff", "calculateDiff_indexed", "simulate", "simulate_indexed", "getLogProb", "getLogProb_indexed", - "getBound", "getBound_indexed", "getParam", "getParam_indexed", "getParamID" + "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")) { @@ -966,6 +1085,14 @@ modelBase_nClass <- nClass( 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. diff --git a/nimbleModel/R/nimbleModel.R b/nimbleModel/R/nimbleModel.R index 4ce6991..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] @@ -672,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) { 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/include/nimbleModel/predef/declFunBase_nC/declFunBase_nClass_cppContent.cpp b/nimbleModel/inst/include/nimbleModel/predef/declFunBase_nC/declFunBase_nClass_cppContent.cpp index 86558fd..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 @@ -80,6 +80,30 @@ return(ETaccessPtr(0.0)); 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 @@ -112,7 +136,12 @@ method("calculate_singleton_cpp", &declFunBase_nClass::calculate_singleton_cpp, 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("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 ccaa8a6..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 @@ -26,6 +26,11 @@ class declFunBase_nClass : public interface_resolver< genericInterfaceC 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 1d10da4..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(1793304473.07212, 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 f99e61d..eca1920 100644 --- a/nimbleModel/inst/include/nimbleModel/predef/declFunClass_/declFunClass_.h +++ b/nimbleModel/inst/include/nimbleModel/predef/declFunClass_/declFunClass_.h @@ -2,6 +2,30 @@ // Add "#include " the compileInfo$HpostCode // so it appears after the declaration of nodeFxnBase_nClass. +// 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 { @@ -74,6 +98,87 @@ class declFunClass_ : public declFunBase_nClass { } } + // 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); } 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 faffc8a..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 @@ -91,6 +91,24 @@ Rprintf("modelBase_nClass getParam_indexed_impl (should not see this)\n"); retur 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 @@ -127,7 +145,11 @@ method("calculateDiff_indexed_impl", &modelBase_nClass::calculateDiff_indexed_im 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("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 da05c69..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 @@ -30,6 +30,10 @@ class modelBase_nClass : public interface_resolver< genericInterfaceC 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 8167e63..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(1793304537.50563, 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 0eceef9..79b6446 100644 --- a/nimbleModel/inst/include/nimbleModel/predef/modelClass_/modelClass_.h +++ b/nimbleModel/inst/include/nimbleModel/predef/modelClass_/modelClass_.h @@ -103,6 +103,41 @@ class modelClass_ : public modelBase_nClass { 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 @@ -234,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-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]")) +})