From 5b7b90ee349df26fac2fcc54d99f4b94d372211f Mon Sep 17 00:00:00 2001 From: Aiyu Zheng <46060642+aiyuz@users.noreply.github.com> Date: Fri, 1 May 2020 21:27:56 -0400 Subject: [PATCH] Add files via upload --- DESCRIPTION | 31 +++ NAMESPACE | 2 + R/methods.R | 403 ++++++++++++++++++++++++++++++ R/model_clonal.R | 121 +++++++++ R/model_light.R | 209 ++++++++++++++++ R/utility_functions.R | 553 ++++++++++++++++++++++++++++++++++++++++++ man/calcPPA.Rd | 14 ++ vignettes/example.R | 54 +++++ 8 files changed, 1387 insertions(+) create mode 100644 DESCRIPTION create mode 100644 NAMESPACE create mode 100644 R/methods.R create mode 100644 R/model_clonal.R create mode 100644 R/model_light.R create mode 100644 R/utility_functions.R create mode 100644 man/calcPPA.Rd create mode 100644 vignettes/example.R diff --git a/DESCRIPTION b/DESCRIPTION new file mode 100644 index 0000000..1713b00 --- /dev/null +++ b/DESCRIPTION @@ -0,0 +1,31 @@ +Package: PPA +Title: Modeling Plant Dynamics With The Perfect Plasticity Approximation +Version: 0.0.0.9000 +Authors@R: c( + person(given = "Jacob", + family = "Levine", + role = c("aut", "cre"), + email = "jacoblevine@princeton.edu"), + person(given = "Aiyu", + family = "Zheng", + role = "aut", + email = "aiyuz@princeton.edu"), + person(given = "Hannes", + family = "De Deurwaerder", + role = "aut", + email = "hannesd@princeton.edu"), + person(given = "Marco", + family = "Visser", + role = "aut", + email = "mvisser@princeton.edu")) +Description: Implementation of the Perfect Plasticity Approximation for light competition in plants. Models are predictive, and we include extensions + for clonal plants, annual plants and lianas. Water competition will soon be incorporated as well. +License: GPL-2 +Encoding: UTF-8 +LazyData: true +Imports: + ggplot2, + gridExtra, + grid, + gtable, + mgcv diff --git a/NAMESPACE b/NAMESPACE new file mode 100644 index 0000000..884a631 --- /dev/null +++ b/NAMESPACE @@ -0,0 +1,2 @@ +# Generated by roxygen2: fake comment so roxygen2 overwrites silently. +exportPattern("^[^\\.]") diff --git a/R/methods.R b/R/methods.R new file mode 100644 index 0000000..043316b --- /dev/null +++ b/R/methods.R @@ -0,0 +1,403 @@ + +## Methods file + +## this file contains functions which define methods for the S3 classes defined in this package + +## --------------------- IS. METHODS -------------------------- + +## methods which define class inheritance, each class requires one + +is.model_light <- function(x) inherits(x, "model_light") + +is.allometry <- function(x) inherits(x, "allometry") + +## --------------------- SUMMARY. METHODS -------------------------- + +## methods which print summary output for S3 class objects: + +#'SUMMARY.model_light +#' +#'summary method for model_light object +#' +#'@param simulation a simulation generate by model_light() +#' +#'@return a summary of relevant data from the input simulation + +summary.model_light <- function(simulation) { + + ## generate table of means and sds + + varnames <- c("diameter", "diameter_growth", "height", "crown_area", "struct_biomass")[c("diameter", "diameter_growth", "height", "crown_area", "struct_biomass") %in% names(simulation)] + table_text <- matrix(0, nrow = length(varnames), ncol = 2) + colnames(table_text) <- c("mean", "std. deviation") + rownames(table_text) <- varnames + for(i in 1:length(varnames)) { + if (varnames[i] == "diameter_growth") t <- simulation$max_t + else t <- simulation$max_t+1 + + table_text[i,1] <- mean(simulation[[varnames[i]]][, t], na.rm = TRUE) + table_text[i,2] <- sd(simulation[[varnames[i]]][, t], na.rm = TRUE) + } + + cat("Simulation:", "\n", + "Max time: ", as.character(simulation[["max_t"]]), "\n", + "Time Step: ", as.character(simulation[["timestep"]]), "\n", + "Iterations: ", as.character(simulation[["max_t"]]/simulation[["timestep"]]), "\n", + "Plot Area:", as.character(simulation[["plot_area"]]), "\n", + "\n", + "Means and standard deviations for variables at last time step:", "\n") + print(table_text) +} + + + +summary.model_clonal <- function(simulation) { + + + +} + + +## --------------------- PLOT. METHODS -------------------------- + +## methods which plot S3 class objects: + + +#'PLOT_ONERESULT +#' +#'This function plots the time trajectory of a single dimensional quantity from a simulation. +#' +#'@param data full.data object from simulation +#'@param y.val The value to plot on the y axis +#' +#'@return returns a plot of y.val vs time + +## not an S3 method. This gets wrapped into the plot.simulation method +plot_oneresult <- function(data, y.val) { + p <- ggplot(aes_string(x="timestep", y.val, color = "individual", linetype = "spp"), data = data) + + geom_line(size = 0.5) + + xlab("Time Step") + + ylab(y.val) + + scale_color_discrete(guide = FALSE) + + scale_linetype_discrete(name = "species") + + theme_bw() + + legend <- gtable_filter(ggplotGrob(p), "guide-box") + + p <- p + theme(legend.position = "none") + + return(list(p, legend)) +} + +#'plot.model_light +#' +#'plot method for model_light objects +#' +#'@param simulation an object of class "model_light" +#'@param y.values which dimensional qualities would you like to plot +#' +#'@return returns a grid of time trajectory plots +#' + +## plot method for simulation class +plot.model_light <- function(simulation, y.values) { + ind <- "full_data" + + if (ind == "full_data" & is.null(simulation[["full_data"]])) { + print("Error: full data must be defined in simulation") + return() + } + + data <- as.data.frame(simulation[[ind]]) + data[, "crown_class"] <- factor(data[, "crown_class"]) + data[, "spp"] <- factor(data[, "spp"]) + data[, "individual"] <- factor(data[, "individual"]) + + if (missing(y.values)) { + y.values <- as.list(colnames(data)[!(colnames(data) %in% c("timestep", "individual","crown_class", "spp", "diameter_growth"))]) + } + + plot.list <- lapply(y.values, plot_oneresult, data = data) + eval.list <- character() + + legend <- plot.list[[1]][[2]] + + for (i in 1:length(y.values)) { + if (i == length(y.values)) eval.list <- paste0(eval.list, "plot.list[[", i, "]][[1]]") + else eval.list <- paste0(eval.list, "plot.list[[", i, "]][[1]], ") + } + + eval(parse(text = paste0("grid.arrange(arrangeGrob(", eval.list, ", nrow = 2, top = textGrob('Stand Simulation Output')), legend, widths=unit.c(unit(1, 'npc') - legend$width, legend$width), nrow = 1)"))) + +} + + + +## plot methods for model_clonal + +plot.model_clonal <- function() { + + + +} + + +## --------------------- INDEXING. METHODS -------------------------- + +## define indexing method for allometry class. +"[.allometry" <- function(x, i, j) { + + output <- x$allometry[i, j] + return(output) + +} + + +## --------------------- ALLOMETRY. METHODS -------------------------- + +## allometry methods for the various simulation class types: + + +## eventually this function should have methods for each allometry type supported by the package: +## i.e. create_allometry.clonal, create_allometry.annual, create_allometry.forest (idk how specific we want to get? spp?) + +## this should eventually have different methods for the different model classes +## (i.e. create_allometry.forest, create_allometry.clonal, create_allometry.hydro) + +## some of the parameters in the resulting output are not, strictly speaking, allometric parameters. +## Perhaps this should be called "defaults", as Marco suggested, or those non-allometric parameters could +## be specified in a seperate defaults method. + + +#'CREATE_ALLOMTRY.MODEL_LIGHT +#' +#'creates an allometry object to be used in the model_light simulator +#' +#'@param n.spp the number of spp to create allometry for +#'@param custom a named list supplying custom values for any or all allometric parameters +#'@param stochastic logical, if true inject random variation into the parameters +#' +#'@return returns an object of class "allometry" to be used in model_light simulator + + +create_allometry.model_light <- function(n.spp, custom, stochastic = TRUE) { + ## allometry will be a dataframe of allometric variables, with rows for vars and columns for spp + ## this function also assigns values for non-allometric, but relevant variables. Not sure if should be in separate file + + ## for now, almost all numbers will be same for all species, some will be drawn from random dist. + out <- as.data.frame(matrix(1, nrow = 25, ncol = n.spp)) + colnames(out) <- c(as.character(seq(1, n.spp, by = 1))) + rownames(out) <- c("l_c", "l_u", "r_c", "r_u", "mu_c", "mu_u", "r", "H", "alpha_s", "alpha_w", "gamma", "c_lb", + "c_rb", "c_bg", "tau_l", "tau_r", "p_r", "p_l", "p_sw", "c_l", "c_r", "V", "k", "L_0", "a_f") + + ## function accepts a list of named vectors for each species with custom allometric constants + if(!missing(custom)) { + list_classes <- lapply(custom, class) + list_names <- lapply(custom, names) + if(class(custom) != "list") { + print("Error: custom variables must be a list of named numeric vectors") + return() + } + else if(any(list_classes != "numeric")) { + print("Error: custom variables must be a list of named numeric vectors") + return() + } + else if(any(is.null(list_names))) { + print("Error: custom variables must be a list of named numeric vectors") + return() + } + } + + for (spp in as.character(1:n.spp)) { + + #### allometric variables + out["l_c", spp] <- 4 ## leaf area index of tree canopy tree + out["l_u", spp] <- 1 ## leaf are index of understory tree + out["r_c", spp] <- 3 + out["r_u", spp] <- 2 + out["mu_c", spp] <- 0.005 + out["mu_u", spp] <- 0.015 + out["r", spp] <- 6.5 ## fine-root surface area per unit crown area + out["H", spp] <- 3.6 ## allometric constant for height + out["alpha_s", spp] <- 0.0815 ## allometric constant for sapwood + out["alpha_w", spp] <- 0.20 ## allometric constant for crown area + out["gamma", spp] <- 1.5 ## allometric exponent + + #### carbon accumulation and allocation variables + out["c_lb", spp] <- 0.07656 ## cost of building leaf biomass per LAI + out["c_rb", spp] <- 0.02933 ## cost of builing root biomass per LAI + out["c_bg", spp] <- 0.2 ## building cost of structural biomass + out["tau_l", spp] <- 1 ## average lifetime of a unit carbon in leaf + out["tau_r", spp] <- 2 ## average lifetime of a unit carbon in roots + out["p_r", spp] <- 0.02933 ## respiration rate of roots per unit surface area + out["p_l", spp] <- 0.0638 ## respiration rate of leaves per unit surface area + out["p_sw", spp] <- 0.0466 ## respiration rate of sapwood + out["c_l", spp] <- out["c_lb", spp] / out["tau_l", spp] + out["p_l", spp] + out["p_r", spp] ## cost of building and maintaining leaf biomass per unit LAI + out["c_r", spp] <- out["c_rb", spp] / out["tau_r", spp] + out["p_r", spp] ## cost of building and maintaining root biomass per unit RAI + out["V", spp] <- 0.6 ## maximum rate of carbon fixation - kgC/m^2/day + out["k", spp] <- 0.33 ## light extinction coefficient + out["L_0", spp] <- 1200 ## light level above highest canopy + out["a_f", spp] <- 0.001 ## conversion rate from photons to carbohydrates + } + + ## reassign custom variables. Perhaps there is a way to do his that doesnt require overwriting? + if(!missing(custom)) { + for (i in 1:length(custom)) { + for (j in 1:length(custom[[i]])) { + out[var = names(custom[[i]])[j], i] <- custom[[i]][j] + } + } + } + + + ## crude implementation of random variation for now + ## needs major refinement + if (stochastic) { + for (spp in as.character(1:n.spp)) { + out[, spp] <- rnorm(nrow(out), mean = out[, spp], sd = 0.05*out[, spp]) + + } + } + outlist <- list(allometry = out) + class(outlist) <- "allometry" + + return(outlist) +} + +#'CREATE_ALLOMTRY.MODEL_CLONAL +#' +#'creates an allometry object to be used in the model_clonal simulator +#' +#'@param n.spp the number of spp to create allometry for +#'@param custom a named light supplying custom values for any or all allometric parameters +#'@param stochastic logical, if true inject random variation into the parameters +#' +#'@return returns an object of class "allometry" to be used in model_clonal simulator + +##eventually, I want this create_allometry.model_clonal to be able to generate clonal allometries based on the general create_allometry.model_light +##function. But this requires altering some parameters in the current structure of model_light. Ideally, this function should generate n species' clonal versions parallel +##to their non-clonal ones, both for annuals and biennals. But we don't have a allometry for trees (as there is no annual reproduction cost now) at +## this point. And also since time step is different for trees, all the rates differ. I just assign values and a seperate set of parameters for now. + +create_allometry.model_clonal <- function(n.spp, custom, stochastic = TRUE) { + ## allometry will be a dataframe of allometric variables, with rows for vars and columns for spp + ## for now, all the species-specific parameter values are drawn from the create_allometry.model_light based on indicated species + allometry<-create_allometry.model_light(n.spp = n.spp) + out <- as.data.frame(matrix(1, nrow = 31, ncol = n.spp)) + colnames(out) <- c(as.character(seq(1, n.spp, by = 1))) + rownames(out) <- c("a_f","L_0","V","k","v","g","cf","f","F","l_c","l_u","delta_rl","sigma","w","omega", + "epsilon_l","epsilon_r","r_l","r_r","lambda","rho_w","phi_z","phi_l","phi_s","theta_s","theta_l","theta_r","theta_z", + "mu_c","mu_u","mu_cl") + + ## function accepts a list of named vectors for each species with custom allometric constants + if(!missing(custom)) { + list_classes <- lapply(custom, class) + list_names <- lapply(custom, names) + if(class(custom) != "list") { + print("Error: custom variables must be a list of named numeric vectors") + return() + } + else if(any(list_classes != "numeric")) { + print("Error: custom variables must be a list of named numeric vectors") + return() + } + else if(any(is.null(list_names))) { + print("Error: custom variables must be a list of named numeric vectors") + return() + } + } + + for (spp in as.character(1:n.spp)) { + #### light calculations + out["a_f",spp]<-0.01 ## relationship between carbon fixation and light intensity + out["L_0",spp]<-2600 ## PAR at the canopy top, MJ PAR/m^2/s + out["V",spp]<-out["L_0",spp]*out["a_f",spp] ## per unit leaf area productivity in full sun,kgC/m^2/yr + out["k",spp]<-0.33 ## beer-lambert extinction coefficient + + #### carbon accumulation and allocation variables + out["v", spp] <- 0.5 ## allocation to vegetative propogation (clonal growth) per unit crown area, kgC/m^2 + out["g", spp] <- 1 ## germination probability for seeds + out["cf", spp] <- 20 ## cost of investing in one new seedling, kgC/m^2 + out["f", spp] <-1.4 ## cost of fecundity per unit crown area,kgC/m^2 + out["F", spp] <-out["f", spp]/out["cf", spp] ## new individuals per year per unit m2 + out["l_c", spp] <- allometry["l_c", spp] ## leaf area index of a canopy tree, m^2/m^2 + out["l_u", spp] <- allometry["l_u",spp] ## leaf area index of an understory tree, m^2/m^2 + out["delta_rl",spp]<-0.6 ## ratio of total root surface area to total leaf area + out["sigma", spp] <- 0.0245 ## LMA, kgC/m^2 + out["w",spp] <-43900 ## specific root length, m/kgC + out["omega",spp]<-35 ## specific root area, m^2/kgC + out["epsilon_l",spp]<-0.3 ## maintenance rate for the tree crown, kgC/m^2 + out["epsilon_r",spp]<-0.24 ## maintenance rate for the roots, kgC/m^2 + out["r_l",spp]<-0.06 ## respiration rate for the tree crown, kgC/m^2 + out["r_r",spp]<-0.04 ## respiration rate for the roots, kgC/m^2 + out["lambda",spp]<-0.75 ## scaling parameter of cylinder volume to woody biomass + out["rho_w",spp]<-230 ## wood density, kg*C/m^3 + + ##allometric scaling parameters + out["phi_z",spp]<-31 ## scaling parameter of height with diameter + out["phi_l",spp]<-130 ## scaling parameter of crown area with diameter + out["phi_s",spp]=0.25*pi*out["lambda",spp]*out["rho_w",spp]*out["phi_z",spp] ## scaling parameter of bole wood with dbh + out["theta_s",spp]<-allometry["gamma",spp]+1 ## allometric exponent for scaling D to bolewood mass + out["theta_l",spp]<-allometry["gamma",spp] ## allometric exponent for scaling D to crown area + out["theta_r",spp]<-allometry["gamma",spp] ## allometric exponent for scaling D to root area + out["theta_z",spp]<-allometry["gamma",spp]-1 ## allometric exponent for scaling D to height + + ##death rates (need to be somewhere else) + out["mu_c", spp] <- 0.01 ## (exponential) death rate in the canopy + out["mu_u", spp] <- 0.12 ## (exponential) death rate in the understory + out["mu_cl", spp] <- 0.1 ## (exponential) death rate for any ramet (root sucker) of a clonal plant + } + + ## reassign custom variables. Perhaps there is a way to do his that doesnt require overwriting? + if(!missing(custom)) { + for (i in 1:length(custom)) { + for (j in 1:length(custom[[i]])) { + out[var = names(custom[[i]])[j], i] <- custom[[i]][j] + } + } + } + ## crude implementation of random variation for now + ## needs major refinement + if (stochastic) { + for (spp in as.character(1:n.spp)) { + out[, spp] <- rnorm(nrow(out), mean = out[, spp], sd = 0.05*out[, spp]) + + } + } + outlist <- list(allometry = out) + class(outlist) <- "allometry" + return(outlist) +} + + +#'CALCULATE_ALLOMETRY.model_light +#' +#'function to calculate full allometry from object of class simulation +#' +#'@param simulation a simulation generated by model_light +#'@param allometry an object of class "allometry" +#'@param spp a named list of species +#' +#'@return returns datasets for height, crown_area and biomass + +calculate_allometry.model_light <- function(simulation, allometry, spp) { + + if(!is.model_light(simulation)) { + print("Error: input must be of class simulation") + return() + } + + simulation[["height"]] <- matrix(1, nrow = nrow(simulation[["diameter"]]), ncol = ncol(simulation[["diameter"]])) + simulation[["crown_area"]] <- matrix(1, nrow = nrow(simulation[["diameter"]]), ncol = ncol(simulation[["diameter"]])) + simulation[["struct_biomass"]] <- matrix(1, nrow = nrow(simulation[["diameter"]]), ncol = ncol(simulation[["diameter"]])) + + for (s in spp) { + + simulation[["height"]][simulation[["spp"]] == s] <- allometry["H", s] * (simulation[["diameter"]][simulation[["spp"]] == s] ^ (allometry["gamma", s] - 1)) + simulation[["crown_area"]][simulation[["spp"]] == s] <- allometry["alpha_w", s] * (simulation[["diameter"]][simulation[["spp"]] == s] ^ allometry["gamma", s]) + simulation[["struct_biomass"]][simulation[["spp"]] == s] <- allometry["alpha_s", s] * (simulation[["diameter"]][simulation[["spp"]] == s] ^ (allometry["gamma", s] + 1)) + + } + + return(simulation) +} diff --git a/R/model_clonal.R b/R/model_clonal.R new file mode 100644 index 0000000..01d125b --- /dev/null +++ b/R/model_clonal.R @@ -0,0 +1,121 @@ +##devtools::load_all("/Users/aiyuzheng/PPA/") +## This file contain code for the function model_clonal. + +##model_clonal is not an explicit simulator at present. It assumes infinite space. It's for simulating the LRS for a clonal tree when the mother/sucker are killed randomly. +##but it outputs growth rates, allometries, and lifetime reproductive sucess. I hope we can eventually incorporate clonal allometries into model_light +##or model_water, but I am still working on ranking clonal crowns. So how far we want to go will depend on collective decision. I will put in the codes for invasion analysis and visualization +##after my general exam :) + +#' This function simulates a clonal tree invading a population of non-clonal trees at equilibrium or vice versa using the perfect plasticity approximation for light competition. +#' +#'@param allometry an object of class: clonal or non-clonal allometries created using the create_allometry method +#'@param max_t the last timestep in the simulation +#'@param timestep the length of each timestep +#'@param N the initial population size of the clonal trees' understory seedlings +#'@param dnot the initial diameter for each tree seedling +#'@param dnot_ramet the initial size of a ramet's diameter +#'@param nstem the number of stems a clonal plant develops +#'@param full.lifehistory TRUE/FALSE specifying whether you would like the model output to include life histories of a clonal tree and a non-clonal tree +#'@param recip TRUE/FALSE specifying the direction of invasion. TRUE indicates the clonal tree being the resident, and FALSE indicates the clonal tree being the invader. +#' +#'@return a table containing life time reproductive success for each simulated tree or the final average LRS + + +## simulator for light competition model. +model_clonal <- function(allometry = create_allometry.model_clonal(n.spp = 1,stochastic = FALSE), + spp, + max_t, + timestep, + N, + dnot, + dnot_ramet, + nstem=2, + full.lifehistory = TRUE, + recip = FALSE) { + ## stop on error + stopifnot(is.numeric(max_t), + is.numeric(timestep), + is.numeric(dnot), + is.numeric(dnot_ramet), + is.numeric(nstem), + timestep <= max_t, + is.logical(full.lifehistory), + is.logical(recip)) + + ## ensure spp name is character not numeric + spp <- as.character(spp) + + ##generate life history data for a non-clonal tree and a clonal tree + + NonclonalLH<-get.growth_nonclonal(dnot=dnot,allometry=allometry,spp=spp,max_t=max_t,timestep=timestep) + + ClonalLH<-get.growth_clonal(dnot=dnot,dnot_ramet=dnot_ramet,allometry = allometry,spp=spp,max_t=max_t,timestep=timestep) + + ##get tstar + tstar<-equil_nonclonal(dnot=dnot,allometry=allometry,spp = spp,max_t=max_t,timestep = timestep)[["tstar"]] + ##get tbreak + tbreak<-tbreak_clonal(ClonalLH=ClonalLH) + + ##N trees to start with + + ####kill some trees in the understory + N_understory<-rexp(N, rate = allometry["mu_u",spp]) + N_canopy<-round(N_understory[N_understory>(tstar*timestep)],digits=1) + + ####assign years of death to a clonal tree's primary stem and secondary stem in the canopy + Death<-matrix(ncol = 4,nrow = length(N_canopy)) + colnames(Death)<-c("YMD","YBD","theme","LRS") + rownames(Death)<-as.character(seq(1,length(N_canopy)),by=1) + Death[,"YMD"]<-tstar+round(rexp(length(N_canopy), rate = allometry["mu_c",spp]),digits = 1)/timestep + Death[,"YBD"]<-tstar+round(rexp(length(N_canopy), rate = allometry["mu_cl",spp]),digits=1)/timestep + # Death[Death[,"YBD"]>tbreak*timestep]<-tbreak+round(rexp(length(Death[Death[,"YBD"]>tbreak*timestep]), rate = allometry["mu_c",spp]),digits = 1)/timestep + Death[Death > max_t] <- max_t ##make sure trees don't live longer than max_t + Death[,"theme"]<-NA + Death[,"LRS"]<-NA + + ##based on data frame death, sum up the fecundity for each tree + for (i in 1:nrow(Death)){ + if(Death[i,"YMD"]<=tbreak){ + if(Death[i,"YBD"]<=tbreak){ + if(Death[i,"YMD"]>Death[i,"YBD"]){ + Death[i,"theme"]<-1 + Death[i,"LRS"]<-Nursing.prim(ClonalLH = ClonalLH,tbreak=Death[i,"YBD"])+Switching.prim(NonclonalLH = NonclonalLH, ClonalLH=ClonalLH, tbreak = Death[i,"YBD"], YMD=Death[i,"YMD"]) + }else{ + Death[i,"theme"]<-2 + Death[i,"LRS"]<-Nursing.prim(ClonalLH = ClonalLH,tbreak=Death[i,"YMD"]) + } + }else{ + if(tbreak_switching.ramet(NonclonalLH = NonclonalLH,ClonalLH = ClonalLH, YMD=Death[i,"YMD"],tstar=tstar)>=Death[i,"YBD"]){ + Death[i,"theme"]<-2 + Death[i,"LRS"]<-Nursing.prim(ClonalLH = ClonalLH,tbreak=Death[i,"YMD"]) + }else{ + Death[i,"theme"]<-4 + Death[i,"LRS"]<-Nursing.prim(ClonalLH = ClonalLH,tbreak=Death[i,"YMD"])+Switching.ramet(NonclonalLH = NonclonalLH,ClonalLH = ClonalLH,YMD = Death[i,"YMD"],YBD=Death[i,"YBD"],tstar=tstar) + } + } + }else{ + if(Death[i,"YBD"]<=tbreak){ + Death[i,"theme"]<-3 + Death[i,"LRS"]<-Nursing.prim(ClonalLH = ClonalLH,tbreak=Death[i,"YBD"])+Switching.prim(NonclonalLH = NonclonalLH, ClonalLH= ClonalLH,tbreak = Death[i,"YBD"], YMD=Death[i,"YMD"]) + }else{ + Death[i,"theme"]<-5 + Death[i,"LRS"]<-Nursing.prim(ClonalLH = ClonalLH,tbreak=Death[i,"YMD"])+Success.ramet(ClonalLH = ClonalLH,YBD=Death[i,"YBD"]) + } + } + } + + ## if user specifies that full data should be reported in single matrix, initialize and populate matrix + if(full.lifehistory){ + outlist<-list("clonal_lifehistory"=ClonalLH,"non-clonal_lifehistory"=NonclonalLH,"simulated_lifehistory"=Death) + class(outlist) <- "model_clonal" + return(outlist) + }else{ + output <- Death + class(output) <- "model_clonal" + return(output) + } +} + + + + diff --git a/R/model_light.R b/R/model_light.R new file mode 100644 index 0000000..87bfdde --- /dev/null +++ b/R/model_light.R @@ -0,0 +1,209 @@ + + +#' Model Light Competition +#' +#' This function simulates a stand of trees using the perfect plasticity approximation for light competition. +#' +#'@param allometry an object of class: allometry created using the create_allometry method +#'@param spp a vector of species names to simulate +#'@param max_t the last timestep in the simulation +#'@param timestep the length of each timestep +#'@param dnot the initial diameter for each tree +#'@param plot.area the total area of the plot to be simulated +#'@param n_start a named list of the starting abundances. Names must correspond to the spp names specified in spp. +#'@param full.allom TRUE/FALSE specifying whether you would like the model output to include height and structural biomass +#'@param full.data TRUE/FALSE specifying whether you would like the model output to include data from each timestep or just the last timestep +#' +#'@return a list object of class "model_light" which includes simulated growth and allometric data + + +## simulator for light competition model. +model_light <- function(allometry = create_allometry.model_light(n.spp = 1), + spp = 1, + max_t, + timestep, + dnot = rnorm(n = n_start, mean = 0.0004, sd = 0.00005), + plot.area, + n_start, + full.allom = TRUE, + full.data = TRUE) { + + ## stop on error + stopifnot(ncol(allometry) >= length(spp), + is.numeric(max_t), + is.numeric(timestep), + is.numeric(dnot), + is.numeric(plot.area), + timestep <= max_t, + is.logical(full.allom), + !is.null(names(n_start))) + + ## ensure spp name is character not numeric + spp <- as.character(spp) + + ## initialize data list + data <- list(spp = matrix(1, nrow = sum(unlist(n_start)), ncol = max_t/timestep + 1), + diameter = matrix(1, nrow = sum(unlist(n_start)), ncol = max_t/timestep+1), + crown_class = matrix(1, nrow = sum(unlist(n_start)), ncol = max_t/timestep + 1), + diameter_growth = matrix(1, nrow = sum(unlist(n_start)), ncol = max_t/timestep + 1)) + + data[["diameter"]][,1] <- dnot + data[["spp"]][,1] <- c(rep(spp[1], times = n_start[spp[1]]), rep(spp[2], times = n_start[spp[2]])) + + A <- calc.A(allometry = allometry, n.crown.class = 2, spp = spp) + + ## loop over days in simulation + for (i in seq(1, max_t, by = timestep)) { + + ## ------------------------- KILL PLANTS ----------------------------- ## + + ## pre-calculate # of each spp in each crown class + spp_pops <- list("1" = sapply(FUN = spp_populations, spp, i = i, crown_class = 1, data = data, simplify = "vector"), + "2" = sapply(FUN = spp_populations, spp, i = i, crown_class = 1, data = data, simplify = "vector")) + + ## create lists of death outcomes based on each spp's individual death probability for canopy and understor + mu_c <- mapply(rbinom, spp_pops[["1"]], matrix(allometry["mu_c", 2:length(spp)], nrow = 1), size = 1) + + if(length(spp_pops[["2"]]) > 0) { + mu_u <- mapply(rbinom, spp_pops[["2"]], matrix(allometry["mu_u", 2:length(spp)], nrow = 1), size = 1) + mu <- list(mu_c) + } + else { + mu <- list(mu_c, mu_u) ## combine into one list + } + ## reorder data to easily index and align death vector and data matrices + + data <- order_simulation(data = data, by = c("crown_class", "spp"), decreasing = TRUE, i = i) + + ## kill trees by removing all trees with mu == 0 from data. In future could include tracking of dead trees, or perhaps mortality rates? + ## Though not particularly interesting with "density independent" death process + data[["spp"]] <- data[["spp"]][unlist(mu) == 0, ] + data[["diameter"]] <- data[["diameter"]][unlist(mu) == 0, ] + data[["crown_class"]] <- data[["crown_class"]][unlist(mu) == 0, ] + data[["diameter_growth"]] <- data[["diameter_growth"]][unlist(mu) == 0, ] + + ## ------------------------- ASSIGN NEW SPP ID ------------------------## + + data[["spp"]][, i + 1] <- data[["spp"]][, i] + + ## ------------------------- GROW PLANTS ------------------------------ ## + + ## grow plants + ## calculate diameter_growth for overstory and understory + ## currently goint to put this into a loop. There should be a more efficient way to do this using apply and/or lookup table, but my head hurts + for (j in 1:length(spp)) { + + if (sum(data[["crown_class"]][, i] == 1 & data[["spp"]][, i] == spp[j]) > 0) { + + data[["diameter_growth"]][data[["crown_class"]][, i] == 1 & data[["spp"]][, i] == spp[j], i] <- sapply(data[["diameter"]][data[["crown_class"]][, i] == 1 & data[["spp"]][, i] == spp[j], i], + calc.dd, allometry = allometry, + l = allometry["l_c", spp[j]], + r <- allometry["r_c", spp[j]], + spp = spp[j], A = A[spp[j], "A_c"], + simplify = "vector") + + } + + if (sum(data[["crown_class"]][, i] == 2 & data[["spp"]][, i] == spp[j]) > 0) { + + if (class(data[["diameter_growth"]]) != "matrix") print(c(i, j)) + + data[["diameter_growth"]][data[["crown_class"]][, i] == 2 & data[["spp"]][, i] == spp[j], i] <- sapply(data[["diameter"]][data[["crown_class"]][, i] == 2 & data[["spp"]][, i] == spp[j], i], + calc.dd, allometry = allometry, + l = allometry["l_u", spp[j]], + r <- allometry["r_u", spp[j]], + spp = spp[j], A = A[spp[j], "A_u"], + simplify = "vector") + + } + + } + + ## this is just to improve interpretability for summary functions, likely there is a cleaner way to do this? + if (i == max_t) { + + data[["diameter_growth"]][, i+1] <- NA + + } + + ## calculate diameter for next timestep + data[["diameter"]][, i+1] <- data[["diameter"]][, i] + data[["diameter_growth"]][, i] + + ## ------------------- CANOPY CLASS REASSIGNMENT ----------------- ## + + CA <- sum(mapply(FUN = calc.CA, + alpha_w = allometry["alpha_w", spp], + gamma = allometry["gamma", spp], + spp = spp, + MoreArgs = list(diam_data = data[["diameter"]], + spp_data = data[["spp"]], + i = i), + SIMPLIFY = "vector")) + + ## check if plants need to added to the understory + if (CA <= plot.area) { + + data[["crown_class"]][, i+1] <- 1 + + } + + else { + ## order data from largest to smallest + data <- order_simulation(data = data, by = c("diameter"), decreasing = TRUE, i = i) + + ## create vector of cumulative CA sums, used to determine which trees to add/remove from canopy + ca.sum <- vector(length = nrow(data[["diameter"]]), mode = "numeric") + + for (j in 1:nrow(data[["diameter"]])) { + + if (j == 1) { + + ca.sum[j] <- allometry["alpha_w", data[["spp"]][j, i+1]] * data[["diameter"]][j, i+1]^allometry["gamma", data[["spp"]][j, i+1]] + + } + + else { + + ca.sum[j] <- sum((allometry["alpha_w", data[["spp"]][j, i+1]] * data[["diameter"]][j, i+1]^allometry["gamma", data[["spp"]][j, i+1]]), ca.sum[j - 1]) + + } + + ## if gap in canopy trees, add understory trees to canopy and vis versa + data[["crown_class"]][ca.sum <= plot.area, i+1] <- 1 + data[["crown_class"]][ca.sum > plot.area, i+1] <- 2 + } + } + + } + + ## assign informative values to list items + data[["max_t"]] <- max_t + data[["timestep"]] <- timestep + data[["plot_area"]] <- plot.area + + ## assign class of return object to class simulation + class(data) <- "model_light" + + ## if user specifies full allometry, calculate and add to the data + if (full.allom) { + + data <- calculate_allometry.model_light(data, allometry = allometry, spp = spp) + + } + + ## if user specifies that full data should be reporeted in single matrix, initialize and populate matrix + if (full.data) { + + colnames <- names(data)[!(names(data) %in% c("avgs", "max_t", "timestep", "plot_area"))] + columns <- lapply(colnames, extract_column, data = data) + full_data <- matrix(unlist(columns), ncol = length(colnames)) + time_column <- rep(seq(1, (max_t+timestep), by = timestep), each = nrow(data[["diameter"]])) + ind_column <- rep(seq(1, nrow(data[["diameter"]]), by = 1), times = (max_t/timestep)+1) + full_data <- cbind(time_column, ind_column, full_data) + colnames(full_data) <- c("timestep", "individual", colnames) + data[["full_data"]] <- full_data + + } + + return(data) +} diff --git a/R/utility_functions.R b/R/utility_functions.R new file mode 100644 index 0000000..c78409d --- /dev/null +++ b/R/utility_functions.R @@ -0,0 +1,553 @@ +## Utility Functions +## anything that doesn't really belong somewhere else, or is generally useful + +#' EXTRACT_COLUMN +#' +#' function to extract a single column from a simulation data matrix +#' +#'@param data the simulation dataset that contains the column you want to extract +#'@param colname the name of the column to extract from data +#' +#'@return returns the column specified + +extract_column <- function(colname, data) { + + data_matrix <- data[[colname]] + + if(!is.matrix(data_matrix)) { + print("Error: column matrix must be of class matrix") + return() + } + + column <- as.numeric(data_matrix) + return(column) + +} + +#' SPP_POPULATIONS +#' +#' function to calculate the number of species of a given crown_class at a given timepoint +#' +#'@param data the simulation dataset +#'@param spp the name of the species whose population you want to extract +#'@param i the timestep at which you want to extract the populatioon +#'@param crown_class which crown class's population do you wish to extract +#' +#'@return returns the population for the species, timestep and crown class you specified +spp_populations <- function(data, spp, i, crown_class) { + if(missing(crown_class)) { + out <- sum(data[["spp"]][,i] == spp) + } + else { + out <- sum(data[["crown_class"]][, i] == crown_class & data[["spp"]][, i] == spp) + } + return(out) +} + +#'CALC.CA +#' +#'A function to calculate the crown area of a plant given its allometric constants and diameter. Works specifically within +#'the context of stand simulators +#' +#'@param alpha_w allometric constant governing the relationship between diameter and crown area +#'@param gamma allometric exponent governing the relationship between diameter and crown area +#'@param spp the name of the species for which you are calculating CA +#'@param diam_data diameter data for the species, specifically from within a model_ call +#'@param spp_data spp data from model_ call +#'@param i current timestep +#' +#'@return returns a numeric value, the crown area for the specied plant or plants + +## need to separate diameter and spp data in order to work with mapply funciton - this is annoying, would rather be able to past data as arg, but mapply tries to interpret it as multiple args +calc.CA <- function(alpha_w, gamma, spp, diam_data, spp_data, i) { + + out <- sum(alpha_w * (diam_data[spp_data[, i + 1] == spp, i + 1]^gamma)) + return(out) + +} + + +#'CALC.A +#' +#' this function calculates max photosynthetic rates given allometric constants, crown classes and the spp info +#' +#'@param allometry an object of class "allometry" +#'@param n.crown.class The number of crown classes to calculate A for, currently only 2 supported +#'@param spp a vector of species names, similar to the input for model_ calls +#' +#'@return returns the maximum photosynthetic rate + +calc.A <- function(allometry, n.crown.class = 2, spp) { + + spp <- as.character(spp) + out <- data.frame(matrix(1, nrow = length(spp), ncol = n.crown.class+1)) + colnames(out) <- c("spp","A_c", "A_u")[1:(n.crown.class+1)] + out[, "spp"] <- spp + + ## calculate maximum light level for understory trees + L_u <- allometry["L_0", spp] * (exp(-(allometry["k", spp]) * allometry["l_c", spp])) + + out[, "A_c"] <- as.numeric((allometry["V", spp] / allometry["k", spp]) * (1 + log((allometry["a_f", spp] * allometry["L_0", spp])/allometry["V", spp]) - ((allometry["a_f", spp] * allometry["L_0", spp]) / allometry["V", spp]) * exp(-allometry["k", spp] * allometry["l_c", spp]))) + + if(n.crown.class == 2) { + out[, "A_u"] <- as.numeric((allometry["a_f", spp] * L_u) / allometry["k", spp] * (1 - exp(-allometry["k", spp] * allometry["l_u", spp]))) + } + if(n.crown.class > 2) { + print("error: only 2 crown classes currently supported.") + return() + } + return(out) +} + +#'CALC.A_CLONAL +#' +#' this function calculates max photosynthetic rates given allometric constants, crown classes and the spp info for model_clonal +#' +#'@param allometry an object of class "allometry" +#'@param n.crown.class The number of crown classes to calculate A for, currently only 2 supported +#'@param spp a vector of species names, similar to the input for model_ calls +#' +#'@return returns the maximum photosynthetic rate + +calc.A_clonal <- function(allometry, n.crown.class = 2, spp) { + + spp <- as.character(spp) + out <- data.frame(matrix(1, nrow = length(spp), ncol = n.crown.class+1)) + colnames(out) <- c("spp","A_c", "A_u")[1:(n.crown.class+1)] + out[, "spp"] <- spp + + ## calculate maximum light level for understory trees + L_u <- allometry["L_0", spp] * (exp(-(allometry["k", spp]) * allometry["l_c", spp])) + + out[, "A_c"] <- as.numeric((allometry["V", spp] / allometry["k", spp]) * (1 - exp(-allometry["k", spp] * allometry["l_c", spp]))) + + if(n.crown.class == 2) { + out[, "A_u"] <- as.numeric((allometry["a_f", spp] * L_u) / allometry["k", spp] * (1 - exp(-allometry["k", spp] * allometry["l_u", spp]))) + } + if(n.crown.class > 2) { + print("error: only 2 crown classes currently supported.") + return() + } + return(out) +} + +#'CALC.DD +#'function to calculate diameter growth +#' +#'@param diameter the current diameter of the plant +#'@param allometry an object of class "allometry" +#'@param spp the spp name +#'@param l the leaf area of the plant +#'@param r the root area of the plant' +#' +#'@return returns a diameter increment + +calc.dd <- function(diameter, allometry, spp, l, r, A) { + dd <- 1/((allometry["alpha_s", spp] * (allometry["gamma", spp]+1) * (1+allometry["c_bg", spp]) * (1/allometry["alpha_w", spp])) + ((allometry["gamma", spp]/diameter) * (l * allometry["c_lb", spp] + r * allometry["c_rb", spp]))) * (A - (l*allometry["c_l", spp]) - (r*allometry["c_r", spp])) + return(dd) +} + +#'CALC.DD_NONCLONAL +#'function to calculate diameter growth rate, i.e., dD/dt for a non-clonal plant in model_clonal +#' +#'@param diameter the current diameter of the plant +#'@param allometry an object of class "allometry" +#'@param spp the spp name +#'@param understory logical, if true calculate the diameter growth rate for an understory plant; if false calculate the diameter growth rate for a canopy tree. +#' +#'@return returns a diameter growth rate at a particular diameter size + +calc.dd_nonclonal <- function(diameter, allometry, spp, understory=TRUE) { + spp <- as.character(spp) + CarbonGain<-calc.A_clonal(allometry = allometry,spp=spp) + A<-CarbonGain[spp==spp,"A_c"] + dDdt<-(allometry["phi_l",spp]*diameter^allometry["theta_l",spp]* + (A-(allometry["epsilon_l",spp]+allometry["r_l",spp])* + allometry["l_c",spp]-(allometry["epsilon_r",spp]+allometry["r_r",spp])* + allometry["delta_rl",spp]*allometry["l_c",spp]/allometry["omega",spp]-allometry["f", spp]))/ + (allometry["theta_s",spp]*allometry["phi_s",spp]*diameter^(allometry["theta_s",spp]-1)+(allometry["sigma",spp]+ + allometry["delta_rl",spp]*allometry["l_c",spp]*allometry["omega",spp])*allometry["phi_l",spp]*allometry["theta_l",spp]*diameter^(-1)) + + if(understory){ + A<-CarbonGain[spp==spp,"A_u"] + dDdt=(allometry["phi_l",spp]*diameter^allometry["theta_l",spp]* + (A-(allometry["epsilon_l",spp]+allometry["r_l",spp])* + allometry["l_u",spp]-(allometry["epsilon_r",spp]+allometry["r_r",spp])* + allometry["delta_rl",spp]*allometry["l_u",spp]/allometry["omega",spp]))/ + (allometry["theta_s",spp]*allometry["phi_s",spp]*diameter^(allometry["theta_s",spp]-1)+(allometry["sigma",spp]+ + allometry["delta_rl",spp]*allometry["l_u",spp]*allometry["omega",spp])*allometry["phi_l",spp]*allometry["theta_l",spp]*diameter^(-1)) + } + ifelse(dDdt<0,print("Error: diameter growth cannot be negative"),return(dDdt)) +} + +####This function below really needs to be faster. Currently I am using a loop because I don't know a better way... + +#'EQUIL_NONCLONAL +#'function to pre-caculate the equilibrium diameter size for a non-clonal plant in the understory over time in model_clonal. +#' +#'@param dnot the initial diameter size of a seedling to begin with, in meters +#'@param allometry an object of class "allometry" +#'@param spp the spp name +#'@param max_t the total number of time steps of the simulation +#'@param timestep the length of each timestep, in years +#' +#'@return returns the equilibrium values of diameter and the number of timesteps to reach the canopy at population equilibrium + +equil_nonclonal<-function(dnot,allometry,spp,max_t,timestep){ + out <- as.data.frame(matrix(1, nrow =max_t, ncol = 2)) + colnames(out) <- c("diameter","LRS_test") + rownames(out) <- c(as.character(seq(1,max_t, by = 1))) + + ##intergrand is the equation for the number of seeds a tree reproduces during its time in the canopy + out[1,"diameter"]<-dnot + for (i in 2:max_t){ + out[i,"diameter"]<-out[i-1,"diameter"]+calc.dd_nonclonal(diameter=out[i-1,"diameter"],allometry=allometry,spp = spp)*timestep + intergrand<-function(t){ + value<-exp(-(allometry["mu_c",spp])*t)*out[i,"diameter"]^allometry["theta_l",spp]*allometry["F",spp]*allometry["g",spp]*allometry["phi_l",spp] + return(value) + } + out[i,"LRS_test"]<-exp(-(allometry["mu_u",spp])*timestep*i)*integrate(intergrand,lower = 0,upper=Inf)$value + if(out[i,"LRS_test"]>1){ + }else{ + break + } + } + ##dstar is the diameter size when a tree reaches canopy in its population at equilibrium, and tstar is the corresponding time it needs + tstar<-i-1 + dstar<-out[tstar,"diameter"] + understory_diameter<-as.matrix(out[1:tstar,1]) + colnames(understory_diameter)<-"diameter" + outlist<-list("understory_diameter"=understory_diameter,"tstar"=tstar,"dstar"=dstar) + return(outlist) +} + + +#'GET.GROWTH_NONCLONAL +#'function to calculate the diameter growth rates dD/dt in the understory and in the canopy based on equilibrium diameter size for a non-clonal plant in model_clonal +#' +#'@param dnot the initial diameter size of a seedling to begin with, in meters +#'@param allometry an object of class "allometry" +#'@param spp the spp name +#'@param max_t the total number of time steps of the simulation +#'@param timestep the length of each timestep, in years +#'@param full.allom logical, if true returns allometric scaling results along with diameter sizes +#' +#'@return returns a data frame containing diameter, bolewood mass, crown area, height based on defined allometry. + +get.growth_nonclonal<-function(dnot,allometry,spp,max_t,timestep,full.allom=TRUE){ + out <- as.data.frame(matrix(2, nrow =max_t, ncol = 2)) + colnames(out) <- c("diameter","crown_class") ##class 2 is understory, 1 is canopy + rownames(out) <- c(as.character(seq(1,max_t, by = 1))) + + ##use equil_nonclonal to generate dstar and understory diameter sizes + u.growth.list<-equil_nonclonal(dnot=dnot,allometry=allometry,spp = spp,max_t=max_t,timestep = timestep) + dstar<-u.growth.list[["dstar"]] + tstar<-u.growth.list[["tstar"]] + out[1:tstar,"diameter"]<-u.growth.list[["understory_diameter"]] + out[((tstar+1):max_t),"crown_class"]<-1 + for(i in (tstar+1):max_t){ + out[i,"diameter"]<-out[i-1,"diameter"]+calc.dd_nonclonal(diameter=out[i-1,"diameter"],allometry=allometry, spp=spp, understory = FALSE)*timestep + } + if(full.allom){ + out[,"time"]<-NA + out[(tstar+1):max_t,"time"]<-seq(1,max_t-tstar)*timestep + out[,"height"]<-allometry["phi_z",spp]*out[,"diameter"]^allometry["theta_z",spp] + out[,"crown_area"]<-allometry["phi_l",spp]*out[,"diameter"]^allometry["theta_l",spp] + out[,"bolewood_mass"]<-allometry["phi_s",spp]*out[,"diameter"]^allometry["theta_s",spp] + out[1:tstar,"step_fecundity"]<-0 + out[(tstar+1):max_t,"step_fecundity"]<- exp(-allometry["mu_u",spp]*tstar*timestep)*exp(-allometry["mu_c",spp]*out[(tstar+1):max_t,"time"])*allometry["F",spp]*allometry["g",spp]*out[(tstar+1):max_t,"crown_area"]*timestep + out[,"cumsum"]<-cumsum(out[,"step_fecundity"]) + } + return(out) +} + +####eventually, the function below should be able to calculate dD/dt for each stem for n-stemmed plants. I'm still working on the math, and so far +####no.stem=2 for a simple root suckering tree. There can be more variants of this function in the future, but now a clonal tree starts building +####its second trunk when the first trunk reaches size D*. At present, cloning only happens when the primary stem is in full light. + +#'CALC.DD_CLONAL_NURSING +#'function to calculate diameter growth rate dD/dt for a clonal plant's nth trunk/ramet before disconnection happens. +#' +#'@param diameter.prim the diameter size of the primary stem +#'@param diameter.ramet a vector of the diameter sizes of the secondary, third... nth stem +#'@param no.stem total number of stems/ramets in the plant's clone +#'@param allometry an object of class "allometry" +#'@param spp the spp name +#'here the plant's status depends on whether the first ramet/stem reaches dstar. +#' +#'@return returns diameter growth rates for each stem at a particular time point + +calc.dd_clonal_nursing<- function(diameter.prim,diameter.ramet,no.stem=2,allometry,spp){ + spp <- as.character(spp) + CarbonGain<-calc.A_clonal(allometry = allometry,spp=spp) + A<-CarbonGain[spp==spp,"A_c"] + ##canopy diameter growth rate for primary stem + dDdt<-(allometry["phi_l",spp]*diameter.prim^allometry["theta_l",spp]* + (A-(allometry["epsilon_l",spp]+allometry["r_l",spp])* + allometry["l_c",spp]-(allometry["epsilon_r",spp]+allometry["r_r",spp])* + allometry["delta_rl",spp]*allometry["l_c",spp]/allometry["omega",spp]-allometry["f", spp]-allometry["v", spp]))/ + (allometry["theta_s",spp]*allometry["phi_s",spp]*diameter.prim^(allometry["theta_s",spp]-1)+(allometry["sigma",spp]+ + allometry["delta_rl",spp]*allometry["l_c",spp]*allometry["omega",spp])*allometry["phi_l",spp]*allometry["theta_l",spp]*diameter.prim^(-1)) + ##nursed growth rate for secondary stem + A_ramet<-CarbonGain[spp==spp,"A_u"] + dDdt_ramet<-(allometry["v",spp]*allometry["phi_l",spp]*diameter.prim^allometry["theta_l",spp]+allometry["phi_l",spp]*diameter.ramet^allometry["theta_l",spp]* + (A_ramet-(allometry["epsilon_l",spp]+allometry["r_l",spp])* + allometry["l_u",spp]-(allometry["epsilon_r",spp]+allometry["r_r",spp])* + allometry["delta_rl",spp]*allometry["l_u",spp]/allometry["omega",spp]))/ + (allometry["theta_s",spp]*allometry["phi_s",spp]*diameter.ramet^(allometry["theta_s",spp]-1)+(allometry["sigma",spp]+ + allometry["delta_rl",spp]*allometry["l_u",spp]*allometry["omega",spp])*allometry["phi_l",spp]*allometry["theta_l",spp]*diameter.ramet^(-1)) + if(dDdt<=0) { + print("Error: diameter growth cannot be negative") + } + if (dDdt_ramet<=0){ + print("Error: diameter growth cannot be negative") + } + vector<-c(dDdt,dDdt_ramet) + names(vector)<-c("prim","second") + return(vector) +} + +#'GET.GROWTH_CLONAL +#'function to calculate the diameter growth rates dD/dt in the understory and in the canopy based on the threshold diameter size in model_clonal. When a ramet reaches the threshold size, it disconnects from the primary stem. +#' +#'@param dnot the initial diameter size of the primary stem +#'@param dnot_ramet the initial size of a ramet's diameter +#'@param allometry an object of class "allometry" +#'@param spp the spp name +#'@param max_t the total number of time steps of the simulation +#'@param timestep the length of each timestep, in years +#'@param custom a vector supplying custom values for dstar and tstar--the size and time to start cloning as well as end nursing +#'@param full.allom logical, if true returns allometric scaling results along with diameter sizes +#' +#'@return returns a data frame containing diameter, bolewood mass, crown area, height based on defined allometry for the clone. + +get.growth_clonal<-function(dnot,dnot_ramet, allometry,spp,max_t,timestep,custom,full.allom=TRUE){ + out <- as.data.frame(matrix(1, nrow =max_t, ncol = 2)) + colnames(out) <- c("diameter.prim","diameter.ramet") + rownames(out) <- c(as.character(seq(1,max_t, by = 1))) + + ##use equil_nonclonal or customized values to generate dstar and understory diameter sizes + if(!missing(custom)) { + dstar<-custom[1] + tstar<-custom[2] + if(any(custom != "numeric")) { + print("Error: custom variables must be a vector of named numeric vectors") + } + } + u.growth.list<-equil_nonclonal(dnot=dnot,allometry=allometry,spp = spp,max_t=max_t,timestep = timestep) + dstar<-u.growth.list[["dstar"]] + tstar<-u.growth.list[["tstar"]] + out[1:tstar,"diameter.prim"]<-u.growth.list[["understory_diameter"]] + out[1:tstar,"diameter.ramet"]<-0 + out[tstar,"diameter.ramet"]<-dnot_ramet + + ##calculate diameter sizes for the primary and secondary stems during the connected nursing period until the ramet reaches dclonal + for(i in (tstar+1):max_t){ + out[i,]<-out[i-1,]+calc.dd_clonal_nursing(diameter.prim=out[i-1,"diameter.prim"],diameter.ramet=out[i-1,"diameter.ramet"],allometry=allometry, spp=spp)*timestep + if(out[i,"diameter.ramet"]>dstar){ + break + } + } + + ##calculate diameter sizes for the primary and secondary stems during the disconnected period in the canopy + for(t in (i+1):max_t){ + out[t,]<-out[t-1,]+calc.dd_nonclonal(diameter=out[t-1,],allometry=allometry, spp=spp, understory = FALSE)*timestep + } + + ##assign statuses: U is before cloning, N is nursing, C is both stems in canopy + out[,"time"]<-NA ##time is time since canopy + out[(tstar+1):max_t,"time"]<-seq(1,max_t-tstar)*timestep + out[1:tstar,"status"]<-"U" + out[(tstar+1):i,"status"]<-"N" + out[(i+1):max_t,"status"]<-"C" + out[,"step_fecundity.prim"]<-0 + out[,"step_fecundity.ramet"]<-0 + + ##get tbreak + tbreak<-nrow(out[out[,"status"]!="C",]) + + if(full.allom){ + out[,"height.prim"]<-allometry["phi_z",spp]*out[,"diameter.prim"]^allometry["theta_z",spp] + out[,"height.ramet"]<-allometry["phi_z",spp]*out[,"diameter.ramet"]^allometry["theta_z",spp] + out[,"crown_area.prim"]<-allometry["phi_l",spp]*out[,"diameter.prim"]^allometry["theta_l",spp] + out[,"crown_area.ramet"]<-allometry["phi_l",spp]*out[,"diameter.ramet"]^allometry["theta_l",spp] + out[,"bolewood_mass.total"]<-allometry["phi_s",spp]*out[,"diameter.prim"]^allometry["theta_s",spp]+allometry["phi_s",spp]*out[,"diameter.ramet"]^allometry["theta_s",spp] + out[out[,"status"]!="U","step_fecundity.prim"]<- exp(-allometry["mu_u",spp]*tstar*timestep)*exp(-allometry["mu_c",spp]*out[out[,"status"]!="U","time"])*allometry["F",spp]*allometry["g",spp]*out[out[,"status"]!="U","crown_area.prim"]*timestep + out[out[,"status"]=="C","step_fecundity.ramet"]<- exp(-allometry["mu_u",spp]*tstar*timestep)*exp(-allometry["mu_cl",spp]*(tbreak-tstar)*timestep)*exp(-allometry["mu_c",spp]*out[out[,"status"]=="C","time"])*allometry["F",spp]*allometry["g",spp]*out[out[,"status"]=="C","crown_area.ramet"]*timestep + out[,"cumsum.prim"]<-cumsum(out[,"step_fecundity.prim"]) + out[,"cumsum.ramet"]<-cumsum(out[,"step_fecundity.ramet"]) + } + return(out) +} + +##eventually this function will return a vector with all the disconnection times for a multi-stemmed plant if the primary stem did not die in nursing. +##But now it only supports the 2-stem root suckering tree. +#'TBREAK_CLONAL +#'function to calculate the number of timesteps when the disconnection between stems happen. +#' +#'@data ClonalLH a data frame of a clonal plant's life history of three statuses, U,N, and C. Can be generated using get.growth_clonal +#' +#'@return returns the developmentally programmed breaking time. +tbreak_clonal<-function(ClonalLH){ + max(as.integer(rownames(ClonalLH[which(ClonalLH$status=="N"),]))) +} + +##Below are a group of functions to calculate fecundity of the clonal tree during each of its life stage, understory/nursing/canopy depending on +##the time when the primary stem dies and the time when the secondary stem dies. They are by no means efficient, but I couldn't figure out +## a better way to do this besides looking through tables. Please let me know if you have any idea on making this cumbersome calculation of +##LRS easier. Also these functions might belong to methods.R by forming a new class? + +#'NURSING.PRIM +#' +#' function to calculate the fecundity of a clonal tree during its primary stem nursing period +#' +#'@data ClonalLH a data frame of a clonal plant's life history of three statuses, U,N, and C. Can be generated using get.growth_clonal +#'@param tbreak a time point, in number of timesteps when the primary stem disconnects from the secondary trunk and stops nursing it with carbohydrates. It can be the time when the primary stem +#'dies during nursing, the time when the secondary stem dies during nursing, YBD. +#' +#'@return returns the fecundity of a clonal tree during the nursing stage + +Nursing.prim<-function(ClonalLH,tbreak){ + return(ClonalLH$cumsum.prim[tbreak]) +} + +#'SWITCHING.PRIM +#' +#' function to calculate the fecundity of a clonal tree's primary stem during the "switching" stage. If a ramet dies before reaching threshold size dstar, the primary stem stops nursing immediately and switches back to normal reproduction. +#' +#'@data NonclonalLH a data frame of a nonclonal plant's life history +#'@data ClonalLH a data frame of a clonal plant's life history +#'@param tbreak a time point, in number of timesteps, when the primary stem disconnects from the secondary trunk and stops nursing it with carbohydrates +#'@param YMD year mom dies---a time point, in number of timesteps, when the primary stem dies +#' +#'@return returns the fecundity of a clonal tree's primary stem during the switching stage + +Switching.prim<-function(NonclonalLH,ClonalLH,tbreak,YMD){ + #find the diameter size of primary stem at tbreak in the clonal plant's life history data to switch a clonal plant back to non-clonal reproduction in the canopy + breakt<-which(abs(NonclonalLH[,"diameter"]-ClonalLH[tbreak,"diameter.prim"])==min(abs(NonclonalLH[,"diameter"]-ClonalLH[tbreak,"diameter.prim"]))) + ifelse(YMD-breakt<2,return(0),return(NonclonalLH$cumsum[YMD-(tbreak-breakt)])) +} + +#'TBREAK_SWITCHING.RAMET +#' +#' function to calculate the expected disconnection time, tbreak, when a ramet is left alone to grow in the understory without support +#' +#'@data NonclonalLH a data frame of a nonclonal plant's life history +#'@data ClonalLH a data frame of a clonal plant's life history +#'@param YBD year baby dies---a time point, in number of timesteps, when the secondary stem dies +#' +#'@return returns the expected tbreak, actual time needed to reach the canopy +tbreak_switching.ramet<-function(NonclonalLH,ClonalLH,YMD,tstar=tstar){ + breakt<-which(abs(NonclonalLH[,"diameter"]-ClonalLH[YMD,"diameter.ramet"])==min(abs(NonclonalLH[,"diameter"]-ClonalLH[YMD,"diameter.ramet"]))) + return((tstar-breakt)+YMD) +} + +#'SWITCHING.RAMET +#' +#' function to calculate the fecundity of a clonal tree's secondary stem during the "switching" stage. If a primary stem dies during nursing while the ramet survives into the canopy, its fecundity can be calculated. +#' +#'@data NonclonalLH a data frame of a nonclonal plant's life history +#'@data ClonalLH a data frame of a clonal plant's life history +#'@param tbreak a time point, in number of timesteps, when the primary stem disconnects from the secondary trunk and stops nursing it with carbohydrates +#'@param YBD year baby dies---a time point, in number of timesteps, when the secondary stem dies +#' +#'@return returns the fecundity of a clonal tree's ramet during the switching stage + +Switching.ramet<-function(NonclonalLH,ClonalLH,YMD,YBD,tstar=tstar){ + #find the diameter size of ramet at break in the clonal plant's life history data to switch a clonal plant back to non-clonal understory + breaktime<-tbreak_switching.ramet(NonclonalLH = NonclonalLH,ClonalLH = ClonalLH,YMD=YMD,tstar=tstar) + ifelse(YBD-(tstar-breaktime)>nrow(NonclonalLH),return(NonclonalLH$cumsum[max_t]),return(NonclonalLH$cumsum[YBD-(tstar-breaktime)])) +} + +#'SUCCESS.RAMET +#' +#' function to calculate the fecundity of a clonal tree's secondary stem if it successfully get nursed into an independent canopy tree. +#' +#'@data ClonalLH a data frame of a clonal plant's life history +#'@param YBD year baby dies---a time point, in number of timesteps, when the secondary stem dies +#' +#'@return returns the fecundity of a clonal tree's ramet + +Success.ramet<-function(ClonalLH,YBD){ + return(ClonalLH$cumsum.ramet[YBD]) +} + +#'ORDER_SIMULATION +#' +#' function to order the simulation data object +#' currently written like a method, could be moved to a methods file if relevant for multiple classes +#' +#'@param data dataset to reorder +#'@param by dimensional value to reorder by +#'@param decreasing if TRUE, ordered from large to small +#'@param i current timestep +#' +#'@return returns a reordered set of data matrices + +order_simulation <- function(data, by, decreasing = TRUE, i) { + ## currently an inelegant solution with if and, could be better? + if (length(by) == 1) { + data[["spp"]] <- data[["spp"]][order(data[[by]][,i], decreasing = decreasing), ] + data[["diameter"]] <- data[["diameter"]][order(data[[by]][,i], decreasing = decreasing), ] + data[["diameter_growth"]] <- data[["diameter_growth"]][order(data[[by]][,i], decreasing = decreasing), ] + data[["crown_class"]] <- data[["crown_class"]][order(data[[by]][,i], decreasing = decreasing), ] + } + else if (length(by) == 2) { + data[["spp"]] <- data[["spp"]][order(data[[by[1]]][,i], data[[by[2]]][,i], decreasing = decreasing), ] + data[["diameter"]] <- data[["diameter"]][order(data[[by[1]]][,i], data[[by[2]]][,i], decreasing = decreasing), ] + data[["diameter_growth"]] <- data[["diameter_growth"]][order(data[[by[1]]][,i], data[[by[2]]][,i], decreasing = decreasing), ] + data[["crown_class"]] <- data[["crown_class"]][order(data[[by[1]]][,i], data[[by[2]]][,i], decreasing = decreasing), ] + } + else { + print("Error: function currently only supports two by arguments") + return() + } + + return(data) +} + + +#'SIMULATION_AVERAGES +#' +#'this function calculates average values for parameters of interest at each time step. +#' +#'@param simulation an object that is a result of a model_ call +#' +#'@return returns the average dimensional values of the plants at each timestep in the simulation + +## function is of limited utility and we might consider getting rid of it +simulation_averages <- function(simulation) { + + ## check that input is of the class simulation + if(!is.simulation(simulation)) { + print("Error: input must be class simulations") + return() + } + + ## create list of potential variable names, initialize matrix to populate with averages + varnames <- c("diameter", "diameter.growth", "height", "crown_area", "struct_biomass") + cnames <- c("timestep","crown.class", varnames[varnames %in% names(simulation)]) + avgs <- matrix(1, nrow = 2*((simulation[["max_t"]]/simulation[["timestep"]])+1), ncol = length(cnames)) + colnames(avgs) <- cnames + + avgs[,2] <- c(rep(1, times = (simulation[["max_t"]]/simulation[["timestep"]])+1), rep(2, times = (simulation[["max_t"]]/simulation[["timestep"]])+1)) + avgs[,1] <- rep(seq(1, simulation[["max_t"]]+1, by = simulation[["timestep"]]), times = 2) + for (i in 3:length(cnames)) { + for(j in 1:(nrow(avgs)/2)) { + ## no diameter growth in last timestep + if(cnames[i]=="diameter.growth" & j == nrow(avgs)/2) { + avgs[j,i] <- NA + avgs[j+(nrow(avgs)/2),i] <- NA + } + else { + avgs[j,i] <- mean(simulation[[cnames[i]]][simulation[["crown.class"]][,j]==1,j]) + + if (any(simulation[["crown.class"]][,j]==2)) { + avgs[j+(nrow(avgs)/2),i] <- mean(simulation[[cnames[i]]][simulation[["crown.class"]][,j]==2,j]) + } + else avgs[j+(nrow(avgs)/2),i] <- NA + } + } + } + return(avgs) +} diff --git a/man/calcPPA.Rd b/man/calcPPA.Rd new file mode 100644 index 0000000..b7edc8c --- /dev/null +++ b/man/calcPPA.Rd @@ -0,0 +1,14 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/calculatePPA.R +\name{calcPPA} +\alias{calcPPA} +\title{calculate PPA} +\usage{ +calcPPA(dat = NULL) +} +\arguments{ +\item{dat}{a census data.frame with x,y and dbh} +} +\description{ +Calculated canopy layers for census data +} diff --git a/vignettes/example.R b/vignettes/example.R new file mode 100644 index 0000000..aae77a0 --- /dev/null +++ b/vignettes/example.R @@ -0,0 +1,54 @@ + + + +## This script is a sandbox for writing out demo's for team members to easily +## examine functionality. As an example I will write out a basic script to demonstrate +## how the light competition model I wrote works. + +## First, load the package using: + +devtools::load_all(##path/to/your/package/PPA) + + +## Then create an allometry object, which will contain all the parameters needed to simulate +## a stand with two species. +## The only argument you need to supply is the number of species. The +## stochastic argument, which is TRUE by default, automatically injects +## some variation into the allometries (right now in a very crude way). +## There is an additional parameter, custom, that can be used to provide +## custom allometries. Eventually, it could be nice to build in some set +## species allometries +test_allometry <- create_allometry.model_light(n.spp = 2, stochastic = TRUE) + +## print test allometry to look inside. Note that its just a named matrix containing +## the relevant parameters. +test_allometry + + +## Next, we can run a simulation using the function: model_light. +test_simulation <- model_light(allometry = test_allometry[["allometry"]], ## supply the allometry object we just created + spp = c(1, 2), ## spp must be a vector of spp names that matches the names in the allometry object + max_t = 100, ## last timestep in simulation + timestep = 1, ## length of a timestep + plot.area = 1000, ## how big is the simulation area, in square meters + n_start = list("1" = 25, "2" = 25), ## how many individuals do we start with + full.data = TRUE, ## do we want the model output to include all the steps on the way? + full.allom = TRUE) ## dow we want the model output to include all of the allometric components? (i.e. biomass, height etc) + +## now we can take a peak at the model_light output. As you can see it is a list of multiple datasets +## This is it in its most verbose form, as we specified both full.data and full.allom +str(test_simulation) + + +## now we can use the summary method to check out what we got: +summary(test_simulation) + +## currently the summary.model_light method prints some basic attributes about the simulation as well as the average values +## for some of the dimensional attributes of the individuals. I don't necessarily think this is the best stuff to display, but +## I think deciding that should be a collaborative effort. + +## finally we can plot the results of the our simulation using plot() +plot(test_simulation) + +## This plot output plots the trajectories of the individual trees through time. Again, the plot output +## should probably be revised collaboratively