migrate caching to get_ZB18a, only do it upon request

CRAN asked me not to write to the user's directory by default. We now have: - get_solution() - which calls get_ZB18a() (which now documents the OS) - which downloads and caches the OS if needed, or reads from cache (unless force == TRUE). - snvec() the main function
japhir · May 9, 2023 · c66f727 · c66f727
1 parent 7dc49f1
commit c66f727
Show file tree

Hide file tree

Showing 18 changed files with 243 additions and 148 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -31,7 +31,8 @@ Imports:
     rlang (>= 0.4.11),
     withr,
     tidyr,
-    glue
+    glue,
+    curl
 Config/testthat/edition: 3
 Depends:
     R (>= 3.5.0)

diff --git a/NAMESPACE b/NAMESPACE
@@ -1,6 +1,8 @@
 # Generated by roxygen2: do not edit by hand
 
 export(.data)
+export(get_ZB18a)
+export(get_solution)
 export(snvec)
 importFrom(cli,cli_abort)
 importFrom(cli,cli_inform)

diff --git a/R/data.R b/R/data.R
@@ -5,8 +5,11 @@
 #'
 #' The HNBody output of Zeebe & Lourens (2019).
 #'
-#' The wikipedia page on [Orbital elements](https://en.wikipedia.org/wiki/Orbital_elements)
-#' describes what the components relate to in order to uniquely specify an orbital plane.
+#' The wikipedia page on [Orbital
+#' elements](https://en.wikipedia.org/wiki/Orbital_elements) describes what the
+#' components relate to in order to uniquely specify an orbital plane. The
+#' function asks to download the files to the user's cache directory so that they
+#' can be accessed more quickly in the future.
 #'
 #' @format ## `ZB18a`
 #' A data frame with 250,001 rows and 20 columns:

diff --git a/R/get_solution.R b/R/get_solution.R
@@ -3,9 +3,16 @@
 #' @param orbital_solution Character vector with the name of the orbital
 #'   solution to use. One of `"ZB18a"` (default) from Zeebe and Lourens (2019),
 #'   or `"La11"` (not yet implemented!).
+# quiet and force are documented in get_ZB18a
+#' @inheritParams get_ZB18a
 #' @seealso [get_ZB18a()]
 #' @inherit get_ZB18a references
