inst/ModelBuild/svrDAEstimateModel.R

In waternumbers/FloodForT: Flood Forecasting Toolbox

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

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

                if(length(input$selDAest)==0){
                    str <- "Please select at least one model"
                    session$sendCustomMessage("messageBox", str)
                    return(NULL)
                }

                ## read the inputs to local variables
                estCodes <- input$selDAest

                ## read existing results into local variables
                mdlTbl <- analysisRecord$mdlTbl
                mdlrec <- analysisRecord$mdl
 
                ## make all models to be estimated
                idx <- !(estCodes %in% rownames(mdlTbl[mdlTbl[,'hasDA']==TRUE,]))
                estCodes <- estCodes[idx]
                
                ## specify number of models to estimate and leave if zero
                nmdl <- length(estCodes)
                if(nmdl==0){
                    return(NULL)
                }

                ## ###### 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
                }

                ## ######## Estimation Loop #########
                ## initialise the progressbar
                withProgress(message = paste('Estimating',nmdl,'models'),
                             detail = 'This will take a while...have a cup of tea (or two)', value = 0,
                             {for (ucode in estCodes){
                                  
                                  ## get the existing varialbes
                                  param <- mdlrec[[ucode]]$param
                                  cal <- mdlrec[[ucode]]$cal
                                  val <- mdlrec[[ucode]]$val
                                  cal.perf <- mdlrec[[ucode]]$cal.perf
                                  val.perf <- mdlrec[[ucode]]$val.perf

                                  ## set initial parameters
                                  thZ <- rep(0,param$mdl[1,'np'])
                                  
                                  lb <- rep(-Inf,param$mdl[1,'np'])
                                  ub <- rep(13,param$mdl[1,'np'])

                                  ## call optimisation routine
                                  out <- suppressWarnings(
                                      tryCatch(nls(~fssest(par,X,param,stype),
                                                   data = list(X=Xcal,param=param,stype="error"),
                                                   start = list(par = thZ),
                                                   control = nls.control(warnOnly = TRUE,tol=1e-3),trace=TRUE,
                                                   lower=lb,upper-ub,algorithm='port'),
                                               error = function(e) e)
                                      )
                                  
                                  ## message box if it fails
                                  if (inherits(out, "error")){
                                      str <- paste("Optimisation Error : ", out$message)
                                      session$sendCustomMessage("messageBox", str)
                                      next
                                  }
                                  
                                  ## make parameter structure
                                  param <- fssest(coef(out),Xcal,param,"param")

                                  ## CALIBRATION
                                  ## run model to get the calibration forecasts
                                  
                                  cal <- cbind(cal,fss(param,Xcal))

                                  ## estimate the uncertainty terms
                                  for(ii in 0:param$mdl[1,'lag']){

                                      ## make names
                                      nm <- paste("f",ii,sep=".")

                                      ## extract and alter time stamp
                                      xx <- cal[,nm]
                                      index(xx) <- index(xx) + param$mdl[1,'dt']*ii
                                      
                                      ## bind with Xcal
                                      e <- cbind(Xcal[,'output'],xx)

                                      ## make uncertainty estimate
                                      param[[nm]] <- qdks(e[,1]-e[,2],e[,2])
                                      
                                      ## work out crps
                                      crps <- cbind(e[,1]-e[,2],qdks(e[,1]-e[,2],e[,2],xout=e[,2])$y)
                                      crps <- apply(crps,1,FUN=function(x){ mean( abs(x[1]-x[-1])) - 0.5*mean(abs(outer(x[-1],x[-1],"-"))) })
                                      crps <- as.xts(crps,index=index(e))
                                      e <- cbind(e,crps=crps)
                                      e[,2] <- (e[,1]-e[,2])^2
                                      
                                      idx <- which(is.finite(e[,1]))
                                      t0 <- ftinit(e,param$mdl)
                                      t0 <- which(index(e)==t0)+1+ii
                                      idx <- idx[idx>t0]
                                      
                                      ## compute thresholded values
                                      thold <- cal.perf$mse[,'thold']
                                      for(tt in 1:length(thold)){
                                          jdx <- idx[ e[idx,1] >= thold[tt] ]
                                          cal.perf$mse[tt,nm] <- mean(e[jdx,2])
                                          cal.perf$vmse[tt,nm] <- var(e[jdx,2])/(length(jdx)-1)
                                          cal.perf$crps[tt,nm] <- mean( e[jdx,3]) 
                                          cal.perf$vcrps[tt,nm] <- var( e[jdx,3]) /(length(jdx)-1)
                                      }
                                  }
                                      
                                  ## VALIDATION ##
                                  if(runVal){
                                      ## run model to get the calibration forecasts
                                      val <- cbind(val,fss(param,Xval))
                                      
                                      ## estimate the uncertainty terms
                                      for(ii in 0:param$mdl[1,'lag']){
                                          ## make names
                                          nm <- paste("f",ii,sep=".")
                                          
                                          ## extract and alter time stamp
                                          xx <- val[,nm]
                                          index(xx) <- index(xx) + param$mdl[1,'dt']*ii
                                          
                                          ## bind with Xcal
                                          e <- cbind(Xval[,'output'],xx)
                                          
                                          ## work out crps
                                          crps <- cbind(e[,1]-e[,2],qdks(e[,1]-e[,2],e[,2],xout=e[,2])$y)
                                          crps <- apply(crps,1,FUN=function(x){ mean( abs(x[1]-x[-1])) - 0.5*mean(abs(outer(x[-1],x[-1],"-"))) })
                                          crps <- as.xts(crps,index=index(e))
                                          e <- cbind(e,crps=crps)
                                          e[,2] <- (e[,1]-e[,2])^2
                                      
                                          idx <- which(is.finite(e[,1]))
                                          t0 <- ftinit(e,param$mdl)
                                          t0 <- which(index(e)==t0)+1+ii
                                          idx <- idx[idx>t0]
                                          
                                      ## compute thresholded values
                                          thold <- val.perf$mse[,'thold']
                                          for(tt in 1:length(thold)){
                                              jdx <- idx[ e[idx,1] >= thold[tt] ]
                                              val.perf$mse[tt,nm] <- mean(e[jdx,2])
                                              val.perf$vmse[tt,nm] <- var(e[jdx,2])/(length(jdx)-1)
                                              val.perf$crps[tt,nm] <- mean( e[jdx,3]) 
                                              val.perf$vcrps[tt,nm] <- var( e[jdx,3]) /(length(jdx)-1)
                                          }
                                      }
                                  }
                              

                                  ## add model output to record
                                  mdlrec[[ucode]]$param <- param
                                  mdlrec[[ucode]]$cal <- cal
                                  mdlrec[[ucode]]$val <- val
                                  mdlrec[[ucode]]$cal.perf <- cal.perf
                                  mdlrec[[ucode]]$val.perf <- val.perf
                                  mdlTbl[ucode,'hasDA'] <- TRUE                                  
                                  ## update progressbar
                                  incProgress(1/nmdl)
                              }

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