Select Git revision
-
Daniel Scheffler authoredRevised tests. Added parameterized as test requirement. Replaced deprecated numpy data types with builtin types. Fixed dead link in the docs. Signed-off-by:Daniel Scheffler <danschef@gfz-potsdam.de>
Daniel Scheffler authoredRevised tests. Added parameterized as test requirement. Replaced deprecated numpy data types with builtin types. Fixed dead link in the docs. Signed-off-by:
Rphree_SurrogateUtils.R 10.96 KiB
## Functions for dealing with surrogate simulations
### Marco De Lucia, delucia@gfz-potsdam.de, 2009-2018
### Time-stamp: "Last modified 2018-05-06 20:16:47 delucia"
##' Computes the average of absolute values of a vector
##' @title Average of absolute values
##' @param x numeric vector
##' @return mean(abs(x))
##' @author MDL
##' @export
mae <- function(x) mean(abs(x))
## ##' .. content for \description{} (no empty lines) ..
## ##'
## ##' .. content for \details{} ..
## ##' @title SoftMax transformation
## ##' @param obj numeric vector
## ##' @return
## ##' @author
## SoftMax <- function(obj) {
## eo <- exp(obj-max(obj))
## ss <- sum(eo)
## scaled <- eo/ss
## attr(res,"back") <- ss
## return(res)
## }
## ## SoftMaxBackTransf <- function(vec, back) return(log(vec)+log(back))
## ##' @title Flatten out a Rphree list
## ##' @param list list of Rphree lists
## ##' @param strip logical, should we delete "ListInfo"? Defaults to
## ##' true
## ##' @return a list whose elements are single Rphree simulations
## ##' @author MDL
## ##' @export
## FlattenList <- function(list) {
## tot <- vector(mode="list", length=sum(sapply(stripped, length)))
## k <- 1
## for (i in seq_along(stripped)) {
## for (j in seq_along(stripped[[i]])) {
## tot[[k]] <- stripped[[i]][[j]]
## k=k+1
## }
## }
## return(tot)
## }
##' @title Find all species occurring in the ensemble
##' @param flatlist the flatted list of Rphree simulations
##' @return vector with unique occurrences of species names
##' @author MDL
##' @export
FindAllSpeciesNames <- function(flatlist)
unique(sort(unlist(sapply(flatlist, function(sol) rownames(sol$species)))))
##' @title Find all pphases occurring in the ensemble
##' @param flatlist the flatted list of Rphree simulations
##' @return vector with unique occurrences of pphases names
##' @author MDL
##' @export
FindAllMinNames <- function(flatlist)
unique(sort(unlist(sapply(flatlist, function(sol) rownames(sol$pphases)))))
##' @title Find all totals occurring in the ensemble
##' @param flatlist the flatted list of Rphree simulations
##' @return vector with unique occurrences of element names##' @author MDL
##' @export
FindAllTotNames <- function(flatlist)
unique(sort(unlist(sapply(flatlist, function(sol) rownames(sol$tot)))))
##' @title Extract species from a flatted Rphree list
##' @param flatlist the list as input to search in
##' @param species vector containing the names of the species to
##' search for
##' @return matix with one species per column
##' @author MDL
##' @export
ExtractSpecies <- function(flatlist, species) {
tmp <- matrix(NA, len <- length(flatlist), speclen <- length(species))
colnames(tmp) <- species
for (i in seq(speclen)) {
tmp[,i] <- RPinfo(flatlist, "species", species[i], flex=TRUE)
}
return(tmp)
}
##' @title Extract all phases from a flatted Rphree list
##' @param flatlist the list as input to search in
##' @param totals vector containing the names of the phases to search
##' for
##' @return matix with one pphases per column
##' @author MDL
##' @export
ExtractPphases <- function(flatlist, pphases) {
tmp <- matrix(NA, len <- length(flatlist), speclen <- length(pphases))
colnames(tmp) <- pphases
for (i in seq(speclen)) {
tmp[,i] <- RPinfo(flatlist, "pphase", pphases[i], flex=TRUE)
}
return(tmp)
}
##' @title Extract total element concentrations from a flatted Rphree
##' list
##' @param flatlist the list as input to search in
##' @param totals vector containing the names of the elements to
##' search for
##' @return Matrix with one element per column
##' @author MDL
##' @export
ExtractTotals <- function(flatlist, totals) {
tmp <- matrix(NA, len <- length(flatlist), speclen <- length(totals))
colnames(tmp) <- totals
for (i in seq(speclen)) {
tmp[,i] <- RPinfo(flatlist, "tot", totals[i], flex=TRUE)
}
return(tmp)
}
##' @title Scan a list of Rphree simulations to form stoichiometric
##' matrix
##' @param flatlist a flatted list of Rphree simulations
##' @param db database (as read by RReadFile(..., is.db=TRUE))
##' @param species vector containing all species to be included in the
##' stoichiometric matrix. If omitted, the flat list is searched
##' for all possible species
##' @param pphases vector containing all EQUILIBRIUM minerals to be
##' included in the stoichiometric matrix. If omitted, the flat
##' list is searched for all possible minerals.
##' @return A stoichiometric matrix
##' @author MDL
##' @export
StoichiometricMatrix <- function(flatlist, db, species, pphases){
if (missing(species))
species <- FindAllSpeciesNames(flatlist)
if(missing(pphases))
minerals <- FindAllMinNames(flatlist)
pos_solsp <- grep("SOLUTION_SPECIES", db, fixed=TRUE)
pos_phases <- grep("PHASES", db, fixed=TRUE)
dbtmp <- db[{pos_solsp+1}:{pos_phases -1}]
tots <- lapply(species, function(x) grep(x,dbtmp,fixed=TRUE))
equations <- dbtmp[grep("^[[:alnum:]]", dbtmp, perl=TRUE)]
parsed <- lapply(equations, ParseFormula)
names(parsed) <- paste("e",seq_along(equations), sep="")
matched <- lapply(parsed, function(x) match(x$comp, species))
matched2 <- sapply(matched, function(x) NA %in% x)
indeq <- which(!matched2)
equations <- equations[indeq]
matched <- matched[indeq]
parsed <- parsed[indeq]
indeq <- which(sapply(matched, length) > 2)
equations <- equations[indeq]
matched <- matched[indeq]
parsed <- parsed[indeq]
parsedmin <- lapply(minerals, FindLogK,db=db, type="rhs")
names(parsedmin) <- minerals
matchedmin <- lapply(parsedmin, function(x) match(x$comp, species))
totparsed <- c(parsed, parsedmin)
totmatch <- c(matched, matchedmin)
mat <- matrix(0, ncol=length(species)+length(parsedmin), nrow=length(totmatch))
colnames(mat) <- c(species,minerals)
rownames(mat) <- names(totmatch)
for (i in seq(nrow(mat))) {
cat(c("Equation: ",i, paste( totparsed[[i]]$comp, collapse=" ")),"\n")
mat[i,totmatch[[i]]] <- totparsed[[i]]$coeff
}
return(list(mat=mat, parsed=totparsed, matched=totmatch, minerals=minerals, species=species))
}
##' @title Returns the composition of a mineral in terms of chemical
##' elements (not of species)
##' @param pphase name of the mineral to parse
##' @param stoichmat full stoichiometric matrix returned by
##' @param elements elements to include in the mass balance equations,
##' defaults to c("C","Ca","Mg","Cl")
##' @return named vector with stoichiometric coefficients of chemical
##' elements for the pphase mineral
##' @author MDL
##' @export
ElementalBalanceMin <- function(pphase, stoichmat, elements = c("C","Ca","Mg","Cl"))
{
ind <- match(pphase[1], rownames(stoichmat))
eq <- stoichmat[ ind[1],]
eq <- eq[eq!=0]
Coefs <- rep(0, length(elements))
names(Coefs) <- elements
for (i in seq_along(elements)) {
Shift <- nchar(elements[i])-1
expr <- paste0(elements[i],'[A-Z0-9\\+\\-\\ ]') ## Element followed by UpperCase, number, "+" space or "-"
indm <- grep(expr, names(eq), perl=TRUE)
if (length(indm)!=0) {
## Now find the coefficient
SpeciesHavingElement <- names(eq)[indm] Multip <- eq[indm]
Pos <- regexpr(expr, SpeciesHavingElement)
NextChar <- substr(SpeciesHavingElement, Pos+1+Shift, Pos+1+Shift)
if (NextChar == "+")
Coefs[i] <- Multip
if (NextChar == "-")
Coefs[i] <- Multip
if (length(grep('[A-Z]', NextChar))>0)
Coefs[i] <- Multip
if (length(grep('[0-9]', NextChar))>0)
Coefs[i] <- as.integer(NextChar)*Multip
}
}
Coefs <- Coefs[Coefs!=0]
return(Coefs)
}
##' @title Form the balance equations from a full stoichiometric
##' matrix
##' @param stmat a stoichiometric matrix returned by
##' @param elements elements to include in the mass balance equations,
##' defaults to c("C","Ca","Mg","Cl")
##' @return a structure (list) containing the balance equations
##' @author MDL
##' @export
BalanceEquations <- function(stmat, elements = c("C","Ca","Mg","Cl"))
{
balance <- lapply(stmat$minerals, ElementalBalanceMin, stoichmat=stmat$mat, elements = elements)
names(balance) <- stmat$minerals
return(balance)
}
##' @title Figures out the formula from a baleq structure
##' @param baleq the data structure expressing the mass balance of
##' minerals
##' @return a named vector with the stoichiometric formula of the
##' minerals
##' @author MDL
##' @export
FormulaFromBal <- function(baleq)
{
sapply(names(baleq), function(x) paste0(baleq[[x]],"*", names(baleq[[x]]), collapse="+" ))
}
##' @title Check mass balance between design and results for
##' surrogates
##' @param baleq structure (list) containing elemental mass balance
##' @param des design (state_T under our conventions)
##' @param res results of the surrogates (corresponds to state_C)
##' @param elements optional, vector of elements on which to perform
##' mass balance check
##' @return a matrix containing a mass balance error, one column for
##' each tested equation (=element)
##' @author MDL
##' @export
CheckBalance <- function(baleq, des, res, elements)
{
pphases <- names(baleq)
desm <- data.matrix(des)
resm <- data.matrix(res)
## deltas <- data.matrix(as.data.frame(lapply(pphases, function(x) resm[,x] - desm[,x]),optional = TRUE))
## colnames(deltas) <- paste0("D",pphases)
if (missing(elements))
elements <- unique(sort(unlist(lapply(baleq, function(x) names(x)))))
diffs <- matrix(NA, ncol=length(elements), nrow=nrow(desm)) colnames(diffs) <- elements
for (i in seq_along(elements)){
## lhs == desmign
## matchinmin <- sapply(baleq, function(x) match(elements[i], names(x)))
## matchinmin <- matchinmin[!is.na(matchinmin)]
coefs <- sapply(baleq, function(x) return(x[elements[i]]))
names(coefs) <- names(baleq)
indcoefs <- which(!is.na(coefs))
coefs <- coefs[indcoefs]
minleft <- sapply(names(coefs), function(x) desm[,x]*coefs[x])
lhs <- desm[,elements[i]] + rowSums(minleft)
minright <- sapply(names(coefs), function(x) resm[,x]*coefs[x])
rhs <- resm[,elements[i]] + rowSums(minright)
diffs[,i] <- rhs - lhs
}
return(diffs)
}
##' From a list containing results of a Reactive Transport simulation,
##' this function extracts the "transport" as a "perturbation" of
##' chemistry. This simply means computing the differences between the
##' state_T at current time and state_C at previous time step.
##' @title Extract the changes in stateC due to transport
##' @param sim a result of Reactive Transport (for each time step, a T
##' and a C data.frame)
##' @return a list
##' @author MDL
##' @export
TranspAsPert <- function (sim) {
out <- lapply(2:length(sim), function(x) sim[[x]]$T - sim[[x-1]]$C)
names(out) <- names(sim)[2:length(sim)]
return(out)
}