```
# Start up a parallel cluster
parallelCluster <- parallel::makeCluster(parallel::detectCores())
print(parallelCluster)
```

`## socket cluster with 4 nodes on host 'localhost'`

```
d <- iris # let "d" refer to one of R's built in data sets
vars <- c('Sepal.Length','Sepal.Width','Petal.Length')
yName <- 'Species'
yLevels <- sort(unique(as.character(d[[yName]])))
fitOneTargetModel <- function(yName,yLevel,vars,data) {
formula <- paste('(',yName,'=="',yLevel,'") ~ ',
paste(vars,collapse=' + '),sep='')
glm(as.formula(formula),family=binomial,data=data)
}
```

```
for(yLevel in yLevels) {
print("*****")
print(yLevel)
print(fitOneTargetModel(yName,yLevel,vars,d))
}
```

```
## [1] "*****"
## [1] "setosa"
```

`## Warning: glm.fit: algorithm did not converge`

`## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred`

```
##
## Call: glm(formula = as.formula(formula), family = binomial, data = data)
##
## Coefficients:
## (Intercept) Sepal.Length Sepal.Width Petal.Length
## -71.80 23.91 13.51 -34.95
##
## Degrees of Freedom: 149 Total (i.e. Null); 146 Residual
## Null Deviance: 191
## Residual Deviance: 3.523e-09 AIC: 8
## [1] "*****"
## [1] "versicolor"
##
## Call: glm(formula = as.formula(formula), family = binomial, data = data)
##
## Coefficients:
## (Intercept) Sepal.Length Sepal.Width Petal.Length
## 7.8238 0.3475 -3.3792 -0.1302
##
## Degrees of Freedom: 149 Total (i.e. Null); 146 Residual
## Null Deviance: 191
## Residual Deviance: 151.5 AIC: 159.5
## [1] "*****"
## [1] "virginica"
```

`## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred`

```
##
## Call: glm(formula = as.formula(formula), family = binomial, data = data)
##
## Coefficients:
## (Intercept) Sepal.Length Sepal.Width Petal.Length
## -38.2154 -3.8524 -0.6385 13.1475
##
## Degrees of Freedom: 149 Total (i.e. Null); 146 Residual
## Null Deviance: 191
## Residual Deviance: 23.77 AIC: 31.77
```

```
worker <- function(yLevel) {
fitOneTargetModel(yName,yLevel,vars,d)
}
models <- lapply(yLevels,worker)
```

`## Warning: glm.fit: algorithm did not converge`

`## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred`

`## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred`

```
names(models) <- yLevels
print(models)
```

```
## $setosa
##
## Call: glm(formula = as.formula(formula), family = binomial, data = data)
##
## Coefficients:
## (Intercept) Sepal.Length Sepal.Width Petal.Length
## -71.80 23.91 13.51 -34.95
##
## Degrees of Freedom: 149 Total (i.e. Null); 146 Residual
## Null Deviance: 191
## Residual Deviance: 3.523e-09 AIC: 8
##
## $versicolor
##
## Call: glm(formula = as.formula(formula), family = binomial, data = data)
##
## Coefficients:
## (Intercept) Sepal.Length Sepal.Width Petal.Length
## 7.8238 0.3475 -3.3792 -0.1302
##
## Degrees of Freedom: 149 Total (i.e. Null); 146 Residual
## Null Deviance: 191
## Residual Deviance: 151.5 AIC: 159.5
##
## $virginica
##
## Call: glm(formula = as.formula(formula), family = binomial, data = data)
##
## Coefficients:
## (Intercept) Sepal.Length Sepal.Width Petal.Length
## -38.2154 -3.8524 -0.6385 13.1475
##
## Degrees of Freedom: 149 Total (i.e. Null); 146 Residual
## Null Deviance: 191
## Residual Deviance: 23.77 AIC: 31.77
```

```
tryCatch(
models <- parallel::parLapply(parallelCluster,
yLevels,worker),
error = function(e) print(e)
)
```

`## <simpleError in checkForRemoteErrors(val): 3 nodes produced errors; first error: could not find function "fitOneTargetModel">`

