Commit 6a5bf6fd authored by pbac's avatar pbac
Browse files

Copied from old repo

parent ff43ac44
^.*\.Rproj$
^\.Rproj\.user$
^vignettes/onlineforecasting_pdf_source
^vignettes/building-heat-load-forecasting_cache-rls
^vignettes/building-heat-load-forecasting_files
^vignettes/setup-and-use-models_cache
^vignettes/setup-and-use-models_files
^vignettes/setup-and-use-models.html
^vignettes/setup-and-use-models.Rmd
^vignettes/setup-and-use-models.R
^vignettes/setup-data_cache
^vignettes/setup-data_files
^vignettes/setup-data.html
^vignettes/setup-data.Rmd
^vignettes/setup-data.R
^vignettes/cache
^vignettes/tmp-output/
^tests
^vignettes/literature.bib
^misc-R$
^.*\.Rhistory$
\ No newline at end of file
.Rproj.user
.Rhistory
.RData
.Ruserdata
NAMESPACE
*.o
src/onlineforecast\.so
inst/doc
modifications_old_notstaged/
cache/
man/
misc-R/*cache*
vignettes/*cache*
vignettes/*genfig*
vignettes/*_files*
vignettes/tmp-output/
vignettes/setup-data_cache/
vignettes/solar-forecasting_cache-rls/
vignettes/building-heat-load-forecasting_cache/
vignettes/onlineforecasting_pdf_source/onlineforecasting\.tex
vignettes/onlineforecasting_pdf_source/*cache*
vignettes/onlineforecasting_pdf_source/*genfig*
vignettes/onlineforecasting_pdf_source/onlineforecasting-tikzDictionary
vignettes/onlineforecasting_pdf_source/onlineforecasting.log
vignettes/onlineforecasting_pdf_source/onlineforecasting.pdf
Package: onlineforecast
Type: Package
Title: Forecast modelling for online application
Version: 1.0.0
Description: A package 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.
License: GPL-3
Encoding: UTF-8
LazyData: true
Authors@R: c(
person("Peder", "Bacher", email="pbac@dtu.dk", role = c("aut", "cre")),
person("Hjorleifur", "Bergsteinsson", "G", email="hgbe@dtu.dk", role = c("aut","cre"))
Depends: R (>= 3.0)
Imports:
Rcpp (>= 0.12.18),
R6 (>= 2.2.2),
splines (>= 3.1.1),
plotly,
digest
LinkingTo: Rcpp, RcppArmadillo
Suggests:
knitr,
rmarkdown,
R.rsp,
testthat (>= 2.1.0)
VignetteBuilder:
knitr,
R.rsp
RoxygenNote: 7.1.0
URL: http://onlineforecasting.org
BugReports: https://lab.compute.dtu.dk/pbac/onlineforecasting/-/issues
\ No newline at end of file
## Do this in a separate file to see the generated help:
#library(devtools)
#document()
#load_all(as.package("../../onlineforecast"))
#?AR
#' Generate auto-regressive (AR) inputs in a model
#'
#' The AR function can be used in an onlineforecast model formulation. It
#' creates the input matrices for including AR inputs in a model during the
#' transformation stage. It takes the values from the model output in the provided data
#' does the needed lagging.
#'
#' The lags must be given according to the one-step ahead model, e.g.:
#'
#' \code{AR(lags=c(0,1))} will give: Y_{t+1|t} = \eqn{\phi_1} y_{t-0} + \eqn{\phi_2} y_{t-1} + \eqn{\epsilon}_{t+1}
#'
#' and:
#'
#' \code{AR(lags=c(0,3,12))} will give: Y_{t+1|t} = \eqn{\phi}_1 y_{t-0} + \eqn{\phi}_2 y_{t-3} + \eqn{\phi}_3 y_{t-12} + \eqn{\epsilon}_{t+1}
#'
#' Note, that
#'
#' For k>1 the coefficients will be fitted individually for each horizon, e.g.:
#'
#' \code{AR(lags=c(0,1))} will be the multi-step AR: Y_{t+k|t} = \eqn{\phi}_{1,k} y_{t-0} + \eqn{\phi}_{2,k} y_{t-1} + \eqn{\epsilon}_{t+k|t}
#'
#' See the details in ??(ref til vignette).
#'
#' @title Auto-Regressive (AR) input
#' @param lags integer vector: The lags of the AR to include.
#' @return A list of matrices, one for each lag in lags, each with columns according to model$kseq.
#' @examples
#'
#' # Setup data and a model for the example
#'
#' model <- forecastmodel$new()
#' model$output = "heatload"
#' # Use the AR in the transformation stage
#' model$add_inputs(AR = "AR(c(0,1))")
#' # Regression parameters
#' model$add_regprm("rls_prm(lambda=0.9)")
#' # kseq must be added
#' model$kseq <- 1:4
#' # In the transformation stage the AR input will be generated
#' # See that it generates two input matrices, simply with the lagged heat load at t for every k
#' model$transform_data(subset(D, 1:10))
#' # Fit with recursive least squares (no parameters prm in the model)
#' fit <- rls_fit(c(lambda=0.99), model, D, returnanalysis=TRUE)
#' # Plot the result, see "?plot_ts.rls_fit"
#' plot(fit, xlim=c(asct("2010-12-20"),max(D$t)))
#' # Plot for a short period with peaks
#' plot(fit, xlim=c("2011-01-05","2011-01-07"))
#' # For online updating, see ??ref{vignette}:
#' # the needed lagged output values are stored in the model for next time new data is available
#' model$yAR
#' # The maximum lag needed is also kept
#' model$maxlagAR
#'
#' @export
AR <- function(lags){
# Just sort them first
lags <- sort(lags)
# Get the data and the model from an environment above (this way has worked until now, not exactly sure why the environments are in this order)
data <- parent.env(parent.frame())$data
model <- parent.env(parent.frame())$self$model
# Remember the max lag for later, only if bigger than current (should make set function doing this check)
if(is.na(model$maxlagAR) | max(lags) > model$maxlagAR){
model$maxlagAR <- max(lags)
}
# Setup the AR inputs, one matrix for each lag
retval <- lapply(lags, function(lag){
# Check if saved output values for AR exists
if(is.na(model$yAR[1])){
# First time its called, so just use output values from data
val <- matrix(lag(data[[model$output]], lag), nrow=length(data$t), ncol=length(model$kseq))
}else{
y <- c(model$yAR, data$y)
# Find the seq for the new y lagged vector
iend <- (length(y)-lag)
istart <- iend - length(data$y) + 1
# Take the sequence
y <- y[istart:iend]
# Insert in a matrix with column for each k
val <- matrix(y, nrow=length(data$t), ncol=length(model$kseq))
}
# Name the columns and return
nams(val) <- pst("k", model$kseq)
return(val)
})
names(retval) <- pst("lag", lags)
return(retval)
}
# Generated by using Rcpp::compileAttributes() -> do not edit by hand
# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
#' Low pass filtering of a vector.
#'
#' This function returns a vector which is x through a unity gain first-order low-pass filter.
#'
#' @name lp_vector_cpp
#' @param x A numeric vector
#' @param a1 the first order low-pass filter coefficient
NULL
lp_vector_cpp <- function(x, a1) {
.Call('_onlineforecast_lp_vector_cpp', PACKAGE = 'onlineforecast', x, a1)
}
#' Calculating k-step recursive least squares estimates
#'
#' This function applies the k-step recursive least squares scheme to estimate
#' parameters in a linear regression model.
#'
#' @name rls_update_cpp
#' @param y Vector of observation
#' @param X Matrix of input variables (design matrix)
#' @param theta Vector of parameters (initial value)
#' @param P Covariance matrix (initial value)
#' @param lambda Forgetting factor
#' @param k Forecast horizon
#' @param n Length of the input
#' @param np Dimension of P (np x np)
#' @param istart Start index
#' @param kmax Keep only the last kmax rows for next time
NULL
rls_update_cpp <- function(y, X, theta, P, lambda, k, n, np, istart, kmax) {
.Call('_onlineforecast_rls_update_cpp', PACKAGE = 'onlineforecast', y, X, theta, P, lambda, k, n, np, istart, kmax)
}
# Do this in a separate file to see the generated help:
# library(devtools)
# document()
# load_all(as.package("../../onlineforecast"))
# ?as.data.list
# ?as.data.list.data.frame
#' These functions will convert the object into a data.list.
#'
#' A data.list is simply a list of vectors and data.frames. For the use in the
#' onlineforecast package the following format must be kept:
#'
#' - t: A vector of time.
#'
#' - vectors with same length as t: Holds observations and values synced to time t.
#'
#' - data.frames with number of rows as time t: Holds forecasts in each column named by \code{kxx} where \code{xx} is the
#' horizon, e.g. \code{k0} is synced as observations, and \code{k1} is one-step ahead.
#'
#' @title Convert to data.list class
#' @param object The object to be converted into a data.list
#' @param ...
#' @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
#'
#' @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
#' \code{xx} is the horizon. See the examples.
#'
#' @title Convertion of data.frame into data.list
#' @param object
#' @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)
#' as.data.list(X)
#'
#' # Convert a dataframe with time, forecast and an observed variable
#' X <- data.frame(t=1:10, x.k1=1:10, x.k2=10:1, y.obs=1:10, y.k1=1:10, y.k2=1:10)
#' as.data.list(X)
#'
#' # Can be converted back and forth
#' X
#' as.data.frame(as.data.list(X))
#'
#' @export
as.data.list.data.frame <- function(object) {
X <- object
#TEST
#grep("\\.[hk][[:digit:]]+$", c("Ta.k1","Ta.k2","I.h1"))
# Check which columns hold forecasts and must be returned as data.frames in the data.list
inmsfor <- grep("\\.[hk][[:digit:]]+$", names(X))
#
if(length(inmsfor) > 0){
# Find the names of them
nmsfor <- unique(unlist(
getse(strsplit(names(X)[inmsfor], "\\."), 1)
))
# Group all in a list
# Note that "u.k2y" is matched, but it maybe shouldn't be: grep("\\.k[:digit:]*", c("ij.k1","i","u.k2y"))
L <- lapply(nmsfor, function(nm){
return(inmsfor[grep(pst("^",nm), names(X)[inmsfor])])
})
names(L) <- nmsfor
# The vectors (time t, and observations)
Lobs <- as.list((1:ncol(X))[-inmsfor])
names(Lobs) <- names(X)[-inmsfor]
}else{
# No forecasts found
L <- list()
# The vectors (time t, and observations)
Lobs <- as.list((1:ncol(X)))
names(Lobs) <- names(X)
}
#
# Combine and sort like the order they came in
L <- c(L, Lobs)
L <- L[order(unlist(getse(L, 1)))]
#
# Make the data.list
val <- lapply(L, function(i) {
tmp <- X[ ,i]
if(!is.null(dim(tmp))){
# Its a forecast, hence a data.frame, remove "name." from the column names
names(tmp) <- getse(strsplit(names(tmp), "\\."), 2)
}
return(tmp)
})
class(val) <- "data.list"
return(val)
}
# Do this in a separate file to see the generated help:
#library(devtools)
#document()
#load_all(as.package("../../onlineforecast"))
#?asct
#?asct.default
#' The object is converted into POSIXct with tz="GMT".
#'
#' A simple helper, which wraps \code{\link{as.POSIXct}}` and sets the time zone to "GMT" per default.
#'
#' @title Convertion to POSIXct
#' @param object The object to convert can be: character, numeric, POSIXct or POSIXlt
#' @param tz Timezone. If set, then the time zone will be changed of the object.
#' @param ...
#' @return An object of class POSIXct
#' @section Methods:
#' @examples
#'
#'
#' # Create a POSIXct with tz="GMT"
#' asct("2019-01-01")
#' class(asct("2019-01-01"))
#' asct("2019-01-01 01:00:05")
#' # Convert to POSIXct
#' class(asct(as.POSIXlt(x)))
#' # To seconds and back again
#' asct(as.numeric(x, units="sec"))
#' # --------
#' # Convert character of time which has summer time leaps
#' # Example from CET (with CEST which is winter time)
#' #
#' # The point of shifting to and from summer time:
#' # DST Start (Clock Forward) DST End (Clock Backward)
#' # Sunday, March 31, 02:00 Sunday, October 27, 03:00
#' # --------
#' # From to winter time to summer time
#' txt <- c("2019-03-31 01:00",
#' "2019-03-31 01:30",
#' "2019-03-31 03:00",
#' "2019-03-31 03:30")
#' x <- asct(txt, tz="CET")
#' x
#' asct(x, tz="GMT")
#' # BE AWARE of this conversion of the 02:00: to 02:59:59 (exact time of shift) will lead to a wrong conversion
#' txt <- c("2019-03-31 01:30",
#' "2019-03-31 02:00",
#' "2019-03-31 03:30")
#' x <- asct(txt, tz="CET")
#' x
#' asct(x, tz="GMT")
#' # Which a diff on the time can detect, since all steps are not equal
#' plot(diff(asct(x, tz="GMT")))
#' # --------
#' # Shift to winter time is more problematic
#' # It works like this
#' txt <- c("2019-10-27 01:30",
#' "2019-10-27 02:00",
#' "2019-10-27 02:30",
#' "2019-10-27 03:00",
#' "2019-10-27 03:30")
#' x <- asct(txt, tz="CET")
#' x
#' asct(x, tz="GMT")
#' # however, timestamps can be given like this
#' txt <- c("2019-10-27 01:30",
#' "2019-10-27 02:00",
#' "2019-10-27 02:30",
#' "2019-10-27 02:00",
#' "2019-10-27 02:30",
#' "2019-10-27 03:00",
#' "2019-10-27 03:30")
#' x <- asct(txt, tz="CET")
#' x
#' asct(x, tz="GMT")
#' # Again can be detected, since all steps are not equal
#' plot(diff(asct(x, tz="GMT")))
#' # This can be fixed by (note that it can go wrong, e.g. with gaps around convertion etc.)
#' asct(x, tz="GMT", duplicatedadd=3600)
#'
#' @export
asct <- function(object, tz, ...){
UseMethod("asct")
}
#' @rdname asct
#' @section Methods:
#' - asct.character: Simply a wrapper for \code{as.POSIXct} with default \code{tz}
#' @export
asct.character <- function(object, tz = "GMT", ...){
as.POSIXct(object, tz=tz, ...)
}
#' @rdname asct
#' @section Methods:
#' - asct.POSIXct: Changes the time zone of the object if \code{tz} is given.
#' @export
asct.POSIXct <- function(object, tz = NA, duplicatedadd = NA){
if(!is.na(tz)){
attr(object, "tzone") <- tz
}
if(!is.na(duplicatedadd)){
# To mitigate the problem of duplicated timestamps at the shift to winter time
# then shift the time of duplicated time stamps with the seconds given
object[which(duplicated(object))] <- object[which(duplicated(object))] + duplicatedadd
}
return(object)
}
#' @rdname asct
#' @section Methods:
#' - asct.POSIXlt: Converts to POSIXct.
#' @export
asct.POSIXlt <- function(object, tz = NA, duplicatedadd = NA){
as.POSIXct(asct.POSIXct(object, tz, duplicatedadd))
}
#' @rdname asct
#' @section Methods:
#' - asct.numeric: Converts from UNIX time in seconds to POSIXct with \code{tz} as GMT.
#' @export
asct.numeric <- function(object){
ISOdate(1970, 1, 1, 0) + object
}
# Do this in a separate file to see the generated help:
#library(devtools)
#document()
#load_all(as.package("../../onlineforecast"))
#?aslt
#?aslt.default
#' The argument is converted into POSIXlt with tz="GMT".
#'
#'
#'
#' @title Convertion to POSIXlt
#' @param object
#' @param tz Timezone. If set, then the time zone will be changed of the object.
#' @param ...
#' @return An object of class POSIXlt
#' @section Methods:
#' #' @examples
#'
#' # Create a POSIXlt with tz="GMT"
#' aslt("2019-01-01")
#' class(aslt("2019-01-01"))
#' aslt("2019-01-01 01:00:05")
#'
#' # Convert between time zones
#' x <- aslt("2019-01-01", tz="CET")
#' aslt(x,tz="GMT")
#'
#' # To seconds and back again
#' aslt(as.numeric(x, units="sec"))
#'
#' @export
aslt <- function(object, tz, ...){
UseMethod("aslt")
}
#' @rdname aslt
#' @section Methods:
#' - aslt.character: Simply a wrapper for \code{as.POSIXlt}
#' @export
aslt.character <- function(object, tz = "GMT", ...){
as.POSIXlt(object, tz = tz, ...)
}
#' @rdname aslt
#' @section Methods:
#' - aslt.POSIXct: Converts to POSIXct.
#' @export
aslt.POSIXct <- function(object, tz = NA){
if(!is.na(tz)){
attr(object, "tzone") <- tz
}
as.POSIXlt(object)
}
#' @rdname aslt
#' @section Methods:
#' - aslt.POSIXlt: Changes the time zone of the object if tz is given.
#' @export
aslt.POSIXlt <- function(object, tz = NA){
if(!is.na(tz)){
attr(object, "tzone") <- tz
}
return(object)
}
#' @rdname aslt
#' @section Methods:
#' - aslt.numeric: Converts from UNIX time in seconds to POSIXlt.
#' @export
aslt.numeric <- function(object){
as.POSIXlt(ISOdate(1970, 1, 1, 0) + object)
}
# # Do this in a separate file to see the generated help:
# library(devtools)
# document()
# load_all(as.package("../../onlineforecast"))
# ?bspline
#' Compute base splines of a variable using the R function \code{splines::bs}, use in the transform stage.
#'
#' Simply wraps the \code{splines::bs}, such that it can be used in the transformation stage.
#'
#' See the help for all arguments with \code{?splines::bs}. NOTE that two arguments have different default values.
#'
#' For more examples of use see ??ref(solar forecast vignette).
#'
#' @family Transform stage functions
#'
#' @param X data.frame (as part of data.list) with horizons as columns named \code{kxx} (i.e. one for each horizon)
#' @param Boundary.knots The value is NA: then the boundaries are set to the range of each horizons (columns in X). See \code{?splines::bs}
#' @param intercept Default value is TRUE: in an onlineforecast model there is no intercept per defauls (set by \code{ones()}. See \code{?splines::bs}
#' @param df See \code{?splines::bs}
#' @param knots See \code{?splines::bs}
#' @param degree See \code{?splines::bs}
#' @return List of data frames with the computed base splines, each with columns for the same horizons as in X
#' @examples
#'
#' # How to make a diurnal curve using splines
#' # Select first 54 hours from the load data
#' D <- subset(Dbuildingheatload, 1:54, kseq=1:4)
#' # Make the hour of the day as a forecast input
#' D$tday <- make_tday(D$t, kseq=1:4)
#' D$tday
#'
#' # Calculate the base splines for each column in tday
#' L <- bspline(D$tday)
#'
#' # Now L holds a data.frame for each base spline
#' str(L)
#' # Hence this will result in four inputs for the regression model
#'
#' # Plot (note that the splines period starts at tday=0)
#' plot(D$t, L$bs1$k1, type="s")
#' for(i in 2:length(L)){
#' lines(D$t, L[[i]]$k1, col=i, type="s")
#' }
#'
#' # In a model formulation it will be:
#' model <- forecastmodel$new()
#' model$add_inputs(mutday = "bspline(tday)")
#' # Such that at the transform stage will give the same as above
#' model$transform_data(D)
#'
#'
#' @export
bspline <- function(X, Boundary.knots = NA, intercept = TRUE, df = NULL, knots = NULL, degree = 3) {
# If a list, then call on each element
if (class(X) == "list") {
# Call again for each element
val <- lapply(1:length(X), function(i) {
bspline(X[[i]], df = df, knots = knits, degree = degree, intercept = intercept,
Boundary.knots = Boundary.knots)
})
nams(val) <- nams(X)
return(val)
}
# X is a data.frame or matrix
# First find the horizons, they are used in the end
nms <- nams(X)
# Run for each horizon and calculate the basis splines
L <- lapply(1:ncol(X), function(i) {
if (is.na(Boundary.knots[1])) {
Boundary.knots <- range(X[, i], na.rm=TRUE)
}
spls <- splines::bs(X[, i], Boundary.knots = Boundary.knots, degree = degree, df = df,
knots = knots, intercept = intercept)
return(spls)
})
# Now we have a bs value in L for each horizon
# Separate each basespline in a data.frame with all horizons
L <- lapply(1:ncol(L[[1]]), function(i) {
tmp <- lapply(L, function(x) {
x[ ,i]
})
tmp <- data.frame(do.call("cbind", tmp))
nams(tmp) <- nms
return(tmp)
})
# Set the extra name
nams(L) <- pst("bs", 1:length(L))
return(L)
}