inst/ModelBuild/svrEstimate.R

In waternumbers/FloodForT: Flood Forecasting Toolbox

## ####################################
## server code for Estimate Model tab
## exception is the reactive function to
## estimate the models which is in
## svrEstimateModel.R
## ####################################

## ####################################
## plots for selecting the minima
output$plotMinima <- renderDygraph(
    {
    tmp <- analysisRecord$mdlData$Variables['output']
    
    if(is.null(tmp)){
        return(NULL)
    }else{
         xtmp <- analysisRecord$baseData[,tmp]
         mtmp <- xtmp
         names(mtmp) <- "Minima"
         mtmp[] <- input$sliderMinima
         return(
                    dygraph(cbind(xtmp,mtmp)) %>% 
                        dyOptions(useDataTimezone = TRUE) )# %>%
                #    dySeries(mtmp) )
     }
})

## ######################################
## table of results
output$tblEstMdl <- renderDataTable(
    {
    z <- analysisRecord$mdlTbl
    if(is.null(z)){return(NULL)}

    ## make the table more compact
    z <- z[,c('Nonlinearity','Paths','Lag','Minima','Rt2','Rt2val','MAE','MAEval','AIC','BIC','YIC')]

    ## change non-linearity names
    z[z[,"Nonlinearity"]=="none","Nonlinearity"] <- "None"
    z[z[,"Nonlinearity"]=="plaw","Nonlinearity"] <- "Power Law"
    z[z[,"Nonlinearity"]=="nexp","Nonlinearity"] <- "Negative Exponential"
    z[z[,"Nonlinearity"]=="sig","Nonlinearity"] <- "Sigmoid"

    ## alter the number formats
    for(ii in c('Rt2','Rt2val','MAE','MAEval','AIC','BIC','YIC')){
        z[,ii] <- round(as.numeric(z[,ii]),3)
    }
    colnames(z) <- c('Nonlinearity','Paths','Lag','Minima','Rt2','Rt2\nvalidation','MAE','MAE\nvalidation','AIC','BIC','YIC')    
    return(z)
})

## plots of hydrographs
output$plotEstMdl <- renderDygraph({
    ## check if any names
    nms <- input$selMdl
    if( is.null(nms) ){
        return(NULL)
    }

    ## make observed data
    if(input$ckValData){
        idx <- paste(
            format(analysisRecord$mdlData$Periods['Valid','start'],tz='GMT'),
            format(analysisRecord$mdlData$Periods['Valid','finish'],tz='GMT'),
            sep="::")
    }else{
        idx <- paste(
            format(analysisRecord$mdlData$Periods['Calib','start'],tz='GMT'),
            format(analysisRecord$mdlData$Periods['Calib','finish'],tz='GMT'),
            sep="::")
    }
    tmp <- analysisRecord$mdlData$Variables['output']
    X <- analysisRecord$baseData[idx,tmp]

    ## handle the case of no observations
    if(all(!is.finite(X))){
        return(NULL)
    }
    
    for(ii in nms){
        tmp <- analysisRecord$mdl[[ii]]
        if(input$ckValData){
            tmp <- tmp$val[,'sim']
        }else{
             tmp <- tmp$cal[,'sim']
         }
        names(tmp) <- ii
        X <- cbind(X,tmp)
    }

    return( dygraph(X) %>% 
               dyOptions(useDataTimezone = TRUE) )# %>%
                #    dySeries(mtmp) )
})
                    
