Benchmark for loading HDF5 arrays into memory

# Define test file
#
suppressPackageStartupMessages({
    library(HDF5Array)
    library(data.table)
})
d <- '/s/project/patient_report/tidy_results/prokisch-full/savedObjects/Full_Prokisch_Analysis'
f <- 'zscore_psiSite.h5'
fp <- file.path(d,f)
h5ls(fp)

##   group           name       otype dclass           dim
## 0     / zscore_psiSite H5I_DATASET  FLOAT 3269955 x 146

da <- HDF5Array(fp, sub('\\.h5','',f))

Subsetting by rows

loading does not scale linearly

log10(nrow(da))

## [1] 6.514542

system.time(as.matrix(da[1:1e3,]))

##    user  system elapsed 
##   0.145   0.005   0.151

system.time(as.matrix(da[1:1e4,]))

##    user  system elapsed 
##   0.669   0.006   0.676

system.time(as.matrix(da[1:1e5,]))

##    user  system elapsed 
##   1.514   0.059   1.572

random access is slower than by blocks

system.time(as.matrix(da[sample(1:nrow(da), 1e3),]))

##    user  system elapsed 
##  14.531   0.284  14.815

Subsetting by columns

ncol(da)

## [1] 146

system.time(as.matrix(da[,1]))

##    user  system elapsed 
##   1.270   0.020   1.291

system.time(as.matrix(da[,1:10]))

##    user  system elapsed 
##   3.021   0.097   3.119

system.time(as.matrix(da[,1:25]))

##    user  system elapsed 
##   5.816   0.479   6.295

system.time(as.matrix(da[,1:50]))

##    user  system elapsed 
##   9.761   1.421  11.184

random access is slower than by blocks

system.time(as.matrix(da[,sample(1:ncol(da), 10)]))

##    user  system elapsed 
##   4.817   0.255   5.072

system.time(as.matrix(da[,sample(1:ncol(da), 25)]))

##    user  system elapsed 
##  11.094   0.539  11.632

Use case I: Particular rows for 1 column

• faster loading time with less but random rows?

system.time(as.matrix(da[,1]))

##    user  system elapsed 
##   1.283   0.037   1.320

system.time(as.matrix(da[sample(1:nrow(da), 1e3),1]))

##    user  system elapsed 
##   0.181   0.003   0.184

system.time(as.matrix(da[sample(1:nrow(da), 1e4),1]))

##    user  system elapsed 
##   7.999   0.002   8.009

If the row subset is faster, you could spend more time in a clever row selection. For us this could be no missing values, no Pvalues>0.5, …, to reduce the number of points we plot in the end. Typical size would be 5k to 10k points to plot.

Use case II: Particular rows for multiple columns

Extract multiple features (organized in different HDF5 files) for 1 patient (1 column per file), but only the relevant rows (random indices).

• faster loading time with less but random rows?

system.time(as.matrix(da[,1:10]))

##    user  system elapsed 
##    2.35    0.16    2.51

system.time(as.matrix(da[sample(1:nrow(da), 1e3), 1:10]))

##    user  system elapsed 
##   1.267   0.020   1.287

system.time(as.matrix(da[sample(1:nrow(da), 1e3), sample(1:ncol(da), 10)]))

##    user  system elapsed 
##   1.850   0.008   1.859

system.time(as.matrix(da[sample(1:nrow(da), 1e4), 1:10]))

##    user  system elapsed 
##  14.308   0.019  14.327

# very slow ~2min
# system.time(as.matrix(da[sample(1:nrow(da), 1e4), sample(1:ncol(da), 10)]))