[CHNOSZ-commits] r519 - in pkg/CHNOSZ: . R inst tests/testthat

noreply at r-forge.r-project.org noreply at r-forge.r-project.org
Tue Dec 10 03:01:32 CET 2019 

Author: jedick
Date: 2019-12-10 03:01:26 +0100 (Tue, 10 Dec 2019)
New Revision: 519

Modified:
   pkg/CHNOSZ/DESCRIPTION
   pkg/CHNOSZ/R/NaCl.R
   pkg/CHNOSZ/R/nonideal.R
   pkg/CHNOSZ/R/util.affinity.R
   pkg/CHNOSZ/R/util.list.R
   pkg/CHNOSZ/inst/NEWS
   pkg/CHNOSZ/tests/testthat/test-nonideal.R
Log:
affinity(): fix bug in IS-T calculations (workaround subcrt() limitation for grid = "IS")


Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION	2019-12-08 15:50:00 UTC (rev 518)
+++ pkg/CHNOSZ/DESCRIPTION	2019-12-10 02:01:26 UTC (rev 519)
@@ -1,6 +1,6 @@
-Date: 2019-12-05
+Date: 2019-12-10
 Package: CHNOSZ
-Version: 1.3.3-15
+Version: 1.3.3-16
 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/NaCl.R
===================================================================
--- pkg/CHNOSZ/R/NaCl.R	2019-12-08 15:50:00 UTC (rev 518)
+++ pkg/CHNOSZ/R/NaCl.R	2019-12-10 02:01:26 UTC (rev 519)
@@ -38,6 +38,8 @@
     gam_Na <- 10^gammas[[1]]$loggam
     gam_Cl <- 10^gammas[[2]]$loggam
     gam_NaCl <- 10^gammas[[3]]$loggam
+    # in case Setchenow calculations are turned off 20191209
+    if(length(gam_NaCl) == 0) gam_NaCl <- rep(1, length(T))
     # solve for m_Cl
     for(i in which(doit)) m_Cl[i] <- uniroot(A, c(0, m_tot), gam_NaCl=gam_NaCl[i], gam_Na=gam_Na[i], gam_Cl=gam_Cl[i], logK=logK[i])$root
     # calculate new total molality

Modified: pkg/CHNOSZ/R/nonideal.R
===================================================================
--- pkg/CHNOSZ/R/nonideal.R	2019-12-08 15:50:00 UTC (rev 518)
+++ pkg/CHNOSZ/R/nonideal.R	2019-12-10 02:01:26 UTC (rev 519)
@@ -142,10 +142,10 @@
       for(j in 1:ncol(myprops)) {
         pname <- colnames(myprops)[j]
         if(!pname %in% c("G", "H", "S", "Cp")) next
-        if(get("thermo", CHNOSZ)$opt$Setchenow == "bgamma") {
+        if(identical(get("thermo", CHNOSZ)$opt$Setchenow, "bgamma")) {
           myprops[, j] <- myprops[, j] + Setchenow(pname, IS, T, m_star, bgamma)
           didneutral <- TRUE
-        } else if(get("thermo", CHNOSZ)$opt$Setchenow == "bgamma0") {
+        } else if(identical(get("thermo", CHNOSZ)$opt$Setchenow, "bgamma0")) {
           myprops[, j] <- myprops[, j] + Setchenow(pname, IS, T, m_star, bgamma = 0)
           didneutral <- TRUE
         }

Modified: pkg/CHNOSZ/R/util.affinity.R
===================================================================
--- pkg/CHNOSZ/R/util.affinity.R	2019-12-08 15:50:00 UTC (rev 518)
+++ pkg/CHNOSZ/R/util.affinity.R	2019-12-10 02:01:26 UTC (rev 519)
@@ -62,7 +62,7 @@
   ibasisvar <- ibasisvar[!is.na(ibasisvar)]
   ## which vars are in P,T,IS
   varissubcrt <- vars %in% c("P","T","IS")
-  if(length(which(varissubcrt)) > 2) stop("sorry, currently only up to 2 of P,T,IS are supported")
+  if(length(which(varissubcrt)) > 2) stop("only up to 2 of P,T,IS are supported")
   ## categorize the basis species:
   # 0 - not in the vars; 1 - one of the vars
   ibasis <- 1:nbasis
@@ -70,9 +70,13 @@
   ibasis1 <- ibasis[ibasis %in% ibasisvar]
   if(identical(what,"logact.basis")) ispecies <- ibasis
   ## what subcrt variable is used to make a 2-D grid?
+  workaround.IS <- FALSE
   if(length(which(varissubcrt)) > 1 & !transect) {
-    if("IS" %in% vars) grid <- "IS"
-    else grid <- vars[varissubcrt][1]
+    grid <- vars[varissubcrt][1]
+    if("IS" %in% vars) {
+      if(grid != "IS") workaround.IS <- TRUE
+      grid <- "IS"
+    }
   } else grid <- NULL
   ### done argument processing
 
@@ -148,6 +152,16 @@
       if("IS" %in% vars) IS <- vals[[which(vars=="IS")]]
       s.args <- list(species=species,property=property,T=T,P=P,IS=IS,grid=grid,convert=FALSE,exceed.Ttr=exceed.Ttr,exceed.rhomin=exceed.rhomin)
       sout <- do.call("subcrt",s.args)
+      # 20191210 workaround a limitation in subcrt(): only IS (if present) can be the 'grid' variable
+      if(workaround.IS) {
+        lenIS <- length(IS)
+        # only T or P has length > 1
+        lenTP <- max(length(T), length(P))
+        # make an indexing matrix "byrow" then vectorize it (not byrow) to get permutation indices
+        indmat <- matrix(1:(lenIS * lenTP), nrow = lenIS, byrow = TRUE)
+        ind <- as.vector(indmat)
+        sout$out <- lapply(sout$out, function(x) x[indmat, ])
+      }
       # species indices are updated by subcrt() for minerals with phase transitions
       # e.g. i <- info("chalcocite"); subcrt(i, T=200)$species$ispecies == i + 1
       # so we should keep the original species index to be able to find the species in a provided 'sout'
@@ -177,7 +191,7 @@
         # two of T,P,IS
         ivar1 <- which(varissubcrt)[1]
         ivar2 <- which(varissubcrt)[2]
-        idim <- ivars(ivar2,ivars(ivar1))
+        idim <- ivars(ivar2, ivars(ivar1))
       }
     }
     return(aperm(array(x,mydim[idim]),idim))
