Skip to content

Commit

Permalink
Merge branch 'develop' into develop
Browse files Browse the repository at this point in the history
  • Loading branch information
istfer authored Jul 12, 2024
2 parents 7f1ffbc + ef1a44f commit 8abc90f
Show file tree
Hide file tree
Showing 14 changed files with 129 additions and 173 deletions.
13 changes: 13 additions & 0 deletions CHANGELOG.md
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,19 @@ For more information about this file see also [Keep a Changelog](http://keepacha

## [Unreleased]

### Added

### Fixed

### Changed
- The following components have changed their licensing. With approval of all their contributors, we now provide them under a BSD 3-clause license rather than the previously used NCSA Open Source license. As a reminder, we intend to relicense the entire system and this list will expand as we gather permission from the relevant copyright owners.
* Shiny apps `dbsync` and `Elicitation`

### Removed


## [1.8.0] - 2024-07-12

### Added
- Created a new soilgrids function to extract the mean soil organic carbon profile with associated undertainty values at each depth for any lat/long location (#3040). Function was created for the CMS SDA workflow

Expand Down
2 changes: 1 addition & 1 deletion modules/assim.sequential/inst/covariates.R
Original file line number Diff line number Diff line change
Expand Up @@ -99,7 +99,7 @@ NA_SoilGrids <- terra::crop(SoilGrids, NA_extent)
names(NA_SoilGrids) <- c("phh2o", "nitrogen", "soc", "sand")

# Resample SoilGrids to match WorldClim maps
NA_SoilGrids <- terra::resample(NA_WorldClim, NA_SoilGrids, method = 'bilinear') # made change
NA_SoilGrids <- terra::resample(NA_SoilGrids, NA_WorldClim, method = 'bilinear')

# Resample land cover to match WorldClim maps (~25 min)
land_cover <- terra::resample(land_cover, NA_WorldClim, method = 'near')
Expand Down
22 changes: 16 additions & 6 deletions modules/data.atmosphere/R/read.register.R
Original file line number Diff line number Diff line change
Expand Up @@ -35,14 +35,24 @@ read.register <- function(register.xml, con) {
} else if ((!is.null(register$format$id) & is.null(register$format$name))
|
(!is.null(register$format$id) & is.null(register$format$mimetype))) {
register$format$name <- PEcAn.DB::db.query(
paste("SELECT name from formats where id = ", register$format$id), con)[[1]]
register$format$mimetype <- PEcAn.DB::db.query(
paste("SELECT mime_type from formats where id = ", register$format$id), con)[[1]]
# Retrieve format name and mimetype from the database
query.format.info <- PEcAn.DB::db.query(
paste(
"SELECT name, type_string AS mimetype",
"FROM formats JOIN mimetypes ON formats.mimetype_id = mimetypes.id",
"WHERE formats.id = ", register$format$id),
con
)

register$format$name <- query.format.info$name
register$format$mimetype <- query.format.info$mimetype

} else if (is.null(register$format$id) & !is.null(register$format$name) & !is.null(register$format$mimetype)) {
register$format$id <- PEcAn.DB::db.query(
paste0("SELECT id from formats where name = '", register$format$name,
"' and mime_type = '", register$format$mimetype, "'"), con)[[1]]
paste0(
"SELECT formats.id FROM formats JOIN mimetypes ON formats.mimetype_id = mimetypes.id ",
"WHERE name = '", register$format$name,
"' AND type_string = '", register$format$mimetype, "'"), con)[[1]]
}
}
return(invisible(register))
Expand Down
2 changes: 1 addition & 1 deletion modules/data.mining/DESCRIPTION
Original file line number Diff line number Diff line change
Expand Up @@ -21,4 +21,4 @@ LazyLoad: yes
LazyData: FALSE
Collate:
Encoding: UTF-8
RoxygenNote: 7.1.2
RoxygenNote: 7.3.1
126 changes: 70 additions & 56 deletions modules/uncertainty/R/sensitivity.analysis.R
Original file line number Diff line number Diff line change
Expand Up @@ -10,23 +10,20 @@
##' Spline estimate of univariate relationship between parameter value and model output
##'
##' Creates a spline function using the splinefun function that estimates univariate response of parameter input to model output
##' @name sa.splinefun
##' @title Sensitivity spline function
##' @param quantiles.input
##' @param quantiles.output
##'
##' @param quantiles.input passed to `x` argument of `stats::splinefun`
##' @param quantiles.output passed to `y` argument of `stats::splinefun`
##' @export
##' @return function
##' @return function
sa.splinefun <- function(quantiles.input, quantiles.output) {
return(stats::splinefun(quantiles.input, quantiles.output, method = "monoH.FC"))
return(stats::splinefun(x = quantiles.input, y = quantiles.output, method = "monoH.FC"))
} # sa.splinefun


