Two station static model: Stan file, data, data prep, implementation

DESCRIPTION

@@ -39,9 +39,12 @@ Imports:
     lubridate,
     magrittr,
     methods,
+    parallel,
+    stats,
     tibble (>= 1.1.0),
     tidyr,
-    unitted (>= 0.2.8)
+    unitted (>= 0.2.8),
+    utils
 Suggests:
     chron,
     devtools,
@@ -80,6 +83,7 @@ Collate:
     'create_calc_DO.R'
     'create_calc_NLL.R'
     'create_calc_dDOdt.R'
+    'data.R'
     'data_metab.R'
     'deprecated.R'
     'load_french_creek.R'
@@ -92,7 +96,9 @@ Collate:
     'specs-class.R'
     'metab_model-class.R'
     'metab_Kmodel.R'
+    'mm_time_by_date_matrix.R'
     'metab_bayes.R'
+    'metab_bayes_2s.R'
     'metab_inputs.R'
     'metab_mle.R'
     'metab_model.get_param_names.R'
@@ -104,6 +110,7 @@ Collate:
     'metab_sim.R'
     'mm_check_mcmc_file.R'
     'mm_data.R'
+    'mm_determine_cores.R'
     'mm_filter_dates.R'
     'mm_filter_hours.R'
     'mm_filter_valid_days.R'
@@ -129,5 +136,6 @@ Collate:
     'streamMetabolizer-deprecated.R'
     'streamMetabolizer.R'
     'zz_build_docs.R'
-RoxygenNote: 7.1.1
+RoxygenNote: 7.3.3
 Encoding: UTF-8
+LazyData: true

NAMESPACE

@@ -12,16 +12,19 @@ S3method(get_param_names,character)
 S3method(get_param_names,metab_model)
 S3method(get_params,metab_Kmodel)
 S3method(get_params,metab_bayes)
+S3method(get_params,metab_bayes_2s)
 S3method(get_params,metab_model)
 S3method(get_params,metab_sim)
 S3method(get_specs,metab_model)
 S3method(get_version,metab_model)
 S3method(predict_DO,metab_Kmodel)
+S3method(predict_DO,metab_bayes_2s)
 S3method(predict_DO,metab_model)
 S3method(predict_DO,metab_night)
 S3method(predict_DO,metab_sim)
 S3method(predict_metab,metab_Kmodel)
 S3method(predict_metab,metab_bayes)
+S3method(predict_metab,metab_bayes_2s)
 S3method(predict_metab,metab_model)
 S3method(print,logs_metab)
 S3method(print,specs)
@@ -69,6 +72,7 @@ export(lookup_usgs_elevation)
 export(metab)
 export(metab_Kmodel)
 export(metab_bayes)
+export(metab_bayes_2s)
 export(metab_inputs)
 export(metab_mle)
 export(metab_model)
@@ -96,6 +100,7 @@ export(sim_pred_Kb)
 export(specs)
 exportClasses(metab_Kmodel)
 exportClasses(metab_bayes)
+exportClasses(metab_bayes_2s)
 exportClasses(metab_mle)
 exportClasses(metab_model)
 exportClasses(metab_night)
@@ -114,7 +119,6 @@ importFrom(LakeMetabolizer,sw.to.par.base)
 importFrom(graphics,abline)
 importFrom(graphics,plot)
 importFrom(graphics,points)
-importFrom(lazyeval,lazy_dots)
 importFrom(lifecycle,deprecate_warn)
 importFrom(lifecycle,deprecated)
 importFrom(lifecycle,is_present)
@@ -124,6 +128,9 @@ importFrom(lubridate,is.Date)
 importFrom(lubridate,is.POSIXct)
 importFrom(lubridate,tz)
 importFrom(lubridate,with_tz)
+importFrom(rlang,as_name)
+importFrom(rlang,enquos)
+importFrom(rlang,quo_is_null)
 importFrom(stats,approx)
 importFrom(stats,approxfun)
 importFrom(stats,coef)
@@ -162,6 +169,7 @@ importFrom(utils,available.packages)
 importFrom(utils,capture.output)
 importFrom(utils,contrib.url)
 importFrom(utils,head)
+importFrom(utils,modifyList)
 importFrom(utils,packageVersion)
 importFrom(utils,read.csv)
 importFrom(utils,tail)

R/data.R

