library(plant)Strategies & traits
guide
A strategy contains the biological settings for one type of plant. This guide shows two ways to create one: change a parameter directly, or provide measured trait values and let plant update the linked parameters. For the equations behind each strategy, see Models.
Inspect or directly change a strategy
Start with the default FF16 strategy. Printing its names shows that the object contains biological parameters, numerical controls, recruitment settings, and other model information:
s <- FF16_Strategy()
names(s)[1] "pars" "control" "collect_all_auxiliary"
[4] "birth_rate_x" "birth_rate_y" "is_variable_birth_rate"The physiological parameters live in the nested pars element. They are ordinary decimal numbers and can be inspected with familiar R tools:
unlist(s$pars) lma rho hmat omega
1.978791e-01 6.080000e+02 1.659587e+01 3.800000e-05
eta theta a_l1 a_l2
1.200000e+01 2.141786e-04 5.440000e+00 3.060000e-01
a_r1 a_b1 r_s r_b
7.000000e-02 1.700000e-01 6.598684e+00 1.319737e+01
r_r r_l a_y a_bio
2.170000e+02 1.984545e+02 7.000000e-01 2.450000e-02
k_l k_b k_s k_r
4.565855e-01 2.000000e-01 2.000000e-01 1.000000e+00
a_p1 a_p2 a_f3 a_f1
1.511778e+02 2.047162e-01 1.140000e-04 1.000000e+00
a_f2 S_D a_d0 d_I
5.000000e+01 2.500000e-01 1.000000e-01 1.000000e-02
a_dG1 a_dG2 k_I recruitment_decay
5.500000e+00 2.000000e+01 5.000000e-01 0.000000e+00 For a controlled experiment, you can change one value directly. Here we double leaf mass per unit leaf area (lma):
s$pars$lma <- s$pars$lma * 2An individual constructed from this strategy uses the revised value:
pl <- FF16_Individual(s)
pl$strategy$pars$lma[1] 0.3957582This direct edit changes only lma. That is useful when deliberately varying one model parameter, but biological traits are often linked. In FF16, LMA affects leaf mass and seedling allometry directly and is also associated with leaf turnover and respiration.
Create strategies from traits
plant represents these links with a hyperparameterisation: an R function that translates a small set of input traits into a consistent set of model parameters. FF16 uses FF16_hyperpar():
FF16_hyperparfunction (m, s, filter = TRUE)
{
with_default <- function(name, default_value = s$pars[[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$pars$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$pars[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: 0x139ff9c08>
<environment: 0x139ff14e8>You will rarely call this function directly. Looking at its output is useful, however, because it makes the parameter links explicit. Supplying LMA updates leaf turnover (k_l) and leaf respiration per mass (r_l). The output also includes photosynthesis per leaf area (a_p1), which is derived from leaf nitrogen per area (narea) rather than LMA itself.
FF16_hyperpar(trait_matrix(0.1, "lma"), s) lma k_l a_p1 r_l
p1 0.1 1.466783 151.1778 392.7The three returned quantities are:
k_l: Turnover rate for leavesa_p1: Leaf photosynthesis per arear_l: Leaf respiration per mass
Because these linking rules live in R, advanced users can define a different trait trade-off without changing the underlying physiological model in C++.
trait_matrix() puts one or more trait values into the row-and-column format expected by the strategy-building functions:
trait_matrix(0.1, "lma") lma
[1,] 0.1trait_matrix(c(0.1, 500), "rho") rho
[1,] 0.1
[2,] 500.0The preferred workflow uses generate_strategy() or add_strategies(). Both start from a Parameters object, which collects the settings needed for a patch run and stores the default strategy from which new strategies are built:
p <- FF16_Parameters()Its top-level names show the main groups of settings:
names(p) [1] "patch_area" "n_patches"
[3] "patch_type" "max_patch_lifetime"
[5] "strategies" "strategy_default"
[7] "node_schedule_times_default" "node_schedule_times"
[9] "ode_times" "initial_state"
[11] "n_initial_cohorts" "initial_node_times"
[13] "initial_patch_density" "initial_pr_patch_survival"
[15] "initial_time" param_hyperpar() retrieves the linking function associated with the model:
identical(param_hyperpar(p), FF16_hyperpar)[1] TRUEThe strategy_default member is the template copied when new strategies are generated. Physiological parameters, including the light extinction coefficient k_I, live under its pars element:
class(p$strategy_default)[1] "FF16_Strategy"Generating a strategy copies that template, sets LMA, and updates the linked parameters:
s <- generate_strategy(p, trait_matrix(0.1, "lma"), birth_rate = 1)[[1]]The same function can generate several strategies at once. Each row of the trait matrix becomes one strategy, so the birth-rate vector also has one value per row:
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 <- generate_strategy(p, lma, birth_rate = rep(1, 5))
length(ss)[1] 5Standard R functions can compare the resulting parameter values:
sapply(ss, function(x) x$pars$lma)[1] 0.1 0.2 0.3 0.4 0.5sapply(ss, function(x) x$pars$k_l)[1] 1.46678341 0.44833712 0.22412414 0.13703875 0.09356799sapply(ss, function(x) x$pars$a_p1)[1] 151.1778 151.1778 151.1778 151.1778 151.1778sapply(ss, function(x) x$pars$r_l)[1] 392.700 196.350 130.900 98.175 78.540Choose numerical accuracy
Strategies describe plants; Control settings describe how the numerical solver runs. The controls include time-step, error-tolerance, and node-schedule settings, and may be numeric or logical:
Control()$function_integration_rule
[1] 21
$shading_model
[1] ""
$ppa_layer_optical_depth
[1] 0.5
$ppa_layer_smoothing
[1] 0.3
$offspring_production_tol
[1] 1e-08
$offspring_production_iterations
[1] 1000
$node_gradient_eps
[1] 1e-06
$node_gradient_direction
[1] -1
$node_gradient_richardson
[1] FALSE
$node_gradient_richardson_depth
[1] 4
$ode_step_size_initial
[1] 1e-06
$ode_step_size_min
[1] 1e-06
$ode_step_size_max
[1] 5
$ode_tol_rel
[1] 1e-04
$ode_tol_abs
[1] 1e-04
$ode_a_y
[1] 1
$ode_a_dydt
[1] 0
$fixed_time_step
[1] 0
$schedule_nsteps
[1] 20
$schedule_eps
[1] 0.02
$schedule_verbose
[1] FALSE
$save_RK45_cache
[1] FALSE
$GSS_tol_abs
[1] 0.001
$vulnerability_curve_ncontrol
[1] 100
$ci_abs_tol
[1] 0.001
$ci_niter
[1] 1000
attr(,"class")
[1] "Control"The Control() defaults are appropriate for routine exploration; control() is a lowercase alias for the same constructor. Use control_accurate() when a final result needs tighter numerical tolerances, accepting the additional run time. The ODE stepping and node-spacing pages explain these trade-offs. The following comparison shows which defaults change:
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