[CHNOSZ-commits] r378 - in pkg/CHNOSZ: . R demo inst man tests/testthat
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Feb 1 06:42:25 CET 2019
Author: jedick
Date: 2019-02-01 06:42:24 +0100 (Fri, 01 Feb 2019)
New Revision: 378
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/mosaic.R
pkg/CHNOSZ/demo/mosaic.R
pkg/CHNOSZ/inst/NEWS
pkg/CHNOSZ/man/mosaic.Rd
pkg/CHNOSZ/man/solubility.Rd
pkg/CHNOSZ/tests/testthat/test-mosaic.R
Log:
mosaic(): rewrite for performance and to handle more than two sets of changing basis species
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2019-02-01 01:45:00 UTC (rev 377)
+++ pkg/CHNOSZ/DESCRIPTION 2019-02-01 05:42:24 UTC (rev 378)
@@ -1,6 +1,6 @@
Date: 2019-02-01
Package: CHNOSZ
-Version: 1.1.3-85
+Version: 1.1.3-86
Title: Thermodynamic Calculations and Diagrams for Geochemistry
Authors at R: c(
person("Jeffrey", "Dick", , "j3ffdick at gmail.com", role = c("aut", "cre"),
Modified: pkg/CHNOSZ/R/mosaic.R
===================================================================
--- pkg/CHNOSZ/R/mosaic.R 2019-02-01 01:45:00 UTC (rev 377)
+++ pkg/CHNOSZ/R/mosaic.R 2019-02-01 05:42:24 UTC (rev 378)
@@ -1,15 +1,22 @@
# CHNOSZ/mosaic.R
# calculate affinities with changing basis species
-# 20141220 jmd
+# 20141220 jmd initial version
+# 20190129 complete rewrite to use any number of groups of changing basis species
+# and improve speed by pre-calculating subcrt values (sout)
+## if this file is interactively sourced, the following are also needed to provide unexported functions:
+#source("basis.R")
+#source("util.character.R")
+#source("util.args.R")
+
# function to calculate affinities with mosaic of basis species
-mosaic <- function(bases, bases2=NULL, blend=FALSE, ...) {
+mosaic <- function(bases, bases2 = NULL, blend = FALSE, mixing = FALSE, ...) {
# argument recall 20190120
# if the first argument is the result from a previous mosaic() calculation,
# just update the remaining arguments
if(is.list(bases)) {
- if(identical(bases[1], list(fun="mosaic"))) {
+ if(identical(bases[1], list(fun = "mosaic"))) {
aargs <- bases$args
# we can only update arguments given in ...
ddd <- list(...)
@@ -23,100 +30,131 @@
}
}
- if(is.null(bases2)) {
- # the arguments for affinity()
- myargs <- list(...)
- } else {
- # the arguments for affinity (first set of basis species; outer loop)
- myargs1 <- list(...)
- # the arguments for mosaic() (second set of basis species; inner loop)
- myargs <- list(bases=bases2, blend=blend, ...)
+ # backward compatibility 20190131:
+ # bases can be a vector instead of a list
+ # bases2 can be present
+ if(!is.list(bases)) {
+ bases <- list(bases)
+ hasbases2 <- FALSE
+ if(!is.null(bases2)) {
+ bases <- c(bases, list(bases2))
+ hasbases2 <- TRUE
+ }
+ otherargs <- list(...)
+ allargs <- c(list(bases = bases, blend = blend, mixing = mixing), otherargs)
+ out <- do.call(mosaic, allargs)
+ # replace A.bases (affinity calculations for all groups of basis species) with backwards-compatbile A.bases and A.bases2
+ if(hasbases2) A.bases2 <- out$A.bases[[2]]
+ A.bases <- out$A.bases[[1]]
+ out$A.bases <- A.bases
+ if(hasbases2) out <- c(out, list(A.bases2 = A.bases2))
+ return(out)
}
- # are the swapped basis species on the plot?
- # (the first one should be present in the starting basis set)
- iswap <- match(bases[1], names(myargs))
- # the log activity of the starting basis species
- logact.swap <- basis()$logact[ibasis(bases[1])]
+ # save starting basis and species definition
+ basis0 <- get("thermo")$basis
+ species0 <- get("thermo")$species
+ # get species indices of requested basis species
+ ispecies <- lapply(bases, info)
+ if(any(is.na(unlist(ispecies)))) stop("one or more of the requested basis species is unavailable")
+ # identify starting basis species
+ ispecies0 <- sapply(ispecies, "[", 1)
+ ibasis0 <- match(ispecies0, basis0$ispecies)
+ # quit if starting basis species are not present
+ ina <- is.na(ibasis0)
+ if(any(ina)) stop("the starting basis does not have ", paste(bases[ina], collapse = ", "))
- # a list where we'll keep the affinity calculations
- affs <- list()
- for(i in seq_along(bases)) {
- message(paste("mosaic: current basis species is", bases[i], sep=" "))
- # set up argument list: name of swapped-in basis species
- if(!is.na(iswap)) names(myargs)[iswap] <- bases[i]
- # calculate affinities
- if(is.null(bases2)) {
- affs[[i]] <- do.call(affinity, myargs)
- } else {
- mcall <- do.call(mosaic, myargs)
- affs[[i]] <- mcall$A.species
- A.bases2 <- mcall$A.bases
- }
- # change the basis species; restore the original at the end of the loop
- if(can.be.numeric(logact.swap)) logact.swap <- as.numeric(logact.swap)
- if(i < length(bases)) {
- swap.basis(bases[i], bases[i+1])
- # TODO: basis() requires the formula to identify the basis species,
- # would be nicer to just use the ibasis here
- bformula <- rownames(basis())[ibasis(bases[i+1])]
- basis(bformula, logact.swap)
- } else {
- swap.basis(bases[i], bases[1])
- bformula <- rownames(basis())[ibasis(bases[1])]
- basis(bformula, logact.swap)
- }
+ # run subcrt() calculations for all basis species and formed species 20190131
+ # this avoids repeating the calculations in different calls to affinity()
+ # add all the basis species here - the formed species are already present
+ lapply(bases, species)
+ sout <- affinity(..., return.sout = TRUE)
+
+ # calculate affinities of the basis species themselves
+ A.bases <- list()
+ for(i in 1:length(bases)) {
+ message("mosaic: calculating affinities of basis species group ", i, ": ", paste(bases[[i]], collapse=" "))
+ species(delete = TRUE)
+ species(bases[[i]])
+ A.bases[[i]] <- suppressMessages(affinity(..., sout = sout))
}
- # calculate affinities of formation of basis species
- message(paste("mosaic: combining diagrams for", paste(bases, collapse=" "), sep=" "))
- ispecies <- species()$ispecies
- species.logact <- species()$logact
- species(delete=TRUE)
- species(bases)
- if(is.null(bases2)) A.bases <- do.call(affinity, myargs)
- else A.bases <- do.call(affinity, myargs1)
- # restore original species with original activities
- species(delete=TRUE)
- species(ispecies, species.logact)
+ # get all combinations of basis species
+ newbases <- as.matrix(expand.grid(ispecies))
+ allbases <- matrix(basis0$ispecies, nrow = 1)[rep(1, nrow(newbases)), , drop = FALSE]
+ allbases[, ibasis0] <- newbases
- # affinities calculated using the first basis species
- A.species <- affs[[1]]
+ # calculate affinities of species for all combinations of basis species
+ aff.species <- list()
+ message("mosaic: calculating affinities of species for all ", nrow(allbases), " combinations of the basis species")
+ # run backwards so that we put the starting basis species back at the end
+ for(i in nrow(allbases):1) {
+ put.basis(allbases[i, ], basis0$logact)
+ # we have to define the species using the current basis
+ species(species0$ispecies, species0$logact)
+ aff.species[[i]] <- suppressMessages(affinity(..., sout = sout))
+ }
+
+ # calculate equilibrium mole fractions for each group of basis species
+ group.fraction <- list()
if(blend) {
- # calculate affinities using relative abundances of basis species
- # this isn't needed (and doesn't work) if all the affinities are NA 20180925
- if(any(!sapply(A.species$values, is.na))) {
- e <- equilibrate(A.bases)
- # what is the total activity of the basis species?
- a.tot <- Reduce("+", lapply(e$loga.equil, function(x) 10^x))
- for(j in seq_along(affs)) {
- for(i in seq_along(A.species$values)) {
- # start with zero affinity
- if(j==1) A.species$values[[i]][] <- 0
- # add affinity scaled by relative abundance of this basis species
- # and include mixing term (-x*log10(x)) 20190121
- x <- 10^e$loga.equil[[j]]/a.tot
- A.species$values[[i]] <- A.species$values[[i]] + affs[[j]]$values[[i]] * x - x * log10(x)
- }
+ for(i in 1:length(A.bases)) {
+ # this isn't needed (and doesn't work) if all the affinities are NA 20180925
+ if(any(!sapply(A.bases[[1]]$values, is.na))) {
+ e <- equilibrate(A.bases[[i]])
+ # exponentiate to get activities then divide by total activity
+ a.equil <- lapply(e$loga.equil, function(x) 10^x)
+ a.tot <- Reduce("+", a.equil)
+ group.fraction[[i]] <- lapply(a.equil, function(x) x / a.tot)
+ } else {
+ group.fraction[[i]] <- A.bases[[i]]$values
}
}
} else {
- # use affinities from the single predominant basis species
- d <- diagram(A.bases, plot.it=FALSE)
- # merge affinities using the second, third, ... basis species
- for(j in tail(seq_along(affs), -1)) {
- is.predominant <- d$predominant==j
- # diagram() produces NA beyond water limits on Eh-pH diagrams (but we can't use NA for indexing, below)
- is.predominant[is.na(is.predominant)] <- FALSE
- for(i in seq_along(A.species$values)) {
- A.species$values[[i]][is.predominant] <- affs[[j]]$values[[i]][is.predominant]
+ # for blend = FALSE, we just look at whether
+ # a basis species predominates within its group
+ for(i in 1:length(A.bases)) {
+ d <- diagram(A.bases[[i]], plot.it = FALSE, limit.water = FALSE)
+ group.fraction[[i]] <- list()
+ for(j in 1:length(bases[[i]])) {
+ # if a basis species predominates, it has a mole fraction of 1, or 0 otherwise
+ yesno <- d$predominant
+ yesno[yesno != j] <- 0
+ yesno[yesno == j] <- 1
+ group.fraction[[i]][[j]] <- yesno
}
}
}
+ # make an indexing matrix for all combinations of basis species
+ ind.mat <- list()
+ for(i in 1:length(ispecies)) ind.mat[[i]] <- 1:length(ispecies[[i]])
+ ind.mat <- as.matrix(expand.grid(ind.mat))
+
+ # calculate mole fractions for each combination of basis species
+ for(i in 1:nrow(ind.mat)) {
+ # multiply fractions from each group
+ for(j in 1:ncol(ind.mat)) {
+ if(j==1) x <- group.fraction[[j]][[ind.mat[i, j]]]
+ else x <- x * group.fraction[[j]][[ind.mat[i, j]]]
+ }
+ # multiply affinities by the mole fractions of basis species
+ # include mixing term (-x*log10(x)) 20190121
+ if(blend & mixing) aff.species[[i]]$values <- lapply(aff.species[[i]]$values, function(values) values * x - x * log10(x))
+ else aff.species[[i]]$values <- lapply(aff.species[[i]]$values, function(values) values * x)
+ }
+
+ # get total affinities for the species
+ A.species <- aff.species[[1]]
+ for(i in 1:length(A.species$values)) {
+ # extract the affinity contributions from each basis species
+ A.values <- lapply(lapply(aff.species, "[[", "values"), "[[", i)
+ # sum them to get total affinities for this species
+ A.species$values[[i]] <- Reduce("+", A.values)
+ }
+
# for argument recall, include all arguments in output 20190120
- allargs <- c(list(bases=bases, bases2=bases2, blend=blend), list(...))
+ allargs <- c(list(bases = bases, blend = blend, mixing = mixing), list(...))
# return the affinities for the species and basis species
- if(is.null(bases2)) return(list(fun="mosaic", args=allargs, A.species=A.species, A.bases=A.bases))
- else return(list(fun="mosaic", args=allargs, A.species=A.species, A.bases=A.bases, A.bases2=A.bases2))
+ return(list(fun = "mosaic", args = allargs, A.species = A.species, A.bases = A.bases))
}
Modified: pkg/CHNOSZ/demo/mosaic.R
===================================================================
--- pkg/CHNOSZ/demo/mosaic.R 2019-02-01 01:45:00 UTC (rev 377)
+++ pkg/CHNOSZ/demo/mosaic.R 2019-02-01 05:42:24 UTC (rev 378)
@@ -1,10 +1,13 @@
+# CHNOSZ/demo/mosaic.R
+# 20141221 first version
+
# Fe-minerals and aqueous species in Fe-S-O-H-C system
# after Garrels and Christ, 1965 Figure 7.21
# to reproduce their diagram as closely as posssible, use their thermodynamic data (from Appendix 2)
-mod.obigt(c("Fe+2", "Fe+3"), G=c(-20300, -2520))
-mod.obigt(c("hematite", "magnetite", "pyrrhotite", "pyrite", "siderite"), G=c(-177100, -242400, -23320, -36000, -161060))
-mod.obigt(c("SO4-2", "HS-", "H2S", "HSO4-"), G=c(-177340, 3010, -6540, -179940))
-mod.obigt(c("CO2", "HCO3-", "CO3-2"), G=c(-92310, -140310, -126220))
+mod.obigt(c("Fe+2", "Fe+3"), G = c(-20300, -2520))
+mod.obigt(c("hematite", "magnetite", "pyrrhotite", "pyrite", "siderite"), G = c(-177100, -242400, -23320, -36000, -161060))
+mod.obigt(c("SO4-2", "HS-", "H2S", "HSO4-"), G = c(-177340, 3010, -6540, -179940))
+mod.obigt(c("CO2", "HCO3-", "CO3-2"), G = c(-92310, -140310, -126220))
# conditions and system definition
pH <- c(0, 14, 400)
Eh <- c(-1, 1, 400)
@@ -22,18 +25,18 @@
# calculate affinities using the predominant basis species
# using blend=TRUE we get curvy lines, particularly at the boundaries with siderite
# compare with the plot in Garrels and Christ, 1965
-m1 <- mosaic(bases, bases2, blend=TRUE, pH=pH, Eh=Eh, T=T)
+m1 <- mosaic(bases, bases2, blend = TRUE, pH = pH, Eh = Eh, T = T)
# make a diagram and add water stability lines
-diagram(m1$A.species, lwd=2)
-water.lines(m1$A.species, col="seagreen", lwd=1.5)
+diagram(m1$A.species, lwd = 2)
+water.lines(m1$A.species, col = "seagreen", lwd = 1.5)
# show lines for Fe(aq) = 10^-4 M
species(c("Fe+2", "Fe+3"), -4)
-m2 <- mosaic(bases, bases2, blend=TRUE, pH=pH, Eh=Eh, T=T)
-diagram(m2$A.species, add=TRUE, names=NULL)
+m2 <- mosaic(bases, bases2, blend = TRUE, pH = pH, Eh = Eh, T = T)
+diagram(m2$A.species, add = TRUE, names = NULL)
title(main=paste("Iron oxides, sulfides and carbonate in water, log(total S) = -6,",
- "log(total C)=0, after Garrels and Christ, 1965", sep="\n"))
+ "log(total C)=0, after Garrels and Christ, 1965", sep = "\n"))
# overlay the carbonate basis species predominance fields
-d <- diagram(m1$A.bases2, add=TRUE, col="blue", names=NULL, lty=3, limit.water=FALSE)
-text(d$namesx, -0.8, as.expression(sapply(m1$A.bases2$species$name, expr.species)), col="blue")
+d <- diagram(m1$A.bases2, add = TRUE, col = "blue", names = NULL, lty = 3, limit.water = FALSE)
+text(d$namesx, -0.8, as.expression(sapply(m1$A.bases2$species$name, expr.species)), col = "blue")
# reset the database, as it was changed in this example
data(thermo)
Modified: pkg/CHNOSZ/inst/NEWS
===================================================================
--- pkg/CHNOSZ/inst/NEWS 2019-02-01 01:45:00 UTC (rev 377)
+++ pkg/CHNOSZ/inst/NEWS 2019-02-01 05:42:24 UTC (rev 378)
@@ -1,4 +1,4 @@
-CHANGES IN CHNOSZ 1.1.3-85 (2019-02-01)
+CHANGES IN CHNOSZ 1.1.3-86 (2019-02-01)
---------------------------------------
BUG FIXES
@@ -69,6 +69,15 @@
of NaCl in water, taking account of activity coefficients and the
reaction Na+ + Cl- = NaCl(aq).
+OTHER NEW FEATURES
+
+- Add dumpdata() for returning/writing all packaged thermodynamic data
+ (including default database and optional data files). The file is
+ also available on the website (chnosz.net/download/alldata.csv).
+
+- mosaic() has been rewritten to handle more than two changing groups
+ of basis species.
+
DOCUMENTATION
- Add demo/gold.R for calculations of Au solubility in hydrothermal
@@ -92,10 +101,6 @@
THERMODYNAMIC DATA
-- Add dumpdata() for returning/writing all packaged thermodynamic data
- (including default database and optional data files). The file is
- also available on the website (chnosz.net/download/alldata.csv).
-
- The Berman data (Berman, 1988 and later additions) have replaced the
SUPCRT92 data (based on Helgeson et al., 1978) for most minerals in
the default database (i.e. the one loaded by data(thermo)). Only
Modified: pkg/CHNOSZ/man/mosaic.Rd
===================================================================
--- pkg/CHNOSZ/man/mosaic.Rd 2019-02-01 01:45:00 UTC (rev 377)
+++ pkg/CHNOSZ/man/mosaic.Rd 2019-02-01 05:42:24 UTC (rev 378)
@@ -7,13 +7,14 @@
}
\usage{
- mosaic(bases, bases2=NULL, blend=FALSE, ...)
+ mosaic(bases, bases2 = NULL, blend = FALSE, mixing = FALSE, ...)
}
\arguments{
- \item{bases}{character, basis species to be changed in the calculation}
+ \item{bases}{character, basis species to be changed in the calculation, or list, containing vectors for each group of changing basis species}
\item{bases2}{character, second set of changing basis species}
\item{blend}{logical, use relative abundances of basis species?}
+ \item{mixing}{logical, include a term for the Gibbs energy of mixing?}
\item{...}{additional arguments to be passed to \code{\link{affinity}}}
}
@@ -28,14 +29,21 @@
The first species listed in \code{bases} should be in the current basis definition.
The arguments in \code{...} are passed to \code{affinity} to specify the conditions.
If \code{blend} is FALSE (the default), the function returns the affinities calculated using the single predominant basis species in \code{bases} at each condition.
-If \code{blend} is TRUE, the function combines the affinities of the formation reactions in proportion to the relative abundances of the basis species at each condition, including a term to account for the Gibbs energy of mixing.
-See the second example in \code{\link{solubility}} for a numerical test of the calculations using \code{blend}.
+If \code{blend} is TRUE, the function combines the affinities of the formation reactions in proportion to the relative abundances of the basis species at each condition.
+Additionally, if \code{mixing} is TRUE, a term is included to account for the Gibbs energy of mixing.
+See the second example in \code{\link{solubility}} for a numerical test of the calculations using \code{blend} and \code{mixing}.
The basis species listed in \code{bases} should all be related to the first basis species there (i.e. all share the same element).
A second, independent set of basis species can be provided in \code{bases2} (for example \samp{CO3-2}, \samp{HCO3-}, \samp{CO2}, if the first set of basis species are the sulfur-bearing ones listed above).
The function then works recursively, by calling itself instead of \code{affinity}, so that the inner loop changes the basis species in \code{bases2}.
In this way, all possible combinations of the two sets of basis species are used in the calculation.
+A more flexible method of specifying multiple sets of basis species is now available.
+Instead of using \code{bases} and \code{bases2}, supply a list for just the \code{bases} argument.
+The list should contain any number of vectors specifying the groups of basis species.
+All combinations of basis species in these groups are used for the calculations.
+This overcomes the prior limitation of only having two changing groups of basis species.
+
}
\value{
Modified: pkg/CHNOSZ/man/solubility.Rd
===================================================================
--- pkg/CHNOSZ/man/solubility.Rd 2019-02-01 01:45:00 UTC (rev 377)
+++ pkg/CHNOSZ/man/solubility.Rd 2019-02-01 05:42:24 UTC (rev 378)
@@ -121,7 +121,7 @@
## method 2: CO2 and carbonate species as basis species
basis(c("calcite", "CO2", "H2O", "O2", "H+"))
species(c("Ca+2"))
-m <- mosaic(c("CO2", "HCO3-", "CO3-2"), pH = c(3, 14), blend = TRUE)
+m <- mosaic(c("CO2", "HCO3-", "CO3-2"), pH = c(3, 14), blend = TRUE, mixing = TRUE)
sm0 <- solubility(m)
smI <- solubility(m, find.IS = TRUE)
## plot the results
Modified: pkg/CHNOSZ/tests/testthat/test-mosaic.R
===================================================================
--- pkg/CHNOSZ/tests/testthat/test-mosaic.R 2019-02-01 01:45:00 UTC (rev 377)
+++ pkg/CHNOSZ/tests/testthat/test-mosaic.R 2019-02-01 05:42:24 UTC (rev 378)
@@ -34,8 +34,9 @@
bases <- c("SO4-2", "HSO4-", "HS-", "H2S")
# calculate affinities using the predominant basis species
pH <- c(0, 14, 29)
- m1 <- mosaic(bases, pH=pH)
- m2 <- mosaic(bases, pH=pH, blend=TRUE)
+ m1 <- mosaic(bases, pH = pH)
+ # calculate affinities with smooth transitions between basis species, including a mixing energy
+ m2 <- mosaic(bases, pH = pH, blend = TRUE, mixing = TRUE)
# these species have no S so the results should be similar,
# 20190121 except for a negative free energy of mixing (positive affinity)
expect_true(all(m2$A.species$values[[1]] - m1$A.species$values[[1]] > 0))
More information about the CHNOSZ-commits
mailing list