diff --git a/nimbleModel/.DS_Store b/nimbleModel/.DS_Store new file mode 100644 index 0000000..59c8ece Binary files /dev/null and b/nimbleModel/.DS_Store differ diff --git a/nimbleModel/DESCRIPTION b/nimbleModel/DESCRIPTION index 0f1462c..770b628 100644 --- a/nimbleModel/DESCRIPTION +++ b/nimbleModel/DESCRIPTION @@ -3,9 +3,9 @@ Title: What the Package Does (one line, title case) Version: 0.0.0.9000 Authors@R: person("First", "Last", email = "first.last@example.com", role = c("aut", "cre")) Description: What the package does (one paragraph). -Depends: R (>= 3.3.3) +Depends: R (>= 4.1.0) Imports: - methods,R6 + methods,R6,nCompiler License: What license is it under? Encoding: UTF-8 LazyData: true @@ -43,6 +43,10 @@ Collate: varRange.R varRules.R varStore.R + declFunBaseClass.R + modelBaseClass.R + instructions.R + nimbleModel.R diff --git a/nimbleModel/NAMESPACE b/nimbleModel/NAMESPACE index 9f9aabf..a7cc763 100644 --- a/nimbleModel/NAMESPACE +++ b/nimbleModel/NAMESPACE @@ -1,5 +1,6 @@ import(methods) import(R6) +import(nCompiler) ## Some exports are for testing convenience and may not need to be exported. ## Classes: R6 classes are exported as regular objects @@ -68,3 +69,6 @@ export(messageIfVerbose) export(calc_dmnormAltParams) export(calc_dwishAltParams) + +export(nimbleModel) +export(makeInstrList) diff --git a/nimbleModel/R/declFunBaseClass.R b/nimbleModel/R/declFunBaseClass.R new file mode 100644 index 0000000..adae7a6 --- /dev/null +++ b/nimbleModel/R/declFunBaseClass.R @@ -0,0 +1,258 @@ +#' @export +declFunBase_nClass <- nClass( + classname = "declFunBase_nClass", + Rpublic = list( + calculate = function(instr) { + calc_op(instr, "calc_one") + }, + calculateDiff = function(instr) { + calc_op(instr, "calcDiff_one") + }, + getLogProb = function(instr) { + calc_op(instr, "getLogProb_one") + }, + calc_op = function(instr, fn) { + if(instr$type == 0) return(calc_0(instr, fn)) + if(instr$type == 1) return(calc_1_seq(instr, fn)) + if(instr$type == 2) return(calc_1_mat(instr, fn)) + if(instr$type == 3) return(calc_1_matp(instr, fn)) + if(instr$type == 4) return(calc_1_matp_ord(instr, fn)) + if(instr$type == 5) return(calc_2_seq_seq(instr, fn)) + if(instr$type == 6) return(calc_2_seq_seq_ord(instr, fn)) + return(0) + }, + calc_0 = function(instr, fn) { + return(self[[fn]](0)) + }, + calc_1_seq = + function(instr, fn) { + logProb = 0 + iStart <- instr$values[[1]][1]-1 + for(i in 1:instr$lens[1]) + logProb <- logProb + self[[fn]](iStart + 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 + for(i in 1:instr$lens[1]) + logProb <- logProb + self[[fn]](instr$values[[1]][(instr$dims[1]*(i-1) + 1):(instr$dims[1]*i)][instr$slots]) + return(logProb) + }, + calc_2_seq_seq = + function(instr, fn) { + logProb <- 0 + idx <- rep(0, 2) + iStart1 <- instr$values[[1]][1]-1 + for(i in 1:instr$lens[1]) { + idx[1] <- iStart1 + i + iStart2 <- instr$values[[2]][1]-1 + for(j in 1:instr$lens[2]) { + idx[2] <- iStart2 + j + logProb <- logProb + self[[fn]](idx) + } + } + return(logProb) + }, + calc_2_seq_seq_ord = + function(instr, fn) { + if(!identical(instr$slots, 1:2)) + stop("Slots not equal to 2,1 in calc_2_seq_seq_ord") + logProb <- 0 + idx <- rep(0, 2) + iStart1 <- instr$values[[1]][1]-1 + for(i in 1:instr$lens[1]) { + idx[2] <- iStart1 + i + iStart2 <- instr$values[[2]][1]-1 + for(j in 1:instr$lens[2]) { + idx[1] <- iStart2 + j + logProb <- logProb + self[[fn]](idx) + } + } + return(logProb) + }, + simulate = function(instr) { + 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_1_matp_ord(instr)) + if(instr$type == 5) return(sim_2_seq_seq(instr)) + if(instr$type == 6) return(sim_2_seq_seq_ord(instr)) + }, + sim_0 = function(instr) { + sim_one(0) + }, + sim_1_seq = function(instr) { + for(i in 1:instr$lens[1]) + sim_one(instr$values[[1]][1]+i) + }, + sim_1_mat = function(instr) { + for(i in 1:instr$lens[1]) + sim_one(instr$values[[1]][i]) + }, + 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)]) + }, + sim_1_matp_ord = function(instr) { + for(i in 1:instr$lens[1]) + sim_one(instr$values[[1]][(instr$dims[1]*(i-1) + 1):(instr$dims[1]*i)][instr$slots]) + }, + sim_2_seq_seq = + function(instr, fn) { + idx <- rep(0, 2) + iStart1 <- instr$values[[1]][1]-1 + for(i in 1:instr$lens[1]) { + idx[1] <- iStart1 + i + iStart2 <- instr$values[[2]][1]-1 + for(j in 1:instr$lens[2]) { + idx[2] <- iStart2 + j + sim_one(idx) + } + } + }, + sim_2_seq_seq_ord = + function(instr, fn) { + if(!identical(instr$slots, 1:2)) + stop("Slots not equal to 2,1 in calc_2_seq_seq_ord") + idx <- rep(0, 2) + iStart1 <- instr$values[[1]][1]-1 + for(i in 1:instr$lens[1]) { + idx[2] <- iStart1 + i + iStart2 <- instr$values[[2]][1]-1 + for(j in 1:instr$lens[2]) { + idx[1] <- iStart2 + j + sim_one(idx) + } + } + } + + + ), + + Cpublic = list( + ## model = 'modelBase_nClass', + ping = nFunction( + name = "ping", + function() {return(TRUE); returnType(logical())}, + compileInfo = list(virtual=TRUE) + ), + calculate_cpp = nFunction( + name = "calculate_cpp", + function(instr) { + stop("Uncompiled version of calculate_cpp should not be called.") + }, + returnType = 'numericScalar', + compileInfo = list(virtual=TRUE, + C_fun = function(instr = 'instr_nClass') { + cppLiteral('Rprintf("declFunBase_nClass virtual base calculate_cpp should never be called (something is wrong)\\n");') + return(0) + }) + ), + calculateDiff_cpp = nFunction( + name = "calculateDiff_cpp", + function(instr) { + stop("Uncompiled version of calculateDiff_cpp should not be called.") + }, + returnType = 'numericScalar', + compileInfo = list(virtual=TRUE, + C_fun = function(instr = 'instr_nClass') { + cppLiteral('Rprintf("declFunBase_nClass virtual base calculateDiff_cpp should never be called (something is wrong)\\n");') + return(0) + }) + ), + getLogProb_cpp = nFunction( + name = "getLogProb_cpp", + function(instr) { + stop("Uncompiled version of getLogProb_cpp should not be called.") + }, + returnType = 'numericScalar', + compileInfo = list(virtual=TRUE, + C_fun = function(instr = 'instr_nClass') { + cppLiteral('Rprintf("declFunBase_nClass virtual base getLogProb_cpp should never be called (something is wrong)\\n");') + return(0) + }) + ), + simulate_cpp = nFunction( + name = "simulate_cpp", + function(instr) { + stop("Uncompiled version of simulate_cpp should not be called.") + }, + returnType = 'void', + compileInfo = list(virtual=TRUE, + C_fun = function(instr = 'instr_nClass') { + cppLiteral('Rprintf("declFunBase_nClass virtual base simulate_cpp should never be called (something is wrong)\\n");') + }) + ) + + + + + + # simulate = nFunction( + # name = "simulate", + # fun = function(instr = 'instr_nClass') { + # ## TODO: how embed determination of vec and parallel cases here? + # if(instr$type == 0) sim_0(instr) + # if(instr$type == 1) sim_1_seq(instr) + # if(instr$type == 2) sim_1_mat(instr) + # if(instr$type == 3) sim_1_matp(instr) + # }, + # compileInfo = list(virtual=TRUE) + # ), + # sim_0 = nFunction( + # name = 'sim_0', + # function(instr = 'instr_nClass') { + # sim_one(0) ## sim_one will always has `idx` as arg? + # } + # ), + # sim_1_seq = nFunction( + # name = 'sim_1_seq', + # function(instr = 'instr_nClass') { + # for(i in 1:instr$lens[1]) + # sim_one(instr$values[[1]][1]+i) + # } + # ), + # sim_1_mat = nFunction( + # name = 'sim_1_mat', + # function(instr = 'instr_nClass') { + # for(i in 1:instr$lens[1]) + # sim_one(instr$values[[1]][i]) + # } + # ), + # sim_1_matp = nFunction( + # name = 'sim_1_mat', + # function(instr = 'instr_nClass') { + # for(i in 1:instr$lens[1]) + # sim_one(instr$values[[1]][(instr$dims[1]*(i-1) + 1):(instr$dims[1]*i)]) ## Ok to call with a vector? + # } + # ), + + + ), + ## We haven't dealt with ensuring a virtual destructor when any method is virtual + ## For now I did it manually by editing the .h and .cpp + predefined = quote(system.file(file.path("include","nimbleModel", "predef"), package="nimbleModel") |> + file.path("declFunBase_nC")), + compileInfo=list(interface="full", + createFromR = FALSE, + exportName = "declFunBase_nClass_new", + needed_units = list("instr_nClass"), + packageNames = c(uncompiled="declFunBase_nClass_R", compiled="declFunBase_nClass") + ) +) diff --git a/nimbleModel/R/instructions.R b/nimbleModel/R/instructions.R new file mode 100644 index 0000000..f614280 --- /dev/null +++ b/nimbleModel/R/instructions.R @@ -0,0 +1,194 @@ +type2itype <- list( + '0' = 0, + "1_seq" = 1, + "1_mat" = 2, + "1_matp" = 3, + "2_seq_seq" = 4, # done + "2_seq_mat" = 5, + "2_mat_seq" = 6, + "2_mat_mat" = 7, + "2_seq_matp" = 8, + "2_matp_seq" = 9, + "2_mat_matp" = 10, + "2_matp_mat" = 11, + "2_matp_matp" = 12, + "3_seq_seq_seq" = 13, + "3_generic" = 14, + "1_matp_ord" = 15, # done + ### Probably not needed given reordering when apply graph rules to create instrList. ### + "2_seq_seq_ord" = 16, # done + "2_seq_mat_ord" = 17, + "2_mat_seq_ord" = 18, + "2_mat_mat_ord" = 19, + ######################################################################################### + "2_seq_matp_ord" = 20, + "2_matp_seq_ord" = 21, + "2_mat_matp_ord" = 22, + "2_matp_mat_ord" = 23, + "2_matp_matp_ord" = 24, + "3_generic_ord" = 25 +) + +## Stand-alone function for setting up inputs to instrClass constructor. +## This could be embedded in the constructor. +range2instr <- function(range) { + instr <- list() + if(!length(range$indexingRange$indexRanges)) { # No indexing + instr$lens <- 1 + instr$index_types <- 0 + instr$nDim <- 0 + } else { + instr$lens <- sapply(range$indexingRange$indexRanges, function(x) x$numElements) + instr$dims <- sapply(range$indexingRange$rangeToIndexSlot, length) + instr$nDim <- sum(instr$dims) + instr$slots <- unlist(range$indexingRange$rangeToIndexSlot) + instr$index_types <- sapply(range$indexingRange$indexRanges, function(x) + switch(class(x)[1], + "indexRangeScalarClass" = 2, + "indexRangeSequenceClass" = 1, + "indexRangeMatrixClass" = 2)) + instr$values <- lapply(range$indexingRange$indexRanges, function(x) + switch(class(x)[1], + "indexRangeScalarClass" = x$value, + "indexRangeSequenceClass" = x$start, + "indexRangeMatrixClass" = c(t(matrix(x$values, nc = x$numColumns))))) # in calcRange, column major; need row major here for simpler/more efficient determination of indices + } + instr$type <- determineInstrType(instr) + instr$sortID <- range$sortID + instr$declID <- range$declID + 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. +## Open question of when to determine if to use parallel calculate. +determineInstrType <- function(instr, use_vec = FALSE) { + type <- NULL + if(!length(instr$dims)) + type <- "0" + if(length(instr$dims) == 1) + if(instr$index_types[1] == 1) { + type <- "1_seq" + } else { + if(instr$dims[1] == 1) type <- "1_mat" else { + if(identical(instr$slots, 1:length(instr$slots))) type <- "1_matp" else type <- "1_matp_ord" + } + } + if(length(instr$dims) == 2) + if(identical(instr$dims, c(1L,1L))) { + ## Some of these not yet written. + if(identical(instr$index_types, c(1,1))) + if(identical(instr$slots, 1:2)) type <- "2_seq_seq" else type <- "2_seq_seq_ord" + if(identical(instr$index_types, c(1,2))) + if(instr$dims[2] == 1) type <- "2_seq_mat" else type <- "2_seq_matp" + if(identical(instr$index_types, c(2,1))) + if(instr$dims[1] == 1) type <- "2_mat_seq" else type <- "2_matp_seq" + if(identical(instr$index_types, c(2,2))) { + if(all(instr$dims == 1)) type <- "2_mat_mat" + if(all(instr$dims == 2)) type <- "2_matp_matp" + if(instr$dims[[1]] == 2) type <- "2_matp_mat" + if(instr$dims[[2]] == 2) type <- "2_mat_matp" + } + } else type <- "2_generic" + if(length(instr$dims) == 3) type <- "3_generic" + if(is.null(type)) stop("no available specific instruction type") + ## TODO: determine how much about slots will be pre-baked. + if(length(instr$dims) && !identical(instr$slots, 1:instr$nDim)) # Non-canonical slot ordering + type <- paste(type, "slots", sep = "_") + return(type2itype[[type]]) +} + +## TODO: document this since it may be user-facing. +#' @export +makeInstrList <- function(model, input, use_vec = FALSE) { + ## `model` simply must contain `modelDef`, so it can be a modelClass or modelBase_nClass object. + ## This works with: + ## (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 instr_nClass objects (not assumed to be in sort order) + if(is(input, 'nList')) + if(!inherits(input[[1]], 'instr_nClass')) { + stop("nList input to `makeInstrList` should contain `instr_nClass` objects") + } else return(input) # Idempotent case. + if(is(input, 'instr_nClass')) + input <- list(input) + if(is.list(input) && all(sapply(input, function(x) inherits(x, 'instr_nClass')))) { + ## Create sort-ordered nList. + instrList <- nList(instr_nClass)$new() + numInstrs <- length(input) + instrList$setLength(numInstrs) + ord <- order(unlist(lapply(input, function(x) x$sortID))) + for(i in 1:numInstrs) + instrList[[i]] <- input[[ord[i]]] + return(instrList) + } + ## At this point we presumably are working with varRange(s). + if(is(input, 'varRangeClass')) input <- list(input) + ## First apply calcRule to get overlap between input and the rule. + ## Then make the calcRange to convert to loop indexing. + ## Note that `calcRule$apply` handles converting char to varRange and handling full variable extent. + ranges <- unlist(lapply(input, function(vr) + lapply(model$modelDef$calcRules[[nimbleModel:::getVarName(vr)]]$rules, function(rule) + rule$makeCalcRange(rule$apply(vr)) + ))) + instrList <- nList(instr_nClass)$new() + numRanges <- length(ranges) + instrList$setLength(numRanges) + ord <- order(unlist(lapply(ranges, function(x) x$sortID))) + for(i in 1:numRanges) + instrList[[i]] <- instr_nClass$new(ranges[[ord[i]]]) + return(instrList) +} + +instr_nClass <- nClass( + classname = "instr_nClass", + Rpublic = list( + initialize = function(calcRange, ...) { + super$initialize(...) + if(!missing(calcRange)) { + instr <- range2instr(calcRange) # This processing could simply be included here in `initialize`. + self$lens <- instr$lens %||% integer() + self$index_types <- instr$index_types %||% integer() + self$nDim <- instr$nDim %||% 0L + self$dims <- instr$dims %||% integer() + self$slots <- instr$slots %||% integer() + self$values <- nList(integerVector)$new() + self$values$setLength(length(self$dims)) + if(self$nDim) + for(i in 1:length(self$dims)) + self$values[[i]] <- instr$values[[i]] + self$type <- instr$type %||% 0L # Use integer for compilation (would char be ok?). + self$sortID <- instr$sortID %||% integer() + self$declID <- instr$declID %||% 0L + } + }), + Cpublic = list( + lens = 'integerVector', + index_types = 'integerVector', + nDim = 'integerScalar', + dims = 'integerVector', + slots = 'integerVector', + values = 'nList(integerVector)', + type = 'integerScalar', + sortID = 'integerVector', + declID = 'integerScalar', + instr_nClass = nFunction( + function() { + values <- nList(integerVector)$new() + }, + compileInfo = list(constructor=TRUE) + ) + ), + predefined = quote(system.file(file.path("include","nimbleModel", "predef"), package="nimbleModel") |> file.path("instr_nClass")), + compileInfo = list(interface = "full", + createFromR = TRUE, + exportName = "instr_nClass_new", + needed_units = list("nList(integerVector)"), + packageNames = c(uncompiled="instr_nClass_R", compiled="instr_nClass") + ) +) + + + diff --git a/nimbleModel/R/model.R b/nimbleModel/R/model.R index 4ad3205..4bdc65d 100644 --- a/nimbleModel/R/model.R +++ b/nimbleModel/R/model.R @@ -2,6 +2,12 @@ ## static information, data information, and methods for querying the model ## structure. +## For the moment this overlaps with the custom model class we create +## for each model using nCompiler. That model class takes some of the +## fields and calls some of the methods set up here. + +## TODO: possibly work all of this into modelBase_nClass. + ## Will need to do some work to extend this to get full current behavior ## where the model is created by `nimbleModel` and the custom model class ## contains fields for the different variables. @@ -42,7 +48,9 @@ modelClass <- R6Class( } } - + origInits <<- inits + origData <<- data + if(length(data) && sum(names(data) == "")) stop("modelClass: 'data' must be a named list") if(any(!sapply(data, function(x) { @@ -65,11 +73,6 @@ modelClass <- R6Class( makeDataRules(data) makePredictiveRules() - - ## Do this once we have a custom model class with fields we can assign into. - ## setData(data) - ## setInits(inits) - }, makeDataRules = function(data) { @@ -124,7 +127,7 @@ modelClass <- R6Class( nonpredictiveRules <<- candidateRules }, - + ## TODO: should this be a standalone function like getNodes, getDependencies? getVarNames = function(includeLogProb = FALSE, nodeRanges) { if(missing(nodeRanges)){ if(includeLogProb) return(modelDef$varNames) diff --git a/nimbleModel/R/modelBaseClass.R b/nimbleModel/R/modelBaseClass.R new file mode 100644 index 0000000..15d3ee9 --- /dev/null +++ b/nimbleModel/R/modelBaseClass.R @@ -0,0 +1,228 @@ +#' @export +modelBase_nClass <- nClass( + classname = "modelBase_nClass", + Rpublic = list( + modelDef = NULL, + dataRules = NULL, + nondataRules = NULL, + predictiveRules = NULL, + nonpredictiveRules = NULL, + initialize = function(sizes = list(), inits = list(), data = list(), ...) { + # It is not very easy to set debug onto the initialize function, so + # here is a magic flag. + if(isTRUE(.GlobalEnv$.debugModelInit)) browser() + super$initialize(...) + + ## TODO: is there a better way to populate declFunNameToIndex in Cpublic? + declFunNameToIndex <- self$declFunNameToIndex_ + + declFunNames <- names(declFunNameToIndex) + if(isCompiled()) { + # self$setup_decl_mgmt_from_names(declFunNames) + # setting up the canonically indexed vector of node functions + # now happens in the C++ constructor. + } else { + self$declFunList <- list() + length(self$declFunList) <- length(declFunNames) + names(self$declFunList) <- declFunNames + for(declFunName in declFunNames) { + self[[declFunName]] <- eval(as.name(self$CpublicDeclFuns[[declFunName]]))$new() + self[[declFunName]]$setModel(self) + self$declFunList[[declFunNameToIndex[[declFunName]]]] <- self[[declFunName]] + } + } + + ## TODO: create a merge_and_set function that handles all three of the following. + allSizes <- self$defaultSizes + if(!missing(sizes)) + for(nm in names(sizes)) + allSizes[[nm]] <- sizes[[nm]] + ## TODO: should we handle 0-dim sizes elsewhere? + allSizes <- allSizes[sapply(allSizes, length) > 0] + if(length(allSizes)) resize_from_list(allSizes[sapply(allSizes, length) > 0]) + + allInits <- self$defaultInits + if(!missing(inits)) + for(nm in names(inits)) + allInits[[nm]] <- inits[[nm]] + if(length(allInits)) set_from_list(allInits) + + if(missing(inits)) { + allInits <- self$defaultInits + } else + if(length(inits)) set_from_list(inits) + + ## TODO: do we want to handle data differently? + ## TODO: need to work through not setting as 'data' if values are NA; + ## check back against how dataRules work in nimbleModel work. + allData <- self$defaultData + if(!missing(inits)) + for(nm in names(inits)) + allData[[nm]] <- inits[[nm]] + if(length(allData)) set_from_list(allData) + }, + getVarNames = function(includeLogProb = FALSE, nodeRanges) { + if(missing(nodeRanges)){ + if(includeLogProb) return(modelDef$varNames) + else return(names(modelDef$varInfo)) + } else { + if(!is.list(nodeRanges)) + nodeRanges <- list(nodeRanges) + return(unique(sapply(nodeRanges, `[[`, 'varName'))) + } + }, + getDependencies = function(nodes, self = TRUE, downstream = FALSE, immediateOnly = FALSE) { + nimbleModel::getDependencies(modelDef, nodes, self, downstream, immediateOnly) + }, + getParents = function(nodes, self = TRUE, upstream = FALSE, immediateOnly = FALSE) { + nimbleModel::getParents(modelDef, nodes, self, upstream, immediateOnly) + }, + ## TODO: not working because `nimbleModel::getNodes` needs the model not just modelDef. + ## Once we integrate modelClass with modelBase_nClass, we should be able to pass `self`. + getNodes = function(nodes = NULL, stochOnly = FALSE, determOnly = FALSE, + includeData = TRUE, dataOnly = FALSE, + includePredictive = TRUE, predictiveOnly = FALSE, + includeRHSonly = FALSE, + topOnly = FALSE, latentOnly = FALSE, endOnly = FALSE) { + nimbleModel::getNodes(modelDef, nodes, stochOnly, determOnly, includeData, dataOnly, + includePredictive, predictiveOnly, includeRHSonly, + topOnly, latentOnly, endOnly) + }, + calc_op = function(instr, fn, fn_cpp) { + if(missing(instr)) + instr <- getVarNames() + instrList <- makeInstrList(self,instr) + if(isCompiled()) { + if(!instrList$isCompiled()) instrList <- makeCompiledInstrList(instrList) + return(self[[fn_cpp]](instrList)) + } + logProb <- 0 + for(i in 1:length(instrList)) { + logProb <- logProb + declFunList[[instrList[[i]]$declID]][[fn]](instrList[[i]]) + } + return(logProb) + }, + calculate = function(instr) { + logProb <- calc_op(instr, "calculate", "calculate_impl") + return(logProb) + }, + calculateDiff = function(instr) { + logProb <- calc_op(instr, "calculateDiff", "calculateDiff_impl") + return(logProb) + }, + getLogProb = function(instr) { + logProb <- calc_op(instr, "getLogProb", "getLogProb_impl") + return(logProb) + }, + simulate = function(instr) { + if(missing(instr)) + instr <- getVarNames() + instrList <- makeInstrList(self,instr) + if(isCompiled()) { + if(!instrList$isCompiled()) instrList <- makeCompiledInstrList(instrList) + self$simulate_impl(instrList) + } else { + for(i in 1:length(instrList)) { + declFunList[[instrList[[i]]$declID]]$simulate(instrList[[i]]) + } + } + return(invisible(NULL)) + } + ), + Cpublic = list( + ## TODO: using 'RcppObject' was resulting in a symbolTBD error - probably nCompiler issue 186. + declFunList = 'numericScalar', # 'RcppObject', # This won't actually be used in C++, but needs to be in Cpublic for accessibility. + declFunNameToIndex = 'RcppList', # Not sure what type this should be for use in C++. + ping = nFunction( + name = "ping", + function() {return(TRUE); returnType(logical())}, + compileInfo = list(virtual=TRUE) + ), + makeCompiledInstrList = nFunction( + name = "makeCompiledInstrList", + function(input = 'SEXP') { + ans <- nList(instr_nClass)$new() + cppLiteral("ans->set_all_values(input);") + return(ans) + }, + returnType = 'nList(instr_nClass)' + ), + calculate_impl = nFunction( + name = "calculate_impl", + function(instrList) { + cat("Uncompiled `calculate_impl` should never be called.\n") + return(0) + }, + returnType = 'numericScalar', + compileInfo = list( + C_fun = function(instrList = 'nList(instr_nClass)') { + ## NOTE: instrList input will be ordered. + cppLiteral('Rprintf("modelBase_nClass calculate_impl (should not see this)\\n");'); return(0) + }, + virtual=TRUE + ) + ), + calculateDiff_impl = nFunction( + name = "calculateDiff_impl", + function(instrList) { + cat("Uncompiled `calculateDiff_impl` should never be called.\n") + return(0) + }, + returnType = 'numericScalar', + compileInfo = list( + C_fun = function(instrList = 'nList(instr_nClass)') { + ## NOTE: instrList input will be ordered. + cppLiteral('Rprintf("modelBase_nClass calculateDiff_impl (should not see this)\\n");'); return(0) + }, + virtual=TRUE + ) + ), + getLogProb_impl = nFunction( + name = "getLogProb_impl", + function(instrList) { + cat("Uncompiled `getLogProb_impl` should never be called.\n") + return(0) + }, + returnType = 'numericScalar', + compileInfo = list( + C_fun = function(instrList = 'nList(instr_nClass)') { + ## NOTE: instrList input will be ordered. + cppLiteral('Rprintf("modelBase_nClass getLogProb_impl (should not see this)\\n");'); return(0) + }, + virtual=TRUE + ) + ), + simulate_impl = nFunction( + name = "simulate_impl", + function(instrList) { + cat("Uncompiled `simulate_impl` should never be called.\n") + return(invisible(NULL)) + }, + returnType = 'void', + compileInfo = list( + C_fun = function(instrList = 'nList(instr_nClass)') { + ## NOTE: instrList input will be ordered. + cppLiteral('Rprintf("modelBase_nClass simulate_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(interface="full", + createFromR = FALSE, + Hincludes = c('"declFunBase_nClass_c_.h"','"instr_nClass_c_.h"'), + needed_units = list("declFunBase_nClass","instr_nClass", "nList(instr_nClass)"), + exportName = "modelBase_nClass_new", + packageNames = c(uncompiled="modelBase_nClass_R", compiled="modelBase_nClass") + ) +) + +# Manually add +# # "#include " to that file, +# after the header content. + + +# nCompile(modelBase_nClass, control=list(generate_predefined=TRUE)) + diff --git a/nimbleModel/R/modelDecl.R b/nimbleModel/R/modelDecl.R index 685279a..3ab444c 100644 --- a/nimbleModel/R/modelDecl.R +++ b/nimbleModel/R/modelDecl.R @@ -8,7 +8,7 @@ modelDeclClass <- R6Class( public = list( context = NULL, # FUTURE: might just use declRule$context sourceLineNumber = NULL, - stoch = FALSE, # Need this here as used before declRule created. + stoch = FALSE, code = NULL, distributionName = NA, valueExpr = NULL, @@ -118,7 +118,7 @@ modelDeclClass <- R6Class( ## Create declRule and symbolic RHS pieces. processDecl = function(nimFunNames, constants = list(), envir) { - declRule <<- declRuleClass$new(self, sourceLineNumber, context, constants) + declRule <<- declRuleClass$new(self, 0, context, constants) makeSymbolicParentNodes(nimFunNames, constants, envir) invisible(NULL) }, @@ -283,10 +283,6 @@ modelDeclClass <- R6Class( invisible(NULL) }, - buildFunctions = function() { - declRule$buildFunctions(calculateCode, genLogProbExpr()) - }, - genLogProbExpr = function() { if(declRule$decl$stoch) { logProbExpr <- code[[2]] diff --git a/nimbleModel/R/modelDef.R b/nimbleModel/R/modelDef.R index ea3d231..48cc6ec 100644 --- a/nimbleModel/R/modelDef.R +++ b/nimbleModel/R/modelDef.R @@ -16,6 +16,8 @@ modelDefClass <- R6Class( contexts = list(), constants = list(), declInfo = list(), + declFunNameToIndex = list(), + declFunIndexToName = NULL, downstreamRules = NULL, upstreamRules = NULL, calcRules = NULL, @@ -53,6 +55,7 @@ modelDefClass <- R6Class( addRemainingDotParams() ## Add additional altParams as needed. replaceAllConstants() ## Simplify expressions introduced in `addRemainingDotParams`. processDecls(userEnv) ## Create declRules and set up symbolicParentNodes (and flags dynamic indexing). + assignDeclIDs() ## Set sequential declID values and declFun mapping. genAltParams() ## Create altParam expressions and create `calculateCode` (without altParams). genBounds() ## Create bound expressions (modifying `calculateCode`). @@ -73,7 +76,6 @@ modelDefClass <- R6Class( warnRHSonlyDynamicIndexing() - buildFunctions() ## Generate calculate and other functions. invisible(NULL) }, @@ -546,6 +548,13 @@ modelDefClass <- R6Class( invisible(NULL) }, + assignDeclIDs = function() { + for(i in seq_along(declInfo)) + declInfo[[i]]$declRule$ID <- as.character(i) + declFunNameToIndex <<- as.list(1:length(declInfo)) + names(declFunNameToIndex) <<- paste("declFun", 1:length(declInfo), sep = "_") + }, + ## Add additional altParams not already addressed in getting canonical params. addRemainingDotParams = function() { for(iDecl in seq_along(declInfo)) { @@ -664,7 +673,8 @@ modelDefClass <- R6Class( } } } - for(iDI in seq_along(declInfo)) { # Do RHS after all LHS so that check for overlap only concerns LHS + for(iDI in seq_along(declInfo)) { # Do RHS after all LHS so that check for overlap only concerns LHS + decl <- declInfo[[iDI]] for(iRHR in seq_along(decl$rhsOriginalRules)) { rhsRule <- decl$rhsOriginalRules[[iRHR]] rhsVar <- rhsRule$varName @@ -867,13 +877,6 @@ modelDefClass <- R6Class( } } invisible(NULL) - }, - - buildFunctions = function() { - for(i in seq_along(declInfo)) { - declInfo[[i]]$buildFunctions() - } - invisible(NULL) } ) ) @@ -893,6 +896,10 @@ modelDefClass <- R6Class( ## Note: data-related flags not handled as that relates to flags on a model ## and not part of modelDef. +## TODO: these should presumably take the model not modelDef as the first arg. +## Once we integrate modelClass with modelBase_nClass, we should be able to +## pass `self` from the getDeps and getParents methods to these functions. + getDependencies <- function(modelDef, nodes, self = TRUE, downstream = FALSE, immediateOnly = FALSE) { diff --git a/nimbleModel/R/nimbleModel.R b/nimbleModel/R/nimbleModel.R new file mode 100644 index 0000000..6ec1889 --- /dev/null +++ b/nimbleModel/R/nimbleModel.R @@ -0,0 +1,496 @@ +## It's unclear what this should return since in nimble one gets a model object +## but with nCompiler, we need a modelClass to compile and then we can create an instance. +## If we create an instance as the output here, one can't then compile that with an algorithm via `nCompile`. +## Need to think more about the workflow for nimble 2.0. +#' @export +nimbleModel <- function(code, + constants = list(), + data = list(), + inits = list(), + dimensions = list(), + compile = FALSE, + returnClass = TRUE, + where = globalenv(), + debug = FALSE, + check = getNimbleOption('checkModel'), + calculate = TRUE, + name = NULL, + buildDerivs = getNimbleOption('buildModelDerivs'), + userEnv = parent.frame()) { + ## TODO: arg list taken from `nimble`. Revisit which options are needed. + ## For the moment this goes through (original) nimbleModel R6 class and then nimbleModel nClass. Clean that up once ideas are in place. + ## Presumably everything would be in Rpublic initialize for modelBaseClass, so this function will just call modelBase_nClass$new(). + nm <- modelClass$new(name = name, code = code, constants = constants, data = data, inits = inits, dimensions = dimensions, userEnv = userEnv) + specificModelClass <- make_modelClass_from_nimbleModel(nm) + if(compile) specificModelClass <- nCompile(specificModelClass) + if(returnClass) return(specificModelClass) # Standard use for when compiling a model(class) and algo(class) together. + model <- specificModelClass$new() # Otherwise return model object for manipulation from R. +} + +make_modelClass_from_nimbleModel <- function(m, compile=FALSE) { + mDef <- m$modelDef + modelVarInfo <- get_varInfo_from_nimbleModel(m) + declInfoList <- list() + declFunClassList <- list() + declFunNames <- names(mDef$declFunNameToIndex) + for(i in seq_along(mDef$declInfo)) { + declInfo <- mDef$declInfo[[i]] + decl_methods <- make_decl_methods_from_declInfo(declInfo) + declVars <- decl_methods |> lapply(\(x) all.vars(body(x))) |> unlist() |> unique() |> setdiff(c("idx", "LocalNewLogProb_", "LocalAns_", "model")) %||% character() + declVarInfo <- modelVarInfo$vars[declVars] + declID <- as.numeric(declInfo$declRule$ID) # Formerly `sourceLineNumber`, which may not be unique. + declFun_membername <- declFunNames[i] + declFun_classname <- sub("declFun", "declFunClass", declFun_membername) # name of an nClass generator + 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)) + 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. + # so we check here. + ordered_decl_names <- lapply(declInfoList, function(x) x$membername) |> unlist() + if(!identical(ordered_decl_names, names(mDef$declFunNameToIndex))) + stop("declaration ordering in declInfoList does not matchdeclFunNameToIndex") + modelClassInstance <- makeModel_nClass(modelVarInfo, declInfoList, inits = m$origInits, data = m$origData, + model = m, classname = "my_model", env = environment()) +} + +## The two "addModelDollarSign" functions are borrowed directly from nimble. +## This should add model$ in front of any names that are not already part of a '$' expression + +nm_addModelDollarSign <- function(expr, exceptionNames = character(0)) { + if(is.numeric(expr)) return(expr) + if(is(expr, 'srcref')) return(expr) + if(is.name(expr)) { + if((as.character(expr) %in% exceptionNames) || (as.character(expr) == '')) return(expr) + proto <- quote(model$a) + proto[[3]] <- expr + return(proto) + } + if(is.call(expr)) { + if(expr[[1]] == '$'){ + expr[2] <- lapply(expr[2], function(listElement) nm_addModelDollarSign(listElement, exceptionNames)) + return(expr) + } + if(expr[[1]] == 'returnType') + return(expr) + if(length(expr) > 1) { + expr[2:length(expr)] <- lapply(expr[-1], function(listElement) nm_addModelDollarSign(listElement, exceptionNames)) + return(expr) + } + } + return(expr) +} + +# Turn variables and methods into a declFun nClass +make_declFun_nClass <- function(varInfo = list(), + methods = list(), + classname, + declID) { + # 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. + CpublicVars <- varInfo |> lapply(\(x) paste0("ref(double(", x$nDim ,", interface=FALSE))")) + names(CpublicVars) <- varInfo |> lapply(\(x) x$name) |> unlist() + + +# varInfo_2_symbol <- \(x) nCompiler:::symbolBasic$new( +# type="double", nDim=x$nDim, name="", isRef=TRUE, isConst=FALSE, interface=FALSE) # In future maybe isConst=TRUE, but it might not matter much +# symbolList <- varInfo |> lapply(varInfo_2_symbol) +# names(symbolList) <- varInfo |> lapply(\(x) x$name) |> unlist() + numVars <- length(varInfo) + +# CpublicVars <- names(symbolList) |> lapply(\(x) eval(substitute(quote(T(symbolList$NAME)), +# list(NAME=as.name(x))))) +# names(CpublicVars) <- names(symbolList) + # This is a kluge to have a model field in the Cpublic_obj, + # needed for uncompiled purposes, and for compiled purposes + # we instead use references to model variables. So + # the declared type here is arbitrary. + initFun <- function(){} + + if(numVars > 0) { + # ctorArgNames <- paste0(names(symbolList), '_') + ctorArgNames <- paste0(names(CpublicVars), '_') + # List used when generating C++ constructor code to allow direct initializers, necessary for references. + # initializersList <- paste0(names(symbolList), '(', ctorArgNames ,')') + initializersList <- paste0(names(CpublicVars), '(', ctorArgNames ,')') + formals(initFun) <- structure(as.pairlist(CpublicVars), names = ctorArgNames) + } else { + initializersList <- character() + } + ## TODO: I don't think this labelCreator (or the one for the model) exist (though they shouldn't be used...) + if(missing(classname)) + classname <- declLabelCreator() + + baseclass <- paste0("declFunClass_<", classname, ">") + + # Rpublic method to set the model pointer/reference. + setModel <- function(model) { + if(!isCompiled()) { + self$model <- model + #private$Cpublic_obj$model <- model + } + else + warning("setModel called on compiled object; no action taken") + } + +# This was a prototype + # Actually, we are using this. Ok? // CJP + declFun_nClass <- substitute( + nClass( + inherit = declFunBase_nClass, + classname = CLASSNAME, + Rpublic = RPUBLIC, + Cpublic = CPUBLIC, + compileInfo = list( + createFromR = FALSE, # Without a default constructor (which we've disabled here), createFromR is impossible + nClass_inherit = list(base = BASECLASS)) # Ideally this line would be obtained from a base nClass, but we insert it directly for now + ), + list( + CPUBLIC = c( + declID = declID, + list( + nFunction( + initFun, + compileInfo = list(constructor=TRUE, initializers = initializersList) + ) + ) |> structure(names = classname), + CpublicVars, + list(model = "RcppList"), + methods + ), + RPUBLIC = list(#model = NULL, + setModel = setModel), + CLASSNAME = classname, + BASECLASS = baseclass + )) + eval(declFun_nClass) +} +#test <- nCompiler:::type2symbol('CppVar(baseType = type2cpp("numericVector"), ref=TRUE, const=TRUE)') + +# Make all the info needed to include a decl in a model class. +# The decl_nClass should be created first. +# Currently it needs to have a name to include in nCompile(). Later we might be able to pass the object itself +# At first drafting this is fairly trivial but could grow in complexity. + +make_decl_info_for_model_nClass <- function(membername, + declFunName, + classname, + varInfo = list() + ) { + ctorArgs <- varInfo |> lapply(\(x) x$name) |> unlist() + + list(declFunName = declFunName, + membername = membername, + classname = classname, + ctorArgs = ctorArgs) +} + +makeModel_nClass <- function(modelVarInfo, + decls = list(), + classname, + inits = list(), + data = list(), + model = NULL, + env = parent.frame() + ) { + ## varInfo will be a list (names not used) of name, nDim, sizes. + CpublicModelVars <- modelVarInfo$vars |> lapply(\(x) paste0("numericArray(nDim=",x$nDim,")")) + names(CpublicModelVars) <- modelVarInfo$vars |> lapply(\(x) x$name) |> unlist() + opDefs <- list( + base_ping = nCompiler:::getOperatorDef("custom_call"), + setup_auto_decl_mgmt = nCompiler:::getOperatorDef("custom_call"), + do_setup_decl_mgmt_from_names = nCompiler:::getOperatorDef("custom_call") + ) + opDefs$base_ping$returnType <- nCompiler:::type2symbol(quote(void())) # How can this be passed into nClass? + opDefs$base_ping$labelAbstractTypes$recurse <- FALSE + opDefs$setup_auto_decl_mgmt$returnType <- nCompiler:::type2symbol(quote(void())) + opDefs$setup_auto_decl_mgmt$labelAbstractTypes$recurse <- FALSE + opDefs$do_setup_decl_mgmt_from_names$returnType <- nCompiler:::type2symbol(quote(void())) + opDefs$do_setup_decl_mgmt_from_names$labelAbstractTypes$recurse <- FALSE + + if(missing(classname)) + classname <- modelLabelCreator() + + CpublicMethods <- list( + do_setup_auto_decl_mgmt = nFunction( + name = "call_setup_auto_decl_mgmt", + function() {}, + compileInfo=list( + C_fun = function() {setup_auto_decl_mgmt()}) + ), + setup_decl_mgmt_from_names = nFunction( + name = "call_setup_decl_mgmt_from_names", + function(declNames) {}, + compileInfo=list( + C_fun = function(declNames="RcppCharacterVector") {do_setup_decl_mgmt_from_names(declNames)}) + ), + # print_decls = nFunction( + # name = "print_decls", + # function() {}, + # compileInfo=list( + # C_fun = function() {cppLiteral('modelClass_::c_print_decls();')}) + # ), + set_from_list = nFunction( + name = "set_from_list", + function(Rlist) {for(v in names(Rlist)) + if(exists(v, self, inherits=FALSE)) self[[v]] <- Rlist[[v]]}, + compileInfo=list( + C_fun=function(Rlist = 'RcppList') {cppLiteral('modelClass_::set_from_list(Rlist);')}) + ), + resize_from_list = nFunction( + name = "resize_from_list", + function(Rlist) {for(v in names(Rlist)) + if(exists(v, self, inherits=FALSE)) self[[v]] <- nArray(value = NA, dim=Rlist[[v]])}, + compileInfo = list( + C_fun=function(Rlist = 'RcppList') {cppLiteral('modelClass_::resize_from_list(Rlist);')}) + ) + ) + # decls will be a list of membername, declName, (decl) classname, ctorArgs (list) + decl_pieces <- decls |> lapply(\(x) { + #nClass_type <- paste0(x$declFxnName, "()") + init_string <- paste0('nCpp("', x$membername, '( new ', x$classname, '(', + paste0(x$ctorArgs, collapse=","), '))")') + list(nClass_type = x$declFunName, + init_string = init_string) + }) + + declFunNameToIndex <- model$modelDef$declFunNameToIndex + + CpublicDeclFuns <- decl_pieces |> lapply(\(x) x$nClass_type) |> setNames(names(declFunNameToIndex)) + # CpublicDeclFuns <- list( + # beta_decl = 'decl_dnorm()' + # ) + CpublicCtor <- list( + nFunction( + function(){ + cppLiteral("setup_decl_mgmt();") # This will be the default but can be overridden by decls that need to do something special. We could also have a version that takes decl names as input and only sets up those. + }, + compileInfo = list(constructor=TRUE, + #initializers = c('nCpp("beta_decl(new decl_dnorm(mu, beta, 1))")')) + initializers = decl_pieces |> lapply(\(x) x$init_string) |> unlist()) + ) + ) |> structure(names = classname) + + declFunPtrsSetupLiterals <- paste0("declFunPtrs[(", as.integer(declFunNameToIndex) , ")-1] = ", names(declFunNameToIndex)) + declFunPtrsResizeLiteral <- paste0("declFunPtrs.resize(", length(declFunNameToIndex) , ")") + setup_decl_mgmt_body <- as.list(c(declFunPtrsResizeLiteral, declFunPtrsSetupLiterals)) |> + lapply(\(x) substitute(nCpp(X), list(X = x))) + setup_decl_mgmt_fun <- function() {} + for(i in seq_along(setup_decl_mgmt_body)) + body(setup_decl_mgmt_fun)[[i+1]] <- setup_decl_mgmt_body[[i]] + Cpublic_setup_decl_mgmt <- list(setup_decl_mgmt = nFunction(name = "setup_decl_mgmt", fun = setup_decl_mgmt_fun)) + + baseclass <- paste0("modelClass_<", classname, ">") + # CpublicDeclFuns has elements like "decl_1 = quote(declFxn_1())" + # We provide it in Cpublic to declare C++ member variables with types. + # We also place the list itself in the class so that we can look up for uncompiled execution + # the objects that need to be created in initialize. + # If we someday make type declarations and initializations more automatic, we can avoid this duplication. + ans <- substitute( + nClass( + classname = CLASSNAME, + inherit = modelBase_nClass, + compileInfo = list(opDefs = OPDEFS, + nClass_inherit = list(base=BASECLASS) #, +# needed_units = list("declFxnBase_nClass"), # needed for package=TRUE +# Hincludes = '"declFxnBase_nClass_c_.h"' # needed for package=TRUE + ), + Rpublic = RPUBLIC, + Cpublic = CPUBLIC, + env = env + ), + list(OPDEFS = opDefs, + # A list of individual elements + RPUBLIC = list( + declFunNameToIndex_ = model$modelDef$declFunNameToIndex, + defaultSizes = modelVarInfo$sizes, + defaultInits = inits, + defaultData = data, + modelDef = model$modelDef, + dataRules = model$dataRules, + nondataRules = model$nondataRules, + predictiveRules = model$predictiveRules, + nonpredictiveRules = model$nonpredictiveRules, + CpublicDeclFuns = CpublicDeclFuns), + # A concatenation of lists + CPUBLIC = c(CpublicDeclFuns, Cpublic_setup_decl_mgmt, CpublicModelVars, CpublicCtor, CpublicMethods), + CLASSNAME = classname, + BASECLASS = baseclass) + ) + eval(ans) +} + +## Get varInfo from new nimbleModel +get_varInfo_from_nimbleModel <- function(model) { + mDef <- model$modelDef + extract <- \(x) x |> lapply(\(x) list(name = x$varName, nDim = x$nDim)) + vars <- mDef$varInfo |> extract() + logProbVars <- mDef$logProbVarInfo |> extract() + # The resize_from_list method will error out if a scalar is included. + # The maxs is empty for scalars, so they are automatically omitted from the sizes result here. + # TODO: CJP sees scalars included as numeric(0) in sizes, so not omitted. Will this be a problem for resize_from_list? + # TODO: If ok, put sizes info into the same list as vars info. + extract_sizes <- \(x) x|> lapply(\(x) x$maxs) + sizes <- mDef$varInfo |> extract_sizes() + logProb_sizes <- mDef$logProbVarInfo |> extract_sizes() + list( + vars = c(vars, logProbVars), + sizes = c(sizes, logProb_sizes) + ) +} + +# make_stoch_calculate <- function(LHSrep, RHSrep, logProbExprRep) { +# lenRHS <- length(RHSrep) +# if(length(RHS) > 1) { +# RHSrep[3:(lenRHS+1)] <- RHSrep[2:lenRHS] +# names(RHSrep)[3:(lenRHS+1)] <- names(RHSrep)[2:lenRHS] +# } +# RHSrep[[2]] <- LHSrep +# names(RHSrep)[2] <- "" +# RHSrep[[lenRHS+2]] <- 1 +# names(RHSrep)[lenRHS+2] <- "log" +# # We create separate code for R and C execution. +# calc1Cfun <- substitute( +# function(idx) {LHS <- RHS; return(LHS)}, +# list(LHS = logProbExprRep, RHS = RHSrep) +# ) |> eval() +# make_calculate_from_Cfun(calc1Cfun) +# } + +make_stoch_sim_line <- function(LHSrep, RHSrep) { + BUGSdistName <- nCompiler:::safeDeparse(RHSrep[[1]]) + distInfo <- getDistributionInfo(BUGSdistName) + sim_code <- as.name(distInfo$simulateName) + if(is.null(sim_code)) stop("Could not find simulation ('r') function for ", BUGSdistName) + RHSrep[[1]] <- sim_code + # scoot all named arguments right 1 position + if(length(RHSrep) > 1) { + for(i in (length(RHSrep)+1):3) { + RHSrep[i] <- RHSrep[i-1] + names(RHSrep)[i] <- names(RHSrep)[i-1] + } + } + RHSrep[[2]] <- 1 + names(RHSrep)[2] <- '' + sim_line <- substitute( + LHS <- RHS, + list(LHS = LHSrep, RHS = RHSrep)) + sim_line +} + +make_stoch_calc_line <- function(LHSrep, RHSrep, logProbExprRep, diff = FALSE) { + lenRHS <- length(RHSrep) + if(length(RHSrep) > 1) { + RHSrep[3:(lenRHS+1)] <- RHSrep[2:lenRHS] + names(RHSrep)[3:(lenRHS+1)] <- names(RHSrep)[2:lenRHS] + } + RHSrep[[2]] <- LHSrep + names(RHSrep)[2] <- "" + RHSrep[[lenRHS+2]] <- 1 + names(RHSrep)[lenRHS+2] <- "log" + # We create separate code for R and C execution. + if(!diff) { + calc_line <- substitute( + LHS <- RHS, + list(LHS = logProbExprRep, RHS = RHSrep)) + } else { + calc_line <- substitute( + LocalNewLogProb_ <- RHS, + list(RHS = RHSrep)) + } + calc_line +} + +make_determ_calc_line <- function(LHSrep, RHSrep) { + calc_line <- substitute( + LHS <- RHS, + list(LHS = LHSrep, RHS = 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_decl_method_nFxn <- function(f, name, returnType='numericScalar') { + Cfun <- f + Rfun <- f + body(Rfun) <- nm_addModelDollarSign(body(f), exceptionNames = c("idx", "LocalNewLogProb_", "LocalAns_")) + if(is.null(returnType)) returnType <- 'void' + nFxn <- nFunction( + name = name, + fun = Rfun, + argTypes = list(idx = 'integerVector'), + returnType = T(returnType), + compileInfo=list(C_fun=Cfun), + ) + nFxn +} + +make_decl_methods_from_declInfo <- function(declInfo) { + # pieces are adapted from Chris' code in nimbleModel and/or old nimble. + # + # This function creates a calc_one nFunction that calculates single index case. + # This will then be used by generic iterator over indices. + # Vectorized cases can be added in this basic framework later. + modelCode <- declInfo$calculateCode + LHS <- modelCode[[2]] + RHS <- modelCode[[3]] + type <- if(modelCode[[1]]=="~") "stoch" else "determ" # or use declInfo$stoch (logical) + context <- declInfo$declRule$context + replacements <- sapply(seq_along(context$singleContexts), + function(i) parse(text = paste0('idx[',i,']'))[[1]]) + names(replacements) <- context$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 = (function(idx) {calc_one(idx)}) |> + make_decl_method_nFxn("sim_one", NULL), + calc_one = (function(idx) {DETERMCALC; return(invisible(0))}) |> + make_decl_method_nFxn("calc_one"), + calcDiff_one = (function(idx) {calc_one(idx);return(invisible(0))}) |> + make_decl_method_nFxn("calcDiff_one"), + getLogProb_one = (function(idx) {return(0)}) |> + make_decl_method_nFxn("getLogProb_one") + ), + list(DETERMCALC = make_determ_calc_line(LHSrep, RHSrep)) + )) + } + if(type == 'stoch') { + logProbExpr <- declInfo$genLogProbExpr() + logProbExprRep <- eval(substitute(substitute(e, replacements), list(e = logProbExpr))) + methodList <- eval(substitute( + list( + sim_one = (function(idx) { STOCHSIM }) |> + make_decl_method_nFxn("sim_one", NULL), + calc_one = (function(idx) { STOCHCALC; return(invisible(LOGPROB)) }) |> + make_decl_method_nFxn("calc_one"), + calcDiff_one = (function(idx) {STOCHCALC_DIFF; LocalAns_ <- LocalNewLogProb_ - LOGPROB; + LOGPROB <- LocalNewLogProb_; return(invisible(LocalAns_))}) |> + make_decl_method_nFxn("calcDiff_one"), + getLogProb_one = (function(idx) { return(LOGPROB) }) |> + make_decl_method_nFxn("getLogProb_one") + ), + list( LOGPROB = logProbExprRep, + STOCHSIM = make_stoch_sim_line(LHSrep, RHSrep), + STOCHCALC = make_stoch_calc_line(LHSrep, RHSrep, logProbExprRep), + STOCHCALC_DIFF = make_stoch_calc_line(LHSrep, RHSrep, logProbExprRep, diff=TRUE)) + )) + } + methodList +} diff --git a/nimbleModel/R/nodeRules.R b/nimbleModel/R/nodeRules.R index c062875..f1f51f8 100644 --- a/nimbleModel/R/nodeRules.R +++ b/nimbleModel/R/nodeRules.R @@ -7,18 +7,12 @@ ## rather it simply represents the valid index values for a variable. ## There is one declRule for each declaration in a model, representing the indexing for the LHS variable. -## A declRule contains the `calculate` function (i.e., the nodeFunction) that operates on a -## single set of index values. ## calcRules are generated by starting with declRules and then fracturing (with `fracture`) ## based on top-down processing. This produces a calcRule for each set of nodes ## from a declaration that can be calculated together (same sortID) ## (as well as the special case of state-space model type formulations). -## TODO: work on `simulate`, `getParam` and other related model methods; will these be part or -## `calcRules`? - - ## Base class for all node-related classes nodeRuleClass <- R6Class( classname = "nodeRuleClass", @@ -117,8 +111,7 @@ nodeRuleClass <- R6Class( ) ) -## Class for representing nodes at the declaration level, containing -## calculate and other nodeFunctions. +## Class for representing nodes at the declaration level. declRuleClass <- R6Class( classname = "declRuleClass", portable = FALSE, @@ -126,65 +119,13 @@ declRuleClass <- R6Class( public = list( originalIndexingRule = NULL, # Determines original indexing (based on context). decl = NULL, # Full declInfo; nodeRuleClass$expr is just LHS. - calculate = NULL, # Generic function for calculation. - ## TODO: remove this: - test = rep(0, 10), ## test element used for testing calculation while modelDef doesn't have acccess to a model - test2 = matrix(0, 3, 5), - initialize = function(decl, ID, context = modelContextClass$new(), constants = list()) { decl <<- decl super$initialize(decl$code[[2]], ID, context = context, constants = constants) ## `expr` in is parent class. originalIndexingRule <<- originalIndexingRuleClass$new(expr, context, constants) - }, - - buildFunctions = function(code, logProbExpr) { - buildCalculateFun(code, logProbExpr, context) - }, - - buildCalculateFun = function(code, logProbExpr, context) { - newCode <- code - if(newCode[[1]] == '<-') { - newCode <- quote(A <<- B) - newCode[2:3] <- code[2:3] - } - replacements <- sapply(seq_along(context$singleContexts), - function(i) parse(text = paste0('idx[',i,']'))[[1]]) - names(replacements) <- context$indexVarNames - - for(i in seq_along(context$singleContexts)) { - newCode <- eval(substitute(substitute(e, replacements), list(e = newCode))) - logProbExpr <- eval(substitute(substitute(e, replacements), list(e = logProbExpr))) - } - - if(decl$stoch) { - ## Insert 'logProb_' and change to assignment, moving LHS in as first argument. - finalCode <- quote(A <<- B) - finalCode[[2]] <- logProbExpr - finalCode[[3]] <- newCode[[3]] - len <- length(finalCode[[3]]) - if(len > 1) { - finalCode[[3]][3:(len+1)] <- finalCode[[3]][2:len] - names(finalCode[[3]])[3:(len+1)] <- names(finalCode[[3]])[2:len] - } - finalCode[[3]][[2]] <- newCode[[2]] - names(finalCode[[3]])[2] <- "" - finalCode[[3]][[len+2]] <- 1 - names(finalCode[[3]])[len+2] <- "log" - calculate <<- function(idx) { - ## logProb_y <- array(0, rep(100, nvals)) # TODO: placeholder so logProb storage exists for testing - } - ## body(calculate)[[length(body(calculate))+1]] <<- finalCode - body(calculate) <<- finalCode - } else { - calculate <<- function(idx) {} - body(calculate) <<- newCode - } - ## TODO: will need to deal with logProb for mv node having single value inserted. - ## TODO: will need to deal with the various complexities we currently deal with - alt params, truncation, etc. } - ) ) @@ -239,7 +180,7 @@ calcRuleClass <- R6Class( indexingRange <- declRule$originalIndexingRule$apply(inputRange) if(is.null(indexingRange)) return(NULL) - result <- calcRangeClass$new(varName, indexingRange, declRule$calculate, sortID, multiSortIDindex) + result <- calcRangeClass$new(varName, indexingRange, as.numeric(declRule$ID), sortID, multiSortIDindex) return(result) }, @@ -385,7 +326,7 @@ calcRuleClass <- R6Class( ) ## calcRanges manage the calculation for one or more nodes, handling the indexing, and -## calling out to the declRule `calculate` function. +## calling out to the declFun `calculate` function. calcRangeClass <- R6Class( classname = "calcRangeClass", portable = FALSE, @@ -395,75 +336,15 @@ calcRangeClass <- R6Class( sortID = NULL, calcFun = NULL, multiSortIDindex = NULL, - initialize = function(varName, indexingRange, calcFun, sortID, multiSortIDindex) { + declID = NULL, + initialize = function(varName, indexingRange, declID, sortID, multiSortIDindex) { varName <<- varName indexingRange <<- indexingRange calcFun <<- calcFun ## note that calcFun itself is not vectorized sortID <<- sortID + declID <<- declID multiSortIDindex <<- multiSortIDindex - }, - - ## Generic calculate function that crosses the indexRanges in the indexingRange (a varRange) - ## and extracts the original indices to feed into calculate nodeFunction - ## that operates on set of scalar indices. - - ## Keep indexing internal to the indexRange to avoid complicated and possibly repetitive - ## calculation of internal indexing. - - ## Will need to figure out how this is going to get compiled. - ## This will rely on nCompiler indexing of eigen tensors for static block indexes, e.g. '3:5' in x[i, 3:5] - ## which will presumably use the information in the symbolicParentNodes. - - calculate = function() { - numRanges <- length(indexingRange$indexRanges) - if(!numRanges) { # no indexing - result <- calcFun(NULL) - } else { - indexRange_lengths <- sapply(indexingRange$indexRanges, - function(x) x$numElements) - indexPositions <- indexingRange$rangeToIndexSlot - len <- prod(indexRange_lengths) - - ## TODO: This is a placeholder so we can test numerical results - ## once fuller workflow is in place, remove this and assignment of output of calcFun() - result <- rep(0, len) - - ## Set up information so `getNext` repeats index values for outer loop indices. - delay <- 1 - for(irIndex in rev(seq_len(numRanges))) { # work from inner-most loop outwards - indexingRange$indexRanges[[irIndex]]$setDelay(delay) - delay <- delay * indexingRange$indexRanges[[irIndex]]$numElements - } - - if(length(sortID) == 1 || len == 1) { - index <- numeric(length(indexingRange$indexSlotToRange)) ## vector to hold the original index values - for(item in seq_len(len)) { - for(irIndex in seq_len(numRanges)) - index[indexPositions[[irIndex]]] <- indexingRange$indexRanges[[irIndex]]$getNext() - result[item] <- calcFun(index) ## scalar calculation - } - } else { - index <- matrix(0, len, length(indexingRange$indexSlotToRange)) ## vector to hold the original index values - sortIDvals <- rep(0, len) - for(item in seq_len(len)) { - for(irIndex in seq_len(numRanges)) { - index[item, indexPositions[[irIndex]]] <- indexingRange$indexRanges[[irIndex]]$getNext() - } - sortIDvals[item] <- sortID[index[item, multiSortIDindex]] - } - for(item in order(sortIDvals)) - result[item] <- calcFun(index[item, ]) - } - } - ## Note that result will be vectorized based on the loop ordering so in simple rectangular - ## settings such as `y[i,j]` with `i` the outer index, it will be row-major. - return(result) - } - - ## simulate() should be very similar to calculate() - ## need to think more about getParam, getBound, etc., but probably also very similar, - ## just swapping out calcFun() for appropriate function. ) ) diff --git a/nimbleModel/R/originalIndexingRules.R b/nimbleModel/R/originalIndexingRules.R index dfba406..01034df 100644 --- a/nimbleModel/R/originalIndexingRules.R +++ b/nimbleModel/R/originalIndexingRules.R @@ -16,12 +16,13 @@ originalIndexingRuleClass <- R6Class( constants = list()) { varName <<- getVarName(LHS) if(length(context$indexVarNames)) { - ## Exclude indices not used in lifted expression, e.g., `i` in `y[i,j] ~ dnorm(mu[i], sigma[j])` + ## 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)] - dummyLHS <- parse(text = paste0(varName, "[", - paste(indexVarNames, collapse = ","), - "]"))[[1]] + indexing <- if(length(indexVarNames)) + paste0("[", paste(indexVarNames, collapse = ","), "]") else "" + dummyLHS <- parse(text = paste0(varName, indexing))[[1]] + ## Unused singleContexts will be removed in graphRuleClass$new(). } else dummyLHS <- as.name(varName) graphRule <<- graphRuleClass$new(dummyLHS, LHS, diff --git a/nimbleModel/inst/include/nimbleModel/.DS_Store b/nimbleModel/inst/include/nimbleModel/.DS_Store new file mode 100644 index 0000000..c4c1caa Binary files /dev/null and b/nimbleModel/inst/include/nimbleModel/.DS_Store differ diff --git a/nimbleModel/inst/include/nimbleModel/predef/.DS_Store b/nimbleModel/inst/include/nimbleModel/predef/.DS_Store new file mode 100644 index 0000000..dead564 Binary files /dev/null and b/nimbleModel/inst/include/nimbleModel/predef/.DS_Store differ diff --git a/nimbleModel/inst/include/nimbleModel/predef/declFunBase_nC/declFunBase_nClass_copyFiles.txt b/nimbleModel/inst/include/nimbleModel/predef/declFunBase_nC/declFunBase_nClass_copyFiles.txt new file mode 100644 index 0000000..e69de29 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 new file mode 100644 index 0000000..fdc948f --- /dev/null +++ b/nimbleModel/inst/include/nimbleModel/predef/declFunBase_nC/declFunBase_nClass_cppContent.cpp @@ -0,0 +1,67 @@ +/* OPENER (Do not edit this comment) */ +#ifndef __declFunBase_nClass_CPP +#define __declFunBase_nClass_CPP +/* BODY (Do not edit this comment) */ +#ifndef R_NO_REMAP +#define R_NO_REMAP +#endif +#include +#include "declFunBase_nClass_c_.h" +using namespace Rcpp; +// [[Rcpp::plugins(nCompiler_Eigen_plugin)]] +// [[Rcpp::depends(RcppParallel)]] +// [[Rcpp::depends(nCompiler)]] +// [[Rcpp::depends(Rcereal)]] +// [[Rcpp::depends(nimbleModel)]] + + bool declFunBase_nClass::ping ( ) { +RESET_EIGEN_ERRORS +return(true); +} + double declFunBase_nClass::calculate_cpp ( std::shared_ptr instr ) { +RESET_EIGEN_ERRORS +Rprintf("declFunBase_nClass virtual base calculate_cpp should never be called (something is wrong)\n");; +return(0.0); +} + double declFunBase_nClass::calculateDiff_cpp ( std::shared_ptr instr ) { +RESET_EIGEN_ERRORS +Rprintf("declFunBase_nClass virtual base calculateDiff_cpp should never be called (something is wrong)\n");; +return(0.0); +} + double declFunBase_nClass::getLogProb_cpp ( std::shared_ptr instr ) { +RESET_EIGEN_ERRORS +Rprintf("declFunBase_nClass virtual base getLogProb_cpp should never be called (something is wrong)\n");; +return(0.0); +} + void declFunBase_nClass::simulate_cpp ( std::shared_ptr instr ) { +RESET_EIGEN_ERRORS +Rprintf("declFunBase_nClass virtual base simulate_cpp should never be called (something is wrong)\n");; +} + declFunBase_nClass::declFunBase_nClass ( ) { +RESET_EIGEN_ERRORS +} + +// [[Rcpp::export(name = "set_CnClass_env_declFunBase_nClass_new")]] + void set_CnClass_env_declFunBase_nClass ( SEXP env ) { +RESET_EIGEN_ERRORS +SET_CNCLASS_ENV(declFunBase_nClass, env);; +} + +// [[Rcpp::export(name = "get_CnClass_env_declFunBase_nClass_new")]] + Rcpp::Environment get_CnClass_env_declFunBase_nClass ( ) { +RESET_EIGEN_ERRORS +return GET_CNCLASS_ENV(declFunBase_nClass);; +} + +NCOMPILER_INTERFACE( +declFunBase_nClass, +NCOMPILER_FIELDS(), +NCOMPILER_METHODS( +method("ping", &declFunBase_nClass::ping, args({{}})), +method("calculate_cpp", &declFunBase_nClass::calculate_cpp, args({{arg("instr",copy)}})), +method("calculateDiff_cpp", &declFunBase_nClass::calculateDiff_cpp, args({{arg("instr",copy)}})), +method("getLogProb_cpp", &declFunBase_nClass::getLogProb_cpp, args({{arg("instr",copy)}})), +method("simulate_cpp", &declFunBase_nClass::simulate_cpp, args({{arg("instr",copy)}})) +) +) +#endif diff --git a/nimbleModel/inst/include/nimbleModel/predef/declFunBase_nC/declFunBase_nClass_filebase.txt b/nimbleModel/inst/include/nimbleModel/predef/declFunBase_nC/declFunBase_nClass_filebase.txt new file mode 100644 index 0000000..75d0892 --- /dev/null +++ b/nimbleModel/inst/include/nimbleModel/predef/declFunBase_nC/declFunBase_nClass_filebase.txt @@ -0,0 +1 @@ +declFunBase_nClass_c_ 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 new file mode 100644 index 0000000..fe06f27 --- /dev/null +++ b/nimbleModel/inst/include/nimbleModel/predef/declFunBase_nC/declFunBase_nClass_hContent.h @@ -0,0 +1,27 @@ +/* OPENER (Do not edit this comment) */ +#ifndef __declFunBase_nClass_H +#define __declFunBase_nClass_H +/* BODY (Do not edit this comment) */ +#ifndef R_NO_REMAP +#define R_NO_REMAP +#endif +#include +#include "instr_nClass_c_.h" + +class declFunBase_nClass : public interface_resolver< genericInterfaceC >, public loadedObjectHookC { +public: + virtual bool ping ( ) ; + virtual double calculate_cpp ( std::shared_ptr instr ) ; + virtual double calculateDiff_cpp ( std::shared_ptr instr ) ; + virtual double getLogProb_cpp ( std::shared_ptr instr ) ; + virtual void simulate_cpp ( std::shared_ptr instr ) ; + declFunBase_nClass ( ) ; +}; + + void set_CnClass_env_declFunBase_nClass ( SEXP env ) ; + + Rcpp::Environment get_CnClass_env_declFunBase_nClass ( ) ; + +#include + +#endif 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 new file mode 100644 index 0000000..72394be --- /dev/null +++ b/nimbleModel/inst/include/nimbleModel/predef/declFunBase_nC/declFunBase_nClass_manifest.txt @@ -0,0 +1,7 @@ +list(saved_at = structure(1780285145.92737, class = c("POSIXct", +"POSIXt")), packet_name = "declFunBase_nClass", elements = c("preamble", +"cppContent", "hContent", "filebase", "post_cpp_compiler", "copyFiles" +), files = list(preamble = "declFunBase_nClass_preamble.cpp", + cppContent = "declFunBase_nClass_cppContent.cpp", hContent = "declFunBase_nClass_hContent.h", + filebase = "declFunBase_nClass_filebase.txt", post_cpp_compiler = "declFunBase_nClass_post_cpp_compiler.txt", + copyFiles = "declFunBase_nClass_copyFiles.txt")) diff --git a/nimbleModel/inst/include/nimbleModel/predef/declFunBase_nC/declFunBase_nClass_post_cpp_compiler.txt b/nimbleModel/inst/include/nimbleModel/predef/declFunBase_nC/declFunBase_nClass_post_cpp_compiler.txt new file mode 100644 index 0000000..e69de29 diff --git a/nimbleModel/inst/include/nimbleModel/predef/declFunBase_nC/declFunBase_nClass_preamble.cpp b/nimbleModel/inst/include/nimbleModel/predef/declFunBase_nC/declFunBase_nClass_preamble.cpp new file mode 100644 index 0000000..1a92ce3 --- /dev/null +++ b/nimbleModel/inst/include/nimbleModel/predef/declFunBase_nC/declFunBase_nClass_preamble.cpp @@ -0,0 +1,5 @@ +#define NCOMPILER_HANDLE_EIGEN_ERRORS +#define NCOMPILER_USES_EIGEN +#define NCOMPILER_USES_NCPPVEC +#define USES_NCOMPILER +#define NCOMPILER_USES_NCLASS_INTERFACE diff --git a/nimbleModel/inst/include/nimbleModel/predef/declFunClass_/declFunClass_.h b/nimbleModel/inst/include/nimbleModel/predef/declFunClass_/declFunClass_.h new file mode 100644 index 0000000..72fa4f2 --- /dev/null +++ b/nimbleModel/inst/include/nimbleModel/predef/declFunClass_/declFunClass_.h @@ -0,0 +1,249 @@ +// to be included from the predefined nodeFxnBase_nClass. +// Add "#include " to that file, +// after the declaration of nodeFxnBase_nClass. + +template +class declFunClass_ : public declFunBase_nClass { +public: + double v; + declFunClass_() {}; + + double calculate_cpp( std::shared_ptr instr) override { + return calc_op_< &Derived::calc_one >(instr); + } + double calculateDiff_cpp( std::shared_ptr instr) override { + return calc_op_< &Derived::calcDiff_one >(instr); + } + double getLogProb_cpp( std::shared_ptr instr) override { + return calc_op_< &Derived::getLogProb_one >(instr); + } + template + double calc_op_ ( std::shared_ptr instr ) { + RESET_EIGEN_ERRORS; + int instr_type = 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_1_matp_ord_< Method >(instr); + if(instr_type == 5) return calc_2_seq_seq_< Method >(instr); + if(instr_type == 6) return calc_2_seq_seq_ord_< Method >(instr); + return(0); + } + template + double calc_0_ (std::shared_ptr instr) { + return( (static_cast(this)->*Method)(instr->lens) ); // lens serves as a dummy here, to have the right type to pass + } + template + double calc_1_seq_(std::shared_ptr instr) { + int len = instr->lens[0]; + if(len < 1) return(0); + int iStart = instr->values->operator[](0)[0]; + int iEnd = iStart + len; + Eigen::Tensor idx(1); + double logProb(0.); + for(int i = iStart; i < iEnd; ++i) { + idx[0] = i; + logProb += (static_cast(this)->*Method)(idx); + } + return(logProb); + } + template + double calc_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 calc_1_mat_"< idx(1); + double logProb(0.); + for(int i = 0; i < len; ++i) { + idx[0] = vals[i]; + logProb += (static_cast(this)->*Method)(idx); + } + return(logProb); + } + template + double calc_1_matp_(std::shared_ptr instr) { + int len = instr->lens[0]; + int dm = instr->dims[0]; + const auto& vals = instr->values->operator[](0); + if(len*dm != vals.size()) std::cout<<"len*dm != vals.size() in calc_1_matp_"< idx(dm); + double logProb(0.); + for(int i = 0; i < len; ++i) { + for(int p = 0; p < dm; ++p) + idx[p] = vals[i*dm+p]; + logProb += (static_cast(this)->*Method)(idx); + } + return(logProb); + } + template + double calc_1_matp_ord_(std::shared_ptr instr) { + int len = instr->lens[0]; + int dm = instr->dims[0]; + const auto& vals = instr->values->operator[](0); + if(len*dm != vals.size()) std::cout<<"len*dm != vals.size() in calc_1_matp_"< idx(dm); + Eigen::Tensor slots(dm); + for(int p = 0; p < dm; ++p) + slots[p] = instr->slots[p]-1; + double logProb(0.); + for(int i = 0; i < len; ++i) { + for(int p = 0; p < dm; ++p) + idx[slots[p]] = vals[i*dm+p]; + logProb += (static_cast(this)->*Method)(idx); + } + return(logProb); + } + template + double calc_2_seq_seq_(std::shared_ptr instr) { + int len1 = instr->lens[0]; + int len2 = instr->lens[1]; + if(len1 < 1) return(0); + if(len2 < 1) return(0); + 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); + double logProb(0.); + for(int i = iStart1; i < iEnd1; ++i) { + idx[0] = i; + for(int j = iStart2; j < iEnd2; ++j) { + idx[1] = j; + logProb += (static_cast(this)->*Method)(idx); + } + } + return(logProb); + } + template + double calc_2_seq_seq_ord_(std::shared_ptr instr) { + if(instr->slots[0] != 2 || instr->slots[1] != 1) + std::cout<<"slots not equal to 2,1 in calc_2_seq_seq_ord_"<lens[0]; + int len2 = instr->lens[1]; + if(len1 < 1) return(0); + if(len2 < 1) return(0); + 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); + double logProb(0.); + for(int i = iStart1; i < iEnd1; ++i) { + idx[1] = i; + for(int j = iStart2; j < iEnd2; ++j) { + idx[0] = j; + logProb += (static_cast(this)->*Method)(idx); + } + } + return(logProb); + } + // simulate + void simulate_cpp ( std::shared_ptr instr ) { + RESET_EIGEN_ERRORS; + int instr_type = 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_1_matp_ord_(instr); + if(instr_type == 5) return sim_2_seq_seq_(instr); + if(instr_type == 6) return sim_2_seq_seq_ord_(instr); + } + 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] + 1; + 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 dm = instr->dims[0]; + const auto& vals = instr->values->operator[](0); + if(len*dm != vals.size()) std::cout<<"len*dm != vals.size() in sim_1_matp_"< idx(dm); + for(int i = 0; i < len; ++i) { + for(int p = 0; p < dm; ++p) + idx[p] = vals[i*dm+p]; + static_cast(this)->sim_one(idx); + } + } + void sim_1_matp_ord_(std::shared_ptr instr) { + int len = instr->lens[0]; + int dm = instr->dims[0]; + const auto& vals = instr->values->operator[](0); + if(len*dm != vals.size()) std::cout<<"len*dm != vals.size() in sim_1_matp_"< idx(dm); + Eigen::Tensor slots(dm); + for(int p = 0; p < dm; ++p) + slots[p] = instr->slots[p]-1; + for(int i = 0; i < len; ++i) { + for(int p = 0; p < dm; ++p) + idx[slots[p]] = vals[i*dm+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_seq_ord_(std::shared_ptr instr) { + if(instr->slots[0] != 2 || instr->slots[1] != 1) + std::cout<<"slots not equal to 2,1 in calc_2_seq_seq_ord_"<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[1] = i; + for(int j = iStart2; j < iEnd2; ++j) { + idx[0] = j; + static_cast(this)->sim_one(idx); + } + } + } + virtual ~declFunClass_() {}; +}; diff --git a/nimbleModel/inst/include/nimbleModel/predef/instr_nClass/instr_nClass_copyFiles.txt b/nimbleModel/inst/include/nimbleModel/predef/instr_nClass/instr_nClass_copyFiles.txt new file mode 100644 index 0000000..e69de29 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 new file mode 100644 index 0000000..656aa1c --- /dev/null +++ b/nimbleModel/inst/include/nimbleModel/predef/instr_nClass/instr_nClass_cppContent.cpp @@ -0,0 +1,54 @@ +/* OPENER (Do not edit this comment) */ +#ifndef __instr_nClass_CPP +#define __instr_nClass_CPP +/* BODY (Do not edit this comment) */ +#ifndef R_NO_REMAP +#define R_NO_REMAP +#endif +#include +#include "instr_nClass_c_.h" +using namespace Rcpp; +// [[Rcpp::plugins(nCompiler_Eigen_plugin)]] +// [[Rcpp::depends(RcppParallel)]] +// [[Rcpp::depends(nCompiler)]] +// [[Rcpp::depends(Rcereal)]] + + instr_nClass::instr_nClass ( ) { +RESET_EIGEN_ERRORS +values = nClass_builder()(); +} + +// [[Rcpp::export(name = "instr_nClass_new")]] + SEXP new_instr_nClass ( ) { +RESET_EIGEN_ERRORS +return CREATE_NEW_NCOMP_OBJECT(instr_nClass);; +} + +// [[Rcpp::export(name = "set_CnClass_env_instr_nClass_new")]] + void set_CnClass_env_instr_nClass ( SEXP env ) { +RESET_EIGEN_ERRORS +SET_CNCLASS_ENV(instr_nClass, env);; +} + +// [[Rcpp::export(name = "get_CnClass_env_instr_nClass_new")]] + Rcpp::Environment get_CnClass_env_instr_nClass ( ) { +RESET_EIGEN_ERRORS +return GET_CNCLASS_ENV(instr_nClass);; +} + +NCOMPILER_INTERFACE( +instr_nClass, +NCOMPILER_FIELDS( +field("lens", &instr_nClass::lens), +field("index_types", &instr_nClass::index_types), +field("nDim", &instr_nClass::nDim), +field("dims", &instr_nClass::dims), +field("slots", &instr_nClass::slots), +field("values", &instr_nClass::values), +field("type", &instr_nClass::type), +field("sortID", &instr_nClass::sortID), +field("declID", &instr_nClass::declID) +), +NCOMPILER_METHODS() +) +#endif diff --git a/nimbleModel/inst/include/nimbleModel/predef/instr_nClass/instr_nClass_filebase.txt b/nimbleModel/inst/include/nimbleModel/predef/instr_nClass/instr_nClass_filebase.txt new file mode 100644 index 0000000..ba6031f --- /dev/null +++ b/nimbleModel/inst/include/nimbleModel/predef/instr_nClass/instr_nClass_filebase.txt @@ -0,0 +1 @@ +instr_nClass_c_ 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 new file mode 100644 index 0000000..f964e2c --- /dev/null +++ b/nimbleModel/inst/include/nimbleModel/predef/instr_nClass/instr_nClass_hContent.h @@ -0,0 +1,32 @@ +/* OPENER (Do not edit this comment) */ +#ifndef __instr_nClass_H +#define __instr_nClass_H +/* BODY (Do not edit this comment) */ +#ifndef R_NO_REMAP +#define R_NO_REMAP +#endif +#include +#include "nList_I1_c_.h" + +class instr_nClass : public interface_resolver< genericInterfaceC >, public loadedObjectHookC { +public: + instr_nClass ( ) ; + Eigen::Tensor lens; + Eigen::Tensor index_types; + int nDim; + Eigen::Tensor dims; + Eigen::Tensor slots; + std::shared_ptr values; + int type; + Eigen::Tensor sortID; + int declID; +}; + + SEXP new_instr_nClass ( ) ; + + void set_CnClass_env_instr_nClass ( SEXP env ) ; + + Rcpp::Environment get_CnClass_env_instr_nClass ( ) ; + + +#endif diff --git a/nimbleModel/inst/include/nimbleModel/predef/instr_nClass/instr_nClass_manifest.txt b/nimbleModel/inst/include/nimbleModel/predef/instr_nClass/instr_nClass_manifest.txt new file mode 100644 index 0000000..a3aa3c2 --- /dev/null +++ b/nimbleModel/inst/include/nimbleModel/predef/instr_nClass/instr_nClass_manifest.txt @@ -0,0 +1,7 @@ +list(saved_at = structure(1780285129.29081, class = c("POSIXct", +"POSIXt")), packet_name = "instr_nClass", elements = c("preamble", +"cppContent", "hContent", "filebase", "post_cpp_compiler", "copyFiles" +), files = list(preamble = "instr_nClass_preamble.cpp", cppContent = "instr_nClass_cppContent.cpp", + hContent = "instr_nClass_hContent.h", filebase = "instr_nClass_filebase.txt", + post_cpp_compiler = "instr_nClass_post_cpp_compiler.txt", + copyFiles = "instr_nClass_copyFiles.txt")) diff --git a/nimbleModel/inst/include/nimbleModel/predef/instr_nClass/instr_nClass_post_cpp_compiler.txt b/nimbleModel/inst/include/nimbleModel/predef/instr_nClass/instr_nClass_post_cpp_compiler.txt new file mode 100644 index 0000000..e69de29 diff --git a/nimbleModel/inst/include/nimbleModel/predef/instr_nClass/instr_nClass_preamble.cpp b/nimbleModel/inst/include/nimbleModel/predef/instr_nClass/instr_nClass_preamble.cpp new file mode 100644 index 0000000..1a92ce3 --- /dev/null +++ b/nimbleModel/inst/include/nimbleModel/predef/instr_nClass/instr_nClass_preamble.cpp @@ -0,0 +1,5 @@ +#define NCOMPILER_HANDLE_EIGEN_ERRORS +#define NCOMPILER_USES_EIGEN +#define NCOMPILER_USES_NCPPVEC +#define USES_NCOMPILER +#define NCOMPILER_USES_NCLASS_INTERFACE diff --git a/nimbleModel/inst/include/nimbleModel/predef/modelBase_nC/modelBase_nClass_copyFiles.txt b/nimbleModel/inst/include/nimbleModel/predef/modelBase_nC/modelBase_nClass_copyFiles.txt new file mode 100644 index 0000000..e69de29 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 new file mode 100644 index 0000000..d545f78 --- /dev/null +++ b/nimbleModel/inst/include/nimbleModel/predef/modelBase_nC/modelBase_nClass_cppContent.cpp @@ -0,0 +1,78 @@ +/* OPENER (Do not edit this comment) */ +#ifndef __modelBase_nClass_CPP +#define __modelBase_nClass_CPP +/* BODY (Do not edit this comment) */ +#ifndef R_NO_REMAP +#define R_NO_REMAP +#endif +#include +#include "modelBase_nClass_c_.h" +using namespace Rcpp; +// [[Rcpp::plugins(nCompiler_Eigen_plugin)]] +// [[Rcpp::depends(RcppParallel)]] +// [[Rcpp::depends(nCompiler)]] +// [[Rcpp::depends(Rcereal)]] +// [[Rcpp::depends(nimbleModel)]] + + bool modelBase_nClass::ping ( ) { +RESET_EIGEN_ERRORS +return(true); +} + std::shared_ptr modelBase_nClass::makeCompiledInstrList ( SEXP input ) { +RESET_EIGEN_ERRORS +std::shared_ptr ans; +ans = nClass_builder()(); +ans->set_all_values(input);; +return(ans); +} + 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 ) { +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 ) { +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 ) { +RESET_EIGEN_ERRORS +Rprintf("modelBase_nClass simulate_impl (should not see this)\n");; +} + modelBase_nClass::modelBase_nClass ( ) { +RESET_EIGEN_ERRORS +} + +// [[Rcpp::export(name = "set_CnClass_env_modelBase_nClass_new")]] + void set_CnClass_env_modelBase_nClass ( SEXP env ) { +RESET_EIGEN_ERRORS +SET_CNCLASS_ENV(modelBase_nClass, env);; +} + +// [[Rcpp::export(name = "get_CnClass_env_modelBase_nClass_new")]] + Rcpp::Environment get_CnClass_env_modelBase_nClass ( ) { +RESET_EIGEN_ERRORS +return GET_CNCLASS_ENV(modelBase_nClass);; +} + +NCOMPILER_INTERFACE( +modelBase_nClass, +NCOMPILER_FIELDS( +field("declFunList", &modelBase_nClass::declFunList), +field("declFunNameToIndex", &modelBase_nClass::declFunNameToIndex) +), +NCOMPILER_METHODS( +method("ping", &modelBase_nClass::ping, args({{}})), +method("makeCompiledInstrList", &modelBase_nClass::makeCompiledInstrList, args({{arg("input",copy)}})), +method("calculate_impl", &modelBase_nClass::calculate_impl, args({{arg("instrList",copy)}})), +method("calculateDiff_impl", &modelBase_nClass::calculateDiff_impl, args({{arg("instrList",copy)}})), +method("getLogProb_impl", &modelBase_nClass::getLogProb_impl, args({{arg("instrList",copy)}})), +method("simulate_impl", &modelBase_nClass::simulate_impl, args({{arg("instrList",copy)}})) +) +) +#endif diff --git a/nimbleModel/inst/include/nimbleModel/predef/modelBase_nC/modelBase_nClass_filebase.txt b/nimbleModel/inst/include/nimbleModel/predef/modelBase_nC/modelBase_nClass_filebase.txt new file mode 100644 index 0000000..e8994f8 --- /dev/null +++ b/nimbleModel/inst/include/nimbleModel/predef/modelBase_nC/modelBase_nClass_filebase.txt @@ -0,0 +1 @@ +modelBase_nClass_c_ 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 new file mode 100644 index 0000000..6422f2a --- /dev/null +++ b/nimbleModel/inst/include/nimbleModel/predef/modelBase_nC/modelBase_nClass_hContent.h @@ -0,0 +1,32 @@ +/* OPENER (Do not edit this comment) */ +#ifndef __modelBase_nClass_H +#define __modelBase_nClass_H +/* BODY (Do not edit this comment) */ +#ifndef R_NO_REMAP +#define R_NO_REMAP +#endif +#include +#include "declFunBase_nClass_c_.h" +#include "instr_nClass_c_.h" +#include "nList_instr_nClass_c_.h" + +class modelBase_nClass : public interface_resolver< genericInterfaceC >, public loadedObjectHookC { +public: + virtual bool ping ( ) ; + std::shared_ptr makeCompiledInstrList ( 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 ) ; + modelBase_nClass ( ) ; + double declFunList; + Rcpp::List declFunNameToIndex; +}; + + void set_CnClass_env_modelBase_nClass ( SEXP env ) ; + + Rcpp::Environment get_CnClass_env_modelBase_nClass ( ) ; + +#include + +#endif 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 new file mode 100644 index 0000000..01c4ca6 --- /dev/null +++ b/nimbleModel/inst/include/nimbleModel/predef/modelBase_nC/modelBase_nClass_manifest.txt @@ -0,0 +1,7 @@ +list(saved_at = structure(1780285516.41106, class = c("POSIXct", +"POSIXt")), packet_name = "modelBase_nClass", elements = c("preamble", +"cppContent", "hContent", "filebase", "post_cpp_compiler", "copyFiles" +), files = list(preamble = "modelBase_nClass_preamble.cpp", cppContent = "modelBase_nClass_cppContent.cpp", + hContent = "modelBase_nClass_hContent.h", filebase = "modelBase_nClass_filebase.txt", + post_cpp_compiler = "modelBase_nClass_post_cpp_compiler.txt", + copyFiles = "modelBase_nClass_copyFiles.txt")) diff --git a/nimbleModel/inst/include/nimbleModel/predef/modelBase_nC/modelBase_nClass_post_cpp_compiler.txt b/nimbleModel/inst/include/nimbleModel/predef/modelBase_nC/modelBase_nClass_post_cpp_compiler.txt new file mode 100644 index 0000000..e69de29 diff --git a/nimbleModel/inst/include/nimbleModel/predef/modelBase_nC/modelBase_nClass_preamble.cpp b/nimbleModel/inst/include/nimbleModel/predef/modelBase_nC/modelBase_nClass_preamble.cpp new file mode 100644 index 0000000..1a92ce3 --- /dev/null +++ b/nimbleModel/inst/include/nimbleModel/predef/modelBase_nC/modelBase_nClass_preamble.cpp @@ -0,0 +1,5 @@ +#define NCOMPILER_HANDLE_EIGEN_ERRORS +#define NCOMPILER_USES_EIGEN +#define NCOMPILER_USES_NCPPVEC +#define USES_NCOMPILER +#define NCOMPILER_USES_NCLASS_INTERFACE diff --git a/nimbleModel/inst/include/nimbleModel/predef/modelClass_/modelClass_.h b/nimbleModel/inst/include/nimbleModel/predef/modelClass_/modelClass_.h new file mode 100644 index 0000000..04db667 --- /dev/null +++ b/nimbleModel/inst/include/nimbleModel/predef/modelClass_/modelClass_.h @@ -0,0 +1,174 @@ +// to be included from the predefined modelBase_nClass. +// Add "#include " to that file, +// after the declaration of modelBase_nClass. + + +template +class modelClass_ : public modelBase_nClass { +public: + modelClass_() {}; + std::vector< std::shared_ptr > declFunPtrs; + std::map name2index_map; + double calculate_impl(std::shared_ptr instrList) override { + // double logProb(0.0); + // const auto& instrVec = instrList->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 { + return calc_op_impl< &declFunBase_nClass::calculateDiff_cpp >(instrList); + } + 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 logProb(0.0); + const auto& instrVec = instrList->contents(); + for (const auto& instr : instrVec) { + logProb += ((*(declFunPtrs[instr->declID - 1 ])).*Method)(instr); + } + return(logProb); + } + + void simulate_impl(std::shared_ptr instrList) { + const auto& instrVec = instrList->contents(); + for (const auto& instr : instrVec) { + declFunPtrs[instr->declID - 1 ]->simulate_cpp(instr); + } + } + + // 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. + // This may become rarely used because we will generate into a derived + // model class a canonical ordering + void do_setup_decl_mgmt_from_names(Rcpp::CharacterVector names) { + // Rprintf("Attempting setup_decl_mgmt_from_names with %d names\n", (int)names.length()); + Derived *self = static_cast(this); + const auto& name2access = self->get_name2access(); + declFunPtrs.clear(); + name2index_map.clear(); + size_t n = names.length(); + for(size_t i = 0; i < n; ++i) { + std::string name = Rcpp::as(names[i]); + auto it = name2access.find(name); + if(it != name2access.end()) { + std::shared_ptr ptr = it->second->getInterfacePtr(dynamic_cast(self)); + // When looking up this way, we do expect always to find objects (ptr valid) and that they are nodeFxn ptrs (ptr2 valid). + // So we can turn these messages into errors once things are working. + bool got_one = (ptr != nullptr); + if(got_one) { + // Rprintf("HOORAY: field %s is genericInterfaceBaseC\n", name.c_str()); + std::shared_ptr ptr2 = std::dynamic_pointer_cast(ptr); + bool step_two = (ptr2 != nullptr); + if(step_two) { + // Rprintf("AND IT IS A NODEFXN PTR!\n"); + name2index_map.emplace(name, declFunPtrs.size()); + declFunPtrs.push_back(ptr2); + } else { + // Rprintf("but it is not a nodefxn ptr\n"); + } + } else { + // Rprintf("field %s is NOT a genericInterfaceBaseC\n", name.c_str()); + } + } + } + } + + // This version scans all members to find nodeFxns. + // The resulting ordering comes from the order of the name2access map, + // and so may not match R. This was written first but may fall out of common use. + // This may become rarely used because we will generate into a derived + // model class a canonical ordering + void setup_auto_decl_mgmt() { + Derived *self = static_cast(this); + const auto& name2access = self->get_name2access(); + size_t n = name2access.size(); + //Rprintf("There are %d member variables indexed:\n", (int)n); + auto i_n2a = name2access.begin(); + auto end_n2a = name2access.end(); + declFunPtrs.clear(); + name2index_map.clear(); + size_t index = 0; + for(; i_n2a != end_n2a; ++i_n2a) { + std::shared_ptr ptr = i_n2a->second->getInterfacePtr(dynamic_cast(self)); + bool got_one = (ptr != nullptr); + if(got_one) { + // Rprintf("HOORAY: field %s is genericInterfaceBaseC\n", i_n2a->first.c_str()); + std::shared_ptr ptr2 = std::dynamic_pointer_cast(ptr); + bool step_two = (ptr2 != nullptr); + if(step_two) { + // Rprintf("AND IT IS A NODEFXN PTR!\n"); + declFunPtrs.push_back(ptr2); + name2index_map.emplace(i_n2a->first, index++); + } else { + // Rprintf("but it is not a nodefxn ptr\n"); + } + } + else { + // Rprintf("field %s is NOT a genericInterfaceBaseC\n", i_n2a->first.c_str()); + } + } + } + void c_print_nodes() { + auto i_n2i = name2index_map.begin(); + auto end_n2i = name2index_map.end(); + Rprintf("0-based index: name\n"); + for(; i_n2i != end_n2i; ++i_n2i) { + Rprintf("%d: %s\n", i_n2i->first.c_str(), (int)i_n2i->second); + } + } + void set_from_list(Rcpp::List Rlist) { + Rcpp::CharacterVector Rnames = Rlist.names(); + size_t len = Rnames.length(); + for(size_t i = 0; i < len; ++i) { + // explicit cast is needed because even though Rnames[i] can cast to a string, + // set_value takes a const string& so we need an object in place here. + // set_value fails safely if a name is not found. + static_cast(this)->set_value(std::string(Rnames[i]), Rlist[i]); + } + } + void resize_from_list(Rcpp::List Rlist) { + Rcpp::CharacterVector Rnames = Rlist.names(); + size_t len = Rnames.length(); + size_t vec_len; + Rcpp::IntegerVector vs; + for(size_t i = 0; i < len; ++i) { + // explicit cast is needed because even though Rnames[i] can cast to a string, + // set_value takes a const string& so we need an object in place here. + vs = Rlist[i]; + vec_len = vs.length(); + std::unique_ptr ETA = static_cast(this)->access(std::string(Rnames[i])); + // if the name was not found, a "Problem:" message was emitted, and we skip using it here. + if(ETA) { + switch(vec_len) { + case 0 : + break; + case 1 : + ETA->template ref<1>().resize(vs[0]); + break; + case 2 : + ETA->template ref<2>().resize(vs[0], vs[1]); + break; + case 3 : + ETA->template ref<3>().resize(vs[0], vs[1], vs[2]); + break; + case 4 : + ETA->template ref<4>().resize(vs[0], vs[1], vs[2], vs[3]); + break; + case 5 : + ETA->template ref<5>().resize(vs[0], vs[1], vs[2], vs[3], vs[4]); + break; + case 6 : + ETA->template ref<6>().resize(vs[0], vs[1], vs[2], vs[3], vs[4], vs[5]); + break; + } + } + } + } +}; diff --git a/nimbleModel/tests/testthat/test-nimbleModel.R b/nimbleModel/tests/testthat/test-nimbleModel.R new file mode 100644 index 0000000..901abe8 --- /dev/null +++ b/nimbleModel/tests/testthat/test-nimbleModel.R @@ -0,0 +1,564 @@ +# Test code needed for new nimbleModel system. + +library(nCompiler) +library(nimbleModel) +library(testthat) + +## TODO: will location and access to predefined nClasses be as described below given they will live +## in `nimbleModel` package? How will dependence on nCompiler work? + +## # To update the set of predefined nClasses +## # generate new predef/instr_nC. Move that directly to package code inst/nimbleModel/predef/instr_nC +## nCompile(instr_nClass = nimbleModel:::instr_nClass, control=list(generate_predefined=TRUE)) +## test <- nCompile(instr_nClass = nimbleModel:::instr_nClass) +## # +## # generate new predef/declFunBase_nC. Move to package and add +## # "#include " in the hContent +## # And add "// [[Rcpp::depends(nimbleModel)]]" to the cppContent +## # after declaration of declFunBase_nClass +## nCompile(nimbleModel:::declFunBase_nClass, control=list(generate_predefined=TRUE)) +## test <- nCompile(nimbleModel:::declFunBase_nClass) +## # +## # generate new predef/modelBase_nC. Move to package and add +## # "#include " to the hContent +## # And add "// [[Rcpp::depends(nimbleModel)]]" to the cppContent +## # after the declaration of modelBase_nClass. +## nCompile(modelBase_nClass = nimbleModel:::modelBase_nClass, control=list(generate_predefined=TRUE)) +## test <- nCompile(nimbleModel:::modelBase_nClass) +## #nCompile(instr_nClass, modelBase_nClass, declFunBase_nClass, control=list(generate_predefined=TRUE)) + +## TODO: revise these tests for instrClass (flattened approach) + +test_that("initial test of compiled model", { + code <- quote({ + tau ~ dunif(0, 100) + mu ~ dnorm(0,1) + for(i in 1:5) { + y[i] ~ dnorm(mu, var = tau) + } + }) + + inits <- list(tau = 25, mu = 0) + data <- list(y = rnorm(5)) + + ## "Manual" workflow not using `nimbleModel()`. + nm <- modelClass$new(code, inits = inits, data = data) + mclass <- nimbleModel:::make_modelClass_from_nimbleModel(nm) + + Cmclass <- nCompile(mclass) + Cobj <- Cmclass$new() + obj <- mclass$new() + + # Check a first calculation on a simple node + Cans <- Cobj$calculate('tau') + ans <- obj$calculate('tau') + check <- dunif(Cobj$tau, 0, 100, log = TRUE) + expect_equal(Cans, ans) + expect_equal(Cans, check) + + # Check entire model, also getting lifted sd node computed + Cans <- Cobj$calculate() + ans <- obj$calculate() + expect_equal(Cans, ans) + + # Check a sequence + Cans <- Cobj$calculate('y[1:3]') + ans <- obj$calculate('y[1:3]') + check <- dnorm(Cobj$y[1:3], Cobj$mu, sqrt(Cobj$tau), log=TRUE) |> sum() + expect_equal(Cans, ans) + expect_equal(Cans, check) + + # Check a non-contiguous pair of nodes (a mat case) + nodes <- c('y[2]','y[4]') + Cans <- Cobj$calculate(nodes) + ans <- obj$calculate(nodes) + check <- dnorm(Cobj$y[c(2, 4)], Cobj$mu, sqrt(Cobj$tau), log=TRUE) |> sum() + expect_equal(Cans, ans) + expect_equal(Cans, check) + + # Check getLogProb + Cans <- Cobj$getLogProb('y[1:4]') + ans <- obj$calculate('y[1:4]') + check <- dnorm(Cobj$y[1:4], Cobj$mu, sqrt(Cobj$tau), log=TRUE) |> sum() + expect_equal(Cans, ans) + expect_equal(Cans, check) + + # Prepare for calculateDiff test below + old_logProb <- dnorm(Cobj$y[3:4], Cobj$mu, sqrt(Cobj$tau), log=TRUE) |> sum() + + # Check simulate + set.seed(1) + Cobj$simulate('y[3:4]') + set.seed(1) + obj$simulate('y[3:4]') + expect_equal(Cobj$y, obj$y) + + # Check getLogProb + # Do this assignment in case the previous test of repeatability fails + obj$y[3:4] <- Cobj$y[3:4] + Cans <- Cobj$calculateDiff('y[3:4]') + ans <- obj$calculateDiff('y[3:4]') + new_logProb <- dnorm(Cobj$y[3:4], Cobj$mu, sqrt(Cobj$tau), log=TRUE) |> sum() + check <- new_logProb - old_logProb + expect_equal(Cans, ans) + expect_equal(Cans, check) + + # Always end compiled tests with removing and garbage collecting + # to ensure gc() happens while the DLL is still in place. + rm(Cobj, obj); gc() +}) + +test_that("initial tests/examples of nimble model using flattened approach", { + + code <- quote({ + tau ~ dunif(0, 100) + mu ~ dnorm(0,1) + for(i in 1:5) { + y[i] ~ dnorm(mu, var = tau) + } + }) + + inits <- list(tau = 25, mu = 0) + data <- list(y = rnorm(5)) + + ## "Manual" workflow not using `nimbleModel()`. + nm <- modelClass$new(code, inits = inits, data = data) + #debug(nimbleModel:::make_modelClass_from_nimbleModel) + #debug(nimbleModel:::makeModel_nClass) + mclass <- nimbleModel:::make_modelClass_from_nimbleModel(nm) + + # Begin Perry + Cmclass <- nCompile(mclass) + Cobj <- Cmclass$new() + #Cobj$calculate_impl + #Cobj$calculate + #debug(Cobj$calculate) + Cobj$calculate('tau') + Cobj$calculate() + Cobj$calculate('y[1]') + dnorm(Cobj$y[1], Cobj$mu, sqrt(Cobj$tau), log=TRUE) + Cobj$calculate('y[1:3]') + dnorm(Cobj$y[1:3], Cobj$mu, sqrt(Cobj$tau), log=TRUE) |> sum() + + + NULL + + obj <- mclass$new() + obj$calculate() + #debug(obj$calculate) + obj$calculate('y[1]') + obj$calculate('y[1:3]') + NULL + # PROBLEM, in nList_<>::set_from_list for uncompiled list input. + # I guess set_all_values should skip NULLs? Or maybe only for non-R targets? + # Give a better message than "Bad type". Pass the name? Check for NULL? + + # Next steps + # initialize the instrs from uncompiled if needed + # + + # check technique of building and copying nList(instr_nClass) as a method: + inC <- nimbleModel:::instr_nClass + test1 <- nFunction( + function(Robj = 'SEXP') { + ans <- inC$new() + cppLiteral("ans->set_all_values(Robj);") + cppLiteral("std::cout<dim<set_all_values(Robj);") +# cppLiteral("std::cout<dim< 95)) + + ## Use of nimbleModel + mclass <- nimbleModel(code, data = data, inits = inits) + m <- mclass$new() + expect_identical(m$calculate('tau'), dunif(m$tau, 0, 100, log = TRUE)) + + m <- nimbleModel(code, data = data, inits = inits, returnClass = FALSE) + expect_identical(m$calculate('tau'), dunif(m$tau, 0, 100, log = TRUE)) + + ## Override init value when creating model instance. + mclass <- nimbleModel(code, data = data, inits = inits) + m <- mclass$new(inits = list(tau = 7)) + expect_identical(m$tau, 7) + +}) + +test_that("multiple index slots, single indexRange case", { + code <- quote({ + for(i in 1:5) + for(j in 1:3) + y[i,j] ~ dnorm(0,1) + }) + data <- list(y = matrix(rnorm(15),5)) + mclass <- nimbleModel(code, data = data) + m <- mclass$new() + vr <- varRangeClass$new(list(newIndexRange(matrix(c(2,4,3,1), ncol=2))), varName='y') + expect_equal(m$calculate(vr), dnorm(data$y[2,3],log=TRUE) + dnorm(data$y[4,1],log=TRUE)) + cmclass <- nCompile(mclass) + cm <- cmclass$new() + expect_equal(cm$calculate(vr), dnorm(data$y[2,3],log=TRUE) + dnorm(data$y[4,1],log=TRUE)) + + set.seed(1) + m$simulate(vr) + set.seed(1) + cm$simulate(vr) + expect_equal(m$y, cm$y) + + ## Now with slot reordering. + vr <- varRangeClass$new(list(newIndexRange(matrix(c(2,4,3,1), ncol=2))), + rangeToIndexSlot = list(2:1), varName='y') + expect_equal(m$calculate(vr), dnorm(data$y[3,2],log=TRUE) + dnorm(data$y[1,4],log=TRUE)) + expect_equal(cm$calculate(vr), dnorm(data$y[3,2],log=TRUE) + dnorm(data$y[1,4],log=TRUE)) + + ## 3-d case for more robust ordering check + code <- quote({ + for(i in 1:5) + for(j in 1:4) + for(k in 1:3) + y[i,j,k] ~ dnorm(0,1) + }) + data <- list(y = array(rnorm(60),c(5,4,3))) + mclass <- nimbleModel(code, data = data) + m <- mclass$new() + vr <- varRangeClass$new(list(newIndexRange(matrix(c(2,3,1,3,5,2,1,2,4), ncol=3))), + rangeToIndexSlot = list(c(3,1,2)), varName='y') + truth <- dnorm(data$y[3,1,2],log=TRUE) + dnorm(data$y[5,2,3],log=TRUE) + dnorm(data$y[2,4,1],log=TRUE) + expect_equal(m$calculate(vr), truth) + cmclass <- nCompile(mclass) + cm <- cmclass$new() + expect_equal(cm$calculate(vr), truth) + + +}) + +test_that("two sequences case", { + code <- quote({ + for(i in 1:5) + for(j in 1:3) + y[i,j] ~ dnorm(0,1) + }) + data <- list(y = matrix(rnorm(15),5)) + mclass <- nimbleModel(code, data = data) + truth <- sum(dnorm(data$y[2:4,1:3],0,1,log=TRUE)) + m <- mclass$new() + expect_equal(m$calculate('y[2:4,1:3]'), truth) + cmclass <- nCompile(mclass) + cm <- cmclass$new() + expect_equal(cm$calculate('y[2:4,1:3]'), truth) + + set.seed(1) + m$simulate('y[2:4,1:3]') + set.seed(1) + cm$simulate('y[2:4,1:3]') + expect_equal(m$y, cm$y) + + ## 2-d case for ordering check. + ## This does not test calc_2_seq_seq_ord because rule application in creating instr + ## already re-sorts the indexRanges. + code <- quote({ + for(i in 1:5) + for(j in 1:2) + y[i,j] ~ dnorm(0,1) + }) + data <- list(y = matrix(rnorm(10),5)) + mclass <- nimbleModel(code, data = data) + m <- mclass$new() + vr <- varRangeClass$new(list(newIndexRange(quote(1:2)), newIndexRange(quote(1:5))), + rangeToIndexSlot = list(2,1), varName='y') + truth <- sum(dnorm(data$y, log = TRUE)) + expect_equal(m$calculate(vr), truth) + cmclass <- nCompile(mclass) + cm <- cmclass$new() + expect_equal(cm$calculate(vr), truth) + +}) + +test_that("basic creation of list of instr_nClass objects", { + + code <- quote({ + for(i in 1:5) { + mu ~ dnorm(0, 1) + y[i] ~ dnorm(mu, 1) + } + }) + + data <- list(y = rnorm(5)) + + m <- nimbleModel(code, data = data, returnClass = FALSE) + + instr0 <- makeInstrList(m, 'mu')[[1]] + expect_identical(instr0$lens, 1) + expect_identical(length(instr0$values), 0L) + expect_identical(instr0$index_types, 0) + expect_identical(instr0$type, 0) + + instr1 <- makeInstrList(m, 'y[3:4]')[[1]] + expect_identical(instr1$lens, 2) + expect_identical(instr1$values[[1]], 2) # offset + expect_identical(instr1$index_types, 1) + expect_identical(instr1$type, 1) + + instr2 <- makeInstrList(m, c('y[c(2,5)]'))[[1]] + expect_identical(instr2$lens, 2) + expect_identical(instr2$values[[1]], c(2,5)) + expect_identical(instr2$index_types, 2) + expect_identical(instr2$type, 2) + + instr2 <- makeInstrList(m, varRangeClass$new(list(newIndexRange(matrix(c(2,5), ncol=1))), varName='y'))[[1]] + expect_identical(instr2$lens, 2) + expect_identical(instr2$values[[1]], c(2,5)) + expect_identical(instr2$index_types, 2) + expect_identical(instr2$type, 2) + + ## TODO: flesh this out with multiple index cases. +}) + +## TODO: modify tests below in light of flattened approach. + +test_that("nimble model prototype works", { + nodeVarInfo <- list(list(name = "x", nDim = 1), list(name = "mu", nDim = 1), + list(name = "sd", nDim = 0)) + calc_one <- nFunction( + name = "calc_one", + fun = function(inds) { + ans <- model$x[inds[1]] + return(ans) + }, + compileInfo = list( + C_fun = function(inds = 'integerVector') { + returnType('numericScalar') + ans <- x[inds[1]] + return(ans) + } + ) + ) + my_nodeFxn <- make_node_nClass(nodeVarInfo, list(calc_one=calc_one), "test_node") + my_nodeInfo <- make_node_info_for_model_nClass("beta_NF1", "my_nodeFxn", "test_node", nodeVarInfo) + + modelVarInfo <- list(list(name="x", nDim = 1), + list(name = "mu", nDim = 1), + list(name = "sd", nDim = 0), + list(name = "gamma", nDim = 2)) + #debug(makeModel_nClass) + ncm1 <- makeModel_nClass(modelVarInfo, list(my_nodeInfo), classname = "my_model", env=environment()) + #undebug(addGenericInterface_impl) + #undebug(nCompile_finish_nonpackage) + for(package in c(FALSE, TRUE)) { + Cncm1 <- nCompile(ncm1, returnList=TRUE, package=package) + #Cncm1 <- nCompile(modelBase_nClass, nodeFxnBase_nClass, calcInstrList_nClass, calcInstr_nClass, nodeInstr_nClass, ncm1, my_nodeFxn) + for(mode in c("uncompiled", "compiled")) { + if(mode=="compiled") { + obj <- Cncm1$ncm1$new() + } else { + obj <- ncm1$new() + } + # obj$do_setup_node_mgmt() + nodeObj <- obj$beta_NF1 + obj$x <- 1:3 + expect_equal(obj$x, 1:3) + + obj$set_from_list(list(x = 10:11)) + # expect Problem msg: (alpha is not a field in the class) + obj$set_from_list(list(mu = 110, x = 11:20, alpha = 101)) + obj$mu + + obj$resize_from_list(list(x = 7)) + # expect Problem msg: + obj$resize_from_list(list(alpha = 5, mu = 3, gamma = c(2, 4))) + expect_equal(length(obj$mu), 3) + expect_equal(dim(obj$gamma), c(2, 4)) + obj$resize_from_list(list(x = 5, gamma = c(3, 5))) + expect_equal(length(obj$x), 5) + expect_equal(dim(obj$gamma), c(3, 5)) + + obj$x <- 11:15 + expect_equal(nodeObj$calc_one(c(3)), 13) + rm(obj, nodeObj); gc() + } + } +}) + +test_that("nodeInstr_nClass and calcInstr_nClass basics work", { + for(package in c(FALSE, TRUE)) { + test <- nCompile(nodeInstr_nClass, calcInstr_nClass, calcInstrList_nClass, control=list(generate_predefined=FALSE), package = package) + calcInstrList <- test$nList_calcInstr_nClass$new() + calcInstr <- test$calcInstr_nClass$new() + expect_equal(calcInstr$nodeInstrVec, NULL) + ni1 <- test$nodeInstr_nClass$new() + ni2 <- test$nodeInstr_nClass$new() + ni1$methodInstr <- 1 + ni2$methodInstr <- 2 +# nList("integerVector")$new() +# ni1$indsInstrVec <- nList("integerVector")$new() + ni1$indsInstrVec[1:2] <- list(1:2, 3:4) + ni2$indsInstrVec + ni2$indsInstrVec[1:2] <- list(11:12, 13:14) + calcInstr$nodeInstrVec + calcInstr$nodeInstrVec[1:2] <- list(ni1, ni2) + + expect_true(length(calcInstr$nodeInstrVec)==2) + expect_identical(calcInstr$nodeInstrVec[[1]]$indsInstrVec |> as.list(), list(1:2, 3:4)) + expect_identical(calcInstr$nodeInstrVec[[2]]$indsInstrVec |> as.list(), list(11:12, 13:14)) + calcInstrList[1] <- list(calcInstr) + expect_equal(calcInstrList |> as.list(), list(calcInstr)) + rm(calcInstrList, calcInstr, ni1, ni2); gc() + } +}) + +###### + +## This is somewhat redundant with the first test +test_that("nimble model variables are set up", { + library(nimbleModel) + code <- quote({ + sd ~ dunif(0, 10) + for(i in 1:5) { + y[i] ~ dnorm(x[i+1], sd = sd) + } + }) + m <- modelClass$new(code) + varInfo <- get_varInfo_from_nimbleModel(m) + modelVars <- varInfo$vars + # Try making a model with no nodeFxns + ncm1 <- makeModel_nClass(modelVars, list(), classname = "my_model", env = environment()) + Cncm1 <- nCompile(ncm1, returnList=TRUE) + #Cncm1 <- nCompile(modelBase_nClass, nodeFxnBase_nClass, calcInstrList_nClass, calcInstr_nClass, nodeInstr_nClass, ncm1) + obj <- Cncm1$ncm1$new() + obj$resize_from_list(varInfo$sizes) + expect_equal(length(obj$x), 6) + expect_equal(length(obj$y), 5) + expect_equal(length(obj$logProb_y), 5) +}) + +######## +# nOptions(pause_after_writing_files=TRUE) +# Try automating the whole model creation including nodeFxns +# Ditto: this works but relies on nimbleModel +test_that("nimble model with stochastic and deterministic nodes is created and compiles", { + library(nimbleModel) + code <- quote({ + sd ~ dunif(0, 10) + for(i in 1:5) { + z[i] <- x[i+1] + 10 + y[i] ~ dnorm(x[i+1], sd = sd) + } + }) + m <- modelClass$new(code) + + ## Check that a separate R implementation was created + mDef_ <- m$modelDef + dI <- mDef_$declInfo[[2]] + nFxn <- make_node_methods_from_declInfo(dI) + expect_true(!is.null(NFinternals(nFxn[[1]])$R_fun)) + dI <- mDef_$declInfo[[3]] + nFxn <- make_node_methods_from_declInfo(dI) + expect_true(!is.null(NFinternals(nFxn[[1]])$R_fun)) + + for(mode in c("uncompiled", "compiled")) { + package_options <- if(mode=="compiled") c(FALSE, TRUE) else TRUE + for(package in package_options) { + nMod <- make_model_from_nimbleModel(m, compile=FALSE) + if(mode=="compiled") { + expect_no_error(CnMod <- nCompile(nMod, package = package)) + nMod <- CnMod + } + expect_no_error(obj <- nMod$new()) + obj$y <- 1:5 + expect_equal(obj$y, 1:5) + vals <- list(x = 2:7, y = 11:15, sd = 8) + obj$set_from_list(vals) + expect_equal(obj$x, vals$x) + rm(obj); gc() + } + } +}) + +message("test-nimbleModel does not have tests of calculate etc.") + +if(FALSE) { + nodeFxn_2_nodeIndex <- c(nodeFxn_1 = 1, nodeFxn_3 = 2) + + calcInputList <- list(list(nodeFxn="nodeFxn_1", # which declaration (nodeFxn) + nodeInputVec = list(list(methodInput=1, # which index iteration method + indsInputVec=list(1))))) # input(s) to index iterations + + calcInstrList <- calcInputList_to_calcInstrList(calcInputList, test) + + obj$calculate(calcInstrList) +}