Newton-Raphson
[openmx:openmx.git] / R / MxCompute.R
1 #
2 #   Copyright 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 setClass(Class = "MxBaseCompute", 
17          representation = representation(
18            "VIRTUAL"),
19          contains = "MxBaseNamed")
20
21 setClassUnion("MxCompute", c("NULL", "MxBaseCompute"))
22
23 setGeneric("convertForBackend",
24         function(.Object, flatModel, model) {
25                 return(standardGeneric("convertForBackend"))
26         })
27
28 #----------------------------------------------------
29
30 setClass(Class = "MxComputeOperation",
31          contains = "MxBaseCompute",
32          representation = representation(
33            free.group = "MxCharOrNumber"))
34
35 setMethod("qualifyNames", signature("MxComputeOperation"),
36         function(.Object, modelname, namespace) {
37                 .Object@name <- imxIdentifier(modelname, .Object@name)
38                 .Object
39         })
40
41 setMethod("convertForBackend", signature("MxComputeOperation"),
42         function(.Object, flatModel, model) {
43                 name <- .Object@name
44                 fg <- match(.Object@free.group, flatModel@freeGroupNames)
45                 if (is.na(fg)) {
46                         stop(paste("Cannot find free group", .Object@free.group,
47                                    "in list of free groups:",
48                                    omxQuotes(flatModel@freeGroupNames)))
49                 } else {
50                         .Object@free.group <- fg - 1L
51                 }
52                 .Object
53         })
54
55 #----------------------------------------------------
56
57 setClass(Class = "MxComputeOnce",
58          contains = "MxComputeOperation",
59          representation = representation(
60            what = "MxOptionalCharOrNumber",
61            context = "character",
62            gradient = "logical",
63            hessian = "logical"))
64
65 setMethod("qualifyNames", signature("MxComputeOnce"),
66         function(.Object, modelname, namespace) {
67                 .Object@name <- imxIdentifier(modelname, .Object@name)
68                 .Object@what <- imxConvertIdentifier(.Object@what, modelname, namespace)
69                 .Object
70         })
71
72 setMethod("convertForBackend", signature("MxComputeOnce"),
73         function(.Object, flatModel, model) {
74                 .Object <- callNextMethod();
75                 name <- .Object@name
76                 if (length(.Object@what) != 1) stop("Can only apply MxComputeOnce to one object")
77                 if (!is.integer(.Object@what)) {
78                         expNum <- match(.Object@what, names(flatModel@expectations))
79                         algNum <- match(.Object@what, append(names(flatModel@algebras),
80                                                              names(flatModel@fitfunctions)))
81                         if (is.na(expNum) && is.na(algNum)) {
82                                 stop("Can only apply MxComputeOnce to MxAlgebra or MxExpectation")
83                         }
84                         if (!is.na(expNum)) {
85                                 .Object@what <- - expNum  # usually negative numbers indicate matrices
86                         } else {
87                                 .Object@what <- algNum - 1L
88                         }
89                 }
90                 .Object
91         })
92
93 setMethod("initialize", "MxComputeOnce",
94           function(.Object, what, free.group, context, gradient, hessian) {
95                   .Object@name <- 'compute'
96                   .Object@what <- what
97                   .Object@free.group <- free.group
98                   .Object@context <- context
99                   .Object@gradient <- gradient
100                   .Object@hessian <- hessian
101                   .Object
102           })
103
104 mxComputeOnce <- function(what, free.group='default', context=character(0), gradient=FALSE, hessian=FALSE) {
105         new("MxComputeOnce", what, free.group, context, gradient, hessian)
106 }
107
108 #----------------------------------------------------
109
110 setClass(Class = "MxComputeGradientDescent",
111          contains = "MxComputeOperation",
112          representation = representation(
113            fitfunction = "MxCharOrNumber",
114            engine = "character"))
115
116 setMethod("qualifyNames", signature("MxComputeGradientDescent"),
117         function(.Object, modelname, namespace) {
118                 .Object@name <- imxIdentifier(modelname, .Object@name)
119                 .Object@fitfunction <- imxConvertIdentifier(.Object@fitfunction, modelname, namespace)
120                 .Object
121         })
122
123 setMethod("convertForBackend", signature("MxComputeGradientDescent"),
124         function(.Object, flatModel, model) {
125                 .Object <- callNextMethod();
126                 name <- .Object@name
127                 if (is.character(.Object@fitfunction)) {
128                         .Object@fitfunction <- imxLocateIndex(flatModel, .Object@fitfunction, name)
129                 }
130                 .Object
131         })
132
133 setMethod("initialize", "MxComputeGradientDescent",
134           function(.Object, free.group, engine, fit) {
135                   .Object@name <- 'compute'
136                   .Object@free.group <- free.group
137                   .Object@fitfunction <- fit
138                   .Object@engine <- engine
139                   .Object
140           })
141
142 mxComputeGradientDescent <- function(type, free.group='default',
143                                      engine=NULL, fitfunction='fitfunction') {
144 # What to do with 'type'?
145 #       if (length(type) != 1) stop("Specific 1 compute type")
146
147         if (is.null(engine)) engine <- as.character(NA)
148
149         new("MxComputeGradientDescent", free.group, engine, fitfunction)
150 }
151
152 #----------------------------------------------------
153
154 setClass(Class = "MxComputeNewtonRaphson",
155          contains = "MxComputeOperation",
156          representation = representation(
157            fitfunction = "MxCharOrNumber",
158            maxIter = "integer",
159            tolerance = "numeric"))
160
161 setMethod("qualifyNames", signature("MxComputeNewtonRaphson"),
162         function(.Object, modelname, namespace) {
163                 .Object@name <- imxIdentifier(modelname, .Object@name)
164                 .Object@fitfunction <- imxConvertIdentifier(.Object@fitfunction, modelname, namespace)
165                 .Object
166         })
167
168 setMethod("convertForBackend", signature("MxComputeNewtonRaphson"),
169         function(.Object, flatModel, model) {
170                 .Object <- callNextMethod();
171                 name <- .Object@name
172                 if (is.character(.Object@fitfunction)) {
173                         .Object@fitfunction <- imxLocateIndex(flatModel, .Object@fitfunction, name)
174                 }
175                 .Object
176         })
177
178 setMethod("initialize", "MxComputeNewtonRaphson",
179           function(.Object, free.group, fit, maxIter, tolerance) {
180                   .Object@name <- 'compute'
181                   .Object@free.group <- free.group
182                   .Object@fitfunction <- fit
183                   .Object@maxIter <- maxIter
184                   .Object@tolerance <- tolerance
185                   .Object
186           })
187
188 mxComputeNewtonRaphson <- function(type, free.group='default',
189                                    fitfunction='fitfunction', maxIter = 100L, tolerance=1e-7) {
190
191         new("MxComputeNewtonRaphson", free.group, fitfunction, maxIter, tolerance)
192 }
193
194 #----------------------------------------------------
195
196 setClass(Class = "MxComputeIterate",
197          contains = "MxBaseCompute",
198          representation = representation(
199            steps = "list",
200            maxIter = "integer",
201            tolerance = "numeric",
202            verbose = "logical"))
203
204 setMethod("initialize", "MxComputeIterate",
205           function(.Object, steps, maxIter, tolerance, verbose) {
206                   .Object@name <- 'compute'
207                   .Object@steps <- steps
208                   .Object@maxIter <- maxIter
209                   .Object@tolerance <- tolerance
210                   .Object@verbose <- verbose
211                   .Object
212           })
213
214 setMethod("qualifyNames", signature("MxComputeIterate"),
215         function(.Object, modelname, namespace) {
216                 .Object@name <- imxIdentifier(modelname, .Object@name)
217                 .Object@steps <- lapply(.Object@steps, function (c) qualifyNames(c, modelname, namespace))
218                 .Object
219         })
220
221 setMethod("convertForBackend", signature("MxComputeIterate"),
222         function(.Object, flatModel, model) {
223                 .Object@steps <- lapply(.Object@steps, function (c) convertForBackend(c, flatModel, model))
224                 .Object
225         })
226
227 mxComputeIterate <- function(steps, maxIter=500L, tolerance=1e-4, verbose=FALSE) {
228         new("MxComputeIterate", steps=steps, maxIter=maxIter, tolerance=tolerance, verbose)
229 }
230
231 displayMxComputeIterate <- function(opt) {
232         cat(class(opt), omxQuotes(opt@name), '\n')
233         cat("@tolerance :", omxQuotes(opt@tolerance), '\n')
234         cat("@maxIter :", omxQuotes(opt@maxIter), '\n')
235         for (step in 1:length(opt@steps)) {
236                 cat("[[", step, "]] :", class(opt@steps[[step]]), '\n')
237         }
238         invisible(opt)
239 }
240
241 setMethod("print", "MxComputeIterate", function(x, ...) displayMxComputeIterate(x))
242 setMethod("show",  "MxComputeIterate", function(object) displayMxComputeIterate(object))
243
244 #----------------------------------------------------
245
246 setClass(Class = "MxComputeEstimatedHessian",
247          contains = "MxComputeOperation",
248          representation = representation(
249            fitfunction = "MxCharOrNumber",
250            se = "logical"))
251
252 setMethod("qualifyNames", signature("MxComputeEstimatedHessian"),
253         function(.Object, modelname, namespace) {
254                 .Object@name <- imxIdentifier(modelname, .Object@name)
255                 .Object@fitfunction <- imxConvertIdentifier(.Object@fitfunction, modelname, namespace)
256                 .Object
257         })
258
259 setMethod("convertForBackend", signature("MxComputeEstimatedHessian"),
260         function(.Object, flatModel, model) {
261                 .Object <- callNextMethod();
262                 name <- .Object@name
263                 if (is.character(.Object@fitfunction)) {
264                         .Object@fitfunction <- imxLocateIndex(flatModel, .Object@fitfunction, name)
265                 }
266                 .Object
267         })
268
269 setMethod("initialize", "MxComputeEstimatedHessian",
270           function(.Object, free.group, fit, want.se) {
271                   .Object@name <- 'compute'
272                   .Object@free.group <- free.group
273                   .Object@fitfunction <- fit
274                   .Object@se <- want.se
275                   .Object
276           })
277
278 mxComputeEstimatedHessian <- function(free.group='default', fitfunction='fitfunction', want.se=TRUE) {
279         new("MxComputeEstimatedHessian", free.group, fitfunction, want.se)
280 }
281
282 #----------------------------------------------------
283
284 setClass(Class = "MxComputeSequence",
285          contains = "MxBaseCompute",
286          representation = representation(
287            steps = "list"))
288
289 setMethod("initialize", "MxComputeSequence",
290           function(.Object, steps) {
291                   .Object@name <- 'compute'
292                   .Object@steps <- steps
293                   .Object
294           })
295
296 setMethod("qualifyNames", signature("MxComputeSequence"),
297         function(.Object, modelname, namespace) {
298                 .Object@name <- imxIdentifier(modelname, .Object@name)
299                 .Object@steps <- lapply(.Object@steps, function (c) qualifyNames(c, modelname, namespace))
300                 .Object
301         })
302
303 setMethod("convertForBackend", signature("MxComputeSequence"),
304         function(.Object, flatModel, model) {
305                 .Object@steps <- lapply(.Object@steps, function (c) convertForBackend(c, flatModel, model))
306                 .Object
307         })
308
309 mxComputeSequence <- function(steps) {
310         new("MxComputeSequence", steps=steps)
311 }
312
313 displayMxComputeSequence <- function(opt) {
314         cat(class(opt), omxQuotes(opt@name), '\n')
315         for (step in 1:length(opt@steps)) {
316                 cat("[[", step, "]] :", class(opt@steps[[step]]), '\n')
317         }
318         invisible(opt)
319 }
320
321 setMethod("print", "MxComputeSequence", function(x, ...) displayMxComputeSequence(x))
322 setMethod("show",  "MxComputeSequence", function(object) displayMxComputeSequence(object))
323
324 #----------------------------------------------------
325
326 displayMxComputeOperation <- function(opt) {
327         cat(class(opt), omxQuotes(opt@name), '\n')
328         cat("@free.group :", omxQuotes(opt@free.group), '\n')
329         invisible(opt)
330 }
331
332 setMethod("print", "MxComputeOperation", function(x, ...) displayMxComputeOperation(x))
333 setMethod("show",  "MxComputeOperation", function(object) displayMxComputeOperation(object))
334
335 displayMxComputeGradientDescent <- function(opt) {
336         cat("@type :", omxQuotes(opt@type), '\n')
337         cat("@engine :", omxQuotes(opt@engine), '\n')
338         cat("@fitfunction :", omxQuotes(opt@fitfunction), '\n')
339         invisible(opt)
340 }
341
342 setMethod("print", "MxComputeGradientDescent",
343           function(x, ...) { callNextMethod(); displayMxComputeGradientDescent(x) })
344 setMethod("show",  "MxComputeGradientDescent",
345           function(object) { callNextMethod(); displayMxComputeGradientDescent(object) })
346
347 convertComputes <- function(flatModel, model) {
348         retval <- lapply(flatModel@computes, function(opt) {
349                 convertForBackend(opt, flatModel, model)
350         })
351         retval
352 }