[CHNOSZ-commits] r422 - in pkg/CHNOSZ: . R inst man tests/testthat
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Mar 5 02:48:00 CET 2019
Author: jedick
Date: 2019-03-05 02:47:59 +0100 (Tue, 05 Mar 2019)
New Revision: 422
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/retrieve.R
pkg/CHNOSZ/inst/NEWS
pkg/CHNOSZ/man/retrieve.Rd
pkg/CHNOSZ/tests/testthat/test-retrieve.R
Log:
retrieve(): improve argument structure
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2019-03-04 08:13:33 UTC (rev 421)
+++ pkg/CHNOSZ/DESCRIPTION 2019-03-05 01:47:59 UTC (rev 422)
@@ -1,6 +1,6 @@
-Date: 2019-03-04
+Date: 2019-03-05
Package: CHNOSZ
-Version: 1.3.1-2
+Version: 1.3.1-3
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/retrieve.R
===================================================================
--- pkg/CHNOSZ/R/retrieve.R 2019-03-04 08:13:33 UTC (rev 421)
+++ pkg/CHNOSZ/R/retrieve.R 2019-03-05 01:47:59 UTC (rev 422)
@@ -1,10 +1,15 @@
# CHNOSZ/retrieve.R
# retrieve species with given elements
# 20190214 initial version
-# 20190224 use ... for multiple arguments (define a chemical system)
+# 20190224 define a chemical system using multiple arguments [defunct; use list() for chemical system]
# 20190304 update the stoichiometric matrix instead of doing a full recalculation when the database changes
+# 20190305 use c() for combination of elements, list() for chemical system,
+# and add 'ligands' argument to retrieve element-bearing species
-retrieve <- function(..., state = NULL, add.charge = TRUE, hide.groups = TRUE, req1 = FALSE) {
+retrieve <- function(elements = NULL, ligands = NULL, state = NULL, add.charge = TRUE, hide.groups = TRUE) {
+ ## empty argument handling
+ if(is.null(elements)) return(integer())
+
## stoichiometric matrix
# what are the formulas of species in the current database?
formula <- thermo()$obigt$formula
@@ -38,65 +43,73 @@
thermo("stoich" = stoich)
}
- ## species identification
- args <- list(...)
- ispecies <- numeric()
- # automatically add charge to a system 20190225
- if(add.charge & length(args) > 1) {
- if(!"Z" %in% unlist(args)) args <- c(args, "Z")
- }
- # for a numeric first argument, limit the result to only those species 20190225
- for(elements in args) {
- if(identical(elements, "all")) {
- ispecies <- 1:nrow(thermo()$obigt)
- names(ispecies) <- thermo()$obigt$formula
+ ## handle 'ligands' argument
+ if(!is.null(ligands)) {
+
+ # if 'ligands' is cr, liq, gas, or aq, use that as the state
+ if(any(ligands %in% c("cr", "liq", "gas", "aq")) & length(ligands)==1) {
+ state <- ligands
+ ispecies <- retrieve(elements, add.charge = add.charge)
} else {
- not.present <- ! elements %in% colnames(stoich)
- if(any(not.present)) {
- if(sum(not.present)==1) stop('"', elements[not.present], '" is not an element that is present in any species')
- else stop('"', paste(elements[not.present], collapse='", "'), '" are not elements that are present in any species')
+ # include the element in the system defined by the ligands list
+ ligands <- c(elements, as.list(ligands))
+ # call retrieve() for each argument and take the intersection
+ r1 <- retrieve(elements, add.charge = add.charge)
+ r2 <- retrieve(ligands, add.charge = add.charge)
+ ispecies <- intersect(r1, r2)
+ }
+
+ } else {
+
+ ## species identification
+ ispecies <- list()
+ # automatically add charge to a system 20190225
+ if(add.charge & is.list(elements) & !"Z" %in% elements) elements <- c(elements, "Z")
+ # proceed element-by-element
+ for(i in seq_along(elements)) {
+ element <- unlist(elements[i])
+ if(identical(element, "all")) {
+ ispecies[[i]] <- 1:nrow(thermo()$obigt)
+ } else {
+ if(! element %in% colnames(stoich)) stop('"', element, '" is not an element that is present in any species')
+ # identify the species that have the element
+ has.element <- rowSums(stoich[, element, drop = FALSE] != 0) == 1
+ ispecies[[i]] <- which(has.element)
}
- # identify the species that have the elements
- has.elements <- rowSums(stoich[, elements, drop = FALSE] != 0) == length(elements)
- # which species are these (i.e. the species index)
- # for req1, remember the species containing the first element 20190225
- if(length(ispecies)==0) ispecies1 <- which(has.elements)
- ispecies <- c(ispecies, which(has.elements))
- ispecies <- ispecies[!duplicated(ispecies)]
}
+
+ # now we have an list of ispecies (one vector for each element)
+ # what we do next depends on whether the argument is a list() or c()
+ if(is.list(elements)) {
+ # for a chemical system, all species are included that do not contain any other elements
+ ispecies <- unique(unlist(ispecies))
+ ielements <- colnames(thermo()$stoich) %in% elements
+ otherstoich <- thermo()$stoich[, !ielements]
+ iother <- rowSums(otherstoich[ispecies, ] != 0) > 0
+ ispecies <- ispecies[!iother]
+ } else {
+ # get species that have all the elements; the species must be present in each vector
+ # Reduce() hint from https://stackoverflow.com/questions/27520310/union-of-intersecting-vectors-in-a-list-in-r
+ ispecies <- Reduce(intersect, ispecies)
+ }
+
}
- # for a chemical system, defined by multiple arguments, the species can not contain any _other_ elements
- if(length(args) > 1) {
- syselements <- unlist(args)
- isyselements <- colnames(thermo()$stoich) %in% syselements
- notsysstoich <- thermo()$stoich[, !isyselements]
- iother <- rowSums(notsysstoich[ispecies, ] != 0) > 0
- ispecies <- ispecies[!iother]
- }
- # keep only species that contain the first element
- if(req1) {
- ispecies <- intersect(ispecies1, ispecies)
- names(ispecies) <- thermo()$obigt$name[ispecies]
- }
+
# exclude groups
if(hide.groups) {
igroup <- grepl("^\\[.*\\]$", thermo()$obigt$name[ispecies])
ispecies <- ispecies[!igroup]
}
- #if(hide.electron) {
- # ielectron <- names(ispecies) == "(Z-1)"
- # ispecies <- ispecies[!ielectron]
- #}
- #if(hide.proton) {
- # iproton <- names(ispecies) == "H+"
- # ispecies <- ispecies[!iproton]
- #}
- # filter states
+ # filter on state
if(!is.null(state)) {
istate <- thermo()$obigt$state[ispecies] %in% state
ispecies <- ispecies[istate]
}
- # for names, use e- instead of (Z-1)
+
+ # assign names; use e- instead of (Z-1)
+ names(ispecies) <- thermo()$obigt$formula[ispecies]
names(ispecies)[names(ispecies)=="(Z-1)"] <- "e-"
+ # if there's nothing, don't give it a name
+ if(length(ispecies)==0) ispecies <- integer()
ispecies
}
Modified: pkg/CHNOSZ/inst/NEWS
===================================================================
--- pkg/CHNOSZ/inst/NEWS 2019-03-04 08:13:33 UTC (rev 421)
+++ pkg/CHNOSZ/inst/NEWS 2019-03-05 01:47:59 UTC (rev 422)
@@ -1,4 +1,4 @@
-CHANGES IN CHNOSZ 1.3.1-2 (2019-03-04)
+CHANGES IN CHNOSZ 1.3.1-3 (2019-03-05)
--------------------------------------
- Add thermo/stoich.csv.xz (loaded as thermo()$stoich), containing a
@@ -8,6 +8,9 @@
- retrieve() now updates the stoichiometric matrix when the database
changes, instead of performing a full recalculation.
+- Add 'ligands' argument to retrieve(), for getting metal-bearing
+ species with a range of possible elements in the ligands.
+
CHANGES IN CHNOSZ 1.3.1 (2019-03-02)
------------------------------------
Modified: pkg/CHNOSZ/man/retrieve.Rd
===================================================================
--- pkg/CHNOSZ/man/retrieve.Rd 2019-03-04 08:13:33 UTC (rev 421)
+++ pkg/CHNOSZ/man/retrieve.Rd 2019-03-05 01:47:59 UTC (rev 422)
@@ -7,37 +7,35 @@
}
\usage{
- retrieve(..., state = NULL, add.charge = TRUE, hide.groups = TRUE,
- req1 = FALSE)
+ retrieve(elements = NULL, ligands = NULL, state = NULL,
+ add.charge = TRUE, hide.groups = TRUE)
}
\arguments{
- \item{...}{list, one or more arguments, each of which is a character vector with the names of one or more chemical elements}
+ \item{elements}{character, combination of elements, or list, elements in a chemical system}
+ \item{ligands}{character, elements present in any ligands}
\item{state}{character, filter the result on these state(s).}
\item{add.charge}{logical, add charge to the system?}
\item{hide.groups}{logical, exclude groups from the result?}
- \item{req1}{logical, include only species that contain the first element in the system?}
}
\details{
-This function retrieves the species in the thermodynamic database (see \code{\link{thermo}}) that have all of the elements specified in the arguments.
-A single argument is interpreted as a combination of one or more elements that must be present in each species.
+This function retrieves the species in the thermodynamic database (see \code{\link{thermo}}) that have the indicated \code{elements}.
+A character value of \code{elements} is interpreted as a combination of one or more elements that must be present in each species.
+A list value of \code{elements} is used to specify a chemical system -- the species must contain one or more of the indicated elements, but no other elements.
+\code{ligands}, if present, gives the elements that may be present in any ligands; this can be used to retrieve all species in a system bearing the \code{elements} (usually a single metal).
+
+The result includes charged species if \code{add.charge} is TRUE (the default) or the user supplies the \dQuote{element} of charge (\samp{Z}).
Results can be filtered on physical state by setting the \code{state} argument.
-
-If more than one argument is present, all of the species identified by each argument are combined, and all species containing any other elements are excluded.
-This can be used to retrieve all of the species in the database within a given chemical system.
-A chemical system includes charged species if \code{add.charge} is TRUE (the default) or the user supplies the \dQuote{element} of charge (\samp{Z}).
-If \code{req1} is TRUE, the result corresponds to the intersection of all of the species in the system with those identified by the first argument (i.e. those bearing the first element).
-
Groups used in group-additivity calculations, which have names with square brackets (e.g. [-CH2-]), are excluded unless \code{hide.groups} is FALSE.
-
-The return value is a named numeric vector giving the species index (i.e. rownumber(s) of \code{thermo()$obigt}) with names corresponding to the chemical formulas of the species.
-However, if the electron is in the result, its name (\samp{e-}) is used instead of its chemical formula (\samp{(Z-1)}).
-If the argument list is empty, then the function returns an empty (length 0) numeric value.
A special argument value \samp{all} can be used to retrieve all species in the thermodynamic database, including filtering on state and hiding of the groups.
+The return value is a named integer vector giving the species index (i.e. rownumber(s) of \code{thermo()$obigt}) with names corresponding to the chemical formulas of the species.
+If the electron is in the result, its name (\samp{e-}) is used instead of its chemical formula (\samp{(Z-1)}).
+An empty (length 0) integer value is returned if no \code{elements} are specified or no species are retrieved.
+
To speed up operation, the function uses a precalculated stoichiometric matrix for the default database, which is loaded with the package (see \code{\link{thermo}}).
-If the function detects a change to any chemical formulas in database, it updates the stoichiometric matrix using \code{\link{i2A}}.
+If the function detects a change to any chemical formulas in the database, it updates the stoichiometric matrix using \code{\link{i2A}}.
}
\seealso{
@@ -54,19 +52,22 @@
retrieve("Au")
# all species that have both Au and Cl
retrieve(c("Au", "Cl"))
-# all species that have Au and/or Cl, including charged species, but no other elements
+# Au-Cl system: species that have Au and/or Cl,
+# including charged species, but no other elements
+retrieve(list("Au", "Cl"))
+# all Au-bearing species in the Au-Cl system
retrieve("Au", "Cl")
-# all uncharged species that have Au and/or Cl
+# all uncharged Au-bearing species in the Au-Cl system
retrieve("Au", "Cl", add.charge = FALSE)
# minerals in the system SiO2-MgO-CaO-CO2
-retrieve("Si", "Mg", "Ca", "C", "O", state="cr")
+retrieve(list("Si", "Mg", "Ca", "C", "O"), state="cr")
-# make an Eh-pH diagram for aqueous S-bearing species
-basis("CHNOSe")
-ispecies <- retrieve("S", "O", "H", req1 = TRUE, state = "aq")
-species(ispecies)
-a <- affinity(pH = c(0, 14), Eh = c(-1, 1))
+# an Eh-pH diagram for Mn-bearing aqueous species
+basis(c("Mn+2", "H2O", "H+", "e-"))
+iMn <- retrieve("Mn", c("O", "H"), "aq")
+species(iMn)
+a <- affinity(pH = c(6, 14), Eh = c(-1, 1))
diagram(a, fill = "terrain")
}
Modified: pkg/CHNOSZ/tests/testthat/test-retrieve.R
===================================================================
--- pkg/CHNOSZ/tests/testthat/test-retrieve.R 2019-03-04 08:13:33 UTC (rev 421)
+++ pkg/CHNOSZ/tests/testthat/test-retrieve.R 2019-03-05 01:47:59 UTC (rev 422)
@@ -12,11 +12,7 @@
})
test_that("errors and recalculations produce expected messages", {
- # this should give an error about one non-element
expect_error(retrieve(c("A", "B", "C")), '"A" is not an element')
- # this should give an error about two non-elements
- expect_error(retrieve(c("A", "B", "C", "D")), '"A", "D" are not elements')
- # this should recalculate the stoichiometric matrix
add.obigt("SUPCRT92")
expect_message(retrieve("Ti"), "updating stoichiometric matrix")
reset()
More information about the CHNOSZ-commits
mailing list