[CHNOSZ-commits] r98 - in pkg/CHNOSZ: . R inst man tests/testthat vignettes
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Nov 9 05:45:01 CET 2015
Author: jedick
Date: 2015-11-09 05:45:00 +0100 (Mon, 09 Nov 2015)
New Revision: 98
Added:
pkg/CHNOSZ/inst/tests/
Removed:
pkg/CHNOSZ/inst/tests/
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/diagram.R
pkg/CHNOSZ/R/equilibrate.R
pkg/CHNOSZ/inst/NEWS
pkg/CHNOSZ/man/equilibrate.Rd
pkg/CHNOSZ/tests/testthat/test-diagram.R
pkg/CHNOSZ/tests/testthat/test-equilibrate.R
pkg/CHNOSZ/vignettes/equilibrium.Rnw
pkg/CHNOSZ/vignettes/equilibrium.lyx
Log:
add 'method' argument to equilibrate()
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2015-11-08 10:54:55 UTC (rev 97)
+++ pkg/CHNOSZ/DESCRIPTION 2015-11-09 04:45:00 UTC (rev 98)
@@ -1,6 +1,6 @@
-Date: 2015-11-08
+Date: 2015-11-09
Package: CHNOSZ
-Version: 1.0.6-3
+Version: 1.0.6-4
Title: Chemical Thermodynamics and Activity Diagrams
Author: Jeffrey Dick
Maintainer: Jeffrey Dick <j3ffdick at gmail.com>
Modified: pkg/CHNOSZ/R/diagram.R
===================================================================
--- pkg/CHNOSZ/R/diagram.R 2015-11-08 10:54:55 UTC (rev 97)
+++ pkg/CHNOSZ/R/diagram.R 2015-11-09 04:45:00 UTC (rev 98)
@@ -42,7 +42,7 @@
if(!"loga.equil" %in% names(eout)) {
eout.is.aout <- TRUE
# get the balancing coefficients
- n.balance <- balance(eout, balance)$n
+ n.balance <- balance(eout, balance)
}
} else if(what %in% rownames(eout$basis)) {
# to calculate the loga of basis species at equilibrium
Modified: pkg/CHNOSZ/R/equilibrate.R
===================================================================
--- pkg/CHNOSZ/R/equilibrate.R 2015-11-08 10:54:55 UTC (rev 97)
+++ pkg/CHNOSZ/R/equilibrate.R 2015-11-09 04:45:00 UTC (rev 98)
@@ -3,15 +3,15 @@
# of species in (metastable) equilibrium
equilibrate <- function(aout, balance=NULL, loga.balance=NULL,
- ispecies=1:length(aout$values), normalize=FALSE, as.residue=FALSE) {
+ ispecies=1:length(aout$values), normalize=FALSE, as.residue=FALSE,
+ method=c("boltzmann", "reaction")) {
### set up calculation of equilibrium activities of species from the affinities
### of their formation reactions from basis species at known activities
### split from diagram() 20120925 jmd
## number of possible species
nspecies <- length(aout$values)
## get the balancing coefficients
- balance <- balance(aout, balance)
- n.balance <- balance$n
+ n.balance <- balance(aout, balance)
## take selected species in 'ispecies'
if(length(ispecies)==0) stop("the length of ispecies is zero")
# take out species that have NA affinities
@@ -27,17 +27,15 @@
## number of species that are left
nspecies <- length(aout$values)
## say what the balancing coefficients are
- if(length(n.balance) < 100) msgout(paste("equilibrate: balancing coefficients are", c2s(n.balance), "\n"))
+ if(length(n.balance) < 100) msgout(paste("equilibrate: n.balance is", c2s(n.balance), "\n"))
## logarithm of total activity of the balance
- if(is.numeric(loga.balance))
- msgout(paste("equilibrate: logarithm of total", balance$description, "(from loga.balance) is", loga.balance, "\n"))
- else {
+ if(is.null(loga.balance)) {
# sum up the activities, then take absolute value
# in case n.balance is negative
sumact <- abs(sum(10^aout$species$logact * n.balance))
loga.balance <- log10(sumact)
- msgout(paste("equilibrate: logarithm of total", balance$description, "is", loga.balance, "\n"))
}
+ msgout(paste0("equilibrate: loga.balance is ", loga.balance, "\n"))
## normalize -- normalize the molar formula by the balance coefficients
m.balance <- n.balance
isprotein <- grepl("_", as.character(aout$species$name))
@@ -54,9 +52,14 @@
# and normalize the value by the nor-molar ratio
(aout$values[[i]] + aout$species$logact[i]) / m.balance[i]
})
- ## compute the equilibrium activities of species
- if(all(n.balance==1)) loga.equil <- equil.boltzmann(Astar, n.balance, loga.balance)
- else loga.equil <- equil.reaction(Astar, n.balance, loga.balance)
+ ## chose a method and compute the equilibrium activities of species
+ if(missing(method)) {
+ if(all(n.balance==1)) method <- method[1]
+ else method <- method[2]
+ }
+ msgout(paste("equilibrate: using", method[1], "method\n"))
+ if(method[1]=="boltzmann") loga.equil <- equil.boltzmann(Astar, n.balance, loga.balance)
+ else if(method[1]=="reaction") loga.equil <- equil.reaction(Astar, n.balance, loga.balance)
## if we normalized the formulas, get back to activities to species
if(normalize & !as.residue) {
loga.equil <- lapply(1:nspecies, function(i) {
@@ -83,7 +86,7 @@
# disadvantage:
# 1) only works for per-residue reactions
# 2) can create NaN logacts if the Astars are huge/small
-
+ if(any(n.balance!=1)) stop("won't run equil.boltzmann for balance <> 1")
# initialize output object
A <- Astar
# remember the dimensions of elements of Astar (could be vector,matrix)
@@ -150,6 +153,11 @@
Abarrange <- function(i) {
# starting guess of Abar (min/max) from range of Astar / n.balance
Abar.range <- range(Astar[i, ] / n.balance)
+ # diff(Abar.range) can't be 0 (dlogadiff.dAbar becomes NaN)
+ if(diff(Abar.range)==0) {
+ Abar.range[1] <- Abar.range[1] - 0.1
+ Abar.range[2] <- Abar.range[2] + 0.1
+ }
# the range of logadiff
logadiff.min <- logadiff(Abar.range[1], i)
logadiff.max <- logadiff(Abar.range[2], i)
@@ -161,7 +169,6 @@
# if one of them is infinite we might have a chance
if(is.infinite(logadiff.min)) {
# decrease the Abar range by increasing the minimum
- # the 0.99 is necessary at least for the residue=FALSE plot in protactiv.Rnw; 0.9 isn't enough!
Abar.range[1] <- Abar.range[1] + 0.99 * diff(Abar.range)
logadiff.min <- logadiff(Abar.range[1], i)
if(is.infinite(logadiff.min)) stop("FIXME: the second initial guess for Abar.min failed")
@@ -246,21 +253,16 @@
if(is.numeric(balance[1])) {
# a numeric vector
n.balance <- rep(balance, length.out=length(aout$values))
- if(all(n.balance==1)) msgout("balance: coefficients are unity\n")
- # use 'balance' below as a name
- if(all(n.balance==1)) balance.description <- "moles of species"
- else balance <- "moles of user-specified coefficients"
+ msgout("balance: from numeric argument value\n")
} else {
# "length" for balancing on protein length
if(identical(balance, "length")) {
if(!all(isprotein)) stop("length was the requested balance, but some species are not proteins")
n.balance <- protein.length(aout$species$name)
- balance.description <- "protein length"
- msgout(paste("balance: coefficients are", balance.description, "\n"))
+ msgout("balance: from protein length\n")
} else if(identical(balance, "volume")) {
n.balance <- info(aout$species$ispecies, check.it=FALSE)$V
- balance.description <- "volume"
- msgout(paste("balance: coefficients are", balance.description, "\n"))
+ msgout("balance: from volume")
} else {
# is the balance the name of a basis species?
if(length(ibalance)==0) {
@@ -269,13 +271,12 @@
}
# the name of the basis species (need this if we got ibalance which which.balance, above)
balance <- colnames(aout$species)[ibalance[1]]
- balance.description <- paste("moles of", balance)
- msgout(paste("balance: coefficients are", balance.description, "in formation reactions\n"))
+ msgout(paste("balance: from moles of", balance, "in formation reactions\n"))
# the balance vector
n.balance <- aout$species[, ibalance[1]]
# we check if that all formation reactions contain this basis species
if(any(n.balance==0)) stop("some species have no ", balance, " in the formation reaction")
}
}
- return(list(n=n.balance, description=balance.description))
+ return(n.balance)
}
Modified: pkg/CHNOSZ/inst/NEWS
===================================================================
--- pkg/CHNOSZ/inst/NEWS 2015-11-08 10:54:55 UTC (rev 97)
+++ pkg/CHNOSZ/inst/NEWS 2015-11-09 04:45:00 UTC (rev 98)
@@ -1,4 +1,4 @@
-CHANGES IN CHNOSZ 1.0.6-3 (2015-11-08)
+CHANGES IN CHNOSZ 1.0.6-4 (2015-11-09)
--------------------------------------
- Add functions usrfig() (get figure limits in user coordinates) and
@@ -31,6 +31,8 @@
nonideal() to set activity coefficients of proton and electron to
zero.
+- Add 'method' argument to equilibrate().
+
CHANGES IN CHNOSZ 1.0.6 (2015-10-19)
------------------------------------
Modified: pkg/CHNOSZ/man/equilibrate.Rd
===================================================================
--- pkg/CHNOSZ/man/equilibrate.Rd 2015-11-08 10:54:55 UTC (rev 97)
+++ pkg/CHNOSZ/man/equilibrate.Rd 2015-11-09 04:45:00 UTC (rev 98)
@@ -10,7 +10,8 @@
\usage{
equilibrate(aout, balance=NULL, loga.balance=NULL,
- ispecies=1:length(aout$values), normalize=FALSE, as.residue=FALSE)
+ ispecies=1:length(aout$values), normalize=FALSE, as.residue=FALSE,
+ method=c("boltzmann", "reaction"))
balance(aout, balance=NULL)
equil.boltzmann(Astar, n.balance, loga.balance)
equil.reaction(Astar, n.balance, loga.balance)
@@ -25,6 +26,7 @@
\item{Astar}{numeric, affinities of formation reactions excluding species contribution}
\item{n.balance}{numeric, number of moles of conserved component in the formation reactions of the species of interest}
\item{loga.balance}{numeric, logarithm of total activity of balanced quantity}
+ \item{method}{character, equilibration method to use}
}
\details{
@@ -58,7 +60,9 @@
\code{ispecies} can be supplied to identify a subset of the species to include in the calculation.
-\code{equil.boltzmann} is used to calculate the equilibrium activities if \code{balance} is \samp{1} (including the normalized result when \code{normalize} is TRUE), otherwise \code{equil.reaction} is called.
+\code{equil.boltzmann} is used to calculate the equilibrium activities if \code{balance} is \samp{1} (or when \code{normalize} or \code{as.residue} is TRUE), otherwise \code{equil.reaction} is called.
+The default behavior can be overriden by specifying either \samp{boltzmann} or \samp{reaction} in \code{method}.
+Using \code{equil.reaction} may be needed for systems with huge (negative or positive) affinities, where \code{equil.boltzmann} produces a NaN result.
}
@@ -76,6 +80,7 @@
In \code{equil.boltzmann} (algorithm available beginning with CHNOSZ-0.8), the chemical activities of species are calculated using the Boltzmann distribution.
This calculation is faster than the algorithm of \code{equil.reaction}, but is limited to systems where the transformations are all balanced on one mole of species.
+If \code{equil.boltzmann} is called with \code{balance} other than \samp{1}, it stops with an error.
}
Modified: pkg/CHNOSZ/tests/testthat/test-diagram.R
===================================================================
--- pkg/CHNOSZ/tests/testthat/test-diagram.R 2015-11-08 10:54:55 UTC (rev 97)
+++ pkg/CHNOSZ/tests/testthat/test-diagram.R 2015-11-09 04:45:00 UTC (rev 98)
@@ -5,7 +5,7 @@
basis("CHNOS")
species(c("glycine", "alanine"))
a <- affinity()
- expect_message(diagram(a, plot.it=FALSE), "coefficients are moles of CO2 in formation reactions")
+ 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")
})
@@ -22,7 +22,7 @@
# we can't calculate the equilibrium activity of a basis species if it's externally buffered
expect_error(diagram(a, "O2"), "is not numeric - was a buffer selected\\?")
# this one works - a barplot of A/2.303RT
- expect_message(diagram(a, plot.it=FALSE), "coefficients are moles of CO2 in formation reactions")
+ expect_message(diagram(a, plot.it=FALSE), "balance: from moles of CO2 in formation reactions")
# if we're plotting A/2.303RT the values can be divided by balancing coefficient or not
d.1 <- diagram(a, balance=1, plot.it=FALSE)
d.CO2 <- diagram(a, plot.it=FALSE)
Modified: pkg/CHNOSZ/tests/testthat/test-equilibrate.R
===================================================================
--- pkg/CHNOSZ/tests/testthat/test-equilibrate.R 2015-11-08 10:54:55 UTC (rev 97)
+++ pkg/CHNOSZ/tests/testthat/test-equilibrate.R 2015-11-09 04:45:00 UTC (rev 98)
@@ -19,18 +19,18 @@
test_that("equilibrate() gives expected messages and errors for balance calculation", {
# the following error is triggered by equil.react, not equil.boltzmann
expect_error(equilibrate(aone), "at least two species needed")
- expect_message(equilibrate(aacid), "coefficients are moles of CO2")
- expect_message(equilibrate(aacid), "balancing coefficients are 2 1 1 2 ")
- expect_message(equilibrate(aacid), "logarithm of total moles of CO2 is -2.221848")
- expect_message(equilibrate(aacid, loga.balance=-3), "logarithm of total moles of CO2 \\(from loga.balance\\) is -3")
+ expect_message(equilibrate(aacid), "balance: from moles of CO2")
+ expect_message(equilibrate(aacid), "n.balance is 2 1 1 2 ")
+ expect_message(equilibrate(aacid), "loga.balance is -2.221848")
+ expect_message(equilibrate(aacid, loga.balance=-3), "loga.balance is -3")
expect_error(equilibrate(aacid, balance="length"), "some species are not proteins")
- expect_message(equilibrate(aacidS), "coefficients are unity")
- expect_message(equilibrate(aacidS), "balancing coefficients are 1 1 1 1 1")
- expect_message(equilibrate(aacidS), "logarithm of total moles of species is -2.301029")
+ expect_message(equilibrate(aacidS), "balance: from numeric argument value")
+ expect_message(equilibrate(aacidS), "n.balance is 1 1 1 1 1")
+ expect_message(equilibrate(aacidS), "loga.balance is -2.301029")
expect_error(equilibrate(aacidS, balance="CO2"), "some species have no CO2 in the formation reaction")
- expect_message(equilibrate(aprot), "coefficients are protein length")
- expect_message(equilibrate(aprot), "balancing coefficients are 129 153 124 104")
- expect_message(equilibrate(aprot), "logarithm of total protein length is -0.292429")
+ expect_message(equilibrate(aprot), "balance: from protein length")
+ expect_message(equilibrate(aprot), "n.balance is 129 153 124 104")
+ expect_message(equilibrate(aprot), "loga.balance is -0.292429")
expect_message(equilibrate(aprot, normalize=TRUE), "using 'normalize' for molar formulas")
})
@@ -78,8 +78,51 @@
expect_true(all(abs(loga.equil-loga.ref) < 0.36))
})
-#}
+test_that("equilibrate() can be used for huge values of Astar", {
+ ## working out some bugs and testing new 'method' argument 20151109
+ ## first, demonstrate that equil.reaction works where equil.boltzmann doesn't
+ # minimal example: Astar=c(0, 0), n.balance=c(1, 1), loga.balance=0
+ # results in equal activities of two species
+ eb0 <- equil.boltzmann(c(0, 0), c(1, 1), 0)
+ expect_equal(unlist(eb0), rep(log10(0.5), 2))
+ # Astar=c(-330, -330)
+ # result is NaN (we probably get an Inf-Inf somewhere)
+ eb330 <- equil.boltzmann(c(-330, -330), c(1, 1), 0)
+ expect_equal(unlist(eb330), rep(NaN, 2))
+ # (fixed bug: while loop in equil.reaction tested a NaN value)
+ # (dlogadiff.dAbar <- 0 / 0)
+ er330 <- equil.reaction(c(-330, -330), c(1, 1), 0)
+ expect_equal(er330, eb0)
+
+ ## second, set up extreme test case and show boltzmann method produces NaN (is.na)
+ basis("CHNOS")
+ basis("O2", 200)
+ species(c("glycine", "alanine", "proline"))
+ a <- affinity()
+ expect_message(eb <- equilibrate(a, balance=1), "using boltzmann method")
+ expect_true(all(is.na(unlist(eb$loga.equil))))
+
+ ## third, check we can use method="reaction"
+ expect_message(er1 <- equilibrate(a, balance=1, method="reaction"), "using reaction method")
+ expect_false(any(is.na(unlist(er1$loga.equil))))
+ # is it an equilibrium solution?
+ species(1:3, unlist(er1$loga.equil))
+ a1 <- affinity()
+ expect_equal(diff(range(unlist(a1$values))), 0)
+
+ ## third, check that we can use arbitrary numeric balance specification
+ # (balance <> 1 here means equilibrate will call equil.reaction)
+ expect_message(er11 <- equilibrate(a, balance=1.000001), "using reaction method")
+ species(1:3, unlist(er11$loga.equil))
+ a11 <- affinity()
+ expect_equal(unlist(a1$values), unlist(a11$values))
+
+ ## fourth, check that equil.boltzmann won't run for balance <> 1
+ expect_error(equilibrate(a, balance=1.000001, method="boltzmann"), "won't run equil.boltzmann")
+})
+
+
# references
# Seewald, J. S. (2001)
Modified: pkg/CHNOSZ/vignettes/equilibrium.Rnw
===================================================================
--- pkg/CHNOSZ/vignettes/equilibrium.Rnw 2015-11-08 10:54:55 UTC (rev 97)
+++ pkg/CHNOSZ/vignettes/equilibrium.Rnw 2015-11-09 04:45:00 UTC (rev 98)
@@ -610,7 +610,7 @@
if(loga==-1) d <- diagram(a)
else d <- diagram(a, add=TRUE, names=NULL)
iCu <- which(d$predominant == 1, arr.ind=TRUE)
- text(a$vals[[1]][max(iCu[, 1])] - 0.2, a$vals[[2]][min(iCu[, 2])] + 0.2, loga)
+ text(a$vals[[1]][max(iCu[, 1])] - 0.03, a$vals[[2]][min(iCu[, 2])] + 0.2, adj=1, loga)
}
water.lines()
@@ -679,7 +679,7 @@
equilibrate the basis species rather than compose the diagram using
the predominant basis species).
\begin{itemize}
-\item See the examples in \texttt{?mosaic} and \texttt{demo(``mosaic'')}
+\item See the examples in \texttt{?mosaic} and \texttt{demo(\textquotedbl{}mosaic\textquotedbl{})}
for calculations of mineral stabilities in the Fe-S-O-$\mathrm{H_{2}O}$
system (after \citealp{GC65}).
\item A calculation of copper solubility limits and speciation with aqueous
Modified: pkg/CHNOSZ/vignettes/equilibrium.lyx
===================================================================
--- pkg/CHNOSZ/vignettes/equilibrium.lyx 2015-11-08 10:54:55 UTC (rev 97)
+++ pkg/CHNOSZ/vignettes/equilibrium.lyx 2015-11-09 04:45:00 UTC (rev 98)
@@ -2638,8 +2638,8 @@
\begin_layout Plain Layout
- text(a$vals[[1]][max(iCu[, 1])] - 0.2, a$vals[[2]][min(iCu[, 2])] + 0.2,
- loga)
+ text(a$vals[[1]][max(iCu[, 1])] - 0.03, a$vals[[2]][min(iCu[, 2])] + 0.2,
+ adj=1, loga)
\end_layout
\begin_layout Plain Layout
@@ -2919,15 +2919,7 @@
\family default
and
\family typewriter
-demo(
-\begin_inset Quotes eld
-\end_inset
-
-mosaic
-\begin_inset Quotes erd
-\end_inset
-
-)
+demo("mosaic")
\family default
for calculations of mineral stabilities in the Fe-S-O-
\begin_inset Formula $\mathrm{H_{2}O}$
More information about the CHNOSZ-commits
mailing list