inst/ModelBuild/svrEstimateModel.R

In waternumbers/FloodForT: Flood Forecasting Toolbox

###################################################
# reactive function to estimate models when button pressed
###################################################

observe({
    if (input$bttnEstMdl == 0) {return(NULL)}
    isolate({
        ## check model data specified
        if(length(analysisRecord$mdlData)==0){
            str <- "Please select calibration data"
            session$sendCustomMessage("messageBox", str)
            return(NULL)
        }

        ## read the inputs to local variables
        numPath <- input$sliderNumPath[1]:input$sliderNumPath[2]
        timeDelay <- input$sliderTimeDelay[1]:input$sliderTimeDelay[2]
        nonLin <- input$selNonLin
        minima <- input$sliderMinima

        if(length(nonLin)==0){
            str <- "Please specify a non-linearity"
            session$sendCustomMessage("messageBox", str)
            return(NULL)
        }
            

        ## read existing results into local variables
        mdlTbl <- analysisRecord$mdlTbl
        mdlrec <- analysisRecord$mdl
                
        ## make unique codes for all models to be estimated
        np <- rep(numPath,length(timeDelay)*length(nonLin))
        td <- rep(timeDelay,length(nonLin),each=length(numPath))
        nl <- rep(nonLin,each=length(timeDelay)*length(numPath))
        minimaStr <- paste(minima)
        estCodes <- paste(nl,np,td,minimaStr,sep="_")

        ## trim codes to those not already estimated
        if(!is.null(mdlTbl)){
            idx <- !(estCodes %in% rownames(mdlTbl))
            estCodes <- estCodes[idx]
            np <- np[idx]
            td <- td[idx]
            nl <- nl[idx]
        }

        ## specify number of models to estimate and leave if zero
        nmdl <- length(estCodes)
        if(nmdl==0){
            return(NULL)
        }

        ## add additonal rows to mdltbl
        nmdl <- length(estCodes)
        tmp <- data.frame(Paths=np,
                          Lag=td,
                          Nonlinearity=nl,
                          Minima = rep(minima,nmdl),
                          Rt2 = rep(NA,nmdl),
                          Rt2val = rep(NA,nmdl),
                          MAE = rep(NA,nmdl),
                          MAEval = rep(NA,nmdl),
                          CRPS = rep(NA,nmdl),
                          AIC = rep(NA,nmdl),
                          BIC = rep(NA,nmdl),
                          YIC = rep(NA,nmdl),
                          hasDA= rep(FALSE,nmdl),
                          stringsAsFactors = FALSE,
                          row.names = estCodes)
        mdlTbl <- rbind(mdlTbl,tmp[estCodes,])
        
        ## ###### MAKE DATA SERIES #############
        ## calibration
        tmp <- analysisRecord$mdlData$Variables[c('output','input')]
        idx <- paste(
            format(analysisRecord$mdlData$Periods['Calib','start'],tz='GMT'),
            format(analysisRecord$mdlData$Periods['Calib','finish'],tz='GMT'),
            sep="::")
        Xcal <- analysisRecord$baseData[idx,tmp]
        names(Xcal) <- c('output','input')

        if(all(!is.finite(Xcal[,'output'])) | all(!is.finite(Xcal[,'input']))){
            str <- "Please select a calibration period with data"
            session$sendCustomMessage("messageBox", str)
            return(NULL)
        }

        
        ## validation
        tmp <- analysisRecord$mdlData$Variables[c('output','input')]
        idx <- paste(
            format(analysisRecord$mdlData$Periods['Valid','start'],tz='GMT'),
            format(analysisRecord$mdlData$Periods['Valid','finish'],tz='GMT'),
            sep="::")
        Xval <- analysisRecord$baseData[idx,tmp]
        names(Xval) <- c('output','input')
        
        runVal <- TRUE
        if(all(!is.finite(Xval[,'output'])) | all(!is.finite(Xval[,'input']))){
            str <- "Validation period is not run. It has no observations"
            session$sendCustomMessage("messageBox", str)
            runVal <- FALSE
        }

        ## specify the thresholds for the crps and rmse
        thold<- quantile(Xcal[,"output"],seq(0,1,0.01),na.rm=TRUE,names=FALSE)

        ## ## make a minimal varTbl for storage with model
        ## tmp <- analysisRecord$mdlData$Variables[c('output','input')]
        ## vTbl <- analysisRecord$varTbl[tmp,]
        ## names(Xval) <- c('output','input')
        
        ## runVal <- TRUE
        ## if(all(!is.finite(Xval[,'output'])) | all(!is.finite(Xval[,'input']))){
        ##     str <- "Validation period is note run. It has no observations"
        ##     session$sendCustomMessage("messageBox", str)
        ##     runVal <- FALSE
        ## }

        ## ## specify the thresholds for the crps and rmse
        ## thold<- quantile(Xcal[,"output"],seq(0,1,0.01),na.rm=TRUE,names=FALSE)

        ## make a minimal varTbl for storage with model
        tmp <- analysisRecord$mdlData$Variables[c('output','input')]
        vTbl <- analysisRecord$varTbl[tmp,]
        rownames(vTbl) <- c('output','input')
        vTbl <- vTbl[,colSums(vTbl)>0]
        
        ## ######## Estimation Loop #########
        ## initialise the progressbar
        withProgress(message = paste('Estimating',nmdl,'models'),
                     detail = 'This may take a while...', value = 0,{ 
                         for (ucode in estCodes){
                             
                             ## define model structure
                             mdl <- data.frame(
                                 np = mdlTbl[ucode,'Paths'],
                                 lag = mdlTbl[ucode,'Lag'],
                                 nl = mdlTbl[ucode,'Nonlinearity'],
                                 minima = mdlTbl[ucode,'Minima'])
                                   
                             ## set initial parameters
                             tmp <- seq(10,100,length=mdl[1,'np'])
                             thpp <- c(log(tmp),
                                       log( rep(1,mdl[1,'np']) ) )
                             names(thpp) <- c(paste("T",1:mdl[1,'np'],sep="."),
                                              paste("SSG",1:mdl[1,'np'],sep="."))
                             
                             lb <- c(log(rep(1/24,mdl[1,'np'])),rep(-Inf,mdl[1,'np']))
                             ub <- c(log(rep(1e5,mdl[1,'np'])),rep(300,mdl[1,'np']))

                             if( mdl[1,'nl'] =="none" ){
                                 thZ <- thpp
                             }
                             if( mdl[1,'nl'] =="plaw" ){
                                 thZ <- c(thpp,phi.1=log(0.3))
                                 lb <- c(lb,-15)
                                 ub <- c(ub,0)
                             }
                             if( mdl[1,'nl'] =="nexp" ){
                                 thZ <- c(thpp,phi.1=log(1))
                                 lb <- c(lb,-Inf)
                                 ub <- c(ub,300)
                             }                             
                             if( mdl[1,'nl'] =="sig" ){
                                 rng <- range(Xcal[,'output'],na.rm=TRUE)
                                 p2 <- mean(rng)
                                 p1 <- -( log((1-0.01)/0.01) / (rng[1]-p2) )
                                 thZ <- c(thpp,phi.1=log(p1),phi.2=log(p2))
                                 lb <- c(lb,-Inf,-Inf)
                                 ub <- c(ub,300,300)
                             }                             

                             ## call optimisation routine
                             
                             out <- suppressWarnings(
                                 tryCatch(nls(~fppest(par,X,mdl,stype,pNms),
                                              data = list(X=Xcal,mdl=mdl,stype="error",pNms=names(thZ)),
                                              start = list(par = thZ),
                                              control = nls.control(warnOnly = TRUE),trace=TRUE,
                                              lower=lb,upper=ub,algorithm='port'),
                                          error = function(e) e)
                                 )
                             
##                             out <- NA

                             ## message box if it fails
                             if (inherits(out, "error")){
                                 str <- paste("Optimisation Error for",ucode,":", out$message)
                                 session$sendCustomMessage("messageBox", str)
                                 next
                             }

                             ## make parameter structure
##                             param <- fppest(thZ,Xcal,mdl,"param",names(thZ))
                             param <- fppest(coef(out),Xcal,mdl,"param",names(thZ))
                             
                             ## compute varaiance
                             ##Vest <- cov(matrix(rnorm(1000*length(thZ)),1000,length(thZ)))
                             Vest <- suppressWarnings( tryCatch(vcovHAC(out),error=function(e)e) )
                             if(inherits(Vest, "error")){
                                 str <- paste("Error in variance computation for",ucode,":", Vest$message)
                                 session$sendCustomMessage("messageBox", str)
                                 Vest <- matrix(NA,length(thZ),length(thZ)) # default value so we don't have to handle it as missing
                             }
                             colnames(Vest) <- rownames(Vest) <- names(thZ)
                             param$Sigma <- Vest
                             

                             
                             ## run model to get the calibration summaries
                             cal <- fpp(param,Xcal)
                             e <- cbind(Xcal[,'output'],cal)
                             e[,2] <- e[,1]-e[,2]
                             idx <- which(is.finite(e[,1]))
                             t0 <- ftinit(e,param$mdl)
                             t0 <- which(index(e)==t0)
                             idx <- idx[idx>t0]
                                                                
                             mse <- mean(e[idx,2]^2)
                             mae <- mean(abs(e[idx,2]))
                             den <- mean( (e[idx,'output'] - mean(e[idx,'output']))^2 )
                             aic <- 2*(length(param$trans)+1) + length(idx)*(log(2*pi)+log(mse)+1)
                             bic <- log(length(idx))*(length(param$trans)+1) + length(idx)*(log(2*pi)+log(mse)+1)
                             yic <- sum(log( diag(param$Sigma) / (param$raw^2) )) + log(mse/den)
                             ## populate the model output table
                             mdlTbl[ucode,c('Rt2','MAE','AIC','BIC','YIC')] <-
                                 c( 1-(mse/den) , mae , aic , bic , yic )
                                   
                             ## populate mse and crps threshold table
                             mse <- data.frame(thold = thold,
                                               sim = rep(NA,length(thold)))
                             vmse <- crps <- vcrps <- mse
                             for(tt in 1:length(thold)){
                                 jdx <- idx[ e[idx,1] >= thold[tt] ]
                                 mse[tt,'sim'] <- mean(e[jdx,2]^2)
                                 vmse[tt,'sim'] <- var(e[jdx,2]^2)/(length(jdx)-1)
                                 crps[tt,'sim'] <- mean( abs(e[jdx,2]) )
                                 vcrps[tt,'sim'] <- var( abs(e[jdx,2]) )/(length(jdx)-1)
                             }

                             cal.perf <- list(mse=mse,
                                              vmse=vmse,
                                              crps=crps,
                                              vcrps=vcrps)

                             ## VALIDATION ###
                             if(runVal){
                                 val <- fpp(param,Xval)
                                 e <- cbind(Xval[,'output'],val)
                                 e[,2] <- e[,1]-e[,2]
                                 idx <- which(is.finite(e[,1]))
                                 t0 <- ftinit(e,param$mdl)
                                 t0 <- which(index(e)==t0)
                                 idx <- idx[idx>t0]
                                 
                                 mse <- mean(e[idx,2]^2)
                                 mae <- mean(abs(e[idx,2]))
                                 den <- mean( (e[idx,'output'] - mean(e[idx,'output']))^2 )
                                 ## populate the model output table
                                 mdlTbl[ucode,c('Rt2val','MAEval')] <-
                                     c( 1- (mse/den) , mae)
                                 
                                 ## populate mse and crps threshold table
                                 mse <- data.frame(thold = thold,
                                                   sim = rep(NA,length(thold)))
                                 vmse <- crps <- vcrps <- mse
                                 for(tt in 1:length(thold)){
                                     jdx <- idx[ e[idx,1] >= thold[tt] ]
                                     mse[tt,'sim'] <- mean(e[jdx,2]^2)
                                     vmse[tt,'sim'] <- var(e[jdx,2]^2)/(length(jdx)-1)
                                     crps[tt,'sim'] <- mean( abs(e[jdx,2]) )
                                     vcrps[tt,'sim'] <- var( abs(e[jdx,2]) )/(length(jdx)-1)
                                 }
                                 
                                 val.perf <- list(mse=mse,
                                                  vmse=vmse,
                                                  crps=crps,
                                                  vcrps=vcrps)
                             }else{
                                  val <- NULL
                                  val.perf <- NULL
                              }
                                      

                             param$varTbl <- vTbl
                             param$lvl <- analysisRecord$mdlData$lvl
                             ## add model output to record
                             mdlrec[[ucode]] <- list(param = param,
                                                     data = analysisRecord$mdlData,
                                                     cal=cal,
                                                     val=val,
                                                     cal.perf=cal.perf,
                                                     val.perf=val.perf)
                             ## update progressbar
                             incProgress(1/nmdl)
                         }
                     })
 


        ## put the results back into analysisRecord
        analysisRecord$mdlTbl <- mdlTbl
        analysisRecord$mdl <- mdlrec
    })
})