Starting changes to LISREL to type='LISREL' and to expectations for mxEval(model...
[openmx:openmx.git] / R / MxEval.R
1 #
2 #   Copyright 2007-2013 The OpenMx Project
3 #
4 #   Licensed under the Apache License, Version 2.0 (the "License");
5 #   you may not use this file except in compliance with the License.
6 #   You may obtain a copy of the License at
7
8 #        http://www.apache.org/licenses/LICENSE-2.0
9
10 #   Unless required by applicable law or agreed to in writing, software
11 #   distributed under the License is distributed on an "AS IS" BASIS,
12 #   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13 #   See the License for the specific language governing permissions and
14 #   limitations under the License.
15
16 mxEval <- function(expression, model, compute = FALSE, show = FALSE, defvar.row = 1,
17                         cache = new.env(parent = emptyenv()), cacheBack = FALSE) {
18         if (missing(expression)) {
19                 stop("'expression' argument is mandatory in call to mxEval function")
20         } else if (missing(model)) {
21                 stop("'model' argument is mandatory in call to mxEval function")
22         } else if (!is.numeric(defvar.row) || length(defvar.row) != 1 || is.na(defvar.row)) {
23                 stop("'defvar.row' argument must be a single numeric value")
24         }       
25         expression <- match.call()$expression
26         modelvariable <- match.call()$model
27         labelsData <- imxGenerateLabels(model)
28         env <- parent.frame()
29         if (compute) {
30                 expression <- formulaEliminateObjectiveFunctions(expression)
31                 namespace <- imxGenerateNamespace(model)
32                 model <- imxFlattenModel(model, namespace)
33                 expression <- qualifyNamesFormula(expression, model@name, namespace)
34         }
35         if (show) {
36                 tuple <- evaluateExpression(expression, deparse(expression), model, 
37                         labelsData, env, compute, show = TRUE, outsideAlgebra = TRUE,
38                         cache = new.env(parent = emptyenv()))
39                 showResult <- tuple[[1]]
40                 showResult <- eval(substitute(substitute(x, list(.zzz = modelvariable)), list(x = showResult)))
41                 print(deparse(showResult), width.cutoff = 450L)
42         }
43         tuple <- evaluateExpression(expression, deparse(expression), model, 
44                 labelsData, env, compute, show = FALSE, outsideAlgebra = TRUE, defvar.row,
45                 cache = cache)
46         result <- tuple[[1]]
47         cache <- tuple[[2]]
48         if (!is.vector(result)) {
49                 result <- eval(result, envir = env)
50         }
51         if (cacheBack) {
52                 return(list(result, cache))
53         } else {
54                 return(result)
55         }
56 }
57
58 translateErrorSymbol <- function(symbol, model) {
59         charsymbol <- as.character(symbol)
60         object <- model[[charsymbol]]
61         if (!is.null(object) && is(object, "MxMatrix") && length(object@display) == 1) {
62                 return(as.symbol(object@display))
63         } else {
64                 return(symbol)
65         }
66 }
67
68 translateErrorFormula <- function(formula, model) {
69         if(length(formula) == 1) {
70                 formula <- translateErrorSymbol(formula, model)
71         } else {
72                 formula <- lapply(formula, translateErrorFormula, model)
73                 formula <- as.call(formula)
74         }
75         return(formula)
76 }
77
78 evaluateExpression <- function(formula, contextString, model, labelsData, 
79         env, compute, show, outsideAlgebra, defvar.row = 1, cache = new.env(parent = emptyenv())) {
80         len <- length(formula)
81         if (len == 0) {
82                 stop("mxEval has reached an invalid state")
83         } else if (len == 1) {
84                 if (!identical(as.character(formula), "")) {
85                         tuple <- evaluateSymbol(formula, contextString, model, 
86                                 labelsData, env, compute, show, outsideAlgebra, defvar.row,
87                                 cache)
88                         return(tuple)
89                 } else {
90                         return(list(formula, cache))
91                 }
92         }
93         originalFormula <- formula
94         for(i in 2:length(formula)) {
95                 tuple <- evaluateExpression(formula[[i]], contextString, 
96                         model, labelsData, env, compute, show, 
97                         outsideAlgebra, defvar.row, cache)
98                 formula[[i]] <- tuple[[1]]
99                 cache <- tuple[[2]]
100         }
101         if (len == 4 && identical(as.character(formula[1]), '[')) {
102                 formula$drop <- FALSE
103                 if(is.matrix(formula[[3]]) && nrow(formula[[3]]) == 0 && ncol(formula[[3]]) == 0) {
104                         formula[[3]] <- TRUE
105                 }
106                 if(is.matrix(formula[[4]]) && nrow(formula[[4]]) == 0 && ncol(formula[[4]]) == 0) {
107                         formula[[4]] <- TRUE
108                 }
109         }
110         operator <- as.character(formula[1])
111         if (len == 3 && operator %in% c('*', '^', '/')) {
112                 formula[[2]] <- substitute(as.matrix(x), list(x = formula[[2]]))
113                 formula[[3]] <- substitute(as.matrix(x), list(x = formula[[3]]))
114         }
115         if (show) {
116                 return(list(result, cache))
117         }
118         result <- tryCatch(do.call(as.character(formula[[1]]), as.list(formula[-1]), envir = env),
119                         error = function(x) {
120                                 originalFormula <- translateErrorFormula(originalFormula, model)
121                                 formulaString <- deparse(originalFormula, width.cutoff = 500L)
122                                 if(identical(formulaString, contextString)) {
123                                         msg <- paste("The following error occurred while",
124                                                 "evaluating the expression", omxQuotes(formulaString), 
125                                                 "in model", omxQuotes(model@name),
126                                                 ":", x$message)
127                                 } else {
128                                         msg <- paste("The following error occurred while",
129                                                 "evaluating the subexpression",
130                                                 omxQuotes(formulaString), "during the evaluation of",
131                                                 omxQuotes(contextString), "in model", omxQuotes(model@name),
132                                                 ":", x$message)
133                                 }
134                                 stop(msg, call. = FALSE)
135                 })
136         return(list(result, cache))
137 }
138
139 evaluateSymbol <- function(symbol, contextString, model, labelsData,
140                         env, compute, show, outsideAlgebra, defvar.row = 1, 
141                         cache) {
142         key <- deparse(symbol)
143         if (exists(key, envir = cache)) {
144                 result <- get(key, envir = cache)
145         } else {
146                 index <- match(key, dimnames(labelsData)[[1]])
147                 if (!is.na(index)) {
148                         targetmodel <- labelsData[index, "model"]
149                         matrix <- labelsData[index, "matrix"]
150                         fullname <- imxIdentifier(targetmodel, matrix)
151                         row <- labelsData[index, "row"]
152                         col <- labelsData[index, "col"]
153                         value <- model[[fullname]]@values[row,col]
154                         if (show) {
155                                 result <- substitute(.zzz[[x]]@values[y,z], list(x = fullname, y = row, z = col))
156                         } else {
157                                 result <- value
158                         }
159                 } else {
160                         lookup <- model[[key]]
161                         if (is.list(lookup) && length(lookup) == 2 && is(lookup[[2]], "MxFitFunction")) {
162                                 lookup <- lookup[[2]]
163                         }
164                         if (is.list(lookup) && length(lookup) == 2 && is(lookup[[2]], "MxExpectation")) {
165                                 lookup <- lookup[[2]]
166                         }
167                         if (imxIsDefinitionVariable(key)) {
168                                 result <- definitionStartingValue(key, contextString, model, defvar.row)
169                         } else if (is.null(lookup)) {
170                                 if (!show && !outsideAlgebra && exists(key, envir = env)) {
171                                         result <- as.matrix(get(key, envir = env))
172                                 } else {
173                                         result <- symbol
174                                 }
175                         } else if (is(lookup, "MxMatrix")) {
176                                 tuple <- evaluateMatrix(lookup, model, labelsData, env, show, compute, defvar.row, cache)
177                                 result <- tuple[[1]]
178                                 cache <- tuple[[2]]
179                         } else if (is(lookup, "MxAlgebra")) {
180                                 tuple <- evaluateAlgebra(lookup, model, labelsData, env, show, compute, defvar.row, cache)
181                                 result <- tuple[[1]]
182                                 cache <- tuple[[2]]
183                         } else if (is(lookup, "MxData")) {
184                                 if (show) {
185                                         result <- substitute(.zzz[[x]]@observed, list(x = key))
186                                 } else {
187                                         result <- lookup@observed
188                                 }
189                         } else if (is(lookup, "MxFitFunction")) {
190                                 if (length(lookup@result) == 0) {
191                                         if (compute) {
192                                                 result <- genericFitInitialMatrix(lookup, model)
193                                         } else {
194                                                 result <- matrix(0, 0, 0)
195                                         }
196                                 } else if (show) {
197                                         result <- substitute(.zzz[[x]]@result, list(x = key))
198                                 } else {
199                                         result <- lookup@result
200                                 }
201                         } else if (is(lookup, "MxExpectation")) {
202                                 stop("mxEval for Expectations is not yet implemented.")
203                                 #if (length(lookup@result) == 0) {
204                                 #       if (compute) {
205                                 #               result <- genericFitInitialMatrix(lookup, model)
206                                 #       } else {
207                                 #               result <- matrix(0, 0, 0)
208                                 #       }
209                                 #} else if (show) {
210                                 #       result <- substitute(.zzz[[x]]@result, list(x = key))
211                                 #} else {
212                                 #       result <- lookup@result
213                                 #}
214                         } else {
215                                 stop(paste("Cannot evaluate the object",
216                                         omxQuotes(key), "in the model",
217                                         omxQuotes(model@name)))
218                         }
219                 }
220                 assign(key, result, envir = cache)
221         }
222         return(list(result, cache))
223 }
224
225 evaluateMatrix <- function(matrix, model, labelsData, env, show, compute, defvar.row, cache) {
226         if (show) {
227                 return(list(substitute(.zzz[[x]]@values, list(x = matrix@name)), cache))
228         } else if (compute) {
229                 tuple <- computeMatrix(matrix, model, labelsData, defvar.row, env, cache)
230         } else {
231                 tuple <- list(matrix@values, cache)
232         }
233         tuple[[1]] <- assignDimnames(matrix, tuple[[1]])
234         return(tuple)
235 }
236
237 computeMatrix <- function(matrix, model, labelsData, defvar.row, env, cache) {
238         matrix <- populateDefVarMatrix(matrix, model, defvar.row)
239         values <- matrix@values
240         labels <- matrix@labels
241         select <- matrix@.squareBrackets
242         if (all(!select)) {
243                 return(list(values, cache))
244         }
245         subs <- labels[select]
246         rows <- row(labels)[select]
247         cols <- col(labels)[select]
248         for (i in 1:length(subs)) {
249                 substitution <- subs[[i]]
250                 row <- rows[[i]]
251                 col <- cols[[i]]
252                 components <- splitSubstitution(substitution)
253                 expression <- parse(text = paste('`', components[[1]], '`[',
254                         as.integer(components[[2]]), ',', as.integer(components[[3]]), ']', sep=''))[[1]]
255                 contextString <- paste("label at row ", row, " and column ", col, " of matrix ", 
256                         omxQuotes(simplifyName(matrix@name, model@name)), sep = '')
257                 tuple <- evaluateExpression(expression, contextString, model, labelsData, env, 
258                         compute=TRUE, show=FALSE, outsideAlgebra=FALSE, cache=cache)
259                 result <- tuple[[1]]
260                 cache <- tuple[[2]]
261                 result <- as.matrix(result)
262                 if (nrow(result) != 1 || ncol(result) != 1) {
263                         stop(paste("The label", 
264                                 omxQuotes(substitution),
265                                 "of matrix", omxQuotes(simplifyName(matrix@name, model@name)),
266                                 "in model", omxQuotes(model@name), 
267                                 "does not evaluate to a (1 x 1) matrix."), call. = FALSE)
268                 }
269                 values[row, col] <- result[1,1]
270         }
271         return(list(values, cache))
272 }
273
274 evaluateAlgebra <- function(algebra, model, labelsData, env, show, compute, defvar.row, cache) {
275         if (compute && show) {
276                 return(computeAlgebra(algebra, model, labelsData, show = TRUE, defvar.row, env, cache))
277         } else if (compute) {
278                 tuple <- computeAlgebra(algebra, model, labelsData, show = FALSE, defvar.row, env, cache)
279                 tuple[[1]] <- as.matrix(tuple[[1]])
280         } else if (show) {
281                 result <- substitute(.zzz[[x]]@result, list(x = algebra@name))
282                 return(list(result, cache))
283         } else {
284                 tuple <- list(algebra@result, cache)
285         }
286         if (!is.null(dimnames(algebra))) {
287                 tuple[[1]] <- assignDimnames(algebra, tuple[[1]])
288         }
289         return(tuple)
290 }
291
292 computeAlgebra <- function(algebra, model, labelsData, show, defvar.row, env, cache) {
293         contextString <- simplifyName(algebra@name, model@name)
294         tuple <- evaluateExpression(algebra@formula, contextString, model, labelsData, env, 
295                 compute = TRUE, show, outsideAlgebra = FALSE, defvar.row, cache)
296         return(tuple)
297 }
298
299 evaluateMxObject <- function(objname, flatModel, labelsData, cache) {
300         return(eval(substitute(evaluateSymbol(x, objname, flatModel, 
301                         labelsData, globalenv(), compute = TRUE, 
302                         show = FALSE, outsideAlgebra = FALSE, defvar.row = 1, 
303                         cache = cache),
304                         list(x = quote(as.symbol(objname))))))
305 }
306
307 evaluateAlgebraWithContext <- function(algebra, context, flatModel, labelsData, cache) {
308         return(evaluateExpression(algebra@formula, context, flatModel, labelsData, globalenv(), 
309                 compute = TRUE, show = FALSE, outsideAlgebra = FALSE, cache = cache))
310 }
311
312 getEntityType <- function(object) {
313         if(is(object, "MxMatrix")) {
314                 entity <- "MxMatrix"
315         } else if(is(object, "MxAlgebra")) {
316                 entity <- "MxAlgebra"
317         } else if(is(object, "MxExpectation")) {
318                 entity <- "MxExpectation"
319         } else if(is(object, "MxFitFunction")) {
320                 entity <- "MxFitFunction"
321         } else {
322                 entity <- "entity"
323         }
324         return(entity)
325 }
326
327 assignDimnames <- function(object, values) {
328         newnames <- dimnames(object)
329         if (!is.null(newnames)) {
330                 if (length(newnames) != 2) {
331                         entity <- getEntityType(object)
332                         msg <- paste("The 'dimnames' argument to", entity,
333                                 omxQuotes(object@name), "must have a length of 2")
334                         stop(msg, call. = FALSE)
335                 }
336                 if (!is.null(newnames[[1]]) && length(newnames[[1]]) != nrow(values) || 
337                         !is.null(newnames[[2]]) && length(newnames[[2]]) != ncol(values)) {
338                         entity <- getEntityType(object)
339                         msg <- paste("The", entity, omxQuotes(object@name), 
340                                 "has specified dimnames with dimensions",
341                                 length(newnames[[1]]), "x", length(newnames[[2]]), "but the", entity,
342                                         "is of dimensions", nrow(values), "x", ncol(values))
343                         stop(msg, call. = FALSE)
344                 }
345         }
346         dimnames(values) <- dimnames(object)
347         return(values)
348 }
349
350 imxGenerateLabels <- function(model) {
351         return(generateLabelsHelper(model, data.frame()))
352 }
353
354 generateLabelsHelper <- function(model, labelsData) {
355         if (length(model@matrices) > 0) {
356                 for(i in 1:length(model@matrices)) {
357                         labelsData <- generateLabelsMatrix(model@name, model@matrices[[i]], labelsData)
358                 }
359         }
360         if (length(model@submodels) > 0) {
361                 for(i in 1:length(model@submodels)) {
362                         labelsData <- generateLabelsHelper(model@submodels[[i]], labelsData)
363                 }
364         }
365         return(labelsData)
366 }
367
368 retainLabel <- function(label, existing) {
369     return(!imxIsDefinitionVariable(label) && 
370                         !hasSquareBrackets(label) &&
371                         !(label %in% existing))
372 }
373
374 generateLabelsMatrix <- function(modelName, matrix, labelsData) {
375         labels <- matrix@labels
376         test <- !is.na(labels)
377         select <- labels[test]
378         rows <- row(labels)[test]
379         cols <- col(labels)[test]
380         if (length(select) > 0) {
381                 retain <- !duplicated(select) & sapply(select, 
382                         retainLabel, rownames(labelsData))
383                 select <- select[retain]
384                 rows <- rows[retain]
385                 cols <- cols[retain]
386                 if (length(select) > 0) {
387                         newLabels <- data.frame(model=modelName, 
388                                 matrix=matrix@name, row=rows, 
389                                 col=cols, stringsAsFactors=FALSE)
390                         rownames(newLabels) <- select
391                         labelsData <- rbind(labelsData, newLabels)
392                 }
393         }
394         return(labelsData)
395 }