Rphree_Utils.R 12.3 KB
Newer Older
Marco De Lucia's avatar
Marco De Lucia committed
1
2
### Utility functions for RedModRphree

3
### Marco De Lucia, delucia@gfz-potsdam.de, 2009-2021
4
### Time-stamp: "Last modified 2021-04-28 17:06:47 delucia"
Marco De Lucia's avatar
Marco De Lucia committed
5
6
7

##' Replicates an input buffer containing only one SOLUTION, taking
##' care of SOLUTION/KINETICS/PURE identifiers. Eventually insert a
8
##' block (e.g: PUNCH) just once under the first SOLUTION.
Marco De Lucia's avatar
Marco De Lucia committed
9
10
11
##'
##' Analogous to standard \code{rep} function.
##' @title Replicate a SOLUTION 
12
13
##' @param sol The initial phreeqc buffer that needs to be replicated
##' @param n Number of replicates to generate
Marco De Lucia's avatar
Marco De Lucia committed
14
##' @param first A new block which needs to be inserted only in the
15
16
##' first solution (typically, PUNCH and alike)
##' @return A buffer upon which RunPQC can be called
Marco De Lucia's avatar
Marco De Lucia committed
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
##' @author MDL
##' @export
RepSol <- function(sol, n, first=NULL)
{
    if (length(dum <- grep("^TITLE|^DATABASE",sol)) > 0 )
        sol <- sol[-dum]
    
    if (is.character(first))
        {
            firstend <- grep("^END",sol)[1]
            newsol <- c(sol[-firstend],first,"END",rep(sol,n-1))
        } else {
            newsol <- rep(sol,n)
        }
    
    if (length(linesol <- grep("^SOLUTION",newsol)) != n)
        {
            stop("too many or no SOLUTION defined\n")
        } else {
            newsol[linesol] <- paste("SOLUTION",1:n)
        }
    
    linepure <- grep("^PURE|^EQUIL",newsol)
    if (length(linepure) > 0 )
        {
            if (length(linepure) != n)
                {
                    stop("too many or no PURE defined\n")
                } else {
                    newsol[linepure] <- paste("PURE",1:n)
                }
        }
    
    linekin <- grep("^KINET",newsol)
    if (length(linekin) > 0)
        {
            if (length(linekin) != n)
                {
                    stop("too many or no KINET defined\n")
                } else {
                    newsol[linekin] <- paste("KINETICS",1:n)
                }
        }
    return(newsol)  
}

##' Function that adds a new property (species, pure phase, ...) to an
##' input buffer.
##' @title Add a new property (species, pure phase, ...) to an input
##' buffer
##' @param input The input buffer.
##' @param name Name of the new property to add (element, PURE
##' mineral,...)
##' @param values The value(s) of the property.
##' @param cat The category of the property. Must be one of
##' \code{c("tot","pphases","kin")}
73
##' @param kinpar Kinetics parameters for the KINETICS case
Marco De Lucia's avatar
Marco De Lucia committed
74
##' @param first A block to be added in the first SOLUTION, as in
75
76
##' \code{RepSol} and \code{Distribute}
##' @return The new input buffer upon which RunPQC can be called
Marco De Lucia's avatar
Marco De Lucia committed
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
##' @author MDL
##' @export
AddProp <- function(input, name, values, cat, kinpar=NULL, first=NULL)
{
    cat <- match.arg(cat,c("tot","pphases","kin"))
    
    ## correct quoting of parenthesis
    reg <- gsub("\\(","\\\\(",name)
    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))

        ## 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 (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)
                ## 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]
            tmp[rest] <- input
            ## use Distribute to create the new input
            newinp <- Distribute(tmp, reg, values, first)      
        }
    }
    return(newinp)  
}

##' Function to distribute different values of one property
##' (concentration, temperature, mineral phase,\dots) properties in an
##' input buffer.
##'
##' The number of SOLUTIONs in the input buffer must be 1 or equal to
##' length(values). If the input buffer contains only one SOLUTION,
##' and length(values)>1, then the new input gets automatically
##' replicated length(values) times by RepSol. For KINETICS use
##' \code{DistributeKin}.
##' 
150
151
152
153
154
##' @title Distribute properties in an input buffer
##' @param input The initial input buffer
##' @param prop The property whose values need to be distributed
##' @param values The numerical value(s) to distribute across the
##'     SOLUTIONS
Marco De Lucia's avatar
Marco De Lucia committed
155
##' @param newname Optional name of a property which is not initially
156
##'     present in the input buffer
Marco De Lucia's avatar
Marco De Lucia committed
157
##' @param first The \code{first} block which will be passed to
158
159
160
161
##'     \code{\link{RepSol}} if needed
##' @param wholeline logical. if TRUE, comments after the properties
##'     (i.e. "as HCO3") get also distributed.
##' @return A new input buffer upon which RunPQC can be run
Marco De Lucia's avatar
Marco De Lucia committed
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
##' @author MDL
##' @export
Distribute <- function(input, prop, values, newname=NULL, first=NULL, wholeline=TRUE)
{

    ## correct quoting of parenthesis
    reg <- gsub("\\(","\\\\(",prop)
    reg <- gsub("\\)","\\\\)",reg)

    ## coerce values to character
    n <- length(val <- as.character(values))
    nsim <- length(proplines <- grep(paste("^", reg," ", sep=""),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")
    }
    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))," "))
            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))
    } else {
        newinp <- input
    }
    if (is.null(newname))
        newname <- prop

    newinp[proplines] <- paste(newname, val, linesave)
    return(newinp)  
}


