diff --git a/.github/workflows/R-CMD-check-Dev.yaml b/.github/workflows/R-CMD-check-Dev.yaml deleted file mode 100644 index 763fcf2..0000000 --- a/.github/workflows/R-CMD-check-Dev.yaml +++ /dev/null @@ -1,53 +0,0 @@ -# For help debugging build failures open an issue on the RStudio community with the 'github-actions' tag. -# https://community.rstudio.com/new-topic?category=Package%20development&tags=github-actions -on: [push, pull_request] - -name: R-CMD-check-Dev - -jobs: - R-CMD-check: - runs-on: ${{ matrix.config.os }} - - name: ${{ matrix.config.os }} (${{ matrix.config.r }}) - - strategy: - fail-fast: false - matrix: - config: - - {os: macOS-latest, r: 'release'} - - {os: windows-latest, r: 'release'} - - {os: ubuntu-latest, r: 'devel', http-user-agent: 'release'} - - {os: ubuntu-latest, r: 'release'} - env: - GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }} - R_KEEP_PKG_SOURCE: yes - - steps: - - uses: actions/checkout@v4 - - - uses: r-lib/actions/setup-pandoc@v2 - - - uses: r-lib/actions/setup-r@v2 - with: - r-version: ${{ matrix.config.r }} - http-user-agent: ${{ matrix.config.http-user-agent }} - use-public-rspm: true - - - uses: r-lib/actions/setup-r-dependencies@v2 - with: - extra-packages: any::rcmdcheck - needs: check - - - uses: r-lib/actions/check-r-package@v2 - - - name: Show testthat output - if: always() - run: find check -name 'testthat.Rout*' -exec cat '{}' \; || true - shell: bash - - - name: Upload check results - if: failure() - uses: actions/upload-artifact@main - with: - name: ${{ runner.os }}-r${{ matrix.config.r }}-results - path: check \ No newline at end of file diff --git a/.github/workflows/R-CMD-check.yaml b/.github/workflows/R-CMD-check.yaml deleted file mode 100644 index 68716aa..0000000 --- a/.github/workflows/R-CMD-check.yaml +++ /dev/null @@ -1,54 +0,0 @@ -# For help debugging build failures open an issue on the RStudio community with the 'github-actions' tag. -# https://community.rstudio.com/new-topic?category=Package%20development&tags=github-actions -on: - pull_request - -name: R-CMD-check - -jobs: - R-CMD-check: - runs-on: ${{ matrix.config.os }} - - name: ${{ matrix.config.os }} (${{ matrix.config.r }}) - - strategy: - fail-fast: false - matrix: - config: - - {os: macOS-latest, r: 'release'} - - {os: windows-latest, r: 'release'} - - {os: ubuntu-latest, r: 'release'} - - env: - GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }} - R_KEEP_PKG_SOURCE: yes - - steps: - - uses: actions/checkout@v4 - - - uses: r-lib/actions/setup-pandoc@v2 - - - uses: r-lib/actions/setup-r@v2 - with: - r-version: ${{ matrix.config.r }} - http-user-agent: ${{ matrix.config.http-user-agent }} - use-public-rspm: true - - - uses: r-lib/actions/setup-r-dependencies@v2 - with: - extra-packages: any::rcmdcheck - needs: check - - - uses: r-lib/actions/check-r-package@v2 - - - name: Show testthat output - if: always() - run: find check -name 'testthat.Rout*' -exec cat '{}' \; || true - shell: bash - - - name: Upload check results - if: failure() - uses: actions/upload-artifact@main - with: - name: ${{ runner.os }}-r${{ matrix.config.r }}-results - path: check \ No newline at end of file diff --git a/.github/workflows/R-hub-PR.yaml b/.github/workflows/R-hub-PR.yaml new file mode 100644 index 0000000..19935bb --- /dev/null +++ b/.github/workflows/R-hub-PR.yaml @@ -0,0 +1,86 @@ +# R-hub's generic GitHub Actions workflow file. It's canonical location is at +# https://github.com/r-hub/actions/blob/v1/workflows/rhub.yaml +# You can update this file to a newer version using the rhub2 package: +# +# rhub::rhub_setup() +# +# It is unlikely that you need to modify this file manually. +on: + pull_request + +name: R-hub-PR + +jobs: + + setup: + runs-on: ubuntu-latest + outputs: + containers: ${{ steps.rhub-setup.outputs.containers }} + platforms: ${{ steps.rhub-setup.outputs.platforms }} + + steps: + # NO NEED TO CHECKOUT HERE + - uses: r-hub/actions/setup@v1 + with: + config: 'linux, macos, windows, macos-arm64, clang-asan' + id: rhub-setup + + linux-containers: + needs: setup + if: ${{ needs.setup.outputs.containers != '[]' }} + runs-on: ubuntu-latest + name: 'R-hub' + strategy: + fail-fast: false + matrix: + config: ${{ fromJson(needs.setup.outputs.containers) }} + container: + image: ${{ matrix.config.container }} + + steps: + - uses: actions/checkout@v4 + with: + submodules: 'true' + - uses: r-hub/actions/platform-info@v1 + with: + token: ${{ secrets.RHUB_TOKEN }} + job-config: ${{ matrix.config.job-config }} + - uses: r-hub/actions/setup-deps@v1 + with: + token: ${{ secrets.RHUB_TOKEN }} + job-config: ${{ matrix.config.job-config }} + - uses: r-hub/actions/run-check@v1 + with: + token: ${{ secrets.RHUB_TOKEN }} + job-config: ${{ matrix.config.job-config }} + + other-platforms: + needs: setup + if: ${{ needs.setup.outputs.platforms != '[]' }} + runs-on: ${{ matrix.config.os }} + name: 'R-hub' + strategy: + fail-fast: false + matrix: + config: ${{ fromJson(needs.setup.outputs.platforms) }} + + steps: + - uses: actions/checkout@v4 + with: + submodules: 'true' + - uses: r-hub/actions/setup-r@v1 + with: + job-config: ${{ matrix.config.job-config }} + token: ${{ secrets.RHUB_TOKEN }} + - uses: r-hub/actions/platform-info@v1 + with: + token: ${{ secrets.RHUB_TOKEN }} + job-config: ${{ matrix.config.job-config }} + - uses: r-hub/actions/setup-deps@v1 + with: + job-config: ${{ matrix.config.job-config }} + token: ${{ secrets.RHUB_TOKEN }} + - uses: r-hub/actions/run-check@v1 + with: + job-config: ${{ matrix.config.job-config }} + token: ${{ secrets.RHUB_TOKEN }} diff --git a/.github/workflows/R-hub-main.yaml b/.github/workflows/R-hub-main.yaml new file mode 100644 index 0000000..bbb036f --- /dev/null +++ b/.github/workflows/R-hub-main.yaml @@ -0,0 +1,90 @@ +# R-hub's generic GitHub Actions workflow file. It's canonical location is at +# https://github.com/r-hub/actions/blob/v1/workflows/rhub.yaml +# You can update this file to a newer version using the rhub2 package: +# +# rhub::rhub_setup() +# +# It is unlikely that you need to modify this file manually. +on: + push: + branches: + - main + schedule: + - cron: '32 4 1 * *' + +name: R-hub-main + +jobs: + + setup: + runs-on: ubuntu-latest + outputs: + containers: ${{ steps.rhub-setup.outputs.containers }} + platforms: ${{ steps.rhub-setup.outputs.platforms }} + + steps: + # NO NEED TO CHECKOUT HERE + - uses: r-hub/actions/setup@v1 + with: + config: 'linux, macos, windows, macos-arm64, clang-asan' + id: rhub-setup + + linux-containers: + needs: setup + if: ${{ needs.setup.outputs.containers != '[]' }} + runs-on: ubuntu-latest + name: 'R-hub' + strategy: + fail-fast: false + matrix: + config: ${{ fromJson(needs.setup.outputs.containers) }} + container: + image: ${{ matrix.config.container }} + + steps: + - uses: actions/checkout@v4 + with: + submodules: 'true' + - uses: r-hub/actions/platform-info@v1 + with: + token: ${{ secrets.RHUB_TOKEN }} + job-config: ${{ matrix.config.job-config }} + - uses: r-hub/actions/setup-deps@v1 + with: + token: ${{ secrets.RHUB_TOKEN }} + job-config: ${{ matrix.config.job-config }} + - uses: r-hub/actions/run-check@v1 + with: + token: ${{ secrets.RHUB_TOKEN }} + job-config: ${{ matrix.config.job-config }} + + other-platforms: + needs: setup + if: ${{ needs.setup.outputs.platforms != '[]' }} + runs-on: ${{ matrix.config.os }} + name: 'R-hub' + strategy: + fail-fast: false + matrix: + config: ${{ fromJson(needs.setup.outputs.platforms) }} + + steps: + - uses: actions/checkout@v4 + with: + submodules: 'true' + - uses: r-hub/actions/setup-r@v1 + with: + job-config: ${{ matrix.config.job-config }} + token: ${{ secrets.RHUB_TOKEN }} + - uses: r-hub/actions/platform-info@v1 + with: + token: ${{ secrets.RHUB_TOKEN }} + job-config: ${{ matrix.config.job-config }} + - uses: r-hub/actions/setup-deps@v1 + with: + job-config: ${{ matrix.config.job-config }} + token: ${{ secrets.RHUB_TOKEN }} + - uses: r-hub/actions/run-check@v1 + with: + job-config: ${{ matrix.config.job-config }} + token: ${{ secrets.RHUB_TOKEN }} diff --git a/.github/workflows/recheck.yaml b/.github/workflows/recheck.yaml new file mode 100644 index 0000000..983ea80 --- /dev/null +++ b/.github/workflows/recheck.yaml @@ -0,0 +1,147 @@ +on: + pull_request: + branches: + - main + paths: + - 'NEWS.md' + +name: Reverse dependency check + +env: + R_LIBS_USER: ${{github.workspace}}/pkglib + _R_CHECK_FORCE_SUGGESTS_: 'FALSE' + _R_CHECK_CRAN_INCOMING_: 'FALSE' + _R_CHECK_ELAPSED_TIMEOUT_: 3600 + GITHUB_ACTIONS: '' + CI: '' + +jobs: + prepare: + name: Prepare dependencies + runs-on: ubuntu-latest + container: ghcr.io/r-devel/recheck + outputs: + oldfile: ${{ steps.filenames.outputs.oldfile }} + newfile: ${{ steps.filenames.outputs.newfile }} + steps: + - name: prepare + run: | + mkdir -p $R_LIBS_USER + R -e ".libPaths()" + + - name: checkout + uses: actions/checkout@v4 + with: + path: source + submodules: 'true' + + - name: Set package source directory + run: | + if [ "${{ inputs.subdirectory }}" ]; then + echo "PKG_SOURCE_DIR=source/${{ inputs.subdirectory }}" >> $GITHUB_ENV + else + echo "PKG_SOURCE_DIR=source" >> $GITHUB_ENV + fi + + - name: download dependencies + run: rechecktools::install_recheck_deps('${{env.PKG_SOURCE_DIR}}', 'strong') + shell: Rscript {0} + + - name: build source package + run: | + mkdir newpkg + R CMD build ${{env.PKG_SOURCE_DIR}} + mv *.tar.gz newpkg/ + rm -Rf source + + - name: Get old version of package + shell: Rscript {0} + run: | + dir.create("oldpkg") + pkg <- sub("_.*", "", list.files("newpkg")) + download.packages(pkg, "oldpkg", repos = "https://cran.r-project.org") + + - name: Get file names + id: filenames + run: | + echo "oldfile=$(cd oldpkg; echo *.tar.gz)" >> "$GITHUB_OUTPUT" + echo "newfile=$(cd newpkg; echo *.tar.gz)" >> "$GITHUB_OUTPUT" + + - name: Save package and library + uses: actions/cache/save@v4 + with: + path: | + pkglib + newpkg + oldpkg + key: ${{ runner.os }}-${{ github.run_id }}-${{github.run_attempt}} + + checks: + needs: prepare + runs-on: ubuntu-latest + name: Recheck ${{matrix.check == 'oldpkg' && needs.prepare.outputs.oldfile || needs.prepare.outputs.newfile}} (${{matrix.check}}) + container: ghcr.io/r-devel/recheck + timeout-minutes: 600 + strategy: + matrix: + check: [ 'oldpkg', 'newpkg' ] + steps: + + - name: Download package and library + uses: actions/cache/restore@v4 + with: + path: | + pkglib + newpkg + oldpkg + key: ${{ runner.os }}-${{ github.run_id }}-${{github.run_attempt}} + + - name: Run checks + shell: Rscript {0} + run: | + checkdir <- "${{matrix.check}}" + options(repos = c(CRAN = 'https://cloud.r-project.org')) + reverse <- list(repos = 'https://cloud.r-project.org', which = "strong") + Sys.setenv(R_PROFILE="nothing") # do not set binary p3m mirror in child proc + tools::check_packages_in_dir(checkdir, reverse = reverse, Ncpus = parallel::detectCores()) + details <- tools::check_packages_in_dir_details(checkdir) + write.csv(details, file.path(checkdir, 'check-details.csv')) + writeLines(paste(format(details), collapse = "\n\n"), file.path(checkdir, 'check-details.txt')) + tools::summarize_check_packages_in_dir_results(checkdir) + + - uses: actions/upload-artifact@v4 + with: + name: ${{matrix.check}}-checklogs + path: | + ${{matrix.check}}/*/*.log + ${{matrix.check}}/*/*.out + ${{matrix.check}}/*/*.Rout + ${{matrix.check}}/*/*.fail + ${{matrix.check}}/Outputs + ${{matrix.check}}/check-details.* + + results: + name: Show results + needs: checks + runs-on: ubuntu-latest + container: ghcr.io/r-hub/r-minimal/r-minimal:latest + steps: + - name: Download log files + uses: actions/download-artifact@v4 + + - name: Get results + shell: Rscript {0} + run: | + cat("\n------- Check results summary ------\n") + tools::summarize_check_packages_in_dir_results("newpkg-checklogs") + + # Check for regressions + cat("\n------- Check for regressions ------\n") + changes <- tools::check_packages_in_dir_changes("newpkg-checklogs", "oldpkg-checklogs") + if(nrow(changes)){ + cat("Changes:\n") + print(changes) + stop("Found potential problems") + } else { + cat("No changes between old and new version\n") + } \ No newline at end of file diff --git a/CRAN-SUBMISSION b/CRAN-SUBMISSION index be99a60..36e1553 100644 --- a/CRAN-SUBMISSION +++ b/CRAN-SUBMISSION @@ -1,3 +1,3 @@ -Version: 1.3.2 -Date: 2024-03-25 11:09:04 UTC -SHA: 9927d26dfeb95922b6f15aae7b303b90a329d75e +Version: 1.3.3 +Date: 2024-09-21 21:57:06 UTC +SHA: e3dd1646c90053393342636b543e96f65eb8dec8 diff --git a/DESCRIPTION b/DESCRIPTION index a7829ad..05abc67 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,8 +1,8 @@ Package: profoc Type: Package Title: Probabilistic Forecast Combination Using CRPS Learning -Version: 1.3.2 -Date: 2024-03-25 +Version: 1.3.3 +Date: 2024-09-21 Authors@R: c( person(given = "Jonathan", family = "Berrisch", @@ -14,7 +14,7 @@ Authors@R: c( role = "aut", comment = c(ORCID = "0000-0002-2974-2660")) ) -Description: Combine probabilistic forecasts using CRPS learning algorithms proposed in Berrisch, Ziel (2021) . The package implements multiple online learning algorithms like Bernstein online aggregation; see Wintenberger (2014) . Quantile regression is also implemented for comparison purposes. Model parameters can be tuned automatically with respect to the loss of the forecast combination. Methods like predict(), update(), plot() and print() are available for convenience. This package utilizes the optim C++ library for numeric optimization . +Description: Combine probabilistic forecasts using CRPS learning algorithms proposed in Berrisch, Ziel (2021) . The package implements multiple online learning algorithms like Bernstein online aggregation; see Wintenberger (2014) . Quantile regression is also implemented for comparison purposes. Model parameters can be tuned automatically with respect to the loss of the forecast combination. Methods like predict(), update(), plot() and print() are available for convenience. This package utilizes the optim C++ library for numeric optimization . License: GPL (>= 3) Encoding: UTF-8 Depends: R (>= 4.3.0) @@ -27,17 +27,18 @@ Imports: generics, tibble, ggplot2 -LinkingTo: Rcpp, RcppArmadillo (>= 0.10.7.5.0), RcppProgress, splines2 (>= 0.4.4), rcpptimer (>= 1.1.0) +LinkingTo: Rcpp, RcppArmadillo (>= 0.10.7.5.0), RcppProgress, splines2 (>= 0.4.4), rcpptimer (>= 1.2.0) URL: https://profoc.berrisch.biz, https://github.com/BerriJ/profoc BugReports: https://github.com/BerriJ/profoc/issues -RoxygenNote: 7.3.1 +RoxygenNote: 7.3.2 Language: en-US Suggests: testthat (>= 3.0.0), gamlss.dist, knitr, rmarkdown, - dplyr + dplyr, + rcpptimer (>= 1.2.0), Config/testthat/edition: 3 Roxygen: list(markdown = TRUE) VignetteBuilder: knitr diff --git a/NEWS.md b/NEWS.md index 3a98835..f03a55f 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,10 @@ +profoc 1.3.3 +============== + +## Improvements + +* We adjusted the integration of `rcpptimer`. This simplifies the code and makes use of the API of `rcpptimer` 1.2.0 which is expected to be stable. + profoc 1.3.2 ============== diff --git a/R/online.R b/R/online.R index f1b56b5..f51ed6c 100644 --- a/R/online.R +++ b/R/online.R @@ -81,15 +81,15 @@ #' y <- rnorm(n = T) # Realized #' experts <- array(dim = c(T, P, N)) # Predictions #' for (t in 1:T) { -#' experts[t, , 1] <- qnorm(prob_grid, mean = -1, sd = 1) -#' experts[t, , 2] <- qnorm(prob_grid, mean = 3, sd = sqrt(4)) +#' experts[t, , 1] <- qnorm(prob_grid, mean = -1, sd = 1) +#' experts[t, , 2] <- qnorm(prob_grid, mean = 3, sd = sqrt(4)) #' } #' #' model <- online( -#' y = matrix(y), -#' experts = experts, -#' tau = prob_grid, -#' p_smooth_pr = list(lambda = 10) +#' y = matrix(y), +#' experts = experts, +#' tau = prob_grid, +#' p_smooth_pr = list(lambda = 10) #' ) #' #' print(model) @@ -115,44 +115,44 @@ online <- function(y, experts, tau, loss_gradient = TRUE, method = "bewa", b_smooth_pr = list( - knots = P, - mu = 0.5, - sigma = 1, - nonc = 0, - tailweight = 1, - deg = 1, - periodic = FALSE + knots = P, + mu = 0.5, + sigma = 1, + nonc = 0, + tailweight = 1, + deg = 1, + periodic = FALSE ), p_smooth_pr = list( - knots = P, - mu = 0.5, - sigma = 1, - nonc = 0, - tailweight = 1, - deg = 1, - ndiff = 1.5, - lambda = -Inf, - periodic = FALSE + knots = P, + mu = 0.5, + sigma = 1, + nonc = 0, + tailweight = 1, + deg = 1, + ndiff = 1.5, + lambda = -Inf, + periodic = FALSE ), b_smooth_mv = list( - knots = D, - mu = 0.5, - sigma = 1, - nonc = 0, - tailweight = 1, - deg = 1, - periodic = FALSE + knots = D, + mu = 0.5, + sigma = 1, + nonc = 0, + tailweight = 1, + deg = 1, + periodic = FALSE ), p_smooth_mv = list( - knots = D, - mu = 0.5, - sigma = 1, - nonc = 0, - tailweight = 1, - deg = 1, - ndiff = 1.5, - lambda = -Inf, - periodic = FALSE + knots = D, + mu = 0.5, + sigma = 1, + nonc = 0, + tailweight = 1, + deg = 1, + ndiff = 1.5, + lambda = -Inf, + periodic = FALSE ), forget_regret = 0, soft_threshold = -Inf, @@ -161,11 +161,11 @@ online <- function(y, experts, tau, gamma = 1, parametergrid_max_combinations = 100, parametergrids = list( - general = NULL, - b_smooth_pr = NULL, - p_smooth_pr = NULL, - b_smooth_mv = NULL, - p_smooth_mv = NULL + general = NULL, + b_smooth_pr = NULL, + p_smooth_pr = NULL, + b_smooth_mv = NULL, + p_smooth_mv = NULL ), forget_past_performance = 0, save_past_performance = FALSE, @@ -176,292 +176,296 @@ online <- function(y, experts, tau, regret = NULL, trace = TRUE, get_timings = FALSE) { - model_instance <- new(conline) - model_instance$trace <- trace - model_instance$get_timings <- get_timings - model_instance$forget_past_performance <- forget_past_performance - model_instance$save_past_performance <- save_past_performance - model_instance$save_predictions_grid <- save_predictions_grid - - if (is.vector(y)) { - y <- matrix(y) + if (requireNamespace("rcpptimer", quietly = TRUE)) { + registerS3method("print", "rcpptimer", rcpptimer:::print.rcpptimer) + } + + model_instance <- new(conline) + model_instance$trace <- trace + model_instance$get_timings <- get_timings + model_instance$forget_past_performance <- forget_past_performance + model_instance$save_past_performance <- save_past_performance + model_instance$save_predictions_grid <- save_predictions_grid + + if (is.vector(y)) { + y <- matrix(y) + } + + model_instance$y <- y + + # preserve names + names <- list(y = dimnames(y)) + + # Prepare experts + + e_list <- init_experts_list( + experts = experts, + y = y, + output_with_names = TRUE + ) + model_instance$experts <- e_list$experts + + names$experts <- list(NULL) + names$experts[[2]] <- e_list$dnames + names$experts[[3]] <- tau + names$experts[[4]] <- e_list$enames + + model_instance$tau <- tau + + # Define dimensions for convenience + T <- dim(e_list$experts)[1] + D <- dim(e_list$experts[[1]])[1] + P <- dim(e_list$experts[[1]])[2] + K <- dim(e_list$experts[[1]])[3] + + if (nrow(e_list$experts) - nrow(y) < 0) { + stop("Number of provided expert predictions has to match or exceed observations.") + } + + if (nrow(y) <= lead_time) { + stop("Number of expert predictions need to exceed lead_time.") + } + + if (is.null(names$y[[2]])) names$y[[2]] <- paste0(1:D, "D") + + model_instance$lead_time <- lead_time + + if (is.null(loss)) { + loss_share <- 0 + } else if (is.array(loss)) { + loss_array <- list() + for (i in 1:T) { + loss_array[[i]] <- array(loss[i, , , ], dim = c(D, P, K)) } - - model_instance$y <- y - - # preserve names - names <- list(y = dimnames(y)) - - # Prepare experts - - e_list <- init_experts_list( - experts = experts, - y = y, - output_with_names = TRUE - ) - model_instance$experts <- e_list$experts - - names$experts <- list(NULL) - names$experts[[2]] <- e_list$dnames - names$experts[[3]] <- tau - names$experts[[4]] <- e_list$enames - - model_instance$tau <- tau - - # Define dimensions for convenience - T <- dim(e_list$experts)[1] - D <- dim(e_list$experts[[1]])[1] - P <- dim(e_list$experts[[1]])[2] - K <- dim(e_list$experts[[1]])[3] - - if (nrow(e_list$experts) - nrow(y) < 0) { - stop("Number of provided expert predictions has to match or exceed observations.") + dim(loss_array) <- c(T, 1) + loss_share <- 1 + model_instance$loss_array <- loss_array + } else if (is.list(loss)) { + loss_array <- list() + for (i in 1:T) { + loss_array[[i]] <- array(loss$loss[i, , , ], dim = c(D, P, K)) } - - if (nrow(y) <= lead_time) { - stop("Number of expert predictions need to exceed lead_time.") - } - - if (is.null(names$y[[2]])) names$y[[2]] <- paste0(1:D, "D") - - model_instance$lead_time <- lead_time - - if (is.null(loss)) { - loss_share <- 0 - } else if (is.array(loss)) { - loss_array <- list() - for (i in 1:T) { - loss_array[[i]] <- array(loss[i, , , ], dim = c(D, P, K)) - } - dim(loss_array) <- c(T, 1) - loss_share <- 1 - model_instance$loss_array <- loss_array - } else if (is.list(loss)) { - loss_array <- list() - for (i in 1:T) { - loss_array[[i]] <- array(loss$loss[i, , , ], dim = c(D, P, K)) - } - dim(loss_array) <- c(T, 1) - loss_share <- loss$share - model_instance$loss_array <- loss_array + dim(loss_array) <- c(T, 1) + loss_share <- loss$share + model_instance$loss_array <- loss_array + } + + if (is.null(regret)) { + regret_share <- 0 + } else if (is.array(regret)) { + regret_array <- list() + for (i in 1:T) { + regret_array[[i]] <- array(regret[i, , , ], dim = c(D, P, K)) } - - if (is.null(regret)) { - regret_share <- 0 - } else if (is.array(regret)) { - regret_array <- list() - for (i in 1:T) { - regret_array[[i]] <- array(regret[i, , , ], dim = c(D, P, K)) - } - dim(regret_array) <- c(T, 1) - regret_share <- 1 - model_instance$regret_array <- regret_array - } else if (is.list(regret)) { - regret_array <- list() - for (i in 1:T) { - regret_array[[i]] <- array(regret$regret[i, , , ], - dim = c(D, P, K) - ) - } - dim(regret_array) <- c(T, 1) - regret_share <- regret$share - model_instance$regret_array <- regret_array + dim(regret_array) <- c(T, 1) + regret_share <- 1 + model_instance$regret_array <- regret_array + } else if (is.list(regret)) { + regret_array <- list() + for (i in 1:T) { + regret_array[[i]] <- array(regret$regret[i, , , ], + dim = c(D, P, K) + ) } - - # Create parameter grid - pars_basis_pr <- list( - n = val_or_def(b_smooth_pr$knots, P), - mu = val_or_def(b_smooth_pr$mu, 0.5), - sigma = val_or_def(b_smooth_pr$sigma, 1), - nonc = val_or_def(b_smooth_pr$nonc, 0), - tailw = val_or_def(b_smooth_pr$tailweight, 1), - deg = val_or_def(b_smooth_pr$deg, 1), - periodic = val_or_def(b_smooth_pr$periodic, FALSE) - ) - pars_basis_pr_n <- prod(sapply(pars_basis_pr, length)) - - pars_basis_mv <- list( - n = val_or_def(b_smooth_mv$knots, D), - mu = val_or_def(b_smooth_mv$mu, 0.5), - sigma = val_or_def(b_smooth_mv$sigma, 1), - nonc = val_or_def(b_smooth_mv$nonc, 0), - tailw = val_or_def(b_smooth_mv$tailweight, 1), - deg = val_or_def(b_smooth_mv$deg, 1), - periodic = val_or_def(b_smooth_mv$periodic, FALSE) - ) - pars_basis_mv_n <- prod(sapply(pars_basis_mv, length)) - - pars_hat_pr <- list( - n = val_or_def(p_smooth_pr$knots, P), - mu = val_or_def(p_smooth_pr$mu, 0.5), - sigma = val_or_def(p_smooth_pr$sigma, 1), - nonc = val_or_def(p_smooth_pr$nonc, 0), - tailw = val_or_def(p_smooth_pr$tailweight, 1), - deg = val_or_def(p_smooth_pr$deg, 1), - ndiff = val_or_def(p_smooth_pr$ndiff, 1.5), - lambda = val_or_def(p_smooth_pr$lambda, -Inf), - periodic = val_or_def(p_smooth_pr$periodic, FALSE) - ) - pars_hat_pr_n <- prod(sapply(pars_hat_pr, length)) - - pars_hat_mv <- list( - n = val_or_def(p_smooth_mv$knots, D), - mu = val_or_def(p_smooth_mv$mu, 0.5), - sigma = val_or_def(p_smooth_mv$sigma, 1), - nonc = val_or_def(p_smooth_mv$nonc, 0), - tailw = val_or_def(p_smooth_mv$tailweight, 1), - deg = val_or_def(p_smooth_mv$deg, 1), - ndiff = val_or_def(p_smooth_mv$ndiff, 1.5), - lambda = val_or_def(p_smooth_mv$lambda, -Inf), - periodic = val_or_def(p_smooth_mv$periodic, FALSE) + dim(regret_array) <- c(T, 1) + regret_share <- regret$share + model_instance$regret_array <- regret_array + } + + # Create parameter grid + pars_basis_pr <- list( + n = val_or_def(b_smooth_pr$knots, P), + mu = val_or_def(b_smooth_pr$mu, 0.5), + sigma = val_or_def(b_smooth_pr$sigma, 1), + nonc = val_or_def(b_smooth_pr$nonc, 0), + tailw = val_or_def(b_smooth_pr$tailweight, 1), + deg = val_or_def(b_smooth_pr$deg, 1), + periodic = val_or_def(b_smooth_pr$periodic, FALSE) + ) + pars_basis_pr_n <- prod(sapply(pars_basis_pr, length)) + + pars_basis_mv <- list( + n = val_or_def(b_smooth_mv$knots, D), + mu = val_or_def(b_smooth_mv$mu, 0.5), + sigma = val_or_def(b_smooth_mv$sigma, 1), + nonc = val_or_def(b_smooth_mv$nonc, 0), + tailw = val_or_def(b_smooth_mv$tailweight, 1), + deg = val_or_def(b_smooth_mv$deg, 1), + periodic = val_or_def(b_smooth_mv$periodic, FALSE) + ) + pars_basis_mv_n <- prod(sapply(pars_basis_mv, length)) + + pars_hat_pr <- list( + n = val_or_def(p_smooth_pr$knots, P), + mu = val_or_def(p_smooth_pr$mu, 0.5), + sigma = val_or_def(p_smooth_pr$sigma, 1), + nonc = val_or_def(p_smooth_pr$nonc, 0), + tailw = val_or_def(p_smooth_pr$tailweight, 1), + deg = val_or_def(p_smooth_pr$deg, 1), + ndiff = val_or_def(p_smooth_pr$ndiff, 1.5), + lambda = val_or_def(p_smooth_pr$lambda, -Inf), + periodic = val_or_def(p_smooth_pr$periodic, FALSE) + ) + pars_hat_pr_n <- prod(sapply(pars_hat_pr, length)) + + pars_hat_mv <- list( + n = val_or_def(p_smooth_mv$knots, D), + mu = val_or_def(p_smooth_mv$mu, 0.5), + sigma = val_or_def(p_smooth_mv$sigma, 1), + nonc = val_or_def(p_smooth_mv$nonc, 0), + tailw = val_or_def(p_smooth_mv$tailweight, 1), + deg = val_or_def(p_smooth_mv$deg, 1), + ndiff = val_or_def(p_smooth_mv$ndiff, 1.5), + lambda = val_or_def(p_smooth_mv$lambda, -Inf), + periodic = val_or_def(p_smooth_mv$periodic, FALSE) + ) + pars_hat_mv_n <- prod(sapply(pars_hat_mv, length)) + + if (is.null(parametergrids$general)) { + parametergrid <- expand_grid_sample( + list( + forget_regret = forget_regret, + soft_threshold = soft_threshold, + hard_threshold = hard_threshold, + fixed_share = fixed_share, + basis_pr_idx = seq_len(pars_basis_pr_n), + basis_mv_idx = seq_len(pars_basis_mv_n), + hat_pr_idx = seq_len(pars_hat_pr_n), + hat_mv_idx = seq_len(pars_hat_mv_n), + gamma = gamma, + loss_share = loss_share, + regret_share = regret_share + ), + n = parametergrid_max_combinations, + verbose = TRUE ) - pars_hat_mv_n <- prod(sapply(pars_hat_mv, length)) - - if (is.null(parametergrids$general)) { - parametergrid <- expand_grid_sample( - list( - forget_regret = forget_regret, - soft_threshold = soft_threshold, - hard_threshold = hard_threshold, - fixed_share = fixed_share, - basis_pr_idx = seq_len(pars_basis_pr_n), - basis_mv_idx = seq_len(pars_basis_mv_n), - hat_pr_idx = seq_len(pars_hat_pr_n), - hat_mv_idx = seq_len(pars_hat_mv_n), - gamma = gamma, - loss_share = loss_share, - regret_share = regret_share - ), - n = parametergrid_max_combinations, - verbose = TRUE - ) - } else { - parametergrid <- parametergrids$general + } else { + parametergrid <- parametergrids$general + } + + # Create basis and hat matix lists + + # # Basis matrices for probabilistic smoothing + tmp <- make_basis_mats( + x = 1:P / (P + 1), + n = val_or_def(b_smooth_pr$knots, P), + mu = val_or_def(b_smooth_pr$mu, 0.5), + sigma = val_or_def(b_smooth_pr$sigma, 1), + nonc = val_or_def(b_smooth_pr$nonc, 0), + tailw = val_or_def(b_smooth_pr$tailweight, 1), + deg = val_or_def(b_smooth_pr$deg, 1), + periodic = val_or_def(b_smooth_pr$periodic, FALSE), + idx = sort(unique(parametergrid[, "basis_pr_idx"])), + params = parametergrids$b_smooth_pr + ) + model_instance$basis_pr <- tmp$basis + model_instance$params_basis_pr <- tmp$params + + # Basis matrices for multivariate smoothing + tmp <- make_basis_mats( + x = 1:D / (D + 1), + n = val_or_def(b_smooth_mv$knots, D), + mu = val_or_def(b_smooth_mv$mu, 0.5), + sigma = val_or_def(b_smooth_mv$sigma, 1), + nonc = val_or_def(b_smooth_mv$nonc, 0), + tailw = val_or_def(b_smooth_mv$tailweight, 1), + deg = val_or_def(b_smooth_mv$deg, 1), + periodic = val_or_def(b_smooth_mv$periodic, FALSE), + idx = sort(unique(parametergrid[, "basis_mv_idx"])), + params = parametergrids$b_smooth_mv + ) + model_instance$basis_mv <- tmp$basis + model_instance$params_basis_mv <- tmp$params + + tmp <- make_hat_mats( + x = 1:P / (P + 1), + n = val_or_def(p_smooth_pr$knots, P), + mu = val_or_def(p_smooth_pr$mu, 0.5), + sigma = val_or_def(p_smooth_pr$sigma, 1), + nonc = val_or_def(p_smooth_pr$nonc, 0), + tailw = val_or_def(p_smooth_pr$tailweight, 1), + deg = val_or_def(p_smooth_pr$deg, 1), + ndiff = val_or_def(p_smooth_pr$ndiff, 1.5), + lambda = val_or_def(p_smooth_pr$lambda, -Inf), + periodic = val_or_def(p_smooth_pr$periodic, FALSE), + idx = sort(unique(parametergrid[, "hat_pr_idx"])), + params = parametergrids$p_smooth_pr + ) + model_instance$hat_pr <- tmp$hat + model_instance$params_hat_pr <- tmp$params + + tmp <- make_hat_mats( + x = 1:D / (D + 1), + n = val_or_def(p_smooth_mv$knots, D), + mu = val_or_def(p_smooth_mv$mu, 0.5), + sigma = val_or_def(p_smooth_mv$sigma, 1), + nonc = val_or_def(p_smooth_mv$nonc, 0), + tailw = val_or_def(p_smooth_mv$tailweight, 1), + deg = val_or_def(p_smooth_mv$deg, 1), + ndiff = val_or_def(p_smooth_mv$ndiff, 1.5), + lambda = val_or_def(p_smooth_mv$lambda, -Inf), + periodic = val_or_def(p_smooth_mv$periodic, FALSE), + idx = sort(unique(parametergrid[, "hat_mv_idx"])), + params = parametergrids$p_smooth_mv + ) + model_instance$hat_mv <- tmp$hat + model_instance$params_hat_mv <- tmp$params + + # Fix basis / hat idx while retaining original order + parametergrid[, "basis_pr_idx"] <- match( + parametergrid[, "basis_pr_idx"], + sort(unique(parametergrid[, "basis_pr_idx"])) + ) + + parametergrid[, "basis_mv_idx"] <- match( + parametergrid[, "basis_mv_idx"], + sort(unique(parametergrid[, "basis_mv_idx"])) + ) + + parametergrid[, "hat_pr_idx"] <- match( + parametergrid[, "hat_pr_idx"], + sort(unique(parametergrid[, "hat_pr_idx"])) + ) + + parametergrid[, "hat_mv_idx"] <- match( + parametergrid[, "hat_mv_idx"], + sort(unique(parametergrid[, "hat_mv_idx"])) + ) + + model_instance$params <- parametergrid + model_instance$allow_quantile_crossing <- allow_quantile_crossing + model_instance$loss_function <- loss_function + model_instance$loss_gradient <- loss_gradient + model_instance$loss_parameter <- loss_parameter + model_instance$method <- method + + # Populate default values for w0, R0 etc. + model_instance$set_defaults() + + # Overwrite defaults if user specified explicitly + if (!is.null(init$init_weights)) { + init$init_weights <- pmax(init$init_weights, exp(-350)) + for (d in 1:D) { + for (p in 1:P) { + init$init_weights[d, p, ] <- + init$init_weights[d, p, ] / sum(init$init_weights[d, p, ]) + } } + model_instance$w0 <- init$init_weights + } - # Create basis and hat matix lists - - # # Basis matrices for probabilistic smoothing - tmp <- make_basis_mats( - x = 1:P / (P + 1), - n = val_or_def(b_smooth_pr$knots, P), - mu = val_or_def(b_smooth_pr$mu, 0.5), - sigma = val_or_def(b_smooth_pr$sigma, 1), - nonc = val_or_def(b_smooth_pr$nonc, 0), - tailw = val_or_def(b_smooth_pr$tailweight, 1), - deg = val_or_def(b_smooth_pr$deg, 1), - periodic = val_or_def(b_smooth_pr$periodic, FALSE), - idx = sort(unique(parametergrid[, "basis_pr_idx"])), - params = parametergrids$b_smooth_pr - ) - model_instance$basis_pr <- tmp$basis - model_instance$params_basis_pr <- tmp$params - - # Basis matrices for multivariate smoothing - tmp <- make_basis_mats( - x = 1:D / (D + 1), - n = val_or_def(b_smooth_mv$knots, D), - mu = val_or_def(b_smooth_mv$mu, 0.5), - sigma = val_or_def(b_smooth_mv$sigma, 1), - nonc = val_or_def(b_smooth_mv$nonc, 0), - tailw = val_or_def(b_smooth_mv$tailweight, 1), - deg = val_or_def(b_smooth_mv$deg, 1), - periodic = val_or_def(b_smooth_mv$periodic, FALSE), - idx = sort(unique(parametergrid[, "basis_mv_idx"])), - params = parametergrids$b_smooth_mv - ) - model_instance$basis_mv <- tmp$basis - model_instance$params_basis_mv <- tmp$params - - tmp <- make_hat_mats( - x = 1:P / (P + 1), - n = val_or_def(p_smooth_pr$knots, P), - mu = val_or_def(p_smooth_pr$mu, 0.5), - sigma = val_or_def(p_smooth_pr$sigma, 1), - nonc = val_or_def(p_smooth_pr$nonc, 0), - tailw = val_or_def(p_smooth_pr$tailweight, 1), - deg = val_or_def(p_smooth_pr$deg, 1), - ndiff = val_or_def(p_smooth_pr$ndiff, 1.5), - lambda = val_or_def(p_smooth_pr$lambda, -Inf), - periodic = val_or_def(p_smooth_pr$periodic, FALSE), - idx = sort(unique(parametergrid[, "hat_pr_idx"])), - params = parametergrids$p_smooth_pr - ) - model_instance$hat_pr <- tmp$hat - model_instance$params_hat_pr <- tmp$params - - tmp <- make_hat_mats( - x = 1:D / (D + 1), - n = val_or_def(p_smooth_mv$knots, D), - mu = val_or_def(p_smooth_mv$mu, 0.5), - sigma = val_or_def(p_smooth_mv$sigma, 1), - nonc = val_or_def(p_smooth_mv$nonc, 0), - tailw = val_or_def(p_smooth_mv$tailweight, 1), - deg = val_or_def(p_smooth_mv$deg, 1), - ndiff = val_or_def(p_smooth_mv$ndiff, 1.5), - lambda = val_or_def(p_smooth_mv$lambda, -Inf), - periodic = val_or_def(p_smooth_mv$periodic, FALSE), - idx = sort(unique(parametergrid[, "hat_mv_idx"])), - params = parametergrids$p_smooth_mv - ) - model_instance$hat_mv <- tmp$hat - model_instance$params_hat_mv <- tmp$params - - # Fix basis / hat idx while retaining original order - parametergrid[, "basis_pr_idx"] <- match( - parametergrid[, "basis_pr_idx"], - sort(unique(parametergrid[, "basis_pr_idx"])) - ) - - parametergrid[, "basis_mv_idx"] <- match( - parametergrid[, "basis_mv_idx"], - sort(unique(parametergrid[, "basis_mv_idx"])) - ) - - parametergrid[, "hat_pr_idx"] <- match( - parametergrid[, "hat_pr_idx"], - sort(unique(parametergrid[, "hat_pr_idx"])) - ) - - parametergrid[, "hat_mv_idx"] <- match( - parametergrid[, "hat_mv_idx"], - sort(unique(parametergrid[, "hat_mv_idx"])) - ) - - model_instance$params <- parametergrid - model_instance$allow_quantile_crossing <- allow_quantile_crossing - model_instance$loss_function <- loss_function - model_instance$loss_gradient <- loss_gradient - model_instance$loss_parameter <- loss_parameter - model_instance$method <- method - - # Populate default values for w0, R0 etc. - model_instance$set_defaults() - - # Overwrite defaults if user specified explicitly - if (!is.null(init$init_weights)) { - init$init_weights <- pmax(init$init_weights, exp(-350)) - for (d in 1:D) { - for (p in 1:P) { - init$init_weights[d, p, ] <- - init$init_weights[d, p, ] / sum(init$init_weights[d, p, ]) - } - } - model_instance$w0 <- init$init_weights - } - - if (!is.null(init$R0)) model_instance$R0 <- init$R0 - - # Populates default values to grid-sized objects - model_instance$set_grid_objects() + if (!is.null(init$R0)) model_instance$R0 <- init$R0 - # Execute online learning - model_instance$learn() + # Populates default values to grid-sized objects + model_instance$set_grid_objects() - model <- post_process_model(model_instance, names) + # Execute online learning + model_instance$learn() - model_instance$get_times() - rm(model_instance) + model <- post_process_model(model_instance, names) - return(model) + model_instance$get_times() + rm(model_instance) + gc() + return(model) } diff --git a/README.md b/README.md index 0e5059d..61f3914 100644 --- a/README.md +++ b/README.md @@ -5,12 +5,12 @@ An R package for probabilistic forecast combination ------------ -[![R-CMD-check](https://img.shields.io/github/actions/workflow/status/berrij/profoc/R-CMD-check.yaml?&style=for-the-badge)](https://github.com/BerriJ/profoc/actions/workflows/R-CMD-check.yaml) +[![R-hub](https://img.shields.io/github/actions/workflow/status/berrij/profoc/R-hub-main.yaml?&style=for-the-badge)](https://github.com/BerriJ/profoc/actions/workflows/R-hub-main.yaml) [![GitHub Workflow Status (branch)](https://img.shields.io/github/actions/workflow/status/berrij/profoc/pkgdown.yaml?label=Documentation&style=for-the-badge)](https://profoc.berrisch.biz/) [![Lifecycle: stable](https://img.shields.io/badge/Lifecycle-stable-brightgreen?style=for-the-badge)](https://lifecycle.r-lib.org/articles/stages.html#stable) -The primary function `online` can be used to combine probabilistic forecasts using the CRPS learning algorithm introduced in Berrisch, Ziel (2021): [Pre-Print](https://arxiv.org/pdf/2102.00968.pdf), [Publication](https://doi.org/10.1016/j.jeconom.2021.11.008). +The primary function `online` can be used to combine probabilistic forecasts using the CRPS learning algorithm introduced in Berrisch, Ziel (2021): [Pre-Print](https://arxiv.org/pdf/2102.00968), [Publication](https://doi.org/10.1016/j.jeconom.2021.11.008). The function `batch` can be used in a similar way for batch optimization. Common methods like `summary`, `print`, `plot`, `update`, and `predict` are available. Installation diff --git a/man/online.Rd b/man/online.Rd index fdc5744..1747b39 100644 --- a/man/online.Rd +++ b/man/online.Rd @@ -214,15 +214,15 @@ prob_grid <- 1:P / (P + 1) y <- rnorm(n = T) # Realized experts <- array(dim = c(T, P, N)) # Predictions for (t in 1:T) { - experts[t, , 1] <- qnorm(prob_grid, mean = -1, sd = 1) - experts[t, , 2] <- qnorm(prob_grid, mean = 3, sd = sqrt(4)) + experts[t, , 1] <- qnorm(prob_grid, mean = -1, sd = 1) + experts[t, , 2] <- qnorm(prob_grid, mean = 3, sd = sqrt(4)) } model <- online( - y = matrix(y), - experts = experts, - tau = prob_grid, - p_smooth_pr = list(lambda = 10) + y = matrix(y), + experts = experts, + tau = prob_grid, + p_smooth_pr = list(lambda = 10) ) print(model) diff --git a/src/conline.cpp b/src/conline.cpp index da3102a..cf9a558 100644 --- a/src/conline.cpp +++ b/src/conline.cpp @@ -17,15 +17,15 @@ using namespace arma; // [[Rcpp::export]] bool test_class_input(conline &obj) { - return obj.trace; + return obj.trace; } // [[Rcpp::export]] conline test_class_output() { - conline instance; - instance.trace = false; - return instance; + conline instance; + instance.trace = false; + return instance; } // Constructor @@ -36,536 +36,533 @@ conline test_class_output() void conline::set_defaults() { - // Expand tau if necessary - if (tau.n_elem == 1) - { - tau.resize(P); - tau.fill(tau(0)); - } + timer.autoreturn = get_timings; - // Initial params - w0.ones(D, P, K); - w0 /= K; + // Expand tau if necessary + if (tau.n_elem == 1) + { + tau.resize(P); + tau.fill(tau(0)); + } - R0.zeros(D, P, K); + // Initial params + w0.ones(D, P, K); + w0 /= K; - predictions_got_sorted.zeros(T + T_E_Y, D); + R0.zeros(D, P, K); + + predictions_got_sorted.zeros(T + T_E_Y, D); } void conline::set_grid_objects() { - - if (get_timings) - timer.tic("init"); - - opt_index.zeros(T + 1); - if (save_past_performance) - { - past_performance.set_size(T); - } - else + timer.tic("init"); + + opt_index.zeros(T + 1); + if (save_past_performance) + { + past_performance.set_size(T); + } + else + { + past_performance.set_size(1); + } + tmp_performance.zeros(D, P); + cum_performance.zeros(X); + + predictions.zeros(T + T_E_Y, D, P); + if (save_predictions_grid) + { + predictions_grid.set_size(T + T_E_Y); + } + else + { + predictions_grid.set_size(1 + lead_time); + } + weights_tmp.set_size(X); + weights.set_size(T + 1); + + loss_for.zeros(T, D, P); + loss_exp.set_size(T); + + V.set_size(X); + E.set_size(X); + eta.set_size(X); + R.set_size(X); + beta.set_size(X); + beta0field.set_size(X); + + for (unsigned int x = 0; x < X; x++) + { + unsigned int Pr = basis_pr(params["basis_pr_idx"](x) - 1).n_cols; + unsigned int Dr = basis_mv(params["basis_mv_idx"](x) - 1).n_cols; + + // Learning parameters + V(x).zeros(Dr, Pr, K); + E(x).zeros(Dr, Pr, K); + + arma::cube eta_(Dr, Pr, K, fill::zeros); + eta(x) = eta_; + if (method == "ml_poly") { - past_performance.set_size(1); + eta_.fill(exp(350)); + eta(x) = eta_; } - tmp_performance.zeros(D, P); - cum_performance.zeros(X); - predictions.zeros(T + T_E_Y, D, P); - if (save_predictions_grid) + R(x).set_size(Dr, Pr, K); + beta(x).set_size(Dr, Pr, K); + weights_tmp(x).set_size(D, P, K); + + for (unsigned int d = 0; d < D; d++) { - predictions_grid.set_size(T + T_E_Y); + weights_tmp(x).row(d) = w0.row(d); } - else + + for (unsigned int k = 0; k < K; k++) { - predictions_grid.set_size(1 + lead_time); + R(x).slice(k) = basis_mv(params["basis_mv_idx"](x) - 1).t() * + R0.slice(k) * + basis_pr(params["basis_pr_idx"](x) - 1); + R(x).slice(k) = basis_mv(params["basis_mv_idx"](x) - 1).t() * + R0.slice(k) * + basis_pr(params["basis_pr_idx"](x) - 1); + beta(x).slice(k) = pinv( + mat(basis_mv(params["basis_mv_idx"](x) - 1))) * + w0.slice(k) * + pinv(mat(basis_pr(params["basis_pr_idx"](x) - 1))).t(); } - weights_tmp.set_size(X); - weights.set_size(T + 1); + beta0field(x) = beta(x); + } - loss_for.zeros(T, D, P); - loss_exp.set_size(T); + // Predictions at t < lead_time using initial weights + for (unsigned int t = 0; t < lead_time; t++) + { - V.set_size(X); - E.set_size(X); - eta.set_size(X); - R.set_size(X); - beta.set_size(X); - beta0field.set_size(X); + weights(t).set_size(D, P, K); - for (unsigned int x = 0; x < X; x++) + if (save_past_performance) { - unsigned int Pr = basis_pr(params["basis_pr_idx"](x) - 1).n_cols; - unsigned int Dr = basis_mv(params["basis_mv_idx"](x) - 1).n_cols; + past_performance(t).set_size(D, P, X); + past_performance(t).fill(datum::nan); + } + predictions_grid(t).set_size(D, P, X); - // Learning parameters - V(x).zeros(Dr, Pr, K); - E(x).zeros(Dr, Pr, K); + // Store predictions w.r.t. grid for time t + cube tmp_preds_cube(D, P, X); - arma::cube eta_(Dr, Pr, K, fill::zeros); - eta(x) = eta_; - if (method == "ml_poly") - { - eta_.fill(exp(350)); - eta(x) = eta_; - } + for (unsigned int d = 0; d < D; d++) + { + // Save final weights weights_tmp + weights(t).row(d) = weights_tmp(opt_index(t)).row(d); - R(x).set_size(Dr, Pr, K); - beta(x).set_size(Dr, Pr, K); - weights_tmp(x).set_size(D, P, K); + // Store expert predictions temporarily + mat experts_mat = experts(t).row(d); - for (unsigned int d = 0; d < D; d++) - { - weights_tmp(x).row(d) = w0.row(d); - } + for (unsigned int x = 0; x < X; x++) + { - for (unsigned int k = 0; k < K; k++) - { - R(x).slice(k) = basis_mv(params["basis_mv_idx"](x) - 1).t() * - R0.slice(k) * - basis_pr(params["basis_pr_idx"](x) - 1); - R(x).slice(k) = basis_mv(params["basis_mv_idx"](x) - 1).t() * - R0.slice(k) * - basis_pr(params["basis_pr_idx"](x) - 1); - beta(x).slice(k) = pinv( - mat(basis_mv(params["basis_mv_idx"](x) - 1))) * - w0.slice(k) * - pinv(mat(basis_pr(params["basis_pr_idx"](x) - 1))).t(); - } - beta0field(x) = beta(x); - } + mat weights_temp = weights_tmp(x).row(d); - // Predictions at t < lead_time using initial weights - for (unsigned int t = 0; t < lead_time; t++) - { + // Forecasters prediction + vec tmp_preds_vec = sum(weights_temp % experts_mat, 1); - weights(t).set_size(D, P, K); + tmp_preds_cube(span(d), span::all, span(x)) = tmp_preds_vec; + } + } + predictions_grid(t) = tmp_preds_cube; + // Final Prediction + predictions.row(t) = tmp_preds_cube.slice(opt_index(t)); + } - if (save_past_performance) - { - past_performance(t).set_size(D, P, X); - past_performance(t).fill(datum::nan); - } - predictions_grid(t).set_size(D, P, X); + start = lead_time; - // Store predictions w.r.t. grid for time t - cube tmp_preds_cube(D, P, X); + timer.toc("init"); +} - for (unsigned int d = 0; d < D; d++) - { - // Save final weights weights_tmp - weights(t).row(d) = weights_tmp(opt_index(t)).row(d); +void conline::learn() +{ + Progress prog(T, trace); - // Store expert predictions temporarily - mat experts_mat = experts(t).row(d); + timer.tic("core"); - for (unsigned int x = 0; x < X; x++) - { + for (unsigned int tp = 0; tp < predictions_grid.n_rows; tp++) + { + predictions_grid(tp).set_size(D, P, X); + } - mat weights_temp = weights_tmp(x).row(d); + for (unsigned int t = start; t < T; t++) + { - // Forecasters prediction - vec tmp_preds_vec = sum(weights_temp % experts_mat, 1); + weights(t).set_size(D, P, K); - tmp_preds_cube(span(d), span::all, span(x)) = tmp_preds_vec; - } - } - predictions_grid(t) = tmp_preds_cube; - // Final Prediction - predictions.row(t) = tmp_preds_cube.slice(opt_index(t)); - } + if (save_past_performance) + past_performance(t).set_size(D, P, X); - start = lead_time; - if (get_timings) - timer.toc("init"); -} + timer.tic("loss"); -void conline::learn() -{ - Progress prog(T, trace); - if (get_timings) - timer.tic("core"); + // Store predictions w.r.t. grid for time t + cube tmp_preds_cube(D, P, X); - for (unsigned int tp = 0; tp < predictions_grid.n_rows; tp++) + for (unsigned int d = 0; d < D; d++) { - predictions_grid(tp).set_size(D, P, X); - } + // Save final weights weights_tmp + weights(t).row(d) = weights_tmp(opt_index(t)).row(d); - for (unsigned int t = start; t < T; t++) - { + // Store expert predictions temporarily + mat experts_mat = experts(t).row(d); - weights(t).set_size(D, P, K); + // Predictions using different parameter values + for (unsigned int x = 0; x < X; x++) + { - if (save_past_performance) - past_performance(t).set_size(D, P, X); + mat weights_temp = weights_tmp(x).row(d); + vec tmp_preds_vec = sum(weights_temp % experts_mat, 1); - if (get_timings) - timer.tic("loss"); + // Sort predictions if quantile_crossing is prohibited + if (!allow_quantile_crossing) + { + if ((x == opt_index(t)) && (!tmp_preds_vec.is_sorted())) + { + predictions_got_sorted(t, d) = 1; + } + tmp_preds_vec = arma::sort(tmp_preds_vec, "ascend", 0); + } + tmp_preds_cube(span(d), span::all, span(x)) = tmp_preds_vec; + } + } - // Store predictions w.r.t. grid for time t - cube tmp_preds_cube(D, P, X); + // Define which observation to base the weight update on. That is the + // most recent prediction unless lead_time is greater than 0. + arma::uword predictions_tmp_idx; + if (save_predictions_grid) + { + predictions_tmp_idx = t - lead_time; + predictions_grid(t) = tmp_preds_cube; + } + else + { + // If save_predictions_grid is false, the first observation is + // allways the one we base the weight update on. This is because + // the size of predictions_grid is only 1 + lead_time. + predictions_tmp_idx = 0; + for (unsigned int tp = 0; tp < predictions_grid.n_rows - 1; tp++) + { + predictions_grid(tp) = predictions_grid(tp + 1); + } + predictions_grid(predictions_grid.n_rows - 1) = tmp_preds_cube; + } - for (unsigned int d = 0; d < D; d++) - { - // Save final weights weights_tmp - weights(t).row(d) = weights_tmp(opt_index(t)).row(d); + // Final prediction + predictions.row(t) = tmp_preds_cube.slice(opt_index(t)); - // Store expert predictions temporarily - mat experts_mat = experts(t).row(d); + timer.toc("loss"); + for (unsigned int x = 0; x < X; x++) + { - // Predictions using different parameter values - for (unsigned int x = 0; x < X; x++) + timer.tic("regret"); + mat lexp_int(P, K); // Experts loss + mat lexp_ext(P, K); // Experts loss + mat lexp(P, K); // Experts loss + vec lfor(P); // Forecasters loss + cube regret_tmp(D, P, K); + cube regret(basis_mv(params["basis_mv_idx"](x) - 1).n_cols, + basis_pr(params["basis_pr_idx"](x) - 1).n_cols, + K); // Dr x Pr x K + + for (unsigned int d = 0; d < D; d++) + { +#pragma omp parallel for + for (unsigned int p = 0; p < P; p++) + { + tmp_performance(d, p) = loss(y(t, d), + predictions_grid(predictions_tmp_idx)(d, p, x), + 9999, // where evaluate loss_gradient + loss_function, // method + tau(p), // tau + loss_parameter, // alpha + false); + + if (params["loss_share"](x) != 1) + { + for (unsigned int k = 0; k < K; k++) { + lexp_int(p, k) = loss(y(t, d), + experts(t)(d, p, k), + predictions_grid(predictions_tmp_idx)(d, p, x), // where evaluate loss_gradient + loss_function, // method + tau(p), // tau + loss_parameter, // alpha + loss_gradient); + } - mat weights_temp = weights_tmp(x).row(d); - vec tmp_preds_vec = sum(weights_temp % experts_mat, 1); - - // Sort predictions if quantile_crossing is prohibited - if (!allow_quantile_crossing) - { - if ((x == opt_index(t)) && (!tmp_preds_vec.is_sorted())) - { - predictions_got_sorted(t, d) = 1; - } - tmp_preds_vec = arma::sort(tmp_preds_vec, "ascend", 0); - } - tmp_preds_cube(span(d), span::all, span(x)) = tmp_preds_vec; + if (params["loss_share"](x) == 0) + { + lexp.row(p) = lexp_int.row(p); } + else + { + lexp_ext.row(p) = arma::vectorise(loss_array(t).tube(d, p)).t(); + lexp.row(p) = (1 - params["loss_share"](x)) * lexp_int.row(p) + params["loss_share"](x) * lexp_ext.row(p); + } + } + else + { + lexp_ext.row(p) = arma::vectorise(loss_array(t).tube(d, p)).t(); + lexp.row(p) = lexp_ext.row(p); + } + lfor(p) = loss(y(t, d), + predictions_grid(predictions_tmp_idx)(d, p, x), + predictions_grid(predictions_tmp_idx)(d, p, x), // where to evaluate loss_gradient + loss_function, // method + tau(p), // tau + loss_parameter, // alpha + loss_gradient); } - // Define which observation to base the weight update on. That is the - // most recent prediction unless lead_time is greater than 0. - arma::uword predictions_tmp_idx; - if (save_predictions_grid) + mat regret_int(P, K); + mat regret_ext(P, K); + + if (params["regret_share"](x) != 1) { - predictions_tmp_idx = t - lead_time; - predictions_grid(t) = tmp_preds_cube; + regret_int = (lfor - lexp_int.each_col()).t(); + + if (params["regret_share"](x) == 0) + { + regret_tmp.row(d) = regret_int.t(); + } + else + { + regret_ext = regret_array(t).row(d); + regret_ext = regret_ext.t(); + regret_tmp.row(d) = ((1 - params["regret_share"](x)) * regret_int + params["regret_share"](x) * regret_ext).t(); + } } else { - // If save_predictions_grid is false, the first observation is - // allways the one we base the weight update on. This is because - // the size of predictions_grid is only 1 + lead_time. - predictions_tmp_idx = 0; - for (unsigned int tp = 0; tp < predictions_grid.n_rows - 1; tp++) - { - predictions_grid(tp) = predictions_grid(tp + 1); - } - predictions_grid(predictions_grid.n_rows - 1) = tmp_preds_cube; + regret_ext = regret_array(t).row(d); + regret_ext = regret_ext.t(); + regret_tmp.row(d) = regret_ext.t(); } + } - // Final prediction - predictions.row(t) = tmp_preds_cube.slice(opt_index(t)); +#pragma omp parallel for + for (unsigned int k = 0; k < K; k++) + { + regret.slice(k) = basis_mv(params["basis_mv_idx"](x) - 1).t() * regret_tmp.slice(k) * basis_pr(params["basis_pr_idx"](x) - 1); + regret.slice(k) *= double(basis_pr(params["basis_pr_idx"](x) - 1).n_cols) / double(P); + regret.slice(k) *= double(basis_mv(params["basis_mv_idx"](x) - 1).n_cols) / double(D); + } + + timer.toc("regret"); - if (get_timings) - timer.toc("loss"); - for (unsigned int x = 0; x < X; x++) + timer.tic("learning"); +#pragma omp parallel for collapse(2) + for (unsigned int dr = 0; dr < regret.n_rows; dr++) + { + for (unsigned int pr = 0; pr < regret.n_cols; pr++) { - if (get_timings) - timer.tic("regret"); - mat lexp_int(P, K); // Experts loss - mat lexp_ext(P, K); // Experts loss - mat lexp(P, K); // Experts loss - vec lfor(P); // Forecasters loss - cube regret_tmp(D, P, K); - cube regret(basis_mv(params["basis_mv_idx"](x) - 1).n_cols, - basis_pr(params["basis_pr_idx"](x) - 1).n_cols, - K); // Dr x Pr x K - - for (unsigned int d = 0; d < D; d++) + + vec r = regret.tube(dr, pr); + + if (method == "ewa") + { + // Update the cumulative regret used by eta + R(x).tube(dr, pr) = vectorise(R(x).tube(dr, pr) * (1 - params["forget_regret"](x))) + r; + eta(x).tube(dr, pr).fill(params["gamma"](x)); + beta(x).tube(dr, pr) = vectorise(beta0field(x).tube(dr, pr)).t() * K % softmax_r(params["gamma"](x) * vectorise(R(x).tube(dr, pr)).t()); + } + else if (method == "ml_poly") + { + // Update the cumulative regret used by ML_Poly + R(x) + .tube(dr, pr) = vectorise(R(x).tube(dr, pr) * (1 - params["forget_regret"](x))) + r; + + // Update the learning rate + eta(x).tube(dr, pr) = 1 / (1 / vectorise(eta(x).tube(dr, pr)).t() + square(r.t())); + + beta(x).tube(dr, pr) = vectorise(beta0field(x).tube(dr, pr)).t() * K * params["gamma"](x) % vectorise(eta(x).tube(dr, pr)).t() % pmax_arma(vectorise(R(x).tube(dr, pr)).t(), exp(-700)); + beta(x).tube(dr, pr) /= accu(beta(x).tube(dr, pr)); + } + else if (method == "boa" || method == "bewa") + { + V(x).tube(dr, pr) = vectorise(V(x).tube(dr, pr)).t() * (1 - params["forget_regret"](x)) + square(r.t()); + + E(x).tube(dr, pr) = pmax_arma(max(vectorise(E(x).tube(dr, pr)).t() * (1 - params["forget_regret"](x)), abs(r.t())), exp(-350)); + + eta(x) + .tube(dr, pr) = + pmin_arma( + min(1 / (2 * vectorise(E(x).tube(dr, pr))), + sqrt(-log(vectorise(beta0field(x).tube(dr, pr))) / pmax_arma(vectorise(V(x).tube(dr, pr)), exp(-350)))), + exp(350)); + + vec r_reg = r - vectorise(eta(x).tube(dr, pr)) % square(r); + + R(x).tube(dr, pr) *= (1 - params["forget_regret"](x)); // forget + R(x).tube(dr, pr) += + 0.5 * (r_reg + conv_to::from(vectorise(eta(x).tube(dr, pr)) % r > 0.5) % (2 * vectorise(E(x).tube(dr, pr)))); + + if (method == "boa") { -#pragma omp parallel for - for (unsigned int p = 0; p < P; p++) - { - tmp_performance(d, p) = loss(y(t, d), - predictions_grid(predictions_tmp_idx)(d, p, x), - 9999, // where evaluate loss_gradient - loss_function, // method - tau(p), // tau - loss_parameter, // alpha - false); - - if (params["loss_share"](x) != 1) - { - for (unsigned int k = 0; k < K; k++) - { - lexp_int(p, k) = loss(y(t, d), - experts(t)(d, p, k), - predictions_grid(predictions_tmp_idx)(d, p, x), // where evaluate loss_gradient - loss_function, // method - tau(p), // tau - loss_parameter, // alpha - loss_gradient); - } - - if (params["loss_share"](x) == 0) - { - lexp.row(p) = lexp_int.row(p); - } - else - { - lexp_ext.row(p) = arma::vectorise(loss_array(t).tube(d, p)).t(); - lexp.row(p) = (1 - params["loss_share"](x)) * lexp_int.row(p) + params["loss_share"](x) * lexp_ext.row(p); - } - } - else - { - lexp_ext.row(p) = arma::vectorise(loss_array(t).tube(d, p)).t(); - lexp.row(p) = lexp_ext.row(p); - } - lfor(p) = loss(y(t, d), - predictions_grid(predictions_tmp_idx)(d, p, x), - predictions_grid(predictions_tmp_idx)(d, p, x), // where to evaluate loss_gradient - loss_function, // method - tau(p), // tau - loss_parameter, // alpha - loss_gradient); - } - - mat regret_int(P, K); - mat regret_ext(P, K); - - if (params["regret_share"](x) != 1) - { - regret_int = (lfor - lexp_int.each_col()).t(); - - if (params["regret_share"](x) == 0) - { - regret_tmp.row(d) = regret_int.t(); - } - else - { - regret_ext = regret_array(t).row(d); - regret_ext = regret_ext.t(); - regret_tmp.row(d) = ((1 - params["regret_share"](x)) * regret_int + params["regret_share"](x) * regret_ext).t(); - } - } - else - { - regret_ext = regret_array(t).row(d); - regret_ext = regret_ext.t(); - regret_tmp.row(d) = regret_ext.t(); - } + // Wintenberger + beta(x).tube(dr, pr) = vectorise(beta0field(x).tube(dr, pr)).t() * K % softmax_r(log(params["gamma"](x) * vectorise(eta(x).tube(dr, pr)).t()) + params["gamma"](x) * vectorise(eta(x).tube(dr, pr)).t() % vectorise(R(x).tube(dr, pr)).t()); } - -#pragma omp parallel for - for (unsigned int k = 0; k < K; k++) + else { - regret.slice(k) = basis_mv(params["basis_mv_idx"](x) - 1).t() * regret_tmp.slice(k) * basis_pr(params["basis_pr_idx"](x) - 1); - regret.slice(k) *= double(basis_pr(params["basis_pr_idx"](x) - 1).n_cols) / double(P); - regret.slice(k) *= double(basis_mv(params["basis_mv_idx"](x) - 1).n_cols) / double(D); + // Gaillard + beta(x).tube(dr, pr) = vectorise(beta0field(x).tube(dr, pr)).t() * K % softmax_r(params["gamma"](x) * vectorise(eta(x).tube(dr, pr)).t() % vectorise(R(x).tube(dr, pr)).t()); } - - if (get_timings) - timer.toc("regret"); - if (get_timings) - timer.tic("learning"); -#pragma omp parallel for collapse(2) - for (unsigned int dr = 0; dr < regret.n_rows; dr++) + } + else + { + Rcpp::stop("Choose 'boa', 'bewa', 'ml_poly' or 'ewa' as method."); + } + + // // Apply thresholds + if (params["soft_threshold"](x) > 0) + { + int best_k = beta(x).tube(dr, pr).index_max(); + + for (double &e : beta(x).tube(dr, pr)) { - for (unsigned int pr = 0; pr < regret.n_cols; pr++) - { - - vec r = regret.tube(dr, pr); - - if (method == "ewa") - { - // Update the cumulative regret used by eta - R(x).tube(dr, pr) = vectorise(R(x).tube(dr, pr) * (1 - params["forget_regret"](x))) + r; - eta(x).tube(dr, pr).fill(params["gamma"](x)); - beta(x).tube(dr, pr) = vectorise(beta0field(x).tube(dr, pr)).t() * K % softmax_r(params["gamma"](x) * vectorise(R(x).tube(dr, pr)).t()); - } - else if (method == "ml_poly") - { - // Update the cumulative regret used by ML_Poly - R(x) - .tube(dr, pr) = vectorise(R(x).tube(dr, pr) * (1 - params["forget_regret"](x))) + r; - - // Update the learning rate - eta(x).tube(dr, pr) = 1 / (1 / vectorise(eta(x).tube(dr, pr)).t() + square(r.t())); - - beta(x).tube(dr, pr) = vectorise(beta0field(x).tube(dr, pr)).t() * K * params["gamma"](x) % vectorise(eta(x).tube(dr, pr)).t() % pmax_arma(vectorise(R(x).tube(dr, pr)).t(), exp(-700)); - beta(x).tube(dr, pr) /= accu(beta(x).tube(dr, pr)); - } - else if (method == "boa" || method == "bewa") - { - V(x).tube(dr, pr) = vectorise(V(x).tube(dr, pr)).t() * (1 - params["forget_regret"](x)) + square(r.t()); - - E(x).tube(dr, pr) = pmax_arma(max(vectorise(E(x).tube(dr, pr)).t() * (1 - params["forget_regret"](x)), abs(r.t())), exp(-350)); - - eta(x) - .tube(dr, pr) = - pmin_arma( - min(1 / (2 * vectorise(E(x).tube(dr, pr))), - sqrt(-log(vectorise(beta0field(x).tube(dr, pr))) / pmax_arma(vectorise(V(x).tube(dr, pr)), exp(-350)))), - exp(350)); - - vec r_reg = r - vectorise(eta(x).tube(dr, pr)) % square(r); - - R(x).tube(dr, pr) *= (1 - params["forget_regret"](x)); // forget - R(x).tube(dr, pr) += - 0.5 * (r_reg + conv_to::from(vectorise(eta(x).tube(dr, pr)) % r > 0.5) % (2 * vectorise(E(x).tube(dr, pr)))); - - if (method == "boa") - { - // Wintenberger - beta(x).tube(dr, pr) = vectorise(beta0field(x).tube(dr, pr)).t() * K % softmax_r(log(params["gamma"](x) * vectorise(eta(x).tube(dr, pr)).t()) + params["gamma"](x) * vectorise(eta(x).tube(dr, pr)).t() % vectorise(R(x).tube(dr, pr)).t()); - } - else - { - // Gaillard - beta(x).tube(dr, pr) = vectorise(beta0field(x).tube(dr, pr)).t() * K % softmax_r(params["gamma"](x) * vectorise(eta(x).tube(dr, pr)).t() % vectorise(R(x).tube(dr, pr)).t()); - } - } - else - { - Rcpp::stop("Choose 'boa', 'bewa', 'ml_poly' or 'ewa' as method."); - } - - // // Apply thresholds - if (params["soft_threshold"](x) > 0) - { - int best_k = beta(x).tube(dr, pr).index_max(); - - for (double &e : beta(x).tube(dr, pr)) - { - threshold_soft(e, params["soft_threshold"](x)); - } - if (accu(beta(x).tube(dr, pr)) == 0) - { - beta(x)(dr, pr, best_k) = 1; - } - } - - if (params["hard_threshold"](x) > 0) - { - int best_k = beta(x).tube(dr, pr).index_max(); - for (double &e : beta(x).tube(dr, pr)) - { - threshold_hard(e, params["hard_threshold"](x)); - } - if (accu(beta(x).tube(dr, pr)) == 0) - { - beta(x)(dr, pr, best_k) = 1; - } - } - - // Add fixed_share - beta(x).tube(dr, pr) = - (1 - params["fixed_share"](x)) * vectorise(beta(x).tube(dr, pr)) + - (params["fixed_share"](x) / K); - } // pr - } // dr - if (get_timings) - timer.toc("learning"); - -#pragma omp parallel for - // Smoothing - for (unsigned int k = 0; k < K; k++) + threshold_soft(e, params["soft_threshold"](x)); + } + if (accu(beta(x).tube(dr, pr)) == 0) { - if (get_timings) - timer.tic("smoothing"); - weights_tmp(x).slice(k) = hat_mv(params["hat_mv_idx"](x) - 1) * - basis_mv(params["basis_mv_idx"](x) - 1) * - beta(x).slice(k) * - basis_pr(params["basis_pr_idx"](x) - 1).t() * - hat_pr(params["hat_pr_idx"](x) - 1); - if (get_timings) - timer.toc("smoothing"); + beta(x)(dr, pr, best_k) = 1; } + } -#pragma omp parallel for - // Enshure that constraints hold - for (unsigned int p = 0; p < P; p++) + if (params["hard_threshold"](x) > 0) + { + int best_k = beta(x).tube(dr, pr).index_max(); + for (double &e : beta(x).tube(dr, pr)) + { + threshold_hard(e, params["hard_threshold"](x)); + } + if (accu(beta(x).tube(dr, pr)) == 0) { - for (unsigned int d = 0; d < D; d++) - { - // Positivity - weights_tmp(x)(span(d), span(p), span::all) = - pmax_arma(weights_tmp(x)(span(d), span(p), span::all), exp(-700)); - - // // Affinity - weights_tmp(x)(span(d), span(p), span::all) /= - accu(weights_tmp(x)(span(d), span(p), span::all)); - } + beta(x)(dr, pr, best_k) = 1; } - if (save_past_performance) - past_performance(t).slice(x) = tmp_performance; - R_CheckUserInterrupt(); + } - // Apply forget - cum_performance(x) *= (1 - forget_past_performance); - // Add new loss - cum_performance(x) += accu(tmp_performance) / double(D * P); + // Add fixed_share + beta(x).tube(dr, pr) = + (1 - params["fixed_share"](x)) * vectorise(beta(x).tube(dr, pr)) + + (params["fixed_share"](x) / K); + } // pr + } // dr - } // x + timer.toc("learning"); - opt_index(t + 1) = cum_performance.index_min(); - prog.increment(); // Update progress +#pragma omp parallel for + // Smoothing + for (unsigned int k = 0; k < K; k++) + { - } // t + timer.tic("smoothing"); + weights_tmp(x).slice(k) = hat_mv(params["hat_mv_idx"](x) - 1) * + basis_mv(params["basis_mv_idx"](x) - 1) * + beta(x).slice(k) * + basis_pr(params["basis_pr_idx"](x) - 1).t() * + hat_pr(params["hat_pr_idx"](x) - 1); - // Save Final Weights and Prediction - weights(T) = weights_tmp(opt_index(T)); + timer.toc("smoothing"); + } - // Predict residual expert forecasts if any are available - for (unsigned int t = T; t < T + T_E_Y; t++) - { +#pragma omp parallel for + // Enshure that constraints hold + for (unsigned int p = 0; p < P; p++) + { for (unsigned int d = 0; d < D; d++) { - mat experts_mat = experts(t).row(d); - mat weights_temp = weights(T).row(d); - vec tmp_preds_vec = sum(weights_temp % experts_mat, 1); + // Positivity + weights_tmp(x)(span(d), span(p), span::all) = + pmax_arma(weights_tmp(x)(span(d), span(p), span::all), exp(-700)); - // Sort predictions if quantile_crossing is prohibited - if (!allow_quantile_crossing) - { - if (!tmp_preds_vec.is_sorted()) - { - predictions_got_sorted(t, d) = 1; - } - tmp_preds_vec = arma::sort(tmp_preds_vec, "ascend", 0); - } - predictions.tube(t, d) = tmp_preds_vec; + // // Affinity + weights_tmp(x)(span(d), span(p), span::all) /= + accu(weights_tmp(x)(span(d), span(p), span::all)); } + } + if (save_past_performance) + past_performance(t).slice(x) = tmp_performance; + R_CheckUserInterrupt(); + + // Apply forget + cum_performance(x) *= (1 - forget_past_performance); + // Add new loss + cum_performance(x) += accu(tmp_performance) / double(D * P); + + } // x + + opt_index(t + 1) = cum_performance.index_min(); + prog.increment(); // Update progress + + } // t + + // Save Final Weights and Prediction + weights(T) = weights_tmp(opt_index(T)); + + // Predict residual expert forecasts if any are available + for (unsigned int t = T; t < T + T_E_Y; t++) + { + for (unsigned int d = 0; d < D; d++) + { + mat experts_mat = experts(t).row(d); + mat weights_temp = weights(T).row(d); + vec tmp_preds_vec = sum(weights_temp % experts_mat, 1); + + // Sort predictions if quantile_crossing is prohibited + if (!allow_quantile_crossing) + { + if (!tmp_preds_vec.is_sorted()) + { + predictions_got_sorted(t, d) = 1; + } + tmp_preds_vec = arma::sort(tmp_preds_vec, "ascend", 0); + } + predictions.tube(t, d) = tmp_preds_vec; } + } + + // Save losses suffered by forecaster and experts - // Save losses suffered by forecaster and experts - if (get_timings) - timer.tic("loss_for_exp"); + timer.tic("loss_for_exp"); #pragma omp parallel for - for (unsigned int t = 0; t < T; t++) - { - loss_exp(t).set_size(D, P, K); + for (unsigned int t = 0; t < T; t++) + { + loss_exp(t).set_size(D, P, K); - for (unsigned int d = 0; d < D; d++) + for (unsigned int d = 0; d < D; d++) + { + for (unsigned int p = 0; p < P; p++) + { + for (unsigned int k = 0; k < K; k++) { - for (unsigned int p = 0; p < P; p++) - { - for (unsigned int k = 0; k < K; k++) - { - loss_exp(t)(d, p, k) = - loss(y(t, d), - experts(t)(d, p, k), - 9999, // where to evaluate the loss_gradient - loss_function, // method - tau(p), // tau - loss_parameter, // alpha - false); // loss_gradient - } - loss_for(t, d, p) = loss(y(t, d), - predictions(t, d, p), - 9999, // where to evaluate the loss_gradient - loss_function, // method - tau(p), // tau - loss_parameter, // alpha - false); // loss_gradient; - } + loss_exp(t)(d, p, k) = + loss(y(t, d), + experts(t)(d, p, k), + 9999, // where to evaluate the loss_gradient + loss_function, // method + tau(p), // tau + loss_parameter, // alpha + false); // loss_gradient } + loss_for(t, d, p) = loss(y(t, d), + predictions(t, d, p), + 9999, // where to evaluate the loss_gradient + loss_function, // method + tau(p), // tau + loss_parameter, // alpha + false); // loss_gradient; + } } - if (get_timings) - timer.toc("loss_for_exp"); - if (get_timings) - timer.toc("core"); + } + + timer.toc("loss_for_exp"); + + timer.toc("core"); } void conline::init_update( @@ -573,115 +570,114 @@ void conline::init_update( arma::mat &new_y, arma::field &new_experts) { - if (get_timings) - timer.tic("init update"); - - // This creates a references not copies - Rcpp::List specification = object["specification"]; - Rcpp::List model_parameters = specification["parameters"]; - Rcpp::List model_data = specification["data"]; - Rcpp::List model_objects = specification["objects"]; - - // Data - - // Join old and new expert_predictions - arma::field old_experts = model_data["experts"]; - experts.set_size(old_experts.n_rows + new_experts.n_rows); - experts.rows(0, old_experts.n_rows - 1) = old_experts; - if (new_experts.n_rows > 0) - { - experts.rows(old_experts.n_rows, experts.n_rows - 1) = new_experts; - } - - y = Rcpp::as(model_data["y"]); - y.insert_rows(y.n_rows, new_y); - start = T - new_y.n_rows; - - if (T_E_Y < 0) - Rcpp::stop("Number of provided expert predictions has to match or exceed observations."); - - tau = Rcpp::as(model_data["tau"]); - - params = Rcpp::as>(object["parametergrid"]); - params_basis_pr = Rcpp::as>(object["params_basis_pr"]); - params_basis_mv = Rcpp::as>(object["params_basis_mv"]); - params_hat_pr = Rcpp::as>(object["params_hat_pr"]); - params_hat_mv = Rcpp::as>(object["params_hat_mv"]); - - opt_index = Rcpp::as(object["opt_index"]); - // Zero indexing in C++ - opt_index -= 1; - opt_index.resize(T + 1); - - tmp_performance.zeros(D, P); - cum_performance = Rcpp::as(model_objects["cum_performance"]); - - weights_tmp = - Rcpp::as>(model_objects["weights_tmp"]); - - // // Output Objects - predictions = Rcpp::as(object["predictions"]); - predictions.resize(T + T_E_Y, D, P); - predictions_got_sorted = Rcpp::as(object["predictions_got_sorted"]); - predictions_got_sorted.resize(T + T_E_Y, D); - weights.set_size(T + 1); - weights.rows(0, start) = Rcpp::as>(object["weights"]); - - basis_pr = Rcpp::as>(model_objects["basis_pr"]); - basis_mv = Rcpp::as>(model_objects["basis_mv"]); - hat_pr = Rcpp::as>(model_objects["hat_pr"]); - hat_mv = Rcpp::as>(model_objects["hat_mv"]); - - V = Rcpp::as>(model_objects["V"]); - E = Rcpp::as>(model_objects["E"]); - eta = Rcpp::as>(model_objects["eta"]); - R = Rcpp::as>(model_objects["R"]); - beta = Rcpp::as>(model_objects["beta"]); - beta0field = Rcpp::as>(model_objects["beta0field"]); - - // // // Misc parameters - lead_time = model_parameters["lead_time"]; - loss_function = Rcpp::as(model_parameters["loss_function"]); - loss_parameter = model_parameters["loss_parameter"]; - loss_gradient = model_parameters["loss_gradient"]; - method = Rcpp::as(model_parameters["method"]); - - forget_past_performance = model_parameters["forget_past_performance"]; - allow_quantile_crossing = model_parameters["allow_quantile_crossing"]; - - save_past_performance = model_parameters["save_past_performance"]; - save_predictions_grid = model_parameters["save_predictions_grid"]; - - if (save_past_performance) - { - past_performance.set_size(T); - past_performance.rows(0, start - 1) = - Rcpp::as>(object["past_performance"]); - } - else - { - past_performance.set_size(1); - } - - if (save_predictions_grid) - { - predictions_grid.set_size(T + T_E_Y); - predictions_grid.rows(0, old_experts.n_rows - 1) = Rcpp::as>(model_objects["predictions_grid"]); - } - else - { - predictions_grid.set_size(1 + lead_time); - predictions_grid = Rcpp::as>(model_objects["predictions_grid"]); - } - - loss_for.zeros(T, D, P); - loss_for.rows(0, start - 1) = - Rcpp::as(object["forecaster_loss"]); - - loss_exp.set_size(T); - loss_exp.rows(0, start - 1) = - Rcpp::as>(object["experts_loss"]); - if (get_timings) - timer.toc("init update"); + timer.tic("init update"); + + // This creates a references not copies + Rcpp::List specification = object["specification"]; + Rcpp::List model_parameters = specification["parameters"]; + Rcpp::List model_data = specification["data"]; + Rcpp::List model_objects = specification["objects"]; + + // Data + + // Join old and new expert_predictions + arma::field old_experts = model_data["experts"]; + experts.set_size(old_experts.n_rows + new_experts.n_rows); + experts.rows(0, old_experts.n_rows - 1) = old_experts; + if (new_experts.n_rows > 0) + { + experts.rows(old_experts.n_rows, experts.n_rows - 1) = new_experts; + } + + y = Rcpp::as(model_data["y"]); + y.insert_rows(y.n_rows, new_y); + start = T - new_y.n_rows; + + if (T_E_Y < 0) + Rcpp::stop("Number of provided expert predictions has to match or exceed observations."); + + tau = Rcpp::as(model_data["tau"]); + + params = Rcpp::as>(object["parametergrid"]); + params_basis_pr = Rcpp::as>(object["params_basis_pr"]); + params_basis_mv = Rcpp::as>(object["params_basis_mv"]); + params_hat_pr = Rcpp::as>(object["params_hat_pr"]); + params_hat_mv = Rcpp::as>(object["params_hat_mv"]); + + opt_index = Rcpp::as(object["opt_index"]); + // Zero indexing in C++ + opt_index -= 1; + opt_index.resize(T + 1); + + tmp_performance.zeros(D, P); + cum_performance = Rcpp::as(model_objects["cum_performance"]); + + weights_tmp = + Rcpp::as>(model_objects["weights_tmp"]); + + // // Output Objects + predictions = Rcpp::as(object["predictions"]); + predictions.resize(T + T_E_Y, D, P); + predictions_got_sorted = Rcpp::as(object["predictions_got_sorted"]); + predictions_got_sorted.resize(T + T_E_Y, D); + weights.set_size(T + 1); + weights.rows(0, start) = Rcpp::as>(object["weights"]); + + basis_pr = Rcpp::as>(model_objects["basis_pr"]); + basis_mv = Rcpp::as>(model_objects["basis_mv"]); + hat_pr = Rcpp::as>(model_objects["hat_pr"]); + hat_mv = Rcpp::as>(model_objects["hat_mv"]); + + V = Rcpp::as>(model_objects["V"]); + E = Rcpp::as>(model_objects["E"]); + eta = Rcpp::as>(model_objects["eta"]); + R = Rcpp::as>(model_objects["R"]); + beta = Rcpp::as>(model_objects["beta"]); + beta0field = Rcpp::as>(model_objects["beta0field"]); + + // // // Misc parameters + lead_time = model_parameters["lead_time"]; + loss_function = Rcpp::as(model_parameters["loss_function"]); + loss_parameter = model_parameters["loss_parameter"]; + loss_gradient = model_parameters["loss_gradient"]; + method = Rcpp::as(model_parameters["method"]); + + forget_past_performance = model_parameters["forget_past_performance"]; + allow_quantile_crossing = model_parameters["allow_quantile_crossing"]; + + save_past_performance = model_parameters["save_past_performance"]; + save_predictions_grid = model_parameters["save_predictions_grid"]; + + if (save_past_performance) + { + past_performance.set_size(T); + past_performance.rows(0, start - 1) = + Rcpp::as>(object["past_performance"]); + } + else + { + past_performance.set_size(1); + } + + if (save_predictions_grid) + { + predictions_grid.set_size(T + T_E_Y); + predictions_grid.rows(0, old_experts.n_rows - 1) = Rcpp::as>(model_objects["predictions_grid"]); + } + else + { + predictions_grid.set_size(1 + lead_time); + predictions_grid = Rcpp::as>(model_objects["predictions_grid"]); + } + + loss_for.zeros(T, D, P); + loss_for.rows(0, start - 1) = + Rcpp::as(object["forecaster_loss"]); + + loss_exp.set_size(T); + loss_exp.rows(0, start - 1) = + Rcpp::as>(object["experts_loss"]); + + timer.toc("init update"); } diff --git a/tests/testthat/test-output.R b/tests/testthat/test-output.R index 541a454..601a027 100644 --- a/tests/testthat/test-output.R +++ b/tests/testthat/test-output.R @@ -16,15 +16,16 @@ y <- rnorm(n = T) # Expert predictions experts <- array(dim = c(T, P, N)) for (t in 1:T) { - experts[t, , 1] <- qnorm(prob_grid, mean = -5, sd = 2) - experts[t, , 2] <- qnorm(prob_grid, mean = 5, sd = 2) + experts[t, , 1] <- qnorm(prob_grid, mean = -5, sd = 2) + experts[t, , 2] <- qnorm(prob_grid, mean = 5, sd = 2) } model <- online( - y = matrix(y), - experts = experts, - tau = prob_grid, - trace = FALSE + y = matrix(y), + experts = experts, + tau = prob_grid, + trace = FALSE, + get_timings = TRUE ) # Return Type @@ -41,4 +42,4 @@ expect_true(all(!is.na(model$predictions))) expect_true(all(!is.na(model$past_perf_wrt_params))) expect_true(all(!is.na(model$opt_index))) expect_true(all(!is.na(model$forecaster_loss))) -expect_true(all(!is.na(model$experts_loss))) \ No newline at end of file +expect_true(all(!is.na(model$experts_loss)))