Add time varying growth and mortality

R/ff16drivers.R

@@ -0,0 +1,308 @@
+# Built from  R/ff16.R on Wed Aug 12 11:12:34 2020 using the scaffolder, from the strategy:  FF16
+
+##' Create a FF16drivers Plant or Node
+##' @title Create a FF16drivers Plant or Node
+##' @param s A \code{\link{FF16drivers_Strategy}} object
+##' @export
+##' @rdname FF16drivers
+##' @examples
+##' pl <- FF16drivers_Individual()
+##' pl$height
+FF16drivers_Individual <- function(s=FF16drivers_Strategy()) {
+  Individual("FF16drivers", "FF16_Env")(s)
+}
+
+##' @export
+##' @rdname FF16drivers
+FF16drivers_Node <- function(s=FF16drivers_Strategy()) {
+  Node("FF16drivers", "FF16_Env")(s)
+}
+
+##' @export
+##' @rdname FF16drivers
+FF16drivers_Species <- function(s=FF16drivers_Strategy()) {
+  Species("FF16drivers", "FF16_Env")(s)
+}
+
+##' @export
+##' @rdname FF16drivers
+FF16drivers_Parameters <- function() {
+  Parameters("FF16drivers","FF16_Env")()
+}
+
+##' @export
+##' @rdname FF16drivers
+##' @param p A \code{Parameters<FF16drivers,FF16_Env>} object
+FF16drivers_Patch <- function(p) {
+  Patch("FF16drivers", "FF16_Env")(p)
+}
+
+##' @export
+##' @rdname FF16drivers
+FF16drivers_SCM <- function(p) {
+  SCM("FF16drivers", "FF16_Env")(p)
+}
+
+##' @export
+##' @rdname FF16drivers
+FF16drivers_StochasticSpecies <- function(s=FF16drivers_Strategy()) {
+  StochasticSpecies("FF16drivers", "FF16_Env")(s)
+}
+
+##' @export
+##' @rdname FF16drivers
+FF16drivers_StochasticPatch <- function(p) {
+  StochasticPatch("FF16drivers", "FF16_Env")(p)
+}
+
+##' @export
+##' @rdname FF16drivers
+FF16drivers_StochasticPatchRunner <- function(p) {
+  StochasticPatchRunner("FF16drivers", "FF16_Env")(p)
+}
+
+## Helper:
+##' @export
+##' @rdname FF16_Environment
+##' @param canopy_rescale_usually turns on environment rescaling (default = TRUE)
+FF16drivers_make_environment <- function(canopy_rescale_usually = TRUE) {
+
+  FF16_make_environment(canopy_rescale_usually)
+}
+
+
+##' This makes a pretend light environment over the plant height,
+##' slightly concave up, whatever.
+##' @title Create a test environment for FF16drivers startegy
+##' @param height top height of environment object
+##' @param n number of points
+##' @param light_env function for light environment in test object
+##' @param n_strategies number of strategies for test environment
+##' @export
+##' @rdname FF16drivers_test_environment
+##' @examples
+##' environment <- FF16drivers_test_environment(10)
+FF16drivers_test_environment <- function(height, n=101, light_env=NULL,
+                             n_strategies=1) {
+  hh <- seq(0, height, length.out=n)
+  if (is.null(light_env)) {
+    light_env <- function(x) {
+      exp(x/(height*2)) - 1 + (1 - (exp(.5) - 1))/2
+    }
+  }
+  ee <- light_env(hh)
+  interpolator <- Interpolator()
+  interpolator$init(hh, ee)
+
+  ret <- FF16drivers_make_environment()
+  ret$canopy$canopy_interpolator <- interpolator
+  attr(ret, "light_env") <- light_env
+  ret
+}
+
+
+
+##' Hyperparameters for FF16drivers physiological model
+##' @title Hyperparameters for FF16drivers physiological model
+##' @param lma_0 Central (mean) value for leaf mass per area [kg /m2]
+##' @param B_kl1 Rate of leaf turnover at lma_0 [/yr]
+##' @param B_kl2 Scaling slope for phi in leaf turnover [dimensionless]
+##' @param rho_0 Central (mean) value for wood density [kg /m3]
+##' @param B_dI1 Rate of instantaneous mortality at rho_0 [/yr]
+##' @param B_dI2 Scaling slope for wood density in intrinsic mortality [dimensionless]
+##' @param B_ks1 Rate of sapwood turnover at rho_0 [/yr]
+##' @param B_ks2 Scaling slope for rho in sapwood turnover [dimensionless]
+##' @param B_rs1 CO_2 respiration per unit sapwood volume [mol / yr / m3 ]
+##' @param B_rb1 CO_2 respiration per unit sapwood volume [mol / yr / m3 ]
+##' @param B_f1 Cost of seed accessories per unit seed mass [dimensionless]
+##' @param narea nitrogen per leaf area [kg / m2]
+##' @param narea_0 central (mean) value for nitrogen per leaf area [kg / m2]
+##' @param B_lf1 Potential CO_2 photosynthesis at average leaf nitrogen [mol / d / m2]
+##' @param B_lf2 Curvature of leaf photosynthetic light response curve [dimensionless]
+##' @param B_lf3 Quantum yield of leaf photosynthetic light response curve [dimensionless]
+##' @param B_lf4 CO_2 respiration per unit leaf nitrogen [mol / yr / kg]
+##' @param B_lf5 Scaling exponent for leaf nitrogen in maximum leaf photosynthesis [dimensionless]
+##' @param k_I light extinction coefficient [dimensionless]
+##' @param latitude degrees from equator (0-90), used in solar model [deg]
+##' @importFrom stats coef nls
+##' @export
+##' @rdname FF16drivers_hyperpar
+make_FF16drivers_hyperpar <- function(
+                                lma_0=0.1978791,
+                                B_kl1=0.4565855,
+                                B_kl2=1.71,
+                                rho_0=608.0,
+                                B_dI1=0.01,
+                                B_dI2=0.0,
+                                B_ks1=0.2,
+                                B_ks2=0.0,
+                                B_rs1=4012.0,
+                                B_rb1=2.0*4012.0,
+                                B_f1 =3.0,
+                                narea=1.87e-3,
+                                narea_0=1.87e-3,
+                                B_lf1=5120.738 * 1.87e-3 * 24 * 3600 / 1e+06,
+                                B_lf2=0.5,
+                                B_lf3=0.04,
+                                B_lf4=21000,
+                                B_lf5=1,
+                                k_I=0.5,
+                                latitude=0) {
+  assert_scalar <- function(x, name=deparse(substitute(x))) {
+    if (length(x) != 1L) {
+      stop(sprintf("%s must be a scalar", name), call. = FALSE)
+    }
+  }
+  assert_scalar(lma_0)
+  assert_scalar(B_kl1)
+  assert_scalar(B_kl2)
+  assert_scalar(rho_0)
+  assert_scalar(B_dI1)
+  assert_scalar(B_dI2)
+  assert_scalar(B_ks1)
+  assert_scalar(B_ks2)
+  assert_scalar(B_rs1)
+  assert_scalar(B_rb1)
+  assert_scalar(B_f1)
+  assert_scalar(narea)
+  assert_scalar(narea_0)
+  assert_scalar(B_lf1)
+  assert_scalar(B_lf2)
+  assert_scalar(B_lf3)
+  assert_scalar(B_lf4)
+  assert_scalar(B_lf5)
+  assert_scalar(k_I)
+  assert_scalar(latitude)
+
+  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)
+
+    ## lma / leaf turnover relationship:
+    k_l   <- B_kl1 * (lma / lma_0) ^ (-B_kl2)
+
+    ## rho / mortality relationship:
+    d_I  <- B_dI1 * (rho / rho_0) ^ (-B_dI2)
+
+    ## rho / wood turnover relationship:
+    k_s  <- B_ks1 *  (rho / rho_0) ^ (-B_ks2)
+
+    ## rho / sapwood respiration relationship:
+
+    ## Respiration rates are per unit mass, so this next line has the
+    ## effect of holding constant the respiration rate per unit volume.
+    ## So respiration rates per unit mass vary with rho, respiration
+    ## rates per unit volume don't.
+    r_s <- B_rs1 / rho
+    # bark respiration follows from sapwood
+    r_b <- B_rb1 / rho
+
+    ## omega / accessory cost relationship
+    a_f3 <- B_f1 * omega
+
+    ## Narea, photosynthesis, respiration
+
+    assimilation_rectangular_hyperbolae <- function(I, Amax, theta, QY) {
+      x <- QY * I + Amax
+      (x - sqrt(x^2 - 4 * theta * QY * I * Amax)) / (2 * theta)
+    }
+
+    ## Photosynthesis  [mol CO2 / m2 / yr]
+    approximate_annual_assimilation <- function(narea, latitude) {
+      E <- seq(0, 1, by=0.02)
+      ## Only integrate over half year, as solar path is symmetrical
+      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-8)) {
+        # line fitting will fail if all have are zero, or potentially same value
+        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
+    }
+
+    # This needed in case narea has length zero, in which case trapezium fails
+    a_p1 <- a_p2 <- 0 * narea
+    ## TODO: Remove the 0.5 hardcoded default for k_I here, and deal
+    ## with this more nicely.
+    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]
+    }
+
+    ## Respiration rates are per unit mass, so convert to mass-based
+    ## rate by dividing with lma
+    ## So respiration rates per unit mass vary with lma, while
+    ## respiration rates per unit area don't.
+    r_l  <- B_lf4 * narea / lma
+
+    extra <- cbind(k_l,                # lma
+                   d_I, k_s, r_s, r_b, # rho
+                   a_f3,               # omega
+                   a_p1, a_p2,         # narea
+                   r_l)                # lma, narea
+
+    overlap <- intersect(colnames(m), colnames(extra))
+    if (length(overlap) > 0L) {
+      stop("Attempt to overwrite generated parameters: ",
+           paste(overlap, collapse=", "))
+    }
+
+    ## Filter extra so that any column where all numbers are with eps
+    ## of the default strategy are not replaced:
+    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
+  }
+}
+
+##' Hyperparameter function for FF16drivers physiological model
+##' @title Hyperparameter function for FF16drivers physiological model
+##' @param m A matrix of trait values, as returned by \code{trait_matrix}
+##' @param s A strategy object
+##' @param filter A flag indicating whether to filter columns. If TRUE, any numbers
+##' that are within eps of the default strategy are not replaced.
+##' @export
+FF16drivers_hyperpar <- make_FF16drivers_hyperpar()

