More progress on LISREL type models. Changed imxExpectationLISREL to mx prefix....
[openmx:openmx.git] / R / MxExpectationStateSpace.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
17 #--------------------------------------------------------------------
18 # Author: Michael D. Hunter
19 # Filename: MxExpectationStateSpace.R
20 # Date: 2012.11.14
21 # Purpose: Define classes and methods for the state space model (SSM)
22 #  expectations.
23 #--------------------------------------------------------------------
24
25 #--------------------------------------------------------------------
26 # Revision History
27 #   Wed Nov 14 13:34:01 Central Standard Time 2012 -- Michael Hunter created file
28 #   Sat Nov 17 16:14:56 Central Standard Time 2012 -- Michael Hunter change names to ExpectationStateSpace instead of ExpectationSSM
29 #   
30
31 #--------------------------------------------------------------------
32 # **DONE**
33 setClass(Class = "MxExpectationStateSpace",
34         representation = representation(
35                 A = "MxCharOrNumber",
36                 B = "MxCharOrNumber",
37                 C = "MxCharOrNumber",
38                 D = "MxCharOrNumber",
39                 Q = "MxCharOrNumber",
40                 R = "MxCharOrNumber",
41                 x = "MxCharOrNumber",
42                 P = "MxCharOrNumber",
43                 thresholds = "MxCharOrNumber",
44                 dims = "character",
45                 definitionVars = "list",
46                 dataColumns = "numeric",
47                 thresholdColumns = "numeric",
48                 thresholdLevels = "numeric",
49                 threshnames = "character"),
50         contains = "MxBaseExpectation")
51
52
53 #--------------------------------------------------------------------
54 # **DONE**
55 setMethod("initialize", "MxExpectationStateSpace",
56         function(.Object, A, B, C, D, Q, R, x, P, dims, thresholds, threshnames,
57                 data = as.integer(NA), name = 'expectation') {
58                 .Object@name <- name
59                 .Object@A <- A
60                 .Object@B <- B
61                 .Object@C <- C
62                 .Object@D <- D
63                 .Object@Q <- Q
64                 .Object@R <- R
65                 .Object@x <- x
66                 .Object@P <- P
67                 .Object@data <- data
68                 .Object@dims <- dims
69                 .Object@thresholds <- thresholds
70                 .Object@definitionVars <- list()
71                 .Object@threshnames <- threshnames
72                 return(.Object)
73         }
74 )
75
76
77 #--------------------------------------------------------------------
78 # TODO: Ask Spiegel and Brick what this is really supposed to do.
79 setMethod("genericExpConvertEntities", "MxExpectationStateSpace",
80         function(.Object, flatModel, namespace, labelsData) {
81                 cache <- new.env(parent = emptyenv())
82                 if(is.na(.Object@data)) {
83                         msg <- paste("The SSM expectation function",
84                                 "does not have a dataset associated with it in model",
85                                 omxQuotes(model@name))
86                         stop(msg, call.=FALSE)
87                 }
88                 flatModel <- updateThresholdDimnames(.Object, flatModel, labelsData)
89                 
90                 return(flatModel)
91         }
92 )
93
94
95 #--------------------------------------------------------------------
96 # **DONE**
97 setMethod("qualifyNames", signature("MxExpectationStateSpace"), 
98         function(.Object, modelname, namespace) {
99                 .Object@name <- imxIdentifier(modelname, .Object@name)
100                 .Object@A <- imxConvertIdentifier(.Object@A, modelname, namespace)
101                 .Object@B <- imxConvertIdentifier(.Object@B, modelname, namespace)
102                 .Object@C <- imxConvertIdentifier(.Object@C, modelname, namespace)
103                 .Object@D <- imxConvertIdentifier(.Object@D, modelname, namespace)
104                 .Object@Q <- imxConvertIdentifier(.Object@Q, modelname, namespace)
105                 .Object@R <- imxConvertIdentifier(.Object@R, modelname, namespace)
106                 .Object@x <- imxConvertIdentifier(.Object@x, modelname, namespace)
107                 .Object@P <- imxConvertIdentifier(.Object@P, modelname, namespace)
108                 .Object@data <- imxConvertIdentifier(.Object@data, modelname, namespace)
109                 .Object@thresholds <- sapply(.Object@thresholds, imxConvertIdentifier, modelname, namespace)
110                 return(.Object)
111         }
112 )
113
114
115 #--------------------------------------------------------------------
116 # TODO: Add lots of error checking.
117 setMethod("genericExpFunConvert", signature("MxExpectationStateSpace"), 
118         function(.Object, flatModel, model, labelsData, defVars, dependencies) {
119                 modelname <- imxReverseIdentifier(model, .Object@name)[[1]]     
120                 name <- .Object@name
121                 aMatrix <- .Object@A
122                 bMatrix <- .Object@B
123                 cMatrix <- .Object@C
124                 dMatrix <- .Object@D
125                 qMatrix <- .Object@Q
126                 rMatrix <- .Object@R
127                 xMatrix <- .Object@x
128                 pMatrix <- .Object@P
129                 data <- .Object@data
130                 if(is.na(data)) {
131                         msg <- paste("The SSM expectation function",
132                                 "does not have a dataset associated with it in model",
133                                 omxQuotes(modelname))
134                         stop(msg, call. = FALSE)
135                 }
136                 mxDataObject <- flatModel@datasets[[.Object@data]]
137                 checkNumericData(mxDataObject)
138                 .Object@A <- imxLocateIndex(flatModel, aMatrix, name)
139                 .Object@B <- imxLocateIndex(flatModel, bMatrix, name)
140                 .Object@C <- imxLocateIndex(flatModel, cMatrix, name)
141                 .Object@D <- imxLocateIndex(flatModel, dMatrix, name)
142                 .Object@Q <- imxLocateIndex(flatModel, qMatrix, name)
143                 .Object@R <- imxLocateIndex(flatModel, rMatrix, name)
144                 .Object@x <- imxLocateIndex(flatModel, xMatrix, name)
145                 .Object@P <- imxLocateIndex(flatModel, pMatrix, name)
146                 .Object@data <- as.integer(imxLocateIndex(flatModel, data, name))
147                 aMatrix <- flatModel[[aMatrix]]
148                 bMatrix <- flatModel[[bMatrix]]
149                 cMatrix <- flatModel[[cMatrix]]
150                 dMatrix <- flatModel[[dMatrix]]
151                 qMatrix <- flatModel[[qMatrix]]
152                 rMatrix <- flatModel[[rMatrix]]
153                 xMatrix <- flatModel[[xMatrix]]
154                 pMatrix <- flatModel[[pMatrix]]
155                 translatedNames <- c(dimnames(cMatrix)[[1]])
156                 if (mxDataObject@type == 'raw') {
157                         threshName <- .Object@thresholds
158                         checkNumberOrdinalColumns(mxDataObject)
159                         .Object@definitionVars <- imxFilterDefinitionVariables(defVars, data)
160                         .Object@dataColumns <- generateDataColumns(flatModel, translatedNames, data)
161                         verifyThresholds(flatModel, model, labelsData, data, translatedNames, threshName)
162                         .Object@thresholds <- imxLocateIndex(flatModel, threshName, name)
163                         retval <- generateThresholdColumns(flatModel, model, labelsData, translatedNames, data, threshName)
164                         .Object@thresholdColumns <- retval[[1]]
165                         .Object@thresholdLevels <- retval[[2]]
166                         if (length(mxDataObject@observed) == 0) {
167                                 .Object@data <- as.integer(NA)
168                         }
169                         if (single.na(.Object@dims)) {
170                                 .Object@dims <- translatedNames
171                         }
172                 }
173                 return(.Object)
174         }
175 )
176
177
178 #--------------------------------------------------------------------
179 # **DONE**
180 setMethod("genericExpDependencies", signature("MxExpectationStateSpace"),
181         function(.Object, dependencies) {
182                 sources <- c(.Object@A, .Object@B, .Object@C, .Object@D, .Object@Q, .Object@R, .Object@x, .Object@P, .Object@thresholds)
183                 sources <- sources[!is.na(sources)]
184                 dependencies <- imxAddDependency(sources, .Object@name, dependencies)
185                 return(dependencies)
186         }
187 )
188
189
190 #--------------------------------------------------------------------
191 # **DONE**
192 setMethod("genericExpRename", signature("MxExpectationStateSpace"),
193         function(.Object, oldname, newname) {
194                 .Object@A <- renameReference(.Object@A, oldname, newname)
195                 .Object@B <- renameReference(.Object@B, oldname, newname)
196                 .Object@C <- renameReference(.Object@C, oldname, newname)
197                 .Object@D <- renameReference(.Object@D, oldname, newname)
198                 .Object@Q <- renameReference(.Object@Q, oldname, newname)
199                 .Object@R <- renameReference(.Object@R, oldname, newname)
200                 .Object@x <- renameReference(.Object@x, oldname, newname)
201                 .Object@P <- renameReference(.Object@P, oldname, newname)
202                 .Object@data <- renameReference(.Object@data, oldname, newname)
203                 .Object@thresholds <- sapply(.Object@thresholds, renameReference, oldname, newname)             
204                 return(.Object)
205         }
206 )
207
208
209 #--------------------------------------------------------------------
210 # Note: this turns off data sorting for the State Space expectation
211 setMethod("genericExpAddEntities", "MxExpectationStateSpace",
212         function(.Object, job, flatJob, labelsData) {
213                 # Do NOT sort data for state space models
214                 if(is.null(job@options$`No Sort Data`)){
215                     # Only add the option if it wasn't there before.
216                     # Do not clobber existing option
217                     key <- "No Sort Data"
218                     value <- getModelName(.Object)
219                     job <- mxOption(job, key, value)
220                 }
221                 
222                 # Run state space models single threaded
223                 key <- "Number of Threads"
224                 job <- mxOption(job, key, 1)
225                 return(job)
226         }
227 )
228
229
230 #--------------------------------------------------------------------
231 checkSSMargument <- function(x, xname) {
232         if (!(single.na(x) || typeof(x) == "character")) {
233                 msg <- paste("argument ", xname, " is not a string ",
234                         "(the name of the '", xname, "' matrix)", sep="")
235                 stop(msg)
236         }
237         if (is.na(x)) x <- as.integer(NA)
238         return(x)
239 }
240
241 #--------------------------------------------------------------------
242 # **DONE**
243 imxExpectationStateSpace <- function(A, B, C, D, Q, R, x, P, dimnames = NA, thresholds = NA, threshnames = dimnames){
244         A <- checkSSMargument(A, "A")
245         B <- checkSSMargument(B, "B")
246         C <- checkSSMargument(C, "C")
247         D <- checkSSMargument(D, "D")
248         Q <- checkSSMargument(Q, "Q")
249         R <- checkSSMargument(R, "R")
250         x <- checkSSMargument(x, "x")
251         P <- checkSSMargument(P, "P")
252         if (single.na(thresholds)) thresholds <- as.character(NA)
253         if (single.na(dimnames)) dimnames <- as.character(NA)
254         if (single.na(threshnames)) threshnames <- as.character(NA)
255         if (!is.vector(dimnames) || typeof(dimnames) != 'character') {
256                 stop("Dimnames argument is not a character vector")
257         }
258         if (!is.vector(threshnames) || typeof(threshnames) != 'character') {
259                 stop("'threshnames' argument is not a character vector")
260         }
261         if (length(thresholds) != 1) {
262                 stop("Thresholds argument must be a single matrix or algebra name")
263         }
264         if (length(dimnames) == 0) {
265                 stop("Dimnames argument cannot be an empty vector")
266         }
267         if (length(threshnames) == 0) {
268                 stop("'threshnames' argument cannot be an empty vector")
269         }
270         if (length(dimnames) > 1 && any(is.na(dimnames))) {
271                 stop("NA values are not allowed for dimnames vector")
272         }
273         if (length(threshnames) > 1 && any(is.na(threshnames))) {
274                 stop("NA values are not allowed for 'threshnames' vector")
275         }
276         return(new("MxExpectationStateSpace", A, B, C, D, Q, R, x, P, dimnames, thresholds, threshnames))
277 }
278
279
280 #--------------------------------------------------------------------
281 # TODO: Add expected mean and cov printouts
282 displayMxExpectationStateSpace <- function(expectation) {
283         cat("MxExpectationStateSpace", omxQuotes(expectation@name), '\n')
284         cat("@A :", omxQuotes(expectation@A), '\n')
285         cat("@B :", omxQuotes(expectation@B), '\n')
286         cat("@C :", omxQuotes(expectation@C), '\n')
287         cat("@D :", omxQuotes(expectation@D), '\n')
288         cat("@Q :", omxQuotes(expectation@Q), '\n')
289         cat("@R :", omxQuotes(expectation@R), '\n')
290         cat("@x :", omxQuotes(expectation@x), '\n')
291         cat("@P :", omxQuotes(expectation@P), '\n')
292         if (single.na(expectation@dims)) {
293                 cat("@dims : NA \n")
294         } else {
295                 cat("@dims :", omxQuotes(expectation@dims), '\n')
296         }               
297         if (single.na(expectation@thresholds)) {
298                 cat("@thresholds : NA \n")
299         } else {
300                 cat("@thresholds :", omxQuotes(expectation@thresholds), '\n')
301         }
302         invisible(expectation)
303 }
304
305 setMethod("print", "MxExpectationStateSpace", function(x,...) { 
306         displayMxExpectationStateSpace(x) 
307 })
308
309 setMethod("show", "MxExpectationStateSpace", function(object) { 
310         displayMxExpectationStateSpace(object) 
311 })
312
313