Segmento vettoriale in base al fatto che i valori superino o meno una soglia in R

Question

Ho un vettore lungo e ho bisogno di dividerlo in segmenti in base a una soglia. Un segmento è un valore consecutivo oltre la soglia. Quando i valori scendono al di sotto della soglia, il segmento termina e il segmento successivo inizia dove i valori ancora una volta superano la soglia. Devo registrare gli indici di inizio e fine di ogni segmento.

Di seguito è una implementazione inefficiente. Qual è il modo più veloce e più appropriato per scrivere questo? Questo è abbastanza brutto, devo supporre che ci sia un'implementazione più pulita.

set.seed(10) 
test.vec <- rnorm(100, 8, 10) 
threshold <- 0 
segments <- list() 
in.segment <- FALSE 
for(i in 1:length(test.vec)){ 

    # If we're in a segment 
    if(in.segment){ 
     if(test.vec[i] > threshold){ 
      next 
     }else{ 
      end.ind <- i - 1 
      in.segment <- FALSE 
      segments[[length(segments) + 1]] <- c(start.ind, end.ind) 
     } 
    } 

    # if not in segment 
    else{ 
     if(test.vec[i] > threshold){   
      start.ind <- i 
      in.segment <- TRUE 
     } 
    } 
} 

EDIT: Runtime di tutte le soluzioni

