Generic ComputeEM with Ramsay (1975) acceleration
[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            id = "integer",
19            "VIRTUAL"),
20          contains = "MxBaseNamed")
21
22 setClassUnion("MxCompute", c("NULL", "MxBaseCompute"))
23
24 setGeneric("convertForBackend",
25         function(.Object, flatModel, model) {
26                 return(standardGeneric("convertForBackend"))
27         })
28
29 setGeneric("assignId",
30         function(.Object, id) {
31                 return(standardGeneric("assignId"))
32         })
33
34 setMethod("assignId", signature("MxBaseCompute"),
35         function(.Object, id) {
36                 .Object@id <- id
37                 .Object
38         })
39
40 setGeneric("getFreeVarGroup",
41         function(.Object) {
42                 return(standardGeneric("getFreeVarGroup"))
43         })
44
45 setMethod("getFreeVarGroup", signature("MxBaseCompute"),
46         function(.Object) {
47                 list()
48         })
49
50 #----------------------------------------------------
51
52 setClass(Class = "MxComputeOperation",
53          contains = "MxBaseCompute",
54          representation = representation(
55            free.set = "MxOptionalChar"))
56
57 setMethod("qualifyNames", signature("MxComputeOperation"),
58         function(.Object, modelname, namespace) {
59                 .Object@name <- imxIdentifier(modelname, .Object@name)
60                 .Object@free.set <- imxConvertIdentifier(.Object@free.set, modelname, namespace)
61                 .Object
62         })
63
64 setMethod("getFreeVarGroup", signature("MxComputeOperation"),
65         function(.Object) {
66                 if (length(.Object@free.set)) {
67                         list(.Object@id, .Object@free.set)
68                 } else {
69                         list()
70                 }
71         })
72
73 setMethod("convertForBackend", signature("MxComputeOperation"),
74         function(.Object, flatModel, model) {
75                 name <- .Object@name
76                 .Object
77         })
78
79 #----------------------------------------------------
80
81 setClass(Class = "MxComputeOnce",
82          contains = "MxComputeOperation",
83          representation = representation(
84            what = "MxCharOrNumber",
85            verbose = "integer",
86            context = "character",
87            maxAbsChange = "logical",
88            fit = "logical",
89            gradient = "logical",
90            hessian = "logical",
91              information = "logical",
92            ihessian = "logical",
93            hgprod = "logical"))
94
95 setMethod("qualifyNames", signature("MxComputeOnce"),
96         function(.Object, modelname, namespace) {
97                 .Object <- callNextMethod();
98                 for (sl in c('what')) {
99                         slot(.Object, sl) <- imxConvertIdentifier(slot(.Object, sl), modelname, namespace)
100                 }
101                 .Object
102         })
103
104 setMethod("convertForBackend", signature("MxComputeOnce"),
105         function(.Object, flatModel, model) {
106                 .Object <- callNextMethod();
107                 name <- .Object@name
108                 if (any(!is.integer(.Object@what))) {
109                         expNum <- match(.Object@what, names(flatModel@expectations))
110                         algNum <- match(.Object@what, append(names(flatModel@algebras),
111                                                              names(flatModel@fitfunctions)))
112                         if (any(is.na(expNum)) && any(is.na(algNum))) {
113                                 stop(paste("Can only apply MxComputeOnce to MxAlgebra or MxExpectation not",
114                                            deparse(.Object@what)))
115                         }
116                         if (!any(is.na(expNum))) {
117                                         # Usually negative numbers indicate matrices; not here
118                                 .Object@what <- - expNum
119                         } else {
120                                 if (any(algNum > length(flatModel@algebras)) && length(algNum) > 1) {
121                                         stop("MxComputeOnce cannot evaluate more than 1 fit function")
122                                 }
123                                 .Object@what <- algNum - 1L
124                         }
125                 }
126                 if (length(.Object@what) == 0) warning("MxComputeOnce with nothing will have no effect")
127                 if (all(.Object@what >= 0) && !.Object@maxAbsChange && !.Object@fit && !.Object@gradient &&
128                             !.Object@hessian && !.Object@information && !.Object@ihessian && !.Object@hgprod) {
129                         warning("MxComputeOnce with no action")
130                 }
131                 .Object
132         })
133
134 setMethod("initialize", "MxComputeOnce",
135           function(.Object, what, free.set, context, maxAbsChange, fit, gradient,
136                    hessian, information, ihessian, hgprod, verbose) {
137                   .Object@name <- 'compute'
138                   .Object@what <- what
139                   .Object@verbose = verbose
140                   .Object@free.set <- free.set
141                   .Object@context <- context
142                   .Object@maxAbsChange <- maxAbsChange
143                   .Object@fit <- fit
144                   .Object@gradient <- gradient
145                   .Object@hessian <- hessian
146                   .Object@information <- information
147                   .Object@ihessian <- ihessian
148                   .Object@hgprod <- hgprod
149                   .Object
150           })
151
152 mxComputeOnce <- function(what, free.set=NULL, context=character(0),
153                           maxAbsChange=FALSE, fit=FALSE, gradient=FALSE,
154                           hessian=FALSE, information=FALSE, ihessian=FALSE, hgprod=FALSE, verbose=0L) {
155         new("MxComputeOnce", what, free.set, context, maxAbsChange, fit, gradient,
156             hessian, information, ihessian, hgprod, verbose)
157 }
158
159 #----------------------------------------------------
160
161 setClass(Class = "MxComputeGradientDescent",
162          contains = "MxComputeOperation",
163          representation = representation(
164            useGradient = "MxOptionalLogical",
165            fitfunction = "MxCharOrNumber",
166            engine = "character",
167            verbose = "integer"))
168
169 setMethod("qualifyNames", signature("MxComputeGradientDescent"),
170         function(.Object, modelname, namespace) {
171                 .Object <- callNextMethod();
172                 for (sl in c('fitfunction')) {
173                         slot(.Object, sl) <- imxConvertIdentifier(slot(.Object, sl), modelname, namespace)
174                 }
175                 .Object
176         })
177
178 setMethod("convertForBackend", signature("MxComputeGradientDescent"),
179         function(.Object, flatModel, model) {
180                 .Object <- callNextMethod();
181                 name <- .Object@name
182                 if (is.character(.Object@fitfunction)) {
183                         .Object@fitfunction <- imxLocateIndex(flatModel, .Object@fitfunction, name)
184                 }
185                 .Object
186         })
187
188 setMethod("initialize", "MxComputeGradientDescent",
189           function(.Object, free.set, engine, fit, useGradient, verbose) {
190                   .Object@name <- 'compute'
191                   .Object@free.set <- free.set
192                   .Object@fitfunction <- fit
193                   .Object@engine <- engine
194                   .Object@useGradient <- useGradient
195                   .Object@verbose <- verbose
196                   .Object
197           })
198
199 ##' imxHasNPSOL
200 ##'
201 ##' @return
202 ##' Returns TRUE if the NPSOL proprietary optimizer is compiled and
203 ##' linked with OpenMx. Otherwise FALSE.
204 imxHasNPSOL <- function() .Call(hasNPSOL_wrapper)
205
206 mxComputeGradientDescent <- function(type=NULL, free.set=NULL, useGradient=NULL,
207                                      engine=NULL, fitfunction='fitfunction', verbose=0L) {
208 # What to do with 'type'?
209 #       if (length(type) != 1) stop("Specific 1 compute type")
210
211         if (missing(engine)) {
212                 engine <- Sys.getenv("IMX_OPT_ENGINE")
213                 if (!nchar(engine)) {
214                         if (imxHasNPSOL()) {
215                                 engine <- "NPSOL"
216                         } else {
217                                 engine <- "CSOLNP"
218                         }
219                 }
220         }
221
222         new("MxComputeGradientDescent", free.set, engine, fitfunction, useGradient, verbose)
223 }
224
225 #----------------------------------------------------
226
227 setClass(Class = "MxComputeNewtonRaphson",
228          contains = "MxComputeOperation",
229          representation = representation(
230            fitfunction = "MxCharOrNumber",
231            maxIter = "integer",
232            tolerance = "numeric",
233            verbose = "integer",
234            carefully = "logical"))
235
236 setMethod("qualifyNames", signature("MxComputeNewtonRaphson"),
237         function(.Object, modelname, namespace) {
238                 .Object <- callNextMethod();
239                 for (sl in c('fitfunction')) {
240                         slot(.Object, sl) <- imxConvertIdentifier(slot(.Object, sl), modelname, namespace)
241                 }
242                 .Object
243         })
244
245 setMethod("convertForBackend", signature("MxComputeNewtonRaphson"),
246         function(.Object, flatModel, model) {
247                 .Object <- callNextMethod();
248                 name <- .Object@name
249                 if (is.character(.Object@fitfunction)) {
250                         .Object@fitfunction <- imxLocateIndex(flatModel, .Object@fitfunction, name)
251                 }
252                 .Object
253         })
254
255 setMethod("initialize", "MxComputeNewtonRaphson",
256           function(.Object, free.set, fit, maxIter, tolerance, verbose, carefully) {
257                   .Object@name <- 'compute'
258                   .Object@free.set <- free.set
259                   .Object@fitfunction <- fit
260                   .Object@maxIter <- maxIter
261                   .Object@tolerance <- tolerance
262                   .Object@verbose <- verbose
263                   .Object@carefully <- carefully
264                   .Object
265           })
266
267 mxComputeNewtonRaphson <- function(type, free.set=NULL,
268                                    fitfunction='fitfunction', maxIter = 100L, tolerance=1e-7,
269                                    verbose=0L, carefully=FALSE) {
270
271         new("MxComputeNewtonRaphson", free.set, fitfunction, maxIter, tolerance, verbose, carefully)
272 }
273
274 #----------------------------------------------------
275
276 setClass(Class = "MxComputeSteps",
277          contains = "MxBaseCompute",
278          representation = representation(
279            steps = "list"))
280
281 setMethod("getFreeVarGroup", signature("MxComputeSteps"),
282         function(.Object) {
283                 result <- list()
284                 for (step in .Object@steps) {
285                         got <- getFreeVarGroup(step)
286                         if (length(got)) result <- append(result, got)
287                 }
288                 result
289         })
290
291 setMethod("assignId", signature("MxComputeSteps"),
292         function(.Object, id) {
293                 steps <- .Object@steps
294                 for (sx in 1:length(steps)) {
295                         steps[[sx]] <- assignId(steps[[sx]], id)
296                         id <- steps[[sx]]@id + 1L
297                 }
298                 .Object@steps <- steps
299                 .Object@id <- id
300                 .Object
301         })
302
303 setMethod("qualifyNames", signature("MxComputeSteps"),
304         function(.Object, modelname, namespace) {
305                 .Object@name <- imxIdentifier(modelname, .Object@name)
306                 .Object@steps <- lapply(.Object@steps, function (c) qualifyNames(c, modelname, namespace))
307                 .Object
308         })
309
310 setMethod("convertForBackend", signature("MxComputeSteps"),
311         function(.Object, flatModel, model) {
312                 .Object@steps <- lapply(.Object@steps, function (c) convertForBackend(c, flatModel, model))
313                 .Object
314         })
315
316 #----------------------------------------------------
317
318 setClass(Class = "MxComputeIterate",
319          contains = "MxComputeSteps",
320          representation = representation(
321            maxIter = "integer",
322            tolerance = "numeric",
323            verbose = "integer"))
324
325 setMethod("initialize", "MxComputeIterate",
326           function(.Object, steps, maxIter, tolerance, verbose) {
327                   .Object@name <- 'compute'
328                   .Object@steps <- steps
329                   .Object@maxIter <- maxIter
330                   .Object@tolerance <- tolerance
331                   .Object@verbose <- verbose
332                   .Object
333           })
334
335 mxComputeIterate <- function(steps, maxIter=500L, tolerance=1e-4, verbose=0L) {
336         new("MxComputeIterate", steps=steps, maxIter=maxIter, tolerance=tolerance, verbose)
337 }
338
339 displayMxComputeIterate <- function(opt) {
340         cat(class(opt), omxQuotes(opt@name), '\n')
341         cat("@tolerance :", omxQuotes(opt@tolerance), '\n')
342         cat("@maxIter :", omxQuotes(opt@maxIter), '\n')
343         for (step in 1:length(opt@steps)) {
344                 cat("[[", step, "]] :", class(opt@steps[[step]]), '\n')
345         }
346         invisible(opt)
347 }
348
349 setMethod("print", "MxComputeIterate", function(x, ...) displayMxComputeIterate(x))
350 setMethod("show",  "MxComputeIterate", function(object) displayMxComputeIterate(object))
351
352 #----------------------------------------------------
353
354 setClass(Class = "MxComputeEM",
355          contains = "MxComputeOperation",
356          representation = representation(
357              what = "MxCharOrNumber",
358              mstep.fit = "MxCompute",
359              fit = "MxCompute",
360              maxIter = "integer",
361              tolerance = "numeric",
362              verbose = "integer",
363              ramsay="logical",
364              information="logical"))
365
366 setMethod("assignId", signature("MxComputeEM"),
367         function(.Object, id) {
368                 .Object@mstep.fit <- assignId(.Object@mstep.fit, id)
369                 .Object@fit <- assignId(.Object@fit, id + 1L)
370                 .Object@id <- id + 2L
371                 .Object
372         })
373
374 setMethod("getFreeVarGroup", signature("MxComputeEM"),
375         function(.Object) {
376                 result <- callNextMethod();
377                 for (step in c(.Object@mstep.fit, .Object@fit)) {
378                         got <- getFreeVarGroup(step)
379                         if (length(got)) result <- append(result, got)
380                 }
381                 result
382         })
383
384 setMethod("qualifyNames", signature("MxComputeEM"),
385         function(.Object, modelname, namespace) {
386                 .Object <- callNextMethod();
387                 for (sl in c('what')) {
388                         slot(.Object, sl) <- imxConvertIdentifier(slot(.Object, sl), modelname, namespace)
389                 }
390                 for (sl in c('mstep.fit', 'fit')) {
391                         slot(.Object, sl) <- qualifyNames(slot(.Object, sl), modelname, namespace)
392                 }
393                 .Object
394         })
395
396 setMethod("convertForBackend", signature("MxComputeEM"),
397         function(.Object, flatModel, model) {
398                 .Object <- callNextMethod();
399                 name <- .Object@name
400                 if (any(!is.integer(.Object@what))) {
401                         expNum <- match(.Object@what, names(flatModel@expectations))
402                         if (any(is.na(expNum))) {
403                                 stop(paste("Can only apply MxComputeEM to MxExpectation",
404                                            deparse(.Object@what)))
405                         }
406                         .Object@what <- expNum - 1L
407                 }
408                 if (length(.Object@what) == 0) warning("MxComputeEM with nothing will have no effect")
409                 for (sl in c('mstep.fit', 'fit')) {
410                         slot(.Object, sl) <- convertForBackend(slot(.Object, sl), flatModel, model)
411                 }
412                 .Object
413         })
414
415 setMethod("initialize", "MxComputeEM",
416           function(.Object, what, mstep.fit, fit, maxIter, tolerance, verbose, ramsay, information) {
417                   .Object@name <- 'compute'
418                   .Object@what <- what
419                   .Object@mstep.fit <- mstep.fit
420                   .Object@fit <- fit
421                   .Object@maxIter <- maxIter
422                   .Object@tolerance <- tolerance
423                   .Object@verbose <- verbose
424                   .Object@ramsay <- ramsay
425                   .Object@information <- information
426                   .Object
427           })
428
429 mxComputeEM <- function(expectation, mstep.fit, fit, maxIter=500L, tolerance=1e-4,
430                         verbose=0L, ramsay=TRUE, information=FALSE) {
431         new("MxComputeEM", what=expectation, mstep.fit, fit, maxIter=maxIter,
432                 tolerance=tolerance, verbose, ramsay, information)
433 }
434
435 displayMxComputeEM <- function(opt) {
436         cat(class(opt), omxQuotes(opt@name), '\n')
437         cat("@tolerance :", omxQuotes(opt@tolerance), '\n')
438         cat("@maxIter :", omxQuotes(opt@maxIter), '\n')
439         invisible(opt)
440 }
441
442 setMethod("print", "MxComputeEM", function(x, ...) displayMxComputeEM(x))
443 setMethod("show",  "MxComputeEM", function(object) displayMxComputeEM(object))
444
445 #----------------------------------------------------
446
447 setClass(Class = "MxComputeEstimatedHessian",
448          contains = "MxComputeOperation",
449          representation = representation(
450            fitfunction = "MxCharOrNumber",
451            se = "logical"))
452
453 setMethod("qualifyNames", signature("MxComputeEstimatedHessian"),
454         function(.Object, modelname, namespace) {
455                 .Object <- callNextMethod();
456                 .Object@fitfunction <- imxConvertIdentifier(.Object@fitfunction, modelname, namespace)
457                 .Object
458         })
459
460 setMethod("convertForBackend", signature("MxComputeEstimatedHessian"),
461         function(.Object, flatModel, model) {
462                 .Object <- callNextMethod();
463                 name <- .Object@name
464                 if (is.character(.Object@fitfunction)) {
465                         .Object@fitfunction <- imxLocateIndex(flatModel, .Object@fitfunction, name)
466                 }
467                 .Object
468         })
469
470 setMethod("initialize", "MxComputeEstimatedHessian",
471           function(.Object, free.set, fit, want.se) {
472                   .Object@name <- 'compute'
473                   .Object@free.set <- free.set
474                   .Object@fitfunction <- fit
475                   .Object@se <- want.se
476                   .Object
477           })
478
479 mxComputeEstimatedHessian <- function(free.set=NULL, fitfunction='fitfunction', want.se=TRUE) {
480         new("MxComputeEstimatedHessian", free.set, fitfunction, want.se)
481 }
482
483 #----------------------------------------------------
484
485 setClass(Class = "MxComputeSequence",
486          contains = "MxComputeSteps")
487
488 setMethod("initialize", "MxComputeSequence",
489           function(.Object, steps) {
490                   .Object@name <- 'compute'
491                   .Object@steps <- steps
492                   .Object
493           })
494
495 mxComputeSequence <- function(steps) {
496         new("MxComputeSequence", steps=steps)
497 }
498
499 displayMxComputeSequence <- function(opt) {
500         cat(class(opt), omxQuotes(opt@name), '\n')
501         for (step in 1:length(opt@steps)) {
502                 cat("[[", step, "]] :", class(opt@steps[[step]]), '\n')
503         }
504         invisible(opt)
505 }
506
507 setMethod("print", "MxComputeSequence", function(x, ...) displayMxComputeSequence(x))
508 setMethod("show",  "MxComputeSequence", function(object) displayMxComputeSequence(object))
509
510 #----------------------------------------------------
511
512 displayMxComputeOperation <- function(opt) {
513         cat(class(opt), omxQuotes(opt@name), '\n')
514         cat("@id :", opt@id, '\n')
515         cat("@free.set :", omxQuotes(opt@free.set), '\n')
516         invisible(opt)
517 }
518
519 setMethod("print", "MxComputeOperation", function(x, ...) displayMxComputeOperation(x))
520 setMethod("show",  "MxComputeOperation", function(object) displayMxComputeOperation(object))
521
522 displayMxComputeGradientDescent <- function(opt) {
523         cat("@type :", omxQuotes(opt@type), '\n')
524         cat("@engine :", omxQuotes(opt@engine), '\n')
525         cat("@fitfunction :", omxQuotes(opt@fitfunction), '\n')
526         invisible(opt)
527 }
528
529 setMethod("print", "MxComputeGradientDescent",
530           function(x, ...) { callNextMethod(); displayMxComputeGradientDescent(x) })
531 setMethod("show",  "MxComputeGradientDescent",
532           function(object) { callNextMethod(); displayMxComputeGradientDescent(object) })
533
534 convertComputes <- function(flatModel, model) {
535         retval <- lapply(flatModel@computes, function(opt) {
536                 convertForBackend(opt, flatModel, model)
537         })
538         retval
539 }