## -----------------------------------------------------------------------------
# Load the package
library(onlineforecast)
# Set the data in D to simplify notation
D <- Dbuilding

## -----------------------------------------------------------------------------
# Print the first time point
D$t[1]
# Set the score period 
D$scoreperiod <- in_range("2010-12-22", D$t)
# Plot to see it
plot(D$t, D$scoreperiod, xlab="Time", ylab="Scoreperiod")

## -----------------------------------------------------------------------------
# Exclude other points example
scoreperiod2 <- D$scoreperiod
scoreperiod2[in_range("2010-12-30",D$t,"2011-01-02")] <- FALSE

## -----------------------------------------------------------------------------
# Generate new object (R6 class)
model <- forecastmodel$new()
# Set the model output
model$output = "heatload"
# Inputs (transformation step)
model$add_inputs(Ta = "Ta",
                 mu = "one()")
# Regression step parameters
model$add_regprm("rls_prm(lambda=0.9)")
# Optimization bounds for parameters
model$add_prmbounds(lambda = c(0.9, 0.99, 0.9999))
# Set the horizons for which the model will be fitted
model$kseq <- 1:36

## -----------------------------------------------------------------------------
# Generate new object
model <- forecastmodel$new()
# Set the model output
model$output = "heatload"

## -----------------------------------------------------------------------------
# Inputs (transformation step)
model$add_inputs(Ta = "Ta",
                 mu = "one()")

## -----------------------------------------------------------------------------
# Regression step parameters
model$add_regprm("rls_prm(lambda=0.9)")

## -----------------------------------------------------------------------------
# The evaluation happens with
eval(parse(text="rls_prm(lambda=0.9)"))
# and the result is stored in
model$regprm 

## -----------------------------------------------------------------------------
# Optimization bounds for parameters
model$add_prmbounds(lambda = c(0.9, 0.99, 0.9999))

## -----------------------------------------------------------------------------
# Set the horizons for which the model will be fitted
model$kseq <- 1:36

## ----output.lines=15----------------------------------------------------------
# Call the optim() wrapper
rls_optim(model, D, kseq = c(3,18))

## -----------------------------------------------------------------------------
# Optimized lambda
model$prm

## -----------------------------------------------------------------------------
# Fit for all on entire period in D
fit1 <- rls_fit(model$prm, model, D)

## -----------------------------------------------------------------------------
# See the summary of the fit
summary(fit1)

## -----------------------------------------------------------------------------
# Put the forecasts in D
D$Yhat1 <- fit1$Yhat
# Plot them for selected horizons
plot_ts(D, c("^heatload$|^Y"), kseq = c(1,6,18,36))

## ----fig.height=4-------------------------------------------------------------
# Select a point
i <- 996-48
# and kseq steps ahead
iseq <- i+model$kseq
# The observations ahead in time
plot(D$t[iseq], D$heatload[iseq], type = "b", xlab = "t", ylab = "y")
title(main=pst("Forecast available at ",D$t[i]))
# The forecasts
lines(D$t[iseq], D$Yhat1[i, ], type = "b", col = 2)
legend("topright", c("Observations",pst("Predictions (",min(model$kseq)," to ",max(model$kseq)," steps ahead)")), lty = 1, col = 1:2)

## ----eval=FALSE---------------------------------------------------------------
#  # This will give error
#  one()

## -----------------------------------------------------------------------------
# Evaluate input expressions
datatr <- model$transform_data(D)
# See what came out
summary.default(datatr)
# In particular for the mu = "one()"
head(datatr$mu)

## ----eval=FALSE---------------------------------------------------------------
#  # Set the function to debug (uncomment the line)
#  #debug(one)
#  # Run the input transformation now and it will stop in one()
#  datatr <- model$transform_data(D)
#  # Set to undebug
#  #undebug(one)

## -----------------------------------------------------------------------------
# Just update the Ta input by
model$add_inputs(Ta = "lp(Ta, a1=0.9)")

## -----------------------------------------------------------------------------
# Define a new model with low-pass filtering of the Ta input
model <- forecastmodel$new()
model$output = "heatload"
model$add_inputs(Ta = "lp(Ta, a1=0.9)",
                 mu = "one()")
model$add_regprm("rls_prm(lambda=0.99)")
model$add_prmbounds(Ta__a1 = c(0.5, 0.9, 0.9999),
                    lambda = c(0.9, 0.99, 0.9999))
model$kseq <- c(3,18)

## -----------------------------------------------------------------------------
model$prmbounds

## -----------------------------------------------------------------------------
# Low-pass filter Ta (with a1=0.9 as defined above)
datatr <- model$transform_data(D)
# Actually, lp() can be called directly (although two warnings are thrown)
Talp <- lp(D$Ta, a1=0.99)

## -----------------------------------------------------------------------------
# Plot the Ta$k1 forecasts
plot(D$t, D$Ta$k1, type="l")
# Add the filtered with a1=0.9
lines(D$t, datatr$Ta[ ,"k1"], col=2)
# Add the filtered with a1=0.99
lines(D$t, Talp[ ,"k1"], col=3)

## ----output.lines=15----------------------------------------------------------
# Optimize the parameters
model$prm <- rls_optim(model, D)$par

## ----fig.height=4-------------------------------------------------------------
# Fit for all horizons
model$kseq <- 1:36
# Fit with RLS
fit2 <- rls_fit(model$prm, model, D)
# Take the forecasts
D$Yhat2 <- fit2$Yhat
# Plot all
plot_ts(D, c("^heatload$|^Y"), kseq = c(1,18))

## -----------------------------------------------------------------------------
summary(fit2)

## -----------------------------------------------------------------------------
# Calculate the score
RMSE1 <- summary(fit1, printit=FALSE)$scoreval
RMSE2 <- summary(fit2, printit=FALSE)$scoreval

## -----------------------------------------------------------------------------
# Check that all NAs in the scoreperiod are at the same positions
all(is.na(fit1$Yhat[fit1$data$scoreperiod, ]) == is.na(fit2$Yhat[fit2$data$scoreperiod, ]))

## -----------------------------------------------------------------------------
# Plot the score for the two models
plot(RMSE1, xlab="Horizon k", ylab="RMSE", type="b", ylim=range(RMSE1,RMSE2))
lines(RMSE2, type="b", col=2)
legend("topleft", c("Input: Ta","Input: Low-pass Ta"), lty=1, col=1:2)

## ----output.lines=28----------------------------------------------------------
make_tday(D$t, kseq=1:3)

## -----------------------------------------------------------------------------
D$tday <- make_tday(D$t, 1:36)

## -----------------------------------------------------------------------------
D$Tao <- make_input(D$Taobs, kseq=1:36)
model$add_inputs(Tao = "lp(Tao, a1=0.99)")

## -----------------------------------------------------------------------------
m1 <- model$clone_deep()