Event Counter - RECOVER salinity envelope

example_code_S79Threshold.r

library(zoo)

# Read data and set up ----------------------------------------------------
dat=read.csv("ex_dat_S79_NA25.csv")
dat$Date=as.Date(dat$Date)
dat$S79.14d=with(dat,rollapply(S79,width=14,FUN=function(x)mean(x,na.rm=T),fill=NA,align="right"))
dat$CRE.low=with(dat,ifelse(S79.14d<750,1,0)) # RECOVER Low flow envelope

# Helper functions for 14-day periods
dat$period.14=c(0,rep(1:(nrow(dat)-1))%/%14)+1
dat$CY=as.numeric(format(dat$Date,"%Y"))
dat$biwk=as.numeric(format(dat$Date,"%j"))%/%15L+1L


# Standard to compare -----------------------------------------------------
## Translated from Excel spreadsheet version (i.e. standard to compare)
dat$CRE.low.count=0 # empty variable filled with zeros

# if the threshold is exceeded during a 14 day period
# if an exceedance occurs then the 14 day period "resets"
for( i in 14:nrow(dat)){
  dat$CRE.low.count[i]=with(dat,ifelse(CRE.low[i]==1&sum(CRE.low.count[(i-13):(i-1)],na.rm=T)==0,1,0))
}
sum(dat$CRE.low.count); #Should be 556


# optimization attempts ---------------------------------------------------
## Attempt at using roll apply
dat$test=with(dat,rollapply(CRE.low,width=13,FUN=function(x) sum(x,na.rm=T),by=13,fill=0,align="right"))
dat$CRE.low.count.roll=with(dat,ifelse(rollapply(CRE.low,width=14,FUN=function(x) sum(x,na.rm=T),by=14,fill=0,align="right")>0,1,0))
sum(dat$CRE.low.count.roll,na.rm=T); #Should be 556
# extra counts
# 
head(dat,28L)
head(dat,56L)

## attempted with sapply
dat$apply.count=0
test.fun=function(i){
  dat$apply.count[i]=with(dat,ifelse(CRE.low[i]==1&sum(apply.count[(i-13):(i-1)],na.rm=T)==0,1,0))
}

dat$apply.count=c(rep(0,13),sapply(14:nrow(dat),FUN=test.fun))
sum(dat$apply.count,na.rm=T)#Should be 556
# extra counts


## Alternate approach
library(plyr)
test=ddply(dat,c("period.14"),summarise,N.period=length(period.14),N.count=ifelse(sum(CRE.low)>0,1,0))
head(test)
sum(test$N.count,na.rm=T)#Should be 556
# extra counts

test=ddply(dat,c("CY","biwk"),summarise,N.period=length(biwk),N.count=ifelse(sum(CRE.low)>0,1,0))
head(test)
sum(test$N.count,na.rm=T)#Should be 556
# extra counts



# benchmarking ------------------------------------------------------------
test=dat
test$loop.count=0
f1=function(test){
  for( i in 14:nrow(test)){
    test$loop.count[i]=with(test,ifelse(CRE.low[i]==1&sum(loop.count[(i-13):(i-1)],na.rm=T)==0,1,0))
  }
}

test$apply.count=0
test.fun=function(i){
  test$loop.count[i]=with(test,ifelse(CRE.low[i]==1&sum(loop.count[(i-13):(i-1)],na.rm=T)==0,1,0))
}
f2=function(test){
  test$apply.count=c(rep(0,13),sapply(14:nrow(test),FUN=test.fun))
}
library(rbenchmark)
benchmark(f1(test))
benchmark(f2(test))

library(microbenchmark)
mbm<-microbenchmark("for loop"=f1(test),
                    "sapply"=f2(test))
library(ggplot2)
autoplot(mbm)