Rphree_Utils.R 7.67 KB
Newer Older
Marco De Lucia's avatar
Marco De Lucia committed
1
2
3
4
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
### Utility functions for RedModRphree

### Marco De Lucia, delucia@gfz-potsdam.de, 2009-2018
### Time-stamp: "Last modified 2018-05-03 11:03:54 delucia"

##' 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)
}