Here is a sketch of design ideas for BUGS modules

This is implemented on branch BUGSmodules

library(nimble)
nimbleOptions(enableBUGSmodules = TRUE) ## feature toggle

Here is the basic idea. A declaration with BUGS module lmPred would get expanded in place:

mc <- nimbleCode({
    a ~ dnorm(0,1)
    Ypred[1:n] <- lmPred(X[1:n])
})

expands to:

nimble:::codeProcessBUGSmodules(mc)

## {
##     a ~ dnorm(0, 1)
##     {
##         X_coeff ~ dnorm(0, sd = 1000)
##         for (i in 1:n) {
##             predY[i] <- X_coeff * X[i]
##         }
##     }
## }

This provides the linear predictors from X[1:n] to predict Ypred[1:n].

When option enableBUGSmodules == TRUE, this will be done early in processing BUGS code:

X <- 1:5
m <- nimbleModel(mc, constants = list(n = 5), data = list(X = X))

## defining model...

## building model...

## setting data and initial values...

## 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: a, X_coeff, predY. This is not an error, but some or all
## variables may need to be initialized for certain algorithms to operate
## properly.

## 

## model building finished.

Note that the functionality in lmPred is currently a sham: It replaces the code but not by truly processing it.

One BUGS module can use another. Say we want a linear model:

mc <- nimbleCode({
    a ~ dnorm(0,1)
    Ypred[1:n] ~ nim_lm(X[1:n])
})

The non-recursed expansion uses the lmPred template:

nimble:::nim_lm$process(Ypred[1:n], nim_lm(X[1:n])) ## The processing is fake, to show what the real result would look like.  It generates the response distributions and uses lmPred to expand linear predictions...

## $code
## {
##     for (i in 1:n) {
##         Y[i] ~ dnorm(predY[i], sd = sigmga)
##     }
##     predY[1:n] ~ lmPred(X[1:n], priors = priors)
## }

This means the recursed expansion has everything:

nimble:::codeProcessBUGSmodules(mc)  ## ... which are then replaced by recursion

## {
##     a ~ dnorm(0, 1)
##     {
##         for (i in 1:n) {
##             Y[i] ~ dnorm(predY[i], sd = sigmga)
##         }
##         {
##             X_coeff ~ dnorm(0, sd = 1000)
##             for (i in 1:n) {
##                 predY[i] <- X_coeff * X[i]
##             }
##         }
##     }
## }

X <- 1:5
## The full expansion works in a model
m <- nimbleModel(mc, constants = list(n = 5), data = list(X = X))

## defining model...

## building model...

## setting data and initial values...

## 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: a, Y, X_coeff, predY, sigmga. This is not an error, but some or
## all variables may need to be initialized for certain algorithms to operate
## properly.

## 

## model building finished.

It would be easy to remove extra {s and collapse identical for loops into one.

Similarly, a glm module could use lmPred but provide an argument to insert a link:

mc <- nimbleCode({
    a ~ dnorm(0,1)
    Ypred[1:n] ~ nim_glm(X[1:n])
})
X <- 1:5
## non-recursed expansion, shows link argument
nimble:::nim_glm$process(Ypred[1:n], nim_glm(X[1:n])) ## In this case a link is requested in the linear expansion.

## $code
## {
##     for (i in 1:n) {
##         Y[i] ~ dexp(predY[i])
##     }
##     predY[1:n] ~ lmPred(X[1:n], priors = priors, link = "log")
## }

## recursed expansion shows full model
nimble:::codeProcessBUGSmodules(mc)  ## ... which is then replaced by recursion

## {
##     a ~ dnorm(0, 1)
##     {
##         for (i in 1:n) {
##             Y[i] ~ dexp(predY[i])
##         }
##         {
##             X_coeff ~ dnorm(0, sd = 1000)
##             for (i in 1:n) {
##                 log(predY[i]) <- X_coeff * X[i]
##             }
##         }
##     }
## }

## It works in a fully processed model
m <- nimbleModel(mc, constants = list(n = 5), data = list(X = X))

## defining model...

## building model...

## setting data and initial values...

## 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: a, Y, X_coeff, predY. This is not an error, but some or all
## variables may need to be initialized for certain algorithms to operate
## properly.

## 

## model building finished.