Как эффективно найти самую длинную последовательность значений выше порога в R

Question

Я работаю над пространственно-временными наблюдениями температуры, хранящимися в массивах размером 100 * 100 * 504 (сетка 100 * 100, для 504 различных часов, представляющих 21 день). Я вычисляю различные индикаторы из этих наблюдений для разных периодов (от 3 до 21 дня), что, очевидно, требует некоторого времени, и я ищу возможности повышения эффективности вычислений. Я не очень хорошо знаком с R, поэтому я не уверен, что то, что я делаю, является наиболее эффективным способом.

Одна из вещей, которую я хочу сделать, это найти (для каждой ячейки) самый длинный непрерывный период время, когда температура выше определенного порога. Вот что я делаю сейчас:

1. Сначала я вычисляю логический массив на основе порога, используя следующую функцию.

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)
}

Затем я преобразую каждый вектор, представляющий почасовое значение для каждой ячейки, в объект rle и возвращаю максимальную длину последовательностей единиц (представляющих непрерывный период сверх порога).

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(), которая работает, потому что вы знаете, что хотите запускать единицу и проводите немного меньшее сравнение. Это дает вам примерно в 2 раза быстрее, вплоть до среднего времени около 250 миллисекунд или около того на моей машине (общий c macbook).

Если вам придется сделать это 8000 раз, вы сэкономите в большинстве случаев путем распараллеливания кода для работы на многоядерной машине, что легко сделать в R (посмотрите, например, пакет parallel).

Ниже кода для ускорения.

# 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

^{Создано 07.05.2020 с помощью пакета REPEX (v0.3.0)}