快速滑动窗口算法的给定坐标问题

Question

我有一个数据表，行数大约在一百万或两百万左右，列数约为200。

每行中的每个条目都有一个相关联的坐标。

以下是数据的一小部分：

[1,] -2.80331471  -0.8874522 -2.34401863   -3.811584   -2.1292443
[2,]  0.03177716   0.2588624  0.82877467    1.955099    0.6321881
[3,] -1.32954665  -0.5433407 -2.19211837   -2.342554   -2.2142461
[4,] -0.60771429  -0.9758734  0.01558774    1.651459   -0.8137684

前4行的坐标：

9928202 9928251 9928288 9928319

我希望您可以提供一个函数，根据给定的数据和窗口大小，在每列上应用平均滑动窗口，并返回相同大小的数据表。或者换句话说 - 对于每个行条目 i ，它将查找介于coords [i] - windsize和coords [i] + windsize之间的坐标，并用该区间内值的平均值替换初始值（对每一列分别执行）。

速度是主要问题。

以下是我编写的这种函数的第一次尝试。

doSlidingWindow <- function(intensities, coords, windsize) {
windHalfSize <- ceiling(windsize/2)
### whole range inds
RANGE <- integer(max(coords)+windsize)
RANGE[coords] <- c(1:length(coords)[1])

### get indeces of rows falling in each window
COORDS <- as.list(coords)
WINDOWINDS <- sapply(COORDS, function(crds){ unique(RANGE[(crds-windHalfSize):
    (crds+windHalfSize)]) })

### do windowing

wind_ints <- intensities
wind_ints[] <- 0
for(i in 1:length(coords)) {
    wind_ints[i,] <- apply(as.matrix(intensities[WINDOWINDS[[i]],]), 2, mean)
}
return(wind_ints)
}

最后一个for循环之前的代码非常快，可以为我获取每个条目所需使用的索引列表。但是，一切都崩溃了，因为我需要磨损for循环一百万次，从我的数据表中取子集，并确保我有多行才能够在apply内同时处理所有列。

我的第二种方法是将实际值粘贴到RANGE列表中，用零填充间隙，并对每个列重复使用zoo包中的rollmean。但这是冗余的，因为rollmean将浏览所有间隙，而我最终只会使用原始坐标的值。

如果不使用C，任何使其更快的帮助都将非常感激。

- Karolis Koncevičius

---

我不是zoo的专家，但你确定使用rollmean(data,fill=NA)不够快吗？ - Carl Witthoft

如果您无论如何都要将数据存储在数据库中：使用PostgreSQL的sqldf可以进行运行窗口统计。 - Dieter Menne

对Carl：rollmean足够快。但是它无法处理任意坐标上的间隔。它只在时间序列上使用固定窗口大小，并且时间序列具有规则间隔。在这种情况下，间隔不规则，两个点之间的间距可以是任意的。因此，如果我为zoo包中的所有间隙填充零，我将得到长度约为5亿的向量。在DataFrame上使用rollmean很痛苦，尤其是当我只需要从那500个计算出的数百万个中获取几个时。 - Karolis Koncevičius

在最后一个循环中，最好将代码更改为：wind_ints[i,] <- apply(matrix(intensities[WINDOWINDS[[i]],], ncol=ncol(intensities)), 2, mean)。当窗口中只有一行时，您的代码会导致错误的结果。 - redmode

Answer 1

数据生成：

N <- 1e5 # rows
M <- 200 # columns
W <- 10  # window size

set.seed(1)
intensities <- matrix(rnorm(N*M), nrow=N, ncol=M)
coords <- 8000000 + sort(sample(1:(5*N), N))

我用于基准测试的带有轻微修改的原始函数：

doSlidingWindow <- function(intensities, coords, windsize) {
  windHalfSize <- ceiling(windsize/2)
  ### whole range inds
  RANGE <- integer(max(coords)+windsize)
  RANGE[coords] <- c(1:length(coords)[1])

  ### get indices of rows falling in each window
  ### NOTE: Each elements of WINDOWINDS holds zero. Not a big problem though.
  WINDOWINDS <- sapply(coords, function(crds) ret <- unique(RANGE[(crds-windHalfSize):(crds+windHalfSize)]))

  ### do windowing
  wind_ints <- intensities
  wind_ints[] <- 0
  for(i in 1:length(coords)) {
    # CORRECTION: When it's only one row in window there was a trouble
    wind_ints[i,] <- apply(matrix(intensities[WINDOWINDS[[i]],], ncol=ncol(intensities)), 2, mean)
  }
  return(wind_ints)
}

---

可能的解决方案：

---

1) data.table

data.table 以子集为特点，速度较快，但是这个页面（和其他与滑动窗口相关的页面）表明，事实并非如此。实际上，data.table 的代码很优雅，但是不幸的是非常缓慢：

require(data.table)
require(plyr)
dt <- data.table(coords, intensities)
setkey(dt, coords)
aaply(1:N, 1, function(i) dt[WINDOWINDS[[i]], sapply(.SD,mean), .SDcols=2:(M+1)])

