A very welcome feature would be the ability to perform the computations of st_apply in parallel. Does anybody have any suggestion on how to do this?
adrfantini
All 24 comments
Are you reading my mind somehow?
edzer
on 2 Oct 2018
I do my best.
adrfantini
on 2 Oct 2018
Integration with package future would be very nice.
adrfantini
on 2 Oct 2018
This is very elementary: if there is a cluster, the (inner) apply loop uses parApply instead of apply.
st_apply has a double loop:
- for all arrays in x (lapply)
- apply function to it (apply)
I now parallelized the inner (apply) one, since the outer one has often length 1, which will not help much. I didn't test even basic functioning, but in case it works I'm looking forward to see crude benchmarks.
edzer
on 2 Oct 2018
I'll do some tests ASAP.
adrfantini
on 2 Oct 2018
Error : object 'getDefaultCluster' is not exported by 'namespace:parallel'
ERROR: lazy loading failed for package 'stars'
adrfantini
on 2 Oct 2018
sessionInfo()?
edzer
on 2 Oct 2018
Before installing latest version:
R version 3.4.4 (2018-03-15)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Linux Mint 18
Matrix products: default
BLAS: /usr/lib/libblas/libblas.so.3.6.0
LAPACK: /usr/lib/lapack/liblapack.so.3.6.0
locale:
[1] C
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] pacman_0.4.6 colorout_1.1-2
loaded via a namespace (and not attached):
[1] compiler_3.4.4
After the (tentative) install.
R version 3.4.4 (2018-03-15)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Linux Mint 18
Matrix products: default
BLAS: /usr/lib/libblas/libblas.so.3.6.0
LAPACK: /usr/lib/lapack/liblapack.so.3.6.0
locale:
[1] C
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] pacman_0.4.6 colorout_1.1-2
loaded via a namespace (and not attached):
[1] httr_1.3.1 compiler_3.4.4 R6_2.2.2 tools_3.4.4 withr_2.1.2 curl_3.2
[7] memoise_1.1.0 knitr_1.20 git2r_0.23.0 digest_0.6.17 devtools_1.13.6
Strange. For now you can use NULL as value for argument CLUSTER, or the output of, e.g.,
cl = makeCluster(4) # use 4 cores
No, sorry I was not clear. The error is during the install process for the latest stars. Not during runtime. Sorry, I thought it was obvious, but clearly it was not. My bad.
adrfantini
on 2 Oct 2018
Could you give the output you see when you do in a fresh R session:
library(parallel)
ls(2)
sessionInfo()
> library(parallel)
> ls(2)
[1] "clusterApply" "clusterApplyLB" "clusterCall" "clusterEvalQ" "clusterExport"
[6] "clusterMap" "clusterSetRNGStream" "clusterSplit" "detectCores" "makeCluster"
[11] "makeForkCluster" "makePSOCKcluster" "mc.reset.stream" "mcMap" "mcaffinity"
[16] "mccollect" "mclapply" "mcmapply" "mcparallel" "nextRNGStream"
[21] "nextRNGSubStream" "parApply" "parCapply" "parLapply" "parLapplyLB"
[26] "parRapply" "parSapply" "parSapplyLB" "pvec" "setDefaultCluster"
[31] "splitIndices" "stopCluster"
> sessionInfo()
R version 3.4.4 (2018-03-15)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Linux Mint 18
Matrix products: default
BLAS: /usr/lib/libblas/libblas.so.3.6.0
LAPACK: /usr/lib/lapack/liblapack.so.3.6.0
locale:
[1] C
attached base packages:
[1] parallel stats graphics grDevices utils datasets methods base
other attached packages:
[1] pacman_0.4.6 colorout_1.1-2
loaded via a namespace (and not attached):
[1] compiler_3.4.4
Here's a test run computing band means from S2 data in starsdata; looks like the overhead is larger than the gain of using 3 cores. In this case.
library(parallel)
cl = makeCluster(3)
library(stars)
# Loading required package: abind
# Loading required package: sf
# Linking to GEOS 3.6.2, GDAL 2.2.3, PROJ 4.9.3
granule = system.file("sentinel/S2A_MSIL1C_20180220T105051_N0206_R051_T32ULE_20180221T134037.zip", package = "starsdata")
base_name = strsplit(basename(granule), ".zip")[[1]]
s2 = paste0("SENTINEL2_L1C:/vsizip/", granule, "/", base_name, ".SAFE/MTD_MSIL1C.xml:10m:EPSG_32632")
(p = read_stars(s2, proxy = TRUE))
# stars_proxy object with 1 attribute in file:
# $`MTD_MSIL1C.xml:10m:EPSG_32632`
# [1] "SENTINEL2_L1C:/vsizip//home/edzer/R/x86_64-pc-linux-gnu-library/3.5/starsdata/sentinel/S2A_MSIL1C_20180220T105051_N0206_R051_T32ULE_20180221T134037.zip/S2A_MSIL1C_20180220T105051_N0206_R051_T32ULE_20180221T134037.SAFE/MTD_MSIL1C.xml:10m:EPSG_32632"
# dimension(s):
# from to offset delta refsys point values
# x 1 10980 3e+05 10 +proj=utm +zone=32 +datum... NA NULL [x]
# y 1 10980 6e+06 -10 +proj=utm +zone=32 +datum... NA NULL [y]
# band 1 4 NA NA NA NA NULL
s2 = st_as_stars(p, downsample = 2)
system.time(n1 <- st_apply(s2, 3, mean, CLUSTER = NULL))
# user system elapsed
# 0.563 0.367 0.930
system.time(n2 <- st_apply(s2, 3, mean, CLUSTER = cl))
# user system elapsed
# 0.848 0.414 1.689
n1
# stars object with 1 dimensions and 1 attribute
# attribute(s):
# MTD_MSIL1C.xml:10m:EPSG_32632
# Min. :2497
# 1st Qu.:2508
# Median :2613
# Mean :2650
# 3rd Qu.:2755
# Max. :2876
# dimension(s):
# from to offset delta refsys point values
# band 1 4 NA NA NA NA NULL
all.equal(n1, n2)
# [1] TRUE
How can I export objects/parameters to be loaded on the nodes?
Using your code, but changing to a different function that requires parameters:
addc_external = 7
f = function(v, addc = addc_external) v + addc
system.time(n2 <- st_apply(s2, 3, f, CLUSTER = cl))
#Error in checkForRemoteErrors(val) :
# 3 nodes produced errors; first error: object 'addc_external' not found
My bad, just pass them explicitly to st_apply: system.time(n2 <- st_apply(s2, 3, f, CLUSTER = cl, addc = addc_external))
adrfantini
on 2 Oct 2018
Works for me, my code went from 36 minutes to 6 minutes using 8 cores. Will test more.
adrfantini
on 2 Oct 2018
st_apply seems to be swapping dimension indexes in some cases, see for example using your data (more downsampled), with 1:2 as margins and a function (rev) that returns a vector:
before applying:
stars object with 3 dimensions and 1 attribute
attribute(s), summary of first 1e+05 cells:
MTD_MSIL1C.xml:10m:EPSG_32632
Min. : 184
1st Qu.: 1742
Median : 2948
Mean : 2955
3rd Qu.: 4004
Max. :13115
dimension(s):
from to offset delta refsys point values
x 1 522 3e+05 210.345 +proj=utm +zone=32 +datum... NA NULL [x]
y 1 522 6e+06 -210.345 +proj=utm +zone=32 +datum... NA NULL [y]
band 1 4 NA NA NA NA NULL
after st_apply(s2, 1:2, rev, CLUSTER = NULL):
stars object with 3 dimensions and 1 attribute
attribute(s), summary of first 1e+05 cells:
MTD_MSIL1C.xml:10m:EPSG_32632
Min. : 31
1st Qu.: 1724
Median : 2707
Mean : 2914
3rd Qu.: 3863
Max. :14030
dimension(s):
from to offset delta refsys point values
rev 1 4 NA NA NA NA NULL
x 1 522 3e+05 210.345 +proj=utm +zone=32 +datum... NA NULL [x]
y 1 522 6e+06 -210.345 +proj=utm +zone=32 +datum... NA NULL [y]
This is not ideal in some cases. If x and y are the first dimensions, they should stay the first dimensions.
adrfantini
on 2 Oct 2018
This happens because apply does this:
> dim(apply(array(1, dim = c(3,4,5)), 1:2, rev))
[1] 5 3 4
Why does dimension order matter?
edzer
on 2 Oct 2018
Consistency. You do not expect the order to change, and you do not want it to change in some cases.
What happened to me in this case is that I read a netCDF file with stars, did my computing, and wanted to write back to the same netCDF file (with ncdf4). It took me a while to figure out I was writing an array with swapped dimensions, which caused of course the resulting data to be scrambled.
aperm fixed this, but it's just not what you would naturally expect.
If the permutation of dimensions is always the same, why not reorder them in st_apply every time? (maybe as an option, since it adds a little bit of computation time / RAM usage, I guess?)
adrfantini
on 2 Oct 2018
I see your point.
Current st_apply is consistent with apply. This saves documentation and managing expectations of those familiar with apply.
Maybe a function that writes stars objects to NetCDF files should pick a sensible default array order?
edzer
on 2 Oct 2018
The CF-Conventions do not enforce any ordering, but suggest:
If any or all of the dimensions of a variable have the interpretations of "date or time" (T), "height or depth" (Z), "latitude" (Y), or "longitude" (X) then we recommend, but do not require (see Section 1.4, "Relationship to the COARDS Conventions" ), those dimensions to appear in the relative order T, then Z, then Y, then X in the CDL definition corresponding to the file. All other dimensions should, whenever possible, be placed to the left of the spatiotemporal dimensions.
http://cfconventions.org/Data/cf-conventions/cf-conventions-1.7/cf-conventions.html#dimensions
adrfantini
on 2 Oct 2018
What do you think about optionally including a progress bar, e.g. with:
https://cran.r-project.org/web/packages/pbapply/pbapply.pdf
?
It is relatively low overhead, and sometimes being able to know if a job is at 5% or 50% can be very helpful.
adrfantini
on 4 Oct 2018
I guess that works only for non-parallel st_apply-applications.
edzer
on 4 Oct 2018
Parallel processing can be enabled through the cl argument. parLapply is called when cl is a
’cluster’ object, mclapply is called when cl is an integer.
adrfantini
on 4 Oct 2018
Related issues
noamross
·
10Comments
courtiol
·
15Comments
michaeldorman
·
10Comments
kendonB
·
9Comments
rgzn
·
7Comments
Most helpful comment
Are you reading my mind somehow?