Nothing
R/make_data.R
In TTR.PGM: Thornley Transport Resistance Plant Growth Model
Defines functions
make_data
required_parameters
add_requirements
parameter_groups
replace.na
Documented in
make_data
replace.na <- function(row_bounds, group_name){
ind.group <- grep(group_name, names(row_bounds))
if(length(ind.group) == 0) return(row_bounds)
if(!any(is.na(row_bounds[ind.group]))) return(row_bounds)
if(length(ind.group) != 4 & length(ind.group) != 2) {
vecprt <- paste0("c(", paste0(names(row_bounds), "=", row_bounds, collapse=", "), ")")
warning(paste0("group ", group_name, " in parameters ", vecprt, " indicates neither trapezoid nor step function."))
}
# trapezoid
if(length(ind.group) == 4) {
row_bounds[ind.group] = c(-1e20,10000,2e20-10000,10000)
return(row_bounds)
}
#loss rate has to be 0, other effects are neutral at 1.
if(group_name == "r_tloss"){
row_bounds[ind.group] = c(1e20,10000)
} else {
row_bounds[ind.group] = c(-1e20, 10000)
}
return(row_bounds)
}
#this is a list which contains all on/off options in the options list used for TTR_Model
#used for checking input for required optimisation parameters, see below
option_switches <- c("pe","initial_mass_par")
#these define the required optimisation parameters for the different versions of the model
#used for checking if bounds for these are defined in input
#Parameters are defined in the enumeration structure in enum struct BetaParameters inst/include/definitions.hpp
required_parameters_list <- list(
base = list( # in all models
beta = c(
"CU_swc_1",
"CU_swc_2",
"CU_Ns_1",
"CU_Ns_2",
"NU_tnur_1",
"NU_tnur_2",
"NU_swc_1",
"NU_swc_2",
"NU_swc_3",
"NU_swc_4",
"NU_Nsoil_1",
"NU_Nsoil_2",
"g_tgrowth_1",
"g_tgrowth_2",
"g_tgrowth_3",
"g_tgrowth_4"
),
alpha = c()
),
std = list(
beta = c(
"CU_tcur_1",
"CU_tcur_2",
"CU_tcur_3",
"CU_tcur_4",
"CU_ppfd_1",
"CU_ppfd_2",
"L_mix_1",
"r_tloss_1",
"r_tloss_2"
),
alpha = c()
),
fqr = list(
beta = c(
"L_mix_1",
"r_tloss_1",
"r_tloss_2"
),
alpha = c()
),
red = list(
beta = c(
"g_swc_1",
"g_swc_2",
"L_mix_1",
"r_tloss_1",
"r_tloss_2"
),
alpha = c()
),
oak = list(
beta = c(
"g_swc_1",
"g_swc_2",
"L_swc_1",
"L_swc_2",
"L_tloss_1",
"L_tloss_2",
"L_mix_1",
"L_mix_2",
"L_mix_3"
),
alpha = c()
),
swc_online = list(
beta = c(),
alpha = c(
"AET_bsum_1",
"AET_bsum_2"
)
),
lc = list(
beta = c(),
alpha = c(
"LC_bsum_1",
"LC_bsum_2"
)
),
fire = list(
beta = c(
"L_fire_1",
"L_fire_2"
),
alpha = c()
),
phen = list(
beta = c(
"phen_t_0",
"phen_t_1",
"phen_t_2",
"phen_t_3",
"phen_t_4"
),
alpha = c()
),
pe = list(
beta = c(
"pe_Ms",
"pe_Mr",
"pe_Cs",
"pe_Cr",
"pe_Ns",
"pe_Nr"
)
),
initial_mass_par = list(
alpha = c(),
beta = c("iM")
)
)
#These entries define which parameters should be interpreted as a group.
#The parameter names should be "<Group name>_<number>"
standard_groups <- c(
"CU_tcur",
"CU_ppfd",
"CU_swc",
"CU_Ns",
"NU_tnur",
"NU_Nsoil",
"NU_swc",
"g_tgrowth",
"r_tloss",
"L_tloss",
"L_swc",
"phen_t",
"L_fire",
"LC_bsum",
"AET_bsum",
"g_swc"
)
parameter_groups <- function(bounds, traps){
out <- list(alpha=list(), beta=list)
for(c in c("alpha","beta")){
bnames <- names(bounds[[c]])
groups <- lapply(traps, function(t){
grep(paste0("^",t,"_[1-9][0-9]*"),bnames)
})
groups <- groups[lapply(groups, length) > 0]
out[[c]] <- groups
}
out
}
add_requirements <- function(req, additional_req){
req$alpha <- c(req$alpha, additional_req$alpha)
req$beta <- c(req$beta, additional_req$beta)
return(req)
}
#this is the function that returns all required parameters for a model configuration as given in options
#used for checking if all required bounds are given
required_parameters <- function(options){
#add required bounds according to version
req <- required_parameters_list[["base"]]
req <- add_requirements(req, required_parameters_list[[options[["ve"]]]])
#add required bounds for each option
for(opt in option_switches){
if(options[[opt]]) {
req <- add_requirements(req, required_parameters_list[[opt]])
}
}
req
}
#' Standard values for global parameters
#'
#' This vector contains the global parameter values for the process model,
#' which must be supplied to the \code{make_data()} function.
#' The time components are the units are per modelled time step and the length of these time
#' steps are defined by the parameter seconds in \code{make_data()}.
#' In the descriptions below M indicates
#' dry mass of structural biomass (either shoot mass MS or root mass MR).
#'
#' @section Units:
#' \describe{
#' \item{mmax}{Loss coefficient, per timestep (kg / kg M per timestep)}
#' \item{gmax}{Growth coefficient, in (kg C kg N kg M^-2)^-1 per timestep}
#' \item{KM}{Size dependency of mass-loss (kg M)}
#' \item{CUmax}{Carbon uptake rate (kg C / kg shoot M per timestep)}
#' \item{NUmax}{Nitrogen uptake rate (kg N / kg shoot M per timestep)}
#' \item{KA}{Size dependency of uptake rates (kg M)}
#' \item{Jc}{Substrate inhibition coefficient for carbon (kg C / kg M)}
#' \item{Jn}{Substrate inhibition coefficient for nitrogen (in kg N / kg M)}
#' \item{q}{Scaling coefficient for substrate transport, no units}
#' \item{RHOc}{Specific carbon transport resistance, per timestep}
#' \item{RHOn}{Specific nitrogen transport resistance, per timestep}
#' \item{Fc}{Fraction of carbon in M (kg C / kg M)}
#' \item{Fn}{Fraction of nitrogen in M (kg N / kg M)}
#' }
#' @export
standard_globals <- c(
mmax = 0.1, # 0.05 per day (proportion)
gmax = 20, # 200 [kg C kg N kg M^-2]^-1 per day
KM = 1.5, # Michaelis-Menten constant for litter (kg M)
CUmax = 0.10, # Carbon uptake (kg C / kg shoot M / day)
NUmax = 0.02, # Nitrogen uptake (kg N / kg shoot M / day)
KA = 3, # Substrate Michaelis constant (kg M)
Jc = 0.1, # Carbon cost (kg C / kg M)
Jn = 0.01, # Nitrogen cost (kg N / kg M)
q = 1.0, # No units
RHOc = 1.0, # Specific carbon transport resistance, per day
RHOn = 1.0, # Specific nirogen transport resistance, per day
Fc = 0.5, # Fraction of carbon (kg C / kg M)
Fn = 0.025 # Fraction of nitrogen (kg N / kg M)
)
#' Standard values for bounds
#' @description This list contains the standard bounds of parameters used in the Model.
#' See \code{make_data()} for more information on setting bounds.
#' @export
standard_bounds = list(
alpha = list(
LC_bsum_1 = c( 0, 1*10 ),
LC_bsum_2 = c( 0, 1*10 ),
AET_bsum_1 = c( 0, 1*10 ),
AET_bsum_2 = c( 0, 1*10 ),
scale = c(0, 50)
),
beta = list(
CU_tcur_1 = c( -50, 50 ),
CU_tcur_2 = c( 0, 50 ),
CU_tcur_3 = c( 0, 50 ),
CU_tcur_4 = c( 0, 50 ),
CU_ppfd_1 = c( 0, 2000),
CU_ppfd_2 = c( 0, 2000),
CU_swc_1 = c( 0, 100 ),
CU_swc_2 = c( 0, 100 ),
CU_Ns_1 = c( 0, 1 ),
CU_Ns_2 = c( 0, 1 ),
NU_tnur_1 = c( -50, 50 ),
NU_tnur_2 = c( 0, 50 ),
NU_Nsoil_1 = c( 0, 1 ),
NU_Nsoil_2 = c( 0, 1 ),
NU_swc_1 = c( 0, 100 ),
NU_swc_2 = c( 0, 100 ),
NU_swc_3 = c( 0, 100 ),
NU_swc_4 = c( 0, 100 ),
g_tgrowth_1 =c( -50, 50 ),
g_tgrowth_2 =c( 0, 50 ),
g_tgrowth_3 =c( 0, 50 ),
g_tgrowth_4 =c( 0, 50 ),
g_swc_1 = c( 0, 100 ),
g_swc_2 = c( 0, 100 ),
r_tloss_1 = c( -50, 50 ),
r_tloss_2 = c( 0, 50 ),
L_tloss_1 = c( -50, 50 ),
L_tloss_2 = c( 0, 50 ),
L_swc_1 = c( 0, 100 ),
L_swc_2 = c( 0, 100 ),
phen_t_0 = c( 0, 1 ),
phen_t_1 = c( 0, 365 ),
phen_t_2 = c( 0, 365 ),
phen_t_3 = c( 0, 365 ),
phen_t_4 = c( 0, 365 ),
L_fire_1 = c( 0, 100 ),
L_fire_2 = c( 0, 100 ),
L_mix_1 = c( 0, 100 ),
L_mix_2 = c( 0, 100 ),
L_mix_3 = c( 0, 100 ),
iM = c( 0, 100 ),
pe_Ms = c( 0, 100 ),
pe_Mr = c( 0, 100 ),
pe_Cs = c( 0, 100 ),
pe_Cr = c( 0, 100 ),
pe_Ns = c( 0, 100 ),
pe_Nr = c( 0, 100 )
)
)
#' Standard options for the model
#' @description This list contains the standard options used by \code{make_data()}.
#' It is provided here as an example. See \code{options} for more information.
#' @export
standard_options <- list(
optim = "DEoptim",
ve = "red",
#output_steps
steps = c(50,100),
initial_mass = 0.5,
initial_mass_par = FALSE, # set initial mass as estimated parameter TRUE or FALSE
pe = FALSE, # consider process error TRUE or FALSE
species = "species_1",
photo = "c3", #for multi-species a vector of length SPECIES
pe_scale = 1e-04
)
#' Options to be used in the \code{make_data()} function
#'
#' @name options
#' @title Option list used in the \code{make_data()} function
#'
#' @param optim
#' String specifying the data calculated by the user defined Model function, e.g.
#' one of \code{c("DEoptim", "LaplacesDemon")}.
#'
#' @param ve
#' String specifying the version of the process model. Must be one of \code{c("std", "fqr", "red", "oak")}.
#'
#' @param steps
#' Integer numeric vector specifying the timesteps to output. The process model
#' will run \code{max(steps)} iterations and the output will contain
#' \code{length(steps)} output times.
#'
#' @param initial_mass_par
#' Boolean specifying if the initial biomass used in the process model should be
#' estimated as a parameter.
#'
#' @param initial_mass
#' The initial biomass to use for the process model, when \code{initial_mass_par} is set
#' to \code{FALSE}.
#'
#' @param species
#' A vector of names of species.
#'
#' @param photo
#' A vector specifying the photosynthesis system used by the species.
#' Must be either \code{c("c3", "c4")} for each species.
#'
#' @param swc_online
#' Boolean specifying if soil water content should be calculated online as part of the
#' process model (\code{TRUE}) or if it is given as a forcing parameter. Currently not supported.
#'
#' @param lc
#' Boolean turning light competition simulation in the process model on/off. Currently not supported.
#'
#'
#' @param fire
#' Boolean turning fire in the process model on/off. Currently not supported.
#'
#' @param pe
#' Boolean turning process error simulation on/off.
#'
#' @param pe_scale
#' The process error scaling factor.
#'
#' @description
#' The options list handed to \code{make_data()}. It specifies the process model options.
#' An example is provided by \code{standard_options}.
NULL
#' Create data object for use with the process model
#' @description This function creates an object representing a configuration of the process model including forcing data.
#' This object is used by \code{run_ttr()}.
#' @param input The forcing data object created by a call to \code{get_input()}, see specification in its help page
#' @param options An options list for running the model, see \code{options} help page
#' @param globals Optional parameter specifying global variables used in the simulation, see \code{standard_globals} help page
#' @param bounds Optional parameter specifying the upper and lower bounds of the model's parameters, see \code{make_data()}. Bounds for additional parameters included in the statistical model need to be specified here. Bounds of the process model parameters can also be modified from their default values.
#' @param groups Optional parameter specifying the groups of process model parameters
#' @param verbose Print out a log detailing checks
#' @return A ttr data object, a list with the following elements:
#' \item{PGF}{Function used by LaplacesDemon for generating initial values.
#' }
#' \item{options}{Options as specified in the \code{options} and \code{standard_options} help pages
#' }
#' \item{globals}{Globals as specified in the \code{standard_globals} help page
#' }
#' \item{n.sites}{Number of sites
#' }
#' \item{n.species}{Number of species
#' }
#' \item{n.time}{Number of timesteps
#' }
#' \item{n.parm}{Number of parameters to be fitted
#' }
#' \item{n.parm.a}{Number of per-process parameters
#' }
#' \item{n.parm.b}{Number of per-species parameters
#' }
#' \item{parm.names}{Names of all parameters for this model configuration
#' }
#' \item{parm.names.a}{Names of per-process parameter names
#' }
#' \item{parm.names.b}{Names of per-species parameter names
#' }
#' \item{mon.names}{Required by LaplacesDemon, always set to "LP"
#' }
#' \item{out}{Names of the response variables returned by \code{run_ttr()}
#' }
#' \item{out.dimnames}{Dimension names for the output returned by \code{run_ttr()}
#' }
#' \item{out.dim}{Dimensions for the output returned by \code{run_ttr()}
#' }
#' \item{dimnames.beta}{Dimension names for the ttr parameters object produced by \code{get_parms()}
#' }
#' \item{bounds}{Bounds for all process model parameters
#' }
#' \item{timeseries}{See \code{get_input()}
#' }
#' \item{timeinvariant}{See \code{get_input()}
#' }
#' \item{y}{observations from \code{get_input()} cast to a matrix
#' }
#' \item{lonlat}{a matrix with lon and lat columns for the nsites rows
#' }
#' \item{pos.oe}{Indices of observation error parameters
#' }
#' \item{pos.pe}{Indices of process error parameters
#' }
#' \item{parm.trap.groups}{Representation of parameters that form a group for the trapezoid functions
#' }
#' \item{arguments}{The unprocessed arguments to this function, excluding data
#' }
#' @details The data structure produced in this function can be handed to an external parameter
#' estimation algorithm e.g. LaplacesDemon or DEoptim to fit the model.
#'
#' The bounds of the process model's parameters need to be defined.
#' The user can supply these to \code{make_data()} or accept the default
#' values returned by a call to \code{make_data()}. Different model variants use
#' different parameters and a call to \code{make_data()} returns the default bounds for the
#' model variant specified in \code{options}.
#'
#' bounds is a list of two lists; one for alpha (per-process) parameters and
#' one for beta (per-species) parameters. The list bounds is formatted as:
#' bounds = list(alpha = list(alpha_parameter = (lower, upper)),
#' beta = list(beta_parameter = c(lower, upper))), where `upper` specifies
#' the maximum value a process model parameter can take on, and
#' `lower` is the corresponding minimum.
#'
#' The external parameter estimation algorithm should then take Data$bounds[,1] as
#' lower and Data$bounds[,2] as the upper bounds.
#'
#' In case the user does not wish to estimate some of the parameters in the
#' external optimisation process, the upper and lower bound of a parameter
#' should be set to the same numeric value. In this case, the parameter will
#' not be part of the Data$bounds matrix passed to the external parameter estimation
#' algorithm, and its value is set internally in the \code{get_parms()} function.
#'
#' For any group of parameters, e.g., CU_Ns_1 and CU_Ns_2, setting any one
#' of the parameters to c(NA, NA) will set the whole group of parameters to
#' a constant value that will negate the influence of this parameter group on the
#' process model. This effectively switches off the effect of these parameters
#' on the model's processes and removes them from Data$bounds. The constant
#' values the parameters are set to can be observed by inspecting the output
#' of \code{get_parms()} for a Data object.
#'
#' The names of the parameters are written as YY_XX_i, where YY specifies the
#' process being considered and XX the environmental driver and i the parameter
#' number in the step or trapezoid function, e.g. CU_swc_1 and CU_swc_2 specify
#' how CU (carbon uptake) is influenced by swc (soil water content) and 1 and 2
#' specify the lower and upper parameters of the step function.
#'
#' An example for specifying all process bounds for a Model run with the
#' standard options is given in the object `standard_bounds`.
#'
#' @importFrom stats rnorm
#' @examples
#' #some dummy data as input
#' input_data <- get_input(
#' observations = c(1),
#' tcur = c(1),
#' tnur = c(1),
#' tgrowth = c(1),
#' tloss = c(1),
#' seconds = c(1),
#' lon = c(1),
#' lat = c(1),
#' rsds = c(1),
#' catm = c(1),
#' pet = c(1),
#' rain = c(1),
#' wp = c(1),
#' fc = c(1)
#' )
#' data <- make_data(
#' input = input_data,
#' options = standard_options,
#' globals = standard_globals,
#' bounds = list(
#' alpha = list(
#' #defining a new parameter, e.g. as part of a model function
#' my_parameter = c(1,2),
#' #overwriting a standard bound
#' scale = c(0,500)
#' ),
#' #no more parameters than needed for the process model in beta, for these standard values.
#' beta = list()
#' )
#' )
#' #all standard values
#' data <- make_data(input_data, standard_options)
#' @export
make_data <- function(input, options=standard_options, globals=standard_globals, bounds=list(alpha = NULL, beta=NULL), groups=c(), verbose=FALSE){
arguments <- list(options=options, globals=globals, bounds=bounds, groups=groups)
ts <- input$timeseries
time.invariant <- input$time.invariant
observations <- input$observations
in_data_dimensions <- input$dimensions
#options
assertList(options)
assertChoice(options[["optim"]], c("DEoptim", "LaplacesDemon", "no"))
assertIntegerish(
options[["steps"]],
any.missing=FALSE,
min.len=1,
lower=0
)
options[["steps"]] <- unique(sort(options[["steps"]]))
assertChoice(options[["ve"]],varnames("TTRVariant"))
assertNumber(
options[["initial_mass"]],
lower=0.0,
finite=T,
)
for(boolprop in option_switches){
assertFlag(options[[boolprop]], .var.name=paste0("options$",boolprop))
}
assertCharacter(
options[["species"]],
min.len = 1,
any.missing = F,
unique = T
)
if(!testFactor(options[["photo"]],
levels=varnames("PhotoType"))){
assertCharacter(
options[["photo"]],
len = length(options[["species"]]),
any.missing = F
)
for(photo_type in options[["photo"]]){
assertChoice(
photo_type,
choices = varnames("PhotoType")
)
}
options$photo <- factor(options$photo, levels = varnames("PhotoType"))
}
assertFactor(
options[["photo"]],
levels=varnames("PhotoType"),
len=length(options[["species"]]),
any.missing=FALSE
)
assertNumber(
options[["pe_scale"]],
finite=T
)
assertList(
bounds
)
assertNames(
names(bounds),
permutation.of = c("alpha","beta")
)
if(!is.null(bounds$alpha)){
assertList(bounds$alpha)
}
if(!is.null(bounds$beta)){
assertList(bounds$beta)
}
if(length(bounds$alpha) > 0 && length(bounds$beta) > 0){
assertNames(
names(bounds$alpha),
disjunct.from = names(bounds$beta)
)
}
req <- required_parameters(options)
missing_alpha_ind <- which(!(req$alpha %in% names(bounds$alpha)))
missing_beta_ind <- which(!(req$beta %in% names(bounds$beta)))
if(length(missing_alpha_ind) > 0 || length(missing_beta_ind) > 0){
if(verbose){
cat(paste0("Missing parameter bounds for required parameters, using standard bounds for: ",
paste(req$beta[missing_beta_ind], collapse = ", "), ", ",
paste(req$alpha[missing_alpha_ind], collapse = ", "), "\n"))
}
}
for(missing_beta in req$beta[missing_beta_ind]){
bounds$beta[missing_beta] <- standard_bounds$beta[missing_beta]
}
for(missing_alpha in req$alpha[missing_alpha_ind]){
bounds$alpha[missing_alpha] <- standard_bounds$alpha[missing_alpha]
}
bounds$beta <- bounds$beta[sort(names(bounds$beta))]
bounds$alpha <- bounds$alpha[sort(names(bounds$alpha))]
#parameter groups
if(!is.null(groups)){
assertCharacter(
groups,
any.missing=FALSE
)
}
groups <- unique(c(groups, standard_groups))
#globals
assertNumeric(
globals,
any.missing=FALSE,
names="unique"
)
assertNames(
names(globals),
type="unique",
permutation.of=varnames("Global")
)
globals <- globals[varnames("Global")]
#Optimisation: q has to be 1, a different value requires recompilation
assertNumber(
globals[["q"]],
lower = 1.0,
upper = 1.0
)
#time series
assertArray(
ts,
mode="double",
any.missing=FALSE,
d=3,
.var.name="time dependent forcing variables"
)
assertNames(
dimnames(ts)[[1]],
type="unique",
identical.to=varnames("TimeSeries"),
.var.name="forcing variable names"
)
#time invariant
assertArray(
time.invariant,
mode="double",
any.missing=FALSE,
d=2,
.var.name="time-independent forcing variables"
)
assertNames(
dimnames(time.invariant)[[1]],
type="unique",
identical.to=varnames("TimeInvariant")
)
if(dim(ts)[[3]] != in_data_dimensions[["nsites"]]){
stop(paste0("Time series array has an incorrect number of sites: ",dim(ts)[[3]], " instead of ", in_data_dimensions[["nsites"]]))
}
if(dim(time.invariant)[[2]] != in_data_dimensions[["nsites"]]){
stop(paste0("Time independent variable array has an incorrect number of sites: ",dim(time.invariant)[[2]], " instead of ", in_data_dimensions[["nsites"]]))
}
#cat("\nchecks complete\n")
#From here, assume all bounds are specifying a parameter in the model
pn <- names(bounds)
nspecies <- length(options$species)
lower_of <- function(bounds){
vapply(bounds, function(entry){entry[[1]]}, FUN.VALUE=c(0.0))
}
upper_of <- function(bounds){
vapply(bounds, function(entry){entry[[2]]}, FUN.VALUE=c(0.0))
}
#replace NAs in bounds with constant values
lower.a <- lower_of(bounds$alpha)
upper.a <- upper_of(bounds$alpha)
for(group in groups){
lower.a <- replace.na(lower.a, group)
upper.a <- replace.na(upper.a, group)
}
lower.b.single <- lower_of(bounds$beta)
upper.b.single <- upper_of(bounds$beta)
for(group in groups){
lower.b.single <- replace.na(lower.b.single, group)
upper.b.single <- replace.na(upper.b.single, group)
}
lower.b <- rep(lower.b.single, nspecies)
upper.b <- rep(upper.b.single, nspecies)
lower <- c(lower.a, lower.b)
upper <- c(upper.a, upper.b)
#remove constant parameters from bounds structure
ind.parm.constant <- which(lower == upper)
parm.constant <- c()
if(length(ind.parm.constant) > 0){
parm.constant <- lower[ind.parm.constant]
lower <- lower[-ind.parm.constant]
upper <- upper[-ind.parm.constant]
}
out.dimnames <- list(
varnames("Biotic"),
options$species,
as.character(options$steps),
as.character(1:in_data_dimensions[["nsites"]])
)
out.dim <- c(
length(varnames("Biotic")),
length(options$species),
length(options$steps),
in_data_dimensions[["nsites"]]
)
names(out.dim) <- c(
"output",
"species",
"steps",
"sites"
)
PGF <- function(Data){
lo <- Data$bounds[,1]
up <- Data$bounds[,2]
mu <- (lo+up)/2
pmin(pmax(rnorm(Data$n.parm, mu, 0.1), lo), up)
}
opt.bounds <- matrix(data=c(lower,upper),ncol=2,dimnames=list(names(lower),c("lower","upper")))
n.constant <- length(ind.parm.constant)
n.constant.a <- length(which(ind.parm.constant <= length(bounds$alpha)))
n.constant.b <- length(which(ind.parm.constant > length(bounds$alpha)))
parm.names <- names(lower)
data <- structure(
list(
PGF=PGF,
options=options,
globals=globals,
n.sites = in_data_dimensions[["nsites"]],
n.species = length(options$species),
n.time = in_data_dimensions[["nsteps"]],
n.parm = dim(opt.bounds)[1],
n.parm.proc = length(bounds$alpha) + nspecies*length(bounds$beta),
n.parm.a = length(bounds$alpha) - n.constant.a,
n.parm.proc.a = length(bounds$alpha),
n.parm.b = length(bounds$beta) - n.constant.b,
n.parm.proc.b = length(bounds$beta),
ind.parm.constant = ind.parm.constant,
parm.constant = parm.constant,
out.dimnames = out.dimnames,
out.dim = out.dim,
parm.names.a=names(bounds$alpha),
parm.names.b=names(bounds$beta),
parm.names=parm.names,
dimnames.beta = list(names(bounds$beta),options$species),
parm.trap.groups = parameter_groups(bounds, groups),
pos.oe = grep("^oe", bounds$alpha),
pos.pe = grep("^pe", bounds$beta),
mon.names = "LP",
out=varnames("Biotic"),
bounds=opt.bounds,
timeseries=ts,
time.invariant=time.invariant,
lonlat=input$lonlat,
y=observations,
arguments=arguments),
class = "ttr.data")
return(data)
}
Try the TTR.PGM package in your browser
Any scripts or data that you put into this service are public.
TTR.PGM documentation built on June 8, 2025, 9:32 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions make_data required_parameters add_requirements parameter_groups replace.na
Documented in make_data
replace.na <- function(row_bounds, group_name){
ind.group <- grep(group_name, names(row_bounds))
if(length(ind.group) == 0) return(row_bounds)
if(!any(is.na(row_bounds[ind.group]))) return(row_bounds)
if(length(ind.group) != 4 & length(ind.group) != 2) {
vecprt <- paste0("c(", paste0(names(row_bounds), "=", row_bounds, collapse=", "), ")")
warning(paste0("group ", group_name, " in parameters ", vecprt, " indicates neither trapezoid nor step function."))
}
# trapezoid
if(length(ind.group) == 4) {
row_bounds[ind.group] = c(-1e20,10000,2e20-10000,10000)
return(row_bounds)
}
#loss rate has to be 0, other effects are neutral at 1.
if(group_name == "r_tloss"){
row_bounds[ind.group] = c(1e20,10000)
} else {
row_bounds[ind.group] = c(-1e20, 10000)
}
return(row_bounds)
}
#this is a list which contains all on/off options in the options list used for TTR_Model
#used for checking input for required optimisation parameters, see below
option_switches <- c("pe","initial_mass_par")
#these define the required optimisation parameters for the different versions of the model
#used for checking if bounds for these are defined in input
#Parameters are defined in the enumeration structure in enum struct BetaParameters inst/include/definitions.hpp
required_parameters_list <- list(
base = list( # in all models
beta = c(
"CU_swc_1",
"CU_swc_2",
"CU_Ns_1",
"CU_Ns_2",
"NU_tnur_1",
"NU_tnur_2",
"NU_swc_1",
"NU_swc_2",
"NU_swc_3",
"NU_swc_4",
"NU_Nsoil_1",
"NU_Nsoil_2",
"g_tgrowth_1",
"g_tgrowth_2",
"g_tgrowth_3",
"g_tgrowth_4"
),
alpha = c()
),
std = list(
beta = c(
"CU_tcur_1",
"CU_tcur_2",
"CU_tcur_3",
"CU_tcur_4",
"CU_ppfd_1",
"CU_ppfd_2",
"L_mix_1",
"r_tloss_1",
"r_tloss_2"
),
alpha = c()
),
fqr = list(
beta = c(
"L_mix_1",
"r_tloss_1",
"r_tloss_2"
),
alpha = c()
),
red = list(
beta = c(
"g_swc_1",
"g_swc_2",
"L_mix_1",
"r_tloss_1",
"r_tloss_2"
),
alpha = c()
),
oak = list(
beta = c(
"g_swc_1",
"g_swc_2",
"L_swc_1",
"L_swc_2",
"L_tloss_1",
"L_tloss_2",
"L_mix_1",
"L_mix_2",
"L_mix_3"
),
alpha = c()
),
swc_online = list(
beta = c(),
alpha = c(
"AET_bsum_1",
"AET_bsum_2"
)
),
lc = list(
beta = c(),
alpha = c(
"LC_bsum_1",
"LC_bsum_2"
)
),
fire = list(
beta = c(
"L_fire_1",
"L_fire_2"
),
alpha = c()
),
phen = list(
beta = c(
"phen_t_0",
"phen_t_1",
"phen_t_2",
"phen_t_3",
"phen_t_4"
),
alpha = c()
),
pe = list(
beta = c(
"pe_Ms",
"pe_Mr",
"pe_Cs",
"pe_Cr",
"pe_Ns",
"pe_Nr"
)
),
initial_mass_par = list(
alpha = c(),
beta = c("iM")
)
)
#These entries define which parameters should be interpreted as a group.
#The parameter names should be "<Group name>_<number>"
standard_groups <- c(
"CU_tcur",
"CU_ppfd",
"CU_swc",
"CU_Ns",
"NU_tnur",
"NU_Nsoil",
"NU_swc",
"g_tgrowth",
"r_tloss",
"L_tloss",
"L_swc",
"phen_t",
"L_fire",
"LC_bsum",
"AET_bsum",
"g_swc"
)
parameter_groups <- function(bounds, traps){
out <- list(alpha=list(), beta=list)
for(c in c("alpha","beta")){
bnames <- names(bounds[[c]])
groups <- lapply(traps, function(t){
grep(paste0("^",t,"_[1-9][0-9]*"),bnames)
})
groups <- groups[lapply(groups, length) > 0]
out[[c]] <- groups
}
out
}
add_requirements <- function(req, additional_req){
req$alpha <- c(req$alpha, additional_req$alpha)
req$beta <- c(req$beta, additional_req$beta)
return(req)
}
#this is the function that returns all required parameters for a model configuration as given in options
#used for checking if all required bounds are given
required_parameters <- function(options){
#add required bounds according to version
req <- required_parameters_list[["base"]]
req <- add_requirements(req, required_parameters_list[[options[["ve"]]]])
#add required bounds for each option
for(opt in option_switches){
if(options[[opt]]) {
req <- add_requirements(req, required_parameters_list[[opt]])
}
}
req
}
#' Standard values for global parameters
#'
#' This vector contains the global parameter values for the process model,
#' which must be supplied to the \code{make_data()} function.
#' The time components are the units are per modelled time step and the length of these time
#' steps are defined by the parameter seconds in \code{make_data()}.
#' In the descriptions below M indicates
#' dry mass of structural biomass (either shoot mass MS or root mass MR).
#'
#' @section Units:
#' \describe{
#' \item{mmax}{Loss coefficient, per timestep (kg / kg M per timestep)}
#' \item{gmax}{Growth coefficient, in (kg C kg N kg M^-2)^-1 per timestep}
#' \item{KM}{Size dependency of mass-loss (kg M)}
#' \item{CUmax}{Carbon uptake rate (kg C / kg shoot M per timestep)}
#' \item{NUmax}{Nitrogen uptake rate (kg N / kg shoot M per timestep)}
#' \item{KA}{Size dependency of uptake rates (kg M)}
#' \item{Jc}{Substrate inhibition coefficient for carbon (kg C / kg M)}
#' \item{Jn}{Substrate inhibition coefficient for nitrogen (in kg N / kg M)}
#' \item{q}{Scaling coefficient for substrate transport, no units}
#' \item{RHOc}{Specific carbon transport resistance, per timestep}
#' \item{RHOn}{Specific nitrogen transport resistance, per timestep}
#' \item{Fc}{Fraction of carbon in M (kg C / kg M)}
#' \item{Fn}{Fraction of nitrogen in M (kg N / kg M)}
#' }
#' @export
standard_globals <- c(
mmax = 0.1, # 0.05 per day (proportion)
gmax = 20, # 200 [kg C kg N kg M^-2]^-1 per day
KM = 1.5, # Michaelis-Menten constant for litter (kg M)
CUmax = 0.10, # Carbon uptake (kg C / kg shoot M / day)
NUmax = 0.02, # Nitrogen uptake (kg N / kg shoot M / day)
KA = 3, # Substrate Michaelis constant (kg M)
Jc = 0.1, # Carbon cost (kg C / kg M)
Jn = 0.01, # Nitrogen cost (kg N / kg M)
q = 1.0, # No units
RHOc = 1.0, # Specific carbon transport resistance, per day
RHOn = 1.0, # Specific nirogen transport resistance, per day
Fc = 0.5, # Fraction of carbon (kg C / kg M)
Fn = 0.025 # Fraction of nitrogen (kg N / kg M)
)
#' Standard values for bounds
#' @description This list contains the standard bounds of parameters used in the Model.
#' See \code{make_data()} for more information on setting bounds.
#' @export
standard_bounds = list(
alpha = list(
LC_bsum_1 = c( 0, 1*10 ),
LC_bsum_2 = c( 0, 1*10 ),
AET_bsum_1 = c( 0, 1*10 ),
AET_bsum_2 = c( 0, 1*10 ),
scale = c(0, 50)
),
beta = list(
CU_tcur_1 = c( -50, 50 ),
CU_tcur_2 = c( 0, 50 ),
CU_tcur_3 = c( 0, 50 ),
CU_tcur_4 = c( 0, 50 ),
CU_ppfd_1 = c( 0, 2000),
CU_ppfd_2 = c( 0, 2000),
CU_swc_1 = c( 0, 100 ),
CU_swc_2 = c( 0, 100 ),
CU_Ns_1 = c( 0, 1 ),
CU_Ns_2 = c( 0, 1 ),
NU_tnur_1 = c( -50, 50 ),
NU_tnur_2 = c( 0, 50 ),
NU_Nsoil_1 = c( 0, 1 ),
NU_Nsoil_2 = c( 0, 1 ),
NU_swc_1 = c( 0, 100 ),
NU_swc_2 = c( 0, 100 ),
NU_swc_3 = c( 0, 100 ),
NU_swc_4 = c( 0, 100 ),
g_tgrowth_1 =c( -50, 50 ),
g_tgrowth_2 =c( 0, 50 ),
g_tgrowth_3 =c( 0, 50 ),
g_tgrowth_4 =c( 0, 50 ),
g_swc_1 = c( 0, 100 ),
g_swc_2 = c( 0, 100 ),
r_tloss_1 = c( -50, 50 ),
r_tloss_2 = c( 0, 50 ),
L_tloss_1 = c( -50, 50 ),
L_tloss_2 = c( 0, 50 ),
L_swc_1 = c( 0, 100 ),
L_swc_2 = c( 0, 100 ),
phen_t_0 = c( 0, 1 ),
phen_t_1 = c( 0, 365 ),
phen_t_2 = c( 0, 365 ),
phen_t_3 = c( 0, 365 ),
phen_t_4 = c( 0, 365 ),
L_fire_1 = c( 0, 100 ),
L_fire_2 = c( 0, 100 ),
L_mix_1 = c( 0, 100 ),
L_mix_2 = c( 0, 100 ),
L_mix_3 = c( 0, 100 ),
iM = c( 0, 100 ),
pe_Ms = c( 0, 100 ),
pe_Mr = c( 0, 100 ),
pe_Cs = c( 0, 100 ),
pe_Cr = c( 0, 100 ),
pe_Ns = c( 0, 100 ),
pe_Nr = c( 0, 100 )
)
)
#' Standard options for the model
#' @description This list contains the standard options used by \code{make_data()}.
#' It is provided here as an example. See \code{options} for more information.
#' @export
standard_options <- list(
optim = "DEoptim",
ve = "red",
#output_steps
steps = c(50,100),
initial_mass = 0.5,
initial_mass_par = FALSE, # set initial mass as estimated parameter TRUE or FALSE
pe = FALSE, # consider process error TRUE or FALSE
species = "species_1",
photo = "c3", #for multi-species a vector of length SPECIES
pe_scale = 1e-04
)
#' Options to be used in the \code{make_data()} function
#'
#' @name options
#' @title Option list used in the \code{make_data()} function
#'
#' @param optim
#' String specifying the data calculated by the user defined Model function, e.g.
#' one of \code{c("DEoptim", "LaplacesDemon")}.
#'
#' @param ve
#' String specifying the version of the process model. Must be one of \code{c("std", "fqr", "red", "oak")}.
#'
#' @param steps
#' Integer numeric vector specifying the timesteps to output. The process model
#' will run \code{max(steps)} iterations and the output will contain
#' \code{length(steps)} output times.
#'
#' @param initial_mass_par
#' Boolean specifying if the initial biomass used in the process model should be
#' estimated as a parameter.
#'
#' @param initial_mass
#' The initial biomass to use for the process model, when \code{initial_mass_par} is set
#' to \code{FALSE}.
#'
#' @param species
#' A vector of names of species.
#'
#' @param photo
#' A vector specifying the photosynthesis system used by the species.
#' Must be either \code{c("c3", "c4")} for each species.
#'
#' @param swc_online
#' Boolean specifying if soil water content should be calculated online as part of the
#' process model (\code{TRUE}) or if it is given as a forcing parameter. Currently not supported.
#'
#' @param lc
#' Boolean turning light competition simulation in the process model on/off. Currently not supported.
#'
#'
#' @param fire
#' Boolean turning fire in the process model on/off. Currently not supported.
#'
#' @param pe
#' Boolean turning process error simulation on/off.
#'
#' @param pe_scale
#' The process error scaling factor.
#'
#' @description
#' The options list handed to \code{make_data()}. It specifies the process model options.
#' An example is provided by \code{standard_options}.
NULL
#' Create data object for use with the process model
#' @description This function creates an object representing a configuration of the process model including forcing data.
#' This object is used by \code{run_ttr()}.
#' @param input The forcing data object created by a call to \code{get_input()}, see specification in its help page
#' @param options An options list for running the model, see \code{options} help page
#' @param globals Optional parameter specifying global variables used in the simulation, see \code{standard_globals} help page
#' @param bounds Optional parameter specifying the upper and lower bounds of the model's parameters, see \code{make_data()}. Bounds for additional parameters included in the statistical model need to be specified here. Bounds of the process model parameters can also be modified from their default values.
#' @param groups Optional parameter specifying the groups of process model parameters
#' @param verbose Print out a log detailing checks
#' @return A ttr data object, a list with the following elements:
#' \item{PGF}{Function used by LaplacesDemon for generating initial values.
#' }
#' \item{options}{Options as specified in the \code{options} and \code{standard_options} help pages
#' }
#' \item{globals}{Globals as specified in the \code{standard_globals} help page
#' }
#' \item{n.sites}{Number of sites
#' }
#' \item{n.species}{Number of species
#' }
#' \item{n.time}{Number of timesteps
#' }
#' \item{n.parm}{Number of parameters to be fitted
#' }
#' \item{n.parm.a}{Number of per-process parameters
#' }
#' \item{n.parm.b}{Number of per-species parameters
#' }
#' \item{parm.names}{Names of all parameters for this model configuration
#' }
#' \item{parm.names.a}{Names of per-process parameter names
#' }
#' \item{parm.names.b}{Names of per-species parameter names
#' }
#' \item{mon.names}{Required by LaplacesDemon, always set to "LP"
#' }
#' \item{out}{Names of the response variables returned by \code{run_ttr()}
#' }
#' \item{out.dimnames}{Dimension names for the output returned by \code{run_ttr()}
#' }
#' \item{out.dim}{Dimensions for the output returned by \code{run_ttr()}
#' }
#' \item{dimnames.beta}{Dimension names for the ttr parameters object produced by \code{get_parms()}
#' }
#' \item{bounds}{Bounds for all process model parameters
#' }
#' \item{timeseries}{See \code{get_input()}
#' }
#' \item{timeinvariant}{See \code{get_input()}
#' }
#' \item{y}{observations from \code{get_input()} cast to a matrix
#' }
#' \item{lonlat}{a matrix with lon and lat columns for the nsites rows
#' }
#' \item{pos.oe}{Indices of observation error parameters
#' }
#' \item{pos.pe}{Indices of process error parameters
#' }
#' \item{parm.trap.groups}{Representation of parameters that form a group for the trapezoid functions
#' }
#' \item{arguments}{The unprocessed arguments to this function, excluding data
#' }
#' @details The data structure produced in this function can be handed to an external parameter
#' estimation algorithm e.g. LaplacesDemon or DEoptim to fit the model.
#'
#' The bounds of the process model's parameters need to be defined.
#' The user can supply these to \code{make_data()} or accept the default
#' values returned by a call to \code{make_data()}. Different model variants use
#' different parameters and a call to \code{make_data()} returns the default bounds for the
#' model variant specified in \code{options}.
#'
#' bounds is a list of two lists; one for alpha (per-process) parameters and
#' one for beta (per-species) parameters. The list bounds is formatted as:
#' bounds = list(alpha = list(alpha_parameter = (lower, upper)),
#' beta = list(beta_parameter = c(lower, upper))), where `upper` specifies
#' the maximum value a process model parameter can take on, and
#' `lower` is the corresponding minimum.
#'
#' The external parameter estimation algorithm should then take Data$bounds[,1] as
#' lower and Data$bounds[,2] as the upper bounds.
#'
#' In case the user does not wish to estimate some of the parameters in the
#' external optimisation process, the upper and lower bound of a parameter
#' should be set to the same numeric value. In this case, the parameter will
#' not be part of the Data$bounds matrix passed to the external parameter estimation
#' algorithm, and its value is set internally in the \code{get_parms()} function.
#'
#' For any group of parameters, e.g., CU_Ns_1 and CU_Ns_2, setting any one
#' of the parameters to c(NA, NA) will set the whole group of parameters to
#' a constant value that will negate the influence of this parameter group on the
#' process model. This effectively switches off the effect of these parameters
#' on the model's processes and removes them from Data$bounds. The constant
#' values the parameters are set to can be observed by inspecting the output
#' of \code{get_parms()} for a Data object.
#'
#' The names of the parameters are written as YY_XX_i, where YY specifies the
#' process being considered and XX the environmental driver and i the parameter
#' number in the step or trapezoid function, e.g. CU_swc_1 and CU_swc_2 specify
#' how CU (carbon uptake) is influenced by swc (soil water content) and 1 and 2
#' specify the lower and upper parameters of the step function.
#'
#' An example for specifying all process bounds for a Model run with the
#' standard options is given in the object `standard_bounds`.
#'
#' @importFrom stats rnorm
#' @examples
#' #some dummy data as input
#' input_data <- get_input(
#' observations = c(1),
#' tcur = c(1),
#' tnur = c(1),
#' tgrowth = c(1),
#' tloss = c(1),
#' seconds = c(1),
#' lon = c(1),
#' lat = c(1),
#' rsds = c(1),
#' catm = c(1),
#' pet = c(1),
#' rain = c(1),
#' wp = c(1),
#' fc = c(1)
#' )
#' data <- make_data(
#' input = input_data,
#' options = standard_options,
#' globals = standard_globals,
#' bounds = list(
#' alpha = list(
#' #defining a new parameter, e.g. as part of a model function
#' my_parameter = c(1,2),
#' #overwriting a standard bound
#' scale = c(0,500)
#' ),
#' #no more parameters than needed for the process model in beta, for these standard values.
#' beta = list()
#' )
#' )
#' #all standard values
#' data <- make_data(input_data, standard_options)
#' @export
make_data <- function(input, options=standard_options, globals=standard_globals, bounds=list(alpha = NULL, beta=NULL), groups=c(), verbose=FALSE){
arguments <- list(options=options, globals=globals, bounds=bounds, groups=groups)
ts <- input$timeseries
time.invariant <- input$time.invariant
observations <- input$observations
in_data_dimensions <- input$dimensions
#options
assertList(options)
assertChoice(options[["optim"]], c("DEoptim", "LaplacesDemon", "no"))
assertIntegerish(
options[["steps"]],
any.missing=FALSE,
min.len=1,
lower=0
)
options[["steps"]] <- unique(sort(options[["steps"]]))
assertChoice(options[["ve"]],varnames("TTRVariant"))
assertNumber(
options[["initial_mass"]],
lower=0.0,
finite=T,
)
for(boolprop in option_switches){
assertFlag(options[[boolprop]], .var.name=paste0("options$",boolprop))
}
assertCharacter(
options[["species"]],
min.len = 1,
any.missing = F,
unique = T
)
if(!testFactor(options[["photo"]],
levels=varnames("PhotoType"))){
assertCharacter(
options[["photo"]],
len = length(options[["species"]]),
any.missing = F
)
for(photo_type in options[["photo"]]){
assertChoice(
photo_type,
choices = varnames("PhotoType")
)
}
options$photo <- factor(options$photo, levels = varnames("PhotoType"))
}
assertFactor(
options[["photo"]],
levels=varnames("PhotoType"),
len=length(options[["species"]]),
any.missing=FALSE
)
assertNumber(
options[["pe_scale"]],
finite=T
)
assertList(
bounds
)
assertNames(
names(bounds),
permutation.of = c("alpha","beta")
)
if(!is.null(bounds$alpha)){
assertList(bounds$alpha)
}
if(!is.null(bounds$beta)){
assertList(bounds$beta)
}
if(length(bounds$alpha) > 0 && length(bounds$beta) > 0){
assertNames(
names(bounds$alpha),
disjunct.from = names(bounds$beta)
)
}
req <- required_parameters(options)
missing_alpha_ind <- which(!(req$alpha %in% names(bounds$alpha)))
missing_beta_ind <- which(!(req$beta %in% names(bounds$beta)))
if(length(missing_alpha_ind) > 0 || length(missing_beta_ind) > 0){
if(verbose){
cat(paste0("Missing parameter bounds for required parameters, using standard bounds for: ",
paste(req$beta[missing_beta_ind], collapse = ", "), ", ",
paste(req$alpha[missing_alpha_ind], collapse = ", "), "\n"))
}
}
for(missing_beta in req$beta[missing_beta_ind]){
bounds$beta[missing_beta] <- standard_bounds$beta[missing_beta]
}
for(missing_alpha in req$alpha[missing_alpha_ind]){
bounds$alpha[missing_alpha] <- standard_bounds$alpha[missing_alpha]
}
bounds$beta <- bounds$beta[sort(names(bounds$beta))]
bounds$alpha <- bounds$alpha[sort(names(bounds$alpha))]
#parameter groups
if(!is.null(groups)){
assertCharacter(
groups,
any.missing=FALSE
)
}
groups <- unique(c(groups, standard_groups))
#globals
assertNumeric(
globals,
any.missing=FALSE,
names="unique"
)
assertNames(
names(globals),
type="unique",
permutation.of=varnames("Global")
)
globals <- globals[varnames("Global")]
#Optimisation: q has to be 1, a different value requires recompilation
assertNumber(
globals[["q"]],
lower = 1.0,
upper = 1.0
)
#time series
assertArray(
ts,
mode="double",
any.missing=FALSE,
d=3,
.var.name="time dependent forcing variables"
)
assertNames(
dimnames(ts)[[1]],
type="unique",
identical.to=varnames("TimeSeries"),
.var.name="forcing variable names"
)
#time invariant
assertArray(
time.invariant,
mode="double",
any.missing=FALSE,
d=2,
.var.name="time-independent forcing variables"
)
assertNames(
dimnames(time.invariant)[[1]],
type="unique",
identical.to=varnames("TimeInvariant")
)
if(dim(ts)[[3]] != in_data_dimensions[["nsites"]]){
stop(paste0("Time series array has an incorrect number of sites: ",dim(ts)[[3]], " instead of ", in_data_dimensions[["nsites"]]))
}
if(dim(time.invariant)[[2]] != in_data_dimensions[["nsites"]]){
stop(paste0("Time independent variable array has an incorrect number of sites: ",dim(time.invariant)[[2]], " instead of ", in_data_dimensions[["nsites"]]))
}
#cat("\nchecks complete\n")
#From here, assume all bounds are specifying a parameter in the model
pn <- names(bounds)
nspecies <- length(options$species)
lower_of <- function(bounds){
vapply(bounds, function(entry){entry[[1]]}, FUN.VALUE=c(0.0))
}
upper_of <- function(bounds){
vapply(bounds, function(entry){entry[[2]]}, FUN.VALUE=c(0.0))
}
#replace NAs in bounds with constant values
lower.a <- lower_of(bounds$alpha)
upper.a <- upper_of(bounds$alpha)
for(group in groups){
lower.a <- replace.na(lower.a, group)
upper.a <- replace.na(upper.a, group)
}
lower.b.single <- lower_of(bounds$beta)
upper.b.single <- upper_of(bounds$beta)
for(group in groups){
lower.b.single <- replace.na(lower.b.single, group)
upper.b.single <- replace.na(upper.b.single, group)
}
lower.b <- rep(lower.b.single, nspecies)
upper.b <- rep(upper.b.single, nspecies)
lower <- c(lower.a, lower.b)
upper <- c(upper.a, upper.b)
#remove constant parameters from bounds structure
ind.parm.constant <- which(lower == upper)
parm.constant <- c()
if(length(ind.parm.constant) > 0){
parm.constant <- lower[ind.parm.constant]
lower <- lower[-ind.parm.constant]
upper <- upper[-ind.parm.constant]
}
out.dimnames <- list(
varnames("Biotic"),
options$species,
as.character(options$steps),
as.character(1:in_data_dimensions[["nsites"]])
)
out.dim <- c(
length(varnames("Biotic")),
length(options$species),
length(options$steps),
in_data_dimensions[["nsites"]]
)
names(out.dim) <- c(
"output",
"species",
"steps",
"sites"
)
PGF <- function(Data){
lo <- Data$bounds[,1]
up <- Data$bounds[,2]
mu <- (lo+up)/2
pmin(pmax(rnorm(Data$n.parm, mu, 0.1), lo), up)
}
opt.bounds <- matrix(data=c(lower,upper),ncol=2,dimnames=list(names(lower),c("lower","upper")))
n.constant <- length(ind.parm.constant)
n.constant.a <- length(which(ind.parm.constant <= length(bounds$alpha)))
n.constant.b <- length(which(ind.parm.constant > length(bounds$alpha)))
parm.names <- names(lower)
data <- structure(
list(
PGF=PGF,
options=options,
globals=globals,
n.sites = in_data_dimensions[["nsites"]],
n.species = length(options$species),
n.time = in_data_dimensions[["nsteps"]],
n.parm = dim(opt.bounds)[1],
n.parm.proc = length(bounds$alpha) + nspecies*length(bounds$beta),
n.parm.a = length(bounds$alpha) - n.constant.a,
n.parm.proc.a = length(bounds$alpha),
n.parm.b = length(bounds$beta) - n.constant.b,
n.parm.proc.b = length(bounds$beta),
ind.parm.constant = ind.parm.constant,
parm.constant = parm.constant,
out.dimnames = out.dimnames,
out.dim = out.dim,
parm.names.a=names(bounds$alpha),
parm.names.b=names(bounds$beta),
parm.names=parm.names,
dimnames.beta = list(names(bounds$beta),options$species),
parm.trap.groups = parameter_groups(bounds, groups),
pos.oe = grep("^oe", bounds$alpha),
pos.pe = grep("^pe", bounds$beta),
mon.names = "LP",
out=varnames("Biotic"),
bounds=opt.bounds,
timeseries=ts,
time.invariant=time.invariant,
lonlat=input$lonlat,
y=observations,
arguments=arguments),
class = "ttr.data")
return(data)
}
Try the TTR.PGM package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.