Lifting square-bracket restrictions in MxMatrix labels.
[openmx:openmx.git] / R / MxAlgebra.R
1 #
2 #   Copyright 2007-2009 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 mxAlgebra <- function(expression, name = NA, dimnames = NA) {
56         if (is.na(name)) {
57                 name <- omxUntitledName()
58         }
59         omxVerifyName(name)
60         retval <- new("MxAlgebra", NA, name)
61         retval@formula <- match.call()$expression
62         if(!(length(dimnames) == 1 && is.na(dimnames))) {
63                 dimnames(retval) <- dimnames
64         }
65         return(retval)  
66 }
67
68 defStringsAsFactors <- getOption('stringsAsFactors')
69 options(stringsAsFactors = FALSE)
70 data(omxSymbolTable)
71 options(stringsAsFactors = defStringsAsFactors)
72
73
74 formulaList <- function(x) {
75         retval <- as.list(x)
76         retval <- lapply(retval, function(x) {
77                 if(is.call(x)) {formulaList(x)} else {x}
78         })
79         return(retval)
80 }
81
82 algebraSymbolCheck <- function(formula, name) {
83         formula <- unlist(formulaList(formula))
84         test <- sapply(formula, function(x) {!is.numeric(x) && !is.character(x)})
85         if(length(formula[test]) == 1) {
86                 msg <- paste("The reference", omxQuotes(formula[test]),
87                         "is unknown in the algebra named", omxQuotes(name))
88                 stop(msg, call. = FALSE)
89         } else if (length(formula[test]) > 1) {
90                 msg <- paste("The references", omxQuotes(formula[test]),
91                         "are unknown in the algebra named", omxQuotes(name))
92                 stop(msg, call. = FALSE)                
93         }
94 }
95
96 generateAlgebraHelper <- function(algebra, matrixNames, algebraNames) {
97         retval <- algebra@formula
98         matrixNumbers <- as.list(as.integer(-1 : (-length(matrixNames))))
99         algebraNumbers <- as.list(as.integer(0 : (length(algebraNames) - 1)))
100         names(matrixNumbers) <- matrixNames
101         names(algebraNumbers) <- algebraNames
102         retval <- eval(substitute(substitute(e, matrixNumbers), list(e = retval)))
103         retval <- eval(substitute(substitute(e, algebraNumbers), list(e = retval)))
104         retval <- substituteOperators(as.list(retval))
105         algebraSymbolCheck(retval, algebra@name)
106         return(list(algebra@initial, retval))
107 }
108
109 substituteOperators <- function(algebra) {
110         if ((length(algebra) == 1) && (is.list(algebra))) {
111                 algebra <- list(0L, algebra[[1]])
112         } else if ((length(algebra) > 1) && (!is.numeric(algebra[[1]]))) {
113                 if (as.character(algebra[[1]]) == '(') {
114                         return(substituteOperators(as.list(algebra[[2]])))
115                 }
116                 names <- omxSymbolTable["R.name"] == as.character(algebra[[1]])
117         variableSymbols <- omxSymbolTable["Number.of.arguments"] == -1
118                 result <- omxSymbolTable[names & variableSymbols, "Num"]
119                 if (length(result) > 1) {
120                                 stop(paste("Ambiguous function with name", algebra[[1]],
121                                         "and", (length(algebra) - 1), "arguments"))
122                 } else if(length(result) == 1) {
123                         head <- as.integer(result[[1]])
124                         tail <- lapply(algebra[-1], substituteOperators)
125                         result <- append(tail, head, after=0)
126                         return(result)
127                 } else {
128                         length <- omxSymbolTable["Number.of.arguments"] == (length(algebra) - 1)
129                         result <- omxSymbolTable[names & length, "Num"]
130                         if (length(result) == 0) {
131                                 stop(paste("Could not find function with name", algebra[[1]],
132                                         "and", (length(algebra) - 1), "arguments"))
133                         } else if (length(result) > 1) {
134                                 stop(paste("Ambiguous function with name", algebra[[1]],
135                                         "and", (length(algebra) - 1), "arguments"))
136                         } else {
137                                 head <- as.integer(result[[1]])
138                             tail <- lapply(algebra[-1], substituteOperators)
139                                 result <- append(tail, head, after=0)
140                                 return(result)
141                         }
142         }
143         }
144         return(algebra)
145 }
146
147 # Any time mxEval() is called by one of the helper functions,
148 # we make sure to use oldFlatModel which does not have the autoboxing
149 # transformations.  Any error messages generated by oldFlatModel
150 # will not show the autoboxed matrices.
151 checkEvaluation <- function(model, flatModel, oldFlatModel) {
152         labelsData <- omxGenerateLabels(model)
153         checkMatrixEvaluation(model, oldFlatModel, labelsData)
154         flatModel <- checkAlgebraEvaluation(model, flatModel, oldFlatModel, labelsData)
155         checkConstraintEvaluation(model, oldFlatModel, labelsData)
156         return(flatModel)
157 }
158
159 checkMatrixEvaluation <- function(model, flatModel, labelsData) {
160         if(length(flatModel@matrices) == 0) { return() }
161         for(i in 1:length(flatModel@matrices)) {
162                 matrix <- flatModel@matrices[[i]]
163                 eval(computeSymbol(as.symbol(matrix@name), flatModel, labelsData))
164         }
165 }
166
167 checkAlgebraEvaluation <- function(model, retval, flatModel, labelsData) {
168         if(length(flatModel@algebras) == 0) { return(retval) }
169         for(i in 1:length(flatModel@algebras)) {
170                 algebra <- flatModel@algebras[[i]]
171                 result <- tryCatch(eval(computeSymbol(as.symbol(algebra@name), flatModel, labelsData)), 
172                         error = function(x) {
173                                 stop(paste("The algebra", 
174                                         omxQuotes(simplifyName(algebra@name, model@name)), 
175                                         "in model", omxQuotes(model@name), 
176                                         "generated the error message:",
177                                         x$message), call. = FALSE)
178                 })
179                 algebra <- retval[[algebra@name]]
180                 algebra@initial <- result
181                 retval[[algebra@name]] <- algebra
182         }
183         return(retval)
184 }
185
186 checkConstraintEvaluation <- function(model, flatModel, labelsData) {
187         if(length(flatModel@constraints) == 0) { return() }
188         for(i in 1:length(flatModel@constraints)) {
189                 constraint <- flatModel@constraints[[i]]
190                 lhs <- eval(computeSymbol(as.symbol(constraint@alg1), flatModel, labelsData))
191                 rhs <- eval(computeSymbol(as.symbol(constraint@alg2), flatModel, labelsData))
192                 if (!all(dim(lhs) == dim(rhs))) {
193                         stop(paste("The algebras/matrices", 
194                                 omxQuotes(c(simplifyName(constraint@alg1, model@name), 
195                                         simplifyName(constraint@alg2, model@name))), 
196                                 "in model", omxQuotes(model@name),
197                                 "are in constraint", omxQuotes(simplifyName(constraint@name, model@name)),
198                                 "and are not of identical dimensions"), call. = FALSE)
199                 }
200         }               
201 }
202
203 displayAlgebra <- function(mxAlgebra) {
204         cat("mxAlgebra", omxQuotes(mxAlgebra@name), '\n')
205         cat("@formula: ", deparse(mxAlgebra@formula, width.cutoff=500L), '\n')
206         if (length(mxAlgebra@result) == 0) {
207                 cat("@result: (not yet computed) ")
208         } else {
209                 cat("@result:\n")
210         }
211         print(mxAlgebra@result)
212         if (is.null(dimnames(mxAlgebra))) {
213                         cat("dimnames: NULL\n")
214         } else {
215                 cat("dimnames:\n")
216                 print(dimnames(mxAlgebra))
217         }
218 }
219
220 setMethod("print", "MxAlgebra", function(x,...) { displayAlgebra(x) })
221 setMethod("show", "MxAlgebra", function(object) { displayAlgebra(object) })