library(plant)There are a large number of parameters to the physiological model, but all are changeable. The default strategy is detailed in the “physiology” vignette:
s <- FF16_Strategy()
names(s)## [1] "lma" "rho" "hmat" "omega" "eta" "theta" "a_l1"
## [8] "a_l2" "a_r1" "a_b1" "r_s" "r_b" "r_r" "r_l"
## [15] "a_y" "a_bio" "k_l" "k_b" "k_s" "k_r" "a_p1"
## [22] "a_p2" "a_f3" "a_f1" "a_f2" "S_D" "a_d0" "d_I"
## [29] "a_dG1" "a_dG2" "control"The strategy object here is a special object of class FF16_Strategy. In contrast with most of the “reference” objects used by plant, this is simply a list with a class attribute. Some validation will be done on the parameters every time it is passed into the C++ bits of the model.
All parameters except for control are floating point (decimal) numbers:
unlist(s[names(s) != "control"])## lma rho hmat omega eta
## 1.978791e-01 6.080000e+02 1.659587e+01 3.800000e-05 1.200000e+01
## theta a_l1 a_l2 a_r1 a_b1
## 2.141786e-04 5.440000e+00 3.060000e-01 7.000000e-02 1.700000e-01
## r_s r_b r_r r_l a_y
## 6.598684e+00 1.319737e+01 2.170000e+02 1.984545e+02 7.000000e-01
## a_bio k_l k_b k_s k_r
## 2.450000e-02 4.565855e-01 2.000000e-01 2.000000e-01 1.000000e+00
## a_p1 a_p2 a_f3 a_f1 a_f2
## 1.511778e+02 2.047162e-01 1.140000e-04 1.000000e+00 5.000000e+01
## S_D a_d0 d_I a_dG1 a_dG2
## 2.500000e-01 1.000000e-01 1.000000e-02 5.500000e+00 2.000000e+01All of these parameters can be directly changed. For example, we can double the leaf mass per unit leaf area (lma) value:
s$lma <- s$lma * 2and then from this construct a plant:
pl <- FF16_Plant(s)
pl$strategy$lma## [1] 0.3957582lma affects a few places in the model; see the source code, but only really for converting from leaf area to leaf mass (leaf mass being leaf area multiplied by leaf mass per unit leaf area). However, as a component of the leaf economic spectrum we imagine LMA as affecting a number of other components of the model.
We capture this through what we call a “hyper-parameterisation”; additional parameters and functions that mean that changing one parameter might affect a number of other lower-level parameters. Our default hyper-parameterisation is via the FF16_hyperpar function:
FF16_hyperpar## function (m, s, filter = TRUE)
## {
## with_default <- function(name, default_value = s[[name]]) {
## rep_len(if (name %in% colnames(m))
## m[, name]
## else default_value, nrow(m))
## }
## lma <- with_default("lma")
## rho <- with_default("rho")
## omega <- with_default("omega")
## narea <- with_default("narea", narea)
## k_l <- B_kl1 * (lma/lma_0)^(-B_kl2)
## d_I <- B_dI1 * (rho/rho_0)^(-B_dI2)
## k_s <- B_ks1 * (rho/rho_0)^(-B_ks2)
## r_s <- B_rs1/rho
## r_b <- B_rb1/rho
## a_f3 <- B_f1 * omega
## assimilation_rectangular_hyperbolae <- function(I, Amax,
## theta, QY) {
## x <- QY * I + Amax
## (x - sqrt(x^2 - 4 * theta * QY * I * Amax))/(2 * theta)
## }
## approximate_annual_assimilation <- function(narea, latitude) {
## E <- seq(0, 1, by = 0.02)
## D <- seq(0, 365/2, length.out = 10000)
## I <- PAR_given_solar_angle(solar_angle(D, latitude = abs(latitude)))
## Amax <- B_lf1 * (narea/narea_0)^B_lf5
## theta <- B_lf2
## QY <- B_lf3
## AA <- NA * E
## for (i in seq_len(length(E))) {
## AA[i] <- 2 * trapezium(D, assimilation_rectangular_hyperbolae(k_I *
## I * E[i], Amax, theta, QY))
## }
## if (all(diff(AA) < 1e-08)) {
## ret <- c(last(AA), 0)
## names(ret) <- c("p1", "p2")
## }
## else {
## fit <- nls(AA ~ p1 * E/(p2 + E), data.frame(E = E,
## AA = AA), start = list(p1 = 100, p2 = 0.2))
## ret <- coef(fit)
## }
## ret
## }
## a_p1 <- a_p2 <- 0 * narea
## if (length(narea) > 0 || k_I != 0.5) {
## i <- match(narea, unique(narea))
## y <- vapply(unique(narea), approximate_annual_assimilation,
## numeric(2), latitude)
## a_p1 <- y["p1", i]
## a_p2 <- y["p2", i]
## }
## r_l <- B_lf4 * narea/lma
## extra <- cbind(k_l, d_I, k_s, r_s, r_b, a_f3, a_p1, a_p2,
## r_l)
## overlap <- intersect(colnames(m), colnames(extra))
## if (length(overlap) > 0L) {
## stop("Attempt to overwrite generated parameters: ", paste(overlap,
## collapse = ", "))
## }
## if (filter) {
## if (nrow(extra) == 0L) {
## extra <- NULL
## }
## else {
## pos <- diff(apply(extra, 2, range)) == 0
## if (any(pos)) {
## eps <- sqrt(.Machine$double.eps)
## x1 <- extra[1, pos]
## x2 <- unlist(s[names(x1)])
## drop <- abs(x1 - x2) < eps & abs(1 - x1/x2) <
## eps
## if (any(drop)) {
## keep <- setdiff(colnames(extra), names(drop)[drop])
## extra <- extra[, keep, drop = FALSE]
## }
## }
## }
## }
## if (!is.null(extra)) {
## m <- cbind(m, extra)
## }
## m
## }
## <bytecode: 0x7fc970721590>
## <environment: 0x7fc97060dc38>You will rarely need to call this function directly (see below) but note how setting of lma affects parameters k_l, a_p1 and r_l
FF16_hyperpar(trait_matrix(0.1, "lma"), s)## lma k_l r_l
## p1 0.1 1.466783 392.7These are: * k_l: Turnover rate for leaves * a_p1: Leaf photosynthesis per area * r_l: Leaf respiration per mass
This means that it is possible to implement different trade-offs between parameters relatively easily by modifying the hyper-parameterisation function in R, rather than having to modify the underlying physiological model in C++.
The trait_matrix function is a simple wrapper that just makes sure the trait matrix has the right format:
trait_matrix(0.1, "lma")## lma
## [1,] 0.1trait_matrix(c(0.1, 500), "rho")## rho
## [1,] 0.1
## [2,] 500.0To make use of the hyper-parameterisation, the preferred way of setting parameters is through the utility functions strategy and strategy_list. These take a Parameters object:
p <- FF16_Parameters()The Parameters object mostly contains information about the patch:
names(p)## [1] "k_I" "patch_area"
## [3] "n_patches" "disturbance_mean_interval"
## [5] "strategies" "seed_rain"
## [7] "is_resident" "control"
## [9] "strategy_default" "cohort_schedule_max_time"
## [11] "cohort_schedule_times_default" "cohort_schedule_times"
## [13] "cohort_schedule_ode_times" "hyperpar"and it comes pre-set with the hyper-parameterisation function:
identical(p$hyperpar, FF16_hyperpar)## [1] TRUEk_I is the light extinction coefficient, disturbance_mean_interval is the mean disturbance interval. The strategy_default member is a Strategy:
class(p$strategy_default)## [1] "FF16_Strategy"This is the Strategy object that all others will be built from by difference. Running
s <- strategy(trait_matrix(0.1, "lma"), p)will create a strategy s where lma is set but also all the parameters that depend on lma.
The function strategy_list can be used to create a list of strategies:
lma <- trait_matrix(seq(0.1, 0.5, length.out=5), "lma")
FF16_hyperpar(trait_matrix(lma, "lma"), s)## lma k_l r_l
## [1,] 0.1 1.46678341 392.700
## [2,] 0.2 0.44833712 196.350
## [3,] 0.3 0.22412414 130.900
## [4,] 0.4 0.13703875 98.175
## [5,] 0.5 0.09356799 78.540ss <- strategy_list(lma, p)
length(ss)## [1] 5We can then use standard R commands to extract variable from this list
sapply(ss, function(x) x$lma)## [1] 0.1 0.2 0.3 0.4 0.5sapply(ss, function(x) x$k_l)## [1] 1.46678341 0.44833712 0.22412414 0.13703875 0.09356799sapply(ss, function(x) x$a_p1)## [1] 151.1778 151.1778 151.1778 151.1778 151.1778sapply(ss, function(x) x$r_l)## [1] 392.700 196.350 130.900 98.175 78.540There’s a convenience function plant_list that returns a set of plants based on a vector of traits:
pp <- plant_list(lma, p)
sapply(pp, function(p) p$area_leaf_above(0))## [1] 1.713949e-04 1.201090e-04 9.206433e-05 7.449914e-05 6.250239e-05In addition to the physiological parameters there are large number of “control” parameters that affect the behaviour of the various numerical algorithms used (note that in contrast to the physiological parameters these have a variety of types)
p$control## $plant_assimilation_adaptive
## [1] TRUE
##
## $plant_assimilation_over_distribution
## [1] FALSE
##
## $plant_assimilation_tol
## [1] 1e-06
##
## $plant_assimilation_iterations
## [1] 1000
##
## $plant_assimilation_rule
## [1] 21
##
## $plant_seed_tol
## [1] 1e-08
##
## $plant_seed_iterations
## [1] 1000
##
## $cohort_gradient_eps
## [1] 1e-06
##
## $cohort_gradient_direction
## [1] 1
##
## $cohort_gradient_richardson
## [1] FALSE
##
## $cohort_gradient_richardson_depth
## [1] 4
##
## $environment_light_tol
## [1] 1e-06
##
## $environment_light_nbase
## [1] 17
##
## $environment_light_max_depth
## [1] 16
##
## $environment_light_rescale_usually
## [1] FALSE
##
## $ode_step_size_initial
## [1] 1e-06
##
## $ode_step_size_min
## [1] 1e-06
##
## $ode_step_size_max
## [1] 0.1
##
## $ode_tol_rel
## [1] 1e-06
##
## $ode_tol_abs
## [1] 1e-06
##
## $ode_a_y
## [1] 1
##
## $ode_a_dydt
## [1] 0
##
## $schedule_nsteps
## [1] 20
##
## $schedule_eps
## [1] 0.001
##
## $schedule_verbose
## [1] FALSE
##
## $schedule_patch_survival
## [1] 6.253026e-05
##
## $equilibrium_nsteps
## [1] 20
##
## $equilibrium_eps
## [1] 1e-05
##
## $equilibrium_large_seed_rain_change
## [1] 10
##
## $equilibrium_verbose
## [1] TRUE
##
## $equilibrium_solver_name
## [1] "iteration"
##
## $equilibrium_extinct_seed_rain
## [1] 0.001
##
## $equilibrium_nattempts
## [1] 5
##
## $equilibrium_solver_logN
## [1] TRUE
##
## $equilibrium_solver_try_keep
## [1] TRUE
##
## attr(,"class")
## [1] "Control"The defaults are rather too slow for many uses, so scm_base_parameters provides a faster set by using fast_control to set many of these to less accurate values.
p2 <- scm_base_parameters("FF16")
p2$control[unlist(p$control) != unlist(p2$control)]## $plant_assimilation_adaptive
## [1] FALSE
##
## $plant_assimilation_tol
## [1] 1e-04
##
## $cohort_gradient_direction
## [1] -1
##
## $environment_light_tol
## [1] 1e-04
##
## $environment_light_rescale_usually
## [1] TRUE
##
## $ode_step_size_max
## [1] 5
##
## $ode_tol_rel
## [1] 1e-04
##
## $ode_tol_abs
## [1] 1e-04
##
## $schedule_eps
## [1] 0.005
##
## $schedule_verbose
## [1] TRUE
##
## $equilibrium_eps
## [1] 0.001