[CHNOSZ-commits] r758 - in pkg/CHNOSZ: . R inst inst/tinytest man vignettes
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Dec 1 07:51:13 CET 2022
Author: jedick
Date: 2022-12-01 07:51:13 +0100 (Thu, 01 Dec 2022)
New Revision: 758
Removed:
pkg/CHNOSZ/R/objective.R
pkg/CHNOSZ/R/revisit.R
pkg/CHNOSZ/inst/tinytest/test-objective.R
pkg/CHNOSZ/inst/tinytest/test-revisit.R
pkg/CHNOSZ/man/objective.Rd
pkg/CHNOSZ/man/revisit.Rd
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/NAMESPACE
pkg/CHNOSZ/R/equilibrate.R
pkg/CHNOSZ/R/examples.R
pkg/CHNOSZ/inst/NEWS.Rd
pkg/CHNOSZ/man/CHNOSZ-package.Rd
pkg/CHNOSZ/man/diagram.Rd
pkg/CHNOSZ/man/equilibrate.Rd
pkg/CHNOSZ/man/util.plot.Rd
pkg/CHNOSZ/vignettes/CHNOSZ.dia
pkg/CHNOSZ/vignettes/CHNOSZ.png
pkg/CHNOSZ/vignettes/anintro.Rmd
pkg/CHNOSZ/vignettes/mklinks.sh
Log:
Remove revisit()
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2022-12-01 05:43:36 UTC (rev 757)
+++ pkg/CHNOSZ/DESCRIPTION 2022-12-01 06:51:13 UTC (rev 758)
@@ -1,6 +1,6 @@
Date: 2022-12-01
Package: CHNOSZ
-Version: 1.9.9-49
+Version: 1.9.9-50
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/NAMESPACE
===================================================================
--- pkg/CHNOSZ/NAMESPACE 2022-12-01 05:43:36 UTC (rev 757)
+++ pkg/CHNOSZ/NAMESPACE 2022-12-01 06:51:13 UTC (rev 758)
@@ -10,7 +10,7 @@
"aminoacids", "ZC.col",
"pinfo", "protein.length", "protein.formula",
"read.fasta", "protein.basis", "add.protein",
- "unitize", "revisit", "seq2aa",
+ "unitize", "seq2aa",
"thermo.refs", "mod.OBIGT",
# examples
"examples", "demos", "mtitle",
@@ -26,7 +26,6 @@
"protein.OBIGT", "which.pmax",
"equil.boltzmann", "equil.reaction", "find.tp",
"ionize.aa", "MP90.cp", "aasum",
- "qqr", "RMSD", "CVRMSD", "spearman", "DGmix", "DDGmix", "DGtr",
"ratlab",
# demos
"protein.equil", "palply",
@@ -40,7 +39,6 @@
# (no other functions are used in the tests)
# other exported functions that are not used above
"checkEOS", "checkGHS", "check.OBIGT",
- "DGinf", "SD", "pearson", "shannon", "CV", "logact",
"basis.elements", "element.mu", "ibasis",
"water.SUPCRT92",
"nonideal", "uniprot.aa",
Modified: pkg/CHNOSZ/R/equilibrate.R
===================================================================
--- pkg/CHNOSZ/R/equilibrate.R 2022-12-01 05:43:36 UTC (rev 757)
+++ pkg/CHNOSZ/R/equilibrate.R 2022-12-01 06:51:13 UTC (rev 758)
@@ -274,7 +274,7 @@
Abar.range <- Abar.range - logadiff.mean / dlogadiff.dAbar
# one iteration is enough for the examples in the package
# but there might be a case where the range of logadiff doesn't cross zero
- # (e.g. for the carboxylic acid example in revisit.Rd)
+ # (e.g. for the carboxylic acid example previously in revisit.Rd)
logadiff.min <- logadiff(Abar.range[1], i)
logadiff.max <- logadiff(Abar.range[2], i)
iter <- 1
Modified: pkg/CHNOSZ/R/examples.R
===================================================================
--- pkg/CHNOSZ/R/examples.R 2022-12-01 05:43:36 UTC (rev 757)
+++ pkg/CHNOSZ/R/examples.R 2022-12-01 06:51:13 UTC (rev 758)
@@ -14,8 +14,7 @@
"solubility", "equilibrate",
"diagram", "mosaic", "mix",
"buffer", "nonideal", "NaCl",
- "add.protein", "ionize.aa",
- "objective", "revisit")
+ "add.protein", "ionize.aa")
plot.it <- FALSE
if(is.character(save.png))
png(paste(save.png,"%d.png",sep=""),width=500,height=500,pointsize=12)
Deleted: pkg/CHNOSZ/R/objective.R
===================================================================
--- pkg/CHNOSZ/R/objective.R 2022-12-01 05:43:36 UTC (rev 757)
+++ pkg/CHNOSZ/R/objective.R 2022-12-01 06:51:13 UTC (rev 758)
@@ -1,198 +0,0 @@
-# CHNOSZ/objective.R
-# objective functions for revisit(), findit()
-# all objective functions in one file, added attributes 20121009 jmd
-
-# input: a1 or loga1 is matrix with a column for each species and a row for each condition
-# output: vector with length equal to number of rows of a1 or loga1
-# attributes:
-# optimum - the target for optimization (minimal or maximal)
-## absolute - is the absolute value of the function used in the calculation? (TRUE or FALSE)
-
-SD <- structure(
- function(a1) {
- # calculate standard deviation
- SD <- apply(a1, 1, sd)
- return(SD)
- },
- optimum="minimal"
-)
-
-CV <- structure(
- function(a1) {
- # get the coefficient of variation
- SD <- SD(a1)
- CV <- SD / rowMeans(a1)
- return(CV)
- },
- optimum="minimal"
-)
-
-shannon <- structure(
- function(a1) {
- # shannon entropy
- p <- a1/rowSums(a1)
- H <- rowSums(-p*log(p))
- return(H)
- },
- optimum="maximal"
-)
-
-DGmix <- structure(
- function(loga1) {
- # Gibbs energy/2.303RT of ideal mixing
- a1 <- 10^loga1
- DGmix <- rowSums(a1*loga1)
- return(DGmix)
- },
- optimum="minimal"
-)
-
-qqr <- structure(
- function(loga1) {
- # correlation coefficient of a q-q plot (for log-normality comparisons) 20100901
- # first, a function to calculate qqr
- qqrfun <- function(x) {
- # this is to catch errors from qqnorm
- qqn <- tryCatch(qqnorm(x, plot.it=FALSE), error=identity)
- qqr <- if(inherits(qqn, "error")) NA else cor(qqn[[1]], qqn[[2]])
- }
- # apply the function to the rows of loga1
- qqr <- apply(loga1, 1, qqrfun)
- return(qqr)
- },
- optimum="maximal"
-)
-
-logact <- structure(
- function(loga1, loga2) {
- # the logarithm of activity of the species identified in loga2 (species number)
- return(loga1[, loga2[1]])
- },
- optimum="maximal"
-)
-
-spearman <- structure(
- function(loga1, loga2) {
- # Spearman's rho (rank correlation coefficient) 20100912
- spearfun <- function(a, b) {
- # based on help(dSpearman) in package SuppDists
- if(length(a)!=length(b)) stop("a and b must be same length")
- if(any(is.na(a)) | any(is.na(b))) return(NA)
- ra <- rank(a)
- rb <- rank(b)
- dr <- rb - ra
- d <- sum(dr^2)
- r <- length(a)
- x <- 1-6*d/(r*(r^2-1))
- return(x)
- }
- rho <- apply(loga1, 1, spearfun, b=loga2)
- return(rho)
- },
- optimum="maximal"
-)
-
-pearson <- structure(
- function(loga1, loga2) {
- # pearson correlation coefficient
- H <- apply(loga1, 1, cor, y=loga2)
- return(H)
- },
- optimum="maximal"
-)
-
-RMSD <- structure(
- function(loga1, loga2) {
- rmsd <- function(a, b) {
- # to calculate the root mean square deviation
- d <- b - a
- sd <- d^2
- msd <- mean(sd)
- rmsd <- sqrt(msd)
- return(rmsd)
- }
- RMSD <- apply(loga1, 1, rmsd, b=loga2)
- return(RMSD)
- },
- optimum="minimal"
-)
-
-CVRMSD <- structure(
- function(loga1, loga2) {
- # coefficient of variation of the root mean squared deviation
- RMSD <- RMSD(loga1, loga2)
- CVRMSD <- RMSD/rowMeans(loga1)
- return(CVRMSD)
- },
- optimum="minimal"
-)
-
-DDGmix <- structure(
- function(loga1, loga2) {
- # difference in Gibbs energy/2.303RT of ideal mixing
- a2 <- 10^loga2
- DGmix2 <- sum(a2*loga2)
- DGmix1 <- DGmix(loga1)
- DDGmix <- DGmix2 - DGmix1
- return(DDGmix)
- },
- optimum="minimal"
-)
-
-DGinf <- structure(
- function(a1, a2) {
- dginf <- function(a1, a2) {
- # informatic Gibbs energy/2.303RT difference between assemblages
- p1 <- a1/sum(a1)
- p2 <- a2/sum(a2)
- sum(p2 * log10(p2/p1))
- }
- DGinf <- apply(a1, 1, dginf, a2=a2)
- return(DGinf)
- },
- optimum="minimal"
-)
-
-DGtr <- structure(
- function(loga1, loga2, Astar) {
- dgtr <- function(loga1, loga2, Astar) {
- # calculate the Gibbs energy/2.303RT of transformation
- # from an initial to a final assemblage at constant T, P and
- # chemical potentials of basis species 20120917
- # loga1 - logarithms of activity of species in the initial assemblage
- # loga2 - logarithms of activity of species in the final assemblage
- # Astar - starved (of activity of the species of interest) values of chemical affinity/2.303RT
- # remove the logarithms
- a1 <- 10^loga1
- a2 <- 10^loga2
- # the molal Gibbs energy in the initial and final states
- # (derived from integrating A = Astar - RTln(a) from a1 to a2)
- G1 <- -a1*Astar + a1*loga1 - a1/log(10)
- G2 <- -a2*Astar + a2*loga2 - a2/log(10)
- # calculate the change in molal Gibbs energy for each species
- DG <- G2 - G1
- # return the sum
- return(sum(DG))
- }
- # we need to index both loga1 and Astar
- DGtr <- unlist(lapply(seq(nrow(loga1)), function(i) {
- dgtr(loga1[i, ], loga2, Astar[i, ])
- }))
- return(DGtr)
- },
- optimum="minimal"
-)
-
-### unexported functions ###
-
-# return the objective function named in objective,
-# or produce an error if the function has no optimum attribute
-get.objfun <- function(objective) {
- # get the function with the name given in 'objective'
- objfun <- get(objective)
- # perform a check on its usability
- if(!"optimum" %in% names(attributes(objfun)))
- stop(paste(objective, "is not a function with an attribute named 'optimum'"))
- # return the function
- return(objfun)
-}
Deleted: pkg/CHNOSZ/R/revisit.R
===================================================================
--- pkg/CHNOSZ/R/revisit.R 2022-12-01 05:43:36 UTC (rev 757)
+++ pkg/CHNOSZ/R/revisit.R 2022-12-01 06:51:13 UTC (rev 758)
@@ -1,266 +0,0 @@
-# CHNOSZ/revisit.R
-# 20090415 functions related to diversity calculations
-# 20100929 merged draw.diversity and revisit
-
-revisit <- function(eout, objective = "CV", loga2 = NULL, loga0 = NULL, ispecies = NULL,
- col = par("fg"), yline = 2, ylim = NULL, cex = par("cex"),
- lwd = par("lwd"), mar = NULL, side = 1:4, xlim = NULL, labcex = 0.6,
- pch = 1, main = NULL, plot.it = NULL, add = FALSE, plot.optval = TRUE,
- style.2D = "contour", bg = par("bg")) {
- # given a list of logarithms of activities of species
- # (as vectors or matrices or higher dimensional arrays)
- # calculate a diversity index or thermodynamic property
- # with the same dimensions 20090316 jmd v1
- # eout can be the output from equilibrate (enables plotting)
- # or simply a list of logarithms of activity
- # (each list entry must have the same dimensions)
-
- # if the entries have the same lengths they are assumed to be logarithms of activity
- ud <- unique(sapply(eout, length))
- if(length(ud)==1) {
- # eout is list of logarithms of activity
- if(missing(plot.it)) plot.it <- FALSE
- if(plot.it) stop("can't make a plot if 'eout' is not the output from equilibrate()")
- loga1 <- eout
- eout.is.eout <- FALSE
- } else {
- # test if eout is the output from equilibrate()
- if(!"loga.equil" %in% names(eout))
- stop(paste("the list provided in 'eout' is not the output from equilibrate()"))
- if(missing(plot.it)) plot.it <- TRUE
- loga1 <- eout$loga.equil
- eout.is.eout <- TRUE
- }
-
- # take a subset (or all) of the species
- if(is.null(ispecies)) ispecies <- 1:length(loga1)
- loga1 <- loga1[ispecies]
- # the dimensions
- dim1 <- dim(as.array(loga1[[1]]))
- # the number of dimensions
- nd <- ifelse(identical(dim1, 1L), 0, length(dim1))
- message(paste("revisit: calculating", objective, "in", nd, "dimensions"))
-
- # get the objective function
- objfun <- get.objfun(objective)
- # the arguments to the function (a1, [a2]) or (loga1, [loga2], [Astar])
- objargs <- names(formals(objfun))
-
- # these objectives only depend on the activities (a1/loga1):
- # shannon, SD, CV, QQR
- if(any(grepl("a1", objargs))) {
- # vectorize the list entries: a1/loga1
- loga1 <- sapply(loga1, c)
- # for 0-D case we need to keep loga1 as a 1-row matrix (sapply simplifies to vector)
- if(nd==0) loga1 <- t(loga1)
- # convert infinite values to NA
- loga1[is.infinite(loga1)] <- NA
- # if we have loga0, calculate the base-2 log ratio (loga1/loga0)
- base <- 10
- if(!is.null(loga0)) {
- loga1 <- t(t(loga1) - loga0) * log2(10)
- base <- 2
- }
- # remove logarithms if necessary
- if(any(grepl("loga1", objargs))) a1 <- loga1
- else a1 <- base^loga1
- }
-
- # these objectives also depend on reference activities (a2/loga2):
- # RMSD, CVRMSD, spearman, pearson, DGtr
- if(any(grepl("a2", objargs))) {
- # check that all needed arguments are present
- if(is.null(loga2)) stop(paste("loga2 must be supplied for", objective))
- # if loga2 is a single value, expand it into the dimensions of loga1
- if(length(loga2)==1) loga2 <- rep(loga2, length.out=ncol(loga1))
- # check that loga2 has the same length as loga1
- if(!identical(ncol(loga1), length(loga2))) stop(paste("loga2 has different length (",
- length(loga2), ") than list in eout (", ncol(loga1), ")", sep=""))
- # remove logarithms if necessary
- if(any(grepl("loga2", objargs))) a2 <- loga2
- else a2 <- base^loga2
- }
-
-
- # construct array of values: Astar (for DGtr)
- if(any(grepl("Astar", objargs))) {
- Astar <- eout$Astar[ispecies]
- # one row for each condition
- Astar <- sapply(Astar, as.vector)
- # for 0-D case we want a 1-row matrix (sapply simplifies to vector)
- if(nd==0) Astar <- t(Astar)
- }
-
- # calculation of the objective function
- # the symbol "H" is reminiscent of the first implemented target, shannon entropy
- if(length(objargs) == 1) H <- objfun(a1)
- else if(length(objargs) == 2) H <- objfun(a1, a2)
- else if(length(objargs) == 3) H <- objfun(a1, a2, Astar)
-
- # replace dims
- dim(H) <- dim1
-
- ## now on to plotting + assembling return values
- # for zero or more than two dimensions we'll just return the values
- # of the objective function and the index of the optimal value
- iopt <- optimal.index(H, objective)
- ret.val <- list(H=H, iopt=iopt)
- # get information about the x-axis
- if(eout.is.eout & nd > 0) {
- xname <- eout$vars[1]
- # the x-values
- xs <- eout$vals[[1]]
- xrange <- range(xs)
- } else xs <- NA
-
- # format the objective name if it's DGtr or DGinf
- if(objective=="DGtr") objexpr <- expr.property("DGtr/2.303RT")
- if(objective=="DGinf") objexpr <- expr.property("DGinf/2.303RT")
- else objexpr <- objective
-
- # make plots and assemble return values
- if(nd==0) {
- # a 0-D plot
- if(plot.it) {
- if(objective=="qqr") {
- # make a q-q plot for qqr
- qqnorm(loga1, col=col, pch=pch, main=NA)
- qqline(loga1)
- } else if(any(grepl("a2", objargs))) {
- # plot the points for a referenced objective
- xlab <- substitute(log~italic(a)[expt])
- ylab <- substitute(log~italic(a)[calc])
- if(is.null(xlim)) xlim <- extendrange(loga2)
- if(is.null(ylim)) ylim <- extendrange(loga1)
- # just initialize the plot here; add points below so that appear above lines
- plot(loga2, loga1, xlab=xlab, ylab=ylab, xlim=xlim, ylim=ylim, type="n")
- # add a 1:1 line
- lines(range(loga2), range(loga2), col="grey")
- # add a loess line
- if(!is.null(lwd)) {
- ls <- loess.smooth(loga2, loga1, family="gaussian")
- lines(ls$x, ls$y, col="red", lwd=lwd)
- }
- # add points
- points(loga2, loga1, pch=pch, col=col, cex=cex, bg=bg)
- } else plot.it <- FALSE
- # add a title
- if(missing(main)) main <- substitute(obj==H, list(obj=objexpr, H=round(H,3)))
- if(plot.it) title(main=main)
- }
- } else if(nd==1) {
- # locate the optimal value
- ixopt <- c(iopt)
- xopt <- xs[ixopt]
- optimum <- H[ixopt]
- ret.val <- list(H=H, ixopt=ixopt, xopt=xopt, optimum=optimum)
- # a 1-D plot
- if(plot.it) {
- if(is.null(ylim)) ylim <- extendrange(H, f=0.075)
- if(is.null(xlim)) xlim <- xrange
- if(!add) thermo.plot.new(xlim=xlim, ylim=ylim, xlab=axis.label(xname),
- ylab=objexpr, yline=yline, cex=cex, lwd=lwd, mar=mar, side=side)
- # plot the values
- lines(xs, c(H), col=col)
- # indicate the optimal value
- if(plot.optval) abline(v=xopt, lty=2)
- }
- } else if(nd==2) {
- # a 2-D plot
- # information about the y-axis
- if(eout.is.eout) {
- yname <- eout$vars[2]
- ys <- eout$vals[[2]]
- yrange <- range(ys)
- } else ys <- NA
- # locate the optimal value
- ixopt <- iopt[, 1]
- iyopt <- iopt[, 2]
- optimum <- H[ixopt, iyopt]
- ret.val <- list(H=H, ixopt=ixopt, iyopt=iyopt, xopt=xs[ixopt], yopt=ys[iyopt], optimum=optimum)
- if(plot.it) {
- # start the plot
- if(is.null(xlim)) xlim <- xrange
- if(is.null(ylim)) ylim <- yrange
- if(!add) thermo.plot.new(xlim=xlim, ylim=ylim, xlab=axis.label(xname),
- ylab=axis.label(yname), yline=yline, side=side, cex=cex, mar=mar)
- if(style.2D=="contour") contour(xs, ys, H, add=TRUE, labcex=labcex)
- else if(style.2D=="image") {
- # get colors based on number of species
- nspecies <- length(loga1)
- if(missing(col)) col <- heat.colors(nspecies)
- image(xs, ys, H, add=TRUE, col=col)
- } else stop(paste("2D plot style", style.2D, "not one of 'contour' or 'image'"))
- if(plot.optval) {
- # plot the location(s) of the extremum
- points(xs[ixopt], ys[iyopt], pch=8, cex=2)
- # show trajectories of the extrema
- iexts <- extremes(H, objective)
- # take out large jumps
- yext <- ys[iexts$y]
- yext.1 <- c(yext[2:length(yext)], yext[length(yext)])
- yext.2 <- c(yext[1], yext[1:length(yext)-1])
- yext[abs(yext.1-yext)/abs(diff(range(ys))) > 0.1] <- NA
- yext[abs(yext.2-yext)/abs(diff(range(ys))) > 0.1] <- NA
- lines(xs, yext, lty=3, col="blue")
- xext <- xs[iexts$x]
- xext.1 <- c(xext[2:length(xext)], xext[length(xext)])
- xext.2 <- c(xext[1], xext[1:length(xext)-1])
- xext[abs(xext.1-xext)/abs(diff(range(xs))) > 0.1] <- NA
- xext[abs(xext.2-xext)/abs(diff(range(xs))) > 0.1] <- NA
- lines(xext, ys, lty=3, col="seagreen")
- }
- }
- }
-
- # return the results
- return(invisible(ret.val))
-}
-
-### unexported functions ###
-
-optimal.index <- function(z, objective) {
- # for a vector, returns the index of the optimum value
- # for a matrix, returns the x, y coordinates of the optimum
- # find the minimum or maximum? look at attribute of objfun
- objfun <- get.objfun(objective)
- optimum <- attributes(objfun)$optimum
-# # do we care about the sign of the index?
-# if(tolower(target) %in% c("sd", "sd.log", "cv", "cv.log", "rmsd", "cvrmsd"))
-# doabs <- TRUE else doabs <- FALSE
-# if(doabs) z <- abs(z)
- # the value of the optimum
- if(optimum=="minimal") optval <- z[which.min(z)]
- else optval <- z[which.max(z)]
- # the index of the optimum
- # (or indices if multiple instances of the optimum)
- ret.val <- which(z==optval, arr.ind=TRUE)
- return(ret.val)
-}
-
-extremes <- function(z, objective) {
- # are we interested in a maximum or minimum?
- objfun <- get.objfun(objective)
- optimum <- attributes(objfun)$optimum
-# # do we care about the sign of the index?
-# if(tolower(target) %in% c("sd", "sd.log", "cv", "cv.log", "rmsd", "cvrmsd"))
-# doabs <- TRUE else doabs <- FALSE
-# if(doabs) z <- abs(z)
- # takes a matrix, returns the y as f(x) and x as f(y)
- # trajectories of the optimum
- y <- x <- numeric()
- xres <- ncol(z)
- yres <- nrow(z)
- if(optimum=="minimal") {
- for(i in 1:xres) y <- c(y, which.min(z[i,]))
- for(i in 1:yres) x <- c(x, which.min(z[,i]))
- } else {
- for(i in 1:xres) y <- c(y, which.max(z[i,]))
- for(i in 1:yres) x <- c(x, which.max(z[,i]))
- }
- # stop if we missed some
- if(length(x)!=xres) stop("optima not found for all y")
- if(length(y)!=yres) stop("optima not found for all x")
- return(list(x=x, y=y))
-}
-
Modified: pkg/CHNOSZ/inst/NEWS.Rd
===================================================================
--- pkg/CHNOSZ/inst/NEWS.Rd 2022-12-01 05:43:36 UTC (rev 757)
+++ pkg/CHNOSZ/inst/NEWS.Rd 2022-12-01 06:51:13 UTC (rev 758)
@@ -12,7 +12,7 @@
% links to vignettes 20220723
\newcommand{\viglink}{\ifelse{html}{\out{<a href="../CHNOSZ/doc/#1.html"><strong>#1.Rmd</strong></a>}}{\bold{#1.Rmd}}}
-\section{Changes in CHNOSZ version 1.9.9-49 (2022-12-01)}{
+\section{Changes in CHNOSZ version 1.9.9-50 (2022-12-01)}{
\subsection{MAJOR USER-VISIBLE CHANGES}{
\itemize{
@@ -141,8 +141,8 @@
\item Remove parallel calculations in \code{read.fasta()} and
\code{count.aa()} (they just made things slower in my tests).
- \item Remove \code{findit()} (functions for \dQuote{computations on
- chemical activities} in the extended workflow).
+ \item Remove \code{revisit()} and \code{findit()} (functions for
+ \dQuote{computations on chemical activities} in the extended workflow).
}
}
Deleted: pkg/CHNOSZ/inst/tinytest/test-objective.R
===================================================================
--- pkg/CHNOSZ/inst/tinytest/test-objective.R 2022-12-01 05:43:36 UTC (rev 757)
+++ pkg/CHNOSZ/inst/tinytest/test-objective.R 2022-12-01 06:51:13 UTC (rev 758)
@@ -1,16 +0,0 @@
-# Load default settings for CHNOSZ
-reset()
-
-# 20200728 moved from objective.Rd
-info <- "Calculations give expected results"
-loga1 <- t(-4:-1)
-loga2 <- loga1 + 1
-expect_true(qqr(loga1) < 1, info = info)
-expect_equal(RMSD(loga1, loga1), 0, info = info)
-expect_equal(RMSD(loga1, loga2), 1, info = info)
-expect_equal(CVRMSD(loga1, loga2), -0.4, info = info) # 1 / mean(-4:-1)
-expect_equal(spearman(loga1, loga2), 1, info = info)
-expect_equal(spearman(loga1, rev(loga2)), -1, info = info)
-# less statistical, more thermodynamical...
-expect_equal(DGmix(loga1), -0.1234, info = info) # as expected for decimal logarithms
-expect_equal(DDGmix(loga1, loga2), 0.0004, info = info)
Deleted: pkg/CHNOSZ/inst/tinytest/test-revisit.R
===================================================================
--- pkg/CHNOSZ/inst/tinytest/test-revisit.R 2022-12-01 05:43:36 UTC (rev 757)
+++ pkg/CHNOSZ/inst/tinytest/test-revisit.R 2022-12-01 06:51:13 UTC (rev 758)
@@ -1,92 +0,0 @@
-# Load default settings for CHNOSZ
-reset()
-
-# initial setup
-suppressMessages({
- # for numerical reproducibility, use the now-superseded properties of glycine 20190208
- mod.OBIGT("glycine", G = -90950, H = -124780, S = 39.29)
- basis("CHNOS")
- basis("O2", -65)
- species(c("leucine", "glycine", "glutamic acid"))
- # a 0-D case
- a <- affinity()
- e0 <- equilibrate(a)
- # a 1-D case
- a <- affinity(O2 = c(-75, -65))
- e1 <- equilibrate(a)
- # a 2-D case
- a <- affinity(O2 = c(-75, -65), NH3 = c(-4, 100, 3))
- e2 <- equilibrate(a)
- # a 3-D case
- a <- affinity(O2 = c(-75, -65, 4), H2O = c(-8, 0, 3), NH3 = c(-6, -4, 2))
- e3 <- equilibrate(a)
-})
-
-info <- "Inconsistent arguments produce an error"
-expect_error(CHNOSZ:::get.objfun("affinity"), "affinity is not a function with an attribute named 'optimum'", info = info)
-expect_error(revisit(list(1, 1), plot.it = TRUE), "can't make a plot if 'eout' is not the output from equilibrate\\(\\)", info = info)
-expect_error(revisit(list(1, c(1, 2))), "the list provided in 'eout' is not the output from equilibrate\\(\\)", info = info)
-expect_error(revisit(e1, "RMSD"), "loga2 must be supplied for RMSD", info = info)
-expect_error(revisit(e1, "RMSD", list(1, 1)), "loga2 has different length \\(2\\) than list in eout \\(3\\)", info = info)
-# commented because a plot is still initialized ...
-#expect_error(revisit(e2, "CV", style.2D = "xxx"), "2D plot style xxx not one of 'contour' or 'image'", info = info)
-
-info <- "0-D, 1-D, 2-D and 3-D calculations give identical results at the same conditions"
-r0.qqr <- revisit(e0, "qqr", plot.it = FALSE)
-r1.qqr <- revisit(e1, "qqr", plot.it = FALSE)
-r2.qqr <- revisit(e2, "qqr", plot.it = FALSE)
-r3.qqr <- revisit(e3, "qqr", plot.it = FALSE)
-# check that we get the same values
-expect_equal(c(r0.qqr$H), c(tail(r1.qqr$H, 1)), info = info)
-expect_equal(c(r0.qqr$H), r3.qqr$H[4, 3, 2], info = info)
-# check that we get the same index and same optimum
-expect_equal(r1.qqr$ixopt, r2.qqr$ixopt, check.attributes = FALSE, info = info)
-expect_equal(r1.qqr$optimum, r2.qqr$optimum, info = info)
-
-info <- "Non-referenced objectives give expected results"
-# the non-referenced objectives use only the logarithms of activities in eout (loga1)
-r1.cv <- revisit(e1, "CV", plot.it = FALSE)
-r1.sd <- revisit(e1, "SD", plot.it = FALSE)
-r1.shannon <- revisit(e1, "shannon", plot.it = FALSE)
-r1.qqr <- revisit(e1, "qqr", plot.it = FALSE)
-# the tests will alert us to significant numerical changes
-# but so far haven't been independently verified
-expect_equal(r1.cv$optimum, 0.29927, tolerance = 1e-5, info = info)
-expect_equal(r1.sd$optimum, 0.00027671, tolerance = 1e-5, info = info)
-expect_equal(r1.shannon$optimum, 1.067228, tolerance = 1e-5, info = info)
-expect_equal(r1.qqr$optimum, 0.9999584, tolerance = 1e-5, info = info)
-
-info <- "Referenced objectives give expected results"
-# the referenced objectives compare the logarithms of activities (loga1) to reference values (loga2)
-# the spearman correlation coefficient
-r1.spearman <- revisit(e1, "spearman", c(1, 2, 3), plot.it = FALSE)
-expect_equal(c(head(r1.spearman$H, 1)), -1, info = info) # perfect anti-rank correlation
-expect_equal(max(r1.spearman$H), 1, info = info) # perfect rank correlation
-# where logarithm of activity of the 3rd species (glutamic acid) maximizes
-r1.logact <- revisit(e1, "logact", 3, plot.it = FALSE)
-expect_equal(r1.logact$ixopt, 141, info = info)
-
-info <- "DGtr objective gives zero at equilibrium and >0 not at equilibrium"
-# let's use n-alkanes
-basis(c("CH4", "H2"), c("gas", "gas"))
-species(c("methane", "ethane", "propane", "butane"), "liq")
-# calculate equilibrium distribution over a range of logaH2
-a1 <- affinity(H2 = c(-10, -5, 101), exceed.Ttr = TRUE)
-e1 <- equilibrate(a1)
-# take the equilibrium distribution at logfH2 = -7.5 as the reference distribution
-loga2 <- list2array(e1$loga.equil)[51, ]
-# calculate the DGtr/RT relative to the reference distribution
-r1 <- revisit(e1, "DGtr", loga2 = loga2, plot.it = FALSE)
-# we should find a minimum of zero at logfH2 = -7.5
-expect_equal(min(r1$H), 0, info = info)
-expect_equal(r1$xopt, -7.5, info = info)
-
-# we can even go into 2 dimensions
-# (it's a slightly longer test, so don't run it on CRAN)
-if(!at_home()) exit_file("Skipping long test")
-a2 <- affinity(H2 = c(-10, -5, 101), T = c(0, 100, 101), exceed.Ttr = TRUE)
-e2 <- equilibrate(a2)
-r2 <- revisit(e2, "DGtr", loga2 = loga2, plot.it = FALSE)
-# we should DGtr = 0 at the temperature of the reference distribution (25 degC)
-expect_equal(min(r2$H), 0, info = info)
-expect_equal(r2$yopt, 25, info = info)
Modified: pkg/CHNOSZ/man/CHNOSZ-package.Rd
===================================================================
--- pkg/CHNOSZ/man/CHNOSZ-package.Rd 2022-12-01 05:43:36 UTC (rev 757)
+++ pkg/CHNOSZ/man/CHNOSZ-package.Rd 2022-12-01 06:51:13 UTC (rev 758)
@@ -20,7 +20,7 @@
Each help page (other than this one) has been given one of the following \dQuote{concept index entries}:
\itemize{
\item Main workflow: \code{\link{info}}, \code{\link{subcrt}}, \code{\link{basis}}, \code{\link{species}}, \code{\link{affinity}}, \code{\link{equilibrate}}, \code{\link{diagram}}
- \item Extended workflow: \code{\link{swap.basis}}, \code{\link{buffer}}, \code{\link{mosaic}}, \code{\link{objective}}, \code{\link{revisit}}
+ \item Extended workflow: \code{\link{swap.basis}}, \code{\link{buffer}}, \code{\link{mosaic}}
\item Thermodynamic data: \code{\link{data}}, \code{\link{extdata}}, \code{\link{add.OBIGT}}, \code{\link{util.data}}
\item Thermodynamic calculations: \code{\link{util.formula}}, \code{\link{makeup}}, \code{\link{util.units}}, \code{\link{Berman}}, \code{\link{nonideal}}, \code{\link{util.misc}}
\item Water properties: \code{\link{water}}, \code{\link{util.water}}, \code{\link{DEW}}, \code{\link{IAPWS95}}
Modified: pkg/CHNOSZ/man/diagram.Rd
===================================================================
--- pkg/CHNOSZ/man/diagram.Rd 2022-12-01 05:43:36 UTC (rev 757)
+++ pkg/CHNOSZ/man/diagram.Rd 2022-12-01 06:51:13 UTC (rev 758)
@@ -210,7 +210,7 @@
}
\seealso{
-\code{\link{Berman}}, \code{\link{mix}}, \code{\link{mosaic}}, \code{\link{nonideal}}, \code{\link{revisit}}, \code{\link{solubility}}, and \code{\link{util.plot}} are other help topics that use \code{diagram} in their examples.
+\code{\link{Berman}}, \code{\link{mix}}, \code{\link{mosaic}}, \code{\link{nonideal}}, \code{\link{solubility}}, and \code{\link{util.plot}} are other help topics that use \code{diagram} in their examples.
See the \code{\link{demos}} for even more examples.
}
Modified: pkg/CHNOSZ/man/equilibrate.Rd
===================================================================
--- pkg/CHNOSZ/man/equilibrate.Rd 2022-12-01 05:43:36 UTC (rev 757)
+++ pkg/CHNOSZ/man/equilibrate.Rd 2022-12-01 06:51:13 UTC (rev 758)
@@ -102,7 +102,7 @@
}
\seealso{
-\code{\link{diagram}} has examples of using \code{equilibrate} to make equilibrium activity diagrams. \code{\link{revisit}} can be used to perform further analysis of the equilibrium activities.
+\code{\link{diagram}} has examples of using \code{equilibrate} to make equilibrium activity diagrams.
\code{\link{palply}} is used by both \code{equil.reaction} and \code{equil.boltzmann} to parallelize intensive parts of the calculations.
See the vignette \viglink{multi-metal} for an example of balancing on two elements (N in the basis species, C in the formed species).
Deleted: pkg/CHNOSZ/man/objective.Rd
===================================================================
--- pkg/CHNOSZ/man/objective.Rd 2022-12-01 05:43:36 UTC (rev 757)
+++ pkg/CHNOSZ/man/objective.Rd 2022-12-01 06:51:13 UTC (rev 758)
@@ -1,124 +0,0 @@
-\encoding{UTF-8}
-\name{objective}
-\alias{objective}
-\alias{SD}
-\alias{CV}
-\alias{shannon}
-\alias{DGmix}
-\alias{qqr}
-\alias{logact}
-\alias{spearman}
-\alias{pearson}
-\alias{RMSD}
-\alias{CVRMSD}
-\alias{DDGmix}
-\alias{DGinf}
-\alias{DGtr}
-\title{Objective Functions}
-\description{
-Calculate statistical and thermodynamic quantities for activities of species.
-These functions can be specified as objectives in \code{\link{revisit}} in order to identify optimal chemical conditions.
-}
-
-\usage{
- SD(a1)
- CV(a1)
- shannon(a1)
- DGmix(loga1)
- qqr(loga1)
- logact(loga1, loga2)
- spearman(loga1, loga2)
- pearson(loga1, loga2)
- RMSD(loga1, loga2)
- CVRMSD(loga1, loga2)
- DDGmix(loga1, loga2)
- DGinf(a1, a2)
- DGtr(loga1, loga2, Astar)
-}
-
-\arguments{
- \item{a1}{numeric matrix, chemical activities of species}
- \item{loga1}{numeric matrix, logarithms of activity}
- \item{loga2}{numeric, reference values of logarithms of activity}
- \item{a2}{numeric, reference values of activity}
- \item{Astar}{numeric, reference values of chemical affinity}
-}
-
-\details{
-
-The value in \code{a1} or \code{loga1} is a matrix of chemical activities or logarithms of activity with a column for each species, and a row for each chemical condition.
-Except for calculations of the Shannon entropy, all logarithmic bases (including in the equations below) are decimal.
-
-\code{SD}, \code{CV} and \code{shannon} calculate the standard deviation, coefficient of variation, and Shannon entropy for the values in each row of \code{a1}. The Shannon entropy is calculated from the fractional abundances: H = sum(-p * log(p)) (natural logarithm), where p=a1/sum(a1).
-
-\code{DGmix} calculates the Gibbs energy/2.303RT of ideal mixing from pure components corresponding to one molal (unit activity) solutions: DGmix/2.303RT = sum(a1 * loga1) (cf. Eq. 7.20 of Anderson, 2005).
-
-\code{qqr} calculates the correlation coefficient on a quantile-quantile (Q-Q) plot (see \code{\link{qqnorm}}) for each row of \code{loga1}, giving some indication of the resemblance of the chemical activities to a log-normal distribution.
-
-\code{logact} returns the logarithm of activity of a single species identified by index in \code{loga2} (which of the species in the system).
-
-\code{spearman}, \code{pearson}, \code{RMSD} and \code{CVRMSD} calculate Spearman's rank correlation coefficient, the Pearson correlation coefficient, the root mean sqaured deviation (RMSD) and the coefficient of variation of the RMSD between each row of \code{loga1} and the values in \code{loga2}.
-The CVRMSD is computed as the RMSD divided by the mean of the values in \code{loga1}.
-
-\code{DDGmix} calculates the difference in Gibbs energy/2.303RT of ideal mixing between the assemblages with logarithms of activity \code{loga1} and \code{loga2}.
-
-\code{DGinf} calculates the difference in Gibbs energy/2.303RT attributed to relative informatic entropy between an initial assemblage with activities \code{a2} and final assemblage(s) with activities with activities in each row of \code{a1}.
-The equation used is DGinf/2.303RT = sum(p2 * (logp2 - logp1)), where p1 and p2 are the proportions, i.e. p1 = a1 / sum(a1) and p2 = a2 / sum(a2).
-This equation has the form of the Kullback-Leibler divergence, sometimes known as relative entropy (Ludovisi and Taticchi, 2006).
-In specific cases (systems where formulas of species are normalized by the balancing coefficients), the values of \code{DGinf} and \code{DGtr} are equal.
-
-\code{DGtr} calculates the change in Gibbs energy/2.303RT of a system in which species with initial logarithms of activitiy (\code{loga1}) are transformed to the same species with different final logarithms of activity (\code{loga2}) at constant temperature, pressure and chemical potentials of basis species.
-It is calculated as the sum over species of (G2-G1) where G1/RT = -a1*Astar + a1*loga1 - a1 + a constant (where a1 is 10^loga1), likewise for G2, and where \code{Astar} is the starved affinity, that is the affinity of the reaction to form one mole of the species at unit activity from the basis species in their defined activities.
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/chnosz -r 758
More information about the CHNOSZ-commits
mailing list