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).

Toy model

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.

Basics: what is the final graph?

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.

The 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.

Nodes, RHSonly nodes, vertices, and elements

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:

  1. 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].

  2. 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].

  3. 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.

  4. 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.)

The user-facing concept of “nodes”

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]"

various IDs: elementIDs, vertexIDs, and graphIDs

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.

Storing IDs in environments with variables matching the model’s variables

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.

Mapping from one kind of ID to another.

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

Other content to be explained

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