Teach MxExpectation to store its relationship in the model tree
[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) == 2) {
53                 retval <- list(pList[[1]], pList[[2]], integer())
54                 return(retval)
55         }
56         locations <- pList[3:length(pList)]
57         deps <- lapply(locations, findDependencies, flatModel, dependencies)
58         depnames <- Reduce(union, deps, character())
59         depnames <- Filter(Negate(isExpectation), depnames)
60         depnumbers <- sapply(depnames, doLocateIndex, flatModel, flatModel@name, USE.NAMES=FALSE)
61         depnumbers <- as.integer(depnumbers)
62         retval <- list(pList[[1]], pList[[2]], depnumbers)
63         return(append(retval, locations))
64 }
65
66 generateParameterList <- function(flatModel, dependencies) {
67         result <- list()
68         if (length(flatModel@matrices) == 0) {
69                 return(result)
70         }
71         for(i in 1:length(flatModel@matrices)) {
72                 matrix <- flatModel@matrices[[i]]
73                 result <- generateParameterListHelper(
74                         matrix, result, i - 1L)
75         }
76         result <- lapply(result, parameterDependencyList, flatModel, dependencies)
77         return(result)
78 }
79
80 generateDefinitionList <- function(flatModel, dependencies) {
81         result <- list()
82         defLocations <- generateDefinitionLocations(flatModel@datasets)
83         if (length(flatModel@matrices) == 0) {
84                 return(result)
85         }
86         for(i in 1:length(flatModel@matrices)) {
87                 result <- generateDefinitionListHelper(
88                         flatModel@matrices[[i]], 
89                         result, defLocations, i - 1L)
90         }
91         result <- lapply(result, parameterDependencyList, flatModel, dependencies)
92         return(result)
93 }
94
95 generateValueList <- function(mList, pList) {
96         mList <- lapply(mList, function(x) { x[[1]] })
97         retval <- vector()
98         if (length(pList) == 0) {
99                 return(retval)
100         }
101         for(i in 1:length(pList)) {
102                 parameter <- pList[[i]]
103                 parameter <- parameter[4:length(parameter)] # Remove (min, max, dependencies)
104                 if (length(parameter) > 1) {
105                         values <- sapply(parameter, generateValueHelper, mList)
106                         if (!all(values == values[[1]])) {
107                                 warning(paste('Parameter',i,'has multiple start values.',
108                                         'Selecting', values[[1]]))
109                         }
110                         retval[i] <- values[[1]]
111                 } else {
112                         retval[i] <- generateValueHelper(parameter[[1]], mList)
113                 }
114         }
115         return(retval)  
116 }
117
118 generateValueHelper <- function(triple, mList) {
119         mat <- triple[1] + 1
120         row <- triple[2] + 1
121         col <- triple[3] + 1
122         val <- mList[[mat]][row,col]
123         if (is.na(val)) {
124                 stop(paste("Starting value in ",names(mList)[[mat]],
125                            "[",row,",",col,"] is missing", sep=""))
126         }
127         return(val)
128 }
129
130 getFitFunctionIndex <- function(flatModel) {
131         fitfunction <- flatModel@fitfunction
132         if(is.null(fitfunction)) {
133                 return(NULL)
134         } else {
135                 return(imxLocateIndex(flatModel, fitfunction@name, flatModel@name))
136         }
137 }
138
139 imxUpdateModelValues <- function(model, flatModel, pList, values) {
140         if(length(pList) != length(values)) {
141                 stop(paste("This model has", length(pList), 
142                         "parameters, but you have given me", length(values),
143                         "values"))
144         }
145         if (length(pList) == 0) {
146                 return(model)
147         }
148         for(i in 1:length(pList)) {
149                 parameters <- pList[[i]]
150                 parameters <- parameters[4:length(parameters)] # Remove min, max, and dependencies
151                 model <- updateModelValuesHelper(parameters, values[[i]], flatModel, model)
152     }
153         return(model)
154 }
155
156 updateModelValuesHelper <- function(triples, value, flatModel, model) {
157         for(i in 1:length(triples)) {
158                 triple <- triples[[i]]
159                 mat <- triple[1] + 1
160                 row <- triple[2] + 1
161                 col <- triple[3] + 1
162                 name <- flatModel@matrices[[mat]]@name
163                 model[[name]]@values[row,col] <- value
164         }
165         return(model)
166 }
167
168 removeTail <- function(lst, tailSize) {
169     newEnd <- length(lst) - tailSize
170     if (newEnd == 0) {
171         return(list())
172     } else {
173         return(lst[1 : newEnd])
174     }
175 }
176
177 updateModelMatrices <- function(model, flatModel, values) {
178         mList <- names(flatModel@matrices)
179         if (length(mList) != length(values)) {
180                 stop(paste("This model has", length(mList), 
181                         "matrices, but the backend has returned", length(values),
182                         "values"))
183         }
184         if (length(values) == 0) {
185                 return(model)
186         }
187         model <- updateModelEntitiesHelper(mList, values, model)
188         return(model)
189 }
190
191
192 updateModelAlgebras <- function(model, flatModel, values) {
193         aNames <- names(flatModel@algebras)
194         oNames <- names(flatModel@fitfunctions)
195         aList <- append(aNames, oNames)
196         if(length(aList) != length(values)) {
197                 stop(paste("This model has", length(aList), 
198                         "algebras, but the backend has returned", length(values),
199                         "values"))
200         }
201         if (length(aList) == 0) {
202                 return(model)
203         }
204         model <- updateModelEntitiesHelper(aList, values, model)
205         return(model)
206 }
207
208
209 updateModelEntitiesTargetModel <- function(model, entNames, values, modelNameMapping) {
210     nextName <- model@name
211     selectEnt <- entNames[modelNameMapping == nextName]
212     selectVal <- values[modelNameMapping == nextName]
213     if (length(selectEnt) > 0) {
214                 for(i in 1:length(selectEnt)) {
215                         name <- selectEnt[[i]]
216                         candidate <- model[[name]]
217                         value <- selectVal[[i]]
218                         if (!is.null(candidate) && (length(value) > 0)
219                                 && !is.nan(value)) {
220                                 if (is(candidate, "MxAlgebra")) {
221                                         dimnames(value) <- dimnames(candidate)
222                                         candidate@result <- value
223                                 } else if (is(candidate,"MxFitFunction")) {
224                                         candidate <- fitFunctionReadAttributes(candidate, value)
225                                 } else if(is(candidate, "MxMatrix")) {
226                                         dimnames(value) <- dimnames(candidate)
227                                         candidate@values <- value
228                                 }
229                                 model[[name]] <- candidate
230                         }
231                 }
232         }
233     if (length(model@submodels) > 0) {
234         model@submodels <- lapply(model@submodels, 
235             updateModelEntitiesTargetModel, entNames, values, modelNameMapping)
236     }
237         return(model)
238 }
239
240 updateModelEntitiesHelper <- function(entNames, values, model) {
241     modelNameMapping <- sapply(entNames, getModelNameString)
242     model <- updateModelEntitiesTargetModel(model, entNames, 
243                 values, modelNameMapping)
244     return(model)
245 }
246
247 imxLocateIndex <- function(model, name, referant) {
248         if (length(name) == 0) return(name)
249         if (is.na(name)) { return(as.integer(name)) }
250         mNames <- names(model@matrices)
251         aNames <- names(model@algebras)
252         fNames <- names(model@fitfunctions)
253         eNames <- names(model@expectations)
254         dNames <- names(model@datasets)         
255         matrixNumber <- match(name, mNames)
256         algebraNumber <- match(name, append(aNames, fNames))
257         dataNumber <- match(name, dNames)
258         expectationNumber <- match(name, eNames)
259         if (is.na(matrixNumber) && is.na(algebraNumber) 
260                 && is.na(dataNumber) && is.na(expectationNumber)) {
261                 msg <- paste("The reference", omxQuotes(name),
262                         "does not exist.  It is used by the named entity",
263                         omxQuotes(referant),".")
264                 stop(msg, call.=FALSE)
265         } else if (!is.na(matrixNumber)) {
266                 return(- matrixNumber)
267         } else if (!is.na(dataNumber)) {
268                 return(dataNumber - 1L)
269         } else if (!is.na(expectationNumber)) {
270                 return(expectationNumber - 1L)
271         } else {
272                 return(algebraNumber - 1L)
273         }
274 }
275
276 imxPreprocessModel <- function(model) {
277         model@matrices <- lapply(model@matrices, findSquareBrackets)
278         model@submodels <- lapply(model@submodels, imxPreprocessModel)
279         return(model)
280 }
281
282
283 imxCheckMatrices <- function(model) {
284         matrices <- model@matrices
285         lapply(matrices, imxVerifyMatrix)
286         submodels <- imxDependentModels(model)
287         lapply(submodels, imxCheckMatrices)
288 }