Refactor (#3)

* refactor in progress
* vignette start + pkgdown

Author: AdrienLeGuillou

.Rbuildignore

@@ -0,0 +1,8 @@
+^renv$
+^renv\.lock$
+^doc$
+^Meta$
+^_pkgdown\.yml$
+^docs$
+^pkgdown$
+^\.github$

.github/.gitignore

@@ -0,0 +1 @@
+*.html

.github/workflows/pkgdown.yaml

@@ -0,0 +1,48 @@
+# Workflow derived from https://github.com/r-lib/actions/tree/v2/examples
+# Need help debugging build failures? Start at https://github.com/r-lib/actions#where-to-find-help
+on:
+  push:
+    branches: [main, master]
+  pull_request:
+    branches: [main, master]
+  release:
+    types: [published]
+  workflow_dispatch:
+
+name: pkgdown
+
+jobs:
+  pkgdown:
+    runs-on: ubuntu-latest
+    # Only restrict concurrency for non-PR jobs
+    concurrency:
+      group: pkgdown-${{ github.event_name != 'pull_request' || github.run_id }}
+    env:
+      GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }}
+    permissions:
+      contents: write
+    steps:
+      - uses: actions/checkout@v3
+
+      - uses: r-lib/actions/setup-pandoc@v2
+
+      - uses: r-lib/actions/setup-r@v2
+        with:
+          use-public-rspm: true
+
+      - uses: r-lib/actions/setup-r-dependencies@v2
+        with:
+          extra-packages: any::pkgdown, local::.
+          needs: website
+
+      - name: Build site
+        run: pkgdown::build_site_github_pages(new_process = FALSE, install = FALSE)
+        shell: Rscript {0}
+
+      - name: Deploy to GitHub pages 🚀
+        if: github.event_name != 'pull_request'
+        uses: JamesIves/[email protected]
+        with:
+          clean: false
+          branch: gh-pages
+          folder: docs

.gitignore

@@ -0,0 +1,12 @@
+inst/doc
+/doc/
+/Meta/
+/renv/
+.Rproj.user
+.Rhistory
+.Rprofile
+.RData
+.Ruserdata
+.DS_Store
+renv/
+docs

DESCRIPTION

@@ -17,9 +17,16 @@ Imports:
    future.apply,
    slurmworkflow,
    fs,
+   lhs,
    dplyr
 Remotes:
     github::EpiModel/slurmworkflow
 Encoding: UTF-8
 Roxygen: list(markdown = TRUE)
 RoxygenNote: 7.2.3
+Suggests:
+    ggplot2,
+    knitr,
+    rmarkdown
+VignetteBuilder: knitr
+URL: https://epimodel.github.io/swfcalib/

NAMESPACE

@@ -3,5 +3,11 @@
 export(calibration_step1)
 export(calibration_step2)
 export(calibration_step3)
+export(determ_end_thresh)
+export(determ_poly_end)
 export(load_sideload)
+export(make_proposer_se_range)
+export(make_shrink_proposer)
+export(render_assessment)
 export(save_sideload)
+importFrom(dplyr,.data)

R/assessment_plots.R