##' Erases the SOLUTION number n (label, not position) from an input
##' buffer containing many.
##' @title Suppress a simulation from an input buffer
##' @param biginp The input buffer
##' @param n The number of the simulation to erase
##' @return The new input buffer 
##' @author MDL
##' @export
SuppressSim <- function(biginp, n=1L)
{
    startsols <- grep("^SOLUTION",biginp)
    newinp <- biginp[ -seq(startsols[n], startsols[n+1]-1) ]
    return(newinp)
}
217
218


219
220
221
222
##' This functions extracts informations extracted from standard
##' PHREEQC output files
##' @title Extract specific values from a phreeqc output
##' @param lin The list containg the parsed phreeqc solution(s)
223
##' @param cat A string indicating the category (or output block) to
224
225
##'     search in, must be one of:
##'     "tot, desc, pphases, master, species, kin, SI"
226
##' @param prop The name of the inquired property (element, species,
227
228
229
230
231
##'     mineral phase,\dots) whose value(s)
##' @param force Logical. If TRUE (default if left unspecified), a
##'     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
232
##' @param flex Logical. If TRUE, expects no "ListInfo" in the
233
234
##'     formatted solution list and performs heuristics to circumvent
##'     this absence
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
##' @return A numeric vector containing the inquired properties
##' @author MDL
##' @examples
##' \dontrun{
##'   RPinfo(SOL,"tot","Fe")
##'   RPinfo(SOL,"pphases","Anhydrite")
##' }
##' @export
RPinfo <- function(lin, cat=c("tot","desc","pphases","master","species","kin","SI"),
                   prop, force=TRUE, flex=FALSE)
{
    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 )
        ## 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)
            } else {
                if (force)
                    {
                        lout <- rep(0,length(lin))
                    } else {
                        stop("RPinfo:: Sure that the right output was selected in the simulation?")
                    }
            }
        } else {
            ## just 1 solution
            if ( cat %in% names(lin)) {
                lout <- RPhreeExt(lin, cat, prop)
            } else {
                if (force)
                    {
                        lout <- 0
                    } else {
                        stop("RPinfo:: Sure that the right output was selected in the simulation?")
                    }
            } 
        }
    return(lout)
}

##' Workhorse function for extraction of results from a solution,
##' primarily intended to be used by \code{\link{RPinfo}}, where it is
##' called inside a \code{sapply} statement.
##' @title Workhorse function for extraction of informations from a
292
##'     solution
293
294
##' @param sol The solution outputted by Rphree
##' @param cat String: the category or output block where to look for
295
##'     the property \code{prop}
296
297
298
299
300
301
302
303
304
305
306
307
308
##' @param prop The name of the property to look for
##' @return Numeric of length 1, the inquired value for the property
##' @author MDL
##' @export
RPhreeExt <- function(sol, cat, prop)
{
    if (is.data.frame(sol[[cat]]))
        return(sol[[cat]][prop,1])
    else
        return(0)
}

##' Workhorse function for extraction of results from a solution,
309
##' specifically equilibrium phases.
310
##'
311
312
##' @title Workhorse function for extraction of equilibrium phases
##'     from a solution list
313
314
##' @param lin The solution outputted by Rphree
##' @param which String containing the name of the specific pphase to
315
##'     look for. If omitted, all pphases are returned.
316
##' @param delta Logical. If TRUE, delta min are returned, if FALSE,
317
##'     total amount
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
##' @return A data.frame containing the results of equilibrium phases
##' @author MDL
##' @export
RGetPhases <- function(lin, which=NULL, delta=TRUE)
{
    i <- 1
    if (delta) i <- 2
    
    tab <- lin$pphases
    if (!is.character(which)) {
        return(tab[,i])
    } else {
        inde <- match(which,rownames(tab))
        return(tab[inde,i])
    }
}


336
337
338
339
340
341
342
343
344
345
346
347
348
## ##' Workhorse function for extraction of results from a solution,
## ##' specifically equilibrium phases. 
## ##'
## ##' @title extract the phases from a PHREEQC output
## ##' @param lin The solution outputted by phree
## ##' @return A data.frame containing the results of equilibrium phases
## ##' @author MDL
## ##' @export
## RReadPhases <- function(lin)
## {
##     ns <- length( phas <- unique(unlist(lapply(lin,function (x) return(rownames(x$pphases))))) )
##     if (ns==0)
##         return(NULL)
349
    
350
351
352
353
354
355
356
357
358
##     mat <- matrix(0,ncol=ns,nrow=n)
##     for (i in seq_along(phas))
##         {
##             mat[,i] <- RPinfo(lin,"pphases",phas[i])
##         }
##     colnames(mat) <- phas
##     mat[is.na(mat)] <- 0
##     return(mat)
## }