The above fails because the worker function can not find the function “fitOneTargetModel.” In fact it can not find any of the code or data we supplied (it will also fail to find yName, vars, and d). This is because in the common default implementation of the parallel interface (“snow”): the function to be evaluated is serialized and unpacked into different processes and the Global environment is not copied over. So the running processes do not see the variable bindings we have in the Global environment.

So we need to make sure all of our variable references are neatly packed together so parallel/snow can transport them to where we need them. The cleanest way to do this is through the use of constructing our own lexical closure. If we define our worker inside another function then anything also in that function activation environment will be available to our worker. We do this using a pattern we will call a worker factory.

```
# build the single argument function we are going to pass to parallel
mkWorker <- function(yName,vars,d) {
# make sure each of the three values we need passed
# are available in this environment
force(yName)
force(vars)
force(d)
# define any and every function our worker function
# needs in this environment
fitOneTargetModel <- function(yName,yLevel,vars,data) {
formula <- paste('(',yName,'=="',yLevel,'") ~ ',
paste(vars,collapse=' + '),sep='')
glm(as.formula(formula),family=binomial,data=data)
}
# finally: define and return our worker function
# the point is worker's "lexical closure"
# (where it looks for unbound variables)
# is mkWorker's activation/execution environment
# and not the usual Global environment.
# the parallel library is willing to transport
# this environment (which it does not
# do for the Global environment).
worker <- function(yLevel) {
fitOneTargetModel(yName,yLevel,vars,d)
}
return(worker)
}
models <- parallel::parLapply(parallelCluster,yLevels,
mkWorker(yName,vars,d))
names(models) <- yLevels
print(models)
```

```
## $setosa
##
## Call: glm(formula = as.formula(formula), family = binomial, data = data)
##
## Coefficients:
## (Intercept) Sepal.Length Sepal.Width Petal.Length
## -71.80 23.91 13.51 -34.95
##
## Degrees of Freedom: 149 Total (i.e. Null); 146 Residual
## Null Deviance: 191
## Residual Deviance: 3.523e-09 AIC: 8
##
## $versicolor
##
## Call: glm(formula = as.formula(formula), family = binomial, data = data)
##
## Coefficients:
## (Intercept) Sepal.Length Sepal.Width Petal.Length
## 7.8238 0.3475 -3.3792 -0.1302
##
## Degrees of Freedom: 149 Total (i.e. Null); 146 Residual
## Null Deviance: 191
## Residual Deviance: 151.5 AIC: 159.5
##
## $virginica
##
## Call: glm(formula = as.formula(formula), family = binomial, data = data)
##
## Coefficients:
## (Intercept) Sepal.Length Sepal.Width Petal.Length
## -38.2154 -3.8524 -0.6385 13.1475
##
## Degrees of Freedom: 149 Total (i.e. Null); 146 Residual
## Null Deviance: 191
## Residual Deviance: 23.77 AIC: 31.77
```

The above works, but it is incredibly tedious and wasteful to have to re-define every function we need every time we need it. Any easier to use way is to use a helper function we supply called “bindToEnv” to do some of the work. With bindToEnv the code looks like the following.

```
source('bindToEnv.R') # http://winvector.github.io/Parallel/bindToEnv.R
# build the single argument function we are going to pass to parallel
mkWorker <- function() {
bindToEnv(objNames=c('yName','vars','d','fitOneTargetModel'))
worker <- function(yLevel) {
fitOneTargetModel(yName,yLevel,vars,d)
}
worker
}
models <- parallel::parLapply(parallelCluster,yLevels,
mkWorker())
names(models) <- yLevels
print(models)
```

