Add Documentation (#17)

* Tweak the FME family of functions to work for inital pool sizes and paramters

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, ...)
 
 }
 

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"]])
 
   }
 

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)
 

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            <NA>   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 

_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

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

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", "")

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
   )
 

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))
 })

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)