[CHNOSZ-commits] r571 - in pkg/CHNOSZ: . R inst man vignettes
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Jul 22 05:11:32 CEST 2020
Author: jedick
Date: 2020-07-22 05:11:31 +0200 (Wed, 22 Jul 2020)
New Revision: 571
Added:
pkg/CHNOSZ/R/mix.R
pkg/CHNOSZ/man/mix.Rd
Removed:
pkg/CHNOSZ/R/flatten.R
pkg/CHNOSZ/man/flatten.Rd
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/NAMESPACE
pkg/CHNOSZ/R/diagram.R
pkg/CHNOSZ/inst/NEWS.Rd
pkg/CHNOSZ/vignettes/mklinks.sh
pkg/CHNOSZ/vignettes/multi-metal.Rmd
Log:
Add mix() for mixing diagrams with bimetallic species
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2020-07-21 02:13:09 UTC (rev 570)
+++ pkg/CHNOSZ/DESCRIPTION 2020-07-22 03:11:31 UTC (rev 571)
@@ -1,6 +1,6 @@
-Date: 2020-07-21
+Date: 2020-07-22
Package: CHNOSZ
-Version: 1.3.6-44
+Version: 1.3.6-45
Title: Thermodynamic Calculations and Diagrams for Geochemistry
Authors at R: c(
person("Jeffrey", "Dick", , "j3ffdick at gmail.com", role = c("aut", "cre"),
Modified: pkg/CHNOSZ/NAMESPACE
===================================================================
--- pkg/CHNOSZ/NAMESPACE 2020-07-21 02:13:09 UTC (rev 570)
+++ pkg/CHNOSZ/NAMESPACE 2020-07-22 03:11:31 UTC (rev 571)
@@ -59,7 +59,7 @@
"CHNOSZ", "thermo", "reset", "OBIGT", "retrieve", "AkDi", "moles",
"lNaCl", "lS", "lT", "lP", "lTP", "lex",
# added 20200716 or later
- "duplex", "flatten"
+ "flatten", "mix", "rebalance"
)
# Load shared objects
Modified: pkg/CHNOSZ/R/diagram.R
===================================================================
--- pkg/CHNOSZ/R/diagram.R 2020-07-21 02:13:09 UTC (rev 570)
+++ pkg/CHNOSZ/R/diagram.R 2020-07-22 03:11:31 UTC (rev 571)
@@ -219,6 +219,68 @@
}
}
+ ## create some names for lines/fields if they are missing
+ is.pname <- FALSE
+ onames <- names
+ if(identical(names, FALSE) | identical(names, NA)) names <- ""
+ else if(!is.character(names)) {
+ # properties of basis species or reactions?
+ if(eout$property %in% c("G.basis", "logact.basis")) names <- rownames(eout$basis)
+ else {
+ if(!missing(groups)) {
+ if(is.null(names(groups))) names <- paste("group", 1:length(groups), sep="")
+ else names <- names(groups)
+ }
+ else names <- as.character(eout$species$name)
+ # remove non-unique organism or protein names
+ if(all(grepl("_", names))) {
+ is.pname <- TRUE
+ # 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 state to distinguish ambiguous species names
+ isdup <- names %in% names[duplicated(names)]
+ if(any(isdup)) names[isdup] <- paste(names[isdup],
+ " (", eout$species$state[isdup], ")", sep="")
+ }
+ }
+ # numeric values indicate a subset 20181007
+ if(all(is.numeric(onames))) {
+ if(isTRUE(all(onames > 0))) names[-onames] <- ""
+ else if(isTRUE(all(onames < 0))) names[-onames] <- ""
+ else stop("numeric 'names' should be all positive or all negative")
+ }
+
+ ## apply formatting to chemical formulas 20170204
+ if(all(grepl("_", names))) is.pname <- TRUE
+ if(format.names & !is.pname) {
+ # check if names are a deparsed expression (used in combine()) 20200718
+ parsed <- FALSE
+ if(any(grepl("paste\\(", names))) {
+ exprnames <- parse(text = names)
+ if(length(exprnames) != length(names)) stop("parse()-ing names gives length not equal to number of names")
+ parsed <- TRUE
+ } else {
+ exprnames <- as.expression(names)
+ # get formatted chemical formulas
+ for(i in seq_along(exprnames)) exprnames[[i]] <- expr.species(exprnames[[i]])
+ }
+ # apply bold or italic
+ bold <- rep(bold, length.out = length(exprnames))
+ italic <- rep(italic, length.out = length(exprnames))
+ for(i in seq_along(exprnames)) {
+ if(bold[i]) exprnames[[i]] <- substitute(bold(a), list(a=exprnames[[i]]))
+ if(italic[i]) exprnames[[i]] <- substitute(italic(a), list(a=exprnames[[i]]))
+ }
+ # only use the expression if it's different from the unformatted names
+ if(parsed | !identical(as.character(exprnames), names)) names <- exprnames
+ }
+
## where we'll put extra output for predominance diagrams (namesx, namesy)
out2D <- list()
@@ -238,68 +300,6 @@
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
- is.pname <- FALSE
- onames <- names
- if(identical(names, FALSE) | identical(names, NA)) names <- ""
- else if(!is.character(names)) {
- # properties of basis species or reactions?
- if(eout$property %in% c("G.basis", "logact.basis")) names <- rownames(eout$basis)
- else {
- if(!missing(groups)) {
- if(is.null(names(groups))) names <- paste("group", 1:length(groups), sep="")
- else names <- names(groups)
- }
- else names <- as.character(eout$species$name)
- # remove non-unique organism or protein names
- if(all(grepl("_", names))) {
- is.pname <- TRUE
- # 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 state to distinguish ambiguous species names
- isdup <- names %in% names[duplicated(names)]
- if(any(isdup)) names[isdup] <- paste(names[isdup],
- " (", eout$species$state[isdup], ")", sep="")
- }
- }
- # numeric values indicate a subset 20181007
- if(all(is.numeric(onames))) {
- if(isTRUE(all(onames > 0))) names[-onames] <- ""
- else if(isTRUE(all(onames < 0))) names[-onames] <- ""
- else stop("numeric 'names' should be all positive or all negative")
- }
-
- ## apply formatting to chemical formulas 20170204
- if(all(grepl("_", names))) is.pname <- TRUE
- if(format.names & !is.pname) {
- # check if names are a deparsed expression (used in combine()) 20200718
- parsed <- FALSE
- if(any(grepl("paste\\(", names))) {
- exprnames <- parse(text = names)
- if(length(exprnames) != length(names)) stop("parse()-ing names gives length not equal to number of names")
- parsed <- TRUE
- } else {
- exprnames <- as.expression(names)
- # get formatted chemical formulas
- for(i in seq_along(exprnames)) exprnames[[i]] <- expr.species(exprnames[[i]])
- }
- # apply bold or italic
- bold <- rep(bold, length.out = length(exprnames))
- italic <- rep(italic, length.out = length(exprnames))
- for(i in seq_along(exprnames)) {
- if(bold[i]) exprnames[[i]] <- substitute(bold(a), list(a=exprnames[[i]]))
- if(italic[i]) exprnames[[i]] <- substitute(italic(a), list(a=exprnames[[i]]))
- }
- # only use the expression if it's different from the unformatted names
- if(parsed | !identical(as.character(exprnames), names)) names <- exprnames
- }
-
if(nd==0) {
### 0-D diagram - bar graph of properties of species or reactions
Deleted: pkg/CHNOSZ/R/flatten.R
===================================================================
--- pkg/CHNOSZ/R/flatten.R 2020-07-21 02:13:09 UTC (rev 570)
+++ pkg/CHNOSZ/R/flatten.R 2020-07-22 03:11:31 UTC (rev 571)
@@ -1,138 +0,0 @@
-# CHNOSZ/flatten.R
-# Combine diagrams for two metals
-# 20200713 first version jmd
-
-# Function to combine two diagrams (simple overlay, no interaction) 20200717
-# -- makes new "species" from all combinations of those in d1 and d2
-flatten <- function(d1, d2) {
- check_d1_d2(d1, d2)
-
- # Index all combinations of species in d1 and d2
- i1 <- 1:nrow(d1$species)
- i2 <- 1:nrow(d2$species)
- combs <- expand.grid(i1, i2)
-
- # Get species rows for each combination
- s1 <- d1$species[combs[, 1], ]
- s2 <- d2$species[combs[, 2], ]
- # Make a new species data frame
- nbasis <- nrow(d1$basis)
- species <- s1[, 1:nbasis] + s2[, 1:nbasis]
- ispecies <- paste(s1$ispecies, s2$ispecies, sep = ",")
- logact <- paste(s1$logact, s2$logact, sep = ",")
- state <- paste(s1$state, s2$state, sep = ",")
- # Use names from diagram()
- if(is.expression(d1$names) & is.expression(d2$names)) {
- name <- lapply(1:nrow(combs), function(i) bquote(.(d1$names[[combs[i, 1]]])+.(d2$names[[combs[i, 2]]])))
- name <- unlist(lapply(name, deparse, width.cutoff = 500, control = NULL))
- } else if(is.expression(d1$names)) {
- name <- lapply(1:nrow(combs), function(i) bquote(.(d1$names[[combs[i, 1]]])+.(d2$names[combs[i, 2]])))
- name <- unlist(lapply(name, deparse, width.cutoff = 500, control = NULL))
- } else if(is.expression(d2$names)) {
- name <- lapply(1:nrow(combs), function(i) bquote(.(d1$names[combs[i, 1]])+.(d2$names[[combs[i, 2]]])))
- name <- unlist(lapply(name, deparse, width.cutoff = 500, control = NULL))
- } else name <- paste(d1$names[combs[, 1]], d2$names[combs[, 2]], sep="+")
- if(length(name) != nrow(combs)) stop("deparse()-ing expressions gives unequal length; try diagram(., format.names = FALSE)")
- species <- cbind(species, ispecies, logact, state, name)
-
- # Get affinities for each combination
- v1 <- d1$values[combs[, 1]]
- v2 <- d2$values[combs[, 2]]
- values <- Map("+", v1, v2)
- # Assign -Inf affinity where a species isn't predominant
- for(i in seq_along(values)) {
- i1 <- combs[i, 1]
- i2 <- combs[i, 2]
- ip1 <- d1$predominant == i1
- ip2 <- d2$predominant == i2
- ip12 <- ip1 & ip2
- values[[i]][!ip12] <- -Inf
- }
-
- # Use d1 as a template for the new affinity object
- anew <- d1[1:11]
- # Insert combined results
- anew$species <- species
- anew$values <- values
- # We don't have sout (results from subcrt()) for the combined "species"
- anew$sout <- NULL
- anew
-}
-
-# Function to make a new "affinity" object from two diagrams 20200713
-# -- uses *secondary* balancing coefficients to combine the diagrams
-duplex <- function(d1, d2, balance = NULL) {
- check_d1_d2(d1, d2)
-
- # Combine the species data frames
- species <- rbind(d1$species, d2$species)
- # Combine the sout objects (results from subcrt())
- only2 <- !d2$sout$species$ispecies %in% d1$sout$species$ispecies
- sout <- d1$sout
- sout$species <- rbind(sout$species, d2$sout$species[only2, ])
- sout$out <- c(sout$out, d2$sout$out[only2])
- # Combine the affinity values divided by the *primary*
- # balancing coefficients ("plotvals" from diagram())
- values <- c(d1$plotvals, d2$plotvals)
-
- # Use d1 as a template for the new affinity object
- anew <- d1[1:11]
- # Insert combined results
- anew$species <- species
- anew$sout <- sout
- anew$values <- values
-
- # Figure out the *secondary* balancing coefficients
- n.balance <- balance(anew, balance = balance)$n.balance
- # In the Fe-Cu-S-O-H example all the coefficients on H+ are negative
- if(all(n.balance < 0)) n.balance <- -n.balance
- n1 <- nrow(d1$species)
- n.balance.1 <- n.balance[1:n1]
- n.balance.2 <- n.balance[(n1+1):length(n.balance)]
-
- # Make empty matrices to hold affinities and balancing coefficients
- a1 <- d1$values[[1]]
- a1[] <- NA
- b2 <- a2 <- b1 <- a1
- # Get the affinities (per mole of species, not divided by any balancing coefficients)
- # and the secondary balancing coefficients for the predominant species in each diagram
- p1 <- d1$predominant
- for(ip in unique(as.vector(p1))) {
- a1[p1 == ip] <- d1$values[[ip]][p1 == ip]
- b1[p1 == ip] <- n.balance.1[ip]
- }
- p2 <- d2$predominant
- for(ip in unique(as.vector(p2))) {
- a2[p2 == ip] <- d2$values[[ip]][p2 == ip]
- b2[p2 == ip] <- n.balance.2[ip]
- }
- # Divide the affinities by the secondary balancing coefficients
- ab1 <- a1 / b1
- ab2 <- a2 / b2
- # Identify the species with the highest affinity (predominant in the *secondary* reactions)
- i1 <- ab1 > ab2
- # Suppress non-predominant species at each grid point
- for(i in 1:n1) anew$values[[i]][!i1] <- -Inf
- for(i in (n1+1):length(n.balance)) anew$values[[i]][i1] <- -Inf
-
- anew
-
-}
-
-### unexported function ###
-
-# Check that d1 and d2 can be combined
-# Extracted from duplex() 20200717
-check_d1_d2 <- function(d1, d2) {
- # Check that the basis species are the same
- if(!identical(d1$basis, d2$basis)) stop("basis species in objects 'd1' and 'd2' are not identical")
- # Check that the variables and their values are the same
- if(!identical(d1$vars, d2$vars)) stop("variable names in objects 'd1' and 'd2' are not identical")
- if(!identical(d1$vals, d2$vals)) stop("variable values in objects 'd1' and 'd2' are not identical")
- # Check that T and P are the same
- if(!identical(d1$T, d2$T)) stop("temperatures in objects 'd1' and 'd2' are not identical")
- if(!identical(d1$P, d2$P)) stop("pressures in objects 'd1' and 'd2' are not identical")
- # Check that we have plotvals and predominant (from diagram())
- if(is.null(d1$plotvals) | is.null(d1$predominant)) stop("object 'd1' is missing 'plotvals' or 'predominant' components (not made by diagram()?)")
- if(is.null(d2$plotvals) | is.null(d2$predominant)) stop("object 'd2' is missing 'plotvals' or 'predominant' components (not made by diagram()?)")
-}
Copied: pkg/CHNOSZ/R/mix.R (from rev 570, pkg/CHNOSZ/R/flatten.R)
===================================================================
--- pkg/CHNOSZ/R/mix.R (rev 0)
+++ pkg/CHNOSZ/R/mix.R 2020-07-22 03:11:31 UTC (rev 571)
@@ -0,0 +1,254 @@
+# CHNOSZ/flatten.R
+# Combine diagrams for two metals
+# 20200713 first version jmd
+
+# Function to combine two diagrams (simple overlay, no interaction) 20200717
+# -- makes new "species" from all combinations of those in d1 and d2
+flatten <- function(d1, d2) {
+ # It's just mixing d1 and d2 (two single-metal diagrams) without adding d3 (bimetal) 20200722
+ mix(d1, d2)
+}
+
+# mix the systems to include bimetal species 20200721
+# 'parts' gives number of moles of each metal in 1 mole of the mixture
+mix <- function(d1, d2, d3 = NULL, parts = c(1, 1), .balance = NULL) {
+ check_d1_d2(d1, d2, .balance)
+
+ # Handle mixing here
+ if(!is.null(d3)) {
+ # Mix d1 (e.g. Fe only) and d2 (e.g. V only)
+ m12 <- mix(d1, d2, parts = parts)
+ # Mix d1 (Fe) and d3 (FeV bimetallic species)
+ # - '.balance' is defined to add d3 to get required amount of V
+ m13 <- mix(d1, d3, parts = parts, .balance = d2$balance)
+ # Mix d2 (V) and d3 (FeV)
+ # - '.balance' is defined to add d3 to get required amount of Fe
+ # - d2 (V) is first, so we need to reverse the 'parts' values
+ m23 <- mix(d2, d3, parts = rev(parts), .balance = d1$balance)
+ # Mix d3 with itself (combinations of bimetallic species)
+ m33 <- mix(d3, d3, parts = parts, .balance = c(d1$balance, d2$balance))
+ # Merge all the species and affinity values
+ species <- rbind(m12$species, m13$species, m23$species, m33$species)
+ values <- c(m12$values, m13$values, m23$values, m33$values)
+ # Remove duplicates
+ # (i.e. bimetallic species that exactly match the composition in 'parts'
+ # and therefore appear in multiple combinations with
+ # mono-metallic species that have zero mole fractions)
+ idup <- duplicated(species$name)
+ species <- species[!idup, ]
+ values <- values[!idup]
+ # Put together the output
+ anew <- m12
+ anew$species <- species
+ anew$values <- values
+ return(anew)
+ }
+
+ # Index all combinations of species in d1 and d2
+ i1 <- 1:nrow(d1$species)
+ i2 <- 1:nrow(d2$species)
+ combs <- expand.grid(i1, i2)
+ # For bimetallic species (m33)
+ if(length(.balance)==2) {
+ # Remove combinations of a species with itself
+ combs <- combs[combs[, 1] != combs[, 2], ]
+ # Remove duplicated combinations
+ isdup <- duplicated(lapply(apply(combs, 1, sort, simplify = FALSE), paste))
+ combs <- combs[!isdup, ]
+ }
+ # Get species rows for each combination
+ s1 <- d1$species[combs[, 1], ]
+ s2 <- d2$species[combs[, 2], ]
+ # Get balancing coefficients for each combination
+ b1 <- d1$n.balance[combs[, 1]]
+ b2 <- d2$n.balance[combs[, 2]]
+ # Locate the columns for the two metals in the formation reactions
+ ibal <- match(c(d1$balance, d2$balance), colnames(s1))
+ # With non-null '.balance', d2 is really d3 so we calculate the needed balancing coefficients
+ if(!is.null(.balance)) {
+ if(length(.balance)==2) {
+ # For combinations of bimetallic species (m33)
+ b1 <- balance(d1, .balance[1])$n.balance[combs[, 1]]
+ b2 <- balance(d2, .balance[2])$n.balance[combs[, 2]]
+ ibal <- match(balance, colnames(s1))
+ } else {
+ b2 <- balance(d2, .balance)$n.balance[combs[, 2]]
+ ibal <- match(c(d1$balance, .balance), colnames(s1))
+ }
+ }
+ # Solve for the mole fractions of each species that give the required mixture
+ frac1 <- frac2 <- numeric()
+ for(i in 1:nrow(combs)) {
+ # The stoichiometric matrix
+ A <- matrix(as.numeric(c(s1[i, ibal], s2[i, ibal])), nrow = 2, byrow = TRUE)
+ x <- solve(t(A), parts)
+ frac1 <- c(frac1, x[1])
+ frac2 <- c(frac2, x[2])
+ }
+ # Note that some of frac1 might be < 0 ... we remove these combinations below
+
+ # Make a new species data frame starting with sum of scaled formation reactions
+ nbasis <- nrow(d1$basis)
+ species <- s1[, 1:nbasis] * frac1 + s2[, 1:nbasis] * frac2
+ ispecies <- paste(s1$ispecies, s2$ispecies, sep = ",")
+ logact <- paste(s1$logact, s2$logact, sep = ",")
+ state <- paste(s1$state, s2$state, sep = ",")
+ # Use names from diagram()
+ if(is.expression(d1$names) | is.expression(d2$names)) {
+ # Convert non-expressions to lists so we can use [[ indexing below
+ if(!is.expression(d1$names)) d1$names <- as.list(d1$names)
+ if(!is.expression(d2$names)) d2$names <- as.list(d2$names)
+ name <- lapply(1:nrow(combs), function(i) {
+ # Don't include names of species that are not present (zero or negative mole fractions)
+ if(frac1[i] <= 0) bquote(.(d2$names[[combs[i, 2]]]))
+ else if(frac2[i] <= 0) bquote(.(d1$names[[combs[i, 1]]]))
+ else bquote(.(d1$names[[combs[i, 1]]])+.(d2$names[[combs[i, 2]]]))
+ })
+ name <- unlist(lapply(name, deparse, width.cutoff = 500, control = NULL))
+ if(length(name) != nrow(combs)) stop("deparse()-ing expressions gives unequal length; try diagram(., format.names = FALSE)")
+ } else {
+ # Plain text names
+ name <- sapply(1:nrow(combs), function(i) {
+ if(frac1[i] <= 0) d2$names[combs[i, 2]]
+ else if(frac2[i] <= 0) d1$names[combs[i, 1]]
+ else paste(d1$names[combs[i, 1]], d2$names[combs[i, 2]], sep="+")
+ })
+ }
+ species <- cbind(species, ispecies, logact, state, name)
+
+ # Get affinities for each combination of species
+ v1 <- d1$values[combs[, 1]]
+ v2 <- d2$values[combs[, 2]]
+ # Scale affinities by mole fractions computed for the mixture
+ v1 <- Map("*", v1, as.list(frac1))
+ v2 <- Map("*", v2, as.list(frac2))
+ # Add together the scaled affinities
+ values <- Map("+", v1, v2)
+ ipredominant <- logical(length(values))
+ if(length(.balance)==2) {
+ # For combinations of bimetallic species (m33), don't do predominance masking
+ # (there are not mono-metallic species to look at)
+ ipredominant[] <- TRUE
+ } else {
+ # Loop over combinations to find predominant species in the single-metal diagram(s)
+ for(i in seq_along(values)) {
+ # Get predominant species in first diagram
+ i1 <- combs[i, 1]
+ ip1 <- d1$predominant == i1
+ # If the mole fraction is zero, it is predominant by definition
+ # (this allows a single bimetallic species to be formed)
+ if(frac1[i] == 0) ip1[] <- TRUE
+ if(is.null(.balance)) {
+ # Get predominant species in second diagram
+ i2 <- combs[i, 2]
+ ip2 <- d2$predominant == i2
+ # If the mole fraction is zero, it is predominant by definition
+ # (this allows us to recover the single-element diagram with frac = c(0, 1) or c(1, 0))
+ if(frac2[i] == 0) ip2[] <- TRUE
+ # Assign -Inf affinity where any species isn't predominant in the corresponding single-metal diagram
+ ip12 <- ip1 & ip2
+ values[[i]][!ip12] <- -Inf
+ if(any(ip12)) ipredominant[i] <- TRUE
+ } else {
+ # For non-null '.balance', only the first diagram is for a single metal
+ values[[i]][!ip1] <- -Inf
+ if(any(ip1)) ipredominant[i] <- TRUE
+ }
+ }
+ }
+ # Remove combinations that:
+ # 1) involve a species with no predominance field in the single-metal diagram or
+ # 2) have a negative mole fraction of a single-metal species
+ inotnegative <- frac1 >= 0
+ values <- values[ipredominant & inotnegative]
+ species <- species[ipredominant & inotnegative, ]
+
+ # Use d1 as a template for the new affinity object
+ anew <- d1[1:11]
+ # Insert combined results
+ anew$species <- species
+ anew$values <- values
+ # We don't have sout (results from subcrt()) for the combined "species"
+ anew$sout <- NULL
+ anew
+}
+
+# Function to make a new "affinity" object from two diagrams 20200713
+# -- uses *secondary* balancing coefficients to combine the diagrams
+rebalance <- function(d1, d2, balance = NULL) {
+ check_d1_d2(d1, d2)
+
+ # Combine the species data frames
+ species <- rbind(d1$species, d2$species)
+ # Combine the sout objects (results from subcrt())
+ only2 <- !d2$sout$species$ispecies %in% d1$sout$species$ispecies
+ sout <- d1$sout
+ sout$species <- rbind(sout$species, d2$sout$species[only2, ])
+ sout$out <- c(sout$out, d2$sout$out[only2])
+ # Combine the affinity values divided by the *primary*
+ # balancing coefficients ("plotvals" from diagram())
+ values <- c(d1$plotvals, d2$plotvals)
+
+ # Use d1 as a template for the new affinity object
+ anew <- d1[1:11]
+ # Insert combined results
+ anew$species <- species
+ anew$sout <- sout
+ anew$values <- values
+
+ # Figure out the *secondary* balancing coefficients
+ n.balance <- balance(anew, balance = balance)$n.balance
+ # In the Fe-Cu-S-O-H example all the coefficients on H+ are negative
+ if(all(n.balance < 0)) n.balance <- -n.balance
+ n1 <- nrow(d1$species)
+ n.balance.1 <- n.balance[1:n1]
+ n.balance.2 <- n.balance[(n1+1):length(n.balance)]
+
+ # Make empty matrices to hold affinities and balancing coefficients
+ a1 <- d1$values[[1]]
+ a1[] <- NA
+ b2 <- a2 <- b1 <- a1
+ # Get the affinities (per mole of species, not divided by any balancing coefficients)
+ # and the secondary balancing coefficients for the predominant species in each diagram
+ p1 <- d1$predominant
+ for(ip in unique(as.vector(p1))) {
+ a1[p1 == ip] <- d1$values[[ip]][p1 == ip]
+ b1[p1 == ip] <- n.balance.1[ip]
+ }
+ p2 <- d2$predominant
+ for(ip in unique(as.vector(p2))) {
+ a2[p2 == ip] <- d2$values[[ip]][p2 == ip]
+ b2[p2 == ip] <- n.balance.2[ip]
+ }
+ # Divide the affinities by the secondary balancing coefficients
+ ab1 <- a1 / b1
+ ab2 <- a2 / b2
+ # Identify the species with the highest affinity (predominant in the *secondary* reactions)
+ i1 <- ab1 > ab2
+ # Suppress non-predominant species at each grid point
+ for(i in 1:n1) anew$values[[i]][!i1] <- -Inf
+ for(i in (n1+1):length(n.balance)) anew$values[[i]][i1] <- -Inf
+
+ anew
+
+}
+
+
+### unexported function ###
+
+# Check that d1 and d2 can be combined
+# Extracted from duplex() (now rebalance()) 20200717
+check_d1_d2 <- function(d1, d2, .balance = NULL) {
+ if(is.null(.balance)) d2txt <- "d2" else d2txt <- "d3"
+ # Check that the basis species are the same
+ if(!identical(d1$basis, d2$basis)) stop(paste0("basis species in objects 'd1' and '", d2txt, "' are not identical"))
+ # Check that the variables and their values are the same
+ if(!identical(d1$vars, d2$vars)) stop(paste0("variable names in objects 'd1' and '", d2txt, "' are not identical"))
+ if(!identical(d1$vals, d2$vals)) stop(paste0("variable values in objects 'd1' and '", d2txt, "' are not identical"))
+ # Check that T and P are the same
+ if(!identical(d1$T, d2$T)) stop(paste0("temperatures in objects 'd1' and '", d2txt, "' are not identical"))
+ if(!identical(d1$P, d2$P)) stop(paste0("pressures in objects 'd1' and '", d2txt, "' are not identical"))
+ # Check that we have plotvals and predominant (from diagram())
+ if(is.null(d1$plotvals) | is.null(d1$predominant)) stop("object 'd1' is missing 'plotvals' or 'predominant' components (not made by diagram()?)")
+ if(is.null(d2$plotvals) | is.null(d2$predominant)) stop(paste0("object '", d2txt, "' is missing 'plotvals' or 'predominant' components (not made by diagram()?)"))
+}
Modified: pkg/CHNOSZ/inst/NEWS.Rd
===================================================================
--- pkg/CHNOSZ/inst/NEWS.Rd 2020-07-21 02:13:09 UTC (rev 570)
+++ pkg/CHNOSZ/inst/NEWS.Rd 2020-07-22 03:11:31 UTC (rev 571)
@@ -6,7 +6,7 @@
\newcommand{\H2O}{\ifelse{latex}{\eqn{\mathrm{H_{2}O}}}{\ifelse{html}{\out{H<sub>2</sub>O}}{H2O}}}
\newcommand{\Hplus}{\ifelse{latex}{\eqn{\mathrm{H^{+}}}}{\ifelse{html}{\out{H<sup>+</sup>}}{H+}}}
-\section{Changes in CHNOSZ version 1.3.6-43 (2020-07-21)}{
+\section{Changes in CHNOSZ version 1.3.6-45 (2020-07-22)}{
\subsection{MAJOR CHANGES}{
\itemize{
@@ -49,12 +49,16 @@
\itemize{
\item Add function \strong{flatten()} for combining two diagrams for
- different systems (i.e., simple overlay with labels for species from both
+ different systems (i.e., simple overlay of diagrams for two single-metal
systems).
- \item Add function \strong{duplex()} for making a new diagram by secondary
- balancing between two systems.
+ \item Add function \strong{mix()} for combining two single-metal diagrams
+ with a third diagram for bimetallic species. This can be used to produce
+ diagrams for a binary system with fixed composition of the metals.
+ \item Add function \strong{rebalance()} for making a new diagram by
+ secondary balancing between two systems.
+
\item Add a \strong{predominant} argument to \code{mosaic()} to use
previously calculated predominances of species (e.g. minerals) for the
changing basis species. This allows \code{mosaic()} calculations to be
@@ -218,6 +222,8 @@
\item TODO: for balance = 1, add or change message text to "formula units".
+ \item TODO: remove species combinations with no stability field in flatten().
+
}
}
Deleted: pkg/CHNOSZ/man/flatten.Rd
===================================================================
--- pkg/CHNOSZ/man/flatten.Rd 2020-07-21 02:13:09 UTC (rev 570)
+++ pkg/CHNOSZ/man/flatten.Rd 2020-07-22 03:11:31 UTC (rev 571)
@@ -1,74 +0,0 @@
-\encoding{UTF-8}
-\name{flatten}
-\alias{flatten}
-\alias{duplex}
-\title{Combine Diagrams}
-\description{
- Combine two diagrams for different systems by flattening them or using a secondary balancing constraint.
-}
-
-\usage{
- flatten(d1, d2)
- duplex(d1, d2, balance = NULL)
-}
-
-\arguments{
- \item{d1}{list, output of \code{\link{diagram}} for first system}
- \item{d2}{list, output of \code{diagram} for second system}
- \item{balance}{character or numeric, specification of secondary balancing coefficients}
-}
-
-\details{
-
-These functions both make a new \code{\link{affinity}} object from two diagrams.
-\code{flatten} identifies the intersection of predominance fields for all possible combinations of species (without interaction between the systems), while \code{duplex} creates a new diagram by comparing the affinities of reactions between species in both systems.
-Both functions mask the non-predominant species by assigning them -Inf values of affinity, so the result can be used to make a new diagram that shows the combined system.
-
-\code{flatten} makes a simple overlay of two diagrams using new "species" generated from all combinations of those in \code{d1} and \code{d2}.
-The new names are formed from the \code{names} used in the source diagrams; for example if "Cp" and "Py" are predominant minerals at the same position in diagrams 1 and 2, the field for the flattened diagram will be labeled "Cp+Py".
-The resulting affinities are simply the sum of affinities of the two species; they are assigned values of -Inf wherever one of the species is not predominant in either of the diagrams.
-
-Diagrams for different systems likely use different \emph{primary} balancing coefficients, such as balancing on different metals.
-\code{duplex} uses \emph{secondary} balancing coefficients, specified acording to \code{balance} (see \code{\link{equilibrate}} for a description of this argument), to determine the reactions between the species in the two systems.
-The affinities of these reactions are then used \emph{only} to identify the predominant species at each grid point.
-The \emph{returned} value of affinity are carried forward from those used to make the source diagrams (\samp{plotvals} in \code{d1} and \code{d2}), and therefore reflect the primary balancing coefficients.
-The returned values are assigned -Inf wherever that species is determined to not predominate according to the secondary balancing.
-
-Because \code{flatten} yields finite values of affinity for only a single species at any grid point, the final diagram can be made with any setting of \code{balance}.
-However, for \code{duplex}, \code{balance} in the final diagram should be set to \samp{1} in order to preserve the primary balancing coefficients.
-
-}
-
-\value{
-A list object with the same structure as the output from \code{\link{affinity}}, so it can be used as input to \code{diagram}.
-}
-
-\seealso{
-A longer example is in the vignette \viglink{multi-metal}.
-}
-
-\examples{\dontshow{opar <- par(no.readonly = TRUE)}
-par(mfrow = c(2, 2))
-# Define basis species with Fe and Cu
-basis(c("Fe+2", "Cu+", "hydrogen sulfide", "oxygen", "H2O", "H+"))
-xlab <- ratlab("Fe+2", "Cu+")
-# Calculate diagram for only Fe-bearing minerals
-species(c("pyrite", "pyrrhotite", "magnetite", "hematite"))
-aFe <- affinity("Fe+2" = c(0, 12), O2 = c(-40, -16), T = 400, P = 2000)
-dFe <- diagram(aFe, xlab = xlab, main = "Fe-S-O-H")
-# Calculate diagram for only Cu-bearing minerals
-species(c("covellite", "chalcocite", "tenorite", "cuprite"))
-aCu <- affinity(aFe) # argument recall
-dCu <- diagram(aCu, xlab = xlab, main = "Cu-S-O-H")
-### flatten() diagram
-ac <- flatten(dFe, dCu)
-diagram(ac, xlab = xlab, main = "Cu-Fe-S-O-H with flatten()")
-### duplex() diagram
-ad <- duplex(dFe, dCu)
-diagram(ad, xlab = xlab, balance = 1, main = "Cu-Fe-S-O-H with duplex()")
-db <- describe.basis(ibasis = 3)
-leg <- lex(lTP(400, 2000), db)
-legend("bottomleft", legend = leg, bty = "n")
-\dontshow{par(opar)}}
-
-\concept{Extended workflow}
Copied: pkg/CHNOSZ/man/mix.Rd (from rev 570, pkg/CHNOSZ/man/flatten.Rd)
===================================================================
--- pkg/CHNOSZ/man/mix.Rd (rev 0)
+++ pkg/CHNOSZ/man/mix.Rd 2020-07-22 03:11:31 UTC (rev 571)
@@ -0,0 +1,96 @@
+\encoding{UTF-8}
+\name{mix}
+\alias{flatten}
+\alias{rebalance}
+\alias{mix}
+\title{Combine Diagrams for Multi-Metal Systems}
+\description{
+ Combine diagrams for different systems by flattening or rebalancing two diagrams or mixing two diagrams with a third.
+}
+
+\usage{
+ flatten(d1, d2)
+ rebalance(d1, d2, balance = NULL)
+ mix(d1, d2, d3, parts = c(1, 1), .balance = NULL)
+}
+
+\arguments{
+ \item{d1}{list, output of \code{\link{diagram}} for the first mono-metallic system}
+ \item{d2}{list, output of \code{diagram} for the second mono-metallic system}
+ \item{balance}{character or numeric, specification of secondary balancing coefficients}
+ \item{d3}{list, output of \code{diagram} for the bimetallic system}
+ \item{parts}{numeric, amount of each metal (i.e. fixed composition) for the mixed system}
+ \item{.balance}{\emph{argument for internal use only}}
+}
+
+\details{
+
+These functions make a new \code{\link{affinity}} object from the output of \code{\link{diagram}}.
+The result can be used to make a new diagram that shows the combined system.
+
+\code{flatten} creates a set of intersecting predominance fields for all possible combinations of species in \code{d1} and \code{d2}.
+The new names are formed from the \code{names} used in the source diagrams; for example if "Cp" and "Py" are predominant minerals at the same position in diagrams 1 and 2, the field for the flattened diagram will be labeled "Cp+Py".
+The affinities are calculated by summing the formation reactions from the two diagrams to give equal parts of the balancing coefficients in \code{d1} and \code{d2} (that is, equal parts of two different metals).
+Note that the actual values of the affinities (and therefore the ratio between the metals) doesn't affect the resulting diagram because the affinities are assigned values of -Inf wherever one of the species is not predominant in the respective single-metal diagram.
+
+\code{mix} is an expanded form of \code{flatten} that allows combinations not only between two single-metal diagrams (\code{d1} and \code{d2}) but also between each of those diagrams and third diagram for bimetallic species (\code{d3}).
+All combinations of species in all crosses between the diagrams (\code{d1-d2}, \code{d1-d3}, \code{d2-d3}, \code{d3-d3}) are identified.
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/chnosz -r 571
More information about the CHNOSZ-commits
mailing list