0

I am trying to achieve parallel Kriging in R on several variables using a loop. Here is a reproducible example using data meuse and a code for parallel kriging that I found here. In the code below, each kriging is overwritten, but that's not the issue, since I can't even run the first kriging.

library(parallel)
library(sp)
library(gstat)
data(meuse)
    names <- colnames(meuse)[3:6]
    coordinates(meuse) = ~x+y
    data(meuse.grid)
    gridded(meuse.grid) = ~x+y
    m <- vgm(.59, "Sph", 874, .04)
    no_cores <- 7
    # ordinary kriging:
    for (i in 1:length(names)) {
      parts <- split(x = 1:length(meuse.grid), f = 1:no_cores)
cl <- makeCluster(no_cores)          
clusterExport(cl = cl, varlist = c("meuse", "meuse.grid", "m", "parts"), envir = .GlobalEnv)
          clusterEvalQ(cl = cl, expr = c(library('sp'), library('gstat')))
          parallelX <- parLapply(cl = cl, X = 1:no_cores, fun = function(x) krige(formula = log(get(names[i]))~1, locations = meuse, newdata = meuse.grid[parts[[x]],], model = m))
          stopCluster(cl)
          # Merge all the predictions    
          mergeParallelX <- maptools::spRbind(parallelX[[1]], parallelX[[2]])
          mergeParallelX <- maptools::spRbind(mergeParallelX, parallelX[[3]])
          mergeParallelX <- maptools::spRbind(mergeParallelX, parallelX[[4]])
          mergeParallelX <- maptools::spRbind(mergeParallelX, parallelX[[5]])
          mergeParallelX <- maptools::spRbind(mergeParallelX, parallelX[[6]])
          mergeParallelX <- maptools::spRbind(mergeParallelX, parallelX[[7]])
          # Create SpatialPixelsDataFrame from mergeParallelX
          mergeParallelX <- SpatialPixelsDataFrame(points = mergeParallelX, data = mergeParallelX@data)
        }

I keep having issues which are, I think, related to the dynamic formula. I tried a lot of other functions, such as paste, paste0, as.formula, objects, etc. Nothing works, I can't paste the variable names into looped parallel kriging, no matter how I code the dynamic formula. Maybe that's because of how the parallel function is coded?

Any idea about a dynamic formula that works with the above code?

5
  • What does "not working" mean here? Is this really only a Q about how to dynamically construct a formula?
    – Spacedman
    Jul 11, 2022 at 14:44
  • Thanks, I didn't know about the vocabulary I had to use for this issue. So, the answer is yes and no, apparently a formula constructed dynamically don't work with the package parallel. For example the structure you gave below doesn't work. I will edit the question and the title to reflect on this issue.
    – ePoQ
    Jul 12, 2022 at 9:14
  • As written you are trying to do n_cores replications of kriging of each of the columns? Because your loop over the columns is outside the cluster creation and execution. This seems a bit off. Normally you'd throw all the replicated work inside the cluster and let the cluster's n_cores specification optimise the number of cores working. That then puts i in side the cluster loop. Also as written you're only going to return the value from the last column's kriging.
    – Spacedman
    Jul 12, 2022 at 9:37
  • As a matter of fact, I have a double loop to record all kringing, and I use detectCores for setting the number of clusters. I wrote the code this way in the Q to have a simple and reproducible example. What doesn't work is the dynamic formula here, but if you have a suggestion for a better structure I am ears. That's the first time I try parallel computing so I am quite a newbie here.
    – ePoQ
    Jul 12, 2022 at 9:45
  • And yes, since Kriging predictions are independent of each other, I am just trying to speed up the processing time by making all my cores work at the same time for each variable successively, rather than allocating different variables to different cores.
    – ePoQ
    Jul 12, 2022 at 9:56

1 Answer 1

1

Guessing that your problem is constructing a formula dynamically, you can build one as character and convert to a formula with as.formula:

So this:

> k1 = krige(log(copper)~1, meuse, meuse.grid, model = m)
[using ordinary kriging]

is the same as:

> k2 = krige(as.formula("log(copper)~1"), meuse, meuse.grid, model = m)
[using ordinary kriging]

which you can construct

> names
[1] "cadmium" "copper"  "lead"    "zinc"   
> i
[1] 2
> k3 = krige(as.formula(paste0("log(",names[i],")~1")), meuse, meuse.grid, model = m)
[using ordinary kriging]

And then:

> identical(k1, k2)
[1] TRUE
> identical(k1, k3)
[1] TRUE
> identical(k2, k3)
[1] TRUE
1
  • I just edited the question/title. as.formula(paste0("log(",names[i],")~1")) doesn't work either with the package parallel
    – ePoQ
    Jul 12, 2022 at 9:24

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.