## plots of thresholded RMSE
output$plotTrmse <- renderPlot({
    
    nms <- input$selMdl
    if( is.null(input$selMdl)){
        return(NULL)
    }

    ## make array of values
    y <- NULL
    
    for(ii in nms){
        tmp <- analysisRecord$mdl[[ii]]
        if(input$ckValData){
            mse <- tmp$val.perf$mse
            vmse <- tmp$val.perf$vmse
        }else{
            mse <- tmp$cal.perf$mse
            vmse <- tmp$cal.perf$vmse
        }

        ## if values don't exist
        if(is.null(mse) |is.null(vmse)){
            next
        }
        
        if(is.null(y)){
            y <- array(NA,c(nrow(mse),3,length(nms)))
            dimnames(y) <- list(NULL,NULL,nms)
        }
        y[,1,ii] <- sqrt( mse[,'sim'] - 2*sqrt(vmse[,'sim']) )
        y[,2,ii] <- sqrt( mse[,'sim'] )
        y[,3,ii] <- sqrt( mse[,'sim'] + 2*sqrt(vmse[,'sim']) )
    }

    ## if nothing to plot
    if(is.null(y)){return(NULL)}
     
    yrng <- c(0.99*min(y,na.rm=TRUE),1.01*max(y,na.rm=TRUE))
    lty <- rep(1,length(nms))
    names(lty) <- nms
    col <- 1:length(nms)
    names(col) <- nms

    matplot(mse[,'thold'],y[,2,],
            type="l",
            lty=lty,
            col=col,
            ylim=yrng,
            xlab="Threshold",
            ylab="RMSE"            
    )

    for(ii in nms){
        for(jj in 1:(dim(y)[1])){
            segments(mse[jj,'thold'], y[jj,1,ii], y1 = y[jj,3,ii],
                     col = col[ii],lty=lty[ii])
        }
    }

    legendStr <- nms
    legend("topleft",legendStr,
           lty = lty, lwd = 1,
           pch = NULL,
           col = col)

})
   
