Implemented mxEval() caching system
[openmx:openmx.git] / R / MxAlgebra.R
1 #
2 #   Copyright 2007-2012 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 setClassUnion("MxListOrNull", c("list", "NULL"))
17 setClassUnion("MxAlgebraFormula", c("call", "name", "language", "logical", "numeric"))
18
19 setClass(Class = "MxAlgebra",
20         representation = representation(
21                 formula = "MxAlgebraFormula",
22                 name = "character",
23                 dirty = "logical",
24                 .dimnames = "MxListOrNull",
25                 initial = "matrix",
26                 result = "matrix"))
27                 
28 setMethod("initialize", "MxAlgebra",
29         function(.Object, formula, name) {
30                 .Object@formula <- sys.call(which=-3)[[3]]
31                 .Object@name <- name
32                 .Object@dirty <- FALSE
33                 .Object@.dimnames <- NULL
34                 return(.Object)
35         }
36 )
37
38 setMethod("dimnames", "MxAlgebra",
39         function(x) { x@.dimnames }
40 )
41
42 setReplaceMethod("dimnames", "MxAlgebra",
43         function(x, value) {
44                 if (is.null(value)) {
45                 } else if (!is.list(value)) {
46                         stop("dimnames of MxAlgebra object must be either NULL or list of length 2")
47                 } else if (length(value) != 2) {
48                         stop("dimnames of MxAlgebra object must be either NULL or list of length 2")
49                 }
50                 x@.dimnames <- value
51                 return(x)
52         }
53 )
54
55
56
57 mxAlgebra <- function(expression, name = NA, dimnames = NA) {
58         if (single.na(name)) {
59                 name <- imxUntitledName()
60         }
61         imxVerifyName(name, 0)
62         retval <- new("MxAlgebra", NA, name)
63         retval@formula <- match.call()$expression
64     algebraErrorChecking(retval@formula, "mxAlgebra")
65         if(!(length(dimnames) == 1 && is.na(dimnames))) {
66                 dimnames(retval) <- dimnames
67         }
68         return(retval)  
69 }
70
71 algebraErrorChecking <- function(formula, context) {
72         if(length(formula) < 2) {
73                 return()
74         }
75         operator <- as.character(formula[[1]])
76         if (identical(operator, "(")) {
77         } else if (operator %in% omxSymbolTable[,'R.name']) {
78         } else {
79                 msg <- paste("Unknown matrix operator or function",
80                         omxQuotes(operator), "in",
81                         deparse(width.cutoff = 400L, imxLocateFunction(context)))
82                 stop(msg, call. = FALSE)
83         }
84         lapply(formula[-1], algebraErrorChecking)
85 }
86
87 defStringsAsFactors <- getOption('stringsAsFactors')
88 options(stringsAsFactors = FALSE)
89 data(omxSymbolTable)
90 options(stringsAsFactors = defStringsAsFactors)
91
92
93 formulaList <- function(x) {
94         retval <- as.list(x)
95         retval <- lapply(retval, function(x) {
96                 if(is.call(x)) {formulaList(x)} else {x}
97         })
98         return(retval)
99 }
100
101 algebraSymbolCheck <- function(formula, name) {
102         formula <- unlist(formulaList(formula))
103         test <- sapply(formula, function(x) {!is.numeric(x) && !is.character(x)})
104         if(length(formula[test]) == 1) {
105                 msg <- paste("The reference", omxQuotes(formula[test]),
106                         "is unknown in the algebra named", omxQuotes(name))
107                 stop(msg, call. = FALSE)
108         } else if (length(formula[test]) > 1) {
109                 msg <- paste("The references", omxQuotes(formula[test]),
110                         "are unknown in the algebra named", omxQuotes(name))
111                 stop(msg, call. = FALSE)                
112         }
113 }
114
115 generateAlgebraHelper <- function(algebra, matrixNames, algebraNames) {
116         retval <- algebra@formula
117         matrixNumbers <- as.list(as.integer(-1 : (-length(matrixNames))))
118         algebraNumbers <- as.list(as.integer(0 : (length(algebraNames) - 1)))
119         names(matrixNumbers) <- matrixNames
120         names(algebraNumbers) <- algebraNames
121         retval <- eval(substitute(substitute(e, matrixNumbers), list(e = retval)))
122         retval <- eval(substitute(substitute(e, algebraNumbers), list(e = retval)))
123         retval <- substituteOperators(as.list(retval), algebra@name)
124         algebraSymbolCheck(retval, algebra@name)
125         return(list(algebra@initial, retval))
126 }
127
128 substituteOperators <- function(algebra, name) {
129         if ((length(algebra) == 1) && (is.list(algebra))) {
130                 algebra <- list(0L, algebra[[1]])
131         } else if ((length(algebra) > 1) && (!is.numeric(algebra[[1]]))) {
132                 if (as.character(algebra[[1]]) == '(') {
133                         return(substituteOperators(as.list(algebra[[2]]), name))
134                 }
135                 names <- omxSymbolTable["R.name"] == as.character(algebra[[1]])
136         variableSymbols <- omxSymbolTable["Number.of.arguments"] == -1
137                 result <- omxSymbolTable[names & variableSymbols, "Num"]
138                 if (length(result) > 1) {
139                                 msg <- paste("In algebra", omxQuotes(name), 
140                                         "ambiguous function with name", algebra[[1]],
141                                         "and", (length(algebra) - 1), "arguments.")
142                                 stop(msg, call. = FALSE)
143                 } else if(length(result) == 1) {
144                         head <- as.integer(result[[1]])
145                         tail <- lapply(algebra[-1], substituteOperators, name)
146                         result <- append(tail, head, after=0)
147                         return(result)
148                 } else {
149                         length <- omxSymbolTable["Number.of.arguments"] == (length(algebra) - 1)
150                         result <- omxSymbolTable[names & length, "Num"]
151                         if (length(result) == 0) {
152                                 msg <- paste("In algebra", omxQuotes(name),
153                                         "could not find function with name", algebra[[1]],
154                                         "and", (length(algebra) - 1), "arguments.")
155                                 stop(msg, call. = FALSE)
156                         } else if (length(result) > 1) {
157                                 msg <- paste("In algebra", omxQuotes(name),
158                                         "ambiguous function with name", algebra[[1]],
159                                         "and", (length(algebra) - 1), "arguments.")
160                                 stop(msg, call. = FALSE)
161                         } else {
162                                 head <- as.integer(result[[1]])
163                             tail <- lapply(algebra[-1], substituteOperators, name)
164                                 result <- append(tail, head, after=0)
165                                 return(result)
166                         }
167         }
168         }
169         return(algebra)
170 }
171
172 checkEvaluation <- function(model, flatModel) {
173         labelsData <- imxGenerateLabels(model)
174         cache <- list()
175         cache <- checkMatrixEvaluation(model, flatModel, labelsData, cache)
176         tuple <- checkAlgebraEvaluation(model, flatModel, labelsData, cache)
177         flatModel <- tuple[[1]]
178         cache <- tuple[[2]]
179         checkConstraintEvaluation(model, flatModel, labelsData, cache)
180         return(flatModel)
181 }
182
183 checkMatrixEvaluation <- function(model, flatModel, labelsData, cache) {
184         if(length(flatModel@matrices) == 0) { return(cache) }
185         for(i in 1:length(flatModel@matrices)) {
186                 matrix <- flatModel@matrices[[i]]
187                 tuple <- evaluateMxObject(matrix@name, flatModel, labelsData, cache)
188                 cache <- tuple[[2]]
189         }
190         return(cache)
191 }
192
193 checkAlgebraEvaluation <- function(model, flatModel, labelsData, cache) {
194         if(length(flatModel@algebras) == 0) { return(list(flatModel, cache)) }
195         for(i in 1:length(flatModel@algebras)) {
196                 algebra <- flatModel@algebras[[i]]
197                 tuple <- evaluateMxObject(algebra@name, flatModel, labelsData, cache)
198                 result <- tuple[[1]]
199                 cache <- tuple[[2]]
200                 algebra <- flatModel[[algebra@name]]
201                 algebra@initial <- as.matrix(result)
202                 flatModel[[algebra@name]] <- algebra
203         }
204         return(list(flatModel, cache))
205 }
206
207 checkConstraintEvaluation <- function(model, flatModel, labelsData, cache) {
208         if(length(flatModel@constraints) == 0) { return() }
209         for(i in 1:length(flatModel@constraints)) {
210                 constraint <- flatModel@constraints[[i]]                
211                 lhsContext <- paste("the left-hand side of constraint", omxQuotes(constraint@name))
212                 rhsContext <- paste("the right-hand side of constraint", omxQuotes(constraint@name))            
213                 tuple <- evaluateAlgebraWithContext(flatModel[[constraint@alg1]], lhsContext, flatModel, labelsData, cache)
214                 lhs <- tuple[[1]]
215                 cache <- tuple[[2]]
216                 tuple <- evaluateAlgebraWithContext(flatModel[[constraint@alg2]], rhsContext, flatModel, labelsData, cache)
217                 rhs <- tuple[[1]]
218                 cache <- tuple[[2]]
219                 if (!all(dim(lhs) == dim(rhs))) {
220                         lhsName <- constraint@formula[[2]]
221                         rhsName <- constraint@formula[[3]]
222                         if (length(lhsName) == 1) {
223                                 lhsName <- simplifyName(deparse(lhsName), model@name)
224                         } else {
225                                 lhsName <- deparse(lhsName)
226                         }
227                         if (length(rhsName) == 1) {
228                                 rhsName <- simplifyName(deparse(rhsName), model@name)
229                         } else {
230                                 rhsName <- deparse(rhsName)
231                         }
232                         stop(paste("The algebras/matrices", 
233                                 omxQuotes(c(lhsName, rhsName)),
234                                 "in model", omxQuotes(model@name),
235                                 "are in constraint", omxQuotes(simplifyName(constraint@name, model@name)),
236                                 "and are not of identical dimensions. The left-hand side is",
237                                 nrow(lhs), "x", ncol(lhs), "and the right-hand side is",
238                                 nrow(rhs), "x", paste(ncol(rhs), ".", sep = '')), call. = FALSE)
239                 }
240         }               
241 }
242
243 displayAlgebra <- function(mxAlgebra) {
244         cat("mxAlgebra", omxQuotes(mxAlgebra@name), '\n')
245         cat("@formula: ", deparse(mxAlgebra@formula, width.cutoff=500L), '\n')
246         if (length(mxAlgebra@result) == 0) {
247                 cat("@result: (not yet computed) ")
248         } else {
249                 cat("@result:\n")
250         }
251         print(mxAlgebra@result)
252         if (is.null(dimnames(mxAlgebra))) {
253                         cat("dimnames: NULL\n")
254         } else {
255                 cat("dimnames:\n")
256                 print(dimnames(mxAlgebra))
257         }
258 }
259
260 setMethod("print", "MxAlgebra", function(x,...) { displayAlgebra(x) })
261 setMethod("show", "MxAlgebra", function(object) { displayAlgebra(object) })