R/strategy_support.R

@@ -11,6 +11,7 @@ make_hyperpar <- function(type) {
          FF16=make_FF16_hyperpar,
          FF16w=make_FF16w_hyperpar,
          FF16r=make_FF16r_hyperpar,
+         FF16drivers=make_FF16drivers_hyperpar,
          K93=make_K93_hyperpar,
          stop("Unknown type ", type))
 }
@@ -23,6 +24,7 @@ param_hyperpar <- function(parameters) {
          FF16_Strategy=FF16_hyperpar,
          FF16w_Strategy=FF16w_hyperpar,
          FF16r_Strategy=FF16r_hyperpar,
+         FF16drivers_Strategy=FF16drivers_hyperpar,
          K93_Strategy=K93_hyperpar,
          stop("Unknown type ", type))
 }
@@ -37,6 +39,7 @@ hyperpar <- function(type) {
          FF16=FF16_hyperpar,
          FF16w=FF16w_hyperpar,
          FF16r=FF16r_hyperpar,
+         FF16drivers=FF16drivers_hyperpar,
          K93=K93_hyperpar,
          stop("Unknown type ", type))
 }
@@ -48,6 +51,7 @@ environment_type <- function(type) {
          FF16=sprintf("FF16_Env"),
          FF16w=sprintf("FF16_Env"),
          FF16r=sprintf("FF16_Env"),
+         FF16drivers=sprintf("FF16_Env"),
          K93=sprintf("K93_Env"),
          stop("Unknown type ", type))
 }
