R - Как перебрать каждый кусочек 4D матрицы

Question

Я использую doMPI в R для распараллеливания сохранения климатических данных 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 матрица для каждого работника. Это занимает чрезмерно большой объем памяти, и мне нужно уменьшить объем копируемых данных. Из этого ответа я подумал о том, что я мог бы выполнять итерацию по каждому временному срезу матрицы, поэтому только данные для временного среза копировались каждому рабочему на каждой итерации. foreach Конструкция позволяет нескольким объектам выполнять итерацию одновременно, поэтому я мог бы без проблем иметь индекс времени вместе с временным интервалом матрицы. К сожалению, я не знаю ни одного способа перебора матрицы по временному интервалу. Есть ли способ сделать так, чтобы на каждой итерации t из foreach цикл для переменной var Я могу иметь переменную data которая содержит 2-мерную матрицу m[,,t,var]?

Я уже попробовал интуитивный подход ниже, но он повторяется по каждому отдельному элементу, а не по целому временному интервалу за раз.

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

Answer 1

Если вы можете массировать ваши данные в основном процессе R, вы можете попробовать преобразовать каждый 2-мерный срез в big.matrix от bigmemory пакет, и используйте это в ваших параллельных рабочих. Это было бы полезно, только если время, необходимое для обработки каждого среза в подчиненном процессе, является значительным.

Посмотрите этот пример и обратите внимание, что вы можете вкладывать 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или если обработка каждого 2-мерного среза слишком быстрая (да, это может быть проблематично для мультиобработки, см. этот ответ), возможно, вы можете извлечь 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,