Увеличение скорости векторной переменной переменной длины R

Question

У меня есть следующий код, целью которого является взять один столбец с числовым столбцом данных и создать список st каждые два элемента вектора, ссылаясь на начальный и конечный индекс кадра данных где среднее значение составляет 0,032.

Пример:

Input: [0.012,0.02,0.032,0.045,0.026,0.06,0.01] 
Output [3,5,6,6] 

, как mean(input(3:5))>0.032 и mean(input(6:6))>0.032

Немного более сложный пример ввода [0,0.08,0.08,0.031,0.031, -0,1] Выход [2,5]

Поэтому я не могу просто идентифицировать элементы выше 0,032, и насколько я вижу, мне нужно перебирать каждый индекс. (отсюда цикл while)

Он работает очень хорошо для «небольших кадров данных», но я пытаюсь заставить его работать на кадры данных с 2 000 000 строк, если не больше.

Моя проблема в том, что он работает очень медленно, когда я встаю на большое количество строк. В частности, она стреляет через значения 0-100000, но резко замедляется после

activityduration<-function(input) 
{ 
datum<-as.matrix(input) 
len=length(datum) 
times <-c() 
i<-1 
while (i <len) 
    { 
    if (i>=len) 
    { 
     break 
    } 
    i<-i+1 
     if (datum[i]<0.032) 
     { 
      next 
     } 
     else 
     { 
     vect = c(datum[i]) 
     x<-i 
     while ((mean(vect)>=0.032)){ 
      print(i) 
      if (i==len) 
       { 
       break 
       } 
      i<-i+1 
      boolean <- TRUE 
      vect <- c(datum[x:i]) 
     } 
     if (i==len) 
       { 
       break 
       } 
     if (boolean) 
      { 
      times <- c(times, c(x,i-1)) 
      boolean<-FALSE 
      } 
     } 
    } 
return(times) 
} 

Что я предполагаю, что это вопрос: Я постоянно растущий вектор vect внутри второго цикла While. (по некоторым моим данным vect может достигать длины = 10000). Это означает, что я обновляю vect's размер, повторяюще вызывающий замедление.

Исправления, которые я пробовал: Первоначально вход (кадр данных) был только что получен как кадр данных, я изменил его на матрицу для существенного увеличения скорости.

я заменил еще с:

