diff --git a/R/fme.R b/R/fme.R index a5f2abb..e016f95 100644 --- a/R/fme.R +++ b/R/fme.R @@ -11,21 +11,15 @@ make_memc_objective <- function(comp_data, x, config) { assert_that("time" %in% names(comp_data), msg = "comp_data must contain time column") comp_data_vars <- names(comp_data)[names(comp_data) != "time"] - assert_that(comp_data_vars %in% names(MEMC::default_initial), msg = "comp_data must contain a MEMC variable") + assert_that(all(comp_data_vars %in% names(MEMC::default_initial)), msg = "comp_data must contain a MEMC variable") fxn <- function(x) { - # split up up the input into model parameters and state value - xx <- split_param_state(x) - new_p <- xx$params - new_s <- xx$state # the time steps to evaluate the model at t <- seq(0, max(comp_data$time)) # solve the model - new_config <- update_config(mod = config, - params = new_p, - state = new_s) + new_config <- update_config(mod = config, new = x) out <- sm_internal(new_config, t) # make sure that the model solved for all time steps @@ -49,9 +43,10 @@ make_memc_objective <- function(comp_data, x, config) { #' Fit a MEMC model to a comparison data #' -#' @param x memc model parameters or initial conditions that will be fit to the data +#' @param x memc model parameters or initial conditions that will be fit to the data, users will need to provide an initial guess for these values. #' @param config memc model configuration object, either one of the pre-built configurations listed in \code{model_configs} or created using \code{configure_model} -#' @param comp_data data frame containing the comparison data that the model will be fit to +#' @param comp_data data frame containing the comparison data that the model will +#' be fit this data frame must contain a column for time, the other columns must be named for the MEMC model variables. #' @param lower lower bounds on the parameters; if unbounded set equal to -Inf #' @param upper bounds on the parameters; if unbounded set equal to Inf #' @param ... addition arguments that may be passed to FME::modFit @@ -87,7 +82,7 @@ memc_modfit <- #' #' @param config a memc model configuration object, either one of the pre-built configurations listed in \code{model_configs} or created using \code{configure_model} #' @param t vector of the time steps to run the model at -#' @param pars vector of the parameters that will be varied +#' @param x vector of the parameters or initial model pool sizes that will be varied #' @param parRange data frame of the min/max parameter values #' @param dist str for the distribution according to which the parameters will be sampled from, options" "unif" (uniformly random samples), "norm", (normally distributed random samples), "latin" (latin hypercube distribution), and "grid" (parameters arranged on a grid). #' @param ... additional arguments passed to FME::sensRange @@ -101,7 +96,7 @@ memc_modfit <- #' prange <- data.frame(min = pars - pars * 0.75, #' max = pars + pars * 0.75) #' t <- floor(seq(0, 365, length.out = 10)) -#' out <- memc_sensrange(config = MEND_model, t = t, pars = pars, +#' out <- memc_sensrange(config = MEND_model, t = t, x = pars, #' parRange = prange, dist = "latin", num = 10) #' plot(summary(out)) #' # Using the helper functions. @@ -111,14 +106,14 @@ memc_modfit <- #' geom_ribbon(aes(time, ymin = Min, ymax = Max), alpha = 0.5) + #' facet_wrap("variable", scales = "free") #'} -memc_sensrange <- function(config, t, pars, parRange, dist, ...){ +memc_sensrange <- function(config, t, x, parRange, dist, ...){ func <- function(pars){ - new_mod <- update_config(mod = config, params = pars) + new_mod <- update_config(mod = config, new = pars) sm_internal(mod = new_mod, time = t) } - out <- FME::sensRange(func, parms = pars, dist = dist, parRange = parRange, ...) + out <- FME::sensRange(func, parms = x, dist = dist, parRange = parRange, ...) return(out) } @@ -130,7 +125,7 @@ memc_sensrange <- function(config, t, pars, parRange, dist, ...){ #' #' @param config a memc model configuration object, either one of the pre-built configurations listed in \code{model_configs} or created using \code{configure_model} #' @param t vector of the time steps to run the model at -#' @param pars vector of the parameters to test +#' @param x vector of the parameters or initial state values to test #' @param ... additional arguments passed to FME::sensFun #' @return the results of the FME::sensFun #' @export @@ -139,7 +134,7 @@ memc_sensrange <- function(config, t, pars, parRange, dist, ...){ #'\dontrun{ #' # Test the sensitivity of the MEND output for V.p, K.p, V.m #' pars <- c("V.d" = 3.0e+00,"V.p" = 1.4e+01,"V.m" = 2.5e-01) -#' out <- memc_sensfunc(config = MEND_model, t = t, pars = pars) +#' out <- memc_sensfunc(config = MEND_model, t = t, x = pars) #' pairs(out) #' plot(out) #' # Using the helper functions to make nice ggplots @@ -148,14 +143,14 @@ memc_sensrange <- function(config, t, pars, parRange, dist, ...){ #' geom_line(aes(time, value, color = parameter)) + #' facet_wrap("variable", scales = "free") #'} -memc_sensfunc <- function(config, t, pars, ...){ +memc_sensfunc <- function(config, t, x, ...){ - func <- function(pars){ - new_mod <- update_config(mod = config, params = pars) + func <- function(x){ + new_mod <- update_config(mod = config, new = x) sm_internal(mod = new_mod, time = t) } - out <- FME::sensFun(func, pars, ...) + out <- FME::sensFun(func, x, ...) } diff --git a/R/helper.R b/R/helper.R index d710ab3..584534b 100644 --- a/R/helper.R +++ b/R/helper.R @@ -62,23 +62,29 @@ update_state <- function(new_vals, state) { #' Update a model configuration this is for internal function use #' -#' @param mod TODO -#' @param params TODO -#' @param state TODO +#' @param mod a MEMC model configuration object created by \code{configure_model} +#' @param new vector containing the parameters and or initial pool values #' @import assertthat #' @family helper functions #' @noRd -update_config <- function(mod, params = NULL, state = NULL) { +update_config <- function(mod, new = NULL) { + assert_that(is_memc_config(mod)) - - if (!is.null(params)) { + + if(is.null(new)){ + return(mod) + } + + x <- split_param_state(new) + + if (length(x$params) >= 1) { mod[["params"]] <- - update_params(new_params = params, param_table = mod[["params"]]) + update_params(new_params = x$params, param_table = mod[["params"]]) } - if (!is.null(state)) { + if (length(x$state) >= 1) { mod[["state"]] <- - update_state(new_vals = state, state = mod[["state"]]) + update_state(new_vals = x$state, state = mod[["state"]]) } diff --git a/R/model_componet.R b/R/model_componet.R index 7f47e52..03fd592 100644 --- a/R/model_componet.R +++ b/R/model_componet.R @@ -130,7 +130,7 @@ carbon_pool_derivs <- # POM = particulate organic carbon dPOM <- (1 - p[["g_d"]]) * F6 - F2 + p[["Input_POM"]] # MOM = mineral-associated organic carbon (MOC) - dMOM <- (1 - p[["f_d"]]) * F2 - F3 + p[["Input_MOM"]] + dMOM <- (1 - p[["f_d"]]) * F2 - F3 # QOMO = active layer of MOC dQOM <- F4 - F5 # MB = microbial biomass carbon @@ -145,7 +145,7 @@ carbon_pool_derivs <- # IC = inorganic carbon (CO2) dIC <- F1 * (1 - p[["CUE"]]) # Tot = the total carbon pool - dTot <- -F1 * (1 - p[["CUE"]]) + (p[["Input_POM"]] + p[["Input_DOM"]] + p[["Input_MOM"]]) + dTot <- -F1 * (1 - p[["CUE"]]) + (p[["Input_POM"]] + p[["Input_DOM"]]) # Return outputs return(list(c(dPOM, dMOM, dQOM, dMB, dDOM, dEP, dEM, dIC, dTot))) @@ -166,6 +166,12 @@ sm_internal <- function(mod, time, ...) { p <- mod[["params"]][["value"]] names(p) <- mod[["params"]][["parameter"]] + + # Check that all the parameters that are fractions are less than 1 + frac_params <- c("f_d", "g_d", "p_ep", "p_em") + frac_params_vals <- p[names(p) %in% frac_params] + assert_that(all(0 < frac_params_vals & frac_params_vals < 1), + msg = "parameters f_d, g_d, p_ep, and p_em must be between 0 and 1") rslt <- deSolve::ode( y = mod[["state"]], @@ -229,14 +235,13 @@ solve_model <- # Update the model configuration with new parameter and initial state values mod <- update_config(mod = mod, - params = params, - state = state) + new = c(params, state)) # Check the arguments assert_that(is_memc_config(obj = mod)) assert_that(is_param_table(table = mod$params)) assert_that(is_state_vector(state = mod$state)) - + results <- sm_internal(mod = mod, time = time, ...) out <- sm_format_out(rslt = results, mod = mod) diff --git a/README.md b/README.md index e8fdb5c..c4aced3 100644 --- a/README.md +++ b/README.md @@ -94,8 +94,7 @@ print(MEND_model) #> 16 dd_beta strength of density dependent microbial decay #> 17 Input_POM POM input #> 18 Input_DOM DOM input -#> 19 Input_MOM MOM input -#> 20 CUE carbon use efficiency +#> 19 CUE carbon use efficiency #> units value #> 1 mgC mgC^-1 h^-1 14.000 #> 2 mgC / g soil 50.000 @@ -115,8 +114,7 @@ print(MEND_model) #> 16 1.000 #> 17 mg C 0.000 #> 18 mg C 0.000 -#> 19 mg C 0.000 -#> 20 0.400 +#> 19 0.400 #> #> $state #> POM MOM QOM MB DOM EP EM IC diff --git a/_pkgdown.yml b/_pkgdown.yml index ddeff9a..fe5cdeb 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -1,9 +1,16 @@ url: https://microbial-explicit-model.github.io/MEMC/ template: bootstrap: 5 +development: + mode: auto navbar: - title: "MEMC" + structure: + right: [search, github] + right: + - icon: fa-github fa-lg + href: "https://github.com/JGCRI/hector" + left: - icon: fa-home fa-lg href: index.html @@ -18,4 +25,4 @@ navbar: - text: "Functions" href: reference/index.html - text: "News" - href: news/index.html + href: news/index.html \ No newline at end of file diff --git a/data-raw/default_params.csv b/data-raw/default_params.csv index d5ecace..073625b 100644 --- a/data-raw/default_params.csv +++ b/data-raw/default_params.csv @@ -1,21 +1,20 @@ -parameter,description,units,value -V_p,maximum specific decomposition rate for POM by EP,mgC mgC^-1 h^-1,14 -K_p,half-saturation constant for decomposition of POM,mgC / g soil,50 -V_m,maximum specific decomposition rate for MOM by EM,mgC mgC^-1 h^-1,0.25 -K_m,half-saturation constant for decomposition of MOM,mg C/g soil,250 -V_d,maximum specific uptake rate of D for growth of MB,mgC mgC^-1 h^-1,3.00E+00 -K_d,half-saturation constant of uptake of D for growth of MB,mg C/g soil,0.25 -f_d,fraction of decomposed P allocated to DOM,NA,0.5 -g_d,fraction of dead B allocated to DOM,NA,0.5 -p_ep,fraction of mR for production of EP,NA,0.01 -p_em,fraction of mR for production of EM,NA,0.01 -r_ep,turnover rate of EP,mgC mgC^-1 h^-1,0.001 -r_em,turnover rate of EM,mgC mgC^-1 h^-1,0.001 -Q_max,maximum DOC sorption capacity,mgC / g soil,3.4 -K_ads,specific adsorption rate,mgC mgC^-1 h^-1,0.006 -K_des,desorption rate,mgC mgC^-1 h^-1,0.001 -dd_beta,strength of density dependent microbial decay,NA,1 -Input_POM,POM input ,mg C,0 -Input_DOM,DOM input,mg C,0 -Input_MOM,MOM input ,mg C,0 +parameter,description,units,value +V_p,maximum specific decomposition rate for POM by EP,mgC mgC^-1 h^-1,14 +K_p,half-saturation constant for decomposition of POM,mgC / g soil,50 +V_m,maximum specific decomposition rate for MOM by EM,mgC mgC^-1 h^-1,0.25 +K_m,half-saturation constant for decomposition of MOM,mg C/g soil,250 +V_d,maximum specific uptake rate of D for growth of MB,mgC mgC^-1 h^-1,3.00E+00 +K_d,half-saturation constant of uptake of D for growth of MB,mg C/g soil,0.25 +f_d,fraction of decomposed P allocated to DOM,NA,0.5 +g_d,fraction of dead B allocated to DOM,NA,0.5 +p_ep,fraction of mR for production of EP,NA,0.01 +p_em,fraction of mR for production of EM,NA,0.01 +r_ep,turnover rate of EP,mgC mgC^-1 h^-1,0.001 +r_em,turnover rate of EM,mgC mgC^-1 h^-1,0.001 +Q_max,maximum DOC sorption capacity,mgC / g soil,3.4 +K_ads,specific adsorption rate,mgC mgC^-1 h^-1,0.006 +K_des,desorption rate,mgC mgC^-1 h^-1,0.001 +dd_beta,strength of density dependent microbial decay,NA,1 +Input_POM,POM input ,mg C,0 +Input_DOM,DOM input,mg C,0 CUE,carbon use efficiency,,0.4 \ No newline at end of file diff --git a/data/BAMS_model.rda b/data/BAMS_model.rda index ce884eb..53de116 100644 Binary files a/data/BAMS_model.rda and b/data/BAMS_model.rda differ diff --git a/data/COMISSION_model.rda b/data/COMISSION_model.rda index 6012d83..1d79fbe 100644 Binary files a/data/COMISSION_model.rda and b/data/COMISSION_model.rda differ diff --git a/data/CORPSE_model.rda b/data/CORPSE_model.rda index 5cb5007..5c4eafb 100644 Binary files a/data/CORPSE_model.rda and b/data/CORPSE_model.rda differ diff --git a/data/MEMS_model.rda b/data/MEMS_model.rda index bf94272..ab32e10 100644 Binary files a/data/MEMS_model.rda and b/data/MEMS_model.rda differ diff --git a/data/MEND_model.rda b/data/MEND_model.rda index 5d2db63..1fd44e1 100644 Binary files a/data/MEND_model.rda and b/data/MEND_model.rda differ diff --git a/data/MIMCS_model.rda b/data/MIMCS_model.rda index e9e1857..756c055 100644 Binary files a/data/MIMCS_model.rda and b/data/MIMCS_model.rda differ diff --git a/data/default_initial.rda b/data/default_initial.rda index d573a1a..6ea9078 100644 Binary files a/data/default_initial.rda and b/data/default_initial.rda differ diff --git a/data/default_params.rda b/data/default_params.rda index e0409b7..dd82a71 100644 Binary files a/data/default_params.rda and b/data/default_params.rda differ diff --git a/data/model_configs.rda b/data/model_configs.rda index 514154b..3f1d027 100644 Binary files a/data/model_configs.rda and b/data/model_configs.rda differ diff --git a/legacy/recreating_MEMC_Dec7.R b/legacy/recreating_MEMC_Dec7.R index de2ea53..ff90f43 100644 --- a/legacy/recreating_MEMC_Dec7.R +++ b/legacy/recreating_MEMC_Dec7.R @@ -15,7 +15,7 @@ desc <- c("maximum specific decomposition rate for P by EP", "half-saturation co "turnover rate of EP", "turnover rate of EM", "maximum DOC sorption capacity", "specific adsorption rate", "desorption rate", "strength of density dependent microbial decay", "POM input", "DOM input", "M input", "carbon use efficiency") value <- c("V_p"=14, "K_p"=50, "V_m"=0.25, "K_m"=250, "V_d"=3, "K_d"=0.25, "f_d"=0.5, "g_d"=0.5, "p_ep"=0.01, "p_em"=0.01, "r_ep"=1e-3, "r_em"=1e-3 , "Q_max"=3.4, "K_ads"=0.006, - "K_des"= 0.001, "dd_beta" = 2, "Input_POM" = 0, "Input_DOM"=0, "Input_MOM"=0, "CUE"=0.4) + "K_des"= 0.001, "dd_beta" = 2, "Input_POM" = 0, "Input_DOM"=0, "CUE"=0.4) units <- c("mgC mgC^-1 h^-1", "mgC / g soil", "mgC mgC^-1 h^-1", "mg C/g soil", "mgC mgC^-1 h^-1", "mg C/g soil", NA, NA, NA, NA, "mgC mgC^-1 h^-1", "mgC mgC^-1 h^-1", "mgC / g soil", "mgC mgC^-1 h^-1", "mgC mgC^-1 h^-1", NA, "mg C", "mg C", "mg C", "") diff --git a/man/figures/README-unnamed-chunk-6-1.png b/man/figures/README-unnamed-chunk-6-1.png index de9a4ec..15125b8 100644 Binary files a/man/figures/README-unnamed-chunk-6-1.png and b/man/figures/README-unnamed-chunk-6-1.png differ diff --git a/man/figures/README-unnamed-chunk-9-1.png b/man/figures/README-unnamed-chunk-9-1.png index d51dcef..9bc6731 100644 Binary files a/man/figures/README-unnamed-chunk-9-1.png and b/man/figures/README-unnamed-chunk-9-1.png differ diff --git a/man/memc_modfit.Rd b/man/memc_modfit.Rd index cdad38e..8b839b8 100644 --- a/man/memc_modfit.Rd +++ b/man/memc_modfit.Rd @@ -9,7 +9,10 @@ memc_modfit(config, x, comp_data, lower = -Inf, upper = Inf, ...) \arguments{ \item{config}{memc model configuration object, either one of the pre-built configurations listed in \code{model_configs} or created using \code{configure_model}} -\item{comp_data}{data frame containing the comparison data thata the model will be fit to} +\item{x}{memc model parameters or initial conditions that will be fit to the data, users will need to provide an initial guess for these values.} + +\item{comp_data}{data frame containing the comparison data that the model will +be fit this data frame must contain a column for time, the other columns must be named for the MEMC model variables.} \item{lower}{lower bounds on the parameters; if unbounded set equal to -Inf} diff --git a/man/memc_sensfunc.Rd b/man/memc_sensfunc.Rd index e0c9b34..f50b3a0 100644 --- a/man/memc_sensfunc.Rd +++ b/man/memc_sensfunc.Rd @@ -4,14 +4,14 @@ \alias{memc_sensfunc} \title{Local sensitivity for a MEMC model} \usage{ -memc_sensfunc(config, t, pars, ...) +memc_sensfunc(config, t, x, ...) } \arguments{ \item{config}{a memc model configuration object, either one of the pre-built configurations listed in \code{model_configs} or created using \code{configure_model}} \item{t}{vector of the time steps to run the model at} -\item{pars}{vector of the parameters to test} +\item{x}{vector of the parameters or initial state values to test} \item{...}{additional arguments passed to FME::sensFun} } @@ -25,7 +25,7 @@ Estimate the local effect of a parameter on a memc model output \dontrun{ # Test the sensitivity of the MEND output for V.p, K.p, V.m pars <- c("V.d" = 3.0e+00,"V.p" = 1.4e+01,"V.m" = 2.5e-01) -out <- memc_sensfunc(config = MEND_model, t = t, pars = pars) +out <- memc_sensfunc(config = MEND_model, t = t, x = pars) pairs(out) plot(out) # Using the helper functions to make nice ggplots diff --git a/man/memc_sensrange.Rd b/man/memc_sensrange.Rd index 514543c..5c35f48 100644 --- a/man/memc_sensrange.Rd +++ b/man/memc_sensrange.Rd @@ -4,14 +4,14 @@ \alias{memc_sensrange} \title{Global sensitivity ranges for a memc model} \usage{ -memc_sensrange(config, t, pars, parRange, dist, ...) +memc_sensrange(config, t, x, parRange, dist, ...) } \arguments{ \item{config}{a memc model configuration object, either one of the pre-built configurations listed in \code{model_configs} or created using \code{configure_model}} \item{t}{vector of the time steps to run the model at} -\item{pars}{vector of the parameters that will be varied} +\item{x}{vector of the parameters or initial model pool sizes that will be varied} \item{parRange}{data frame of the min/max parameter values} @@ -32,7 +32,7 @@ pars <- c("V.d" = 3.0e+00,"V.p" = 1.4e+01,"V.m" = 2.5e-01) prange <- data.frame(min = pars - pars * 0.75, max = pars + pars * 0.75) t <- floor(seq(0, 365, length.out = 10)) -out <- memc_sensrange(config = MEND_model, t = t, pars = pars, +out <- memc_sensrange(config = MEND_model, t = t, x = pars, parRange = prange, dist = "latin", num = 10) plot(summary(out)) # Using the helper functions. diff --git a/tests/testthat/test-argument.R b/tests/testthat/test-argument.R index 243ec6b..a5b5279 100644 --- a/tests/testthat/test-argument.R +++ b/tests/testthat/test-argument.R @@ -10,7 +10,7 @@ test_that("is_param_table works", { # A number of errors should be thrown if the the parameter table does not meet the correct conditions. expect_error( is_param_table(params[1:10, ]), - "param_table is missing a parameter value(s) for: r_ep, r_em, Q_max, K_ads, K_des, dd_beta, Input_POM, Input_DOM, Input_MOM, CUE", + "param_table is missing a parameter value(s) for: r_ep, r_em, Q_max, K_ads, K_des, dd_beta, Input_POM, Input_DOM, CUE", fixed = TRUE ) diff --git a/tests/testthat/test-fme.R b/tests/testthat/test-fme.R index 4f9ce3b..63bfe8a 100644 --- a/tests/testthat/test-fme.R +++ b/tests/testthat/test-fme.R @@ -65,9 +65,9 @@ test_that("memc_sensrange", { prange <- data.frame(min = p - p * frac, max = p + p * frac) # Run the sense range and check the output. - out <- memc_sensrange(MEND_model, t = time, pars = p, n = n, dist = "latin", parRange = prange) + out <- memc_sensrange(MEND_model, t = time, x = p, n = n, dist = "latin", parRange = prange) expect_true(all(dim(out) == c(n, (length(time) * 9) + 1))) - out2 <- memc_sensrange(MEND_model, t = time, pars = p, num = 6, dist = "latin", parRange = prange) + out2 <- memc_sensrange(MEND_model, t = time, x = p, num = 6, dist = "latin", parRange = prange) expect_true(nrow(out2) > nrow(out)) }) diff --git a/tests/testthat/test-helper.R b/tests/testthat/test-helper.R index 04f571c..c227f2d 100644 --- a/tests/testthat/test-helper.R +++ b/tests/testthat/test-helper.R @@ -54,7 +54,7 @@ test_that("bad model configuration will fail", { # Errors should be thrown expect_error( configure_model(ptable[1:9,], state), - "param_table is missing a parameter value(s) for: p_em, r_ep, r_em, Q_max, K_ads, K_des, dd_beta, Input_POM, Input_DOM, Input_MOM, CUE", + "param_table is missing a parameter value(s) for: p_em, r_ep, r_em, Q_max, K_ads, K_des, dd_beta, Input_POM, Input_DOM, CUE", fixed = TRUE ) expect_error( @@ -74,10 +74,7 @@ test_that("bad model configuration will fail", { ) # Only change one parameter value and one state value - new_out <- - update_config(out1, - params = c("V_d" = 50), - state = c("MB" = 50)) + new_out <- update_config(mod = out1, new = c("V_d" = 50, "MB" = 50)) expect_equal(sum(new_out$params$value != out1$params$value), 1) expect_equal(sum(new_out$state != out1$state), 1) diff --git a/tests/testthat/test-jz.R b/tests/testthat/test-jz.R index f536aad..d82be0d 100644 --- a/tests/testthat/test-jz.R +++ b/tests/testthat/test-jz.R @@ -12,7 +12,7 @@ # # params <- c("V_p", "K_p", "V_m", "K_m", "V_d", "K_d", "f_d", "g_d", "p_ep", # "p_em", "r_ep", "r_em", "Q_max", "K_ads", -# "K_des", "dd_beta", "Input_POM", "Input_DOM", "Input_MOM", "CUE") +# "K_des", "dd_beta", "Input_POM", "Input_DOM", "CUE") # desc <- c("maximum specific decomposition rate for P by EP", "half-saturation constant for decomposition of P", # "maximum specific decomposition rate for M by EM", # "half-saturation constant for decomposition of M", @@ -26,7 +26,7 @@ # value <- c("V_p" = 14, "K_p" = 50, "V_m" = 0.25, "K_m" = 250, "V_d" = 3, "K_d"= 0.25, # "f_d" = 0.5, "g_d" = 0.5, "p_ep" = 0.01, "p_em" = 0.01, "r_ep" = 1e-3, # "r_em" = 1e-3 , "Q_max" = 3.4, "K_ads" = 0.006, "K_des"= 0.001, "dd_beta" = 2, -# "Input_POM" = 0, "Input_DOM" = 0, "Input_MOM" = 0, "CUE" = 0.4) +# "Input_POM" = 0, "Input_DOM" = 0, "CUE" = 0.4) # units <- c("mgC mgC^-1 h^-1", "mgC / g soil", "mgC mgC^-1 h^-1", "mg C/g soil", # "mgC mgC^-1 h^-1", "mg C/g soil", NA, # NA, NA, NA, "mgC mgC^-1 h^-1", "mgC mgC^-1 h^-1", "mgC / g soil", diff --git a/tests/testthat/test-model_componet.R b/tests/testthat/test-model_componet.R index 260be97..5296fff 100644 --- a/tests/testthat/test-model_componet.R +++ b/tests/testthat/test-model_componet.R @@ -56,7 +56,7 @@ test_that("configure_model", { test_that("change param", { out1 <- solve_model(mod = mod, time) - out2 <- solve_model(mod = mod, time, params = c("V_d" = 10)) + out2 <- solve_model(mod = mod, time, params = c("K_m" = 10)) expect_gte(mean((out1$value - out2$value)^2), zero) }) diff --git a/vignettes/articles/Background.Rmd b/vignettes/articles/Background.Rmd index 9794718..7f51e9a 100644 --- a/vignettes/articles/Background.Rmd +++ b/vignettes/articles/Background.Rmd @@ -14,38 +14,45 @@ library(ggplot2) theme_set(theme_bw(base_size = 14)) ``` +# Motivation -### Model Structure +Microbial explicit carbon models include process-based representations of microbial activity in soil carbon models. However, microbial-explicit models are diverse in terms of pool structure, representations of carbon flux dynamics, and environmental drivers, as illustrated in **Figure 1**. -The MEMC package provides users with the ability to explore alternative configurations of a soil organic matter (SOM) model with alternative structural formulations (i.e. representations of fluxes between carbon pools). MEMC uses the common model structure visualized in Figure 1. In the first version of this package, users have the ability to decide which implementation of DOM uptake (F1), POM decomposition (F2), and microbial biomass decay (F8) should be used with the MEMC pool structure (Figure 2). This flexible framework allows users to isolate the effects of flux dynamics on model behavior. +![**Figure 1:** Conceptrual diagram illustrating the variations in pool strucutre, carbon flux dynamics, and environmental drivers of various popular microbial-explicit carbon models.](som_diversity.png) -![Figure 1: Conceptual diagram of the soil organic matter pool structure used by MEMC. The pools include particulate organic matter (POM), mineral-associated organic carbon matter (MOM), dissolved organic carbon matter (DOM), adsorbed phase of DOM (QOM), microbial biomass (MB), POM-degrading enzymes (EP), and MOM-degrading enzymes (EM). External inputs for POM, DOM, and MOM are directly are fed into POM, DOM, and MOM, respectively. The fluxes between the pools are detonated with solid arrows, they include DOM uptake by microbial biomass (F1), POM decomposition (F2), mineral-associated carbon decomposition (F3), adsorption of DOC (F4), microbial biomass decay (F6), and enzyme production used in POM and DOM decomposition (F7) and enzyme turn over (F8). The dashed arrows indicate pools that affect flux rates.](MEMC_conceptual_diagram.png) +\ +The variations in carbon model structures and dynamics contribute to diverging long-term model projections, making model inter-comparisons difficult. In order to explore, attribute, and better understand the implications of varying the representations of the C dynamics, the Microbial Explicit Model Consortium (MEMC) adopts a uniform pool structure to enable easy model comparison exercises. By implementing MEMC as an R package, we hope to make soil carbon modeling accessible to a wider audience. +\ +# Model Structure -![Figure 2: Illustrates the different combinations of soil pool fluxes dynamics for POM decomposition, DOM uptake, and microbial biomass decay supported by MEMC V1. Where MM denotes Michaelis–Menten and Reverse MM refers to reverse Michaelis–Menten.](MEMC_configs.png) +MEMC uses the common model structure visualized in **Figure 2**. In the first version of this package, users have the ability to decide which implementation POM decomposition (1), DOM uptake (2), and microbial biomass decay (3) as indicated in (**Figure 3**). This flexible framework allows users to isolate the effects of flux dynamics on model behavior. +![**Figure 2:** Conceptual diagram of the soil organic matter pool structure used by MEMC. The pools include particulate organic matter (POM), mineral-associated organic carbon matter (MAOM), dissolved organic carbon matter (DOM), adsorbed phase of DOM (QOM), microbial biomass (MB), POM-degrading enzymes (EP), and MOM-degrading enzymes (EM). External inputs for POM, DOM, and MOM are directly are fed into POM, DOM, and MOM, respectively. The fluxes between the pools are detonated with solid arrows, they include POM decomposition (1), DOM uptake by microbial biomass (2), mineral-associated carbon decomposition, adsorption of DOC, microbial biomass decay (3), and enzyme production used in POM and DOM decomposition and enzyme turn over. The dashed arrows indicate enzymes that affect flux rate.](pool_structure.png) -TODO?? +# Model Configurations +The MEMC R package ships with six model pre-configured models (MICMS, MEND, CORPSE, BAMS, COMISSION, and MEMS). All of these model implementations use the uniform pool structure from **Figure 2** but use a different combination of C flux dynamics as indicated in **Figure 3**. -### MM -### Reverse MM +![**Figure 3:** Illustrates the different combinations of soil pool fluxes dynamics for POM decomposition, DOM uptake, and microbial biomass decay supported by MEMC V1. The numbers correspond to the arrows labeled in **Figure 2**. MM denotes Michaelis–Menten and Reverse MM refers to reverse Michaelis–Menten.](configs.png) -### Linear -### Density Dependent +# Dynamics +* **MM:** Michaelis Menten model, a commonly used equation in biogeochemical modelings which describe the rate of a reaction in response to a change in substrate. +* **Reverse MM:** Reverse Michaelis Menten model (RMM), a modified version of Michaelis Menten model, however now the modeled reaction rate changes varies in response to enzyme availability. +* **Linear:** a simple linear model (LM) is used. +* **Density Dependent:** in the Density Dependent (DD) model the density of the population (the microbial biomass) affects the rate of a given process. +* **ECA:** Equilibrium Chemistry Approximation kinetics an alternative representation of biogeochemical processes that is more appropriate over a wide range of substrate availability. -TODO add descipritons of the MM, RMM, Linear, Den - diff --git a/vignettes/articles/ChangeParams.Rmd b/vignettes/articles/ChangeParams.Rmd new file mode 100644 index 0000000..2013668 --- /dev/null +++ b/vignettes/articles/ChangeParams.Rmd @@ -0,0 +1,103 @@ +--- +title: "Change Model Parameters" +--- + +```{r, include = FALSE, echo = FALSE, warning = FALSE, message = FALSE} +knitr::opts_chunk$set( + collapse = TRUE, + comment = "#>", + fig.width = 8, + fig.height = 5 +) +library(MEMC) +library(ggplot2) +theme_set(theme_bw(base_size = 14)) +``` + +This example will show how to update the model parameter values for any MEMC model configuration. + + +```{r, eval=FALSE} +# Assumes MEMC has already been installed. +library(MEMC) + +# Use to visualize results. +library(ggplot2) +``` + +The `MEMC` package includes default parameter and initial pool values that are included as package data that are based on Wang et al. 2015^[Wang, G., Jagadamma, S., Mayes, M. et al. Microbial dormancy improves development and experimental validation of ecosystem model. ISME J 9, 226–237 (2015). https://doi.org/10.1038/ismej.2014.120]. This example will use these values and show users how to update these values. We will also be running the MEMC MEND model configuration. \ + +# Default Parameters + +Start by taking a look at the default model parameter values. + +```{r} +print(default_initial) +# Alternatively can use the :: to access the MEMC package data. +# print(MEMC::default_initial) +``` + + +```{r} +print(default_params) +# Alternatively can use the :: to access the MEMC package data. +# print(MEMC::default_params) +``` + +# Run deafult MEMC MEND + +Run the pre-configured MEMC MEND model with the default parameter values. + +```{r} +time <- 0:600 +mod1 <- MEND_model +out1 <- solve_model(mod = mod1, time = time) +out1$name <- "Km = 250 (default)" +``` + +# Update parameter value + +Users can update the parameters directly via `solve_model`. Any parameter(s) value(s) fed into the `solve_model` params argument will overwrite the parameter values defined the model configuration. + +```{r} +out2 <- solve_model(mod = MEND_model, time = time, params = c("K_m" = 10)) +out2$name <- "Km = 10 (lower)" +``` + +```{r} +ggplot(data = rbind(out1, out2)) + + geom_line(aes(time, value, color = name, linetype = name), linewidth = 0.75) + + facet_wrap("variable", scales = "free") + + labs(y = "mg C/g soil", + title = "Change Km Values") + + theme(legend.title = element_blank()) +``` + +# Update parameter table using `update_params` + +`update_params` is a helper function that can change the parameter values in the parameter table. You will need to set up a new model configuration with the new parameter table. + +```{r} +new_km <- update_params(new_params = c("K_m" = 100), param_table = default_params) + +# Set up the new model configuration with the +mod_km <- configure_model(params = new_km, + state = default_initial, + name = "Km = 100", + DOMuptake = "MM", + POMdecomp = "MM", + MBdecay = "LM") + +out3 <- solve_model(mod = mod_km, time = time) +``` + + + +```{r} +ggplot(data = rbind(out1, out2, out3)) + + geom_line(aes(time, value, color = name, linetype = name), linewidth = 0.75) + + facet_wrap("variable", scales = "free") + + labs(y = "mg C/g soil", + title = "Change Km Values") + + theme(legend.title = element_blank()) +``` \ No newline at end of file diff --git a/vignettes/articles/Fit-to-Data.Rmd b/vignettes/articles/Fit-to-Data.Rmd index 43dd3aa..f2a24ac 100644 --- a/vignettes/articles/Fit-to-Data.Rmd +++ b/vignettes/articles/Fit-to-Data.Rmd @@ -89,6 +89,10 @@ ggplot(data = out[out$variable == "IC", ]) + geom_point(data = comp_data, aes(time, IC)) + facet_wrap("variable", scales = "free") + labs(x = "Time", y = "mg C/g soil", title = "MIMCS") - ``` + + + + + diff --git a/vignettes/articles/FutureDynamics.png b/vignettes/articles/FutureDynamics.png new file mode 100644 index 0000000..a073fac Binary files /dev/null and b/vignettes/articles/FutureDynamics.png differ diff --git a/vignettes/articles/MEMC_conceptual_diagram.png b/vignettes/articles/MEMC_conceptual_diagram.png deleted file mode 100644 index e6a2f0f..0000000 Binary files a/vignettes/articles/MEMC_conceptual_diagram.png and /dev/null differ diff --git a/vignettes/articles/MEMC_configs.png b/vignettes/articles/MEMC_configs.png deleted file mode 100644 index cd27bf9..0000000 Binary files a/vignettes/articles/MEMC_configs.png and /dev/null differ diff --git a/vignettes/articles/Manual.Rmd b/vignettes/articles/Manual.Rmd index 8839dde..4166989 100644 --- a/vignettes/articles/Manual.Rmd +++ b/vignettes/articles/Manual.Rmd @@ -13,11 +13,15 @@ knitr::opts_chunk$set( # Science * [Background](Background.html) on MEMC model structure and flexible flux dynamics +* Take a look at our tentative [future development plans](RoadMap.html) # Getting Started * [Build & Install MEMC](Install-Build.html) -* +* [How to run all the pre-configured models](ModelComparison.html) +* [How to set up a MEMC model](ToyModel.html) +* [How to change model parameters](ChangeParams.html) +* More examples can be found under the examples tab # Contributing to MEMC @@ -27,4 +31,4 @@ knitr::opts_chunk$set( **** -[Complete index of online documenation](index.html) +[Complete index of online documentation](index.html) diff --git a/vignettes/articles/RoadMap.Rmd b/vignettes/articles/RoadMap.Rmd new file mode 100644 index 0000000..36c81d2 --- /dev/null +++ b/vignettes/articles/RoadMap.Rmd @@ -0,0 +1,9 @@ +--- +title: "Devlopment Plan" +--- + +![](roadmap.png) + +Version 2 of the MEMC Package will support more flexible dynamics, this will affect the model configurations and results. + +![](FutureDynamics.png) diff --git a/vignettes/articles/ToyModel.Rmd b/vignettes/articles/ToyModel.Rmd new file mode 100644 index 0000000..7d906ec --- /dev/null +++ b/vignettes/articles/ToyModel.Rmd @@ -0,0 +1,92 @@ +--- +title: "Configure a MEMC Model" +--- + +```{r, include = FALSE, echo = FALSE, warning = FALSE, message = FALSE} +knitr::opts_chunk$set( + collapse = TRUE, + comment = "#>", + fig.width = 8, + fig.height = 5 +) +library(MEMC) +library(ggplot2) +theme_set(theme_bw(base_size = 14)) +``` + + +The `configure_model` function allows users to build their own SOM model. All of the MEMC models use the same SOM pool structure but users will be able to decide which of the following dynamics. + +![](pool_structure.png) + +The `configure_model` function allows users to build their own SOM model. All of the MEMC models use the same SOM pool structure but users will be able to decide which of the formulations (LM, MM, RMM, ECA, or DD) are used for POM decomposition (1), DOM uptake (2), and microbial biomass decay (3). The following example will demonstrate how changing the DOM uptake flux impacts the model projections. + +Set up the model configuration using the default paramaters and inital carbon pool values. + + + + +```{r} +# Running configure_model will print a table describing the model configuration. +mod1 <- configure_model(params = default_params, + state = default_initial, + name = "Model 1", + DOMuptake = "MM", + POMdecomp = "LM", + MBdecay = "LM") +``` + + +Use `solve_model` to run our model. + +```{r, warning=FALSE} +mod1_out <- solve_model(mod = mod1, time = 0:50) +head(mod1_out) +``` + + +Now switch the DOM uptake to use the RMM dynamics. + +```{r} +# Running configure_model will print a table describing the model configuration. +mod2 <- configure_model(params = default_params, + state = default_initial, + name = "Model 2", + DOMuptake = "RMM", + POMdecomp = "LM", + MBdecay = "LM") +mod2_out <- solve_model(mod = mod2, time = 0:50) +``` + + +Now try setting the DOM uptake to ECA dynamics. + + +```{r} +# Running configure_model will print a table describing the model configuration. +mod3 <- configure_model(params = default_params, + state = default_initial, + name = "Model 3", + DOMuptake = "ECA", + POMdecomp = "LM", + MBdecay = "LM") +mod3_out <- solve_model(mod = mod3, time = 0:50) +``` + + +Compare the model results with one another. + + +```{r, echo = FALSE} +ggplot(data = rbind(mod1_out, mod2_out, mod3_out)) + + geom_line(aes(time, value, color = name, linetype = name), linewidth = 0.75) + + facet_wrap("variable", scales = "free") + + labs(y = "mg C/g soil", + title = "Comparing Model Results") + + theme(legend.title = element_blank()) +``` + + +For more examples and tutorials please see our [online documentation](https://microbial-explicit-model.github.io/MEMC/). + + diff --git a/vignettes/articles/configs.png b/vignettes/articles/configs.png new file mode 100644 index 0000000..3c824ab Binary files /dev/null and b/vignettes/articles/configs.png differ diff --git a/vignettes/articles/pool_structure.png b/vignettes/articles/pool_structure.png new file mode 100644 index 0000000..393cb38 Binary files /dev/null and b/vignettes/articles/pool_structure.png differ diff --git a/vignettes/articles/roadmap.png b/vignettes/articles/roadmap.png new file mode 100644 index 0000000..7a0239b Binary files /dev/null and b/vignettes/articles/roadmap.png differ diff --git a/vignettes/articles/sens.Rmd b/vignettes/articles/sens.Rmd index 0790998..cfaa96e 100644 --- a/vignettes/articles/sens.Rmd +++ b/vignettes/articles/sens.Rmd @@ -22,9 +22,9 @@ theme_set(theme_bw()) -## Define the paaramter ranges +## Define the paramter ranges -Set up a data frame defining the parameter ranges to sample from. Any model parameters or initial pool sizes can be used, in this example we will explore the effects of maximum specific decomposition constants for POM, DOM, and MOM have on the model behvaior. Users should use more informed parameter ranges from the literature, here we sample a range of broad range defined as ± 50% of the parameter value. +Set up a data frame defining the parameter ranges to sample from. Any model parameters or initial pool sizes can be used, in this example we will explore the effects of maximum specific decomposition constants for POM, DOM, and MOM have on the model behavior. Users should use more informed parameter ranges from the literature, here we sample a range of broad range defined as ± 50% of the parameter value. ```{r} @@ -52,7 +52,7 @@ Description of inputs to MEMC function `memc_sensrange` are as follows: ```{r, fig.width=6} MENDsens_out <- memc_sensrange(config = MEND_model, t = seq(0, 365), - pars = pars, + x = pars, parRange = prange, dist = "latin", n = 50) @@ -67,11 +67,11 @@ head(summary(MENDsens_out)) Quickly visualize results using R's base `plot` function, this is a quick way to look at results however for more user control over the figures we recommend using ggplot2. ```{r, fig.width=8, fig.height=8} -plot(summary(MENDsens_out)) +plot(MENDsens_out) ``` ```{r, fig.width=8, fig.height=8} -plot(MENDsens_out) +plot(summary(MENDsens_out)) ``` @@ -88,7 +88,7 @@ Description of inputs to MEMC function `memc_sensfun` are as follows: ```{r} sens_out <- memc_sensfunc(config = MEND_model, t = seq(0, 365), - pars = c("V_d" = 3.0e+00,"V_p" = 1.4e+01,"V_m" = 2.5e-01)) + x = c("V_d" = 3.0e+00,"V_p" = 1.4e+01,"V_m" = 2.5e-01)) ``` Take a look at the results. @@ -124,37 +124,37 @@ Use the same parameter ranges with each of the different MEMC model configuratio time <- seq(0, 100) MENDsens_out <- memc_sensrange(config = MEND_model, t = time, - pars = pars, + x = pars, parRange = prange, dist = "latin", n = 50) COMISSIONsens_out <- memc_sensrange(config = COMISSION_model, t = time, - pars = pars, + x = pars, parRange = prange, dist = "latin", n = 50) CORPSEsens_out <- memc_sensrange(config = CORPSE_model, t = time, - pars = pars, + x = pars, parRange = prange, dist = "latin", n = 50) MEMSens_out <- memc_sensrange(config = MEMS_model, t = time, - pars = pars, + x = pars, parRange = prange, dist = "latin", n = 50) BAMSens_out <- memc_sensrange(config = BAMS_model, t = time, - pars = pars, + x = pars, parRange = prange, dist = "latin", n = 50) MIMCSens_out <- memc_sensrange(config = MIMCS_model, t = time, - pars = pars, + x = pars, parRange = prange, dist = "latin", n = 50) @@ -196,12 +196,12 @@ Use the `memc_sensfun` to compare the sensitivity of the DOM pool is to V_d and time <- 0:100 MENDsens_out <- memc_sensfunc(config = MEND_model, t = time, - pars = c("V_d" = 3.0e+00,"K_d" = 0.250), + x = c("V_d" = 3.0e+00,"K_d" = 0.250), sensvar = "DOM") CORPSEsens_out <- memc_sensfunc(config = CORPSE_model, t = time, - pars = c("V_d" = 3.0e+00,"K_d" = 0.250), + x = c("V_d" = 3.0e+00,"K_d" = 0.250), sensvar = "DOM") ``` @@ -216,9 +216,39 @@ ggplot(data = out) + geom_line(aes(time, value, color = name)) + facet_wrap("parameter", scales = "free") + labs(title = "Local Senstivity to DOM", - x = "Time", y = "" ) + x = "Time", y = "" ) + + scale_color_manual(values = MEMC::colorMEMCPalette(out$model)) +``` + + + + + + +# Model Sensitivity to Inital Pool Sizes + +The MEMC models are also sensitive to the initial pool size! Here we will demonstrate how model results are sensitive to uncertainty within a single pool! This being said the initial carbon pool sizes matter and users will need to make sure that they are using the appropriate values for their applications. + +```{r} +pools <- c("DOM" = 1) +prange <- data.frame(pools = pools, + min = pools * 0.50, + max = pools * 1.50) ``` +```{r, fig.width=6} +inital_pool_uncertainty <- memc_sensrange(config = MEND_model, + t = seq(0, 365), + x = pools, + parRange = prange, + dist = "latin", + n = 50) +``` + + +```{r} +plot(summary(inital_pool_uncertainty)) +``` diff --git a/vignettes/articles/som_diversity.png b/vignettes/articles/som_diversity.png new file mode 100644 index 0000000..5b77796 Binary files /dev/null and b/vignettes/articles/som_diversity.png differ