Fixed merge conflict from develop

.gitignore

@@ -3,19 +3,14 @@
 .RData
 .Ruserdata
 
-NAMESPACE
-
 *.o
 src/onlineforecast\.so
-inst/doc
 
 modifications_old_notstaged/
 
 
 cache/
 
-man/
-
 misc-R/*cache*
 
 vignettes/*cache*

DESCRIPTION

 Package: onlineforecast
 Type: Package
 Title: Forecast Modelling for Online Applications
-Version: 0.9.3
+Version: 0.10.1
 Description: A framework for fitting adaptive forecasting models. Provides a way to use forecasts as input to models, e.g. weather forecasts for energy related forecasting. The models can be fitted recursively and can easily be setup for updating parameters when new data arrives. See the included vignettes, the website <https://onlineforecasting.org> and the paper "Short-term heat load forecasting for single family houses" <doi:10.1016/j.enbuild.2013.04.022>.
 License: GPL-3
 Encoding: UTF-8
@@ -21,10 +21,11 @@ Suggests:
     knitr,
     rmarkdown,
     R.rsp,
-    testthat (>= 2.1.0),
+    testthat (>= 3.0.0),
     data.table,
     plotly
 VignetteBuilder: knitr
 RoxygenNote: 7.1.1
 URL: https://onlineforecasting.org
 BugReports: https://lab.compute.dtu.dk/packages/onlineforecast/-/issues
+Config/testthat/edition: 3

NAMESPACE

+# Generated by roxygen2: do not edit by hand
+
+S3method("==",data.list)
+S3method(as.data.frame,data.list)
+S3method(as.data.list,data.frame)
+S3method(aslt,POSIXct)
+S3method(aslt,POSIXlt)
+S3method(aslt,character)
+S3method(aslt,numeric)
+S3method(check,data.list)
+S3method(complete_cases,data.frame)
+S3method(complete_cases,list)
+S3method(ct,POSIXct)
+S3method(ct,POSIXlt)
+S3method(ct,character)
+S3method(ct,numeric)
+S3method(lagdf,character)
+S3method(lagdf,data.frame)
+S3method(lagdf,factor)
+S3method(lagdf,logical)
+S3method(lagdf,matrix)
+S3method(lagdf,numeric)
+S3method(pairs,data.list)
+S3method(plot_ts,data.frame)
+S3method(plot_ts,data.list)
+S3method(plot_ts,matrix)
+S3method(plot_ts,rls_fit)
+S3method(plotly_ts,data.frame)
+S3method(plotly_ts,data.list)
+S3method(print,forecastmodel)
+S3method(resample,data.frame)
+S3method(residuals,data.frame)
+S3method(residuals,forecastmodel_fit)
+S3method(residuals,list)
+S3method(residuals,matrix)
+S3method(score,data.frame)
+S3method(score,list)
+S3method(subset,data.list)
+S3method(summary,rls_fit)
+export("%**%")
+export("nams<-")
+export(AR)
+export(as.data.list)
+export(aslt)
+export(bspline)
+export(cache_name)
+export(cache_save)
+export(check)
+export(complete_cases)
+export(ct)
+export(data.list)
+export(forecastmodel)
+export(fs)
+export(getse)
+export(gof)
+export(in_range)
+export(lagdf)
+export(lagdl)
+export(lagvec)
+export(lapply_cbind)
+export(lapply_cbind_df)
+export(lapply_rbind)
+export(lapply_rbind_df)
+export(lm_fit)
+export(lm_optim)
+export(lm_predict)
+export(long_format)
+export(lp)
+export(make_input)
+export(make_periodic)
+export(make_tday)
+export(nams)
+export(one)
+export(par_ts)
+export(pbspline)
+export(persistence)
+export(plot_ts)
+export(plotly_ts)
+export(pst)
+export(resample)
+export(rls_fit)
+export(rls_optim)
+export(rls_predict)
+export(rls_prm)
+export(rls_summary)
+export(rls_update)
+export(rmse)
+export(score)
+export(setpar)
+export(stairs)
+export(step_optim)
+importFrom(Rcpp,sourceCpp)
+importFrom(grDevices,colorRampPalette)
+importFrom(grDevices,graphics.off)
+importFrom(graphics,axis)
+importFrom(graphics,axis.POSIXct)
+importFrom(graphics,legend)
+importFrom(graphics,lines)
+importFrom(graphics,mtext)
+importFrom(graphics,pairs)
+importFrom(graphics,panel.smooth)
+importFrom(graphics,par)
+importFrom(graphics,plot)
+importFrom(graphics,points)
+importFrom(graphics,title)
+importFrom(parallel,mclapply)
+importFrom(stats,aggregate)
+importFrom(stats,complete.cases)
+importFrom(stats,lm)
+importFrom(stats,optim)
+importFrom(stats,predict)
+importFrom(stats,residuals)
+importFrom(stats,sd)
+useDynLib(onlineforecast)

R/as.data.list.R

@@ -22,16 +22,14 @@
 #' @param object The object to be converted into a data.list
 #' @return a value of class data.list
 #' @seealso \code{For specific detailed info see the children, e.g. \link{as.data.list.data.frame} }
-#' @family as.data.list
 #' 
+#' @rdname as.data.list
 #' @export
 as.data.list <- function(object){
     UseMethod("as.data.list")
 }
 
 
-
-
 #' Convert a data.frame into a data.list
 #'
 #' The convention is that columns with forecasts are postfixed with \code{.kxx} where
@@ -41,7 +39,6 @@ as.data.list <- function(object){
 #' @param object The data.frame to be converted.
 #' @return a data.list
 #' @seealso as.data.list
-#' @family as.data.list
 #' @examples
 #' # Convert a dataframe with time and two observed variables
 #' X <- data.frame(t=1:10, x=1:10, y=1:10)
@@ -55,7 +52,9 @@ as.data.list <- function(object){
 #' X
 #' as.data.frame(as.data.list(X))
 #'
+#' @rdname as.data.list
 #' @export
+
 as.data.list.data.frame <- function(object) {
     X <- object
     #TEST

R/bspline.R

@@ -68,7 +68,7 @@
 #' @export
 bspline <- function(X, Boundary.knots = NA, intercept = FALSE, df = NULL, knots = NULL, degree = 3, bknots = NA, periodic = FALSE) {
     # bknots is just a short for Boundary.knots and replace if Boundary.knots are not given.
-    if(is.na(Boundary.knots)){
+    if(is.na(Boundary.knots[1])){
         Boundary.knots <- bknots
     }
     # If a list, then call on each element
@@ -81,9 +81,12 @@ bspline <- function(X, Boundary.knots = NA, intercept = FALSE, df = NULL, knots
         nams(val) <- nams(X)
         return(val)
     }
-    # X is a data.frame or matrix
+    # X must be a data.frame or matrix
+    if(!(class(X) %in% c("data.frame","matrix"))){ stop("X must be a data.frame or matrix") }
     # First find the horizons, they are used in the end
     nms <- nams(X)
+    # All columns must be like "k12"
+    #if(length(grep("[k|h][[:digit:]]+$", nms)) != length(nms)){ stop("All column names must indicate a horizon, so start with k and then an integer") }
     # Run for each horizon and calculate the basis splines
     L <- lapply(1:ncol(X), function(i) {
       if (is.na(Boundary.knots[1])) {

R/cache_name.R

@@ -40,16 +40,17 @@
 #' # In real use this should not be temporary! (changes between R sessions with tempdir())
 #' cachedir <- tempdir()
 #'
+#' # Uncomment to run:
 #' # First time it takes at least 1 sec.
-#' fun(x=2,y=2)
+#' #fun(x=2,y=2)
 #' # Second time it loads the cache and is much faster
-#' fun(x=2,y=2)
+#' #fun(x=2,y=2)
 #' # Try changing the arguments (x,y) and run again
 #'
 #' # See the cache file(s)
-#' dir(cachedir)
+#' #dir(cachedir)
 #' # Delete the cache folder
-#' unlink(cachedir, recursive=TRUE)
+#' #unlink(cachedir, recursive=TRUE)
 #'
 #' # Demonstrate how cache_name() is functioning
 #' # Cache using the all objects given in the function calling, i.e. both x and y
@@ -81,7 +82,18 @@
 cache_name <- function(..., cachedir = "cache"){
     # Get the name, definition and arguments of the function from which cache_name was called
     funname <- strsplit(deparse(sys.calls()[[sys.nframe()-1]]), "\\(")[[1]][1]
-    fundef <- digest::digest(attr(eval(parse(text=funname)), "srcref"))
+    # Find the function in the nearest environment in the stack (i.e. parent calls)
+    ## for(i in rev(sys.parents())){
+    ##     browser()
+    ##     if(funname %in% ls(parent.frame(i+1))){
+    ##         val <- mget(funname, parent.frame(i+1))
+    ##         break
+    ##     }
+    ## }
+    ## fundef <- digest::digest(attr(eval(val[[funname]]), "srcref"))
+    # Somehow the above stopped working, don't know why! just take it, this should do the same I guess
+    fundef <- try(get(funname), silent=TRUE)
+    fundef <- digest::digest(fundef)
     # if no arguments were given, then use the arguments function from which cache_name was called
     if(length(list(...)) == 0){
         funargs <- digest::digest(as.list( match.call(definition = sys.function( -1 ), call = sys.call(-1)))[-1])

R/complete_cases.R

+#' Returns a logical vector indicating the time points which 
+#'
+#' Given a forecast matrix the forecasts are lagged "+k" steps to align them and
+#' then 'complete.cases()' is run on that .
+#'
+#' Gieven a list of forecast matrices the points where all are complete (also all horizons) are complete are TRUE.
+#' 
+#' @title Find complete cases in forecast matrices
+#' @param object A data.frame (with columns named 'kxx') or a list of
+#'     data.frames.
+#' @param kseq integer vector: If given then only these horizons are processed.
+#' @return A logical vector specifying if there is no missing
+#'     values across all horizonsd.
+#' @author Peder Bacher
+#'
+#' @examples
+#' # Take a small data set
+#' D <- subset(Dbuilding, 1:20, kseq=1:5)
+#' # Check the forecast matrix of ambient temperature
+#' D$Ta
+#' # Which are complete over all horizons? The first are not since not all horizons
+#' # have a value there (after lagging)
+#' complete_cases(D$Ta)
+#' # Same goes if given as a list
+#' complete_cases(D["Ta"])
+#' # and if more than one is given
+#' complete_cases(D[c("Ta","I")])
+#' 
+#' # Set some NA of some horizon
+#' D$I$k3[8:9] <- NA
+#' # Now they are recognized as not complete
+#' complete_cases(D[c("Ta","I")])
+#'
+#' # If we deal with residuals, which are observations and there for have column names "hxx"
+#' Resid <- residuals(D$Ta, D$Taobs)
+#' names(Resid)
+#' # With columns with "h" instead of "k" no lagging occurs in complete_cases
+#' complete_cases(Resid)
+#' #
+#' Resid2 <- Resid
+#' Resid$h3[8:9] <- NA
+#' complete_cases(list(Resid,Resid2))
+#' 
+#' @rdname complete_cases
+#' @importFrom stats complete.cases
+#' @export
+complete_cases <- function(object, kseq = NA){
+    UseMethod("complete_cases")
+}
+
+
+#' @rdname complete_cases
+#' @importFrom stats complete.cases
+#' @export
+complete_cases.list <- function(object, kseq=NA){
+    if(length(grep("[k][[:digit:]]+$", nams(object[[1]])))){
+        prefix <- "k"
+    }else{
+        prefix <- "h"
+    }
+    # Do it only for the given kseq horizons, if kseq is not set, then take it from the first
+    if(is.na(kseq[1])){
+        kseq <- as.integer(gsub(prefix, "", nams(object[[1]])))
+    }
+    # Check that they all have kseq horizons
+    # Init a logical matrix with each point
+    ok <- matrix(TRUE, nrow(object[[1]]), length(kseq))
+    #
+    for(i in 1:length(object)){
+        if(!all(kseq %in% as.integer(gsub(prefix, "", nams(object[[i]]))))){
+            stop("Not all forecast matrices have kseq horizons: ",kseq)
+        }
+        ok <- ok & !is.na(object[[i]][ ,pst(prefix,kseq)])
+    }
+    #
+    ok <- as.data.frame(ok)
+    if(prefix == "k"){
+        # Lag to match resiuduals in time
+        names(ok) <- pst("k",kseq)
+        ok <- lagdf(ok, "+k")
+    }
+    # Finally, the vector with TRUE for all points with no NAs for any forecast
+    ok <- apply(ok, 1, all)
+    # Just set NAs to false
+    ok[is.na(ok)] <- FALSE
+    #
+    return(ok)
+}
+
+
+#' @rdname complete_cases
+#' @importFrom stats complete.cases
+#' @export
+complete_cases.data.frame <- function(object, kseq=NA){
+    # Make a distinction between "kxx" which are forecasts or "hxx", which are observations
+    if(length(grep("[k][[:digit:]]+$", nams(object)[1]))){
+        # Forecasts, so lag them
+        object <- lagdf(object, "+k")        
+    }
+    complete.cases(object)
+}

R/data.list.R

@@ -117,10 +117,10 @@ subset.data.list <- function(x, subset = NA, nms = NA, kseq = NA, lagforecasts =
             X <- D[[nms[i]]]
             if(class(X)[1] == "data.frame" ){
                 # Check if holds forecasts by checking if any name is "kxx"
-                if(length(grep("^k[[:digit:]]+$", names(X))) > 0){
+                if(length(grep("k[[:digit:]]+$", names(X))) > 0){
                     # If it holds forecasts, check that they are all there
                     if( !all(pst("k",kseq) %in% names(X)) ){
-                        warning(pst("The variable ",nms[i]," contain ",pst(names(X),collapse=",")," hence doesn't contain all k in kseq = ",pst(kseq,collapse=",")))
+                        warning(pst("The variable ",nms[i]," contains ",pst(names(X),collapse=",")," hence doesn't contain all k in kseq = ",pst(kseq,collapse=",")))
                     }
                 }
             }
@@ -160,7 +160,7 @@ subset.data.list <- function(x, subset = NA, nms = NA, kseq = NA, lagforecasts =
         val <- lapply(D[nms], function(X) {
             if (any(class(X) == "data.frame")) {
                 # Check if holds forecasts by checking if any name is "kxx"
-                if(length(grep("^k[[:digit:]]+$", names(X))) > 0){
+                if(length(grep("k[[:digit:]]+$", names(X))) > 0){
                     return(X[subset,pst("k",kseq), drop=FALSE])
                 }else{
                     return(X[subset, , drop=FALSE])
@@ -173,7 +173,7 @@ subset.data.list <- function(x, subset = NA, nms = NA, kseq = NA, lagforecasts =
     # Lag the forecasts k if specified
     if(lagforecasts){
         val <- lapply(val, function(X){
-            if(any(class(X) == "data.frame") & length(grep("^k[[:digit:]]+$",names(X))) > 0) {
+            if(any(class(X) == "data.frame") & length(grep("k[[:digit:]]+$",names(X))) > 0) {
                 return(lagdf.data.frame(X, lagseq="+k"))
             }else{
                 return(X)
@@ -213,7 +213,7 @@ as.data.frame.data.list <- function(x, row.names=NULL, optional=FALSE, ...){
         val <- as.data.frame(val)
     }
     # Fix names of data.frames (i.e. forecasts, their names are now "kxx", but should be X.kxx)
-    i <- grep("^k[[:digit:]]+$", names(val))
+    i <- grep("k[[:digit:]]+$", names(val))
     if(length(i) > 0){
         names(val)[i] <- pst(names(x)[i],".",names(val)[i])
     }
@@ -330,8 +330,11 @@ check.data.list <- function(object){
     # Which is data.frame or matrix?
     dfOrMat <- sapply(D, function(x){ (class(x) %in% c("matrix","data.frame"))[1] })
     # Vectors check
-    vecchecks <- c("ok","NAs","length","class")
     vecseq <- which(!dfOrMat & names(dfOrMat) != "t")
+    Observations <- NA
+    if(length(vecseq) > 0){
+        cat("Observation vectors:\n")
+        vecchecks <- c("ok","NAs","length","class")
         Observations <- data.frame(matrix("", nrow=length(vecseq), ncol=length(vecchecks), dimnames=list(names(vecseq),vecchecks)), stringsAsFactors=FALSE)
         Observations$ok <- "V"
         #
@@ -352,9 +355,14 @@ check.data.list <- function(object){
                 Observations$ok[i] <- ""
             }
         }
+        print(Observations)
+    }
     #
     # For forecasts
     dfseq <- which(dfOrMat)
+    Forecasts <- NA
+    if(length(dfseq) > 0){
+        cat("\nForecast data.frames or matrices:\n")
         dfchecks <- c("ok","maxNAs","meanNAs","nrow","colnames","sameclass","class")
         Forecasts <- data.frame(matrix("", nrow=length(dfseq), ncol=length(dfchecks), dimnames=list(names(dfseq),dfchecks)), stringsAsFactors=FALSE)
         Forecasts$ok <- "V"
@@ -374,7 +382,7 @@ check.data.list <- function(object){
                 Forecasts$nrow[i] <- nrow(D[[nm]])
             }
             # Check the colnames, are they unique and all k+integer?
-        if(!length(unique(grep("^k[[:digit:]]+$",colnms,value=TRUE))) == length(colnms)){
+            if(!length(unique(grep("k[[:digit:]]+$",colnms,value=TRUE))) == length(colnms)){
                 Forecasts$colnames[i] <- "X"
             }
             if(!length(unique(sapply(colnms, function(colnm){ class(D[[nm]][ ,colnm]) }))) == 1){
@@ -387,11 +395,8 @@ check.data.list <- function(object){
                 Forecasts$ok[i] <- ""
             }
         }
-    #
-    message("Observation vectors:")
-    print(Observations)
-    message("\nForecast data.frames or matrices:")
         print(Forecasts)
+    }
 
     invisible(list(Observations=Observations, Forecasts=Forecasts))
 }

R/depth.R

+#' Depth of a list
+#'
+#' Returns the depth of a list
+#' @title Depth of a list
+#' @param this list
+#' @return integer
+depth <- function(this) ifelse(is.list(this), 1L + max(sapply(this, depth)), 0L)

R/flattenlist.R

+#' Flattens list in a single list of data.frames
+#'
+#' Flattens list. Can maybe be made better. It might end up copying data in
+#' memory!? It might change the order of the elements.
+#' @title Flattens list
+#' @param x List to flatten.
+#' @return A flatten list
+flattenlist <- function(x){
+    (n <- depth(x))
+    if(n == 2){
+        # Its fine
+        return(x)
+    }else if(n ==3){
+        unlist(x, recursive=FALSE)
+    }else{
+        morelists <- sapply(x, function(xprime) class(xprime)[1]=="list")
+        out <- c(x[!morelists], unlist(x[morelists], recursive=FALSE))
+        if(sum(morelists)){ 
+            Recall(out)
+        }else{
+            return(out)
+        }
+    }
+}

R/forecastmodel.R

@@ -23,7 +23,9 @@ forecastmodel <- R6::R6Class("forecastmodel", public = list(
     #
     # The horizons to fit for
     kseq = NA,
-    # The (transformation stage) parameters used for the fit
+    # The horizons to optimize for
+    kseqopt = NA,
+    # The (transformation stage) parameters (only the ones set in last call of insert_prm())
     prm = NA,
     # Stores the maximum lag for AR terms
     maxlagAR = NA,
@@ -83,7 +85,7 @@ forecastmodel <- R6::R6Class("forecastmodel", public = list(
 
 
     #----------------------------------------------------------------
-    # Get the transformation parameters
+    # Get the transformation parameters (set for optimization)
     get_prmbounds = function(nm){
         if(nm == "init"){
             if(is.null(dim(self$prmbounds))){
@@ -130,8 +132,15 @@ forecastmodel <- R6::R6Class("forecastmodel", public = list(
         }
         # MUST INCLUDE SOME checks here and print useful messages if something is not right
         if(any(is.na(prm))){ stop(pst("None of the parameters (in prm) must be NA: prm=",prm)) }
-
-        # Keep the prm
+        # If given without names, then set in same order as in the prm_bounds
+        if(is.null(nams(prm))){
+            if(length(prm) != nrow(self$prmbounds)){
+                stop("prm was given without names and length not same as prmbounds, so don't know what to do")
+            }else{
+                nams(prm) <- row.names(self$prmbounds)
+            }
+        }
+        # Keep the prm given
         self$prm <- prm
         # Find if any opt parameters, first the one with "__" hence for the inputs
         pinputs <- prm[grep("__",nams(prm))]
@@ -152,7 +161,7 @@ forecastmodel <- R6::R6Class("forecastmodel", public = list(
                     # Find if the input i have prefix match with the opt. parameter ii
                     if(pnms[ii]==nams(self$inputs)[i]){
                         # if the opt. parameter is in the expr, then replace
-                        self$inputs[[i]]$expr <- private$replace_value(name = pprm[ii],
+                        self$inputs[[i]]$expr <- private$replace_prmvalue(name = pprm[ii],
                                                                        value = pinputs[ii],
                                                                        expr = self$inputs[[i]]$expr)
                     }
@@ -160,12 +169,12 @@ forecastmodel <- R6::R6Class("forecastmodel", public = list(
             }
         }
         # ################
-        # For the fit parameters, insert from prm if any found
+        # For the regression parameters, insert from prm if any found
         if (length(preg) & any(!is.na(self$regprmexpr))) {
             nams(preg)
             for(i in 1:length(preg)){
                 # if the opt. parameter is in the expr, then replace
-                self$regprmexpr <- private$replace_value(name = nams(preg)[i],
+                self$regprmexpr <- private$replace_prmvalue(name = nams(preg)[i],
                                                          value = preg[i],
                                                          expr = self$regprmexpr)
             }
@@ -175,6 +184,32 @@ forecastmodel <- R6::R6Class("forecastmodel", public = list(
     #----------------------------------------------------------------
 
 
+    #----------------------------------------------------------------
+    # Return the values of the parameter names given
+    get_prmvalues = function(prmnames){
+        #
+        regprm <- eval(parse(text = self$regprmexpr))
+        # From the input parameters
+        val <- sapply(prmnames, function(nm){
+            if(length(grep("__",nm))){
+                tmp <- strsplit(nm, "__")[[1]]
+                if(tmp[1] %in% names(self$inputs)){
+                    return(as.numeric(private$get_exprprmvalue(tmp[2], self$inputs[[tmp[1]]]$expr)))
+                }else{
+                    return(NA)
+                }
+            }else{
+                if(nm %in% names(regprm)){
+                    return(as.numeric(regprm[nm]))
+                }else{
+                    return(NA)
+                }
+            }
+        })
+        return(val)
+    },
+    #----------------------------------------------------------------
+
     #----------------------------------------------------------------
     # Function for transforming the input data to the regression data
     transform_data = function(data){
@@ -187,11 +222,11 @@ forecastmodel <- R6::R6Class("forecastmodel", public = list(
             if(class(L)[1]=="data.frame"){ return(list(L)) }
             if(class(L)[1]!="list"){ stop(pst("The value returned from evaluating: ",input$expr,", was not a matrix, data.frame or a list of them."))}
             if(class(L[[1]])[1]=="matrix"){ return(lapply(L, function(mat){ return(as.data.frame(mat)) })) }
-            return(L)
+            return(flattenlist(L))
         })
-        # Put together in one data.list
-        L <- structure(do.call(c, L), class="data.list")
-        #
+        # Make it a data.list with no subsubelements (it's maybe not a data.list, since it miss "t", however to take subsets etc., it must be a data.list)
+        L <- flattenlist(L)
+        class(L) <- "data.list"
         return(L)
     },
     #----------------------------------------------------------------
@@ -257,10 +292,12 @@ forecastmodel <- R6::R6Class("forecastmodel", public = list(
                         stop("The input variable '",nm,"' doesn't have the same number of observations as time vector 't' in the data. It has ",length(data[[nm]]),", but 't' has ",length(data$t))
                     }
                 }else{
+                    if(!nm == "pi"){
                         stop("The variable '",nm,"' is missing in data, or it has the wrong class.\nIt must be class: data.frame, matrix or vector.\nIt is needed for the input expression '",self$inputs[[i]]$expr[[1]],"'")
                     }
                 }
             }
+        }
     },
 
     #----------------------------------------------------------------
@@ -289,7 +326,7 @@ forecastmodel <- R6::R6Class("forecastmodel", public = list(
 
     #----------------------------------------------------------------
     # Replace the value in "name=value" in expr
-    replace_value = function(name, value, expr){
+    replace_prmvalue = function(name, value, expr){
         # First make regex
         pattern <- gsub("\\.", ".*", name)
         # Try to find it in the input
@@ -298,7 +335,7 @@ forecastmodel <- R6::R6Class("forecastmodel", public = list(
         if(pos>0){
             pos <- c(pos+attr(pos,"match.length"))
             # Find the substr to replace with the prm value
-            (tmp <- substr(expr, pos, nchar(expr)))
+            tmp <- substr(expr, pos, nchar(expr))
             pos2 <- regexpr(",|)", tmp)
             # Insert the prm value and return
             expr <- pst(substr(expr,1,pos-1), "=", value, substr(expr,pos+pos2-1,nchar(expr)))
@@ -309,6 +346,30 @@ forecastmodel <- R6::R6Class("forecastmodel", public = list(
     },
     #----------------------------------------------------------------
 
+    #----------------------------------------------------------------
+    get_exprprmvalue = function(name, expr){
+        #name <- "degree"
+        #expr <- "bspline(tday, Boundary.knots = c(start=6,18), degree = 5, intercept=TRUE) %**% ones() + 2 + ones()"
+        #expr <- "one()"
+        expr <- gsub(" ", "", expr)
+        
+        # First make regex
+        pattern <- gsub("\\.", ".*", name)
+        # Try to find it in the input
+        pos <- regexpr(pattern, expr)
+        # Only replace if prm was found
+        if(pos>0){
+            pos <- c(pos+attr(pos,"match.length"))
+            # Find the substr to replace with the prm value
+            (tmp <- substr(expr, pos, nchar(expr)))
+            pos2 <- regexpr(",|)", tmp)
+            return(substr(tmp, 2, pos2-1))
+        }else{
+            return(NA)
+        }
+    },
+    #----------------------------------------------------------------
+    
     #----------------------------------------------------------------
     # For deep cloning, in order to get the inputs list of R6 objects copied
     deep_clone = function(name, value) {
@@ -344,14 +405,16 @@ print.forecastmodel <- function(x, ...){
     model <- x
     #    cat("\nObject of class forecastmodel (R6::class)\n\n")
     cat("\nOutput:",model$output)
-    cat("Inputs: ")
+    cat("\nInputs: ")
     if(length(model$inputs) == 0 ){
-        cat("No inputs\n")
+        cat("\nNo inputs\n\n")
     }else{
         cat(names(model$inputs)[1],"=",model$inputs[[1]]$expr,"\n")
+        if(length(model$inputs) > 1){
             for(i in 2:length(model$inputs)){
                 cat("       ",names(model$inputs)[i],"=",model$inputs[[i]]$expr,"\n")
             }
+        }
         cat("\n")
     }
 }

R/forecastmodel.R-documentation.R

@@ -48,6 +48,8 @@
 #'
 #'     - kseq = NA: The horizons to fit for.
 #'
+#'     - kseqopt = NA: The horizons to fit for when optimizing.
+#' 
 #'     - p = NA: The (transformation stage) parameters used for the fit.
 #'
 #'     - Lfits = list(): The regression fits, one for each k in kseq (simply a list with the latest fit).
@@ -135,7 +137,7 @@
 
 #----------------------------------------------------------------
 #' @section \code{$insert_prm(prm)}:
-#' Insert the transformation parameters prm in the input expressions and regression expressions, and keep them (simply string manipulation).
+#' Insert the transformation parameters prm in the input expressions and regression expressions, and keep them in $prm (simply string manipulation).
 #'
 #' @examples
 #' 

R/in_range.R

@@ -43,15 +43,30 @@
 #' # since it's covering to "2010-12-25 23:00:00" to "2010-12-26 00:00:00"
 #' D$t[in_range("2010-12-26", D$t, "2010-12-27")]
 #'
+#' # When characters are given, then they are translated to the time zone of the time vector
+#' D <- subset(Dbuilding, c("2010-12-15", "2010-12-16"))
+#' D$t <- ct(D$t, tz="CET")
+#' D$t[in_range("2010-12-15 02:00", D$t, "2010-12-15 05:00")]
 #'
 #' @export
 
 in_range <- function(tstart, time, tend=NA) {
-    if (class(tstart)[1] == "character") 
-        tstart <- ct(tstart)
-    if (is.na(tend))
+    if (is.na(tend)){
         tend <- time[length(time)]
-    if (class(tend)[1] == "character") 
-        tend <- ct(tend)
-    ct(tstart) < time & time <= ct(tend)
+    }
+    if(is.numeric(time)){
+        return(tstart < time & time <= tend)
+    }else{
+        # Use the time zone of the time when converting from character
+        tz <- attr(time, "tzone")
+        #
+        if (class(tstart)[1] == "character"){
+            tstart <- ct(tstart, tz=tz)
+        }
+        if (class(tend)[1] == "character"){
+            tend <- ct(tend, tz=tz)
+        }
+        #timezone <- attr(time, "tzone")
+        return(tstart < time & time <= tend)
+    }
 }

R/lagdf.R

@@ -22,9 +22,9 @@
 #' lagdf(1:10, 3)
 #' # Back in time
 #' lagdf(1:10, -3)
-#' # Works but returns a numric
+#' # Works but returns a numeric column
 #' lagdf(as.factor(1:10), 3)
-#' # Works and returns a character
+#' # Works and returns a character column
 #' lagdf(as.character(1:10), 3)
 #' # Giving several lag values
 #' lagdf(1:10, c(1:3))
@@ -114,12 +114,12 @@ lagdf.data.frame <- function(x, lagseq) {
         if (lagseq %in% c("+k","+h")) {
             lagseq <- rep(0, length(nms))
             ## lagseq according to the k value of the columnnames
-            i <- grep("^[k|h][[:digit:]]+$", nms)
+            i <- grep("[k|h][[:digit:]]+$", nms)
             lagseq[i] <- as.integer(sapply(strsplit(nms[i], "[k|h]"), function(x){ x[length(x)] }))
         } else if (lagseq %in% c("-k","-h")) {
             lagseq <- rep(0, length(nms))
             ## lagseq according to the negative k value of the columnnames
-            i <- grep("^[k|h][[:digit:]]+$", nms)
+            i <- grep("[k|h][[:digit:]]+$", nms)
             lagseq[i] <- -as.integer(sapply(strsplit(nms[i], "[k|h]"), function(x){ x[length(x)] }))
         }
     }

R/lagdl.R

+## Do this in a separate file to see the generated help:
+#library(devtools)
+#document()
+#load_all(as.package("../../onlineforecast"))
+#?lagdl
+
+#' Lagging by shifting the values back or fourth always returning a data.list.
+#'
+#' This function lags (shifts) the values of the vector. A data.list is always returned with each data.frame lagged with \code{lagdf}.
+#'
+#' 
+#' @title Lagging which returns a data.list
+#' @param DL The data.list to be lagged.
+#' @param lagseq The integer(s) setting the lag steps.
+#' @return A data.list.
+#' @examples
+#' # The values are simply shifted in each data.frame with lagdf
+#'
+#'@export
+
+lagdl <- function(DL, lagseq){
+    iseq <- which(sapply(DL,class) %in% c("data.frame","matrix"))
+    DL[iseq] <- lapply(iseq, function(i){
+        lagdf(DL[[i]], lagseq)
+    })
+    return(DL)
+}

R/lapply.R

@@ -7,36 +7,40 @@
 #' Helper which does lapply and then cbind
 #' @param X object to apply on
 #' @param FUN function to apply
+#' @param ... passed on to lapply
 #' @export
-lapply_cbind <- function(X, FUN){
-  val <- lapply(X, FUN)
+lapply_cbind <- function(X, FUN, ...){
+  val <- lapply(X, FUN, ...)
   return(do.call("cbind", val))
 }
 
 #' Helper which does lapply and then rbind
 #' @param X object to apply on
 #' @param FUN function to apply
+#' @param ... passed on to lapply
 #' @export
-lapply_rbind <- function(X, FUN){
-  val <- lapply(X, FUN)
+lapply_rbind <- function(X, FUN, ...){
+  val <- lapply(X, FUN, ...)
   return(do.call("rbind", val))
 }
 
 #' Helper which does lapply, cbind and then as.data.frame
 #' @param X object to apply on
 #' @param FUN function to apply
+#' @param ... passed on to lapply
 #' @export
-lapply_cbind_df <- function(X, FUN){
-  val <- lapply(X, FUN)
+lapply_cbind_df <- function(X, FUN, ...){
+  val <- lapply(X, FUN, ...)
   return(as.data.frame(do.call("cbind", val)))
 }
 
 #' Helper which does lapply, rbind and then as.data.frame
 #' @param X object to apply on
 #' @param FUN function to apply
+#' @param ... passed on to lapply
 #' @export
-lapply_rbind_df <- function(X, FUN){
-  val <- lapply(X, FUN)
+lapply_rbind_df <- function(X, FUN, ...){
+  val <- lapply(X, FUN, ...)
   return(as.data.frame(do.call("rbind", val)))
 }
 

R/lm_fit.R

@@ -37,9 +37,8 @@
 #' # Define a simple model 
 #' model <- forecastmodel$new()
 #' model$output <- "y"
-#' model$add_inputs(Ta = "Ta",
+#' model$add_inputs(Ta = "lp(Ta, a1=0.9)",
 #'                  mu = "one()")
-#' model$add_regprm("rls_prm(lambda=0.99)")
 #'
 #' # Before fitting the model, define which points to include in the evaluation of the score function
 #' D$scoreperiod <- in_range("2010-12-20", D$t)
@@ -53,13 +52,13 @@
 #' class(fit)
 #' # The one-step forecast
 #' plot(D$y, type="l")
-#' lines(fit$Yhat$k1, col=2)
+#' lines(lagvec(fit$Yhat$k1,-1), col=2)
 #' # Get the residuals
 #' plot(residuals(fit)$h1)
 #' # Score for each horizon
-#' score_fit(fit)
+#' score(residuals(fit))
 #'
-#' # The lm_fit don't put anything in
+#' # The lm_fit don't put anything in this field
 #' fit$Lfitval
 #' # Find the lm fits here
 #' model$Lfits
@@ -68,6 +67,10 @@
 #' # Some diurnal pattern is present
 #' acf(residuals(fit)$h1, na.action=na.pass, lag.max=96)
 #'
+#' # Run with other parameters and return the RMSE
+#' lm_fit(c(Ta__a1=0.8), model, D, scorefun=rmse, returnanalysis=FALSE)
+#' lm_fit(c(Ta__a1=0.9), model, D, scorefun=rmse, returnanalysis=FALSE)
+#' 
 #' @importFrom stats lm residuals
 #' @export
 lm_fit <- function(prm=NA, model, data, scorefun = NA, returnanalysis = TRUE, printout = TRUE){

R/lm_optim.R

@@ -12,8 +12,10 @@
 #' @title Optimize parameters for onlineforecast model fitted with LM
 #' @param model The onlineforecast model, including inputs, output, kseq, p
 #' @param data The data.list including the variables used in the model.
+#' @param kseq The horizons to fit for (if not set, then model$kseq is used)
 #' @param scorefun The function to be score used for calculating the score to be optimized.
 #' @param cachedir A character specifying the path (and prefix) of the cache file name. If set to \code{""}, then no cache will be loaded or written. See \url{https://onlineforecasting.org/vignettes/nice-tricks.html} for examples.
+#' @param cachererun A logical controlling whether to run the optimization even if the cache exists.
 #' @param printout A logical determining if the score function is printed out in each iteration of the optimization.
 #' @param method The method argument for \code{\link{optim}}.
 #' @param ... Additional parameters to \code{\link{optim}}
@@ -55,7 +57,7 @@
 #'
 #' @importFrom stats optim
 #' @export
-lm_optim <- function(model, data, scorefun = rmse, cachedir="", printout=TRUE, method="L-BFGS-B", ...){
+lm_optim <- function(model, data, kseq = NA, scorefun = rmse, cachedir="", cachererun=FALSE, printout=TRUE, method="L-BFGS-B", ...){
     ## Take the parameters bounds from the parameter bounds set in the model
     init <- model$get_prmbounds("init")
     lower <- model$get_prmbounds("lower")
@@ -64,21 +66,36 @@ lm_optim <- function(model, data, scorefun = rmse, cachedir="", printout=TRUE, m
     if(any(is.na(lower))){ lower[is.na(lower)] <- -Inf}
     if(any(is.na(upper))){ lower[is.na(upper)] <- Inf}
 
+    # Clone the model no matter what (at least model$kseq should not be changed no matter if optimization is stopped)
+    m <- model$clone_deep()
+    if(!is.na(kseq[1])){
+        m$kseq <- kseq
+    }else if(!is.na(m$kseqopt[1])){
+        m$kseq <- m$kseqopt
+    }
+
     ## Caching the results based on some of the function arguments
     if(cachedir != ""){
-        ## Have to insert the parameters in the expressions
-        model$insert_prm(init)
+        # Have to insert the parameters in the expressions to get the right state of the model for unique checksum
+        m$insert_prm(init)
+        # Have to reset the state first to remove dependency of previous calls
+        m$reset_state()
         ## Give all the elements to calculate the unique cache name
-        cnm <- cache_name(lm_fit, lm_optim, model$outputrange, model$regprm, model$transform_data(data),
-                          data[[model$output]], scorefun, init, lower, upper, cachedir = cachedir)
-        ## Maybe load the cached result
-        if(file.exists(cnm)){ return(readRDS(cnm)) }
+        cnm <- cache_name(lm_fit, lm_optim, m$outputrange, m$regprm, m$transform_data(data),
+                          data[[m$output]], scorefun, init, lower, upper, cachedir = cachedir)
+        # Load the cached result if it exists
+        if(file.exists(cnm) & !cachererun){
+            res <- readRDS(cnm)
+            # Set the optimized parameters into the model
+            model$insert_prm(res$par)
+            return(res)
+        }
     }
 
-    ## Run the optimization
+    # Run the optimization
     res <- optim(par = init,
                  fn = lm_fit,
-                 model = model,
+                 model = m,
                  data = data,
                  scorefun = scorefun,
                  printout = printout,
@@ -87,8 +104,9 @@ lm_optim <- function(model, data, scorefun = rmse, cachedir="", printout=TRUE, m
                  upper = upper,
                  method = method,
                  ...)
-    ## Save the result in the cachedir
+    # Save the result in the cachedir
     if(cachedir != ""){ cache_save(res, cnm) }
-    ## Return the result
+    # Set the optimized parameters into the model
+    model$insert_prm(res$par)
     return(res)
 }

R/make_input.R

@@ -28,7 +28,7 @@ make_input <- function(observations, kseq){
     val <- sapply(kseq, function(k){
         observations
     })
-    ## set row and column names
+    # set row and column names
     nams(val) <- paste0('k', kseq)
     return( as.data.frame(val) )
 }

R/make_periodic.R

+## Do this in a separate tmp.R file to check the documentation
+#library(devtools)
+#document()
+#load_all(as.package("../../onlineforecast"))
+#?make_periodic
+
+#' Make an forecast matrix with a periodic time signal.
+#'
+#' This function creates a data.frame with k-steps-ahead values of a periodic time signal,
+#' such that it can be added to a data.list and used inputs to a forecast model.
+#' 
+#' @param time vector of times of class "POSIXct" "POSIXt".
+#' @param kseq vector of integers, representing the desired "k-steps ahead".
+#' @param period a numeric setting the length of the period in seconds.
+#' @param offset a numeric setting an offset in the period start in seconds.
+#' @param tstep step time of k in seconds.
+#' @return Returns a forecast matrix (data.frame) with rownames = times, colnames = k1, k2, k3, ...
+#' The content of the data frame is the hour of day.
+#' @keywords periodic lags data.frame
+#' @seealso make_tday
+#' @examples
+#' # Create a time sequence of 30 min sample period
+#' tseq <- seq(ct("2019-01-01"), ct("2019-02-01 12:00"), by=1800)
+#' 
+#' # Make the three hourly periodic sequence
+#' make_periodic(tseq, 1:10, 3*3600)
+#'
+#' # With an offset of one hour
+#' make_periodic(tseq, 1:10, 3*3600, 3600)
+#' 
+#' @export
+
+make_periodic <- function(time, kseq, period, offset=0, tstep=NA){
+    # If tstep not given, then take it from time assuming same k-step is the sampling period
+    if(is.na(tstep)){
+        tstep <- as.numeric(time[2] - time[1], units="secs")
+    }
+    # The time of day (in the specified units)
+    tday <- sapply(kseq, function(k){
+        tk <- time + k * tstep - offset
+        as.numeric(tk,units="secs") %% period
+    })
+    # set row and column names
+    nams(tday) <- paste0('k', kseq)
+    return( as.data.frame(tday) )
+}