R/fpp.R

In waternumbers/FloodForT: Flood Forecasting Toolbox

#' Simulates a parrallel path model
#' 
#' Return xts object of model output
#' 
#' @param param parameter and model structure
#' @param X an xts object of data
#' 
#' @keywords FloorForT
#' @export
#' @examples
#' # Not Run
#' # ModelBuild()
fpp <- function(param,X){

    ## time step of data
    dt <- unique( diff( as.numeric(index(X)) ) )

    ## fail if time steps are different
    if(dt != param$mdl[1,'dt']){
        stop("Incorrect time step in data")
    }

    ## switch depending on nl
    fnl <- switch(as.character(param$mdl[1,'nl']),
                  none = function(x,p){rep(1)},
                  plaw = function(x,p){x^p},
                  nexp = function(x,p){ 1 - exp(p*x) },
                  sig = function(x,p){1 / (1+exp(-p[1]*(x-p[2])))}
                  )

    ## compute the input
    nms <- names(param$trans)
    idx <- which(substr(nms,1,3)=="phi")
    s <- pmax(0,(X[,'output']-param$mdl[1,'minima']),na.rm=TRUE)
    u <- fnl(s,param$trans[idx])*pmax(0,X[,'input'],na.rm=TRUE)
    u <- filter(u,c(rep(0,param$mdl[1,'lag']),1),'convolution',sides=1)

    
    ## find initial time and states
    t0 <- ftinit(X,param$mdl)
    t0 <- which(index(X)==t0)
    y0 <- max(0,(X[t0,'output']-param$mdl[1,'minima']))

    ## initialise the output
    yhat <- X[,'output']
    yhat[] <- param$mdl[1,'minima']
    names(yhat) <- 'sim'
    yhat[1:t0] <- NA

    denSSG <- sum(param$trans[which(substr(nms,1,3)=="SSG")])
    idx <- 1:param$mdl[1,'np']
    a <- exp(-dt/param$trans[paste('T',idx,sep=".")])
    b <- param$trans[paste('SSG',idx,sep=".")]*(1-a)
    yinit <- y0 * param$trans[paste('SSG',idx,sep=".")] / denSSG

    for(ii in 1:param$mdl[1,'np']){
        uu <- u*b[ii]
        uu <- uu[-(1:t0)]
        yy <- filter(uu,a[ii],'recursive',init = yinit[ii])
        yhat[-(1:t0)] <- yhat[-(1:t0)] + as.numeric(yy)
    }
    
    return(yhat)
}