Bumped to 0.3.3; fixed AddProp, removed all references to "ListInfo",...

Bumped to 0.3.3; fixed AddProp, removed all references to "ListInfo", introduced some checks in Pourbaix

.gitignore

+
+### R ###
+# History files
+.Rhistory
+.Rapp.history
+
+# Session Data files
+.RData
+
+# User-specific files
+.Ruserdata
+
+# Example code in package build process
+*-Ex.R
+
+# Output files from R CMD build
+/*.tar.gz
+
+# Output files from R CMD check
+/*.Rcheck/
+
+# RStudio files
+.Rproj.user/
+
+# produced vignettes
+vignettes/*.html
+vignettes/*.pdf
+
+# OAuth2 token, see https://github.com/hadley/httr/releases/tag/v0.3
+.httr-oauth
+
+# knitr and R markdown default cache directories
+*_cache/
+/cache/
+
+# Temporary files created by R markdown
+*.utf8.md
+*.knit.md
+
+# R Environment Variables
+.Renviron

DESCRIPTION

 Package: RedModRphree
-Title: Programmable Interface to the PHREEQC Geochemical Solver 
-Version: 0.3.2
+Title: Utilities Leveraging the R Interface to the PHREEQC Geochemical Solver 
+Version: 0.3.3
 Authors@R: c(person(given = "Marco",
                     family = "De Lucia",
                     email = "delucia@gfz-potsdam.de",
@@ -10,10 +10,11 @@ Authors@R: c(person(given = "Marco",
 	            family = "Jatnieks",
 		    email = "deltaxzz@gmail.com",
 		    role = c("ctb")))
-Description: Programmable interface to the phreeqc geochemical solver
-    adding features such as surrogate models, reactive transport and
-    Pourbaix diagrams.
-Depends: R (>= 3.2.0), doParallel, phreeqc, mgcv, graphics, methods, stats, utils, plyr, foreach
+Description: Utilities and building blocks for programming with the
+    PHREEQC geochemical solver. Includes features such as 1D reactive
+    transport with surrogate models, database manipulation, parsing of
+    PHREEQC output files and Pourbaix diagrams.
+Depends: R (>= 4.0.1), doParallel, phreeqc, mgcv, graphics, methods, stats, utils, plyr, foreach
 License: LGPL-2.1
 Encoding: UTF-8
 RoxygenNote: 7.1.1

NAMESPACE

@@ -20,6 +20,7 @@ export(FindAllSpeciesNames)
 export(FindAllTotNames)
 export(FindLogK)
 export(FindPhase)
+export(FlattenList)
 export(FormSelectedOutput)
 export(FormulaFromBal)
 export(FormulaToExpression)

R/Rphree_Files.R

 ### File-related utility functions for RedModRphree
 
-### Marco De Lucia, delucia@gfz-potsdam.de, 2009-2018
-### Time-stamp: "Last modified 2021-04-28 16:52:34 delucia"
+### Marco De Lucia, delucia@gfz-potsdam.de, 2009-2021
+### Time-stamp: "Last modified 2021-04-29 18:02:02 delucia"
 
 
 ##' Reads a normal PHREEQC input file and prepares it for
@@ -245,19 +245,18 @@ RReadOut <- function(out)
     
 }
 
-##' Kinetics simulation: reads a phreeqc output file or buffer and
+##' Kinetics simulation: reads a PHREEQC output file or buffer and
 ##' forms an output list.
-##' @title RReadOutKin, import the output file of a kinetic simulation
-##'     into R
+##' @title RReadOutKin, imports the output file of a kinetic
+##'     simulation into R
 ##' @param out The PHREEQC output file or buffer
-##' @param strip logical. If TRUE, the "ListInfo" element will be
-##'     appended to the list
 ##' @param verbose logical. If TRUE more output is given to the
 ##'     console (for debugging)
-##' @return An output list
+##' @return An output list containing one data.frame per
+##'     "logical block": desc, tot, pphases, species, kin
 ##' @author MDL
 ##' @export
-RReadOutKin <- function(out, strip=TRUE, verbose=FALSE)
+RReadOutKin <- function(out, verbose=FALSE)
 {
     if (length(out)==1) {## is it a buffer or a file?
         ## note that empty lines are removed!
@@ -301,10 +300,7 @@ RReadOutKin <- function(out, strip=TRUE, verbose=FALSE)
     endtot <- grep('End of simulation',tot,fixed=TRUE)
     endsim <- c(endsim, endtot-1)
 
-    if (strip) ## do you want "ListInfo"?
-        reslen <- ntot
-    else
-        reslen <- ntot + 1
+    reslen <- ntot
     
     ## create the container
     res <- vector(mode="list",length=reslen)
@@ -407,13 +403,9 @@ RReadOutKin <- function(out, strip=TRUE, verbose=FALSE)
         ## end loop over time steps
     }
 
-    if (!strip) {
-        res[[reslen]] <- list(n=ntot,format=TRUE,years=years)
-        names(res) <- c(paste("z",1:ntot,sep=""),"ListInfo")
-    } else {
-        names(res) <- c(paste("z",1:ntot,sep=""))
-        attr(res,"time") <- years
-    }
+    names(res) <- c(paste("z",1:ntot,sep=""))
+    attr(res,"time") <- years
+    
     return(res)
 }
 
@@ -426,9 +418,8 @@ RReadOutKin <- function(out, strip=TRUE, verbose=FALSE)
 ##'     of specific valence state of an element
 ##' @param out The PHREEQC output file or buffer
 ##' @param valence the specified element valence, e.g., Fe(2)
-##' @return An output list, as if the simulation would have being run
-##'     through Rphree (the same blocks and the same names are
-##'     returned)
+##' @return An output list containing one data.frame per
+##'     "logical block": desc, tot, pphases, species
 ##' @author MDL
 ##' @export
 RReadOutVal <- function(out, valence)

R/Rphree_Pourbaix.R

 ### Licence: LGPL version 2.1
-## Time-stamp: "Last modified 2021-04-28 17:20:11 delucia"
+## Time-stamp: "Last modified 2021-04-29 15:56:47 delucia"
 
 ##' @title Revert a string
 ##' @param x character vector. All element of the vector will be
@@ -196,16 +196,39 @@ FormSelectedOutput <- function(sim, element) {
 ##'               Patm=1, base=base, aqonly=TRUE,  db=llnl.dat)
 ##' b <- Pourbaix(element="Cu(2)", pe=seq(-10,12, length=51), pH=seq(2,10, length=51),
 ##'               Patm=1, base=base, aqonly=TRUE,  db=llnl.dat)
-Pourbaix <- function(element, pe=seq(-5,8), pH=seq(6,12), Tc=25, Patm=1, base, first, db,
-                     aqonly=FALSE, suppress, ann.title, ann.sub, colors, procs=1, plot=TRUE,
+Pourbaix <- function(element,
+                     pe=seq(-5, 8, length.out = 31),
+                     pH=seq(6, 12, length.out = 31),
+                     Tc=25,
+                     Patm=1,
+                     base,
+                     first,
+                     db,
+                     aqonly=FALSE,
+                     suppress,
+                     ann.title,
+                     ann.sub,
+                     colors,
+                     procs=1,
+                     plot=TRUE,
                      ...)
 {  
     phreeqc::phrLoadDatabaseString(db)
     phreeqc::phrSetOutputStringsOn(TRUE)
 
-    ## treat "first" as a really optional argument
+    ## "first" is an optional argument
     if (missing(first))
         first <- ""
+
+    ## some checks in "base" - pe, pH, pressure and temp must be
+    ## there, if not, add them to base script
+    for (prop in c("temp", "pe", "pH", "pressure")) {
+        if (length(grep(prop,base))==0) {
+            base <- AddProp(base, name=prop, values = 1, cat="tot")
+            msg("Adding ", prop, "to base script")
+        }
+    }
+    
     ## testrun of base script to get the output
     aa  <- phreeqc::phrRunString(c(first, base))
     con <- phreeqc::phrGetOutputStrings()
@@ -215,7 +238,7 @@ Pourbaix <- function(element, pe=seq(-5,8), pH=seq(6,12), Tc=25, Patm=1, base, f
         utils::capture.output(res <- RReadOutVal(con)[[1]])
         selout <- FormSelectedOutput(res)
     } else {
-        ## first we distinguish if is total amount of element or if it is a given valence
+        ## first we distinguish if we want total amount of element or a specific valence
         if (length(grep("(", element, fixed=TRUE))>0) {
             ## now we call "valence" the full expression (including the parentheses)
             valence <- element
@@ -232,6 +255,7 @@ Pourbaix <- function(element, pe=seq(-5,8), pH=seq(6,12), Tc=25, Patm=1, base, f
     }
 
     ## do we only want aqueous species?
+    ## if yes, just suppress SIs from selected output
     if (aqonly) {
         selout <- selout[-grep("-saturation_indices", selout, fixed=TRUE)]
     }
@@ -284,9 +308,9 @@ Pourbaix <- function(element, pe=seq(-5,8), pH=seq(6,12), Tc=25, Patm=1, base, f
         }
 
         doParallel::registerDoParallel(ThisRunCluster)
-        msg("Registered default doParallel cluster with ", procs, "nodes\n")
+        msg("Registered default doParallel cluster with ", procs, "nodes")
         parallel::clusterCall(cl=ThisRunCluster, phreeqc::phrLoadDatabaseString, db)
-        msg("Database loaded on each worker\n")
+        msg("Database loaded on each worker")
 
         items <- nrow(fcombs)
         if (items %/% procs > 0)
@@ -472,8 +496,14 @@ Pourbaix <- function(element, pe=seq(-5,8), pH=seq(6,12), Tc=25, Patm=1, base, f
         box()
     } ## plotting
     
-    invisible(list(mat=mat, displayed=todisplay, formula=parsednames, is.aq=is.aqueous,
+    invisible(list(mat=mat,
+                   displayed=todisplay,
+                   formula=parsednames,
+                   is.aq=is.aqueous,
                    phreeqc=list(input=c(first, selout, biginp), pqc=bigres, db=db, res=res),
-                   dat=dat, vec=vec, aq=aq, eq=eq))
+                   dat=dat,
+                   vec=vec,
+                   aq=aq,
+                   eq=eq))
 }
     

R/Rphree_SurrogateUtils.R

 ##  Functions for dealing with surrogate simulations
 
 ### Marco De Lucia, delucia@gfz-potsdam.de, 2009-2021
-### Time-stamp: "Last modified 2021-04-28 16:09:22 delucia"
+### Time-stamp: "Last modified 2021-04-29 18:05:53 delucia"
 
 ##' Computes the average of absolute values of a vector
 ##' @title Average of absolute values
@@ -11,41 +11,24 @@
 ##' @export
 mae <- function(x) mean(abs(x))
 
-## ##' .. content for \description{} (no empty lines) ..
-## ##'
-## ##' .. content for \details{} ..
-## ##' @title SoftMax transformation
-## ##' @param obj numeric vector
-## ##' @return 
-## ##' @author 
-## SoftMax <- function(obj) {
-##     eo <- exp(obj-max(obj))
-##     ss <- sum(eo)
-##     scaled <- eo/ss
-##     attr(res,"back") <- ss
-##     return(res)
-## }
-
-## ## SoftMaxBackTransf <- function(vec, back) return(log(vec)+log(back))
-
-## ##' @title Flatten out a Rphree list
-## ##' @param list list of Rphree lists
-## ##' @param strip logical, should we delete "ListInfo"? Defaults to
-## ##'     true
-## ##' @return a list whose elements are single Rphree simulations
-## ##' @author MDL
-## ##' @export
-## FlattenList <- function(list) {
-##    tot <- vector(mode="list", length=sum(sapply(stripped, length)))
-##    k <- 1
-##    for (i in seq_along(stripped)) {
-##       for (j in seq_along(stripped[[i]])) {
-##          tot[[k]] <- stripped[[i]][[j]]
-##          k=k+1
-##       }
-##    }
-##    return(tot)
-## }
+
+##' This function is useful when parsing many output files in parallel.
+##' @title Flatten a list of Rphree lists
+##' @param list list of Rphree lists
+##' @return a list whose elements are single Rphree simulations
+##' @author MDL
+##' @export
+FlattenList <- function(list) {
+   tot <- vector(mode="list", length=sum(sapply(list, length)))
+   k <- 1
+   for (i in seq_along(list)) {
+      for (j in seq_along(list[[i]])) {
+         tot[[k]] <- list[[i]][[j]]
+         k=k+1
+      }
+   }
+   return(tot)
+}
 
 
 ##' @title Find all species occurring in the ensemble
@@ -85,7 +68,7 @@ ExtractSpecies <- function(flatlist, species) {
    tmp <- matrix(NA, len <- length(flatlist), speclen <- length(species))
    colnames(tmp) <- species
    for (i in seq(speclen)) {
-      tmp[,i] <- RPinfo(flatlist, "species", species[i], flex=TRUE)
+      tmp[,i] <- RPinfo(flatlist, "species", species[i])
    }
    return(tmp)
 }
@@ -101,7 +84,7 @@ ExtractPphases <- function(flatlist, pphases) {
    tmp <- matrix(NA, len <- length(flatlist), speclen <- length(pphases))
    colnames(tmp) <- pphases
    for (i in seq(speclen)) {
-      tmp[,i] <- RPinfo(flatlist, "pphase", pphases[i], flex=TRUE)
+      tmp[,i] <- RPinfo(flatlist, "pphase", pphases[i])
    }
    return(tmp)
 }
@@ -119,7 +102,7 @@ ExtractTotals <- function(flatlist, totals) {
    tmp <- matrix(NA, len <- length(flatlist), speclen <- length(totals))
    colnames(tmp) <- totals
    for (i in seq(speclen)) {
-      tmp[,i] <- RPinfo(flatlist, "tot", totals[i], flex=TRUE)
+      tmp[,i] <- RPinfo(flatlist, "tot", totals[i])
    }
    return(tmp)
 }

R/Rphree_Utils.R

 ### Utility functions for RedModRphree
 
 ### Marco De Lucia, delucia@gfz-potsdam.de, 2009-2021
-### Time-stamp: "Last modified 2021-04-28 17:06:47 delucia"
+### Time-stamp: "Last modified 2021-04-29 17:36:29 delucia"
 
 ##' Replicates an input buffer containing only one SOLUTION, taking
 ##' care of SOLUTION/KINETICS/PURE identifiers. Eventually insert a
@@ -31,7 +31,7 @@ RepSol <- function(sol, n, first=NULL)
     
     if (length(linesol <- grep("^SOLUTION",newsol)) != n)
         {
-            stop("too many or no SOLUTION defined\n")
+            stopmsg("too many or no SOLUTION defined")
         } else {
             newsol[linesol] <- paste("SOLUTION",1:n)
         }
@@ -41,7 +41,7 @@ RepSol <- function(sol, n, first=NULL)
         {
             if (length(linepure) != n)
                 {
-                    stop("too many or no PURE defined\n")
+                    stopmsg("too many or no PURE defined")
                 } else {
                     newsol[linepure] <- paste("PURE",1:n)
                 }
@@ -52,7 +52,7 @@ RepSol <- function(sol, n, first=NULL)
         {
             if (length(linekin) != n)
                 {
-                    stop("too many or no KINET defined\n")
+                    stopmsg("too many or no KINET defined")
                 } else {
                     newsol[linekin] <- paste("KINETICS",1:n)
                 }
@@ -85,56 +85,67 @@ AddProp <- function(input, name, values, cat, kinpar=NULL, first=NULL)
     reg <- gsub("\\)","\\\\)",reg)
     
     ## coerce values to character and find lengths
-    n <- length(val <- as.character(values))
-    nsim <- length(grep("^SOLUTION",input))
-
-    if (n != nsim && nsim != 1)
-        stop("AddProp :: Lengths of simulation and values to add do not agree!\n")
-    
-    if (cat == "kin") {
-        if (nsim != 1)
-            stop(":: AddProp with KINETICS is able to handle only one simul at a time!")
-
-        markkin <- max(grep("^KINETIC|step",input))
+    val <- as.character(values)
+    nval <- length(val)
+    nsim <- length(marksol <- grep("^SOLUTION",input))
 
-        ## add the lines
-        tmp <- c(input[1:markkin],name,paste("-m0", val), paste("-parms", kinpar), input[(markkin+1):length(input)])
-        if (n!=1)
-            ## use Distribute to create the new input
-            newinp <- Distribute(tmp, reg, val, first)
-        else
-            newinp <- tmp
-    } else { ## from here pphases and tot
-        ## for now, use "PURE" as delimiter 
-        markpure <- grep("PURE",input)
-
-        ## one line before if is a "tot" component, one line after if is
-        ## a "pphases"
-        if (cat == "tot") 
-            markpure <- markpure - 1
+    if (nval != nsim && nsim != 1)
+        stopmsg("Lengths of simulation and values to add do not agree!")
 
+    if (cat == "tot") {
         if (nsim == 1) {
-            ## if input is just 1 simulation, then first add the line and then
-            ## repeat it n times (adding "first")
-            ## add the line
-            tmp <- c(input[1:markpure],paste(name,val),input[(markpure+1):length(input)])
-            if (n!=1)
+            ## if input is just 1 simulation, then just insert the line after "SOLUTION"
+            tmp <- c(input[1:marksol],paste(name,val),input[(marksol+1):length(input)])
+            
+            if (nval>1)
                 ## use Distribute to create the new input
                 newinp <- Distribute(tmp, reg, val, first)
             else
                 newinp <- tmp
         } else {
-            ## beautiful indexes arithmetic :)
-            markpure <- markpure+seq_along(markpure) - ifelse(cat=="tot",1,0)
-            tmp <- character(length(input)+n)
-            tmp[markpure] <- paste(reg,val)
-            rest <- seq_along(tmp)[-markpure]
+            marksol <- marksol + seq_along(marksol) ## the line after "SOLUTION"
+            tmp <- character(length(input)+nval)
+            tmp[marksol] <- paste(reg,val)
+            rest <- seq_along(tmp)[-marksol]
             tmp[rest] <- input
             ## use Distribute to create the new input
             newinp <- Distribute(tmp, reg, values, first)      
         }
     }
-    return(newinp)  
+
+    if (cat == "pphases") {
+        markpure <- grep("PURE", input)
+        if (nsim == 1) {
+            tmp <- c(input[1:markpure], paste(name, val[1]), input[(markpure + 1):length(input)])
+            if (nval != 1) 
+                newinp <- Distribute(tmp, reg, val, first)
+            else newinp <- tmp
+        } else {
+            markpure <- markpure + seq_along(markpure)
+            tmp <- character(length(input) + nval)
+            tmp[markpure] <- paste(reg, val)
+            rest <- seq_along(tmp)[-markpure]
+            tmp[rest] <- input
+            newinp <- Distribute(tmp, reg, values, first)
+        }
+    }    
+
+    if (cat == "kin") {
+        if (nsim != 1)
+            stopmsg(" with KINETICS is able to handle only one simul at a time!")
+
+        markkin <- max(grep("^KINETIC|step",input))
+
+        ## add the lines
+        tmp <- c(input[1:markkin],name,paste("-m0", val), paste("-parms", kinpar), input[(markkin+1):length(input)])
+        if (nval != 1)
+            ## use Distribute to create the new input
+            newinp <- Distribute(tmp, reg, val, first)
+        else
+            newinp <- tmp
+    }
+
+    return(newinp)
 }
 
 ##' Function to distribute different values of one property
@@ -170,25 +181,24 @@ Distribute <- function(input, prop, values, newname=NULL, first=NULL, wholeline=
 
     ## coerce values to character
     n <- length(val <- as.character(values))
-    nsim <- length(proplines <- grep(paste("^", reg," ", sep=""),input))
+    nsim <- length(proplines <- grep(paste0("^", reg),input))
     
     if (n != nsim && nsim != 1) {
-        ##  cat("Lengths of simulation and values do not agree, n=",n,", nsim=",nsim,"!\n")
-        stop("Distribute :: Lengths of simulation and values do not agree")
+        stopmsg("Lengths of simulation and values do not agree")
     }
     linesave <- NULL
     
     if (nsim == 1) {
         ## deal with comments after the properties (i.e. "as HCO3")
         if (!wholeline) {
-            dum <- unlist(strsplit(gsub(' +',' ',grep(paste("^", reg," ",sep=""),input,value=TRUE))," "))
+            dum <- unlist(strsplit(gsub(' +',' ',grep(paste0("^", reg), input, value=TRUE))," "))
             if (length(dum)>2)
                 linesave <- paste(dum[3:length(dum)],collapse=" ")
         }
         ## if the simulation is just 1, then repeat it n times (adding "first")
         ## grep(paste("^", reg, " ", sep=""), input, value=TRUE)
         newinp <- RepSol(input, n, first)
-        nsim <- length(proplines <- grep(paste("^", reg, " ", sep=""), newinp))
+        nsim <- length(proplines <- grep(paste0("^", reg), newinp))
     } else {
         newinp <- input
     }
@@ -229,9 +239,6 @@ SuppressSim <- function(biginp, n=1L)
 ##'     valid numeric value (0) is returned even if the property is
 ##'     not found value is returned if the inquired property is not
 ##'     found in the solution
-##' @param flex Logical. If TRUE, expects no "ListInfo" in the
-##'     formatted solution list and performs heuristics to circumvent
-##'     this absence
 ##' @return A numeric vector containing the inquired properties
 ##' @author MDL
 ##' @examples
@@ -241,47 +248,35 @@ SuppressSim <- function(biginp, n=1L)
 ##' }
 ##' @export
 RPinfo <- function(lin, cat=c("tot","desc","pphases","master","species","kin","SI"),
-                   prop, force=TRUE, flex=FALSE)
+                   prop, force=TRUE)
 {
     cat <- match.arg(cat)
-    
-    if ("ListInfo" %in% names(lin))
-        nsim <- lin$ListInfo$n
-    else {
-        if (!flex)
-            stop("RPinfo:: no ListInfo! Perhaps not a valid list produced by Rphree?
-  Specify flex=TRUE if you still want to try")
-        else
-            nsim <- length(lin)
-    }
-    
-    if ( nsim > 1 )
+    nsim <- length(lin)
+       
+    if ( nsim > 1 ) {
         ## we have a list of solutions, check if the name exists
-        {
-            if ( cat %in% names(lin[[1]])) {
-                ## everything seems fine, let's (s)apply
-                lout <- sapply(lin[1:nsim], RPhreeExt, cat=cat, prop=prop)
+        if ( cat %in% names(lin[[1]])) {
+            ## everything seems fine, let's (s)apply
+            lout <- sapply(lin[1:nsim], RPhreeExt, cat=cat, prop=prop)
+        } else {
+            if (force) {
+                lout <- rep(0,length(lin))
             } else {
-                if (force)
-                    {
-                        lout <- rep(0,length(lin))
-                    } else {
-                        stop("RPinfo:: Sure that the right output was selected in the simulation?")
-                    }
+                stopmsg("Sure that the correct output was selected in the simulation?")
             }
+        }
+    } else {
+        ## just 1 solution
+        if ( cat %in% names(lin)) {
+            lout <- RPhreeExt(lin, cat, prop)
         } else {
-            ## just 1 solution
-            if ( cat %in% names(lin)) {
-                lout <- RPhreeExt(lin, cat, prop)
+            if (force) {
+                lout <- 0
             } else {
-                if (force)
-                    {
-                        lout <- 0
-                    } else {
-                        stop("RPinfo:: Sure that the right output was selected in the simulation?")
-                    }
-            } 
-        }
+                stopmsg("Sure that the right output was selected in the simulation?")
+            }
+        } 
+    }
     return(lout)
 }
 

demo/00Index

@@ -3,4 +3,5 @@ demo-kinetics         1D reactive transport with kinetic chemistry, full-physics
 demo-eq-generate-data Generate the training data set from 3 equilibrium simulations
 demo-eq-surr-RF       Generate the training data set for equilibrium simulations, train randomForest surrogates and perform Reactive Transport simulations using them
 demo-kin-surr-RF      Generate the training data set for kinetic simulations, train randomForest surrogates and perform Reactive Transport simulations using them
-demo-gmd2020-445      Pointer to the actual gitlab repository for the code
+demo-Pourbaix         Example of Pourbaix diagram
+demo-gmd2020-445      Just a pointer to the actual gitlab repo for the code

demo/demo-Pourbaix.R

+### RedModRphree, demo for Pourbaix
+
+### Licence: LGPL version 2.1
+### Time-stamp: "Last modified 2021-04-29 15:47:44 delucia"
+library(RedModRphree)
+
+base <- c("SOLUTION 1", 
+          "units mol/kgw", 
+          "temp 25",
+          "pressure 1",
+          "pH 7", 
+          "pe 4.0", 
+          "Ca 0.25",
+