[CHNOSZ-commits] r25 - in pkg/CHNOSZ: . R inst inst/doc inst/tests man vignettes
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Oct 2 16:15:54 CEST 2012
Author: jedick
Date: 2012-10-02 16:15:53 +0200 (Tue, 02 Oct 2012)
New Revision: 25
Added:
pkg/CHNOSZ/R/equilibrate.R
pkg/CHNOSZ/inst/doc/equilibrium.pdf
pkg/CHNOSZ/inst/doc/hotspring.pdf
pkg/CHNOSZ/inst/tests/test-equilibrate.R
pkg/CHNOSZ/man/equilibrate.Rd
pkg/CHNOSZ/vignettes/equilibrium.Rnw
pkg/CHNOSZ/vignettes/equilibrium.lyx
pkg/CHNOSZ/vignettes/flowchart.dia
pkg/CHNOSZ/vignettes/flowchart.pdf
pkg/CHNOSZ/vignettes/hotspring.Rnw
pkg/CHNOSZ/vignettes/hotspring.lyx
Removed:
pkg/CHNOSZ/R/equil.R
pkg/CHNOSZ/inst/doc/hs-chemistry.pdf
pkg/CHNOSZ/man/equil.Rd
pkg/CHNOSZ/vignettes/chnosz_new.pdf
pkg/CHNOSZ/vignettes/hs-chemistry.Rnw
pkg/CHNOSZ/vignettes/hs-chemistry.lyx
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/affinity.R
pkg/CHNOSZ/R/anim.R
pkg/CHNOSZ/R/basis.R
pkg/CHNOSZ/R/buffer.R
pkg/CHNOSZ/R/diagram.R
pkg/CHNOSZ/R/examples.R
pkg/CHNOSZ/R/findit.R
pkg/CHNOSZ/R/protein.info.R
pkg/CHNOSZ/R/revisit.R
pkg/CHNOSZ/R/transfer.R
pkg/CHNOSZ/R/util.affinity.R
pkg/CHNOSZ/R/util.expression.R
pkg/CHNOSZ/R/util.list.R
pkg/CHNOSZ/R/util.program.R
pkg/CHNOSZ/inst/NEWS
pkg/CHNOSZ/inst/doc/anintro.pdf
pkg/CHNOSZ/inst/doc/wjd.pdf
pkg/CHNOSZ/inst/tests/test-affinity.R
pkg/CHNOSZ/inst/tests/test-basis.R
pkg/CHNOSZ/inst/tests/test-diagram.R
pkg/CHNOSZ/inst/tests/test-findit.R
pkg/CHNOSZ/inst/tests/test-revisit.R
pkg/CHNOSZ/man/affinity.Rd
pkg/CHNOSZ/man/diagram.Rd
pkg/CHNOSZ/man/findit.Rd
pkg/CHNOSZ/man/protein.Rd
pkg/CHNOSZ/man/protein.info.Rd
pkg/CHNOSZ/man/read.expr.Rd
pkg/CHNOSZ/man/revisit.Rd
pkg/CHNOSZ/man/sideeffects.Rd
pkg/CHNOSZ/man/subcrt.Rd
pkg/CHNOSZ/man/util.list.Rd
pkg/CHNOSZ/vignettes/anintro.Rnw
pkg/CHNOSZ/vignettes/anintro.lyx
pkg/CHNOSZ/vignettes/wjd.Rnw
pkg/CHNOSZ/vignettes/wjd.lyx
Log:
equilibrate() is liberated from diagram()
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2012-09-26 02:05:06 UTC (rev 24)
+++ pkg/CHNOSZ/DESCRIPTION 2012-10-02 14:15:53 UTC (rev 25)
@@ -1,6 +1,6 @@
-Date: 2012-09-25
+Date: 2012-10-02
Package: CHNOSZ
-Version: 0.9-7.98
+Version: 0.9.8
Title: Chemical Thermodynamics and Activity Diagrams
Author: Jeffrey M. Dick
Maintainer: Jeffrey M. Dick <jmdick at asu.edu>
Modified: pkg/CHNOSZ/R/affinity.R
===================================================================
--- pkg/CHNOSZ/R/affinity.R 2012-09-26 02:05:06 UTC (rev 24)
+++ pkg/CHNOSZ/R/affinity.R 2012-10-02 14:15:53 UTC (rev 25)
@@ -144,7 +144,7 @@
protein.fun <- function(ip) {
if(ip %% 50 == 0) msgout(paste(ip,"..",sep=""))
tpext <- as.numeric(thermo$protein[iprotein[ip],5:25])
- return(psum(pprod(a[ires],tpext)) - loga.protein[ip])
+ return(Reduce("+", pprod(a[ires],tpext)) - loga.protein[ip])
}
# use another level of indexing to let the function
# report on its progress
@@ -177,29 +177,55 @@
}
# put together return values
+ # constant T and P, in Kelvin and bar
T <- args$T
P <- args$P
- xname <- yname <- ""
- xlim <- ylim <- ""
- xvals <- yvals <- ""
- if(length(args$vars) > 0) {
- xname <- names(args$lims)[1]
- xlim <- args$lims[[1]]
- xvals <- args$vals[[1]]
- if(xname=="T") xvals <- outvert(xvals,"K")
+ # the variable names and values
+ vars <- args$vars
+ vals <- args$vals
+ # if T or P is variable it's not constant;
+ # and set it to user's units
+ iT <- match("T", vars)
+ if(!is.na(iT)) {
+ T <- numeric()
+ vals[[iT]] <- outvert(vals[[iT]], "K")
}
- if(length(args$vars) > 1 & !args$transect) {
- #yname <- args$vars[2]
- yname <- names(args$lims)[2]
- ylim <- args$lims[[2]]
- yvals <- args$vals[[2]]
- if(yname=="T") yvals <- outvert(yvals,"K")
+ iP <- match("P", vars)
+ if(!is.na(iP) > 0) {
+ P <- numeric()
+ vals[[iP]] <- outvert(vals[[iP]], "bar")
}
+ # get Eh
+ args.orig <- list(...)
+ iEh <- match("Eh", names(args.orig))
+ if(!is.na(iEh)) {
+ vars[[iEh]] <- "Eh"
+ # we have to reconstruct the values used because
+ # they got converted to log_a(e-) at an unknown temperature
+ Eharg <- args.orig[[iEh]]
+ if(length(Eharg) > 3) Ehvals <- Eharg
+ else if(length(Eharg) == 3) Ehvals <- seq(Eharg[1], Eharg[2], length.out=Eharg[3])
+ else if(length(Eharg) == 2) Ehvals <- seq(Eharg[1], Eharg[2], length.out=128)
+ vals[[iEh]] <- Ehvals
+ }
+ # get pe and pH and Eh
+ ipe <- match("pe", names(args.orig))
+ if(!is.na(ipe)) {
+ ie <- match("e-", names(args$lims))
+ vars[[ie]] <- "pe"
+ vals[[ie]] <- -args$vals[[ie]]
+ }
+ ipH <- match("pH", names(args.orig))
+ if(!is.na(ipH)) {
+ iH <- match("H+", names(args$lims))
+ vars[[iH]] <- "pH"
+ vals[[iH]] <- -args$vals[[iH]]
+ }
- if(return.buffer) return(c(tb, list(xname=xname, xlim=xlim, xvals=xvals, yname=yname, ylim=ylim)))
- a <- list(sout=sout,property=property,basis=mybasis,species=myspecies,T=T,P=P,
- xname=xname,xlim=xlim,xvals=xvals,yname=yname,ylim=ylim,yvals=yvals,values=a)
- if(buffer) a <- c(a,list(buffer=tb))
+ # content of return value depends on buffer request
+ if(return.buffer) return(c(tb, list(vars=vars, vals=vals)))
+ a <- list(sout=sout, property=property, basis=mybasis, species=myspecies, T=T, P=P, vars=vars, vals=vals, values=a)
+ if(buffer) a <- c(a, list(buffer=tb))
return(a)
}
Modified: pkg/CHNOSZ/R/anim.R
===================================================================
--- pkg/CHNOSZ/R/anim.R 2012-09-26 02:05:06 UTC (rev 24)
+++ pkg/CHNOSZ/R/anim.R 2012-10-02 14:15:53 UTC (rev 25)
@@ -50,7 +50,7 @@
# don't fill fields with heat colors if we're doing a high-T overlay
if(high.T) color <- "lightgrey" else color <- "heat"
# take a single O2-H2O slice along pH
- diagram(slice.affinity(a.cold,3,iframes[i]),color=color,cex=1.5)
+ diagram(slice.affinity(a.cold,3,iframes[i]),fill=color,cex=1.5)
ltext <- substitute(paste(italic(T)==x~degree*C),list(x=25))
lcol <- "black"
# overlay the high-T diagram
@@ -111,17 +111,17 @@
np <- c(rep(head(np,1),8),np,rep(tail(np,1),8))
# now loop until we get to two proteins
for(i in 1:length(np)) {
- d <- diagram(a, ispecies=ispecies, cex=1.5)
+ d <- diagram(a, groups=as.list(ispecies), names=pname[ispecies], normalize=TRUE, cex=1.5)
# note that the darker colors go with higher abundances
# as reported by Anderson and Anderson, 2003
# how many show up on the diagram?
- nshow <- length(unique(as.numeric(d$out)))
+ nshow <- length(unique(as.numeric(d$predominant)))
title(main=paste(np[i]," human plasma proteins (",nshow,
" showing)",sep=""),cex.main=1)
if(c(0,diff(np))[i] != 0) {
# identify the protein with the greatest
# area on the diagram
- imost <- which.max(tabulate(as.numeric(d$out)))
+ imost <- which.max(tabulate(as.numeric(d$predominant)))
# take out that proteins
ispecies <- ispecies[-imost]
}
@@ -195,15 +195,15 @@
a <- affinity(T=T)
} else a <- affinity(T=T,H2=H2)
# calculate activities
- d <- diagram(a,residue=TRUE,legend.x=NULL,ylim=c(-6,-2),plot.it=FALSE)
+ e <- equilibrate(a, normalize=TRUE)
# for each point make a rank plot
- rank <- 1:length(d$logact)
+ rank <- 1:length(e$loga.equil)
for(i in 1:length(ido)) {
# print some progress
if(i%%20 == 0) cat("\n") else cat(".")
# keep track of positions of previous points
loga <- numeric()
- for(j in 1:length(d$logact)) loga <- c(loga,d$logact[[j]][ido[i]])
+ for(j in 1:length(e$loga.equil)) loga <- c(loga, e$loga.equil[[j]][ido[i]])
if(i > 4) myrank4 <- myrank3
if(i > 3) myrank3 <- myrank2
if(i > 2) myrank2 <- myrank1
Modified: pkg/CHNOSZ/R/basis.R
===================================================================
--- pkg/CHNOSZ/R/basis.R 2012-09-26 02:05:06 UTC (rev 24)
+++ pkg/CHNOSZ/R/basis.R 2012-10-02 14:15:53 UTC (rev 25)
@@ -72,10 +72,10 @@
}
thermo$basis$logact[ib] <<- state[i]
} else {
- # look for this state of the species
- ispecies <- suppressMessages(info(species[i], state[i], check.it=FALSE))
+ # look for a species with this name in the requested state
+ ispecies <- suppressMessages(info(thermo$obigt$name[thermo$basis$ispecies[ib]], state[i], check.it=FALSE))
if(is.na(ispecies) | is.list(ispecies))
- stop(paste("state or buffer '",state[i],"' not found for ",species[i],"\n",sep=""))
+ stop(paste("state or buffer '", state[i], "' not found for ", thermo$obigt$name[thermo$basis$ispecies[ib]], "\n", sep=""))
thermo$basis$ispecies[ib] <<- ispecies
thermo$basis$state[ib] <<- state[i]
}
Modified: pkg/CHNOSZ/R/buffer.R
===================================================================
--- pkg/CHNOSZ/R/buffer.R 2012-09-26 02:05:06 UTC (rev 24)
+++ pkg/CHNOSZ/R/buffer.R 2012-10-02 14:15:53 UTC (rev 25)
@@ -20,7 +20,7 @@
if(length(imod)>0) {
if(state[1]=='') {
thermo$buffers <<- thermo$buffers[-imod,]
- cat(paste('mod.buffer: removed ',c2s(species),' in ',
+ msgout(paste('mod.buffer: removed ',c2s(species),' in ',
c2s(unique(name)),' buffer.\n',sep=''))
} else {
if(missing(state)) state <- thermo$buffers$state[imod]
@@ -32,10 +32,10 @@
thermo$buffers$state[imod] <<- state
thermo$buffers$logact[imod] <<- logact
if(identical(state.old,state) & identical(logact.old,logact)) {
- cat(paste('mod.buffer: nothing changed for ',
+ msgout(paste('mod.buffer: nothing changed for ',
c2s(species),' in ',c2s(unique(name)),' buffer.\n',sep=''))
} else {
- cat(paste('mod.buffer: changed state and/or logact of ',
+ msgout(paste('mod.buffer: changed state and/or logact of ',
c2s(species),' in ',c2s(unique(name)),' buffer.\n',sep=''))
}
}
@@ -47,7 +47,7 @@
if(state[1]=='') state <- rep(thermo$opt$state,length.out=ls)
t <- data.frame(name=name,species=species,state=state,logact=logact)
thermo$buffers <<- rbind(thermo$buffers,t)
- cat(paste('mod.buffer: added ',c2s(unique(name)),'.\n',sep=''))
+ msgout(paste('mod.buffer: added ',c2s(unique(name)),'.\n',sep=''))
}
return(invisible(thermo$buffers[thermo$buffers$name %in% name,]))
}
@@ -158,7 +158,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='')
- cat(paste('buffer: ( ',bufname,' ) for activity of ',reqtext,addtext,context,'\n',sep=''))
+ msgout(paste('buffer: ( ',bufname,' ) for activity of ',reqtext,addtext,context,'\n',sep=''))
#print(bufbasis)
# there could still be stuff here (over-defined system?)
xx <- bufbasis[,-ibasis,drop=FALSE]
Modified: pkg/CHNOSZ/R/diagram.R
===================================================================
--- pkg/CHNOSZ/R/diagram.R 2012-09-26 02:05:06 UTC (rev 24)
+++ pkg/CHNOSZ/R/diagram.R 2012-10-02 14:15:53 UTC (rev 25)
@@ -1,648 +1,426 @@
# CHNOSZ/diagram.R
-# plot predominance or activity diagrams
-# from affinities of formation reactions
-# 20061023 jmd
+# plot equilibrium chemical activity and predominance diagrams
+# 20061023 jmd v1
+# 20120927 work with output from either equil() or affinity(),
+# gather plotvals independently of plot parameters (including nd),
+# single return statement
-diagram <- function(affinity,what="logact",ispecies=NULL,balance=NULL,
- names=NA,color=NA,add=FALSE,dotted=0,cex=par('cex'),col=par('col'),
- pe=TRUE,pH=TRUE,ylog=TRUE,main=NULL,cex.names=1,legend.x='topright',
- lty=NULL,col.names=par('fg'),cex.axis=par('cex'),loga.balance=NA,lwd=par('lwd'),
- alpha=FALSE,mar=NULL,residue=FALSE,yline=par('mgp')[1]+1,xrange=NULL,
- ylab=NULL,xlab=NULL,plot.it=TRUE,as.residue=FALSE,mam=TRUE,group=NULL,
- bg=par("bg"),side=1:4,xlim=NULL,ylim=NULL) {
+diagram <- function(
+ # primary input
+ eout,
+ # what to plot
+ what="loga.equil", alpha=FALSE, normalize=FALSE, balance=NULL,
+ groups=as.list(1:length(eout$values)), xrange=NULL,
+ # plot dimensions
+ mar=NULL, yline=par("mgp")[1]+1, side=1:4,
+ # axes
+ ylog=TRUE, xlim=NULL, ylim=NULL, xlab=NULL, ylab=NULL,
+ # sizes
+ cex=par("cex"), cex.names=1, cex.axis=par("cex"),
+ # line styles
+ lty=NULL, lwd=par("lwd"), dotted=0,
+ # colors
+ bg=par("bg"), col=par("col"), col.names=par("col"), fill=NULL,
+ # labels
+ names=NULL, main=NULL, legend.x="topright",
+ # plotting controls
+ add=FALSE, plot.it=TRUE
+) {
- # store input values
- aval <- affinity$values
- asl <- affinity$species$logact
- # number of species possible
- nspecies.orig <- nspecies <- length(aval)
- # number of dimensions (T, P or chemical potentials that are varied)
- mydim <- dim(aval[[1]])
- nd <- length(mydim)
- if(nd==1) if(mydim==1) nd <- 0
- # 'ispecies' which species to consider
- if(is.null(ispecies)) {
- ispecies <- 1:nspecies
- # take out species that have NA affinities
- for(i in 1:nspecies) if(all(is.na(aval[[i]]))) ispecies[i] <- NA
- ispecies <- ispecies[!is.na(ispecies)]
- if(length(ispecies)==0) stop("all species have NA affinities")
- }
- if(!identical(ispecies,1:nspecies)) {
- if(length(ispecies)==0) {
- cn <- caller.name()
- stop(paste("the length of ispecies is zero ( from",cn,")"))
+ ### argument handling ###
+
+ ## check that eout is valid input
+ if(!"sout" %in% names(eout)) stop("'eout' does not look like output from equil() or affinity()")
+
+ ## 'what' can be:
+ # loga.equil - equilibrium activities of species of interest (eout)
+ # basis species - equilibrium activity of a basis species (aout)
+ # missing - property from affinity() or predominances of species (aout)
+ eout.is.aout <- FALSE
+ plot.loga.basis <- FALSE
+ if(missing(what)) {
+ if(!"loga.equil" %in% names(eout)) {
+ eout.is.aout <- TRUE
+ # get the balancing coefficients
+ n.balance <- balance.coeffs(eout, balance)
}
- msgout(paste('diagram: using',length(ispecies),'of',nspecies,'species\n'))
- affinity$species <- affinity$species[ispecies,]
- aval <- aval[ispecies]
+ } else if(what %in% rownames(eout$basis)) {
+ # to calculate the loga of basis species at equilibrium
+ if(!missing(groups)) stop("can't plot equilibrium activities of basis species for grouped species")
+ if(alpha) stop("equilibrium activities of basis species not available with alpha=TRUE")
+ plot.loga.basis <- TRUE
+ } 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")
}
- # number of species that are left
- nspecies <- length(aval)
- # what property we're plotting (jmd 20101014)
- do.loga.equil <- FALSE
- if(affinity$property=="A") {
- if(what %in% rownames(affinity$basis)) {
- # to calculate the loga of basis species at equilibrium
- do.loga.equil <- TRUE
- residue <- FALSE
- balance <- 1
- }
- } else {
- # it might be the values returned by affinity()
- if(missing(what)) what <- affinity$property
- else if(what != affinity$property) stop(paste("can't plot",what,"since the property in 'a' is",affinity$property))
+
+ ## consider a different number of species if we're grouping them together
+ ngroups <- length(groups)
+
+ ## keep the values we plot in one place so they can be modified, plotted and eventually returned
+ # unless something happens below, we'll plot the loga.equil from equilibrate()
+ plotvals <- eout$loga.equil
+ plotvar <- "loga.equil"
+
+ ## deal with output from affinity()
+ if(eout.is.aout) {
+ # plot property from affinity(), divided by balancing coefficients
+ plotvals <- lapply(1:length(eout$values), function(i) {
+ # we divide by the balancing coefficients if we're working with affinities
+ # this is not normalizing the formulas! it's balancing the reactions...
+ # normalizing the formulas is done below
+ eout$values[[i]] / n.balance[i]
+ })
+ plotvar <- eout$property
+ # we change 'A' to 'A/2.303RT' so the axis label is made correctly
+ if(plotvar=="A") plotvar <- "A/2.303RT"
+ msgout(paste("diagram: plotting", plotvar, "from affinity(), divided by balancing coefficients\n"))
}
- # consider a different number of species if we're grouping them together
- if(what=='affinity' & (!mam | nd==0 | nd==1) & !is.null(group)) {
- nspecies.orig <- ngroup <- length(group)
- ispecies <- 1:ngroup
- } else {
- ngroup <- nspecies
+ ## sum activities of species together in groups 20090524
+ # using lapply/Reduce 20120927
+ if(!missing(groups)) {
+ # loop over the groups
+ plotvals <- lapply(groups, function(ispecies) {
+ # remove the logarithms
+ act <- lapply(plotvals[ispecies], function(x) 10^x)
+ # sum the activities
+ return(Reduce("+", act))
+ })
+ # restore the logarithms
+ plotvals <- lapply(plotvals, function(x) log10(x))
+ # we also combine the balancing coefficients for calculations using affinity
+ if(eout.is.aout) n.balance <- sapply(groups, function(ispecies) sum(n.balance[ispecies]))
}
- # handle line type/width/color arguments
- if(plot.it) {
- if(is.null(lty)) lty <- 1:ngroup
- if(is.null(lwd)) lwd <- rep(1,ngroup)
- if(length(lty)!=ngroup) lty <- rep(lty,length.out=ngroup)
- if(length(lwd)!=ngroup) lwd <- rep(lwd,length.out=ngroup)
- if(missing(col.names)) col.names <- col
- if(length(col.names) != ngroup) col.names <- rep(col.names,length.out=ngroup)
- # figure out colors
- if(missing(color)) {
- if(nd==2) {
- if( add | names(dev.cur())=='postscript' ) color <- NULL
- else color <- 'heat'
- } else color <- 'black'
- }
- if(!is.null(color)) {
- color <- rep(color,length.out=ngroup)
- if(is.character(color[1])) {
- # 20091027: take colors as a subset of the original number of species
- if(color[1]=='rainbow') color <- rainbow(nspecies.orig)[ispecies]
- else if(color[1]=='heat') color <- heat.colors(nspecies.orig)[ispecies]
- }
- }
+
+ ## calculate the equilibrium logarithm of activity of a basis species
+ ## (such that affinities of formation reactions are zero)
+ if(plot.loga.basis) {
+ ibasis <- match(what, rownames(eout$basis))
+ # the logarithm of activity used in the affinity calculation
+ is.loga.basis <- can.be.numeric(eout$basis$logact[ibasis])
+ if(!is.loga.basis) stop(paste("the logarithm of activity for basis species", what, "is not numeric - was a buffer selected?"))
+ loga.basis <- as.numeric(eout$basis$logact[ibasis])
+ # the reaction coefficients for this basis species
+ nu.basis <- eout$species[, ibasis]
+ # the logarithm of activity where affinity = 0
+ plotvals <- lapply(1:length(eout$values), function(x) {
+ # eout$values is a strange name for affinity ... should be named something like eout$affinity ...
+ loga.basis - eout$values[[x]]/nu.basis[x]
+ })
+ plotvar <- what
}
- # handle legend argument (turn off legend
- # in 1D diagrams if too many species)
- if(nd==1 & ngroup > 10 & missing(legend.x)) legend.x <- NULL
- # covert from activities of proton or electron to pH or pe
- if(affinity$xname=='H+' & pH) { affinity$xname <- 'pH'; affinity$xlim[1:2] <- -affinity$xlim[1:2]; affinity$xvals <- -affinity$xvals }
- if(affinity$xname=='e-' & pe) { affinity$xname <- 'pe'; affinity$xlim[1:2] <- -affinity$xlim[1:2]; affinity$xvals <- -affinity$xvals }
- if(affinity$yname=='H+' & pH) { affinity$yname <- 'pH'; affinity$ylim[1:2] <- -affinity$ylim[1:2]; affinity$yvals <- -affinity$yvals }
- if(affinity$yname=='e-' & pe) { affinity$yname <- 'pe'; affinity$ylim[1:2] <- -affinity$ylim[1:2]; affinity$yvals <- -affinity$yvals }
- # T/P conversions
- if(T.units()=='C') {
- if(affinity$xname=='T') { affinity$xlim[1:2] <- convert(affinity$xlim[1:2],T.units()) }
- if(affinity$yname=='T') { affinity$ylim[1:2] <- convert(affinity$ylim[1:2],T.units()) }
+
+ ## alpha: plot fractional degree of formation instead of logarithms of activities
+ ## scale the activities to sum=1 ... 20091017
+ if(alpha) {
+ # remove the logarithms
+ act <- lapply(eout$loga.equil, function(x) 10^x)
+ # sum the activities
+ sumact <- Reduce("+", act)
+ # divide activities by the total
+ alpha <- lapply(act, function(x) x/sumact)
+ plotvals <- alpha
+ plotvar <- "alpha"
}
- if(P.units()=='MPa') {
- if(affinity$xname=='P') { affinity$xlim[1:2] <- convert(affinity$xlim[1:2],P.units()) }
- if(affinity$yname=='P') { affinity$ylim[1:2] <- convert(affinity$ylim[1:2],P.units()) }
- }
- # 'balance' how to balance reactions
- # (name of basis species, 'PBB', 1, or NULL for auto-select)
- # (or vector to set nbalance)
- # 'nbalance' vector of coefficients of the balanced component
- isprotein <- grep('_',as.character(affinity$species$name))
- if(is.null(balance) & length(isprotein)==nspecies)
- balance <- 'PBB'
- # 20091119 make residue=TRUE the default for protein reactions
- if(missing(residue) & length(isprotein)==nspecies)
- residue <- TRUE
- if(is.null(balance)) {
- ib <- which.balance(affinity$species)
- if(length(ib)==0) stop("a shared basis species is not present in all formation reactions")
- if(!is.null(ib)) {
- balance <- rownames(basis())[ib[1]]
- nbalance <- affinity$species[,ib[1]]
- msgout(paste('diagram: balanced quantity is moles of',balance,'\n'))
- } else {
- balance <- 1
- nbalance <- rep(1,nspecies)
- msgout('diagram: balanced quantity is moles of chemical formulas as written\n')
+
+ ## identify predominant species
+ predominant <- NA
+ if(plotvar %in% c("loga.equil", "alpha", "A/2.303RT")) {
+ pv <- plotvals
+ for(i in 1:length(pv)) {
+ # change any NAs in the plotvals to -Inf, so that
+ # they don't get on the plot, but permit others to
+ pv[[i]][is.na(pv[[i]])] <- -Inf
+ # 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(is.character(balance[1])) {
- # is the balance is "PBB" for protein backbone
- if(identical(balance,'PBB')) {
- if(length(isprotein) != nspecies)
- stop('diagram: PBB was the requested balance, but ',
- affinity$species$name[-isprotein],' is not protein.')
- nbalance <- protein.length(affinity$species$name)
- msgout('diagram: balanced quantity is moles of protein backbone group\n')
- } else {
- # is the balance the name of a basis species
- ib <- match(balance,rownames(basis()))
- if(!is.na(ib)) {
- nbalance <- affinity$species[,ib]
- if(TRUE %in% (nbalance==0)) {
- inotbalance <- which(nbalance==0)
- # take out species that can't be balanced
- if(length(inotbalance)==1) hh <- "has" else hh <- "have"
- msgout(paste('diagram: removing ',c2s(affinity$species$name[inotbalance]),
- ', which ',hh,' no ',balance,'\n',sep=''))
- aval <- aval[-inotbalance]
- affinity$species <- affinity$species[-inotbalance,]
- nbalance <- nbalance[-inotbalance]
- nspecies <- nspecies - length(inotbalance)
- }
- } else {
- stop('requested balance (',balance,') is not a the name of a basis species or other value I understand.')
- }
- }
- } else if(is.numeric(balance[1])) {
- # user-defined balance
- nbalance <- rep(balance,length.out=nspecies)
- if(all(nbalance==1)) balance <- "species" else balance <- "user-defined"
- }
+ predominant <- which.pmax(pv)
+ dim(predominant) <- dim(pv[[1]])
}
- if(length(nbalance) < 100) msgout(paste('diagram: balancing coefficients are ',c2s(nbalance),'\n',sep=''))
- # rewrite the balance vector for residue reactions
- if(residue) {
- if(any(nbalance < 0)) stop("negative balance prohibits using residue equivalents")
- # affinities divided by balance
- for(i in 1:nspecies) aval[[i]] <- aval[[i]]/nbalance[i]
- oldlogact <- affinity$species$logact
- affinity$species$logact <- affinity$species$logact + log10(nbalance)
- oldbalance <- nbalance
- nbalance <- rep(1,length(nbalance))
- msgout(paste('diagram: using residue equivalents\n'))
- }
- # get labels for the plot
- if(!is.null(names)) {
- if(is.na(names[1])) {
- # if names are missing ... make them up
+
+ ### now on to the plotting ###
+
+ if(plot.it) {
+
+ ### 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)
+ lwd <- rep(lwd, length.out=ngroups)
+ col <- rep(col, length.out=ngroups)
+ col.names <- rep(col.names, length.out=ngroups)
+
+ ## make up some names for lines/fields if they are missing
+ if(missing(names)) {
# properties of basis species or reactions?
- if(affinity$property %in% c('G.basis','logact.basis')) names <- rownames(affinity$basis)
+ if(eout$property %in% c("G.basis", "logact.basis")) names <- rownames(eout$basis)
else {
- if(!is.null(group)) {
- if(is.null(names(group))) names(group) <- paste("group",1:length(group),sep="")
- names <- names(group)
+ if(!missing(groups)) {
+ if(is.null(names(groups))) names <- paste("group", 1:length(groups), sep="")
+ else names <- names(groups)
}
- else names <- as.character(affinity$species$name)
- # remove common organism label for proteins
- if(all(grep("_",names)>0)) {
- pname <- oname <- character()
- for(i in 1:length(names)) {
- mynames <- s2c(names[i],sep="_",keep.sep=FALSE)
- pname <- c(pname,mynames[1])
- oname <- c(oname,mynames[2])
- }
- if(length(unique(oname))==1) for(j in 1:length(names)) names[j] <- pname[j]
+ else names <- as.character(eout$species$name)
+ # remove non-unique organism or protein names
+ if(all(grepl("_", names))) {
+ # everything before the underscore (the protein)
+ pname <- gsub("_.*$", "", names)
+ # everything after the underscore (the organism)
+ oname <- gsub("^.*_", "", names)
+ # if the pname or oname are all the same, use the other one as identifying name
+ if(length(unique(pname))==1) names <- oname
+ if(length(unique(oname))==1) names <- pname
}
- # append species label to distinguish ambiguous names
- #if(length(unique(names))!=length(names))
- # names <- paste(names,' (',affinity$species$state,')',sep='')
+ # append state to distinguish ambiguous species names
isdup <- names %in% names[duplicated(names)]
if(any(isdup)) names[isdup] <- paste(names[isdup],
- " (",affinity$species$state[isdup],")",sep="")
+ " (", eout$species$state[isdup], ")", sep="")
}
}
- }
- ### end variable assignment
- ### various plotting functions
- ### color plot function
- plot.color <- function(xs,ys,out,color,nspecies) {
- # handle min/max reversal
- if(xs[1] > xs[length(xs)]) {
- tc <- out; t <- numeric()
- for(i in 1:length(xs)) {
- t[i] <- xs[length(xs)+1-i]
- tc[,i] <- out[,length(xs)+1-i]
- }
- out <- tc; xs <- t
- }
- if(ys[1] > ys[length(ys)]) {
- tc <- out; t <- numeric()
- for(i in 1:length(ys)) {
- t[i] <- ys[length(ys)+1-i]
- tc[i,] <- out[length(ys)+1-i,]
- }
- out <- tc; ys <- t
- }
- # the z values
- zs <- out
- for(i in 1:nrow(zs)) zs[i,] <- out[nrow(zs)+1-i,]
- zs <- t(zs)
- breaks <- c(0,1:nspecies) + 0.5
- image(x=xs,y=ys,z=zs,col=color,add=TRUE,breaks=breaks)
- }
- ### curve plot function
- # 20091116 replaced plot.curve with plot.line; different
- # name, same functionality, *much* faster
- plot.line <- function(out,xlim,ylim,dotted,col,lwd,xrange) {
- # plot boundary lines between predominance fields
- vline <- function(out,ix) {
- ny <- nrow(out)
- xs <- rep(ix,ny*2+1)
- ys <- c(rep(ny:1,each=2),0)
- y1 <- out[,ix]
- y2 <- out[,ix+1]
- # no line segment inside a stability field
- iy <- which(y1==y2)
- ys[iy*2] <- NA
- # no line segment at a dotted position
- iyd <- which(ys%%dotted==0)
- ys[iyd] <- NA
- return(list(xs=xs,ys=ys))
- }
- hline <- function(out,iy) {
- nx <- nrow(out)
- ny <- ncol(out)
- ys <- rep(ny-iy,nx*2+1)
- xs <- c(0,rep(1:nx,each=2))
- x1 <- out[iy,]
- x2 <- out[iy+1,]
- # no line segment inside a stability field
- ix <- which(x1==x2)
- xs[ix*2] <- NA
- # no line segment at a dotted position
- ixd <- which(xs%%dotted==0)
- xs[ixd] <- NA
- return(list(xs=xs,ys=ys))
- }
- clipfun <- function(z,zlim) {
- if(zlim[2] > zlim[1]) {
- z[z>zlim[2]] <- zlim[2]
- z[z<zlim[1]] <- zlim[1]
- } else {
- z[z>zlim[1]] <- zlim[1]
- z[z<zlim[2]] <- zlim[2]
- }
- return(z)
- }
- rx <- (xlim[2] - xlim[1]) / (ncol(out) - 1)
- ry <- (ylim[2] - ylim[1]) / (nrow(out) - 1)
- # vertical lines
- xs <- ys <- NA
- for(ix in 1:(ncol(out)-1)) {
- vl <- vline(out,ix)
- xs <- c(xs,vl$xs,NA)
- ys <- c(ys,vl$ys,NA)
- }
- xs <- xlim[1] + (xs - 0.5) * rx
- ys <- ylim[1] + (ys - 0.5) * ry
- ys <- clipfun(ys,ylim)
- lines(xs,ys,col=col,lwd=lwd)
- # horizontal lines
- xs <- ys <-NA
- for(iy in 1:(nrow(out)-1)) {
- hl <- hline(out,iy)
- xs <- c(xs,hl$xs,NA)
- ys <- c(ys,hl$ys,NA)
- }
- xs <- xlim[1] + (xs - 0.5) * rx
- ys <- ylim[1] + (ys - 0.5) * ry
- xs <- clipfun(xs,xlim)
- if(!is.null(xrange)) xs <- clipfun(xs,xrange)
- lines(xs,ys,col=col,lwd=lwd)
- }
- ### label plot function
- # calculate coordinates for field labels
- plot.names <- function(out,xs,ys,names) {
- ll <- ngroup
- lx <- numeric(ll); ly <- numeric(ll); n <- numeric(ll)
- for(j in nrow(out):1) {
- # 20091116 for speed, loop over ngroup instead of k (columns)
- for(i in 1:ngroup) {
- k <- which(out[j,]==i)
- if(length(k)==0) next
- lx[i] <- lx[i] + sum(xs[k])
- ly[i] <- ly[i] + length(k)*ys[nrow(out)+1-j]
- n[i] <- n[i] + length(k)
- }
- }
- lx <- lx[n!=0]
- ly <- ly[n!=0]
- is <- n!=0
- n <- n[n!=0]
- lx <- lx/n
- ly <- ly/n
- # plot field labels
- # the cex argument in this function specifies the character
- # expansion of the labels relative to the current
- text(lx,ly,labels=names[is],cex=cex.names,col=col.names[is])
- }
- ### property description
- axis.title <- function(what,suffix="") {
- if(what=="A") return(as.expression(substitute(italic(bold("A"))/2.303*italic(R)*italic(T)~~x,list(x=suffix))))
- else if(what=="logact") return(as.expression(substitute(log*italic(a)~~x,list(x=suffix))))
- else if(what=='logact.basis') return(paste('logQ*',suffix))
- else if(!what %in% c("logK","logQ")) return(paste("Delta",what,suffix))
- else return(what)
- }
- balance.title <- function(balance) {
- if(balance==1) return('per mole of product')
- else if(is.numeric(balance)) return(paste('per',balance,'moles of product'))
- else return(paste('per mole of',balance))
- }
- ### done with function definitions
+ if(nd==0) {
- ### now on to the calculation and plots
- # do speciation calculations
- # unless using mam (maximum affinity method) for 2-D diagram
- if(what=="logact" & (!mam | nd==0 | nd==1) ) {
- # compute the activities of species
- # logarithm of total activity of the balanced component
- ib <- match(balance,colnames(affinity$species))
- if(!is.numeric(loga.balance)) {
- thisa <- sum(10^affinity$species$logact * nbalance)
- if(thisa < 0) thisa <- -thisa
- loga.balance <- log10(thisa)
- if(can.be.numeric(balance)) msgout('diagram: log total activity of species is ',loga.balance,'\n',sep='')
- else msgout('diagram: log total activity of ',balance,' (from species) is ',loga.balance,'\n',sep='')
- } else {
- msgout('diagram: log total activity of ',balance,' (from argument) is ',loga.balance,'\n',sep='')
- }
- # Astar: the affinities/2.303RT of formation of species
- # in their standard-state activities
- Astar <- aval
- if(residue) for(j in 1:length(Astar)) Astar[[j]] <- Astar[[j]] + oldlogact[j]/oldbalance[j]
- else for(j in 1:length(Astar)) Astar[[j]] <- Astar[[j]] + affinity$species$logact[j]
- # now get the equilibrium activities
- if(all(nbalance==1)) A <- equil.boltzmann(Astar, nbalance, loga.balance)
- else A <- equil.reaction(Astar, nbalance, loga.balance)
- # if we rewrote the formation reactions per residue,
- # get back to activities of species
- if(residue & !as.residue) for(j in 1:length(A)) A[[j]] <- A[[j]] - log10(oldbalance[j])
- # groups the species together into meta-species 20090524
- if(!is.null(group)) {
- # store our dimensions
- mydim <- dim(A[[1]])
- # set up the output
- B <- A[1:length(group)]
- for(i in 1:length(B)) {
- B[[i]] <- as.numeric(A[[i]])
- B[[i]][] <- 0
- }
- # add up the activities
- for(i in 1:length(group)) {
- for(j in 1:length(group[[i]])) {
- B[[i]] <- B[[i]] + 10^as.numeric(A[[group[[i]][j]]])
- }
- }
- # return to logarithms and replace dimensions
- for(i in 1:length(B)) {
- # do we want to divide by the number of representative for
- # each group? probably not 20091017
- #B[[i]] <- log10(B[[i]]/length(group[[i]]))
- B[[i]] <- log10(B[[i]])
- dim(B[[i]]) <- mydim
- }
- A <- B
- }
- } else if(do.loga.equil) {
- # calculate the logarithm of activity of a basis species
- # at equilibrium
- ibasis <- match(what,rownames(affinity$basis))
- # the reference logact
- AV <- aval
- loga.basis <- as.numeric(affinity$basis$logact[ibasis])
- if(is.na(loga.basis)) stop(paste("the logact for basis species",what,"is not numeric"))
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/chnosz -r 25
More information about the CHNOSZ-commits
mailing list