Grazie per tutte le risposte, questo è stato utile e molto istruttivo. Di seguito è riportato un piccolo test di tutte e cinque le soluzioni (le quattro fornite più l'esempio originale). Come potete vedere, tutti e quattro sono un enorme miglioramento rispetto alla soluzione originale, ma la soluzione di Khashaa è di gran lunga la più veloce.

set.seed(1) 
test.vec <- rnorm(1e6, 8, 10);threshold <- 0 

originalFunction <- function(x, threshold){ 
    segments <- list() 
    in.segment <- FALSE 
    for(i in 1:length(test.vec)){ 

    # If we're in a segment 
     if(in.segment){ 
      if(test.vec[i] > threshold){ 
       next 
      }else{ 
       end.ind <- i - 1 
       in.segment <- FALSE 
       segments[[length(segments) + 1]] <- c(start.ind, end.ind) 
      } 
     } 

    # if not in segment 
     else{ 
      if(test.vec[i] > threshold){   
       start.ind <- i 
       in.segment <- TRUE 
      } 
     } 
    } 
    segments 
} 

SimonG <- function(x, threshold){ 

    hit <- which(x > threshold) 
    n <- length(hit) 

    ind <- which(hit[-1] - hit[-n] > 1) 

    starts <- c(hit[1], hit[ ind+1 ]) 
    ends <- c(hit[ ind ], hit[n]) 

    cbind(starts,ends) 
} 

Rcpp::cppFunction('DataFrame Khashaa(NumericVector x, double threshold) { 
    x.push_back(-1); 
    int n = x.size(), startind, endind; 
    std::vector<int> startinds, endinds; 
    bool insegment = false; 
    for(int i=0; i<n; i++){ 
    if(!insegment){ 
     if(x[i] > threshold){   
     startind = i + 1; 
     insegment = true;   } 
    }else{ 
     if(x[i] < threshold){ 
     endind = i; 
     insegment = false; 
     startinds.push_back(startind); 
     endinds.push_back(endind); 
     } 
    } 
    } 
    return DataFrame::create(_["start"]= startinds, _["end"]= endinds); 
}') 

bgoldst <- function(x, threshold){ 
    with(rle(x>threshold), 
     t(matrix(c(0L,rep(cumsum(lengths),2L)[-length(lengths)]),2L,byrow=T)+1:0)[values,]) 
} 

ClausWilke <- function(x, threshold){ 
    suppressMessages(require(dplyr, quietly = TRUE)) 
    in.segment <- (x > threshold) 
    start <- which(c(FALSE, in.segment) == TRUE & lag(c(FALSE, in.segment) == FALSE)) - 1 
    end <- which(c(in.segment, FALSE) == TRUE & lead(c(in.segment, FALSE) == FALSE)) 
    data.frame(start, end)  
} 

system.time({ originalFunction(test.vec, threshold); }) 
## user system elapsed 
## 66.539 1.232 67.770 
system.time({ SimonG(test.vec, threshold); }) 
## user system elapsed 
## 0.028 0.008 0.036 
system.time({ Khashaa(test.vec, threshold); }) 
## user system elapsed 
## 0.008 0.000 0.008 
system.time({ bgoldst(test.vec, threshold); }) 
## user system elapsed 
## 0.065 0.000 0.065 
system.time({ ClausWilke(test.vec, threshold); }) 
## user system elapsed 
## 0.274 0.012 0.285 

Answer 1

Mi piace for loops per la traduzione in Rcpp è semplice.

Rcpp::cppFunction('DataFrame findSegment(NumericVector x, double threshold) { 
    x.push_back(-1); 
    int n = x.size(), startind, endind; 
    std::vector<int> startinds, endinds; 
    bool insegment = false; 
    for(int i=0; i<n; i++){ 
    if(!insegment){ 
     if(x[i] > threshold){   
     startind = i + 1; 
     insegment = true;   } 
    }else{ 
     if(x[i] < threshold){ 
     endind = i; 
     insegment = false; 
     startinds.push_back(startind); 
     endinds.push_back(endind); 
     } 
    } 
    } 
    return DataFrame::create(_["start"]= startinds, _["end"]= endinds); 
}') 
set.seed(1); test.vec <- rnorm(1e7,8,10); threshold <- 0; 
system.time(findSegment(test.vec, threshold)) 

# user system elapsed 
# 0.045 0.000 0.045 

# @SimonG's solution 
system.time(findSegments(test.vec, threshold)) 
# user system elapsed 
# 0.533 0.012 0.548 

Answer 2

with(rle(test.vec>threshold),t(matrix(c(0L,rep(cumsum(lengths),2L)[-length(lengths)]),2L,byrow=T)+1:0)[values,]); 
##  [,1] [,2] 
## [1,] 1 8 
## [2,] 10 13 
## [3,] 16 17 
## [4,] 20 26 
## [5,] 28 28 
## [6,] 30 34 
## [7,] 36 38 
## [8,] 41 46 
## [9,] 48 49 
## [10,] 51 53 
## [11,] 55 81 
## [12,] 84 90 
## [13,] 92 100 

---

Spiegazione

test.vec>threshold 
## [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE FALSE TRUE TRUE TRUE TRUE FALSE FALSE TRUE TRUE FALSE FALSE TRUE TRUE TRUE TRUE TRUE TRUE TRUE FALSE TRUE FALSE TRUE TRUE TRUE TRUE TRUE FALSE TRUE TRUE TRUE FALSE FALSE TRUE TRUE TRUE TRUE TRUE TRUE FALSE TRUE TRUE FALSE TRUE TRUE TRUE FALSE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE FALSE FALSE TRUE TRUE TRUE TRUE TRUE TRUE TRUE FALSE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE 

Calcola quali elementi del vettore di ingresso sono al di sopra della soglia mediante confronto vettorializzare.

---

rle(...) 
## Run Length Encoding 
## lengths: int [1:25] 8 1 4 2 2 2 7 1 1 1 ... 
## values : logi [1:25] TRUE FALSE TRUE FALSE TRUE FALSE ... 

Calcola il run-length encoding del vettore logico. Ciò restituisce un elenco classificato come 'rle' che contiene due componenti denominati: lengths, contenente le lunghezze di ciascuna lunghezza di esecuzione e values, contenente il valore di quella lunghezza, che in questo caso sarà TRUE o FALSE, con il precedente che rappresenta un segmento di interesse e quest'ultimo rappresenta una durata della corsa non di segmento.

---

with(...,...) 

Il primo argomento è la codifica run-length come descritto sopra. Questo valuterà il secondo argomento in un ambiente virtuale costituito dall'elenco 'rle' -classed, rendendo quindi i componenti lengths e values accessibili come variabili lessicali.

Di seguito mi tuffo nel contenuto del secondo argomento.

---

cumsum(lengths) 
## [1] 8 9 13 15 17 19 26 27 28 29 34 35 38 40 46 47 49 50 53 54 81 83 90 91 100 

Calcola la somma cumulativa del lengths. Ciò costituirà la base per calcolare sia gli indici iniziali che gli indici finali di ogni run-length. Punto critico: ogni elemento della cumsum rappresenta l'indice finale di quel run-length.

---

rep(...,2L) 
## [1] 8 9 13 15 17 19 26 27 28 29 34 35 38 40 46 47 49 50 53 54 81 83 90 91 100 8 9 13 15 17 19 26 27 28 29 34 35 38 40 46 47 49 50 53 54 81 83 90 91 100 

duplicare la somma cumulativa. La prima ripetizione servirà come base per gli indici di partenza, la seconda per la fine. D'ora in poi mi riferirò a queste ripetizioni come "ripetizione dell'indice iniziale" e "ripetizione dell'indice finale".

---

c(0L,...[-length(lengths)]) 
## [1] 0 8 9 13 15 17 19 26 27 28 29 34 35 38 40 46 47 49 50 53 54 81 83 90 91 8 9 13 15 17 19 26 27 28 29 34 35 38 40 46 47 49 50 53 54 81 83 90 91 100 

Ciò rimuove l'ultimo elemento al termine della ripetizione start-index, e antepone uno zero all'inizio di esso. Questo in effetti rallenta la ripetizione dell'indice iniziale di un elemento. Questo è necessario perché abbiamo bisogno di calcolare ogni indice iniziale aggiungendo uno all'indice finale precedente alla lunghezza dell'esecuzione, prendendo come valore zero l'indice finale della run-length inesistente prima del primo.

---

matrix(...,2L,byrow=T) 
##  [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25] 
## [1,] 0 8 9 13 15 17 19 26 27 28 29 34 35 38 40 46 47 49 50 53 54 81 83 90 91 
## [2,] 8 9 13 15 17 19 26 27 28 29 34 35 38 40 46 47 49 50 53 54 81 83 90 91 100 

Questo costruisce una matrice a due ranghi di risultato precedente. La ripetizione ritardata dell'indice iniziale è la riga superiore, la ripetizione dell'indice finale è la riga inferiore.

---

...+1:0 
##  [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25] 
## [1,] 1 9 10 14 16 18 20 27 28 29 30 35 36 39 41 47 48 50 51 54 55 82 84 91 92 
## [2,] 8 9 13 15 17 19 26 27 28 29 34 35 38 40 46 47 49 50 53 54 81 83 90 91 100 

R cicli questi due elementi addendo in righe prima, poi attraverso le colonne, quindi questo aggiunge uno alla riga superiore. Questo completa il calcolo degli indici di partenza.

---

t(...) 
##  [,1] [,2] 
## [1,] 1 8 
## [2,] 9 9 
## [3,] 10 13 
## [4,] 14 15 
## [5,] 16 17 
## [6,] 18 19 
## [7,] 20 26 
## [8,] 27 27 
## [9,] 28 28 
## [10,] 29 29 
## [11,] 30 34 
## [12,] 35 35 
## [13,] 36 38 
## [14,] 39 40 
## [15,] 41 46 
## [16,] 47 47 
## [17,] 48 49 
## [18,] 50 50 
## [19,] 51 53 
## [20,] 54 54 
## [21,] 55 81 
## [22,] 82 83 
## [23,] 84 90 
## [24,] 91 91 
## [25,] 92 100 

Trasposizione di una matrice a due colonne. Questo non è del tutto necessario, se si sta bene ottenendo il risultato come una matrice a due righe.

---

...[values,] 
##  [,1] [,2] 
## [1,] 1 8 
## [2,] 10 13 
## [3,] 16 17 
## [4,] 20 26 
## [5,] 28 28 
## [6,] 30 34 
## [7,] 36 38 
## [8,] 41 46 
## [9,] 48 49 
## [10,] 51 53 
## [11,] 55 81 
## [12,] 84 90 
## [13,] 92 100 

sottoinsieme solo i segmenti di interesse. Poiché values è un vettore logico che rappresenta quali lunghezze di esecuzione hanno superato la soglia, possiamo utilizzarlo direttamente come vettore di indice di riga.

---

prestazioni

Credo di essere l'avvitamento me stesso qui, ma la soluzione di SimonG esegue circa due volte così come la mia:

bgoldst <- function() with(rle(test.vec>threshold),t(matrix(c(0L,rep(cumsum(lengths),2L)[-length(lengths)]),2L,byrow=T)+1:0)[values,]); 
simong <- function() findSegments(test.vec,threshold); 
set.seed(1); test.vec <- rnorm(1e7,8,10); threshold <- 0; 
identical(bgoldst(),unname(simong())); 
## [1] TRUE 
system.time({ bgoldst(); }) 
## user system elapsed 
## 1.344 0.204 1.551 
system.time({ simong(); }) 
## user system elapsed 
## 0.656 0.109 0.762 

+1 da me ...

Answer 3

Ecco un'altra opzione, principalmente utilizzando which. I punti iniziale e finale sono determinati individuando gli elementi non consecutivi della sequenza hit.

test.vec <- rnorm(100, 8, 10) 
threshold <- 0 


findSegments <- function(x, threshold){ 

    hit <- which(x > threshold) 
    n <- length(hit) 

    ind <- which(hit[-1] - hit[-n] > 1) 

    starts <- c(hit[1], hit[ ind+1 ]) 
    ends <- c(hit[ ind ], hit[n]) 

    cbind(starts,ends) 

} 

findSegments(test.vec, threshold=0) 

Questo dà qualcosa di simile:

> findSegments(test.vec, threshold=0) 
     starts ends 
[1,]  1 3 
[2,]  5 7 
[3,]  9 11 
[4,]  13 28 
[5,]  30 30 
[6,]  32 32 
[7,]  34 36 
[8,]  38 39 
[9,]  41 41 
[10,]  43 43 
[11,]  46 51 
[12,]  54 54 
[13,]  56 61 
[14,]  63 67 
[15,]  69 72 
[16,]  76 77 
[17,]  80 81 
[18,]  83 84 
[19,]  86 88 
[20,]  90 92 
[21,]  94 95 
[22,]  97 97 
[23,] 100 100 

Confronta che alla sequenza originale:

> round(test.vec,1) 
    [1] 20.7 15.7 4.3 -15.1 24.6 9.4 23.2 -4.5 16.9 20.9 13.2 -1.2 
[13] 22.6 7.7 6.0 6.6 4.1 21.3 5.3 16.7 11.4 16.7 19.6 16.7 
[25] 11.6 7.3 3.7 8.4 -4.5 11.7 -7.1 8.4 -18.5 12.8 22.5 11.0 
[37] -3.3 11.1 6.9 -7.9 22.9 -3.7 3.5 -7.1 -5.9 3.5 13.2 20.0 
[49] 13.2 23.4 15.9 -5.0 -6.3 10.0 -6.2 4.7 2.1 26.4 5.9 27.3 
[61] 14.3 -12.4 28.4 30.9 18.2 11.4 5.7 -4.5 6.2 12.0 10.9 11.1 
[73] -2.0 -9.0 -1.4 15.4 19.1 -1.6 -5.4 5.4 7.8 -5.6 15.2 13.8 
[85] -18.8 7.1 17.1 9.3 -3.9 22.6 1.7 28.9 -21.3 21.2 8.2 -15.4 
[97] 3.2 -10.2 -6.2 14.1 

Answer 4

Ecco un'altra soluzione che secondo me è più semplice. Si noti che è necessario utilizzare set.seed(10), non set.seed <- 10, per impostare il seme del generatore di numeri casuali.

require(dplyr) # for lead() and lag() 

set.seed(10) 
test.vec <- rnorm(100, 8, 10) 
threshold <- 0 

in.segment <- (test.vec > threshold) 
start <- which(c(FALSE, in.segment) == TRUE & lag(c(FALSE, in.segment) == FALSE)) - 1 
end <- which(c(in.segment, FALSE) == TRUE & lead(c(in.segment, FALSE) == FALSE)) 
segments <- data.frame(start, end) 

head(segments) 
## start end 
## 1  1 2 
## 2  4 6 
## 3  8 8 
## 4 10 16 
## 5 18 21 
## 6 23 23 

In generale, in R, se vi trovate a scrivere loop complicate e se le dichiarazioni il gioco è probabilmente facendo male. La maggior parte dei problemi può essere risolta in una forma vettoriale.