foreach supplied object is not found within dofuture

[m-yday]

Describe the bug
in my code below: I am experiencing a strange "object 'nu' not found" coming from the traceback soon after the foreach function is run. The thing is, nu is defined within that foreach which is going directly into the %dofuture%!

Judging by the fact that I simply merged the structure of supplied code from both progressr and dofuture to fit my needs, I am really shocked that this didn't work, and that I'm unable to figure out precisely why.

This file is found on my partner's GitHub repository under Ass2.qmd. This should make this easier to see how the error actually came up.
I've now spent a couple of hours debugging this without coming to any reasonable conclusions about how the nu value is not available to the dofuture section

The line of code in question begins "Val_error<-" within the "with_progress" block

#| label: validation-analysis

# the lengthy-execution part
n_val <- 10
nu_seq <- exp(seq(-6,2,length= n_val))

#should give us a progress bar, while computing the validation error for this
#neural net, for the given nu values.
# instead, somehow the nu value—defined in foreach—is not found? 
# it may have just been faster to process this sequentially, but i am (perhaps
# foolishly) committed to learning these. 
# I feel so close to a breakthrough I SWEAR

with_progress({
  p <- progressor(along=nu_seq)
  Val_error <- foreach(nu = nu_seq,
                       .combine = rbind) %dofuture%{
    
    #Validation analysis, per nu
    res_opt = nlm(obj_pen,theta_rand,iterlim=max_iter)
    
    res_val = af_net(valid_X,valid_Y,res_opt$estimate,m,0) 
    #we are done with the penalised values now
    error = res_val$E1
    
    p(sprintf("nu=%g",nu))
    data.frame(nu=nu, 
               error=error, 
               pid=Sys.getpid() # curiosity may kill this cat
               )
  } %seed% FALSE #no random numbers are generated within this as far as i'm aware, but double-check for consistency
},handlers=handler_pbcol())

This is a chunk within a qmd file.

Reproduce example

If I paste in the other user defined sections that will get this to run, I get this:

library(tidyverse,attach.required = T)
library(future)
library(progressr)
library(doFuture)
library(parallelly)

options(parallelly.availableCores.omit = 1)

if(supportsMulticore()) # multicore = mac, linux, etc. (outside of RStudio)
  plan(multicore) else plan(multisession) #multisession = windows and other

# to ensure that the workers are correctly shut down
on.exit(
  plan(sequential)
) 


#dummy data
dat <- data.frame(X1=runif(90,-4,4),X2=runif(90,-4,4),X3=sample(c(1,0),90,replace=T),Y1=c(rep(1,30),rep(0,60)),Y2=c(rep(0,30),rep(1,30),rep(0,30)),Y3=c(rep(0,60),rep(1,30))) |> as_tibble()


# ------------------------------- model data ----------------------------------
# matrices for use in the NN functions

m <- 4 #hidden layer size
p <- 3 #predictors
q <- 3 #responses
N <- dat |> nrow()
X <- dat |> select(  1:p) |> as.matrix(ncol=p, nrow=N)
Y <- dat |> select(p+1:q) |> as.matrix(ncol=q, nrow=N)

npars <- p*p +2*p*m +m*m +m*q +p+m+m+q


g <- function(Yhat, Y){
  N <- nrow(Y)
  C <- numeric(N)
  for(i in 1:N){
    C[i] <- -log(Yhat[i, which.max(Y[i,])])
  }
  return(sum(C)) 
}

# --- activation functions ----------------------------------------------------

# augmentation layer
sig1 <- function(z)
{
  tanh(z)
}
# first m hidden layer
sig2 <- function(z)
{
  tanh(z)
}
# second m hidden layer
sig3 <- function(z)
{
  tanh(z)
}
# output layer
sig4 <- function(z)
{
  ez <- exp(z) 
  ez/matrix(1,q,q)%*%ez
}

# --- neural network ----------------------------------------------------------

