Defines functions do_acirun acifun_wrap do_fit_method_bilinear_bestcitrans do_fit_method_bilinear do_fit_method1 set_Rdmeas set_Tleaf set_PPFD guess_Rd guess_Vcmax guess_Jmax rmse_acifit fitaci

Documented in fitaci

#' Fit the Farquhar-Berry-von Caemmerer model of leaf photosynthesis
#' @description Fits the Farquhar-Berry-von Caemmerer model of photosynthesis to measurements of 
#' photosynthesis and intercellular \eqn{CO_2}{CO2} concentration (Ci). Estimates Jmax, Vcmax, Rd 
#' and their standard errors. A simple plotting method is also included, as well as the function 
#' \code{\link{fitacis}} which quickly fits multiple A-Ci curves (see its help page). Temperature 
#' dependencies of the parameters are taken into account following Medlyn et al. (2002), see 
#' \code{\link{Photosyn}} for more details.
#' @param data Dataframe with Ci, Photo, Tleaf, PPFD (the last two are optional). For \code{fitacis}, 
#' also requires a grouping variable.
#' @param varnames List of names of variables in the dataset (see Details).
#' @param Tcorrect If TRUE, Vcmax and Jmax are corrected to 25C. Otherwise, Vcmax and Jmax 
#' are estimated at measurement temperature. \strong{Warning} : since package version 1.4, the default parameters have been adjusted (see Details).
#' @param Patm Atmospheric pressure (kPa)
#' @param citransition If provided, fits the Vcmax and Jmax limited regions separately (see Details).
#' @param quiet If TRUE, no messages are written to the screen.
#' @param startValgrid If TRUE (the default), uses a fine grid of starting values to increase 
#' the chance of finding a solution.
#' @param fitmethod Method to fit the A-Ci curve. Either 'default' (Duursma 2015), 'bilinear' (See Details), or 'onepoint' (De Kauwe et al. 2016).
#' @param algorithm Passed to \code{\link{nls}}, sets the algorithm for finding parameter values.
#' @param fitTPU Logical (default FALSE). Attempt to fit TPU limitation (fitmethod set to 'bilinear' 
#' automatically if used). See Details.
#' @param alphag When estimating TPU limitation (with \code{fitTPU}), an 
#' additional parameter (see Details).
#' @param useRd If Rd provided in data, and useRd=TRUE (default is FALSE), uses measured Rd in 
#' fit. Otherwise it is estimated from the fit to the A-Ci curve.
#' @param PPFD Photosynthetic photon flux density ('PAR') (mu mol m-2 s-1)
#' @param Tleaf Leaf temperature (degrees C)
#' @param alpha Quantum yield of electron transport (mol mol-1)
#' @param theta Shape of light response curve.
#' @param gmeso Mesophyll conductance (mol m-2 s-1 bar-1). If not NULL (the default), Vcmax 
#' and Jmax are chloroplastic rates.
#' @param EaV,EdVC,delsC Vcmax temperature response parameters
#' @param EaJ,EdVJ,delsJ Jmax temperature response parameters
#' @param Km,GammaStar Optionally, provide Michaelis-Menten coefficient for Farquhar model, and 
#' Gammastar. If not provided, they are calculated with a built-in function of leaf temperature.
#' @param id Names of variables (quoted, can be a vector) in the original dataset to be stored
#'  in the result. Most useful when using \code{\link{fitacis}}, see there for examples of its use.
#' @param \dots Further arguments (ignored at the moment).
#' @details 
#' \subsection{Fitting method}{The default method to fit A-Ci curves (set by 
#' \code{fitmethod="default"}) uses non-linear regression to fit the A-Ci curve. No assumptions 
#' are made on which part of the curve is Vcmax or Jmax limited. Normally, all three parameters 
#' are estimated: Jmax, Vcmax and Rd, unless Rd is provided as measured (when \code{useRd=TRUE}, 
#' and Rd is contained in the data). This is the method as described by Duursma (2015, Plos One).
#' The 'bilinear' method to fit A-Ci curves (set by \code{fitmethod="bilinear"}) linearizes 
#' the Vcmax and Jmax-limited regions, and applies linear regression twice to estimate first 
#' Vcmax and Rd, and then Jmax (using Rd estimated from the Vcmax-limited region). The transition 
#' point is found as the one which gives the best overall fit to the data (i.e. all possible 
#' transitions are tried out, similar to Gu et al. 2010, PCE). The advantage of this method is that 
#' it \emph{always} returns parameter estimates, so it should be used in cases where the default 
#' method fails. Be aware, though, that the default method fails mostly when the curve contains 
#' bad data (so check your data before believing the fitted parameters).
#' When \code{citransition} is set, it splits the data into a Vcmax-limited (where Ci < citransition), 
#' and Jmax-limited region (Ci > citransition). Both parameters are then estimated separately for 
#' each region (Rd is estimated only for the Vcmax-limited region). \bold{Note} that the actual 
#' transition point as shown in the standard plot of the fitted A-Ci curve may be quite different 
#' from that provided, since the fitting method simply decides which part of the dataset to use 
#' for which limitation, it does not constrain the actual estimated transition point directly. 
#' See the example below. If \code{fitmethod="default"}, it applies non-linear regression to 
#' both parts of the data, and when fitmethod="bilinear", it uses linear regression on the 
#' linearized photosynthesis rate. Results will differ between the two methods (slightly).}
#' The 'onepoint' fitting method is a very simple estimation of Vcmax and Jmax for every point 
#' in the dataset, simply by inverting the photosynthesis equation. See De Kauwe et al. (2016) 
#' for details. The output will give the original data with Vcmax and Jmax added (note you can 
#' set \code{Tcorrect} as usual!). For increased reliability, this method only works if 
#' dark respiration (Rd) is included in the data (\code{useRd} is set automatically when 
#' setting \code{fitmethod='one-point'}). This method is not recommended for full A-Ci curves, 
#' but rather for spot gas exchange measurements, when a simple estimate of Vcmax or Jmax 
#' is needed, for example when separating stomatal and non-stomatal drought effects on 
#' photosynthesis (Zhou et al. 2013, AgForMet). The user will have to decide whether the Vcmax 
#' or Jmax rates are used in further analyses. This fitting method can not be used in \code{fitacis}, 
#' because Vcmax and Jmax are already estimated for each point in the dataset.
#' \subsection{TPU limitation}{Optionally, the \code{fitaci} function estimates the triose-phosphate 
#' utilization (TPU) rate. The TPU can act as another limitation on photosynthesis, and can be 
#' recognized by a 'flattening out' of the A-Ci curve at high Ci. When \code{fitTPU=TRUE}, the 
#' fitting method used will always be 'bilinear'. The TPU is estimated by trying out whether the 
#' fit improves when the last n points of the curve are TPU-limited (where n=1,2,...). When TPU is 
#' estimated, it is possible (though rare) that no points are Jmax-limited (in which case estimated 
#' Jmax will be NA). A minimum of two points is always reserved for the estimate of Vcmax and Rd. 
#' An additional parameter (\code{alphag}) can be set that affects the behaviour at high Ci (see 
#' Ellsworth et al. 2015 for details, and also \code{\link{Photosyn}}). See examples.}
#' \subsection{Temperature correction}{When \code{Tcorrect=TRUE} (the default), Jmax and Vcmax 
#' are re-scaled to 25C, using the temperature response parameters provided (but Rd is always 
#' at measurement temperature). When \code{Tcorrect=FALSE}, estimates of all parameters are at 
#' measurement temperature. If TPU is fit, it is never corrected for temperature. Important 
#' parameters to the fit are GammaStar and Km, both of which are calculated from leaf temperature 
#' using standard formulations. Alternatively, they can be provided as known inputs. \strong{Warning} : 
#' since package version 1.4, the default parameters have been adjusted. The new parameter values 
#' (EaV, EdVJ, delSJ, etc.) were based on a comprehensive literature review. See 
#' \code{vignette("new_T_responses")} or the article on remkoduursma.github.io/plantecophys.}
#' \subsection{Mesophyll conductance}{It is possible to provide an estimate of the mesophyll 
#' conductance as input (\code{gmeso}), in which case the fitted Vcmax and Jmax are to be interpreted 
#' as chloroplastic rates. When using gmeso, it is recommended to use the 'default' fitting 
#' method (which will use the Ethier&Livingston equations inside \code{Photosyn}). It is also 
#' implemented with the 'bilinear' method but it requires more testing (and seems to give 
#' some strange results). When gmeso is set to a relatively low value, the resulting fit 
#' may be quite strange.}
#' \subsection{Other parameters}{The A-Ci curve parameters depend on the values of a number 
#' of other parameters. For Jmax, PPFD is needed in order to express it as the asymptote. If 
#' PPFD is not provided in the dataset, it is assumed to equal 1800 mu mol m-2 s-1 (in which 
#' case a warning is printed). It is possible to either provide PPFD as a variable in the 
#' dataset (with the default name 'PARi', which can be changed), or as an argument to 
#' the \code{fitaci} directly.}
#' \subsection{Plotting and summarizing}{The default \strong{plot} of the fit is constructed 
#' with \code{\link{plot.acifit}}, see Examples below. When plotting the fit, the A-Ci curve 
#' is simulated using the \code{\link{Aci}} function, with leaf temperature (Tleaf) and PPFD 
#' set to the mean value for the dataset. The \strong{coefficients} estimated in the fit 
#' (Vcmax, Jmax, and usually Rd) are extracted with \code{coef}. The summary of the fit is 
#' the same as the 'print' method, that is \code{myfit} will give the same output as 
#' \code{summary(myfit)} (where \code{myfit} is an object returned by \code{fitaci}).
#' Because fitaci returns the fitted \code{\link{nls}} object, more details on statistics of 
#' the fit can be extracted with standard tools. The Examples below shows the use of the 
#' \pkg{nlstools} to extract many details of the fit at once. The fit also includes the 
#' \strong{root mean squared error} (RMSE), which can be extracted as \code{myfit$RMSE}. This 
#' is a useful metric to compare the different fitting methods.}
#' \subsection{Predicting and the CO2 compensation point}{The fitted object contains two functions 
#' that reproduce the fitted curve exactly. Suppose your object is called 'myfit', then 
#' \code{myfit$Photosyn(200)} will give the fitted rate of photosynthesis at a Ci of 200. 
#' The inverse, calculating the Ci where some rate of photosynthesis is achieved, can be done with
#' \code{myfit$Ci(10)} (find the Ci where net photosynthesis is ten). The (fitted!) CO2 
#' compensation point can then be calculated with : \code{myfit$Ci(0)}}.
#' \subsection{Atmospheric pressure correction}{Note that atmospheric pressure (Patm) is taken 
#' into account, assuming the original data are in molar units (Ci in mu mol mol-1, or ppm). 
#' During the fit, Ci is converted to mu bar, and Km and Gammastar are recalculated accounting 
#' for the effects of Patm on the partial pressure of oxygen. When plotting the fit, though,
#'  molar units are shown on the X-axis. Thus, you should get (nearly) the same fitted curve when
#'   Patm was set to a value lower than 100kPa, but the fitted Vcmax and Jmax will be higher. 
#'   This is because at low Patm, photosynthetic capacity has to be higher to achieve the same 
#'   measured photosynthesis rate.}
#' @troubleshooting From time to time, the \code{fitaci} function returns an error, indicating 
#' that the A-Ci curve could not be fit. In the majority of cases, this indicates a bad curve 
#' that probably could not (and should not) be fit in any case. Inspect the raw data to check 
#' if the curve does not include severe outliers, or large deviations from the expected, typical 
#' A-Ci curve. If the curve looks fine, refit using the option \code{fitmethod="bilinear"}, 
#' which will always return estimated parameters. Whatever you do, do not trust fitted curves 
#' without inspecting the raw data, and the fit of the model to the data.
#' @references
#' Duursma, R.A., 2015. Plantecophys - An R Package for Analysing and Modelling Leaf Gas 
#' Exchange Data. PLoS ONE 10, e0143346. doi:10.1371/journal.pone.0143346
#' De Kauwe, M. G. et al. 2016. A test of the 'one-point method' for estimating maximum 
#' carboxylation capacity from field-measured, light-saturated photosynthesis. 
#' New Phytol 210, 1130-1144.
#' @return A list of class 'acifit', with the following components:
#' \describe{
#' \item{df}{A dataframe with the original data, including the measured photosynthetic 
#' rate (Ameas), the fitted photosynthetic rate (Amodel), Jmax and Vcmax-limited gross 
#' rates (Aj, Ac), TPU-limited rate (Ap), dark respiration (Rd), leaf temperature (Tleaf), 
#' chloroplastic CO2 (Cc), PPFD, atmospheric pressure (Patm), and 'original Ci, i.e. the 
#' Ci used as input (which is different from the Ci used in fitting if Patm was not set to 100kPa)}
#' \item{pars}{Contains the parameter estimates and their approximate standard errors}
#' \item{nlsfit}{The object returned by \code{\link{nls}}, and contains more detail on 
#' the quality of the fit}
#' \item{Tcorrect}{whether the temperature correction was applied (logical)}
#' \item{Photosyn}{A copy of the \code{\link{Photosyn}} function with the arguments adjusted for 
#' the current fit. That is, Vcmax, Jmax and Rd are set to those estimated in the fit, and Tleaf and 
#' PPFD are set to the mean value in the dataset. All other parameters that were set in fitaci are 
#' also used (e.g. temperature dependency parameters, TPU, etc.).}
#' \item{Ci}{As Photosyn, except the opposite: calculate the Ci where some rate of net photosynthesis 
#' is achieved.}
#' \item{Ci_transition}{The Ci at which photosynthesis transitions from Vcmax 
#' to Jmax limited photosynthesis.}
#' \item{Ci_transition2}{The Ci at which photosynthesis transitions from Jmax to
#' TPU limitation. Set to NA is either TPU was not estimated, or it could not be 
#' estimated from the data.}
#' \item{Rd_measured}{Logical - was Rd provided as measured input?}
#' \item{GammaStar}{The value for GammaStar, either calculated or provided to the fit.}
#' \item{Km}{he value for Km, either calculated or provided to the fit.}
#' \item{kminput}{Was Km provided as input? (If FALSE, it was calculated from Tleaf)}
#' \item{gstarinput}{Was GammaStar provided as input? (If FALSE, it was calculated from Tleaf)}
#' \item{fitmethod}{The fitmethod uses, either default or bilinear}
#' \item{citransition}{The input citransition (NA if it was not provided as input)}
#' \item{gmeso}{The mesophyll conductance used in the fit (NA if it was not set)}
#' \item{fitTPU}{Was TPU fit?}
#' \item{alphag}{The value of alphag used in estimating TPU.}
#' \item{RMSE}{The Root-mean squared error, calculated as \code{sqrt(sum((Ameas-Amodel)^2))}.}
#' \item{runorder}{The data returned in the 'df' slot are ordered by Ci, but in rare cases the 
#' original order of the data contains information; 'runorder' is 
#' the order in which the data were provided.}
#' }
#' @examples
#' \dontrun{
#' # Fit an A-Ci curve on a dataframe that contains Ci, Photo and optionally Tleaf and PPFD. 
#' # Here, we use the built-in example dataset 'acidata1'.
#' f <- fitaci(acidata1)
#' # Note that the default behaviour is to correct Vcmax and Jmax for temperature, 
#' # so the estimated values are at 25C. To turn this off:
#' f2 <- fitaci(acidata1, Tcorrect=FALSE)
#' # To use different T response parameters (see ?Photosyn),
#' f3 <- fitaci(acidata1, Tcorrect=TRUE, EaV=25000)
#' # Make a standard plot
#' plot(f)
#' # Look at a summary of the fit
#' summary(f)
#' # Extract coefficients only
#' coef(f)
#' # The object 'f' also contains the original data with predictions.
#' # Here, Amodel are the modelled (fitted) values, Ameas are the measured values.
#' with(f$df, plot(Amodel, Ameas))
#' abline(0,1)
#' # The fitted values can also be extracted with the fitted() function:
#' fitted(f)
#' # The non-linear regression (nls) fit is stored as well,
#' summary(f$nlsfit)
#' # Many more details can be extracted with the nlstools package
#' library(nlstools)
#' overview(f$nlsfit)
#' # The curve generator is stored as f$Photosyn:
#' # Calculate photosynthesis at some value for Ci, using estimated 
#' # parameters and mean Tleaf, PPFD for the dataset.
#' f$Photosyn(Ci=820)
#' # Photosynthetic rate at the transition point:
#' f$Photosyn(Ci=f$Ci_transition)$ALEAF
#' # Set the transition point; this will fit Vcmax and Jmax separately. Note that the *actual* 
#' # transition is quite different from that provided, this is perfectly fine : 
#' # in this case Jmax is estimated from the latter 3 points only (Ci>800), but the actual 
#' # transition point is at ca. 400ppm.
#' g <- fitaci(acidata1, citransition=800)
#' plot(g)
#' g$Ci_transition
#' # Use measured Rd instead of estimating it from the A-Ci curve. 
#' # The Rd measurement must be added to the dataset used in fitting, 
#' # and you must set useRd=TRUE.
#' acidata1$Rd <- 2
#' f2 <- fitaci(acidata1, useRd=TRUE)
#' f2
#' # Fit TPU limitation
#' ftpu <- fitaci(acidata1, fitTPU=TRUE, PPFD=1800, Tcorrect=TRUE)
#' plot(ftpu)
#' }
#' @export
#' @rdname fitaci
#' @importFrom stats lm
fitaci <- function(data, 
                   varnames=list(ALEAF="Photo", Tleaf="Tleaf", 
                                 Ci="Ci", PPFD="PARi", Rd="Rd"),
                   quiet=FALSE, startValgrid=TRUE, 
                   EaV = 82620.87,
                   EdVC = 0,
                   delsC = 645.1013,
                   EaJ = 39676.89,
                   EdVJ = 200000,
                   delsJ = 641.3615,
                   GammaStar = NULL,
                   Km = NULL,
  fitmethod <- match.arg(fitmethod)
  if(fitTPU & fitmethod == "default")fitmethod <- "bilinear"
  if(fitTPU & fitmethod == "onepoint"){
    Stop("TPU limitation not implemented with one-point method.")
  if(fitmethod == "onepoint")useRd <- TRUE
  gstarinput <- !is.null(GammaStar)
  kminput <- !is.null(Km)
  if(is.null(gmeso))gmeso <- -999  # cannot pass NULL value to nls
  # Check data 
  if(nrow(data) == 0){
    Stop("No rows in data - check observations.")
  # Make sure data is a proper dataframe, not some tibble or other nonsense.
  data <- as.data.frame(data)
  # Extra parameters not used at the moment; warn when provided.
  # Add PPFD and Tleaf to data, if needed (uses default values, or input values)
  data <- set_PPFD(varnames, data, PPFD, quiet)
  data <- set_Tleaf(varnames, data, Tleaf, quiet)
  # Set measured Rd if provided (or warn when provided but not used)
  Rd_meas <- set_Rdmeas(varnames, data, useRd, citransition, quiet)
  haveRd <- !is.na(Rd_meas)
  # Must have Rd with one-point method.
  if(!haveRd && fitmethod == "onepoint"){
    Stop("Must add Rd to the dataset (and check that varnames is set correctly).")
  # Extract Ci and apply pressure correction
  data$Ci_original <- data[,varnames$Ci]
  data$Ci <- data[,varnames$Ci] * Patm/100
  # Extract measured net leaf photosynthesis
  data$ALEAF <- data[,varnames$ALEAF]
  # Calculate Km and GammaStar, if not input 
  # (Photosyn does this as well, but have to repeat it here 
  # since we cannot have NULL in call to nls)
  Km_v <- if(!kminput) TKm(data$Tleaf, Patm) else rep(Km, nrow(data))
  GammaStar_v <- if(!gstarinput){
    TGammaStar(data$Tleaf, Patm)
  } else {
    rep(GammaStar, nrow(data))
  # Citransition not defined:
  # default method = full non-linear model
  # bilinear = optimize transition point in bilinear fit
    if(fitmethod == "default"){
      f <- do_fit_method1(data, haveRd, Rd_meas, Patm, startValgrid, 
                          Tcorrect, algorithm,alpha,theta,gmeso,EaV,
    if(fitmethod == "bilinear"){
      f <- do_fit_method_bilinear_bestcitrans(data, haveRd, fitTPU, alphag, Rd_meas, 
                                              Patm, Tcorrect, algorithm,
                                              GammaStar_v, Km_v)
    if(fitmethod == "onepoint"){
      f <- do_fit_method_bilinear(data, haveRd, alphag, Rd_meas, Patm, 
                                  NA,NA,  # don't matter
                                  Tcorrect=Tcorrect, algorithm,
                                  GammaStar, Km, onepoint=TRUE)
  # Ci transition defined. 
  # default - two nonlinear fits (first Vcmax, then Jmax region)
  # bilinear - two linear fits
    # NOTE-- bug in default method when citransition specified. Switch off for now.
    fitmethod <- "bilinear"
    # if(fitmethod == "default"){
    #   f <- do_fit_method2(data, haveRd, Rd_meas, Patm, citransition, Tcorrect, 
    #                       algorithm,alpha,theta,gmeso,EaV,EdVC,
    #                       delsC,EaJ,EdVJ,delsJ,GammaStar_v,Km_v)
    # }
    if(fitmethod == "bilinear"){
      f <- do_fit_method_bilinear(data, haveRd, alphag, Rd_meas, Patm, citransition, 
                                  NULL, Tcorrect, algorithm,
                                GammaStar_v, Km_v)
  # TPU. If it is not estimated, put something in because we need it later.
  if(!("TPU" %in% names(f)))f$TPU <- 1000
  # Only used to add 'Amodel' to the output
  acirun <- do_acirun(data,f,Patm,Tcorrect,  # Note several parameters are stored in 'f'
                      delsJ=delsJ,GammaStar=GammaStar_v, Km=Km_v)

  # If Ap is never actually limiting, set estimated TPU to 1000
  # This means there is no evidence for TPU-limitation.
  if(!any(acirun$Ap < acirun$Aj))f$TPU <- 1000
  # Organize output
  l <- list()  
  runorder <- order(acirun$Ci)
  l$df <- acirun[runorder,]
    if(any(!id %in% names(data))){
      Warning("id ignored - not all variables are in dataset provided")
    } else {
      l$df <- cbind(l$df, data[id])
  l$id <- id
  l$pars <- f$pars
    val <- ifelse(f$TPU < 1000, f$TPU, NA)
    tpu <- matrix(c(val,NA), ncol=2)
    colnames(tpu) <- colnames(l$pars)
    rownames(tpu) <- "TPU"
    l$pars <- rbind(l$pars, tpu)
  l$nlsfit <- f$fit
  l$Tcorrect <- Tcorrect
  # Save function itself, the formals contain the parameters 
  # used to fit the A-Ci curve.
  # First save Tleaf, PPFD in the formals (as the mean of the dataset)
  formals(Photosyn)$Tleaf <- mean(data$Tleaf)
  formals(Photosyn)$Patm <- Patm
  formals(Photosyn)$PPFD <- mean(data$PPFD)
  formals(Photosyn)$Vcmax <- l$pars[1]
  formals(Photosyn)$Jmax <- l$pars[2]
  formals(Photosyn)$Rd <- l$pars[3]
  formals(Photosyn)$TPU <- if(fitTPU)l$pars["TPU","Estimate"] else 1000
  formals(Photosyn)$alphag <- alphag
  formals(Photosyn)$Tcorrect <- Tcorrect
  formals(Photosyn)$alpha <- alpha
  formals(Photosyn)$theta <- theta
  formals(Photosyn)$gmeso <- gmeso
  formals(Photosyn)$EaV <- EaV
  formals(Photosyn)$EdVC <- EdVC
  formals(Photosyn)$delsC <- delsC
  formals(Photosyn)$EaJ <- EaJ
  formals(Photosyn)$EdVJ <- EdVJ
  formals(Photosyn)$delsJ <- delsJ
  if(gstarinput)formals(Photosyn)$GammaStar <- GammaStar
  if(kminput)formals(Photosyn)$Km <- Km
  l$Photosyn <- Photosyn
  # The inverse - find Ci given a rate of photosyntheiss
  l$Ci <- function(ALEAF){
    O <- function(ci, photo){
      (l$Photosyn(Ci=ci)$ALEAF - photo)^2
    o <- optimize(O, interval=c(1,10^5), photo =ALEAF)
    if(o$objective > 1e-02){
    } else {
  # Transition points.
  trans <- findCiTransition(l$Photosyn)
  l$Ci_transition <- trans[1]
  l$Ci_transition2 <- trans[2]
  l$Rd_measured <- haveRd
  # Save GammaStar and Km (either evaluated at mean temperature, 
  # or input if provided)
  l$GammaStar <- mean(GammaStar_v)
  l$Km <- mean(Km_v)
  l$kminput <- kminput
  l$gstarinput <- gstarinput
  l$fitmethod <- fitmethod
  l$citransition <- ifelse(is.null(citransition), NA, citransition)
  l$gmeso <- ifelse(is.null(gmeso) || gmeso < 0, NA, gmeso)
  l$fitTPU <- fitTPU
  l$alphag <- alphag
  l$RMSE <- rmse_acifit(l)
  l$runorder <- runorder
  class(l) <- "acifit"

#-- Component functions of fitaci

rmse_acifit <- function(x)sqrt(sum((x$df$Ameas - x$df$Amodel)^2))

guess_Jmax <- function(data, Rd_guess, Patm, Tcorrect){
  # Guess Jmax from max A, T-corrected gammastar (starting value)
  maxCi <- max(data$Ci)
  mi <- which.max(data$Ci)

  maxPhoto <- data$ALEAF[mi]

  Tl <- data$Tleaf[mi]
  gammastar <- TGammaStar(Tl,Patm)
  VJ <- (maxPhoto+Rd_guess) / ((maxCi - gammastar) / (maxCi + 2*gammastar))
  Jmax_guess <- VJ*4
    Teffect <- TJmax(Tl,  EaJ=39676.89, delsJ=641.3615, EdVJ=200000)
    Jmax_guess <- Jmax_guess / Teffect

guess_Vcmax <- function(data, Jmax_guess, Rd_guess, Patm, Tcorrect){
  # Guess Vcmax, from section of curve that is definitely Vcmax-limited
  dato <- data[data$Ci < 150 & data$Ci > 60 & data$ALEAF > 0,]
  if(nrow(dato) > 0){
    Km <- TKm(dato$Tleaf,Patm)
    gammastar <- TGammaStar(dato$Tleaf,Patm)
    vcmax <- with(dato, (ALEAF + Rd_guess) / ((Ci - gammastar)/(Ci + Km)))
    Vcmax_guess <- median(vcmax)
  } else {
    Vcmax_guess <- Jmax_guess/1.8 
    Teffect <- TVcmax(mean(data$Tleaf), EaV=82620.87, delsC=645.1013, EdVC=0)
    Vcmax_guess <- Vcmax_guess / Teffect

guess_Rd <- function(haveRd, Rd_meas, default=1.5){
    Rd_guess <- Rd_meas
  } else {
    Rd_guess <- default

set_PPFD <- function(varnames, data, PPFD, quiet){
  if(!varnames$PPFD %in% names(data) & is.null(PPFD)){
    data$PPFD <- 1800
    if(!quiet)Warning("PARi not in dataset; assumed PARi = 1800.")
  } else {
      data$PPFD <- PPFD
      data$PPFD <- data[,varnames$PPFD]

set_Tleaf <- function(varnames, data, Tleaf, quiet){
  # Check if Tleaf is provided
  if(!varnames$Tleaf %in% names(data) & is.null(Tleaf)){
    data$Tleaf <- 25
    if(!quiet)Warning("Tleaf not in dataset; assumed Tleaf = 25.")
  } else {
      data$Tleaf <- Tleaf
      data$Tleaf <- data[,varnames$Tleaf]

# Check if Rd is provided and whether we want to set it as known.
set_Rdmeas <- function(varnames, data, useRd, citransition, quiet){
  Rd_meas <- NA
    if(varnames$Rd %in% names(data) && useRd){
      # Has to be a single unique value for this dataset
      Rd_meas <- data[,varnames$Rd]
      Rd_meas <- unique(Rd_meas)
      if(length(Rd_meas) > 1){
        Stop("If Rd provided as measured, it must be a single",
             "\nunique value for an A-Ci curve.")
      # Use positive value throughout.
      Rd_meas <- abs(Rd_meas)
      haveRd <- TRUE
        Stop("At the moment cannot provide citransition as well as measured Rd.")
    if(varnames$Rd %in% names(data) && !useRd){
        message(paste("Rd found in dataset but useRd set to FALSE.",
                      "Set to TRUE to use measured Rd."))

do_fit_method1 <- function(data, haveRd, Rd_meas, Patm, startValgrid, 
                           Tcorrect, algorithm,
  # Guess Rd (starting value)
  Rd_guess <- guess_Rd(haveRd, Rd_meas)
  Jmax_guess <- guess_Jmax(data, Rd_guess, Patm, Tcorrect)
  Vcmax_guess <- guess_Vcmax(data, Jmax_guess, Rd_guess, Patm, Tcorrect)
  # Fine-tune starting values; try grid of values around initial estimates.
      aciSS <- function(Vcmax, Jmax, Rd){
        Photo_mod <- acifun_wrap(data$Ci, PPFD=data$PPFD, 
                                 Vcmax=Vcmax, Jmax=Jmax, 
                                 Rd=Rd, Tleaf=data$Tleaf, Patm=Patm,
        SS <- sum((data$ALEAF - Photo_mod)^2)
      d <- 0.4
      n <- 25
      gg <- expand.grid(Vcmax=seq(Vcmax_guess*(1-d),Vcmax_guess*(1+d),length=n),
      m <- with(gg, mapply(aciSS, Vcmax=Vcmax, Jmax=Jmax_guess, Rd=Rd))
      ii <- which.min(m)
      Vcmax_guess <- gg$Vcmax[ii]
      Rd_guess <- gg$Rd[ii]
    } else {
      aciSS <- function(Vcmax, Jmax, Rd){
        Photo_mod <- acifun_wrap(data$Ci, PPFD=data$PPFD, 
                                 Vcmax=Vcmax, Jmax=Jmax, 
                                 Rd=Rd, Tleaf=data$Tleaf, Patm=Patm,
        SS <- sum((data$ALEAF - Photo_mod)^2)
      d <- 0.3
      n <- 20
      gg <- data.frame(Vcmax=seq(Vcmax_guess*(1-d),Vcmax_guess*(1+d),length=n))
      m <- with(gg, mapply(aciSS, Vcmax=Vcmax, Jmax=Jmax_guess, Rd=Rd_meas))
      ii <- which.min(m)
      Vcmax_guess <- gg$Vcmax[ii]
    # Fit Vcmax, Jmax and Rd
    nlsfit <- try(nls(ALEAF ~ acifun_wrap(Ci, PPFD=PPFD, Vcmax=Vcmax, 
                                      Jmax=Jmax, Rd=Rd, Tleaf=Tleaf, Patm=Patm,
                  data=data, control=nls.control(maxiter=500, minFactor=1/10000),
                  start=list(Vcmax=Vcmax_guess, Jmax=Jmax_guess, Rd=Rd_guess)), silent=TRUE)
    if(inherits(nlsfit, "try-error")){
      Stop("Could not fit curve - check quality of data or fit using fitmethod='bilinear'.")
    p <- coef(nlsfit)
    pars <- summary(nlsfit)$coefficients[,1:2]
  } else {
    # Fit Vcmax and Jmax
    nlsfit <- try(nls(ALEAF ~ acifun_wrap(Ci, PPFD=PPFD, Vcmax=Vcmax, 
                                      Jmax=Jmax, Rd=Rd_meas, Tleaf=Tleaf, Patm=Patm,
                  data=data, control=nls.control(maxiter=500, minFactor=1/10000),
                  start=list(Vcmax=Vcmax_guess, Jmax=Jmax_guess)), silent=TRUE)
    if(inherits(nlsfit, "try-error")){
      Stop("Could not fit curve -",
           "\ncheck quality of data or fit using fitmethod='bilinear'.")
    p <- coef(nlsfit)
    p[[3]] <- Rd_meas
    names(p)[3] <- "Rd"
    pars <- summary(nlsfit)$coefficients[,1:2]
    pars <- rbind(pars, c(Rd_meas, NA))
    rownames(pars)[3] <- "Rd"
  return(list(pars=pars, fit=nlsfit))

# do_fit_method2 <- function(data, haveRd, Rd_meas, Patm, citransition, 
#                            Tcorrect,algorithm,alpha,theta,gmeso,EaV,EdVC,
#                            delsC,EaJ,EdVJ,delsJ,GammaStar,Km){
#   # Guess Rd (starting value)
#   Rd_guess <- guess_Rd(haveRd, Rd_meas)
#   Jmax_guess <- guess_Jmax(data, Rd_guess, Patm, Tcorrect)
#   Vcmax_guess <- guess_Vcmax(data, Jmax_guess, Rd_guess, Patm, Tcorrect)
#   # If citransition provided, fit twice.
#   dat_vcmax <- data[data$Ci < citransition,]
#   dat_jmax <- data[data$Ci >= citransition,]
#   if(nrow(dat_vcmax) > 0){
#     nlsfit_vcmax <- nls(ALEAF ~ acifun_wrap(Ci, PPFD=PPFD, Vcmax=Vcmax, 
#                                             Jmax=10^6, Rd=Rd, Tleaf=mean(Tleaf), 
#                                             Patm=Patm, returnwhat="Ac",
#                                             TcorrectVJ=Tcorrect,
#                                             alpha=alpha,theta=theta,
#                                             gmeso=gmeso,EaV=EaV,
#                                             EdVC=EdVC,delsC=delsC,
#                                             EaJ=EaJ,EdVJ=EdVJ,
#                                             delsJ=delsJ,Km=Km,GammaStar=GammaStar),
#                         algorithm=algorithm,
#                         data=dat_vcmax, 
#                         control=nls.control(maxiter=500, minFactor=1/10000),
#                         start=list(Vcmax=Vcmax_guess, Rd=Rd_guess))
#     p1 <- coef(nlsfit_vcmax)
#   } else {
#     nlsfit_vcmax <- NULL
#     p1 <- c(Vcmax=NA, Rd=NA)
#   }
#   # Don't re-estimate Rd; it is better estimated from the Vcmax-limited region.
#   Rd_vcmaxguess <- if(!is.null(nlsfit_vcmax))coef(nlsfit_vcmax)[["Rd"]] else 1.5
#   if(nrow(dat_jmax) > 0){
#     nlsfit_jmax <- nls(ALEAF ~ acifun_wrap(Ci, PPFD=PPFD, Vcmax=10000, 
#                                            Jmax=Jmax, Rd=Rd_vcmaxguess, 
#                                            Tleaf=mean(Tleaf), Patm=Patm, returnwhat="Aj",
#                                            TcorrectVJ=Tcorrect,
#                                            alpha=alpha,theta=theta,
#                                            gmeso=gmeso,EaV=EaV,
#                                            EdVC=EdVC,delsC=delsC,
#                                            EaJ=EaJ,EdVJ=EdVJ,
#                                            delsJ=delsJ,Km=Km,GammaStar=GammaStar),
#                        algorithm=algorithm,
#                        data=dat_jmax, control=nls.control(maxiter=500, minFactor=1/10000),
#                        start=list(Jmax=Jmax_guess))
#     p2 <- coef(nlsfit_jmax)
#   } else { 
#     nlsfit_jmax <- NULL
#     p2 <- c(Jmax=NA)
#   }
#   p <- c(p1[1],p2,p1[2])
#   pars1 <- if(!is.null(nlsfit_vcmax)){
#     summary(nlsfit_vcmax)$coefficients[,1:2] 
#   } else {
#     matrix(rep(NA,4),ncol=2)
#   }
#   pars2 <- if(!is.null(nlsfit_jmax)){
#     summary(nlsfit_jmax)$coefficients[,1:2] 
#   } else {
#     matrix(rep(NA,2),ncol=2)
#   }
#   pars <- rbind(pars1[1,],pars2,pars1[2,])
#   rownames(pars) <- c("Vcmax","Jmax","Rd")
#   nlsfit <- list(nlsfit_vcmax=nlsfit_vcmax,nlsfit_jmax=nlsfit_jmax)
#   return(list(pars=pars, fit=nlsfit))      
# }

do_fit_method_bilinear <- function(data, haveRd, alphag, Rd_meas, 
                                   Patm, citransition, citransition2=NULL,
                                   Tcorrect, algorithm,
                                   GammaStar, Km, onepoint=FALSE){
  # Calculate T-dependent parameters
  ppar <- Photosyn(Tleaf=data$Tleaf, Patm=Patm, Tcorrect=Tcorrect,
  if(!is.null(GammaStar))ppar$GammaStar <- GammaStar
  if(!is.null(Km))ppar$Km <- Km
  if(is.null(citransition2))citransition2 <- 10^6
  # Calculate Cc if gmeso included
  if(!is.null(gmeso) && gmeso > 0){
    Conc <- data$Ci - data$ALEAF/gmeso
  } else {
    Conc <- data$Ci
  # Linearize
  data$vcmax_pred <- (Conc - ppar$GammaStar)/(Conc + ppar$Km)
  data$Jmax_pred <- (Conc - ppar$GammaStar)/(Conc + 2*ppar$GammaStar)
  data$TPU_part <- (Conc - ppar$GammaStar)/(Conc - (1+3*alphag)*ppar$GammaStar)
    orig_dat <- data[, 1:(match("Ci_original", names(data))-1)]
    orig_dat$Vcmax <- (data$ALEAF + Rd_meas)  / data$vcmax_pred
    J <- (data$ALEAF + Rd_meas)  / data$Jmax_pred
    orig_dat$Jmax <- inverseJfun(mean(data$PPFD), alpha, 4 * J, theta)
      orig_dat$Jmax <- orig_dat$Jmax / TJmax(data$Tleaf, EaJ, delsJ, EdVJ)
      orig_dat$Vcmax <- orig_dat$Vcmax / TVcmax(data$Tleaf,EaV, delsC, EdVC)
  # Fit Vcmax and Rd from linearized portion
  datv <- data[data$Ci < citransition & data$Ci < citransition2,]
  if(nrow(datv) == 0){
    return(list(pars=NA, fit=NA, TPU=NA, success=FALSE))
    fitv <- lm(ALEAF ~ vcmax_pred, data=datv)
    Rd_fit <- coef(fitv)[[1]]
    Vcmax_fit <- coef(fitv)[[2]]
  } else {
    # If using measured Rd, add to Anet, and remove intercept from fit.
    datv$ALEAFg <- datv$ALEAF + Rd_meas
    fitv <- lm(ALEAFg ~ vcmax_pred-1, data=datv)
    Rd_fit <- -Rd_meas
    Vcmax_fit <- coef(fitv)[[1]]
  # Fit Jmax from linearized portion
  datj <- data[data$Ci >= citransition & data$Ci < citransition2,]
  datp <- data[data$Ci >= citransition2,]

  # Manual fix: if only one point for TPU, and none for Jmax, abandon fit.
  # In this case it would be more defensible to use the single point for Jmax.
  if(nrow(datp) == 1 && nrow(datj) == 0){
    return(list(pars=NA, fit=NA, TPU=NA, success=FALSE))
  if(nrow(datj) > 0){
    # Fit gross photo using fitted Rd
    # (Rd_fit is negative)
    datj$Agross <- datj$ALEAF - Rd_fit
    # One point, calculate directly
    if(nrow(datj) == 1){
      J_fit <- with(datj, 4 * Agross / Jmax_pred)
    } else {
      fitj <- lm(Agross ~ Jmax_pred-1, data=datj)
      J_fit <- 4 * coef(fitj)[[1]]
    # And solve for Jmax from inverse non-rect. hyperbola
    Jmax_fit <- inverseJfun(mean(data$PPFD), alpha, J_fit, theta)
  } else {
    Jmax_fit <- 10^6  # not elegant but will do for now 
                      # (avoids trouble elsewhere)
  # TPU
  if(nrow(datp) > 0){
    datp$Agross <- datp$ALEAF - Rd_fit
    tpu_vals <- (1/3) * datp$Agross / datp$TPU_part
    TPU <- mean(tpu_vals)
  } else {
    TPU <- 1000  # same as default in Photosyn
  # The above estimates are at the measured Tleaf.
  # Express at 25C?
    Jmax_fit <- Jmax_fit / TJmax(mean(data$Tleaf), EaJ, delsJ, EdVJ)
    Vcmax_fit <- Vcmax_fit / TVcmax(mean(data$Tleaf),EaV, delsC, EdVC)
  ses <- summary(fitv)$coefficients[,2]
    pars <- matrix(c(Vcmax_fit, Jmax_fit, -Rd_fit,
                   ses[2],NA,ses[1]), ncol=2)
  } else {
    pars <- matrix(c(Vcmax_fit, Jmax_fit, -Rd_fit,
                     ses[1],NA,NA), ncol=2)
  rownames(pars) <- c("Vcmax","Jmax","Rd")
  colnames(pars) <- c("Estimate","Std. Error")
  return(list(pars=pars, fit=fitv, TPU=TPU, success=TRUE))

do_fit_method_bilinear_bestcitrans <- function(data, haveRd, fitTPU, alphag, 
                                               Rd_meas, Patm, Tcorrect,
                                               delsJ,GammaStar, Km){
  # Possible Ci transitions
  ci <- data$Ci
  nci <- length(ci)
  citransitions <- diff(ci)/2 + ci[-nci]
  # at least two Ci values to estimate Vcmax and Rd, so delete first
  citransitions1 <- citransitions[-1]
  # Possible transitions to TPU
    citransitions2 <- max(ci) + 1  # outside range, on purpose
  } else {
    # start at top, all the way down, leave lowest 2 points alone
    citransitions2 <- c(max(ci) + 1, rev(citransitions1))
  citransdf <- expand.grid(ci1=citransitions1, ci2=citransitions2)
  citransdf <- citransdf[citransdf$ci1 <= citransdf$ci2,]
  SS <- c()
  # Note that Tcorrect is set to FALSE inside the loop. If Tcorrect is needed, it is done
  # after the loop finishes (avoids a bug).
  for(i in seq_len(nrow(citransdf))){
    fit <- do_fit_method_bilinear(data, haveRd, alphag, Rd_meas, Patm, 
                                  citransdf$ci1[i], citransdf$ci2[i], 
                                  Tcorrect=FALSE, algorithm,
                                  GammaStar, Km)
    if(fit$success && !any(is.na(fit$pars[,"Estimate"]))){
        run <- do_acirun(data,fit,Patm,Tcorrect=FALSE,
        SS[i] <- sum((run$Ameas - run$Amodel)^2)  
    } else {
      SS[i] <- 10^6
  # Best Ci transitions
  bestcis <- citransdf[which.min(SS),]
  f <- do_fit_method_bilinear(data, haveRd, alphag, Rd_meas, Patm, 
                              bestcis$ci1, bestcis$ci2, Tcorrect, algorithm,
                              GammaStar, Km)
  if(f$pars["Jmax","Estimate"] < 0){
    Stop("Cannot invert light response curve to estimate Jmax - increase alpha or theta.")


# Wrapper around Photosyn; this wrapper will be sent to nls. 
acifun_wrap <- function(Ci,..., TcorrectVJ, returnwhat="ALEAF"){
  r <- Photosyn(Ci=Ci,Tcorrect=TcorrectVJ,...)
  if(returnwhat == "ALEAF")return(r$ALEAF)
  if(returnwhat == "Ac")return(r$Ac - r$Rd)
  if(returnwhat == "Aj")return(r$Aj - r$Rd)

# Using fitted coefficients, get predictions from model.
do_acirun <- function(data,f,Patm,Tcorrect,...){
  acirun <- Photosyn(Ci=data$Ci, 
                     Vcmax=f$pars[1,1], Jmax=f$pars[2,1], 
  acirun$Ameas <- data$ALEAF
  acirun$ELEAF <- NULL
  acirun$GS <- NULL
  acirun$Ca <- NULL
  acirun$Ci_original <- data$Ci_original
  names(acirun)[names(acirun) == "ALEAF"] <- "Amodel"
  # shuffle
  avars <- match(c("Ci","Ameas","Amodel"),names(acirun))
  acirun <- acirun[,c(avars, setdiff(1:ncol(acirun), avars))]

Try the plantecophys package in your browser

Any scripts or data that you put into this service are public.

plantecophys documentation built on May 2, 2019, 6:36 a.m.