##############################################
## Parameter table output
output$paramRes <- renderUI(
    {
    if ( is.null(input$selMdl) ){
        return(NULL)
    }
    isolate({
                estCodes <- input$selMdl
                out <- list()
                cnt <- 0
                for(ii in estCodes){
                    ## model title
                    cnt <- cnt + 1
                    out[[cnt]] <- h3(fprettyNames(ii))

                    ## make matrix of parameter values
                    tmp <- analysisRecord$mdl[[ii]]$param
                    nms <- names(tmp$trans)
                    nms <- nms[substr(nms,1,1) != 'Q']
                    pmat <- matrix(NA,length(nms),3)
                    rownames(pmat) <- nms
                    pmat[nms,1] <- tmp$trans[nms]
                    pmat[nms,2] <- ftrans(tmp$raw[nms] - 2*sqrt(diag(tmp$Sigma[nms,nms])),tmp$mdl)
                    pmat[nms,3] <- ftrans(tmp$raw[nms] + 2*sqrt(diag(tmp$Sigma[nms,nms])),tmp$mdl)
                    
                    ## loop paths
                    for(jj in 1:tmp$mdl[1,'np']){
                        cnt <- cnt+1
                        out[[cnt]] <- p(paste("Pathway",jj))
                        cnt <- cnt + 1
                        str <- paste("T",jj,sep=".")
                        out[[cnt]] <- p(paste0("Time constant [h]: ",
                                               round(pmat[str,1]/3600,3),
                                               " (",
                                               round(pmat[str,2]/3600,3),
                                               ",",
                                               round(pmat[str,3]/3600,3),
                                               ")"
                                               ))
                        cnt <- cnt + 1
                        str <- paste("SSG",jj,sep=".")
                        out[[cnt]] <- p(paste0("Steady State Gain: ",
                                               round(pmat[str,1],3),
                                               " (",
                                               round(pmat[str,2],3),
                                               ",",
                                               round(pmat[str,3],3),
                                               ")"
                                               ))
                        
                    }
                    ## do non-linear parameters is any
                    nms <- nms[substr(nms,1,3)=="phi"]
                    if(length(nms)>0){
                        cnt <- cnt+1
                        out[[cnt]] <- p("Non-linearity parameters")
                        for(jj in nms){
                            cnt <- cnt+1
                            out[[cnt]] <- p(paste0(jj,": ",
                                                   round(pmat[jj,1],3),
                                                   " (",
                                                   round(pmat[jj,2],3),
                                                   ",",
                                                   round(pmat[jj,3],3),
                                                   ")"
                                                   ))
                        }
                    }
                        
                                       

## round(
##                         tmp <- analysisRecord$mdl[[ii]]$param
##                         nms <- names(tmp$trans)
##                         nms <- nms[substr(nms,1,1) != 'Q']

                    
##                     pmat <- matrix(NA,length(nms),3)
##                     rownames(pmat) <- nms
                    
##                     pmat[nms,1] <- tmp$trans[nms]
##                     print(rownames(pmat))
##                     pmat[nms,2] <- ftrans(tmp$raw[nms] - 2*sqrt(diag(tmp$Sigma[nms,nms])),tmp$mdl)
##                     pmat[nms,3] <- ftrans(tmp$raw[nms] + 2*sqrt(diag(tmp$Sigma[nms,nms])),tmp$mdl)
##                     cnt <- cnt + 1
##                     out[[cnt]] <- h3(fprettyNames(ii))
##                     for(jj in 1:nrow(pmat)){
##                         cnt <- cnt+1
##                         str <- paste(rownames(pmat)[jj]," : ",pmat[jj,1]," (",pmat[jj,2],",",pmat[jj,3],")")
##                         out[[cnt]] <- p(str)
##                     }
                    ## print("variables in model")
                    ## print(names(tmp))
                    ## print(names(tmp$param))
                    ## pmat <- matrix(NA,
                    ## ttmp <- tmp$param$trans
                    ## print("transformed parameter")
                    ## print(ttmp)
                    ## jj <- which(estCodes %in% ii)
                    ## out[[(2*jj)-1]] <- p(ii)
                    ## out[[2*jj]] <- p(paste( paste(names(ttmp),ttmp),collapse="<br>"))
                }
                return(out)
            })
})
##                 out <- list()
##                 cnt <- 0
##                 for(ii in 1:mdl$mdl[1,'np']){
##                     Tstr <- paste("T",ii,sep=".")
##                     Sstr <- paste("SSG",ii,sep=".")
##                     tmp <- c(
##                         Tlw = param$raw[Tstr] - 2*sqrt(param$Vest[Tstr,Tstr]),
##                         Tup = param$raw[Tstr] + 2*sqrt(param$Vest[Tstr,Tstr])
##                         T = param$raw[Tstr],
##                         SSGlw = param$raw[Sstr] - 2*sqrt(param$Vest[Sstr,Sstr]),
##                         SSgup = param$raw[Sstr] + 2*sqrt(param$Vest[Sstr,Sstr])
##                         SSG = param$raw[Sstr])
##                     tmp <- ftrans(tmp)
##                     cnt <- cnt + 1
##                     out[[cnt]] <- p(paste0("Time Constant: ",
##                                           formatC(tmp['T'],digits=2,format="f"),
##                                           " (",
##                                           formatC(tmp['Tlw'],digits=2,format="f"),
##                                           ",",
##                                           formatC(tmp['Tup'],digits=2,format="f"),
##                                           ")"))
##                     cnt <- cnt + 1
##                     out[[cnt]] <- p(paste0("Steady State Gain: ",
##                                           formatC(tmp['SSG'],digits=2,format="f"),
##                                           " (",
##                                           formatC(tmp['SSGlw'],digits=2,format="f"),
##                                           ",",
##                                           formatC(tmp['SSGup'],digits=2,format="f"),
##                                           ")"))

##                     idx <- which( substr(names(param$raw),1,3)=="phi" )
##                     if(length(idx) > 0){
                        
                        
##                         cnt <- cnt + 1
##                         out[[cnt]] <- h3("Non-linearity parameters")
                        
##                         tmp <- c(
##                             Tlw = param$raw[Tstr] - 2*sqrt(param$Vest[Tstr,Tstr]),
##                         Tup = param$raw[Tstr] + 2*sqrt(param$Vest[Tstr,Tstr])
##                         T = param$raw[Tstr],
##                         SSGlw = param$raw[Sstr] - 2*sqrt(param$Vest[Sstr,Sstr]),
##                         SSgup = param$raw[Sstr] + 2*sqrt(param$Vest[Sstr,Sstr])
##                         SSG = param$raw[Sstr])
##                     tmp <- ftrans(tmp)
                    
##                     ssg <- c( = param$raw[paste("T",ii,sep=".")]
                
##                 out[[1]] <- h3("Model Summaries")

                
##                 ## steady state gains
                    
                
    
