Skip to content

Commit

Permalink
format docs
Browse files Browse the repository at this point in the history
  • Loading branch information
@kongdd
kongdd committed Jan 1, 2025
1 parent 2a10064 commit b4a18c8
Show file tree
Hide file tree
Showing 48 changed files with 2,936 additions and 630 deletions.
28 changes: 14 additions & 14 deletions NAMESPACE
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,21 +2,21 @@

export("%>%")
export(BerkeleyJulianDateToPOSIXct)
export(LRC_computeCost)
export(LRC_computeLRCGradient)
export(LRC_fitLRC)
export(LRC_getOptimizedParameterPositions)
export(LRC_getParameterInitials)
export(LRC_getParameterNames)
export(LRC_getPriorLocation)
export(LRC_getPriorScale)
export(LRC_isParameterInBounds)
export(LRC_optimLRC)
export(LRC_optimLRCBounds)
export(LRC_optimLRCOnAdjustedPrior)
export(LRC_predictGPP)
export(LRC_predictLRC)
export(LightResponseCurveFitter)
export(LightResponseCurveFitter_computeCost)
export(LightResponseCurveFitter_computeLRCGradient)
export(LightResponseCurveFitter_fitLRC)
export(LightResponseCurveFitter_getOptimizedParameterPositions)
export(LightResponseCurveFitter_getParameterInitials)
export(LightResponseCurveFitter_getParameterNames)
export(LightResponseCurveFitter_getPriorLocation)
export(LightResponseCurveFitter_getPriorScale)
export(LightResponseCurveFitter_isParameterInBounds)
export(LightResponseCurveFitter_optimLRC)
export(LightResponseCurveFitter_optimLRCBounds)
export(LightResponseCurveFitter_optimLRCOnAdjustedPrior)
export(LightResponseCurveFitter_predictGPP)
export(LightResponseCurveFitter_predictLRC)
export(LogisticSigmoidLRCFitter)
export(LogisticSigmoidLRCFitter_computeGPPGradient)
export(LogisticSigmoidLRCFitter_predictGPP)
Expand Down
51 changes: 19 additions & 32 deletions R/EddyGapfilling.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -441,6 +441,11 @@ sEddyProc$methods(sFillMDC = sEddyProc_sFillMDC)

#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

