R - 如何迭代 4D 矩阵的每个切片

Question

我在 R 中使用 doMPI 来并行保存 netCDF 气候数据。该数据存储在 R 中的 4 维矩阵 m 中，包含 6 个变量的数据，在 20000 个时间点，在一个纬度和经度网格上。 m 因此被索引为m[lon,lat,time,variable]。根据 netCDF 将其数据存储在磁盘上的方式，将数据写入磁盘的最有效方式是按时间片。因此，我想为每个变量一次迭代m 一个时间片。目前，我的代码如下所示：

    ntime <- 20000
    output.vars <- list("rainfall", "snowfallwateq", "snowmelt", "newsnow", "snowdepth", "swe") 
    for (var.index in seq_along(output.vars)) {
        ncout <- nc_open(output.files[var.index], write=TRUE)

        val <- foreach(time.index=1:ntime, .packages=c("ncdf4")) %dopar%
        {
            ncvar_put(ncout, output.vars[[var.index]], 
                      vals=m[,,time.index,var.index],
                      start=c(1, 1, time.index),
                      count=c(nlon, nlat, 1))
        }

        nc_close(ncout)
    }

这不必要地将整个m 矩阵复制到每个工作人员。这占用了大量的内存，我需要减少复制的数据量。我对from this answer 的想法是我可以迭代矩阵的每个时间片，因此每次迭代时只有时间片的数据被复制到每个工作人员。 foreach 构造允许多个对象同时迭代，因此我什至可以将时间索引与矩阵时间片放在一起而不会出现问题。不幸的是，我不知道有任何方法可以按时间片迭代矩阵。有没有办法做到这一点，这样在foreach 循环变量var 的每次迭代t 中，我可以有一个变量data 保存二维矩阵m[,,t,var]？

我已经尝试过下面的直观方法，但它会遍历每个单独的元素，而不是一次遍历整个时间片。

val <- foreach(time.index=1:ntime, slice=m[,,,var], ...

Answer 1

如果您可以在 R 主进程中处理您的数据， 您可以尝试将每个二维切片从bigmemory 包中转换为big.matrix， 并在您的并行工作人员中使用它。 这只有在从属进程中处理每个切片所需的时间很长时才有用。

看这个例子，注意你可以用%:%嵌套2个foreach循环

m <- as.numeric(1:16)
dim(m) <- rep(2L, 4L)

# use %do% for sequential processing, without copying the data to parallel workers
big_m <- foreach(i=1L:2L, .combine=c) %:% foreach(j=1L:2L, .combine=list) %do% {
  as.big.matrix(m[,,i,j], type="double")
}

descriptors <- lapply(big_m, describe)

# specify .noexport to avoid copying the data to each worker
foreach(m_slice_desc=descriptors, .packages=c("bigmemory"), .noexport=ls(all.names=TRUE)) %dopar% {
  # you could even modify the slices in parallel if you wanted
  m_slice <- attach.big.matrix(m_slice_desc)
  for (i in 1L:2L) {
    for (j in 1L:2L) {
      m_slice[i,j] <- m_slice[i,j] * 2
    }
  }
  # return nothing
  NULL
}

# just to show that the values were modified in place
for (bm in big_m) { print(bm[,]) }
     [,1] [,2]
[1,]    2    6
[2,]    4    8
     [,1] [,2]
[1,]   18   22
[2,]   20   24
     [,1] [,2]
[1,]   10   14
[2,]   12   16
     [,1] [,2]
[1,]   26   30
[2,]   28   32

如果你不能/不会使用bigmemory， 或者如果处理每个二维切片太快 （是的，这对于多处理可能有问题， 见this answer), 也许您可以从数据中提取 3 维切片并使用 .noexport 一次只复制一个， 类似：

slices_3d <- lapply(1L:2L, function(i) { m[,,,i] })
foreach(slice_3d=slices_3d, .noexport=ls(all.names=TRUE)) %dopar% {
  for (j in 1L:2L) {
    slice_2d <- slice_3d[,,j]
    # do something
  }
  # return nothing
  NULL
}

我实际上不是 100% 确定以上内容会阻止复制整个 slices_3d， 如果没有，您可能必须在主 R 进程中手动提取块中的子集 （例如slices_3d[1L:num_parallel_workers] 等每次）， 并确保每次调用foreach时只导出一个块。