Альтернатива для цикла for с "динамическими" переменными с R

Question

Я новичок в Stackru, хотя уже давно играю с R. Я борюсь с проблемой, из-за которой мне не удалось найти ответ на сайте. Пожалуйста, поправьте меня, если мой квест был недостаточно точным.

У меня есть два 3D-массива, в этом упрощенном случае 256x256x200. Первый из них представляет собой поле, второй состоит из индексов, охватывающих от 1 до 8. Я хочу вычислить среднее значение для каждого вертикального уровня в соответствии со значениями и количеством индексов, то есть среднее значение для поля для 200 уровней. для каждого индекса (от 1 до 8). Это должно быть сделано только в том случае, если имеется достаточное количество индексов (то есть условие if в цикле). Мой вывод должен быть матрицей 8х200.

Для примера я создаю два случайных массива. Здесь ниже есть основной код, который я использую:

nz=200
lev=1:nz
indices=8
var0=array(rnorm(256*256*nz),dim=c(256,256,nz))
#octo=array(sample(1:indices),dim=c(256,256,nz)) 
octo=array(sample(1:indices,size=256*256*nz,replace=T),dim=c(256,256,nz))
counts=apply(octo,3,function(x) table(factor(x,levels=1:indices)))
#thr=0.1
thr=0.125
np=length(var0[,1,1])*length(var0[1,,1])
profile=array(NA,dim=c(nz,indices))


t0=proc.time()
for (i in 1:indices)
{
    for (z in 1:length(lev)) 
    {
       if (counts[i,z]/np>thr) 
       {v0=var0[,,z];  profile[z,i]=counts[i,z]/np*mean(v0[octo[,,z]==i],na.rm=T)} 
    }
}
print(proc.time()-t0)

user  system elapsed 
5.169   0.001   5.170 

Я попытался применить семейство функций apply, но я не могу записать его разумным и эффективным способом, учитывая, что мне нужно, чтобы каждое вычисление учитывало "динамическую" переменную, которая изменяет свой уровень (то есть octo и count vars). Мой реальный случай сделан с помощью больших матриц, и это должно быть сделано на десятках полей, поэтому время довольно актуально. Вам известны какие-либо более быстрые альтернативы? Большое спасибо за любую помощь!

РЕДАКТИРОВАТЬ: я исправил первоначальное определение octo и я скорректировал порог thr. Таким образом, условие if имеет смысл, поскольку оно не всегда соблюдается.

Answer 1

Вот data.table изменить решение, которое избегает циклов и / или применять операторы:

nz=200
lev=1:nz
indices=8
var0=array(rnorm(256*256*nz),dim=c(256,256,nz))
octo=array(sample(1:indices),dim=c(256,256,nz))
counts=apply(octo,3,function(x) table(factor(x,levels=1:indices)))
thr=0.1
np=length(var0[,1,1])*length(var0[1,,1])
profile=array(NA,dim=c(nz,indices))


# From here load data.table to do the manipulation
# reshape2 to convert back into a matrix at the end
library(data.table)
library(reshape2)

# Take the data long and convert to data.table
var01 <- setDT(melt(var0))
octo1 <- setDT(melt(octo))

# Join the data to get corresponding data
# EDIT, it currently works, but I think that's because all data is defined
# adding nomatch in case of missing data
octo1 <- octo1[var01, on = c('Var1','Var2','Var3'), nomatch = NA] 

# Make our calculation grouping by the vertical dimension and the value
profile <- octo1[,if(.N/np > thr) .N / np * mean(i.value, na.rm = TRUE) else NA, by = .(value,Var3)]

# Recast to matrix
profile <- acast(profile, value ~ Var3, mean, value.var = 'V1')

Answer 2

Похоже, это быстрее на моей машине:

profile2 <- sapply(lev, function(i){
    v0 <- var0[,,i]
    mV <- sapply(1:indices, function(j){
        mean(v0[octo[,,i] == j], na.rm = TRUE)
    })
    counts[,i]/np*mV
})
profile2[counts/np > thr] <- NA
profile2<- t(profile2)

all.equal(profile, profile2)
## TRUE

Я пытался сравнить их с microbenchmark пакет, но это занимает довольно много времени... Вот быстрое сравнение, которое я сделал с rbenchmark пакет

f1 <- function(){
    for (i in 1:indices){
        for (z in 1:length(lev)) {
            if (counts[i,z]/np>thr){
                v0=var0[,,z];  profile[z,i]=counts[i,z]/np*mean(v0[octo[,,z]==i],na.rm=T)
            } 
        }
    }
}

f2 <- function(){
    prof <- sapply(lev, function(i){
        v0 <- var0[,,i]
        mV <- sapply(1:indices, function(j){
            mean(v0[octo[,,i] == j], na.rm = TRUE)
        })
        counts[,i]/np*mV
    })
    profile2[counts/np > thr] <- NA
    profile2<- t(profile2)
}

library(rbenchmark)
benchmark(f1(), f2(), replications = 10)

Я поместил оба кода в функцию и проверил. Вот результат:

##   test replications elapsed relative user.self sys.self
## 1 f1()           10   89.03    1.342     85.15     1.72
## 2 f2()           10   66.34    1.000     61.50     0.75

Answer 3

Я думаю, что я нахожу хорошее решение с sapply в том числе ч

f1<-function()
{   
for (i in 1:indices)
{
for (z in 1:length(lev)) {if (counts[i,z]/np>thr) {v0=var0[,,z]; profile[z,i]=counts[i,z]/np*mean(v0[octo[,,z]==i],na.rm=T) } }
}
return(profile)
}

f2<-function()
{
profile=sapply(lev, function(i) {
            v0=var0[,,i];
            mV=sapply(1:indices, function(j) {mean(v0[octo[,,i] == j], na.rm = TRUE)})
            counts[,i]/np*mV
    })

profile[counts/np <= thr]=NA
profile<-matrix(profile, nz, indices, byrow = TRUE)
return(profile)
}

f3<-function()
{
profile=sapply(lev, function(i) {
            v0=var0[,,i];
            mV=sapply(1:indices, function(j) {if (counts[j,i]/np>thr) {mean(v0[octo[,,i] == j], na.rm = TRUE)} else {NA}})
            counts[,i]/np*mV
    })

profile<-matrix(profile, nz, indices, byrow = TRUE)
return(profile)
}

На самом деле f1() - оригинал, f2() - @parksw3, а f3() - моя версия немного улучшена.

benchmark(f1(),f2(),f3(),replications=10)

test   replications elapsed relative user.self sys.self user.child  sys.child
1 f1()           10  27.382    1.411    27.375        0          0         0
2 f2()           10  35.195    1.814    35.186        0          0         0
3 f3()           10  19.403    1.000    19.392        0          0         0

Таким образом, это всегда быстрее, чем стандартный цикл. data.table скорее всего быстрее, но требует полного изменения структуры данных, которую я пока не могу выполнить. Надеюсь это поможет!