@@ -0,0 +1,44 @@
+#' Example two-station (VFTS) input data
+#'
+#' A 30-day example dataset for fitting a two-station (upstream/downstream,
+#' Variable Flow Two-Station) metabolism model with
+#' \code{\link{metab_bayes_2s}}. It is a subset of the \code{VFTS-2}
+#' (variable-travel-time) run from the published two-station metabolism modeling
+#' dataset for a reach of the Colorado River in Glen Canyon, covering 2011-07-31
+#' through 2011-08-29 at the source data's native 15-minute timestep, plus a
+#' lead-in block of upstream DO observations (exactly as long as the longest
+#' travel time in the dataset requires -- see \code{\link{metab_bayes_2s}}'s
+#' "Two-station data requirements" section) immediately before 2011-07-31.
+#'
+#' @format A data.frame with 2904 rows and the 9 columns expected by
+#'   \code{\link{metab_bayes_2s}}'s \code{data} argument, each carrying
+#'   \code{\link[unitted]{unitted}} units matching \code{\link{mm_data}}:
+#'   \describe{
+#'     \item{solar.time}{POSIXct timestamp, UTC}
+#'     \item{DO.obs.up}{dissolved oxygen observed at the upstream station,
+#'       mgO2 L^-1}
+#'     \item{DO.sat.up}{dissolved oxygen at equilibrium saturation at the
+#'       upstream station, mgO2 L^-1}
+#'     \item{DO.obs.down}{dissolved oxygen observed at the downstream
+#'       station, mgO2 L^-1}
+#'     \item{DO.sat.down}{dissolved oxygen at equilibrium saturation at the
+#'       downstream station, mgO2 L^-1}
+#'     \item{light}{photosynthetically active radiation, umol m^-2 s^-1}
+#'     \item{depth}{reach depth, m}
+#'     \item{temp.water}{water temperature at the downstream station, degC}
+#'     \item{travel.time}{reach travel time between the upstream and
+#'       downstream stations, d}
+#'   }
+#'
+#' @source Filtered to \code{model_run == 'VFTS-2'}; see
+#'   \code{data-raw/two_station_example.R} for the extraction/renaming code.
+#'
+#'   Bishop, I.W., Deemer, B.R., Kennedy, T.A., Payn, R.A., Hall Jr, R.O. and
+#'   Yackulic, C.B., 2026. A simplified two-station approach for modeling
+#'   metabolism in dam tailwaters subject to diel flow variation. Limnology
+#'   and Oceanography: Methods, p.e70066.
+#'   \url{https://aslopubs.onlinelibrary.wiley.com/doi/pdf/10.1002/lom3.70066}
+#'
+#'   Data archived at ScienceBase:
+#'   \url{https://www.sciencebase.gov/catalog/item/6887d457d4be024722b4aae2}
+"two_station_example"

R/metab.R

@@ -59,11 +59,12 @@ metab <- function(specs=specs(mm_name()), data=v(mm_data(NULL)), data_daily=v(mm
   model_type <- mm_parse_name(specs$model_name)$type
   metab_fun <- switch(
     model_type,
-    bayes  = metab_bayes,
-    Kmodel = metab_Kmodel,
-    mle    = metab_mle,
-    night  = metab_night,
-    sim    = metab_sim)
+    bayes           = metab_bayes,
+    bayes_2s        = metab_bayes_2s,
+    Kmodel          = metab_Kmodel,
+    mle             = metab_mle,
+    night           = metab_night,
+    sim             = metab_sim)
 
   # run the model
   metab_fun(specs=specs, data=data, data_daily=data_daily, info=info)

R/metab_bayes.R

@@ -1,4 +1,4 @@
-#' @include metab_model-class.R
+#' @include metab_model-class.R mm_time_by_date_matrix.R
 NULL
 
 #' Basic Bayesian metabolism model fitting function
@@ -330,9 +330,11 @@ bayes_1ply <- function(
 #' Called from metab_bayes().
 #'
 #' @param data_all data.frame of the form \code{mm_data(solar.time, DO.obs,
-#'   DO.sat, depth, temp.water, light)} and containing data for just one
-#'   estimation-day (this may be >24 hours but only yields estimates for one
-#'   24-hour period)
+#'   DO.sat, depth, temp.water, light)} containing the full (possibly
+#'   multi-day) filtered dataset for this model - unlike \code{bayes_1ply()}'s
+#'   \code{data_ply}, which receives one estimation-day at a time,
+#'   \code{bayes_allply()} is called once with all valid dates together (used
+#'   when \code{specs$split_dates==FALSE})
 #' @param data_daily_all data.frame of daily priors, if appropriate to the given
 #'   model_path
 #' @param removed data.frame of dates that were removed and why
@@ -479,15 +481,13 @@ prepdata_bayes <- function(
     )
     stop("dates have differing numbers of rows; observations cannot be combined in matrix")
   }
-  time_by_date_matrix <- function(vec) {
-    matrix(data=vec, nrow=num_daily_obs, ncol=num_dates, byrow=FALSE)
-  }
+  time_by_date_matrix <- mm_time_by_date_matrix(num_daily_obs, num_dates)
 
   # double-check that our dates are going to line up with the input dates. this
   # should be redundant w/ above date_table checks, so just being extra careful
-  obs_dates <- time_by_date_matrix(format(data$date, format="%Y-%m-%d"))
-  unique_dates <- apply(obs_dates, MARGIN=2, FUN=function(timevec) unique(timevec))
-  if(!all.equal(unique_dates, names(date_table))) stop("couldn't fit given dates into matrix")
+  mm_check_dates_contiguous(
+    time_by_date_matrix(format(data$date, format="%Y-%m-%d")), date_table,
+    "couldn't fit given dates into matrix")
 
   # confirm that every day has the same modal timestep and put a value on that
   # timestep. the tolerance for uniqueness within each day is set by the default
@@ -625,10 +625,7 @@ runstan_bayes <- function(
   verbose=FALSE, ...) {
 
   # determine how many cores to use
-  tot_cores <- detectCores()
-  if (!is.finite(tot_cores)) { tot_cores <- 1 }
-  n_cores <- min(tot_cores, n_cores)
-  if(verbose) message(paste0("MCMC (","Stan","): requesting ",n_chains," chains on ",n_cores," of ",tot_cores," available cores"))
+  n_cores <- mm_determine_cores(n_cores, n_chains=n_chains, verbose=verbose)
 
   # stan() can't find its own function cpp_object_initializer() unless the
   # namespace is loaded. requireNamespace is somehow not doing this. Thoughts