revisit: Plots and Optima of Objective Functions

Description

Calculate values of an objective function from logarithms of activities of chemical species and (for some objectives) reference logarithms of activity. Make line or contour plots showing the values of the objective function and the positions of the optima (minimum or maximum).

Usage

revisit(eout, objective = "CV", loga2 = 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")
  optimal.index(z, objective)
  extremes(z, objective)

Arguments

eout

list, output from equilibrate, containing logarithms of activities of species

objective

character, name of objective function

loga2

numeric array, reference values of logarithm of activities

ispecies

numeric, which species to consider

col

character, color to use for points or lines

yline

numeric, margin line for y-axis label

ylim

numeric, limits of y axis

cex

numeric, character expansion factor

lwd

numeric, line width

mar

numeric, plot margin specifications

side

numeric, which sides of plot to draw axes

xlim

numeric, limits of x axis

labcex

numeric, character expansion factor for species labels

pch

numeric, plotting symbol(s) to use for points

main

character, main title for plot

plot.it

logical, make a plot?

add

logical, add to an existing plot?

plot.optval

logical, show the location of the optimal value(s)?

style.2D

character, type of 2-D plot

numeric, matrix of values (H element of output from revisit)

Value

revisit returns a list containing at least an element named H, giving the calculated values of the objective function. For 1 or 2 dimensions of varibility of chemical conditions, the output also contains the elements ixopt and iyopt (1-D and 2-D) and iyopt and yopt (2-D) indicating the positions and values of the optimum. For calculations in more than two dimensions, the output contains iopt which is a matrix
optimum.index returns the index (or indices) of the optimal values in z using array(..., arr.ind=TRUE). extremes returns $x$ values that are the column numbers where the optimum is found for each row, and $y$ values that are the row numbers where the optimum is found for each column. The optimum attribute of the objective function indicates whether minimal or maximal values are used.

Details

Most commonly, the output of equilibrate is used as the value for eout. The type of calculation is indicated by objective, giving the name of an objective function. The list of logarithms of chemical activities in eout$loga.equil is passed as loga1 to the objective function. If the objective function has an argument a1 instead of loga1, the activities instead of their logarithms are passed. Reference logarithms of activities, for those objective functions that require them, are supplied in the loga2 argument. Generally, loga2 must be a numeric vector with length equal to that of loga1 (i.e., number of species). If a single numeric value is supplied for loga2, it is recycled to the length of loga1.

For calculations at a single condition (0-D, no variation), with the qqr objective, a quantile-quantile plot (qqnorm) is shown. For rmsd and other objective functions having reference values (loga2), a scatter plot is shown with a smooth line calculated using loess.smooth. The line can be suppressed using lwd=NULL. Otherwise, no plot is made for 0-D calculations for the other objective functions.

If plot.it is TRUE, and eout is the output from equilibrate, and the number of variables is 1 or 2, the results are plotted --- a line diagram in 1 dimension or a contour plot in 2 dimensions. style.2D can be set to image to fill the plot with colors instead of the contour plot that is the default.

If plot.optval is TRUE, the location of the optimum (or optima) is indicated by a dashed vertical line(s) on a 1-D plot or a point(s) marked by an asterisk on a 2-D plot. Also, on 2-D plots, the locations of the optima at each grid line perpendicular to the $x$ and $y$ axes are plotted. These points follow major ridges or valleys, and are plotted as dashed lines colored green for the $x$ and blue for the $y$ values returned by extremes.

An alternative source for the eout argument is any list of numeric values, each element of which corresponds to a different observation (such as a single species), all having the same dimensions (as vectors, matrices or higher-dimensional arrays) In this case, plotting is disabled, since the names of the variables are not in the input.

revisit is a partial anagram of diversity, which was the provisional name of the function but was changed in CHNOSZ-0.9. While the diversity function (in vegan) operates on a matrix with (biological) species on the columns, revisit operates on a list with (chemical) species as the elements of the list. The name of the H output value is the conventional symbol for the Shannon diversity index, which was the first target statistic to be implemented in revisit.

Examples

Run this code

data(thermo)
## example of defining a new objective function
# count the species with logarithms of activity greater than loga2
count <- function(loga1, loga2) rowSums(loga1 > loga2)
# set the attribute indicating the type of optimum
attr(count, "optimum") <- "maximal"
# equilibrate a system of amino acids
basis("CHNOS")
species(aminoacids(""))
a <- affinity(O2=c(-80, -60))
e <- equilibrate(a)
# make a plot
r <- revisit(e, "count", -5)
title(main="amino acids with metastable log activities > -5")

# can also make a 2-D plot
a <- affinity(O2=c(-74, -60), H2O=c(-3, 3))
e <- equilibrate(a)
r <- revisit(e, "count", -5, style.2D="image", plot.optval=FALSE)
title(main="amino acids with metastable log activities > -5")

## 'revisit' calculations for amino acids
opar <- par(mfrow=c(2, 2))
basis("CHNOS+")
species(aminoacids(""))
# chemical affinities as a function of logarithm of oxygen fugacity
a <- affinity(O2=c(-85, -60))
# shows the equilibrium abundances of the amino acids
e <- equilibrate(a)
diagram(e)
mtitle(c("20 amino acids", "balanced on CO2"))
# show a legend with input constraints
db <- describe.basis(ibasis=3)
dp <- describe.property("T", 25)
legend("bottomright", c(dp, db))
# default is to plot coefficient of variation
r <- revisit(e)
# show a title with the optimal conditions
mincv <- format(r$optimum, digits=3)
t1 <- paste("minimum coeff of variation,", mincv, "at:")
# the logfO2 that minimized the C.V.
basis("O2", r$x)
t2 <- describe.basis(ibasis=5)
mtitle(c(t1, as.expression(t2)))
# chemical affinities as a function of two other variables
a <- affinity(NH3=c(-10, 10, 40), T=c(0, 80, 40))
diagram(a, fill="heat")
# show a legend with input constraints
db <- describe.basis(ibasis=5)
legend("bottomright", as.expression(db))
# contour plot of the CV
e <- equilibrate(a)
r <- revisit(e)
# show a title with the optimal conditions
mincv <- format(r$optimum, digits=3)
t1 <- paste("minimum coeff of variation,", mincv, "at:")
# the logaNH3 and T that minimized the C.V.
basis("NH3", r$x)
db <- describe.basis(ibasis=3)
dp <- describe.property("T", r$y)
t2 <- substitute(list(dp, db), list(dp=dp[[1]], db=db[[1]]))
mtitle(c(t1, as.expression(t2)))
par(opar)

## calculations for proteins in Pelagibacter ubique
## using grep.file, read.fasta, add.protein
f <- system.file("extdata/fasta/HTCC1062.faa.xz", package="CHNOSZ")
# what proteins to select (set to "" for all proteins)
w <- "ribosomal"
# locate entries whose names contain w
j <- grep.file(f, w)
# get the amino acid compositions of these protein
aa <- read.fasta(f, j)
# add these proteins to CHNOSZ's inventory
ip <- add.protein(aa)
# set up a the chemical system
basis("CHNOS+")
# calculate affinities of formation in logfO2 space
a <- affinity(O2=c(-85, -60), iprotein=ip)
# show the equilibrium activities
opar <- par(mfrow=c(2, 2))
e <- equilibrate(a, loga.balance=0, normalize=TRUE)
diagram(e, names=NULL)
# make a title
expr <- as.expression(substitute(x~y~"proteins in"~
  italic("P. ubique"), list(x=length(j), y=w)))
mtitle(c("Equilibrium activities of", expr))
# show the coefficient of variation
revisit(e, "CV")
mtitle(c("CV of equilibrium activities of", expr))
# calculate affinities in logfO2-logaH2O space
a <- affinity(O2=c(-85,-65), H2O=c(-10,0), iprotein=ip)
# show the predominances
diagram(a, normalize=TRUE, fill="heat")
# calculate the equilibrium activities
e <- equilibrate(a, loga.balance=0, normalize=TRUE)
# show the coefficient of variation
r <- revisit(e, "CV")
stopifnot(r$ix==37)
stopifnot(r$iy==53)
mtitle(c("CV of equilibrium activities of", expr))
par(opar)