Rework specification of free param groups per Tim's suggestion
[openmx:openmx.git] / R / MxModelFunctions.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 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         result <- list()
67         for (step in flatModel@computes) {
68                 got <- getFreeVarGroup(step)
69                 if (length(got)) result <- append(result, got)
70         }
71         if (length(result)) {
72                 members <- unlist(result[seq(2,length(result),2)])
73                 recog <- match(members, names(flatModel@matrices))
74                 if (any(is.na(recog))) {
75                         stop(paste("free.set component", omxQuotes(members[is.na(recog)]), "not found"))
76                 }
77         }
78         result
79 }
80
81 generateParameterList <- function(flatModel, dependencies, freeVarGroups) {
82         mList <- flatModel@matrices
83         pList <- list()
84         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         result <- list()
127         defLocations <- generateDefinitionLocations(flatModel@datasets)
128         if (length(flatModel@matrices) == 0) {
129                 return(result)
130         }
131         for(i in 1:length(flatModel@matrices)) {
132                 result <- generateDefinitionListHelper(
133                         flatModel@matrices[[i]], 
134                         result, defLocations, i - 1L)
135         }
136         result <- lapply(result, definitionDependencyList, flatModel, dependencies)
137         return(result)
138 }
139
140 generateValueHelper <- function(triple, mList) {
141         mat <- triple[1] + 1
142         row <- triple[2] + 1
143         col <- triple[3] + 1
144         val <- mList[[mat]]@values[row,col]
145         if (is.na(val)) {
146                 stop(paste("Starting value in ",names(mList)[[mat]],
147                            "[",row,",",col,"] is missing", sep=""))
148         }
149         return(val)
150 }
151
152 imxUpdateModelValues <- function(model, flatModel, values) {
153         pList <- flatModel@parameters
154         if(length(pList) != length(values)) {
155                 stop(paste("This model has", length(pList), 
156                         "parameters, but you have given me", length(values),
157                         "values"))
158         }
159         if (length(pList) == 0) {
160                 return(model)
161         }
162         for(i in 1:length(pList)) {
163                 parameters <- pList[[i]]
164                 parameters <- parameters[5:(length(parameters)-1)]
165                 model <- updateModelValuesHelper(parameters, values[[i]], flatModel, model)
166     }
167         return(model)
168 }
169
170 updateModelValuesHelper <- function(triples, value, flatModel, model) {
171         for(i in 1:length(triples)) {
172                 triple <- triples[[i]]
173                 mat <- triple[1] + 1
174                 row <- triple[2] + 1
175                 col <- triple[3] + 1
176                 name <- flatModel@matrices[[mat]]@name
177                 model[[name]]@values[row,col] <- value
178         }
179         return(model)
180 }
181
182 removeTail <- function(lst, tailSize) {
183     newEnd <- length(lst) - tailSize
184     if (newEnd == 0) {
185         return(list())
186     } else {
187         return(lst[1 : newEnd])
188     }
189 }
190
191 updateModelMatrices <- function(model, flatModel, values) {
192         mList <- names(flatModel@matrices)
193         if (length(mList) != length(values)) {
194                 stop(paste("This model has", length(mList), 
195                         "matrices, but the backend has returned", length(values),
196                         "values"))
197         }
198         if (length(values) == 0) {
199                 return(model)
200         }
201         model <- updateModelEntitiesHelper(mList, values, model)
202         return(model)
203 }
204
205
206 updateModelAlgebras <- function(model, flatModel, values) {
207         aNames <- names(flatModel@algebras)
208         oNames <- names(flatModel@fitfunctions)
209         aList <- append(aNames, oNames)
210         if(length(aList) != length(values)) {
211                 stop(paste("This model has", length(aList), 
212                         "algebras, but the backend has returned", length(values),
213                         "values"))
214         }
215         if (length(aList) == 0) {
216                 return(model)
217         }
218         model <- updateModelEntitiesHelper(aList, values, model)
219         return(model)
220 }
221
222 updateModelExpectations <- function(model, flatModel, values) {
223         eNames <- names(flatModel@expectations)
224         if (length(eNames) != length(values)) {
225                 stop(paste("This model has", length(eNames),
226                            "expectations, but the backend has returned", length(values),
227                            "values"))
228         }
229         if (length(eNames) == 0) return(model)
230         updateModelEntitiesHelper(eNames, values, model)
231 }
232
233 updateModelEntitiesTargetModel <- function(model, entNames, values, modelNameMapping) {
234     nextName <- model@name
235     selectEnt <- entNames[modelNameMapping == nextName]
236     selectVal <- values[modelNameMapping == nextName]
237     if (length(selectEnt) > 0) {
238                 for(i in 1:length(selectEnt)) {
239                         name <- selectEnt[[i]]
240                         candidate <- model[[name]]
241                         value <- selectVal[[i]]
242                         if (!is.null(candidate) && (length(value) > 0)
243                                 && !is.nan(value)) {
244                                 if (is(candidate, "MxAlgebra")) {
245                                         dimnames(value) <- dimnames(candidate)
246                                         candidate@result <- value
247                                 } else if (is(candidate,"MxFitFunction")) {
248                                         candidate@result <- as.matrix(value[1,1])
249                                         attr <- attributes(value)
250                                         attr$dim <- NULL
251                                         candidate@info <- attr
252                                 } else if(is(candidate, "MxMatrix")) {
253                                         dimnames(value) <- dimnames(candidate)
254                                         candidate@values <- value
255                                 } else if (is(candidate, "MxExpectation")) {
256                                         for (sl in names(attributes(value))) {
257                                                 slot(candidate, sl) <- attr(value, sl)
258                                         }
259                                 }
260                                 model[[name]] <- candidate
261                         }
262                 }
263         }
264     if (length(model@submodels) > 0) {
265         model@submodels <- lapply(model@submodels, 
266             updateModelEntitiesTargetModel, entNames, values, modelNameMapping)
267     }
268         return(model)
269 }
270
271 updateModelEntitiesHelper <- function(entNames, values, model) {
272     modelNameMapping <- sapply(entNames, getModelNameString)
273     model <- updateModelEntitiesTargetModel(model, entNames, 
274                 values, modelNameMapping)
275     return(model)
276 }
277
278 imxLocateIndex <- function(model, name, referant) {
279         if (length(name) == 0) return(name)
280         if (is.na(name)) { return(as.integer(name)) }
281         mNames <- names(model@matrices)
282         aNames <- names(model@algebras)
283         fNames <- names(model@fitfunctions)
284         eNames <- names(model@expectations)
285         oNames <- names(model@computes)
286         dNames <- names(model@datasets)         
287         matrixNumber <- match(name, mNames)
288         algebraNumber <- match(name, append(aNames, fNames))
289         dataNumber <- match(name, dNames)
290         expectationNumber <- match(name, eNames)
291         computeNumber <- match(name, oNames)
292         if (is.na(matrixNumber) && is.na(algebraNumber) 
293                 && is.na(dataNumber) && is.na(expectationNumber) &&
294             is.na(computeNumber)) {
295                 msg <- paste("The reference", omxQuotes(name),
296                         "does not exist.  It is used by the named entity",
297                         omxQuotes(referant),".")
298                 stop(msg, call.=FALSE)
299         } else if (!is.na(matrixNumber)) {
300                 return(- matrixNumber)
301         } else if (!is.na(dataNumber)) {
302                 return(dataNumber - 1L)
303         } else if (!is.na(expectationNumber)) {
304                 return(expectationNumber - 1L)
305         } else if (!is.na(computeNumber)) {
306                 return(computeNumber - 1L)
307         } else {
308                 return(algebraNumber - 1L)
309         }
310 }
311
312 imxPreprocessModel <- function(model) {
313         model@matrices <- lapply(model@matrices, findSquareBrackets)
314         model@submodels <- lapply(model@submodels, imxPreprocessModel)
315         return(model)
316 }
317
318
319 imxCheckMatrices <- function(model) {
320         matrices <- model@matrices
321         lapply(matrices, imxVerifyMatrix)
322         submodels <- imxDependentModels(model)
323         lapply(submodels, imxCheckMatrices)
324 }