#' sEddyProc_sMDSGapFill
#'
#' MDS gap filling algorithm adapted after the PV-Wave code and paper by Markus
#' Reichstein.
#'
#' @param Var Variable to be filled
#' @param QFVar Quality flag of variable to be filled
#' @param QFValue Value of quality flag for _good_ (original) data, other data is set to missing
Expand All @@ -464,38 +469,20 @@ sEddyProc$methods(sFillMDC = sEddyProc_sFillMDC)
#'
#' @export
sEddyProc_sMDSGapFill <- function(
### MDS gap filling algorithm adapted after the PV-Wave code and paper by Markus Reichstein.
Var = Var.s ##<< Variable to be filled
, QFVar = if (!missing(QFVar.s)) QFVar.s else 'none' ##<< Quality flag
## of variable to be filled
, QFValue = if (!missing(QFValue.n)) QFValue.n else NA_real_ ##<< Value of
## quality flag for _good_ (original) data, other data is set to missing
, V1 = if (!missing(V1.s)) V1.s else 'Rg' ##<< Condition variable 1
## (default: Global radiation 'Rg' in W m-2)
, T1 = if (!missing(T1.n)) T1.n else 50 ##<< Tolerance interval 1
## (default: 50 W m-2)
, V2 = if (!missing(V2.s)) V2.s else 'VPD' ##<< Condition variable 2
## (default: Vapour pressure deficit 'VPD' in hPa)
, T2 = if (!missing(T2.n)) T2.n else 5 ##<< Tolerance interval 2
## (default: 5 hPa)
, V3 = if (!missing(V3.s)) V3.s else 'Tair' ##<< Condition variable 3
## (default: Air temperature 'Tair' in degC)
, T3 = if (!missing(T3.n)) T3.n else 2.5 ##<< Tolerance interval 3
## (default: 2.5 degC)
, FillAll = if (!missing(FillAll.b)) FillAll.b else TRUE ##<< Fill
## all values to estimate uncertainties
, isVerbose = if (!missing(Verbose.b)) Verbose.b else TRUE ##<< Print
## status information to screen
, suffix = if (!missing(Suffix.s)) Suffix.s else '' ##<< String
## suffix needed for different processing setups on the same dataset
## (for explanations see below)
, minNWarnRunLength = ##<< scalar integer:
## warn if number of subsequent
## numerically equal values exceeds this number.
## Set to Inf or NA for no warnings.
## defaults for "NEE" to records across 4 hours and no warning for others.
if (Var == "NEE") 4 * .self$sINFO$DTS/24 else NA_integer_
, Var.s ##<< deprecated
Var = Var.s,
QFVar = if (!missing(QFVar.s)) QFVar.s else 'none',
QFValue = if (!missing(QFValue.n)) QFValue.n else NA_real_,
V1 = if (!missing(V1.s)) V1.s else 'Rg',
T1 = if (!missing(T1.n)) T1.n else 50,
V2 = if (!missing(V2.s)) V2.s else 'VPD',
T2 = if (!missing(T2.n)) T2.n else 5,
V3 = if (!missing(V3.s)) V3.s else 'Tair',
T3 = if (!missing(T3.n)) T3.n else 2.5,
FillAll = if (!missing(FillAll.b)) FillAll.b else TRUE,
isVerbose = if (!missing(Verbose.b)) Verbose.b else TRUE,
suffix = if (!missing(Suffix.s)) Suffix.s else '',
minNWarnRunLength = if (Var == "NEE") 4 * .self$sINFO$DTS/24 else NA_integer_,
Var.s ##<< deprecated
, QFVar.s ##<< deprecated
, QFValue.n ##<< deprecated
, V1.s ##<< deprecated
Expand Down
97 changes: 48 additions & 49 deletions R/LRC_base.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,10 +5,10 @@
#' @exportClass LightResponseCurveFitter
LightResponseCurveFitter <- setRefClass("LightResponseCurveFitter")

# if this method is adjusted in subclass, then also adjust getPriorLocation,
# if this method is adjusted in subclass, then also need to adjust getPriorLocation,
# getPriorScale.

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

LightResponseCurveFitter$methods(
getParameterNames = LightResponseCurveFitter_getParameterNames)
getParameterNames = LRC_getParameterNames)