{ 
newVal = c(datum[i]) 
x<-i 
n<-0 
meanValue<-0 
while (((meanValue*n+newVal)>=(0.032*(n+1))){ 
    print(i) 
    if (i==len) 
     { 
     break 
     } 
    meanValue<-(meanValue*n+newVal)/n+1 
    n<n+1 
    i<-i+1 

    } 

который устранил необходимость вектора при сохранении той же операции, однако это вызывает еще большее замедляться. скорее всего, из-за большого количества выполненных операций.

Я также попытался: Инициирование вектора vect с 700000 элементов, так что никогда не должен расти, но для того, чтобы сделать это мне нужно изменить:

mean(vect)>=0.032 к любому sum(vect)/n >=0.032 или mean(vect[!vect==0]) и этого результата в еще большем замедлении.

Кто-нибудь знает, как я могу увеличить скорость?

Answer 1

Вот еще один алгоритм для получения результатов, идентичных тому, что получает @Joseph Wood.

activityduration <- function(input, th) { 
    epsilon <- 2*.Machine$double.eps 
    a <- input - th 
    s <- 0; l <- 1; r <- 1; f <- F; 
    n <- length(input) 
    res <- vector(mode = "integer", length = 2 * n) 
    j <- 0 
    for (i in 1:n) { 
     s <- s + a[i] 
     if (s < 0 - epsilon) { 
      if (f) { 
       j <- j + 1 
       res[c(2 * j - 1, 2 * j)] <- c(l, r) 
       f <- F 
      } else { 
       l <- i + 1 
      } 
      s <- 0 
     } else { 
      r <- i 
      if (!f) { 
       f <- T 
       l <- i 
      } 
     } 
    } 
    if (f) { 
     j <- j + 1 
     res[c(2 * j - 1, 2 * j)] <- c(l, r) 
    } 
    return(res[res > 0]) 
} 

Испытания на оригинальных примерах

print(activityduration(c(0.012,0.02,0.032,0.045,0.026,0.06,0.01), 0.032)) 
[1] 3 7 
print(activityduration(c(0,0.08,0.08,0.031,0.031,-0.1), 0.032)) 
[1] 2 5 

испытания по данным @ Джозеф Вуд

set.seed(1313) 
options(scipen = 999) 
HighSamp <- sample(51:75, 10, replace = TRUE) 
MidSamp <- sample(36:50, 25, replace = TRUE) 
LowSamp <- sample(11:35, 30, replace = TRUE) 
MinSamp <- sample(1:10, 35, replace = TRUE) 
Samp1 <- sample(c(MinSamp, LowSamp, MidSamp, HighSamp), 20000, replace=TRUE)/1000 
Samp2 <- sample(c(MinSamp, LowSamp, MidSamp, HighSamp), 100000, replace=TRUE)/1000 
Samp3 <- sample(c(MinSamp, LowSamp, MidSamp, HighSamp), 1000000, replace=TRUE)/1000 


JoeTest <- VariableMean(Samp1, 0.032) 
SomeTest <- activityduration(Samp1, 0.032) 

all(JoeTest == SomeTest) 
[1] TRUE 

тесты производительности

library("microbenchmark") 
microbenchmark(Joseph=VariableMean(Samp1, 0.032), SomeAlgo=activityduration(Samp1, 0.032), times = 10) 
Unit: milliseconds 
    expr  min  lq  mean median  uq  max neval 
    Joseph 38.94056 39.54052 40.59358 40.41387 41.83913 42.14377 10 
SomeAlgo 38.14466 38.53188 39.47474 38.91653 40.24965 41.72669 10 
microbenchmark(Joseph=VariableMean(Samp2, 0.032), SomeAlgo=activityduration(Samp2, 0.032), times = 10) 
Unit: milliseconds 
    expr  min  lq  mean median  uq  max neval 
    Joseph 201.9639 212.5006 226.1548 217.6033 238.1169 266.1831 10 
SomeAlgo 194.1691 200.7253 203.0191 203.6269 205.4802 211.1224 10 

system.time(VariableMean(Samp3, 0.032)) 
    user system elapsed 
    2.12 0.01 2.16 
system.time(activityduration(Samp3, 0.032)) 
    user system elapsed 
    2.08 0.02 2.10 

Обсуждения
1. Этот алгоритм имеет коэффициент усиления скорости, хотя и очень умеренный;
2. Ядро алгоритма состоит в том, чтобы избежать прямого вычисления среднего значения, вместо этого оно вычисляет, если совокупная сумма меняет свой знак.

Answer 2

Find_indexes <- function(input,target = 0.032) { 
    l <- length(input) 
    starts <- which(input >= target) 
    seqs <- c(0,diff(starts)) 
    contiguousIdx <- which(seqs == 1) 
    MStarts <- starts[-contiguousIdx] 
    ranges <- vector('list',length(MStarts)) 
    current <- 1 
    it <- 0 
    for (i in MStarts) { 
    fidx <- i 
    i <- i + 1 
    if (i > l) i <- l 
    while (mean(input[fidx:i]) >= target) { 
     i <- i + 1 
     if (i > l) break 
    } 
    ranges[[current]] <- fidx:(i - 1) 
    current <- current + 1 
    } 
    ranges 
} 

set.seed(123) 
qinput <- c(0.012,0.02,0.032,0.045,0.026,0.06,0.01) 
largeinput <- sample(qinput,1e6,replace = TRUE) 

library(microbenchmark) 
microbenchmark(Find_indexes(largeinput,0.032),times=3) 

Идея заключается в том, чтобы ограничить как можно больше петли, поэтому сначала мы ищем запись на входе равной или выше 0,032, рядом с diff и wich мы ищем «смежную» Запись (индекс только +1 или предшественника) и построить вектор только начальных точек.

Далее мы перебираем эти начальные точки и строим список индексов, в то время как среднее от начальной точки до фактической позиции все еще> = до target (0.032 по умолчанию)

Функция возвращает список индексов, если вы хотите только первый и последний indexe вы можете передать результат в lapply и используя функцию ranges

результаты Benchmarck на векторе 1E6:

Unit: seconds 
             expr  min  lq  mean median  uq  max neval 
result <- Find_indexes(largeinput, 0.032) 14.24063 14.25262 14.40224 14.2646 14.48304 14.70147  3 

Это еще 14 секунд на моей машине, но это звучит лучше, чем у вас на самом деле. Извините).

Есть один недостаток, он записывает перекрывающиеся диапазоны.

Выходы:

> head(largeinput,10) 
[1] 0.032 0.060 0.032 0.010 0.010 0.012 0.045 0.010 0.045 0.045 


> head(result) 
[[1]] 
[1] 1 2 3 4 

[[2]] 
[1] 7 

[[3]] 
[1] 9 10 11 12 13 14 

[[4]] 
[1] 12 13 14 

[[5]] 
[1] 19 

[[6]] 
[1] 27 28 29 

> head(lapply(result,range)) 
[[1]] 
[1] 1 4 

[[2]] 
[1] 7 7 

[[3]] 
[1] 9 14 

[[4]] 
[1] 12 14 

[[5]] 
[1] 19 19 

[[6]] 
[1] 27 29 

Answer 3

Try This:

VariableMean <- function(v, Lim) {options(scipen = 999) 
    s <- which(v >= Lim) 
    Len <- length(v) 
    stInd <- c(s[1L], s[which(diff(s) > 1L)+1L]) 
    size <- length(stInd) 
    myIndex <- vector(mode="integer", length = 2*size) 
    bContinue <- FALSE 
    epsilon <- 2*.Machine$double.eps ## added to account for double precision 

    i <- r <- 1L; j <- stInd[i] 
    while (i < size) { 
     k <- stInd[i+1L]-1L 
     temp <- j:k 
     myMeans <- cumsum(v[temp])/(1:length(temp)) 
     myEnd <- temp[max(which(myMeans >= (Lim-epsilon)))] 
     i <- i+1L 
     if (myEnd+1L < stInd[i]) { 
      myIndex[2L*(r-1L)+1L] <- j; myIndex[2L*r] <- myEnd 
      j <- stInd[i] 
      r <- r+1L 
      bContinue <- FALSE 
     } else { 
      bContinue <- TRUE 
     } 
    } 

    if (!bContinue) {j <- stInd[size]} 
    temp <- j:Len 
    mySums <- cumsum(v[temp]) 
    myEnd <- temp[max(which(mySums >= Lim*(1:length(temp))))] 
    myIndex[2L*(r-1L)+1L] <- j; myIndex[2L*r] <- myEnd 

    myIndex[which(myIndex > 0L)] 
} 

VariableMean(c(0.012,0.02,0.032,0.045,0.026,0.06,0.01), 0.032) 
[1] 3 7 

VariableMean(c(0,0.08,0.08,0.031,0.031,-0.1), 0.032) 
[1] 2 5 

Ниже приведены некоторые тесты и тесты (не сравнивал равенство с алгоритмом поставляемого @Tensibai, так как они не делать то же самое (т.е. есть перекрываться @ Tensibai в алгоритме)):

тестовых данных

set.seed(1313) 
options(scipen = 999) 
HighSamp <- sample(51:75, 10, replace = TRUE) 
MidSamp <- sample(36:50, 25, replace = TRUE) 
LowSamp <- sample(11:35, 30, replace = TRUE) 
MinSamp <- sample(1:10, 35, replace = TRUE) 
Samp1 <- sample(c(MinSamp, LowSamp, MidSamp, HighSamp), 20000, replace=TRUE)/1000 
Samp2 <- sample(c(MinSamp, LowSamp, MidSamp, HighSamp), 100000, replace=TRUE)/1000 
Samp3 <- sample(c(MinSamp, LowSamp, MidSamp, HighSamp), 1000000, replace=TRUE)/1000 

Проверки Равенство/Проверка

JoeTest <- VariableMean(Samp1, 0.032) 
OPTest <- activityduration(Samp1) 
FoehnTest <- activityduration2(Samp1, 0.032) 

length(JoeTest) 
[1] 5466 
length(OPTest) 
[1] 5464 

tail(JoeTest) 
[1] 19966 19967 19971 19993 19999 20000 
tail(OPTest)  ## OP's algo doesn't handle the end case 
[1] 19960 19961 19966 19967 19971 19993 

mean(Samp1[19999:20000]) 
[1] 0.065 ## > 0.032 

all(JoeTest[1:length(OPTest)]==OPTest) ## testing equality expect for the end 
[1] TRUE 

all(JoeTest==FoehnTest) 
[1] TRUE 

## Ensuring mean of intervals is greater than 0.032 
TestMean <- sapply(seq.int(1,length(JoeTest),2), function(x) mean(Samp1[JoeTest[x]:JoeTest[x+1L]])) 
all(TestMean >= 0.032) 
[1] TRUE 

Бенчмаркинг

microbenchmark(Joseph=VariableMean(Samp1, 0.032), 
       Foehn=activityduration2(Samp1, 0.032), 
       Tensibai=Find_indexes(Samp1, 0.032), 
       OPAlgo=activityduration(Samp1), times = 10) 
Unit: milliseconds 
    expr  min   lq  mean  median   uq  max neval 
    Joseph 18.191671 19.027055 20.362151 20.917034 21.325900 22.214652 10 
    Foehn 6.848098 7.238491 8.079705 7.829212 9.083315 9.794315 10 
Tensibai 140.924588 142.171712 149.936844 143.188952 148.294031 198.626850 10 
    OPAlgo 122.381933 123.829385 129.934586 128.347027 136.496846 143.782135 10 

microbenchmark(Joseph=VariableMean(Samp2, 0.032), 
       Foehn=activityduration2(Samp2, 0.032), 
       Tensibai=Find_indexes(Samp2, 0.032), 
       OPAlgo=activityduration(Samp2), times = 10) 
Unit: milliseconds 
    expr  min  lq  mean  median   uq  max neval 
    Joseph 95.38979 99.82943 106.67638 101.45689 102.99117 154.21767 10 
    Foehn 36.63334 37.75115 39.00842 38.97406 39.97898 41.26387 10 
Tensibai 709.57490 725.15861 740.39442 737.45620 747.31374 803.22536 10 
    OPAlgo 994.43310 996.61208 1025.54683 1030.84784 1046.03234 1063.52655 10 

system.time(VariableMean(Samp3, 0.032)) 
user system elapsed 
0.98 0.00 1.00 

system.time(activityduration2(Samp3, 0.032)) 
user system elapsed 
0.37 0.00 0.37 

system.time(activityduration(Samp3)) 
user system elapsed 
51.37 0.42 51.82 

system.time(Find_indexes(Samp3, 0.032)) 
user system elapsed 
7.69 0.00 7.72 

На моей машине, алгоритм поставляется @Foehn является самым быстрым и хорошим битом (~ В 3 раза быстрее, чем у меня). @Tensibai, @Foehn и мои алгоритмы кажутся хорошо масштабируемыми и стабильными с большими наборами данных (как видно из диапазона в разы (т. Е. Разница между минимальными и максимальными временами) из тестов).