Commit e23f1366 authored by Marco De Lucia's avatar Marco De Lucia
Browse files

Merge branch 'cran' into 'master'

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

See merge request !2
parents fd398bb7 0e74fb2f
### 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
Package: RedModRphree
Title: Leverage geochemical modelling with phreeqc
Version: 0.1.1
Authors@R: c(person("Marco", "De Lucia", email = "delucia@gfz-potsdam.de", role = c("aut", "cre"), comment=c(ORCID = "0000-0002-1186-4491")),
person("Janis", "Jatnieks", email = "deltaxzz@gmail.com", role = c("ctb")))
Description: Utilities to program algorithms involving geochemical models
Depends: R (>= 3.2.0), doParallel, plyr, data.table, phreeqc, caret, mgcv, randomForest, graphics, methods
Author: Marco De Lucia [aut, cre]
Maintainer: Marco De Lucia <delucia@gfz-potsdam.de>
Title: Utilities Leveraging the R Interface to the PHREEQC Geochemical Solver
Version: 0.3.4
Authors@R: c(person(given = "Marco",
family = "De Lucia",
email = "delucia@gfz-potsdam.de",
role = c("aut", "cre"),
comment = c(ORCID = "0000-0003-0918-3766")),
person(given = "Janis",
family = "Jatnieks",
email = "deltaxzz@gmail.com",
role = c("ctb")))
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
LLazyData: true
RoxygenNote: 7.1.1
LazyData: true
This diff is collapsed.
# Generated by roxygen2: do not edit by hand
export(AME)
export(Act2pH)
export(AddProp)
export(AdvectionPQC)
export(AllowSomeNegativeColumns)
export(AnalyticalLogK)
export(BalanceEquations)
export(BatchPreProc)
export(BatchPreProc_)
export(BestModelCluster)
export(CHT)
export(CheckBalance)
export(CheckColsNumValues)
export(Cols2Fact)
export(CompareTimesteps)
export(Deltas)
export(Distribute)
export(DistributeKin)
export(DistributeKinMatrix)
export(DistributeMatrix)
export(DownSampler)
export(ElementalBalanceMin)
export(ExtractPphases)
export(ExtractSamples)
export(ExtractSpecies)
export(ExtractTotals)
export(FastClust)
export(Filtering)
export(FindAllMinNames)
export(FindAllSpeciesNames)
export(FindAllTotNames)
export(FindLogK)
export(FindPhase)
export(FitSurrogates)
export(FitWithDice)
export(FlattenList)
export(FormSelectedOutput)
export(FormulaFromBal)
export(GetModelNames)
export(GetRanges)
export(InitPreProc)
export(Loadata)
export(MAD)
export(MASE)
export(FormulaToExpression)
export(Matplot)
export(MatplotSingle)
export(ModelSelector)
export(MultiRound)
export(NPSC)
export(NumericallyCompareTables)
export(OverrideValueRange)
export(ParseFormula)
export(PlotModsInSample)
export(Pourbaix)
export(RGetPhases)
export(RPhreeExt)
export(RPhreeFile)
export(RPhreeWriteInp)
export(RPinfo)
export(RReadOut)
export(RReadOutKin)
export(RSS)
export(RangeTableCreator)
export(RReadOutVal)
export(ReactTranspBalanceEq)
export(ReactTranspBalanceKin)
export(RecomposeState)
export(ReduceState)
export(Refit)
export(Relcal)
export(RepSol)
export(Reshaping)
export(RunModelList)
export(RunPQC)
export(RunSurrogates)
export(SAD)
export(SDrift)
export(SelectActiveColumns)
export(SelectColsByPredix)
export(SelectMinActiveColumns)
export(SelectPreProc)
export(ShowTopProfiler)
export(SparseChangeCluster)
export(StoichiometricMatrix)
export(SubstractCommonTableCols)
export(SuppressSim)
export(Surrogate)
export(TByElem)
export(Tminus)
export(TrailSpaces)
export(Train)
export(TranspAsPert)
export(WriteModel)
export(WriteModelResiduals)
export(all_mda_mars_impvars)
export(allcoluniquevals)
export(cf)
export(cfa)
export(coltypes)
export(coluniquevals)
export(coluniratio)
export(colwise_RSS)
export(colwise_rescale)
export(createTuneGrids)
export(custom_scaler)
export(dt2v)
export(exc)
export(exc_type)
export(f2i)
export(hb)
export(mae)
export(mda_mars_importance)
export(mrm)
export(msg)
export(na.to.zero)
export(namewithin)
export(nan.to.zero)
export(outminus)
export(pH2Act)
export(paralap)
export(rescaling)
export(roundbycol)
export(rsi)
export(safe_get_cores)
export(seekcoltypeidx)
export(seekcoltypename)
export(sharedvalues)
export(sliding_join)
export(smartround)
export(splitMultiFix)
export(splitMultiKin)
export(start_up)
export(stopmsg)
export(striplast)
export(striplastx)
export(take_top_pct_cols)
export(to.zero)
export(uvc)
export(uvcm)
export(uvt)
export(write_evals)
import(caret)
import(data.table)
export(strrev)
import(doParallel)
import(foreach)
import(graphics)
import(methods)
import(mgcv)
import(phreeqc)
import(plyr)
import(randomForest)
## Functions for manipulating PHREEQC databases
### Marco De Lucia, delucia@gfz-potsdam.de, 2009-2018
### Time-stamp: "Last modified 2018-05-03 17:23:09 delucia"
### Marco De Lucia, delucia@gfz-potsdam.de, 2009-2021
### Time-stamp: "Last modified 2021-04-28 18:39:07 delucia"
##' @title Finds the expression of a logK of a given compound/mineral
##' in a PHREEQC database
......@@ -155,7 +155,7 @@ ParseFormula <- function(line, type="all")
return(out)
}
##' @title Find a PURE_PHASE in a database and parses its formula
##' @title Find a PURE_PHASE in a database and parse its formula
##' @param species name of phase to search for
##' @param db a PHREEQC database (buffer)
##' @return a list with the components of the formula and the
......
### File-related utility functions for RedModRphree
### Marco De Lucia, delucia@gfz-potsdam.de, 2009-2018
### Time-stamp: "Last modified 2018-05-03 11:00:16 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
......@@ -88,12 +88,12 @@ RReadOut <- function(out)
{
if (length(out)==1) {## is it a buffer or a file?
## note that empty lines are removed!
cat(paste("RReadOut:: opening file",out,"\n"))
msg("opening file",out)
tot <- RPhreeFile(out, is.db=FALSE, tabs=TRUE)
} else {
## remove empty lines and tabs as in RPhreeFile, store
## everything in "tot"
cat("RReadOut:: scanning the buffer... \n")
msg("scanning the buffer...")
tot <- sub(' +$', '', out) ## remove spaces at the end of the lines
tot <- sub("#.*$","", tot) ## remove everything after the #'s
tot <- gsub('\t', ' ', tot) ## substitute tabs with spaces
......@@ -113,8 +113,8 @@ RReadOut <- function(out)
## This is unique (only for calculated solution)
endsim <- grep('End of simulation.',tot,fixed=TRUE)
cat(paste("RReadOut::",ntot," simulation in the given output"))
msg(ntot," simulations in the given output")
## "phase assemblage" is unique, but could not be there
has_pphases <- TRUE
endkin <- grep('-Phase assemblage-',tot,fixed=TRUE)
......@@ -145,7 +145,7 @@ RReadOut <- function(out)
error <- grep("^ERROR",tot[ solutions[n] : endsim[n] ])
if (length(error) > 0) {
res[[n]] <- "error"
cat(" ERROR!!\n")
msg(" ERROR in simulation",n,", skipping")
next
}
......@@ -156,7 +156,7 @@ RReadOut <- function(out)
npure <- endpure[n] - startpure - 1
pure_conn <- textConnection( tot[startpure:(startpure+npure)])
pure <- read.table( pure_conn, row.names=1, fill=TRUE, as.is=TRUE)
pure <- utils::read.table( pure_conn, row.names=1, fill=TRUE, as.is=TRUE)
close(pure_conn)
indreactants <- which(pure[,2]=="reactant")
if (length(indreactants) > 0)
......@@ -180,7 +180,7 @@ RReadOut <- function(out)
startcomp <- endpure[n] + 2
ncomp <- endcomp[n]-startcomp
comp_conn <- textConnection(tot[startcomp:(startcomp+ncomp-1)])
comp <- read.table(comp_conn,row.names=1,fill=TRUE)
comp <- utils::read.table(comp_conn,row.names=1,fill=TRUE)
close(comp_conn)
colnames(comp) <- c("molal","moles")
## if (verbose)
......@@ -218,7 +218,7 @@ RReadOut <- function(out)
block_conn <- textConnection(block)
## Now we can read.table it
tmp <- read.table( block_conn, fill=TRUE, as.is=TRUE)[,c(1,2,3)]
tmp <- utils::read.table( block_conn, fill=TRUE, as.is=TRUE)[,c(1,2,3)]
close(block_conn)
## remove duplicated
rnames <- tmp$V1
......@@ -230,7 +230,7 @@ RReadOut <- function(out)
## Now saturation indexes
block <- tot[(endspec[n] + 2) :( endsim[n] - 4)]
block_conn <- textConnection(block)
SI <- read.table(block_conn, fill=TRUE, row.names=1, as.is=TRUE)
SI <- utils::read.table(block_conn, fill=TRUE, row.names=1, as.is=TRUE)
close(block_conn)
names(SI) <- c("SI","IAP","logK","formula")
......@@ -239,37 +239,33 @@ RReadOut <- function(out)
} ## end of loop over simulations
cat(" OK\n")
msg(" Finished, bye!")
return(res)
}
##' Kinetics simulation: reads a phreeqc output file and forms an
##' output list - as if the calculation was made through
##' Rphree.
##'
##'
##' @title RReadOutKin, import the output file of a kinetic simulation
##' into R
##' @param out The PHREEQC output file.
##' @param strip logical. If TRUE, the "ListInfo" element will be
##' appended to the list.
##' Kinetics simulation: reads a PHREEQC output file or buffer and
##' forms an output list.
##' @title RReadOutKin, imports the output file of a kinetic
##' simulation into R
##' @param out The PHREEQC output file or buffer
##' @param verbose logical. If TRUE more output is given to the
##' console (for debugging).
##' @return An output list as per Rphree call.
##' console (for debugging)
##' @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!
cat(paste("RReadOutKin:: opening file",out,"\n"))
msg("opening file",out)
tot <- RPhreeFile(out, is.db=FALSE, tabs=TRUE)
} else {
## remove empty lines and tabs as in RPhreeFile, store
## everything in "tot"
cat("RReadOutKin:: scanning the buffer... \n")
msg("scanning the buffer...")
tot <- sub(' +$', '', out) ## remove spaces at the end of the lines
tot <- sub("#.*$","", tot) ## remove everything after the #'s
tot <- gsub('\t', ' ', tot) ## substitute tabs with spaces
......@@ -281,7 +277,8 @@ RReadOutKin <- function(out, strip=TRUE, verbose=FALSE)
if (length(toremove)>0) {
toremove <- sort(c(toremove, toremove+1))
tot <- tot[-toremove]
if (verbose) cat("RReadOutKin:: Appears to be PHREEQC version 3\n")
if (verbose)
msg("Appears to be PHREEQC version 3")
}
times <- grep("Time step:", tot, fixed=TRUE)
......@@ -290,7 +287,7 @@ RReadOutKin <- function(out, strip=TRUE, verbose=FALSE)
years <- round(as.numeric(sub('\ .*$','',gsub('.*:\ ','',tot[times]))),1)
if (verbose) {
cat(out,":\n")
cat(paste("RReadOutKin:: Found ",ntot,"time steps, the last at time",years[ntot],"\n"))
msg("Found ",ntot,"time steps, the last at time",years[ntot])
}
endkin <- grep('-Phase assemblage-',tot,fixed=TRUE)
endpure <- grep('-Solution composition-',tot,fixed=TRUE)[-1]
......@@ -303,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)
......@@ -314,7 +308,7 @@ RReadOutKin <- function(out, strip=TRUE, verbose=FALSE)
## loop over all steps
for (n in seq_along(times)) {
if (verbose)
cat(paste("RReadOutKin:: Reading", n, "solution, time ",years[n],"\n"))
msg("Reading", n, "solution, time ",years[n])
## find the last solution
start <- times[n]+2
......@@ -322,35 +316,35 @@ RReadOutKin <- function(out, strip=TRUE, verbose=FALSE)
nkin <- endkin[n] - start
kconn <- textConnection(tot[(start):(start+nkin-1)])
kin <- read.table(kconn, row.names=1,fill=TRUE)[,c(2,1)]
kin <- utils::read.table(kconn, row.names=1,fill=TRUE)[,c(2,1)]
close(kconn)
colnames(kin) <- c("moles","delta")
if (verbose)
cat(paste("RReadOutKin:: Read kinetic block ", n, "\n"))
msg("Read kinetic block ", n)
## next block in output file is EQUILIBRIUM_PHASES
startpure <- endkin[n]+3
npure <- endpure[n] - startpure - 1
pconn <- textConnection(tot[startpure:(startpure+npure)])
pure <- read.table(pconn, row.names=1,fill=TRUE)[,c(5,6)]
pure <- utils::read.table(pconn, row.names=1,fill=TRUE)[,c(5,6)]
close(pconn)
colnames(pure) <- c("moles","delta")
if (verbose)
cat(paste("RReadOutKin:: Read pphases block ", n, "of length",npure+1," \n"))
msg("Read pphases block ", n, "of length",npure+1)
## now solutes
startcomp <- endpure[n] + 2
ncomp <- endcomp[n]-startcomp
sconn <- textConnection(tot[startcomp:(startcomp+ncomp-1)])
comp <- read.table(sconn, row.names=1,fill=TRUE)
comp <- utils::read.table(sconn, row.names=1,fill=TRUE)
close(sconn)
colnames(comp) <- c("molal","moles")
if (verbose)
cat(paste("RReadOutKin:: Read total solutes block ", n, "\n"))
msg("Read total solutes block ", n)
## desc: pH, pe, ecc
block <- tot[(endcomp[n]+1):(enddesc[n]-1)]
......@@ -384,7 +378,7 @@ RReadOutKin <- function(out, strip=TRUE, verbose=FALSE)
block_conn <- textConnection(block)
## Now we can read.table it
tmp <- read.table( block_conn, fill=TRUE, as.is=TRUE)[,c(1,2,3)]
tmp <- utils::read.table( block_conn, fill=TRUE, as.is=TRUE)[,c(1,2,3)]
close(block_conn)
## remove duplicated
rnames <- tmp$V1
......@@ -394,12 +388,12 @@ RReadOutKin <- function(out, strip=TRUE, verbose=FALSE)
colnames(species) <- c("molal","act")
if (verbose)
cat(paste("RReadOutKin:: Read speciation ", n, "\n"))
msg("Read speciation ", n)
## Now saturation indexes
block <- tot[(endspec[n] + 2) :( endsim[n]-1)]
block_conn <- textConnection(block)
SI <- read.table(block_conn, fill=TRUE, row.names=1, as.is=TRUE)
SI <- utils::read.table(block_conn, fill=TRUE, row.names=1, as.is=TRUE)
close(block_conn)
names(SI) <- c("SI","IAP","logK","formula")
......@@ -409,13 +403,203 @@ 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")
names(res) <- c(paste("z",1:ntot,sep=""))
attr(res,"time") <- years
return(res)
}
##' Reads a phreeqc output file or buffer and forms a results list.
##' This version looks for a precise valence of an element. It is used
##' primarily by \code{Pourbaix} to track the species and the SIs
##' which must be retained for the diagram.
##' @title RReadOutVal, parse a PHREEQC output file/buffer in search
##' 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 containing one data.frame per
##' "logical block": desc, tot, pphases, species
##' @author MDL
##' @export
RReadOutVal <- function(out, valence)
{
if (length(out)==1) {## is it a buffer or a file?
## note that empty lines are removed!
msg("opening file",out)
tot <- RPhreeFile(out, is.db=FALSE, tabs=TRUE)
} else {
## remove empty lines and tabs as in RPhreeFile, store
## everything in "tot"
msg("scanning the buffer...")
tot <- sub(' +$', '', out) ## remove spaces at the end of the lines
tot <- sub("#.*$","", tot) ## remove everything after the #'s
tot <- gsub('\t', ' ', tot) ## substitute tabs with spaces
tot <- tot[tot!=""] ## remove empty elements
}
## some extra lines in PHREEQC version 3 we need to take care of
toremove <- grep("^\\*\\*For", tot)
if (length(toremove)>0) {
toremove <- sort(c(toremove, toremove+1))
tot <- tot[-toremove]
}
solutions <- grep("Beginning of initial solution calculations", tot)
ntot <- length(solutions)
## This is unique (only for calculated solution)
endsim <- grep('End of simulation.',tot,fixed=TRUE)
msg(ntot," simulations in the buffer/file")
## "phase assemblage" is unique, but could not be there
has_pphases <- TRUE
endkin <- grep('-Phase assemblage-',tot,fixed=TRUE)
if (length(endkin)==0) {
## This means we don't have any pure phase in the simulations,
## and that the only solutions to read are the initial ones!
has_pphases <- FALSE
initials <- seq(1,ntot)
} else {
names(res) <- c(paste("z",1:ntot,sep=""))
attr(res,"time") <- years
## this index stores the lines marking initial solutions
initials <- -(seq(1,ntot)*2 -1)
}
## for these, don't consider the "initial solution" instead
endpure <- grep('-Solution composition-',tot,fixed=TRUE)[initials]
endcomp <- grep('-Description of solution-',tot,fixed=TRUE)[initials]
enddesc <- grep('-Distribution of species-',tot,fixed=TRUE)[initials]
endspec <- grep('-Saturation indices-',tot,fixed=TRUE)[initials]
## create the container
res <- vector(mode="list",length=ntot)
## loop over all steps
for (n in seq_along(solutions)) {
## check if there is a fatal error
error <- grep("^ERROR",tot[ solutions[n] : endsim[n] ])
if (length(error) > 0) {
res[[n]] <- "error"
stopmsg(" ERROR!!")
next
}
## next block in output file is the "optional"
## EQUILIBRIUM_PHASES
if (has_pphases) {
startpure <- endkin[n]+3
npure <- endpure[n] - startpure - 1
pure_conn <- textConnection( tot[startpure:(startpure+npure)])
pure <- utils::read.table( pure_conn, row.names=1, fill=TRUE, as.is=TRUE)
close(pure_conn)
indreactants <- which(pure[,2]=="reactant")
if (length(indreactants) > 0)
{
for (ir in indreactants) {
rownames(pure)[ir-1] <- paste(rownames(pure)[ir-1], rownames(pure)[ir],sep="_")
pure[(ir-1),c(4,5,6)] <- pure[ir,c(3,4,5)]
}
pure <- pure[-indreactants,]
}
pure <- pure[,c(5,6)]
colnames(pure) <- c("moles","delta")
## if (verbose)
## cat(paste(":: Read pphases block ", n, "of length",npure+1," \n"))
} else {
pure <- NA
}
## now total solutes
startcomp <- endpure[n] + 2
ncomp <- endcomp[n]-startcomp
comp_conn <- textConnection(tot[startcomp:(startcomp+ncomp-1)])
comp <- utils::read.table(comp_conn,row.names=1,fill=TRUE)
close(comp_conn)
colnames(comp) <- c("molal","moles")
## if (verbose)
## cat(paste(":: Read total solutes block ", n, "\n"))
## desc: pH, pe, ecc
block <- tot[(endcomp[n]+1):(enddesc[n]-1)]
pH <- as.numeric(unlist(strsplit(grep('^pH',block, value=TRUE)," +"))[3])
pe <- as.numeric(unlist(strsplit(grep('^pe',block, value=TRUE)," +"))[3])
temp <- as.numeric(unlist(strsplit(grep('^Temperature ',block,value=TRUE)," = "))[2])
water <- as.numeric(unlist(strsplit(grep('Mass of water',block,fixed=TRUE,value=TRUE)," = "))[2])
WaterActivity <- as.numeric(unlist(strsplit(grep('Activity of water',block,fixed=TRUE,value=TRUE)," = "))[2])
IonStr <- as.numeric(