af_net <- function(X, Y, theta, m, nu)
{
  N <- X |> nrow()
  p <- X |> ncol()
  q <- Y |> ncol()
  
  # allocating theta to weight and biases
  index <- 1:(p*p)
  # to auto-feature layer
  W1 <- matrix(theta[index],p,p)    # p*p
  index <- max(index)+1:(2*p*m)
  # from auto-feature (as well as from regular input) layer
  W2 <- matrix(theta[index],2*p,m)  # 2p*m
  index <- max(index)+1:(m*m)
  W3 <- matrix(theta[index],m,m)    # m*m
  index <- max(index)+1:(m*q)
  W4 <- matrix(theta[index],m,q)    # m*q
    
  index <- max(index)+1:p
  b1 <- matrix(theta[index],p,1)    # p*1
  index <- max(index)+1:m
  b2 <- matrix(theta[index],m,1)    # m*1
  index <- max(index)+1:m
  b3 <- matrix(theta[index],m,1)    # m*1
  index <- max(index)+1:q
  b4 <- matrix(theta[index],q,1)    # q*1
  
  
  # Forward Evaluation of network
  
  ones <- matrix(1,1,N) # vector of 1s transposed 1*N
  
  # auto-features inputs  
  A0   <- t(X)                        # p*N
  # regular inputs, followed by auto-features
  A1   <- sig1(t(W1)%*%A0 +b1%*%ones) # p*N
  # from auto-feature and regular inputs
  A1.  <- rbind(t(X),A1)              # 2p*N
  A2   <- sig2(t(W2)%*%A1.+b2%*%ones) # m*N
  A3   <- sig3(t(W3)%*%A2 +b3%*%ones) # m*N
  A4   <- sig4(t(W4)%*%A3 +b4%*%ones) # q*N
  Yhat <- t(A4)                       # N*q
  
  # Calculating error
  E1 <- g(Yhat,Y)
  # Penalised error
  E2 <- E1/N -nu/N*(sum(W1^2,W2^2,W3^2,W4^2)) #MSE_penalised
  
  # returning predictions and error
  return(list(Yhat = Yhat, E1 = E1, E2 = E2))
}

set.seed(2025)

valid_split <- 0.8

selected <- sample(1:N,valid_split*N,replace = F)

train_X  <- X[ selected,,drop=F]
train_Y  <- Y[ selected,,drop=F]
valid_X  <- X[-selected,,drop=F]
valid_Y  <- Y[-selected,,drop=F]

nu <- 0

obj_pen <- function(pars)
{
  res_model <- af_net(train_X,train_Y,pars,m,nu)
  return(res_model$E2)
}

max_iter   <- 1000
theta_rand <- runif(npars,-1,1)

res_opt <- nlm(obj_pen,theta_rand,iterlim = max_iter)

res_fitted <- af_net(X,Y,res_opt$estimate,m,0)

# the lengthy-execution part
n_val <- 10
nu_seq <- exp(seq(-6,2,length= n_val))

#should give us a progress bar, while computing the validation error for this
#neural net, for the given nu values.
# instead, somehow the nu value—defined in foreach—is not found? 
# it may have just been faster to process this sequentially, but i am (perhaps
# foolishly) committed to learning these. 
# I feel so close to a breakthrough I SWEAR

with_progress({
  p <- progressor(along=nu_seq)
  Val_error <- foreach(nu = nu_seq,
                       .combine = rbind) %dofuture%{
    
    #Validation analysis, per nu
    res_opt = nlm(obj_pen,theta_rand,iterlim=max_iter)
    
    res_val = af_net(valid_X,valid_Y,res_opt$estimate,m,0) 
    #we are done with the penalised values now
    error = res_val$E1
    
    p(sprintf("nu=%g",nu))
    data.frame(nu=nu, 
               error=error, 
               pid=Sys.getpid() # curiosity may kill this cat
               )
  } %seed% FALSE #no random numbers are generated within this as far as i'm aware, but double-check for consistency
},handlers=handler_pbcol())

Expected behavior

The nu value should be one of the most accessible values to the parallelised function, so this error really does feel shocking.

Session information

R version 4.4.2 (2024-10-31)
Platform: aarch64-apple-darwin20
Running under: macOS Sequoia 15.3.1

Matrix products: default
BLAS:   /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib 
LAPACK: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.0

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

time zone: Africa/Johannesburg
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] doFuture_1.0.2    future_1.40.0     foreach_1.5.2     progressr_0.15.1  plotly_4.10.4    
 [6] parallelly_1.43.0 futureverse_0.1.0 lubridate_1.9.4   forcats_1.0.0     stringr_1.5.1    
