[CHNOSZ-commits] r307 - in pkg/CHNOSZ: . R demo inst man tests/testthat vignettes
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Mar 14 18:54:46 CET 2018
Author: jedick
Date: 2018-03-14 18:54:46 +0100 (Wed, 14 Mar 2018)
New Revision: 307
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/diagram.R
pkg/CHNOSZ/demo/buffer.R
pkg/CHNOSZ/inst/NEWS
pkg/CHNOSZ/man/buffer.Rd
pkg/CHNOSZ/man/diagram.Rd
pkg/CHNOSZ/tests/testthat/test-diagram.R
pkg/CHNOSZ/vignettes/anintro.Rmd
Log:
diagram(): allow 'saturation' to be specified for 'type' of diagram
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2018-03-07 15:46:30 UTC (rev 306)
+++ pkg/CHNOSZ/DESCRIPTION 2018-03-14 17:54:46 UTC (rev 307)
@@ -1,6 +1,6 @@
-Date: 2018-03-07
+Date: 2018-03-15
Package: CHNOSZ
-Version: 1.1.3-14
+Version: 1.1.3-15
Title: Thermodynamic Calculations for Geobiochemistry
Authors at R: c(
person("Jeffrey", "Dick", , "j3ffdick at gmail.com", role = c("aut", "cre"),
Modified: pkg/CHNOSZ/R/diagram.R
===================================================================
--- pkg/CHNOSZ/R/diagram.R 2018-03-07 15:46:30 UTC (rev 306)
+++ pkg/CHNOSZ/R/diagram.R 2018-03-14 17:54:46 UTC (rev 307)
@@ -14,8 +14,8 @@
diagram <- function(
# species affinities or activities
eout,
- # what to plot
- what="loga.equil", alpha=FALSE, normalize=FALSE, as.residue=FALSE, balance=NULL,
+ # type of plot
+ type="auto", alpha=FALSE, normalize=FALSE, as.residue=FALSE, balance=NULL,
groups=as.list(1:length(eout$values)), xrange=NULL,
# figure size and sides for axis tick marks
mar=NULL, yline=par("mgp")[1]+0.3, side=1:4,
@@ -24,7 +24,7 @@
# sizes
cex=par("cex"), cex.names=1, cex.axis=par("cex"),
# line styles
- lty=NULL, lwd=par("lwd"), dotted=NULL, spline.method=NULL,
+ lty=NULL, lwd=par("lwd"), dotted=NULL, spline.method=NULL, contour.method="edge",
# colors
col=par("col"), col.names=par("col"), fill=NULL,
fill.NA="slategray1", limit.water=TRUE,
@@ -41,25 +41,27 @@
## check that eout is valid input
if(!"sout" %in% names(eout)) stop("'eout' does not look like output from equil() or affinity()")
- ## 'what' can be:
- # loga.equil - equilibrium activities of species of interest (eout)
- # basis species - equilibrium activity of a basis species (aout)
- # missing - property from affinity() or predominances of species (aout)
+ ## 'type' can be:
+ # 'auto' - property from affinity() (1D) or maximum affinity (affinity 2D) (aout) or loga.equil (eout)
+ # 'loga.equil' - equilibrium activities of species of interest (eout)
+ # name of basis species - equilibrium activity of a basis species (aout)
+ # 'saturation' - affinity=0 line for each species (2D)
eout.is.aout <- FALSE
plot.loga.basis <- FALSE
- if(missing(what)) {
+ if(type %in% c("auto", "saturation")) {
if(!"loga.equil" %in% names(eout)) {
eout.is.aout <- TRUE
# get the balancing coefficients
- n.balance <- balance(eout, balance)
+ if(type=="auto") n.balance <- balance(eout, balance)
+ else n.balance <- rep(1, length(eout$values))
}
- } else if(what %in% rownames(eout$basis)) {
+ } else if(type %in% rownames(eout$basis)) {
# to calculate the loga of basis species at equilibrium
if(!missing(groups)) stop("can't plot equilibrium activities of basis species for grouped species")
if(isTRUE(alpha) | is.character(alpha)) stop("equilibrium activities of basis species not available with alpha=TRUE")
plot.loga.basis <- TRUE
- } else if(what=="loga.equil" & !"loga.equil" %in% names(eout)) stop("'eout' is not the output from equil()")
- else if(what!="loga.equil") stop(what, " is not a basis species or 'loga.equil'")
+ } else if(type=="loga.equil" & !"loga.equil" %in% names(eout)) stop("'eout' is not the output from equil()")
+ else if(type!="loga.equil") stop(type, " is not a valid diagram type")
## consider a different number of species if we're grouping them together
ngroups <- length(groups)
@@ -69,6 +71,11 @@
plotvals <- eout$loga.equil
plotvar <- "loga.equil"
+ ## number of dimensions (T, P or chemical potentials that are varied)
+ # length(eout$vars) - the number of variables = the maximum number of dimensions
+ # length(dim(eout$values[[1]])) - nd=1 if it was a transect along multiple variables
+ nd <- min(length(eout$vars), length(dim(eout$values[[1]])))
+
## deal with output from affinity()
if(eout.is.aout) {
# plot property from affinity(), divided by balancing coefficients
@@ -83,15 +90,12 @@
# 20171027 use parentheses to avoid ambiguity about order of operations
if(plotvar=="A") {
plotvar <- "A/(2.303RT)"
- message("diagram: plotting A/(2.303RT) / n.balance (maximum affinity method for 2-D diagrams)")
+ if(nd==2 & type=="auto") message("diagram: using maximum affinity method for 2-D diagram")
+ else if(nd==2 & type=="saturation") message("diagram: plotting saturation lines for 2-D diagram")
+ else message("diagram: plotting A/(2.303RT) / n.balance")
} else message(paste("diagram: plotting", plotvar, " / n.balance"))
}
- ## number of dimensions (T, P or chemical potentials that are varied)
- # length(eout$vars) - the number of variables = the maximum number of dimensions
- # length(dim(eout$values[[1]])) - nd=1 if it was a transect along multiple variables
- nd <- min(length(eout$vars), length(dim(eout$values[[1]])))
-
## use molality instead of activity if the affinity calculation include ionic strength 20171101
use.molality <- "IS" %in% names(eout)
@@ -125,10 +129,10 @@
## calculate the equilibrium logarithm of activity of a basis species
## (such that affinities of formation reactions are zero)
if(plot.loga.basis) {
- ibasis <- match(what, rownames(eout$basis))
+ ibasis <- match(type, rownames(eout$basis))
# the logarithm of activity used in the affinity calculation
is.loga.basis <- can.be.numeric(eout$basis$logact[ibasis])
- if(!is.loga.basis) stop(paste("the logarithm of activity for basis species", what, "is not numeric - was a buffer selected?"))
+ if(!is.loga.basis) stop(paste("the logarithm of activity for basis species", type, "is not numeric - was a buffer selected?"))
loga.basis <- as.numeric(eout$basis$logact[ibasis])
# the reaction coefficients for this basis species
nu.basis <- eout$species[, ibasis]
@@ -137,7 +141,7 @@
# eout$values is a strange name for affinity ... should be named something like eout$affinity ...
loga.basis - eout$values[[x]]/nu.basis[x]
})
- plotvar <- what
+ plotvar <- type
}
## alpha: plot fractional degree of formation
@@ -157,7 +161,7 @@
## identify predominant species
predominant <- NA
- if(plotvar %in% c("loga.equil", "alpha", "A/(2.303RT)")) {
+ if(plotvar %in% c("loga.equil", "alpha", "A/(2.303RT)") & type!="saturation") {
pv <- plotvals
# some additional steps for affinity values, but not for equilibrated activities
if(eout.is.aout) {
@@ -197,9 +201,6 @@
}
}
- # a warning about that we can only show properties of the first species on a 2-D diagram
- if(nd==2 & length(plotvals) > 1 & identical(predominant, NA)) warning("showing only first species in 2-D property diagram")
-
## where we'll put extra output for predominance diagrams (lx, ly, is)
out2D <- list()
@@ -534,8 +535,31 @@
}
# colors and curves (predominance), or contours (properties)
if(identical(predominant, NA)) {
- zs <- plotvals[[1]]
- contour(xs, ys, zs, add=TRUE, col=col, lty=lty, lwd=lwd, labcex=cex)
+ if(type=="saturation") {
+ # for saturation plot, contour affinity=0 for all species
+ for(i in 1:length(plotvals)) {
+ zs <- plotvals[[i]]
+ # skip plotting if this species has no possible saturation line, or a line outside the plot range
+ if(length(unique(as.numeric(zs)))==1) {
+ message("diagram: no saturation line possible for ", names[i])
+ next
+ }
+ if(all(zs < 0) | all(zs > 0)) {
+ message("diagram: beyond range for saturation line of ", names[i])
+ next
+ }
+ if(identical(contour.method, NULL) | identical(contour.method, NA) | identical(contour.method, ""))
+ contour(xs, ys, zs, add=TRUE, col=col, lty=lty, lwd=lwd, labcex=cex, levels=0, labels=names[i], drawlabels=FALSE)
+ else contour(xs, ys, zs, add=TRUE, col=col, lty=lty, lwd=lwd, labcex=cex, levels=0, labels=names[i], method=contour.method)
+ }
+ } else {
+ # otherwise, make contours of properties using first species only
+ if(length(plotvals) > 1) warning("showing only first species in 2-D property diagram")
+ print('hello')
+ print(length(plotvals))
+ zs <- plotvals[[1]]
+ contour(xs, ys, zs, add=TRUE, col=col, lty=lty, lwd=lwd, labcex=cex, method=contour.method)
+ }
pn <- list(lx=NULL, ly=NULL, is=NULL)
} else {
# put predominance matrix in the right order for image() etc
@@ -551,6 +575,14 @@
} # end if(nd==2)
} # end if(plot.it)
+ # warn if we have a system with both Berman and Helgeson minerals 20180315
+ ref1 <- get("thermo")$obigt$ref1
+ ref2 <- get("thermo")$obigt$ref2
+ ispecies <- eout$species$ispecies
+ hasHelgeson <- any(grepl("HDNB78", ref1[ispecies])) | any(grepl("HDNB78", ref2[ispecies]))
+ hasBerman <- any(get("thermo")$obigt$state[ispecies]=="cr_Berman")
+ if(hasHelgeson & hasBerman) warning("the system has minerals from both the Helgeson and Berman datasets; data may not be internally consistent")
+
out <- c(eout, list(plotvar=plotvar, plotvals=plotvals, names=names, predominant=predominant), out2D)
return(invisible(out))
}
Modified: pkg/CHNOSZ/demo/buffer.R
===================================================================
--- pkg/CHNOSZ/demo/buffer.R 2018-03-07 15:46:30 UTC (rev 306)
+++ pkg/CHNOSZ/demo/buffer.R 2018-03-14 17:54:46 UTC (rev 307)
@@ -19,12 +19,12 @@
bufferline("QFM", 38)
bufferline("PPM", 102)
bufferline("HM", 51)
-# method 2: in diagram(), use the `what` argument
+# method 2: in diagram(), use the `type` argument
basis("H2", 0)
for(logact in c(-6, -10, -15)) {
species(c("formaldehyde", "HCN"), logact)
a <- affinity(T=xlim, P=300)
- d <- diagram(a, what="H2", lty=c(3, 2), add=TRUE)
+ d <- diagram(a, type="H2", lty=c(3, 2), add=TRUE)
text(a$vals[[1]][13], mean(sapply(d$plotvals, c)[13, ]), logact)
}
# add legends and title
Modified: pkg/CHNOSZ/inst/NEWS
===================================================================
--- pkg/CHNOSZ/inst/NEWS 2018-03-07 15:46:30 UTC (rev 306)
+++ pkg/CHNOSZ/inst/NEWS 2018-03-14 17:54:46 UTC (rev 307)
@@ -1,4 +1,4 @@
-CHANGES IN CHNOSZ 1.1.3-14 (2018-03-07)
+CHANGES IN CHNOSZ 1.1.3-15 (2018-03-15)
---------------------------------------
- Lines in 1-D diagram()s can optionally be drawn as splines using the
@@ -37,6 +37,12 @@
- Change files in extdata/Berman to be more like winTWQ format (no
value multipliers, except 10^5, 10^5, 10^5, and 10^8 on v1 to v4).
+- diagram(): rename 'what' argument to 'type'.
+
+- digram(): add new type of diagram, 'saturation', which is used to
+ plot saturation lines for minerals (where their affinity equals
+ zero).
+
CHANGES IN CHNOSZ 1.1.3 (2017-11-13)
------------------------------------
Modified: pkg/CHNOSZ/man/buffer.Rd
===================================================================
--- pkg/CHNOSZ/man/buffer.Rd 2018-03-07 15:46:30 UTC (rev 306)
+++ pkg/CHNOSZ/man/buffer.Rd 2018-03-14 17:54:46 UTC (rev 307)
@@ -71,7 +71,7 @@
mod.buffer("AC",logact=-10)
affinity(O2=c(-85,-70,4),T=c(25,100,4),return.buffer=TRUE)
# see below for a different strategy using the
-# 'what' argument of diagram
+# 'type' argument of diagram
## buffer made of three species
## Pyrite-Pyrrhotite-Magnetite (PPM)
@@ -133,7 +133,7 @@
affinity(return.buffer=TRUE)
# one can solve for the logarithm of activity of a
-# basis species using the 'what' argument of diagram
+# basis species using the 'type' argument of diagram
basis("CHNOS")
basis("CO2", 999)
species("acetic acid", -3)
@@ -142,16 +142,16 @@
lCO2 <- axis.label("CO2")
main <- substitute(a~~b~~c,list(a=lCO2, b="buffered by",
c="acetic acid"))
-d <- diagram(a, what="CO2", main=main)
+d <- diagram(a, type="CO2", main=main)
species(1, -10)
a <- affinity(O2=c(-85, -70, 4), T=c(25, 100, 4))
-d <- diagram(a, what="CO2", add=TRUE, lty=2)
+d <- diagram(a, type="CO2", add=TRUE, lty=2)
# add a legend
lAC <- expr.species("CH3COOH", log="aq")
ltext <- c(as.expression(lAC), -3, -10)
lty <- c(NA, 1, 2)
legend("topright", legend=ltext, lty=lty, bg="white")
-# do return.buffer and diagram(what) give the same results?
+# do return.buffer and diagram(type=...) give the same results?
and <- as.numeric(d$plotvals[[1]])
basis("CO2", "AC")
mod.buffer("AC", logact=-10)
Modified: pkg/CHNOSZ/man/diagram.Rd
===================================================================
--- pkg/CHNOSZ/man/diagram.Rd 2018-03-07 15:46:30 UTC (rev 306)
+++ pkg/CHNOSZ/man/diagram.Rd 2018-03-14 17:54:46 UTC (rev 307)
@@ -12,8 +12,8 @@
diagram(
# species affinities or activities
eout,
- # what to plot
- what = "loga.equil", alpha = FALSE, normalize = FALSE,
+ # type plot
+ type = "auto", alpha = FALSE, normalize = FALSE,
as.residue = FALSE, balance=NULL, groups=as.list(1:length(eout$values)),
# figure size and sides for axis tick marks
xrange=NULL, mar=NULL, yline=par("mgp")[1]+0.3, side=1:4,
@@ -22,7 +22,7 @@
# character sizes
cex=par("cex"), cex.names=1, cex.axis=par("cex"),
# line styles
- lty=NULL, lwd=par("lwd"), dotted=NULL, spline.method = NULL,
+ lty=NULL, lwd=par("lwd"), dotted=NULL, spline.method = NULL, contour.method = "edge",
# colors
col=par("col"), col.names=par("col"), fill=NULL,
fill.NA="slategray1", limit.water=TRUE,
@@ -39,7 +39,7 @@
\arguments{
\item{eout}{list, object returned by \code{\link{equilibrate}} or \code{\link{affinity}}}
- \item{what}{character, what property to calculate and plot}
+ \item{type}{character, type of plot, or name of basis species whose activity to plot}
\item{alpha}{logical or character (\samp{balance}), for speciation diagrams, plot degree of formation instead of activities?}
\item{normalize}{logical, divide chemical affinities by balance coefficients (rescale to whole formulas)?}
\item{as.residue}{logical, divide chemical affinities by balance coefficients (no rescaling)?}
@@ -61,6 +61,7 @@
\item{lwd}{numeric, line width}
\item{dotted}{numeric, how often to skip plotting points on predominance field boundaries (to gain the effect of dotted or dashed boundary lines)}
\item{spline.method}{character, method used in \code{\link{splinefun}}}
+ \item{contour.method}{character, labelling method used in \code{\link{contour}} (use NULL for no labels).}
\item{col}{character, color of activity lines (1D diagram) or predominance field boundaries (2D diagram), or colors of bars in a strip diagram (\code{strip})}
\item{col.names}{character, colors for labels of species}
\item{fill}{character, colors used to fill predominance fields}
@@ -140,14 +141,16 @@
\section{Affinity Diagrams}{
The function behaves differently when the output from \code{affinity} is being used instead of the equilibrium activities from \code{equilibrate}.
-If \code{what} is missing, and the number of dimensions is 0 or 1, a property of a reaction, such as the equilibrium constant (\samp{logK}), is plotted.
+If \code{type} is \samp{auto}, and the number of dimensions is 0 or 1, the property computed by \code{\link{affinity}} for each species is plotted.
+This is usually the affinity of the formation reaction, but can be set to some other property, such as the equilibrium constant (\samp{logK}).
+If \code{type} is \samp{auto}, and the number of dimensions is 2, then equilibrium predominance (maximum affinity) fields are plotted.
+This algorithm is based on a comparison of the affinities of the formation reactions scaled by the balancing coefficients that are determined by the \code{balance} argument.
-If \code{what} is missing, and the number of dimensions is 2, then highest potential (equilibrium predominance) fields are plotted.
-This algorithm is based on comparing the relative magnitudes of the affinities of the formation reactions, also referred to as the maximum affinity method.
-In this case, the balancing coefficients are determined usingi the \code{balance} argument.
-
-If \code{what} is the name of a basis species, it indicates to plot the equilibrium activity of a selected basis species in all of the formation reactions.
-A contour plot is made in the case of 2-D diagrams of the equilibrium activity of a basis species (see the \CO2-acetic acid example in \code{\link{buffer}}, and only the first species of interest is used in the calculation; a warning is produced if there is more than one.
+If \code{type} is \samp{saturation}, the function plots the line for each species where the affinity of formation equals zero.
+If for a given species no saturation line is possible or the range of the diagram is beyond the saturation line, the function prints a message instead.
+If \code{type} is the name of a basis species, the equilibrium activity of the selected basis species in each of the formation reactions is plotted.
+In the case of 2-D diagrams, both of these options use \code{\link{contour}} to draw the lines (see the \CO2-acetic acid example in \code{\link{buffer}}.
+The \samp{saturation} diagram can handle multiple species, but if \code{type} is the name a basis species, then only the first species of interest is used in the calculation, and a warning is produced if there is more than one.
}
\section{Activity Coefficients}{
@@ -187,7 +190,7 @@
basis("CHNOS")
species(c("ethanol", "lactic acid", "deoxyribose", "ribose"))
a <- affinity(T=c(0, 150))
-diagram(a, what="O2", legend.x="topleft", col=rev(rainbow(4)), lwd=2)
+diagram(a, type="O2", legend.x="topleft", col=rev(rainbow(4)), lwd=2)
title(main="Equilibrium logfO2 for 1e-3 mol/kg of CO2 and ... ")
### 1-D diagrams: logarithms of activities
Modified: pkg/CHNOSZ/tests/testthat/test-diagram.R
===================================================================
--- pkg/CHNOSZ/tests/testthat/test-diagram.R 2018-03-07 15:46:30 UTC (rev 306)
+++ pkg/CHNOSZ/tests/testthat/test-diagram.R 2018-03-14 17:54:46 UTC (rev 307)
@@ -7,7 +7,7 @@
a <- affinity()
expect_message(diagram(a, plot.it=FALSE), "balance: from moles of CO2 in formation reactions")
e <- equilibrate(a)
- expect_error(diagram(e, "Z"), "Z is not a basis species")
+ expect_error(diagram(e, "Z"), "Z is not a valid diagram type")
})
test_that("expected messages, errors and results arise using output from affinity()", {
@@ -81,9 +81,8 @@
basis("CHNOS")
species(c("alanine", "glycine", "serine", "methionine"))
a <- affinity(T=c(0, 200, 6), O2=c(-90, -60, 5))
- # TODO: fix plot.line() function in diagram() so that the plot can be made
- #expect_equal(diagram(a), diagram(a, plot.it=FALSE))
- expect_warning(diagram(a, what="CO2", plot.it=FALSE), "showing only first species in 2-D property diagram")
+ # now the warning is invokes next to the actual plotting, so no warning is produced with plot.it=FALSE 20180315
+ #expect_warning(diagram(a, type="CO2", plot.it=FALSE), "showing only first species in 2-D property diagram")
})
test_that("NaN values from equilibrate() are preserved (as NA in predominance calculation)", {
Modified: pkg/CHNOSZ/vignettes/anintro.Rmd
===================================================================
--- pkg/CHNOSZ/vignettes/anintro.Rmd 2018-03-07 15:46:30 UTC (rev 306)
+++ pkg/CHNOSZ/vignettes/anintro.Rmd 2018-03-14 17:54:46 UTC (rev 307)
@@ -807,7 +807,7 @@
* the buffers are active in calculations of affinity of other species;
* use <span style="color:red">`mod.buffer()`</span> to change or add buffers in `thermo$buffer`;
* [<span style="color:blue">`demo(buffer)`</span>](../demo) uses it for mineral buffers (solid lines).
-2. Use the `what` argument of <span style="color:green">`diagram()`</span> to solve for the activity of the indicated basis species:
+2. Use the `type` argument of <span style="color:green">`diagram()`</span> to solve for the activity of the indicated basis species:
* more convenient (the buffers come from the currently defined species of interest), but only a single basis species can be buffered, and it's not used in the calculation of affinity;
* [<span style="color:blue">`demo(buffer)`</span>](../demo) uses it for aqueous organic species as buffers (dotted and dashed lines).
More information about the CHNOSZ-commits
mailing list