@@ -69,6 +73,7 @@ make_environment <- function(type = NULL, parameters = NULL, ...) {
          FF16=FF16_make_environment(...),
          FF16w=FF16w_make_environment(...),
          FF16r=FF16r_make_environment(...),
+         FF16drivers=FF16drivers_make_environment(...),
          K93=K93_make_environment(...),
          stop("Unknown type ", type))
 }
@@ -79,6 +84,7 @@ node_schedule_default <- function(p) {
          "Parameters<FF16,FF16_Env>"=`node_schedule_default__Parameters___FF16__FF16_Env`,
          "Parameters<FF16w,FF16_Env>"=`node_schedule_default__Parameters___FF16w__FF16_Env`,
          "Parameters<FF16r,FF16_Env>"=`node_schedule_default__Parameters___FF16r__FF16_Env`,
+         "Parameters<FF16drivers,FF16_Env>"=`node_schedule_default__Parameters___FF16drivers__FF16_Env`,
          "Parameters<K93,K93_Env>"=`node_schedule_default__Parameters___K93__K93_Env`,
          stop("Unknown type: ", cl))(p)
 }
@@ -89,6 +95,7 @@ make_node_schedule <- function(p) {
          "Parameters<FF16,FF16_Env>"=`make_node_schedule__Parameters___FF16__FF16_Env`,
          "Parameters<FF16w,FF16_Env>"=`make_node_schedule__Parameters___FF16w__FF16_Env`,
          "Parameters<FF16r,FF16_Env>"=`make_node_schedule__Parameters___FF16r__FF16_Env`,
+         "Parameters<FF16drivers,FF16_Env>"=`make_node_schedule__Parameters___FF16drivers__FF16_Env`,
          "Parameters<K93,K93_Env>"=`make_node_schedule__Parameters___K93__K93_Env`,
                   stop("Unknown type: ", cl))(p)
 }