R中的条件模拟(使用Kriging)是否具有并行性？

Question

我在R中使用gstat包来生成顺序高斯模拟。我的pc有4个内核，我尝试按照Guzmán提供的回答问题How to achieve parallel Kriging in R to speed up the process?的脚本，使用并行包对krige()函数进行并行化。

然而，由此产生的模拟不同于当时只使用一个核心的模拟(没有并行化)。它看起来像是一个几何问题，但我找不到如何解决它。

接下来，我将提供一个生成2个模拟的示例(使用4个内核)。您将看到，在运行代码后，从并行化派生的模拟映射显示了一些伪像(如垂直线)，并且与当时只使用一个内核的映射不同。

代码需要库gstat、sp、raster、parallel和spatstat。如果有任何library()行不起作用，请先运行install.packages()。

library(gstat)
library(sp)
library(raster)
library(parallel)
library(spatstat)

# create a regular grid
nx=100 # number of columns
ny=100 # number of rows
srgr <- expand.grid(1:ny, nx:1)
names(srgr) <- c('x','y')
gridded(srgr)<-~x+y

# generate a spatial process (unconditional simulation)
g<-gstat(formula=z~x+y, locations=~x+y, dummy=T, beta=15, model=vgm(psill=3, range=10, nugget=0,model='Exp'), nmax=20)
sim <- predict(g, newdata=srgr, nsim=1)
r<-raster(sim)

# generate sample data (Poisson process)  
int<-0.02
rpp<-rpoispp(int,win=owin(c(0,nx),c(0,ny)))
df<-as.data.frame(rpp)
coordinates(df)<-~x+y 

# assign raster values to sample data
dfpp <-raster::extract(r,df,df=TRUE)
smp<-cbind(coordinates(df),dfpp)
smp<-smp[complete.cases(smp), ]
coordinates(smp)<-~x+y

# fit variogram to sample data
vs <- variogram(sim1~1, data=smp)
m <- fit.variogram(vs, vgm("Exp"))
plot(vs, model = m)

# generate 2 conditional simulations with one core processor
one <- krige(formula = sim1~1, locations = smp, newdata = srgr, model = m,nmax=12,nsim=2)

# plot simulation 1 and 2: statistics (min, max) are ok, simulations are also ok.
spplot(one["sim1"], main = "conditional simulation")
spplot(one["sim2"], main = "conditional simulation")

# generate 2 conditional with parallel processing
no_cores<-detectCores()
cl<-makeCluster(no_cores)
parts <- split(x = 1:length(srgr), f = 1:no_cores)
clusterExport(cl = cl, varlist = c("smp", "srgr", "parts","m"), envir = .GlobalEnv)
clusterEvalQ(cl = cl, expr = c(library('sp'), library('gstat')))
par <- parLapply(cl = cl, X = 1:no_cores, fun = function(x) krige(formula=sim1~1, locations=smp, model=m, newdata=srgr[parts[[x]],],  nmax=12, nsim=2))
stopCluster(cl)

# merge all parts    
mergep <- maptools::spRbind(par[[1]], par[[2]])
mergep <- maptools::spRbind(mergep, par[[3]])
mergep <- maptools::spRbind(mergep, par[[4]])

# create SpatialPixelsDataFrame from mergep
mergep <- SpatialPixelsDataFrame(points = mergep, data = mergep@data)

# plot mergep: statistics (min, max) are ok, but simulated maps show "vertical lines". i don't understand why.
spplot(mergep[1], main = "conditional simulation")
spplot(mergep[2], main = "conditional simulation")

Answer 1

我试过你的代码，我认为问题出在你划分工作的方式上：

parts <- split(x = 1:length(srgr), f = 1:no_cores)

在我的双核机器上，这意味着srgr中的所有奇数索引都由一个进程处理，所有偶数索引都由另一个进程处理。这可能是您所看到的垂直工件的来源。

更好的方法应该是将数据分成连续的块，如下所示：

parts <- parallel::splitIndices(length(srgr), no_cores)

对您的其余代码使用这种拆分，我得到的结果看起来与顺序的结果相当。至少在我未受过训练的眼里。

原始答案，这只是一个很小的影响。使用用于顺序处理的set.seed和用于并行处理的clusterSetRNGStream来修复种子可能仍然是有意义的。

从我读到的关于克里格法的文章来看，它要求你随机抽取数字。这些随机数在并行处理时会有所不同。有关更多详细信息，请参阅parallel vignette (vignette("parallel"))第6节。