[CHNOSZ-commits] r75 - in pkg/CHNOSZ: . R inst inst/tests man vignettes
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun Feb 15 07:42:56 CET 2015
Author: jedick
Date: 2015-02-15 07:42:50 +0100 (Sun, 15 Feb 2015)
New Revision: 75
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/diagram.R
pkg/CHNOSZ/R/equilibrate.R
pkg/CHNOSZ/R/wjd.R
pkg/CHNOSZ/inst/NEWS
pkg/CHNOSZ/inst/tests/test-diagram.R
pkg/CHNOSZ/inst/tests/test-equilibrate.R
pkg/CHNOSZ/man/diagram.Rd
pkg/CHNOSZ/man/equilibrate.Rd
pkg/CHNOSZ/vignettes/hotspring.Rnw
pkg/CHNOSZ/vignettes/hotspring.lyx
pkg/CHNOSZ/vignettes/wjd.Rnw
pkg/CHNOSZ/vignettes/wjd.lyx
Log:
add 'as.residue' argument to diagram()
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2015-01-07 13:06:31 UTC (rev 74)
+++ pkg/CHNOSZ/DESCRIPTION 2015-02-15 06:42:50 UTC (rev 75)
@@ -1,6 +1,6 @@
-Date: 2015-01-07
+Date: 2015-02-15
Package: CHNOSZ
-Version: 1.0.3-12
+Version: 1.0.3-13
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-01-07 13:06:31 UTC (rev 74)
+++ pkg/CHNOSZ/R/diagram.R 2015-02-15 06:42:50 UTC (rev 75)
@@ -9,7 +9,7 @@
# primary input
eout,
# what to plot
- what="loga.equil", alpha=FALSE, normalize=FALSE, balance=NULL,
+ what="loga.equil", alpha=FALSE, normalize=FALSE, as.residue=FALSE, balance=NULL,
groups=as.list(1:length(eout$values)), xrange=NULL,
# plot dimensions
mar=NULL, yline=par("mgp")[1]+0.7, side=1:4,
@@ -24,7 +24,7 @@
# labels
names=NULL, main=NULL, legend.x="topright",
# plotting controls
- add=FALSE, plot.it=TRUE
+ add=FALSE, plot.it=TRUE, tplot=TRUE
) {
### argument handling ###
@@ -52,12 +52,6 @@
} else if(what=="loga.equil" & !"loga.equil" %in% names(eout)) stop("'eout' is not the output from equil()")
else if(what!="loga.equil") stop(what, " is not a basis species or 'loga.equil'")
- ## we only normalize if eout.is.aout
- if(normalize) {
- if(!eout.is.aout) stop("normalizing formulas is only possible if 'eout' is the output from affinity()")
- else msgout("diagram: normalizing formulas in calculation of predominant species\n")
- }
-
## consider a different number of species if we're grouping them together
ngroups <- length(groups)
@@ -81,6 +75,20 @@
msgout(paste("diagram: plotting", plotvar, "from affinity(), divided by balancing coefficients\n"))
}
+ ## number of dimensions (T, P or chemical potentials that are varied)
+ # length(eout$vars) - the number of variables = the maximum number of dimensions
+ # length(dim(eout$values[[1]])) - nd=1 if it was a transect along multiple variables
+ nd <- min(length(eout$vars), length(dim(eout$values[[1]])))
+
+ ## when can normalize and as.residue be used
+ if(normalize | as.residue) {
+ if(normalize & as.residue) stop("'normalize' and 'as.residue' can not both be TRUE")
+ if(!eout.is.aout) stop("'normalize' or 'as.residue' can be TRUE only if 'eout' is the output from affinity()")
+ if(nd!=2) stop("'normalize' or 'as.residue' can be TRUE only for a 2-D (predominance) diagram")
+ if(normalize) msgout("diagram: using 'normalize' in calculation of predominant species\n")
+ else msgout("diagram: using 'as.residue' in calculation of predominant species\n")
+ }
+
## sum activities of species together in groups 20090524
# using lapply/Reduce 20120927
if(!missing(groups)) {
@@ -139,6 +147,7 @@
# TODO: see vignette for an explanation for how this is normalizing
# the formulas in a predominance calculation
if(normalize & eout.is.aout) pv[[i]] <- (pv[[i]] + eout$species$logact[i] / n.balance[i]) - log10(n.balance[i])
+ else if(as.residue & eout.is.aout) pv[[i]] <- pv[[i]] + eout$species$logact[i] / n.balance[i]
}
predominant <- which.pmax(pv)
dim(predominant) <- dim(pv[[1]])
@@ -153,11 +162,6 @@
### general plot parameters ###
- ## number of dimensions (T, P or chemical potentials that are varied)
- # length(eout$vars) - the number of variables = the maximum number of dimensions
- # length(dim(eout$values[[1]])) - nd=1 if it was a transect along multiple variables
- nd <- min(length(eout$vars), length(dim(eout$values[[1]])))
-
## handle line type/width/color arguments
if(is.null(lty)) lty <- 1:ngroups
lty <- rep(lty, length.out=ngroups)
@@ -216,7 +220,8 @@
myval <- sapply(plotvals, xfun)
ylim <- extendrange(myval)
}
- thermo.plot.new(xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab, cex=cex, mar=mar, yline=yline, side=side)
+ if(tplot) thermo.plot.new(xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab, cex=cex, mar=mar, yline=yline, side=side)
+ else plot(0, 0, type="n", xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)
}
# draw the lines
for(i in 1:length(plotvals)) lines(xvalues, plotvals[[i]], col=col[i], lty=lty[i], lwd=lwd[i])
@@ -397,8 +402,9 @@
if(!add) {
if(is.null(xlab)) xlab <- axis.label(eout$vars[1], basis=eout$basis)
if(is.null(ylab)) ylab <- axis.label(eout$vars[2], basis=eout$basis)
- thermo.plot.new(xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab,
+ if(tplot) thermo.plot.new(xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab,
cex=cex, cex.axis=cex.axis, mar=mar, yline=yline, side=side)
+ else plot(0, 0, type="n", xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)
# add a title
if(!is.null(main)) title(main=main)
}
Modified: pkg/CHNOSZ/R/equilibrate.R
===================================================================
--- pkg/CHNOSZ/R/equilibrate.R 2015-01-07 13:06:31 UTC (rev 74)
+++ pkg/CHNOSZ/R/equilibrate.R 2015-02-15 06:42:50 UTC (rev 75)
@@ -3,7 +3,7 @@
# of species in (metastable) equilibrium
equilibrate <- function(aout, balance=NULL, loga.balance=NULL,
- ispecies=1:length(aout$values), normalize=FALSE, stay.normal=FALSE) {
+ ispecies=1:length(aout$values), normalize=FALSE, as.residue=FALSE) {
### 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
@@ -41,10 +41,11 @@
## normalize -- normalize the molar formula by the balance coefficients
m.balance <- n.balance
isprotein <- grepl("_", as.character(aout$species$name))
- if(normalize) {
+ if(normalize | as.residue) {
if(any(n.balance < 0)) stop("one or more negative balancing coefficients prohibit using normalized molar formulas")
n.balance <- rep(1, nspecies)
- msgout(paste("equilibrate: normalizing molar formulas by the balancing coefficients\n"))
+ if(as.residue) msgout(paste("equilibrate: using 'as.residue' for molar formulas\n"))
+ else msgout(paste("equilibrate: using 'normalize' for molar formulas\n"))
} else m.balance <- rep(1, nspecies)
## Astar: the affinities/2.303RT of formation reactions with
## formed species in their standard-state activities
@@ -57,7 +58,7 @@
if(all(n.balance==1)) loga.equil <- equil.boltzmann(Astar, n.balance, loga.balance)
else loga.equil <- equil.reaction(Astar, n.balance, loga.balance)
## if we normalized the formulas, get back to activities to species
- if(normalize & !stay.normal) {
+ if(normalize & !as.residue) {
loga.equil <- lapply(1:nspecies, function(i) {
loga.equil[[i]] - log10(m.balance[i])
})
Modified: pkg/CHNOSZ/R/wjd.R
===================================================================
--- pkg/CHNOSZ/R/wjd.R 2015-01-07 13:06:31 UTC (rev 74)
+++ pkg/CHNOSZ/R/wjd.R 2015-02-15 06:42:50 UTC (rev 75)
@@ -264,7 +264,6 @@
if(!"limSolve" %in% row.names(installed.packages())) {
msgout("guess: skipping 'central' method as limSolve package is not available\n")
} else {
- require(limSolve)
# the inequality constraints for moles of species
G <- diag(nrow(A))
# minX is the minimum mole number we will accept
Modified: pkg/CHNOSZ/inst/NEWS
===================================================================
--- pkg/CHNOSZ/inst/NEWS 2015-01-07 13:06:31 UTC (rev 74)
+++ pkg/CHNOSZ/inst/NEWS 2015-02-15 06:42:50 UTC (rev 75)
@@ -1,4 +1,4 @@
-CHANGES IN CHNOSZ 1.0.3-12 (2015-01-07)
+CHANGES IN CHNOSZ 1.0.3-13 (2015-02-15)
---------------------------------------
- Add files with average amino acid compositions of proteins from Bison
@@ -32,6 +32,19 @@
some resolutions (i.e. blank lines), and greatly reduces the size
of PDF files.
+- Add 'as.residue' argument to diagram(), to divide reactions by
+ balance coefficients without rescaling to whole species formulas
+ (cf. 'normalize' argument which uses rescaling).
+
+- In equilibrate(), change name of 'stay.normal' argument to
+ 'as.residue'.
+
+- In test-diagram.R, test that same diagrams are produced using
+ 'normalize' in equilibrate() or diagram().
+
+- In test-diagram.R, test that same diagrams are produced using
+ 'as.residue' in equilibrate() or diagram().
+
CHANGES IN CHNOSZ 1.0.3 (2014-01-12)
------------------------------------
Modified: pkg/CHNOSZ/inst/tests/test-diagram.R
===================================================================
--- pkg/CHNOSZ/inst/tests/test-diagram.R 2015-01-07 13:06:31 UTC (rev 74)
+++ pkg/CHNOSZ/inst/tests/test-diagram.R 2015-02-15 06:42:50 UTC (rev 75)
@@ -6,10 +6,8 @@
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, normalize=TRUE, plot.it=FALSE), "normalizing formulas in calculation of predominant species")
e <- equilibrate(a)
expect_error(diagram(e, "Z"), "Z is not a basis species")
- expect_error(diagram(e, normalize=TRUE), "normalizing formulas is only possible if 'eout' is the output from affinity\\(\\)")
})
test_that("expected messages, errors and results arise using output from affinity()", {
@@ -42,7 +40,7 @@
expect_equal(a$values[[4]], array(numeric(128)))
})
-test_that("'groups', 'alpha' and 'normalize' work as expected", {
+test_that("'groups' and 'alpha' work as expected", {
basis("CHNOS+")
species(c("formic acid", "formate", "acetic acid", "acetate"))
# 1-D
@@ -57,15 +55,27 @@
d <- diagram(e, alpha=TRUE, plot.it=FALSE)
# we should find that the sum of alphas is one
expect_equal(Reduce("+", d$plotvals), array(rep(1, 128)))
-# # normalize multiplies the activities by the balancing coefficients
-# d <- diagram(e)
-# d.norm <- diagram(e, normalize=TRUE)
-# # formic acid, 1 CO2 in formation reaction
-# expect_equal(d$plotvals[[1]], d.norm$plotvals[[1]])
-# # acetic acid, 2 CO2 in formation reaction
-# expect_equal(d$plotvals[[3]], d.norm$plotvals[[3]] - log10(2))
})
+test_that("'normalize' and 'as.residue' work as expected", {
+ basis("CHNOS")
+ species(c("LYSC_CHICK", "MYG_PHYCA", "RNAS1_BOVIN", "CYC_BOVIN"))
+ # 1-D
+ a <- affinity(O2=c(-80, -70))
+ expect_error(diagram(a, normalize=TRUE), "can be TRUE only for a 2-D \\(predominance\\) diagram")
+ # 2-D
+ a <- affinity(H2O=c(-10, 0), O2=c(-80, -70))
+ d1 <- diagram(a, normalize=TRUE, plot.it=FALSE)
+ e <- equilibrate(a, normalize=TRUE)
+ d2 <- diagram(e, plot.it=FALSE)
+ expect_equal(d1$predominant, d2$predominant)
+ expect_error(diagram(e, normalize=TRUE), "can be TRUE only if 'eout' is the output from affinity\\(\\)")
+ d3 <- diagram(a, as.residue=TRUE, plot.it=FALSE)
+ e <- equilibrate(a, as.residue=TRUE)
+ d4 <- diagram(e, plot.it=FALSE)
+ expect_equal(d3$predominant, d4$predominant)
+})
+
test_that("diagram() handles 2D plots with different x and y resolution and warns for >1 species in contour plot", {
basis("CHNOS")
species(c("alanine", "glycine", "serine", "methionine"))
Modified: pkg/CHNOSZ/inst/tests/test-equilibrate.R
===================================================================
--- pkg/CHNOSZ/inst/tests/test-equilibrate.R 2015-01-07 13:06:31 UTC (rev 74)
+++ pkg/CHNOSZ/inst/tests/test-equilibrate.R 2015-02-15 06:42:50 UTC (rev 75)
@@ -31,7 +31,7 @@
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, normalize=TRUE), "normalizing molar formulas by the balancing coefficients")
+ expect_message(equilibrate(aprot, normalize=TRUE), "using 'normalize' for molar formulas")
})
test_that("equilibrate() gives expected messages and errors for species selection", {
Modified: pkg/CHNOSZ/man/diagram.Rd
===================================================================
--- pkg/CHNOSZ/man/diagram.Rd 2015-01-07 13:06:31 UTC (rev 74)
+++ pkg/CHNOSZ/man/diagram.Rd 2015-02-15 06:42:50 UTC (rev 75)
@@ -9,13 +9,14 @@
\usage{
diagram(eout, what = "loga.equil", alpha = FALSE, normalize = FALSE,
- balance=NULL, groups=as.list(1:length(eout$values)), xrange=NULL,
- mar=NULL, yline=par("mgp")[1]+0.7, side=1:4,
+ as.residue = FALSE, balance=NULL, groups=as.list(1:length(eout$values)),
+ xrange=NULL, mar=NULL, yline=par("mgp")[1]+0.7, side=1:4,
ylog=TRUE, xlim=NULL, ylim=NULL, xlab=NULL, ylab=NULL,
cex=par("cex"), cex.names=1, cex.axis=par("cex"),
lty=NULL, lwd=par("lwd"), dotted=0,
bg=par("bg"), col=par("col"), col.names=par("col"), fill=NULL,
- names=NULL, main=NULL, legend.x="topright", add=FALSE, plot.it=TRUE)
+ names=NULL, main=NULL, legend.x="topright", add=FALSE, plot.it=TRUE,
+ tplot=TRUE)
strip(affinity, ispecies = NULL, col = NULL, ns = NULL,
xticks = NULL, ymin = -0.2, xpad = 1, cex.names = 0.7)
find.tp(x)
@@ -25,7 +26,8 @@
\item{eout}{list, object returned by \code{\link{equilibrate}}} or \code{\link{affinity}}
\item{what}{character, what property to calculate and plot}
\item{alpha}{logical, for speciation diagrams, plot degree of formation instead of activities?}
- \item{normalize}{logical, normalize chemical formulas or reaction properties by the balance vector?}
+ \item{normalize}{logical, divide chemical affinities by balance coefficients (rescale to whole formulas)?}
+ \item{as.residue}{logical, divide chemical affinities by balance coefficients (no rescaling)?}
\item{balance}{character, balancing constraint; see \code{\link{balance}}}
\item{groups}{list of numeric, groups of species to consider as a single effective species}
\item{xrange}{numeric, range of x-values between which predominance field boundaries are plotted}
@@ -52,6 +54,7 @@
\item{legend.x}{character, description of legend placement passed to \code{\link{legend}}}
\item{add}{logical, add to current plot?}
\item{plot.it}{logical, make a plot?}
+ \item{tplot}{logical, set up plot with \code{\link{thermo.plot.new}}?}
\item{affinity}{list, object returned by \code{\link{affinity}}}
\item{ispecies}{numeric, which species to consider (default of \code{NULL} is to consider all species)}
\item{ns}{numeric, numbers of species, used to make inset plots for strip diagrams}
@@ -78,11 +81,15 @@
If \code{legend.x} is NA, instead of a legend, the lines are labeled with the names of the species near the maximum value.
The line type and line width can be controlled with \code{lty} and \code{lwd}, respectively.
-On 2-D diagrams, the fields represent the species with the highest equilibrium activity, after any \code{normalize} or \code{alpha} operations.
+On 2-D diagrams, the fields represent the species with the highest equilibrium activity.
\code{fill} determines the color of the predominance fields, \code{col} that of the boundary lines.
By default, \code{\link{heat.colors}} are used to fill the predominance fields in diagrams on the screen plot device.
The style of the boundary lines on 2-D diagrams can be altered by supplying one or more non-zero integers in \code{dotted}, which indicates the fraction of line segments to omit; a value of \samp{1} or NULL for \code{dotted} has the effect of not drawing the boundary lines.
+\code{normalize} and \code{as.residue} apply only to the 2-D diagrams, and only when \code{eout} is the output from \code{affinity}.
+With \code{normalize}, the activity boundaries are calculated as between the residues of the species (the species divided by the balance coefficients), then the activities are rescaled to the whole species formulas.
+With \code{as.residue}, the activity boundaries are calculated as between the residues of the species, and no rescaling is performed.
+
For all diagrams, the \code{names} of the species and their colors in \code{col.names} can be changed, as can \code{cex}, \code{cex.names}, and \code{cex.axis} to adjust the overall character expansion factors (see \code{\link{par}}) and those of the species names and axis labels.
The x- and y-axis labels are automatically generated unless they are supplied in \code{xlab} and \code{ylab}.
A new plot is started unless \code{add} is TRUE.
@@ -97,7 +104,7 @@
If \code{what} is missing, option (1) is selected if the number of dimensions is 0 or 1, and option (3) is selected if the number of dimensions is 2.
The latter is referred to as the maximum affinity method.
In cases where it applies (see Warning), the maximum affinity method is much faster than an equilibrium calculation.
-\code{balance} is the option, sent to \code{\link{balance}}, that determines the balancing coefficients used in the normalization (this argument has no effect on the calculations of equilibrium activities.)
+\code{balance} is the option, sent to \code{\link{balance}}, that determines the balancing coefficients (this argument has no effect on the calculations of equilibrium activities.)
If \code{what} is the name of a basis species, it refers to option (2) above.
A contour plot is made in the case of 2-D diagrams of the equilibrium activity of a basis species (see \code{\link{demos}("CO2Ac")}, and only the first species of interest is used in the calculation; a warning is produced if there is more than one.
Modified: pkg/CHNOSZ/man/equilibrate.Rd
===================================================================
--- pkg/CHNOSZ/man/equilibrate.Rd 2015-01-07 13:06:31 UTC (rev 74)
+++ pkg/CHNOSZ/man/equilibrate.Rd 2015-02-15 06:42:50 UTC (rev 75)
@@ -10,7 +10,7 @@
\usage{
equilibrate(aout, balance=NULL, loga.balance=NULL,
- ispecies=1:length(aout$values), normalize=FALSE, stay.normal=FALSE)
+ ispecies=1:length(aout$values), normalize=FALSE, as.residue=FALSE)
balance(aout, balance=NULL)
equil.boltzmann(Astar, n.balance, loga.balance)
equil.reaction(Astar, n.balance, loga.balance)
@@ -21,7 +21,7 @@
\item{balance}{character or numeric, how to balance the transformations}
\item{ispecies}{numeric, which species to include}
\item{normalize}{logical, normalize the molar formulas of species by the balancing coefficients?}
- \item{stay.normal}{logical, report results for the normalized formulas?}
+ \item{as.residue}{logical, report results for the normalized formulas?}
\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}
@@ -54,7 +54,7 @@
This operation is intended for systems of polymers, such as proteins, whose conventional formulas are much larger than the basis speices.
The normalization also applies to the balancing coefficients, which as a result consist of \samp{1}s.
\code{normalize} has the same effect as did \code{diagram(..., residue=TRUE)} in versions of CHNOSZ before 0.9-9.
-After normalization and equilibration, the equilibrium activities are then un-normalized (for the original formulas of the species), unless \code{stay.normal} is TRUE.
+After normalization and equilibration, the equilibrium activities are then re-scaled (for the original formulas of the species), unless \code{as.residue} is TRUE.
\code{ispecies} can be supplied to identify a subset of the species to include in the calculation.
Modified: pkg/CHNOSZ/vignettes/hotspring.Rnw
===================================================================
--- pkg/CHNOSZ/vignettes/hotspring.Rnw 2015-01-07 13:06:31 UTC (rev 74)
+++ pkg/CHNOSZ/vignettes/hotspring.Rnw 2015-02-15 06:42:50 UTC (rev 75)
@@ -1,4 +1,4 @@
-%% LyX 2.1.0beta2 created this file. For more info, see http://www.lyx.org/.
+%% LyX 2.1.3 created this file. For more info, see http://www.lyx.org/.
%% Do not edit unless you really know what you are doing.
\documentclass[english]{article}
\usepackage{mathpazo}
@@ -473,7 +473,7 @@
basis("pH", bison.pH[i])
# calculate metastable equilibrium activities of the residues
a <- affinity(H2=c(-11, -1, 101), T=bison.T[i], iprotein=ip.phyla[iloc])
- e <- equilibrate(a, loga.balance=0, normalize=TRUE, stay.normal=TRUE)
+ e <- equilibrate(a, loga.balance=0, as.residue=TRUE)
# remove the logarithms to get relative abundances
a.residue <- 10^sapply(e$loga.equil, c)
colnames(a.residue) <- aa.phyla$organism[iloc]
Modified: pkg/CHNOSZ/vignettes/hotspring.lyx
===================================================================
--- pkg/CHNOSZ/vignettes/hotspring.lyx 2015-01-07 13:06:31 UTC (rev 74)
+++ pkg/CHNOSZ/vignettes/hotspring.lyx 2015-02-15 06:42:50 UTC (rev 75)
@@ -1550,6 +1550,7 @@
\begin_layout Plain Layout
+
\backslash
clearpage
\end_layout
@@ -1691,6 +1692,7 @@
\begin_layout Plain Layout
+
\backslash
clearpage
\end_layout
@@ -1994,7 +1996,7 @@
\begin_layout Plain Layout
- e <- equilibrate(a, loga.balance=0, normalize=TRUE, stay.normal=TRUE)
+ e <- equilibrate(a, loga.balance=0, as.residue=TRUE)
\end_layout
\begin_layout Plain Layout
@@ -2268,6 +2270,7 @@
\begin_layout Plain Layout
+
\backslash
clearpage
\end_layout
@@ -2953,6 +2956,7 @@
\begin_layout Plain Layout
+
\backslash
clearpage
\end_layout
Modified: pkg/CHNOSZ/vignettes/wjd.Rnw
===================================================================
--- pkg/CHNOSZ/vignettes/wjd.Rnw 2015-01-07 13:06:31 UTC (rev 74)
+++ pkg/CHNOSZ/vignettes/wjd.Rnw 2015-02-15 06:42:50 UTC (rev 75)
@@ -1,4 +1,4 @@
-%% LyX 2.0.6 created this file. For more info, see http://www.lyx.org/.
+%% LyX 2.1.3 created this file. For more info, see http://www.lyx.org/.
%% Do not edit unless you really know what you are doing.
\documentclass[english,noae,round]{article}
\usepackage{mathpazo}
@@ -71,7 +71,6 @@
@
-
\section{Is it winding?}
The ``winding'' in the title refers to the observation that the
@@ -87,6 +86,7 @@
<<wjd>>=
w <- wjd()
@
+
What are the most abundant species, and how many iterations were taken?
<<abundant>>=
@@ -98,6 +98,7 @@
niter <- length(w$lambda)
niter
@
+
<<abundantnames,echo=FALSE>>=
# the formulas of the two most abundant species
# use CHNOSZ's as.chemical.formula function, after
@@ -106,6 +107,7 @@
f2 <- as.chemical.formula(w$A[oX[2],w$A[oX[2],]!=0])
f4 <- as.chemical.formula(w$A[oX[4],w$A[oX[4],]!=0])
@
+
We find that \Sexpr{f1} and \Sexpr{f2} are the most abundant species,
after \Sexpr{niter} iterations. The fourth most-abundant species,
\Sexpr{f4}, is discussed below.
@@ -121,10 +123,13 @@
return(x)
}
@
+
Then apply the function over the different numbers of iterations,
from 0 (initial conditions) to \Sexpr{niter}, and plot the values.
\setkeys{Gin}{width=1.0\textwidth}
+
+<<>>=
<<iterplot, fig=TRUE, height=3>>=
par(mfrow=c(1, 2))
sa <- sapply(0:niter, iterfun, i=oX[1])
@@ -132,9 +137,9 @@
sa <- sapply(0:niter, iterfun, i=oX[4])
plot(0:niter, sa, xlab="iteration", ylab=paste("x", f4))
@
+
\setkeys{Gin}{width=1.0\textwidth}
-
A bit of winding: the mole fraction of \Sexpr{f1} generally increases
(there is a very slight decrease at the fifth iteration), but the
mole fraction of \Sexpr{f4} increases in the first two iterations,
@@ -146,6 +151,7 @@
<<Gdown>>=
all(diff(w$F.Y) < 0)
@
+
The decrease in Gibbs energy becomes smaller with every iteration,
as an equilibrium distribution of species is approached:
@@ -154,7 +160,6 @@
@
-
\section{Is it near the bottom? Testing for equilibrium}
@@ -172,8 +177,9 @@
are the fractional changes with each iteration:
<<Gfrac>>=
-diff(w$F.Y)/w$F.Y[1:6]
+diff(w$F.Y)/w$F.Y[1:7]
@
+
Only the last value is below the default value of \texttt{Gfrac} in
\texttt{wjd()}.
@@ -197,23 +203,27 @@
first look at a plot showing the results after 3 iterations.
\setkeys{Gin}{width=0.6\textwidth}
+
+<<>>=
<<w3_ep_plot, fig=TRUE>>=
w3 <- wjd(imax=3)
ep3 <- element.potentials(w3, plot.it=TRUE)
@
+
\setkeys{Gin}{width=1.0\textwidth}
-
Here's the plot for the default settings of \texttt{wjd()} which takes
6 iterations:
\setkeys{Gin}{width=0.6\textwidth}
+
+<<>>=
<<w_ep_plot, fig=TRUE>>=
ep <- element.potentials(w, plot.it=TRUE)
@
+
\setkeys{Gin}{width=1.0\textwidth}
-
That one shows considerably less deviation than the first plot. The
differences don't become zero, but perhaps they are small enough to
accept as an operational solution. There is a function, \texttt{is.near.equil()},
@@ -232,7 +242,6 @@
@
-
\subsection{Different initial solutions yield similar results}
@@ -268,18 +277,21 @@
#Y2
@
-
Now let's run the minimization using the different guesses, count
the numbers of iterations, and test them for equilibrium.
\setkeys{Gin}{width=0.6\textwidth}
+
+<<>>=
<<stoich_guess>>=
wY1 <- wjd(Y=Y1)
niterY1 <- length(wY1$lambda)
niterY1
is.near.equil(wY1, tol=0.0001)
@
+
\setkeys{Gin}{width=0.6\textwidth}
+
Using the first guess generated using the ``stoich'' method there
were \Sexpr{niterY1-niter} more iterations than using the default
initial solution in \texttt{wjd()}.
@@ -289,14 +301,16 @@
Do the different initial guesses actually give similar results?
\setkeys{Gin}{width=0.6\textwidth}
+
+<<>>=
<<compare_guesses, fig=TRUE>>=
plot(1:10, w$X, xlab="species", ylab="mole fraction")
points(1:10, wY1$X, pch=0)
#points(1:10, wY2$X, pch=2)
@
+
\setkeys{Gin}{width=0.6\textwidth}
-
At that scale it's fairly difficult to tell them apart.
@@ -317,12 +331,14 @@
nguess <- length(iYs)
nguess
@
+
Now let's run the minimization using each of those guesses, and test
if each one \texttt{is.near.equil()}:
<<allguess.equil>>=
sapply(iYs,function(i) is.near.equil(wjd(Y=Ys[[i]])))
@
+
Looks like they are. What if we're really picky though and want to
make sure the chemical potentials of the elements are very well-defined?
We can decrease the tolerance of \texttt{is.near.equil()}:
@@ -330,6 +346,7 @@
<<allguess.equil.2>>=
sapply(iYs,function(i) is.near.equil(wjd(Y=Ys[[i]]),tol=0.0001))
@
+
Well some of the tests failed. But maybe this is because we have stopped
iterating too soon. Let's iterate until a we have smaller change in
Gibbs energy of the system:
@@ -339,6 +356,7 @@
is.near.equil(wjd(Y=Ys[[i]], Gfrac=1e-9), tol=0.0001)
})
@
+
We're back to passing all but one of the equilibrium tests based on
uniformity of chemical potentials of the elements. What are the standard
deviations of the resulting species abundances?
@@ -347,6 +365,7 @@
Xs <- sapply(iYs, function(i) wjd(Y=Ys[[i]], Gfrac=1e-9)$X)
apply(Xs, 1, sd)
@
+
Based on this exercise, starting from different initial species abundances
for the same bulk composition, and arriving at a similar near-equilibrium
solution, it seems likely that we are near the global minimum.
@@ -372,6 +391,7 @@
dlen <- read.csv(file, stringsAsFactors=FALSE, row.names=1)
t(dlen[, 1, drop=FALSE])
@
+
Now let's write some code to define the system. The bulk composition
in C:H:O space has 40 percent oxygen, a percentage of carbon given
by \texttt{xC}, and hydrogen as the remainder; the C:N ratio is taken
@@ -395,6 +415,7 @@
return(w)
}
@
+
Notice the increase in \texttt{imax} and the decrease in \texttt{Gfrac}
from the defaults! These changes are needed to get closer to equilibrium
(and we'll know we're not there if a message is shown). What does
@@ -404,6 +425,7 @@
<<atmos_C15>>=
sort(min.atmos(15)$X, decreasing=TRUE)
@
+
Well that's exciting! The order of abundances of $\mathrm{H_{2}O}$,
$\mathrm{CO_{2}}$, $\mathrm{CH_{4}}$, $\mathrm{NH_{3}}$, $\mathrm{CO}$,
ethane ($\mathrm{C_{2}H_{6}}$), formic acid ($\mathrm{CH_{2}O_{2}}$)
@@ -421,7 +443,6 @@
text(7, log10(Xs[, 1]), dlen$formula, adj=1)
@
-
We triggered some ``not near equilibrium'' messages, but overall
it seems well behaved (note that using the ``central'' method in
\texttt{guess()} is one way that leads to fewer of these messages).
@@ -457,7 +478,6 @@
iaro <- info(aromatics, "liq")
@
-
Let's find the equilibrium distribution of species for a bulk composition
corresponding to a single mole of each species.
@@ -467,7 +487,6 @@
is.near.equil(waa)
@
-
The \texttt{waa\$elements} shows that the bulk chemical composition
in the Gibbs energy minimization was $\mathrm{C_{84}H_{106}}$ (H/C
= 1.261905). That FALSE means the chemical potentials of the elements
@@ -482,7 +501,6 @@
is.near.equil(waa)
@
-
Let's plot the abundances of the species:
<<barplot_aa, fig=TRUE>>=
@@ -490,7 +508,6 @@
text(bp, rep(0.2,8), c(alkanes, aromatics), srt=90, adj=0)
@
-
Now let's calculate the chemical potentials of the elements (dimensionless
values, $\mu/RT$) ...
@@ -499,7 +516,6 @@
print(ep)
@
-
... and the corresponding logarithms of chemical activities of a set
of basis species:
@@ -509,7 +525,6 @@
@
-
\subsection{Bulk composition instead of moles of species}
\texttt{run.wjd()} can be given B (bulk chemical composition), instead
@@ -524,7 +539,6 @@
text(bp, rep(0.2,8), c(alkanes, aromatics), srt=90, adj=0)
@
-
Compared with the previous example, here the lower H/C = 0.7 defines
a more oxidized system, so the increase in aromatic content is expected.
As before, we can calculate the chemical potentials of the elements
@@ -536,7 +550,6 @@
@
-
\section{Document Information}
Revision history
@@ -551,7 +564,6 @@
sessionInfo()
@
-
\bibliographystyle{plainnat}
\bibliography{vig}
Modified: pkg/CHNOSZ/vignettes/wjd.lyx
===================================================================
--- pkg/CHNOSZ/vignettes/wjd.lyx 2015-01-07 13:06:31 UTC (rev 74)
+++ pkg/CHNOSZ/vignettes/wjd.lyx 2015-02-15 06:42:50 UTC (rev 75)
@@ -1,5 +1,5 @@
-#LyX 2.0 created this file. For more info see http://www.lyx.org/
-\lyxformat 413
+#LyX 2.1 created this file. For more info see http://www.lyx.org/
+\lyxformat 474
\begin_document
\begin_header
\textclass article
@@ -22,13 +22,13 @@
\font_roman palatino
\font_sans default
\font_typewriter default
+\font_math auto
\font_default_family default
\use_non_tex_fonts false
\font_sc false
\font_osf false
\font_sf_scale 100
\font_tt_scale 100
-
\graphics default
\default_output_format default
\output_sync 0
@@ -51,15 +51,24 @@
\pdf_quoted_options "citecolor=blue"
\papersize letterpaper
\use_geometry true
-\use_amsmath 1
-\use_esint 1
-\use_mhchem 1
-\use_mathdots 1
-\cite_engine natbib_authoryear
+\use_package amsmath 1
+\use_package amssymb 1
+\use_package cancel 1
+\use_package esint 1
+\use_package mathdots 1
+\use_package mathtools 1
+\use_package mhchem 1
+\use_package stackrel 1
+\use_package stmaryrd 1
+\use_package undertilde 1
+\cite_engine natbib
+\cite_engine_type authoryear
+\biblio_style plainnat
\use_bibtopic false
\use_indices false
\paperorientation portrait
\suppress_date false
+\justification true
\use_refstyle 0
\branch short
\selected 1
@@ -230,24 +239,32 @@
\begin_inset Branch short
status open
-\begin_layout Chunk
+\begin_layout Standard
+\begin_inset Flex Chunk
+status open
-<<libraryCHNOSZ, echo=FALSE>>=
+\begin_layout Plain Layout
+
+\begin_inset Argument 1
+status open
+
+\begin_layout Plain Layout
+libraryCHNOSZ, echo=FALSE
\end_layout
-\begin_layout Chunk
+\end_inset
library(CHNOSZ)
\end_layout
-\begin_layout Chunk
+\begin_layout Plain Layout
data(thermo)
\end_layout
-\begin_layout Chunk
+\end_inset
-@
+
\end_layout
\end_inset
@@ -297,23 +314,35 @@
\begin_inset Branch short
status open
-\begin_layout Chunk
+\begin_layout Standard
+\begin_inset Flex Chunk
+status open
-<<wjd>>=
+\begin_layout Plain Layout
+
+\begin_inset Argument 1
+status open
+
+\begin_layout Plain Layout
+wjd
\end_layout
-\begin_layout Chunk
+\end_inset
w <- wjd()
\end_layout
-\begin_layout Chunk
+\end_inset
-@
+
\end_layout
\end_inset
+
+\end_layout
+
+\begin_layout Standard
What are the most abundant species, and how many iterations were taken?
\end_layout
@@ -321,102 +350,125 @@
\begin_inset Branch short
status open
-\begin_layout Chunk
+\begin_layout Standard
+\begin_inset Flex Chunk
+status open
-<<abundant>>=
+\begin_layout Plain Layout
+
+\begin_inset Argument 1
+status open
+
+\begin_layout Plain Layout
+abundant
\end_layout
-\begin_layout Chunk
+\end_inset
# the order of species abundance
\end_layout
-\begin_layout Chunk
+\begin_layout Plain Layout
oX <- order(w$X, decreasing=TRUE)
\end_layout
-\begin_layout Chunk
+\begin_layout Plain Layout
# the stoichiometries of the two most abundant species
\end_layout
-\begin_layout Chunk
+\begin_layout Plain Layout
w$A[head(oX,2),]
\end_layout
-\begin_layout Chunk
+\begin_layout Plain Layout
# the number of iterations
\end_layout
-\begin_layout Chunk
+\begin_layout Plain Layout
niter <- length(w$lambda)
\end_layout
-\begin_layout Chunk
+\begin_layout Plain Layout
niter
\end_layout
-\begin_layout Chunk
+\end_inset
-@
+
\end_layout
\end_inset
+\end_layout
+
+\begin_layout Standard
\begin_inset Branch short
status open
-\begin_layout Chunk
+\begin_layout Standard
+\begin_inset Flex Chunk
+status open
-<<abundantnames,echo=FALSE>>=
+\begin_layout Plain Layout
+
+\begin_inset Argument 1
+status open
+
+\begin_layout Plain Layout
+abundantnames,echo=FALSE
\end_layout
-\begin_layout Chunk
+\end_inset
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/chnosz -r 75
More information about the CHNOSZ-commits
mailing list