Added $ accessor for intervals in MxModels. model produces the same result as model...
[openmx:openmx.git] / R / MxModel.R
1 #
2 #   Copyright 2007-2014 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 setClass(Class = "MxModel",
17         representation = representation(
18                 name = "character",
19                 matrices = "list",
20                 algebras = "list",
21                 constraints = "list",
22                 intervals = "list",
23                 latentVars = "MxCharOrList",
24                 manifestVars = "MxCharOrList",
25                 data = "MxData",
26                 submodels = "list",
27                 expectation = "MxExpectation",
28                 fitfunction = "MxFitFunction",
29                 compute = "MxCompute",
30                 independent = "logical",
31                 options = "list",
32                 output = "list",
33                 runstate = "list",
34                 .forcesequential = "logical",
35                 .newobjects = "logical",
36                 .resetdata = "logical",
37                 .modifiedSinceRun = "logical"
38 ))
39
40 imxModelTypes[['default']] <- "MxModel"
41
42 setMethod("initialize", "MxModel",
43         function(.Object, name = character()) {
44                 .Object@name <- name
45                 .Object@latentVars <- character()
46                 .Object@manifestVars <- character()
47                 .Object@matrices <- list()
48                 .Object@algebras <- list()
49                 .Object@constraints <- list()
50                 .Object@data <- NULL
51                 .Object@submodels <- list()
52                 .Object@expectation <- NULL
53                 .Object@fitfunction <- NULL
54                 .Object@compute <- NULL
55                 .Object@independent <- FALSE
56                 .Object@options <- list()
57                 .Object@output <- list()
58                 .Object@runstate <- list()
59                 .Object@.newobjects <- FALSE
60                 .Object@.resetdata <- FALSE
61                 .Object@.modifiedSinceRun <- FALSE
62                 .Object <- imxInitModel(.Object)
63                 return(.Object)
64         }
65 )
66
67 ##' imxInitModel
68 ##'
69 ##' This is an internal function exported for those people who know
70 ##' what they are doing.
71 ##'
72 ##' @param model model
73 ##' @aliases
74 ##' imxInitModel,MxModel-method
75 ##' imxInitModel,MxRAMModel-method
76 ##' imxInitModel,MxLISRELModel-method
77 setGeneric("imxInitModel", function(model) {
78         return(standardGeneric("imxInitModel")) } )
79
80 ##' imxModelBuilder
81 ##'
82 ##' This is an internal function exported for those people who know
83 ##' what they are doing.
84 ##' 
85 ##' TODO: It probably makes sense to split this into separate
86 ##' methods. For example, modelAddVariables and modelRemoveVariables
87 ##' could be their own methods. This would reduce some cut&paste
88 ##' duplication.
89 ##'
90 ##' @param model model
91 ##' @param lst lst
92 ##' @param name name
93 ##' @param manifestVars manifestVars
94 ##' @param latentVars latentVars
95 ##' @param submodels submodels
96 ##' @param remove remove
97 ##' @param independent independent
98 ##' @aliases
99 ##' imxModelBuilder,MxLISRELModel-method
100 ##' imxModelBuilder,MxModel-method
101 ##' imxModelBuilder,MxRAMModel-method
102 setGeneric("imxModelBuilder", function(model, lst, name, 
103         manifestVars, latentVars, submodels, remove, independent) {
104         return(standardGeneric("imxModelBuilder")) } )
105
106 ##' imxTypeName
107 ##'
108 ##' This is an internal function exported for those people who know
109 ##' what they are doing.
110 ##'
111 ##' @param model model
112 ##' @aliases
113 ##' imxTypeName,MxLISRELModel-method
114 ##' imxTypeName,MxModel-method
115 ##' imxTypeName,MxRAMModel-method
116 setGeneric("imxTypeName", function(model) { 
117         return(standardGeneric("imxTypeName")) 
118 })
119
120 ##' imxVerifyModel
121 ##' 
122 ##' This is an internal function exported for those people who know
123 ##' what they are doing.
124 ##'
125 ##' @param model model
126 ##' @aliases
127 ##' imxVerifyModel,MxLISRELModel-method
128 ##' imxVerifyModel,MxModel-method
129 ##' imxVerifyModel,MxRAMModel-method
130 setGeneric("imxVerifyModel", function(model) {
131     return(standardGeneric("imxVerifyModel"))
132 })
133
134 # End declaration of generics
135
136 generateParentNames <- function(model) {
137         retval <- generateLocalNames(model)
138         if (length(model@submodels) > 0) {
139                 retval <- union(retval, names(model@submodels))
140                 childNames <- unlist(lapply(model@submodels, generateChildNames))
141                 retval <- union(retval, childNames)
142         }
143         return(retval)
144 }
145
146 generateChildNames <- function(model) {
147         retval <- generateLocalNames(model)     
148         if (!is.null(retval)) {
149                 retval <- paste(model@name, retval, sep = ".")
150         }
151         if (length(model@submodels) > 0) {
152                 retval <- union(retval, names(model@submodels))
153                 childNames <- unlist(lapply(model@submodels, generateChildNames))
154                 retval <- union(retval, childNames)
155         }
156         return(retval)
157 }
158
159 generateLocalNames <- function(model) {
160         matrices <- names(model@matrices)
161         algebras <- names(model@algebras)
162         constraints <- names(model@constraints)
163         retval <- union(matrices, algebras)
164         retval <- union(retval, constraints)
165         if (!is.null(model@fitfunction)) {
166                 retval <- union(retval, model@fitfunction@name)
167         }
168         if (!is.null(model@expectation)) {
169                 retval <- union(retval, model@expectation@name)
170         }
171         if (!is.null(model@data)) {
172                 retval <- union(retval, model@data@name)
173         }
174         return(retval)
175 }
176
177 setMethod("names", "MxModel",
178         function(x) {
179                 generateParentNames(x)
180         }
181 )
182
183 setMethod("[[", "MxModel",
184         function(x, i, j, ..., drop = FALSE) {
185                 return(imxExtractMethod(x, i))
186         }
187 )
188
189 setReplaceMethod("[[", "MxModel",
190         function(x, i, j, value) {
191                 return(imxReplaceMethod(x, i, value))
192         }
193 )
194
195 setMethod("$", "MxModel",
196         function(x, name) {
197         result <- imxExtractMethod(x, name)
198         if(name %in% c("name", "matrices", "algebras", "data", "submodels", "output", "compute", "options", "intervals")) {
199             result <- imxExtractSlot(x, name)
200         }
201                 return(result)
202         }
203 )
204
205 setReplaceMethod("$", "MxModel",
206         function(x, name, value) {
207                 if(name == "output") {
208                         stop("You cannot directly set the output of a model.  Use mxRun() if you want output.")
209                 } 
210                 if(name == "name") {
211                         stop("You cannot directly set the name of a model.  To rename the model, use model<-mxModel(model, name=\"NewName\").")
212                 }
213                 if(name %in% c("matrices", "algebras", "submodels")) {
214                         stop(paste("You cannot directly set the", name,
215                                    "of a model.  To set objects in the model, use the mxModel() function."))
216                 }
217                 return(imxReplaceMethod(x, name, value))
218         }
219 )
220
221
222 ##' imxSameType
223 ##' 
224 ##' This is an internal function exported for those people who know
225 ##' what they are doing.
226 ##'
227 ##' @param a a
228 ##' @param b b
229 imxSameType <- function(a, b) {
230         return( (is(a, "MxModel") && is(b, "MxModel")) ||
231                         (is(a, "MxMatrix") && is(b, "MxMatrix")) ||
232                         (is(a, "MxAlgebra") && is(b, "MxAlgebra")) ||
233                         (is(a, "MxExpectation") && is(b, "MxExpectation")) ||
234                         (is(a, "MxFitFunction") && is(b, "MxFitFunction")) ||
235                         (is(a, "MxCompute") && is(b, "MxCompute")) ||
236                         (is(a, "MxConstraint") && is(b, "MxConstraint")) ||
237                         (is(a, "MxData") && is(b, "MxData")))
238 }
239
240 mxModel <- function(model = NA, ..., manifestVars = NA, latentVars = NA,
241                     remove = FALSE, independent = NA, type = NA, name = NA) {
242         retval <- firstArgument(model, name)
243         first <- retval[[1]]
244         model <- retval[[2]]
245         name  <- retval[[3]]
246         model <- typeArgument(model, type)
247         lst <- c(first, list(...))
248         lst <- unlist(lst)
249         filter <- sapply(lst, is, "MxModel")
250         submodels <- lst[filter]
251         lst <- lst[!filter]
252         model <- imxModelBuilder(model, lst, name, manifestVars,
253                 latentVars, submodels, remove, independent)
254         return(model)
255 }
256
257 firstArgument <- function(model, name) {
258         first <- NULL
259         defaultType <- imxModelTypes[[getOption("mxDefaultType")]]
260         if (is(model, "MxModel")) {
261         } else {
262                 if (single.na(model)) {
263                 } else if (typeof(model) == "character") {
264                         name <- model
265                 } else if (isS4(model)) {
266                         first <- model
267                 } else {
268                         first <- list(model)
269                 }
270                 if (length(name) > 0 && is.na(name)) {
271                         name <- imxUntitledName()
272                 }
273                 imxVerifyName(name, -1)
274                 model <- new(defaultType, name)
275         }
276         return(list(first, model, name))
277 }
278
279 typeArgument <- function(model, type) {
280         if (!is.na(type)) {
281                 if (is.null(imxModelTypes[[type]])) {
282                         stop(paste("The model type", omxQuotes(type), 
283                                 "is not in the the list of acceptable types:",
284                                 omxQuotes(names(imxModelTypes))), call. = FALSE)
285                 }
286                 typename <- imxModelTypes[[type]]
287                 class(model) <- typename
288                 model <- imxInitModel(model)
289         }
290         return(model)
291 }
292
293 ##' imxGenericModelBuilder
294 ##'
295 ##' This is an internal function exported for those people who know
296 ##' what they are doing.
297 ##'
298 ##' @param model model
299 ##' @param lst lst
300 ##' @param name name
301 ##' @param manifestVars manifestVars
302 ##' @param latentVars latentVars
303 ##' @param submodels submodels
304 ##' @param remove remove
305 ##' @param independent independent
306 imxGenericModelBuilder <- function(model, lst, name, 
307         manifestVars, latentVars, submodels, remove, independent) {
308         model <- nameArgument(model, name)
309         model <- variablesArgument(model, manifestVars, latentVars, submodels, remove)
310         model <- listArgument(model, lst, remove)
311         model <- independentArgument(model, independent)
312         return(model)
313 }
314
315 varsToCharacter <- function(vars, vartype) {
316         if (is.list(vars)) {
317                 varnames <- names(vars)
318                 if (length(varnames) == 0) {
319                         return(as.character(unlist(vars)))      
320                 } else {
321                         result <- pmatch(varnames, imxVariableTypes)
322                         illegal <- which(is.na(result))
323                         if (length(illegal) > 0) {
324                                 if (length(illegal) == 1) {
325                                         ctgMsg <- "category"
326                                 } else {
327                                         ctgMsg <- "categories"
328                                 }
329                                 msg <- paste("In the", vartype, "variables",
330                                         "the", ctgMsg,
331                                         omxQuotes(varnames[illegal]), "did not match",
332                                         "to a valid category or two categories matched",
333                                         "to the same string (see 'imxVariableTypes'",
334                                         "for the list of legal categories)")
335                                 stop(msg, call. = FALSE)
336                         }
337                         varnames <- imxVariableTypes[result]
338                         vars <- lapply(vars, as.character)
339                         names(vars) <- varnames
340                         return(vars)
341                 }
342         } else {
343                 return(as.character(vars))
344         }
345 }
346
347 variablesArgument <- function(model, manifestVars, latentVars, submodels, remove) {
348         if (single.na(manifestVars)) {
349                 manifestVars <- character()
350         }
351         if (single.na(latentVars)) {
352                 latentVars <- character()
353         }
354         if (remove == TRUE) {
355                 model <- modelRemoveVariables(model, latentVars, manifestVars)
356                 if (length(submodels)) for(i in 1:length(submodels)) {
357                         model <- removeSingleNamedEntity(model, submodels[[i]])
358                 }
359         } else {
360                 if (length(manifestVars) + length(latentVars) > 0) {
361                         latentVars <- varsToCharacter(latentVars, "latent")
362                         manifestVars <- varsToCharacter(manifestVars, "manifest")
363                         checkVariables(model, latentVars, manifestVars)
364                         model <- modelAddVariables(model, latentVars, manifestVars)
365                 }
366                 if (length(submodels)) for(i in 1:length(submodels)) {
367                         model <- addSingleNamedEntity(model, submodels[[i]])
368                 }
369         }
370         return(model)
371 }
372
373 listArgument <- function(model, lst, remove) {
374         if(remove == TRUE) {
375                 model <- modelRemoveEntries(model, lst)
376         } else {
377                 model <- modelAddEntries(model, lst)
378         }
379         return(model)
380 }
381
382 independentArgument <- function(model, independent) {
383         if(!is.na(independent)) {
384                 model@independent <- independent
385         }
386         return(model)
387 }
388
389 nameArgument <- function(model, name) {
390         if(!is.na(name)) {
391                 model@name <- name
392         }
393         return(model)
394 }
395
396 checkVariables <- function(model, latentVars, manifestVars, submodels) {
397         latentVars   <- unlist(latentVars, use.names = FALSE)
398         manifestVars <- unlist(manifestVars, use.names = FALSE)
399         modelLatent <- unlist(model@latentVars, use.names = FALSE)
400         modelManifest <- unlist(model@manifestVars, use.names = FALSE)
401         common <- intersect(latentVars, manifestVars)    
402         if (length(common) > 0) {
403                 stop(paste("The following variables cannot",
404                         "be both latent and manifest:",
405                         omxQuotes(common)), call. = FALSE)
406         }
407         common <- intersect(modelLatent, manifestVars)
408         if (length(common) > 0) {
409                 stop(paste("The following variables cannot",
410                         "be both latent and manifest:",
411                         omxQuotes(common)), call. = FALSE)
412         }
413         common <- intersect(modelManifest, latentVars)
414         if (length(common) > 0) {
415                 stop(paste("The following variables cannot",
416                         "be both latent and manifest",
417                         omxQuotes(common)), call. = FALSE)
418         }
419         common <- intersect(modelManifest, manifestVars)
420         if (length(common) > 0) {
421                 stop(paste("The following manifest variables",
422                         "have already been declared",
423                         omxQuotes(common)), call. = FALSE)
424         }
425         common <- intersect(modelLatent, latentVars)
426         if (length(common) > 0) {
427                 stop(paste("The following latent variables",
428                         "have already been declared",
429                         omxQuotes(common)), call. = FALSE)
430         }
431         if (any(is.na(latentVars))) {
432                 stop("NA is not allowed as a latent variable", call. = FALSE)
433         }
434         if (any(is.na(manifestVars))) {
435                 stop("NA is not allowed as a manifest variable", call. = FALSE)
436         }
437         if (length(unique(latentVars)) != length(latentVars)) {
438                 stop(paste("The following variables in the latentVars list are duplicated:", 
439                 omxQuotes(latentVars[duplicated(latentVars)])), call. = FALSE)
440         }
441         if (length(unique(manifestVars)) != length(manifestVars)) {
442                 stop(paste("The following variables in the manifestVars list are duplicated:", 
443                 omxQuotes(manifestVars[duplicated(manifestVars)])), call. = FALSE)
444         }
445 }
446
447 # Begin implementation of generics
448
449 setMethod("imxModelBuilder", "MxModel", imxGenericModelBuilder)
450
451 setMethod("imxInitModel", "MxModel", function(model) { 
452         return(model)
453 })
454
455 setMethod("imxTypeName", "MxModel", function(model) { 
456         return("default")
457 })
458
459 setMethod("imxVerifyModel", "MxModel", function(model) {
460     return(TRUE)
461 })
462
463 # End implementation of generics
464
465 addVariablesHelper <- function(model, vartype, vars) {
466         modelvars <- slot(model, vartype)
467
468         if (length(vars) == 0) {
469                 return(model)
470         } else if (length(modelvars) == 0) {
471                 slot(model, vartype) <- vars
472                 return(model)
473         }
474
475         if (is.list(vars) && !is.list(modelvars)) {
476                 msg <- paste("The", vartype, "variables in",
477                         "the call to mxModel() have been separated",
478                         "into categories, and the existing", vartype,
479                         "variables do not have categories.")
480                 stop(msg, call. = FALSE)
481         } else if (!is.list(vars) && is.list(modelvars)) {
482                 msg <- paste("The", vartype, "variables in",
483                         "the call to mxModel() have not been separated",
484                         "into categories, and the existing", vartype,
485                         "variables do have categories.")
486                 stop(msg, call. = FALSE)
487         }
488
489         if (is.character(vars) && is.character(modelvars)) {
490                 modelvars <- c(modelvars, vars)
491                 slot(model, vartype) <- modelvars
492         } else {
493                 varnames <- names(vars)
494                 offsets <- pmatch(varnames, imxVariableTypes)
495                 for(i in 1:length(vars)) {
496                         currOffset <- offsets[[i]]
497                         currName <- imxVariableTypes[[currOffset]]
498                         modelvars[[currName]] <- c(modelvars[[currName]], vars[[i]])
499                 }
500                 slot(model, vartype) <- modelvars
501         }
502
503         return(model)
504 }
505
506 modelAddVariables <- function(model, latent, manifest) {
507         model <- addVariablesHelper(model, "latentVars", latent)
508         model <- addVariablesHelper(model, "manifestVars", manifest)
509         return(model)
510 }
511
512 modelRemoveVariables <- function(model, latent, manifest) {
513         model@latentVars <- setdiff(model@latentVars, latent)
514         model@manifestVars <- setdiff(model@manifestVars, manifest)
515         return(model)
516 }
517         
518 modelAddEntries <- function(model, entries) {
519         if (length(entries) == 0) {
520                 return(model)
521         }
522         tuple <- modelModifyFilter(model, entries, "add")
523         namedEntities <- tuple[[1]]
524         bounds        <- tuple[[2]]
525         intervals     <- tuple[[3]]
526         intervals     <- expandIntervals(intervals)
527         names(intervals) <- sapply(intervals, slot, "reference")
528         if (length(namedEntities) > 0) for(i in 1:length(namedEntities)) {
529                 model <- addSingleNamedEntity(model, namedEntities[[i]])
530         }
531         model <- modelAddBounds(model, bounds)
532         model <- modelAddIntervals(model, intervals)
533         return(model)
534 }
535
536 modelRemoveEntries <- function(model, entries) {
537         if (length(entries) == 0) {
538                 return(model)
539         }
540         tuple <- modelModifyFilter(model, entries, "remove")
541         namedEntities <- tuple[[1]]
542         bounds        <- tuple[[2]]
543         intervals     <- tuple[[3]]
544         intervals     <- expandIntervals(intervals)     
545         names(intervals) <- sapply(intervals, slot, "reference")
546         if (length(namedEntities) > 0) for(i in 1:length(namedEntities)) {
547                 model <- removeSingleNamedEntity(model, namedEntities[[i]])
548         }
549         model <- modelRemoveBounds(model, bounds)
550         model <- modelRemoveIntervals(model, intervals)
551         return(model)
552 }
553
554 actionCorrespondingPredicate <- c('add' = 'into', 'remove' = 'from')
555
556 modelModifyFilter <- function(model, entries, action) {
557         boundsFilter <- sapply(entries, is, "MxBounds")
558         intervalFilter <- sapply(entries, is, "MxInterval")
559         namedEntityFilter <- sapply(entries, function(x) {"name" %in% slotNames(x)})
560         characterFilter <- sapply(entries, is.character)
561         pathFilter <- sapply(entries, is, "MxPath")
562   thresholdFilter <- sapply(entries, is, "MxThreshold")
563         unknownFilter <- !(boundsFilter | namedEntityFilter | intervalFilter | characterFilter | thresholdFilter)
564         if (any(pathFilter)) {
565                 stop(paste("The model type of model",
566                         omxQuotes(model@name), "does not recognize paths."),
567                         call. = FALSE)
568         }
569         if (any(thresholdFilter)) {
570           stop(paste("The model type of model",
571                      omxQuotes(model@name), "does not recognize thresholds."),
572                call. = FALSE)
573         }
574         if (any(unknownFilter)) {
575                 stop(paste("Cannot", action, "the following item(s)", 
576                         actionCorrespondingPredicate[[action]], "the model:", 
577                         omxQuotes(sapply(entries[unknownFilter], deparse))), call. = FALSE)
578         }
579         if (any(namedEntityFilter) && action == 'remove') {
580                 stop(paste("Cannot use named entities when remove = TRUE.",
581                         "Instead give the name of the entity when removing it.",
582                         "See http://openmx.psyc.virginia.edu/wiki/mxmodel-help#Remove_an_object_from_a_model"))
583         }
584         if (any(characterFilter) && action == 'add') {
585                 stop(paste("I don't know what to do with the following strings",
586                         omxQuotes(entries[characterFilter]),
587                         "that have been passed into the function:",
588                         deparse(width.cutoff = 400L, imxLocateFunction("mxModel"))), call. = FALSE)
589         }
590         if (identical(action, 'add')) {
591                 return(list(entries[namedEntityFilter], entries[boundsFilter], entries[intervalFilter]))
592         } else if (identical(action, 'remove')) {
593                 return(list(entries[characterFilter], entries[boundsFilter], entries[intervalFilter]))
594         } else {
595                 stop(paste("Internal error, unidentified action:", omxQuotes(action)))
596         }
597 }
598
599 addSingleNamedEntity <- function(model, entity) {
600         if (!nzchar(entity@name)) {
601                 stop(paste("Entity",
602                         omxQuotes(class(entity)), "in model",
603                         omxQuotes(model@name), "needs a name"), call. = FALSE)
604         }
605         if (model@name == entity@name) {
606                 stop(paste("You cannot insert an entity named",
607                         omxQuotes(entity@name), "into a model named",
608                         omxQuotes(model@name)), call. = FALSE)
609         }
610         model[[entity@name]] <- entity
611         return(model)
612 }
613
614 removeSingleNamedEntity <- function(model, name) {
615         model[[name]] <- NULL
616         return(model)
617 }
618
619 setMethod("imxVerifyModel", "MxModel",
620     function(model) {
621         if (length(model@submodels) > 0) {
622                 return(all(sapply(model@submodels, imxVerifyModel)))
623         }
624         return(TRUE)
625     }
626 )