##     ## sample of parameters
##     parsamp <- mvrnorm(n = 10000, mdl$th, mdl$V, tol = 1e-6, empirical = FALSE, EISPACK = FALSE)
##     parsamp <- rbind(mdl$th,parsamp)
    
##     # number of nonlinear parameters
##     nnlnpar <- length(fparSep(mdl)$nlpar)
    
##     # seperate out sample so only keep the linear parts
##     parsamp <- parsamp[,(nnlnpar+1):ncol(parsamp)] # parameters of linear system
    
##     # seperate parsamp into roots and poles - alter into parmeter ranges
##     ord <- mdl$numPath
##     p <- parsamp[,1:ord] # poles
##     p[p>1-1e-6] <- 1-1e-6
##     p[p<0] <- 0
##     r <- parsamp[,-(1:ord)] # roots
##     r[r<1e-6] <- 1e-6
##     rm(parsamp)
    
##     # make sure they are matrices to simplify later code
##     r <- as.matrix(r)
##     p <- as.matrix(p)
    
##     # compute summaries
##     ssg  <-  r / (1-p) # steady stage gains
##     ssg[p==0] <- 0
##     et <- -1/log(p)
##     #et <- p /(1-p) # based on deiscrete time equivilent to the exponential distribution
##     et[p==0] <- 0
##     frac <- ssg/rowSums(ssg) # fraction to each path
##     frac[rowSums(ssg)==0,] <- 0
##     conctime <- mdl$lag + rowSums(frac*et)
    
##     # apply the timestep to the concentration time and time constant
##     # convert to hours
 
##     lag <- mdl$lag*info$timeStep /(60*60)
##     et <- et*info$timeStep /(60*60)
##     conctime <- conctime*info$timeStep /(60*60)
    
##     # make output of function which is a list
##     out <- list()

##     if(nnlnpar > 0){
##       for(ii in 1:nnlnpar){
##         str <- "Non linear Parameters +/- 2*s.e.:"
##         str <- paste(str,
##                      formatC(mdl$th[ii],digits=2,format="f"),
##                      "+/-",
##                      formatC(2*sqrt(mdl$V[ii,ii]),digits=2,format="f"))
##       }
##       out[[2]] <- p(str)
##     }else{
##       out[[2]] <- NULL
##     }
##     out[[3]] <- p(paste("Advective time delay [h]:",
##                         formatC(as.numeric(lag),digits=2,format="f")))
##     val <- formatC(c(conctime[1],quantile(conctime[-1],c(0.05,0.95))),
##                   digits=2,format="f") #scientific=FALSE)
##     out[[4]] <- p(paste("Time of Concentration [h]:",val[1],"(",val[2],",",val[3],")"))
##     cnt <- 5
##     for(ii in 1:ord){
##       out[[cnt]] <- h3(paste("Pathway",ii))
##       cnt <- cnt + 1
##       val <- formatC(c(r[1,ii],quantile(r[-1,ii],c(0.05,0.95))),
##                     digits=2,format="f")
##       out[[cnt]] <- p(paste("Root:",val[1],"(",val[2],",",val[3],")"))
##       cnt <- cnt + 1
##       val <- formatC(c(p[1,ii],quantile(p[-1,ii],c(0.05,0.95))),
##                     digits=2,format="f")
##       out[[cnt]] <- p(paste("Pole:",val[1],"(",val[2],",",val[3],")"))
##       cnt <- cnt + 1
##       val <- formatC(c(et[1,ii],quantile(et[-1,ii],c(0.05,0.95))),
##                     digits=2,format="f")
##       out[[cnt]] <- p(paste("Approximate time constant [h]:",val[1],"(",val[2],",",val[3],")"))
##       cnt <- cnt + 1
##       val <- formatC(c(frac[1,ii],quantile(frac[-1,ii],c(0.05,0.95))),
##                     digits=2,format="f")
##       out[[cnt]] <- p(paste("Fraction of effective input received:",val[1],"(",val[2],",",val[3],")"))
##       cnt <- cnt + 1
##     }
##     return(out)
##   })