Skip to content

Commit

Permalink
second
Browse files Browse the repository at this point in the history
  • Loading branch information
kongdd committed Jan 1, 2025
1 parent b4a18c8 commit 4b44cd9
Show file tree
Hide file tree
Showing 22 changed files with 237 additions and 226 deletions.
14 changes: 7 additions & 7 deletions NAMESPACE
Original file line number Diff line number Diff line change
Expand Up @@ -2,12 +2,16 @@

export("%>%")
export(BerkeleyJulianDateToPOSIXct)
export(LRC_NonRect_computeLRCGradient)
export(LRC_NonRect_getParamNames)
export(LRC_NonRect_predictGPP)
export(LRC_NonRect_predictLRC)
export(LRC_computeCost)
export(LRC_computeLRCGradient)
export(LRC_fitLRC)
export(LRC_getOptimizedParameterPositions)
export(LRC_getParameterInitials)
export(LRC_getParameterNames)
export(LRC_getOptimizedParamPos)
export(LRC_getParamInitials)
export(LRC_getParamNames)
export(LRC_getPriorLocation)
export(LRC_getPriorScale)
export(LRC_isParameterInBounds)
Expand All @@ -21,10 +25,6 @@ export(LogisticSigmoidLRCFitter)
export(LogisticSigmoidLRCFitter_computeGPPGradient)
export(LogisticSigmoidLRCFitter_predictGPP)
export(NonrectangularLRCFitter)
export(NonrectangularLRCFitter_computeLRCGradient)
export(NonrectangularLRCFitter_getParameterNames)
export(NonrectangularLRCFitter_predictGPP)
export(NonrectangularLRCFitter_predictLRC)
export(POSIXctToBerkeleyJulianDate)
export(REddyProc_defaultunits)
export(RectangularLRCFitter)
Expand Down
21 changes: 10 additions & 11 deletions R/LRC_base.R
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,7 @@ LightResponseCurveFitter <- setRefClass("LightResponseCurveFitter")
# if this method is adjusted in subclass, then also need to adjust getPriorLocation,
# getPriorScale.

