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