[11] dplyr_1.1.4       purrr_1.0.4       readr_2.1.5       tidyr_1.3.1       tibble_3.2.1     
[16] ggplot2_3.5.2     tidyverse_2.0.0  

loaded via a namespace (and not attached):
 [1] generics_0.1.3      stringi_1.8.7       listenv_0.9.1       hms_1.1.3          
 [5] digest_0.6.37       magrittr_2.0.3      evaluate_1.0.3      grid_4.4.2         
 [9] timechange_0.3.0    RColorBrewer_1.1-3  iterators_1.0.14    fastmap_1.2.0      
[13] jsonlite_2.0.0      httr_1.4.7          crosstalk_1.2.1     viridisLite_0.4.2  
[17] scales_1.4.0        codetools_0.2-20    lazyeval_0.2.2      cli_3.6.5          
[21] crayon_1.5.3        rlang_1.1.6         future.apply_1.11.3 yaml_2.3.10        
[25] withr_3.0.2         tools_4.4.2         parallel_4.4.2      tzdb_0.5.0         
[29] globals_0.17.0      vctrs_0.6.5         R6_2.6.1            lifecycle_1.0.4    
[33] htmlwidgets_1.6.4   pkgconfig_2.0.3     pillar_1.10.2       gtable_0.3.6       
[37] glue_1.8.0          data.table_1.17.0   highr_0.11          xfun_0.52          
[41] tidyselect_1.2.1    rstudioapi_0.17.1   knitr_1.50          farver_2.1.2       
[45] htmltools_0.5.8.1   labeling_0.4.3      compiler_4.4.2 

[HenrikBengtsson]

Thanks for reporting. There has been reports on future 1.40.0 introducing regression bugs, where some users have reported on Futureverse failing to pick up some variable, whereas it worked in future 1.34.0 (#778). If this is affecting you too, then:

1. Retry with future 1.34.0;

https://cran.r-project.org/src/contrib/Archive/future/future_1.34.0.tar.gz

If it works with future 1.34.0, then it's likely that this is the same bug. If so, then continue trying with:

2. Retry with the development versions of the globals and the future packages;

remotes::install_github("futureverse/globals", ref = "develop")
remotes::install_github("futureverse/future", ref = "develop")

If that works, then "awesome!", because then I know the bug has been fixed.

If the above does not do the trick, please try to minimize your reproducible example under "Reproduce example", so I don't have to spend time doing it. For example, I doubt the progressr code makes a difference, so you should be able to remove all that code. There's also lots of meaty pieces of code that I think could start peeling off one by one, while making sure the bug is still there. Eventually, you should end up with a minimal reproducible example that I can work with, e.g. check if it is a mistake in the code, or a bug in the futureverse code.

[m-yday]

Thank you for that! Unfortunately both things above that didn't work. So it's not the same issue. hrm. I did want to apologise for the massive code example – I had typed this up at 3am and was trying to get something viable out. Now that I'm decently awake, I'm working on a minimal reproducible example for ya :)

thank you for your massively timely response, I'm very impressed with that, and the help you've already offered!

Minimal Reproducible Example
Code now runs next-to instantly on my laptop! <3

# --- SET UP ------------------------------------------
library(doFuture)
library(future)
#install.packages("~/Downloads/future_1.34.0.tar.gz",repos=NULL,type="source")
#remotes::install_github("futureverse/globals", ref = "develop")
#remotes::install_github("futureverse/future", ref = "develop")

plan(multisession)

set.seed(2025)
m <- p <- q <- 1
npars <- 9
N <- 50

X <- as.matrix(runif(N,-4,4), nrow=N)
Y <- as.matrix(c(rep(1,N/2),rep(0,N-N/2)), nrow=N)

