R:优化尖峰修剪功能

R: Optimise spike pruning function

由于我还没有找到用于分析电生理数据的 R 包,所以我使用了我组的尖峰 p运行ing 函数:

prune.spikes <- function(spikes, min.isi) {
    # copy spike matrix
    prunedspikes <- spikes

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

该函数采用由 0 和 1 值组成的峰值向量或矩阵,如果出现在最小间隔内,则删除任何 1。 由于两个嵌套循环,它需要很长时间才能 运行。为了优化它,我想出了这个解决方案(删除一个循环):

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
    for (i in 1:length(spikes)) {
        if (spikes[i] == 1) {
            if (i - last < min.isi) {
                prunedspikes[i] <- 0;  # remove the spike
            }
            else {
                last <- i
            }
        }
    }
    return(prunedspikes)
}

与原始版本相比,在大型数据集上调用 prune.cols 的速度明显加快(约 60 次)。不过,一个循环仍然存在。到目前为止,我还没有想出一个好的简单的解决方案。如何进一步完善功能?

如果速度差异还不是问题,可能保留循环而不是使用Rcpp会更好

根据 Hadley Wickham 的文章 Loops that should be left as is,使用此循环并不是一个坏主意,因为它可以归类为 递归关系 情况。

一旦速度成为瓶颈,那么求助于 Rcpp 或 this page(文章也建议)可能是解决方案。

像@Khashaa提出的那样,我借助Rcpp实现了功能:

NumericMatrix prunespikes(NumericMatrix spikes, double minisi) {
  NumericMatrix prunedspikes = spikes;
  int ncol = spikes.ncol();
  int nrow = spikes.nrow();
  for (int i = 0; i < ncol; i++) {
    int last = 0;
    while (spikes(last, i) == 0) {
      last++;
    }
    for (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;
}