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