rls_optim.R

# Do this in a separate file to see the generated help:
#library(devtools)
#document()
#load_all(as.package("../../onlineforecast"))
#?rls_optim


#' Optimize parameters (transformation stage) of RLS model
#'
#' This is a wrapper for \code{\link{optim}} to enable easy use of bounds and caching in the optimization.
#'
#' One smart trick, is to cache the optimization results. Caching can be done by providing a path to the
#' \code{cachedir} argument (relative to the current working directory).
#' E.g. \code{rls_optim(model, D, cachedir="cache")} will write a file in the folder 'cache', such that
#' next time the same call is carried out, then the file is read instead of running the optimization again.
#' See the example in \url{https://onlineforecasting.org/vignettes/nice-tricks.html}.
#' 
#' @title Optimize parameters for onlineforecast model fitted with RLS
#' @param model The onlineforecast model, including inputs, output, kseq, p
#' @param data The data.list which holds the data on which the model is fitted.
#' @param kseq The horizons to fit for (if not set, then model$kseq is used)
#' @param scorefun The function to be score used for calculating the score to be optimized.
#' @param cachedir A character specifying the path (and prefix) of the cache file name. If set to \code{""}, then no cache will be loaded or written. See \url{https://onlineforecasting.org/vignettes/nice-tricks.html} for examples.
#' @param cachererun A logical controlling whether to run the optimization even if the cache exists. 
#' @param printout A logical determining if the score function is printed out in each iteration of the optimization.
#' @param method The method argument for \code{\link{optim}}.
#' @param ... Additional parameters to \code{\link{optim}}
#' @return Result object of optim().
#' Parameters resulting from the optimization can be found from \code{result$par}
#' @seealso \code{link{optim}} for how to control the optimization.
#' @examples
#'
#' # Take data
#' D <- subset(Dbuilding, c("2010-12-15", "2011-01-01"))
#' D$y <- D$heatload
#' # Define a simple model 
#' model <- forecastmodel$new()
#' model$add_inputs(Ta = "Ta", mu = "one()")
#' model$add_regprm("rls_prm(lambda=0.99)")
#'
#' # Before fitting the model, define which points to include in the evaluation of the score function
#' D$scoreperiod <- in_range("2010-12-20", D$t)
#' # And the sequence of horizons to fit for
#' model$kseq <- 1:6
#' # Now we can fit the model and get the score, as it is
#' rls_fit(model=model, data=D, scorefun=rmse, returnanalysis=FALSE)
#' # Or we can change the lambda
#' rls_fit(c(lambda=0.9), model, D, rmse, returnanalysis=FALSE)
#'
#' # This could be passed to optim() (or any optimizer, see forecastmodel$insert_prm()).
#' optim(c(lambda=0.98), rls_fit, model=model, data=D, scorefun=rmse, returnanalysis=FALSE,
#'       lower=c(lambda=0.9), upper=c(lambda=0.999), method="L-BFGS-B")
#'
#' # rls_optim is simply a helper, it's makes using bounds easiere and enables caching of the results
#' # First add bounds for lambda (lower, init, upper)
#' model$add_prmbounds(lambda = c(0.9, 0.98, 0.999))
#'
#' # Now the same optimization as above can be done by
#' val <- rls_optim(model, D)
#' val
#' 
#' 
#' @export
rls_optim <- function(model, data, kseq = NA, scorefun = rmse, cachedir="", cachererun=FALSE, printout=TRUE, method="L-BFGS-B", ...){
    # Take the parameters bounds from the parameter bounds set in the model
    init <- model$get_prmbounds("init")
    lower <- model$get_prmbounds("lower")
    upper <- model$get_prmbounds("upper")
    # If bounds are NA, then set
    if(any(is.na(lower))){ lower[is.na(lower)] <- -Inf}
    if(any(is.na(upper))){ lower[is.na(upper)] <- Inf}

    # Clone the model no matter what (at least model$kseq should not be changed no matter if optimization is stopped)
    m <- model$clone_deep()
    if(!is.na(kseq[1])){
        m$kseq <- kseq
    }

    # Caching the results based on some of the function arguments
    if(cachedir != ""){
        # Have to insert the parameters in the expressions to get the right state of the model for unique checksum
        m$insert_prm(init)
        # Have to reset the state first to remove dependency of previous calls
        m$reset_state()
        # Give all the elements needed to calculate the unique cache name
        # This is maybe smarter, don't have to calculate the transformation of the data: cnm <- cache_name(m$regprm, getse(m$inputs, nms="expr"), m$output, m$prmbounds, m$kseq, data, objfun, init, lower, upper, cachedir = cachedir)
        cnm <- cache_name(rls_fit, rls_optim, m$outputrange, m$regprm, m$transform_data(data), data[[m$output]], scorefun, init, lower, upper, kseq, cachedir = cachedir)
        # Load the cached result if it exists
        if(file.exists(cnm) & !cachererun){
            res <- readRDS(cnm)
            # Set the optimized parameters into the model
            model$insert_prm(res$par)
            return(res)
        }
    }

    # Run the optimization
    res <- optim(par = init,
                 fn = rls_fit,
                 # Parameters to pass to rls_fit
                 model = m,
                 data = data,
                 scorefun = scorefun,
                 printout = printout,
                 returnanalysis = FALSE,
                 # Parameters to pass to optim
                 lower = lower,
                 upper = upper,
                 method =  method,
                 ...)        
    # Save the result in the cachedir
    if(cachedir != ""){ cache_save(res, cnm)}
    # Set the optimized parameters into the model
    model$insert_prm(res$par)
    return(res)
}