Rphree_Utils.R 11.8 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-29 17:36:29 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
##' @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)
        {
34
            stopmsg("too many or no SOLUTION defined")
Marco De Lucia's avatar
Marco De Lucia committed
35
36
37
38
39
40
41
42
43
        } else {
            newsol[linesol] <- paste("SOLUTION",1:n)
        }
    
    linepure <- grep("^PURE|^EQUIL",newsol)
    if (length(linepure) > 0 )
        {
            if (length(linepure) != n)
                {
44
                    stopmsg("too many or no PURE defined")
Marco De Lucia's avatar
Marco De Lucia committed
45
46
47
48
49
50
51
52
53
54
                } else {
                    newsol[linepure] <- paste("PURE",1:n)
                }
        }
    
    linekin <- grep("^KINET",newsol)
    if (length(linekin) > 0)
        {
            if (length(linekin) != n)
                {
55
                    stopmsg("too many or no KINET defined")
Marco De Lucia's avatar
Marco De Lucia committed
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
                } 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
##' @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
88
89
90
    val <- as.character(values)
    nval <- length(val)
    nsim <- length(marksol <- grep("^SOLUTION",input))
Marco De Lucia's avatar
Marco De Lucia committed
91

92
93
    if (nval != nsim && nsim != 1)
        stopmsg("Lengths of simulation and values to add do not agree!")
Marco De Lucia's avatar
Marco De Lucia committed
94

95
    if (cat == "tot") {
Marco De Lucia's avatar
Marco De Lucia committed
96
        if (nsim == 1) {
97
98
99
100
            ## 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)
Marco De Lucia's avatar
Marco De Lucia committed
101
102
103
104
105
                ## use Distribute to create the new input
                newinp <- Distribute(tmp, reg, val, first)
            else
                newinp <- tmp
        } else {
106
107
108
109
            marksol <- marksol + seq_along(marksol) ## the line after "SOLUTION"
            tmp <- character(length(input)+nval)
            tmp[marksol] <- paste(reg,val)
            rest <- seq_along(tmp)[-marksol]
Marco De Lucia's avatar
Marco De Lucia committed
110
111
112
113
114
            tmp[rest] <- input
            ## use Distribute to create the new input
            newinp <- Distribute(tmp, reg, values, first)      
        }
    }
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

    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)
Marco De Lucia's avatar
Marco De Lucia committed
149
150
151
152
153
154
155
156
157
158
159
160
}

##' 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}.
##' 
161
162
163
164
165
##' @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
166
##' @param newname Optional name of a property which is not initially
167
##'     present in the input buffer
Marco De Lucia's avatar
Marco De Lucia committed
168
##' @param first The \code{first} block which will be passed to
169
170
171
172
##'     \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
173
174
175
176
177
178
179
180
181
182
183
##' @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))
184
    nsim <- length(proplines <- grep(paste0("^", reg),input))
Marco De Lucia's avatar
Marco De Lucia committed
185
186
    
    if (n != nsim && nsim != 1) {
187
        stopmsg("Lengths of simulation and values do not agree")
Marco De Lucia's avatar
Marco De Lucia committed
188
189
190
191
192
193
    }
    linesave <- NULL
    
    if (nsim == 1) {
        ## deal with comments after the properties (i.e. "as HCO3")
        if (!wholeline) {
194
            dum <- unlist(strsplit(gsub(' +',' ',grep(paste0("^", reg), input, value=TRUE))," "))
Marco De Lucia's avatar
Marco De Lucia committed
195
196
197
198
199
200
            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)
201
        nsim <- length(proplines <- grep(paste0("^", reg), newinp))
Marco De Lucia's avatar
Marco De Lucia committed
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
    } 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)
}
227
228


229
230
231
232
##' 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)
233
##' @param cat A string indicating the category (or output block) to
234
235
##'     search in, must be one of:
##'     "tot, desc, pphases, master, species, kin, SI"
236
##' @param prop The name of the inquired property (element, species,
237
238
239
240
241
##'     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
242
243
244
245
246
247
248
249
250
##' @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"),
251
                   prop, force=TRUE)
252
253
{
    cat <- match.arg(cat)
254
255
256
    nsim <- length(lin)
       
    if ( nsim > 1 ) {
257
        ## we have a list of solutions, check if the name exists
258
259
260
261
262
263
        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))
264
            } else {
265
                stopmsg("Sure that the correct output was selected in the simulation?")
266
            }
267
268
269
270
271
        }
    } else {
        ## just 1 solution
        if ( cat %in% names(lin)) {
            lout <- RPhreeExt(lin, cat, prop)
272
        } else {
273
274
            if (force) {
                lout <- 0
275
            } else {
276
277
278
279
                stopmsg("Sure that the right output was selected in the simulation?")
            }
        } 
    }
280
281
282
283
284
285
286
    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
287
##'     solution
288
289
##' @param sol The solution outputted by Rphree
##' @param cat String: the category or output block where to look for
290
##'     the property \code{prop}
291
292
293
294
295
296
297
298
299
300
301
302
303
##' @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,
304
##' specifically equilibrium phases.
305
##'
306
307
##' @title Workhorse function for extraction of equilibrium phases
##'     from a solution list
308
309
##' @param lin The solution outputted by Rphree
##' @param which String containing the name of the specific pphase to
310
##'     look for. If omitted, all pphases are returned.
311
##' @param delta Logical. If TRUE, delta min are returned, if FALSE,
312
##'     total amount
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
##' @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])
    }
}


331
332
333
334
335
336
337
338
339
340
341
342
343
## ##' 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)
344
    
345
346
347
348
349
350
351
352
353
##     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)
## }