#--------------------------------------------------------------------------------------------------#
##' Calculates the standard deviation of the variance estimate
##'
##' Uses the equation \eqn{\sigma^4\left(\frac{2}{n-1}+\frac{\kappa}{n}\right)}{\sigma^4 (2/(n-1) + \kappa/n)}
##' @name sd.var
##' @title Standard deviation of sample variance
##' @param x sample
##' @return estimate of standard deviation of the sample variance
##' @export
Expand All @@ -42,17 +39,15 @@ sd.var <- function(x) {
##'
##' Note that this calculates the 'excess kurtosis', which is defined as kurtosis - 3.
##' This statistic is used in the calculation of the standard deviation of sample variance
##' in the function \code{\link{sd.var}}.
##' Additional details
##' @name kurtosis
##' @title Calculate excess kurtosis from a vector
##' in the function \code{\link{sd.var}}.
##'
##' @param x vector of values
##' @return numeric value of kurtosis
##' @author David LeBauer
##' @references NIST/SEMATECH e-Handbook of Statistical Methods, \url{http://www.itl.nist.gov/div898/handbook/eda/section3/eda35b.htm}, 2011-06-20.
kurtosis <- function(x) {
kappa <- sum((x - mean(x, na.rm = TRUE)) ^ 4) /
((sum(!is.na(x)) - 1) * stats::sd(x, na.rm = TRUE) ^ 4) - 3
kappa <- sum((x - mean(x, na.rm = TRUE))^4) /
((sum(!is.na(x)) - 1) * stats::sd(x, na.rm = TRUE)^4) - 3
return(kappa)
} # kurtosis
# ==================================================================================================#
Expand All @@ -64,10 +59,9 @@ kurtosis <- function(x) {
##' This function evaluates the sensitivity of a model to a parameter.
##' This is done by evaluating the first derivative of the univariate spline estimate
##' of the model response at the parameter median.
##' @name get.sensitivity
##' @title Calculate Sensitivity
##' @param trait.samples
##' @param sa.splinefun
##'
##' @param trait.samples parameter values to evaluate at their median
##' @param sa.splinefun fitted spline function. Must take two arguments.
##' @export
##' @return numeric estimate of model sensitivity to parameter
get.sensitivity <- function(trait.samples, sa.splinefun) {
Expand All @@ -77,10 +71,10 @@ get.sensitivity <- function(trait.samples, sa.splinefun) {


#--------------------------------------------------------------------------------------------------#
##' Get coefficient of variance
##'
##' Given a set of numbers (a numeric vector), this returns the set's coefficient of variance.
##'
##' @name get.coef.var
##' @title Get coefficient of variance
##' @param set numeric vector of trait values
##' @export
##' @return coeficient of variance
Expand All @@ -93,86 +87,106 @@ get.coef.var <- function(set) {
##' Generic function for the elasticity
##'
##' Given the sensitivity, samples, and outputs for a single trait, return elasticity
##' @name get.elasticity
##' @title Get Elasticity
##'
##' @param sensitivity univariate sensitivity of model to a parameter, can be calculated by \code{\link{get.sensitivity}}
##' @param samples samples from trait distribution
##' @param outputs model output from ensemble runs
##' @export
##' @return elasticity = normalized sensitivity
##' @return elasticity = normalized sensitivity
get.elasticity <- function(sensitivity, samples, outputs) {
return(sensitivity / (stats::median(outputs) / stats::median(samples)))
} # get.elasticity


#--------------------------------------------------------------------------------------------------#
##' Performs univariate sensitivity analysis and variance decomposition
##' Performs univariate sensitivity analysis and variance decomposition
##'
##' This function estimates the univariate responses of a model to a parameter for a set of traits, calculates the model sensitivity at the median, and performs a variance decomposition. This function results in a set of sensitivity plots (one per variable) and plot_variance_decomposition.
##' @name sensitivity.analysis
##' @title Sensitivity Analysis
##' @param trait.samples list of vectors, one per trait, representing samples of the trait value, with length equal to the mcmc chain length. Samples are taken from either the prior distribution or meta-analysis results
##' @param sa.samples data.frame with one column per trait and one row for the set of quantiles used in sensitivity analysis. Each cell contains the value of the trait at the given quantile.
##' @param sa.output list of data.frames, similar to sa.samples, except cells contain the results of a model run with that trait x quantile combination and all other traits held at their median value
##' This function estimates the univariate responses of a model to a parameter for a set of traits, calculates the model sensitivity at the median,
##' and performs a variance decomposition. This function results in a set of sensitivity plots (one per variable) and plot_variance_decomposition.
##'
##' @param trait.samples list of vectors, one per trait, representing samples of the trait value, with length equal to the mcmc chain length.
##' Samples are taken from either the prior distribution or meta-analysis results
##' @param sa.samples data.frame with one column per trait and one row for the set of quantiles used in sensitivity analysis.
## Each cell contains the value of the trait at the given quantile.
##' @param sa.output list of data.frames, similar to sa.samples, except cells contain the results of a model run
##' with that trait x quantile combination and all other traits held at their median value
##' @param outdir directory to which plots are written
##' @return results of sensitivity analysis
##' @export
##' @author David LeBauer
##' @examples
##' \dontrun{
##'sensitivity.analysis(
##' sensitivity.analysis(
##' trait.samples[[pft$name]],
##' sa.samples[[pft$name]],
##' sa.agb[[pft$name]],
##' sa.samples[[pft$name]],
##' sa.agb[[pft$name]],
##' pft$outdir
##')
##' )
##' }
sensitivity.analysis <- function(trait.samples, sa.samples, sa.output, outdir) {
traits <- names(trait.samples)
sa.splines <- sapply(traits,
function(trait) sa.splinefun(sa.samples[[trait]], sa.output[[trait]]))

spline.estimates <- lapply(traits,
function(trait) spline.truncate(sa.splines[[trait]](trait.samples[[trait]])))
sa.splines <- sapply(
traits,
function(trait) sa.splinefun(sa.samples[[trait]], sa.output[[trait]])
)

spline.estimates <- lapply(
traits,
function(trait) spline.truncate(sa.splines[[trait]](trait.samples[[trait]]))
)
names(spline.estimates) <- traits
sensitivities <- sapply(traits,
function(trait) get.sensitivity(trait.samples[[trait]], sa.splines[[trait]]))
elasticities <- sapply(traits,
function(trait) get.elasticity(sensitivities[[trait]],
trait.samples[[trait]],
spline.estimates[[trait]]))
sensitivities <- sapply(
traits,
function(trait) get.sensitivity(trait.samples[[trait]], sa.splines[[trait]])
)
elasticities <- sapply(
traits,
function(trait) {
get.elasticity(
sensitivities[[trait]],
trait.samples[[trait]],
spline.estimates[[trait]]
)
}
)
variances <- sapply(traits, function(trait) stats::var(spline.estimates[[trait]]))
partial.variances <- variances / sum(variances)

coef.vars <- sapply(trait.samples, get.coef.var)
outlist <- list(sensitivity.output = list(sa.samples = sa.samples,
sa.splines = sa.splines),
variance.decomposition.output = list(coef.vars = coef.vars,
elasticities = elasticities,
sensitivities = sensitivities,
variances = variances,
partial.variances = partial.variances))
outlist <- list(
sensitivity.output = list(
sa.samples = sa.samples,
sa.splines = sa.splines
),
variance.decomposition.output = list(
coef.vars = coef.vars,
elasticities = elasticities,
sensitivities = sensitivities,
variances = variances,
partial.variances = partial.variances
)
)
return(outlist)
} # sensitivity.analysis


##' Truncate spline at zero if..
##' Truncate spline at zero if...
##'
##' Truncate spline at zero if P[x<0] < pnorm(-3)
##' pnorm(-3) chosen as default value for min quantile
##' because this is the default low end of range for the
##' sensitivity analysis.
##' This parameter could be determined based on minimum value in
##' settings$sensitivity.analysis$quantiles
##' @title Truncate spline
##'
##' @param x vector
##' @param min.quantile threshold quantile for testing lower bound on variable
##' @return either x or a vector with values < 0 converted to zero
##' @author David LeBauer
##' @export
##' @examples
##' set.seed(0)
##' x <- c(rgamma(998,1,1), rnorm(10))
##' x <- c(rgamma(998,1,1), rnorm(10))
##' min(x) # -0.5238
##' min(PEcAn.uncertainty::spline.truncate(x))
spline.truncate <- function(x, min.quantile = stats::pnorm(-3)) {
Expand Down
5 changes: 1 addition & 4 deletions modules/uncertainty/man/get.elasticity.Rd

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

9 changes: 3 additions & 6 deletions modules/uncertainty/man/get.sensitivity.Rd

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

8 changes: 2 additions & 6 deletions modules/uncertainty/man/kurtosis.Rd

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

9 changes: 3 additions & 6 deletions modules/uncertainty/man/sa.splinefun.Rd

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

5 changes: 1 addition & 4 deletions modules/uncertainty/man/sd.var.Rd

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

Loading

0 comments on commit 8abc90f

Please sign in to comment.