library(plant)Strategies & traits
guide
How to inspect and modify the parameters that define a plant strategy, and how trait values propagate through the hyper-parameterisation. For what the strategies mean physiologically, see the model docs.
There are a large number of parameters to the physiological model, but all are changeable. The default strategy is detailed in the model docs:
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" "k_I" "recruitment_decay" "control" "collect_all_auxiliary" "birth_rate_x"
[36] "birth_rate_y" "is_variable_birth_rate"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 theta a_l1 a_l2
1.978791e-01 6.080000e+02 1.659587e+01 3.800000e-05 1.200000e+01 2.141786e-04 5.440000e+00 3.060000e-01
a_r1 a_b1 r_s r_b r_r r_l a_y a_bio
7.000000e-02 1.700000e-01 6.598684e+00 1.319737e+01 2.170000e+02 1.984545e+02 7.000000e-01 2.450000e-02
k_l k_b k_s k_r a_p1 a_p2 a_f3 a_f1
4.565855e-01 2.000000e-01 2.000000e-01 1.000000e+00 1.511778e+02 2.047162e-01 1.140000e-04 1.000000e+00
a_f2 S_D a_d0 d_I a_dG1 a_dG2 k_I recruitment_decay
5.000000e+01 2.500000e-01 1.000000e-01 1.000000e-02 5.500000e+00 2.000000e+01 5.000000e-01 0.000000e+00
collect_all_auxiliary birth_rate_y is_variable_birth_rate
0.000000e+00 1.000000e+00 0.000000e+00 All 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_Individual(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_hyperparfunction (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_I <- s$k_I
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(dplyr::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) {
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 (any(is.infinite(extra))) {
stop("Attempt to use infinite value in derived parameters: ",
paste(colnames(extra)[is.infinite(extra)], 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: 0x8b95b32d8>
<environment: 0x8b95ae2a8>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 a_p1 r_l
p1 0.1 1.466783 151.1778 392.7These are:
k_l: Turnover rate for leavesa_p1: Leaf photosynthesis per arear_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] "patch_area" "n_patches" "patch_type" "max_patch_lifetime" "strategies" "strategy_default"
[7] "node_schedule_times_default" "node_schedule_times" "ode_times" and it comes pre-set with the hyper-parameterisation function:
identical(p$hyperpar, FF16_hyperpar)[1] FALSEk_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_list(trait_matrix(0.1, "lma"), p, birth_rate_list = 1)[[1]]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, birth_rate_list = rep(1, 5))
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.540In 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$controlNULLThe Control() defaults are the pragmatic, fast-ish settings used for essentially all of plant’s runs (control() is a lowercase alias for the same constructor). When you need a high-accuracy run, control_accurate() tightens the ODE and schedule tolerances. We can see which fields it changes:
ctrl_accurate <- control_accurate()
ctrl_accurate[unlist(control()) != unlist(ctrl_accurate)]$ode_step_size_max
[1] 0.1
$ode_tol_rel
[1] 1e-06
$ode_tol_abs
[1] 1e-06
$schedule_eps
[1] 0.001