diff --git a/DESCRIPTION b/DESCRIPTION index 4af05e3c..b3c27f84 100644 --- a/DESCRIPTION +++ b/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 diff --git a/NAMESPACE b/NAMESPACE index b161194a..eb89725f 100644 --- a/NAMESPACE +++ b/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) diff --git a/R/data.R b/R/data.R new file mode 100644 index 00000000..64c1c255 --- /dev/null +++ b/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" diff --git a/R/metab.R b/R/metab.R index 8a45c07a..4d55efb5 100644 --- a/R/metab.R +++ b/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) diff --git a/R/metab_bayes.R b/R/metab_bayes.R index 9db1c0d0..43d49c3e 100644 --- a/R/metab_bayes.R +++ b/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 diff --git a/R/metab_bayes_2s.R b/R/metab_bayes_2s.R new file mode 100644 index 00000000..d1582703 --- /dev/null +++ b/R/metab_bayes_2s.R @@ -0,0 +1,498 @@ +#' @include metab_model-class.R metab_bayes.R +NULL + +# Suppress R CMD CHECK NOTEs for column names used as unbound globals in +# dplyr NSE calls (e.g. mutate, summarise). These are data frame column +# names resolved at runtime, not missing variable declarations. +utils::globalVariables(c(".", "metab_50pct", "DO.mod.down")) + +#' Two-station Bayesian metabolism model fitting function +#' +#' Fits a two-station (upstream/downstream, Variable Flow Two-Station) Bayesian +#' model to estimate GPP, ER, and K600 from paired upstream and downstream DO, +#' temperature, light, and travel-time data, using the single fixed Stan model +#' in \code{inst/models/b2_np_oi_tr_plrckm.stan}. See \code{\link{mm_name}} to +#' choose a Bayesian model and \code{\link{specs}} for relevant options for the +#' \code{specs} argument. +#' +#' Unlike \code{\link{metab_bayes}}, which supports many model structures via +#' \code{split_dates}/\code{pool_K600}/etc., \code{metab_bayes_2s} always +#' fits every date jointly in a single Stan call (\code{specs$split_dates} is +#' forced to \code{FALSE} by \code{\link{specs}}), because the +#' upstream-downstream lag shift ties each date's first modeled rows to the +#' previous date's last rows. +#' +#' @inheritParams metab +#' @return A metab_bayes_2s object containing the fitted model. This object +#' can be inspected with the functions in the +#' \code{\link{metab_model_interface}} and also \code{\link{get_mcmc}}. +#' +#' @section Two-station data requirements: In addition to the checks +#' performed by \code{\link{mm_validate_data}}, \code{data$travel.time} (the +#' reach travel time between the upstream and downstream stations, in days) +#' must be strictly positive and no greater than 8/24 days (8 hours). +#' Values above this limit either indicate travel time was supplied in the +#' wrong units (e.g., minutes or hours instead of days), or reflect a reach +#' whose actual travel time exceeds the 8-hour limit required to prevent +#' the previous day's light conditions from influencing the following +#' day's metabolism estimate. There must also be enough lead-in +#' observations of upstream DO before the first row of \code{data} to +#' cover the longest travel time in the dataset. +#' +#' @export +#' @family metab_model +#' @importFrom utils modifyList +metab_bayes_2s <- function( + specs=specs(mm_name('bayes_2s')), + data=mm_data(solar.time, DO.obs.up, DO.sat.up, DO.obs.down, DO.sat.down, + light, depth, temp.water, travel.time), + data_daily=mm_data(date, optional='all'), + info=NULL +) { + + stanfit <- NULL + fitting_time <- system.time({ + # Check data for correct column names, units, and travel.time bounds + # (mm_validate_data()), then check lead-in coverage + # (mm_validate_data_2station()), before any data prep begins + dat_list <- mm_validate_data(data, data_daily, 'metab_bayes_2s') + mm_validate_data_2station(dat_list$data) + + data_v <- v(dat_list$data) + travel_time <- data_v$travel.time + solar_time <- data_v$solar.time + + # Reconstruct the same "modeled rows" (post-lead-in-trim) index set that + # prepdata_bayes_2s() computes internally, using the identical + # timestep_days/lag/max_lag formula, so that Stan's date_index/time_index + # can be mapped back to actual dates and solar.times for the daily and + # instantaneous results below. (Duplicated here rather than exposed by + # prepdata_bayes_2s(), which returns only the Stan-ready matrices.) + timestep_days <- stats::median(as.numeric(diff(solar_time), units='days')) + max_lag <- max(round(travel_time / timestep_days)) + n_total <- nrow(dat_list$data) + keep <- seq.int(max_lag + 1, n_total) + modeled_solar_time <- solar_time[keep] + modeled_dates <- as.Date(modeled_solar_time) + date_df <- tibble::tibble(date=unique(modeled_dates), date_index=seq_along(unique(modeled_dates))) + n_days <- nrow(date_df) + if(length(keep) %% n_days != 0) { + stop(paste0( + 'dates have differing numbers of modeled rows after lead-in removal; ', + 'observations cannot be combined into a matrix: ', + paste(sprintf('%s (%d rows)', names(table(modeled_dates)), table(modeled_dates)), collapse=', '))) + } + n_obs <- length(keep) / n_days + obs_index_df <- tibble::tibble( + solar.time=modeled_solar_time, + DO.obs.down=data_v$DO.obs.down[keep], + date_index=rep(date_df$date_index, each=n_obs), + time_index=rep(seq_len(n_obs), times=n_days)) + + # Prepare the Stan data list (matrices from data, plus scalar priors from + # specs). modifyList (not c()) is used because prepdata_bayes_2s() already + # supplies K600_lnorm_meanlog/K600_lnorm_sdlog (read from specs), and + # those two names are also in specs$params_in; a plain c() would create + # duplicate-named list elements instead of overriding + data_list <- prepdata_bayes_2s(dat_list$data, specs=specs) + data_list <- modifyList(data_list, specs[specs$params_in]) + + # Check and parse model file path + specs$model_path <- mm_locate_filename(specs$model_name) + + # determine how many cores to use, as in runstan_bayes() + n_cores <- mm_determine_cores(specs$n_cores, n_chains=specs$n_chains, verbose=specs$verbose) + + # Fit the model, collecting errors/warnings as strings rather than + # letting a bad dataset halt execution without reporting anything back + stop_strs <- character(0) + warn_strs <- character(0) + daily <- NULL + inst <- NULL + withCallingHandlers( + tryCatch({ + if (!requireNamespace("rstan", quietly = TRUE)) stop("rstan is required but not installed. Install it with: install.packages('rstan')") + + consolelog <- utils::capture.output( + stanfit <- rstan::stan( + file=specs$model_path, data=data_list, pars=specs$params_out, + chains=specs$n_chains, cores=n_cores, + iter=specs$burnin_steps + specs$saved_steps, warmup=specs$burnin_steps, + thin=specs$thin_steps, verbose=specs$verbose, open_progress=FALSE), + split=specs$verbose) + + if(stanfit@mode == 2L) { + # mirror runstan_bayes()'s warn-and-continue pattern for a failed + # run: report the diagnostic as a warning (caught below) and skip + # the post-processing steps, which all assume a successful fit + warning(paste(utils::capture.output(print(stanfit)), collapse='\n')) + } else { + # format the Stan summary matrix into per-variable data.frames + stan_mat <- rstan::summary(stanfit)$summary + mcmc_out <- format_mcmc_mat_nosplit( + stan_mat, data_list$n_days, data_list$n_obs, specs$model_name, + keep_mcmc=isTRUE(specs$keep_mcmcs), stanfit) + + # daily GPP/ER/K600 estimates: join Stan's date_index back to dates + date_index <- time_index <- index <- '.dplyr.var' + daily <- mcmc_out$daily %>% + dplyr::left_join(date_df, by='date_index') %>% + dplyr::select(-date_index, -time_index, -index) %>% + dplyr::select(date, dplyr::everything()) + + # instantaneous DO.mod.down estimates come from the 'metab' Stan + # transformed parameter (posterior median), which format_mcmc_mat_nosplit() + # buckets by row count rather than by name since 'metab' isn't in its + # par_homes lookup table; find that bucket by its column names instead + is_metab_bucket <- sapply(mcmc_out, function(df) is.data.frame(df) && any(grepl('^metab_', names(df)))) + metab_bucket_name <- names(mcmc_out)[is_metab_bucket][1] + if(is.na(metab_bucket_name)) { + stop("could not find 'metab' in the Stan output; check that specs$params_out includes 'metab'") + } + inst <- mcmc_out[[metab_bucket_name]] %>% + dplyr::select(date_index, time_index, DO.mod.down=metab_50pct) %>% + dplyr::inner_join(obs_index_df, by=c('date_index','time_index')) %>% + dplyr::select(solar.time, DO.obs.down, DO.mod.down) %>% + dplyr::arrange(solar.time) + } + + }, error=function(err) { + stop_strs <<- c(stop_strs, err$message) + }), warning=function(war) { + warn_strs <<- c(warn_strs, war$message) + invokeRestart("muffleWarning") + }) + + # if fitting failed, fill in NA daily estimates (with real dates) so the + # returned model at least reports which dates were attempted + if(length(stop_strs) > 0 || is.null(daily)) { + na_vec <- rep(as.numeric(NA), nrow(date_df)) + daily <- data.frame( + date=date_df$date, + GPP_daily_2.5pct=na_vec, GPP_daily_50pct=na_vec, GPP_daily_97.5pct=na_vec, + ER_daily_2.5pct=na_vec, ER_daily_50pct=na_vec, ER_daily_97.5pct=na_vec, + K600_daily_2.5pct=na_vec, K600_daily_50pct=na_vec, K600_daily_97.5pct=na_vec) + inst <- NULL + } + daily <- dplyr::mutate(daily, valid_day=TRUE, warnings='', errors='') + + fit <- list( + daily=daily, inst=inst, + warnings=trimws(unique(warn_strs)), errors=trimws(unique(stop_strs))) + }) + + # Package and return results + mm <- metab_model( + model_class="metab_bayes_2s", + info=info, + fit=fit, + log=NULL, + mcmc=if(isTRUE(specs$keep_mcmcs)) stanfit else NULL, + mcmc_data=if(isTRUE(specs$keep_mcmc_data)) data_list else NULL, + fitting_time=fitting_time, + compile_time=system.time({}), # rstan::stan() compiles & samples in one call; not timed separately + specs=specs, + data=dat_list$data, # keep the units if given + data_daily=dat_list$data_daily) + + # Update data with DO predictions + success <- !is.null(fit$inst) && length(fit$errors) == 0 + if(success) { + mm@data <- predict_DO(mm) + } else { + warntxt <- paste0( + 'Modeling failed\n', + if(length(fit$warnings) > 0) paste0(' Warnings:\n', paste0(' ', fit$warnings, collapse='\n')), + if(length(fit$errors) > 0) paste0(' Errors:\n', paste0(' ', fit$errors, collapse='\n'))) + warning(warntxt) + } + + # Return + mm +} + + +#' Reshape long-format two-station data into the list expected by the +#' two-station Stan model +#' +#' Time-shifts the upstream DO series to match the travel time between +#' stations, then pivots the result into the \code{n_obs x n_days} matrices +#' expected by the \code{data} block of \code{inst/models/b2_np_oi_tr_plrckm.stan} +#' (see \code{\link{metab_bayes_2s}}). +#' +#' The upstream observation that "matches" a given downstream observation at +#' row \code{i} was recorded \code{lag[i] <- round(travel.time[i] / +#' timestep_days)} timesteps earlier, where \code{timestep_days} is the +#' median timestep of \code{data$solar.time}. This must be computed the same +#' way as in \code{\link{mm_validate_data_2station}}'s lead-in check, so that the +#' \code{max(lag)} computed here always agrees with the lead-in requirement +#' already validated there. The first \code{max(lag)} rows of \code{data} are +#' lead-in rows: they supply upstream DO for the shift but are never +#' themselves treated as modeled (downstream) observations. +#' +#' @param data data.frame as validated by \code{\link{mm_validate_data}} for +#' \code{\link{metab_bayes_2s}}: must contain \code{solar.time}, +#' \code{DO.obs.up}, \code{DO.sat.up}, \code{DO.obs.down}, +#' \code{DO.sat.down}, \code{light}, \code{depth}, \code{temp.water}, +#' \code{travel.time}, sorted ascending by \code{solar.time}, and must +#' include the lead-in rows required to cover the longest travel time (see +#' \code{\link{metab_bayes_2s}}). +#' @param specs a list of model specs (see \code{\link{specs}}), expected to +#' already contain \code{K600_lnorm_meanlog} and \code{K600_lnorm_sdlog} +#' -- e.g., the object returned by \code{specs(mm_name('bayes_2s'))}, which +#' populates them with sensible defaults. This function does not supply +#' its own fallback values; if \code{specs} is omitted or missing these +#' fields, the resulting Stan data will contain NULL/missing values for +#' them. +#' @return a named list with all variables in the Stan model's data block: +#' \code{n_obs}, \code{n_days}, \code{DO_obs_up}, \code{DO_sat_up}, +#' \code{DO_obs_down}, \code{DO_sat_down}, \code{light}, \code{depth}, +#' \code{temp_water}, \code{travel_time} (each an \code{n_obs x n_days} +#' matrix, unitless), and \code{K600_lnorm_meanlog}/\code{K600_lnorm_sdlog} +#' @importFrom unitted v +#' @keywords internal +prepdata_bayes_2s <- function(data, specs=NULL) { + + # strip units; Stan cannot handle unitted vectors/matrices + data <- v(data) + + # timestep_days must match mm_validate_data_2station()'s lead-in check + # exactly (median timestep, in days), so that max_lag here agrees with + # what was validated there + timestep_days <- stats::median(as.numeric(diff(data$solar.time), units='days')) + + # lag, in timesteps, that the upstream series must be shifted by to line up + # with each row's downstream observation + lag <- round(data$travel.time / timestep_days) + max_lag <- max(lag) + n_total <- nrow(data) + + # the first max_lag rows are lead-in only (upstream DO used for the shift, + # but never modeled themselves). every row i > max_lag is guaranteed to + # have a valid shift target (i - lag[i] >= 1) because lag[i] <= max_lag + keep <- seq.int(max_lag + 1, n_total) + shift_idx <- keep - lag[keep] + if(any(shift_idx < 1)) { + # should be unreachable given max_lag's definition; guards against + # programming errors rather than expected user input + stop('internal error: upstream shift index falls before the first row of data') + } + + modeled <- data.frame( + solar.time = data$solar.time[keep], + DO_obs_up = data$DO.obs.up[shift_idx], + DO_sat_up = data$DO.sat.up[shift_idx], + DO_obs_down = data$DO.obs.down[keep], + DO_sat_down = data$DO.sat.down[keep], + light = data$light[keep], + depth = data$depth[keep], + temp_water = data$temp.water[keep], + travel_time = data$travel.time[keep] + ) + + # pivot into n_obs x n_days matrices, one column per unique date, using the + # same mm_time_by_date_matrix()/mm_check_dates_contiguous() helpers shared + # with prepdata_bayes() (see mm_time_by_date_matrix.R) + date_vec <- as.character(as.Date(modeled$solar.time)) + date_table <- table(date_vec) + n_days <- length(date_table) + n_obs_per_day <- unique(unname(date_table)) + if(length(n_obs_per_day) > 1) { + stop( + 'dates have differing numbers of modeled rows after lead-in removal; ', + 'observations cannot be combined into a matrix: ', + paste(sprintf('%s (%d rows)', names(date_table), date_table), collapse=', ')) + } + n_obs <- n_obs_per_day + + to_matrix <- mm_time_by_date_matrix(n_obs, n_days) + + # confirm each date occupies a contiguous block of rows, i.e., that data + # was sorted by solar.time; otherwise the matrix pivot below would silently + # scramble which rows belong to which date + mm_check_dates_contiguous( + to_matrix(date_vec), date_table, + 'data must be sorted by solar.time so that each date occupies a contiguous block of rows') + + list( + n_obs = n_obs, + n_days = n_days, + DO_obs_up = to_matrix(modeled$DO_obs_up), + DO_sat_up = to_matrix(modeled$DO_sat_up), + DO_obs_down = to_matrix(modeled$DO_obs_down), + DO_sat_down = to_matrix(modeled$DO_sat_down), + light = to_matrix(modeled$light), + depth = to_matrix(modeled$depth), + temp_water = to_matrix(modeled$temp_water), + travel_time = to_matrix(modeled$travel_time), + K600_lnorm_meanlog = specs$K600_lnorm_meanlog, + K600_lnorm_sdlog = specs$K600_lnorm_sdlog + ) +} + + +#### metab_bayes_2s class #### + +#' Metabolism model fitted by two-station (Variable Flow Two-Station) Bayesian +#' MCMC +#' +#' \code{metab_bayes_2s} models use Bayesian MCMC methods to fit values of +#' GPP, ER, and K600 from paired upstream/downstream DO curves. This class +#' inherits from \code{metab_bayes} (same \code{log}/\code{mcmc}/ +#' \code{mcmc_data}/\code{compile_time} slots, and therefore the same +#' \code{\link{get_mcmc}}, \code{\link{get_mcmc_data}}, and +#' \code{\link{get_log}} methods), but \code{predict_metab} and +#' \code{predict_DO} are overridden below because two-station's fitted-value +#' structure and output columns differ from one-station's. +#' +#' @exportClass metab_bayes_2s +#' @family metab.model.classes +setClass("metab_bayes_2s", contains="metab_bayes") + + +#' @describeIn get_params Does the same Stan-output-to-streamMetabolizer +#' renaming as \code{get_params.metab_bayes}, but does not delegate to +#' \code{get_params.metab_model} via \code{NextMethod()}: the two-station +#' steady-state model's daily GPP/ER/K600 parameters don't fit that +#' generic's one-station, ODE-based parameter-name lookup. \code{fixed} +#' column/star annotations (relevant only to models that can take fixed +#' daily parameters from \code{data_daily}) are not supported here. +#' @export +#' @import dplyr +get_params.metab_bayes_2s <- function( + metab_model, date_start=NA, date_end=NA, uncertainty=c('sd','ci','none'), messages=TRUE, ...) { + + # not delegated to get_params.metab_model via NextMethod(): that generic's + # parameter-name lookup builds an ODE-based dDOdt function from one-station's + # ode_method/GPP_fun/ER_fun/deficit_src specs, which don't exist for this + # steady-state model + uncertainty <- match.arg(uncertainty) + + fit <- metab_model@fit$daily + if(is.null(fit)) return(NULL) + + # Stan prohibits '.' in variable names, so convert back from '_' to '.', + # as in get_params.metab_bayes + parnames <- setNames(gsub('_', '\\.', metab_model@specs$params_out), metab_model@specs$params_out) + parnames <- parnames[order(nchar(parnames), decreasing=TRUE)] + for(i in seq_along(parnames)) { + names(fit) <- gsub(names(parnames[i]), parnames[[i]], names(fit)) + } + names(fit) <- gsub('_mean$', '', names(fit)) + names(fit) <- gsub('_sd$', '.sd', names(fit)) + names(fit) <- gsub('_50pct$', '.median', names(fit)) + names(fit) <- gsub('_2.5pct$', '.lower', names(fit)) + names(fit) <- gsub('_97.5pct$', '.upper', names(fit)) + + fit <- mm_filter_dates(fit, date_start=date_start, date_end=date_end) + + metab.vars <- c('GPP.daily', 'ER.daily', 'K600.daily') + for(mv in metab.vars) { + if(paste0(mv, '.median') %in% names(fit)) fit[[mv]] <- fit[[paste0(mv, '.median')]] + } + keep.cols <- c('date', unlist(lapply(metab.vars, function(mv) grep(paste0('^', mv, '($|\\.)'), names(fit), value=TRUE)))) + params <- fit[intersect(keep.cols, names(fit))] + + params <- switch( + uncertainty, + 'none' = params[!grepl('\\.median$|\\.sd$|\\.lower$|\\.upper$', names(params))], + 'sd' = params[!grepl('\\.median$|\\.lower$|\\.upper$', names(params))], + 'ci' = params[!grepl('\\.median$|\\.sd$', names(params))]) + + # attach raw warnings/errors columns (not yet compressed into a single + # column); show()'s pretty_print_ddat()/compress_msgs() does that + # compression itself at print time, as in get_params.metab_model + if(messages && exists('date', fit) && any(c('warnings','errors') %in% names(fit))) { + msgs <- fit[c('date','warnings','errors') %>% { .[. %in% names(fit)] }] + params <- left_join(params, msgs, by='date', copy=TRUE) + } + + params +} + + +#' @describeIn predict_metab Pulls daily GPP, ER, and K600 estimates out of +#' the two-station Stan model results. +#' @export +#' @import dplyr +#' @importFrom lifecycle deprecated is_present +#' @importFrom unitted get_units u +predict_metab.metab_bayes_2s <- function(metab_model, date_start=NA, date_end=NA, ..., attach.units=deprecated()) { + + Var1 <- Var2 <- '.dplyr.var' + + # check units-related arguments + if (lifecycle::is_present(attach.units)) { + unitted_deprecate_warn("predict_metab(attach.units)") + } else { + attach.units <- FALSE + } + + fit.names <- expand.grid(c('50pct','2.5pct','97.5pct'), c('GPP_daily','ER_daily','K600_daily'), stringsAsFactors=FALSE) %>% + select(Var2, Var1) %>% + apply(MARGIN=1, FUN=function(row) do.call(paste, c(as.list(row), list(sep='_')))) + metab.names <- expand.grid(c('','.lower','.upper'), c('GPP','ER','K600'), stringsAsFactors=FALSE) %>% + select(Var2, Var1) %>% + apply(MARGIN=1, FUN=function(row) do.call(paste0, as.list(row))) + + fit <- metab_model@fit$daily %>% + mm_filter_dates(date_start=date_start, date_end=date_end) + if(is.null(fit) || !all(fit.names %in% names(fit))) { + stop('could not find GPP_daily, ER_daily, and K600_daily estimates in the model fit') + } + preds <- fit[c('date', fit.names)] %>% + setNames(c('date', metab.names)) + + # add date-specific fitting warnings/errors, as in predict_metab.metab_bayes + warnings <- errors <- '.dplyr.var' + if(!is.null(fit) && all(c('date','warnings','errors') %in% names(fit))) { + messages <- fit %>% + select(date, warnings, errors) %>% + compress_msgs('msgs.fit', warnings.overall=metab_model@fit$warnings, errors.overall=metab_model@fit$errors) + preds <- full_join(preds, messages, by='date', copy=TRUE) + } else { + preds <- mutate(preds, msgs.fit=NA) + } + + preds <- mutate( + preds, + warnings=if(length(metab_model@fit$errors) > 0) NA else '', + errors=if(length(metab_model@fit$errors) > 0) NA else '') + + # attach.units if requested + if(attach.units) { + pred.units <- get_units(mm_data())[sapply(names(preds), function(x) strsplit(x, '\\.')[[1]][1], USE.NAMES=FALSE)] + preds <- u(preds, pred.units) + } + preds +} + + +#' @describeIn predict_DO Two-station (Variable Flow Two-Station) models. +#' Returns a data.frame with columns \code{solar.time}, \code{DO.obs.down} +#' (the observed downstream DO from the input data), and \code{DO.mod.down} +#' (the posterior median of the two-station Stan model's fitted downstream DO) +#' -- unlike the one-station \code{predict_DO} methods, which return +#' \code{DO.obs}/\code{DO.mod}. The values are those computed once at fitting +#' time (see \code{\link{metab_bayes_2s}}); \code{use_saved=FALSE} (on-demand +#' recomputation from the fitted daily GPP/ER/K600 medians) is not +#' implemented. +#' @export +predict_DO.metab_bayes_2s <- function(metab_model, date_start=NA, date_end=NA, ..., use_saved=TRUE) { + + if(!isTRUE(use_saved)) { + stop("predict_DO(use_saved=FALSE) is not implemented for metab_bayes_2s; only the fitted-time DO.mod.down values are available") + } + + inst <- metab_model@fit$inst + if(is.null(inst)) { + stop("no DO.mod.down predictions are available; the model fit may have failed (see get_fit(metab_model))") + } + + # NOTE: R/plot_DO_preds.R and tests/testthat/helper-rmse_DO.R both + # hard-code the one-station DO.obs/DO.mod column names; they still need to + # branch on (or be parameterized for) DO.obs.down/DO.mod.down before those + # tools will work with two-station predictions -- deferred, out of scope + # for this PR. + mm_filter_dates(inst, date_start=date_start, date_end=date_end) +} diff --git a/R/mm_data.R b/R/mm_data.R index d9ce2398..da229d62 100644 --- a/R/mm_data.R +++ b/R/mm_data.R @@ -24,6 +24,24 @@ #' equilibrium saturation \eqn{mg O[2] L^{-1}}{mg O2 / L}. Calculate using #' \link{calc_DO_sat}} #' +#' \item{ \code{DO.obs.up} dissolved oxygen concentration observations at the +#' upstream station of a two-station reach, \eqn{mg O[2] L^{-1}}{mg O2 / L}} +#' +#' \item{ \code{DO.sat.up} dissolved oxygen concentrations at equilibrium +#' saturation at the upstream station of a two-station reach, \eqn{mg O[2] +#' L^{-1}}{mg O2 / L}} +#' +#' \item{ \code{DO.obs.down} dissolved oxygen concentration observations at +#' the downstream station of a two-station reach, \eqn{mg O[2] L^{-1}}{mg O2 +#' / L}} +#' +#' \item{ \code{DO.sat.down} dissolved oxygen concentrations at equilibrium +#' saturation at the downstream station of a two-station reach, \eqn{mg O[2] +#' L^{-1}}{mg O2 / L}} +#' +#' \item{ \code{travel.time} reach travel time between the upstream and +#' downstream stations of a two-station reach, in days, \eqn{d}{d}} +#' #' \item{ \code{depth} stream depth, \eqn{m}{m}}. #' #' \item{ \code{temp.water} water temperature, \eqn{degC}}. @@ -97,9 +115,14 @@ mm_data <- function(..., optional='none') { solar.time = u(as.POSIXct("2050-03-14 15:10:00", tz="UTC"), NA), DO.obs = u(10.1,"mgO2 L^-1"), DO.sat = u(14.2,"mgO2 L^-1"), + DO.obs.up = u(10.1,"mgO2 L^-1"), + DO.sat.up = u(14.2,"mgO2 L^-1"), + DO.obs.down = u(9.8,"mgO2 L^-1"), + DO.sat.down = u(14.0,"mgO2 L^-1"), depth = u(0.5,"m"), temp.water = u(21.8,"degC"), light = u(300.9,"umol m^-2 s^-1"), + travel.time = u(0.05,"d"), discharge = u(9,"m^3 s^-1"), velocity = u(2,"m s^-1"), date = u(as.Date("2050-03-14", tz="UTC"), NA), @@ -130,6 +153,9 @@ mm_data <- function(..., optional='none') { ER = u(-5,"gO2 m^-2 d^-1"), ER.lower = u(-6,"gO2 m^-2 d^-1"), ER.upper = u(-4,"gO2 m^-2 d^-1"), + K600 = u(10,"d^-1"), + K600.lower = u(4.5,"d^-1"), + K600.upper = u(15.6,"d^-1"), D = u(5,"gO2 m^-3 d^-1"), D.lower = u(5,"gO2 m^-3 d^-1"), D.upper = u(5,"gO2 m^-3 d^-1") diff --git a/R/mm_determine_cores.R b/R/mm_determine_cores.R new file mode 100644 index 00000000..7c9ea33e --- /dev/null +++ b/R/mm_determine_cores.R @@ -0,0 +1,27 @@ +#' Determine how many cores to use for an MCMC run +#' +#' Shared by \code{runstan_bayes} (one-station) and \code{metab_bayes_2s} +#' (two-station): detects the number of cores available on the machine, +#' falls back to 1 if detection fails, and caps the requested core count at +#' whatever is actually available. +#' +#' @param n_cores the number of cores requested for this run +#' @param n_chains the number of chains being requested, used only to +#' reconstruct the verbose status message; if NULL (the default), the +#' chains clause is omitted from the message. Ignored when +#' \code{verbose=FALSE} +#' @param verbose logical. if TRUE, emit a status message reporting the +#' number of cores requested vs. available +#' @return the number of cores to actually use, i.e. +#' \code{min(detected_cores, n_cores)} +#' @keywords internal +mm_determine_cores <- function(n_cores, n_chains=NULL, verbose=FALSE) { + tot_cores <- parallel::detectCores() + if(!is.finite(tot_cores)) tot_cores <- 1 + n_cores <- min(tot_cores, n_cores) + if(verbose) { + chains_clause <- if(is.null(n_chains)) "" else paste0(n_chains," chains on ") + message(paste0("MCMC (","Stan","): requesting ",chains_clause,n_cores," of ",tot_cores," available cores")) + } + n_cores +} diff --git a/R/mm_name.R b/R/mm_name.R index cc41d30d..17079dd4 100644 --- a/R/mm_name.R +++ b/R/mm_name.R @@ -54,13 +54,16 @@ #' @param type character. The model type. Options: \itemize{ \item \code{mle}: #' maximum likelihood estimation (see also \code{\link{metab_mle}}) \item #' \code{bayes}: bayesian hierarchical models \code{\link{metab_bayes}} \item -#' \code{night}: nighttime regression (see also \code{\link{metab_night}}) -#' \item \code{Kmodel}: regression of \emph{daily} estimates of -#' \code{K600.daily} versus discharge, time, etc., usually for 3-phase -#' estimation of K alone (by MLE or nighttime regression), K vs discharge -#' (using this model), and then GPP and ER with fixed K (by MLE) (see also -#' \code{\link{metab_Kmodel}}) \item \code{sim}: simulation of \code{DO.obs} -#' 'data' for testing other models (see also \code{\link{metab_sim}}) } +#' \code{bayes_2s}: two-station (upstream/downstream, Variable Flow +#' Two-Station) Bayesian model with a single fixed structure (see also +#' \code{\link{metab_bayes_2s}}) \item \code{night}: nighttime regression (see +#' also \code{\link{metab_night}}) \item \code{Kmodel}: regression of +#' \emph{daily} estimates of \code{K600.daily} versus discharge, time, etc., +#' usually for 3-phase estimation of K alone (by MLE or nighttime regression), +#' K vs discharge (using this model), and then GPP and ER with fixed K (by +#' MLE) (see also \code{\link{metab_Kmodel}}) \item \code{sim}: simulation of +#' \code{DO.obs} 'data' for testing other models (see also +#' \code{\link{metab_sim}}) } #' @param pool_K600 character. [How] should the model pool information among #' days to get more consistent daily estimates for K600? Options (see Details #' for more): \itemize{ \item \code{none}: no pooling of K600 \item @@ -155,7 +158,7 @@ #' mm_name('sim', err_proc_acor=TRUE) #' mm_name('bayes', pool_K600='binned') mm_name <- function( - type=c('mle','bayes','night','Kmodel','sim'), + type=c('mle','bayes','bayes_2s','night','Kmodel','sim'), #pool_GPP='none', pool_ER='none', pool_eoi='alldays', pool_epc='alldays', pool_epi='alldays', pool_K600=c('none', 'normal','normal_sdzero','normal_sdfixed', @@ -174,17 +177,39 @@ mm_name <- function( deficit_src=c('DO_mod','DO_obs','DO_obs_filter','NA'), engine=c('stan','nlm','lm','mean','loess','rnorm'), check_validity=TRUE) { - - # determine type - type <- match.arg(type) - + + # determine type. 'bayes_2s' is matched exactly, before match.arg's + # partial-prefix matching, because 'b' would otherwise be an ambiguous + # abbreviation between 'bayes' and 'bayes_2s' -- so unlike the other + # types, 'bayes_2s' must be spelled out in full (no abbreviations). + # match.arg's choices are narrowed to exclude it so that pre-existing + # abbreviations like 'b' (-> 'bayes') and 'm' (-> 'mle') stay unambiguous. + if(missing(type)) { + type <- eval(formals(mm_name)$type)[1] + } else if(length(type) == 1 && identical(type, 'bayes_2s')) { + type <- 'bayes_2s' + } else { + type <- match.arg(type, choices=setdiff(eval(formals(mm_name)$type), 'bayes_2s')) + } + + # bayes_2s has a single fixed model structure rather than being built + # from combinations of pool_K600/err_*/ode_method/GPP_fun/ER_fun/ + # deficit_src/engine, so skip the argument-combination machinery below and + # return the one valid name directly + if(type == 'bayes_2s') { + mmname <- 'b2_np_oi_tr_plrckm.stan' + check_validity <- if(!is.logical(check_validity)) stop("need check_validity to be a logical of length 1") else check_validity[1] + if(isTRUE(check_validity)) mm_validate_name(mmname) + return(mmname) + } + # set type-specific defaults where values weren't specified . <- '.dplyr.var' if(type != 'Kmodel') { relevant_args <- names(formals(mm_name)) %>% .[!(. %in% c('type','check_validity'))] } else { # only one argument allowed for Kmodel - relevant_args <- 'engine' + relevant_args <- 'engine' # directly specify all the rest pool_K600='complete' pool_all='complete' @@ -205,9 +230,9 @@ mm_name <- function( assign(ms, default_args[[ms]]) } } - - # check arguments and throw errors as needed. these checks define the names - # that are possible to create; will be supplemented by call to mm_valid_names + + # check arguments and throw errors as needed. these checks define the names + # that are possible to create; will be supplemented by call to mm_valid_names # to see if a specific arg combo is actually implemented if(type != 'Kmodel') { pool_K600 <- match.arg(pool_K600) @@ -229,7 +254,7 @@ mm_name <- function( engine <- match.arg(engine) if(!(engine %in% list(bayes='stan', mle='nlm', night='lm', Kmodel=c('mean','lm','loess'), sim='rnorm')[[type]])) stop("mismatch between type (",type,") and engine (",engine,")") - + # make the name mmname <- paste0( c(bayes='b', mle='m', night='n', Kmodel='K', sim='s')[[type]], '_', @@ -237,19 +262,19 @@ mm_name <- function( c(none_or_fitted='', sdzero='0', sdfixed='x')[[tryCatch(strsplit(pool_K600, '_')[[1]][[2]], error=function(e) 'none_or_fitted')]], c(none='np', partial='', complete='')[[pool_all]], '_', if(err_obs_iid) 'oi', if(err_proc_acor) 'pc', if(err_proc_iid) 'pi', if(err_proc_GPP) 'pp', '_', - c(Euler='Eu', pairmeans='pm', trapezoid='tr', rk2='r2', - lsoda='o1', lsode='o2', lsodes='o3', lsodar='o4', vode='o5', daspk='o6', euler='eu', rk4='o8', + c(Euler='Eu', pairmeans='pm', trapezoid='tr', rk2='r2', + lsoda='o1', lsode='o2', lsodes='o3', lsodar='o4', vode='o5', daspk='o6', euler='eu', rk4='o8', ode23='o9', ode45='o10', radau='o11', bdf='o12', bdf_d='o13', adams='o14', impAdams='o15', impAdams_d='o16', 'NA'='')[[ode_method]], '_', c(linlight='pl', satlight='ps', satlightq10temp='pq', 'NA'='')[[GPP_fun]], c(constant='rc', q10temp='rq', 'NA'='')[[ER_fun]], - c(DO_mod='km', DO_obs='ko', DO_obs_filter='kf', 'NA'='')[[deficit_src]], + c(DO_mod='km', DO_obs='ko', DO_obs_filter='kf', 'NA'='')[[deficit_src]], '.', engine) - + # check validity if requested check_validity <- if(!is.logical(check_validity)) stop("need check_validity to be a logical of length 1") else check_validity[1] if(isTRUE(check_validity)) mm_validate_name(mmname) - + # return mmname } diff --git a/R/mm_parse_name.R b/R/mm_parse_name.R index 43cd18a8..d169369a 100644 --- a/R/mm_parse_name.R +++ b/R/mm_parse_name.R @@ -1,20 +1,20 @@ #' Parse a model name into its features -#' -#' Returns a data.frame with one column per model structure detail and one row -#' per `model_name` supplied to this function. See \code{?\link{mm_name}} for a +#' +#' Returns a data.frame with one column per model structure detail and one row +#' per `model_name` supplied to this function. See \code{?\link{mm_name}} for a #' description of each of the data.frame columns that is returned. -#' -#' Custom model files (for MCMC) may have additional characters after an -#' underscore at the end of the name and before the prefix. For example, +#' +#' Custom model files (for MCMC) may have additional characters after an +#' underscore at the end of the name and before the prefix. For example, #' 'b_np_pcpi_eu_ko.stan' and 'b_np_pcpi_eu_ko_v2.stan' are parsed the same; the #' _v2 is ignored by this function. -#' +#' #' @seealso The converse of this function is \code{\link{mm_name}}. -#' +#' #' @param model_name character: the model name #' @param expand logical: should additional columns such as model_name and #' pool_K600_type be added? If expand=TRUE then the result cannot be passed -#' directly back into mm_name, but the additional columns may be helpful for +#' directly back into mm_name, but the additional columns may be helpful for #' interpreting the model structure. #' @import dplyr #' @importFrom stats na.omit @@ -25,21 +25,24 @@ mm_parse_name <- function(model_name, expand=FALSE) { # define function that gets used to parse prk_terms - match_or_NA <- function(key, pairs) { - matches <- c(unname(na.omit(key[pairs]))) + match_or_NA <- function(key, pairs) { + matches <- c(unname(na.omit(key[pairs]))) if(length(matches) == 0) { 'NA' } else if(length(matches) > 1) { - stop('found too many matches in PRK terms') + stop('found too many matches in PRK terms') } else { matches } } # parse the name parsed <- strsplit(basename(model_name), "_|\\.") sapply(1:length(parsed), function(pnum) if(length(parsed[[pnum]]) <= 5) stop('missing one or more pieces in name: ', model_name[pnum])) - type <- unname(c(b='bayes', m='mle', n='night', K='Kmodel', s='sim')[sapply(parsed, `[`, 1)]) + # the "_|\\." split regex above handles 'b2' correctly. No change + # to the token-extraction logic was needed to support two-station names -- + # only a new lookup entry below, mapping the 'b2' token to its type name. + type <- unname(c(b='bayes', b2='bayes_2s', m='mle', n='night', K='Kmodel', s='sim')[sapply(parsed, `[`, 1)]) pool_K600 <- unname(c( - np='none', + np='none', Kn='normal', Kn0='normal_sdzero', Knx='normal_sdfixed', Kl='linear', Kl0='linear_sdzero', Klx='linear_sdfixed', Kb='binned', Kb0='binned_sdzero', Kbx='binned_sdfixed', @@ -59,8 +62,8 @@ mm_parse_name <- function(model_name, expand=FALSE) { err_proc_iid <- grepl('pi', sapply(parsed, `[`, 3)) err_proc_GPP <- grepl('pp', sapply(parsed, `[`, 3)) ode_method <- unname( - c(Eu='Euler', pm='pairmeans', tr='trapezoid', r2='rk2', o1='lsoda', o2='lsode', o3='lsodes', - o4='lsodar', o5='vode', o6='daspk', o7='euler', eu='euler', o8='rk4', o9='ode23', o10='ode45', o11='radau', + c(Eu='Euler', pm='pairmeans', tr='trapezoid', r2='rk2', o1='lsoda', o2='lsode', o3='lsodes', + o4='lsodar', o5='vode', o6='daspk', o7='euler', eu='euler', o8='rk4', o9='ode23', o10='ode45', o11='radau', o12='bdf', o13='bdf_d', o14='adams', o15='impAdams', o16='impAdams_d')[sapply(parsed, `[`, 4)]) prk_terms <- bind_rows(lapply(parsed, function(parsed1) { prk_term <- parsed1[5] @@ -75,7 +78,7 @@ mm_parse_name <- function(model_name, expand=FALSE) { ER_fun <- prk_terms$ER_fun deficit_src <- prk_terms$deficit_src engine <- sapply(parsed, function(vec) vec[length(vec)]) # the last one - leaves room for custom name endings before the suffix - + # combine the parsed pieces into a data.frame df <- data.frame( model_name=model_name, @@ -91,10 +94,10 @@ mm_parse_name <- function(model_name, expand=FALSE) { GPP_fun=ifelse(is.na(GPP_fun), 'NA', GPP_fun), ER_fun=ifelse(is.na(ER_fun), 'NA', ER_fun), deficit_src=ifelse(is.na(deficit_src), 'NA', deficit_src), - engine=ifelse(is.na(engine), 'NA', engine), + engine=ifelse(is.na(engine), 'NA', engine), stringsAsFactors=FALSE) - + if(!expand) df$model_name <- df$pool_K600_type <- df$pool_K600_sd <- NULL - + df } diff --git a/R/mm_time_by_date_matrix.R b/R/mm_time_by_date_matrix.R new file mode 100644 index 00000000..bbdbde3b --- /dev/null +++ b/R/mm_time_by_date_matrix.R @@ -0,0 +1,43 @@ +#' Build a closure that pivots a per-row vector into a time-by-date matrix +#' +#' Both \code{prepdata_bayes} (one-station) and \code{prepdata_bayes_2s} +#' (two-station) reshape several per-row vectors (DO, depth, light, etc.) +#' into an obs-per-day x num-days matrix, assuming that the input vector is +#' already sorted so that each date occupies a contiguous block of +#' \code{n_per_group} rows. This function returns the reshaping closure; +#' \code{\link{mm_check_dates_contiguous}} verifies the contiguous-block +#' assumption actually holds. +#' +#' @param n_per_group the number of rows per date (must be the same for +#' every date; callers are responsible for having already confirmed this) +#' @param n_groups the number of distinct dates +#' @return a function of one argument, \code{vec}, that reshapes \code{vec} +#' into an \code{n_per_group} x \code{n_groups} matrix, filling by column +#' @keywords internal +mm_time_by_date_matrix <- function(n_per_group, n_groups) { + function(vec) matrix(vec, nrow=n_per_group, ncol=n_groups, byrow=FALSE) +} + +#' Confirm that each date occupies a contiguous block of rows +#' +#' Used immediately after pivoting a date vector with the closure from +#' \code{\link{mm_time_by_date_matrix}}, to catch input that wasn't actually +#' sorted by date/time before pivoting (which would otherwise let the matrix +#' reshape silently scramble which rows belong to which date). Shared by +#' \code{prepdata_bayes} and \code{prepdata_bayes_2s}. +#' +#' @param date_mat the date-identifier vector (as used to build +#' \code{date_table}) already pivoted via +#' \code{\link{mm_time_by_date_matrix}}'s closure +#' @param date_table a table of date counts, as from \code{table(date_vec)}, +#' giving the expected date for each column of \code{date_mat} +#' @param error_message the message to pass to \code{stop()} if the dates +#' are not contiguous; left to the caller so each can keep its own wording +#' @keywords internal +mm_check_dates_contiguous <- function(date_mat, date_table, error_message) { + unique_dates_per_col <- apply(date_mat, MARGIN=2, FUN=unique) + if(is.list(unique_dates_per_col) || !isTRUE(all.equal(unname(unique_dates_per_col), names(date_table)))) { + stop(error_message) + } + invisible(TRUE) +} diff --git a/R/mm_valid_names.R b/R/mm_valid_names.R index 69a7111d..9423984e 100644 --- a/R/mm_valid_names.R +++ b/R/mm_valid_names.R @@ -10,7 +10,7 @@ #' @examples #' mm_valid_names('mle') #' @export -mm_valid_names <- function(type=c('bayes','mle','night','Kmodel','sim')) { +mm_valid_names <- function(type=c('bayes','bayes_2s','mle','night','Kmodel','sim')) { type <- match.arg(type, several.ok=TRUE) @@ -39,6 +39,11 @@ mm_valid_names <- function(type=c('bayes','mle','night','Kmodel','sim')) { mnames <- grep('^b_', dir(system.file('models', package='streamMetabolizer')), value=TRUE) favorites <- c('b_np_oipi_tr_plrckm.stan','b_np_oi_tr_plrckm.stan','b_np_pi_tr_plrckm.stan','b_np_oipp_tr_plrckm.stan') }, + bayes_2s={ + # single fixed model structure; no combinatorial name-building needed + mnames <- 'b2_np_oi_tr_plrckm.stan' + favorites <- mnames + }, mle={ opts <- expand.grid( type='mle', diff --git a/R/mm_validate_data.R b/R/mm_validate_data.R index c9fbea4f..5e6f14f7 100644 --- a/R/mm_validate_data.R +++ b/R/mm_validate_data.R @@ -59,7 +59,10 @@ mm_validate_data <- function( # missing_cols was not among the data_tests or the metab_model data were # specified without a timestamp column if('na_times' %in% data_tests) { - timecol <- grep('date|time', names(dat), value=TRUE) + # match against the known timestamp column names rather than a + # substring grep for 'date'/'time', which would also match non- + # timestamp columns such as 'travel.time' + timecol <- intersect(c('solar.time','date'), names(dat)) if(length(timecol) != 1) stop("in ", data_type, " found ", length(timecol), " possible timestamp columns", call.=FALSE) na.times <- which(is.na(dat[[timecol]])) if(length(na.times) > 0) { @@ -84,16 +87,73 @@ mm_validate_data <- function( data.units <- get_units(dat)[mismatched.units] expected.units <- get_units(expected.data)[mismatched.units] stop(paste0("unexpected units in ", data_type, ": ", paste0( - "(", 1:length(mismatched.units), ") ", + "(", 1:length(mismatched.units), ") ", names(data.units), " = ", data.units, ", expected ", expected.units, collapse="; ")), call.=FALSE) } } - + + # check travel.time bounds, if present. only two-station models supply + # this column, so this is a no-op for other model types. travel.time is + # expected in days; values outside (0, 8/24] either reflect a units + # mistake (e.g., minutes or hours rather than days) or a reach whose + # travel time exceeds the 8-hour limit needed to keep the previous + # day's light from bleeding into the following day's metabolism estimate + if('travel.time' %in% names(dat)) { + travel.time <- v(dat$travel.time) + if(any(travel.time <= 0)) { + stop('travel.time must be > 0', call.=FALSE) + } + if(any(travel.time > 8/24)) { + stop('travel.time must be <= 8/24 days (8 hours); values above this either suggest incorrect units ', + '(expected days, e.g. not minutes or hours) or a reach travel time that exceeds the 8-hour limit ', + "required to prevent the previous day's light conditions from influencing the following day's ", + 'metabolism estimate', call.=FALSE) + } + } + # return the data, whose columns may be reordered/filtered dat }) - + # return the data.frames, which may have had their columns reordered during validation and are packaged as a list return(dat_all) } + + +#' Two-station-specific data validation +#' +#' Checks the lead-in coverage requirement described in +#' \code{\link{metab_bayes_2s}}, using the (median) timestep of +#' \code{data$solar.time} to compute the required lag. Column presence, +#' timestamp validity, and travel.time bounds are expected to have already +#' been checked by \code{\link{mm_validate_data}}. +#' +#' @param data data.frame as returned by \code{\link{mm_validate_data}} for +#' \code{\link{metab_bayes_2s}}: must contain \code{solar.time} and +#' \code{travel.time}, sorted ascending by \code{solar.time}. +#' @keywords internal +mm_validate_data_2station <- function(data) { + + data_v <- v(data) + travel_time <- data_v$travel.time + solar_time <- data_v$solar.time + + # there must be enough lead-in rows of upstream DO before the first + # modeled row to cover the longest travel time in the dataset. timestep_days + # is the median observation interval, in days; max_lag is the number of + # timesteps by which upstream data must lead downstream predictions. The + # first max_lag rows of data serve only as lead-in and cannot themselves be + # modeled, so at least max_lag + 1 rows are required overall. + timestep_days <- stats::median(as.numeric(diff(solar_time), units='days')) + max_lag <- max(round(travel_time / timestep_days)) + if(nrow(data) <= max_lag) { + lead_in_needed <- max_lag - nrow(data) + 1 + stop(paste0( + 'insufficient lead-in data for upstream DO: the longest travel.time implies a lag of ', + max_lag, ' timestep(s), but only ', nrow(data), ' row(s) were supplied; ', + 'need ', lead_in_needed, ' more lead-in timestep(s) of upstream data before the first modeled row')) + } + + invisible(NULL) +} diff --git a/R/specs.R b/R/specs.R index 637c1a5e..5528a639 100644 --- a/R/specs.R +++ b/R/specs.R @@ -250,6 +250,12 @@ #' proportional to light (with noise) and is applied to GPP rather than to #' dDO/dt. #' +#' @param K600_lnorm_meanlog hyperparameter for \code{type='bayes_2s'}. +#' The mean of a lognormal prior distribution for K600_daily. +#' @param K600_lnorm_sdlog hyperparameter for \code{type='bayes_2s'}. The +#' standard deviation parameter of a lognormal prior distribution for +#' K600_daily. +#' #' @param params_in Character vector of hyperparameters to pass from the specs #' list into the data list for the MCMC run. Will be automatically generated #' during the specs() call; need only be revised if you're using a custom @@ -337,39 +343,39 @@ #' specs(mm_name(type='bayes', pool_K600='normal')) #' @export specs <- function( - + ## All or several models - + model_name = mm_name(), engine, - + # inheritParams mm_model_by_ply day_start = 4, day_end = 28, - + # inheritParams mm_is_valid_day day_tests=c('full_day', 'even_timesteps', 'complete_data', 'pos_discharge', 'pos_depth'), required_timestep=NA, - - + + ## MLE - + # initial values - init.GPP.daily = 8, + init.GPP.daily = 8, init.Pmax = 10, init.alpha = 0.0001, - init.ER.daily = -10, + init.ER.daily = -10, init.ER20 = -10, init.K600.daily = 10, - - + + ## Bayes - + # model setup split_dates, keep_mcmcs = TRUE, keep_mcmc_data = TRUE, - + # hyperparameters for non-hierarchical GPP & ER GPP_daily_mu = 3.1, GPP_daily_lower = -Inf, @@ -381,41 +387,41 @@ specs <- function( ER_daily_mu = -7.1, ER_daily_upper = Inf, ER_daily_sigma = 7.1, - + # hyperparameters for non-hierarchical K600 K600_daily_meanlog = log(12), - + # hyperparameters for hierarchical K600 - normal K600_daily_meanlog_meanlog = log(12), K600_daily_meanlog_sdlog = 1.32, - + # hyperparameters for hierarchical K600 - linear. defaults should be # reasonably constrained, not too wide lnK600_lnQ_intercept_mu = 2, lnK600_lnQ_intercept_sigma = 2.4, lnK600_lnQ_slope_mu = 0, lnK600_lnQ_slope_sigma = 0.5, - - # hyperparameters for hierarchical K600 - binned. K600_daily ~ - # lognormal(K600_daily_nodes_meanlog[lnQ_bin], + + # hyperparameters for hierarchical K600 - binned. K600_daily ~ + # lognormal(K600_daily_nodes_meanlog[lnQ_bin], # K600_daily_nodes_sdlog[lnQ_bin]) with linear interpolation among bins before - # exponentiating. nodes_meanlog and nodes_sdlog may be length b = - # length(K600_daily_lnQ_nodes) or length 1 (to be replicated to length b). - # -8:6 covers almost all points in Raymond et al. 2012 and will therefore + # exponentiating. nodes_meanlog and nodes_sdlog may be length b = + # length(K600_daily_lnQ_nodes) or length 1 (to be replicated to length b). + # -8:6 covers almost all points in Raymond et al. 2012 and will therefore # always be too broad a range for a single stream. -3:3 will catch some # streams to rivers as a first cut, though users should still modify K600_lnQ_nodes_centers = -3:3, # the x=lnQ values for the nodes K600_lnQ_nodediffs_sdlog = 0.5, # for centers 1 apart; for centers 0.2 apart, use 1/5 of this K600_lnQ_nodes_meanlog = rep(log(12), length(K600_lnQ_nodes_centers)), # distribs for the y=K600 values of the nodes K600_lnQ_nodes_sdlog = rep(1.32, length(K600_lnQ_nodes_centers)), - + # hyperparameters for any K pooling or non-pooling strategy K600_daily_sdlog = switch(mm_parse_name(model_name)$pool_K600, none=1, normal_sdfixed=0.05, NA), K600_daily_sigma = switch(mm_parse_name(model_name)$pool_K600, linear_sdfixed=10, binned_sdfixed=5, NA), K600_daily_sdlog_sigma = switch(mm_parse_name(model_name)$pool_K600, normal=0.05, NA), K600_daily_sigma_sigma = switch(mm_parse_name(model_name)$pool_K600, linear=1.2, binned=0.24, NA), # normal_sdzero, linear_sdzero, and binned_sdzero all have no parameters for this - + # hyperparameters for error terms err_obs_iid_sigma_scale = 0.03, err_proc_iid_sigma_scale = 5, @@ -423,10 +429,17 @@ specs <- function( err_proc_acor_phi_beta = 1, err_proc_acor_sigma_scale = 1, err_mult_GPP_sdlog_sigma = 1, - + + # hyperparameters for two-station (bayes_2s) K600. GPP_daily_mu, + # GPP_daily_sigma, ER_daily_mu, and ER_daily_sigma above are reused as-is. + # These are the sole source of the K600 lognormal prior defaults -- + # prepdata_bayes_2s() reads them from specs$K600_lnorm_meanlog/sdlog + K600_lnorm_meanlog = 2.484907, + K600_lnorm_sdlog = 1.0, + # vector of hyperparameters to include as MCMC data params_in, - + # inheritParams runstan_bayes params_out, n_chains = 4, @@ -435,22 +448,22 @@ specs <- function( saved_steps = 500, thin_steps = 1, verbose = FALSE, - - + + ## Kmodel - + #inheritParams prepdata_Kmodel weights = c("K600/CI"), # 'K600/CI' is argued for in stream_metab_usa issue #64 filters = c(CI.max=NA, discharge.daily.max=NA, velocity.daily.max=NA), - + #inheritParams Kmodel_allply - predictors = c("discharge.daily"), + predictors = c("discharge.daily"), transforms = c(K600='log', date=NA, velocity.daily="log", discharge.daily="log"), other_args = c(), - - + + ## Sim - + # multi-day simulation parameters. already above for bayes: # K600_lnQ_nodes_centers, K600_lnQ_nodediffs_sdlog K600_lnQ_cnode_meanlog = log(6), # distrib for the y=K600 values of the middle (or just past middle) node @@ -461,7 +474,7 @@ specs <- function( sim_Kb(K600_lnQ_nodes_centers, K600_lnQ_cnode_meanlog, K600_lnQ_cnode_sdlog, K600_lnQ_nodediffs_meanlog, K600_lnQ_nodediffs_sdlog) }, - + # daily simulation parameters discharge_daily = function(n, ...) rnorm(n, 20, 3), DO_mod_1 = NULL, @@ -471,29 +484,29 @@ specs <- function( alpha = function(n, ...) pmax(0, rnorm(n, 0.0001, 0.00002)), ER_daily = function(n, ...) pmin(0, rnorm(n, -10, 5)), ER20 = function(n, ...) pmin(0, rnorm(n, -10, 4)), - + # sub-daily simulation parameters err_obs_sigma = 0.01, err_obs_phi = 0, err_proc_sigma = 0.2, err_proc_phi = 0, err_round = NA, - + # simulation replicability sim_seed = NA - + ) { - + # make it easier to enter custom specs by creating the type-specific default if model_name %in% 'mle', etc. if(model_name %in% eval(formals(mm_name)$type)) model_name <- mm_name(type=model_name) - + # check the validity of the model_name against the list of officially accepted model names mm_validate_name(model_name) - + # parse the model_name features <- mm_parse_name(model_name, expand=TRUE) - + # collect info about the arguments required <- 'model_name' all_possible <- names(formals(specs)) @@ -501,11 +514,11 @@ specs <- function( yes_missing <- all_possible[!(all_possible %in% not_missing)] prefer_missing <- setdiff(all_possible[sapply(formals(specs), is.symbol)], 'params_out') # the arguments w/o defaults, mostly prefer_not_missing <- if(features$type == 'bayes' && features$GPP_fun == 'satlight') { - c('alpha_meanlog', 'alpha_sdlog', 'Pmax_mu', 'Pmax_sigma') + c('alpha_meanlog', 'alpha_sdlog', 'Pmax_mu', 'Pmax_sigma') } else { c() # could be made more extensive } - + # argument checks if(any(required %in% yes_missing)) stop("missing and required argument: ", paste(required[required %in% yes_missing], collapse=", ")) @@ -521,16 +534,16 @@ specs <- function( if(length(redundant) > 0) { warning("argument[s] that should usually be specified in revise() rather than specs(): ", paste(redundant, collapse=", ")) } - + # collect the defaults + directly specified arguments all_specs <- as.list(environment()) - + # copy/calculate arguments as appropriate to the model specs <- list() switch( features$type, 'bayes' = { - + # list the specs that will make it all the way to the Stan model as data all_specs$params_in <- c( switch( @@ -555,27 +568,27 @@ specs <- function( if(features$err_proc_iid) 'err_proc_iid_sigma_scale', if(features$err_proc_GPP) 'err_mult_GPP_sdlog_sigma' ) - + # list all needed arguments included <- c( # model setup 'model_name', 'engine', 'split_dates', 'keep_mcmcs', 'keep_mcmc_data', - + # date ply day_tests 'day_start', 'day_end', 'day_tests', 'required_timestep', - - # discharge binning parameters are not params_in, though they're + + # discharge binning parameters are not params_in, though they're # conceptually related and therefore colocated in formals(specs) if(features$pool_K600_type == 'binned') c('K600_lnQ_nodes_centers'), - + # params_in is both a vector of specs to include and a vector to include in specs all_specs$params_in, 'params_in', - + # inheritParams runstan_bayes - 'params_out', 'n_chains', 'n_cores', + 'params_out', 'n_chains', 'n_cores', 'burnin_steps', 'saved_steps', 'thin_steps', 'verbose' ) - + # compute some arguments if('engine' %in% yes_missing) { all_specs$engine <- features$engine @@ -584,7 +597,7 @@ specs <- function( all_specs$split_dates <- switch( features$pool_K600_type, 'none' = FALSE, # pretty sure FALSE is faster. also allows hierarchical error terms - 'normal'=, 'linear'=, 'binned' = FALSE, + 'normal'=, 'linear'=, 'binned' = FALSE, stop("unknown pool_K600; unsure how to set split_dates")) } if(features$pool_K600_type == 'binned') { @@ -605,7 +618,7 @@ specs <- function( none=c(), normal=c('K600_daily_predlog'), linear=c('K600_daily_predlog', 'lnK600_lnQ_intercept', 'lnK600_lnQ_slope'), - binned=c('K600_daily_predlog', 'lnK600_lnQ_nodes')), + binned=c('K600_daily_predlog', 'lnK600_lnQ_nodes')), if(features$pool_K600_sd == 'fitted') switch( features$pool_K600_type, @@ -616,31 +629,85 @@ specs <- function( if(features$err_proc_iid) c('err_proc_iid_sigma', 'err_proc_iid'), if(features$err_proc_GPP) c('err_proc_GPP', 'GPP_pseudo_R2')) } - + + # check for errors/inconsistencies + model_path <- tryCatch( + mm_locate_filename(model_name), + error=function(e) { + warning(e) + return(model_name) + }) + if(features$engine == "NA") + stop('engine must be specified for Bayesian models') + + }, + 'bayes_2s' = { + + # bayes_2s has a single fixed model structure (see + # inst/models/b2_np_oi_tr_plrckm.stan), so params_in/params_out are + # hardcoded here rather than built up from pool_K600/GPP_fun/ER_fun/ + # err_* toggles as in the 'bayes' case above + + # the six scalar prior hyperparameters spliced into the Stan data list + # by prepdata_bayes_2s() + all_specs$params_in <- c( + 'GPP_daily_mu', 'GPP_daily_sigma', + 'ER_daily_mu', 'ER_daily_sigma', + 'K600_lnorm_meanlog', 'K600_lnorm_sdlog') + + # list all needed arguments + included <- c( + # model setup + 'model_name', 'engine', 'split_dates', 'keep_mcmcs', 'keep_mcmc_data', + + # params_in is both a vector of specs to include and a vector to include in specs + all_specs$params_in, 'params_in', + + # inheritParams runstan_bayes + 'params_out', 'n_chains', 'n_cores', + 'burnin_steps', 'saved_steps', 'thin_steps', 'verbose' + ) + + # compute some arguments + if('engine' %in% yes_missing) { + all_specs$engine <- 'stan' + } + if('split_dates' %in% yes_missing) { + # forced FALSE: the upstream/downstream lag shift ties consecutive + # days together, so days can't be modeled independently + all_specs$split_dates <- FALSE + } + if('params_out' %in% yes_missing) { + # 'metab' (the Stan transformed-parameter matrix of modeled + # downstream DO) is included so that predict_DO() has a fitted value + # to report; without it the model would only yield daily GPP/ER/K600 + all_specs$params_out <- c('GPP_daily', 'ER_daily', 'K600_daily', 'sigma', 'metab') + } + # check for errors/inconsistencies model_path <- tryCatch( - mm_locate_filename(model_name), + mm_locate_filename(model_name), error=function(e) { warning(e) return(model_name) }) - if(features$engine == "NA") + if(features$engine == "NA") stop('engine must be specified for Bayesian models') - + }, 'mle' = { # determine which init values will be needed . <- '.dplyr.var' init.needs <- paste0('init.', get_param_names(model_name)$required) - + # list all needed arguments included <- c('model_name', 'day_start', 'day_end', 'day_tests', 'required_timestep', init.needs) - }, + }, 'night' = { # list all needed arguments included <- c('model_name', 'day_start', 'day_end', 'day_tests', 'required_timestep') - + # some different defaults for night relative to other models if('day_start' %in% yes_missing) { all_specs$day_start <- 12 @@ -651,18 +718,18 @@ specs <- function( if('day_tests' %in% yes_missing) { all_specs$day_tests <- c(day_tests, 'include_sunset') } - - }, + + }, 'Kmodel' = { # list all needed arguments included <- c( 'model_name', 'engine', 'day_start', 'day_end', 'day_tests', 'required_timestep', 'weights', 'filters', 'predictors', 'transforms', 'other_args') - + if('engine' %in% yes_missing) { all_specs$engine <- features$engine } - + # some different defaults for each engine, because no one set of defaults # makes sense for all engines #if('weights' %in% yes_missing) all_specs$weights <- c("K600/CI") # same for all, so use default as in Usage @@ -687,12 +754,12 @@ specs <- function( if('other_args' %in% yes_missing) all_specs$other_args <- list(possible_args=names(formals('loess'))[-which(names(formals('loess')) %in% c('formula','data','weights'))]) } ) - + }, 'sim' = { # determine which daily parameters will be needed par_needs <- gsub('\\.', '_', unlist(get_param_names(model_name)[c('optional','required')])) - + # list all needed arguments included <- c( 'model_name', 'day_start', 'day_end', 'day_tests', 'required_timestep', @@ -701,23 +768,23 @@ specs <- function( none=c(), normal=stop("pool_K600='normal' unavailable for now; try 'binned' instead"), linear=stop("pool_K600='linear' unavailable for now; try 'binned' instead"), # 'discharge_daily', etc. - binned=c('K600_lnQ_nodes_centers', + binned=c('K600_lnQ_nodes_centers', 'K600_lnQ_cnode_meanlog', 'K600_lnQ_cnode_sdlog', 'K600_lnQ_nodediffs_meanlog', 'K600_lnQ_nodediffs_sdlog', 'lnK600_lnQ_nodes')), par_needs, 'err_round', 'sim_seed') - + if(features$pool_K600 == 'binned') { if('K600_lnQ_nodes_centers' %in% yes_missing) # override the default, which is for 'bayes' rather than 'sim' all_specs$K600_lnQ_nodes_centers <- function(discharge.daily, ...) calc_bins(log(discharge.daily), 'width', width=0.2)$bounds } } ) - + # stop if truly irrelevant arguments were given - if(length(irrelevant <- not_missing[!(not_missing %in% included)]) > 0) + if(length(irrelevant <- not_missing[!(not_missing %in% included)]) > 0) stop("irrelevant argument: ", paste(irrelevant, collapse=", ")) - + # return just the arguments we actually need add_specs_class(all_specs[included]) - + } diff --git a/data-raw/two_station_example.R b/data-raw/two_station_example.R new file mode 100644 index 00000000..b2749ebd --- /dev/null +++ b/data-raw/two_station_example.R @@ -0,0 +1,88 @@ +# Builds data/two_station_example.rda from the VFTS (Variable Flow +# Two-Station) paper's published input data. Not run automatically as part +# of the package build/check; run +# manually (with the package root as the working directory) whenever the +# example dataset needs to be regenerated. +# + +# Download from ScienceBase: https://www.sciencebase.gov/catalog/item/6887d457d4be024722b4aae2 + +# Paper URL: https://aslopubs.onlinelibrary.wiley.com/doi/pdf/10.1002/lom3.70066 + +library(dplyr) +library(unitted) + +# --- read & filter ----------------------------------------------------- + +raw <- read.csv( + file.path('..', '2_station', 'Data', '2_VFTS_and_One-station_model_input.csv'), + stringsAsFactors=FALSE) + +# VFTS-2 is the variable-travel-time two-station run (as opposed to the +# fixed-travel-time VFTS-1/VFTS-3 runs or the one-station-only 'OS' run) and +# covers an intermediate date range (2008-03-11 to 2014-02-27) +vfts2 <- raw %>% + filter(model_run == 'VFTS-2') %>% + mutate(datetime = as.POSIXct(datetime, format='%Y-%m-%dT%H:%M:%SZ', tz='UTC')) %>% + arrange(datetime) + +# 30 consecutive days from the middle of the dataset (2011, roughly halfway +# between 2008 and 2014), avoiding the start/end edges of both the overall +# dataset and of this particular gap-free stretch of observations (verified +# separately to run gap-free from 2011-07-28 to 2011-08-30 at the source +# data's native 15-minute timestep). +modeled_start <- as.POSIXct('2011-07-31 00:00:00', tz='UTC') +modeled_end <- as.POSIXct('2011-08-29 23:45:00', tz='UTC') +timestep_days <- 15/(24*60) + +# metab_bayes_2s()'s upstream-DO lag shift (see prepdata_bayes_2s() and +# metab_bayes_2s()'s "Two-station data requirements" section) needs +# max_lag = max(round(travel.time / timestep_days)) rows of lead-in +# immediately before modeled_start -- and because prepdata_bayes_2s() trims +# max_lag rows off the *start of the whole array*, not off each calendar +# day, that lead-in window must be exactly max_lag rows (not e.g. a whole +# extra day) or the first modeled date ends up with a different row count +# than the rest, which prepdata_bayes_2s() rejects. max_lag is computed from a +# generous 2-day candidate lead-in window and then trimmed to size. +candidate_start <- modeled_start - as.difftime(2, units='days') +candidate <- vfts2 %>% filter(datetime >= candidate_start, datetime <= modeled_end) +max_lag <- max(round(candidate$travel_time / timestep_days)) +lead_in_start <- modeled_start - as.difftime(max_lag * 15, units='mins') + +vfts2_window <- vfts2 %>% + filter(datetime >= lead_in_start, datetime <= modeled_end) + +# confirm the window is gap-free at the native 15-min timestep, and that +# trimming the lead-in rows (as prepdata_bayes_2s() does) leaves exactly 30 +# modeled dates with equal row counts +stopifnot(all(abs(diff(as.numeric(vfts2_window$datetime)) - 15*60) < 1e-6)) +modeled_dates <- as.Date(vfts2_window$datetime[seq.int(max_lag+1, nrow(vfts2_window))]) +stopifnot(length(unique(table(modeled_dates))) == 1) +stopifnot(length(unique(modeled_dates)) == 30) + +# --- rename to package conventions & attach units ----------------------- + +renamed <- vfts2_window %>% + transmute( + solar.time = datetime, + DO.obs.up = upstream_DO, + DO.sat.up = upstream_DO_sat, + DO.obs.down = downstream_DO, + DO.sat.down = downstream_DO_sat, + light = light, + depth = reach_depth, + temp.water = downstream_temp, + travel.time = travel_time) +# model_run and lag are intentionally dropped (not package columns) + +template <- mm_data(solar.time, DO.obs.up, DO.sat.up, DO.obs.down, DO.sat.down, + light, depth, temp.water, travel.time) +two_station_example <- renamed +for(col in names(template)) { + two_station_example[[col]] <- u(renamed[[col]], get_units(template[[col]])) +} +two_station_example <- two_station_example[names(template)] + +# --- save ----------------------------------------------------------------- + +usethis::use_data(two_station_example, overwrite=TRUE) diff --git a/data/two_station_example.rda b/data/two_station_example.rda new file mode 100644 index 00000000..d6292ef0 Binary files /dev/null and b/data/two_station_example.rda differ diff --git a/inst/models/b2_np_oi_tr_plrckm.stan b/inst/models/b2_np_oi_tr_plrckm.stan new file mode 100644 index 00000000..bf2ce4b7 --- /dev/null +++ b/inst/models/b2_np_oi_tr_plrckm.stan @@ -0,0 +1,49 @@ +// b2_np_oi_tr_plrckm.stan + +data { + int n_obs; // number of total do observations + int n_days; // number of days + array[n_obs] vector[n_days] DO_obs_up; + array[n_obs] vector[n_days] DO_sat_up; + array[n_obs] vector[n_days] DO_obs_down; + array[n_obs] vector[n_days] DO_sat_down; + array[n_obs] vector[n_days] light; + array[n_obs] vector[n_days] depth; + array[n_obs] vector[n_days] temp_water; + array[n_obs] vector[n_days] travel_time; + real GPP_daily_mu; + real GPP_daily_sigma; + real ER_daily_mu; + real ER_daily_sigma; + real K600_lnorm_meanlog; + real K600_lnorm_sdlog; +} + +parameters { + vector[n_days] GPP_daily; + vector[n_days] ER_daily; + vector[n_days] K600_daily; + real sigma; +} + +transformed parameters { + array[n_obs] vector[n_days] metab; + array[n_obs] vector[n_days] KO2; + for (i in 1:n_obs){ + for (t in 1:n_days){ + KO2[i,t] = (K600_daily[t] / depth[i,t]) / ((600 / (1800.6 - (temp_water[i,t] * 120.1) + (3.7818 * temp_water[i,t]^2) - (0.047608 * temp_water[i,t]^3)))^-0.5); + } + metab[i] = (DO_obs_up[i] + GPP_daily .* light[i] ./ depth[i] + ER_daily .* travel_time[i] ./ depth[i] + (KO2[i] .* travel_time[i] .* (DO_sat_up[i] - DO_obs_up[i] + DO_sat_down[i]) / 2)) ./ (1 + (KO2[i] .* travel_time[i]) / 2); + } +} + +model { + GPP_daily ~ normal(GPP_daily_mu, GPP_daily_sigma); + ER_daily ~ normal(ER_daily_mu, ER_daily_sigma); + for (i in 1:n_days){ + K600_daily[i] ~ lognormal(K600_lnorm_meanlog, K600_lnorm_sdlog); + } + for (i in 1:n_obs){ + DO_obs_down[i] ~ normal(metab[i], sigma); + } +} diff --git a/man/bayes_allply.Rd b/man/bayes_allply.Rd index 46e5c258..c1c77e89 100644 --- a/man/bayes_allply.Rd +++ b/man/bayes_allply.Rd @@ -9,9 +9,11 @@ bayes_allply(data_all, data_daily_all, removed, specs) } \arguments{ \item{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})} \item{data_daily_all}{data.frame of daily priors, if appropriate to the given model_path} diff --git a/man/get_params.Rd b/man/get_params.Rd index 644e2377..69526bd1 100644 --- a/man/get_params.Rd +++ b/man/get_params.Rd @@ -1,10 +1,12 @@ % Generated by roxygen2: do not edit by hand % Please edit documentation in R/metab_model_interface.R, R/metab_Kmodel.R, -% R/metab_bayes.R, R/metab_model.get_params.R, R/metab_sim.R +% R/metab_bayes.R, R/metab_bayes_2s.R, R/metab_model.get_params.R, +% R/metab_sim.R \name{get_params} \alias{get_params} \alias{get_params.metab_Kmodel} \alias{get_params.metab_bayes} +\alias{get_params.metab_bayes_2s} \alias{get_params.metab_model} \alias{get_params.metab_sim} \title{Extract the metabolism parameters (fitted and/or fixed) from a model.} @@ -42,6 +44,15 @@ get_params( attach.units = deprecated() ) +\method{get_params}{metab_bayes_2s}( + metab_model, + date_start = NA, + date_end = NA, + uncertainty = c("sd", "ci", "none"), + messages = TRUE, + ... +) + \method{get_params}{metab_model}( metab_model, date_start = NA, @@ -109,22 +120,30 @@ parameters describing the rates and/or shapes of GPP, ER, or reaeration. } \section{Methods (by class)}{ \itemize{ -\item \code{metab_Kmodel}: Make daily re-predictions of K600.daily based on the +\item \code{get_params(metab_Kmodel)}: Make daily re-predictions of K600.daily based on the across-days model of K600.daily versus predictors. Only returns estimates for K600.daily, not any of the other daily parameters -\item \code{metab_bayes}: Does a little formatting to convert from Stan output +\item \code{get_params(metab_bayes)}: Does a little formatting to convert from Stan output to streamMetabolizer parameter names; otherwise the same as \code{get_params.metab_model} -\item \code{metab_model}: This implementation is shared by many model types +\item \code{get_params(metab_bayes_2s)}: Does the same Stan-output-to-streamMetabolizer +renaming as \code{get_params.metab_bayes}, but does not delegate to +\code{get_params.metab_model} via \code{NextMethod()}: the two-station +steady-state model's daily GPP/ER/K600 parameters don't fit that +generic's one-station, ODE-based parameter-name lookup. \code{fixed} +column/star annotations (relevant only to models that can take fixed +daily parameters from \code{data_daily}) are not supported here. -\item \code{metab_sim}: Generates new simulated values for daily parameters if +\item \code{get_params(metab_model)}: This implementation is shared by many model types + +\item \code{get_params(metab_sim)}: Generates new simulated values for daily parameters if they were described with evaluatable expressions in \code{\link{specs}}, or returns the fixed values for daily parameters if they were set in \code{data_daily} -}} +}} \examples{ dat <- data_metab('3', day_start=12, day_end=36) mm <- metab_night(specs(mm_name('night')), data=dat) @@ -135,10 +154,10 @@ get_params(mm, date_start=get_fit(mm)$date[2]) \code{\link{predict_metab}} for daily average rates of GPP and ER Other metab_model_interface: -\code{\link{get_data_daily}()}, \code{\link{get_data}()}, -\code{\link{get_fitting_time}()}, +\code{\link{get_data_daily}()}, \code{\link{get_fit}()}, +\code{\link{get_fitting_time}()}, \code{\link{get_info}()}, \code{\link{get_param_names}()}, \code{\link{get_specs}()}, diff --git a/man/metab_Kmodel-class.Rd b/man/metab_Kmodel-class.Rd index 74fe446d..4036c68f 100644 --- a/man/metab_Kmodel-class.Rd +++ b/man/metab_Kmodel-class.Rd @@ -12,6 +12,7 @@ available data to reach better, less variable daily estimates of K \seealso{ Other metab.model.classes: \code{\link{metab_bayes-class}}, +\code{\link{metab_bayes_2s-class}}, \code{\link{metab_mle-class}}, \code{\link{metab_model-class}}, \code{\link{metab_night-class}}, diff --git a/man/metab_Kmodel.Rd b/man/metab_Kmodel.Rd index f8368f60..6f2f37b5 100644 --- a/man/metab_Kmodel.Rd +++ b/man/metab_Kmodel.Rd @@ -135,6 +135,7 @@ plot_metab_preds(mm3) \seealso{ Other metab_model: \code{\link{metab_bayes}}, +\code{\link{metab_bayes_2s}}, \code{\link{metab_mle}}, \code{\link{metab_night}}, \code{\link{metab_sim}} diff --git a/man/metab_bayes-class.Rd b/man/metab_bayes-class.Rd index f74b4d3d..7435bf4e 100644 --- a/man/metab_bayes-class.Rd +++ b/man/metab_bayes-class.Rd @@ -11,6 +11,7 @@ and K for a given DO curve. \seealso{ Other metab.model.classes: \code{\link{metab_Kmodel-class}}, +\code{\link{metab_bayes_2s-class}}, \code{\link{metab_mle-class}}, \code{\link{metab_model-class}}, \code{\link{metab_night-class}}, diff --git a/man/metab_bayes.Rd b/man/metab_bayes.Rd index e7727952..ae1b3d1a 100644 --- a/man/metab_bayes.Rd +++ b/man/metab_bayes.Rd @@ -80,6 +80,7 @@ file.edit(get_specs(mm)$model_path) \seealso{ Other metab_model: \code{\link{metab_Kmodel}}, +\code{\link{metab_bayes_2s}}, \code{\link{metab_mle}}, \code{\link{metab_night}}, \code{\link{metab_sim}} diff --git a/man/metab_bayes_2s-class.Rd b/man/metab_bayes_2s-class.Rd new file mode 100644 index 00000000..9051296f --- /dev/null +++ b/man/metab_bayes_2s-class.Rd @@ -0,0 +1,27 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/metab_bayes_2s.R +\docType{class} +\name{metab_bayes_2s-class} +\alias{metab_bayes_2s-class} +\title{Metabolism model fitted by two-station (Variable Flow Two-Station) Bayesian +MCMC} +\description{ +\code{metab_bayes_2s} models use Bayesian MCMC methods to fit values of +GPP, ER, and K600 from paired upstream/downstream DO curves. This class +inherits from \code{metab_bayes} (same \code{log}/\code{mcmc}/ +\code{mcmc_data}/\code{compile_time} slots, and therefore the same +\code{\link{get_mcmc}}, \code{\link{get_mcmc_data}}, and +\code{\link{get_log}} methods), but \code{predict_metab} and +\code{predict_DO} are overridden below because two-station's fitted-value +structure and output columns differ from one-station's. +} +\seealso{ +Other metab.model.classes: +\code{\link{metab_Kmodel-class}}, +\code{\link{metab_bayes-class}}, +\code{\link{metab_mle-class}}, +\code{\link{metab_model-class}}, +\code{\link{metab_night-class}}, +\code{\link{metab_sim-class}} +} +\concept{metab.model.classes} diff --git a/man/metab_bayes_2s.Rd b/man/metab_bayes_2s.Rd new file mode 100644 index 00000000..bab4567d --- /dev/null +++ b/man/metab_bayes_2s.Rd @@ -0,0 +1,80 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/metab_bayes_2s.R +\name{metab_bayes_2s} +\alias{metab_bayes_2s} +\title{Two-station Bayesian metabolism model fitting function} +\usage{ +metab_bayes_2s( + specs = specs(mm_name("bayes_2s")), + data = mm_data(solar.time, DO.obs.up, DO.sat.up, DO.obs.down, DO.sat.down, light, + depth, temp.water, travel.time), + data_daily = mm_data(date, optional = "all"), + info = NULL +) +} +\arguments{ +\item{specs}{a list of model specifications and parameters for a model. +Although this may be specified manually (it's just a list), it is easier +and safer to use \code{\link{specs}} to generate the list, because the set +of required parameters and their defaults depends on the model given in the +\code{model_name} argument to \code{specs}. The help file for +\code{\link{specs}} lists the necessary parameters, describes them in +detail, and gives default values.} + +\item{data}{data.frame (not a tbl_df) of input data at the temporal +resolution of raw observations (unit-value). Columns must have the same +names, units, and format as the default. The solar.time column must also +have a timezone code ('tzone' attribute) of 'UTC'. See the +\strong{'Formatting \code{data}'} section below for a full description.} + +\item{data_daily}{data.frame containing inputs with a daily timestep. See the +\strong{'Formatting \code{data_daily}'} section below for a full +description.} + +\item{info}{any information, in any format, that you would like to store +within the metab_model object} +} +\value{ +A metab_bayes_2s object containing the fitted model. This object + can be inspected with the functions in the + \code{\link{metab_model_interface}} and also \code{\link{get_mcmc}}. +} +\description{ +Fits a two-station (upstream/downstream, Variable Flow Two-Station) Bayesian +model to estimate GPP, ER, and K600 from paired upstream and downstream DO, +temperature, light, and travel-time data, using the single fixed Stan model +in \code{inst/models/b2_np_oi_tr_plrckm.stan}. See \code{\link{mm_name}} to +choose a Bayesian model and \code{\link{specs}} for relevant options for the +\code{specs} argument. +} +\details{ +Unlike \code{\link{metab_bayes}}, which supports many model structures via +\code{split_dates}/\code{pool_K600}/etc., \code{metab_bayes_2s} always +fits every date jointly in a single Stan call (\code{specs$split_dates} is +forced to \code{FALSE} by \code{\link{specs}}), because the +upstream-downstream lag shift ties each date's first modeled rows to the +previous date's last rows. +} +\section{Two-station data requirements}{ + In addition to the checks + performed by \code{\link{mm_validate_data}}, \code{data$travel.time} (the + reach travel time between the upstream and downstream stations, in days) + must be strictly positive and no greater than 8/24 days (8 hours). + Values above this limit either indicate travel time was supplied in the + wrong units (e.g., minutes or hours instead of days), or reflect a reach + whose actual travel time exceeds the 8-hour limit required to prevent + the previous day's light conditions from influencing the following + day's metabolism estimate. There must also be enough lead-in + observations of upstream DO before the first row of \code{data} to + cover the longest travel time in the dataset. +} + +\seealso{ +Other metab_model: +\code{\link{metab_Kmodel}}, +\code{\link{metab_bayes}}, +\code{\link{metab_mle}}, +\code{\link{metab_night}}, +\code{\link{metab_sim}} +} +\concept{metab_model} diff --git a/man/metab_mle-class.Rd b/man/metab_mle-class.Rd index 62f2301c..4017dc1c 100644 --- a/man/metab_mle-class.Rd +++ b/man/metab_mle-class.Rd @@ -12,6 +12,7 @@ likelihood to fit values of GPP, ER, and K for a given DO curve. Other metab.model.classes: \code{\link{metab_Kmodel-class}}, \code{\link{metab_bayes-class}}, +\code{\link{metab_bayes_2s-class}}, \code{\link{metab_model-class}}, \code{\link{metab_night-class}}, \code{\link{metab_sim-class}} diff --git a/man/metab_mle.Rd b/man/metab_mle.Rd index 1c9ab901..de2211d5 100644 --- a/man/metab_mle.Rd +++ b/man/metab_mle.Rd @@ -72,6 +72,7 @@ plot_DO_preds(predict_DO(mm)) Other metab_model: \code{\link{metab_Kmodel}}, \code{\link{metab_bayes}}, +\code{\link{metab_bayes_2s}}, \code{\link{metab_night}}, \code{\link{metab_sim}} } diff --git a/man/metab_model-class.Rd b/man/metab_model-class.Rd index e3070d9e..7103041d 100644 --- a/man/metab_model-class.Rd +++ b/man/metab_model-class.Rd @@ -34,6 +34,7 @@ function.} Other metab.model.classes: \code{\link{metab_Kmodel-class}}, \code{\link{metab_bayes-class}}, +\code{\link{metab_bayes_2s-class}}, \code{\link{metab_mle-class}}, \code{\link{metab_night-class}}, \code{\link{metab_sim-class}} diff --git a/man/metab_night-class.Rd b/man/metab_night-class.Rd index dcf06605..362b14d9 100644 --- a/man/metab_night-class.Rd +++ b/man/metab_night-class.Rd @@ -12,6 +12,7 @@ given DO time series. Other metab.model.classes: \code{\link{metab_Kmodel-class}}, \code{\link{metab_bayes-class}}, +\code{\link{metab_bayes_2s-class}}, \code{\link{metab_mle-class}}, \code{\link{metab_model-class}}, \code{\link{metab_sim-class}} diff --git a/man/metab_night.Rd b/man/metab_night.Rd index e2fa1581..8aa9db3e 100644 --- a/man/metab_night.Rd +++ b/man/metab_night.Rd @@ -57,6 +57,7 @@ plot_DO_preds(predict_DO(mm)) Other metab_model: \code{\link{metab_Kmodel}}, \code{\link{metab_bayes}}, +\code{\link{metab_bayes_2s}}, \code{\link{metab_mle}}, \code{\link{metab_sim}} } diff --git a/man/metab_sim-class.Rd b/man/metab_sim-class.Rd index 59c59225..eb745606 100644 --- a/man/metab_sim-class.Rd +++ b/man/metab_sim-class.Rd @@ -12,6 +12,7 @@ including GPP, ER, and K600 values Other metab.model.classes: \code{\link{metab_Kmodel-class}}, \code{\link{metab_bayes-class}}, +\code{\link{metab_bayes_2s-class}}, \code{\link{metab_mle-class}}, \code{\link{metab_model-class}}, \code{\link{metab_night-class}} diff --git a/man/metab_sim.Rd b/man/metab_sim.Rd index fa6e0bf4..3dde3a6d 100644 --- a/man/metab_sim.Rd +++ b/man/metab_sim.Rd @@ -108,6 +108,7 @@ library(ggplot2) Other metab_model: \code{\link{metab_Kmodel}}, \code{\link{metab_bayes}}, +\code{\link{metab_bayes_2s}}, \code{\link{metab_mle}}, \code{\link{metab_night}} } diff --git a/man/mm_check_dates_contiguous.Rd b/man/mm_check_dates_contiguous.Rd new file mode 100644 index 00000000..8f2821b1 --- /dev/null +++ b/man/mm_check_dates_contiguous.Rd @@ -0,0 +1,27 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/mm_time_by_date_matrix.R +\name{mm_check_dates_contiguous} +\alias{mm_check_dates_contiguous} +\title{Confirm that each date occupies a contiguous block of rows} +\usage{ +mm_check_dates_contiguous(date_mat, date_table, error_message) +} +\arguments{ +\item{date_mat}{the date-identifier vector (as used to build +\code{date_table}) already pivoted via +\code{\link{mm_time_by_date_matrix}}'s closure} + +\item{date_table}{a table of date counts, as from \code{table(date_vec)}, +giving the expected date for each column of \code{date_mat}} + +\item{error_message}{the message to pass to \code{stop()} if the dates +are not contiguous; left to the caller so each can keep its own wording} +} +\description{ +Used immediately after pivoting a date vector with the closure from +\code{\link{mm_time_by_date_matrix}}, to catch input that wasn't actually +sorted by date/time before pivoting (which would otherwise let the matrix +reshape silently scramble which rows belong to which date). Shared by +\code{prepdata_bayes} and \code{prepdata_bayes_2s}. +} +\keyword{internal} diff --git a/man/mm_data.Rd b/man/mm_data.Rd index 7fe0459c..af079a96 100644 --- a/man/mm_data.Rd +++ b/man/mm_data.Rd @@ -43,6 +43,24 @@ Produces a unitted data.frame with the column names, units, and equilibrium saturation \eqn{mg O[2] L^{-1}}{mg O2 / L}. Calculate using \link{calc_DO_sat}} + \item{ \code{DO.obs.up} dissolved oxygen concentration observations at the + upstream station of a two-station reach, \eqn{mg O[2] L^{-1}}{mg O2 / L}} + + \item{ \code{DO.sat.up} dissolved oxygen concentrations at equilibrium + saturation at the upstream station of a two-station reach, \eqn{mg O[2] + L^{-1}}{mg O2 / L}} + + \item{ \code{DO.obs.down} dissolved oxygen concentration observations at + the downstream station of a two-station reach, \eqn{mg O[2] L^{-1}}{mg O2 + / L}} + + \item{ \code{DO.sat.down} dissolved oxygen concentrations at equilibrium + saturation at the downstream station of a two-station reach, \eqn{mg O[2] + L^{-1}}{mg O2 / L}} + + \item{ \code{travel.time} reach travel time between the upstream and + downstream stations of a two-station reach, in days, \eqn{d}{d}} + \item{ \code{depth} stream depth, \eqn{m}{m}}. \item{ \code{temp.water} water temperature, \eqn{degC}}. diff --git a/man/mm_determine_cores.Rd b/man/mm_determine_cores.Rd new file mode 100644 index 00000000..cedda661 --- /dev/null +++ b/man/mm_determine_cores.Rd @@ -0,0 +1,30 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/mm_determine_cores.R +\name{mm_determine_cores} +\alias{mm_determine_cores} +\title{Determine how many cores to use for an MCMC run} +\usage{ +mm_determine_cores(n_cores, n_chains = NULL, verbose = FALSE) +} +\arguments{ +\item{n_cores}{the number of cores requested for this run} + +\item{n_chains}{the number of chains being requested, used only to +reconstruct the verbose status message; if NULL (the default), the +chains clause is omitted from the message. Ignored when +\code{verbose=FALSE}} + +\item{verbose}{logical. if TRUE, emit a status message reporting the +number of cores requested vs. available} +} +\value{ +the number of cores to actually use, i.e. + \code{min(detected_cores, n_cores)} +} +\description{ +Shared by \code{runstan_bayes} (one-station) and \code{metab_bayes_2s} +(two-station): detects the number of cores available on the machine, +falls back to 1 if detection fails, and caps the requested core count at +whatever is actually available. +} +\keyword{internal} diff --git a/man/mm_generate_mcmc_file.Rd b/man/mm_generate_mcmc_file.Rd index 144d9249..1884b17c 100644 --- a/man/mm_generate_mcmc_file.Rd +++ b/man/mm_generate_mcmc_file.Rd @@ -23,13 +23,16 @@ mm_generate_mcmc_file( \item{type}{character. The model type. Options: \itemize{ \item \code{mle}: maximum likelihood estimation (see also \code{\link{metab_mle}}) \item \code{bayes}: bayesian hierarchical models \code{\link{metab_bayes}} \item -\code{night}: nighttime regression (see also \code{\link{metab_night}}) -\item \code{Kmodel}: regression of \emph{daily} estimates of -\code{K600.daily} versus discharge, time, etc., usually for 3-phase -estimation of K alone (by MLE or nighttime regression), K vs discharge -(using this model), and then GPP and ER with fixed K (by MLE) (see also -\code{\link{metab_Kmodel}}) \item \code{sim}: simulation of \code{DO.obs} -'data' for testing other models (see also \code{\link{metab_sim}}) }} +\code{bayes_2s}: two-station (upstream/downstream, Variable Flow +Two-Station) Bayesian model with a single fixed structure (see also +\code{\link{metab_bayes_2s}}) \item \code{night}: nighttime regression (see +also \code{\link{metab_night}}) \item \code{Kmodel}: regression of +\emph{daily} estimates of \code{K600.daily} versus discharge, time, etc., +usually for 3-phase estimation of K alone (by MLE or nighttime regression), +K vs discharge (using this model), and then GPP and ER with fixed K (by +MLE) (see also \code{\link{metab_Kmodel}}) \item \code{sim}: simulation of +\code{DO.obs} 'data' for testing other models (see also +\code{\link{metab_sim}}) }} \item{pool_K600}{character. [How] should the model pool information among days to get more consistent daily estimates for K600? Options (see Details diff --git a/man/mm_name.Rd b/man/mm_name.Rd index abf4a4df..c1e8037e 100644 --- a/man/mm_name.Rd +++ b/man/mm_name.Rd @@ -5,7 +5,7 @@ \title{Find the name of a model by its features} \usage{ mm_name( - type = c("mle", "bayes", "night", "Kmodel", "sim"), + type = c("mle", "bayes", "bayes_2s", "night", "Kmodel", "sim"), pool_K600 = c("none", "normal", "normal_sdzero", "normal_sdfixed", "linear", "linear_sdzero", "linear_sdfixed", "binned", "binned_sdzero", "binned_sdfixed", "complete"), @@ -27,13 +27,16 @@ mm_name( \item{type}{character. The model type. Options: \itemize{ \item \code{mle}: maximum likelihood estimation (see also \code{\link{metab_mle}}) \item \code{bayes}: bayesian hierarchical models \code{\link{metab_bayes}} \item -\code{night}: nighttime regression (see also \code{\link{metab_night}}) -\item \code{Kmodel}: regression of \emph{daily} estimates of -\code{K600.daily} versus discharge, time, etc., usually for 3-phase -estimation of K alone (by MLE or nighttime regression), K vs discharge -(using this model), and then GPP and ER with fixed K (by MLE) (see also -\code{\link{metab_Kmodel}}) \item \code{sim}: simulation of \code{DO.obs} -'data' for testing other models (see also \code{\link{metab_sim}}) }} +\code{bayes_2s}: two-station (upstream/downstream, Variable Flow +Two-Station) Bayesian model with a single fixed structure (see also +\code{\link{metab_bayes_2s}}) \item \code{night}: nighttime regression (see +also \code{\link{metab_night}}) \item \code{Kmodel}: regression of +\emph{daily} estimates of \code{K600.daily} versus discharge, time, etc., +usually for 3-phase estimation of K alone (by MLE or nighttime regression), +K vs discharge (using this model), and then GPP and ER with fixed K (by +MLE) (see also \code{\link{metab_Kmodel}}) \item \code{sim}: simulation of +\code{DO.obs} 'data' for testing other models (see also +\code{\link{metab_sim}}) }} \item{pool_K600}{character. [How] should the model pool information among days to get more consistent daily estimates for K600? Options (see Details diff --git a/man/mm_parse_name.Rd b/man/mm_parse_name.Rd index 670b22f5..cc7a3f52 100644 --- a/man/mm_parse_name.Rd +++ b/man/mm_parse_name.Rd @@ -11,17 +11,17 @@ mm_parse_name(model_name, expand = FALSE) \item{expand}{logical: should additional columns such as model_name and pool_K600_type be added? If expand=TRUE then the result cannot be passed -directly back into mm_name, but the additional columns may be helpful for +directly back into mm_name, but the additional columns may be helpful for interpreting the model structure.} } \description{ -Returns a data.frame with one column per model structure detail and one row -per `model_name` supplied to this function. See \code{?\link{mm_name}} for a +Returns a data.frame with one column per model structure detail and one row +per `model_name` supplied to this function. See \code{?\link{mm_name}} for a description of each of the data.frame columns that is returned. } \details{ -Custom model files (for MCMC) may have additional characters after an -underscore at the end of the name and before the prefix. For example, +Custom model files (for MCMC) may have additional characters after an +underscore at the end of the name and before the prefix. For example, 'b_np_pcpi_eu_ko.stan' and 'b_np_pcpi_eu_ko_v2.stan' are parsed the same; the _v2 is ignored by this function. } diff --git a/man/mm_time_by_date_matrix.Rd b/man/mm_time_by_date_matrix.Rd new file mode 100644 index 00000000..762b5cec --- /dev/null +++ b/man/mm_time_by_date_matrix.Rd @@ -0,0 +1,28 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/mm_time_by_date_matrix.R +\name{mm_time_by_date_matrix} +\alias{mm_time_by_date_matrix} +\title{Build a closure that pivots a per-row vector into a time-by-date matrix} +\usage{ +mm_time_by_date_matrix(n_per_group, n_groups) +} +\arguments{ +\item{n_per_group}{the number of rows per date (must be the same for +every date; callers are responsible for having already confirmed this)} + +\item{n_groups}{the number of distinct dates} +} +\value{ +a function of one argument, \code{vec}, that reshapes \code{vec} + into an \code{n_per_group} x \code{n_groups} matrix, filling by column +} +\description{ +Both \code{prepdata_bayes} (one-station) and \code{prepdata_bayes_2s} +(two-station) reshape several per-row vectors (DO, depth, light, etc.) +into an obs-per-day x num-days matrix, assuming that the input vector is +already sorted so that each date occupies a contiguous block of +\code{n_per_group} rows. This function returns the reshaping closure; +\code{\link{mm_check_dates_contiguous}} verifies the contiguous-block +assumption actually holds. +} +\keyword{internal} diff --git a/man/mm_valid_names.Rd b/man/mm_valid_names.Rd index 0982a1e8..fb578cd7 100644 --- a/man/mm_valid_names.Rd +++ b/man/mm_valid_names.Rd @@ -4,19 +4,22 @@ \alias{mm_valid_names} \title{Get the valid names for a given model type or types} \usage{ -mm_valid_names(type = c("bayes", "mle", "night", "Kmodel", "sim")) +mm_valid_names(type = c("bayes", "bayes_2s", "mle", "night", "Kmodel", "sim")) } \arguments{ \item{type}{character. The model type. Options: \itemize{ \item \code{mle}: maximum likelihood estimation (see also \code{\link{metab_mle}}) \item \code{bayes}: bayesian hierarchical models \code{\link{metab_bayes}} \item -\code{night}: nighttime regression (see also \code{\link{metab_night}}) -\item \code{Kmodel}: regression of \emph{daily} estimates of -\code{K600.daily} versus discharge, time, etc., usually for 3-phase -estimation of K alone (by MLE or nighttime regression), K vs discharge -(using this model), and then GPP and ER with fixed K (by MLE) (see also -\code{\link{metab_Kmodel}}) \item \code{sim}: simulation of \code{DO.obs} -'data' for testing other models (see also \code{\link{metab_sim}}) }} +\code{bayes_2s}: two-station (upstream/downstream, Variable Flow +Two-Station) Bayesian model with a single fixed structure (see also +\code{\link{metab_bayes_2s}}) \item \code{night}: nighttime regression (see +also \code{\link{metab_night}}) \item \code{Kmodel}: regression of +\emph{daily} estimates of \code{K600.daily} versus discharge, time, etc., +usually for 3-phase estimation of K alone (by MLE or nighttime regression), +K vs discharge (using this model), and then GPP and ER with fixed K (by +MLE) (see also \code{\link{metab_Kmodel}}) \item \code{sim}: simulation of +\code{DO.obs} 'data' for testing other models (see also +\code{\link{metab_sim}}) }} } \description{ Returns a vector of the \code{model_name}s for the type[s] indicated. If diff --git a/man/mm_validate_data_2station.Rd b/man/mm_validate_data_2station.Rd new file mode 100644 index 00000000..9b6c0371 --- /dev/null +++ b/man/mm_validate_data_2station.Rd @@ -0,0 +1,21 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/mm_validate_data.R +\name{mm_validate_data_2station} +\alias{mm_validate_data_2station} +\title{Two-station-specific data validation} +\usage{ +mm_validate_data_2station(data) +} +\arguments{ +\item{data}{data.frame as returned by \code{\link{mm_validate_data}} for +\code{\link{metab_bayes_2s}}: must contain \code{solar.time} and +\code{travel.time}, sorted ascending by \code{solar.time}.} +} +\description{ +Checks the lead-in coverage requirement described in +\code{\link{metab_bayes_2s}}, using the (median) timestep of +\code{data$solar.time} to compute the required lag. Column presence, +timestamp validity, and travel.time bounds are expected to have already +been checked by \code{\link{mm_validate_data}}. +} +\keyword{internal} diff --git a/man/predict_DO.Rd b/man/predict_DO.Rd index e7ba5e82..db2c026a 100644 --- a/man/predict_DO.Rd +++ b/man/predict_DO.Rd @@ -1,9 +1,11 @@ % Generated by roxygen2: do not edit by hand % Please edit documentation in R/metab_model_interface.R, R/metab_Kmodel.R, -% R/metab_model.predict_DO.R, R/metab_night.R, R/metab_sim.R +% R/metab_bayes_2s.R, R/metab_model.predict_DO.R, R/metab_night.R, +% R/metab_sim.R \name{predict_DO} \alias{predict_DO} \alias{predict_DO.metab_Kmodel} +\alias{predict_DO.metab_bayes_2s} \alias{predict_DO.metab_model} \alias{predict_DO.metab_night} \alias{predict_DO.metab_sim} @@ -20,6 +22,8 @@ predict_DO( \method{predict_DO}{metab_Kmodel}(metab_model, date_start = NA, date_end = NA, ..., use_saved = TRUE) +\method{predict_DO}{metab_bayes_2s}(metab_model, date_start = NA, date_end = NA, ..., use_saved = TRUE) + \method{predict_DO}{metab_model}( metab_model, date_start = NA, @@ -64,22 +68,32 @@ oxygen. } \section{Methods (by class)}{ \itemize{ -\item \code{metab_Kmodel}: Throws an error because models of type 'Kmodel' can't +\item \code{predict_DO(metab_Kmodel)}: Throws an error because models of type 'Kmodel' can't predict DO. \code{metab_Kmodel} predicts K at daily timesteps and usually knows nothing about GPP or ER. So it's not possible to predict DO from this model. Try passing the output to metab_mle and THEN predicting DO. -\item \code{metab_model}: This implementation is shared by many model types +\item \code{predict_DO(metab_bayes_2s)}: Two-station (Variable Flow Two-Station) models. +Returns a data.frame with columns \code{solar.time}, \code{DO.obs.down} +(the observed downstream DO from the input data), and \code{DO.mod.down} +(the posterior median of the two-station Stan model's fitted downstream DO) +-- unlike the one-station \code{predict_DO} methods, which return +\code{DO.obs}/\code{DO.mod}. The values are those computed once at fitting +time (see \code{\link{metab_bayes_2s}}); \code{use_saved=FALSE} (on-demand +recomputation from the fitted daily GPP/ER/K600 medians) is not +implemented. + +\item \code{predict_DO(metab_model)}: This implementation is shared by many model types -\item \code{metab_night}: Generate nighttime dissolved oxygen predictions from a +\item \code{predict_DO(metab_night)}: Generate nighttime dissolved oxygen predictions from a nighttime regression model. \code{metab_night} only fits ER and K, and only for the darkness hours, so predictions are only generated for those hours. -\item \code{metab_sim}: Simulate values for DO.obs (with process and +\item \code{predict_DO(metab_sim)}: Simulate values for DO.obs (with process and observation error), DO.mod (with process error only), and DO.pure (with no error). The errors are randomly generated on every new call to predict_DO. -}} +}} \examples{ dat <- data_metab('3', day_start=12, day_end=36) mm <- metab_night(specs(mm_name('night')), data=dat) @@ -88,10 +102,10 @@ head(preds) } \seealso{ Other metab_model_interface: -\code{\link{get_data_daily}()}, \code{\link{get_data}()}, -\code{\link{get_fitting_time}()}, +\code{\link{get_data_daily}()}, \code{\link{get_fit}()}, +\code{\link{get_fitting_time}()}, \code{\link{get_info}()}, \code{\link{get_param_names}()}, \code{\link{get_params}()}, diff --git a/man/predict_metab.Rd b/man/predict_metab.Rd index ca653ff7..e4a5df9e 100644 --- a/man/predict_metab.Rd +++ b/man/predict_metab.Rd @@ -1,9 +1,10 @@ % Generated by roxygen2: do not edit by hand % Please edit documentation in R/metab_model_interface.R, R/metab_bayes.R, -% R/metab_model.predict_metab.R +% R/metab_bayes_2s.R, R/metab_model.predict_metab.R \name{predict_metab} \alias{predict_metab} \alias{predict_metab.metab_bayes} +\alias{predict_metab.metab_bayes_2s} \alias{predict_metab.metab_model} \title{Predict metabolism from a fitted model.} \usage{ @@ -26,6 +27,14 @@ predict_metab( attach.units = deprecated() ) +\method{predict_metab}{metab_bayes_2s}( + metab_model, + date_start = NA, + date_end = NA, + ..., + attach.units = deprecated() +) + \method{predict_metab}{metab_model}( metab_model, date_start = NA, @@ -91,16 +100,19 @@ GPP, ER, and K600. } \section{Methods (by class)}{ \itemize{ -\item \code{metab_bayes}: Pulls daily metabolism estimates out of the Stan +\item \code{predict_metab(metab_bayes)}: Pulls daily metabolism estimates out of the Stan model results; looks for \code{GPP} or \code{GPP_daily} and for \code{ER} or \code{ER_daily} among the \code{params_out} (see \code{\link{specs}}), which means you can save just one (or both) of those sets of daily parameters when running the Stan model. Saving fewer parameters can help models run faster and use less RAM. -\item \code{metab_model}: This implementation is shared by many model types -}} +\item \code{predict_metab(metab_bayes_2s)}: Pulls daily GPP, ER, and K600 estimates out of +the two-station Stan model results. + +\item \code{predict_metab(metab_model)}: This implementation is shared by many model types +}} \examples{ dat <- data_metab('3', day_start=12, day_end=36) mm <- metab_night(specs(mm_name('night')), data=dat) @@ -109,10 +121,10 @@ predict_metab(mm, date_start=get_fit(mm)$date[2]) } \seealso{ Other metab_model_interface: -\code{\link{get_data_daily}()}, \code{\link{get_data}()}, -\code{\link{get_fitting_time}()}, +\code{\link{get_data_daily}()}, \code{\link{get_fit}()}, +\code{\link{get_fitting_time}()}, \code{\link{get_info}()}, \code{\link{get_param_names}()}, \code{\link{get_params}()}, diff --git a/man/prepdata_bayes_2s.Rd b/man/prepdata_bayes_2s.Rd new file mode 100644 index 00000000..85557198 --- /dev/null +++ b/man/prepdata_bayes_2s.Rd @@ -0,0 +1,51 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/metab_bayes_2s.R +\name{prepdata_bayes_2s} +\alias{prepdata_bayes_2s} +\title{Reshape long-format two-station data into the list expected by the +two-station Stan model} +\usage{ +prepdata_bayes_2s(data, specs = NULL) +} +\arguments{ +\item{data}{data.frame as validated by \code{\link{mm_validate_data}} for +\code{\link{metab_bayes_2s}}: must contain \code{solar.time}, +\code{DO.obs.up}, \code{DO.sat.up}, \code{DO.obs.down}, +\code{DO.sat.down}, \code{light}, \code{depth}, \code{temp.water}, +\code{travel.time}, sorted ascending by \code{solar.time}, and must +include the lead-in rows required to cover the longest travel time (see +\code{\link{metab_bayes_2s}}).} + +\item{specs}{a list of model specs (see \code{\link{specs}}), expected to +already contain \code{K600_lnorm_meanlog} and \code{K600_lnorm_sdlog} +-- e.g., the object returned by \code{specs(mm_name('bayes_2s'))}, which +populates them with sensible defaults. This function does not supply +its own fallback values; if \code{specs} is omitted or missing these +fields, the resulting Stan data will contain NULL/missing values for +them.} +} +\value{ +a named list with all variables in the Stan model's data block: + \code{n_obs}, \code{n_days}, \code{DO_obs_up}, \code{DO_sat_up}, + \code{DO_obs_down}, \code{DO_sat_down}, \code{light}, \code{depth}, + \code{temp_water}, \code{travel_time} (each an \code{n_obs x n_days} + matrix, unitless), and \code{K600_lnorm_meanlog}/\code{K600_lnorm_sdlog} +} +\description{ +Time-shifts the upstream DO series to match the travel time between +stations, then pivots the result into the \code{n_obs x n_days} matrices +expected by the \code{data} block of \code{inst/models/b2_np_oi_tr_plrckm.stan} +(see \code{\link{metab_bayes_2s}}). +} +\details{ +The upstream observation that "matches" a given downstream observation at +row \code{i} was recorded \code{lag[i] <- round(travel.time[i] / +timestep_days)} timesteps earlier, where \code{timestep_days} is the +median timestep of \code{data$solar.time}. This must be computed the same +way as in \code{\link{mm_validate_data_2station}}'s lead-in check, so that the +\code{max(lag)} computed here always agrees with the lead-in requirement +already validated there. The first \code{max(lag)} rows of \code{data} are +lead-in rows: they supply upstream DO for the shift but are never +themselves treated as modeled (downstream) observations. +} +\keyword{internal} diff --git a/man/specs.Rd b/man/specs.Rd index d0e64fa7..58e65423 100644 --- a/man/specs.Rd +++ b/man/specs.Rd @@ -42,12 +42,11 @@ specs( K600_lnQ_nodediffs_sdlog = 0.5, K600_lnQ_nodes_meanlog = rep(log(12), length(K600_lnQ_nodes_centers)), K600_lnQ_nodes_sdlog = rep(1.32, length(K600_lnQ_nodes_centers)), - K600_daily_sdlog = switch(mm_parse_name(model_name)$pool_K600, none = 1, - normal_sdfixed = 0.05, NA), + K600_daily_sdlog = switch(mm_parse_name(model_name)$pool_K600, none = 1, normal_sdfixed + = 0.05, NA), K600_daily_sigma = switch(mm_parse_name(model_name)$pool_K600, linear_sdfixed = 10, binned_sdfixed = 5, NA), - K600_daily_sdlog_sigma = switch(mm_parse_name(model_name)$pool_K600, normal = 0.05, - NA), + K600_daily_sdlog_sigma = switch(mm_parse_name(model_name)$pool_K600, normal = 0.05, NA), K600_daily_sigma_sigma = switch(mm_parse_name(model_name)$pool_K600, linear = 1.2, binned = 0.24, NA), err_obs_iid_sigma_scale = 0.03, @@ -56,6 +55,8 @@ specs( err_proc_acor_phi_beta = 1, err_proc_acor_sigma_scale = 1, err_mult_GPP_sdlog_sigma = 1, + K600_lnorm_meanlog = 2.484907, + K600_lnorm_sdlog = 1, params_in, params_out, n_chains = 4, @@ -74,9 +75,11 @@ specs( K600_lnQ_cnode_sdlog = 1, K600_lnQ_nodediffs_meanlog = 0.2, lnK600_lnQ_nodes = function(K600_lnQ_nodes_centers, K600_lnQ_cnode_meanlog, - K600_lnQ_cnode_sdlog, K600_lnQ_nodediffs_meanlog, K600_lnQ_nodediffs_sdlog, ...) { - sim_Kb(K600_lnQ_nodes_centers, K600_lnQ_cnode_meanlog, K600_lnQ_cnode_sdlog, - K600_lnQ_nodediffs_meanlog, K600_lnQ_nodediffs_sdlog) }, + K600_lnQ_cnode_sdlog, K600_lnQ_nodediffs_meanlog, K600_lnQ_nodediffs_sdlog, ...) { + + sim_Kb(K600_lnQ_nodes_centers, K600_lnQ_cnode_meanlog, K600_lnQ_cnode_sdlog, + K600_lnQ_nodediffs_meanlog, K600_lnQ_nodediffs_sdlog) + }, discharge_daily = function(n, ...) rnorm(n, 20, 3), DO_mod_1 = NULL, K600_daily = function(n, K600_daily_predlog = log(10), ...) pmax(0, rnorm(n, @@ -337,6 +340,13 @@ estimate GPP_inst. The effect is a special kind of process error that is proportional to light (with noise) and is applied to GPP rather than to dDO/dt.} +\item{K600_lnorm_meanlog}{hyperparameter for \code{type='bayes_2s'}. +The mean of a lognormal prior distribution for K600_daily.} + +\item{K600_lnorm_sdlog}{hyperparameter for \code{type='bayes_2s'}. The +standard deviation parameter of a lognormal prior distribution for +K600_daily.} + \item{params_in}{Character vector of hyperparameters to pass from the specs list into the data list for the MCMC run. Will be automatically generated during the specs() call; need only be revised if you're using a custom diff --git a/man/streamMetabolizer.Rd b/man/streamMetabolizer.Rd index 25a0bf0f..dd330519 100644 --- a/man/streamMetabolizer.Rd +++ b/man/streamMetabolizer.Rd @@ -106,3 +106,23 @@ See http://usgs-r.github.io/streamMetabolizer for vignettes on the web. } } +\seealso{ +Useful links: +\itemize{ + \item \url{https://github.com/USGS-R/streamMetabolizer} + \item \url{http://usgs-r.github.io/streamMetabolizer/} + \item Report bugs at \url{https://github.com/USGS-R/streamMetabolizer/issues} +} + +} +\author{ +\strong{Maintainer}: Alison P. Appling \email{aappling@usgs.gov} + +Authors: +\itemize{ + \item Robert O. Hall + \item Maite Arroita + \item Charles B. Yackulic +} + +} diff --git a/man/two_station_example.Rd b/man/two_station_example.Rd new file mode 100644 index 00000000..7bcafb24 --- /dev/null +++ b/man/two_station_example.Rd @@ -0,0 +1,55 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/data.R +\docType{data} +\name{two_station_example} +\alias{two_station_example} +\title{Example two-station (VFTS) input data} +\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} +} +\usage{ +two_station_example +} +\description{ +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. +} +\keyword{datasets} diff --git a/tests/testthat/test-metab_bayes.R b/tests/testthat/test-metab_bayes.R index 929872e5..d5b527b6 100644 --- a/tests/testthat/test-metab_bayes.R +++ b/tests/testthat/test-metab_bayes.R @@ -6,6 +6,82 @@ # skip_on_appveyor() # skip_if_not_installed('deSolve') + +# prepdata_bayes() --------------------------------------------------------- +# +# Direct, fast (no Stan) unit tests for prepdata_bayes()'s matrix-pivot and +# contiguous-sort-check behavior. Written to pin current behavior ahead of +# extracting shared logic with prepdata_bayes_2s() (see mm_time_by_date_matrix.R). + +# Build a minimal, traceable two-day data.frame: 4 rows/day, DO.obs == +# original row index so the pivoted matrix's values can be checked exactly. +make_bayes_prepdata <- function(n_per_day=4, n_days=2, shuffle=FALSE) { + n_total <- n_per_day * n_days + solar.time <- as.POSIXct("2050-06-01 00:00:00", tz="UTC") + + as.difftime( + rep((seq_len(n_per_day) - 1) * (24 / n_per_day), n_days) + + rep((seq_len(n_days) - 1) * 24, each=n_per_day), + units="hours") + dat <- data.frame( + solar.time = solar.time, + date = as.Date(solar.time), + DO.obs = seq_len(n_total), + DO.sat = rep(10, n_total), + depth = rep(1, n_total), + temp.water = rep(20, n_total), + light = rep(100, n_total) + ) + if(shuffle) { + # interleave day1/day2 rows to break contiguity while keeping an equal + # row count per date (so the earlier date_table-based check still passes + # and only the contiguous-sort check is exercised) + dat <- dat[order(rep(seq_len(n_per_day), n_days)), ] + } + dat +} + +test_that("prepdata_bayes() pivots data into the expected n x d matrix (dims and values)", { + sp <- specs(mm_name('bayes')) + dat <- make_bayes_prepdata(n_per_day=4, n_days=2) + + out <- prepdata_bayes(data=dat, data_daily=NULL, ply_date=NA, specs=sp) + + expect_equal(out$d, 2) + expect_equal(out$n, 4) + expect_equal(dim(out$DO_obs), c(4, 2)) + expect_equal(dim(out$DO_sat), c(4, 2)) + expect_equal(dim(out$depth), c(4, 2)) + # values: matrix(vec, nrow=4, ncol=2, byrow=FALSE) fills column-wise, so + # day 1 = original rows 1:4, day 2 = original rows 5:8 + expect_equal(out$DO_obs[,1], as.numeric(1:4)) + expect_equal(out$DO_obs[,2], as.numeric(5:8)) +}) + +test_that("prepdata_bayes()'s contiguous-sort check passes for properly sorted data", { + sp <- specs(mm_name('bayes')) + dat <- make_bayes_prepdata(n_per_day=4, n_days=2) + + expect_silent(prepdata_bayes(data=dat, data_daily=NULL, ply_date=NA, specs=sp)) +}) + +test_that("prepdata_bayes() errors for data that isn't sorted by date", { + # HISTORY: prior to the mm_time_by_date_matrix()/mm_check_dates_contiguous() + # extraction, this check was `if(!all.equal(unique_dates, names(date_table))) + # stop("couldn't fit given dates into matrix")`, unguarded by isTRUE(). + # Confirmed empirically pre-refactor: whenever the dates were truly + # non-contiguous, all.equal() returned a character vector (not TRUE), and + # `!` on a character vector always errors with "invalid argument type" in + # base R -- so the intended "couldn't fit given dates into matrix" message + # was actually unreachable. The shared mm_check_dates_contiguous() helper + # fixes this as a side effect of the extraction, so this test now + # asserts the originally-intended message. + sp <- specs(mm_name('bayes')) + dat <- make_bayes_prepdata(n_per_day=4, n_days=2, shuffle=TRUE) + + expect_error(prepdata_bayes(data=dat, data_daily=NULL, ply_date=NA, specs=sp), "couldn't fit given dates into matrix") +}) + + manual_test4 <- function() { library(streamMetabolizer) diff --git a/tests/testthat/test-metab_bayes_2s.R b/tests/testthat/test-metab_bayes_2s.R new file mode 100644 index 00000000..34abe0f1 --- /dev/null +++ b/tests/testthat/test-metab_bayes_2s.R @@ -0,0 +1,275 @@ + +# Build a minimal, valid two-station data.frame. Defaults give a 5-minute +# timestep (0.0034722 days) and a 0.01-day travel time, so +# max_lag = round(0.01 / 0.0034722) = 3 timesteps of required upstream lead-in. +make_2station_data <- function(n=10, timestep_min=5, travel_time=0.01) { + data.frame( + solar.time = as.POSIXct("2050-06-01 00:00:00", tz="UTC") + + as.difftime((seq_len(n) - 1) * timestep_min, units="mins"), + DO.obs.up = rep(9, n), + DO.sat.up = rep(10, n), + DO.obs.down = rep(8.8, n), + DO.sat.down = rep(9.9, n), + light = rep(300, n), + depth = rep(0.5, n), + temp.water = rep(20, n), + travel.time = rep(travel_time, n) + ) +} + +test_that("mm_validate_data catches missing required columns", { + dat <- dplyr::select(make_2station_data(), -DO.obs.up) + expect_error(metab_bayes_2s(data=dat), "missing these columns") +}) + +test_that("travel.time <= 0 triggers an error", { + dat <- make_2station_data(travel_time=0) + expect_error(metab_bayes_2s(data=dat), "travel.time must be > 0") + + dat <- make_2station_data(travel_time=-0.01) + expect_error(metab_bayes_2s(data=dat), "travel.time must be > 0") +}) + +test_that("travel.time > 8/24 days (8 hours) triggers an error with a units/limit hint", { + dat <- make_2station_data(travel_time=1) + expect_error(metab_bayes_2s(data=dat), "travel.time must be <= 8/24 days.*incorrect units") + + dat <- make_2station_data(travel_time=1.5) + expect_error(metab_bayes_2s(data=dat), "travel.time must be <= 8/24 days.*incorrect units") +}) + +test_that("insufficient lead-in data triggers an error", { + # 2 rows but max_lag=3 timesteps of upstream lead-in are needed + dat <- make_2station_data(n=2) + expect_error(metab_bayes_2s(data=dat), "insufficient lead-in data") +}) + + +# prepdata_bayes_2s() ----------------------------------------------------- + +# Build a two-day, unit-labeled data.frame with a known, traceable +# DO.obs.up/DO.sat.up series (sequential integers) so the shift can be +# checked by exact value, plus a leading lead-in block. 5-minute timestep +# (0.0034722 days) and 0.01-day travel.time give +# max_lag = round(0.01 / 0.0034722) = 3 lead-in timesteps. +# Day 1 = 10 rows (3 lead-in + 7 modeled), Day 2 = 7 rows (all modeled), so +# both modeled days end up with n_obs = 7 rows. +make_ts_data <- function(n_leadin=3, n_day1=10, n_day2=7, travel_time=0.01, unitted=FALSE) { + n_total <- n_day1 + n_day2 + solar.time <- c( + as.POSIXct("2050-06-01 00:00:00", tz="UTC") + as.difftime((seq_len(n_day1) - 1) * 5, units="mins"), + as.POSIXct("2050-06-02 00:00:00", tz="UTC") + as.difftime((seq_len(n_day2) - 1) * 5, units="mins")) + dat <- data.frame( + solar.time = solar.time, + DO.obs.up = seq_len(n_total), # traceable: value == original row index + DO.sat.up = seq_len(n_total) + 100, # traceable, offset so it's distinguishable from DO.obs.up + DO.obs.down = seq_len(n_total) + 1000, # traceable, offset so it's distinguishable from up/sat values + DO.sat.down = rep(9.9, n_total), + light = rep(300, n_total), + depth = rep(0.5, n_total), + temp.water = rep(20, n_total), + travel.time = rep(travel_time, n_total) + ) + if(unitted) { + units_template <- get_units(mm_data( + solar.time, DO.obs.up, DO.sat.up, DO.obs.down, DO.sat.down, light, depth, temp.water, travel.time)) + dat <- u(dat, unname(units_template[names(dat)])) + } + dat +} + +test_that("upstream DO is shifted by the correct lag", { + dat <- make_ts_data() + out <- prepdata_bayes_2s(dat) + + # max_lag=3, so modeled row i (original index i) uses upstream data from + # original row (i - 3). day 1's 7 modeled rows are original rows 4:10, so + # they pick up DO.obs.up from original rows 1:7; day 2's 7 modeled rows are + # original rows 11:17, picking up DO.obs.up from original rows 8:14. + expect_equal(out$DO_obs_up[,1], as.numeric(1:7)) + expect_equal(out$DO_obs_up[,2], as.numeric(8:14)) + expect_equal(out$DO_sat_up[,1], as.numeric(1:7) + 100) + expect_equal(out$DO_sat_up[,2], as.numeric(8:14) + 100) +}) + +test_that("lead-in rows are excluded from the output matrices", { + dat <- make_ts_data(n_leadin=3, n_day1=10, n_day2=7) + out <- prepdata_bayes_2s(dat) + + # 17 total rows in, 3 are lead-in-only, so 14 modeled rows should remain + expect_equal(out$n_obs * out$n_days, nrow(dat) - 3) + # none of the lead-in DO.obs.up values (1, 2, 3) should appear as a + # DOWNSTREAM-paired value, i.e., the first modeled column should start at + # the shifted value 1, not 1:3 appearing as downstream/lead-in rows + expect_false(any(dat$DO.obs.down[1:3] %in% unlist(out$DO_obs_down))) +}) + +test_that("output matrices have n_obs x n_days dimensions", { + dat <- make_ts_data() + out <- prepdata_bayes_2s(dat) + + expect_equal(out$n_obs, 7) + expect_equal(out$n_days, 2) + for(varname in c('DO_obs_up','DO_sat_up','DO_obs_down','DO_sat_down','light','depth','temp_water','travel_time')) { + expect_equal(dim(out[[varname]]), c(7, 2), info=varname) + } +}) + +test_that("all required Stan data block variables are present", { + dat <- make_ts_data() + # K600_lnorm_meanlog/sdlog are owned by specs() (see PR D-6/I1); pass + # distinctive marker values here to confirm prepdata_bayes_2s() just reads + # them through from specs rather than computing its own defaults + out <- prepdata_bayes_2s(dat, specs=list(K600_lnorm_meanlog=1.23, K600_lnorm_sdlog=4.56)) + + expected_names <- c( + 'n_obs','n_days','DO_obs_up','DO_sat_up','DO_obs_down','DO_sat_down', + 'light','depth','temp_water','travel_time','K600_lnorm_meanlog','K600_lnorm_sdlog') + expect_true(all(expected_names %in% names(out))) + + expect_equal(out$K600_lnorm_meanlog, 1.23) + expect_equal(out$K600_lnorm_sdlog, 4.56) +}) + +test_that("units are stripped from all numeric outputs", { + dat <- make_ts_data(unitted=TRUE) + expect_true(is.unitted(dat)) + + out <- prepdata_bayes_2s(dat, specs=list(K600_lnorm_meanlog=2.484907, K600_lnorm_sdlog=1.0)) + for(varname in c('DO_obs_up','DO_sat_up','DO_obs_down','DO_sat_down','light','depth','temp_water','travel_time')) { + expect_false(is.unitted(out[[varname]]), info=varname) + } + expect_false(is.unitted(out$K600_lnorm_meanlog)) + expect_false(is.unitted(out$K600_lnorm_sdlog)) + expect_false(is.unitted(out$n_obs)) + expect_false(is.unitted(out$n_days)) +}) + + +# mm_parse_name() for two-station models --------------------------------- + +test_that("mm_parse_name recognizes the b2_ prefix for two-station models", { + parsed <- mm_parse_name('b2_np_oi_tr_plrckm.stan') + + expect_equal(parsed$type, 'bayes_2s') + # the rest of the name is shared syntax with one-station bayes models and + # should parse the same way regardless of the b vs. b2 prefix + expect_equal(parsed$pool_K600, 'none') + expect_true(parsed$err_obs_iid) + expect_false(parsed$err_proc_acor) + expect_false(parsed$err_proc_iid) + expect_false(parsed$err_proc_GPP) + expect_equal(parsed$ode_method, 'trapezoid') + expect_equal(parsed$GPP_fun, 'linlight') + expect_equal(parsed$ER_fun, 'constant') + expect_equal(parsed$deficit_src, 'DO_mod') + expect_equal(parsed$engine, 'stan') + + # a one-station name with the same suffix should still parse as plain 'bayes' + expect_equal(mm_parse_name('b_np_oi_tr_plrckm.stan')$type, 'bayes') +}) + + +# mm_name() / mm_valid_names() / specs() for bayes_2s --------------- + +test_that("mm_name(type='bayes_2s') returns the single two-station model name", { + expect_equal(mm_name(type='bayes_2s'), 'b2_np_oi_tr_plrckm.stan') +}) + +test_that("mm_valid_names('bayes_2s') returns the single two-station model name", { + expect_equal(mm_valid_names('bayes_2s'), 'b2_np_oi_tr_plrckm.stan') +}) + +test_that("specs(mm_name('bayes_2s')) has the expected params_in/params_out/split_dates", { + sp <- specs(mm_name('bayes_2s')) + + expect_equal( + sp$params_in, + c('GPP_daily_mu', 'GPP_daily_sigma', 'ER_daily_mu', 'ER_daily_sigma', + 'K600_lnorm_meanlog', 'K600_lnorm_sdlog')) + expect_equal(sp$params_out, c('GPP_daily', 'ER_daily', 'K600_daily', 'sigma', 'metab')) + expect_false(sp$split_dates) + expect_equal(sp$engine, 'stan') +}) + + +# metab_bayes_2s() fitting, predict_metab(), predict_DO() ------------------ + +# Subset two_station_example to just a few modeled days for a faster test +# fit. Naively slicing rows doesn't work: max_lag (the number of upstream +# lead-in rows required) is recomputed from whatever travel.time values are +# present in the slice, so an arbitrary row range can leave a partial first +# date once prepdata_bayes_2s() trims max_lag rows off the front -- the same +# lead-in-sizing logic used in data-raw/two_station_example.R is needed here +# too. +subset_2station_data <- function(full_data, n_modeled_days) { + solar_time <- v(full_data$solar.time) + timestep_days <- stats::median(as.numeric(diff(solar_time), units='days')) + + all_dates <- unique(as.Date(solar_time)) + modeled_dates <- all_dates[2:(1 + n_modeled_days)] + modeled_start <- as.POSIXct(paste0(modeled_dates[1], ' 00:00:00'), tz='UTC') + modeled_end <- as.POSIXct(paste0(modeled_dates[length(modeled_dates)], ' 23:45:00'), tz='UTC') + + candidate_start <- modeled_start - as.difftime(1, units='days') + candidate <- full_data[solar_time >= candidate_start & solar_time <= modeled_end, ] + max_lag <- max(round(v(candidate$travel.time) / timestep_days)) + lead_in_start <- modeled_start - as.difftime(max_lag * timestep_days, units='days') + + full_data[solar_time >= lead_in_start & solar_time <= modeled_end, ] +} + +test_that("metab() fits a two-station model and predict_metab()/predict_DO() work", { + skip_on_cran() + skip_if_not_installed('rstan') + + small_dat <- subset_2station_data(two_station_example, n_modeled_days=3) + + sp <- specs( + mm_name('bayes_2s'), + n_chains=1, n_cores=1, burnin_steps=100, saved_steps=100, verbose=FALSE) + + mm <- metab(specs=sp, data=small_dat) + expect_s4_class(mm, 'metab_bayes_2s') + + pm <- predict_metab(mm) + expect_s3_class(pm, 'data.frame') + expect_true(all(c('GPP','ER','K600') %in% names(pm))) + expect_equal(nrow(pm), 3) + + pdo <- predict_DO(mm) + expect_s3_class(pdo, 'data.frame') + expect_true(all(c('DO.obs.down','DO.mod.down') %in% names(pdo))) +}) + +test_that("a failed Stan run (mode==2L) warns and continues, matching runstan_bayes()'s pattern, rather than erroring out", { + skip_if_not_installed('rstan') + + # stand in for rstan::stan()'s return value on a failed run: only the + # 'mode' slot is inspected by metab_bayes_2s() before deciding to skip + # post-processing, so a minimal S4 object with that slot is sufficient + setClass('fake_failed_stanfit', representation(mode='integer')) + fake_stanfit <- methods::new('fake_failed_stanfit', mode=2L) + testthat::local_mocked_bindings(stan=function(...) fake_stanfit, .package='rstan') + + dat <- make_ts_data() + sp <- specs( + mm_name('bayes_2s'), + n_chains=1, n_cores=1, burnin_steps=10, saved_steps=10, verbose=FALSE) + + expect_warning( + mm <- metab_bayes_2s(specs=sp, data=dat), + 'Modeling failed') + + expect_s4_class(mm, 'metab_bayes_2s') + fit <- mm@fit + expect_true(nrow(fit$daily) > 0) + expect_true(all(is.na(fit$daily$GPP_daily_50pct))) + expect_true(all(is.na(fit$daily$ER_daily_50pct))) + expect_true(all(is.na(fit$daily$K600_daily_50pct))) + expect_true(all(fit$daily$valid_day)) + expect_null(fit$inst) + expect_equal(length(fit$errors), 0) + expect_true(length(fit$warnings) > 0) + expect_true(any(grepl('fake_failed_stanfit', fit$warnings))) +}) diff --git a/tests/testthat/test-mm_determine_cores.R b/tests/testthat/test-mm_determine_cores.R new file mode 100644 index 00000000..4c8c615d --- /dev/null +++ b/tests/testthat/test-mm_determine_cores.R @@ -0,0 +1,46 @@ +context("mm_determine_cores") + +test_that("falls back to 1 core when detectCores() is non-finite", { + local_mocked_bindings(detectCores = function() NA_real_, .package = "parallel") + expect_equal(mm_determine_cores(n_cores=5), 1) + + local_mocked_bindings(detectCores = function() Inf, .package = "parallel") + expect_equal(mm_determine_cores(n_cores=5), 1) +}) + +test_that("caps the requested core count at the number detected", { + local_mocked_bindings(detectCores = function() 4, .package = "parallel") + expect_equal(mm_determine_cores(n_cores=10), 4) + expect_equal(mm_determine_cores(n_cores=2), 2) + expect_equal(mm_determine_cores(n_cores=4), 4) +}) + +test_that("emits a status message only when verbose=TRUE", { + local_mocked_bindings(detectCores = function() 4, .package = "parallel") + + expect_message( + result <- mm_determine_cores(n_cores=10, n_chains=3, verbose=TRUE), + "requesting 3 chains on 4 of 4 available cores", + fixed=TRUE) + expect_equal(result, 4) + + expect_no_message(mm_determine_cores(n_cores=10, n_chains=3, verbose=FALSE)) +}) + +test_that("verbose message matches runstan_bayes()'s original wording exactly", { + local_mocked_bindings(detectCores = function() 8, .package = "parallel") + + expect_message( + mm_determine_cores(n_cores=4, n_chains=4, verbose=TRUE), + "MCMC (Stan): requesting 4 chains on 4 of 8 available cores", + fixed=TRUE) +}) + +test_that("omits the chains clause when n_chains is not supplied", { + local_mocked_bindings(detectCores = function() 4, .package = "parallel") + + expect_message( + mm_determine_cores(n_cores=10, verbose=TRUE), + "MCMC (Stan): requesting 4 of 4 available cores", + fixed=TRUE) +})