Nothing
#' R5 reference parent class for describing the NEP ~ PAR relationship
#'
#' @import methods
#' @export LightResponseCurveFitter
#' @exportClass LightResponseCurveFitter
LightResponseCurveFitter <- setRefClass('LightResponseCurveFitter')
#' @export
LightResponseCurveFitter_getParameterNames <- function(
### return the parameter names used by this Light Response Curve Function
) {
# if this method is adjusted in subclass
# , then also adjust getPriorLocation, getPriorScale
##value<< string vector of parameter names. Positions are important.
c(k = "k" ##<< VPD effect
, beta = "beta" ##<< saturation of GPP at high radiation
, alpha = "alpha" ##<< initial slope
, RRef = "RRef" ##<< basal respiration (units of provided NEE,
## usually mumol CO2 m-^-2 s^-2)
, E0 = "E0") ##<< temperature sensitivity estimated
## from night-time data (K)
}
LightResponseCurveFitter$methods(getParameterNames =
LightResponseCurveFitter_getParameterNames)
#' @export
LightResponseCurveFitter_fitLRC <- function(
### Optimize rectangular hyperbolic light response curve in one window
dsDay ##<< data.frame with columns NEE, Rg, Temp_C, VPD, and no NAs in NEE
, E0 ##<< temperature sensitivity of respiration
, sdE0 ##<< standard deviation of E_0.n
, RRefNight ##<< basal respiration estimated from night time data
, controlGLPart = partGLControl() ##<< further default parameters
## (see \code{\link{partGLControl}})
, lastGoodParameters = rep(NA_real_, 7L) ##<< numeric vector returned
## by last reasonable fit
) {
##author<< TW, MM
##seealso<< \code{partGLFitLRCWindows}
#
# Definition of initial guess theta, theta2 and theta3.
#Three initial guess vectors are defined according to Lasslop et al., 2010
parNames <- .self$getParameterNames() # hook method from derived classes
nPar <- length(parNames)
##details<<
## Optimization is performed for three initial parameter sets that differ
## by \code{beta0 (* 1.3, * 0.8)}.
## From those three, the optimization result is selected that yielded
## the lowest misfit.
## Starting values are: \code{k = 0},
## \code{beta = interpercentileRange(0.03, 0.97)} of
## respiration, \code{alpha = 0.1}, \code{R_ref}
## from nightTime estimate.
## \code{E0} is fixed to the night-time estimate, but varies for estimating
## parameter uncertainty.
parameterPrior <- .self$getPriorLocation(dsDay$NEE, RRefNight = RRefNight
, E0 = E0)
thetaInitials <- .self$getParameterInitials(parameterPrior)
#
##seealso<< \code{\link{LightResponseCurveFitter_optimLRCBounds}}
resOpt3 <- apply(thetaInitials, 1, function(theta0) {
resOpt <- .self$optimLRCBounds(theta0, parameterPrior, dsDay = dsDay
, ctrl = controlGLPart
, lastGoodParameters = lastGoodParameters)
})
iValid <- which(sapply(resOpt3, function(resOpt) { is.finite(resOpt$theta[1]) }))
resOpt3Valid <- resOpt3[iValid]
optSSE <- sapply(resOpt3Valid, "[[", "value")
getNAResult <- function(convergenceCode) {
list(
thetaOpt = structure(rep(NA_real_, nPar), names = colnames(thetaInitials) )
, iOpt = integer(0) ##<< index of parameters that have been optimized
, thetaInitialGuess = ##<< the initial guess from data
thetaInitials[1, ]
, covParms = matrix(NA_real_, nPar, nPar, dimnames = list(
colnames(thetaInitials), colnames(thetaInitials)))
, convergence = convergenceCode
)}
if (sum(!is.na(optSSE)) == 0L) {
return(getNAResult(resOpt3[[1]]$convergence) )
} else {
resOpt <- resOpt3Valid[[iBest <- which.min(optSSE)]] # select lowest cost
thetaOpt<- resOpt$theta
if (controlGLPart$nBootUncertainty == 0L) {
##details<< If \code{controlGLPart$nBootUncertainty == 0L} then the
## covariance matrix of the
## parameters is estimated by the Hessian of the LRC curve at optimum.
## Then, the additional uncertainty and covariance with uncertainty E0
## is neglected.
#seParmsHess <- seParmsHess0 <- sqrt(abs(diag(solve(resOpt$hessian))))
covParmsLRC <- try(if ( (resOpt$hessian[1L, 1L] < 1e-8) ) {
# case where k = 0 and not varying: cov(k, :) = cov(:, ) = 0
covParmsLRC <- structure(diag(0, nrow = nrow(resOpt$hessian))
, dimnames = dimnames(resOpt$hessian))
covParmsLRC[-1L, -1L] <- solve(resOpt$hessian[-1L, -1L])
covParmsLRC
} else {
solve(resOpt$hessian)
}, silent = TRUE)
if (inherits(covParmsLRC, "try-error")) {
return(getNAResult(1006L)) # count not invert the Hessian
}
covParms <- structure(diag(0, nrow = length(resOpt$theta))
, dimnames = list(parNames, parNames))
covParms[5L, 5L] <- sdE0^2
covParms[resOpt$iOpt, resOpt$iOpt] <- covParmsLRC
if (any(diag(covParms) < 0)) {
#ggplot(dsDay, aes(Rg, -NEE, color = Temp)) + geom_point()
return(getNAResult(1005L) )
}
} else {
##details<<
## If \code{controlGLPart.l$nBootUncertainty > 0L} then the
## covariance matrix of the
## parameters is estimated by a bootstrap of the data.
## In each draw, E0 is drawn from N ~ (E_0, sdE_0).
# #seealso<< \code{\link{.bootStrapLRCFit}}
resBoot <- .bootStrapLRCFit(resOpt$theta, resOpt$iOpt, dsDay, sdE0
, parameterPrior, controlGLPart, lrcFitter = .self)
iFiniteRows <- which(is.finite(resBoot[, 1L]))
##details<< If there are no estimates for more than 20% of the
## bootstrapped samples
## The an NA-result with convergence code 1001L is returned.
if (length(iFiniteRows) < 0.8 * nrow(resBoot))
return(getNAResult(1001L) )
covParms <- cov(resBoot[iFiniteRows, ])
if (covParms[2, 2] < 1e-8) getNAResult(1004L)
}
}
# further parameter checking after parameter uncertainty has been computed
# (before here, other parameter checking is done in optimLRCBounds)
sdTheta <- thetaOpt; sdTheta[] <- NA
sdTheta[resOpt$iOpt] <- sqrt(diag(covParms)[resOpt$iOpt])
if (!.self$isParameterInBounds(thetaOpt, sdTheta, RRefNight = RRefNight
, ctrl = controlGLPart) )
return(getNAResult(1002L))
# debugging: tracing specific LRC fits:
##value<< a list, If none of the optimizations from different starting
## conditions converged, the parameters are NA.
ans <- list(
thetaOpt = thetaOpt ##<< numeric vector of optimized parameters
## including the fixed ones and E0
, iOpt = c(resOpt$iOpt, 5L) ##<< index of parameters that have been
## optimized, here including E0, which has been optimized prior to
## this function.
, thetaInitialGuess = ##<< the initial guess from data
thetaInitials[1, ]
, covParms = covParms ##<< numeric matrix of the covariance matrix
## of parameters, including E0
, convergence = resOpt$convergence ##<< integer code specifying
## convergence problems: \\
## 0: good convergence \\
## , 1-1000: see \code{\link{optim}} \\
## , 1001: too few bootstraps converged \\
## , 1002: fitted parameters were outside reasonable bounds \\
## , 1003: too few valid records in window \\
## , 1004: near zero covariance in bootstrap indicating bad fit \\
## , 1005: covariance from curvature of fit yielded negative variances
## indicating bad fit \\
## , 1006: prediction of highest PAR in window was far from saturation
## indicating insufficient data to constrain LRC \\
## , 1010: no temperature-respiration relationship found \\
## , 1011: too few valid records in window
## (from different location: partGLFitLRCOneWindow) \\
)
}
LightResponseCurveFitter$methods(fitLRC = LightResponseCurveFitter_fitLRC)
#' @export
LightResponseCurveFitter_getPriorLocation <- function(
### return the prior distribution of parameters
NEEDay ##<< numeric vector of daytime NEE
, RRefNight ##<< numeric scalar of basal respiration estimated
## from night-time data
, E0 ##<< numeric scalar of night-time estimate of temperature sensitivity
) {
##value<< a numeric vector with prior estimates of the parameters
parameterPrior <- c(
k = 0.05
#, beta = as.vector(abs(quantile(NEEDay, 0.03)-quantile(NEEDay, 0.97)))
, beta = as.vector(abs(quantile(NEEDay, 0.03, na.rm = TRUE) -
quantile(NEEDay, 0.97, na.rm = TRUE)))
, alpha = 0.1
#, RRef = mean(NEENight.V.n, na.rm = T)
, RRef = if (is.finite(RRefNight) ) as.vector(RRefNight) else stop(
"must provide finite RRefNight") #mean(NEENight.V.n, na.rm = T)
, E0 = as.vector(E0)
)
}
LightResponseCurveFitter$methods(getPriorLocation =
LightResponseCurveFitter_getPriorLocation)
#' @export
LightResponseCurveFitter_getPriorScale <- function(
### return the prior distribution of parameters
thetaPrior ##<< numeric vector of location of priors
, medianRelFluxUncertainty ##<< numeric scalar: median across the
## relative uncertainty of the flux values, i.e. sdNEE / NEE
, nRec ##<< integer scalar: number of finite observations
, ctrl ##<< list of further controls, with entry
## \code{isLasslopPriorsApplied}
) {
##value<< a numeric vector with prior estimates of the parameters
sdParameterPrior <- if (ctrl$isLasslopPriorsApplied) {
#twutz: changed to no prior for logitconv
c(k = 50, beta = 600, alpha = 10, RRef = 80, E0 = NA)
} else {
##details<< The beta parameter is quite well defined. Hence use a prior with
##a standard deviation. The specific results are sometimes a bit sensitive
##to the uncertainty of the beta prior. This uncertainty is set
##corresponding to 20 times the median relative flux uncertainty. The prior
##is weighted n times the observations in the cost. Hence, overall it is
##using a weight of 1 / 20 of the weight of all observations.
##
##However, its not well defined if PAR does not reach saturation. Need to
##check before applying this prior
sdBetaPrior <- 20 * medianRelFluxUncertainty * thetaPrior[2] / sqrt(nRec)
c(k = NA, beta = as.vector(sdBetaPrior), alpha = NA, RRef = NA, E0 = NA)
}
}
LightResponseCurveFitter$methods(getPriorScale = LightResponseCurveFitter_getPriorScale)
#' @export
LightResponseCurveFitter_getParameterInitials <- function(
### return the prior distribution of parameters
thetaPrior ##<< numeric vector prior estimate of parameters
) {
# only beta (second parameter) is varied between different initial guesses
##value<< a numeric matrix (3, nPar) of initial values for fitting parameters
thetaInitials <- matrix(rep(thetaPrior, each = 3), 3, length(thetaPrior)
, dimnames = list(NULL, names(thetaPrior)))
thetaInitials [2, 2] <- thetaPrior[2] * 1.3
thetaInitials [3, 2] <- thetaPrior[2] * 0.8
thetaInitials
}
LightResponseCurveFitter$methods(getParameterInitials =
LightResponseCurveFitter_getParameterInitials)
#' @export
LightResponseCurveFitter_optimLRCBounds <- function(
### Optimize parameters with refitting with some fixed parameters if outside bounds
theta0 ##<< initial parameter estimate
, parameterPrior ##<< prior estimate of model parameters
, ... ##<< further parameters to \code{.optimLRC},
, dsDay ##<< argument to \code{.optimLRC}, here checked
## for occurrence of high VPD
, lastGoodParameters ##<< parameters vector of last successful fit
, ctrl ##<< list of further controls, such as
## \code{isNeglectVPDEffect = TRUE}
) {
##author<< TW, MM
##seealso<< \code{\link{LightResponseCurveFitter_fitLRC}}
# twutz 161014: default alpha
if (!is.finite(lastGoodParameters[3L]) ) lastGoodParameters[3L] <- 0.22
isNeglectVPDEffect <- isTRUE(ctrl$isNeglectVPDEffect)
VPD0 = 10 ##<< VPD0 [hPa] -> Parameters VPD0 fixed to 10 hPa
## according to Lasslop et al 2010
isUsingFixedVPD <- isNeglectVPDEffect ||
(sum(dsDay$VPD >= VPD0, na.rm = TRUE) == 0)
isUsingFixedAlpha <- FALSE
getIOpt <- .self$getOptimizedParameterPositions
theta0Adj <- theta0 # initial estimate with some parameters adjusted to bounds
if (isNeglectVPDEffect) theta0Adj[1] <- 0
resOpt <- resOpt0 <- .self$optimLRCOnAdjustedPrior(theta0Adj, iOpt =
getIOpt(isUsingFixedVPD, isUsingFixedAlpha)
, parameterPrior = parameterPrior, ctrl, dsDay = dsDay, ...)
##details<<
## If parameters alpha or k are outside bounds (Table A1 in Lasslop 2010),
## refit with some parameters fixed
## to values from fit of previous window.
#dsDay <- list(ctrl, ...)$dsDay
if (is.na(resOpt$theta[1L]) || (resOpt$theta[1L] < 0) ) { #k
isUsingFixedVPD <- TRUE
theta0Adj[1L] <- 0
resOpt <- .self$optimLRCOnAdjustedPrior(theta0Adj, iOpt = getIOpt(
isUsingFixedVPD, isUsingFixedAlpha)
, parameterPrior = parameterPrior, ctrl, dsDay = dsDay, ...)
# check alpha, in case refit with fixed alpha of last window
if ( (is.na(resOpt$theta[3L]) || (resOpt$theta[3L] > 0.22)) &&
is.finite(lastGoodParameters[3L]) ) {
isUsingFixedAlpha <- TRUE
theta0Adj[3L] <- lastGoodParameters[3L]
resOpt <- .self$optimLRCOnAdjustedPrior(theta0Adj, iOpt = getIOpt(
isUsingFixedVPD, isUsingFixedAlpha)
, parameterPrior = parameterPrior, ctrl, dsDay = dsDay, ...)
}
} else {
# check alpha, if gt 0.22 estimate parameters with fixed alpha of last window
# if not last window exists, let alpha > 0.22
if ( (is.na(resOpt$theta[3L]) || (resOpt$theta[3L] > 0.22))
&& is.finite(lastGoodParameters[3L]) ) {
isUsingFixedAlpha <- TRUE
theta0Adj[3L] <- lastGoodParameters[3L]
resOpt <- .self$optimLRCOnAdjustedPrior(theta0Adj, iOpt = getIOpt(
isUsingFixedVPD, isUsingFixedAlpha)
, parameterPrior = parameterPrior, ctrl, dsDay = dsDay, ...)
# check k, if less than zero estimate parameters without VPD effect
# and with fixed alpha of last window
if (is.na(resOpt$theta[1L]) || (resOpt$theta[1L] < 0)) {
isUsingFixedVPD <- TRUE
theta0Adj[1L] <- 0
resOpt <- .self$optimLRCOnAdjustedPrior(theta0Adj, iOpt = getIOpt(
isUsingFixedVPD, isUsingFixedAlpha)
, parameterPrior = parameterPrior, ctrl, dsDay = dsDay, ...)
}
}
}
##details<<
## No parameters are reported if alpha<0 or RRef < 0 or beta0 < 0
## or beta0 > 250
# positions in theta0: "k" "beta0" "alpha" "RRef" "E0"
if (resOpt$convergence != 0) {
resOpt$theta <- NA
}
if (!is.na(resOpt$theta[1L]) && ((resOpt$theta[3L] < 0) ||
(resOpt$theta[4L] < 0) || (resOpt$theta[2L] < 0) ||
(resOpt$theta[2L] >= 250))
) {
# TODO estimate RRef from daytime data?
#LloydT_E0fix
#stop("case with alpha or beta < 0")
resOpt$theta[] <- NA
resOpt$convergence <- 1002
}
##details<<
## Not parameters are reported if the data did not contain records that
## are near light saturation.
## This is checked by comparing the prediction at highest PAR with the
## beta parameter
if (is.finite(resOpt$theta[1]) & !isTRUE(ctrl$isUsingLasslopQualityConstraints) &
is.finite(ctrl$minPropSaturation)
) {
dsDay <- list(...)$dsDay
iMaxRg <- which.max(dsDay$Rg)
dsDayMax <- dsDay[iMaxRg, ,drop = FALSE]
# compute prediction at maximum observed PAR and and PAR = 2000 and
# compare how close its to saturation
predMaxGPP <- .self$predictLRC(theta = resOpt$theta,
, Rg = c(dsDayMax$Rg, 2000)
, VPD = 0 #dsDayMax$VPD
, Temp = NA#dsDayMax$Temp
#, VPD0 = 10 # TODO: think of providing VPD0 and TRef to function
#, TRef = 15
)$GPP
if (predMaxGPP[1] < ctrl$minPropSaturation * predMaxGPP[2]) {
#plot(-NEE ~ Rg, dsDay)
resOpt$theta[] <- NA
resOpt$convergence <- 1006
}
}
# Further checks are done, after parameter uncertainty has been determined,
# by call from fitLRC to isParameterInBounds
#
##value<< list result of optimization as of
## \code{\link{LightResponseCurveFitter_optimLRCOnAdjustedPrior}} with entries
## \item{theta}{ numeric parameter vector that includes the fixed components}
## \item{iOpt}{ integer vector of indices of the vector that have been optimized}
## \item{convergence}{ scalar integer indicating bad conditions on fitting
## (see \code{\link{LightResponseCurveFitter_fitLRC}}) }
resOpt
}
LightResponseCurveFitter$methods(optimLRCBounds =
LightResponseCurveFitter_optimLRCBounds)
#' @export
LightResponseCurveFitter_getOptimizedParameterPositions <- function(
### get the positions of the parameters to optimize for given Fixed
isUsingFixedVPD ##<< boolean scalar: if TRUE,
## VPD effect set to zero and is not optimized
, isUsingFixedAlpha ##<< boolean scalar: if TRUE,
## initial slope is fixed and is not optimized
) {
##details<< If subclasses extend the parameter vector, they need
## to override this method.
iOpt <-
if (!isUsingFixedVPD & !isUsingFixedAlpha) c(1:4) else
if ( isUsingFixedVPD & !isUsingFixedAlpha) 2:4 else
if (!isUsingFixedVPD & isUsingFixedAlpha) c(1L, 2L, 4L) else
if ( isUsingFixedVPD & isUsingFixedAlpha) c(2L, 4L)
##value<< integer vector of positions in parameter vector
iOpt
}
LightResponseCurveFitter$methods(getOptimizedParameterPositions =
LightResponseCurveFitter_getOptimizedParameterPositions)
#' @export
LightResponseCurveFitter_optimLRCOnAdjustedPrior = function(
### Lower bound flux uncertainty and adjust prior uncertainty before calling optimLRC
theta ##<< numeric vector of starting values
, iOpt ##<< integer vector: positions of subset of parameters
## that are optimized
, dsDay ##<< dataframe of NEE, sdNEE and predictors Rg, VPD and Temp
, parameterPrior ##<< numeric vector of prior parameter estimates
## (corresponding to theta) # TODO rename to thetaPrior
, ctrl ##<< list of further controls
, ... ##<< further arguments to
## \code{\link{LightResponseCurveFitter_optimLRC}}
## (passed to \code{\link{LightResponseCurveFitter_computeCost}})
) {
if (!all(is.finite(theta))) stop("need to provide finite starting values.")
##details<<
## Only those records are used for optimization where both NEE
## and sdNEE are finite.
## In larger settings, already filtered at
dsDayFinite <- dsDay[is.finite(dsDay$NEE) & is.finite(dsDay$sdNEE), ]
if (nrow(dsDayFinite) < ctrl$minNRecInDayWindow) {
stop("inspect too few records, should be already filtered "
,"in partGLFitLRCOneWindow")
return(list(
theta = {theta[] <- NA; theta}
, iOpt = integer(0)
, convergence = 1003L
))
}
##details<< Optimization of LRC parameters takes into account the uncertainty
##of the flux values. In order to avoid very strong leverage, values with a
##very low uncertainty (< a lower quantile) are assigned the lower quantile is
##assigned. This procedure downweighs records with a high uncertainty, but does
##not apply a large leverage for records with a very low uncertainty. Avoid
##this correction by setting \code{ctrl$isBoundLowerNEEUncertainty =
##FALSE}
minUnc <- quantile(dsDayFinite$sdNEE, 0.3)
if (minUnc == 0) stop(
"Too many zeros in uncertainty of NEE.",
" This cannot be handled in daytime partitioning.")
Fc_unc <- if (isTRUE(ctrl$isBoundLowerNEEUncertainty) ) {
#twutz: avoid excessive weights by small uncertainties (of 1 / unc^2)
pmax(dsDayFinite$sdNEE, minUnc )
} else {
dsDayFinite$sdNEE
}
#plot(Fc_unc ~ dsDayFinite$sdNEE)
medianRelFluxUncertainty <- abs(median(Fc_unc / dsDayFinite$NEE))
##details<<
## The uncertainty of the prior, that maybe derived from fluxes) is allowed to
## adapt to the uncertainty of the fluxes.
## This is done in \code{link{LightResponseCurveFitter_getPriorScale}}
sdParameterPrior <- .self$getPriorScale(parameterPrior, medianRelFluxUncertainty
, nrow(dsDayFinite), ctrl = ctrl)
sdParameterPrior[-iOpt] <- NA
isUsingHessian <- (ctrl$nBootUncertainty == 0L)
.self$optimLRC(theta
, iOpt = iOpt
, flux = -dsDayFinite$NEE
, sdFlux = Fc_unc
, sdParameterPrior = sdParameterPrior
, parameterPrior = parameterPrior
, Rg = dsDayFinite$Rg
, VPD = dsDayFinite$VPD
, Temp = dsDayFinite$Temp
, isUsingHessian = isUsingHessian
, ctrl = ctrl
)
##value<< result of \code{\link{LightResponseCurveFitter_optimLRC}} with
## items theta, iOpt and convergence
}
LightResponseCurveFitter$methods(optimLRCOnAdjustedPrior =
LightResponseCurveFitter_optimLRCOnAdjustedPrior)
#' @export
LightResponseCurveFitter_isParameterInBounds <- function(
### Check if estimated parameter vector is within reasonable bounds
theta ##<< estimate of parameter
, sdTheta ##<< estimate of uncertainty of the parameter
, RRefNight ##<< numeric scalar: night-time based estimate of
## basal respiration
, ctrl ##<< list of further controls
) {
##author<< TW, MM
#
##details<< check the Beta bounds that depend on uncertainty:
## outside if (beta > 100 and sdBeta >= beta)
if (!is.finite(theta[2]) ) return(FALSE)
if (isTRUE(as.vector( (theta[2] > 100) && (sdTheta[2] >= theta[2]) )))
return(FALSE)
##value<< FALSE if parameters are outside reasonable bounds, TRUE otherwise
return(TRUE)
}
LightResponseCurveFitter$methods(isParameterInBounds =
LightResponseCurveFitter_isParameterInBounds)
#' @export
LightResponseCurveFitter_optimLRC <- function(
### call the optimization function
theta ##<< numeric vector: starting parameters
, iOpt ##<< integer vector: positions of parameters to optimize
, sdParameterPrior ##<< numeric vector: prior uncertainty
, ... ##<< further arguments to the cost function
, ctrl ##<< list of further controls
, isUsingHessian ##<< scalar boolean: set to TRUE to compute
## Hessian at optimum
) {
thetaOrig <- theta
resOptim <- optim(thetaOrig[iOpt], .self$computeCost
#resOptim <- optim(theta, .partGLRHLightResponseCost
#tmp <- .partGLRHLightResponseCost(theta[iOpt]
, theta = thetaOrig
, iOpt = iOpt
, sdParameterPrior = sdParameterPrior
, ...
, control = list(reltol = ctrl$LRCFitConvergenceTolerance)
, method = "BFGS", hessian = isUsingHessian)
##value<<
## list of result of \code{\link{optim}} amended with list
thetaOpt <- theta; thetaOpt[iOpt] <- resOptim$par
ans <- list(
theta = thetaOpt ##<< numeric vector: optimized parameter
## vector including the fixed components
, iOpt = iOpt ##<< integer vector: position of
## parameters that have been optimized
)
c(resOptim, ans)
}
LightResponseCurveFitter$methods(optimLRC = LightResponseCurveFitter_optimLRC)
#' @export
LightResponseCurveFitter_computeCost <- function(
### Computing residual sum of squares for predictions vs. data of NEE
thetaOpt ##<< parameter vector with components of theta0 that are optimized
, theta ##<< parameter vector with positions as in argument of
## \code{\link{LightResponseCurveFitter_getParameterNames}}
, iOpt ##<< position in theta that are optimized
, flux ##<< numeric: NEP (-NEE) or GPP time series [umolCO2 / m2 / s],
## should not contain NA
, sdFlux ##<< numeric: standard deviation of Flux [umolCO2 / m2 / s],
##should not contain NA
, parameterPrior ##<< numeric vector along theta: prior estimate of
## parameter (range of values)
, sdParameterPrior ##<< standard deviation of parameterPrior
#, weightMisfitPar2000 = NA ##<< weight of misfit of difference between
## saturation and prediction at PAR = 2000
, ... ##<< other arguments to
##\code{\link{LightResponseCurveFitter_predictLRC}}, such as VPD0, fixVPD
) {
theta[iOpt] <- thetaOpt
#print(theta)
resPred <- .self$predictLRC(theta, ...)
NEP_mod <- resPred$NEP
#if (is.na(mean(NEP_mod)) == TRUE) {
# recover()
#}
misFitPrior <- (((theta - parameterPrior)) / (sdParameterPrior))^2
misFitObs <- sum(((NEP_mod-flux) / sdFlux)^2)
RSS <- misFitObs + sum(misFitPrior, na.rm = TRUE)
#if (!is.finite(RSS) ) recover() # debugging the fit
RSS
}
LightResponseCurveFitter$methods(computeCost =
LightResponseCurveFitter_computeCost)
#' @export
LightResponseCurveFitter_predictLRC <- function(
### Light Response Function
theta ##<< numeric vector of parameters
, Rg ##<< ppfd [numeric] -> photosynthetic flux density
## [umol / m2 / s] or Global Radiation
, VPD ##<< VPD [numeric] -> Vapor Pressure Deficit [hPa]
, Temp ##<< Temp [degC] -> Temperature [degC]
, VPD0 = 10 ##<< VPD0 [hPa] -> Parameters VPD0 fixed to 10 hPa
## according to Lasslop et al 2010
, fixVPD = (k == 0) ##<< boolean scalar or vector of nrow theta:
## fixVPD if TRUE the VPD effect is not considered and VPD is not part
## of the computation
, TRef = 15 ##<< numeric scalar of Temperature (degree Celsius) for
## reference respiration RRef
) {
##details<<
## Predict ecosystem fluxes (Reco, GPP, NEP = GPP-Reco) for given parameters
## and environmental conditions.
##
## The VPD effect is included according to Lasslop et al., 2010.
##details<<
## If theta is a matrix, a different row of parameters is used for different
## entries of other inputs
if (is.matrix(theta) ) {
k <- theta[, 1]
beta <- theta[, 2]
alpha <- theta[, 3]
RRef <- theta[, 4]
E0 <- theta[, 5]
} else {
k <- theta[1]
beta <- theta[2]
alpha <- theta[3]
RRef <- theta[4]
E0 <- theta[5]
}
if (length(fixVPD) != length(VPD) )
if (length(fixVPD) == 1L) fixVPD <- rep(fixVPD, length(VPD) ) else
stop("Length of vector argument fixVPD must correspond to rows in theta.")
Amax <- ifelse(fixVPD, beta,
#ifelse(is.finite(VPD) & (VPD > VPD0), beta * exp(-k * (VPD-VPD0)), beta)
#deprecated: better filter for k: twutz: 170927: introduced pmin(1, ...)
# after looking at pvWave code, can happen if k is negative
ifelse((VPD > VPD0), beta * exp(-k * (VPD - VPD0)), beta)
)
Reco <- RRef * exp(E0 * (1 / ((273.15 + TRef) - 227.13) - 1 /
(Temp + 273.15 - 227.13)))
GPP <- .self$predictGPP(Rg, Amax = Amax, alpha = alpha)
NEP <- GPP - Reco
## a data.frame of length of Rg of computed
ans <- list(
NEP = NEP ##<< Net ecosystem production (-NEE), vector of length(Rg)
, Reco = Reco ##<< Ecosystem respiration
, GPP = GPP ##<< Gross primary production
)
}
LightResponseCurveFitter$methods( predictLRC = LightResponseCurveFitter_predictLRC)
#' @export
LightResponseCurveFitter_predictGPP <- function(
### Light Response function for GPP
Rg ##<< ppfd [numeric] -> photosynthetic flux density [mumol / m2 / s]
## or Global Radiation
, ... ##<< further parameters to the LRC
) {
##seealso<< \code{\link{partitionNEEGL}}
#
##details<<
## This method must be be implemented by a specific subclass.
## Currently there are several alternatives: \cr
## \itemize{
## \item Rectangular: \code{\link{RectangularLRCFitter_predictGPP}}
## \item Nonrectangular: \code{\link{NonrectangularLRCFitter_predictGPP}}
## \item Rectangular: \code{\link{LogisticSigmoidLRCFitter_predictGPP}}
## }
#
stop("Abstract method. Need to define in derived LRC class.")
##value<< numeric vector of length(Rg) of GPP
}
LightResponseCurveFitter$methods(predictGPP = LightResponseCurveFitter_predictGPP)
#' @export
LightResponseCurveFitter_computeLRCGradient <- function(
### Gradient of \code{\link{LightResponseCurveFitter_predictLRC}}
theta ##<< theta [numeric] -> parameter vector
## (theta[1] = k (k), theta[2] = beta (beta), theta[3] = alpha,
## theta[4] = RRef (rb), theta[4] = E0)
, Rg ##<< ppfd [numeric] -> photosynthetic flux density
## [umol / m2 / s] or Global Radiation
, VPD ##<< VPD [numeric] -> Vapor Pressure Deficit [hPa]
, Temp ##<< Temp [degC] -> Temperature [degC]
, VPD0 = 10 ##<< VPD0 [hPa] -> Parameters VPD0 fixed to 10
## hPa according to Lasslop et al 2010
, fixVPD = (k == 0) ##<< boolean scalar or vector of nrow(theta):
## fixVPD if TRUE the VPD effect is not considered and VPD is not part
## of the computation
, TRef = 15 ##<< numeric scalar of Temperature (degree Celsius) for
## reference respiration RRef
) {
if (is.matrix(theta) ) {
k <- theta[, 1]
beta <- theta[, 2]
alpha <- theta[, 3]
RRef <- theta[, 4]
E0 <- theta[, 5]
} else {
k <- theta[1]
beta <- theta[2]
alpha <- theta[3]
RRef <- theta[4]
E0 <- theta[5]
}
if (length(fixVPD) != length(VPD) )
if (length(fixVPD) == 1L) fixVPD <- rep(fixVPD, length(VPD) ) else
stop("Length of vector argument fixVPD must correspond to rows in theta.")
Amax <- ifelse(fixVPD, beta,
#ifelse(is.finite(VPD) & (VPD > VPD0), beta * exp(-k * (VPD-VPD0)), beta)
ifelse((VPD > VPD0), beta * exp(-k * (VPD - VPD0)), beta)
)
#ex <- expression(beta * exp(-k * (VPD - VPD0)) ); deriv(ex, c("beta", "k"))
dAmax_dkVPD <- ifelse(fixVPD, 0,
ifelse(VPD > VPD0, beta * -(VPD - VPD0) * exp(-k * (VPD - VPD0)), 0)
)
dAmax_dbeta0 <- ifelse(fixVPD, 0,
ifelse(VPD > VPD0, exp(-k * (VPD - VPD0)), 1)
)
#Reco<- RRef * exp(E0 * (1 / ((273.15 + 10)-227.13)-1 / (Temp + 273.15-227.13)))
#ex <- expression(RRef * exp(E0 * (1 / ((273.15 + TRef)-227.13)-1 /
# (Temp + 273.15-227.13))) ); deriv(ex, c("RRef", "E0"))
.expr7 <- 1 / (273.15 + TRef - 227.13) - 1 / (Temp + 273.15 - 227.13)
.expr9 <- exp(E0 * .expr7)
gradReco <- matrix(0, ncol = 2L, nrow = length(.expr9), dimnames =
list(NULL, c("RRef", "E0")))
gradReco[, "RRef"] <- dReco_dRRef <- .expr9
gradReco[, "E0"] <- dReco_dE0 <- RRef * (.expr9 * .expr7)
#
gradGPP <- array(0, c(nrow(gradReco), 3L), list(NULL, c("k", "beta", "alpha")))
dGPP_dAMax <- .self$computeGPPGradient(Rg, Amax, alpha)
gradGPP[, "beta"] <- dGPP_dAMax[, 1] * dAmax_dbeta0
gradGPP[, "k"] <- dGPP_dAMax[, 1] * dAmax_dkVPD
gradGPP[, "alpha"] <- dGPP_dAMax[, 2]
#NEP <- GPP - Reco
gradNEP <- cbind(gradGPP, -gradReco)
## list with gradient matrices. For each record
## (length(Rg)), c("k", "beta", "alpha", "RRef")
ans <- list(
NEP = gradNEP
, Reco = gradReco
, GPP = gradGPP
)
}
LightResponseCurveFitter$methods(computeLRCGradient =
LightResponseCurveFitter_computeLRCGradient)
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.