[CHNOSZ-commits] r552 - in pkg/CHNOSZ: . R inst tests/testthat vignettes
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Jul 7 10:53:16 CEST 2020
Author: jedick
Date: 2020-07-07 10:53:15 +0200 (Tue, 07 Jul 2020)
New Revision: 552
Added:
pkg/CHNOSZ/tests/testthat/test-buffer.R
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/affinity.R
pkg/CHNOSZ/R/buffer.R
pkg/CHNOSZ/inst/NEWS
pkg/CHNOSZ/tests/testthat/test-berman.R
pkg/CHNOSZ/tests/testthat/test-mosaic.R
pkg/CHNOSZ/vignettes/OBIGT.Rmd
Log:
Make buffer calculations a little more robust and add tests
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2020-07-07 02:09:27 UTC (rev 551)
+++ pkg/CHNOSZ/DESCRIPTION 2020-07-07 08:53:15 UTC (rev 552)
@@ -1,6 +1,6 @@
Date: 2020-07-07
Package: CHNOSZ
-Version: 1.3.6-25
+Version: 1.3.6-26
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/affinity.R
===================================================================
--- pkg/CHNOSZ/R/affinity.R 2020-07-07 02:09:27 UTC (rev 551)
+++ pkg/CHNOSZ/R/affinity.R 2020-07-07 08:53:15 UTC (rev 552)
@@ -127,12 +127,6 @@
logK[[i]] <- logK[[i]] + thermo$species$logact[i]
for(j in 1:length(logact.basis.new)) {
logK[[i]] <- logK[[i]] - logact.basis.new[[j]] * thermo$species[i,j]
- # add ionization correction to proteins
- #if(i %in% is.buffer & length(grep('_',as.character(thermo$species$name[i])))>0 &
- # thermo$opt$ionize & rownames(mybasis)[j]=='H+') {
- # logK[[i]] <- logK[[i]] - logact.basis[[j]] *
- # as.data.frame(charge[[match(thermo$species$ispecies[i],names(charge))]])
- #}
}
}
lbn <- buffer(logK=logK,ibasis=ibasis,logact.basis=logact.basis.new,
@@ -167,9 +161,9 @@
if(nd < 3) {
for(i in 1:length(tb)) {
#tb[[i]] <- as.data.frame(tb[[i]])
- if(nd > 0) colnames(tb[[i]]) <-
+ if(nd > 0) rownames(tb[[i]]) <-
seq(args$lims[[1]][1],args$lims[[1]][2],length.out=args$lims[[1]][3])
- if(nd > 1) rownames(tb[[i]]) <-
+ if(nd > 1) colnames(tb[[i]]) <-
seq(args$lims[[2]][1],args$lims[[2]][2],length.out=args$lims[[2]][3])
}
}
Modified: pkg/CHNOSZ/R/buffer.R
===================================================================
--- pkg/CHNOSZ/R/buffer.R 2020-07-07 02:09:27 UTC (rev 551)
+++ pkg/CHNOSZ/R/buffer.R 2020-07-07 08:53:15 UTC (rev 552)
@@ -165,7 +165,7 @@
if(is.null(newbasis)) context <- '' else context <- paste(', ',basisnames[newbasis],' (conserved)',sep='')
reqtext <- paste(c2s(basisnames[ibasisrequested]),' (active)',sep='')
if(length(ibasisadded)==0) addtext <- '' else addtext <- paste(', ',c2s(basisnames[ibasisadded]),sep='')
- message(paste('buffer: ( ',bufname,' ) for activity of ',reqtext,addtext,context,sep=''))
+ message(paste("buffer: '", bufname, "' for activity of ", reqtext, addtext, context, sep = ""))
#print(bufbasis)
# there could still be stuff here (over-defined system?)
xx <- bufbasis[,-ibasis,drop=FALSE]
@@ -184,18 +184,14 @@
# do nothing
} else {
for(j in 1:ncol(xx)) {
- #if(i %in% are.proteins & colnames(xx)[j]=='H+' & thermo$opt$ionize) {
- # bs <- as.data.frame(charge[[match(ispecies[i],names(charge))]])[i,j] *
- # logact.basis[[match(colnames(xx)[j],rownames(thermo$basis))]]
- #} else bs <- xx[i,j] * logact.basis[[match(colnames(xx)[j],rownames(thermo$basis))]]
bs <- xx[i,j] * logact.basis[[match(colnames(xx)[j],basisnames)]]
if(!is.matrix(bs)) bs <- matrix(bs,byrow=TRUE,nrow=nrow(as.data.frame(logact.basis[[1]])))
- bs <- as.data.frame(bs)
+ if(ncol(bs)==1) b <- matrix(b)
b <- b - bs
}
}
}
- B[[i]] <- as.data.frame(b)
+ B[[i]] <- b
}
# a place to put the results
X <- rep(B[1],length(ibasis))
@@ -204,21 +200,7 @@
b <- numeric()
for(k in 1:length(B)) b <- c(b,B[[k]][i,j])
AAA <- A
- # here we calculate the coefficient on H+ if ionized proteins are present
- #if('H+' %in% colnames(A) & thermo$opt$ionize) {
- # H.coeff <- numeric()
- # for(l in 1:length(ispecies)) {
- # coeff <- as.data.frame(charge[[match(ispecies[l],names(charge))]])[i,j]
- # if(l %in% are.proteins) H.coeff <- c(H.coeff,coeff) else H.coeff <- c(H.coeff,0)
- # }
- # # apply the same type of balance and row-eliminating as above
- # if(length(ispecies)>1) {
- # H.coeff <- H.coeff/nb
- # for(l in 2:length(H.coeff)) H.coeff[l] <- H.coeff[l] - H.coeff[1]
- # H.coeff <- H.coeff[2:length(H.coeff)]
- # }
- # AAA[,match('H+',colnames(AAA))] <- H.coeff
- #}
+ # solve for the activities in the buffer
t <- solve(AAA,b)
for(k in 1:length(ibasis))
X[[k]][i,j] <- t[k]
@@ -226,7 +208,7 @@
}
# store results
for(i in 1:length(ibasis)) {
- if(ncol(X[[i]])==1) X[[i]] <- as.numeric(X[[i]][[1]])
+ if(ncol(X[[i]])==1) X[[i]] <- as.numeric(X[[i]])
else if(nrow(X[[i]])==1) X[[i]] <- as.matrix(X[[i]],nrow=1)
logact.basis[[ibasis[i]]] <- X[[i]]
}
Modified: pkg/CHNOSZ/inst/NEWS
===================================================================
--- pkg/CHNOSZ/inst/NEWS 2020-07-07 02:09:27 UTC (rev 551)
+++ pkg/CHNOSZ/inst/NEWS 2020-07-07 08:53:15 UTC (rev 552)
@@ -1,4 +1,4 @@
-CHANGES IN CHNOSZ 1.3.6-24 (2020-07-07)
+CHANGES IN CHNOSZ 1.3.6-26 (2020-07-07)
---------------------------------------
MAJOR CHANGES
Modified: pkg/CHNOSZ/tests/testthat/test-berman.R
===================================================================
--- pkg/CHNOSZ/tests/testthat/test-berman.R 2020-07-07 02:09:27 UTC (rev 551)
+++ pkg/CHNOSZ/tests/testthat/test-berman.R 2020-07-07 08:53:15 UTC (rev 552)
@@ -31,18 +31,18 @@
# now we can compare Berman and Helgeson G, H, S, Cp, V
# minerals with missing properties are not matched here
-# (i.e. fluorphlogopite, fluortremolite, glaucophane, and pyrope: no G and H in prop_Helgeson data)
+# (i.e. fluortremolite: no G and H in prop_Helgeson data)
test_that("Berman and Helgeson tabulated properties have large differences for few minerals", {
# which minerals differ in DGf by more than 4 kcal/mol?
idiffG <- which(abs(prop_Berman$G - prop_Helgeson$G) > 4000)
expect_match(mineral[idiffG],
- "paragonite|anthophyllite|antigorite|Ca-Al-pyroxene|lawsonite|margarite|merwinite")
+ "paragonite|anthophyllite|antigorite|Ca-Al-pyroxene|lawsonite|margarite|merwinite|fluorphlogopite")
# which minerals differ in DHf by more than 4 kcal/mol?
idiffH <- which(abs(prop_Berman$H - prop_Helgeson$H) > 4000)
expect_match(mineral[idiffH],
- "paragonite|anthophyllite|antigorite|Ca-Al-pyroxene|lawsonite|margarite|merwinite|clinozoisite")
+ "paragonite|anthophyllite|antigorite|Ca-Al-pyroxene|lawsonite|margarite|merwinite|clinozoisite|fluorphlogopite")
# which minerals differ in S by more than 4 cal/K/mol?
idiffS <- which(abs(prop_Berman$S - prop_Helgeson$S) > 4)
@@ -93,7 +93,7 @@
})
test_that("NA values of P are handled", {
- sresult <- subcrt("H2O", T = seq(0, 500, 100))
+ sresult <- suppressWarnings(subcrt("H2O", T = seq(0, 500, 100)))
T <- sresult$out$water$T
P <- sresult$out$water$P
# this stopped with a error prior to version 1.1.3-37
Added: pkg/CHNOSZ/tests/testthat/test-buffer.R
===================================================================
--- pkg/CHNOSZ/tests/testthat/test-buffer.R (rev 0)
+++ pkg/CHNOSZ/tests/testthat/test-buffer.R 2020-07-07 08:53:15 UTC (rev 552)
@@ -0,0 +1,38 @@
+context("buffer")
+
+# clear out any previous basis definition or database alterations
+suppressMessages(reset())
+
+test_that("simple buffer works in 0, 1, and 2 dimensions", {
+
+ add.OBIGT("SUPCRT92")
+ # define 4 temperatures
+ T <- c(200, 300, 400, 500)
+ # calculate SiO2 activity buffered by quartz at 1000 bar
+ logaSiO2 <- subcrt(c("quartz", "SiO2"), c(-1, 1), T = T, P = 1000)$out$logK
+
+ # set up system
+ basis(c("Al+3", "SiO2", "Na+", "K+", "H2O", "O2", "H+"))
+ species(c("K-feldspar", "albite", "paragonite", "dickite", "muscovite"))
+ # calculate logact(SiO2) with quartz buffer
+ basis("SiO2", "quartz")
+
+ # 0 dimensions: constant T
+ a0 <- affinity(T = 200, P = 1000)
+ expect_equal(a0$buffer[[1]], logaSiO2[1])
+
+ # 1 dimension: variable T
+ a1 <- affinity(T = T, P = 1000)
+ expect_equal(a1$buffer[[1]], logaSiO2)
+
+ # 2 dimensions: variable T and Na+ activity (affinity errored in version <= 1.3.6)
+ a2 <- affinity(T = c(200, 500, 4), "Na+" = c(1, 5, 5), P = 1000)
+ expect_equivalent(a2$buffer[[1]][, 1], logaSiO2)
+
+ # 2 dimensions: variable K+ and Na+ activity (diagram errored in version <= 1.3.6)
+ A2 <- affinity("K+" = c(1, 5, 5), "Na+" = c(1, 5, 5), T = 200, P = 1000)
+ #D2 <- diagram(A2)
+ expect_equivalent(unique(as.vector(A2$buffer[[1]])), logaSiO2[1])
+
+})
+
Modified: pkg/CHNOSZ/tests/testthat/test-mosaic.R
===================================================================
--- pkg/CHNOSZ/tests/testthat/test-mosaic.R 2020-07-07 02:09:27 UTC (rev 551)
+++ pkg/CHNOSZ/tests/testthat/test-mosaic.R 2020-07-07 08:53:15 UTC (rev 552)
@@ -16,11 +16,11 @@
m2_25 <- mosaic("NH3", "CO2", blend = FALSE)
expect_equal(a25$values, m2_25$A.species$values)
# make sure the function works when all affinities are NA
- a500 <- affinity(T=500)
+ a500 <- suppressWarnings(affinity(T=500))
# using blend=TRUE was failing prior to version 1.1.3-37
- m1_500 <- mosaic("NH3", "CO2", T=500)
+ m1_500 <- suppressWarnings(mosaic("NH3", "CO2", T=500))
expect_equal(a500$values, m1_500$A.species$values)
- m2_500 <- mosaic("NH3", "CO2", blend = FALSE, T=500)
+ m2_500 <- suppressWarnings(mosaic("NH3", "CO2", blend = FALSE, T=500))
expect_equal(a500$values, m2_500$A.species$values)
})
Modified: pkg/CHNOSZ/vignettes/OBIGT.Rmd
===================================================================
--- pkg/CHNOSZ/vignettes/OBIGT.Rmd 2020-07-07 02:09:27 UTC (rev 551)
+++ pkg/CHNOSZ/vignettes/OBIGT.Rmd 2020-07-07 08:53:15 UTC (rev 552)
@@ -389,7 +389,7 @@
<div id="nothing-open" style="display: block">
*Press a button above to show the citations in that data file.*
-*Or, <a href = "#showall" onclick='OpenAllDivs()'>press here</a> to show all citations.*
+*Or, <a href = "#showall" onclick='OpenAllDivs()'>click here</a> to show all citations.*
<b>Total count of species</b>: References were found for `r length(used)` of `r nrow(thermo()$OBIGT)` species in the default OBIGT database and `r length(optused)` optional species.
</div>
More information about the CHNOSZ-commits
mailing list