Sliding window calculation using R

Question

I have the following function count.ones that takes in a data frame data.l in which the observations have been subdivided into bins defined by indices, a vector of consecutive indices x, and the length of the sliding window width over consecutive indices. The goal of this function is to slide over the indices and count the number of observations labeled as ones.

count.ones = function(x,data.l,width){
  labvec = data.l$lab
  indvec = data.l$index
  return(sapply(1:(length(x)-(width-1)),
                       function(j) sum(labvec[indvec%in%x[j:(j+(width-1))]])))
}

My data is as follows:

library(ks)
library(zoo)
library(plyr)
library(BBmisc)
library(ggplot2)
library(rje)

set.seed(1234)
n1 = n2 = 50000
X1 = rnorm.mixt(n=n1,mus=c(0.4,0.8),sigmas=c(0.04,0.04),props=c(0.97,0.03))
X2 = rnorm.mixt(n=n2,mus=c(0.4,0.8),sigmas=c(0.04,0.04),props=c(0.98,0.02))

n1 = length(X1)
n2 = length(X2)

#pooled observations - label X1's as 1 and X2's as 0
dat = data.frame(X=c(X1,X2),lab=c(rep(1,n1),rep(0,n2)))

N = n1+n2 

#number in each bin
n = N^(1/3)

dat$quant = with(dat, cut_number(dat$X, n = round(N/n)))

#number of bins
M = length(levels(dat$quant))
dat$index = factor(dat$quant, levels = levels(dat$quant), 
                   labels = 1:length(levels(dat$quant)))

Example function call:

indx.contig = seq(M)
count.ones(x=indx.contig,data.l=dat,width=1) #can vary width for any integer >= 1

This function seems to be slower than expected as it's only performing the simple task of counting the number of 1's in each sliding window. This would become problematic as n1 and n2 increase into the millions (which is ultimately what I want to use). Is there a way to speed this up? I initially looked into rollapply in the zoo package, but it is not as flexible as my code.

It is important for me to look at consecutive indices. If, for example, I had:

indx.vec = c(1:10,14:20)

I would first create a list with two elements: one containing 1:10 and the other containing 14:20. I would then run count.ones on each of the elements using sapply.

Answer 1

Your code is inefficient because each iteration within the sapply has to search through the 100k rows of data to find the relevant ones. Without changing your code at all, you could get a significant speed improvement by first pre-aggregating your input data:

dat <- aggregate(lab ~ index, data = dat, FUN = sum)

This way, each iteration will only have to search through $M = 2154$ rows instead of the 100k. But it is not fixing everything and I prefer rewriting your code with your more general goal in mind, where

indx.vec <- c(1:10, 14:20)

First, aggregate your data using tapply:

agg.lab <- with(dat, tapply(lab, list(index), FUN = sum))
length(agg.lab)
# [1] 2154
head(agg.lab)
#  1  2  3  4  5  6 
# 25 23 26 21 22 21

It returns a vector containing sum(lab) for each distinct value of index. In the event that the data has no value for a given index, you could overwrite the output-ed NA with a 0 so as to not affect the upcoming roll sums:

agg.lab[is.na(agg.lab)] <- 0

Next, we subset by only keeping the indices in indx.vec

agg.lab   <- agg.lab[indx.vec]
#  1  2  3  4  5  6  7  8  9 10 14 15 16 17 18 19 20 
# 25 23 26 21 22 21 19 21 24 28 21 18 24 22 23 21 20 

Next, we map each value of indx.vec to the "group number" it belongs to (here, 1:14 belong to group #1 then 14:20 belong to group #2):

indx.grp  <- cumsum(c(TRUE, diff(indx.vec) != 1))
# [1] 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2

Using this vector of group indices, we break agg.lab into a list of vectors, one for each group number:

agg.split <- split(agg.lab, indx.grp)
# $`1`
#  1  2  3  4  5  6  7  8  9 10 
# 25 23 26 21 22 21 19 21 24 28 
# 
# $`2`
# 14 15 16 17 18 19 20 
# 21 18 24 22 23 21 20 

Finally, we rollsum within each vector (where I used width = 4 for illustration)

out <- lapply(agg.split, zoo::rollsum, k = width)
# $`1`
#  2  3  4  5  6  7  8 
# 95 92 90 83 83 85 92 
# 
# $`2`
# 15 16 17 18 
# 85 87 90 86

With indx.vec <- 1:M, this code runs in 0.015 seconds on my machine, versus 12+ seconds for your count.ones. It also scales pretty well as it takes less than 0.4 sec for n1 = n2 = 1e6.