数据生成:
N <- 1e5
M <- 200
W <- 10
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)
RANGE <- integer(max(coords)+windsize)
RANGE[coords] <- c(1:length(coords)[1])
WINDOWINDS <- sapply(coords, function(crds) ret <- unique(RANGE[(crds-windHalfSize):(crds+windHalfSize)]))
wind_ints <- intensities
wind_ints[] <- 0
for(i in 1:length(coords)) {
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
cl <- makeCluster(rep("localhost", NC), type="SOCK")
registerDoSNOW(cl)
N <- ncol(intensities)
chunk <- ceiling(N/NC)
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))
system.time(res2 <- doSlidingWindow2(intensities, coords, W))
all.equal(res, res2, check.attributes=FALSE)
3) Rcpp
是的,我知道您想要“不用C语言”,但请看一下这个。这段代码是内联的,并且相当简单:
require(Rcpp)
require(inline)
doSlidingWindow3 <- cxxfunction(signature(intens="matrix", crds="numeric", wsize="numeric"), plugin="Rcpp", body='
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))
system.time(res3 <- doSlidingWindow3(intensities, coords, W))
all.equal(res, res3, check.attributes=FALSE)
我希望这个结果能够激励你。当数据能够适应内存时,Rcpp
版本的速度非常快。例如:N <- 1e6
和 M <-100
,我获得了以下结果:
user system elapsed
2.873 0.076 2.951
当R开始使用交换空间时,所有内容都会变慢。对于无法放入内存的大型数据,您应该考虑使用sqldf
、ff
或bigmemory
。
zoo
的专家,但你确定使用rollmean(data,fill=NA)
不够快吗? - Carl Witthoftwind_ints[i,] <- apply(matrix(intensities[WINDOWINDS[[i]],], ncol=ncol(intensities)), 2, mean)
。当窗口中只有一行时,您的代码会导致错误的结果。 - redmode