-get_solution <- function(orbital_solution = "ZB18a") {
+#' @returns `get_solution()` returns a [tibble][tibble::tibble-package] with the
+#'   orbital solution input and some preprocessed new columns.
+#' @examples
+#' get_solution()
+#' @export
+get_solution <- function(orbital_solution = "ZB18a", quiet = FALSE) {
   solutions <- c("ZB18a", "La11")
   if (!orbital_solution %in% solutions) {
     cli::cli_abort(c("{.var orbital_solution} must be one of: {.or {.q {solutions}}}",
@@ -14,7 +21,7 @@ get_solution <- function(orbital_solution = "ZB18a") {
 
   if (orbital_solution == "ZB18a") {
     # read in the (new?) cached result
-    dat <- get_ZB18a()
+    dat <- get_ZB18a(quiet = quiet)
   }
   if (orbital_solution == "La11") {
     ## dat <- snvecR::La11
@@ -27,16 +34,99 @@ get_solution <- function(orbital_solution = "ZB18a") {
   return(dat)
 }
 
-#' Get the ZB18a solution from the cache
+#' Get the ZB18a solution
 #'
-# The cache is created in when the package is loaded, in zzz.R.
+#' @param quiet Be quiet?
+#'
+#'   * If `TRUE`, hide info messages.
+#'
+#'   * If `FALSE` (the default) print info messages and timing.
+#' @param force Force re-downloading the results, even if the solution is saved
+#'   to the cache.
+#' @returns `get_ZB18a()` returns a [tibble][tibble::tibble-package] with the
+#'   orbital solution input and some preprocessed new columns.
 #' @rdname ZB18a
-get_ZB18a <- function() {
+#' @seealso [prepare_solution()] Processes orbital solution input to include
+#'   helper columns.
+#' @examples
+#' get_ZB18a()
+#' @export
+get_ZB18a <- function(quiet = FALSE, force = FALSE) {
+  ZB18a_url <- "http://www.soest.hawaii.edu/oceanography/faculty/zeebe_files/Astro/PrecTilt/OS/ZB18a/ems-plan3.dat"
+
   cachedir <- tools::R_user_dir("snvecR", which = "cache")
-  ZB18apath <- paste0(cachedir, "/ZB18a.rds")
+  raw_path <- paste0(cachedir, "/ems-plan3.dat")
+  csv_path <- paste0(cachedir, "/ZB18a.csv")
+  rds_path <- paste0(cachedir, "/ZB18a.rds")
+
+  # if it doesn't exist, or the user wants to override the file
+  if (!file.exists(rds_path) || force) {
+    if (!file.exists(csv_path) || force) {
+      if (!quiet) cli::cli_alert_info("The orbital solution ZB18a has not been downloaded.")
+
+      # default to Yes downloading if not interactive (i.e. GitHub actions)
+      if (!interactive()) {
+        download <- TRUE
+        # default to no caching, so that we're not writing to the user's directory
+        save_cache <- FALSE
+      } else {# we're interactive
+        # a logical, TRUE if Yes, no if otherwise
+        download <- utils::menu(c("Yes", "No"), title = "Would you like to download and process it now?") == 1L
+      }
+
+      # the user would NOT like to download and process the orbital solution
+      if (!download) {
+        cli::cli_abort("Cannot `get_ZB18a()` without downloading the file.")
+      } else {# the user would like to download and process the orbital solution
+        if (interactive()) {
+          save_cache <- utils::menu(c("Yes", "No"), title = "Would you like to save the results to the cache?") == 1L
+        }
+
+        if (!quiet) cli::cli_alert_info("Reading {.file {basename(raw_path)}} from website {.url {ZB18a_url}}.")
+
+        # read the file from the website
+        ZB18a <- readr::read_table(ZB18a_url,
+                                   col_names = c("t", # time in days
+                                                 "aa", # semimajor axis
+                                                 "ee", # eccentricity
+                                                 "inc", # inclination
+                                                 "lph", # long periapse
+                                                 "lan", # long ascending node
+                                                 "arp", # argument of periapse
+                                                 "mna"), # mean anomaly
+                                   skip = 3,
+                                   comment = "#") #|>
+
+        if (save_cache) {
+          dir.create(cachedir, recursive = TRUE, showWarnings = FALSE)
+          if (!quiet) cli::cli_alert_info("The cache directory is {.file {cachedir}}.")
+          # also copy the raw file to disk, even though we've read it in using read_table directly
+          curl::curl_download(ZB18a_url, destfile = raw_path)
+          if (!quiet) cli::cli_alert_info("Saved {.file {basename(raw_path)}} to cache.")
+
+          # write intermediate result to csv
+          readr::write_csv(ZB18a, csv_path)
+          if (!quiet) cli::cli_alert_info("Saved cleaned-up {.file {basename(csv_path)}} to cache.")
+        }
+      }
+    } else {# if we've downloaded the file but haven't prepared the solution somehow
+      ZB18a <- readr::read_csv(csv_path)
+    }
+
+    # prepare the solution (i.e. calculate helper columns)
+    ZB18a <- ZB18a |>
+      prepare_solution()
+
+    if (save_cache) {
+      readr::write_rds(ZB18a, rds_path)
+      if (!quiet) {
+        cli::cli_alert("Saved solution with helper columns {.file {basename(rds_path)}} to cache.",
+                       "i" = "Future runs will read from the cache, unless you specify `force = TRUE`.")
+      }
+    }
+  } else {# if the rds file already exist, read it from the cache
+    ZB18a <- readr::read_rds(rds_path)
+  }
 
-  # if you picked no during the onLoad, I assume read_rds will tell you the
-  # file doesn't exist
-  ZB18a <- readr::read_rds(ZB18apath)
   return(ZB18a)
 }
diff --git a/R/interpolation.R b/R/interpolation.R
@@ -1,47 +1,40 @@
-## using approxfun
-## http://desolve.r-forge.r-project.org/ has an article on time-varying inputs
-## we can use approxfun to generate a function that approximates =col= for timestep t.
+# using approxfun
+# http://desolve.r-forge.r-project.org/ has an article on time-varying inputs
+# we can use approxfun to generate a function that approximates =col= for timestep t.
 
-
-##' Linearly interpolate a dataframe.
-##'
-##' @param dat The dataframe.
-##' @param col A character vector with the column we want to interpolate on.
-##' @returns A function that predicts `col` as a function of `t` in `dat`.
-##'
-##' @seealso [qinterp()] for quick interpolation of a single timestep.
-##' @importFrom stats approxfun
-##' @importFrom dplyr select
-##' @noRd
+#' Linearly interpolate a dataframe.
+#'
+#' @param dat The dataframe.
+#' @param col A character vector with the column we want to interpolate on.
+#' @returns A function that predicts `col` as a function of `t` in `dat`.
+#'
+#' @seealso [qinterp()] for quick interpolation of a single timestep.
+#' @noRd
 approxdat <- function(dat, col) {
   # I'm not putting any input checks because it's an internal function
   dat |>
     dplyr::select(all_of(c("t", col))) |>
     approxfun(rule = 2)
 }
 
-## implement qinterp similar to the C-routine
-## :PROPERTIES:
-## :header-args:R: :tangle R/interpolation.R :comments org :session *R:snvec-R* :exports both :results output :eval no-export
-## :END:
-
-##' Quickly interpolate a single value.
-##'
-##' `qinterp()` is a custom linear interpolation algorithm that is much faster
-##' than using the full vectorized `[approx()]` or `[approxfun()]`, because it
-##' only interpolates the single value of the current timestep.
-##'
-##' @param y The vector to interpolate.
-##' @param ds The difference in timestep in the astronomical solution.
-##' @param dx The difference between the current timestep and the timestep in the astronomical solution.
-##' @param m The index variable of the current position in the astronomical solution.
-##' @returns The vector of interpolated results.
-##' @examples
-##' # interpolate ZB18a$lph[[1:4]]
-##' qinterp(ZB18a$lph[[1:4]], ds = -146100, dx =  -18262.5, m = 2)
-##'
-##' @seealso [approxdat()] for linear interpolation of the full astronomical solution.
-##' @noRd
+# implement qinterp similar to the C-routine
+#' Quickly interpolate a single value.
+#'
+#' `qinterp()` is a custom linear interpolation algorithm that is much faster
+#' than using the full vectorized `[approx()]` or `[approxfun()]`, because it
+#' only interpolates the single value of the current timestep.
+#'
+#' @param y The vector to interpolate.
+#' @param ds The difference in timestep in the astronomical solution.
+#' @param dx The difference between the current timestep and the timestep in the astronomical solution.
+#' @param m The index variable of the current position in the astronomical solution.
+#' @returns The vector of interpolated results.
+#' @examples
+#' # interpolate ZB18a$lph[[1:4]]
+#' qinterp(ZB18a$lph[[1:4]], ds = -146100, dx =  -18262.5, m = 2)
+#'
+#' @seealso [approxdat()] for linear interpolation of the full astronomical solution.
+#' @noRd
 qinterp <- function(y, ds, dx, m) {
   yi <- y[m]
   dy <- 0

diff --git a/R/prepare_solution.R b/R/prepare_solution.R
@@ -9,7 +9,8 @@
 #' * `lph` Longitude of perihelion \eqn{\varpi} (degrees).
 #' * `lan` Longitude of the ascending node \eqn{\Omega} (degrees).
 #' * `inc` Inclination \eqn{I} (degrees).
-#'
+# inherit quiet
+#' @inheritParams get_ZB18a
 #' @returns A [tibble][tibble::tibble-package] with the new columns added.
 #' @seealso [get_ZB18a()] [get_solution()]
 #'
@@ -43,18 +44,18 @@
 #  \item{npx, npy, npz}{The \eqn{x}, \eqn{y}, and \eqn{z}-components of the unit
 #   normal vector \eqn{\vec{n}'}{n'}, relative to ECLIPJ2000.}
 # IOP = instantaneous orbit plane
-prepare_solution <- function(data) {
+prepare_solution <- function(data, quiet = FALSE) {
   if (!all(c("lph", "lan", "t", "inc") %in% colnames(data))) {
     cli::cli_abort("Column names 't', 'lph', 'lan', 'inc' must be in data.")
   }
 
-  cli::cli_alert_info("Calculating helper columns")
+  if (!quiet) cli::cli_alert_info("Calculating helper columns.")
   data |>
     dplyr::mutate(
       lphu = unwrap(.data$lph),
       lanu = unwrap(.data$lan)
     ) |>
-    dplyr::mutate(age = -.data$t / KY2D, .after = .data$t) |>
+    dplyr::mutate(age = -.data$t / KY2D, .after = "t") |>
     dplyr::mutate(
       hh = .data$ee * sin(.data$lph / R2D),
       kk = .data$ee * cos(.data$lph / R2D),

diff --git a/R/snvec.R b/R/snvec.R
@@ -11,6 +11,7 @@
 #'   Defaults to `1.0`.
 #' @param td Tidal dissipation \eqn{T_{d}}{Td}, normalized to modern. Defaults
 #'   to `0.0`.
+# orbital_solution comes from get_solution
 #' @inheritParams get_solution
 #' @param tres Output timestep resolution in thousands of years (kyr). Defaults
 #'   to `0.4`.
@@ -21,12 +22,9 @@
 #' @param solver Character vector specifying the method passed to
 #'   [deSolve::ode()]'s `method`. Defaults to `"vode"` for stiff problems
 #'   with a variable timestep.
-#' @param quiet Be quiet?
-#'
-#'   * If `TRUE`, hide info messages.
-#'
-#'   * If `FALSE` (the default) print info messages and timing.
-#'
+# quiet comes from get_ZB18a. Force does too, but I hope that because we don't
+# use it here, it won't get inherited.
+#' @inheritParams get_ZB18a
 #' @param output Character vector with name of desired output. One of:
 #'
 #'   * `"nice"` (the default) A [tibble][tibble::tibble-package] with the
@@ -45,8 +43,7 @@
 #'   2010) means that the R routine returns an evenly-spaced time grid, whereas
 #'   the C-routine has a variable time-step.
 #'
-#' @returns `snvec()` returns different output depending on the `outputs`
-#'   argument.
+#' @returns `snvec()` returns different output depending on the `outputs` argument.
 #'
 #' If `output = "nice"` (the default), returns a
 #' [tibble][tibble::tibble-package] with the following columns:
@@ -102,8 +99,15 @@
 #' `sy`, and `sz` (see above). This can be useful for i.e.
 #' [deSolve::diagnostics()].
 #'
-#' @seealso [deSolve::ode()] from Soetaert et al., (2010) for the ODE solver
-#'   that we use.
+#' @seealso
+#'
+#' * [deSolve::ode()] from Soetaert et al., (2010) for the ODE solver that we
+#'   use.
+#'
+#' * [get_ZB18a()] Documents the default orbital solution input.
+#'
+#' * [get_solution()] A general function that in the future may be used to get
+#'   other orbital solutions.
 #'
 #' @references
 #'
@@ -144,22 +148,20 @@ snvec <- function(tend = -1e3,
                   quiet = FALSE,
                   output = "nice") {
 
-  dat <- get_solution(orbital_solution = orbital_solution)
-
   outputs <- c("nice", "all", "ode")
   if (!output %in% outputs) {
-    cli::cli_abort(c("{.var output} must be one of {.or {.q {outputs}}}",
-                     "x" = "You've supplied {.q {output}}"))
+    cli::cli_abort(c("{.var output} must be one of {.or {.q {outputs}}}.",
+                     "x" = "You've supplied {.q {output}}."))
   }
   ## tend
   if (tend >= 0) {
-    cli::cli_abort(c("{.var tend} must be < 0",
-      "x" = "You've supplied {tend}"
+    cli::cli_abort(c("{.var tend} must be < 0.",
+      "x" = "You've supplied {tend}."
     ))
   }
-  if (tend < min(dat$t / KY2D)) {
-    cli::cli_abort(c("{.var tend} must be > the orbital solution {min(dat$t)/KY2D}",
-      "x" = "You've supplied {tend}"
+  if (tend < -1e5) {
+    cli::cli_abort(c("{.var tend} must be > the orbital solution {-1e5}, e.g. -100 Ma.",
+      "x" = "You've supplied {tend}."
     ))
   }
 
@@ -182,13 +184,13 @@ snvec <- function(tend = -1e3,
   if (ed < .9 | ed > 1.1) {
     cli::cli_warn(c("!" = "Dynamic ellipticity likely varied between 0.9980 and 1.0005 during the past 45 Ma!",
                     "i" = "{.var ed} = {ed}",
-                    "*" = "See Zeebe & Lourens 2022 Pal&Pal <https://doi.org/10.1029/2021PA004349>"))
+                    "*" = "See Zeebe & Lourens 2022 Pal&Pal <https://doi.org/10.1029/2021PA004349>."))
   }
 
   if (td < 0 | td > 1.2) {
     cli::cli_warn(c("Tidal dissipation likely varied between 0 and 1!",
                     "i" = "{.var td} = {td}",
-                    "*" = "See Zeebe & Lourens 2022 Pal&Pal <https://doi.org/10.1029/2021PA004349>"))
+                    "*" = "See Zeebe & Lourens 2022 Pal&Pal <https://doi.org/10.1029/2021PA004349>."))
   }
 
   if (atol < 1e-12 | atol > 1e-3) {
@@ -198,6 +200,8 @@ snvec <- function(tend = -1e3,
     cli::cli_warn("Input relative tolerance should be smaller than 1e-3.")
   }
 
+  dat <- get_solution(orbital_solution = orbital_solution)
+
   # message user about inputs
   if (!quiet) {
     startdate <- lubridate::now()

diff --git a/R/snvecR-package.R b/R/snvecR-package.R
@@ -10,6 +10,7 @@
 #' @importFrom purrr map2_dbl
 #' @importFrom readr write_rds read_rds read_table
 #' @importFrom rlang .data
+#' @importFrom stats approxfun
 #' @importFrom tibble tibble
 #' @importFrom tidyr unnest
 #' @importFrom tidyselect all_of