#' LRC_getParameterNames
#' LRC_getParamNames
#'
#' @details
#' - `k` : VPD effect
Expand All @@ -19,12 +19,11 @@ LightResponseCurveFitter <- setRefClass("LightResponseCurveFitter")
#'
#' @importFrom stats setNames
#' @export
LRC_getParameterNames <- function() {
LRC_getParamNames <- function() {
vars = c("k", "beta", "alpha", "RRef", "E0")
setNames(vars, vars)
}
LightResponseCurveFitter$methods(
getParameterNames = LRC_getParameterNames)
LightResponseCurveFitter$methods(getParamNames = LRC_getParamNames)

#' Optimize rectangular hyperbolic light response curve in one window
#'
Expand Down Expand Up @@ -75,11 +74,11 @@ LRC_fitLRC <- function(
controlGLPart = partGLControl(), lastGoodParameters = rep(NA_real_, 7L)) {

# Three initial guess vectors are defined according to Lasslop et al., 2010
parNames <- .self$getParameterNames() # hook method from derived classes
parNames <- .self$getParamNames() # hook method from derived classes
nPar <- length(parNames)

parPrior <- .self$getPriorLocation(dsDay$NEE, RRefNight = RRefNight, E0 = E0)
thetaInitials <- .self$getParameterInitials(parPrior)
thetaInitials <- .self$getParamInitials(parPrior)

resOpt3 <- apply(thetaInitials, 1, function(theta0) {
resOpt <- .self$optimLRCBounds(theta0, parPrior,
Expand Down Expand Up @@ -234,7 +233,7 @@ LightResponseCurveFitter$methods(getPriorLocation = LRC_getPriorLocation)

# ' @return A numeric matrix (3, nPar) of initial values for fitting parameters
#' @export
LRC_getParameterInitials <- function(thetaPrior) {
LRC_getParamInitials <- function(thetaPrior) {
theta0 <- matrix(rep(thetaPrior, each = 3), 3, length(thetaPrior),
dimnames = list(NULL, names(thetaPrior))
)
Expand All @@ -243,7 +242,7 @@ LRC_getParameterInitials <- function(thetaPrior) {
theta0 # thetaInitials
}
LightResponseCurveFitter$methods(
getParameterInitials = LRC_getParameterInitials)
getParamInitials = LRC_getParamInitials)

#' Optimize parameters with refitting with some fixed parameters if outside bounds
#'
Expand Down Expand Up @@ -277,7 +276,7 @@ LRC_optimLRCBounds <- function(
isUsingFixedVPD <- isNeglectVPDEffect || (sum(dsDay$VPD >= VPD0, na.rm = TRUE) == 0)
isUsingFixedAlpha <- FALSE

getIOpt <- .self$getOptimizedParameterPositions
getIOpt <- .self$getOptimizedParamPos

theta0Adj <- theta0 # initial estimate with some parameters adjusted to bounds
if (isNeglectVPDEffect) theta0Adj[1] <- 0
Expand Down Expand Up @@ -382,7 +381,7 @@ LightResponseCurveFitter$methods(
#'
#' @return integer vector of positions in parameter vector
#' @export
LRC_getOptimizedParameterPositions <- function(isUsingFixedVPD, isUsingFixedAlpha) {
LRC_getOptimizedParamPos <- function(isUsingFixedVPD, isUsingFixedAlpha) {
if (!isUsingFixedVPD & !isUsingFixedAlpha) {
c(1:4)
} else if (isUsingFixedVPD & !isUsingFixedAlpha) {
Expand All @@ -394,7 +393,7 @@ LRC_getOptimizedParameterPositions <- function(isUsingFixedVPD, isUsingFixedAlph
} # iOpt
}
LightResponseCurveFitter$methods(
getOptimizedParameterPositions = LRC_getOptimizedParameterPositions)
getOptimizedParamPos = LRC_getOptimizedParamPos)

#' optimLRCOnAdjustedPrior
#'
Expand Down
12 changes: 5 additions & 7 deletions R/LRC_base_optim.R
Original file line number Diff line number Diff line change
Expand Up @@ -73,11 +73,11 @@ LightResponseCurveFitter$methods(
#'
#' @param theta numeric vector of parameters. If theta is a matrix, a different
#' row of parameters is used for different entries of other inputs.
#' @param Rg -> photosynthetic flux density [mumol / m2 / s] or
#' @param Rg photosynthetic flux density [mumol / m2 / s] or
#' Global Radiation.
#' @param VPD -> Vapor Pressure Deficit [hPa]
#' @param Temp -> Temperature [degC]
#' @param VPD0 [hPa] -> Parameters VPD0 fixed to 10 hPa according to Lasslop et al 2010
#' @param VPD Vapor Pressure Deficit [hPa]
#' @param Temp Temperature [degC]
#' @param VPD0 [hPa], Parameters VPD0 fixed to 10 hPa according to Lasslop et al 2010
#' @param fixVPD if TRUE the VPD effect is not considered and VPD is not part of the computation
#' @param TRef numeric scalar of Temperature (degree Celsius) for reference respiration RRef
#'
Expand Down Expand Up @@ -215,6 +215,4 @@ LRC_computeLRCGradient <- function(
list(NEP = gradNEP, Reco = gradReco, GPP = gradGPP)
}

LightResponseCurveFitter$methods(
computeLRCGradient = LRC_computeLRCGradient
)
LightResponseCurveFitter$methods(computeLRCGradient = LRC_computeLRCGradient)
124 changes: 48 additions & 76 deletions R/LRC_nonrectangular.R
Original file line number Diff line number Diff line change
Expand Up @@ -7,74 +7,55 @@ NonrectangularLRCFitter <- setRefClass("NonrectangularLRCFitter",
contains = "LightResponseCurveFitter"
)

#' - `NonrectangularLRCFitter_getParameterNames`: return the parameter names
#' - `LRC_NonRect_getParamNames`: return the parameter names
#' used by this Light Response Curve Function
#'
#' @return string vector of parameter names. Positions are important.
#' Adds sixth parameter, `logitconv` to the parameters of
#' [LRC_getParameterNames()]
#' [LRC_getParamNames()]
#'
#' @seealso [NonrectangularLRCFitter_predictGPP()]
#' @seealso [LRC_NonRect_predictGPP()]
#' @export
NonrectangularLRCFitter_getParameterNames <- function() {
LRC_NonRect_getParamNames <- function() {
ans <- callSuper()
c(ans, logitconf = "logitconv") ## << logit-transformed convexity parameter.
## The value at original scale is obtained by
## \code{conv = 1 / (1 + exp(-logitconv))}
}

NonrectangularLRCFitter$methods(
getParameterNames = NonrectangularLRCFitter_getParameterNames)
getParamNames = LRC_NonRect_getParamNames)

# LRC_getPriorLocation
LRC_NonRect_getPriorLocation <- function(NEEDay, RRefNight, E0) {
ans <- callSuper(NEEDay = NEEDay, RRefNight = RRefNight, E0 = E0)
c(ans, logitconv = logit(0.75))
}
NonrectangularLRCFitter$methods(getPriorLocation = LRC_NonRect_getPriorLocation)

NonrectangularLRCFitter$methods(
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
) {
ans <- callSuper(NEEDay = NEEDay, RRefNight = RRefNight, E0 = E0)
## value<< a numeric vector with prior estimates of the parameters
c(ans,
logitconv = logit(0.75)
) ## << logit-transformed convexity parameter.
## The valua at original scale is obtained by
## \code{conv = 1 / (1 + exp(-logitconv))}
}
)

NonrectangularLRCFitter$methods(
#' Get the prior distribution of parameters
#' @param thetaPrior numeric vector of location of priors
#' @param medianRelFluxUncertainty numeric scalar: median across the relative
#' uncertainty of the flux values, i.e. sdNEE / NEE
#' @param nRec integer scalar: number of finite observations
#' @param ctrl list of further controls, with entry `isLasslopPriorsApplied`
getPriorScale = function(thetaPrior, medianRelFluxUncertainty, nRec, ctrl) {
ans <- callSuper(
thetaPrio = thetaPrior,
medianRelFluxUncertainty = medianRelFluxUncertainty, nRec = nRec, ctrl = ctrl)
## value<< adds NA prior for logitconv
c(ans, logitconv = NA)
}
)
#' Get the prior distribution of parameters
#' @keywords internal
#' @inheritParams LRC_getPriorScale
LRC_NonRect_getPriorScale <- function(thetaPrior, medianRelFluxUncertainty, nRec, ctrl) {
ans <- callSuper(thetaPrior, medianRelFluxUncertainty, nRec, ctrl)
c(ans, logitconv = NA)
}
NonrectangularLRCFitter$methods(getPriorScale = LRC_NonRect_getPriorScale)

#' Get the positions of the parameters to optimize for given Fixed
#'
#' @param isUsingFixedVPD boolean scalar: if TRUE, VPD effect set to zero and is not optimized
#' @param isUsingFixedAlpha boolean scalar: if TRUE, initial slope is fixed and is not optimized
#' @param isUsingFixedVPD boolean: if TRUE, VPD effect set to zero and is not optimized
#' @param isUsingFixedAlpha boolean: if TRUE, initial slope is fixed and is not optimized
#'
#' @return integer vector of positions of the parameters to optimize
NonrectangularLRCFitter_getOptimizedParameterPositions <- function(
LRC_NonRect_getOptimizedParamPos <- function(
isUsingFixedVPD, isUsingFixedAlpha) {
iOpt <- callSuper(isUsingFixedVPD = isUsingFixedVPD, isUsingFixedAlpha = isUsingFixedAlpha)
c(iOpt, 6) # add the convexity parameter
}
NonrectangularLRCFitter$methods(
getOptimizedParameterPositions =
NonrectangularLRCFitter_getOptimizedParameterPositions
getOptimizedParamPos =
LRC_NonRect_getOptimizedParamPos
)

#' Nonrectangular Hyperbolic Light Response function: (Gilmanov et al., 2003)
Expand All @@ -88,9 +69,9 @@ NonrectangularLRCFitter$methods(
#' @param fixVPD boolean scalar or vector of nrow theta: fixVPD if TRUE the VPD
#' effect is not considered and VPD is not part of the computation
#' @param TRef numeric scalar of Temperature (degree Celsius) for reference respiration RRef
#'
#'
#' @export
NonrectangularLRCFitter_predictLRC <- function(
LRC_NonRect_predictLRC <- function(
theta, Rg, VPD, Temp, VPD0 = 10, fixVPD = (k == 0), TRef = 15) {
if (is.matrix(theta)) {
k <- theta[, 1]
Expand All @@ -117,41 +98,41 @@ NonrectangularLRCFitter_predictLRC <- function(
}
}
Amax <- ifelse(fixVPD, beta,
ifelse((VPD > VPD0), beta * exp(-k * (VPD - VPD0)), beta)
)
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, conv = conv)
NEP <- GPP - Reco
list(NEP = NEP, Reco = Reco, GPP = GPP)
}
NonrectangularLRCFitter$methods(predictLRC = NonrectangularLRCFitter_predictLRC)
NonrectangularLRCFitter$methods(predictLRC = LRC_NonRect_predictLRC)

#' Nonrectangular Light Response function for GPP
#'
#'
#' @details
#' This function generalizes the \code{\link{RectangularLRCFitter_predictGPP}}
#' by adding the convexity parameter \code{conv}.
#' For \code{conv -> 0 (logitconv -> -Inf)}: approaches the rectangular hyperbolic.
#' For \code{conv -> 1 (logitconv -> + Inf)}: approaches a step function.
#' Expected values of \code{conv} are about 0.7-0.9 (Moffat 2012).
#'
#' @param Rg photosynthetic flux density [umol / m2 / s] or Global Radiation
#' @param Amax saturation (beta parameter) adjusted for effect of VPD
#' @param alpha slope at Rg = 0
#' @param conv convexity parameter (see details)
#'
#'
#' @seealso [LRC_predictGPP()]
#' @return vector of GPP
#' @export
NonrectangularLRCFitter_predictGPP <- function(Rg, Amax, alpha, conv) {
## details<<
## This function generalizes the \code{\link{RectangularLRCFitter_predictGPP}}
## by adding the convexity parameter \code{conv}.
## For \code{conv -> 0 (logitconv -> -Inf)}: approaches the rectangular hyperbolic.
## For \code{conv -> 1 (logitconv -> + Inf)}: approaches a step function.
## Expected values of \code{conv} are about 0.7-0.9 (Moffat 2012).
LRC_NonRect_predictGPP <- function(Rg, Amax, alpha, conv) {
zRoot <- ((alpha * Rg + Amax)^2) - (4 * alpha * Rg * conv * Amax)
zRoot[which(zRoot < 0)] <- 0
## value<< numeric vector of length(Rg) of GPP
GPP <- (1 / (2 * conv)) * (alpha * Rg + Amax - sqrt(zRoot))
}
NonrectangularLRCFitter$methods(predictGPP = NonrectangularLRCFitter_predictGPP)
NonrectangularLRCFitter$methods(predictGPP = LRC_NonRect_predictGPP)


#' Gradient of [NonrectangularLRCFitter_predictLRC()]
#' Gradient of [LRC_NonRect_predictLRC()]
#'
#' @details This function computes the gradient of the Nonrectangular Light Response
#' Curve function with respect to the parameters.
Expand All @@ -165,19 +146,10 @@ NonrectangularLRCFitter$methods(predictGPP = NonrectangularLRCFitter_predictGPP)
#' theta[3] = alpha,
#' theta[4] = RRef (rb), # theta[4] = E0,
#' theta[5] = logitconv)
#' @param Rg photosynthetic flux density [umol / m2 / s] or Global Radiation
#' @param VPD Vapor Pressure Deficit [hPa]
#' @param Temp Temperature [degC]
#' @param VPD0 [hPa] -> Parameters VPD0 fixed to 10 hPa according to
#' Lasslop et al 2010
#' @param fixVPD boolean scalar or vector of nrow(theta):
#' fixVPD if TRUE the VPD effect is not considered and VPD is not part
#' of the computation
#' @param TRef numeric scalar of Temperature (degree Celsius) for reference
#' respiration `RRef`
#' @inheritParams LRC_predictLRC
#'
#' @export
NonrectangularLRCFitter_computeLRCGradient <- function(
LRC_NonRect_computeLRCGradient <- function(
theta, Rg, VPD, Temp, VPD0 = 10, fixVPD = (k == 0), TRef = 15) {
if (is.matrix(theta)) {
k <- theta[, 1]
Expand Down Expand Up @@ -239,17 +211,17 @@ NonrectangularLRCFitter_computeLRCGradient <- function(
}

NonrectangularLRCFitter$methods(
computeLRCGradient = NonrectangularLRCFitter_computeLRCGradient)
computeLRCGradient = LRC_NonRect_computeLRCGradient)

.tmp.f <- function() {
iNonFinite <- which(!is.finite(gradReco[, "E0"]))
}

#' Gradient of [NonrectangularLRCFitter_predictGPP()]
#' Gradient of [LRC_NonRect_predictGPP()]
#'
#' @inheritParams NonrectangularLRCFitter_predictGPP
#' @inheritParams LRC_NonRect_predictGPP
#' @param logitconv -> logit-transformed convexity parameter
NonrectangularLRCFitter_computeGPPGradient <- function(Rg, Amax, alpha, logitconv) {
LRC_NonRect_computeGPPGradient <- function(Rg, Amax, alpha, logitconv) {
zRoot <- ((alpha * Rg + Amax)^2) - (4 * alpha * Rg * invlogit(logitconv) * Amax)
iNegRoot <- which(zRoot < 0)
# GPP<- (1 / (2 * (1 / (1 + exp(-logitconv))) )) * (alpha * Rg + Amax-sqrt(zRoot))
Expand Down Expand Up @@ -297,4 +269,4 @@ NonrectangularLRCFitter_computeGPPGradient <- function(Rg, Amax, alpha, logitcon
}

NonrectangularLRCFitter$methods(
computeGPPGradient = NonrectangularLRCFitter_computeGPPGradient)
computeGPPGradient = LRC_NonRect_computeGPPGradient)
2 changes: 1 addition & 1 deletion R/LRC_rectangular.R
Original file line number Diff line number Diff line change
Expand Up @@ -60,7 +60,7 @@ RectangularLRCFitterCVersion <- setRefClass('RectangularLRCFitterCVersion',
#'
#' @param thetaOpt parameter vector with components of theta0 that are optimized
#' @param theta parameter vector with positions as in argument of
#' [LRC_getParameterNames()]
#' [LRC_getParamNames()]
#' @param iOpt position in theta that are optimized
#' @param flux numeric: NEP (-NEE) or GPP time series [umolCO2 / m2 / s],
#' should not contain NA
Expand Down
Loading

0 comments on commit 4b44cd9

Please sign in to comment.