[R] partial cumsum

[Tony Plate]

William Dunlap wrote:
> 
> Bill Dunlap
> Spotfire, TIBCO Software
> wdunlap tibco.com  
> 
>> -----Original Message-----
>> From: r-help-bounces at r-project.org 
>> [mailto:r-help-bounces at r-project.org] On Behalf Of smu
>> Sent: Wednesday, November 11, 2009 7:58 AM
>> To: r-help at r-project.org
>> Subject: [R] partial cumsum
>>
>> Hello,
>>
>> I am searching for a function to calculate "partial" cumsums.
>>
>> For example it should calculate the cumulative sums until a 
>> NA appears,
>> and restart the cumsum calculation after the NA.
>>
>> this:
>>
>> x <- c(1, 2, 3, NA, 5, 6, 7, 8, 9, 10)
>>
>> should become this:
>>
>> 1 3 6 NA 5  11  18  26  35  45
> 
> Perhaps
>    > ave(x, rev(cumsum(rev(is.na(x)))), FUN=cumsum)
>     [1]  1  3  6 NA  5 11 18 26 35 45
> 

Nice simple function!  Here's a different approach I use that's faster for long vectors with many NA values.  Note however, that this approach can suffer from catastrophic round-off error because it does a cumsum over the whole vector (after replacing NA's with zeros) and then subtracting out the cumsum at the most recent NA values.  Most of the body of this function is devoted to allowing (an unreasonable degree of) flexibility in specification of where to reset.

cumsum.reset <- function(x, reset.at=which(is.na(x)), na.rm=F) {
    # compute the cumsum of x, resetting the cumsum to 0 at each element indexed by reset.at
    if (is.logical(reset.at)) {
        if (length(reset.at)>length(x)) {
            if ((length(reset.at) %% length(x))!=0)
                stop("length of reset.at must be a multiple of length of x")
            x <- rep(x, len=length(reset.at))
        } else if (length(reset.at)<length(x)) {
            if ((length(x) %% length(reset.at))!=0)
                stop("length of x must be a multiple of length of reset.at")
            reset.at <- rep(reset.at, len=length(x))
        }
        reset.at <- which(reset.at)
    } else if (!is.numeric(reset.at)) {
        stop("reset.at must be logical or numeric")
    }
    if (length(i <- which(reset.at<=1)))
        reset.at <- reset.at[-i]
    if (any(i <- is.na(x[reset.at])))
        x[reset.at[i]] <- 0
    if (na.rm && any(i <- which(is.na(x))))
        x[i] <- 0
    y <- cumsum(x)
    d <- diff(c(0, y[reset.at-1]))
    r <- rep(0, length(x))
    r[reset.at] <- d
    return(y - cumsum(r))
}



> x <- c(1, 2, 3, NA, 5, 6, 7, 8, 9, 10)
> cumsum.reset(x)
 [1]  1  3  6  0  5 11 18 26 35 45
> ave(x, rev(cumsum(rev(is.na(x)))), FUN=cumsum)
 [1]  1  3  6 NA  5 11 18 26 35 45
> 

The speedup from not breaking the input vector into smaller vectors is (to me) surprisingly small -- only a factor of 3:

> x <- replace(rnorm(1e6), sample(1e6, 10000), NA)
> all.equal(replace(ave(x, rev(cumsum(rev(is.na(x)))), FUN=cumsum), is.na(x), 0), cumsum.reset(x))
[1] TRUE
> system.time(cumsum.reset(x))
   user  system elapsed 
   0.31    0.03    0.35 
> system.time(ave(x, rev(cumsum(rev(is.na(x)))), FUN=cumsum))
   user  system elapsed 
   0.99    0.05    1.15 
> 

So, I'd go with the ave() approach unless this type of cumsum is the core of a long computationally intensive job.  And if that's the case, it would make sense to code it in C.

-- Tony Plate

> Bill Dunlap
> Spotfire, TIBCO Software
> wdunlap tibco.com 
>  
>> any ideas?
>>
>> thank you and best regards,
>>
>>     stefan
>>
>> ______________________________________________
>> R-help at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide 
>> http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
>>
> 
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>