如何在R中高效地找到阈值以上的最长值序列

Question

我正在进行温度的时空观测，存储在大小为100*100*504（100*100网格，代表21天的504个不同小时）的阵列中。我正在根据这些观察结果计算不同时期（3到21天）的各种指标，这显然需要一些时间，我正在考虑提高计算效率。我不太习惯R，所以我不确定我所做的是不是最有效的方法。

我想做的事情之一是找到(每个细胞)温度高于某个阈值的最长连续时间。这就是我此刻正在做的事情:

< li >首先，我使用以下函数根据阈值计算一个布尔数组。

utci_test = array(runif(100*100*504, min = 18, max = 42), c(100,100,504))
to_hs = function(utci, period=1:length(utci[1,1,]), hs_threshold){
  utci_hs = utci*0
  utci_hs[which(utci > hs_threshold)] = 1
  utci_hs[is.na(utci)] = 0
  return(utci_hs)
}

max_duration_hs = function(utci_hs, period=1:length(utci_hs[1,1,]) ){
  apply(utci_hs, MARGIN=c(1,2), FUN=function(x){
    r = rle(x)
    max(r$lengths[as.logical(r$values)], fill = 0)
  })
}

考虑到所需的时间，我注意到第二步需要一些时间（请记住，我必须重复此操作总共8000次）

system.time(to_hs(utci_test, hs_threshold=32.0))
# utilisateur     système      écoulé 
#      0.051       0.004       0.055 
system.time(to_hs(utci_test, hs_threshold=32.0))
# utilisateur     système      écoulé 
#      0.053       0.000       0.052 
utci_test_sh = to_hs(utci_test, hs_threshold=32.0)
system.time(max_duration_hs(utci_test_sh))
# utilisateur     système      écoulé 
#      0.456       0.012       0.468 

所以，我想知道是否有一种更有效的方法来实现这一点，因为我认为转换为rle对象可能效率低下？

Answer 1

你可以通过编写你自己版本的rle（）函数来获得一点减速带，因为你知道你想要运行1，并且做的比较少一点。这让你快了大约2倍，在我的机器（通用macbook）上下降到大约250毫秒左右的中值时间。

如果你必须这样做 8，000 次，你将通过并行化代码以在多核机器上运行来节省自己的大部分时间，这在 R 中很简单（例如查看并行包）。

在加速的代码下方。

# generate data
set.seed(123)
utci_test <- array(runif(100*100*504, min = 18, max = 42), c(100,100,504))

# original functions
to_hs = function(utci, period=1:length(utci[1,1,]), hs_threshold){
  utci_hs = utci*0
  utci_hs[which(utci > hs_threshold)] = 1
  utci_hs[is.na(utci)] = 0
  return(utci_hs)
}

max_duration_hs = function(utci_hs, period=1:length(utci_hs[1,1,]) ){
  apply(utci_hs, MARGIN=c(1,2), FUN=function(x){
    r = rle(x)
    max(r$lengths[as.logical(r$values)], fill = 0)
  })
}

# helper func for rle
rle_max <- function(v) {
  max(diff(c(0L, which(v==0), length(v)+1))) - 1
}

max_dur_hs_2 <- function(utci_hs) {
  apply(utci_hs, MARGIN=c(1,2), FUN= rle_max)
 }

# Check equivalence
utci_hs <- to_hs(utci = utci_test, hs_threshold = 32)

all.equal(max_dur_hs_2(utci_hs), 
          max_duration_hs(utci_hs))
#> [1] TRUE

# Test speed
library(microbenchmark)

microbenchmark(max_dur_hs_2(utci_hs), 
               max_duration_hs(utci_hs))
#> Unit: milliseconds
#>                      expr      min       lq     mean   median       uq      max
#>     max_dur_hs_2(utci_hs) 216.1481 236.7825 250.9277 247.9918 262.4369 296.0146
#>  max_duration_hs(utci_hs) 454.5740 476.5710 501.5119 489.9536 509.8750 774.9963
#>  neval cld
#>    100  a 
#>    100   b

^{创建于 2020-05-07 由 reprex 软件包 （v0.3.0）}