---

2) foreach+doSNOW

基本例程很容易并行运行，因此我们可以从中受益：

require(doSNOW)
doSlidingWindow2 <- function(intensities, coords, windsize) {
  NC <- 2 # number of nodes in cluster
  cl <- makeCluster(rep("localhost", NC), type="SOCK")
  registerDoSNOW(cl)

  N <- ncol(intensities) # total number of columns
  chunk <- ceiling(N/NC) # number of columns send to the single node

  result <- foreach(i=1:NC, .combine=cbind, .export=c("doSlidingWindow")) %dopar% {
    start <- (i-1)*chunk+1
    end   <- ifelse(i!=NC, i*chunk, N)
    doSlidingWindow(intensities[,start:end], coords, windsize)    
  }

  stopCluster(cl)
  return (result)
}

基准测试显示我的双核处理器速度显著提升：

system.time(res <- doSlidingWindow(intensities, coords, W))
#    user  system elapsed 
# 306.259   0.204 307.770
system.time(res2 <- doSlidingWindow2(intensities, coords, W))
#  user  system elapsed 
# 1.377   1.364 177.223
all.equal(res, res2, check.attributes=FALSE)
# [1] TRUE

---

3) Rcpp

是的，我知道您想要“不用C语言”，但请看一下这个。这段代码是内联的，并且相当简单：

require(Rcpp)
require(inline)
doSlidingWindow3 <- cxxfunction(signature(intens="matrix", crds="numeric", wsize="numeric"), plugin="Rcpp", body='
  #include <vector>
  Rcpp::NumericMatrix intensities(intens);
  const int N = intensities.nrow();
  const int M = intensities.ncol();
  Rcpp::NumericMatrix wind_ints(N, M);

  std::vector<int> coords = as< std::vector<int> >(crds);
  int windsize = ceil(as<double>(wsize)/2);  

  for(int i=0; i<N; i++){
    // Simple search for window range (begin:end in coords)
    // Assumed that coords are non-decreasing
    int begin = (i-windsize)<0?0:(i-windsize);
    while(coords[begin]<(coords[i]-windsize)) ++begin;
    int end = (i+windsize)>(N-1)?(N-1):(i+windsize);
    while(coords[end]>(coords[i]+windsize)) --end;

    for(int j=0; j<M; j++){
      double result = 0.0;
      for(int k=begin; k<=end; k++){
        result += intensities(k,j);
      }
      wind_ints(i,j) = result/(end-begin+1);
    }
  }

  return wind_ints;
')

基准测试：

system.time(res <- doSlidingWindow(intensities, coords, W))
#    user  system elapsed 
# 306.259   0.204 307.770
system.time(res3 <- doSlidingWindow3(intensities, coords, W))
#  user  system elapsed 
# 0.328   0.020   0.351
all.equal(res, res3, check.attributes=FALSE)
# [1] TRUE

我希望这个结果能够激励你。当数据能够适应内存时，Rcpp 版本的速度非常快。例如：N <- 1e6 和 M <-100，我获得了以下结果：

   user  system elapsed 
  2.873   0.076   2.951

当R开始使用交换空间时，所有内容都会变慢。对于无法放入内存的大型数据，您应该考虑使用sqldf、ff或bigmemory。

- redmode

---

没问题，这些评论已经涵盖了它。在网上也有不良使用 data.table 的例子是好的。 - Matt Dowle

@redmode 这里需要多少调整才能处理从0开始的coords？也就是说，如果我执行coords <- sort(sample(1:(5*N), N))，代码会抛出错误。 - Jota

@redmode，我最初使用了doSlidingWindow。现在，我尝试了doSlidingWindow3，它可以工作，但值得注意的是，我需要输入我所需窗口大小的1/2而不是窗口大小本身，并且在定义函数后我收到了一个cygwin警告消息。 - Jota

@Frank，是什么样的消息？虽然我无法复现它。 - redmode

@Frank，这只是CYGWIN的警告，你可以在这里学习如何关闭它：https://cygwin.com/cygwin-ug-net/setup-env.html。这个答案也可能有帮助：https://dev59.com/QWHVa4cB1Zd3GeqPpLBu。 - redmode

Answer 2

Rollapply在小数据集上表现很好。但是，如果您正在处理数百万行（基因组学），它会变得非常缓慢。

以下函数速度超快：

data <- c(runif(100000, min=0, max=.1),runif(100000, min=.05, max=.1),runif(10000, min=.05, max=1), runif(100000, min=0, max=.2))
slideFunct <- function(data, window, step){
  total <- length(data)
  spots <- seq(from=1, to=(total-window), by=step)
  result <- vector(length = length(spots))
  for(i in 1:length(spots)){
    result[i] <- mean(data[spots[i]:(spots[i]+window)])
  }
  return(result)
}

这里详细信息。

- Jacky Smith