[CHNOSZ-commits] r176 - in pkg/CHNOSZ: . R data inst man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sat Feb 25 03:08:43 CET 2017
Author: jedick
Date: 2017-02-25 03:08:42 +0100 (Sat, 25 Feb 2017)
New Revision: 176
Removed:
pkg/CHNOSZ/man/util.affinity.Rd
pkg/CHNOSZ/man/util.args.Rd
pkg/CHNOSZ/man/util.character.Rd
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/NAMESPACE
pkg/CHNOSZ/R/IAPWS95.R
pkg/CHNOSZ/R/basis.R
pkg/CHNOSZ/R/buffer.R
pkg/CHNOSZ/R/diagram.R
pkg/CHNOSZ/R/equilibrate.R
pkg/CHNOSZ/R/examples.R
pkg/CHNOSZ/R/hkf.R
pkg/CHNOSZ/R/info.R
pkg/CHNOSZ/R/makeup.R
pkg/CHNOSZ/R/objective.R
pkg/CHNOSZ/R/revisit.R
pkg/CHNOSZ/R/species.R
pkg/CHNOSZ/R/util.affinity.R
pkg/CHNOSZ/R/util.args.R
pkg/CHNOSZ/R/util.character.R
pkg/CHNOSZ/R/util.data.R
pkg/CHNOSZ/R/util.fasta.R
pkg/CHNOSZ/R/util.formula.R
pkg/CHNOSZ/R/util.list.R
pkg/CHNOSZ/R/util.misc.R
pkg/CHNOSZ/R/util.plot.R
pkg/CHNOSZ/R/util.program.R
pkg/CHNOSZ/R/util.protein.R
pkg/CHNOSZ/R/util.units.R
pkg/CHNOSZ/data/thermo.R
pkg/CHNOSZ/inst/NEWS
pkg/CHNOSZ/man/EOSregress.Rd
pkg/CHNOSZ/man/IAPWS95.Rd
pkg/CHNOSZ/man/affinity.Rd
pkg/CHNOSZ/man/basis.Rd
pkg/CHNOSZ/man/buffer.Rd
pkg/CHNOSZ/man/diagram.Rd
pkg/CHNOSZ/man/eos.Rd
pkg/CHNOSZ/man/equilibrate.Rd
pkg/CHNOSZ/man/findit.Rd
pkg/CHNOSZ/man/info.Rd
pkg/CHNOSZ/man/ionize.aa.Rd
pkg/CHNOSZ/man/makeup.Rd
pkg/CHNOSZ/man/nonideal.Rd
pkg/CHNOSZ/man/objective.Rd
pkg/CHNOSZ/man/revisit.Rd
pkg/CHNOSZ/man/species.Rd
pkg/CHNOSZ/man/swap.basis.Rd
pkg/CHNOSZ/man/util.blast.Rd
pkg/CHNOSZ/man/util.data.Rd
pkg/CHNOSZ/man/util.fasta.Rd
pkg/CHNOSZ/man/util.formula.Rd
pkg/CHNOSZ/man/util.list.Rd
pkg/CHNOSZ/man/util.misc.Rd
pkg/CHNOSZ/man/util.plot.Rd
pkg/CHNOSZ/man/util.program.Rd
pkg/CHNOSZ/man/util.units.Rd
pkg/CHNOSZ/man/water.Rd
pkg/CHNOSZ/man/wjd.Rd
Log:
undocument unexported functions
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2017-02-24 13:00:08 UTC (rev 175)
+++ pkg/CHNOSZ/DESCRIPTION 2017-02-25 02:08:42 UTC (rev 176)
@@ -1,6 +1,6 @@
-Date: 2017-02-24
+Date: 2017-02-25
Package: CHNOSZ
-Version: 1.0.8-65
+Version: 1.0.8-66
Title: Chemical Thermodynamics and Activity Diagrams
Author: Jeffrey Dick
Maintainer: Jeffrey Dick <j3ffdick at gmail.com>
Modified: pkg/CHNOSZ/NAMESPACE
===================================================================
--- pkg/CHNOSZ/NAMESPACE 2017-02-24 13:00:08 UTC (rev 175)
+++ pkg/CHNOSZ/NAMESPACE 2017-02-25 02:08:42 UTC (rev 176)
@@ -22,14 +22,14 @@
"def2gi", "read.blast", "id.blast",
"add.obigt", "RH2obigt",
"expr.property", "expr.units",
- "mass", "entropy", "get.formula", "GHS", "water",
+ "mass", "entropy", "GHS", "water",
"i2A", "invertible.combs",
"dPdTtr", "Ttr",
"count.aa", "nucleic.complement", "nucleic.formula",
"rho.IAPWS95", "IAPWS95", "water.AW90", "WP02.auxiliary", "water.IAPWS95",
"getrank", "parent", "sciname", "allparents", "getnodes", "getnames",
"protein.obigt", "hkf", "cgl", "which.pmax",
- "equil.boltzmann", "equil.reaction", "find.TP",
+ "equil.boltzmann", "equil.reaction", "find.tp",
"ionize.aa", "MP90.cp", "aasum", "read.expr",
"anim.carboxylase",
"qqr", "RMSD", "CVRMSD", "spearman", "DGmix", "DDGmix", "DGtr",
@@ -44,9 +44,18 @@
# equilibrium vignette
"usrfig",
# hotspring vignette
- "strip"
-# no other functions used in tests
-# co-documented with above
+ "strip",
+# ecipex package
+ "count.elements",
+# (no other functions are used in the tests)
+# other exported functions that are not used above
+ "write.blast", "checkEOS", "checkGHS", "check.obigt",
+ "V_s_var", "Cp_s_var",
+ "DGinf", "SD", "pearson", "shannon", "CV", "logact",
+ "feldspar", "apc", "EOSlab", "EOScalc",
+ "basis.elements", "element.mu", "ibasis",
+ "water.SUPCRT92", "water.props", "eqdata", "plot_findit",
+ "nonideal", "anim.TCA", "uniprot.aa", "run.guess"
)
# Load shared objects
Modified: pkg/CHNOSZ/R/IAPWS95.R
===================================================================
--- pkg/CHNOSZ/R/IAPWS95.R 2017-02-24 13:00:08 UTC (rev 175)
+++ pkg/CHNOSZ/R/IAPWS95.R 2017-02-25 02:08:42 UTC (rev 176)
@@ -1,6 +1,116 @@
# functions for properties of water using
# the IAPWS-95 formulation (Wagner and Pruss, 2002)
+IAPWS95 <- function(property,T=298.15,rho=1000) {
+ ## the IAPWS-95 formulation for ordinary water substance
+ ## Wagner and Pruss, 2002
+ property <- tolower(property)
+ # triple point
+ T.triple <- 273.16 # K
+ P.triple <- 611.657 # Pa
+ rho.triple.liquid <- 999.793
+ rho.triple.vapor <- 0.00485458
+ # normal boiling point
+ T.boiling <- 373.124
+ P.boiling <- 0.101325
+ rho.boiling.liquid <- 958.367
+ rho.boiling.vapor <- 0.597657
+ # critical point constants
+ T.critical <- 647.096 # K
+ rho.critical <- 322 # kg m-3
+ # specific and molar gas constants
+ R <- 0.46151805 # kJ kg-1 K-1
+ # R.M <- 8.314472 # J mol-1 K-1
+ # molar mass
+ M <- 18.015268 # g mol-1
+ ## define functions idealgas and residual, supplying arguments delta and tau
+ idealgas <- function(p) IAPWS95.idealgas(p, delta, tau)
+ residual <- function(p) IAPWS95.residual(p, delta, tau)
+ ## relation of thermodynamic properties to Helmholtz free energy
+ a <- function() {
+ x <- idealgas('phi')+residual('phi')
+ return(x*R*T)
+ }
+ # Table 6.3
+ p <- function() {
+ x <- 1 + delta*residual('phi.delta')
+ return(x*rho*R*T/1000) # for MPa
+ }
+ s <- function() {
+ x <- tau * (idealgas('phi.tau')+residual('phi.tau'))-idealgas('phi')-residual('phi')
+ return(x*R)
+ }
+ u <- function() {
+ x <- tau * (idealgas('phi.tau')+residual('phi.tau'))
+ return(x*R*T)
+ }
+ h <- function() {
+ x <- 1 + tau * (idealgas('phi.tau')+residual('phi.tau')) + delta*residual('phi.delta')
+ return(x*R*T)
+ }
+ g <- function() {
+ x <- 1 + idealgas('phi') + residual('phi') + delta*residual('phi.delta')
+ return(x*R*T)
+ }
+ cv <- function() {
+ x <- -tau^2*(idealgas('phi.tau.tau')+residual('phi.tau.tau'))
+ return(x*R)
+ }
+ cp <- function() {
+ x <- -tau^2*(idealgas('phi.tau.tau')+residual('phi.tau.tau')) +
+ (1+delta*residual('phi.delta')-delta*tau*residual('phi.delta.tau'))^2 /
+ (1+2*delta*residual('phi.delta')+delta^2*residual('phi.delta.delta'))
+ return(x*R)
+ }
+# 20090420 speed of sound calculation is incomplete
+# (delta.liquid and drhos.dT not visible)
+# cs <- function() {
+# x <- -tau^2*(idealgas('phi.tau.tau')+residual('phi.tau.tau')) +
+# (1+delta*residual('phi.delta')-delta*tau*residual('phi.delta.tau'))^2 /
+# (1+2*delta*residual('phi.delta')+delta^2*residual('phi.delta.delta')) *
+# ((1+delta.liquid*residual('phi.delta')-delta.liquid*tau*residual('phi.tau.tau'))-rho.critical/(R*delta.liquid)*drhos.dT)
+# return(x*R)
+# }
+ w <- function() {
+ x <- 1 + 2*delta*residual('phi.delta') + delta^2*residual('phi.delta.delta') -
+ (1+delta*residual('phi.delta')-delta*tau*residual('phi.delta.tau'))^2 /
+ tau^2*(idealgas('phi.tau.tau')+residual('phi.tau.tau'))
+ return(sqrt(x*R*T))
+ }
+ mu <- function() {
+ x <- -(delta*residual('phi.delta')+delta^2*residual('phi.delta.delta')+delta*tau*residual('phi.delta.tau')) /
+ ( ( 1+delta*residual('phi.delta')-delta*tau*residual('phi.delta.tau')^2 ) - tau^2 *
+ (idealgas('phi.tau.tau')+residual('phi.tau.tau'))*(1+2*delta*residual('phi.delta')+delta^2*residual('phi.delta.delta')) )
+ return(x/(R*rho))
+ }
+ ## run the calculations
+ ww <- NULL
+ my.T <- T
+ my.rho <- rho
+ for(j in 1:length(property)) {
+ t <- numeric()
+ for(i in 1:length(my.T)) {
+ T <- my.T[i]
+ rho <- my.rho[i]
+ # Equation 6.4
+ delta <- rho / rho.critical
+ tau <- T.critical / T
+ t <- c(t,get(property[j])())
+ }
+ t <- data.frame(t)
+ if(j==1) ww <- t else ww <- cbind(ww,t)
+ }
+ colnames(ww) <- property
+ return(ww)
+}
+
+### unexported functions ###
+
+# IAPWS95.idealgas and IAPWS95.residual are supporting functions to IAPWS95 for calculating
+# the ideal-gas and residual parts in the IAPWS-95 formulation.
+# The value of p can be one of phi, phi.delta, phi.delta.delta, phi.tau, phi.tau.tau, or phi.delta.tau,
+# to calculate the specific dimensionless Helmholtz free energy (phi) or one of its derivatives.
+
IAPWS95.idealgas <- function(p, delta, tau) {
## the ideal gas part in the IAPWS-95 formulation
# from Table 6.1 of Wagner and Pruss, 2002
@@ -142,107 +252,3 @@
}
return(get(p)())
}
-
-IAPWS95 <- function(property,T=298.15,rho=1000) {
- ## the IAPWS-95 formulation for ordinary water substance
- ## Wagner and Pruss, 2002
- property <- tolower(property)
- # triple point
- T.triple <- 273.16 # K
- P.triple <- 611.657 # Pa
- rho.triple.liquid <- 999.793
- rho.triple.vapor <- 0.00485458
- # normal boiling point
- T.boiling <- 373.124
- P.boiling <- 0.101325
- rho.boiling.liquid <- 958.367
- rho.boiling.vapor <- 0.597657
- # critical point constants
- T.critical <- 647.096 # K
- rho.critical <- 322 # kg m-3
- # specific and molar gas constants
- R <- 0.46151805 # kJ kg-1 K-1
- # R.M <- 8.314472 # J mol-1 K-1
- # molar mass
- M <- 18.015268 # g mol-1
- ## define functions idealgas and residual, supplying arguments delta and tau
- idealgas <- function(p) IAPWS95.idealgas(p, delta, tau)
- residual <- function(p) IAPWS95.residual(p, delta, tau)
- ## relation of thermodynamic properties to Helmholtz free energy
- a <- function() {
- x <- idealgas('phi')+residual('phi')
- return(x*R*T)
- }
- # Table 6.3
- p <- function() {
- x <- 1 + delta*residual('phi.delta')
- return(x*rho*R*T/1000) # for MPa
- }
- s <- function() {
- x <- tau * (idealgas('phi.tau')+residual('phi.tau'))-idealgas('phi')-residual('phi')
- return(x*R)
- }
- u <- function() {
- x <- tau * (idealgas('phi.tau')+residual('phi.tau'))
- return(x*R*T)
- }
- h <- function() {
- x <- 1 + tau * (idealgas('phi.tau')+residual('phi.tau')) + delta*residual('phi.delta')
- return(x*R*T)
- }
- g <- function() {
- x <- 1 + idealgas('phi') + residual('phi') + delta*residual('phi.delta')
- return(x*R*T)
- }
- cv <- function() {
- x <- -tau^2*(idealgas('phi.tau.tau')+residual('phi.tau.tau'))
- return(x*R)
- }
- cp <- function() {
- x <- -tau^2*(idealgas('phi.tau.tau')+residual('phi.tau.tau')) +
- (1+delta*residual('phi.delta')-delta*tau*residual('phi.delta.tau'))^2 /
- (1+2*delta*residual('phi.delta')+delta^2*residual('phi.delta.delta'))
- return(x*R)
- }
-# 20090420 speed of sound calculation is incomplete
-# (delta.liquid and drhos.dT not visible)
-# cs <- function() {
-# x <- -tau^2*(idealgas('phi.tau.tau')+residual('phi.tau.tau')) +
-# (1+delta*residual('phi.delta')-delta*tau*residual('phi.delta.tau'))^2 /
-# (1+2*delta*residual('phi.delta')+delta^2*residual('phi.delta.delta')) *
-# ((1+delta.liquid*residual('phi.delta')-delta.liquid*tau*residual('phi.tau.tau'))-rho.critical/(R*delta.liquid)*drhos.dT)
-# return(x*R)
-# }
- w <- function() {
- x <- 1 + 2*delta*residual('phi.delta') + delta^2*residual('phi.delta.delta') -
- (1+delta*residual('phi.delta')-delta*tau*residual('phi.delta.tau'))^2 /
- tau^2*(idealgas('phi.tau.tau')+residual('phi.tau.tau'))
- return(sqrt(x*R*T))
- }
- mu <- function() {
- x <- -(delta*residual('phi.delta')+delta^2*residual('phi.delta.delta')+delta*tau*residual('phi.delta.tau')) /
- ( ( 1+delta*residual('phi.delta')-delta*tau*residual('phi.delta.tau')^2 ) - tau^2 *
- (idealgas('phi.tau.tau')+residual('phi.tau.tau'))*(1+2*delta*residual('phi.delta')+delta^2*residual('phi.delta.delta')) )
- return(x/(R*rho))
- }
- ## run the calculations
- ww <- NULL
- my.T <- T
- my.rho <- rho
- for(j in 1:length(property)) {
- t <- numeric()
- for(i in 1:length(my.T)) {
- T <- my.T[i]
- rho <- my.rho[i]
- # Equation 6.4
- delta <- rho / rho.critical
- tau <- T.critical / T
- t <- c(t,get(property[j])())
- }
- t <- data.frame(t)
- if(j==1) ww <- t else ww <- cbind(ww,t)
- }
- colnames(ww) <- property
- return(ww)
-}
-
Modified: pkg/CHNOSZ/R/basis.R
===================================================================
--- pkg/CHNOSZ/R/basis.R 2017-02-24 13:00:08 UTC (rev 175)
+++ pkg/CHNOSZ/R/basis.R 2017-02-25 02:08:42 UTC (rev 176)
@@ -85,7 +85,7 @@
return(put.basis(ispecies, logact))
}
-## non-exported functions
+### unexported functions ###
# to add the basis to thermo$obigt
put.basis <- function(ispecies, logact = rep(NA, length(ispecies))) {
@@ -212,3 +212,4 @@
logact[is.na(logact)] <- -3
return(logact)
}
+
Modified: pkg/CHNOSZ/R/buffer.R
===================================================================
--- pkg/CHNOSZ/R/buffer.R 2017-02-24 13:00:08 UTC (rev 175)
+++ pkg/CHNOSZ/R/buffer.R 2017-02-25 02:08:42 UTC (rev 176)
@@ -56,6 +56,8 @@
return(invisible(thermo$buffers[thermo$buffers$name %in% name,]))
}
+### unexported functions ###
+
buffer <- function(logK=NULL,ibasis=NULL,logact.basis=NULL,is.buffer=NULL,balance='PBB') {
thermo <- get("thermo")
# if logK is NULL load the buffer species
Modified: pkg/CHNOSZ/R/diagram.R
===================================================================
--- pkg/CHNOSZ/R/diagram.R 2017-02-24 13:00:08 UTC (rev 175)
+++ pkg/CHNOSZ/R/diagram.R 2017-02-25 02:08:42 UTC (rev 176)
@@ -565,7 +565,7 @@
}
}
-find.TP <- function(x) {
+find.tp <- function(x) {
# find triple points in an matrix of integers 20120525 jmd
# these are the locations closest to the greatest number of different values
# rearrange the matrix in the same way that diagram() does for 2-D predominance diagrams
Modified: pkg/CHNOSZ/R/equilibrate.R
===================================================================
--- pkg/CHNOSZ/R/equilibrate.R 2017-02-24 13:00:08 UTC (rev 175)
+++ pkg/CHNOSZ/R/equilibrate.R 2017-02-25 02:08:42 UTC (rev 176)
@@ -228,6 +228,9 @@
return(logact)
}
+### unexported functions ###
+
+# return a list containing the balancing coefficients (n) and a textual description (description)
balance <- function(aout, balance=NULL) {
## generate n.balance from user-given or automatically identified basis species
## extracted from equilibrate() 20120929
Modified: pkg/CHNOSZ/R/examples.R
===================================================================
--- pkg/CHNOSZ/R/examples.R 2017-02-24 13:00:08 UTC (rev 175)
+++ pkg/CHNOSZ/R/examples.R 2017-02-25 02:08:42 UTC (rev 176)
@@ -6,10 +6,10 @@
# run all the examples in CHNOSZ documentation
.ptime <- proc.time()
topics <- c("thermo", "sideeffects", "examples",
- "util.args", "util.array", "util.blast", "util.character", "util.data", "util.expression",
- "util.fasta", "util.formula", "util.matrix", "util.misc", "util.program", "util.seq", "util.units",
+ "util.array", "util.blast", "util.data", "util.expression",
+ "util.fasta", "util.formula", "util.matrix", "util.misc", "util.seq", "util.units",
"util.water", "taxonomy", "info", "protein.info", "hkf", "water", "IAPWS95", "subcrt",
- "makeup", "basis", "swap.basis", "species", "affinity", "util.affinity", "equil.boltzmann",
+ "makeup", "basis", "swap.basis", "species", "affinity", "equil.boltzmann",
"diagram", "buffer", "nonideal", "add.protein", "protein", "ionize.aa", "more.aa", "read.expr",
"anim", "objective", "revisit", "transfer", "EOSregress", "wjd")
plot.it <- FALSE
@@ -46,3 +46,4 @@
}
return(invisible(out))
}
+
Modified: pkg/CHNOSZ/R/hkf.R
===================================================================
--- pkg/CHNOSZ/R/hkf.R 2017-02-24 13:00:08 UTC (rev 175)
+++ pkg/CHNOSZ/R/hkf.R 2017-02-25 02:08:42 UTC (rev 176)
@@ -175,6 +175,8 @@
return(x)
}
+### unexported functions ###
+
gfun <- function(rhohat, Tc, P, alpha, daldT, beta) {
## g and f functions for describing effective electrostatic radii of ions
## split from hkf() 20120123 jmd
Modified: pkg/CHNOSZ/R/info.R
===================================================================
--- pkg/CHNOSZ/R/info.R 2017-02-24 13:00:08 UTC (rev 175)
+++ pkg/CHNOSZ/R/info.R 2017-02-25 02:08:42 UTC (rev 176)
@@ -6,6 +6,62 @@
# these functions expect arguments of length 1;
# info() handles longer arguments
+info <- function(species=NULL, state=NULL, check.it=TRUE) {
+ ## return information for one or more species in thermo$obigt
+ ## if no species are requested, summarize the available data 20101129
+ thermo <- get("thermo")
+ if(is.null(species)) {
+ message("info: 'species' is NULL; summarizing information about thermodynamic data...")
+ message(paste("thermo$obigt has", nrow(thermo$obigt[thermo$obigt$state=="aq", ]), "aqueous,",
+ nrow(thermo$obigt), "total species"))
+ message(paste("number of literature sources: ", nrow(thermo$refs), ", elements: ",
+ nrow(thermo$element), ", buffers: ", length(unique(thermo$buffers$name)), sep=""))
+ message(paste("number of proteins in thermo$protein is", nrow(thermo$protein), "from",
+ length(unique(thermo$protein$organism)), "organisms"))
+ # print information about SGD.csv, ECO.csv, HUM.csv
+ more.aa(organism="Sce")
+ more.aa(organism="Eco")
+ #pdata.aa(organism="HUM")
+ # print information about yeastgfp.csv
+ yeastgfp()
+ return()
+ }
+ ## run info.numeric or info.character depending on the input type
+ if(is.numeric(species)) {
+ out <- lapply(species, info.numeric, check.it)
+ # if we different states the column names could be different
+ if(length(unique(unlist(lapply(out, names)))) > ncol(thermo$obigt)) {
+ # make them the same as thermo$obigt
+ out <- lapply(out, function(row) {
+ colnames(row) <- colnames(thermo$obigt); return(row)
+ })
+ }
+ # turn the list into a data frame
+ out <- do.call(rbind, out)
+ } else {
+ # state and species should be same length
+ if(!is.null(state)) {
+ lmax <- max(length(species), length(state))
+ state <- rep(state, length.out=lmax)
+ species <- rep(species, length.out=lmax)
+ }
+ # loop over the species
+ out <- sapply(seq_along(species), function(i) {
+ # first look for exact match
+ ispecies <- info.character(species[i], state[i])
+ # if no exact match and it's not a protein, show approximate matches (side effect of info.approx)
+ if(identical(ispecies, NA) & !grepl("_", species[i])) ispecies.notused <- info.approx(species[i], state[i])
+ # do not accept multiple matches
+ if(length(ispecies) > 1) ispecies <- NA
+ return(ispecies)
+ })
+ }
+ ## all done!
+ return(out)
+}
+
+### unexported functions ###
+
info.text <- function(ispecies) {
# a textual description of species name, formula, source, e.g.
# CO2 [CO2(aq)] (SSW01, SHS89, 11.Oct.07)
@@ -173,57 +229,3 @@
return(NA)
}
-info <- function(species=NULL, state=NULL, check.it=TRUE) {
- ## return information for one or more species in thermo$obigt
- ## if no species are requested, summarize the available data 20101129
- thermo <- get("thermo")
- if(is.null(species)) {
- message("info: 'species' is NULL; summarizing information about thermodynamic data...")
- message(paste("thermo$obigt has", nrow(thermo$obigt[thermo$obigt$state=="aq", ]), "aqueous,",
- nrow(thermo$obigt), "total species"))
- message(paste("number of literature sources: ", nrow(thermo$refs), ", elements: ",
- nrow(thermo$element), ", buffers: ", length(unique(thermo$buffers$name)), sep=""))
- message(paste("number of proteins in thermo$protein is", nrow(thermo$protein), "from",
- length(unique(thermo$protein$organism)), "organisms"))
- # print information about SGD.csv, ECO.csv, HUM.csv
- more.aa(organism="Sce")
- more.aa(organism="Eco")
- #pdata.aa(organism="HUM")
- # print information about yeastgfp.csv
- yeastgfp()
- return()
- }
- ## run info.numeric or info.character depending on the input type
- if(is.numeric(species)) {
- out <- lapply(species, info.numeric, check.it)
- # if we different states the column names could be different
- if(length(unique(unlist(lapply(out, names)))) > ncol(thermo$obigt)) {
- # make them the same as thermo$obigt
- out <- lapply(out, function(row) {
- colnames(row) <- colnames(thermo$obigt); return(row)
- })
- }
- # turn the list into a data frame
- out <- do.call(rbind, out)
- } else {
- # state and species should be same length
- if(!is.null(state)) {
- lmax <- max(length(species), length(state))
- state <- rep(state, length.out=lmax)
- species <- rep(species, length.out=lmax)
- }
- # loop over the species
- out <- sapply(seq_along(species), function(i) {
- # first look for exact match
- ispecies <- info.character(species[i], state[i])
- # if no exact match and it's not a protein, show approximate matches (side effect of info.approx)
- if(identical(ispecies, NA) & !grepl("_", species[i])) ispecies.notused <- info.approx(species[i], state[i])
- # do not accept multiple matches
- if(length(ispecies) > 1) ispecies <- NA
- return(ispecies)
- })
- }
- ## all done!
- return(out)
-}
-
Modified: pkg/CHNOSZ/R/makeup.R
===================================================================
--- pkg/CHNOSZ/R/makeup.R 2017-02-24 13:00:08 UTC (rev 175)
+++ pkg/CHNOSZ/R/makeup.R 2017-02-25 02:08:42 UTC (rev 176)
@@ -1,3 +1,139 @@
+# CHNOSZ/makeup.R
+
+makeup <- function(formula, multiplier=1, sum=FALSE, count.zero=FALSE) {
+ # return the elemental makeup (counts) of a chemical formula
+ # that may contain suffixed and/or parenthetical subformulas and/or charge
+ # and negative or positive, fractional coefficients
+ # a matrix is processed by rows
+ if(is.matrix(formula))
+ return(lapply(seq_len(nrow(formula)),
+ function(i) makeup(formula[i, ])))
+ # a named object or list of named objects is returned untouched
+ if(!is.null(names(formula))) return(formula)
+ if(is.list(formula) & !is.null(names(formula[[1]]))) return(formula)
+ # prepare to multiply the formula by the multiplier, if given
+ if(length(multiplier) > 1 & length(multiplier) != length(formula))
+ stop("multiplier does not have length = 1 or length = number of formulas")
+ multiplier <- rep(multiplier, length(formula))
+ # if the formula argument has length > 1, apply the function over each formula
+ if(length(formula) > 1) {
+ # get formulas for any species indices in the argument
+ formula <- get.formula(formula)
+ out <- lapply(seq_along(formula), function(i) {
+ makeup(formula[i], multiplier[i])
+ })
+ # if sum is TRUE, take the sum of all formulas
+ if(sum) {
+ out <- unlist(out)
+ out <- tapply(out, names(out), sum)
+ } else if(count.zero) {
+ # if count.zero is TRUE, all elements appearing in any
+ # of the formulas are counted for each species
+ # first construct elemental makeup showing zero of each element
+ em0 <- unlist(out)
+ em0 <- tapply(em0, names(em0), sum)
+ em0[] <- 0
+ # then sum each formula and the zero vector,
+ # using tapply to group the elements
+ out <- lapply(out, function(x) {
+ xem <- c(x, em0)
+ tapply(xem, names(xem), sum)
+ })
+ }
+ return(out)
+ }
+ # if the formula argument is numeric,
+ # and if the thermo object is available,
+ # get the formula of that numbered species from thermo$obigt
+ if("CHNOSZ" %in% search()) {
+ thermo <- get("thermo", "CHNOSZ")
+ if(is.numeric(formula)) formula <- thermo$obigt$formula[formula]
+ }
+ # first deal with charge
+ cc <- count.charge(formula)
+ # count.elements doesn't know about charge so we need
+ # to explicate the elemental symbol for it
+ formula <- cc$uncharged
+ if(cc$Z != 0 ) formula <- paste(formula, "Z", cc$Z, sep="")
+ # now "Z" will be counted
+ # if there are no subformulas, just use count.elements
+ if(length(grep("(\\(|\\)|\\*|\\:)", formula))==0) {
+ out <- count.elements(formula)
+ } else {
+ # count the subformulas
+ cf <- count.formulas(formula)
+ # count the elements in each subformula
+ ce <- lapply(names(cf), count.elements)
+ # multiply elemental counts by respective subformula counts
+ mcc <- lapply(seq_along(cf), function(i) ce[[i]]*cf[i])
+ # unlist the subformula counts and sum them together by element
+ um <- unlist(mcc)
+ out <- tapply(um, names(um), sum)
+ }
+ # all done with the counting, now apply the multiplier
+ out <- out * multiplier
+ # complain if there are any elements that look strange
+ if("CHNOSZ" %in% search()) {
+ are.elements <- names(out) %in% thermo$element$element
+ if(!all(are.elements)) warning(paste("element(s) not in thermo$element:",
+ paste(rownames(out)[!are.elements], collapse=" ") ))
+ }
+ # done!
+ return(out)
+}
+
+count.elements <- function(formula) {
+ # count the elements in a chemical formula 20120111 jmd
+ # this function expects a simple formula,
+ # no charge or parenthetical or suffixed subformulas
+ # regular expressions inspired by an answer on
+ # http://stackoverflow.com/questions/4116786/parsing-a-chemical-formula-from-a-string-in-c
+ #elementRegex <- "([A-Z][a-z]*)([0-9]*)"
+ elementSymbol <- "([A-Z][a-z]*)"
+ # here, element coefficients can be signed (+ or -) and have a decimal point
+ elementCoeff <- "((\\+|-|\\.|[0-9])*)"
+ elementRegex <- paste(elementSymbol, elementCoeff, sep="")
+ # stop if it doesn't look like a chemical formula
+ validateRegex <- paste("^(", elementRegex, ")+$", sep="")
+ if(length(grep(validateRegex, formula)) == 0)
+ stop(paste("'",formula,"' is not a simple chemical formula", sep="", collapse="\n"))
+ # where to put the output
+ element <- character()
+ count <- numeric()
+ # from now use "f" for formula to make writing the code easier
+ f <- formula
+ # we want to find the starting positions of all the elemental symbols
+ # make substrings, starting at every position in the formula
+ fsub <- sapply(1:nchar(f), function(i) substr(f, i, nchar(f)))
+ # get the numbers (positions) that start with an elemental symbol
+ # i.e. an uppercase letter
+ ielem <- grep("^[A-Z]", fsub)
+ # for each elemental symbol, j is the position before the start of the next
+ # symbol (or the position of the last character of the formula)
+ jelem <- c(tail(ielem - 1, -1), nchar(f))
+ # assemble the stuff: each symbol-coefficient combination
+ ec <- sapply(seq_along(ielem), function(i) substr(f, ielem[i], jelem[i]))
+ # get the individual element symbols and coefficients
+ myelement <- gsub(elementCoeff, "", ec)
+ mycount <- as.numeric(gsub(elementSymbol, "", ec))
+ # any missing coefficients are unity
+ mycount[is.na(mycount)] <- 1
+ # append to the output
+ element <- c(element, myelement)
+ count <- c(count, mycount)
+ # in case there are repeated elements, sum all of their counts
+ # (tapply hint from https://stat.ethz.ch/pipermail/r-help/2011-January/265341.html)
+ out <- tapply(count, element, sum)
+ # tapply returns alphabetical sorted list. keep the order appearing in the formula
+ out <- out[match(unique(element), names(out))]
+ return(out)
+}
+
+### unexported functions ###
+
+# returns a list with named elements
+# `Z` - the numeric value of the charge
+# `uncharged` - the original formula string excluding the charge
count.charge <- function(formula) {
# count the charge in a chemical formula 20120113 jmd
# everything else is counted as the uncharged part of the formula
@@ -34,6 +170,7 @@
return(out)
}
+# returns a numeric vector with names refering to each of the subformulas or the whole formula if there are no subformulas
count.formulas <- function(formula) {
# count the subformulas in a chemical formula 20120112 jmd
# the formula may include multiple unnested parenthetical subformulas
@@ -131,131 +268,3 @@
return(out)
}
-count.elements <- function(formula) {
- # count the elements in a chemical formula 20120111 jmd
- # this function expects a simple formula,
- # no charge or parenthetical or suffixed subformulas
- # regular expressions inspired by an answer on
- # http://stackoverflow.com/questions/4116786/parsing-a-chemical-formula-from-a-string-in-c
- #elementRegex <- "([A-Z][a-z]*)([0-9]*)"
- elementSymbol <- "([A-Z][a-z]*)"
- # here, element coefficients can be signed (+ or -) and have a decimal point
- elementCoeff <- "((\\+|-|\\.|[0-9])*)"
- elementRegex <- paste(elementSymbol, elementCoeff, sep="")
- # stop if it doesn't look like a chemical formula
- validateRegex <- paste("^(", elementRegex, ")+$", sep="")
- if(length(grep(validateRegex, formula)) == 0)
- stop(paste("'",formula,"' is not a simple chemical formula", sep="", collapse="\n"))
- # where to put the output
- element <- character()
- count <- numeric()
- # from now use "f" for formula to make writing the code easier
- f <- formula
- # we want to find the starting positions of all the elemental symbols
- # make substrings, starting at every position in the formula
- fsub <- sapply(1:nchar(f), function(i) substr(f, i, nchar(f)))
- # get the numbers (positions) that start with an elemental symbol
- # i.e. an uppercase letter
- ielem <- grep("^[A-Z]", fsub)
- # for each elemental symbol, j is the position before the start of the next
- # symbol (or the position of the last character of the formula)
- jelem <- c(tail(ielem - 1, -1), nchar(f))
- # assemble the stuff: each symbol-coefficient combination
- ec <- sapply(seq_along(ielem), function(i) substr(f, ielem[i], jelem[i]))
- # get the individual element symbols and coefficients
- myelement <- gsub(elementCoeff, "", ec)
- mycount <- as.numeric(gsub(elementSymbol, "", ec))
- # any missing coefficients are unity
- mycount[is.na(mycount)] <- 1
- # append to the output
- element <- c(element, myelement)
- count <- c(count, mycount)
- # in case there are repeated elements, sum all of their counts
- # (tapply hint from https://stat.ethz.ch/pipermail/r-help/2011-January/265341.html)
- out <- tapply(count, element, sum)
- # tapply returns alphabetical sorted list. keep the order appearing in the formula
- out <- out[match(unique(element), names(out))]
- return(out)
-}
-
-makeup <- function(formula, multiplier=1, sum=FALSE, count.zero=FALSE) {
- # return the elemental makeup (counts) of a chemical formula
- # that may contain suffixed and/or parenthetical subformulas and/or charge
- # and negative or positive, fractional coefficients
- # a matrix is processed by rows
- if(is.matrix(formula))
- return(lapply(seq_len(nrow(formula)),
- function(i) makeup(formula[i, ])))
- # a named object or list of named objects is returned untouched
- if(!is.null(names(formula))) return(formula)
- if(is.list(formula) & !is.null(names(formula[[1]]))) return(formula)
- # prepare to multiply the formula by the multiplier, if given
- if(length(multiplier) > 1 & length(multiplier) != length(formula))
- stop("multiplier does not have length = 1 or length = number of formulas")
- multiplier <- rep(multiplier, length(formula))
- # if the formula argument has length > 1, apply the function over each formula
- if(length(formula) > 1) {
- # get formulas for any species indices in the argument
- formula <- get.formula(formula)
- out <- lapply(seq_along(formula), function(i) {
- makeup(formula[i], multiplier[i])
- })
- # if sum is TRUE, take the sum of all formulas
- if(sum) {
- out <- unlist(out)
- out <- tapply(out, names(out), sum)
- } else if(count.zero) {
- # if count.zero is TRUE, all elements appearing in any
- # of the formulas are counted for each species
- # first construct elemental makeup showing zero of each element
- em0 <- unlist(out)
- em0 <- tapply(em0, names(em0), sum)
- em0[] <- 0
- # then sum each formula and the zero vector,
- # using tapply to group the elements
- out <- lapply(out, function(x) {
- xem <- c(x, em0)
- tapply(xem, names(xem), sum)
- })
- }
- return(out)
- }
- # if the formula argument is numeric,
- # and if the thermo object is available,
- # get the formula of that numbered species from thermo$obigt
- if("CHNOSZ" %in% search()) {
- thermo <- get("thermo", "CHNOSZ")
- if(is.numeric(formula)) formula <- thermo$obigt$formula[formula]
- }
- # first deal with charge
- cc <- count.charge(formula)
- # count.elements doesn't know about charge so we need
- # to explicate the elemental symbol for it
- formula <- cc$uncharged
- if(cc$Z != 0 ) formula <- paste(formula, "Z", cc$Z, sep="")
- # now "Z" will be counted
- # if there are no subformulas, just use count.elements
- if(length(grep("(\\(|\\)|\\*|\\:)", formula))==0) {
- out <- count.elements(formula)
- } else {
- # count the subformulas
- cf <- count.formulas(formula)
- # count the elements in each subformula
- ce <- lapply(names(cf), count.elements)
- # multiply elemental counts by respective subformula counts
- mcc <- lapply(seq_along(cf), function(i) ce[[i]]*cf[i])
- # unlist the subformula counts and sum them together by element
- um <- unlist(mcc)
- out <- tapply(um, names(um), sum)
- }
- # all done with the counting, now apply the multiplier
- out <- out * multiplier
- # complain if there are any elements that look strange
- if("CHNOSZ" %in% search()) {
- are.elements <- names(out) %in% thermo$element$element
- if(!all(are.elements)) warning(paste("element(s) not in thermo$element:",
- paste(rownames(out)[!are.elements], collapse=" ") ))
- }
- # done!
- return(out)
-}
Modified: pkg/CHNOSZ/R/objective.R
===================================================================
--- pkg/CHNOSZ/R/objective.R 2017-02-24 13:00:08 UTC (rev 175)
+++ pkg/CHNOSZ/R/objective.R 2017-02-25 02:08:42 UTC (rev 176)
@@ -2,16 +2,6 @@
# objective functions for revisit(), findit()
# all objective functions in one file, added attributes 20121009 jmd
-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)
-}
-
# 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:
@@ -193,3 +183,17 @@
},
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)
+}
Modified: pkg/CHNOSZ/R/revisit.R
===================================================================
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/chnosz -r 176
More information about the CHNOSZ-commits
mailing list