af_net <- function(X, Y, theta, m, nu)
{
  N <- X |> nrow()
  p <- X |> ncol()
  q <- Y |> ncol()
  
  W1 <- matrix(theta[1]  ,1,1)    
  W2 <- matrix(theta[2:3],2,1)
  W3 <- matrix(theta[4]  ,1,1)
  W4 <- matrix(theta[5]  ,1,1)
  b1 <- matrix(theta[6]  ,1,1)
  b2 <- matrix(theta[7]  ,1,1)
  b3 <- matrix(theta[8]  ,1,1)
  b4 <- matrix(theta[9]  ,1,1)
  
  ones <- matrix(1,1,N) 
  
  A0     <- t(X)               
  A1     <- tanh(t(W1)%*%A0 +b1%*%ones)
  A1.    <- rbind(t(X),A1) 
  A2     <- tanh(t(W2)%*%A1.+b2%*%ones) 
  A3     <- tanh(t(W3)%*%A2 +b3%*%ones)
  A4<-1/(1+exp(-(t(W4)%*%A3 +b4%*%ones))) #logistic activation
  Yhat <- t(A4)
  
  E1 <- sum((Yhat-Y)^2)
  E2 <- E1/N -nu/N*(sum(W1^2,W2^2,W3^2,W4^2))
  
  return(list(Yhat = Yhat, E1 = E1, E2 = E2))
}

selected <- sample(1:N,0.8*N,replace = F)

train_X  <- X[ selected,,drop=F]
train_Y  <- Y[ selected,,drop=F]
valid_X  <- X[-selected,,drop=F]
valid_Y  <- Y[-selected,,drop=F]

nu <- 1
obj_pen <- function(pars)
{
  res_model <- af_net(train_X,train_Y,pars,m,nu)
  return(res_model$E2)
}
max_iter   <- 100
theta_rand <- runif(npars,-1,1)
res_opt    <- nlm(obj_pen,theta_rand,iterlim = max_iter)
res_fitted <- af_net(X,Y,res_opt$estimate,m,0)

n_val  <- 5
nu_seq <- exp(seq(-6,2,length= n_val))

# --- CRITICAL COMPONENT  -------------------------------

Val_error <- foreach(nu = nu_seq,.combine = rbind) %dofuture% {
  
  res_opt <- nlm(obj_pen,theta_rand,iterlim=max_iter)
  res_val <- af_net(valid_X,valid_Y,res_opt$estimate,m,0) 
  error   <- res_val$E1
  
  data.frame(nu=nu,  error=error, pid=Sys.getpid())
} %seed% FALSE 

It appears that the issue is present with the multisession plan. Unsure how to test multicore with a qmd file. It appears to me that qmds do not support multicore.
Cannot rule out that it is possible in multicore, but can definitely confirm the problem is within multisession plans.

note: trying the alternative

registerDoFuture()
plan(multisession)

with %dopar% results in the same issue.

And as a final double-check,
using doParallel directly with %dopar%:

library(doParallel)
registerDoParallel()

does in fact allow the function to generate results properly.

[HenrikBengtsson]

You have the objective function:

obj_pen <- function(pars)
{
  res_model <- af_net(train_X,train_Y,pars,m,nu)
  return(res_model$E2)
}

that you minimize via:

res_opt <- nlm(obj_pen,theta_rand,iterlim=max_iter)

If you read help("nlm"), you'll find argument ... described as "additional arguments to be passed to f."

In other words, the best-practices way is to pass all variables that obj_pen() needs (i.e. train_X, train_Y, m nu) as arguments to it. Your current approach is to rely on those variables being available as global variables. That is an okay approach and works, but it's better to be explicit in every step what you variables are used and where.

Then reason your current code complaints about nu is that it is both (i) a global variable (nu <- 1) and (ii) local variable (foreach(nu = nu_seq, ...)). That causes futureverse to get confused.

In the end of the day, do you want to fit:

res_model <- af_net(train_X,train_Y,pars,m,nu)

using the global, constant nu <- 1, or do you want to use the nu values that foreach() iterates over?

[m-yday]

Aha, typically when optimising functions, I don't expect the globals to be those used. I attempted to make it more explicit by adding a nu <- nu in the body of the parallel process, but this too gave an issue. I thus assumed it wasn't an issue of passing params to nlm, since it was erroring out before that (or so I believe it is?)
I want to use the nu values that foreach iterates over, as cross-validation to tune the hyperparameter nu. The global nu is meant for setting the selected nu for the final value, based on what the cross-validation information suggests.

I didn't realise this was not best practice, as this is the way I was guided by my lecturer, but I do massively appreciate the clarification. But my question: would passing the nu values to nlm ensure it is able to be parallelised correctly? That is, one nlm function (with one nu value) per thread?

also: I should clarify I assumed I would not need to pass the nu_seq to the dofuture explicitly, due to it being part of the foreach (each thread requiring one nu only). Is this mistaken?