@@ -0,0 +1,89 @@
+make_rmse_plot <- function(job_assess) {
+  d <- job_assess$measures
+  ggplot2::ggplot(d, ggplot2::aes(x = .data$iteration, y = .data$rmse_mean)) +
+    ggplot2::geom_line() +
+    ggplot2::geom_ribbon(
+      ggplot2::aes(ymin = rmse_mean - rmse_sd, ymax = rmse_mean + rmse_sd),
+      alpha = 0.3
+    ) +
+    ggplot2::scale_y_log10() +
+    ggplot2::theme_light() +
+    ggplot2::labs(
+      title = paste0("RMSE Evolution for: ", job_assess$infos$job_id),
+      x = "Iteration",
+      y = "RMSE \n(log10 scale)"
+  )
+}
+
+make_param_volume_plot <- function(job_assess) {
+  d <- job_assess$measures
+  ggplot2::ggplot(d, ggplot2::aes(x = .data$iteration,
+                                  y = .data$param_volume)) +
+    ggplot2::geom_line() +
+    ggplot2::scale_y_log10() +
+    ggplot2::theme_light() +
+    ggplot2::labs(
+      title =
+        paste0("Parameter Space Volume Evolution for: ",
+               job_assess$infos$job_id),
+      x = "Iteration",
+      y = "Parameter Space Volume \n(log10 scale)"
+    )
+}
+
+make_param_spread_plot <- function(job_assess, param) {
+  d <- job_assess$measures
+  d[["y"]] <- d[[paste0("spread__", param)]]
+  ggplot2::ggplot(d, ggplot2::aes(x = .data$iteration, y = .data$y)) +
+    ggplot2::geom_line() +
+    ggplot2::scale_y_log10() +
+    ggplot2::theme_light() +
+    ggplot2::labs(
+      title = paste0("Spread of Parameter: ", param),
+      x = "Iteration",
+      y = "Spread \n(log10 scale)"
+    )
+}
+
+make_target_err_plot <- function(job_assess, target) {
+  d <- job_assess$measures
+  d[["y"]] <- d[[paste0("mean_err__", target)]]
+  d[["ys"]] <- d[[paste0("sd_err__", target)]]
+  ggplot2::ggplot(d, ggplot2::aes(x = .data$iteration, y = .data$y)) +
+    ggplot2::geom_line() +
+    ggplot2::geom_ribbon(
+     ggplot2::aes(ymin = .data$y - .data$ys, ymax = .data$y + .data$ys),
+      alpha = 0.3
+    ) +
+    ggplot2::geom_hline(yintercept = 0, linetype = "dashed") +
+    ggplot2::theme_light() +
+    ggplot2::labs(
+      title = paste0("Mean Error on: ", target),
+      x = "Iteration",
+      y = "Error \n(mean + sd)"
+    )
+}
+
+make_job_plots <- function(job_assess) {
+  out <- list()
+  infos <- job_assess$infos
+  out$rmse <- make_rmse_plot(job_assess)
+  out$volume <- make_param_volume_plot(job_assess)
+  out$params <- lapply(infos$params, make_param_spread_plot, job = job_assess)
+  names(out$params) <- infos$params
+  out$targets <- lapply(infos$targets, make_target_err_plot, job = job_assess)
+  names(out$targets) <- infos$targets
+  out
+}
+
+make_wave_plots <- function(wave_assess) {
+  out <- lapply(wave_assess, make_job_plots)
+  names(out) <- vapply(wave_assess, function(x) x$infos$job_id, character(1))
+  out
+}
+
+make_assessments_plots <- function(assessments) {
+  out <- lapply(assessments, make_wave_plots)
+  names(out) <- names(assessments)
+  out
+}

R/assessment_rmd.R

@@ -0,0 +1,73 @@
+make_wave_rmd <- function(assessments, wave_num) {
+  wave <- assessments[[paste0("wave", wave_num)]]
+  cat("# Wave", wave_num, "\n\n")
+  for (i in seq_along(wave)) {
+    make_job_rmd(wave[[i]])
+  }
+}
+
+make_job_rmd <- function(job_assess) {
+  cat("##", job_assess$infos$job_id, "\n\n")
+
+  cat("### Targets and Parameters", "\n\n")
+  dplyr::tibble(
+    target_name = job_assess$infos$targets,
+    target_value = job_assess$infos$targets_val
+  ) |> knitr::kable(align = "ll") |> print()
+
+  dplyr::tibble(
+    parameter = job_assess$infos$params,
+    initial_range = vapply(
+      job_assess$infos$params_ranges,
+      \(x) paste0(x[1], " - ", x[2]),
+      ""
+    )
+  ) |> knitr::kable(align = "ll") |> print()
+
+  cat("\n\n")
+
+  cat("### Parameter Space and RMSE Evolution", "\n\n")
+
+  make_param_volume_plot(job_assess) |> print()
+  make_rmse_plot(job_assess) |> print()
+  cat("\n\n")
+
+  cat("### Parameter Spreads", "\n\n")
+  for (p in job_assess$infos$params) {
+    make_param_spread_plot(job_assess, p) |> print()
+    cat("\n\n")
+  }
+
+  cat("### Target Errors", "\n\n")
+  for (t in job_assess$infos$targets) {
+    make_target_err_plot(job_assess, t) |> print()
+    cat("\n\n")
+  }
+  cat("\n\n")
+
+}
+
+#' Generate an html report of the auto-calibration
+#'
+#' The report contains descriptions of the parameters spaces and residual errors
+#' over the duration of the calibration.
+#'
+#' @param path_to_assessments Path to an `assessments.rds` file generated by an
+#'   `swfcalib` process.
+#' @param output_filename Name of the html report (default = "assessment.html")
+#' @param output_dir Directory where to store the report (default = current
+#'   working directory)
+#'
+#' @export
+render_assessment <- function(path_to_assessments,
+                              output_filename = "assessment.html",
+                              output_dir = NULL) {
+  if (is.null(output_dir)) output_dir <- getwd()
+  rmarkdown::render(
+    system.file("rmd/assessment.Rmd", package = "swfcalib"),
+    output_file = output_filename,
+    output_dir = output_dir,
+    knit_root_dir = getwd(),
+    params = list(path_to_assessments = path_to_assessments)
+  )
+}

