This document gives a guided tour of some internals of model objects in R with the goal of illustrating the concepts and data structures of a model defined by BUGS code. It does not cover the C++ representation of the final directed acyclic graph (DAG).
Here is a toy model whose only purpose is to illustrate concepts:
mc <- nimbleCode({
x[1:10] <- seq(.1, .2, length = 10)
y[1] <- sum(x[1:5])
y[2] <- sum(x[6:10])
z[1] <- sum(x[1:2])
z[2] ~ dnorm( mean = x[3] + sum(a[1:2]), sd = 1)
})
m <- nimbleModel(mc)
## defining model...
## building model...
## running calculate on model (any error reports that follow may simply
## reflect missing values in model variables) ...
##
## checking model sizes and dimensions...
## note that missing values (NAs) or non-finite values were found in model
## variables: x, y, z, a. This is not an error, but some or all variables may
## need to be initialized for certain algorithms to operate properly.
##
## model building finished.
The final graph is simple enough that the igraph
output and plot for it are useful:
m$modelDef$graph
## IGRAPH DN-- 10 11 --
## + attr: name (v/c)
## + edges (vertex names):
## [1] x[1:2] ->y[1] x[3] ->y[1] x[4:5] ->y[1] x[6:10]->y[2]
## [5] x[1:2] ->z[1] x[3] ->z[2] a[1:2] ->z[2] x[1:10]->x[1:2]
## [9] x[1:10]->x[3] x[1:10]->x[4:5] x[1:10]->x[6:10]
library(igraph)
plot(m$modelDef$graph)
Notice that some nodes were declared explicitly (x[1:10]
) while others were inferred from right-hand side (RHS) expressions (x[1:2]
). Such “LHSinferred” nodes (yes, that is what they are called) are included as children of the larger node of which they are a part. They require special handling during getDependencies
. They are also the reason for a fair bit of the complexity described below, which is why this toy model includes them.
modelDef
and maps
The modelDef
(model definition) contains all the information about a model’s structure that is not specific to one instance of the model. This is most of the interesting stuff. The exception is that data labeling can be different for difference instances and so is in the model
, not the modelDef
.
Within the modelDef
, the trove of information about nodes and their relationships is in the maps
, so-called because many of its objects provide a way to map between one kind of information and another such as a node name and a graphID
.
Let’s extract the maps
and look at its contents.
maps <- m$modelDef$maps
ls(maps)
## [1] "edgesFrom"
## [2] "edgesFrom2ParentExprID"
## [3] "edgesFrom2To"
## [4] "edgesParentExprID"
## [5] "edgesTo"
## [6] "elementID_2_vertexID"
## [7] "elementNames"
## [8] "end_IDs"
## [9] "graphID_2_declID"
## [10] "graphID_2_logProbName"
## [11] "graphID_2_nodeFunctionName"
## [12] "graphID_2_nodeName"
## [13] "graphID_2_unrolledIndicesMatrixRow"
## [14] "graphIDs"
## [15] "isEndNode_byGID"
## [16] "latent_IDs"
## [17] "nimbleGraph"
## [18] "nodeName_2_graphID"
## [19] "nodeName_2_nodeFunctionName"
## [20] "nodeNames"
## [21] "nodeNamesLHSall"
## [22] "nodeNamesRHSonly"
## [23] "notStoch"
## [24] "setPositions3"
## [25] "top_IDs"
## [26] "types"
## [27] "vars2GraphID_functions"
## [28] "vars2GraphID_functions_and_RHSonly"
## [29] "vars2GraphID_values"
## [30] "vars2ID_elements"
## [31] "vars2LogProbName"
## [32] "vertexID_2_nodeID"
Much of the content of maps
is redundant information organized for fast lookup in different ways for different needs. We’ll walk through pieces of this.
When we talk casually about a graph for nimbleFunction
purposes, we use “node” with multiple meanings. For building and understanding the model, several distinct concepts have evolved:
A “node with a nodeFunction”, or “nodeFunction” for short, is an entity that appeared as a unit on the LHS of a BUGS declaration. This means it will be calculated (simulated, etc.) as a unit via a nodeFunction. Example: x[1:10]
.
A “RHSonly node” (right-hand-side-only node) is an entity that appeared only the RHS of a BUGS declaration and was never part of something in a LHS declaration. Examples: a[1:2]
is a RHSonly node. x[1:5]
is not a RHSonly node because it is part of x[1:10]
.
A “vertex” is anything represented in the graph. This includes nodes and the smallest subsets of their elements that appear as a unit on the RHS (or their complements). Examples: x[1:2]
, x[3]
and x[4:5]
are vertices. x[1:5]
would have been a vertex if x[1:2]
and x[3]
had not been used separately, but we always split into smallest elements; notice that x[1:2]
, x[3]
, and x[4:5]
all have y[1]
as a childe. x[4:5]
never appears in the BUGS code but it is a vertex because it is the complement (what is left over) of splitting x[1:5]
into x[1:2]
and x[3]
. All vertices must use contiguous index blocks, so there can be multiple vertices created from the complement.
An “element” is any scalar element of any variable. Unfortunately, it appeared to be necessary to assign a unique elementID
to every scalar element for construction and lookup purposes.
The names and uses of maps
variables related to these concepts can be confusing, partly because there are ambiguous names (for historical reasons that perhaps we should try to improve) and partly because the concepts need to be combined in particular ways. For example, the nodeNames
includes names of vertices (including nodeFunctions and RHSonly nodes).
Let’s look at some of these:
maps$elementNames ## a name for every scalar element
## [1] "" "a[1]" "a[2]" "x[4]" "x[5]" "x[1]" "x[2]" "x[3]"
## [9] "x[6]" "x[7]" "x[8]" "x[9]" "x[10]" "z[1]" "y[1]" "z[2]"
## [17] "y[2]"
maps$nodeNames ## a name for everything in the graph
## [1] "x[1:10]" "a[1:2]" "x[4:5]" "x[1:2]" "x[3]" "x[6:10]" "z[1]"
## [8] "y[1]" "z[2]" "y[2]"
maps$nodeNamesRHSonly ## which are RHSonly
## [1] "a[1:2]"
maps$types ## labels of node types
## [1] "determ" "RHSonly" "LHSinferred" "LHSinferred" "LHSinferred"
## [6] "LHSinferred" "determ" "determ" "stoch" "determ"
cbind(maps$nodeNames, maps$types) ## types correspond to nodeNames
## [,1] [,2]
## [1,] "x[1:10]" "determ"
## [2,] "a[1:2]" "RHSonly"
## [3,] "x[4:5]" "LHSinferred"
## [4,] "x[1:2]" "LHSinferred"
## [5,] "x[3]" "LHSinferred"
## [6,] "x[6:10]" "LHSinferred"
## [7,] "z[1]" "determ"
## [8,] "y[1]" "determ"
## [9,] "z[2]" "stoch"
## [10,] "y[2]" "determ"
Note that vertices that are in the graph as if they had been on the LHS of a declaration are called “LHSinferred”. (One could argue that a better name would be “RHSinferred”, since they were inferred to be necessary based on appearance on the RHS of a BUGS declaration.)
A user shouldn’t need to think about elements vs. vertices vs. nodes. They might need to think about RHSonly vs. stoch vs. determ nodes (and also data, but that’s a label on some stochastic nodes, not a node type). Hence the following:
args(m$getNodeNames) ## fair bit of control
## function (determOnly = FALSE, stochOnly = FALSE, includeData = TRUE,
## dataOnly = FALSE, includeRHSonly = FALSE, topOnly = FALSE,
## latentOnly = FALSE, endOnly = FALSE, returnType = "names",
## returnScalarComponents = FALSE)
## NULL
m$getNodeNames() ## nodeFunctions!
## [1] "x[1:10]" "z[1]" "y[1]" "z[2]" "y[2]"
m$getNodeNames(includeRHSonly = TRUE) ## The full graph except split vertices
## [1] "x[1:10]" "a[1:2]" "z[1]" "y[1]" "z[2]" "y[2]"
A purpose of splitting vertices as described above is to find dependencies correctly. E.g. the dependencies of x[2]
should include x[1:10]
, y[1]
and z[1]
. At the time of this writing there is a bug with a pending fix that returns all dependencies of x[1:10]
. In typical use this won’t make algorithms wrong but it will make them slow by calculating more than necessary.
m$getDependencies('x[2]')
## [1] "x[1:10]" "z[1]" "y[1]" "z[2]" "y[2]"
The various IDs are important for looking up information and mapping between the different concepts.
The graphIDs
aligns with the nodeNames
and types
data.frame(maps$nodeNames, maps$graphIDs, maps$types)
## maps.nodeNames maps.graphIDs maps.types
## 1 x[1:10] 1 determ
## 2 a[1:2] 2 RHSonly
## 3 x[4:5] 3 LHSinferred
## 4 x[1:2] 4 LHSinferred
## 5 x[3] 5 LHSinferred
## 6 x[6:10] 6 LHSinferred
## 7 z[1] 7 determ
## 8 y[1] 8 determ
## 9 z[2] 9 stoch
## 10 y[2] 10 determ
The elementIDs are simply a unique integer for each scalar element. Let’s use those to illustrate the next concept about how some IDs are stored.
Variables in the model can be matrices, arrays, etc. We often need to go from character node names to various IDs in order to query graph relationships (and back to character names). Early versions of procesing to do this became bottlenecks. Eventually we settled on the following, illustrated for elementIDs
maps$vars2ID_elements ## this is an environment
## <environment: 0x7fa6108578a8>
ls(maps$vars2ID_elements) ## with variables matching those in the model
## [1] "a" "x" "y" "z"
maps$vars2ID_elements$x ## variables are real size but contain IDs
## [1] 6 7 8 4 5 9 10 11 12 13
## which lets us use R's parser and evaluator to extract IDs from character code
## of arbitrary model variable expressions
eval(parse(text = 'x[3:7]', keep.source=FALSE)[[1]], envir = maps$vars2ID_elements)
## [1] 8 4 5 9 10
See parseEvalNumericMany
for how we do this efficiently for many cases in one eval
.
Some elements of maps
are integer vectors that convert one kind of ID to another. For example, if we want to get the vertexIDs of scalar elements of x
:
maps$elementID_2_vertexID[maps$vars2ID_elements$x]
## [1] 4 4 5 3 3 6 6 6 6 6
maps$vars2ID_elements$x
## [1] 6 7 8 4 5 9 10 11 12 13
maps$elementID_2_vertexID
## [1] 1 2 2 3 3 4 4 5 6 6 6 6 6 7 8 9 10
maps$elementID_2_vertexID[maps$vars2ID_elements$x]
## [1] 4 4 5 3 3 6 6 6 6 6
maps$vertexID_2_nodeID
## [1] 1 0 1 1 1 1 7 8 9 10
maps$vertexID_2_nodeID[maps$elementID_2_vertexID[maps$vars2ID_elements$x]]
## [1] 1 1 1 1 1 1 1 1 1 1
maps$graphID_2_nodeName
## [1] "x[1:10]" "a[1:2]" "x[4:5]" "x[1:2]" "x[3]" "x[6:10]" "z[1]"
## [8] "y[1]" "z[2]" "y[2]"
maps$graphID_2_nodeFunctionName
## [1] "x[1:10]" "a[1:2]" "x[1:10]" "x[1:10]" "x[1:10]" "x[1:10]" "z[1]"
## [8] "y[1]" "z[2]" "y[2]"
maps$nodeNamesRHSonly
## [1] "a[1:2]"
maps$vars2GraphID_functions$x
## [1] 1 1 1 1 1 1 1 1 1 1
maps$vars2GraphID_values$x
## [1] 4 4 5 3 3 6 6 6 6 6
maps$vars2GraphID_functions_and_RHSonly$x
## [1] 1 1 1 1 1 1 1 1 1 1
maps$vars2GraphID_functions$a
## [1] NA NA
maps$vars2GraphID_values$a
## [1] 2 2
maps$vars2GraphID_functions_and_RHSonly$a
## [1] 2 2
maps$vars2GraphID_values$x
## [1] 4 4 5 3 3 6 6 6 6 6