```
## $setosa
##
## Call: glm(formula = as.formula(formula), family = binomial, data = data)
##
## Coefficients:
## (Intercept) Sepal.Length Sepal.Width Petal.Length
## -71.80 23.91 13.51 -34.95
##
## Degrees of Freedom: 149 Total (i.e. Null); 146 Residual
## Null Deviance: 191
## Residual Deviance: 3.523e-09 AIC: 8
##
## $versicolor
##
## Call: glm(formula = as.formula(formula), family = binomial, data = data)
##
## Coefficients:
## (Intercept) Sepal.Length Sepal.Width Petal.Length
## 7.8238 0.3475 -3.3792 -0.1302
##
## Degrees of Freedom: 149 Total (i.e. Null); 146 Residual
## Null Deviance: 191
## Residual Deviance: 151.5 AIC: 159.5
##
## $virginica
##
## Call: glm(formula = as.formula(formula), family = binomial, data = data)
##
## Coefficients:
## (Intercept) Sepal.Length Sepal.Width Petal.Length
## -38.2154 -3.8524 -0.6385 13.1475
##
## Degrees of Freedom: 149 Total (i.e. Null); 146 Residual
## Null Deviance: 191
## Residual Deviance: 23.77 AIC: 31.77
```

What bindToEnv does is for each value name listed is copy a reference to it into the current environment.

In addition (and this is the tricky step) bindToEnv also re-binds the lexical environment of any function it is asked to work with to the current target environment. This way two functions can each find each other (as they see each being named in the now common lexical environment). This is somewhat bad (we shouldn’t be doing directly manipulation of environment linkages), but the cat is already out of that bag as the parallel library is clearly also re-wiring environments. It also isn’t as bad as it may seem as (due to R’s immutability rules) only our copy of each function has its environment re-wired, nobody outside of our use sees this change.

What can go wrong is we may break any functions that depend on an functional idea called “Currying” (or any other clever use of closures for that matter). By us getting too clever with closures we defeat other smart uses of closures.

An example of what goes wrong is shown below.

```
# Function taking many arguments
mkPrefixXI <- function(x,i) {
paste(x,i,sep='.')
}
# Wrapper function that curries the value x into a prefix pasting
# This is the problem: our library breaks the results of helpers like CurryFn()
CurryFn <- function(fn,xarg) {
function(yarg) {
fn(xarg,yarg)
}
}
f <- CurryFn(mkPrefixXI,'x')
f(2)
```

`## [1] "x.2"`

```
mkWorker <- function() {
bindToEnv(objNames=c('mkPrefixXI','f'))
function(i) {
f(i)
}
}
tryCatch(parallel::parLapply(parallelCluster,1:3,mkWorker()),
error = function(e) print(e))
```

`## <simpleError in checkForRemoteErrors(val): 3 nodes produced errors; first error: could not find function "fn">`

```
# Our own wrapper also does not work
# Wrapper function that curries the value x into a prefix pasting
# This is the problem: our library breaks the results of helpers like mkPrefixX
mkPrefixX <- function(xarg) {
function(yarg) {
mkPrefixXI(xarg,yarg)
}
}
f <- mkPrefixX('x')
f(2)
```

`## [1] "x.2"`

```
mkWorker <- function() {
bindToEnv(objNames=c('mkPrefixXI','mkPrefixXI','f'))
function(i) {
f(i)
}
}
tryCatch(parallel::parLapply(parallelCluster,1:3,mkWorker()),
error = function(e) print(e))
```

`## <simpleError in checkForRemoteErrors(val): 3 nodes produced errors; first error: object 'xarg' not found>`

This is a problem as Currying is a very common and useful pattern (and fairly similar to what we are doing to prepare functions for use with the parallel library).

The work around is to not combine Currying and function wrapping. Always Curry directly by plugging in known arguments to a function. Do not try to use any sort of syntax shortening wrapper like CurryFn or mkPrefixX.

```
# Curry by plugging in value directly, not by a function wrapper as in last example
# Notice we type in our value x='x' directly, so it isn't coming through an environment
# (which could later stripped).
f <- function(i) { mkPrefixXI('x',i) }
mkWorker <- function() {
bindToEnv(objNames=c('mkPrefixXI','f'))
function(i) {
f(i)
}
}
parallel::parLapply(parallelCluster,1:3,mkWorker())
```

```
## [[1]]
## [1] "x.1"
##
## [[2]]
## [1] "x.2"
##
## [[3]]
## [1] "x.3"
```

Note: this overall methodology does not work with useful structures such as lists of functions, as bindToEnv doesn’t see these functions (and hence does not rebind them).

```
# Shutdown cluster neatly
if(!is.null(parallelCluster)) {
parallel::stopCluster(parallelCluster)
parallelCluster <- c()
}
```