Added cycle detection to algebra expressions.
[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", "logical"))
18
19 setClass(Class = "MxAlgebra",
20         representation = representation(
21                 formula = "MxAlgebraFormula",
22                 name = "character",
23                 dirty = "logical",
24                 .dimnames = "MxListOrNull",
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 mxAlgebra <- function(expression, name = NA, dimnames = NA) {
55         if (is.na(name)) {
56                 name <- omxUntitledName()
57         }
58         omxVerifyName(name)
59         retval <- new("MxAlgebra", NA, name)
60         retval@formula <- match.call()$expression
61         if(!(length(dimnames) == 1 && is.na(dimnames))) {
62                 dimnames(retval) <- dimnames
63         }
64         return(retval)  
65 }
66
67 defStringsAsFactors <- getOption('stringsAsFactors')
68 options(stringsAsFactors = FALSE)
69 data(omxSymbolTable)
70 options(stringsAsFactors = defStringsAsFactors)
71
72
73 formulaList <- function(x) {
74         retval <- as.list(x)
75         retval <- lapply(retval, function(x) {
76                 if(is.call(x)) {formulaList(x)} else {x}
77         })
78         return(retval)
79 }
80
81 algebraSymbolCheck <- function(formula, name) {
82         formula <- unlist(formulaList(formula))
83         test <- sapply(formula, function(x) {!is.numeric(x)})
84         if(length(formula[test]) == 1) {
85                 msg <- paste("The reference", omxQuotes(formula[test]),
86                         "is unknown in the algebra named", omxQuotes(name))
87                 stop(msg, call. = FALSE)
88         } else if (length(formula[test]) > 1) {
89                 msg <- paste("The references", omxQuotes(formula[test]),
90                         "are unknown in the algebra named", omxQuotes(name))
91                 stop(msg, call. = FALSE)                
92         }
93 }
94
95 generateAlgebraHelper <- function(algebra, matrixNames, algebraNames) {
96         retval <- algebra@formula
97         matrixNumbers <- as.list(as.integer(-1 : (-length(matrixNames))))
98         algebraNumbers <- as.list(as.integer(0 : (length(algebraNames) - 1)))
99         names(matrixNumbers) <- matrixNames
100         names(algebraNumbers) <- algebraNames
101         retval <- eval(substitute(substitute(e, matrixNumbers), list(e = retval)))
102         retval <- eval(substitute(substitute(e, algebraNumbers), list(e = retval)))
103         retval <- substituteOperators(as.list(retval))
104         algebraSymbolCheck(retval, algebra@name)
105         return(retval)
106 }
107
108 substituteOperators <- function(algebra) {
109         if ((length(algebra) == 1) && (is.list(algebra))) {
110                 algebra <- list(0L, algebra[[1]])
111         } else if ((length(algebra) > 1) && (!is.numeric(algebra[[1]]))) {
112                 if (as.character(algebra[[1]]) == '(') {
113                         return(substituteOperators(as.list(algebra[[2]])))
114                 }
115                 names <- omxSymbolTable["R.name"] == as.character(algebra[[1]])
116         variableSymbols <- omxSymbolTable["Number.of.arguments"] == -1
117                 result <- omxSymbolTable[names & variableSymbols, "Num"]
118                 if (length(result) > 1) {
119                                 stop(paste("Ambiguous function with name", algebra[[1]],
120                                         "and", (length(algebra) - 1), "arguments"))
121                 } else if(length(result) == 1) {
122                         head <- as.integer(result[[1]])
123                         tail <- lapply(algebra[-1], substituteOperators)
124                         result <- append(tail, head, after=0)
125                         return(result)
126                 } else {
127                         length <- omxSymbolTable["Number.of.arguments"] == (length(algebra) - 1)
128                         result <- omxSymbolTable[names & length, "Num"]
129                         if (length(result) == 0) {
130                                 stop(paste("Could not find function with name", algebra[[1]],
131                                         "and", (length(algebra) - 1), "arguments"))
132                         } else if (length(result) > 1) {
133                                 stop(paste("Ambiguous function with name", algebra[[1]],
134                                         "and", (length(algebra) - 1), "arguments"))
135                         } else {
136                                 head <- as.integer(result[[1]])
137                             tail <- lapply(algebra[-1], substituteOperators)
138                                 result <- append(tail, head, after=0)
139                                 return(result)
140                         }
141         }
142         }
143         return(algebra)
144 }
145
146 checkAlgebras <- function(model, flatModel) {
147         cycleDetection(flatModel)
148         checkAlgebraEvaluation(model, flatModel)
149 }
150
151 checkAlgebraEvaluation <- function(model, flatModel) {
152         if(length(flatModel@algebras) == 0) { return() }
153         for(i in 1:length(flatModel@algebras)) {
154                 algebra <- flatModel@algebras[[i]]
155                 expr <- substitute(mxEval(x, model, compute=TRUE),
156                         list(x = as.symbol(algebra@name)))
157                 tryCatch(eval(expr), error = function(x) {
158                         stop(paste("The algebra", 
159                                 omxQuotes(simplifyName(algebra@name, model@name)), 
160                                 "in model", omxQuotes(model@name), 
161                                 "generated the error message:",
162                                 x$message), call. = FALSE)
163                 })
164         }
165 }
166
167 displayAlgebra <- function(mxAlgebra) {
168         cat("mxAlgebra", omxQuotes(mxAlgebra@name), '\n')
169         cat("@formula: ", deparse(mxAlgebra@formula, width.cutoff=500L), '\n')
170         if (length(mxAlgebra@result) == 0) {
171                 cat("@result: (not yet computed) ")
172         } else {
173                 cat("@result:\n")
174         }
175         print(mxAlgebra@result)
176         if (is.null(dimnames(mxAlgebra))) {
177                         cat("dimnames: NULL\n")
178         } else {
179                 cat("dimnames:\n")
180                 print(dimnames(mxAlgebra))
181         }
182 }
183
184 setMethod("print", "MxAlgebra", function(x,...) { displayAlgebra(x) })
185 setMethod("show", "MxAlgebra", function(object) { displayAlgebra(object) })