diff --git a/.Rbuildignore b/.Rbuildignore
index 1e05c5746714b35fff3b992d50f0e19c30fe7a3f..5a454b84f52f82e9ec73cf7d61732c3ebbb38187 100644
--- a/.Rbuildignore
+++ b/.Rbuildignore
@@ -1,21 +1,12 @@
 ^.*\.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
+^make.R
+^data/all
+^vignettes/tmp-output
+^vignettes/make.R
+^vignettes/init.R
 ^vignettes/cache
-^vignettes/tmp-output/
-^tests
 ^vignettes/literature.bib
+^tests
 ^misc-R$
 ^.*\.Rhistory$
\ No newline at end of file
diff --git a/DESCRIPTION b/DESCRIPTION
index cb28706a6a679af8f34b6a46d2d7c459ae672b24..0611fd22d4e4a6e84c6a0e3d7e30c2254cf30ce5 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,30 +1,30 @@
 Package: onlineforecast
 Type: Package
-Title: Forecast modelling for online application
+Title: Forecast Modelling for Online Applications
 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.
+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.
 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"))
+    person("Hjorleifur G", "Bergsteinsson", email="hgbe@dtu.dk", role = c("aut")))
 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,
+    plotly,
     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
+BugReports: https://lab.compute.dtu.dk/packages/onlineforecast/-/issues
\ No newline at end of file
diff --git a/R/AR.R b/R/AR.R
index 080a158e2c80ab7cd9f0c7313d4aeab521aeb984..d34dfa61e6a569689c62756dfbc672ec517c1839 100644
--- a/R/AR.R
+++ b/R/AR.R
@@ -33,7 +33,7 @@
 #' @examples
 #'
 #' # Setup data and a model for the example
-#'
+#' D <- Dbuilding
 #' model <- forecastmodel$new()
 #' model$output = "heatload"
 #' # Use the AR in the transformation stage
@@ -45,15 +45,15 @@
 #' # 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_ts(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"))
-
+#' plot_ts(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
diff --git a/R/bspline.R b/R/bspline.R
index 92848d62976602fe026da41c5caefdf82fbff67e..72ae30792b41b58a1ce1b94c8e978dee605532fe 100644
--- a/R/bspline.R
+++ b/R/bspline.R
@@ -26,7 +26,7 @@
 #'
 #' # How to make a diurnal curve using splines
 #' # Select first 54 hours from the load data
