Chapter 15 Writing nimbleFunctions to interact with models
15.1 Overview
When you write an R function, you say what the input arguments are,
you provide the code for execution, and in that code you give the
value to be returned29. Using the
function keyword in R triggers the operation of
creating an object that is the function.
Creating nimbleFunctions is similar, but there are two kinds of code and two steps of execution:
Setup code is provided as a regular R function, but the programmer does not control what it returns. Typically the inputs to setup code are objects like a model, a vector of nodes, a modelValues object or a modelValues configuration, or another nimbleFunction. The setup code, as its name implies, sets up information for run-time code. It is executed in R, so it can use any aspect of R.
Run code is provided in the NIMBLE language, which was introduced in Chapter 11. This is similar to a narrow subset of R, but it is important to remember that it is different – defined by what can be compiled – and much more limited. Run code can use the objects created by the setup code. In addition, some information on variable types must be provided for input arguments, the return value, and in some circumstances for local variables. There are two kinds of run code:
- There is always a primary function, given as the argument
run30. - There can optionally be other functions, or ‘methods’ in the language of object-oriented programming, that share the same objects created by the setup function.
- There is always a primary function, given as the argument
Here is a small example to fix ideas:
logProbCalcPlus <- nimbleFunction(
setup = function(model, node) {
dependentNodes <- model$getDependencies(node)
valueToAdd <- 1
},
run = function(P = double(0)) {
model[[node]] <<- P + valueToAdd
return(model$calculate(dependentNodes))
returnType(double(0))
})
code <- nimbleCode({
a ~ dnorm(0, 1)
b ~ dnorm(a, 1)
})
testModel <- nimbleModel(code, check = FALSE)
logProbCalcPlusA <- logProbCalcPlus(testModel, "a")
testModel$b <- 1.5
logProbCalcPlusA$run(0.25) ## [1] -2.650377## [1] -2.650377## [1] 1.25The call to the R function called nimbleFunction returns a
function, similarly to defining a function in R. That function,
logProbCalcPlus, takes arguments for its setup function,
executes it, and returns an object, logProbCalcPlusA, that has a
run member function (method) accessed by $run. In this case, the
setup function obtains the stochastic dependencies of the
node using the getDependencies member function of the model
(see Section 13.1.3) and stores them in
dependentNodes. In this way, logProbCalcPlus can adapt to any
model. It also creates a variable, valueToAdd, that can be used by the nimbleFunction.
The object logProbCalcPlusA, returned by logProbCalcPlus,
is permanently bound to the results of the processed setup
function. In this case, logProbCalcPlusA$run takes a scalar input value, P,
assigns P + valueToAdd to
the given node in the model, and returns the sum of the log
probabilities of that node and its stochastic
dependencies31. We say logProbCalcPlusA is an ‘instance’ of
logProbCalcPlus that is
‘specialized’ or ‘bound’ to a and testModel. Usually, the
setup code will be where information about the model
structure is determined, and then the run code can use that
information without repeatedly, redundantly recomputing it. A
nimbleFunction can
be called repeatedly (one can think of it as a generator), each time returning a specialized
nimbleFunction.
Readers familiar with object-oriented programming may find it useful
to think in terms of class definitions and objects. nimbleFunction
creates a class definition. Each specialized nimbleFunction is one object
in the class. The setup arguments are used to define member data in
the object.
15.2 Using and compiling nimbleFunctions
To compile the nimbleFunction, together with its model, we use compileNimble:
CnfDemo <- compileNimble(testModel, logProbCalcPlusA)
CtestModel <- CnfDemo$testModel
ClogProbCalcPlusA <- CnfDemo$logProbCalcPlusAThese have been initialized with the values from their uncompiled versions and can be used in the same way:
## [1] 1.25## [1] 1.5## [1] -5.462877# verify the answer:
dnorm(CtestModel$b, CtestModel$a, 1, log = TRUE) +
dnorm(CtestModel$a, 0, 1, log = TRUE) ## [1] -5.462877## [1] 2.5## [1] 1.2515.3 Writing setup code
15.3.1 Useful tools for setup functions
The setup function is typically used to determine information on nodes in a model, set up modelValues or nimbleList objects, set up (nested) nimbleFunctions or nimbleFunctionLists, and set up any persistent numeric objects. For example, the setup code of an MCMC nimbleFunction creates the nimbleFunctionList of sampler nimbleFunctions. The values of numeric objects created in setup code can be modified by run code and will persist across calls.
Some of the useful tools and objects to create in setup functions include:
- vectors of node names, often from a model Often these are obtained from the
getNodeNames,getDependencies, and other methods of a model, described in Sections 13.1-13.2. - modelValues objects These are discussed in Sections 14.1 and 15.4.4.
- nimbleList objects New instances of
nimbleLists can be created from a nimbleList definition in either setup or run code. See Section 14.2 for more information. - specializations of other nimbleFunctions A useful NIMBLE programming technique is to have one nimbleFunction contain other nimbleFunctions, which it can use in its run-time code (Section 15.4.7).
- lists of other nimbleFunctions In addition to containing single other nimbleFunctions, a nimbleFunction can contain a list of other nimbleFunctions (Section 15.4.8).
If one wants a nimbleFunction that does get specialized but has
empty setup code, use setup = function() {} or setup = TRUE.
15.3.2 Accessing and modifying numeric values from setup
While models and nodes created during setup cannot be modified32, numeric values and modelValues can be, as illustrated by extending the example from above.
## [1] 1## [1] 1## [1] -16.46288## [1] 4.515.3.3 Determining numeric types in nimbleFunctions
For numeric variables from the setup function that
appear in the run function or other member functions (or
are declared in setupOutputs), the
type is determined from the values created by the setup
code. The types created by setup code must be
consistent across all specializations of the nimbleFunction. For
example if X is created as a matrix (two-dimensional double) in one
specialization but as a vector (one-dimensional double) in another, there
will be a problem during compilation. The sizes may differ in each specialization.
Treatment of vectors of length one presents special challenges because
they could be treated as scalars or vectors. Currently they are
treated as scalars. If you want a vector, ensure that the length is
greater than one in the setup code and then use setSize in the
run-time code.
15.3.4 Control of setup outputs
Sometimes setup code may create variables that are not used in
run code. By default, NIMBLE inspects run code and omits
variables from setup that do not appear in run code from
compilation. However, sometimes a programmer may want to force a
numeric or character variable to be included in compilation, even if it
is not used directly in run code. As shown below, such variables
can be directly accessed in one nimbleFunction from another, which
provides a way of using nimbleFunctions as general data structures.
To force NIMBLE to include variables during compilation, for
example X and Y, simply include
anywhere in the setup code.
15.4 Writing run code
In Chapter 11 we described the functionality of the NIMBLE language that could be used in run code without setup code (typically in cases where no models or modelValues are needed). Next we explain the additional features that allow use of models and modelValues in the run code.
15.4.1 Driving models: calculate, calculateDiff, simulate, getLogProb
These four functions are the primary ways to operate a model. Their
syntax was explained in Section 13.3. Except
for getLogProb, it is usually important for the nodes
vector to be sorted in
topological order. Model member functions such as getDependencies and
expandNodeNames will
always return topoligically sorted node names.
Most R-like indexing of a node vector is allowed within the
argument to calculate, calculateDiff, simulate, and getLogProb. For example, all of the following are allowed:
myModel$calculate(nodes)
myModel$calculate(nodes[i])
myModel$calculate(nodes[1:3])
myModel$calculate(nodes[c(1,3)])
myModel$calculate(nodes[2:i])
myModel$calculate(nodes[ values(model, nodes) + 0.1 < x ])Note that in the last line of code, one must have that the length of nodes is
equal to that of values(model, nodes), which means that all the nodes in nodes
must be scalar nodes.
Also note that one cannot create new vectors of nodes in run code. They
can only be indexed within a call to calculate, calculateDiff,
simulate or getLogProb.
15.4.2 Getting and setting variable and node values
15.4.2.1 Using indexing with nodes
Here is an example that illustrates getting and setting of nodes, subsets of nodes, or variables. Note the following:
- In
model[[v]],vcan only be a single node or variable name, not a vector of multiple nodes nor an element of such a vector (model[[ nodes[i] ]]does not work). The node itself may be a vector, matrix or array node. - In fact,
vcan be a node-name-like character string, even if it is not actually a node in the model. See example 4 in the code below. - One can also use
model$varName, with the caveat thatvarNamemust be a variable name. This usage would only make sense for a nimbleFunction written for models known to have a specific variable. (Note that ifais a scalar node inmodel, thenmodel[['a']]will be a scalar butmodel$awill be a vector of length 1. - one should use the
<<-global assignment operator to assign into model nodes.
Note that NIMBLE does not allow variables to change dimensions. Model nodes are the same, and indeed are more restricted because they can’t change sizes. In addition, NIMBLE distinguishes between scalars and vectors of length 1. These rules, and ways to handle them correctly, are illustrated in the following code as well as in Section 11.3.
code <- nimbleCode({
z ~ dnorm(0, sd = sigma)
sigma ~ dunif(0, 10)
y[1:n] ~ dmnorm(zeroes[1:n], cov = C[1:5, 1:5])
})
n <- 5
m <- nimbleModel(code, constants = list(n = n, zeroes = rep(0, n),
C = diag(n)))
cm <- compileNimble(m)
nfGen <- nimbleFunction(
setup = function(model) {
# node1 and node2 would typically be setup arguments, so they could
# have different values for different models. We are assigning values
# here so the example is clearer.
node1 <- 'sigma' # a scalar node
node2 <- 'y[1:5]' # a vector node
notReallyANode <- 'y[2:4]' # y[2:4] allowed even though not a node!
},
run = function(vals = double(1)) {
tmp0 <- model[[node1]] # 1. tmp0 will be a scalar
tmp1 <- model[[node2]] # 2. tmp1 will be a vector
tmp2 <- model[[node2]][1] # 3. tmp2 will be a scalar
tmp3 <- model[[notReallyANode]] # 4. tmp3 will be a vector
tmp4 <- model$y[3:4] # 5. hard-coded access to a model variable
# 6. node1 is scalar so can be assigned a scalar:
model[[node1]] <<- runif(1)
model[[node2]][1] <<- runif(1)
# 7. an element of node2 can be assigned a scalar
model[[node2]] <<- runif(length(model[[node2]]))
# 8. a vector can be assigned to the vector node2
model[[node2]][1:3] <<- vals[1:3]
# elements of node2 can be indexed as needed
returnType(double(1))
out <- model[[node2]] # we can return a vector
return(out)
}
)
Rnf <- nfGen(m)
Cnf <- compileNimble(Rnf)
Cnf$run(rnorm(10))## [1] -1.4887622 -1.1607519 1.4577382 0.4855725 0.6932443Use of [[ ]] allows one to programmatically access a node based on a character variable containing the node name; this character variable would generally be set in setup code. In contrast, use of $ hard codes the variable name and would not generally be suitable for nimbleFunctions intended for use with arbitrary models.
15.4.2.2 Getting and setting more than one model node or variable at a time using values
Sometimes it is useful to set a collection of nodes or variables at
one time. For example, one might want a nimbleFunction that will
serve as the objective function for an optimizer. The input to the
nimbleFunction would be a vector, which should be used to fill a
collection of nodes in the model before calculating their log
probabilities. This can be done using values:
# get values from a set of model nodes into a vector
P <- values(model, nodes)
# or put values from a vector into a set of model nodes
values(model, nodes) <- Pwhere the first line would assign the collection of values from nodes
into P, and the second would do the inverse. In both cases, values
from nodes with two or more dimensions are flattened into a vector in
column-wise order.
values(model, nodes) may be used as a
vector in other expressions, e.g.,
One can also index elements of nodes in the argument to values, in the same manner as discussed for calculate and related functions in Section 15.4.1.
Note again the potential for confusion between scalars and vectors of
length 1. values returns a vector and expects a vector when used
on the left-hand side of an assignment. If only a single value is
being assigned, it must be a vector of length 1, not a scalar. This
can be achieved by wrapping a scalar in c() when necessary.
For example:
15.4.3 Getting parameter values and node bounds
Sections 13.2.3-13.2.4 describe how to get the parameter values for a node and the range (bounds) of possible values for the node using getParam and getBound. Both of these can be used in run code.
15.4.4 Using modelValues objects
The modelValues structure was introduced in Section
14.1. Inside nimbleFunctions, modelValues are
designed to easily save values from a model object during the running
of a nimbleFunction. A modelValues object used in run code
must always exist in the setup code, either by passing it in as a
setup argument or creating it in the setup code.
To illustrate this, we will create a nimbleFunction for computing
importance weights for importance sampling. This function will use two
modelValues objects. propModelValues will contain a set of
values simulated from the importance sampling distribution and a field propLL
for their log
probabilities (densities). savedWeights will contain the
difference in log probability (density) between the model and the
propLL value provided for each set of values.
# Accepting modelValues as a setup argument
swConf <- modelValuesConf(vars = "w",
types = "double",
sizes = 1)
setup = function(propModelValues, model, savedWeightsConf){
# Building a modelValues in the setup function
savedWeights <- modelValues(conf = savedWeightsConf)
# List of nodes to be used in run function
modelNodes <- model$getNodeNames(stochOnly = TRUE,
includeData = FALSE)
}The simplest way to pass values back and forth between models and
modelValues inside of a nimbleFunction is with copy, which
has the synonym nimCopy. See help(nimCopy) for argument details.
Alternatively, the values may be accessed via indexing of individual
rows, using the notation mv[var, i], where mv is a
modelValues object, var is a variable name (not a node name),
and i is a row number. Likewise, the getsize and
resize functions can be used as discussed in Section 14.1. However the function
as.matrix does not work in run code.
Here is a run function to use these modelValues:
run = function(){
# gets the number of rows of propSamples
m <- getsize(propModelValues)
# resized savedWeights to have the proper rows
resize(savedWeights, m)
for(i in 1:m){
# Copying from propSamples to model.
# Node names of propSamples and model must match!
nimCopy(from = propModelValues, to = model, row = i,
nodes = modelNodes, logProb = FALSE)
# calculates the log likelihood of the model
targLL <- model$calculate()
# retrieves the saved log likelihood from the proposed model
propLL <- propModelValues["propLL",i][1]
# saves the importance weight for the i-th sample
savedWeights["w", i][1] <<- exp(targLL - propLL)
}
# does not return anything
}Once the nimbleFunction is built, the modelValues object can be accessed
using $, which is shown in more detail below. In
fact, since modelValues, like most NIMBLE objects, are reference class
objects, one can get a reference to them before the function is
executed and then use that reference afterwards.
# simple model and modelValues for example use with code above
targetModelCode <- nimbleCode({
x ~ dnorm(0,1)
for(i in 1:4)
y[i] ~ dnorm(0,1)
})
# code for proposal model
propModelCode <- nimbleCode({
x ~ dnorm(0,2)
for(i in 1:4)
y[i] ~ dnorm(0,2)
})
# creating the models
targetModel = nimbleModel(targetModelCode, check = FALSE)
propModel = nimbleModel(propModelCode, check = FALSE)
cTargetModel = compileNimble(targetModel)
cPropModel = compileNimble(propModel)
sampleMVConf = modelValuesConf(vars = c("x", "y", "propLL"),
types = c("double", "double", "double"),
sizes = list(x = 1, y = 4, propLL = 1) )
sampleMV <- modelValues(sampleMVConf)
# nimbleFunction for generating proposal sample
PropSamp_Gen <- nimbleFunction(
setup = function(mv, propModel){
nodeNames <- propModel$getNodeNames()
},
run = function(m = integer() ){
resize(mv, m)
for(i in 1:m){
propModel$simulate()
nimCopy(from = propModel, to = mv, nodes = nodeNames, row = i)
mv["propLL", i][1] <<- propModel$calculate()
}
}
)
# nimbleFunction for calculating importance weights
# uses setup and run functions as defined in previous code chunk
impWeights_Gen <- nimbleFunction(setup = setup,
run = run)
# making instances of nimbleFunctions
# note that both functions share the same modelValues object
RPropSamp <- PropSamp_Gen(sampleMV, propModel)
RImpWeights <- impWeights_Gen(sampleMV, targetModel, swConf)
# compiling
CPropSamp <- compileNimble(RPropSamp, project = propModel)
CImpWeights <- compileNimble(RImpWeights, project = targetModel)
# generating and saving proposal sample of size 10
CPropSamp$run(10)
# calculating the importance weights and saving to mv
CImpWeights$run()
# retrieving the modelValues objects
# extracted objects are C-based modelValues objects
savedPropSamp_1 = CImpWeights$propModelValues
savedPropSamp_2 = CPropSamp$mv
# Subtle note: savedPropSamp_1 and savedPropSamp_2
# both provide interface to the same compiled modelValues objects!
# This is because they were both built from sampleMV.
savedPropSamp_1["x",1]## [1] -0.3919851## [1] -0.3919851## [1] 0# viewing the saved importance weights
savedWeights <- CImpWeights$savedWeights
unlist(savedWeights[["w"]])## [1] 0.5024644 1.7765357 3.5155732 0.3893390 0.5530012 0.5085069 2.3890876
## [8] 0.6548453 0.5940573 0.4943950# viewing first 3 rows -- note that savedPropSsamp_1["x", 1] was altered
as.matrix(savedPropSamp_1)[1:3, ]## propLL[1] x[1] y[1] y[2] y[3] y[4]
## [1,] -4.951100 0.00000000 -0.3951785 -1.0040021 -0.6420246 0.5993671
## [2,] -7.476891 0.35317374 -0.4510107 1.7080201 1.0471602 -0.5225390
## [3,] -8.841966 -0.06913039 -2.2884703 0.6767462 0.3955090 -0.3519276Importance sampling could also be written using simple vectors for the weights, but we illustrated putting them in a modelValues object along with model variables.
15.4.5 Using model variables and modelValues in expressions
Each way of accessing a variable, node, or modelValues can be used amidst mathematical expressions, including with indexing, or passed to another nimbleFunction as an argument. For example, the following two statements would be valid:
if Z is a vector or matrix, and
if B is a vector or matrix.
The NIMBLE language allows scalars, but models defined from BUGS code
are never created as purely
scalar nodes. Instead, a single node such as defined by z ~ dnorm(0, 1) is implemented as a vector of length 1, similar to R.
When using z via model$z or model[["z"]], NIMBLE
will try to do the right thing by treating this as a scalar. In the
event of problems33, a more explicit way to
access z is model$z[1] or model[["z"]][1].
15.4.6 Including other methods in a nimbleFunction
Other methods can be included with the methods argument to
nimbleFunction. These methods can use the objects created in
setup code in just the same ways as the run function. In
fact, the run function is just a default main method name. Any method can then call another method.
methodsDemo <- nimbleFunction(
setup = function() {sharedValue <- 1},
run = function(x = double(1)) {
print("sharedValues = ", sharedValue, "\n")
increment()
print("sharedValues = ", sharedValue, "\n")
A <- times(5)
return(A * x)
returnType(double(1))
},
methods = list(
increment = function() {
sharedValue <<- sharedValue + 1
},
times = function(factor = double()) {
return(factor * sharedValue)
returnType(double())
}))
methodsDemo1 <- methodsDemo()
methodsDemo1$run(1:10)## sharedValues = 1
##
## sharedValues = 2## [1] 10 20 30 40 50 60 70 80 90 100## sharedValues = 1
##
## sharedValues = 2## [1] 10 20 30 40 50 60 70 80 90 10015.4.7 Using other nimbleFunctions
One nimbleFunction can use another nimbleFunction that was passed to
it as a setup argument or was created in the setup function. This can
be an effective way to program. When a nimbleFunction needs to
access a setup variable or method of another nimbleFunction, use
$.
usePreviousDemo <- nimbleFunction(
setup = function(initialSharedValue) {
myMethodsDemo <- methodsDemo()
},
run = function(x = double(1)) {
myMethodsDemo$sharedValue <<- initialSharedValue
print(myMethodsDemo$sharedValue)
A <- myMethodsDemo$run(x[1:5])
print(A)
B <- myMethodsDemo$times(10)
return(B)
returnType(double())
})
usePreviousDemo1 <- usePreviousDemo(2)
usePreviousDemo1$run(1:10)## 2
## sharedValues = 2
##
## sharedValues = 3
##
## 15 30 45 60 75## [1] 30## 2
## sharedValues = 2
##
## sharedValues = 3
##
## 15
## 30
## 45
## 60
## 75## [1] 3015.4.8 Virtual nimbleFunctions and nimbleFunctionLists
Often it is useful for one nimbleFunction to have a list of other nimbleFunctions, all of whose methods have the same arguments and return types. For example, NIMBLE’s MCMC engine contains a list of samplers that are each nimbleFunctions.
To make such a list, NIMBLE provides a way to declare the arguments
and return types of methods: virtual nimbleFunctions created by
nimbleFunctionVirtual. Other nimbleFunctions can inherit from
virtual nimbleFunctions, which in R is called ‘containing’ them.
Readers familiar with object oriented programming will recognize this
as a simple class inheritance system. In Version 1.2.1 it is limited to
simple, single-level inheritance.
Here is how it works:
baseClass <- nimbleFunctionVirtual(
run = function(x = double(1)) {returnType(double())},
methods = list(
foo = function() { returnType(double()) }
))
derived1 <- nimbleFunction(
contains = baseClass,
setup = function(){},
run = function(x = double(1)) {
print("run 1")
return(sum(x))
returnType(double())
},
methods = list(
foo = function() {
print("foo 1")
return(rnorm(1, 0, 1))
returnType(double())
}
)
)
derived2 <- nimbleFunction(
contains = baseClass,
setup = function(){},
run = function(x = double(1)) {
print("run 2")
return(prod(x))
returnType(double())
},
methods = list(
foo = function() {
print("foo 2")
return(runif(1, 100, 200))
returnType(double())
}
)
)
useThem <- nimbleFunction(
setup = function() {
nfl <- nimbleFunctionList(baseClass)
nfl[[1]] <- derived1()
nfl[[2]] <- derived2()
},
run = function(x = double(1)) {
for(i in seq_along(nfl)) {
print( nfl[[i]]$run(x) )
print( nfl[[i]]$foo() )
}
}
)
useThem1 <- useThem()
set.seed(1)
useThem1$run(1:5) ## run 1
## 15
## foo 1
## -0.6264538
## run 2
## 120
## foo 2
## 157.2853## run 1
## 15
## foo 1
## -0.626454
## run 2
## 120
## foo 2
## 157.285One can also use seq_along with nimbleFunctionLists (and only with nimbleFunctionLists). As in R, seq_along(myFunList) is equivalent to
1:length(myFunList) if the length of myFunList is greater
than zero. It is an empty sequence if the length is zero.
Virtual nimbleFunctions cannot define setup values to be inherited.
15.4.9 Character objects
NIMBLE provides limited uses of character objects in run code.
Character vectors created in setup code will be available in
run code, but the only thing you can really do with them is
include them in a print or stop statement.
Note that character vectors of model node and variable names are
processed during compilation. For example, in model[[node]], node
may be a character object, and the NIMBLE compiler processes this
differently than print("The node name was ", node). In the
former, the NIMBLE compiler sets up a C++ pointer directly to the
node in the model, so that the character content of
node is never needed in C++. In the latter, node is used as
a C++ string and therefore is needed in C++.
15.4.10 User-defined data structures
Before the introduction of nimbleLists in Version 0.6-4, NIMBLE did not explicitly have user-defined data structures. An alternative way to create a data structure in NIMBLE is to use nimbleFunctions to achieve a similar effect. To do so, one can define setup code with whatever variables are wanted and ensure
they are compiled using setupOutputs. Here is an example:
dataNF <- nimbleFunction(
setup = function() {
X <- 1
Y <- as.numeric(c(1, 2))
Z <- matrix(as.numeric(1:4), nrow = 2)
setupOutputs(X, Y, Z)
})
useDataNF <- nimbleFunction(
setup = function(myDataNF) {},
run = function(newX = double(), newY = double(1), newZ = double(2)) {
myDataNF$X <<- newX
myDataNF$Y <<- newY
myDataNF$Z <<- newZ
})
myDataNF <- dataNF()
myUseDataNF <- useDataNF(myDataNF)
myUseDataNF$run(as.numeric(100), as.numeric(100:110),
matrix(as.numeric(101:120), nrow = 2))
myDataNF$X## [1] 100## [1] 100 101 102 103 104 105 106 107 108 109 110## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## [1,] 101 103 105 107 109 111 113 115 117 119
## [2,] 102 104 106 108 110 112 114 116 118 120## [1] 100nimbleOptions(buildInterfacesForCompiledNestedNimbleFunctions = TRUE)
CmyUseDataNF <- compileNimble(myUseDataNF)
CmyUseDataNF$run(-100, -(100:110), matrix(-(101:120), nrow = 2))
CmyDataNF <- CmyUseDataNF$myDataNF
CmyDataNF$X## [1] -100## [1] -100 -101 -102 -103 -104 -105 -106 -107 -108 -109 -110## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## [1,] -101 -103 -105 -107 -109 -111 -113 -115 -117 -119
## [2,] -102 -104 -106 -108 -110 -112 -114 -116 -118 -120You’ll notice that:
- After execution of the compiled function, access to
the
X,Y, andZis the same as for the uncompiled case. This occurs becauseCmyUseDataNFis an interface to the compiled version ofmyUseDataNF, and it provides access to member objects and functions. In this case, one member object ismyDataNF, which is an interface to the compiled version ofmyUseDataNF$myDataNF, which in turn provides access toX,Y, andZ. To reduce memory use, NIMBLE defaults to not providing full interfaces to nested nimbleFunctions likemyUseDataNF$myDataNF. In this example we made it provide full interfaces by setting thebuildInterfacesForCompiledNestedNimbleFunctionsoption vianimbleOptionsto TRUE. If we had left that option FALSE (its default value), we could still get to the values of interest usingvalueInCompiledNimbleFunction(CmyDataNF, 'X') - We need to take care that at the time of compilation, the
X,YandZvalues contain doubles viaas.numericso that they are not compiled as integer objects. - The
myDataNFcould be created in the setup code. We just provided it as a setup argument to illustrate that option.
15.5 Example: writing user-defined samplers to extend NIMBLE’s MCMC engine
One important use of nimbleFunctions is to write additional samplers that can be used in NIMBLE’s MCMC engine. This allows a user to write a custom sampler for one or more nodes in a model, as well as for programmers to provide general samplers for use in addition to the library of samplers provided with NIMBLE.
The following code illustrates how a NIMBLE developer would implement and use a Metropolis-Hastings random walk sampler with fixed proposal standard deviation.
my_RW <- nimbleFunction(
contains = sampler_BASE,
setup = function(model, mvSaved, target, control) {
# proposal standard deviation
scale <- if(!is.null(control$scale)) control$scale else 1
calcNodes <- model$getDependencies(target)
},
run = function() {
# initial model logProb
model_lp_initial <- getLogProb(model, calcNodes)
# generate proposal
proposal <- rnorm(1, model[[target]], scale)
# store proposal into model
model[[target]] <<- proposal
# proposal model logProb
model_lp_proposed <- model$calculate(calcNodes)
# log-Metropolis-Hastings ratio
log_MH_ratio <- model_lp_proposed - model_lp_initial
# Metropolis-Hastings step: determine whether or
# not to accept the newly proposed value
u <- runif(1, 0, 1)
if(u < exp(log_MH_ratio)) jump <- TRUE
else jump <- FALSE
# keep the model and mvSaved objects consistent
if(jump) copy(from = model, to = mvSaved, row = 1,
nodes = calcNodes, logProb = TRUE)
else copy(from = mvSaved, to = model, row = 1,
nodes = calcNodes, logProb = TRUE)
},
methods = list( reset = function () {} )
)The name of this sampler function, for the purposes of using it in an
MCMC algorithm, is my_RW. Thus, this sampler can be added
to an exisiting MCMC configuration object conf using:
To be used within the MCMC engine, sampler functions definitions must adhere exactly to the following:
- The nimbleFunction must include the contains statement
contains = sampler_BASE. - The
setupfunction must have the four argumentsmodel, mvSaved, target, control, in that order. - The
runfunction must accept no arguments, and have no return value. Further, after execution it must leave themvSavedmodelValues object as an up-to-date copy of the values and logProb values in the model object. - The nimbleFunction must have a member method called
reset, which takes no arguments and has no return value.
The purpose of the setup function is generally two-fold. First,
to extract control parameters from the control list; in the
example, the proposal standard deviation scale. It is good
practice to specify default values for any control parameters that are
not provided in the control argument, as done in the example. Second, to
generate any sets of nodes needed in the run function. In many
sampling algorithms, as here, calcNodes is used to represent the
target node(s) and dependencies up to the first layer of
stochastic nodes, as this is precisely what is required for
calculating the Metropolis-Hastings acceptance probability. These
probability calculations are done using model$calculate(calcNodes).
The purpose of the mvSaved modelValues object is to store the state of the model, including both node values and log probability values, as it exists before any changes are made by the sampler. This allows restoration of the state when a proposal is rejected, as can be seen in the example above. When a proposal is accepted, one should copy from the model into the mvSaved object. NIMBLE’s MCMC engine expects that mvSaved contains an up-to-date copy of model values and logProb values at the end of the run code of a sampler.
Note that NIMBLE generally expects the user-defined sampler to be defined in the global environment. If you define it in a function (which would generally be the case if you are using it in the context of parallelization), one approach would be to assign the user-defined sampler to the global environment in your function:
15.5.1 User-defined samplers and posterior predictive nodes
As of version 0.13.0, NIMBLE’s handling of posterior predictive nodes in MCMC sampling has changed in order to improve MCMC mixing. Samplers for nodes that are not posterior predictive nodes no longer condition on the values of the posterior predictive nodes. This is done by two changes:
- turning off posterior predictive nodes as dependencies in the use of
getDependencieswhen building an MCMC, and - moving all posterior predictive samplers to be last in the order of samplers used in an MCMC iteration.
The first change calls for careful consideration when writing new samplers. It is done by making buildMCMC set a NIMBLE system option (described below) that tells getDependencies to ignore posterior predictive nodes (defined as nodes that are themselves not data and have no data nodes in their entire downstream (descendant) dependency network) before it builds the samplers (i.e., runs the setup code of each sampler). That way, the setup code of all samplers that use getDependencies will automatically comply with the new system.
User-defined samplers that determine dependencies using getDependencies should automatically work in the new system (although it is worth checking). However if a user-defined sampler manually specifies dependencies, and if those dependencies include posterior predictive nodes, the MCMC results could be incorrect. Whether incorrect results occur could depend on other parts of the model structure connected to node(s) sampled by the user-defined sampler (i.e., the target node(s)) and/or posterior predictive nodes. Therefore, it is strongly recommended that all user-defined samplers rely on getDependencies to determine which nodes depend on the target node(s).
There are several new NIMBLE system options available to take control of the behavior just described.
getDependenciesIncludesPredictiveNodesdetermines whether posterior predictive nodes are included in results ofgetDependencies. This defaults toTRUE, so that outside of building samplers, the behavior ofgetDependenciesshould be identical to previous versions of NIMBLE.MCMCusePredictiveDependenciesInCalculationsgives the value to which
getDependenciesIncludesPredictiveNodeswill be set while building samplers. This defaults toFALSE.MCMCorderPosteriorPredictiveSamplersLastdetermines whether posterior predictive samplers are moved to the end of the sampler list. This default toTRUE. Behavior prior to version 0.13.0 corresponds toFALSE.
Here are some examples of how to use these options:
- If you want to have MCMC samplers condition on posterior predictive nodes, do
nimbleOptions(MCMCusePredictiveDependenciesInCalculations = TRUE). - If you want to get identical MCMC samplers, including their ordering, as in versions prior to 0.13.0, do
nimbleOptions(MCMCusePredictiveDependenciesInCalculations = TRUE)andnimbleOptions(MCMCorderPosteriorPredictiveSamplersLast = FALSE). - If you want to experiment with the behavior of
getDependencieswhen ignoring posterior predictive nodes, donimbleOptions(getDependenciesIncludesPredictiveNodes = FALSE).
15.6 Copying nimbleFunctions (and NIMBLE models)
NIMBLE relies heavily on R’s reference class system. When models, modelValues, and nimbleFunctions with setup code are created, NIMBLE generates a new, customized reference class definition for each. As a result, objects of these types are passed by reference and hence modified in place by most NIMBLE operations. This is necessary to avoid a great deal of copying and returning and having to reassign large objects, both in processing models and nimbleFunctions and in running algorithms.
One cannot generally copy NIMBLE models or nimbleFunctions
(specializations or generators) in a safe fashion, because of the
references to other objects embedded within NIMBLE objects. However,
the model member function newModel will create a new copy of
the model from the same model definition
(Section 6.1.3). This new model can then be used with
newly instantiated nimbleFunctions.
The reliable way to create new copies of nimbleFunctions is to re-run
the R function called nimbleFunction and record the result in a
new object. For example, say you have a nimbleFunction called
foo and 1000 instances of foo are compiled as part of an
algorithm related to a model called model1. If you then need to use foo in
an algorithm for another model, model2, doing so may work without
any problems. However, there are cases where the NIMBLE compiler will
tell you during compilation that the second set of foo instances
cannot be built from the previous compiled version. A solution is to
re-define foo from the beginning – i.e. call nimbleFunction
again – and then proceed with building and compiling the algorithm
for model2.
15.7 Debugging nimbleFunctions
One of the main reasons that NIMBLE provides an R (uncompiled) version
of each nimbleFunction is for debugging. One can call debug on
nimbleFunction methods (in particular the main run method, e.g., debug(mynf$run) and
then step through the code in R using R’s debugger. One can also
insert browser calls into run code and then run the
nimbleFunction from R.
In contrast, directly debugging a compiled nimbleFunction is
difficult, although those familiar with running R through a debugger
and accessing the underlying C code may be able to operate similarly
with NIMBLE code.
We often resort to using print statements for debugging compiled code.
Expert users fluent in C++ may also try setting nimbleOptions(pauseAfterWritingFiles = TRUE) and adding debugging code into the generated C++ files.
15.8 Timing nimbleFunctions with run.time
If your nimbleFunctions are correct but slow to run, you can use benchmarking tools to look for bottlenecks and to compare different implementations.
If your functions are very long-running (say 1ms or more), then standard R benchmarking tools may suffice, e.g. the microbenchmark package
library(microbenchmark)
microbenchmark(myCompiledFunVersion1(1.234),
myCompiledFunVersion2(1.234)) # Beware R <--> C++ overhead!If your nimbleFunctions are very fast, say under 1ms, then microbenchmark will be inaccurate due to R-to-C++ conversion overhead (that won’t happen in your actual functions).
To get timing information in C++, NIMBLE provides a run.time function that avoids the R-to-C++ overhead.
myMicrobenchmark <- compileNimble(nimbleFunction(
run = function(iters = integer(0)){
time1 <- run.time({
for (t in 1:iters) myCompiledFunVersion1(1.234)
})
time2 <- run.time({
for (t in 1:iters) myCompiledFunVersion2(1.234)
})
return(c(time1, time2))
returnType(double(1))
}))
print(myMicroBenchmark(100000))15.9 Clearing and unloading compiled objects
Sometimes it is useful to clear all the compiled objects from a
project and unload the shared library produced by your C++ compiler.
To do so, you can use clearCompiled(obj) where obj is a
compiled object such as a compiled model or nimbleFunction (e.g., a
compiled MCMC algorithm). This will clear all compiled objects
associated with your NIMBLE project. For example, if cModel is a
compiled model, clearCompiled(cModel) will clear both the
model and all associated nimbleFunctions such as compiled MCMCs that
use that model. Be careful, use of clearCompiled can be
dangerous. There is some risk that if you have copies of the R
objects that interfaced to compiled C++ objects that have been
removed, and you attempt to use those R objects after clearing their
compiled counterparts, you will crash R. We have tried to
minimize that risk, but we can’t guarantee safe behavior.
15.10 Reducing memory usage
NIMBLE can create a lot of objects in its processing, and some of them
use R features such as reference classes that are heavy in memory
usage. We have noticed that building large models can use lots of
memory. To help alleviate this, we provide two options, which can be
controlled via nimbleOptions.
As noted above, the option
buildInterfacesForCompiledNestedNimbleFunctions defaults to FALSE,
which means NIMBLE will not build full interfaces to compiled
nimbleFunctions that ony appear within other nimbleFunctions. If you
want access to all such nimbleFunctions, use the option
buildInterfacesForCompiledNestedNimbleFunctions = TRUE. This will
use more memory but can be useful for debugging.
The option clearNimbleFunctionsAfterCompiling is more drastic, and it is
experimental, so ‘buyer beware’. This will clear much of the
contents of an uncompiled nimbleFunction object after it has been
compiled in an effort to free some memory. We expect to be able to
keep making NIMBLE more efficient – faster execution and lower memory
use – in the future.
Normally this is the value of the last evaluated code, or the argument to
return.↩︎This can be omitted if you don’t need it.↩︎
Note the use of the global assignment operator to assign into the model. This is necessary for assigning into variables from the
setupfunction, at least if you want to avoid warnings from R. These warnings come from R’s reference class system.↩︎Actually, they can be, but only for uncompiled nimbleFunctions.↩︎
Please tell us!↩︎