R/assessments.R

@@ -0,0 +1,85 @@
+update_assessments <- function(calib_object, results) {
+  out <- load_assessments(calib_object)
+  if (nrow(results) == 0) {
+    save_assessments(calib_object, out)
+    return(invisible(calib_object))
+  }
+
+  cur_wave <- paste0("wave", get_current_wave(calib_object))
+
+  assessments <- future.apply::future_lapply(
+    get_current_jobs(calib_object),
+    make_job_assessment,
+    calib_object = calib_object,
+    results = results
+  )
+
+  out[[cur_wave]] <- merge_wave_assements(assessments, out[[cur_wave]])
+  save_assessments(calib_object, out)
+  invisible(calib_object)
+}
+
+merge_wave_assements <- function(new, old) {
+  old <- if (is.null(old)) list() else old
+  for (nme in names(new)) {
+    new[[nme]] <- merge_job_assessment(new[[nme]], old[[nme]])
+  }
+  new
+}
+
+merge_job_assessment <- function(new, old) {
+  list(
+    infos = new$infos,
+    measures = dplyr::bind_rows(old$measures, new$measures)
+  )
+}
+
+save_assessments <- function(calib_object, assessments) {
+  saveRDS(assessments, get_assessments_path(calib_object))
+}
+
+get_assessments_path <- function(calib_object) {
+  fs::path(get_root_dir(calib_object), "assessments.rds")
+}
+
+load_assessments <- function(calib_object) {
+  f_path <- get_assessments_path(calib_object)
+  if (fs::file_exists(f_path)) readRDS(f_path) else list()
+}
+
+make_job_assessment <- function(calib_object, job, results) {
+  out <- list()
+
+  out$infos <- job[c("targets", "targets_val", "params")]
+  out$infos$job_id <- get_job_id(job)
+  out$infos$params_ranges <- lapply(job$initial_proposals, range)
+
+  current_iteration <- get_current_iteration(calib_object)
+  current_wave <- get_current_wave(calib_object)
+  d <- dplyr::filter(
+    results,
+    .data$.iteration == current_iteration,
+    .data$.wave == current_wave
+  )
+
+  make_rmse <- function(x, target) sqrt(mean((target - x)^2))
+  iter_rmse <- apply(d[job$targets], 1, make_rmse, target = job$targets_val)
+
+  get_spread <- function(x) diff(range(x))
+  spreads <- vapply(d[job$params], get_spread, numeric(1))
+
+  out$measures <- dplyr::tibble(
+    iteration = current_iteration,
+    rmse_mean = mean(iter_rmse),
+    rmse_sd = sd(iter_rmse),
+    param_volume = prod(spreads)
+  )
+
+  errors <- Map(function(x, target) target - x, d[job$targets], job$targets_val)
+
+  out$measures[paste0("spread__", names(spreads))] <- as.list(spreads)
+  out$measures[paste0("mean_err__", job$targets)] <- lapply(errors, mean)
+  out$measures[paste0("sd_err__", job$targets)] <- lapply(errors, sd)
+
+  out
+}