-#' D <- subset(Dbuildingheatload, 1:54, kseq=1:4)
+#' D <- subset(Dbuilding, 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
diff --git a/R/data.list.R b/R/data.list.R
index e4c0379c4cb414a1e7ffb7b8b094dcc3062224de..dd83eab7985130e4dfe05aebba3aa13489e46b64 100644
--- a/R/data.list.R
+++ b/R/data.list.R
@@ -55,7 +55,7 @@ data.list <- function(...) {
 #' @return a data.list with the subset.
 #' @examples
 #' # Use the data.list with building heat load 
-#' D <- Dbuildingheatload
+#' D <- Dbuilding
 #' # Take a subset for the example
 #' D <- subset(D, 1:10, nms=c("t","Ta.obs","Ta","I.obs","I"), kseq=1:3)
 #' 
@@ -82,7 +82,7 @@ data.list <- function(...) {
 #' subset(D, kseq=2, pattern="^I")
 #' 
 #' # Take data for Ta and lag the forecasts (good for plotting and fitting a model)
-#' X <- subset(Dbuildingheatload, 1:1000, pattern="^Ta", kseq = 10, lagforecasts = TRUE)
+#' X <- subset(Dbuilding, 1:1000, pattern="^Ta", kseq = 10, lagforecasts = TRUE)
 #' 
 #' # A scatter plot between the forecast and the observations (try lagforecasts = FALSE and see the difference)
 #' plot(X$Ta$k10, X$Ta.obs)
@@ -193,7 +193,7 @@ subset.data.list <- function(x, subset = NA, nms = NA, kseq = NA, lagforecasts =
 #' @examples
 #'
 #' #' # Use the data.list with building heat load 
-#' D <- Dbuildingheatload
+#' D <- Dbuilding
 #' # Take a subset
 #' D <- subset(D, 1:5, nms=c("t","Ta.obs","Ta","I.obs","I"), kseq=1:3)
 #'
@@ -231,7 +231,7 @@ as.data.frame.data.list <- function(x){
 #' @param ... Passed to pairs().
 #' @examples
 #' # Take a subset for the example
-#' D <- subset(Dbuildingheatload, c("2010-12-15","2011-01-15"), pattern="^Ta|^I", kseq=1:3)
+#' D <- subset(Dbuilding, c("2010-12-15","2011-01-15"), pattern="^Ta|^I", kseq=1:3)
 #' pairs(D)
 #'
 #' # If the forecasts and the observations are not aligned in time it is easy to see by comparing to the previous plot.
@@ -240,7 +240,7 @@ as.data.frame.data.list <- function(x){
 #' # Hence if the forecasts were not synced properly, then it can be detected using this type of plot.
 #'
 #' # Alternatively, lag when taking the subset
-#' D <- subset(Dbuildingheatload, c("2010-12-15","2011-01-15"), pattern="^Ta|^I", kseq=1:3, lagforecasts=TRUE)
+#' D <- subset(Dbuilding, c("2010-12-15","2011-01-15"), pattern="^Ta|^I", kseq=1:3, lagforecasts=TRUE)
 #' pairs(D, lagforecasts=FALSE)
 #' 
 #' @export
@@ -262,7 +262,7 @@ pairs.data.list <- function(x, subset = NA, nms = NA, kseq = NA, lagforecasts =
 #' @return The tables generated.
 #'
 #' # Check a data.list (see \code{?\link{check.data.list}})
-#' check(Dbuildingheatload)
+#' check(Dbuilding)
 #' 
 #' @export
 check <- function(object){
@@ -293,15 +293,15 @@ check <- function(object){
 #' @return The tables generated.
 #'
 #' # Check a data.list (see \code{?\link{check.data.list}})
-#' check(Dbuildingheatload)
+#' check(Dbuilding)
 #'
 #' # Vector with observations not same length as t
-#' D <- Dbuildingheatload
+#' D <- Dbuilding
 #' D$heatload <- D$heatload[1:10]
 #' check(D)
 #'
 #' # Some NAs in k1 forecast
-#' D <- Dbuildingheatload
+#' D <- Dbuilding
 #' D$Ta$k1[1:1500] <- NA
 #' check(D)
 #'
@@ -400,16 +400,16 @@ check.data.list <- function(object){
 #'
 #' @examples
 #'
-#' Dbuildingheatload == Dbuildingheatload
+#' Dbuilding == Dbuilding
 #'
-#' D <- Dbuildingheatload
+#' D <- Dbuilding
 #' D$Ta$k2[1] <- NA
-#' Dbuildingheatload == D
+#' Dbuilding == D
 #'
-#' D <- Dbuildingheatload
+#' D <- Dbuilding
 #' names(D)[5] <- "I"
 #' names(D)[6] <- "Ta"
-#' Dbuildingheatload == D
+#' Dbuilding == D
 #' 
 #' 
 
diff --git a/R/forecastmodel.R-documentation.R b/R/forecastmodel.R-documentation.R
index abb3c1c1577d74c34434ca6afb27ac611208ea97..b351e9fcec76e2e6c5247dd7848c30dd2a1d4b7a 100644
--- a/R/forecastmodel.R-documentation.R
+++ b/R/forecastmodel.R-documentation.R
@@ -168,13 +168,13 @@
 #' @examples
 #' 
 #' # Check if the model is setup and can be used with a given data.list
-#' model$check(Dbuildingheatload)
+#' model$check(Dbuilding)
 #' # Add the model output
 #' model$output <- "heatload"
-#' model$check(Dbuildingheatload)
+#' model$check(Dbuilding)
 #' # Add the horizons to fit for
 #' model$kseq <- 1:4
 #' # No errors, it's fine :)
-#' model$check(Dbuildingheatload)
+#' model$check(Dbuilding)
 NULL
 
diff --git a/R/in_range.R b/R/in_range.R
index 0f168cf30f8515201584eccc0841ee6bf682fadf..efe14f5c56fd7019388e428ecf3011ebe9d224a3 100644
--- a/R/in_range.R
+++ b/R/in_range.R
@@ -21,7 +21,7 @@
 #' @examples
 #'
 #' # Take a subset
-#' D <- subset(Dbuildingheatload, c("2010-12-15", "2011-01-01"))
+#' D <- subset(Dbuilding, c("2010-12-15", "2011-01-01"))
 #'
 #' # Just a logical returning TRUE in a specified period
 #' in_range("2010-12-20", D$t, "2010-12-22")
diff --git a/R/inputs_class.R-documentation.R b/R/inputs_class.R-documentation.R
index d5770a9fba6b775714a1c77ad58e34874622f59b..06aa65ce75585983b7a4f972db02f6c39ff5451f 100644
--- a/R/inputs_class.R-documentation.R
+++ b/R/inputs_class.R-documentation.R
@@ -52,7 +52,7 @@
 #'
 #' # Now the transformation stage can be carried out to create the regression stage data
 #' # Take a data.list subset for the example
-#' D <- subset(Dbuildingheatload, 1:10, kseq=1:4)
+#' D <- subset(Dbuilding, 1:10, kseq=1:4)
 #' # Transform the data
 #' model$inputs[[1]]$evaluate(D)
 #' # What happens is simply that the expression is evaluated with the data
@@ -73,7 +73,7 @@
 #' # the lp() has saved it's state for next time
 #' model$inputs[[1]]$state_L
 #' # New data arrives
-#' Dnew <- subset(Dbuildingheatload, 11, kseq=1:4)
+#' Dnew <- subset(Dbuilding, 11, kseq=1:4)
 #' # So in lp() the state is read and it continues
 #' model$inputs[[1]]$evaluate(Dnew)
 #'
diff --git a/R/lm_fit.R b/R/lm_fit.R
index 17a54bf72b73dfef48d4f132076640675ae7ae47..aadd593396247ad8dfb6f7e8082a31d03d6755dc 100644
--- a/R/lm_fit.R
+++ b/R/lm_fit.R
@@ -32,7 +32,7 @@
 #' @examples
 #'
 #' # Take data (See vignette ??(ref) for better model and more details)
-#' D <- subset(Dbuildingheatload, c("2010-12-15", "2011-01-01"))
+#' D <- subset(Dbuilding, c("2010-12-15", "2011-01-01"))
 #' D$y <- D$heatload
 #' # Define a model 
 #' model <- forecastmodel$new()
diff --git a/R/lm_optim.R b/R/lm_optim.R
index e27c8182765216d8eff9c72cdf029905d1f5ceb9..fd9f9b6f5cb80f5ce616ab401402a850bc6a0594 100644
--- a/R/lm_optim.R
+++ b/R/lm_optim.R
@@ -23,7 +23,7 @@
 #' @examples
 #'
 #' # Take data (See vignette ??(ref) for better model and more details)
-#' D <- subset(Dbuildingheatload, c("2010-12-15", "2011-01-01"))
+#' D <- subset(Dbuilding, c("2010-12-15", "2011-01-01"))
 #' D$y <- D$heatload
 #' # Define a model 
 #' model <- forecastmodel$new()
diff --git a/R/lm_predict.R b/R/lm_predict.R
index bff7ba22ede9896019a139536de6b674ce57decc..5f573b4b7bc487c11ea860c5ce8fb587cf49d3f3 100644
--- a/R/lm_predict.R
+++ b/R/lm_predict.R
@@ -10,7 +10,7 @@
 #'
 #'
 #' # Take data
-#' D <- subset(Dbuildingheatload, c("2010-12-15", "2011-01-01"))
+#' D <- subset(Dbuilding, c("2010-12-15", "2011-01-01"))
 #' D$y <- D$heatload
 #' # Define a model 
 #' model <- forecastmodel$new()
diff --git a/R/make_input.R b/R/make_input.R
index b19d8b96e62304b05b78b38ffcb3ee7cb40af0ba..ff525d114e1b2d1c60182314a3340b39c9cb1c6c 100644
--- a/R/make_input.R
+++ b/R/make_input.R
@@ -15,7 +15,7 @@
 #' @examples
 #'
 #' # Data for example
-#' D <- subset(Dbuildingheatload, c("2010-12-15","2010-12-20"))
+#' D <- subset(Dbuilding, c("2010-12-15","2010-12-20"))
 #' 
 #' # Generate the input
 #' make_input(D$heatload, 1:4)
diff --git a/R/ones.R b/R/ones.R
index 2f84ba33771d5a46867d255496e08e257a04ff1c..e2d3237e9a2a82c344c9c4193ecf0968803fd533 100644
--- a/R/ones.R
+++ b/R/ones.R
@@ -22,7 +22,7 @@
 #' # Set the forecast horizons
 #' model$kseq <- 1:4
 #' # During the transformation stage the ones will be generated for the horizons
-#' model$transform_data(subset(Dbuildingheatload, 1:7))
+#' model$transform_data(subset(Dbuilding, 1:7))
 #'
 #' @export
 
diff --git a/R/par_ts.R b/R/par_ts.R
index 187db70b4abe4322a717e6556326071d7a448222..ac8f51b116fa945dbf2c7c9fd40cd30426e57622 100644
--- a/R/par_ts.R
+++ b/R/par_ts.R
@@ -31,7 +31,7 @@
 #' @examples
 #'
 #' # Data for plots
-#' D <- subset(Dbuildingheatload, 1:192)
+#' D <- subset(Dbuilding, 1:192)
 #'
 #' # See the parameters which can be set
 #' p <- par_ts()
diff --git a/R/plot_ts.R b/R/plot_ts.R
index 86899aae71172cf42b83178f16e3ac0c04a45b88..f4c134c3adf60bc162c7e291f383e3b59f8ced7f 100644
--- a/R/plot_ts.R
+++ b/R/plot_ts.R
@@ -39,7 +39,7 @@
 #' @examples
 #'
 #' # Time series plots for \code{data.list}, same as for \code{data.frame} except use of \code{kseq}
-#' D <- Dbuildingheatload
+#' D <- Dbuilding
 #' plot_ts(D, c("heatload","Ta"), kseq=c(1,24))
 #' # Make two plots (and set the space for the legend)
 #' plot_ts(D, c("heatload","Ta"), kseq=c(1,24), legendspace=11)
@@ -453,7 +453,7 @@ plot_ts_series <- function(data, pattern, iplot = 1,
 #' @examples
 #'
 #' # Fit a model (see vignette 'setup-and-use-model'
-#' D <- Dbuildingheatload
+#' D <- Dbuilding
 #' D$scoreperiod <- in_range("2010-12-22", D$t)
 #' model <- forecastmodel$new()
 #' model$output = "heatload"
diff --git a/R/plotly_ts.R b/R/plotly_ts.R
index 70629801393122a8df60d43d78f52ccc8f410349..01166f7ee600d11627aa6b7bc782a20e12f5cb93 100644
--- a/R/plotly_ts.R
+++ b/R/plotly_ts.R
@@ -17,7 +17,7 @@
 #' \code{\link{plot_ts}}
 #' @examples
 #'
-#' D <- Dbuildingheatload
+#' D <- Dbuilding
 #' plotly_ts(D, c("heatload","Ta"), kseq=c(1,24))
 #' plotly_ts(D, c("heatload","Ta"), kseq=c(1,24))
 #' plotly_ts(D, c("heatload","Ta$|Ta.obs$"), kseq=c(1,24))
diff --git a/R/rls_fit.R b/R/rls_fit.R
index 7dfb91ebf53f203f967651ca1dc46a7f118dd262..cc5d391709f72d16f8dc683efd60c1bdd00a4611 100644
--- a/R/rls_fit.R
+++ b/R/rls_fit.R
@@ -48,7 +48,7 @@
 #'
 #'
 #' # Take data (See vignette ??(ref) for better model and more details)
-#' D <- subset(Dbuildingheatload, c("2010-12-15", "2011-01-01"))
+#' D <- subset(Dbuilding, c("2010-12-15", "2011-01-01"))
 #' D$y <- D$heatload
 #' # Define a model 
 #' model <- forecastmodel$new()
diff --git a/R/rls_optim.R b/R/rls_optim.R
index 23995fb8d0cc366f828b35935f793773ce1a6167..fb606fa31572d7a8b7c17f3602e91b04ff1fba44 100644
--- a/R/rls_optim.R
+++ b/R/rls_optim.R
@@ -23,7 +23,7 @@
 #' @examples
 #'
 #' # Take data (See vignette ??(ref) for better model and more details)
-#' D <- subset(Dbuildingheatload, c("2010-12-15", "2011-01-01"))
+#' D <- subset(Dbuilding, c("2010-12-15", "2011-01-01"))
 #' D$y <- D$heatload
 #' # Define a model 
 #' model <- forecastmodel$new()
diff --git a/R/rls_predict.R b/R/rls_predict.R
index 8c6213aa24c1ea92e068a93380c3d72a8a51ac08..47f8f608a7d6e8a5cf9598ac21477d2458a7aff4 100644
--- a/R/rls_predict.R
+++ b/R/rls_predict.R
@@ -9,7 +9,7 @@
 #' @examples
 #'
 #' # Take data (See vignette ??(ref) for better model and more details)
-#' D <- subset(Dbuildingheatload, c("2010-12-15", "2011-01-01"))
+#' D <- subset(Dbuilding, c("2010-12-15", "2011-01-01"))
 #' D$y <- D$heatload
 #' # Define a model 
 #' model <- forecastmodel$new()
@@ -44,7 +44,7 @@
 #' plot_ts(D, c("y|Yhat"), kseq=1)
 #'
 #' # Recursive updating and prediction
-#' Dnew <- subset(Dbuildingheatload, c("2011-01-01", "2011-01-02"))
+#' Dnew <- subset(Dbuilding, c("2011-01-01", "2011-01-02"))
 #'
 #' for(i in 1:length(Dnew$t)){
 #'     # New data arrives
diff --git a/R/rls_prm.R b/R/rls_prm.R
index 22a95ca2cc3c50c092ab4fbb505f30f79bd15ad9..ddec0890fdb6060282ae470f77f4773ff3076501 100644
--- a/R/rls_prm.R
+++ b/R/rls_prm.R
@@ -13,7 +13,7 @@
 #' @examples
 #'
 #' # Take data (See vignette ??(ref) for better model and more details)
-#' D <- subset(Dbuildingheatload, c("2010-12-15", "2011-01-01"))
+#' D <- subset(Dbuilding, c("2010-12-15", "2011-01-01"))
 #' D$y <- D$heatload
 #' D$scoreperiod <- in_range("2010-12-20", D$t)
 #' # Define a model 
diff --git a/R/rls_summary.R b/R/rls_summary.R
index 5bdb99ff21d73f01f014d1fe9200db8ea66ef447..fd10077fc9727c20b77f3fc7d0bd1ce5d93c8d48 100644
--- a/R/rls_summary.R
+++ b/R/rls_summary.R
@@ -40,7 +40,7 @@
 #' @examples
 #'
 #' # Take data (See vignette ??(ref) for better model and more details)
-#' D <- subset(Dbuildingheatload, c("2010-12-15", "2011-01-01"))
+#' D <- subset(Dbuilding, c("2010-12-15", "2011-01-01"))
 #' D$y <- D$heatload
 #' D$scoreperiod <- in_range("2010-12-20", D$t)
 #' # Define a model
diff --git a/R/stairs.R b/R/stairs.R
index 47f10d03d5c72565316955c7c0560626ceaddb3d..8df3a7fd5f666008cbf3f0b91a33161136e76cb8 100644
--- a/R/stairs.R
+++ b/R/stairs.R
@@ -27,16 +27,16 @@
 #' stairs(1:10, x)
 #'
 #' # Use for time series plotting
-#' plot_ts(Dbuildingheatload, "heatload", c("2010-12-15","2010-12-16"), plotfun=stairs)
+#' plot_ts(Dbuilding, "heatload", c("2010-12-15","2010-12-16"), plotfun=stairs)
 #'
 #' # Set it globally for all plot_ts
 #' p <- par_ts()
 #' p$plotfun <- stairs
 #' options(par_ts=p)
-#' plot_ts(Dbuildingheatload, "heatload", c("2010-12-15","2010-12-16"))
+#' plot_ts(Dbuilding, "heatload", c("2010-12-15","2010-12-16"))
 #'
 #' # Modify it to only lines
-#' plot_ts(Dbuildingheatload, "heatload", c("2010-12-15","2010-12-16"), plotfun=function(x,y,...){stairs(x,y, type="l")})
+#' plot_ts(Dbuilding, "heatload", c("2010-12-15","2010-12-16"), plotfun=function(x,y,...){stairs(x,y, type="l")})
 #'
 #' @export
 stairs <- function(x, y, type="b", preline=FALSE, pch=19, ...)
diff --git a/data/Dbuilding.rda b/data/Dbuilding.rda
new file mode 100644
index 0000000000000000000000000000000000000000..236ca3811d9dc73d7fd42dddd107493fd48a17a7
Binary files /dev/null and b/data/Dbuilding.rda differ
diff --git a/data/Dbuildingheatload.rda b/data/Dbuildingheatload.rda
deleted file mode 100644
index e1de9397212b07ab7cb7d149fa6159417648b100..0000000000000000000000000000000000000000
Binary files a/data/Dbuildingheatload.rda and /dev/null differ
diff --git a/data/Dsolarpower.rda b/data/Dsolarpower.rda
deleted file mode 100644
index e94a4e1a0732b848f55104e3700baf71899b283d..0000000000000000000000000000000000000000
Binary files a/data/Dsolarpower.rda and /dev/null differ
diff --git a/misc-R/data_soenderborg.csv b/data/all/data_soenderborg.csv
similarity index 100%
rename from misc-R/data_soenderborg.csv
rename to data/all/data_soenderborg.csv
diff --git a/data/all/make.R b/data/all/make.R
new file mode 100644
index 0000000000000000000000000000000000000000..12e787c1482f1284e98e28e76fe91a90b80984d6
--- /dev/null
+++ b/data/all/make.R
@@ -0,0 +1,89 @@
+# setting work directory and libraries #
+rm(list = ls())
+
+# # Packages used
+# require(R6)
+require(data.table)
+# require(Rcpp)
+# require(splines)
+library(devtools)
+library(roxygen2)
+
+pack <- as.package("../../../onlineforecast")
+load_all(pack)
+
+
+# Importing data # First unzip to get the .csv system('unzip
+# ../data/DataSoenderborg.zip')
+data_or <- fread("data_soenderborg.csv", sep = ",", header = TRUE)
+data_or[, `:=`(t, asct(data_or$t))]
+setDF(data_or)
+names(data_or)[names(data_or) == "Ig.obs"] <- "I.obs"
+
+# Make a data.table for each variable
+tmp <- unlist(getse(strsplit(names(data_or)[-1], "\\."), 2))
+colnm <- unlist(getse(strsplit(names(data_or)[-1], "\\."), 1))
+nms <- unique(colnm[grep("^k\\d*$", tmp)])
+
+kmax <- 48
+data <- list()
+data$t <- data_or$t
+for (ii in 1:length(nms)) {
+    # Find the columns
+    i <- grep(pst("^", nms[ii], "$"), unlist(getse(strsplit(names(data_or)[-1],"\\."), 1))) + 1
+    # Take only with kxx
+    i <- i[grep("k[[:digit:]]+$", names(data_or)[i])]
+    # 
+    #
+    data[[nms[ii]]] <- lag(data_or[ ,i], -1:-length(i))
+    names(data[[nms[ii]]]) <- pst("k", 1:length(i))
+    row.names(data[[nms[ii]]]) <- NULL
+    data[[nms[ii]]] <- as.data.frame(data[[nms[ii]]])
+    # Check if observed climate variable is there
+    nm <- pst(nms[[ii]], ".obs")
+    if (nm %in% names(data_or)) {
+        data[[nm]] <- data_or[, nm]
+    }
+}
+# More
+cols <- pst("Heatload.house", 1:16)
+data[["Heatload"]] <- data_or[, cols]
+names(data[["Heatload"]]) <- pst("house", 1:16)
+# 
+data[["cosAoi"]] <- data_or[, "cosAoi.obs"]
+data[["sunElevation"]] <- data_or[, "sunElevation.obs"]
+
+
+# # The time of day
+# ncol <- ncol(data$Ta)
+# tmp <- aslt(data$t)$hour
+# tmp <- matrix(tmp, nrow = length(tmp), ncol = ncol)
+# tmp <- data.frame(t(t(tmp) + (0:(ncol - 1))))
+# names(tmp) <- pst("k", 0:(ncol - 1))
+# data$tday <- tmp%%24
+
+#
+class(data) <- c("data.list","list")
+
+# Save for other scripts to read it
+#saveRDS(data, "data_soenderborg.RDS")
+
+
+data$heatloadtotal <- sapply(1:nrow(data$Heatload), function(i){
+    mean(unlist(data$Heatload[i, ]), na.rm=TRUE)
+})
+data$totaln <- sapply(1:nrow(data$Heatload), function(i){
+    sum(is.na(unlist(data$Heatload[i, ])))
+})
+plot(data$t, data$totaln)
+
+# Write for building heat load forecasting
+#Dbuilding <- subset(data, c("2010-12-15","2011-03-01"), nms=c("t","Heatload","Ta","I","Ws","Wd","Ta.obs","I.obs","Wd.obs","Ws.obs","cosAoi","sunElevation","tday"))
+data$heatload <- data$Heatload$house9
+Dbuilding <- subset(data, c("2010-12-15","2011-03-01"), nms=c("t","heatload","heatloadtotal","Ta.obs","I.obs","Ta","I"))
+rownames(Dbuilding$Ta) <- NULL
+Dbuilding$Ta <- Dbuilding$Ta[ ,1:36]
+rownames(Dbuilding$I) <- NULL
+Dbuilding$I <- Dbuilding$I[ ,1:36]
+#
+usethis::use_data(Dbuilding, overwrite=TRUE)
diff --git a/make.R b/make.R
index b9b2e60552fdba2ed36cf288a372228225f6ccc1..22f3db399cb8560f07e976c6bd51854681aab0ad 100644
--- a/make.R
+++ b/make.R
@@ -51,10 +51,10 @@ docit()
 
 # ----------------------------------------------------------------
 # Build the package (remember to rebuild vignettes for release)
-build(".", vignettes=TRUE)
+build(".", vignettes=FALSE)
 
 # Install it
-install.packages("../onlineforecast_0.1.0.tar.gz")
+install.packages("../onlineforecast_1.0.0.tar.gz")
 
 library(onlineforecast)
 
diff --git a/misc-R/.Rhistory~ b/misc-R/.Rhistory~
deleted file mode 100644
index 5503199be76081abc656fc99a3d14c4f0dbe7f1e..0000000000000000000000000000000000000000
--- a/misc-R/.Rhistory~
+++ /dev/null
@@ -1,512 +0,0 @@
-Dnew_transformed$AR.Order1
-modelAR$yAR
-rls_update(modelAR, Dnew_transformed, Dnew[[model$output]])
-yhatAR <- rls_predict(modelAR, Dnew_transformed)
-modelAR$yAR
-(i <- iseq[length(iseq)] + 2)
-Dnew <- subset(D, i)
-Dnew_transformed <- model$transform_data(Dnew)
-rls_update(model, Dnew_transformed, Dnew[[model$output]])
-yhat <- rls_predict(model, Dnew_transformed)
-iseq <- i+modelAR$kseq
-plot(D$t[iseq], D$y[iseq], type = "b", xlab = "t", ylab = "y")
-lines(D$t[iseq], yhat[1,], type = "b", col = 2)
-lines(D$t[iseq], yhatAR[1,], type = "b", col = 4)
-legend("topright", c("observations",pst("predictions (",min(modelAR$kseq)," to ",max(modelAR$kseq)," steps ahead)"),
-pst("Predictions AR (",min(model$kseq)," to ",max(model$kseq)," steps ahead)")), lty = 1, col = c(1,2,4))
-Ex1data <- data.frame(x = c(71,74,82,76,91,82,82,75,79,82,72,90))
-10:90
-ObsVar <- var(Ex1data$x)
-ObsVar
-1/sqrt(2*pi*ObsVar) * exp(-(71 - 70)^2/(2*ObsVar))
-1/sqrt(2*pi*ObsVar) * exp(-(72 - 70)^2/(2*ObsVar))
-1/sqrt(2*pi*ObsVar) * exp(-(c(71,72) - 70)^2/(2*ObsVar))
-1/sqrt(2*pi*ObsVar) * exp(-(Ex1data$x - 70)^2/(2*ObsVar))
-prod(1/sqrt(2*pi*ObsVar) * exp(-(Ex1data$x - 70)^2/(2*ObsVar)))
-LikelihoodAll
-ObsVar <- var(Ex1data$x)
-LikelihoodAll <- sapply(70:91, function(mean) prod(1/sqrt(2*pi*ObsVar) * exp(-(Ex1data$x - mean)^2/(2*ObsVar))))
-LikelihoodAll
-plot(LikelihoodAll)
-plot(70:91LikelihoodAll)
-plot(70:91, LikelihoodAll)
-ObsVar <- 1#var(Ex1data$x)
-test <- c(4.6,6.3,5.0)
-LikelihoodAll <- sapply(4:6.5, function(mean) prod(1/sqrt(2*pi*ObsVar) * exp(-(test - mean)^2/(2*ObsVar))))
-plot(4:6.5, LikelihoodAll)
-LikelihoodAll <- sapply(seq(4,6.5,0.1), function(mean) prod(1/sqrt(2*pi*ObsVar) * exp(-(test - mean)^2/(2*ObsVar))))
-plot(seq(4,6.5,0.1), LikelihoodAll)
-LikelihoodAll <- sapply(seq(70,95,0.1), function(mean) prod(1/sqrt(2*pi*ObsVar) * exp(-(Ex1data$x - mean)^2/(2*ObsVar))))
-Seq <- seq(70,95,0.1
-LikelihoodAll <- sapply(), function(mean) prod(1/sqrt(2*pi*ObsVar) * exp(-(Ex1data$x - mean)^2/(2*ObsVar))))
-Seq <- seq(70,95,0.1)
-LikelihoodAll <- sapply(Seq, function(mean) prod(1/sqrt(2*pi*ObsVar) * exp(-(Ex1data$x - mean)^2/(2*ObsVar))))
-plot(Seq, LikelihoodAll)
-Seq <- seq(75,85,0.1)
-LikelihoodAll <- sapply(Seq, function(mean) prod(1/sqrt(2*pi*ObsVar) * exp(-(Ex1data$x - mean)^2/(2*ObsVar))))
-plot(Seq, LikelihoodAll)
-L.complete.data <- function(theta) prod(dnorm(Ex1data$x, mean = theta, sd = sqrt(s2)))
-plot(Seq, LikelihoodAll)
-L.complete.data <- function(theta) prod(dnorm(Ex1data$x, mean = theta, sd = sqrt(ObsVar)))
-th <- seq(mean(Ex1data$x) - 3*sqrt(ObsVar), mean(Ex1data$x) + 3* sqrt(ObsVar), length = 200)
-L <- sapply(th, L.complete.data)
-plot(th, L)
-plot(Seq, LikelihoodAll)
-plot(th, L)
-plot(Seq, LikelihoodAll)
-plot(th, L)
-plot(th, L/max(L), ylab = "L", xlab = expression(theta))
-L.ave <- function(theta) dnorm(mean(Ex1data$x), mean = theta, sd = sqrt(ObsVar/length(Ex1data$x)))
-LAverage <- sapply(th, L.ave)
-lines(th, LAverage/max(LAverage), col = "red")
-mean(Ex1data$x)
-abline(v = mean(Ex1data$x), col = "blue")
-mle.estimate <- function(th) -sum(dnorm(Ex1data$x, mean = theta, sd = sqrt(ObsVar), log = TRUE))
-fit <- optiom(par = mean(Ex1data$x),
-fn = mle.estimate,
-hessian = TRUE)
-fit <- optim(par = mean(Ex1data$x),
-fn = mle.estimate,
-hessian = TRUE)
-?optim
-mle.estimate <- function(th) -sum(dnorm(Ex1data$x, mean = th, sd = sqrt(ObsVar), log = TRUE))
-fit <- optim(par = mean(Ex1data$x),
-fn = mle.estimate,
-hessian = TRUE)
-mle.estimate <- function(th) -sum(dnorm(Ex1data$x, mean = th, sd = sqrt(ObsVar), log = TRUE))
-fit <- optim(par = mean(Ex1data$x),
-fn = mle.estimate,
-hessian = TRUE,
-method = "Brent")
-fit <- optim(par = mean(Ex1data$x),
-fn = mle.estimate,
-hessian = TRUE,
-method = "Brent",
-lower = min(Ex1data$x),
-upper = max(Ex1data$x))
-print(fit)
-fit <- optim(par = mean(Ex1data$x),
-fn = mle.estimate,
-hessian = TRUE,
-method = "Brent",
-lower = min(Ex1data$x),
-upper = max(Ex1data$x))
-print(fit)
-fit <- optim(par = mean(Ex1data$x),
-fn = mle.estimate,
-hessian = TRUE,
-method = "Brent",
-lower = min(Ex1data$x),
-upper = max(Ex1data$x))
-print(fit)
-fit$hessian
-1/fit$hessian
-ObsVar
-1
-ObsVar <- var(Ex1data$x)
-mle.estimate <- function(th) -sum(dnorm(Ex1data$x, mean = th, sd = sqrt(ObsVar), log = TRUE))
-fit <- optim(par = mean(Ex1data$x),
-fn = mle.estimate,
-hessian = TRUE,
-method = "Brent",
-lower = min(Ex1data$x),
-upper = max(Ex1data$x))
-print(fit)
-1/fit$hessian
-ObsVar
-print(solve(fit$hessian))
-print(ObsVar)
-print(solve(fit$hessian))
-print(sqrt(ObsVar))
-print(solve(fit$hessian))
-print(ObsVar/length(Ex1data$x))
-print(solve(fit$hessian))
-Ex3data <- data.frame(x = c(4,6,3,7,2,4))
-?dpois
-lambdaseq = seq(-4,0,0.1)
-Posres <- optim(par = lambdaseq,
-fn = PosLogLikelihood)
-PosLogLikelihood <- function(lambda) -sum(dpois(x = Ex3data$x, lambda = lambda, log = TRUE))
-lambdaseq = seq(-4,0,0.1)
-Posres <- optim(par = lambdaseq,
-fn = PosLogLikelihood)
-lambdaseq = seq(-4,0,0.1)
-lambdaseq
-Posres <- optim(par = lambdaseq,
-fn = PosLogLikelihood,
-)
-Posres
-plot(lambdaseq, PosLogLikelihood(lambdaseq))
-dpois(x = Ex3data$x, lambda = 0, log = TRUE)
-Ex3data
-PosLogLikelihood <- function(lambda) sum(dpois(x = Ex3data$x, lambda = lambda, log = TRUE))
-lambdaseq = seq(-4,0,0.1)
-plot(lambdaseq, PosLogLikelihood(lambdaseq))
-lambdaseq = seq(1,8,0.1)
-plot(lambdaseq, PosLogLikelihood(lambdaseq))
-PosLogLikelihood <- function(lambda) sum(dpois(x = Ex3data$x, lambda = lambda, log = TRUE))
-lambdaseq = seq(1,8,0.1)
-plot(lambdaseq, PosLogLikelihood(lambdaseq))
-lambdaseq
-PosLogLikelihood(lambdaseq)
-plot(lambdaseq, sapply(lambdaseq, PosLogLikelihood))
-L <- sapply(lambdaseq, PosLogLikelihood)
-plot(lambdaseq, L - max(L))
-lambdaseq = seq(2,7,0.1)
-L <- sapply(lambdaseq, PosLogLikelihood)
-plot(lambdaseq, L - max(L))
-lambdaseq = seq(2.5,7,0.1)
-L <- sapply(lambdaseq, PosLogLikelihood)
-plot(lambdaseq, L - max(L))
-lambdaseq = seq(2.3,7,0.1)
-L <- sapply(lambdaseq, PosLogLikelihood)
-plot(lambdaseq, L - max(L))
-lambdaseq = seq(2.3,7.4,0.1)
-L <- sapply(lambdaseq, PosLogLikelihood)
-plot(lambdaseq, L - max(L))
-lambdaseq = seq(2.3,7.3,0.1)
-L <- sapply(lambdaseq, PosLogLikelihood)
-plot(lambdaseq, L - max(L))
-L.quad <- sapply(lambdaseq, l.quad)
-l.quad <- function(lambda) dpois(x = mean(Ex3data$x), lambda = lambda, log = TRUE) - 0.5*length(Ex3data$x)/mean(Ex3data$x) * (lambda - mean(Ex3data$x))^2
-L.quad <- sapply(lambdaseq, l.quad)
-warnings()
-L.quad
-l.quad <- function(lambda) dpois(x = Ex3data$x, lambda =  mean(Ex3data$x), log = TRUE) - 0.5*length(Ex3data$x)/mean(Ex3data$x) * (lambda - mean(Ex3data$x))^2
-L.quad <- sapply(lambdaseq, l.quad)
-lines(lambdaseq, l.quad, col = "red")
-lines(lambdaseq, L.quad, col = "red")
-L.quad
-l.quad <- function(lambda) sum(dpois(x = Ex3data$x, lambda =  mean(Ex3data$x), log = TRUE)) - 0.5*length(Ex3data$x)/mean(Ex3data$x) * (lambda - mean(Ex3data$x))^2
-L.quad <- sapply(lambdaseq, l.quad)
-lines(lambdaseq, L.quad, col = "red")
-L.quad
-lines(lambdaseq, L.quad - max(L.quad), col = "red")
-PosLogLikelihood <- function(lambda) sum(dpois(x = Ex3data$x, lambda = lambda, log = TRUE))
-lambdaseq = seq(2,7.3,0.1)
-L <- sapply(lambdaseq, PosLogLikelihood)
-plot(lambdaseq, L - max(L))
-l.quad <- function(lambda) sum(dpois(x = Ex3data$x, lambda =  mean(Ex3data$x), log = TRUE)) - 0.5*length(Ex3data$x)/mean(Ex3data$x) * (lambda - mean(Ex3data$x))^2
-L.quad <- sapply(lambdaseq, l.quad)
-lines(lambdaseq, L.quad - max(L.quad), col = "red")
-PosLogLikelihood <- function(lambda) sum(dpois(x = Ex3data$x, lambda = lambda, log = TRUE))
-lambdaseq = seq(2.3,7.3,0.1)
-L <- sapply(lambdaseq, PosLogLikelihood)
-plot(lambdaseq, L - max(L))
-l.quad <- function(lambda) sum(dpois(x = Ex3data$x, lambda =  mean(Ex3data$x), log = TRUE)) - 0.5*length(Ex3data$x)/mean(Ex3data$x) * (lambda - mean(Ex3data$x))^2
-L.quad <- sapply(lambdaseq, l.quad)
-lines(lambdaseq, L.quad - max(L.quad), col = "red")
-plot(lambdaseq, exp(L - max(L)))
-lines(lambdaseq, exp(L.quad - max(L.quad)), col = "red")
-par(mfrow = c(1,2))
-plot(lambdaseq, L - max(L))
-lines(lambdaseq, L.quad - max(L.quad), col = "red")
-plot(lambdaseq, exp(L - max(L)))
-lines(lambdaseq, exp(L.quad - max(L.quad)), col = "red")
-## Demo on how the AR part works in onlineforecast
-rm(list = ls())
-library(devtools)
-load_all(as.package("../../onlineforecast"))
-class(Dbuildingheatload)
-D <- Dbuildingheatload
-D$y <- D$Heatload$house9
-plot_ts(D, c("^y","Ta"), kseq=c(0,12))
-plot_ts(D, c("^y","Ta"), "2010-12-15", "2010-12-25", kseq=c(0,12))
-Dtrain <- subset(D, c("2010-12-15", "2011-01-01"))
-Dtrain$fiteval <- period("2010-12-20", Dtrain$t)
-model <- forecastmodel$new()
-model$output = "y"
-model$add_inputs(Ta = "lp(Ta, a1=0.9)",
-I = "lp(I, a1=0.7)",
-mu_tday = "fs(tday/24, nharmonics=10)",
-mu = "ones()")
-model$add_regp("rls_prm(lambda=0.9)")
-model$add_pb(Ta__a1 = c(0.8, 0.9, 0.9999),
-I__a1 =  c(0.4, 0.8, 0.9999),
-lambda = c(0.9, 0.99, 0.9999))
-model$kseq <- c(1,18)
-model$p <- rls_optim(model, Dtrain, control=list(maxit=2), cachedir = "cache-building-heat-load-forecasting")$par
-model$kseq <- 1:36
-val <- rls_fit(model$p, model, D, returnanalysis = TRUE)
-D$Yhat <- val$Yhat
-i <- 200
-iseq <- i+model$kseq
-plot(D$t[iseq], D$y[iseq], type = "b", xlab = "t", ylab = "y")
-lines(D$t[iseq], D$Yhat[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)
-iseq <- which(period("2010-12-15",D$t,"2011-01-01"))
-Dfit <- subset(D, iseq)
-rls_fit(model$p, model, Dfit)
-(i <- iseq[length(iseq)] + 1)
-Dnew <- subset(D, i)
-Dnew_transformed <- model$transform_data(Dnew)
-rls_update(model, Dnew_transformed, Dnew[[model$output]])
-yhat <- rls_predict(model, Dnew_transformed)
-## AR part
-#source("AR.R")
-modelAR <- forecastmodel$new()
-modelAR$output = "y"
-modelAR$add_inputs(Ta = "lp(Ta, a1=0.9)",
-I = "lp(I, a1=0.7)",
-mu_tday = "fs(tday/24, nharmonics=10)",
-mu = "ones()",
-AR = "AR(lags=c(0))")
-modelAR$add_regp("rls_prm(lambda=0.9)")
-modelAR$add_pb(Ta__a1 = c(0.8, 0.9, 0.9999),
-I__a1 =  c(0.4, 0.8, 0.9999),
-lambda = c(0.9, 0.99, 0.9999))
-modelAR$kseq <- c(1,18)
-modelAR$p <- rls_optim(modelAR, Dtrain, control=list(maxit=2), cachedir = "cache-building-heat-load-forecasting")$par
-modelAR$kseq <- 1:36
-valAR <- rls_fit(modelAR$p, modelAR, D, returnanalysis = TRUE)
-D$YhatAR <- valAR$Yhat
-plot_ts(D, c("^y|^Y"), "2011-01-01", "2011-02-01", kseq = c(1,18))
-plot_ts(D, c("^y|^Y"), "2011-01-01", "2011-01-10", kseq = c(1,18))
-i <- 200
-iseq <- i+modelAR$kseq
-plot(D$t[iseq], D$y[iseq], type = "b", xlab = "t", ylab = "y")
-lines(D$t[iseq], D$Yhat[i, ], type = "b", col = 2)
-lines(D$t[iseq], D$YhatAR[i, ], type = "b", col = 4)
-legend("topright", c("Observations",pst("Predictions (",min(model$kseq)," to ",max(model$kseq)," steps ahead)"),
-pst("Predictions AR (",min(model$kseq)," to ",max(model$kseq)," steps ahead)")), lty = 1, col = c(1,2,4))
-sqrt(mean( as.numeric((D$y[iseq] -  D$Yhat[i, ])^2), na.rm = T))
-sqrt(mean( as.numeric((D$y[iseq] -  D$YhatAR[i, ])^2), na.rm = T))
-## Update
-## THINKABOUT WITH bigger lag and more new observation!!!
-model <- forecastmodel$new()
-model$output = "y"
-modelAR$add_inputs(Ta = "lp(Ta, a1=0.9)",
-I = "lp(I, a1=0.7)",
-mu_tday = "fs(tday/24, nharmonics=10)",
-mu = "ones()",
-AR = "AR(lags=c(0,1,4))")
-modelAR$add_regp("rls_prm(lambda=0.9)")
-modelAR$add_pb(Ta__a1 = c(0.8, 0.9, 0.9999),
-I__a1 =  c(0.4, 0.8, 0.9999),
-lambda = c(0.9, 0.99, 0.9999))
-modelAR$kseq <- c(1,18)
-iseq <- which(period("2010-12-15",D$t,"2011-01-01"))
-Dfit <- subset(D, iseq)
-rls_fit(modelAR$p, modelAR, Dfit)
-str(modelAR$Lfits[1:2])
-(i <- iseq[length(iseq)] + 1)
-Dnew <- subset(D, i:(i))
-Dnew$y
-modelAR$yAR
-D$y[(i-6):i]
-Dnew_transformed <- modelAR$transform_data(Dnew)
-Dnew_transformed$AR.lag0
-Dnew_transformed$AR.lag1
-Dnew_transformed$AR.lag4
-## modelAR$yAR
-## rls_update(modelAR, Dnew_transformed, Dnew[[model$output]])
-## yhatAR <- rls_predict(modelAR, Dnew_transformed)
-## modelAR$yAR
-## iseq <- i+modelAR$kseq
-## plot(D$t[iseq], D$y[iseq], type = "b", xlab = "t", ylab = "y")
-## lines(D$t[iseq], yhat[1,], type = "b", col = 2)
-## lines(D$t[iseq], yhatAR[1,], type = "b", col = 4)
-## legend("topright", c("observations",pst("predictions (",min(modelAR$kseq)," to ",max(modelAR$kseq)," steps ahead)"),
-##                      pst("Predictions AR (",min(model$kseq)," to ",max(model$kseq)," steps ahead)")), lty = 1, col = c(1,2,4))
-## (i <- iseq[length(iseq)] + 2)
-## Dnew <- subset(D, i:(i))
-## Dnew$y
-## modelAR$yAR
-## tail(Dfit$y)
-## Dnew_transformed <- modelAR$transform_data(Dnew)
-## Dnew_transformed$AR.lag3
-## Dnew_transformed$AR.Order2
-## Dnew_transformed$AR.Order1
-## modelAR$yAR
-## rls_update(modelAR, Dnew_transformed, Dnew[[model$output]])
-## yhatAR <- rls_predict(modelAR, Dnew_transformed)
-## modelAR$yAR
-## (i <- iseq[length(iseq)] + 2)
-## Dnew <- subset(D, i)
-## Dnew_transformed <- model$transform_data(Dnew)
-## rls_update(model, Dnew_transformed, Dnew[[model$output]])
-## yhat <- rls_predict(model, Dnew_transformed)
-## iseq <- i+modelAR$kseq
-## plot(D$t[iseq], D$y[iseq], type = "b", xlab = "t", ylab = "y")
-## lines(D$t[iseq], yhat[1,], type = "b", col = 2)
-## lines(D$t[iseq], yhatAR[1,], type = "b", col = 4)
-## legend("topright", c("observations",pst("predictions (",min(modelAR$kseq)," to ",max(modelAR$kseq)," steps ahead)"),
-##                      pst("Predictions AR (",min(model$kseq)," to ",max(model$kseq)," steps ahead)")), lty = 1, col = c(1,2,4))
-## Demo on how the AR part works in onlineforecast
-rm(list = ls())
-library(devtools)
-load_all(as.package("../../onlineforecast"))
-load_all(as.package("../onlineforecast"))
-library(onlineforecast)
-load_all(as.package("../../onlineforecast"))
-library(devtools)
-load_all(as.package("../../onlineforecast"))
-load_all(as.package("../onlineforecast"))
-setwd("~/Documents/phd/Projects/Git/onlineforecast/misc-R")
-load_all(as.package("../onlineforecast"))
-load_all(as.package("../../onlineforecast"))
-## Demo on how the AR part works in onlineforecast
-rm(list = ls())
-library(devtools)
-load_all(as.package("../../onlineforecast"))
-class(Dbuildingheatload)
-D <- Dbuildingheatload
-D$y <- D$Heatload$house9
-plot_ts(D, c("^y","Ta"), kseq=c(0,12))
-plot_ts(D, c("^y","Ta"), "2010-12-15", "2010-12-25", kseq=c(0,12))
-Dtrain <- subset(D, c("2010-12-15", "2011-01-01"))
-Dtrain$fiteval <- period("2010-12-20", Dtrain$t)
-model <- forecastmodel$new()
-model$output = "y"
-model$add_inputs(Ta = "lp(Ta, a1=0.9)",
-I = "lp(I, a1=0.7)",
-mu_tday = "fs(tday/24, nharmonics=10)",
-mu = "ones()")
-model$add_regp("rls_prm(lambda=0.9)")
-model$add_pb(Ta__a1 = c(0.8, 0.9, 0.9999),
-I__a1 =  c(0.4, 0.8, 0.9999),
-lambda = c(0.9, 0.99, 0.9999))
-model$kseq <- c(1,18)
-model$p <- rls_optim(model, Dtrain, control=list(maxit=2), cachedir = "cache-building-heat-load-forecasting")$par
-model$kseq <- 1:36
-val <- rls_fit(model$p, model, D, returnanalysis = TRUE)
-D$Yhat <- val$Yhat
-i <- 200
-iseq <- i+model$kseq
-plot(D$t[iseq], D$y[iseq], type = "b", xlab = "t", ylab = "y")
-lines(D$t[iseq], D$Yhat[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)
-iseq <- which(period("2010-12-15",D$t,"2011-01-01"))
-Dfit <- subset(D, iseq)
-rls_fit(model$p, model, Dfit)
-(i <- iseq[length(iseq)] + 1)
-Dnew <- subset(D, i)
-Dnew_transformed <- model$transform_data(Dnew)
-rls_update(model, Dnew_transformed, Dnew[[model$output]])
-yhat <- rls_predict(model, Dnew_transformed)
-## AR part
-#source("AR.R")
-modelAR <- forecastmodel$new()
-modelAR$output = "y"
-modelAR$add_inputs(Ta = "lp(Ta, a1=0.9)",
-I = "lp(I, a1=0.7)",
-mu_tday = "fs(tday/24, nharmonics=10)",
-mu = "ones()",
-AR = "AR(lags=c(0))")
-modelAR$add_regp("rls_prm(lambda=0.9)")
-modelAR$add_pb(Ta__a1 = c(0.8, 0.9, 0.9999),
-I__a1 =  c(0.4, 0.8, 0.9999),
-lambda = c(0.9, 0.99, 0.9999))
-modelAR$kseq <- c(1,18)
-modelAR$p <- rls_optim(modelAR, Dtrain, control=list(maxit=2), cachedir = "cache-building-heat-load-forecasting")$par
-modelAR$kseq <- 1:36
-valAR <- rls_fit(modelAR$p, modelAR, D, returnanalysis = TRUE)
-D$YhatAR <- valAR$Yhat
-plot_ts(D, c("^y|^Y"), "2011-01-01", "2011-02-01", kseq = c(1,18))
-plot_ts(D, c("^y|^Y"), "2011-01-01", "2011-01-10", kseq = c(1,18))
-i <- 200
-iseq <- i+modelAR$kseq
-plot(D$t[iseq], D$y[iseq], type = "b", xlab = "t", ylab = "y")
-lines(D$t[iseq], D$Yhat[i, ], type = "b", col = 2)
-lines(D$t[iseq], D$YhatAR[i, ], type = "b", col = 4)
-legend("topright", c("Observations",pst("Predictions (",min(model$kseq)," to ",max(model$kseq)," steps ahead)"),
-pst("Predictions AR (",min(model$kseq)," to ",max(model$kseq)," steps ahead)")), lty = 1, col = c(1,2,4))
-sqrt(mean( as.numeric((D$y[iseq] -  D$Yhat[i, ])^2), na.rm = T))
-sqrt(mean( as.numeric((D$y[iseq] -  D$YhatAR[i, ])^2), na.rm = T))
-## Update
-## THINKABOUT WITH bigger lag and more new observation!!!
-model <- forecastmodel$new()
-model$output = "y"
-modelAR$add_inputs(Ta = "lp(Ta, a1=0.9)",
-I = "lp(I, a1=0.7)",
-mu_tday = "fs(tday/24, nharmonics=10)",
-mu = "ones()",
-AR = "AR(lags=c(0,1,4))")
-modelAR$add_regp("rls_prm(lambda=0.9)")
-modelAR$add_pb(Ta__a1 = c(0.8, 0.9, 0.9999),
-I__a1 =  c(0.4, 0.8, 0.9999),
-lambda = c(0.9, 0.99, 0.9999))
-modelAR$kseq <- c(1,18)
-iseq <- which(period("2010-12-15",D$t,"2011-01-01"))
-Dfit <- subset(D, iseq)
-rls_fit(modelAR$p, modelAR, Dfit)
-str(modelAR$Lfits[1:2])
-(i <- iseq[length(iseq)] + 1)
-Dnew <- subset(D, i:(i))
-Dnew$y
-modelAR$yAR
-D$y[(i-6):i]
-Dnew_transformed <- modelAR$transform_data(Dnew)
-Dnew_transformed$AR.lag0
-Dnew_transformed$AR.lag1
-Dnew_transformed$AR.lag4
-## modelAR$yAR
-## rls_update(modelAR, Dnew_transformed, Dnew[[model$output]])
-## yhatAR <- rls_predict(modelAR, Dnew_transformed)
-## modelAR$yAR
-## iseq <- i+modelAR$kseq
-## plot(D$t[iseq], D$y[iseq], type = "b", xlab = "t", ylab = "y")
-## lines(D$t[iseq], yhat[1,], type = "b", col = 2)
-## lines(D$t[iseq], yhatAR[1,], type = "b", col = 4)
-## legend("topright", c("observations",pst("predictions (",min(modelAR$kseq)," to ",max(modelAR$kseq)," steps ahead)"),
-##                      pst("Predictions AR (",min(model$kseq)," to ",max(model$kseq)," steps ahead)")), lty = 1, col = c(1,2,4))
-## (i <- iseq[length(iseq)] + 2)
-## Dnew <- subset(D, i:(i))
-## Dnew$y
-## modelAR$yAR
-## tail(Dfit$y)
-## Dnew_transformed <- modelAR$transform_data(Dnew)
-## Dnew_transformed$AR.lag3
-## Dnew_transformed$AR.Order2
-## Dnew_transformed$AR.Order1
-## modelAR$yAR
-## rls_update(modelAR, Dnew_transformed, Dnew[[model$output]])
-## yhatAR <- rls_predict(modelAR, Dnew_transformed)
-## modelAR$yAR
-## (i <- iseq[length(iseq)] + 2)
-## Dnew <- subset(D, i)
-## Dnew_transformed <- model$transform_data(Dnew)
-## rls_update(model, Dnew_transformed, Dnew[[model$output]])
-## yhat <- rls_predict(model, Dnew_transformed)
-## iseq <- i+modelAR$kseq
-## plot(D$t[iseq], D$y[iseq], type = "b", xlab = "t", ylab = "y")
-## lines(D$t[iseq], yhat[1,], type = "b", col = 2)
-## lines(D$t[iseq], yhatAR[1,], type = "b", col = 4)
-## legend("topright", c("observations",pst("predictions (",min(modelAR$kseq)," to ",max(modelAR$kseq)," steps ahead)"),
-##                      pst("Predictions AR (",min(model$kseq)," to ",max(model$kseq)," steps ahead)")), lty = 1, col = c(1,2,4))
-valAR$Yhat
-head(D)
-head(D$t)
-dim(valAR$Yhat)
-dummyDF <- data.frame(PredTime = D$t[1],
-Time = D$t[2:dim(valAR$Yhat)[2]])
-dummyDF
-dummyDF <- data.frame(PredTime = D$t[1],
-Time = D$t[2:dim(valAR$Yhat)[2]],
-k = 1:(dim(valAR$Yhat)[2]))
-dummyDF
-predictionDataFrame <- function()
-{
-head(D$t)
-dummyDF <- data.frame(PredTime = D$t[1],
-Time = D$t[2:dim(valAR$Yhat)[2]],
-k = 1:(dim(valAR$Yhat)[2]))
-dummyDF
-}
-dummyDF <- data.frame(PredTime = D$t[1],
-Time = D$t[2:dim(valAR$Yhat)[2]],
-k = 1:(dim(valAR$Yhat)[2]))
-dummyDF <- data.frame(PredTime = D$t[1],
-Time = D$t[2:(dim(valAR$Yhat)[2]+1)],
-k = 1:(dim(valAR$Yhat)[2]))
-dummyDF
-dummyDF <- data.frame(PredTime = D$t[1],
-Time = D$t[2:(dim(valAR$Yhat)[2]+1)],
-k = 1:(dim(valAR$Yhat)[2]),
-Pred = valAR$Yhat[1,])
-valAR$Yhat[1,]
-dummyDF <- data.frame(PredTime = D$t[1],
-Time = D$t[2:(dim(valAR$Yhat)[2]+1)],
-k = 1:(dim(valAR$Yhat)[2]),
-Pred = as.numeric(valAR$Yhat[1,]))
-dummyDF
diff --git a/misc-R/arfunction.R b/misc-R/arfunction.R
deleted file mode 100644
index 296ec2cbed2842390bd79b8edb7f31cf70714739..0000000000000000000000000000000000000000
--- a/misc-R/arfunction.R
+++ /dev/null
@@ -1,190 +0,0 @@
-## Demo on how the AR part works in onlineforecast
-rm(list = ls())
-library(devtools)
-load_all(as.package("../../onlineforecast"))
-
-class(Dbuildingheatload)
-D <- Dbuildingheatload
-D$y <- D$heatload
-
-plot_ts(D, c("^y","Ta"), kseq=c(0,12))
-plot_ts(D, c("^y","Ta"), "2010-12-15", "2010-12-25", kseq=c(0,12))
-
-Dtrain <- subset(D, c("2010-12-15", "2011-01-01"))
-Dtrain$scoreperiod <- in_range("2010-12-20", Dtrain$t)
-model <- forecastmodel$new()
-model$output = "y"
-model$add_inputs(Ta = "lp(Ta, a1=0.9)", 
-                 I = "lp(I, a1=0.7)", 
-                 mu_tday = "fs(tday/24, nharmonics=10)",
-                 mu = "ones()")
-model$add_regprm("rls_prm(lambda=0.9)")
-model$add_prmbounds(Ta__a1 = c(0.8, 0.9, 0.9999),
-                    I__a1 =  c(0.4, 0.8, 0.9999),
-                    lambda = c(0.9, 0.99, 0.9999))
-model$kseq <- c(1,18)
-model$prm <- rls_optim(model, Dtrain, control=list(maxit=2), cachedir = "cache-building-heat-load-forecasting")$par
-model$kseq <- 1:36
-val <- rls_fit(model$prm, model, D, returnanalysis = TRUE)
-D$Yhat <- val$Yhat
-
-i <- 200
-iseq <- i+model$kseq
-plot(D$t[iseq], D$y[iseq], type = "b", xlab = "t", ylab = "y")
-lines(D$t[iseq], D$Yhat[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)
-)
-iseq <- which(in_range("2010-12-15",D$t,"2011-01-01"))
-Dfit <- subset(D, iseq)
-rls_fit(model$prm, model, Dfit, returnanalysis = FALSE)
-(i <- iseq[length(iseq)] + 1)
-Dnew <- subset(D, i)
-Dnew_transformed <- model$transform_data(Dnew)
-rls_update(model, Dnew_transformed, Dnew[[model$output]])
-yhat <- rls_predict(model, Dnew_transformed)
-
-iseq <- i+model$kseq
-plot(D$t[iseq], D$y[iseq], type = "b", xlab = "t", ylab = "y")
-lines(D$t[iseq], yhat, type = "b", col = 2)
-legend("topright", c("observations",pst("predictions (",min(model$kseq)," to ",max(model$kseq)," steps ahead)")), lty = 1, col = 1:2)
-
-## AR part
-modelAR <- forecastmodel$new()
-modelAR$output = "y"
-modelAR$add_inputs(Ta = "lp(Ta, a1=0.9)", 
-                   I = "lp(I, a1=0.7)", 
-                   mu_tday = "fs(tday/24, nharmonics=10)",
-                   mu = "ones()", 
-                   AR = "AR(lags=c(0))")
-modelAR$add_regprm("rls_prm(lambda=0.9)")
-modelAR$add_prmbounds(Ta__a1 = c(0.8, 0.9, 0.9999),
-                      I__a1 =  c(0.4, 0.8, 0.9999),
-                      lambda = c(0.9, 0.99, 0.9999))
-modelAR$kseq <- c(1,18)
-
-modelAR$prm <- rls_optim(modelAR, Dtrain, control=list(maxit=2), cachedir = "cache-building-heat-load-forecasting")$par
-modelAR$kseq <- 1:36
-valAR <- rls_fit(modelAR$prm, modelAR, D, returnanalysis = TRUE)
-names(valAR)
-names(modelAR)
-
-D$YhatAR <- valAR$Yhat
-plot_ts(D, c("^y|^Y"), "2011-01-01", "2011-02-01", kseq = c(1,18))
-plot_ts(D, c("^y|^Y"), "2011-01-01", "2011-01-10", kseq = c(1,18))
-i <- 200
-iseq <- i+modelAR$kseq
-plot(D$t[iseq], D$y[iseq], type = "b", xlab = "t", ylab = "y")
-lines(D$t[iseq], D$Yhat[i, ], type = "b", col = 2)
-lines(D$t[iseq], D$YhatAR[i, ], type = "b", col = 4)
-legend("topright", c("Observations",pst("Predictions (",min(model$kseq)," to ",max(model$kseq)," steps ahead)"),
-                     pst("Predictions AR (",min(model$kseq)," to ",max(model$kseq)," steps ahead)")), lty = 1, col = c(1,2,4))
-
-sqrt(mean( as.numeric((D$y[iseq] -  D$Yhat[i, ])^2), na.rm = T))
-sqrt(mean( as.numeric((D$y[iseq] -  D$YhatAR[i, ])^2), na.rm = T))
-
-## Prediction Data Frame
-PredictionDF <- long_format(fit = valAR, Time = D$t) # long_format(fit = valAR, Time = D$t)
-head(PredictionDF[2000:2500,], 50)
-head(PredictionDF[(36*36):(36*36+36),],40)
-
-
-i <- 200
-iseq <- i+modelAR$kseq
-plot(D$t[iseq], D$y[iseq], type = "b", xlab = "t", ylab = "y")
-lines(D$t[iseq], D$Yhat[i, ], type = "b", col = 2)
-lines(PredictionDF[PredictionDF$PredTime == D$t[i],]$Time, PredictionDF[PredictionDF$PredTime == D$t[i],]$Pred, type = "b", col = 5, cex = 0.4)
-lines(D$t[iseq], D$YhatAR[i, ], type = "b", col = 4)
-legend("topright", c("Observations",pst("Predictions (",min(model$kseq)," to ",max(model$kseq)," steps ahead)"),
-                     pst("Predictions AR (",min(model$kseq)," to ",max(model$kseq)," steps ahead)")), lty = 1, col = c(1,2,4)
-)
-
-## Selecting just the one hour forecasts
-OneHourPred <- PredictionDF[PredictionDF$k == 1,]
-length(OneHourPred$PredTime)
-head(OneHourPred)
-length(valAR$Yhat$k1)
-head(valAR$Yhat$k1)
-
-plot(OneHourPred$Time, OneHourPred$Pred, type = "l")
-lines(D$t, D$y, col = "red")
-
-## Update
-
-## THINKABOUT WITH bigger lag and more new observation!!!
-
-model <- forecastmodel$new()
-model$output = "y"
-modelAR$add_inputs(Ta = "lp(Ta, a1=0.9)", 
-                   I = "lp(I, a1=0.7)", 
-                   mu_tday = "fs(tday/24, nharmonics=10)",
-                   mu = "ones()", 
-                   AR = "AR(lags=c(0,1,4))")
-modelAR$add_regprm("rls_prm(lambda=0.9)")
-modelAR$add_prmbounds(Ta__a1 = c(0.8, 0.9, 0.9999),
-                      I__a1 =  c(0.4, 0.8, 0.9999),
-                      lambda = c(0.9, 0.99, 0.9999))
-modelAR$kseq <- c(1,18)
-
-
-iseq <- which(in_range("2010-12-15",D$t,"2011-01-01"))
-Dfit <- subset(D, iseq)
-rls_fit(modelAR$prm, modelAR, Dfit, returnanalysis = FALSE)
-str(modelAR$Lfits[1:2])
-
-(i <- iseq[length(iseq)] + 1)
-Dnew <- subset(D, i:(i))
-Dnew$y
-
-modelAR$yAR
-D$y[(i-6):i]
-Dnew_transformed <- modelAR$transform_data(Dnew)
-
-Dnew_transformed$AR.lag0
-Dnew_transformed$AR.lag1
-Dnew_transformed$AR.lag4
-
-
-## modelAR$yAR
-## rls_update(modelAR, Dnew_transformed, Dnew[[model$output]])
-## yhatAR <- rls_predict(modelAR, Dnew_transformed)
-
-## modelAR$yAR
-
-## iseq <- i+modelAR$kseq
-## plot(D$t[iseq], D$y[iseq], type = "b", xlab = "t", ylab = "y")
-## lines(D$t[iseq], yhat[1,], type = "b", col = 2)
-## lines(D$t[iseq], yhatAR[1,], type = "b", col = 4)
-## legend("topright", c("observations",pst("predictions (",min(modelAR$kseq)," to ",max(modelAR$kseq)," steps ahead)"),
-##                      pst("Predictions AR (",min(model$kseq)," to ",max(model$kseq)," steps ahead)")), lty = 1, col = c(1,2,4))
-
-## (i <- iseq[length(iseq)] + 2)
-## Dnew <- subset(D, i:(i))
-## Dnew$y
-
-## modelAR$yAR
-## tail(Dfit$y)
-## Dnew_transformed <- modelAR$transform_data(Dnew)
-
-## Dnew_transformed$AR.lag3
-## Dnew_transformed$AR.Order2
-## Dnew_transformed$AR.Order1
-
-
-## modelAR$yAR
-## rls_update(modelAR, Dnew_transformed, Dnew[[model$output]])
-## yhatAR <- rls_predict(modelAR, Dnew_transformed)
-
-## modelAR$yAR
-
-## (i <- iseq[length(iseq)] + 2)
-## Dnew <- subset(D, i)
-## Dnew_transformed <- model$transform_data(Dnew)
-## rls_update(model, Dnew_transformed, Dnew[[model$output]])
-## yhat <- rls_predict(model, Dnew_transformed)
-
-## iseq <- i+modelAR$kseq
-## plot(D$t[iseq], D$y[iseq], type = "b", xlab = "t", ylab = "y")
-## lines(D$t[iseq], yhat[1,], type = "b", col = 2)
-## lines(D$t[iseq], yhatAR[1,], type = "b", col = 4)
-## legend("topright", c("observations",pst("predictions (",min(modelAR$kseq)," to ",max(modelAR$kseq)," steps ahead)"),
-##                      pst("Predictions AR (",min(model$kseq)," to ",max(model$kseq)," steps ahead)")), lty = 1, col = c(1,2,4))
diff --git a/misc-R/building-electricity-load-forecast.Rmd b/misc-R/building-electricity-load-forecast.Rmd
deleted file mode 100644
index fd42e149f438221314a2188b0a0e15ee8f3cd2b3..0000000000000000000000000000000000000000
--- a/misc-R/building-electricity-load-forecast.Rmd
+++ /dev/null
@@ -1,223 +0,0 @@
----
-title: "DRAFT Building electricity load forecasting"
-author: "Peder Bacher"
-date: "`r Sys.Date()`"
-output: rmarkdown::html_vignette
-vignette: >
-  %\VignetteIndexEntry{Building heat load forecasting}
-  %\VignetteEngine{knitr::rmarkdown}
-  %\VignetteEncoding{UTF-8}
-bibliography: literature.bib
----
-
-```{r setup, include = FALSE, purl = FALSE}
-knitr::opts_chunk$set(
-  collapse = TRUE,
-  comment = ">",
-  cache = FALSE, # Somehow it doesn't work
-  fig.height=6.5, 
-  fig.width=13, 
-  out.width="685px"
-)
-options(digits=3)
-```
-
-## Intro
-This vignette presents an example of using the onlineforecasting
-package for fitting a model and calculating forecasts as carried out
-by @Bacher2012.
-
-## Data
-Data for the forecasting examples is taken from a data set collected
-in Sønderborg, Denmark. It comprises heat load measurements for
-sixteen houses, together with local climate observations and weather
-forecasts (NWPs). The houses are generally built in the sixties and
-seventies, with a floor plan in the range of 85 to 170 $\mr{m^2}$, and
-constructed in bricks. For each house only the total heat load, including both space heating and hot tab water heating, is available. The climate observations are measured at the local district heating plant within 10 kilometers from the houses. The NWPs are from the HIRLAM-S05 model and provided by the Danish Meteorological Institute. All times are in UTC and the time stamp for average values are set to the end of the time interval.
-
-Load the package:
-```{r}
-##library(devtools)
-##load_all(as.package("../../onlineforecast"), export_all=FALSE)
-## Remember
-##install.packages("onlineforecast_0.1.0.tar.gz")
-library(onlineforecast)
-```
-
-The "Dbuilding" data is included in the package. Its a data.list (see the vignette onlineforecasting.pdf):
-```{r}
-Dbuilding <- readRDS("Dbuilding.Rda")
-```
-
-Keep it in D and see the content:
-```{r}
-D <- Dbuilding
-names(D)
-```
-
-The time:
-```{r}
-head(D$t)
-```
-
-The observed heat load (in kW) of the different houses are kept in a data.frame:
-```{r}
-head(D$Electricityload$house1)
-```
-
-The Numerical Weather Predictions (NWPs) of ambient temperature steps 0 to 8 hours ahead are:
-```{r}
-head(D$Ta[ ,1:9])
-```
-So at "2008-12-01 01:00:00 GMT" the latest available forecasts is the first row of Ta.
-
-Choose the heat load of House 9 for the example, just keep it as a vector y:
-```{r}
-D$y <- D$Electricityload$house3
-```
-
-A time series plot, see "?plot_ts.data.list" (Note how the forecasts Ta are lagged to be synced with observations (i.e. then also with each other))
-```{r}
-plot_ts(D, c("^y","Ta"), kseq=c(1,12))
-```
-
-A shorter period:
-```{r}
-plot_ts(D, c("^y","Ta"), "2010-12-15", "2010-12-25", kseq=c(1,12))
-
-plot_ts(D, c("^y","Ta"), "2009-12-15", "2011-12-25", kseq=c(1,12))
-
-
-tmp <- subset(D, c("2010-12-15", "2011-12-15"), kseq=1:6, pattern=c("^y|^Ta"))
-pairs(tmp)
-```
-
-Set the index of the training period and which period to evaluate (when fitting the points with scoreperiod==false are not included in the score evaluation)
-```{r}
-Dtrain <- subset(D, c("2010-06-01", "2011-06-01"))
-Dtrain$scoreperiod <- in_range("2010-06-10", Dtrain$t)
-```
-
-Define a model
-```{r}
-model <- forecastmodel$new()
-model$output = "y"
-model$add_inputs(Ta = "lp(Ta, a1=0.9)", 
-                 I = "lp(I, a1=0.7)", 
-                 mu_tday = "fs(tday/24, nharmonics=10)",
-                 mu = "ones()")
-model$add_regprm("rls_prm(lambda=0.9)")
-```
-
-Define the parameters to be optimized offline (their lower, init and upper bound)
-```{r}
-model$add_prmbounds(Ta__a1 = c(0.8, 0.9, 0.9999),
-                    I__a1 =  c(0.4, 0.8, 0.9999),
-                    lambda = c(0.9, 0.99, 0.9999))
-```
-
-Tune the parameters: first set the horizons to run and the 
-```{r, results="hide"}
-model$kseq <- c(1,18)
-model$prm <- rls_optim(model, Dtrain, control=list(maxit=2), cachedir = "cache-building-heat-load-forecasting")$par
-```
-
-Now fit with the optimized parameters on the entire period
-```{r}
-model$kseq <- 1:36
-val <- rls_fit(model$prm, model, D, returnanalysis = TRUE)
-```
-
-Plot the forecasts (Yhat adheres to the forecast matrix format and in plot_ts the forecasts are lagged k steps to sync with the observations)
-```{r, fig.height=4}
-D$Yhat <- val$Yhat
-plot_ts(D, c("^y|^Y"), "2010-06-01", "2011-06-01", kseq = c(1,18))
-plot_ts(D, c("^y|^Y"), "2010-07-01", "2010-07-15", kseq = c(1,18))
-plot_ts(D, c("^y|^Y"), "2010-08-01", "2010-08-05", kseq = c(1,18))
-plot_ts(D, c("^y|^Y"), "2010-12-15", "2010-12-30", kseq = c(1,18))
-```
-
-Plot a forecast for a particular time point
-```{r, fig.height=4}
-i <- 5000
-iseq <- i+model$kseq
-plot(D$t[iseq], D$y[iseq], type = "b", xlab = "t", ylab = "y")
-lines(D$t[iseq], D$Yhat[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)
-```
-
-### Recursive update and prediction
-First fit on a period
-```{r}
-iseq <- which(in_range("2010-12-15",D$t,"2011-01-01"))
-Dfit <- subset(D, iseq)
-rls_fit(model$prm, model, Dfit)
-```
-
-Now the fits are saved in the model object (its an R6 object, hence passed by reference to the functions and can be changed inside the functions). A list of fits with an entry for each horizon is in Lfits, see the two first
-```{r}
-str(model$Lfits[1:2])
-```
-
-Now new data arrives, take the point right after the fit period
-```{r}
-(i <- iseq[length(iseq)] + 1)
-Dnew <- subset(D, i)
-```
-
-First we need to transform the new data (This must only be done once for each new data, since some transform functions, e.g. lp(), actually keep states, see the detailed vignette onlineforecasting.pdf)
-```{r}
-Dnew_transformed <- model$transform_data(Dnew)
-```
-
-Then we can update the parameters using the transformed data
-```{r}
-rls_update(model, Dnew_transformed, Dnew[[model$output]])
-```
-
-Calculate predictions using the new data and the updated fits (rls coefficient estimates in model$Lfits[[k]]$theta)
-```{r}
-yhat <- rls_predict(model, Dnew_transformed)
-```
-
-Plot to see that it fits the observations
-```{r}
-iseq <- i+model$kseq
-plot(D$t[iseq], D$y[iseq], type = "b", xlab = "t", ylab = "y")
-lines(D$t[iseq], yhat, type = "b", col = 2)
-legend("topright", c("observations",pst("predictions (",min(model$kseq)," to ",max(model$kseq)," steps ahead)")), lty = 1, col = 1:2)
-```
-
-Run this for a longer period to verify that the same forecasts are obtained (in one go vs. iteratively)
-
-First in one go
-```{r}
-val <- rls_fit(model$prm, model, D, returnanalysis = TRUE)
-D$Yhat1 <- val$Yhat
-```
-
-and then iteratively
-```{r}
-itrain <- which(in_range("2010-12-15",D$t,"2011-01-01"))
-itest <- which(in_range("2011-01-01",D$t,"2011-01-04"))
-rls_fit(model$prm, model, subset(D, itrain))
-
-D$Yhat2 <- data.frame(matrix(NA, nrow(D$Yhat1), ncol(D$Yhat1)))
-names(D$Yhat2) <- names(D$Yhat1)
-for(i in itest){
-    print(i)
-    Dnew <- subset(D, i)
-    Dnewtr <- model$transform_data(Dnew)
-    rls_update(model, Dnewtr, Dnew[[model$output]])
-    D$Yhat2[i, ] <- as.numeric(rls_predict(model, Dnewtr))
-}
-```
-
-Compare to see the difference between the one step forecasts
-```{r}
-D$Yhat1$k1[itest] - D$Yhat2$k1[itest]
-plot(D$Yhat1$k1[itest], type="b")
-lines(D$Yhat2$k1[itest], type="b", col=2)
-```
-
-## Literature
diff --git a/misc-R/building-heat-load-forecasting-AR-and-error-models.R b/misc-R/building-heat-load-forecasting-AR-and-error-models.R
deleted file mode 100644
index 54c7b57a22ca9bd2a20658564525822b26ec657d..0000000000000000000000000000000000000000
--- a/misc-R/building-heat-load-forecasting-AR-and-error-models.R
+++ /dev/null
@@ -1,277 +0,0 @@
-## ----------------------------------------------------------------
-## Load the current version directly from the folder
-library(devtools)
-load_all(as.package("../../onlineforecast"))
-
-
-## ------------------------------------------------------------------------
-D <- Dbuildingheatload
-D$y <- D$heatload ## Here we are missing this one! so results might not be as interesting: D$Heatload$house15
-plot_ts(D, c("^y","Ta"), kseq=c(1,12))
-plot_ts(D, c("^y","Ta"), "2010-12-15", "2011-01-10", kseq=c(1,12))
-
-
-## ------------------------------------------------------------------------
-D$scoreperiod <- period("2010-12-20", D$t)
-itrain <- period(D$t, "2011-01-01")
-ieval <- period("2011-01-01", D$t)
-
-
-## ------------------------------------------------------------------------
-## Model with no AR
-model <- forecastmodel$new()
-model$output = "y"
-model$add_inputs(
-          Ta = "lp(Ta, a1=0.9)", 
-          I = "lp(I, a1=0.7)", 
-          mu_tday = "fs(tday/24, nharmonics=10)",
-          mu = "ones()")
-model$add_regp("rls_prm(lambda=0.9)")
-##
-model$add_prmbounds(Ta__a1 = c(0.8, 0.9, 0.9999),
-             I__a1 =  c(0.4, 0.8, 0.9999),
-             lambda = c(0.9, 0.99, 0.9999))
-##
-model$kseq <- c(1,18)
-model$prm <- rls_optim(model, subset(D,itrain), control=list(maxit=2))$par
-##
-model$kseq <- 1:36
-val <- rls_fit(model$prm, model, D, returnanalysis = TRUE)
-
-
-## ------------------------------------------------------------------------
-## Investigate the one step residuals, lag them to match the observations
-D$resid <- lag(val$Resid$k1, 1)
-
-plot(D$t[ieval], D$resid[ieval])
-
-acf(D$resid[ieval], na.action=na.pass)
-
-par(mfrow=c(1,2))
-plot(D$y[ieval], D$resid[ieval])
-plot(D$resid[ieval], val$Resid$k1[ieval])
-
-
-
-## ------------------------------------------------------------------------
-## Add an AR part
-model$add_inputs(AR = "AR(0)")
-##
-model$kseq <- c(1,18)
-model$prm <- rls_optim(model, subset(D,itrain), control=list(maxit=2))$par
-##
-model$kseq <- 1:36
-valAR <- rls_fit(model$prm, model, D, returnanalysis = TRUE)
-
-
-## ------------------------------------------------------------------------
-## An MA model for the k=1 ahead (AR on the residuals)
-modelMA <- forecastmodel$new()
-modelMA$output <- "resid"
-modelMA$add_inputs(AR = "AR(0)")
-modelMA$add_regp("rls_prm(lambda=0.9)")
-##
-modelMA$add_prmbounds(lambda = c(0.9, 0.99, 0.9999))
-##
-modelMA$kseq <- 1
-modelMA$prm <- rls_optim(modelMA, subset(D,itrain), control=list(maxit=2), cachedir = "")$par
-##
-valMA <- rls_fit(modelMA$prm, modelMA, D, returnanalysis = TRUE)
-
-
-## ------------------------------------------------------------------------
-## See the k step forecasts
-ieval <- period("2011-01-01", D$t)
-
-plot(D$t[ieval], lag(val$Resid$k1[ieval],1), type = "l", col = 1)
-lines(D$t[ieval], lag(valAR$Resid$k1[ieval],1), type = "l", col = 2)
-lines(D$t[ieval], lag(valMA$Resid$k1[ieval],1), type = "l", col = 3)
-
-par(mfrow=c(1,3))
-acf(val$Resid$k1[ieval], na.action=na.pass)
-acf(valAR$Resid$k1[ieval], na.action=na.pass)
-acf(valMA$Resid$k1[ieval], na.action=na.pass)
-
-rmse(val$Resid$k1[ieval])
-rmse(valAR$Resid$k1[ieval])
-rmse(valMA$Resid$k1[ieval])
-
-
-
-## ------------------------------------------------------------------------
-## Fit the MA (AR error) model for each horizon
-L <- lapply(1:36, function(k){
-    ## The k step residuals from the no AR model
-    D$resid <- lag(valAR$Resid[ ,k], k)
-    ##
-    modelMA$kseq <- k
-    ##
-    modelMA$prm <- rls_optim(modelMA, subset(D,itrain), control=list(maxit=2), cachedir = "")$par
-    ##
-    valMA <- rls_fit(modelMA$prm, modelMA, D, returnanalysis = TRUE)
-    valMA$Resid
-})
-valMAmulti <- list()
-valMAmulti$Resid <- do.call("cbind", L)
-
-
-
-## ------------------------------------------------------------------------
-## RMSE
-## Only the points which are not NA for all horizons, and have values for both models
-iAR <- apply(!is.na(valAR$Resid), 1, all)
-iMA <- apply(!is.na(valMAmulti$Resid), 1, all)
-i <- apply(!is.na(val$Resid), 1, all)
-i <- iAR & iMA & i
-## Only evaluation period
-i <- i & ieval
-##
-plot(val$Resid$k1[i], type="l")
-lines(valAR$Resid$k1[i], col=2)
-lines(valMAmulti$Resid$k1[i], col=3)
-##
-tmpAR <- apply(valAR$Resid[i, ], 2, rmse)
-tmpMA <- apply(valMAmulti$Resid[i, ], 2, rmse)
-tmp <- apply(val$Resid[i, ], 2, rmse)
-plot(tmpAR, type="b", ylim=range(tmp,tmpAR))
-points(tmpMA, type="b", col=2)
-points(tmp, type="b")
-
-acf(valMAmulti$Resid$k1[i], na.action=na.pass)
-pacf(valMAmulti$Resid$k1[i], na.action=na.pass)
-
-plot(lag(val$Resid$k1[i], 1), val$Resid$k1[i])
-plot(lag(valMAmulti$Resid$k1[i], 1), valMAmulti$Resid$k1[i])
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-## ------------------------------------------------------------------------
-## RMSE
-## Only the points which are not NA for all horizons, and have values for both models
-iAR <- apply(!is.na(valAR$Resid), 1, all)
-i <- apply(!is.na(val$Resid), 1, all)
-i <- iAR & i
-## Remove start
-i[1:200] <- FALSE
-##
-tmpAR <- apply(valAR$Resid[i, ], 2, rmse)
-tmp <- apply(val$Resid[i, ], 2, rmse)
-plot(model$kseq, tmpAR, type="b", ylim=range(tmp,tmpAR))
-points(model$kseq, tmp, type="b")
-
-acf(val$Resid$k1[i], na.action=na.pass)
-pacf(val$Resid$k1[i], na.action=na.pass)
-
-plot(lag(val$Resid$k1[i], 1), val$Resid$k1[i])
-plot(lag(val$Resid$k2[i], 2), val$Resid$k2)
-
-y <- val$Resid$k1[i]
-x <- lag(val$Resid$k1[i], 1)
-fit <- lm(y ~ 0 + x)
-plot(fit$residuals)
-rmse(fit$residuals)
-rmse(fit$residuals)
-
-## ##
-## setpar("ts", mfrow=c(ncol(val$Resid),1))
-## apply(val$Resid[i, ], 2, plot, type="l")
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-## Error model
-
-In the applied model there is no auto-regressive part, hence it is almost
-certain that there is auto-correlation in the residuals - which, for the shorter
-horizons, can be used to improve the forecasts.
-
-Check the auto-correlation for the one-step residuals
-
-
-acf(val$Resid$k1[i], na.action=na.pass)
-acf(val$Resid$k2[i], na.action=na.pass)
-acf(val$Resid$k3[i], na.action=na.pass)
-
-## Take the error from the training period
-
-D$residual.k1 <- lag(val$Resid$k1, 1)
-Dtrain$residual.k1 <- D$residual.k1[D$t %in% Dtrain$t]
-
-##
-plot_ts(Dtrain, "residual.k1")
-
-model <- forecastmodel$new()
-model$output <- "residual.k1"
-model$add_inputs(
-          AR = "lp(AR(c(0)), a1=0.2)")
-model$add_regp("rls_prm(lambda=0.9)")
-##
-model$add_prmbounds(AR__a1 =  c(0.4, 0.8, 0.9999),
-                    lambda = c(0.9, 0.99, 0.9999))
-##
-model$kseq <- c(1)
-model$prm <- rls_optim(model, Dtrain, control=list(maxit=2), cachedir = "")$par
-##
-model$kseq <- 1:36
-valMA <- rls_fit(model$prm, model, D, returnanalysis = TRUE)
-
-
-
-## ------------------------------------------------------------------------
-## RMSE
-## Only the points which are not NA for all horizons, and have values for both models
-iAR <- apply(!is.na(valAR$Resid), 1, all)
-iMA <- apply(!is.na(valMA$Resid), 1, all)
-i <- apply(!is.na(val$Resid), 1, all)
-i <- iAR & iMA & i
-## Remove start
-i[1:200] <- FALSE
-##
-plot(val$Resid$k1[i], type="l")
-lines(valAR$Resid$k1[i], col=2)
-lines(valMA$Resid$k1[i], col=3)
-##
-tmpAR <- apply(valAR$Resid[i, ], 2, rmse)
-tmpMA <- apply(valMA$Resid[i, ], 2, rmse)
-tmp <- apply(val$Resid[i, ], 2, rmse)
-plot(model$kseq, tmpAR, type="b", ylim=range(tmp,tmpAR))
-points(model$kseq, tmpMA, type="b")
-points(model$kseq, tmp, type="b")
-
-acf(valMA$Resid$k1[i], na.action=na.pass)
-pacf(valMA$Resid$k1[i], na.action=na.pass)
-
-plot(lag(val$Resid$k1[i], 1), val$Resid$k1[i])
-plot(lag(valMA$Resid$k1[i], 1), valMA$Resid$k1[i])
diff --git a/misc-R/building-heat-load-forecasting-test-AR.R b/misc-R/building-heat-load-forecasting-test-AR.R
deleted file mode 100644
index 79420967e49cc71b5eb302d44556e3dbacb064e1..0000000000000000000000000000000000000000
--- a/misc-R/building-heat-load-forecasting-test-AR.R
+++ /dev/null
@@ -1,310 +0,0 @@
-## ----------------------------------------------------------------
-## Load the current version directly from the folder
-library(devtools)
-load_all(as.package("../../onlineforecast"))
-
-
-## ------------------------------------------------------------------------
-class(Dbuildingheatload)
-
-
-## ------------------------------------------------------------------------
-D <- Dbuildingheatload
-names(D)
-
-
-## ------------------------------------------------------------------------
-head(D$t)
-
-
-## ------------------------------------------------------------------------
-head(D$Heatload[ ,1:9])
-
-
-## ------------------------------------------------------------------------
-head(D$Ta[ ,1:9])
-
-
-## ------------------------------------------------------------------------
-D$y <- D$Heatload$house15
-## ## ------------------------------------------------------------------------
- plot_ts(D, c("^y","Ta"), kseq=c(1,12))
-
-
-## ## ------------------------------------------------------------------------
- plot_ts(D, c("^y","Ta"), "2010-12-15", "2010-12-25", kseq=c(1,12))
-
-
-## ------------------------------------------------------------------------
-Dtrain <- subset(D, c("2010-12-15", "2011-01-01"))
-Dtrain$scoreperiod <- in_range("2010-12-20", Dtrain$t)
-
-
-## ------------------------------------------------------------------------
-model <- forecastmodel$new()
-model$output = "y"
-model$add_inputs(
-          Ta = "lp(Ta, a1=0.9)", 
-          I = "lp(I, a1=0.7)", 
-          mu_tday = "fs(tday/24, nharmonics=10)",
-          mu = "ones()",
-          AR = "lp(AR(c(0,1)), a1=0.2)")
-model$add_regprm("rls_prm(lambda=0.9)")
-
-## ## ------------------------------------------------------------------------
-## model <- forecastmodel$new()
-## model$output = "y"
-## model$add_inputs(Ta = "lp(Ta, a1=0.9)", 
-##                  I = "lp(I, a1=0.7)", 
-##                  mu_tday = "fs(tday/24, nharmonics=10)",
-##                  mu = "ones()",
-##                  AR = "AR(c(1,2))")
-## model$add_regp("rls_prm(lambda=0.9)")
-
-
-
-## ------------------------------------------------------------------------
-model$add_prmbounds(Ta__a1 = c(0.8, 0.9, 0.9999),
-                    I__a1 =  c(0.4, 0.8, 0.9999),
-                    AR__a1 =  c(0.4, 0.8, 0.9999),
-                    lambda = c(0.9, 0.99, 0.9999))
-
-
-model$kseq <- c(1,18)
-#datatr <- model$transform_data(Dtrain)
-#oi
-
-#val <- rls_fit(model$get_prmbounds("init"), model, D, returnanalysis = TRUE)
-
-##model2 <- forecastmodel$new()
-
-## ---- results="hide"-----------------------------------------------------
-model$kseq <- c(1,18)
-model$prm <- rls_optim(model, Dtrain, control=list(maxit=2), cachedir = "building-heat-load-forecasting_cache-rls")$par
-
-
-## ------------------------------------------------------------------------
-model$kseq <- 1:36
-val <- rls_fit(model$prm, model, D, returnanalysis = TRUE)
-
-model$reset_state()
-datatr <- model$transform_data(Dtrain)
-
-## ---- fig.height=4-------------------------------------------------------
-D$Yhat <- val$Yhat
-plot_ts(D, c("^y|^Y"), "2011-01-01", "2011-02-01", kseq = c(1,18))
-
-
-## ---- fig.height=4-------------------------------------------------------
-i <- 200
-iseq <- i+model$kseq
-plot(D$t[iseq], D$y[iseq], type = "b", xlab = "t", ylab = "y")
-lines(D$t[iseq], D$Yhat[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)
-
-
-## ------------------------------------------------------------------------
-iseq <- which(in_range("2010-12-15",D$t,"2011-01-01"))
-Dfit <- subset(D, iseq)
-rls_fit(model$prm, model, Dfit)
-
-
-## ------------------------------------------------------------------------
-str(model$Lfits[1:2])
-
-
-## ------------------------------------------------------------------------
-(i <- iseq[length(iseq)] + 1)
-##i <- i:(i+1)
-Dnew <- subset(D, i)
-
-
-## ------------------------------------------------------------------------
-Dnew_transformed <- model$transform_data(Dnew)
-
-
-## ------------------------------------------------------------------------
-rls_update(model, Dnew_transformed, Dnew[[model$output]])
-
-
-## ------------------------------------------------------------------------
-yhat <- rls_predict(model, Dnew_transformed)
-
-
-## ------------------------------------------------------------------------
-iseq <- i+model$kseq
-plot(D$t[iseq], D$y[iseq], type = "b", xlab = "t", ylab = "y")
-lines(D$t[iseq], yhat, type = "b", col = 2)
-legend("topright", c("observations",pst("predictions (",min(model$kseq)," to ",max(model$kseq)," steps ahead)")), lty = 1, col = 1:2)
-
-
-## ------------------------------------------------------------------------
-val <- rls_fit(model$prm, model, D, returnanalysis = TRUE)
-D$Yhat1 <- val$Yhat
-
-
-## ------------------------------------------------------------------------
-itrain <- which(in_range("2010-12-15",D$t,"2011-01-01"))
-itest <- which(in_range("2011-01-01",D$t,"2011-01-04"))
-rls_fit(model$prm, model, subset(D, itrain))
-
-D$Yhat2 <- data.frame(matrix(NA, nrow(D$Yhat1), ncol(D$Yhat1)))
-names(D$Yhat2) <- names(D$Yhat1)
-for(i in itest){
-    print(i)
-    Dnew <- subset(D, i)
-    Dnewtr <- model$transform_data(Dnew)
-    rls_update(model, Dnewtr, Dnew[[model$output]])
-    D$Yhat2[i, ] <- as.numeric(rls_predict(model, Dnewtr))
-}
-
-
-## ------------------------------------------------------------------------
-D$Yhat1$k1[itest] - D$Yhat2$k1[itest]
-
-
-valAR <- val
-
-
-## ------------------------------------------------------------------------
-## Compare with model with no AR
-model <- forecastmodel$new()
-model$output = "y"
-model$add_inputs(
-          Ta = "lp(Ta, a1=0.9)", 
-          I = "lp(I, a1=0.7)", 
-          mu_tday = "fs(tday/24, nharmonics=10)",
-          mu = "ones()")
-model$add_regp("rls_prm(lambda=0.9)")
-##
-model$add_prmbounds(Ta__a1 = c(0.8, 0.9, 0.9999),
-             I__a1 =  c(0.4, 0.8, 0.9999),
-             lambda = c(0.9, 0.99, 0.9999))
-##
-model$kseq <- c(1,18)
-model$prm <- rls_optim(model, Dtrain, control=list(maxit=2), cachedir = "building-heat-load-forecasting_cache-rls")$par
-##
-model$kseq <- 1:36
-val <- rls_fit(model$prm, model, D, returnanalysis = TRUE)
-
-
-## ------------------------------------------------------------------------
-## See the k step forecasts
-i <- 200
-iseq <- i:(i+7*24)
-plot(D$t[iseq], D$y[iseq], type = "b", xlab = "t", ylab = "y")
-k <- 1
-lines(D$t[iseq], lag(val$Yhat[iseq,k],k), type = "b", col = 2)
-lines(D$t[iseq], lag(valAR$Yhat[iseq,k],k), type = "b", col = 3)
-#legend("topright", c("Observations",pst("Predictions (",min(model$kseq)," to ",max(model$kseq)," steps ahead)")), lty = 1, col = 1:2)
-
-
-
-## ------------------------------------------------------------------------
-## RMSE
-## Only the points which are not NA for all horizons, and have values for both models
-iAR <- apply(!is.na(valAR$Resid), 1, all)
-i <- apply(!is.na(val$Resid), 1, all)
-i <- iAR & i
-## Remove start
-i[1:200] <- FALSE
-##
-tmpAR <- apply(valAR$Resid[i, ], 2, rmse)
-tmp <- apply(val$Resid[i, ], 2, rmse)
-plot(model$kseq, tmpAR, type="b", ylim=range(tmp,tmpAR))
-points(model$kseq, tmp, type="b")
-
-acf(val$Resid$k1[i], na.action=na.pass)
-pacf(val$Resid$k1[i], na.action=na.pass)
-
-plot(lag(val$Resid$k1[i], 1), val$Resid$k1[i])
-plot(lag(val$Resid$k2[i], 2), val$Resid$k2)
-
-y <- val$Resid$k1[i]
-x <- lag(val$Resid$k1[i], 1)
-fit <- lm(y ~ 0 + x)
-plot(fit$residuals)
-rmse(fit$residuals)
-rmse(fit$residuals)
-
-## ##
-## setpar("ts", mfrow=c(ncol(val$Resid),1))
-## apply(val$Resid[i, ], 2, plot, type="l")
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-## Error model
-
-In the applied model there is no auto-regressive part, hence it is almost
-certain that there is auto-correlation in the residuals - which, for the shorter
-horizons, can be used to improve the forecasts.
-
-Check the auto-correlation for the one-step residuals
-
-
-acf(val$Resid$k1[i], na.action=na.pass)
-acf(val$Resid$k2[i], na.action=na.pass)
-acf(val$Resid$k3[i], na.action=na.pass)
-
-## Take the error from the training period
-
-D$residual.k1 <- lag(val$Resid$k1, 1)
-Dtrain$residual.k1 <- D$residual.k1[D$t %in% Dtrain$t]
-
-##
-plot_ts(Dtrain, "residual.k1")
-
-model <- forecastmodel$new()
-model$output <- "residual.k1"
-model$add_inputs(
-          AR = "lp(AR(c(0)), a1=0.2)")
-model$add_regp("rls_prm(lambda=0.9)")
-##
-model$add_prmbounds()(AR__a1 =  c(0.4, 0.8, 0.9999),
-             lambda = c(0.9, 0.99, 0.9999))
-##
-model$kseq <- c(1)
-model$prm <- rls_optim(model, Dtrain, control=list(maxit=2))$par
-##
-model$kseq <- 1:36
-valMA <- rls_fit(model$prm, model, D, returnanalysis = TRUE)
-
-
-
-## ------------------------------------------------------------------------
-## RMSE
-## Only the points which are not NA for all horizons, and have values for both models
-iAR <- apply(!is.na(valAR$Resid), 1, all)
-iMA <- apply(!is.na(valMA$Resid), 1, all)
-i <- apply(!is.na(val$Resid), 1, all)
-i <- iAR & iMA & i
-## Remove start
-i[1:200] <- FALSE
-##
-plot(val$Resid$k1[i], type="l")
-lines(valAR$Resid$k1[i], col=2)
-lines(valMA$Resid$k1[i], col=3)
-##
-tmpAR <- apply(valAR$Resid[i, ], 2, rmse)
-tmpMA <- apply(valMA$Resid[i, ], 2, rmse)
-tmp <- apply(val$Resid[i, ], 2, rmse)
-plot(model$kseq, tmpAR, type="b", ylim=range(tmp,tmpAR))
-points(model$kseq, tmpMA, type="b")
-points(model$kseq, tmp, type="b")
-
-acf(valMA$Resid$k1[i], na.action=na.pass)
-pacf(valMA$Resid$k1[i], na.action=na.pass)
-
-plot(lag(val$Resid$k1[i], 1), val$Resid$k1[i])
-plot(lag(valMA$Resid$k1[i], 1), valMA$Resid$k1[i])
diff --git a/misc-R/data_soenderborg_heatsplit/heatSplit.sp_1.sn_4218597.rda b/misc-R/data_soenderborg_heatsplit/heatSplit.sp_1.sn_4218597.rda
deleted file mode 100644
index 9564ec094a40dcac87af30d58e86b7a73f069308..0000000000000000000000000000000000000000
Binary files a/misc-R/data_soenderborg_heatsplit/heatSplit.sp_1.sn_4218597.rda and /dev/null differ
diff --git a/misc-R/data_soenderborg_heatsplit/heatSplit.sp_1.sn_4218598.rda b/misc-R/data_soenderborg_heatsplit/heatSplit.sp_1.sn_4218598.rda
deleted file mode 100644
index a07a383a2c82367a7ebd96ede563c1937ff69a26..0000000000000000000000000000000000000000
Binary files a/misc-R/data_soenderborg_heatsplit/heatSplit.sp_1.sn_4218598.rda and /dev/null differ
diff --git a/misc-R/data_soenderborg_heatsplit/heatSplit.sp_1.sn_4711176.rda b/misc-R/data_soenderborg_heatsplit/heatSplit.sp_1.sn_4711176.rda
deleted file mode 100644
index 03058f1f0a5c1526e260aeefe7f583db4f4febfa..0000000000000000000000000000000000000000
Binary files a/misc-R/data_soenderborg_heatsplit/heatSplit.sp_1.sn_4711176.rda and /dev/null differ
diff --git a/misc-R/data_soenderborg_heatsplit/heatSplit.sp_1.sn_4724106.rda b/misc-R/data_soenderborg_heatsplit/heatSplit.sp_1.sn_4724106.rda
deleted file mode 100644
index df94d6800acc4ca2627f69066c56ec7041e7c1f9..0000000000000000000000000000000000000000
Binary files a/misc-R/data_soenderborg_heatsplit/heatSplit.sp_1.sn_4724106.rda and /dev/null differ
diff --git a/misc-R/data_soenderborg_heatsplit/heatSplit.sp_1.sn_4836681.rda b/misc-R/data_soenderborg_heatsplit/heatSplit.sp_1.sn_4836681.rda
deleted file mode 100644
index dd1718f1cfe956a7e6fc305f918097d7cc6140c5..0000000000000000000000000000000000000000
Binary files a/misc-R/data_soenderborg_heatsplit/heatSplit.sp_1.sn_4836681.rda and /dev/null differ
diff --git a/misc-R/data_soenderborg_heatsplit/heatSplit.sp_1.sn_4964553.rda b/misc-R/data_soenderborg_heatsplit/heatSplit.sp_1.sn_4964553.rda
deleted file mode 100644
index ada23acbb129bd9b074083b134b4f506357d49fa..0000000000000000000000000000000000000000
Binary files a/misc-R/data_soenderborg_heatsplit/heatSplit.sp_1.sn_4964553.rda and /dev/null differ
diff --git a/misc-R/data_soenderborg_heatsplit/heatSplit.sp_1.sn_5036505.rda b/misc-R/data_soenderborg_heatsplit/heatSplit.sp_1.sn_5036505.rda
deleted file mode 100644
index 5a3d750f93449ea78e821a22605bff05d9343a52..0000000000000000000000000000000000000000
Binary files a/misc-R/data_soenderborg_heatsplit/heatSplit.sp_1.sn_5036505.rda and /dev/null differ
diff --git a/misc-R/data_soenderborg_heatsplit/heatSplit.sp_1.sn_5107720.rda b/misc-R/data_soenderborg_heatsplit/heatSplit.sp_1.sn_5107720.rda
deleted file mode 100644
index 8d309aaa5574a7062dbfbca43ab8ba51120a3834..0000000000000000000000000000000000000000
Binary files a/misc-R/data_soenderborg_heatsplit/heatSplit.sp_1.sn_5107720.rda and /dev/null differ
diff --git a/misc-R/data_soenderborg_heatsplit/heatSplit.sp_1.sn_5159799.rda b/misc-R/data_soenderborg_heatsplit/heatSplit.sp_1.sn_5159799.rda
deleted file mode 100644
index 822d4e986c618223e7422aa9cfb7efc7b6564083..0000000000000000000000000000000000000000
Binary files a/misc-R/data_soenderborg_heatsplit/heatSplit.sp_1.sn_5159799.rda and /dev/null differ
diff --git a/misc-R/data_soenderborg_heatsplit/heatSplit.sp_1.sn_5164534.rda b/misc-R/data_soenderborg_heatsplit/heatSplit.sp_1.sn_5164534.rda
deleted file mode 100644
index ca7f47ad4d4b0bd8168de3eeb9d666a7a31eb09a..0000000000000000000000000000000000000000
Binary files a/misc-R/data_soenderborg_heatsplit/heatSplit.sp_1.sn_5164534.rda and /dev/null differ
diff --git a/misc-R/data_soenderborg_heatsplit/heatSplit.sp_1.sn_5183232.rda b/misc-R/data_soenderborg_heatsplit/heatSplit.sp_1.sn_5183232.rda
deleted file mode 100644
index 2cf1561fa66ac2218e197467f1208295bd62fbf2..0000000000000000000000000000000000000000
Binary files a/misc-R/data_soenderborg_heatsplit/heatSplit.sp_1.sn_5183232.rda and /dev/null differ
diff --git a/misc-R/data_soenderborg_heatsplit/heatSplit.sp_1.sn_5193768.rda b/misc-R/data_soenderborg_heatsplit/heatSplit.sp_1.sn_5193768.rda
deleted file mode 100644
index 837f6bd411c7c83a0e2ed55be09c71eb2d3a38ee..0000000000000000000000000000000000000000
Binary files a/misc-R/data_soenderborg_heatsplit/heatSplit.sp_1.sn_5193768.rda and /dev/null differ
diff --git a/misc-R/data_soenderborg_heatsplit/heatSplit.sp_1.sn_5194732.rda b/misc-R/data_soenderborg_heatsplit/heatSplit.sp_1.sn_5194732.rda
deleted file mode 100644
index be52c3c44a4dc5ca28f0c6aec10af0dbc9f45417..0000000000000000000000000000000000000000
Binary files a/misc-R/data_soenderborg_heatsplit/heatSplit.sp_1.sn_5194732.rda and /dev/null differ
diff --git a/misc-R/data_soenderborg_heatsplit/heatSplit.sp_1.sn_5194965.rda b/misc-R/data_soenderborg_heatsplit/heatSplit.sp_1.sn_5194965.rda
deleted file mode 100644
index e129d9a4ed0e015b7dc7fb7270bdf0f05312289e..0000000000000000000000000000000000000000
Binary files a/misc-R/data_soenderborg_heatsplit/heatSplit.sp_1.sn_5194965.rda and /dev/null differ
diff --git a/misc-R/data_soenderborg_heatsplit/heatSplit.sp_1.sn_5197381.rda b/misc-R/data_soenderborg_heatsplit/heatSplit.sp_1.sn_5197381.rda
deleted file mode 100644
index 6bfeccb293b7b2d3722bf6da8cdb4fe4cf82915a..0000000000000000000000000000000000000000
Binary files a/misc-R/data_soenderborg_heatsplit/heatSplit.sp_1.sn_5197381.rda and /dev/null differ
diff --git a/misc-R/data_soenderborg_heatsplit/heatSplit.sp_1.sn_5223036.rda b/misc-R/data_soenderborg_heatsplit/heatSplit.sp_1.sn_5223036.rda
deleted file mode 100644
index 9d2104bc2b978b594e85ddc46ccb491a45a58400..0000000000000000000000000000000000000000
Binary files a/misc-R/data_soenderborg_heatsplit/heatSplit.sp_1.sn_5223036.rda and /dev/null differ
diff --git a/misc-R/data_soenderborg_load.R b/misc-R/data_soenderborg_load.R
deleted file mode 100644
index b90dc01619d5adaacf3bbc40864eda6fc39ef8bb..0000000000000000000000000000000000000000
--- a/misc-R/data_soenderborg_load.R
+++ /dev/null
@@ -1,328 +0,0 @@
-# setting work directory and libraries #
-rm(list = ls())
-
-# # Packages used
-# require(R6)
-require(data.table)
-# require(Rcpp)
-# require(splines)
-library(devtools)
-library(roxygen2)
-
-pack <- as.package("../../onlineforecast")
-load_all(pack)
-
-
-# Importing data # First unzip to get the .csv system('unzip
-# ../data/DataSoenderborg.zip')
-data_or <- fread("data_soenderborg.csv", sep = ",", header = TRUE)
-data_or[, `:=`(t, asct(data_or$t))]
-setDF(data_or)
-names(data_or)[names(data_or) == "Ig.obs"] <- "I.obs"
-
-# Make a data.table for each variable
-tmp <- unlist(getse(strsplit(names(data_or)[-1], "\\."), 2))
-colnm <- unlist(getse(strsplit(names(data_or)[-1], "\\."), 1))
-nms <- unique(colnm[grep("^k\\d*$", tmp)])
-
-kmax <- 48
-data <- list()
-data$t <- data_or$t
-for (ii in 1:length(nms)) {
-    # Find the columns
-    i <- grep(pst("^", nms[ii], "$"), unlist(getse(strsplit(names(data_or)[-1],"\\."), 1))) + 1
-    # Take only with kxx
-    i <- i[grep("k[[:digit:]]+$", names(data_or)[i])]
-    # 
-    #
-    data[[nms[ii]]] <- lag(data_or[ ,i], -1:-length(i))
-    names(data[[nms[ii]]]) <- pst("k", 1:length(i))
-    row.names(data[[nms[ii]]]) <- NULL
-    data[[nms[ii]]] <- as.data.frame(data[[nms[ii]]])
-    # Check if observed climate variable is there
-    nm <- pst(nms[[ii]], ".obs")
-    if (nm %in% names(data_or)) {
-        data[[nm]] <- data_or[, nm]
-    }
-}
-# More
-cols <- pst("Heatload.house", 1:16)
-data[["Heatload"]] <- data_or[, cols]
-names(data[["Heatload"]]) <- pst("house", 1:16)
-# 
-data[["cosAoi"]] <- data_or[, "cosAoi.obs"]
-data[["sunElevation"]] <- data_or[, "sunElevation.obs"]
-
-
-# # The time of day
-# ncol <- ncol(data$Ta)
-# tmp <- aslt(data$t)$hour
-# tmp <- matrix(tmp, nrow = length(tmp), ncol = ncol)
-# tmp <- data.frame(t(t(tmp) + (0:(ncol - 1))))
-# names(tmp) <- pst("k", 0:(ncol - 1))
-# data$tday <- tmp%%24
-
-#
-class(data) <- c("data.list","list")
-
-# Save for other scripts to read it
-#saveRDS(data, "data_soenderborg.RDS")
-
-
-data$heatloadtotal <- sapply(1:nrow(data$Heatload), function(i){
-    mean(unlist(data$Heatload[i, ]), na.rm=TRUE)
-})
-data$totaln <- sapply(1:nrow(data$Heatload), function(i){
-    sum(is.na(unlist(data$Heatload[i, ])))
-})
-plot(data$t, data$totaln)
-
-# Write for building heat load forecasting
-#Dbuildingheatload <- subset(data, c("2010-12-15","2011-03-01"), nms=c("t","Heatload","Ta","I","Ws","Wd","Ta.obs","I.obs","Wd.obs","Ws.obs","cosAoi","sunElevation","tday"))
-data$heatload <- data$Heatload$house9
-Dbuildingheatload <- subset(data, c("2010-12-15","2011-03-01"), nms=c("t","heatload","heatloadtotal","Ta.obs","I.obs","Ta","I"))
-rownames(Dbuildingheatload$Ta) <- NULL
-Dbuildingheatload$Ta <- Dbuildingheatload$Ta[ ,1:36]
-rownames(Dbuildingheatload$I) <- NULL
-Dbuildingheatload$I <- Dbuildingheatload$I[ ,1:36]
-#
-usethis::use_data(Dbuildingheatload, overwrite=TRUE) 
-
-
-
-
-#----------------------------------------------------------------
-# The electricity was not in the above data, so add it
-# Open the info
-info <- read.table("data_soenderborg_selectedInfo.csv",header=TRUE,sep=",")
-
-# Open the hourly data, and make one data.frame with the power for each house
-load("data_soenderborg_orig.rda")
-#load("../data/heat_sp-1.rda")
-D <- D[D$sn %in% info$sn, ]
-
-# Make the acc. el energy (kWh) into into power (kW)
-L <- split(D,D$sn)
-x <- L[[3]]
-Lre <- lapply(L, function(x){
-  x$el <- c(NA,diff(x$P2) / as.numeric(diff(x$t),unit="hours"))
-  # Remove outliers
-  x$el[x$el>100] <- NA
-  x$el[x$el<0] <- NA
-# plot(x$t,x$el,ylim=c(0,20))
-#   plot(x$t,cumsumWithNA(x$el))
-# plot(x$t,x$P2)  
-  return(x)
-})
-D <- do.call("rbind",Lre)
-
-# Make into a data.frame with each the power for each seperate house as columns
-L <- split(D[ ,c("t","el")],D$sn)
-X <- L[[1]]
-houseid <- info$houseid[ which(info$sn==as.numeric(names(L)[1])) ]
-names(X)[2] <- pst("el",houseid)
-# Rename the columns
-for(i in 2:length(L))
-  {
-    tmp <- L[[i]]
-    houseid <- info$houseid[ which(info$sn==as.numeric(names(L)[i])) ]
-    names(tmp)[2] <- pst("el",houseid)
-    X <- merge(X,tmp,by="t",all=TRUE)
-  }
-
-# To hourly
-X <- resample(X, ts=3600, tstart=asct("2008-12-01"), tend=asct("2011-05-01"))
-
-#---------
-# SOME OF THE FOLLOWING MIGHT BE ERRORFUL (some leftover stuff from a merge was making problems)
-# Read the splitted heat
-load(pst("data_soenderborg_heatsplit/heatSplit.sp_1.sn_",info$sn[1],".rda"))
-tmp <- DH1
-names(tmp)[-1] <- pst("house",1,"_",names(tmp)[-1])
-Xor <- tmp
-#
-for(i in 2:nrow(info)){
-    load(pst("data_soenderborg_heatsplit/heatSplit.sp_1.sn_",info$sn[i],".rda"))
-    #
-    tmp <- DH1
-    names(tmp)[-1] <- pst("house",i,"_",names(tmp)[-1])
-    Xor <- merge(Xor,tmp,by="t")
-}
-
-
-# Put together data and X and Xor
-Dbuilding <- subset(data, nms=c("t","Heatload","Ta","I","Ws","Wd","Ta.obs","I.obs","Wd.obs","Ws.obs","cosAoi","sunElevation","tday"))
-range(Dbuilding$t)
-range(X$t)
-
-tmp <- X#[in_range("2010-12-15 01:00:00 GMT", X$t, "2011-02-01 00:00:00 GMT"), ]
-tmp <- tmp[ ,-1]
-names(tmp) <- pst("house", 1:16)
-Dbuilding$Electricityload <- tmp
-
-
-# Xor
-Dbuilding <- subset(Dbuilding, c("2009-01-01 00:00:00 GMT", "2011-05-01 00:00:00 GMT"))
-range(Dbuilding$t)
-range(Xor$t)
-tmp <- Xor[period("2008-12-01 01:00:00 GMT",Xor$t,"2011-05-01 00:00:00 GMT"), ]
-tmp <- tmp[ ,-1]
-Dbuilding$All <- tmp
-
-
-#
-str(Dbuilding)
-
-#
-plot(Dbuilding$Ta.obs, Dbuilding$Ta$k1)
-plot(Dbuilding$Ta.obs, lag(Dbuilding$Ta$k1, 0))
-plot(Dbuilding$Ta.obs, lag(Dbuilding$Ta$k1, 1))
-plot(Dbuilding$Ta.obs, lag(Dbuilding$Ta$k1, 2))
-
-plot(Dbuilding$All$house1_effekt)
-lines(Dbuilding$All$house1_Pheat)
-lines(Dbuilding$All$house1_Pwater)
-plot(Dbuilding$Electricityload$house1, Dbuilding$All$house1_P1)
-
-
-# Don't save it in this folder
-# saveRDS(Dbuilding, file="XX/Dbuilding.Rda")
-
-
-
-#----------------------------------------------------------------
-# Make for solar power forecasting example
-source("functions/aoi.R")
-data_all <- data
-# Take the observed global radiation
-names(data_all)
-data_all$y <- data_all$I.obs
-
-# Make a transformation into solar power, simply a projection on an inclined surface
-names(data_all)
-
-X <- data.frame(t=data_all$t, G=data_all$I.obs, sunElevation=data_all$sunElevation)
-X$IbNwp <- lag_vector(data_all$Ib$k1, 1)
-X$IdNwp <- lag_vector(data_all$Id$k1, 1)
-#
-# Some small positive morning values
-plot_ts(X[period("2011-04-01",X$t,"2011-04-10"), ], c("^G","sunElevation"))
-X$G[X$sunElevation < 0] <- 0
-plot_ts(X[period("2011-04-01",X$t,"2011-04-10"), ], c("^G","sunElevation"))
-#
-#----------------------------------------------------------------
-# First split into direct and diffuse
-zenith <- (pi/2) - X$sunElevation
-# Clearness index with fitted clear sky radiation
-G0 <- 1367 * (1 + 0.033*cos( 2*pi*as.POSIXlt(X$t)$yday/365)) * cos(zenith)
-kt <- X$G / (G0/1000)
-# air mass
-m <- 1/cos(zenith)     
-m[zenith>(0.95*pi/2)] <- 20
-# Only solar elevation above xx deg
-iNA <- 89<aoiToDeg(zenith) | aoiToDeg(zenith)<0
-iNA[is.na(iNA)] <- FALSE
-# i <- 10000:12800
-# plotTSBeg(2)
-# plot(X$t[i],kt[i],ylim=c(0,1.5))
-# plot(X$t[i],X$Ig[i],ylim=c(0,1000))
-# lines(X$t[i],cl$Surf$I[i],col=2)
-# plotTSXAxis(X$t[i])
-kt[iNA] <- NA
-kt[kt>1] <- 1
-# Split into Diffuse
-# Id <- Ig * (0.944-1.538*exp(-exp(2.808-4.759*kt+2.276*kt^2)))# + 0.125*m + 0.013*m^2)))
-# Id <- Ig * (0.952-1.041*exp(-exp(2.3-4.702*kt)))
-# oldktFunction <- function(kt,a=0.952,b=1.041,c=2.3,d=4.702){ 0.3*(a-b*exp(-exp(c-d*kt)))+0.1 }
-ktSigmoid <- function(kt,minOut,maxOut,offset,slope){ minOut+(maxOut-minOut)*(1-1/(1+exp(-(slope*(kt-offset))))) }
-# See how the sigmoid looks
-# t1 <- seq(0,1,by=0.01)
-# plot(t1,ktSigmoid(t1,minOut=0.12,maxOut=0.85,offset=0.45,slope=10),ylim=c(0,1))
-# Split it
-X$Gdiffuse <- X$G * ktSigmoid(kt,minOut=0.12,maxOut=0.85,offset=0.45,slope=10)
-X$Gdirect <- X$G - X$Gdiffuse
-# For solar evelation below xx deg, set to all diffuse
-X$Gdiffuse[iNA] <- X$G[iNA]
-X$Gdirect[iNA] <- 0
-#
-#plotmulti(X[period("2009-06-01",X$t,"2009-06-10"), ], c("^G","^Ib|^Id"))
-#
-#
-X$sinsunelev <- sin(X$sunElevation)
-# Clip it
-X$sinsunelev[X$sinsunelev<0.01] <- 0.01
-# Project the solar beam plane 
-X$Ib_solarplane <- X$Gdirect / X$sinsunelev
-
-#
-panel_power <- function(pAzimuth){
-    latitude <- 54.909275
-    longitude <- 9.807340
-    slope <- 45
-    X$cosAOI <- aoiCos1(X$t, latitude, longitude, slope, pAzimuth)
-    X$sinAOI <- sqrt(1 - X$cosAOI^2)
-    # Since the sun is behind the plane when cosAOI < 0
-    X$sinAOI[X$cosAOI < 0] <- 0
-    X$Ib_PV_plane <- X$Ib_solarplane * X$sinAOI
-    #
-    #plotmulti(X[period("2009-06-01",X$t,"2009-06-10"), ], c("^I|^G","cos|sin"))
-    #
-    # Apply an angle of incidence modifier on the direct
-    fKtab <- function(b0,cosTheta)
-    {
-        # Calculate the AOI modifier
-        Ktab <- 1 - b0*(1/cosTheta - 1)
-        # set the AOI modifier to 0 for cosTheta<0
-                                        #    Ktab[cosTheta<=0] <- 0
-                                        #    Ktab[Ktab<=0] <- 0
-        cosTheta0 <- 1/(1/b0+1)
-        x <- 1/(1+exp(-1000*(cosTheta-cosTheta0)))
-        Ktab <- x*Ktab
-        return(Ktab)
-    }
-    X$fKtab <- fKtab(0.1, X$cosAOI)
-    # The power is then a mix of direct and diffuse
-    X$P_PV_plane <- 0.6 * (X$fKtab * X$Ib_PV_plane + X$Gdiffuse) + 0.4 * X$G
-    #plotmulti(X[period("2009-07-01",X$t,"2009-07-10"), ], c("^I|^G|^P","sin","fKtab","^G$|^P_PV_plane$"))
-    #
-    #
-    X$y <- 6 * X$P_PV_plane
-    X$y[X$y > 5000] <- 2000
-    #
-    #
-    #    tmp <- as.datalist(data_all, period("2011-04-01", data_all$t, "2011-04-08"))
-    tmp <- X[period("2011-03-07", X$t, "2011-03-14"), ]
-    plot(tmp$t, tmp$G/max(tmp$G, na.rm=TRUE), type="l", ylab="Normalized")
-    title(paste("Panel surface towards",pAzimuthSeq_text[i]), line=-1.2)
-    lines(tmp$t, tmp$y/max(tmp$y, na.rm=TRUE), type="l", col=2)
-    X$P_PV_plane[X$P_PV_plane>1000] <- 998
-    return(X$P_PV_plane)
-}
-
-# Fit for these azimuth angles
-pAzimuthSeq <- c(-90, -20, 0, 90)
-pAzimuthSeq_text <- c("East", "South", "20 deg. East", "West")
-#
-
-setpar("ts", mfrow=c(length(pAzimuthSeq),1))
-for(i in 1:length(pAzimuthSeq)){
-    # The sine of the AOI on an inclined surface
-    panel_power(pAzimuthSeq[i])
-}
-legend("topright", c("Global radiation    ","Solar power"), lty=1, col=c(1,2))
-axis.POSIXct(1, X$t, xaxt="s")
-
-
-pAzimuthSeq <- c(-40, -20, 0, 20, 40)
-setpar("ts", mfrow=c(length(pAzimuthSeq),1))
-x <- lapply_cbind_df(pAzimuthSeq, panel_power)
-names(x) <- pst("azimuth.",gsub("-","m",pAzimuthSeq))
-
-data_all$PVpower <- x * 5
-plot_ts.data.list(data_all, c("^I$","PVpower"), kseq=1)
-
-
-# Write for solar power forecasting
-Dsolarpower <- subset(data_all, c("2010-01-01","2011-01-01"), nms=c("t","PVpower","I","Ta","I.obs","Ta.obs","cosAoi","sunElevation","tday"))
-#
-usethis::use_data(Dsolarpower, overwrite=TRUE) 
diff --git a/misc-R/data_soenderborg_make_fullset_with_electricity.R b/misc-R/data_soenderborg_make_fullset_with_electricity.R
deleted file mode 100644
index 4b841b5e5dc1f0c56652e4386752de30ee34297c..0000000000000000000000000000000000000000
--- a/misc-R/data_soenderborg_make_fullset_with_electricity.R
+++ /dev/null
@@ -1,110 +0,0 @@
-#### setting work directory and libraries ####
-rm(list = ls())
-
-## ## Packages used
-## require(R6)
-require(data.table)
-## require(Rcpp)
-## require(splines)
-library(devtools)
-library(roxygen2)
-
-pack <- as.package("../../onlineforecast")
-load_all(pack)
-
-
-##### Importing data #### First unzip to get the .csv system('unzip
-##### ../data/DataSoenderborg.zip')
-data_or <- fread("data_soenderborg.csv", sep = ",", header = TRUE)
-data_or[, `:=`(t, asct(data_or$t))]
-setDF(data_or)
-
-
-##----------------------------------------------------------------
-## The electricity was not in the above data, so add it
-## Open the info
-info <- read.table("data_soenderborg_selectedInfo.csv",header=TRUE,sep=",")
-
-## Open the hourly data, and make one data.frame with the power for each house
-load("data_soenderborg_orig.rda")
-#load("../data/heat_sp-1.rda")
-D <- D[D$sn %in% info$sn, ]
-
-## Make the acc. el energy (kWh) into into power (kW)
-L <- split(D,D$sn)
-x <- L[[3]]
-Lre <- lapply(L, function(x){
-  x$el <- c(NA,diff(x$P2) / as.numeric(diff(x$t),unit="hours"))
-  ## Remove outliers
-  x$el[x$el>100] <- NA
-  x$el[x$el<0] <- NA
-## plot(x$t,x$el,ylim=c(0,20))
-##   plot(x$t,cumsumWithNA(x$el))
-## plot(x$t,x$P2)  
-  return(x)
-})
-D <- do.call("rbind",Lre)
-
-## Make into a data.frame with each the power for each seperate house as columns
-L <- split(D[ ,c("t","el")],D$sn)
-X <- L[[1]]
-houseid <- info$houseid[ which(info$sn==as.numeric(names(L)[1])) ]
-names(X)[2] <- pst("el",houseid)
-## Rename the columns
-for(i in 2:length(L))
-  {
-    tmp <- L[[i]]
-    houseid <- info$houseid[ which(info$sn==as.numeric(names(L)[i])) ]
-    names(tmp)[2] <- pst("el",houseid)
-    X <- merge(X,tmp,by="t",all=TRUE)
-  }
-
-## To hourly
-X <- resample(X, ts=3600, tstart=asct("2008-12-01"), tend=asct("2011-05-01"))
-
-##---------
-## Read the splitted heat
-load(pst("data_soenderborg_heatsplit/heatSplit.sp_1.sn_",info$sn[1],".rda"))
-Xheat <- DH1[,c("t","Pheat")]
-names(Xheat) <- c("t",pst("heat",1))
-##
-for(i in 2:nrow(info))
-  {
-    load(pst("data_soenderborg_heatsplit/heatSplit.sp_1.sn_",info$sn[i],".rda"))
-    tmp <- DH1[,c("t","effekt")]
-    names(tmp) <- c("t",pst("heat",i))
-    Xheat <- merge(Xheat,tmp,by="t")
-  }
-
-## Put them together
-X <- merge(Xheat,X,by="t")
-
-## Check that they are the same houses:
-## tmp <- merge(X,data_or, by="t")
-## tmp$heat8 == tmp$Heatload.house8
-## plot(tmp$heat2 - tmp$Heatload.house2)
-## ##
-## i <- 1:1000
-## plot(tmp$heat8[i])
-## lines(tmp$Heatload.house8[i])
-
-range(data_or$t)
-names(data_or)
-range(X$t)
-names(X)
-
-## The electricity to data_or
-tmp <- X[ ,c(1,grep("^el",names(X)))]
-names(tmp) <- gsub("el", "Elecload.house", names(tmp))
-data <- merge(data_or,tmp, by="t")
-
-unique(diff(data$t))
-range(data$t)
-
-tmp <- as.data.list(data)
-
-
-
-## Write it
-tmp <- data[in_range("2010-01-01", data$t, "2011-01-01"), ]
-write.table(tmp, "~/tmp/data_soenderborg_2010.csv", sep=",", row.names=FALSE)
diff --git a/misc-R/data_soenderborg_orig.rda b/misc-R/data_soenderborg_orig.rda
deleted file mode 100644
index dea9776305edff74d62062762ffcc9fa6ab5dfcb..0000000000000000000000000000000000000000
Binary files a/misc-R/data_soenderborg_orig.rda and /dev/null differ
diff --git a/misc-R/data_soenderborg_selectedInfo.csv b/misc-R/data_soenderborg_selectedInfo.csv
deleted file mode 100644
index 711b4f5ac29e6cac8ffb4503b286955c7a3f4bda..0000000000000000000000000000000000000000
--- a/misc-R/data_soenderborg_selectedInfo.csv
+++ /dev/null
@@ -1,17 +0,0 @@
-"sn","pdfside","adresse","BBRm2","opfoert","Type","Ydervaeg","Tag","bem","beboere","pejs","bortrejst","elgulvvarme","solvarme","supplVarme","natsaenkning","pdfsidenr","houseid"
-4218597,15,"Skovbrynet 23",151,1970,"Fritliggende enfamilieshus (parcelhus)","Mursten (tegl, kalksandsten, cementsten)","Tegl","",2,0,0,0,0,0,0,15,1
-4218598,36,"Frejasvej 10",163,1969,"Fritliggende enfamilieshus (parcelhus)","Mursten (tegl, kalksandsten, cementsten)","Tegl","",2,0,0,0,0,0,0,36,2
-4711176,40,"Laurids Skaus vej 2",140,1963,"Fritliggende enfamilieshus (parcelhus)","Mursten (tegl, kalksandsten, cementsten)","Fibercement, herunder asbest (bølge- eller skifereternit)","",2,0,0,0,0,0,0,40,3
-4724106,21,"Violvej 1",86,1952,"Fritliggende enfamilieshus (parcelhus)","Mursten (tegl, kalksandsten, cementsten)","Tagpap (med taghældning)","",2,0,0,0,0,0,0,21,4
-4836681,6,"Henrik Ibsens Vej 5",111,1966,"Fritliggende enfamilieshus (parcelhus)","Letbeton (lette bloksten, gasbeton)","Metalplader (bølgeblik, aluminium og lignende)","",1,0,0,0,0,0,0,6,5
-4964553,8,"Dybbøløstenvej 4",119,1963,"","","","",2,0,0,0,0,0,0,8,6
-5036505,10,"Møllegade 47",119,1947,"","","","",3,0,0,0,0,0,0,10,7
-5107720,2,"Agtoftsvej 39",160,1965,"Fritliggende enfamilieshus (parcelhus)","Mursten (tegl, kalksandsten, cementsten)","Fibercement, herunder asbest (bølge- eller skifereternit)","olietank?",4,0,0,0,0,0,0,2,8
-5159799,17,"Peter Graus Vej 5",173,1965,"Fritliggende enfamilieshus (parcelhus)","Mursten (tegl, kalksandsten, cementsten)","Fibercement, herunder asbest (bølge- eller skifereternit)","",1,0,0,0,0,0,0,17,9
-5164534,41,"Brunhoved 2",135,1996,"","","","",2,0,0,0,0,0,1,41,10
-5183232,33,"Aprilvej 4",122,1966,"","","","",2,0,0,0,0,0,0,33,11
-5193768,16,"Brombærhegnet 18",136,1975,"","","","",2,0,0,0,0,0,0,16,12
-5194732,NA,"Parkgade 44",86,1937,"","","","",2,0,1,0,0,0,1,29,13
-5194965,27,"Gammel Aabenraavej 29",123,1965,"","","","",2,0,0,0,0,0,1,27,14
-5197381,26,"Midtkobbel 23",127,1953,"","","","",2,0,0,0,0,0,0,26,15
-5223036,7,"Gyvelvej 3",137,1967,"","","","",5,0,0,0,0,0,1,7,16
diff --git a/misc-R/functions/aoi.R b/misc-R/functions/aoi.R
deleted file mode 100644
index dd5cacbab94ffd2041afe0d7d0591c90d952a334..0000000000000000000000000000000000000000
--- a/misc-R/functions/aoi.R
+++ /dev/null
@@ -1,173 +0,0 @@
-aoiToRad <- function(angle)
-  {
-    angle/180 * pi
-  }                                   
-
-aoiToDeg <- function(rad)
-  {
-    rad/pi * 180
-  }
-
-aoiEquationOfTime <- function(time)
-{
-  day <- as.POSIXlt(time)$yday
-  # Day is the number of days since the start of the year
-  b <- aoiToRad((360/365) * (day-81))
-  # Equation of time in minutes
-  9.87 * sin(2*b) - 7.53 * cos(b) - 1.5 * sin(b)
-}
-
-aoiLocalSolarTime <- function(time, longitude)
-{
-  ## time is in UTC
-  ## longitude is in degrees. Positive east of greenwich. The earth rotates 1 degrees in 4 minutes.
-  time + ( 4 * aoiToDeg(longitude) + aoiEquationOfTime(time) ) * 60
-}
-
-aoiSunHourAngle <- function(time, longitude)
-{
-  ## localSolarTime is in seconds since 1970 given as an POSIXct object
-  t <- as.POSIXlt(aoiLocalSolarTime(time, longitude))
-  LST.tod <- t$hour + t$min/60 + t$sec/3600
-  ## Output hourAngle in rad
-  aoiToRad( 15 * (LST.tod  - 12) )
-}
-
-aoiSunDeclination <- function(time)
-  {
-    #### Calculate declination angle
-    day <- as.POSIXlt(time)$yday + 1
-    ## Return the result in rad
-    asin( sin(aoiToRad(23.45)) * sin(aoiToRad((360/365)*(day-81))) )
-  }
-
-aoiSunElevation <- function(latitude, hourAngle, declination)
-  {
-    asin( cos(hourAngle) * cos(declination) * cos(latitude)  + sin(declination) * sin(latitude) )
-  }
-
-aoiSunElevationDeg <- function(time, latitude, longitude)
-  {
-    ## All input angles are given in degrees, transform them into rad
-    latitude <- aoiToRad(latitude)
-    longitude <- aoiToRad(longitude)
-    ## Calculate the earth declination
-    declination <- aoiSunDeclination(time)
-    ## Calculate the hourAngle
-    hourAngle <- aoiSunHourAngle(time, longitude)
-    ## Calculate the elevation angle of the sun
-    aoiSunElevation(latitude, hourAngle, declination)
-  }
-
-aoiSunAzimuth <- function(latitude, elevation, declination, hourAngle)
-  {
-    ## Works only for latitudes above the max declination of the earth: 23.45 degrees
-    if(abs(aoiToDeg(latitude))<=23.45){ stop("Works only for latitudes above the max declination of the earth: 23.45 degrees") }
-    ##
-    sAzimuth <- acos( (sin(declination) * cos(latitude) - cos(declination) * sin(latitude) * cos(hourAngle))/( cos(elevation) ) )
-    i <- hourAngle > 0
-    hourAngle[i] <- { sAzimuth[i] <- 2*pi - sAzimuth[i] }
-    sAzimuth
-  }
-
-aoiSunAzimuthDeg <- function(time, latitude, longitude)
-  {
-    if(latitude<=23.45){ stop("Works only for latitudes above the max declination of the earth: 23.45 degrees") }
-    ## All input angles are given in degrees, transform them into rad
-    latitude <- aoiToRad(latitude)
-    longitude <- aoiToRad(longitude)
-    ## Calculate the earth declination
-    declination <- aoiSunDeclination(time)
-    ## Calculate the hourAngle
-    hourAngle <- aoiSunHourAngle(time, longitude)
-    ## Calculate the elevation angle of the sun
-    elevation <- aoiSunElevation(latitude, hourAngle, declination)
-    ## Return in Rad
-    aoiSunAzimuth(latitude, elevation, declination, hourAngle)
-  }
-
-aoiCos1 <- function(time, latitude, longitude, slope, pAzimuth)
-{
-  ## All angle are given in degrees.
-  ## slope: is the angle between the normal to the ground surface, and the normal of the panel.
-  ## pAzimuth: azimuth of the panel where 0 degrees is due south. + is toward west, - is toward east.
-  
-  ## All input angles are given in degrees, transform them into rad
-  longitude <- aoiToRad(longitude)
-  latitude <- aoiToRad(latitude)
-  slope <- aoiToRad(slope)
-  pAzimuth <- aoiToRad(pAzimuth)
-
-  ## Calculate the earth declination
-  dcl <- aoiSunDeclination(time)
-
-  ## Calculate the hourAngle
-  hourAngle <- aoiSunHourAngle(time, longitude)
-
-  ## Calculate the angle of incidence
-  cosAOI <- (sin(dcl) * sin(latitude) * cos(slope) - sin(dcl) * cos(latitude) * sin(slope) * cos(pAzimuth)
-               + cos(dcl) * cos(latitude) * cos(slope) * cos(hourAngle)
-               + cos(dcl) * sin(latitude) * sin(slope) * cos(pAzimuth) * cos(hourAngle)
-               + cos(dcl) * sin(slope) * sin(pAzimuth) * sin(hourAngle))
-  
-  ## Return the result 
-  return(cosAOI)
-}
-
-aoiCos2 <- function(time, latitude, longitude, slope, pAzimuth)
-  {
-    ## All angle are given in degrees.
-    ## slope: is the angle between the normal to the ground surface, and the normal of the panel.
-    ## pAzimuth: azimuth of the panel where 0 degrees is due south. + is toward west, - is toward east.
-    pAzimuth <- pAzimuth + 180 # The algorithm uses a panel azimuth where 180 degree is due south
-    
-    ## All input angles are given in degrees, transform them into rad
-    longitude <- aoiToRad(longitude)
-    latitude <- aoiToRad(latitude)
-    slope <- aoiToRad(slope)
-    pAzimuth <- aoiToRad(pAzimuth)
-
-    ## Calculate the earth declination
-    declination <- aoiSunDeclination(time)
-
-    ## Calculate the hourAngle
-    hourAngle <- aoiSunHourAngle(time, longitude)
-
-    ## Calculate the zenith angle of the sun
-    elevation <- aoiSunElevation(latitude, hourAngle, declination)
-    zenith <- (pi/2) - elevation
-
-    ## Calculate the azimuth of the sun
-    sAzimuth <- aoiSunAzimuth(latitude, elevation, declination, hourAngle)
-
-    ## Calculate the angle of incidence
-    cosAOI <-  cos(slope) * cos(zenith) + sin(slope) * sin(zenith) * cos(sAzimuth - pAzimuth)
-
-    ## Return the result
-    return(cosAOI)
-  }
-
-
-
-## ###########################################################################
-## ## Tests and examples
-## ## All times must be in GMT
-## t <- seq(ISOdate(2009,6,1,0),ISOdate(2009,6,3,0),by=60)
-
-## ## This will calculate the elevation at the solar collectors at Byg in radians
-## latitude <- 25
-## longitude <- 12
-## ## latitude <- 37
-## ## longitude <- -2
-
-## ## Solar elevation
-## elevRad <- aoiSunElevation( aoiToRad(latitude), aoiSunHourAngle(t,aoiToRad(longitude)), aoiSunDeclination(t))
-## elevRad2 <- aoiSunElevationDeg(t,latitude,longitude)
-## ##
-## plot(aoiSunAzimuthDeg(t,latitude,longitude),aoiToDeg(elevRad))
-
-## lines(t,aoiToDeg(elevRad2),col=2)
-
-## ## Angle of incidence of some surface
-## par(mfrow=c(2,1))
-## plot(t,aoiCos1(t, longitude, latitude, slope=90, pAzimuth=0))
diff --git a/misc-R/lm-example.R b/misc-R/lm-example.R
deleted file mode 100644
index cfab7c4e80071094e6d6f85abcdac5f023a83322..0000000000000000000000000000000000000000
--- a/misc-R/lm-example.R
+++ /dev/null
@@ -1,179 +0,0 @@
-rm(list = ls())
-library(devtools)
-load_all(as.package("../../onlineforecast"))
-
-D <- Dbuildingheatload
-names(D)
-D$y <- D$heatload
-
-
-plot_ts(D, c("^y","Ta"), kseq=c(1,12))
-
-
-
-
-Dtrain <- subset(D, c("2010-12-15", "2011-01-01"))
-Dtrain$scoreperiod <- in_range("2010-12-20", Dtrain$t)
-
-model <- forecastmodel$new()
-model$output = "y"
-model$add_inputs(Ta = "lp(Ta, a1=0.9)",
-                 AR = "AR(lags=c(0))",     # The ambient temperature
-                 mu = "ones()") # ones() generates a matrix of ones (i.e. an intercept is included)
-
-model$add_prmbounds(Ta__a1 = c(0.8, 0.9, 0.9999))
-
-model$kseq <- c(1,18)
-model$p <- lm_optim(model, Dtrain, control=list(maxit=2))$par
-
-model$kseq <- 1:36
-val <- lm_fit(model$p, model, D, returnanalysis = TRUE)
-
-undebug(lm_fit)
-
-names(val)
-val$Lfitval$k1
-val$scoreperiod
-head(val[[6]])
-
-head(val$Yhat, 20)
-
-names(model)
-D$Yhat <- val$Yhat
-dev.new()
-plot_ts(D, c("^y|^Y"), kseq = c(1,18))
-
-
-
-####
-
-model$insert_p(model$p)
-datatr <- model$transform_data(D)
-names(datatr)
-
-X <- as.data.frame(subset(datatr, kseq = 1, lagforecasts = TRUE))
-inputnms <- names(X)
-## Add the model output to the data.frame for lm()
-X[ ,model$output] <- D[[model$output]]
-
-head(X,10)
-## Generate the formula
-frml <- pst(model$output, " ~ ", pst(inputnms, collapse=" + "), " - 1")
-## Fit the model
-fit <- lm(frml, X)
-
-summary(fit)
-head(fitted(fit))
-
-
-head(predict(fit, X), 10)
-head(predict(fit, Xpred))
-
-
-head(X,10)
-head(Xpred,10)
-#k <- model$kseq[1]
-#fit <- model$Lfits[[1]]
-## Form the regressor matrix, don't lag
-Xpred <- as.data.frame(subset(datatr, kseq = 1))
-pred <- predict(fit, Xpred)
-pred2 <- predict(fit, X)
-head(pred)
-head(pred2,10)
-summary(fit)
-head(Xpred)
-
-X$Ta.k1[2]*fit$coef[1] + X$mu.k1[2]*fit$coef[2]
-
-
-head(fit$residuals)
-head(X$y - pred)
-
-###
-
-head(val$Yhat)
-
-head(pred, 10)
-head(pred2,10)
-
-## Predictions
-x <- rnorm(15)
-y <- x + rnorm(15)
-fitTest <- lm(y ~ x)
-predict(fitTest)
-fitted(fitTest)
-
-new <- data.frame(x = seq(-3, 3, 0.5))
-predict(fitTest, new, se.fit = TRUE)
-
-
-
-# Take data (See vignette "building-heat-load-forecasting" for better model and more details)
-D <- subset(Dbuildingheatload, c("2010-12-15", "2011-01-01"))
-D$y <- D$Heatload[ ,1]
-# Define a model 
-model <- forecastmodel$new()
-model$output <- "y"
-model$add_inputs(Ta = "Ta")
-model$add_regp("rls_prm(lambda=0.9)")
-
-# 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)
-# And the sequence of horizons to fit for
-model$kseq <- 1:6
- 
-# Now we can fit the model and get the model validation analysis data
-L <- rls_fit(p = c(lambda=0.99), model = model, data = D, returnanalysis = TRUE)
-names(L)
-plot(L$Yhat$k1) # The one-step forecast
-plot(L$Lfitval$k1) # The one-step RLS coefficient over time
-# Fitting with lower lambda makes the RLS parameter change faster
-L <- rls_fit(p = c(lambda=0.9), model = model, data = D, returnanalysis = TRUE)
-names(L)
-plot(L$Lfitval$k1) # The one-step RLS coefficient over time
-# It can return a score
-rls_fit(c(lambda=0.99), model, D, scorefun=rmse, returnanalysis = FALSE)
-# Such that it can be passed to an optimzer (see ?rls_optim for a nice wrapper of optim)
-val <- optim(c(lambda=0.99), rls_fit, model = model, data = D, scorefun = rmse, returnanalysis = FALSE, method = "L-BFGS-B", lower = 0.5, upper = 0.9999)
-val$p
-
-# See rmse as a function of horizon
-val <- rls_fit(p = c(lambda=0.9), model = model, data = D, returnanalysis = TRUE, scorefun = rmse)
-names(val)
-head(val$scoreval, 100)
-plot(val$scoreval)
-
-
-
-
-# Take data (See vignette "building-heat-load-forecasting" for better model and more details)
-D <- subset(Dbuildingheatload, c("2010-12-15", "2011-01-01"))
-D$y <- D$Heatload[ ,1]
-# Define a model 
-model <- forecastmodel$new()
-model$output <- "y"
-model$add_inputs(Ta = "lp(Ta, a1=0.9)")
-model$add_prmbounds(Ta__a1 = c(0.8, 0.9, 0.9999))
-# 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)
-# And the sequence of horizons to fit for
-model$kseq <- 1:6
- 
-# Now we can fit the model and get the model validation analysis data
-L <- lm_fit(p = c(Ta__a1 = 0.7), model = model, data = D, returnanalysis = TRUE)
-names(L)
-plot(L$Yhat$k1) # The one-step forecast
-
-# The coefficients for each model
-head(L$Lfitval)
-# It can return a score
-lm_fit(c(Ta__a1=0.7), model, D, scorefun=rmse, returnanalysis = FALSE)
-# Such that it can be passed to an optimzer (see ?rls_optim for a nice wrapper of optim)
-val <- optim(c(Ta__a1=0.7), lm_fit, model = model, data = D, scorefun = rmse, returnanalysis = FALSE, method = "L-BFGS-B", lower = 0.5, upper = 0.9999)
-val$p
-
-# See rmse as a function of horizon
-val <- lm_fit(p = c(Ta__a1 = 0.7), model = model, data = D, returnanalysis = TRUE, scorefun = rmse)
-names(val)
-head(val$scoreval, 100)
-plot(val$scoreval)
diff --git a/misc-R/loadforecast_example_short.R b/misc-R/loadforecast_example_short.R
deleted file mode 100644
index f34cdd1c659d53cb747378cd98698c3f37b98121..0000000000000000000000000000000000000000
--- a/misc-R/loadforecast_example_short.R
+++ /dev/null
@@ -1,269 +0,0 @@
-## setting work directory and libraries ####
-rm(list = ls())
-
-## Packages used
-library(devtools)
-library(roxygen2)
-
-pack <- as.package("../../onlineforecast")
-load_all(pack)
-
-## Read the data
-load("../data/Dbuildingheatload.rda")
-Dall <- Dbuildingheatload
-
-## Set the model output y
-Dall$y <- Dall$heatload
-
-kseq <- 0:49
-Dall$AR0 <- lapply_cbind_df(kseq, function(k){
-    Dall$y
-})
-nams(Dall$AR0) <- pst("k",kseq)
-
-Dall$AR0
-
-
-
-## Plot some time series
-##plot_ts(Dall, patterns=c("Ta","Ws","Wd","Id|I$|Ib","^y$"), kseq=c(1,6,12,18,24,36), tstart="2011-01-01", tend="2011-01-10")
-
-## A pairs plot
-##pairs(subset(Dall, pattern="Ta", kseq=c(1:3,6,12,18)), cex=0.5)
-##pairs(subset(Dall, pattern="Ta", kseq=c(1:3,6,12,18), lagforecasts=TRUE), cex=0.5)
-
-##############################################################
-## Recursive Least Squares fitting the model has inputs (a list of inputobjects)
-## Define the model
-##
-Model <- model_class$new(output = "y",
-                         inputs = list(Ta      = "Ta",#"lp(Ta, a1=0.9)",
-                                       I       = "I",#"lp(I, a1=0.7)",
-                                       AR0     = "AR0"#"lp(I, a1=0.7)",
-##                                     mu_tday = "fs(tday/24, nharmonics=10)",
-                                        #mu      = "ones()"),
-                                       ),
-                         MA     = c(0),
-                         fitprm = "rls_prm(lambda=0.9)")
-
-## Set the parameters which are optimized in the offlines setting
-## Model$prmopt_upper <- c(Ta__a1 = 1,   I__a1 = 1,    lambda = 1) - 1e-04
-## Model$prmopt_init <-  c(Ta__a1 = 0.9, I__a1 = 0.81, lambda = 0.99)
-## Model$prmopt_lower <- c(Ta__a1 = 0.5, I__a1 = 0.5,  lambda = 0.9)
-Model$prmopt_upper <- c(lambda = 1) - 1e-04
-Model$prmopt_init <-  c(lambda = 0.99)
-Model$prmopt_lower <- c(lambda = 0.8)
-
-
-##############################################################
-## Recursive Least Squares fitting the model has inputs (a list of inputobjects)
-## The horizons to fit for
-Model$kseq <- c(1)
-prmopt <- Model$prmopt_init
-
-## The fit function always initialize a new fit, no matter what is in 'Model'
-Dpast <- subset(Dall, in_range("2010-12-01", Dall$t, "2011-02-01"))
-## Define a logical series which sets the fitting period
-Dpast$scoreperiod <- in_range("2010-12-15", Dpast$t, "2011-02-01")
-## From there make a test set
-Dtest <- subset(Dall, in_range("2011-02-01", Dall$t, "2011-02-03"))
-
-## Fit to see it is working
-rls_fit(prmopt = Model$prmopt_init,
-        model = Model,
-        data = Dpast,
-        scorefun = rmse)
-
-
-## This can be passed to an optimizer to optimize the offline parameters (keep them)
-Model$prmopt <- optim(par = Model$prmopt_init,
-                      fn = rls_fit,
-                      model = Model,
-                      data = Dpast,
-                      scorefun = rmse,
-                      lower = Model$prmopt_lower,
-                      upper = Model$prmopt_upper)$par
-
-
-## The model could be copied (useful for saving a particular instance)
-## Model1 <- Model$clone(deep = TRUE)
-
-## Use the optimized parameters and fit for all horizons
-#Model$kseq <- 1:36
-#rls_fit(Model$prmopt, Model, Dpast)
-
-## Calculate predictions on a new data (NOTE: This uses the latest RLS parameter values, it doesn't update recursively)
-#prd <- rls_predict(Model, Model$transform_data(Dtest))
-#prd
-
-
-## ################################################################
-## Plot k step ahead prediction series
-Model$kseq <- c(1)
-Yhat <- rls_fit(Model$prmopt, Model, Dpast, return_analysis=TRUE)$Yhat
-Dpast$yhat <- Yhat$k1
-
-tmp <- subset(Dpast, in_range("2011-01-01", Dpast$t, "2011-01-05"))
-plot(tmp$t, tmp$y, type = "l")
-lines(tmp$t, tmp$yhat, type = "l", col = 2)
-
-
-## ################################################################
-## Plot a prediction for a particular time point
-prd <- rls_predict(Model, Model$transform_data(Dtest))
-##
-plot(1:length(Dtest$y), Dtest$y)
-##
-lines(1:length(prd[1, ]), prd[1, ])
-for (i in 2:nrow(prd)) {
-  lines(i:(length(prd[1, ])+i-1), prd[i, ], col=i)
-}
-
-
-## ################################################################
-## Recursive update example
-
-## First fit on the past data (resets input and parameter states)
-rls_fit(Model$prmopt, Model, Dpast)
-
-## Recursive updating parameters and prediction
-prd <- as.data.frame(matrix(NA, nrow = length(Dtest$t), ncol = length(Model$kseq)))
-names(prd) <- paste0("k",Model$kseq)
-##
-for (it in 1:length(Dtest$t)) {
-  print(paste(it,"of",length(Dtest$t)))
-  ## A new datalist
-  Dnew <- subset(Dtest, it)
-  ## Generate the inputs
-  ## Important: Note this must only be done once for new input data, since it continues from last state when lowpass filtering
-  D <- Model$transform_data(Dnew)
-  ## New output observations
-  y <- Dnew[[Model$output]]
-  ## We can now update the parameters
-  rls_update(Model, D, y)
-  ## and make a prediction
-  prd[it, ] <- unlist(rls_predict(Model, D))
-}
-
-## Lag the predictions to match observations
-prd <- lag(prd, lags = "+k")
-plot(Dtest$y)
-lines(prd$k8)
-
-## Check if we get the same as when we do it in one fit
-Dboth <- subset(Dall, in_range(min(Dpast$t)-3600, Dall$t, max(Dtest$t)))
-Valanalysis <- rls_fit(Model$prmopt, model = Model, data = Dboth, return_analysis = TRUE)
-
-n <- length(Dboth$t)
-prd2 <- Valanalysis$Yhat[(n-nrow(prd)+1):n, ]
-prd2 - prd
-max(prd2 - prd, na.rm = TRUE)
-
-
-## ################################################################
-## But the error is auto-correlated
-tmp <- rls_fit(Model$prmopt, Model, Dpast, return_analysis=TRUE)
-## For 1-step ahead
-Dpast$Yhat <- tmp$Yhat
-Dpast$yhat <- tmp$Yhat$k1
-residuals <- Dpast$y - Dpast$yhat
-##
-plot(residuals)
-acf(residuals, na.action = na.pass)
-
-## For 2-step ahead
-Dpast$yhat <- tmp$Yhat$k2
-residuals <- Dpast$y - Dpast$yhat
-##
-plot(residuals)
-acf(residuals, na.action = na.pass)
-
-plot_ts(Dpast$Yhat[950:1000, ], "k1$|k2$")
-
-## Keep the residuals to use as input to an error model
-## Keep the residuals to use as output
-Dpast$r <- Dpast$y - tmp$Yhat$k1
-## (Now use only the one-step ahead forecast as input (and output))
-##Dpast$R <- Dpast$y - tmp$Yhat
-Dpast$R <- matrix(Dpast$r, nrow = length(residuals), ncol = 36)
-nams(Dpast$R) <- pst("k",1:36)
-
-## See them for different horizonts (they are equal now)
-plot_ts(Dpast$R[900:1000, ], "k[[:digit:]]$")
-
-## Define a model
-Me <- model_class$new(output = "r",
-                      inputs = list(R = "R"),
-                      fitprm = "rls_prm(lambda=0.9)")
-
-## Set the parameters which are optimized in the offlines setting
-Me$prmopt_upper <- c(lambda = 1) - 1e-04
-Me$prmopt_init <-  c(lambda = 0.95)
-Me$prmopt_lower <- c(lambda = 0.9)
-
-## Set the horizons for offline optimization
-Me$kseq <- c(1,2)
-prmopt <- Me$prmopt_init
-
-## Fit to check it is working
-Val <- rls_fit(prmopt = Me$prmopt_init,
-        model = Me,
-        data = Dpast,
-        scorefun = rmse,
-        return_analysis = TRUE)
-
-tmp <- cbind(r = Dpast$r, Val$Yhat)
-plot_ts(tmp[which(Dpast$scoreperiod)[1000:1100], ], "*")
-
-
-Dpast$ResidualsYhat <- Val$Yhat
-
-plot_ts(subset(Dpast,Dpast$scoreperiod), "Residuals", 1:2)
-plot_ts(subset(Dpast,which(Dpast$scoreperiod)[1:200]), "Residuals", 1:2)
-
-## This can be passed to an optimizer to optimize the offline parameters (keep them)
-Me$prmopt <- optim(par = Me$prmopt_init,
-                   fn = rls_fit,
-                   model = Me,
-                   data = Dpast,
-                   scorefun = rmse,
-                   lower = Me$prmopt_lower,
-                   upper = Me$prmopt_upper)$par
-
-
-## Predict on both train and test period
-## First make residuals on the entire period
-Dboth <- subset(Dall, in_range(min(Dpast$t)-3600, Dall$t, max(Dtest$t)))
-Yhat <- rls_fit(Model$prmopt, model = Model, data = Dboth, return_analysis = TRUE)$Yhat
-
-Dboth$r <- Dboth$y - Yhat$k1
-Dboth$R <- matrix(Dboth$r, nrow = length(Dboth$r), ncol = 36)
-nams(Dboth$R) <- pst("k",1:36)
-
-## Use the optimized parameters and fit a model for all horizons
-Me$kseq <- 1:36
-
-Rhat <- rls_fit(Me$prmopt, Me, Dboth, return_analysis = TRUE)$Yhat
-
-scoreperiod <- in_range("2010-12-15", Dboth$t, "2011-02-01")
-#itest <- in_range("2011-02-01", Dboth$t, "2011-02-03")
-
-Yhat_r <- Dboth$y[scoreperiod] - Yhat[scoreperiod, ]
-(rmse_Yhat_r <- apply(Yhat_r, 2, rmse))
-
-## The final forecasts
-Yhat_combined <- Yhat + Rhat
-
-Rhat_r <- Dboth$y[scoreperiod] - Yhat_combined[scoreperiod, ]
-(rmse_Rhat_r <- apply(Rhat_r, 2, rmse))
-#apply(Dboth$r[scoreperiod] - tmp2$Yhat[scoreperiod, ], 2, rmse)
-
-plot(rmse_Yhat_r, ylim = range(rmse_Yhat_r,rmse_Rhat_r))
-points(rmse_Rhat_r, col=2)
-
-
-k <- 1
-plot_ts(cbind(y=Dboth$y, yhat=Yhat[ ,1], yhat_combined=Yhat_combined[ ,1])[1100:1400, ], "*")
-
-k <- 2
-plot_ts(cbind(y=Dboth$y, yhat=Yhat[ ,k], yhat_combined=Yhat_combined[ ,k])[1100:1400, ], "*")
diff --git a/misc-R/make_example_julian.R b/misc-R/make_example_julian.R
deleted file mode 100644
index 48b81b8f50b3dee572180182bd2e446a4f1794d4..0000000000000000000000000000000000000000
--- a/misc-R/make_example_julian.R
+++ /dev/null
@@ -1,11 +0,0 @@
-## Go first and compile the package
-
-##
-library(knitr)
-purl("building-electricity-load-forecast.Rmd")
-
-## Then add files
-files <- c("../../onlineforecast_0.1.0.tar.gz",
-           "./Dbuilding.Rda",
-           "./building-electricity-load-forecast.R")
-zip("~/tmp/julianExample.zip", files, flags="-r9X-j")
diff --git a/misc-R/plotly-test.R b/misc-R/plotly-test.R
index c441871c8136bc64c2e1b802402502bdb45ecdbd..1bd3dabf5ff4f7c4b800b8578c6885e42ebf329d 100644
--- a/misc-R/plotly-test.R
+++ b/misc-R/plotly-test.R
@@ -4,7 +4,7 @@ library(devtools)
 load_all(as.package("../../onlineforecast"))
 
 
-D <- Dbuildingheatload
+D <- Dbuilding
 
 
 
diff --git a/misc-R/temp.R b/misc-R/temp.R
index 2605cc7ebfdd8070107c273b09d48b2998e13d84..e13dd6977747875f1b7117aca1353083f77f9d06 100644
--- a/misc-R/temp.R
+++ b/misc-R/temp.R
@@ -4,7 +4,7 @@ load_all(as.package("../../onlineforecast"))
 ?lm_predict
 
 # Take data
-D <- subset(Dbuildingheatload, c("2010-12-15", "2011-01-01"))
+D <- subset(Dbuilding, c("2010-12-15", "2011-01-01"))
 D$y <- D$heatload
 # Define a model 
 model <- forecastmodel$new()
diff --git a/tests/testthat/test-rls-heat-load.R b/tests/testthat/test-rls-heat-load.R
index aa6b5ccf537258234feb36898424415d740c0d98..fdf9835fb762733e4d47168425cd1efe772d68da 100644
--- a/tests/testthat/test-rls-heat-load.R
+++ b/tests/testthat/test-rls-heat-load.R
@@ -8,7 +8,7 @@ context("running RLS test")
 test_that("run", {
 
     ## ------------------------------------------------------------------------
-    D <- Dbuildingheatload
+    D <- Dbuilding
     D$y <- D$heatload
     D$tday <- make_tday(D$t, kseq=1:36)
 
diff --git a/vignettes/forecast-evaluation.Rmd b/vignettes/forecast-evaluation.Rmd
index 9a9ea8c429b9a974a60d0106aa1abd2187d68e1e..da823e9d4ca86174d44b8b44281a2ef76f4d2c90 100644
--- a/vignettes/forecast-evaluation.Rmd
+++ b/vignettes/forecast-evaluation.Rmd
@@ -15,19 +15,21 @@ bibliography: literature.bib
 library(knitr)
 # This vignettes name
 vignettename <- "forecast-evaluation"
-# Read external code from init.R
-knitr::read_chunk("init.R")
+# Read external code
+knitr::read_chunk("shared-init.R")
 ```
 ```{r init, cache=FALSE, include=FALSE, purl=FALSE}
 ```
+```{r links, cache=FALSE, echo=FALSE, purl=FALSE, results="asis"}
+```
 
 
 ## Intro
 This vignette provides a short overview of the basics of forecast evaluation
 with the functions from the onlineforecast package. It follows up on the
-vignettes ??{ref} and ??{ref} and continues the building load forecast modelling
-presented there. If something is introduced in the present text, but not
-explained, then have a look in the two preceding vignettes to find an explanation.
+vignettes [setup-data](setup-data.html) and [setup-and-use-model](setup-and-use-model.html), and continues the building load forecast
+modelling presented there. If something is introduced in the present text, but
+not explained, then have a look in the two preceding vignettes to find an explanation.
 
 ## Load forecasts
 
@@ -42,7 +44,7 @@ load_all(as.package("../../onlineforecast"))
 Just start by:
 ```{r}
 # Keep the data in D to simplify notation
-D <- Dbuildingheatload
+D <- Dbuilding
 # Keep the model output in y (just easier code later)
 D$y <- D$heatload
 # 
@@ -175,7 +177,7 @@ model. This is however not at all trivial, since the suitable reference model
 depends on the particular case of forecasting, e.g. the suitable reference model
 for wind power forecasting is not the same as for solar power forecasting -
 even within the same application the suitable reference model can be different
-depending on particular conditions etc. ??(referencer)
+depending on particular conditions etc. <!-- ??(referencer) -->
 
 In general the fundamental reference model should be the simplest reasonable
 model not relying on any inputs, hence either a model based on a mean
@@ -381,7 +383,7 @@ The RLS coefficients are in the fit for each horizon:
 str(fit1$Lfitval$k1)
 ```
 
-A pairs plot with residuals and inputs ??comment if patterns are left:
+A pairs plot with residuals and inputs <!-- ??comment --> if patterns are left:
 ```{r plotpairs, fig.height=figwidth}
 kseq <- c(1,36)
 D$Residuals <- residuals(fit2)[ ,pst("h",kseq)]
diff --git a/vignettes/make.R b/vignettes/make.R
index fe04dc7053ed19ab111d8f37f0f47776a9b22b03..97eb472332c785faab9dc5fda61ec2cbdc40771e 100644
--- a/vignettes/make.R
+++ b/vignettes/make.R
@@ -5,9 +5,9 @@ library(rmarkdown)
 dirnam <- "tmp-output/"
 dir.create(dirnam)
 
-makeit <- function(nam, openit=FALSE){
+makeit <- function(nam, openit=FALSE, clean=TRUE){
     namrmd <- paste0(nam,".Rmd")
-    render(namrmd, output_file=paste0(dirnam,nam))
+    render(namrmd, output_file=paste0(dirnam,nam), clean=clean)
     purl(namrmd)
     system(paste0("mv ",nam,".R ",dirnam,nam,".R"))
     if(openit){ system(paste0("chromium-browser ",dirnam,nam,".html &")) }
@@ -20,7 +20,9 @@ makeit("setup-data", openit=FALSE)
 file.remove(dir("cache", full.names=TRUE))
 file.remove("cache")
 file.remove(dir("tmp-output/tmp-setup-and-use-model/", full.names=TRUE))
-makeit("setup-and-use-model", openit=FALSE)
+
+makeit("setup-and-use-model", openit=FALSE, clean=FALSE)
+oi
 
 #
 file.remove(dir("tmp-output/tmp-forecast-evaluation/", full.names=TRUE))
diff --git a/vignettes/online-updating.Rmd b/vignettes/online-updating.Rmd
index 4e08d3d348379d3358f6e19593fd6293e6052196..36f77563ff91b28f7205b7d7654fbc3612341656 100644
--- a/vignettes/online-updating.Rmd
+++ b/vignettes/online-updating.Rmd
@@ -18,11 +18,13 @@ bibliography: literature.bib
 library(knitr)
 # This vignettes name
 vignettename <- "online-updating"
-# Read external code from init.R
-knitr::read_chunk("init.R")
+# Read external code
+knitr::read_chunk("shared-init.R")
 ```
 ```{r init, cache=FALSE, include=FALSE, purl=FALSE}
 ```
+```{r links, cache=FALSE, echo=FALSE, purl=FALSE, results="asis"}
+```
 
 
 ## Intro
@@ -39,7 +41,7 @@ load_all(as.package("../../onlineforecast"))
 Load data, setup and define a model:
 ```{r, output.lines=10}
 # Keep the data in D to simplify notation
-D <- Dbuildingheatload
+D <- Dbuilding
 # Set the score period 
 D$scoreperiod <- in_range("2010-12-20", D$t)
 # Set the training period
@@ -82,7 +84,9 @@ Now new data arrives, take the point right after the fit period
 Dnew <- subset(D, i)
 ```
 
-First we need to transform the new data (This must only be done once for each new data, since some transform functions, e.g. lp(), actually keep states, see the detailed vignette in ??)
+First we need to transform the new data (This must only be done once for each
+new data, since some transform functions, e.g. lp(), actually keep states, see
+the detailed description in )
 ```{r}
 Dnew_transformed <- model$transform_data(Dnew)
 ```
diff --git a/vignettes/setup-and-use-model.Rmd b/vignettes/setup-and-use-model.Rmd
index 037a94760accb3116ad4de291410b6a3fe50e0ab..67ed731aa07dfe1d6c99d44bd28f4f4656ad13d3 100644
--- a/vignettes/setup-and-use-model.Rmd
+++ b/vignettes/setup-and-use-model.Rmd
@@ -18,37 +18,36 @@ bibliography: literature.bib
 library(knitr)
 # This vignettes name
 vignettename <- "setup-and-use-model"
-# Read external code from init.R
-knitr::read_chunk("init.R")
+# Read external code
+knitr::read_chunk("shared-init.R")
 ```
 ```{r init, cache=FALSE, include=FALSE, purl=FALSE}
 ```
+```{r links, cache=FALSE, echo=FALSE, purl=FALSE, results="asis"}
+```
 
 
 ## Intro
 This vignette explains how to setup and use an onlineforecast
-model. This takes offset in the example of building heat load forecasting??(ref)
-and assumes that the data is setup correctly, as explained in
-[setup-data vignette](setup-data.html).
+model. This takes offset in the example of [building heat load
+forecasting] and assumes that the data is setup correctly, as explained in
+[setup-data](setup-data.html) vignette. The R code is available
+[here](setup-and-use-model.R). More information on [onlineforecasting.org].
 
-Load the package:
+Start by loading the package:
 ```{r}
-## Load the package
+# Load the package
 library(onlineforecast)
-```
-
-Just start by:
-```{r}
-# Keep the data in D to simplify notation
-D <- Dbuildingheatload
+# Set the data in D to simplify notation
+D <- Dbuilding
 ```
 
 
 ## Score period
 
 Set the `scoreperiod` as a logical vector with same length as `t`. It controls
-which points are included in score calculations in the package functions when the parameters are being optimized. this
-must be set in the `data.list` used.
+which points are included in score calculations in functions for optimization
+etc. It must be set.
 
 Use it to exclude a burn-in period of one week:
 ```{r}
@@ -114,23 +113,22 @@ model$add_inputs(Ta = "Ta",
                  mu = "ones()")
 ```
 So this is really where the structure of the model is specified. The inputs are
-given a name (`Ta` and `mu`), which each are set as an R expression (in a
+given a name (`Ta` and `mu`), which each are set as an R expression (given as a
 string). The expressions defines the **transformation step**: they will each
 be evaluated in an environment with a given `data.list`. This means that the
-variables from the data (e.g. `D`) can be used in the expressions - below in [Input transformations] we will detail this evaluation.
+variables from the data can be used in the expressions (e.g. `Ta` is in `D`) - below in [Input transformations] we will detail this evaluation.
 
 Next step for setting up the model is to set the parameters for the **regression
-step** by providing an expression, as a string, which returns the regression
+step** by providing an expression, which returns the regression
 parameter values. In the present case we will use the Recursive Least Squares
-(RLS) when regressing and we need to set the forgetting factor `lambda` by:
+(RLS) when regressing, and we need to set the forgetting factor `lambda` by:
 ```{r}
 # Regression step parameters
 model$add_regprm("rls_prm(lambda=0.9)")
 ```
 
 The expression is just of a function, which returns
-a list - in this case with the value of `lambda` ??(see onlineforecast
-vignette). The result of it begin evaluated is kept in:
+a list - in this case with the value of `lambda` (see [onlineforecasting]). The result of it begin evaluated is kept in:
 ```{r}
 # The evaluation happens with
 eval(parse(text="rls_prm(lambda=0.9)"))
@@ -154,8 +152,6 @@ The horizons to fit for is actually not directly related to the model, but
 rather the fitting of the model. In principle, it would be more "clean" if the
 model, data and fit was kept separate, however for recursive fitting this
 becomes un-feasible.
-<!-- see more ??(ref, where to we emphasize the recursive
-fitting, maybe a vignette by it self!?) -->
 
 
 ### Tune the parameters
@@ -201,7 +197,7 @@ Plot the forecasts (`Yhat` adheres to the forecast matrix format and in
 # 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))
+plot_ts(D, c("^heatload$|^Y"), kseq = c(1,6,18,36))
 ```
 We clearly see the burn-in period, where the forecasts vary a lot, 
 
@@ -308,7 +304,7 @@ There are quite a few functions available for input transformations:
 - `bspline()` wraps the `bs()` function for generating base splines.
 - `AR()` generates auto-regressive model inputs.
 
-and they can even be combined, see more details in ??(ref) and in their help
+and they can even be combined, see more details in [onlineforecasting] and in their help
 description, e.g. `?fs`.
 
 
@@ -331,8 +327,6 @@ D$Yhat2 <- fit2$Yhat
 # Plot all
 plot_ts(D, c("^heatload$|^Y"), kseq = c(1,18))
 ```
-See more on how to extend this model even further in ??(ref til
-buildingloadforecast vignette på hjemmeside)
 
 We can see the summary:
 ```{r}
@@ -363,10 +357,9 @@ legend("topleft", c("Input: Ta","Input: Low-pass Ta"), lty=1, col=1:2)
 ```
 We can see, that we obtained improvements. Around 3-4% for the longer horizons.
 
-For more on evaluation, see the vignette ??(ref til forecast-evaluation.html)
-
-For more development of the load forecast model see ??(building load forecast).
+For more on evaluation, see the vignette [forecast-evaluation](forecast-evaluation.html).
 
+See more on how to extend this model even further in [building heat load forecasting].
 
 
 ## Time of day and using observations as input
diff --git a/vignettes/setup-data.Rmd b/vignettes/setup-data.Rmd
index 7f758b547ddea260409babbd6d12751bd5ee7531..8efe8f6f95e865a1db11b8df2f010ec3f296f3fc 100644
--- a/vignettes/setup-data.Rmd
+++ b/vignettes/setup-data.Rmd
@@ -20,18 +20,20 @@ library(knitr)
 ## This vignettes name
 vignettename <- "setup-data"
 ## Read external code
-knitr::read_chunk("init.R")
+knitr::read_chunk("shared-init.R")
 ```
 ```{r init, cache=FALSE, include=FALSE, purl=FALSE}
 ```
+```{r links, cache=FALSE, echo=FALSE, purl=FALSE, results="asis"}
+```
+
 
 
 ## Intro
-This vignette explains how to setup data consisting of
-observations and forecasts, such that it can be used for onlineforecast
-models. A generic introduction and description is in available in the paper
-?(ref) [onlineforecasting vignette](onlineforecasting.pdf). ??(more on the other
-vignettes and website and code from this vignette)
+This vignette explains how to setup data consisting of observations and
+forecasts, such that it can be used for onlineforecast models. A generic
+introduction and description is in available in [onlineforecasting]. The code is
+available [here](setup-data.R). More information on [onlineforecasting.org].
 
 
 ## Data example
@@ -43,14 +45,14 @@ library(onlineforecast)
 ```
 
 In the package different data sets are included. The
-`Dbuildingheatload` holds the data used for the example of
+`Dbuilding` holds the data used for the example of
 heat load forecasting in the building-heat-load-forecasting vignette.
 
 When the package is loaded the data is also loaded, so we can access it
 directly. Let's start out by:
 ```{r}
 ## Keep it in D to simplify notation
-D <- Dbuildingheatload
+D <- Dbuilding
 ```
 
 The class is 'data.ĺist':
@@ -308,9 +310,8 @@ lines(D$t[i], y)
 legend("topright", c("8-step forecasts lagged","Observations"), lty=1, col=2:1)
 ```
 
-Of course that model was very simple, see how to make a better model in the
-??(examples and website)
-[solar forecast vignette](solar-power-forecasting.html).
+Of course that model was very simple, see how to make a better model in
+[building-heat-load-forecasting] and more information on the [website].
 
 
 
diff --git a/vignettes/init.R b/vignettes/shared-init.R
similarity index 84%
rename from vignettes/init.R
rename to vignettes/shared-init.R
index 0195156219f50ae67abe84c599ad0bd8f2a9fe0a..d72bc9b16fa41738826c67b7d66c7d2d3de4191e 100644
--- a/vignettes/init.R
+++ b/vignettes/shared-init.R
@@ -56,3 +56,9 @@ knit_hooks$set(output = function(x, options) {
   x <- paste(c(x, ""), collapse = "\n")
   hook_output(x, options)
 })
+
+
+# ---- links
+cat("[onlineforecasting]: https://onlineforecasting.org/articles/onlineforecasting.pdf\n")
+cat("[building heat load forecasting]: https://onlineforecasting.org/examples/building-heat-load-forecasting.html\n")
+cat("[onlineforecasting.org]: https://onlineforecasting.org\n")