mosaic() can now change basis species that are axis variables on a di…

…agram

git-svn-id: svn://scm.r-forge.r-project.org/svnroot/chnosz/pkg/CHNOSZ@799 edb9625f-4e0d-4859-8d74-9fd3b1da38cb

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", , "[email protected]", role = c("aut", "cre"),

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, ...)

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

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)
 }

inst/NEWS.Rd

@@ -12,7 +12,7 @@
 % links to vignettes 20220723
 \newcommand{\viglink}{\ifelse{html}{\out{<a href="../CHNOSZ/doc/#1.html"><strong>#1.Rmd</strong></a>}}{\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.).
+
     }
 
 }