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

### Marco De Lucia, delucia@gfz-potsdam.de, 2009-2018
4
### Time-stamp: "Last modified 2021-04-12 12:58:36 delucia"
Marco De Lucia's avatar
Marco De Lucia committed
5
6
7
8
9
10
11
12
13
14
15
16
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
73
74
75
76
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
150
151
152
153
154
155
156
157
158
159
160
161
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
217
218

##' Replicates an input buffer containing only one SOLUTION, taking
##' care of SOLUTION/KINETICS/PURE identifiers. Eventually insert a
##' block (e.g: PUNCH) under the first SOLUTION.
##'
##' Analogous to standard \code{rep} function.
##' @title Replicate a SOLUTION 
##' @param sol The initial phreeqc buffer that needs to be replicated.
##' @param n Number of replicates to generate.
##' @param first A new block which needs to be inserted only in the
##' first solution (typically, PUNCH and alike).
##' @return A buffer upon which Rphree can be called.
##' @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.
##'
##' .. content for DETAILS ..
##' @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")}
##' @param kinpar Kinetics parameters for the KINETICS case.
##' @param first A block to be added in the first SOLUTION, as in
##' \code{RepSol} and \code{Distribute}.
##' @return The new input buffer upon which Rphree can be called.
##' @author MDL
##' @export
AddProp <- function(input, name, values, cat, kinpar=NULL, first=NULL)
### Add a property (species, pure phase, ...) to an input buffer
{
    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}.
##' 
##' @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.
##' @param newname Optional name of a property which is not initially
##' present in the input buffer.
##' @param first The \code{first} block which will be passed to
##' \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 Rphree can be run.
##' @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)
}
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
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
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360


##' This functions extracts informations from calculated Rphree
##' formatted solutions.
##'
##' @title Extract specific values from an Rphree output
##' @param lin The list containg the Rphree solution(s)
##' @param cat A string indicating the category (or output block) to
##' search in, must be one of: "tot, desc, pphases, master, species,
##' kin, SI"
##' @param prop The name of the inquired property (element, species,
##' 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
##' @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
##' \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
##' solution
##' @param sol The solution outputted by Rphree
##' @param cat String: the category or output block where to look for
##' the property \code{prop}
##' @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,
##' specifically equilibrium phases. 
##'
##' @title Workhorse function for extraction of 
##' @param lin The solution outputted by Rphree
##' @param which String containing the name of the specific pphase to
##' look for. If omitted, all pphases are returned.
##' @param delta Logical. If TRUE, delta min are returned, if FALSE,
##' total amount
##' @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])
    }
}


##' Workhorse function for extraction of results from a solution,
##' specifically equilibrium phases. Variant.
##'
##' @title extract the phases from a Rphree solution
##' @param lin The solution outputted by Rphree
##' @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)
    
    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)
}