New mxConstraint() interface implementation
[openmx:openmx.git] / R / MxAlgebra.R
1 #\r
2 #   Copyright 2007-2009 The OpenMx Project\r
3 #\r
4 #   Licensed under the Apache License, Version 2.0 (the "License");\r
5 #   you may not use this file except in compliance with the License.\r
6 #   You may obtain a copy of the License at\r
7\r
8 #        http://www.apache.org/licenses/LICENSE-2.0\r
9\r
10 #   Unless required by applicable law or agreed to in writing, software\r
11 #   distributed under the License is distributed on an "AS IS" BASIS,\r
12 #   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.\r
13 #   See the License for the specific language governing permissions and\r
14 #   limitations under the License.\r
15 \r
16 setClassUnion("MxListOrNull", c("list", "NULL"))\r
17 setClassUnion("MxAlgebraFormula", c("call", "name", "language", "logical", "numeric"))\r
18 \r
19 setClass(Class = "MxAlgebra",\r
20         representation = representation(\r
21                 formula = "MxAlgebraFormula",\r
22                 name = "character",\r
23                 dirty = "logical",\r
24                 .dimnames = "MxListOrNull",\r
25                 initial = "matrix",\r
26                 result = "matrix"))\r
27                 \r
28 setMethod("initialize", "MxAlgebra",\r
29         function(.Object, formula, name) {\r
30                 .Object@formula <- sys.call(which=-3)[[3]]\r
31                 .Object@name <- name\r
32                 .Object@dirty <- FALSE\r
33                 .Object@.dimnames <- NULL\r
34                 return(.Object)\r
35         }\r
36 )\r
37 \r
38 setMethod("dimnames", "MxAlgebra",\r
39         function(x) { x@.dimnames }\r
40 )\r
41 \r
42 setReplaceMethod("dimnames", "MxAlgebra",\r
43         function(x, value) {\r
44                 if (is.null(value)) {\r
45                 } else if (!is.list(value)) {\r
46                         stop("dimnames of MxAlgebra object must be either NULL or list of length 2")\r
47                 } else if (length(value) != 2) {\r
48                         stop("dimnames of MxAlgebra object must be either NULL or list of length 2")\r
49                 }\r
50                 x@.dimnames <- value\r
51                 return(x)\r
52         }\r
53 )\r
54 \r
55 mxAlgebra <- function(expression, name = NA, dimnames = NA) {\r
56         if (is.na(name)) {\r
57                 name <- omxUntitledName()\r
58         }\r
59         omxVerifyName(name)\r
60         retval <- new("MxAlgebra", NA, name)\r
61         retval@formula <- match.call()$expression\r
62         if(!(length(dimnames) == 1 && is.na(dimnames))) {\r
63                 dimnames(retval) <- dimnames\r
64         }\r
65         return(retval)  \r
66 }\r
67 \r
68 defStringsAsFactors <- getOption('stringsAsFactors')\r
69 options(stringsAsFactors = FALSE)\r
70 data(omxSymbolTable)\r
71 options(stringsAsFactors = defStringsAsFactors)\r
72 \r
73 \r
74 formulaList <- function(x) {\r
75         retval <- as.list(x)\r
76         retval <- lapply(retval, function(x) {\r
77                 if(is.call(x)) {formulaList(x)} else {x}\r
78         })\r
79         return(retval)\r
80 }\r
81 \r
82 algebraSymbolCheck <- function(formula, name) {\r
83         formula <- unlist(formulaList(formula))\r
84         test <- sapply(formula, function(x) {!is.numeric(x) && !is.character(x)})\r
85         if(length(formula[test]) == 1) {\r
86                 msg <- paste("The reference", omxQuotes(formula[test]),\r
87                         "is unknown in the algebra named", omxQuotes(name))\r
88                 stop(msg, call. = FALSE)\r
89         } else if (length(formula[test]) > 1) {\r
90                 msg <- paste("The references", omxQuotes(formula[test]),\r
91                         "are unknown in the algebra named", omxQuotes(name))\r
92                 stop(msg, call. = FALSE)                \r
93         }\r
94 }\r
95 \r
96 generateAlgebraHelper <- function(algebra, matrixNames, algebraNames) {\r
97         retval <- algebra@formula\r
98         matrixNumbers <- as.list(as.integer(-1 : (-length(matrixNames))))\r
99         algebraNumbers <- as.list(as.integer(0 : (length(algebraNames) - 1)))\r
100         names(matrixNumbers) <- matrixNames\r
101         names(algebraNumbers) <- algebraNames\r
102         retval <- eval(substitute(substitute(e, matrixNumbers), list(e = retval)))\r
103         retval <- eval(substitute(substitute(e, algebraNumbers), list(e = retval)))\r
104         retval <- substituteOperators(as.list(retval))\r
105         algebraSymbolCheck(retval, algebra@name)\r
106         return(list(algebra@initial, retval))\r
107 }\r
108 \r
109 substituteOperators <- function(algebra) {\r
110         if ((length(algebra) == 1) && (is.list(algebra))) {\r
111                 algebra <- list(0L, algebra[[1]])\r
112         } else if ((length(algebra) > 1) && (!is.numeric(algebra[[1]]))) {\r
113                 if (as.character(algebra[[1]]) == '(') {\r
114                         return(substituteOperators(as.list(algebra[[2]])))\r
115                 }\r
116                 names <- omxSymbolTable["R.name"] == as.character(algebra[[1]])\r
117         variableSymbols <- omxSymbolTable["Number.of.arguments"] == -1\r
118                 result <- omxSymbolTable[names & variableSymbols, "Num"]\r
119                 if (length(result) > 1) {\r
120                                 stop(paste("Ambiguous function with name", algebra[[1]],\r
121                                         "and", (length(algebra) - 1), "arguments"))\r
122                 } else if(length(result) == 1) {\r
123                         head <- as.integer(result[[1]])\r
124                         tail <- lapply(algebra[-1], substituteOperators)\r
125                         result <- append(tail, head, after=0)\r
126                         return(result)\r
127                 } else {\r
128                         length <- omxSymbolTable["Number.of.arguments"] == (length(algebra) - 1)\r
129                         result <- omxSymbolTable[names & length, "Num"]\r
130                         if (length(result) == 0) {\r
131                                 stop(paste("Could not find function with name", algebra[[1]],\r
132                                         "and", (length(algebra) - 1), "arguments"))\r
133                         } else if (length(result) > 1) {\r
134                                 stop(paste("Ambiguous function with name", algebra[[1]],\r
135                                         "and", (length(algebra) - 1), "arguments"))\r
136                         } else {\r
137                                 head <- as.integer(result[[1]])\r
138                             tail <- lapply(algebra[-1], substituteOperators)\r
139                                 result <- append(tail, head, after=0)\r
140                                 return(result)\r
141                         }\r
142         }\r
143         }\r
144         return(algebra)\r
145 }\r
146 \r
147 # Any time mxEval() is called by one of the helper functions,\r
148 # we make sure to use oldFlatModel which does not have the autoboxing\r
149 # transformations.  Any error messages generated by oldFlatModel\r
150 # will not show the autoboxed matrices.\r
151 checkEvaluation <- function(model, flatModel, oldFlatModel) {\r
152         labelsData <- omxGenerateLabels(model)\r
153         checkMatrixEvaluation(model, oldFlatModel, labelsData)\r
154         flatModel <- checkAlgebraEvaluation(model, flatModel, oldFlatModel, labelsData)\r
155         checkConstraintEvaluation(model, flatModel, labelsData)\r
156         return(flatModel)\r
157 }\r
158 \r
159 checkMatrixEvaluation <- function(model, flatModel, labelsData) {\r
160         if(length(flatModel@matrices) == 0) { return() }\r
161         for(i in 1:length(flatModel@matrices)) {\r
162                 matrix <- flatModel@matrices[[i]]\r
163                 eval(computeSymbol(as.symbol(matrix@name), flatModel, labelsData))\r
164         }\r
165 }\r
166 \r
167 checkAlgebraEvaluation <- function(model, retval, flatModel, labelsData) {\r
168         if(length(flatModel@algebras) == 0) { return(retval) }\r
169         for(i in 1:length(flatModel@algebras)) {\r
170                 algebra <- flatModel@algebras[[i]]\r
171                 result <- tryCatch(eval(computeSymbol(as.symbol(algebra@name), flatModel, labelsData)), \r
172                         error = function(x) {\r
173                                 stop(paste("The algebra", \r
174                                         omxQuotes(simplifyName(algebra@name, model@name)), \r
175                                         "in model", omxQuotes(model@name), \r
176                                         "generated the error message:",\r
177                                         x$message), call. = FALSE)\r
178                 })\r
179                 algebra <- retval[[algebra@name]]\r
180                 algebra@initial <- result\r
181                 retval[[algebra@name]] <- algebra\r
182         }\r
183         return(retval)\r
184 }\r
185 \r
186 checkConstraintEvaluation <- function(model, flatModel, labelsData) {\r
187         if(length(flatModel@constraints) == 0) { return() }\r
188         for(i in 1:length(flatModel@constraints)) {\r
189                 constraint <- flatModel@constraints[[i]]
190                 lhs <- tryCatch(eval(computeSymbol(as.symbol(constraint@alg1), flatModel, labelsData)), \r
191                         error = function(x) {\r
192                                 stop(paste("The left hand side of constraint", \r
193                                         omxQuotes(simplifyName(constraint@name, model@name)), \r
194                                         "in model", omxQuotes(model@name), \r
195                                         "generated the error message:",\r
196                                         x$message), call. = FALSE)\r
197                 })
198                 rhs <- tryCatch(eval(computeSymbol(as.symbol(constraint@alg2), flatModel, labelsData)), \r
199                         error = function(x) {\r
200                                 stop(paste("The right hand side of constraint", \r
201                                         omxQuotes(simplifyName(constraint@name, model@name)), \r
202                                         "in model", omxQuotes(model@name), \r
203                                         "generated the error message:",\r
204                                         x$message), call. = FALSE)\r
205                 })\r
206                 if (!all(dim(lhs) == dim(rhs))) {
207                         lhsName <- constraint@formula[[2]]
208                         rhsName <- constraint@formula[[3]]
209                         if (length(lhsName) == 1) {
210                                 lhsName <- simplifyName(deparse(lhsName), model@name)
211                         } else {
212                                 lhsName <- deparse(lhsName)
213                         }
214                         if (length(rhsName) == 1) {
215                                 rhsName <- simplifyName(deparse(rhsName), model@name)
216                         } else {
217                                 rhsName <- deparse(rhsName)
218                         }\r
219                         stop(paste("The algebras/matrices", \r
220                                 omxQuotes(c(lhsName, rhsName)),
221                                 "in model", omxQuotes(model@name),\r
222                                 "are in constraint", omxQuotes(simplifyName(constraint@name, model@name)),\r
223                                 "and are not of identical dimensions"), call. = FALSE)\r
224                 }\r
225         }               \r
226 }\r
227 \r
228 displayAlgebra <- function(mxAlgebra) {\r
229         cat("mxAlgebra", omxQuotes(mxAlgebra@name), '\n')\r
230         cat("@formula: ", deparse(mxAlgebra@formula, width.cutoff=500L), '\n')\r
231         if (length(mxAlgebra@result) == 0) {\r
232                 cat("@result: (not yet computed) ")\r
233         } else {\r
234                 cat("@result:\n")\r
235         }\r
236         print(mxAlgebra@result)\r
237         if (is.null(dimnames(mxAlgebra))) {\r
238                         cat("dimnames: NULL\n")\r
239         } else {\r
240                 cat("dimnames:\n")\r
241                 print(dimnames(mxAlgebra))\r
242         }\r
243 }\r
244 \r
245 setMethod("print", "MxAlgebra", function(x,...) { displayAlgebra(x) })\r
246 setMethod("show", "MxAlgebra", function(object) { displayAlgebra(object) })