R/run_ttr.R

In TTR.PGM: Thornley Transport Resistance Plant Growth Model

#' interval_parms: truncate parameters into bounds
#' @description This function takes a numeric vector and constrains it to specified bounds by reflecting values outside 
#' the bounds into the interval. The bounds are defined by a data object returned by \code{make_data()}
#' @param par A numeric vector of parameters as supplied by LaplacesDemon or DEoptim
#' @param Data A data object as returned by \code{make_data()}
#' @return A numeric vector in which each parameter is reflected into bounds.
#' @export
interval_parms <- function(par, Data){
  interval_c(par,Data$bounds)
}

#' get_parms: Pack parameters in an easily accessible list
#' @description This function takes a numeric vector of all process model parameters and returns a 
#' representation that is human readable and can be used by \code{run_ttr()}
#' @param par A numeric vector of parameters as supplied by an external parameter estimation algorithm, e.g. LaplacesDemon or DEoptim
#' @param Data data object as returned by the \code{make_data()} function
#' @param no.structure If 'TRUE', only run \code{interval_parms()} and \code{unname()} the result, this is convenient when using LaplacesDemon.
#' @return A list containing a vector of the alpha parameters and an array of dimension (nbeta, nspecies)
#'         where nbeta is the number of beta parameters and nspecies is the number of species, containing the
#'         beta parameters per species
#' @export
get_parms <- function(par, Data, no.structure=FALSE){
  #bound parameters
  par <- interval_parms(par, Data)

  #Add constant parameters back into the structure
  if(length(Data$ind.parm.constant) > 0){
    ppar <- rep(0,Data$n.parm.proc)
    ppar[Data$ind.parm.constant] <- Data$parm.constant
    ppar[-Data$ind.parm.constant] <- par
    par <- ppar
  }

  par <- unname(par)
  if(no.structure){
    return(par)
  }
  #split alpha/beta parameters
  nalpha <- Data$n.parm.proc.a
  nbeta <- Data$n.parm.proc.b
  nspecies <- Data$n.species
  if(nalpha > 0){
    alpha_pars <- par[1:nalpha]
    beta_pars <- par[-c(1:nalpha)]
    names(alpha_pars) <- Data$parm.names.a
  }
  else{
    alpha_pars <- c()
    beta_pars <- par
  }
  #format parameter output
  attr(beta_pars, "dim") <- c(nbeta, nspecies)
  attr(beta_pars, "dimnames") <- Data$dimnames.beta
  #sum up trapezoid parameters
  groups <- Data$parm.trap.groups
  for(g in groups$alpha){
    partial <- 0
    for(ind in g){
      partial <- partial + alpha_pars[[ind]]
      alpha_pars[[ind]] <- partial
    }
  }
  for(g in groups$beta){
    partial <- rep(0,nspecies)
    for(ind in g){
      partial <- partial + beta_pars[ind,]
      beta_pars[ind,] <- partial
    }
  }
  return(list(alpha=alpha_pars, beta=beta_pars))
}

#' run_ttr: Simulate the TTR model
#'
#' @param parm The parameters as returned by the \code{get_parms()} function
#' @param data The model object as defined in the \code{make_data()} function, see its help page for details
#' @return A four-dimensional matrix-object with the dimensions (output_var,species,time,site)
#' where:
#' \item{output_var}{refers to the array of output values produced by the process model. 
#'   For identifiers, check \code{Data$out} as produced by \code{make_data()} }
#' \item{species}{is the nth species}
#' \item{time}{is the nth output time given in model$options$steps}
#' \item{site}{is the nth site}
#' @export
run_ttr <- function(parm, data){
  if(testNumeric(parm)){
    parm <- get_parms(parm, data)
  }
  res <- run_ttr_cpp(parm, data$timeseries, data$time.invariant, data$options,data$globals)
  dim(res) <- data$out.dim
  dimnames(res) <- data$out.dimnames
  res
}