From 649ff31773c351b921d3b3554131fb1fdc96e646 Mon Sep 17 00:00:00 2001 From: Jeffrey Dick Date: Wed, 9 Aug 2023 04:02:36 +0000 Subject: [PATCH] mosaic() can now change basis species that are axis variables on a diagram git-svn-id: svn://scm.r-forge.r-project.org/svnroot/chnosz/pkg/CHNOSZ@799 edb9625f-4e0d-4859-8d74-9fd3b1da38cb --- DESCRIPTION | 2 +- R/diagram.R | 18 ++++++++++++--- R/mosaic.R | 54 +++++++++++++++++++++++++++++++++++++-------- R/util.expression.R | 9 +++++++- inst/NEWS.Rd | 7 +++++- 5 files changed, 75 insertions(+), 15 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 905f8b7..822f6f4 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Date: 2023-08-09 Package: CHNOSZ -Version: 2.0.0-18 +Version: 2.0.0-19 Title: Thermodynamic Calculations and Diagrams for Geochemistry Authors@R: c( person("Jeffrey", "Dick", , "j3ffdick@gmail.com", role = c("aut", "cre"), diff --git a/R/diagram.R b/R/diagram.R index b0b5b0d..e59ef46 100644 --- a/R/diagram.R +++ b/R/diagram.R @@ -333,6 +333,18 @@ diagram <- function( lty <- rep(lty, length.out = length(plotvals)) lwd <- rep(lwd, length.out = length(plotvals)) col <- rep(col, length.out = length(plotvals)) + + # Function to get label for i'th variable 20230809 + # (uses custom labels from 'labels' list element added by mosaic) + getlabel <- function(ivar) { + label <- eout$vars[ivar] + if(!is.null(eout$labels)) { + if(label %in% names(eout$labels)) { + label <- eout$labels[[label]] + } + } + label + } if(nd == 0) { @@ -349,7 +361,7 @@ diagram <- function( if(missing(xlim)) xlim <- range(xvalues) # TODO: this is backward if the vals are not increasing # Initialize the plot if(!add) { - if(missing(xlab)) xlab <- axis.label(eout$vars[1], basis = eout$basis, molality = molality) + if(missing(xlab)) xlab <- axis.label(getlabel(1), basis = eout$basis, molality = molality) if(missing(ylab)) { ylab <- axis.label(plotvar, units = "", molality = molality) if(plotvar == "rank.affinity") ylab <- "Average affinity ranking" @@ -695,8 +707,8 @@ diagram <- function( } # Initialize the plot if(!add) { - if(is.null(xlab)) xlab <- axis.label(eout$vars[1], basis = eout$basis, molality = molality) - if(is.null(ylab)) ylab <- axis.label(eout$vars[2], basis = eout$basis, molality = molality) + if(is.null(xlab)) xlab <- axis.label(getlabel(1), basis = eout$basis, molality = molality) + if(is.null(ylab)) ylab <- axis.label(getlabel(2), basis = eout$basis, molality = molality) if(tplot) thermo.plot.new(xlim = xlim, ylim = ylim, xlab = xlab, ylab = ylab, cex = cex, cex.axis = cex.axis, mar = mar, yline = yline, side = side, ...) else plot(0, 0, type = "n", xlim = xlim, ylim = ylim, xlab = xlab, ylab = ylab, ...) diff --git a/R/mosaic.R b/R/mosaic.R index 7fc9011..f2f8157 100644 --- a/R/mosaic.R +++ b/R/mosaic.R @@ -69,7 +69,7 @@ mosaic <- function(bases, blend = TRUE, stable = list(), loga_aq = NULL, ...) { ina <- is.na(ibasis0) if(any(ina)) { names0 <- unlist(lapply(bases, "[", 1)) - stop("the starting basis does not have ", paste(names0[ina], collapse = " and ")) + stop("the starting basis species do not have ", paste(names0[ina], collapse = " and ")) } if("sout" %in% names(affinityargs)) { @@ -85,7 +85,6 @@ mosaic <- function(bases, blend = TRUE, stable = list(), loga_aq = NULL, ...) { sout <- do.call(affinity, c(affinityargs, list(return.sout = TRUE))) } - # Calculate affinities of the basis species themselves A.bases <- list() for(i in 1:length(bases)) { @@ -99,15 +98,40 @@ mosaic <- function(bases, blend = TRUE, stable = list(), loga_aq = NULL, ...) { else A.bases[[i]] <- suppressMessages(do.call(affinity, c(affinityargs, list(sout = sout)))) } - # Get all combinations of basis species + # Get all combinations of basis species (species indices in OBIGT) newbases <- as.matrix(expand.grid(ispecies)) allbases <- matrix(basis0$ispecies, nrow = 1)[rep(1, nrow(newbases)), , drop = FALSE] allbases[, ibasis0] <- newbases + # Also get all combinations of names of basis species (for modifying affinityargs) 20230809 + newbnames <- as.matrix(expand.grid(bases)) + allbnames <- matrix(rownames(basis0), nrow = 1)[rep(1, nrow(newbnames)), , drop = FALSE] + allbnames[, ibasis0] <- newbnames + # Look for argument names for affinity() in starting basis species + # (i.e., basis species that are variables on the diagram) + matches.bnames <- names(affinityargs) %in% allbnames[1, ] + # Find the name(s) of the starting basis species that are variables on the diagram + ibnames <- match(names(affinityargs)[matches.bnames], allbnames[1, ]) + # Figure out the element to make labels (total C, total S, etc.) + labels <- NULL + if(any(matches.bnames)) { + element.matrix <- basis0[, 1:nrow(basis0)] + elements.in.basis0 <- colSums(element.matrix) + labelnames <- allbnames[1, ibnames] + labels <- lapply(1:length(labelnames), function(i) { + has.element <- element.matrix[match(labelnames[i], rownames(element.matrix)), ] > 0 + ielement <- has.element & elements.in.basis0 == 1 + # Use the element or fallback to species name if element isn't found + if(any(ielement)) colnames(element.matrix)[ielement][1] + else labelnames[i] + }) + names(labels) <- labelnames + } + # Calculate affinities of species for all combinations of basis species aff.species <- list() message("mosaic: calculating affinities of species for all ", nrow(allbases), " combinations of the basis species") - # Run backwards so that we put the starting basis species back at the end + # Run backwards so that we end up with the starting basis species for(i in nrow(allbases):1) { # Get default loga from starting basis species thislogact <- basis0$logact @@ -127,10 +151,19 @@ mosaic <- function(bases, blend = TRUE, stable = list(), loga_aq = NULL, ...) { } } put.basis(allbases[i, ], thislogact) - # We have to define the species using the current basis + # Load the formed species using the current basis species(species0$ispecies, species0$logact) - if(affinityargs_has_sout) aff.species[[i]] <- suppressMessages(do.call(affinity, affinityargs)) - else aff.species[[i]] <- suppressMessages(do.call(affinity, c(affinityargs, list(sout = sout)))) + + # If mosaic() changes variables on the diagram, argument names for affinity() also have to be changed 20230809 + myaffinityargs <- affinityargs + if(any(matches.bnames)) { + # At least one basis species in 'bases' is a variable on the diagram + # Use the name of the current swapped-in basis species + names(myaffinityargs)[matches.bnames] <- allbnames[i, ibnames] + } + + if(affinityargs_has_sout) aff.species[[i]] <- suppressMessages(do.call(affinity, myaffinityargs)) + else aff.species[[i]] <- suppressMessages(do.call(affinity, c(myaffinityargs, list(sout = sout)))) } # Calculate equilibrium mole fractions for each group of basis species @@ -141,8 +174,8 @@ mosaic <- function(bases, blend = TRUE, stable = list(), loga_aq = NULL, ...) { if(blend[i] & is.null(stable[i][[1]])) { # This isn't needed (and doesn't work) if all the affinities are NA 20180925 if(any(!sapply(A.bases[[1]]$values, is.na))) { - # 20190504: when equilibrating the changing basis species, use a total activity equal to the activity from the basis definition - # 20191111 use equilibrate(loga.balance = ) instead of setting activities in species definition + # When equilibrating the changing basis species, use a total activity equal to the activity from the basis definition 20190504 + # Use equilibrate(loga.balance = ) instead of setting activities in species definition 20191111 e <- equilibrate(A.bases[[i]], loga.balance = as.numeric(basis0$logact[ibasis0[i]])) # Exponentiate to get activities then divide by total activity a.equil <- lapply(e$loga.equil, function(x) 10^x) @@ -211,6 +244,9 @@ mosaic <- function(bases, blend = TRUE, stable = list(), loga_aq = NULL, ...) { A.species$values[[i]] <- Reduce("+", A.values) } + # Insert custom labels 20230809 + A.species$labels <- labels + # For argument recall, include all arguments in output 20190120 allargs <- c(list(bases = bases, blend = blend), affinityargs) # Return the affinities for the species and basis species diff --git a/R/util.expression.R b/R/util.expression.R index 664ebb8..f95842b 100644 --- a/R/util.expression.R +++ b/R/util.expression.R @@ -159,14 +159,20 @@ axis.label <- function(label, units = NULL, basis = thermo()$basis, prefix = "", # Make a formatted axis label from a generic description # It can be a chemical property, condition, or chemical activity in the system; # if the label matches one of the basis species or if the state is specified, it's a chemical activity + # 20090826: Just return the argument if a comma is already present - # (it's good for custom labels that shouldn't be italicized) + # (used for custom labels that shouldn't be italicized) if(grepl(",", label)) return(label) + if(label %in% rownames(basis)) { # 20090215: The state this basis species is in state <- basis$state[match(label, rownames(basis))] # Get the formatted label desc <- expr.species(label, state = state, log = TRUE, molality = molality) + } else if(label %in% colnames(basis)) { + # Make a label for an element (total C, total S, etc.) 20230809 + if(molality) desc <- bquote(log~italic(m)~"(total "*.(label)*")") + else desc <- bquote(log~italic(a)~"(total "*.(label)*")") } else { # The label is for a chemical property or condition # Make the label by putting a comma between the property and the units @@ -176,6 +182,7 @@ axis.label <- function(label, units = NULL, basis = thermo()$basis, prefix = "", if(units == "") desc <- substitute(a, list(a = property)) else desc <- substitute(a~"("*b*")", list(a = property, b = units)) } + # Done! return(desc) } diff --git a/inst/NEWS.Rd b/inst/NEWS.Rd index bd676fa..90a8ec1 100644 --- a/inst/NEWS.Rd +++ b/inst/NEWS.Rd @@ -12,7 +12,7 @@ % links to vignettes 20220723 \newcommand{\viglink}{\ifelse{html}{\out{#1.Rmd}}{\bold{#1.Rmd}}} -\section{Changes in CHNOSZ version 2.0.0-13 (2023-06-24)}{ +\section{Changes in CHNOSZ version 2.0.0-19 (2023-08-09)}{ \itemize{ @@ -28,6 +28,11 @@ \code{check.GHS()} and \code{check.EOS()} and make \code{return.difference} TRUE by default. + \item Where the changing basis species include one of the axis variables + on a diagram, \code{mosaic()} now calls \code{affinity()} with the + appropriate argument names for basis species and adjusts the labels for + the diagram (\dQuote{total C}, \dQuote{total S}, etc.). + } }