Alternative for loop with "dynamic" variables with R

Question

I am new to StackOverflow, although I played with R. I am struggling with a problem for which I could not find an answer on the site. Please correct me if my searches were not accurate enough.

I have two 3D arrays, in this simplified case 256x256x200. The first is the field, the second is the indices, from 1 to 8. I want to calculate the average value at each vertical level in accordance with the values ​​and calculation of the indices, i.e. An average field for 200 levels for each index (from 1 to 8). This should only be done if there are enough indexes (i.e., the if condition inside the loop). My output should be an 8x200 matrix.

In this example, I create two random arrays. Below is the base code that I use:

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 

, , , , "" , (.. octo ). , , . - ? !

EDIT: 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 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, thr

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 one 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 , , . , !