Fixed a bug where could not be accessed from a model object, and adjusted one test...
[openmx:openmx.git] / R / MxModelFunctions.R
1 #
2 #   Copyright 2007-2015 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 generateMatrixList <- function(model) {
18         matvalues <- lapply(model@matrices, generateMatrixValuesHelper)
19         matnames  <- names(model@matrices)
20         names(matvalues) <- matnames
21         references <- generateMatrixReferences(model)
22         retval <- mapply(function(x,y) { c(list(x), y) }, 
23                         matvalues, references, SIMPLIFY = FALSE)
24         return(retval)
25 }
26
27 generateAlgebraList <- function(model) {
28         mNames <- names(model@matrices)
29         aNames <- append(names(model@algebras), names(model@fitfunctions))      
30         mNumbers <- as.list(as.integer(-1 : (-length(mNames))))
31         aNumbers <- as.list(as.integer(0 : (length(aNames) - 1)))
32         names(mNumbers) <- mNames
33         names(aNumbers) <- aNames
34         retval <- lapply(model@algebras, generateAlgebraHelper, 
35                 mNumbers, aNumbers)
36     return(retval)
37 }
38
39 findDependencies <- function(triple, flatModel, dependencies) {
40         mNames <- names(flatModel@matrices)
41         matrixNum <- triple[[1]] + 1
42         matrixName <- mNames[[matrixNum]]
43         return(dependencies[[matrixName]])      
44 }
45
46         
47 isExpectation <- function(name) {
48         return(length(grep("expectation", name, fixed=TRUE)) > 0)
49 }
50
51 parameterDependencyList <- function(pList, flatModel, dependencies) {
52         if (length(pList) == 3) {
53                 return(retval)
54         }
55         locations <- pList[4:length(pList)]
56         deps <- lapply(locations, findDependencies, flatModel, dependencies)
57         depnames <- Reduce(union, deps, character())
58         depnames <- Filter(Negate(isExpectation), depnames)
59         depnumbers <- sapply(depnames, doLocateIndex, flatModel, flatModel@name, USE.NAMES=FALSE)
60         depnumbers <- as.integer(depnumbers)
61         retval <- append(pList[1:3], list(depnumbers))
62         append(retval, locations)
63 }
64
65 buildFreeVarGroupList <- function(flatModel) {
66         if (is.null(flatModel@compute)) return()
67         got <- getFreeVarGroup(flatModel@compute)
68         if (length(got)) {
69                 if (anyDuplicated(unlist(got[seq(1,length(got),2)]))) {
70                         stop("freeSet IDs are not unique") #for debugging
71                 }
72                 members <- unlist(got[seq(2,length(got),2)])
73                 recog <- match(members, names(flatModel@matrices))
74                 if (any(is.na(recog))) {
75                         stop(paste("freeSet component", omxQuotes(members[is.na(recog)]), "not found"))
76                 }
77         }
78         got
79 }
80
81 generateParameterList <- function(flatModel, dependencies, freeVarGroups) {
82         mList <- flatModel@matrices
83         pList <- list()
84         if (length(mList)) for(i in 1:length(mList)) {
85                 matrix <- mList[[i]]
86                 pList <- generateParameterListHelper(matrix, pList, i - 1L, freeVarGroups)
87         }
88         pList <- lapply(pList, parameterDependencyList, flatModel, dependencies)
89
90         if (length(pList)) for(i in 1:length(pList)) {
91                 original <- pList[[i]]
92                 svalues <- original[5:length(original)]
93                 svalue <- NA
94                 if (length(svalues) > 1) {
95                         values <- sapply(svalues, generateValueHelper, mList)
96                         if (!all(values == values[[1]])) {
97                                 warning(paste('Parameter',i,'has multiple start values.',
98                                               'Selecting', values[[1]]))
99                         }
100                         svalue <- values[[1]]
101                 } else {
102                         svalue <- generateValueHelper(svalues[[1]], mList)
103                 }
104                 pList[[i]] <- c(original, svalue)
105         }
106         flatModel@parameters <- pList
107         flatModel
108 }
109
110 definitionDependencyList <- function(pList, flatModel, dependencies) {
111         if (length(pList) == 2) {
112                 retval <- list(pList[[1]], pList[[2]], integer())
113                 return(retval)
114         }
115         locations <- pList[3:length(pList)]
116         deps <- lapply(locations, findDependencies, flatModel, dependencies)
117         depnames <- Reduce(union, deps, character())
118         depnames <- Filter(Negate(isExpectation), depnames)
119         depnumbers <- sapply(depnames, doLocateIndex, flatModel, flatModel@name, USE.NAMES=FALSE)
120         depnumbers <- as.integer(depnumbers)
121         retval <- list(pList[[1]], pList[[2]], depnumbers)
122         return(append(retval, locations))
123 }
124
125 generateDefinitionList <- function(flatModel, dependencies) {
126         if (length(flatModel@matrices) == 0) return(list())
127         result <- list()
128         for(i in 1:length(flatModel@matrices)) {
129                 result <- matrixDefinitions(flatModel,
130                         flatModel@matrices[[i]], 
131                         result, i - 1L)
132         }
133         result <- lapply(result, definitionDependencyList, flatModel, dependencies)
134         return(result)
135 }
136
137 generateValueHelper <- function(triple, mList) {
138         mat <- triple[1] + 1
139         row <- triple[2] + 1
140         col <- triple[3] + 1
141         val <- mList[[mat]]@values[row,col]
142         if (is.na(val)) {
143                 stop(paste("Starting value in ",names(mList)[[mat]],
144                            "[",row,",",col,"] is missing", sep=""))
145         }
146         return(val)
147 }
148
149 imxUpdateModelValues <- function(model, flatModel, values) {
150         pList <- flatModel@parameters
151         if(length(pList) != length(values)) {
152                 stop(paste("This model has", length(pList), 
153                         "parameters, but you have given me", length(values),
154                         "values"))
155         }
156         if (length(pList) == 0) {
157                 return(model)
158         }
159         for(i in 1:length(pList)) {
160                 parameters <- pList[[i]]
161                 parameters <- parameters[5:(length(parameters)-1)]
162                 model <- updateModelValuesHelper(parameters, values[[i]], flatModel, model)
163     }
164         return(model)
165 }
166
167 updateModelValuesHelper <- function(triples, value, flatModel, model) {
168         for(i in 1:length(triples)) {
169                 triple <- triples[[i]]
170                 mat <- triple[1] + 1
171                 row <- triple[2] + 1
172                 col <- triple[3] + 1
173                 name <- flatModel@matrices[[mat]]@name
174                 model[[name]]@values[row,col] <- value
175         }
176         return(model)
177 }
178
179 removeTail <- function(lst, tailSize) {
180     newEnd <- length(lst) - tailSize
181     if (newEnd == 0) {
182         return(list())
183     } else {
184         return(lst[1 : newEnd])
185     }
186 }
187
188 updateModelMatrices <- function(model, flatModel, values) {
189         mList <- names(flatModel@matrices)
190         if (length(mList) != length(values)) {
191                 stop(paste("This model has", length(mList), 
192                         "matrices, but the backend has returned", length(values),
193                         "values"))
194         }
195         if (length(values) == 0) {
196                 return(model)
197         }
198         model <- updateModelEntitiesHelper(mList, values, model)
199         return(model)
200 }
201
202
203 updateModelAlgebras <- function(model, flatModel, values) {
204         aNames <- names(flatModel@algebras)
205         oNames <- names(flatModel@fitfunctions)
206         aList <- append(aNames, oNames)
207         if(length(aList) != length(values)) {
208                 stop(paste("This model has", length(aList), 
209                         "algebras, but the backend has returned", length(values),
210                         "values"))
211         }
212         if (length(aList) == 0) {
213                 return(model)
214         }
215         model <- updateModelEntitiesHelper(aList, values, model)
216         return(model)
217 }
218
219 updateModelExpectations <- function(model, flatModel, values) {
220         eNames <- names(flatModel@expectations)
221         if (length(eNames) != length(values)) {
222                 stop(paste("This model has", length(eNames),
223                            "expectations, but the backend has returned", length(values),
224                            "values"))
225         }
226         if (length(eNames) == 0) return(model)
227         updateModelEntitiesHelper(eNames, values, model)
228 }
229
230 updateModelData <- function(model, flatModel, values) {
231         dNames <- names(flatModel@datasets)
232         if (length(dNames) != length(values)) {
233                 stop(paste("This model has", length(dNames),
234                            "data, but the backend has returned", length(values),
235                            "values"))
236         }
237         if (length(dNames) == 0) return(model)
238         updateModelEntitiesHelper(dNames, values, model)
239 }
240
241 updateModelEntitiesTargetModel <- function(model, entNames, values, modelNameMapping) {
242     nextName <- model@name
243     selectEnt <- entNames[modelNameMapping == nextName]
244     selectVal <- values[modelNameMapping == nextName]
245     if (length(selectEnt) > 0) {
246                 for(i in 1:length(selectEnt)) {
247                         name <- selectEnt[[i]]
248                         candidate <- model[[name]]
249                         value <- selectVal[[i]]
250                         if (!is.null(candidate) && (length(value) > 0)
251                                 && !is.nan(value)) {
252                                 if (is(candidate, "MxAlgebra")) {
253                                         cdim <- sapply(dimnames(candidate), length)
254                                         mask <- cdim != 0
255                                         if (any(cdim[mask] != dim(value)[mask])) {
256                                                 warning(paste(paste(model@name, candidate@name, sep="."),
257                                                               "is dimension", omxQuotes(dim(value)),
258                                                               "but the dimnames imply dimension",
259                                                               omxQuotes(cdim)))
260                                         } else {
261                                                 dimnames(value) <- dimnames(candidate)
262                                         }
263                                         candidate@result <- value
264                                 } else if (is(candidate,"MxFitFunction")) {
265                                         candidate@result <- as.matrix(value)
266                                         attr <- attributes(value)
267                                         attr$dim <- NULL
268                                         candidate@info <- attr
269                                 } else if(is(candidate, "MxMatrix")) {
270                                         if (any(dim(candidate@values) != dim(value))) {
271                                                 msg <- paste("Backend returned a", omxQuotes(dim(value)),
272                                                              "matrix for a", omxQuotes(dim(candidate@values)),
273                                                              "matrix. Not sure how to proceed.",
274                                                              "Please report this error to the OpenMx support team.")
275                                                 stop(msg)
276                                         } else {
277                                                 dimnames(value) <- dimnames(candidate)
278                                                 candidate@values <- value
279                                         }
280                                 } else if (is(candidate, "MxExpectation")) {
281                                         for (sl in names(attributes(value))) {
282                                                 slot(candidate, sl) <- attr(value, sl)
283                                         }
284                                 } else if (is(candidate, "MxDataDynamic")) {
285                                         candidate@numObs <- value
286                                 }
287                                 model[[name]] <- candidate
288                         }
289                 }
290         }
291     if (length(model@submodels) > 0) {
292         model@submodels <- lapply(model@submodels, 
293             updateModelEntitiesTargetModel, entNames, values, modelNameMapping)
294     }
295         return(model)
296 }
297
298 updateModelEntitiesHelper <- function(entNames, values, model) {
299     modelNameMapping <- sapply(entNames, getModelNameString)
300     model <- updateModelEntitiesTargetModel(model, entNames, 
301                 values, modelNameMapping)
302     return(model)
303 }
304
305 clearModifiedSinceRunRecursive <- function(model) {
306         if (length(model@submodels) > 0) {
307                 model@submodels <- lapply(model@submodels, clearModifiedSinceRunRecursive)
308         }
309         model@.wasRun <- TRUE
310         model@.modifiedSinceRun <- FALSE
311         model
312 }
313
314 ##' imxLocateIndex
315 ##'
316 ##' This is an internal function exported for those people who know
317 ##' what they are doing.
318 ##'
319 ##' @param model model
320 ##' @param name name
321 ##' @param referant referant
322 imxLocateIndex <- function(model, name, referant) {
323         if (length(name) == 0) return(name)
324         if (is.na(name)) { return(as.integer(name)) }
325         mNames <- names(model@matrices)
326         aNames <- names(model@algebras)
327         fNames <- names(model@fitfunctions)
328         eNames <- names(model@expectations)
329         dNames <- names(model@datasets)         
330         matrixNumber <- match(name, mNames)
331         algebraNumber <- match(name, append(aNames, fNames))
332         dataNumber <- match(name, dNames)
333         expectationNumber <- match(name, eNames)
334         if (is.na(matrixNumber) && is.na(algebraNumber) 
335                 && is.na(dataNumber) && is.na(expectationNumber)) {
336                 reftype <- "named reference"
337                 if (typeof(referant) == "S4") {
338                         reftype <- 'object'
339                 }
340                 msg <- paste("The reference", omxQuotes(name),
341                         "does not exist.  It is used by",
342                              reftype, omxQuotes(referant), ".")
343                 stop(msg, call.=FALSE)
344         } else if (!is.na(matrixNumber)) {
345                 return(- matrixNumber)
346         } else if (!is.na(dataNumber)) {
347                 return(dataNumber - 1L)
348         } else if (!is.na(expectationNumber)) {
349                 return(expectationNumber - 1L)
350         } else {
351                 return(algebraNumber - 1L)
352         }
353 }
354
355 LocateOptionalMatrix <- function(model, name, referant) {
356         if (is.na(name)) { return(as.integer(name)) }
357         mNames <- names(model@matrices)
358         aNames <- names(model@algebras)
359         fNames <- names(model@fitfunctions)
360         matrixNumber <- match(name, mNames)
361         algebraNumber <- match(name, append(aNames, fNames))
362         if (is.na(matrixNumber) && is.na(algebraNumber)) {
363                 return(NULL)
364         } else if (!is.na(matrixNumber)) {
365                 return(- matrixNumber)
366         } else {
367                 return(algebraNumber - 1L)
368         }
369 }
370
371 zapExtraneousMatrices <- function(model){
372   keepers <- unlist(lapply(model@matrices,function(x){ifelse(".persist" %in% slotNames(x),x@.persist,TRUE)}))
373   #^^^^^For the sake of backward compatibility, treat matrices with no .persist slot as having that slot be TRUE.
374   if(is.logical(keepers)){model@matrices <- model@matrices[which(keepers)]}
375   if(length(model@submodels)>0){model@submodels <- lapply(model@submodels,zapExtraneousMatrices)}
376   return(model)
377 }
378
379 ##' imxPreprocessModel
380 ##'
381 ##' This is an internal function exported for those people who know
382 ##' what they are doing.
383 ##'
384 ##' @param model model
385 imxPreprocessModel <- function(model) {
386         model@matrices <- lapply(model@matrices, findSquareBrackets)
387         model@submodels <- lapply(model@submodels, imxPreprocessModel)
388         return(model)
389 }
390
391
392 ##' imxCheckMatrices
393 ##'
394 ##' This is an internal function exported for those people who know
395 ##' what they are doing.
396 ##'
397 ##' @param model model
398 imxCheckMatrices <- function(model) {
399         matrices <- model@matrices
400         lapply(matrices, imxVerifyMatrix)
401         submodels <- imxDependentModels(model)
402         lapply(submodels, imxCheckMatrices)
403 }