vectorization - R: Optimise spike pruning function -


since have not found r package analysis of electrophysiological data, have used function spike pruning group:

prune.spikes <- function(spikes, min.isi) {     # copy spike matrix     prunedspikes <- spikes      # initialise index of last spike: infinitely before first one.     (i in 1:ncol(spikes)) {         last <- -inf         (j in 1:nrow(spikes)) {             if (spikes[j, i] == 1) {                 if (j - last < min.isi) {                     prunedspikes[j, i] <- 0;  # remove spike                 }                 else {                     last <- j                 }             }         }        }     return(prunedspikes) } 

the function takes spike vector or matrix consisting of 0 , 1 values , removes 1 if occurred within minimum interval. because of 2 nested loops takes ages run. in order optimise have come solution (removes 1 loop):

prune.cols <- function(spikes, min.isi) {     prunedspikes <- apply(spikes, 2, fun = prune.rows, min.isi = min.isi)     return(prunedspikes) }  prune.rows <- function(spikes, min.isi) {     prunedspikes <- spikes     last <- -inf     (i in 1:length(spikes)) {         if (spikes[i] == 1) {             if (i - last < min.isi) {                 prunedspikes[i] <- 0;  # remove spike             }             else {                 last <-             }         }     }     return(prunedspikes) } 

calling prune.cols on large data set noticeable faster compared original version (~60 times). 1 loop remains, though. far not come nice , simple solution. how can function further improved?

like @khashaa proposed, implemented function of rcpp:

numericmatrix prunespikes(numericmatrix spikes, double minisi) {   numericmatrix prunedspikes = spikes;   int ncol = spikes.ncol();   int nrow = spikes.nrow();   (int = 0; < ncol; i++) {     int last = 0;     while (spikes(last, i) == 0) {       last++;     }     (int j = last + 1; j < nrow; j++) {       if (spikes(j, i) == 1) {         if (j - last < minisi) {           prunedspikes(j, i) = 0;         } else {           last = j;         }       }     }   }   return prunedspikes; } 

Comments