R/fss.R

In waternumbers/FloodForT: Flood Forecasting Toolbox

#' Generates forecasts for a parrallel path model run with data assimilation
#' 
#' 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()
fss <- 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)
    
    ## 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
    tmp <- X[,'output']
    tmp[] <- NA
    yhat <- NULL
    for(ii in 0:param$mdl[1,'lag']){
        yhat <- cbind(yhat,tmp)
    }
    names(yhat) <- paste("f",0:param$mdl[1,'lag'],sep=".")

    ## set up the model matrices
    idx <- 1:param$mdl[1,'np']
    a <- exp(-dt/param$trans[paste('T',idx,sep=".")])
    b <- param$trans[paste('SSG',idx,sep=".")]*(1-a)
    q <- param$trans[paste('Q',idx,sep=".")]
    AA <- diag(a,length(a))
    A <- matrix(0,param$mdl[1,'lag']*param$mdl[1,'np'],param$mdl[1,'lag']*param$mdl[1,'np'])
    AQ <- A
    kdx <- idx
    for(ii in 1:param$mdl[1,'lag']){
        Ap <- AA
        jdx <- idx + (ii-1)*param$mdl[1,'np']
        for(jj in 0:(param$mdl[1,'lag']-ii)){
            AQ[kdx,jdx] <- Ap
            if(jj == (param$mdl[1,'lag']-ii)){
                A[kdx,jdx] <- Ap
            }
            Ap <- Ap %*% AA
            jdx <- jdx + param$mdl[1,'np']
        }
        kdx <- kdx + param$mdl[1,'np']
    }

    Q <- AQ %*% diag( rep( q , param$mdl[1,'lag']) ) %*% t(AQ)
    
    B <- matrix(0,param$mdl[1,'lag']*param$mdl[1,'np'],param$mdl[1,'lag'])
    H <- B
    jdx <- idx
    for(ii in 1:param$mdl[1,'lag']){
        Ap <- diag(ncol(AA))
        for(jj in 0:(param$mdl[1,'lag']-ii)){
            B[jdx,ii+jj] <- Ap%*% b
            Ap <- Ap %*% AA
        }
        H[jdx,ii] <- 1
        jdx <- jdx + param$mdl[1,'np']
    }

    H <- t(H)
    H <- H[param$mdl[1,'lag']:1,,drop=FALSE]

    ## initialise the states
    denSSG <- sum(param$trans[which(substr(nms,1,3)=="SSG")])
    x <- matrix(c( rep(0,(param$mdl[1,'lag']-1)*param$mdl[1,'np']),
           y0 * param$trans[paste('SSG',idx,sep=".")] / denSSG ),ncol=1)
    P <- diag(1e6,nrow(x))
    x <- A%*%x # forecast states

    HH <- H[1,,drop=FALSE]#nrow(H),,drop=FALSE]

    for(tt in (t0+1):length(u)){
        if( is.finite(X[tt,'output']) ){
            yy <- param$mdl[1,'minima']  + HH%*% x
            vv <- HH%*%P%*%t(HH) + 1
            K <- (P %*% t(HH))/as.numeric(vv)
            x <- x + K*as.numeric( (X[tt,'output']-yy) )
            P <- P - K%*%HH%*%P
        }

        ## filtered estimate
        yfilt <- HH%*%x + param$mdl[1,'minima']
        ## forecast
        uu <- matrix(u[tt:(tt+1-param$mdl[1,'lag'])],ncol=1)

        x <- A%*%x + B%*%uu
        P <- A%*%P%*%t(A) + Q

        ## observe
        yhat[tt,] <- c(yfilt , param$mdl[1,'minima']  + H%*%x)
    }
        
    return(yhat)
}