R/dynatop.R

In dynatop: An Implementation of Dynamic TOPMODEL Hydrological Model in R

#' R6 Class for Dynamic TOPMODEL
#' @examples
#' ## the vignettes contains further details of the method calls.
#'
#' data("Swindale") ## example data
#' ctch_mdl <- dynatop$new(Swindale$model) ## create with model
#' ctch_mdl$add_data(Swindale$obs) ## add observations
#' ctch_mdl$initialise() ## initialise model
#' ctch_mdl$sim() ## simulate model
#' @export
dynatop <- R6Class(
            "dynatop",
    public = list(
        #' @description Creates a dynatop class object from the a list based model description as generated by dynatopGIS.
        #'
        #' @param model a dynamic TOPMODEL list object
        #' @param use_states logical if states should be imported
        #' @param drop_map logical if the map should be dropped
        #' @param delta error term in checking redistribution sums
        #'
        #' @return invisible(self) suitable for chaining
        #'
        #' @details This function makes some basic consistency checks on a list representing a dynamic TOPMODEL model. The checks performed and basic 'sanity' checks. They do not check for the logic of the parameter values nor the consistency of states and parameters. Sums of the redistribution matrices are checked to be in the range 1 +/- delta.
        initialize = function(model, use_states=FALSE, delta = 1e-13){
            ## digest model with checks - will fail if errors identified
            private$digest_model(model,use_states,delta)
            invisible(self)
        },
        #' @description Adds observed data to a dynatop object
        #'
        #' @param obs_data an xts object of observed data
        #'
        #' @return invisible(self) suitable for chaining
        #'
        #' @details This function makes some basic consistency checks on the observations to ensure they have uniform timestep and all required series are present.
        add_data = function(obs_data){
            self$clear_data()
            ## check input and get model timestep
            private$check_obs(obs_data)
            private$digest_obs(obs_data)
            invisible(self)
        },
        #' @description Clears all forcing and simulation data except current states
        #'
        #' @return invisible(self) suitable for chaining
        clear_data = function(){
            private$time_series <- list()
            private$info$ts <- list()
        },
        #' @description Initialises a dynatop object in the most simple way possible.
        #'
        #' @param tol tolerance for the solution for the saturated zone
        #' @param max_it maximum number of iterations to use in the solution of the saturated zone
        #'
        #' @return invisible(self) suitable for chaining
        initialise = function(tol = 2*.Machine$double.eps, max_it = 1000){

            private$init_hs(tol,max_it)
            private$init_ch()
            invisible(self)
        },
        #' @description Initialises only the channel part of a dynatop object in the most simple way possible.
        #'
        #' @return invisible(self) suitable for chaining
        initialise_channel = function(){
            private$init_ch()
            invisible(self)
        },
        #' @description Simulate the hillslope output of a dynatop object
        #' @param keep_states a vector of POSIXct objects (e.g. from xts) giving the time stamp at which the states should be kept
        #' @param sub_step simulation timestep in seconds, default value of NULL results in data time step
        #' @param tol tolerance on width of bounds in the solution for the saturated zone
        #' @param max_it maximum number of iterations to use in the solution of the saturated zone
        #' @param ftol tolerance in closeness to 0 in the solution for the saturated zone
        #'
        #' @details Both saving the states at every timestep and keeping the mass balance can generate very large data sets!!
        #' While ftol is implemented it is currently set to \code{Inf} to mimic the behaviour of previous versions. This will change in the future.
        #'
        #' @return invisible(self) for chaining
        sim_hillslope = function(keep_states=NULL,sub_step=NULL,tol = 2*.Machine$double.eps, max_it = 1000, ftol= Inf){

            ## check the solver options
            tol <- as.double(tol); ftol <- as.double(ftol)
            if(any( c(tol,ftol) < .Machine$double.eps)){
                stop("A solution tolerance is set lower then machine eps")
            }
            max_it <- as.integer(max_it)
            if(max_it < 10){stop("Please use at least 10 iterations")}

            ## check presence of finite states
            sv <- c("s_sf","s_rz","s_uz","s_sz")
            has_states <- all(sapply(private$model$hillslope[,sv],
                                     function(x){all(!is.na(x))}))
            if( !has_states ){
                stop("Model states are either not initialised or have non-finite values")
            }

            if( !is.null(sub_step) && !is.finite(sub_step[1]) ){
                stop("sub_step should be a single finite value")
            }

            ## check presense of obs
            if( length(private$time_series$index) < 2 ){
                stop("Insufficent data to perform a simulation")
            }

            ## check keep_states is valid
            if( length(keep_states)>0 ){
                if( !("POSIXct" %in% class(keep_states)) ){
                    stop("Times for returning states should be POSIXct object")
                }
            }
            keep_states <- keep_states[keep_states %in% private$time_series$index]

            ## simulate
            private$sim_hs(keep_states,sub_step[1],tol,max_it,ftol)

            invisible(self)
        },
        #' @description Simulate the channel output of a dynatop object
        #' @return invisible(self) for chaining
        sim_channel=function(){

            if(!private$info$can_solve_channel){
                stop("Cannot simulate channel - check connectivity")
            }

            ## check presence of channel_inflow
            if( nrow(private$time_series$channel_inflow$surface) !=
                length(private$time_series$index) ){
                stop("Suitable channel_inflow not available")
            }
            if( length(private$time_series$index) < 2 ){
                stop("Insufficent data to perform a simulation")
            }
            ## check initialised
            if(!("linear_time" %in% names(private$summary$channel))){
                stop("Channel solution is not initialised")
            }

            private$sim_ch()

            invisible(self)
        },
        #' @description Simulate the hillslope and channel components of a dynatop object
        #' @param keep_states a vector of POSIXct objects (e.g. from xts) giving the time stamp at which the states should be kept
        #' @param mass_check Flag indicating is a record of mass balance errors should be kept
        #' @param sub_step simulation timestep in seconds, default value of NULL results in data time step
        #' @param tol tolerance on width of bounds in the solution for the saturated zone
        #' @param max_it maximum number of iterations to use in the solution of the saturated zone
        #' @param ftol tolerance in closeness to 0 in the solution for the saturated zone
        #'
        #' @details Calls the sim_hillslope and sim_channel in sequence. Both saving the states at every timestep and keeping the mass balance can generate very large data sets!!
        #'
        #' @return invisible(self) for chaining
        sim = function(keep_states=NULL,sub_step=NULL,tol=2*.Machine$double.eps,max_it=1000,ftol=Inf){
            self$sim_hillslope(keep_states,sub_step,tol,max_it)
            self$sim_channel()
            invisible(self)
        },
        ## ############
        ## Functions for extracting and plotting data
        #' @description Return channel inflow as an xts series or list of xts series
        #' @param total logical if plot total inflow is to be plotted
        #' @param separate logical if the surface and saturated zone inflows should be returned separately
        get_channel_inflow = function(total=FALSE,separate=FALSE){
            x <- private$time_series$channel_inflow
            if(total){
                x <- lapply(x,rowSums)
            }
            if(separate){
                x$surface <- xts::xts(x$surface,
                                   order.by=private$time_series$index)
                x$saturated <- xts::xts(x$saturated,
                                   order.by=private$time_series$index)
            }else{
                x <- x$surface + x$saturated
                x <- xts::xts(x,
                              order.by=private$time_series$index)
            }
            return(x)
        },
        #' @description Plot the channel inflow
        #' @param total logical if total inflow is to be plotted
        #' @param separate logical logical if the surface and saturated zone inflows should be plotted separately
        plot_channel_inflow = function(total=FALSE,separate=FALSE){
            x <- self$get_channel_inflow(total,separate)
            if(total){
                if(separate){
                    x <- merge(x$surface,x$saturated)
                    names(x) = c("surface","saturated")
                    lloc <- "topright"
                    plot(x,main="Channel Inflow",legend.loc=lloc)
                }else{
                    lloc <- NULL
                    plot(x,main="Channel Inflow",legend.loc=lloc)
                }
            }else{
                if(separate){
                    lloc <- "topright"
                    oldpar <- par(no.readonly = TRUE)
                    on.exit(par(oldpar))
                    par(mfrow=c(2,1))
                    ## print seems to make this work....
                    print(plot(x$surface,main="Channel Inflow: surface",legend.loc=lloc,on=1))
                    print(plot(x$saturated,main="Channel Inflow: saturated",legend.loc=lloc,on=2))
                }else{
                    lloc <- "topright"
                    plot(x,main="Channel Inflow",legend.loc=lloc)
                }

            }
        },
        #' @description Return flow at the gauges as an xts series
        #' @param gauge names of gauges to return (default is all gauges)
        get_gauge_flow = function(gauge=colnames(private$time_series$gauge_flow)){
            gauge <- match.arg(gauge,colnames(private$time_series$gauge_flow),
                               several.ok=TRUE)
            xts::xts(private$time_series$gauge_flow[,gauge,drop=FALSE],
                     order.by=private$time_series$index)
        },
        #' @description Get the flow at gauges
        #' @param gauge names of gauges to return (default is all gauges)
        plot_gauge_flow = function(gauge=colnames(private$time_series$gauge_flow)){
            plot( self$get_gauge_flow(gauge),main="Simulated Flows",legend.loc="topright")
        },
        #' @description Get the observed data
        get_obs_data = function(){
            xts::xts(private$time_series$obs,
                     order.by=private$time_series$index)
        },
        #' @description Return the model
        get_model = function(){
            private$reform_model()
        },
        #' @description Return the model
        get_mass_errors = function(){
            if( !("mass_balance" %in% names(private$time_series)) ){
                stop("Mass errors are not available")
            }
            xts::xts(private$time_series$mass_balance,
                     order.by=private$time_series$index)
        },
        #' @description Return states
        #' @param record logical TRUE if the record should be returned. Otherwise the current states returned
        get_states = function(record=FALSE){
            
            if( record ){
                #browser()
                tmp <- setNames(private$time_series$state_record,
                                private$time_series$index)
                idx <- sapply(tmp, function(x){length(x)>0})
                return(tmp[idx])
##                return( setNames(private$time_series$state_record,
##                                 private$time_series$index) )
            }else{
                nm <- c("s_sf","s_rz","s_uz","s_sz")
                tmp <- private$model$hillslope[,c("id",nm)] + 0## need to modify to force change of memory...
                return( tmp )
            }
       },
       #' @description Plot a current state of the system
       #' @param state the name of the state to be plotted
       #' @param add_channel Logical indicating if the channel should be added to the plot
       plot_state = function(state,add_channel=TRUE){

           if( is.null(private$model$map$hillslope) ){
               stop("The model contains no map of HSU locations")
           }
           if( !file.exists(private$model$map$hillslope) ){ stop("The model map file is missing") }
           if( add_channel & !file.exists(private$model$map$channel) ){ warnings("File containing the channel network does not exist") }

           if(!(state%in%colnames(private$model$hillslope))){
               stop("Model state does not exist")
           }


           if( !requireNamespace("raster",quietly=TRUE) ){
               stop( "The raster package is required for plotting the maps of states - please install or add to libPath" )
           }

           rst <- raster::raster(private$model$map$hillslope)
           rst <- raster::subs(rst, private$model$hillslope[,c("id",state)])

           raster::plot( rst)
           if( add_channel & file.exists(private$model$map$channel) ){
               chn <- raster::shapefile(private$model$map$channel)
               raster::plot(chn,add=TRUE)
           }

       }

    ),
    private = list(
        ## stores of data
        version = "0.2.2",
        model = list(), # storage for model object
        summary = list(), # storage for intermediate computed values used in code
        time_series = list(),
        info = list(can_sim_channel=FALSE),
        ## decribes the model tables in terms of names, data type and role
        ## roles are either
        ## - attribute
        ## - parameter
        ## - data_series
        ## - state
        ## - output_label
        model_description = list(
            ## hillslope data frame
            hillslope = list(
                "id" = list(class = "integer", opt = "all", min=1, max=Inf),
                "atb_bar" = list(class = "numeric", opt = "all", min=0, max=Inf),
                "s_bar" = list(class = "numeric", opt = "all", min=0, max=Inf),
                "area" = list(class = "numeric", opt = "all", min=0, max=Inf),
                "width" = list(class = "numeric", opt = "all", min=0, max=Inf),
                "opt" = list(class = "character", opt = "all", min=0, max=Inf),
                "s_raf" = list(class = "numeric", opt = "all", min=0, max=Inf),
                "t_raf" = list(class = "numeric", opt = "all", min=0, max=Inf),
                "r_sfmax" = list(class = "numeric", opt = "all", min=0, max=Inf),
                "c_sf" = list(class = "numeric", opt = "all", min=0, max=Inf),
                "s_rzmax" = list(class = "numeric", opt = "all", min=0, max=Inf),
                "t_d" = list(class = "numeric", opt = "all", min=0, max=Inf),
                "ln_t0" = list(class = "numeric", opt=c("exp","bexp","dexp"), min=-Inf, max=Inf),
                "c_sz" = list(class = "numeric", opt = "cnst", min=0, max=Inf),
                "m" = list(class = "numeric", opt = c("exp","bexp","dexp"), min=0, max=Inf),
                "D" = list(class = "numeric", opt = c("cnst","bexp"), min=0, max=Inf),
                "m_2" = list(class = "numeric", opt = "dexp", min=0, max=Inf),
                "omega" = list(class = "numeric", opt = "dexp", min=0, max=1),
                "s_rz0" = list(class = "numeric", opt = "all", min=0, max=Inf),
                "r_uz_sz0" = list(class = "numeric", opt = "all", min=0, max=Inf),
                "s_sf" = list(class = "numeric", opt = "all", min=0, max=Inf),
                "s_rz" = list(class = "numeric", opt = "all", min=0, max=Inf),
                "s_uz" = list(class = "numeric", opt = "all", min=0, max=Inf),
                "s_sz" = list(class = "numeric", opt = "all", min=0, max=Inf)
            ),
            ## channel data frame
            channel =  list(
                "id" = list(class = "integer", opt = "all", min=1, max=Inf),
                "area" = list(class = "numeric", opt = "all", min=0, max=Inf),
                "length" = list(class = "numeric", opt = "all", min=0, max=Inf),
                "v_ch" = list(class = "numeric", opt = "all", min=0, max=Inf)
            ),
            ## linkages between HSUs
            flow_direction = list(
                "from" = list(class = "integer", opt = "all", min=1, max=Inf),
                "to" = list(class = "integer", opt = "all", min=1, max=Inf),
                "frc" = list(class = "numeric", opt = "all", min=0, max=1)
            ),
            ## rainfall inputs
            precip_input = list(
                "name" = list(class = "character", opt = "all", min=NA, max=NA),
                "id" = list(class = "integer", opt = "all", min=1, max=Inf),
                "frc" = list(class = "numeric", opt = "all", min=0, max=1)
            ),
            ## pet inputs
            pet_input = list(
                "name" = list(class = "character", opt = "all", min=NA, max=NA),
                "id" = list(class = "integer", opt = "all", min=1, max=Inf),
                "frc" = list(class = "numeric", opt = "all", min=0, max=1)
            ),
            ## point inflow to channels
            point_inflow = list(
                "name" = list(class = "character", opt = "all", min=NA, max=NA),
                "id" = list(class = "integer", opt = "all", min=1, max=Inf)
            ),
            ## diffuse inflow to channels
            diffuse_inflow = list(
                "name" = list(class = "character", opt = "all", min=NA, max=NA),
                "id" = list(class = "integer", opt = "all", min=1, max=Inf)
            ),
            ## location of gauges
            gauge = list(
                "name" = list(class = "character", opt = "all", min=NA, max=NA),
                "id" = list(class = "integer", opt = "all", min=1, max=Inf)
            )
        ),
        ## private function for checking a data frame against description
        ftblCheck = function(tbl,prop,lbl){
            ## check tbl is a data frame
            if(!is.data.frame(tbl)){
                stop(paste("Table",lbl,"should be a data.frame"))
            }
            ## check all columns are present
            idx <- names(prop) %in% names(tbl)
            if( !all( idx ) ){# check it has required columns
                stop( paste("Table",lbl,"is missing columns:",
                            paste(names(prop)[!idx],collapse=",")) )
            }
            ## see if there is an opt column
            tbl_has_opt <- ("opt" %in% names(tbl)) & (class(tbl$opt)=="character")
            desc_opt_val <- unique(unlist(sapply(prop,function(x){x$opt}))) ## values of options in description
            if(tbl_has_opt){
                idx <- tbl$opt %in% desc_opt_val
                if(!all(idx)){
                    stop("Invalid options values in Table",lbl)
                }
            }
            ## catch case where table has no rows but correct structure
            if( nrow(tbl)==0 ){return(tbl)}

            ## loop columns
            str <- ""
            idx <- 1:nrow(tbl)
            for(ii in names(prop)){
                if(class(tbl[[ii]])==prop[[ii]]$class){
                    ## continue checks
                    if(tbl_has_opt){
                        if("all" %in% prop[[ii]]$opt){
                            idx <- rep(TRUE,nrow(tbl))
                        }else{
                            idx <- tbl$opt %in% prop[[ii]]$opt
                        }
                    }
                    if( prop[[ii]]$class == "character" ){
                        if( any(nchar(tbl[[ii]][idx])) == 0 ){
                            str <- paste(str,
                                         paste("All required strings in",ii,
                                               "should have non-zero length"),
                                         sep="\n")
                        }
                    }else{
                        if( any(is.na(tbl[[ii]][idx])) ){
                            str <- paste(str,
                                         paste("All required values in",ii,"should not be missing"),
                                         sep="\n")
                        }
                        if( (!is.na( prop[[ii]]$min ) & any( tbl[[ii]][idx] < prop[[ii]]$min)) |
                            (!is.na( prop[[ii]]$max ) & any( tbl[[ii]][idx] > prop[[ii]]$max)) ){
                            str <- paste(str,
                                         paste("All required values in",ii,"should be between",
                                               prop[[ii]]$min,"and",prop[[ii]]$max),
                                         sep="\n")
                        }
                    }

                    tbl[[ii]][!idx] <- NA

                }else{
                    ## add message to string
                    str <- paste(str,
                                 paste("Class of",ii,"should be",prop[[ii]]$class),
                                 sep="\n")
                }
            }

            if(nchar(str)>0){
                str <- paste(paste("table",lbl,"has the following errors:"),
                             str,sep="\n")
                stop(str)
            }
            return(tbl)
        },
        ## this code checks and digests the model
        digest_model = function(model, use_states, delta=1e-13){

            ## initialise store of required data series names
            req_names <- list(output_label = list(),
                              data_series = list())


            ## check all components of the model exist
            components <- names(private$model_description)
            idx <- components %in% names(model)
            if( !all(idx) ){
                stop(paste("Missing components:",paste(components[!idx],collapse=",")))
            }

            ## ################################################
            ## Digest hillslope table
            ## ################################################
            hsd <- private$model_description$hillslope
            ## adjust options on states so set to NA if required
            if(!use_states){
                for(ii in c("s_sf","s_rz","s_uz","s_sz")){
                    hsd[[ii]]$opt <- ""
                }
            }
            ## check table
            model$hillslope <- private$ftblCheck(model$hillslope,hsd,"hillslope")
            ## sort by id - this is required
            model$hillslope <-
                    model$hillslope[order(model$hillslope$id,decreasing=TRUE),]

            ## ################################################
            ## Digest channel table
            ## ################################################
            ## check table
            model$channel <- private$ftblCheck(model$channel,private$model_description$channel,"channel")
            ## sort by id - this is required
            model$hillslope <-
                    model$hillslope[order(model$hillslope$id,decreasing=TRUE),]

            ## ################################################
            ## check sequential numbering of ids
            ## ################################################
            all_hru_id<- c(model$hillslope$id,model$channel$id)
            if( length(all_hru_id) != length(unique(all_hru_id)) ){
                stop("HRU id values should be unique") }
            if( !all(is.finite(all_hru_id)) ){ stop("HRU id values should be finite") }
            if( !all(range(all_hru_id)==c(1,length(all_hru_id))) ){
                stop("HRU id's should be numbered consecuativly from 1")
            }

            ## ################################################
            ## Check flow direction
            ## ################################################
            ## check table
            model$flow_direction <- private$ftblCheck(model$flow_direction,private$model_description$flow_direction,
                                              "flow_direction")
            ## General checks on flow directions
            if(!all( model$flow_direction$from %in% all_hru_id)){
                stop("Flow link from an unspecified HRU")
            }
            if(!all( model$flow_direction$to %in% all_hru_id)){
                stop("Flow link to unspecified HSU")
            }
            if(!all( model$flow_direction$to < model$flow_direction$from )){
                stop("Flow link going to a HSU with higher id - out of order")
            }
            tmp <- tapply(model$flow_direction$frc,model$flow_direction$from,sum)
            idx <- abs(tmp-1)<delta
            if(!all(idx)){
                stop("The fractions in the flow directions of the followin HRU's do not sum to 1: ",
                     paste(names(tmp[!idx]),collapse=", "))
            }
            ## specific check on hillslope
            if(!all(model$hillslope$id %in% model$flow_direction$from)){
                stop("Hillslope HSU cannot be an outlet or a sink")
            }
            ## specific checks on channel network connectivity
            cfd <- model$flow_direction[ model$flow_direction$from %in% model$channel$id,]
            if(!all(cfd$to %in% model$channel$id)){
                stop("Channel's must flow to channel's")
            }
            ## see if channel routing can be evaluated
            n_chn_con <- table(cfd$from)
            if( any(n_chn_con>1) ){
                warning("Only channel HRUs routing to single channel HRUs are supported by the channel simualtion",
                        "\n",
                        "Channels HRUs with multiple downstream connections are ids: ",
                        paste(names(n_chn_con[n_chn_con>1]),collapse=", "))
                private$info$can_solve_channel <- FALSE
            }else{
                private$info$can_solve_channel <- TRUE
            }
            ## check for circularity in channel network
            is_outlet <- setdiff(paste(model$channel$id),names(n_chn_con)) # no flow from it so outlet
            to_outlet <- setNames(rep(FALSE,length(model$channel$id)),paste(model$channel$id))
            to_outlet[is_outlet] <- TRUE
            chng <- as.integer(is_outlet)
            while(length(chng)>0){
                chng <- unique( cfd$from[cfd$to %in% chng] )
                to_outlet[paste(chng)] <- TRUE
            }
            if(!all(to_outlet)){
                stop("The following channels do not reach an outlet","\n",
                     paste(names(to_outlet)[!to_outlet],serp=", "))
            }
            ## sort flow direction - this is required
            model$flow_direction <-
                model$flow_direction[order(model$flow_direction$from,decreasing=TRUE),]
            
            ## #############################################
            ## Check precip_input
            ## #############################################
            ## check table
            model$precip_input <- private$ftblCheck(model$precip_input,private$model_description$precip_input,
                                            "precip_input")
            ## check it contains onlt valid ids
            if(!all( model$precip_input$id %in% all_hru_id)){
                stop("precip_input table contains id's not in the HRU tables")
            }
            ## check rainfall fractions
            tmp <- tapply(model$precip_input$frc,model$precip_input$id,sum)
            idx <- abs(tmp-1)<delta
            if(!all(idx)){
                stop("PET input fractions for the following HRU id's do not sum to 1: ",
                     paste(names(tmp[!idx]),collapse=", "))
            }
            ## check all non-zero areas receive rainfall
            tmpp <- paste(c(model$channel$id[model$channel$area>0],
                            model$hillslope$id[model$hillslope$area>0]))
            idx <- tmpp %in% names(tmp)
            if(!all(idx)){
                warning(paste0("The following HRUs do not receive PET and have area greater then 0: "),
                        paste(tmpp[!idx],collapse=", "))
            }
            ## sort - not needed but might help speed?
            model$precip_input <-
                model$precip_input[order(model$precip_input$id,decreasing=TRUE),]
            ## set required names
            req_names$data_series[["precip_input"]] <- unique(model$precip_input$name)

            ## #############################################
            ## Check pet_input
            ## #############################################
            ## check table
            model$pet_input <- private$ftblCheck(model$pet_input,private$model_description$pet_input,
                                         "pet_input")
            ## check it contains onlt valid ids
            if(!all( model$pet_input$id %in% all_hru_id)){
                stop("pet_input table contains id's not in the HRU tables")
            }
            ## check series fractions
            tmp <- tapply(model$pet_input$frc,model$pet_input$id,sum)
            idx <- abs(tmp-1)<delta
            if(!all(idx)){
                stop("Precipitation input fractions for the following HRU id's do not sum to 1: ",
                     paste(names(tmp[!idx]),collapse=", "))
            }
            ## check all non-zero areas have some pet
            tmpp <- paste(c(model$channel$id[model$channel$area>0],
                            model$hillslope$id[model$hillslope$area>0]))
            idx <- tmpp %in% names(tmp)
            if(!all(idx)){
                warning(paste0("The following HRUs do not receive precipitation and have area greater then 0: "),
                        paste(ttmp[!idx],collapse=", "))
            }
            ## sort - not needed but might help speed?
            model$pet_input <-
                model$pet_input[order(model$pet_input$id,decreasing=TRUE),]
            ## set required names
            req_names$data_series[["pet_input"]] <- unique(model$pet_input$name)

            ## ##############################################
            ## checks on point_inflow
            ## ##############################################
            ## check table
            model$point_inflow <- private$ftblCheck(model$point_inflow,private$model_description$point_inflow,
                                           "point_inflow")
            ## check goes to channel
            if(nrow(model$point_inflow) > 0){
                idx <- model$point_inflow$id %in% model$channel$id
                if( !all(idx) ){
                    stop(paste("The following point_inflows do not go to channels:",
                               paste(model$point_inflow$name[!idx],collapse=", "),
                               sep="\n"))
                }
            }
            ## set required names
            req_names$data_series[["point_inflow"]] <- unique(model$point_inflow$name)

            ## ##############################################
            ## checks on diffuse_inflow
            ## ##############################################
            ## check table
            model$diffuse_inflow <- private$ftblCheck(model$diffuse_inflow,private$model_description$diffuse_inflow,
                                           "diffuse_inflow")
            ## check goes to channel
            if(nrow(model$diffuse_inflow) > 0){
                idx <- model$diffuse_inflow$id %in% model$channel$id
                if( !all(idx) ){
                    stop(paste("The following diffuse_inflows do not go to channels:",
                               paste(model$diffuse_inflow$name[!idx],collapse=", "),
                               sep="\n"))
                }
            }
            ## set required names
            req_names$data_series[["diffuse_inflow"]] <- unique(model$diffuse_inflow$name)

            ## ##############################################
            ## checks on gauge
            ## ##############################################
            ## check table
            model$gauge <- private$ftblCheck(model$gauge,private$model_description$gauge,
                                           "gauge")
            ## check goes to channel
            if(nrow(model$gauge) > 0){
                idx <- model$gauge$id %in% model$channel$id
                if( !all(idx) ){
                    stop(paste("The following gauges do not go to channels:",
                               paste(model$gauge$name[!idx],collapse=", "),
                               sep="\n"))
                }
            }
            ## set required names
            req_names$output_series[["gauge"]] <- unique(model$gauge$name)

            ## #############################################
            ## check on req_names
            ## #############################################
            ## unpack the required names to vectors
            for(jj in names(req_names)){
                req_names[[jj]] <- do.call(c,req_names[[jj]])
            }
            ## check all output series have unique names
            if( length(req_names$output_names) != length(unique(req_names$output_names)) ){
                stop("All output series should have a unique name")
            }

            ## ########################################################
            ## Checks on maps
            ## ########################################################
            tmp <- sapply(model$map,file.exists)
            if(!all(tmp)){
                warning("The following maps are missing:\n",
                        paste(names(tmp[!tmp]),collapse=", "))
            }

            ## ########################################################
            ## if we have passed all tests and can store the model
            ## But first do things that require reversing in reform model
            ## ########################################################

            ## set D for exp and dexp options
            ## this is used in the c++ as the upper search limit
            model$hillslope$D[ model$hillslope$opt %in% c("exp","dexp") ] <- 1e32

            ## convert hillslope opt into integer values
            tmp <- c("exp"=1,"cnst"=2,"bexp"=3,"dexp"=4)
            model$hillslope$opt <- as.integer(tmp[model$hillslope$opt])
            if(any(is.na(model$hillslope$opt))){
                stop("Failed to convert hillslope options to corresponding integers")
            }

            ## #######################################################
            ## store in private variables
            ## #######################################################
            private$model <- model
            private$info$data_series <- req_names$data_series

            ## return
            invisible( self )

        },
        ## convert the form the internal storage to that input
        ## we presume the model has been checked!!
        reform_model = function(){

            model <- private$model

            ## ######################################################
            ## Udo dirty conversion of the hillslope table values
            ## ######################################################
            ## change class back from integer
            tmp <- c("exp","cnst","bexp","dexp")
            model$hillslope$opt <- tmp[model$hillslope$opt]
            ## change D back to NA
            model$hillslope$D[ model$hillslope$opt %in% c("exp","dexp") ] <- NA

            return(model)
        },
        ## check and add obsservations
        check_obs = function(obs){
            req_series <- private$info$data_series

            ## check types
            if(!is.xts(obs)){ stop("observations should be an xts object") }
            if(!is.vector(req_series) | !all(sapply(req_series,class)=='character') ){ stop("req_series should be a character vector") }

            ## check we have all the series needed
            if( !all( req_series %in% names(obs) ) ){
                stop("Missing input series:",setdiff( req_series , names(obs) ))
            }

            ## check constant time step
            tmp <- diff(as.numeric(index(obs)))
            if( !all( tmp == tmp[1] ) ){
                stop("Time steps in data are not unique")
            }

            ## check all values are finite
            if( !all(is.finite(obs[,req_series])) ){
                stop("There are non finite values in the required time series")
            }
        },
        digest_obs = function(obs){
            ## Assumes all checks passed
            private$time_series$obs_data <- as.matrix(obs)
            private$time_series$index <- index(obs)

            ## set up the column vectors for the rainfall and pet
            private$summary$precip_input <- private$model$precip_input
            private$summary$precip_input$col_idx <- as.integer(
                match(private$model$precip_input$name, colnames(private$time_series$obs)) - 1) # minus 1 to get C++ index
            private$summary$pet_input <- private$model$pet_input
            private$summary$pet_input$col_idx <- as.integer(
                match(private$model$pet_input$name, colnames(private$time_series$obs)) -1 ) # minus 1 to get C++ index
        },
        ## compute the simulation timestep
        comp_ts = function(sub_step=NULL){
            ## work out time steps for use in simulation
            ts <- list()
            ts$step <- diff(as.numeric(private$time_series$index[1:2])) # seconds
            if(is.null(sub_step)){sub_step <- ts$step}
            ts$n_sub_step <- max(1,floor(ts$step/sub_step)) # dimensionless
            ts$sub_step <- ts$step / ts$n_sub_step
            private$info$ts <- ts
        },
        ## ###########################################
        ## Initialise the states
        init_hs = function(tol,max_it){

            dt_init(private$model$hillslope,
                    private$model$channel,
                    private$model$flow_direction,
                    as.double(tol),
                    as.integer(max_it)
                    )
        },
        ## ###############################
        ## function to perform simulations
        sim_hs = function(keep_states,sub_step,tol,max_it,ftol){
          
            ## compute time substep
            if( !is.null(sub_step) && !is.finite(sub_step[1]) ){
                stop("sub_step should be a single finite value")
            }
            ts <- private$comp_ts(sub_step)

            ## work out the courant numbers
            courant <- matrix(as.numeric(NA),length(private$model$hillslope$id),2)
            dt_courant(private$model$hillslope,courant,ts$step,ts$n_sub_step)

            ## Display number of sub steps required
            if( any(courant[,1]>0.7) ){
                warning("Courant number for surface zone is over 0.7\n",
                        "Suggest maximum sub step is: ",
                        round( min( (0.7/courant[,1]) * (ts$step/ts$n_sub_step) ),2 ),
                        "seconds")
            }
            if( any(courant[,2]>0.7) ){
                warning("Courant number for saturated zone is over 0.7\n",
                        "Suggest maximum sub step is: ",
                        round( min( (0.7/courant[,2]) * (ts$step/ts$n_sub_step) ),2 ),
                        "seconds")
            }


            ## Logical if states to be kept and store
            keep_states <- private$time_series$index %in% keep_states
            if(any(keep_states)){
                private$time_series$state_record <- rep(list(as.data.frame(NULL)),length(private$time_series$index))
            }else{
                private$time_series$state_record <- list()
            }

            ## Initialise the mass error store
            private$time_series$mass_balance <- matrix(as.numeric(NA),nrow(private$time_series$obs),6)
            colnames(private$time_series$mass_balance) <-
                c("initial_state","e_t","p","channel_inflow","final_state","error")

            ## set up  channel_inflow
            private$time_series$channel_inflow <- list(
                surface = matrix(as.numeric(NA),
                                 nrow(private$time_series$obs),
                                 length(private$model$channel$id)),
                saturated = matrix(as.numeric(NA),
                                   nrow(private$time_series$obs),
                                   length(private$model$channel$id)))

            colnames(private$time_series$channel_inflow$surface) <-
                colnames(private$time_series$channel_inflow$saturated) <-
                private$model$channel$id
            
            dt_implicit_sim(private$model$hillslope,
                            private$model$channel,
                            private$model$flow_direction,
                            private$summary$precip_input,
                            private$summary$pet_input,
                            private$time_series$obs,
                            private$time_series$channel_inflow$surface,
                            private$time_series$channel_inflow$saturated,
                            private$time_series$mass_balance,
                            as.logical( keep_states ),
                            private$time_series$state_record,
                            ts$step,
                            ts$n_sub_step,
                            as.double(tol),
                            as.integer(max_it),
                            as.double(ftol))
        },
        ## #############################
        init_ch = function(){

            channel <- private$model$channel
            gauge <- private$model$gauge
            point_inflow <- private$model$point_inflow

            ## get the channels upstream of each id
            chn_con <- lapply(channel$id,
                              function(x){
                                  ## work out links going to reach from channels
                                  idx <- (private$model$flow_direction$to==x) &
                                      (private$model$flow_direction$from %in% channel$id)
                                  ## limit to only channels
                                  return( paste(private$model$flow_direction$from[idx]) )
                              }
                              )
            names(chn_con) <- paste(channel$id)


            ## compute the time to travel down each reach
            reach_time <- setNames( channel[,"length"] / channel[,"v_ch"], #private$model$param[channel$v_ch],
                                   channel$id )

            ## storage for output
            linear_time <- setNames(rep(list(NULL),length(gauge$name)),
                                    gauge$name)

            ## Loop gauges
            for(gnm in 1:length(gauge$name)){
                ## initialise diffuse input matrix
                df <- matrix(as.numeric(NA),length(channel$id),2)
                colnames(df) <- c("min_time","max_time")
                rownames(df) <- channel$id
                ## start with location of the gauge
                idx <- paste(gauge$id[gnm])
                df[idx,"min_time"] <- 0
                while(length(idx)>0){
                    ii <- idx[1] ## current reach
                    df[ii,"max_time"] <- df[ii,"min_time"] + reach_time[ii] ## work out max time
                    jdx <- chn_con[[ii]] ## get upstream reaches
                    df[jdx,"min_time"] <- df[ii,"max_time"] ## set min time to max time of downstream
                    idx <- c(idx[-1],jdx)
                }
                ## initialise the point inflow matrix
                ii <- paste(point_inflow$id) ## find id of reach
                pnt <- df[ii,] ## copy rows of df table
                pnt[,"min_time"] <- pnt[,"max_time"] # since point inflow at head of reach
                rownames(pnt) <- point_inflow$name ## rename

                ## trim locations that don't go to that gauge
                df <- df[is.finite(df[,"min_time"])&is.finite(df[,"max_time"]),,drop=FALSE]
                pnt <- pnt[is.finite(pnt[,"min_time"])&is.finite(pnt[,"max_time"]),,drop=FALSE]

                linear_time[[gnm]] <- list(diffuse=df,point=pnt)
            }

            private$summary$channel$linear_time <- linear_time
        },
        sim_ch = function(){
            ## initialise the output
            out <- matrix(NA,length(private$time_series$index),
                          length(private$summary$channel$linear_time))
            colnames(out) <- names(private$summary$channel$linear_time)


            ## ## function to make polynonial representing time delay histogram
            ## ## this works for instananeous flows
            ## fpoly <- function(x){
            ##     if( x["min_time"] == x["max_time"] ){
            ##         ## point input
            ##         ms <- floor(x["max_time"]/private$info$ts$step)+1
            ##         ply <- rep(0,ms)
            ##         ply[ms] <- 1
            ##     }else{
            ##         ms <- ceiling(x["max_time"]/private$info$ts$step)
            ##         #ply <- rep(NA,ms)
            ##         fnsh <- (1:ms)*private$info$ts$step
            ##         strt <- fnsh - private$info$ts$step
            ##         strt <- pmax(strt,x["min_time"])
            ##         fnsh <- pmin(fnsh,x["max_time"])
            ##         ply <- pmax(0,fnsh-strt)
            ##     }
            ##     return(ply/sum(ply))
            ## }
            ## function for point inputs
            fp <- function(x){
                tau0 <- x["min_time"]
                tauL <- x["max_time"]-x["min_time"]
                Dt <- private$info$ts$step
                rL <- floor((tau0+tauL)/Dt)
                irL <- rL+1 ## index in vector since R starts at 1 not zero
                b <- rep(0,irL+1)
                b[irL] <- ((rL+1)*Dt - tau0-tauL)/Dt
                b[irL+1] <- (tau0+tauL - rL*Dt)/Dt
                return(b)
            }
            ## function for diffuse inputs
            fd <- function(x){
                tau0 <- x["min_time"]
                tauL <- x["max_time"]-x["min_time"]
                Dt <- private$info$ts$step
                r0 <- floor(tau0/Dt)
                ir0 <- r0+1 ## index in vector since R starts at 1 not zero
                rL <- floor((tau0+tauL)/Dt)
                irL <- rL+1 ## index in vector since R starts at 1 not zero
                b <- rep(0,irL+1)

                if(rL>r0){
                    b[ir0:(irL-1)] <- Dt # inital values valid unless over written
                    b[ir0] <- ( ((r0+1)*Dt - tau0)^2 ) / (2*Dt)
                    b[ir0+1] <- b[ir0] + ( ((r0+1)*Dt - tau0)*(tau0 - (r0*Dt)) ) / Dt
                    b[irL+1] <- ( (tau0+tauL - (rL*Dt))^2 ) / (2*Dt)
                    b[irL] <- b[irL] + b[irL+1] + ( ((tau0+tauL - (rL*Dt))*((rL+1)*Dt - tau0-tauL)) / Dt ) ## added to self since rL could equal r0+1
                    if( rL > (r0+1) ){
                        b[ir0+1] <- b[ir0+1] + Dt/2
                        b[irL] <- b[irL] + Dt/2
                    }
                }else{
                    b[ir0] <- (tauL/Dt)*( (r0+1)*Dt - tau0 - (tauL/2) )
                    b[ir0+1] <- (tauL/Dt)*(tau0 + (tauL/2) - r0*Dt )
                }

                return(b/sum(b)) #this should really be b*v_ch/L
            }

            ## Loop gauges
            for(gnm in names(private$summary$channel$linear_time)){

                ## initialise the point - set to 0
                out[,gnm] <- 0

                ## diffuse inputs
                df <- private$summary$channel$linear_time[[gnm]]$diffuse

                for(ii in rownames(df)){
                    ply <- fd(df[ii,]) ##fpoly(df[ii,])
                    ## compute input
                    x <- private$time_series$channel_inflow$surface[,ii] +
                        private$time_series$channel_inflow$saturated[,ii]
                    ## add diffuse inputs to x
                    idx <- paste(private$model$diffuse_inflow$id)==ii
                    nm <- private$model$diffuse_inflow$name[idx]
                    x <- x + rowSums(private$time_series$obs[,nm])
                    ## pas and filter
                    np <- length(ply)-1
                    x <- c(rep(x[1],np),x)
                    q <- filter(x,ply,method="conv",sides=1)
                    if(np>0){q <- q[-(1:np)]}
                    out[,gnm] <- out[,gnm] + q
                }

                ## handle point inputs
                pnt <- private$summary$channel$linear_time[[gnm]]$point
                ## loop point inputs upstream
                for(ii in rownames(pnt)){
                    ply <- fp(pnt[ii,]) ##fpoly(pnt[ii,])
                    np <- length(ply)-1
                    ## compute input
                    x <- private$time_series$obs[,ii]
                    x <- c(rep(x[1],np),x)
                    ## apply polynomial
                    q <- filter(x,ply,method="conv",sides=1)
                    if(np>0){q <- q[-(1:np)]}
                    out[,gnm] <- out[,gnm] + q
                }
            }

            private$time_series$gauge_flow <- out
        }
    )
    )