#' Optimize rectangular hyperbolic light response curve in one window
#'
Expand Down Expand Up @@ -69,9 +68,9 @@ LightResponseCurveFitter$methods(
#' + 1011 : too few valid records in window (from different location:
#' partGLFitLRCOneWindow)
#'
#' @seealso [partGLFitLRCWindows()], [LightResponseCurveFitter_optimLRCBounds()]
#' @seealso [partGLFitLRCWindows()], [LRC_optimLRCBounds()]
#' @export
LightResponseCurveFitter_fitLRC <- function(
LRC_fitLRC <- function(
dsDay, E0, sdE0, RRefNight,
controlGLPart = partGLControl(), lastGoodParameters = rep(NA_real_, 7L)) {

Expand All @@ -86,9 +85,10 @@ LightResponseCurveFitter_fitLRC <- function(
resOpt <- .self$optimLRCBounds(theta0, parPrior,
dsDay = dsDay, ctrl = controlGLPart, lastGoodParameters = lastGoodParameters)
})
iValid <- which(sapply(resOpt3, function(resOpt) { is.finite(resOpt$theta[1]) }))
iValid <- which(sapply(resOpt3, \(r) is.finite(r$theta[1]) ))
resOpt3Valid <- resOpt3[iValid]
optSSE <- sapply(resOpt3Valid, "[[", "value")

getNAResult <- function(convergenceCode) {
list(
thetaOpt = structure(rep(NA_real_, nPar), names = colnames(thetaInitials)),
Expand Down Expand Up @@ -172,10 +172,10 @@ LightResponseCurveFitter_fitLRC <- function(
covParms = covParms,
convergence = resOpt$convergence)
}
LightResponseCurveFitter$methods(fitLRC = LightResponseCurveFitter_fitLRC)
LightResponseCurveFitter$methods(fitLRC = LRC_fitLRC)


#' LightResponseCurveFitter_getPriorScale
#' LRC_getPriorScale
#'
#' @details
#' The beta parameter is quite well defined. Hence use a prior with
Expand All @@ -199,7 +199,7 @@ LightResponseCurveFitter$methods(fitLRC = LightResponseCurveFitter_fitLRC)
#'
#' @return a numeric vector with prior estimates of the parameters
#' @export
LightResponseCurveFitter_getPriorScale <- function(
LRC_getPriorScale <- function(
thetaPrior, medianRelFluxUncertainty, nRec, ctrl) {
sdParameterPrior <- if (ctrl$isLasslopPriorsApplied) {
# twutz: changed to no prior for logitconv
Expand All @@ -209,33 +209,32 @@ LightResponseCurveFitter_getPriorScale <- function(
c(k = NA, beta = as.vector(sdBetaPrior), alpha = NA, RRef = NA, E0 = NA)
}
}
LightResponseCurveFitter$methods(getPriorScale = LightResponseCurveFitter_getPriorScale)

LightResponseCurveFitter$methods(getPriorScale = LRC_getPriorScale)

#' LightResponseCurveFitter_getPriorLocation
#' LRC_getPriorLocation
#'
#' @param NEEDay numeric vector of daytime NEE
#' @param RRefNight numeric scalar of basal respiration estimated from night-time data
#' @param E0 numeric scalar of night-time estimate of temperature sensitivity
#'
#'
#' @return a numeric vector with prior estimates of the parameters
#' @export
LightResponseCurveFitter_getPriorLocation <- function(NEEDay, RRefNight, E0) {
LRC_getPriorLocation <- function(NEEDay, RRefNight, E0) {
if (!is.finite(RRefNight)) stop("must provide finite RRefNight")
zs = quantile(NEEDay, c(0.03, 0.97), na.rm = TRUE)
c(k = 0.05,
zs <- quantile(NEEDay, c(0.03, 0.97), na.rm = TRUE)
c(
k = 0.05,
beta = as.vector(abs(diff(zs))),
alpha = 0.1,
RRef = as.vector(RRefNight),
alpha = 0.1,
RRef = as.vector(RRefNight),
E0 = as.vector(E0)
) # parameterPrior
}
LightResponseCurveFitter$methods(
getPriorLocation = LightResponseCurveFitter_getPriorLocation)
LightResponseCurveFitter$methods(getPriorLocation = LRC_getPriorLocation)

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

#' Optimize parameters with refitting with some fixed parameters if outside bounds
#'
Expand All @@ -259,15 +258,15 @@ LightResponseCurveFitter$methods(
#' @param lastGoodParameters numeric vector of last successful fit
#' @param ctrl list of further controls, such as `isNeglectVPDEffect = TRUE`
#'
#' @return list of result of [LightResponseCurveFitter_optimLRCOnAdjustedPrior()]
#' @return list of result of [LRC_optimLRCOnAdjustedPrior()]
#' - `theta`: vector of optimized parameters
#' - `iOpt`: vector of positions of parameters that are optimized
#' - `thetaInitialGuess`: initial guess from data
#' see [LightResponseCurveFitter_fitLRC()]
#' see [LRC_fitLRC()]
#'
#' @seealso [LightResponseCurveFitter_fitLRC()]
#' @seealso [LRC_fitLRC()]
#' @export
LightResponseCurveFitter_optimLRCBounds <- function(
LRC_optimLRCBounds <- function(
theta0, parameterPrior, ..., dsDay, lastGoodParameters, ctrl) {
# twutz 161014: default alpha
if (!is.finite(lastGoodParameters[3L])) lastGoodParameters[3L] <- 0.22 # default alpha
Expand All @@ -284,7 +283,8 @@ LightResponseCurveFitter_optimLRCBounds <- function(
if (isNeglectVPDEffect) theta0Adj[1] <- 0
resOpt <- resOpt0 <- .self$optimLRCOnAdjustedPrior(theta0Adj,
iOpt = getIOpt(isUsingFixedVPD, isUsingFixedAlpha),
parameterPrior = parameterPrior, ctrl, dsDay = dsDay, ...)
parameterPrior = parameterPrior, ctrl, dsDay = dsDay, ...
)

fun_optim <- function() {
.self$optimLRCOnAdjustedPrior(theta0Adj,
Expand Down Expand Up @@ -371,7 +371,7 @@ LightResponseCurveFitter_optimLRCBounds <- function(
resOpt
}
LightResponseCurveFitter$methods(
optimLRCBounds = LightResponseCurveFitter_optimLRCBounds)
optimLRCBounds = LRC_optimLRCBounds)

#' Get the positions of the parameters to optimize for given Fixed
#'
Expand All @@ -382,7 +382,7 @@ LightResponseCurveFitter$methods(
#'
#' @return integer vector of positions in parameter vector
#' @export
LightResponseCurveFitter_getOptimizedParameterPositions <- function(isUsingFixedVPD, isUsingFixedAlpha) {
LRC_getOptimizedParameterPositions <- function(isUsingFixedVPD, isUsingFixedAlpha) {
if (!isUsingFixedVPD & !isUsingFixedAlpha) {
c(1:4)
} else if (isUsingFixedVPD & !isUsingFixedAlpha) {
Expand All @@ -394,7 +394,7 @@ LightResponseCurveFitter_getOptimizedParameterPositions <- function(isUsingFixed
} # iOpt
}
LightResponseCurveFitter$methods(
getOptimizedParameterPositions = LightResponseCurveFitter_getOptimizedParameterPositions)
getOptimizedParameterPositions = LRC_getOptimizedParameterPositions)

#' optimLRCOnAdjustedPrior
#'
Expand All @@ -416,45 +416,44 @@ LightResponseCurveFitter$methods(
#' @param parameterPrior numeric vector of prior parameter estimates (corresponding to theta)
#' @param ctrl list of further controls
#' @param ... further arguments to
#' [LightResponseCurveFitter_optimLRC()] (passed to
#' [LightResponseCurveFitter_computeCost()])
#' [LRC_optimLRC()] (passed to
#' [LRC_computeCost()])
#'
#' @return list of result of [LightResponseCurveFitter_optimLRC()] amended with list
#' @return list of result of [LRC_optimLRC()] amended with list
#' `theta`, `iOpt` and `convergence`
#' @export
LightResponseCurveFitter_optimLRCOnAdjustedPrior <- function(
LRC_optimLRCOnAdjustedPrior <- function(
theta, iOpt, dsDay, parameterPrior, ctrl, ...) {

if (!all(is.finite(theta))) stop("need to provide finite starting values.")
dsDayFinite <- dsDay[is.finite(dsDay$NEE) & is.finite(dsDay$sdNEE), ]
dsDayFinite <- subset(dsDay, is.finite(NEE) & is.finite(sdNEE))

npar = length(theta)
if (nrow(dsDayFinite) < ctrl$minNRecInDayWindow) {
stop(
"inspect too few records, should be already filtered ",
stop("inspect too few records, should be already filtered ",
"in partGLFitLRCOneWindow")
return(list(
theta = {
theta[] <- NA
theta
}, iOpt = integer(0), convergence = 1003L))
theta = rep(NA, npar), iOpt = integer(0), convergence = 1003L))
}

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 [LightResponseCurveFitter_getPriorScale()]
## This is done in [LRC_getPriorScale()]
sdParameterPrior <- .self$getPriorScale(parameterPrior, medianRelFluxUncertainty,
nrow(dsDayFinite), ctrl = ctrl)
sdParameterPrior[-iOpt] <- NA
Expand All @@ -474,7 +473,7 @@ LightResponseCurveFitter_optimLRCOnAdjustedPrior <- function(
}

LightResponseCurveFitter$methods(
optimLRCOnAdjustedPrior = LightResponseCurveFitter_optimLRCOnAdjustedPrior)
optimLRCOnAdjustedPrior = LRC_optimLRCOnAdjustedPrior)

#' isParameterInBounds
#'
Expand All @@ -489,10 +488,10 @@ LightResponseCurveFitter$methods(
#'
#' @return logical scalar: TRUE if parameters are within bounds, FALSE otherwise
#' @export
LightResponseCurveFitter_isParameterInBounds <- function(theta, sdTheta, RRefNight, ctrl) {
LRC_isParameterInBounds <- function(theta, sdTheta, RRefNight, ctrl) {
if (!is.finite(theta[2])) return(FALSE)
if (isTRUE(as.vector((theta[2] > 100) && (sdTheta[2] >= theta[2])))) return(FALSE)
return(TRUE)
}
LightResponseCurveFitter$methods(
isParameterInBounds = LightResponseCurveFitter_isParameterInBounds)
isParameterInBounds = LRC_isParameterInBounds)
Loading

0 comments on commit b4a18c8

Please sign in to comment.