@@ -187,7 +201,7 @@
   logK <- function() lapply(ispecies,X.reaction,"logK")
   # A/2.303RT
   A <- function() {
-    out <- lsub(X.fun("logK"),logQ())
+    out <- mapply(`-`, X.fun("logK"), logQ(), SIMPLIFY = FALSE)
     # deal with affinities of protein ionization here 20120527
     if("H+" %in% rownames(mybasis)) {
       # which species are proteins

Modified: pkg/CHNOSZ/R/util.list.R
===================================================================
--- pkg/CHNOSZ/R/util.list.R	2019-12-08 15:50:00 UTC (rev 518)
+++ pkg/CHNOSZ/R/util.list.R	2019-12-10 02:01:26 UTC (rev 519)
@@ -47,14 +47,6 @@
 
 ### unexported functions ###
 
-lsub <- function(x,y) {
-  # subtract elements of list y from those of
-  # list x to give list z
-  z <- x
-  for(i in 1:length(x)) z[[i]] <- x[[i]] - y[[i]]
-  return(z)
-}
-
 lsum <- function(x,y) {
   # sum up the respective elements of lists
   # x and y to give list z of the same length

Modified: pkg/CHNOSZ/inst/NEWS
===================================================================
--- pkg/CHNOSZ/inst/NEWS	2019-12-08 15:50:00 UTC (rev 518)
+++ pkg/CHNOSZ/inst/NEWS	2019-12-10 02:01:26 UTC (rev 519)
@@ -1,40 +1,25 @@
-CHANGES IN CHNOSZ 1.3.3-14 (2019-12-05)
+CHANGES IN CHNOSZ 1.3.3-16 (2019-12-10)
 ---------------------------------------
 
-- describe.reaction(): remove the Unicode double arrow. A simple equals
-  sign is displayed correctly on all systems.
+THERMODYNAMIC DATA
 
-- OBIGT: Change data for dawsonite, scheelite, and ferberite from
-  calories to Joules (as originally reported).
+- Change data for dawsonite, scheelite, and ferberite from calories to
+  Joules (as originally reported).
 
-- OBIGT: Use correct Cp coefficients for ferberite (2nd term in equation
-  of Wood and Samson, 2000 is off by a factor of 10). Thanks to
-  Xiangchong Liu for the bug report and David Polya for data, and both
-  for helpful advice.
+- Use correct Cp coefficients for ferberite (2nd term in equation of
+  Wood and Samson, 2000 is off by a factor of 10). Thanks to Xiangchong
+  Liu and David Polya for the bug report and assistance.
 
-- Fix bug in unexported function obigt2eos(): lambda (exponent in heat
-  capacity equation) was incorrectly going through a units conversion.
+- Move more carbonates from SUPCRT92 back into the default database:
+  artinite, azurite.
 
-- Remove files Sce.csv.xz and yeastgfp.csv.xz and functions yeastgfp()
-  and yeast.aa(). These have been moved to the JMDplots package
-  (https://github.com/jedick/JMDplots).
+FEATURE ENHANCEMENTS
 
-- Remove protein_refseq.csv.xz (based on RefSeq release 61). An updated
-  version, based on RefSeq release 95, is available in JMDplots.
-
-- Move some carbonates from SUPCRT92 back into default database:
-  artinite, azurite.
-
 - subcrt(): properties of minerals are now output at the listed upper
   T limit (or transition temperature, for e.g. cr1 -> cr2). Previously,
-  properties were set to NA at (and not only above) the T limit.
+  properties were set to NA at, and not only above, the T limit. Thanks
+  to Evgeniy Bastrakov for the request.
 
-- mosaic(): set the activities of only aqueous basis species to the
-  total activity (taken from the incoming basis() definition). This
-  fixes a bug where activities of minerals (particularly, native sulfur)
-  were unexpectedly changed; their logarithms are now set to 0. Thanks
-  to Evgeniy Bastrakov for the bug report.
-
 - Improve the handling of solids in equilibrate(). They are now excluded
   from the equilibrium calculation, but their stability fields are
   calculated using the maximum affinity method (with diagram()). Where
@@ -42,6 +27,29 @@
   to -999, and vice versa. Thanks to Feng Lai for the request and test
   case.
 
+CLEANUP
+
+- Remove files Sce.csv.xz and yeastgfp.csv.xz and functions yeastgfp()
+  and yeast.aa(). These have been moved to the JMDplots package
+  (https://github.com/jedick/JMDplots).
+
+- Remove protein_refseq.csv.xz (based on RefSeq release 61). An updated
+  version, based on RefSeq release 95, is available in JMDplots.
+
+BUG FIXES
+
+- describe.reaction(): remove the Unicode double arrow. Unlike that, a
+  simple equals sign is displayed correctly on all systems.
+
+- Fix bug in unexported function obigt2eos(): lambda (exponent in heat
+  capacity equation) was incorrectly going through a units conversion.
+
+- In mosaic(), set the activities of only aqueous basis species to the
+  total activity, which is taken from the incoming basis() definition.
+  This fixes a bug where activities of minerals (particularly, native
+  sulfur) were unexpectedly changed; their logarithms are now set to 0.
+  Thanks to Evgeniy Bastrakov for the bug report.
+
 - Fix a bug where berman() returned NA values for K-feldspar below
   298.15 K, due to improperly initialized values for the disorder
   properties. Thanks to Kaustubh Hakim for the bug report.
@@ -49,6 +57,10 @@
 - Use inherits() instead of class() for checking try() errors (don't
   use conditions of length > 1 in if() statements -- fix for R-devel).
 
+- Fix a bug in affinity() with calculations in IS-T (but not T-IS)
+  space. This involves a workaround for the limitation in subcrt() that
+  only IS (if present) can be the 'grid' variable.
+
 CHANGES IN CHNOSZ 1.3.3 (2019-08-02)
 ------------------------------------
 

Modified: pkg/CHNOSZ/tests/testthat/test-nonideal.R
===================================================================
--- pkg/CHNOSZ/tests/testthat/test-nonideal.R	2019-12-08 15:50:00 UTC (rev 518)
+++ pkg/CHNOSZ/tests/testthat/test-nonideal.R	2019-12-10 02:01:26 UTC (rev 519)
@@ -162,3 +162,18 @@
   # all done, reset nonideal and E.units options
   reset()
 })
+
+test_that("affinity has IS-T and T-IS calculations", {
+  basis("CHNOS+")
+  species(c("H2S", "HS-", "HSO4-", "SO4-2", "SO2"))
+  # reference values: affinity as a function of IS at T = 250 degC
+  a250 <- affinity(IS = seq(0, 1, 0.2), T = 250)
+  # increasing ionic strength should raise the affinity of HS-
+  expect_true(all(diff(a250$values[[2]]) > 0))
+  # T-IS grid; use different resolutions to help identify indexing bugs
+  a_T.IS <- affinity(IS = c(0, 1, 6), T = c(190, 310, 5))
+  expect_equivalent(a_T.IS$values[[2]][, 3], a250$values[[2]])
+  # IS-T grid
+  a_IS.T <- affinity(T = c(190, 310, 5), IS = c(0, 1, 6))
+  expect_equivalent(a_IS.T$values[[2]][3, ], a250$values[[2]])
+})

More information about the CHNOSZ-commits mailing list