Document optimizers
[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              output = "list",
20              debug = "list",
21            "VIRTUAL"),
22          contains = "MxBaseNamed")
23
24 setClassUnion("MxCompute", c("NULL", "MxBaseCompute"))
25
26 setGeneric("convertForBackend",
27         function(.Object, flatModel, model) {
28                 return(standardGeneric("convertForBackend"))
29         })
30
31 setGeneric("updateFromBackend",
32         function(.Object, computes) {
33                 return(standardGeneric("updateFromBackend"))
34         })
35
36 setGeneric("assignId",
37         function(.Object, id) {
38                 return(standardGeneric("assignId"))
39         })
40
41 setMethod("assignId", signature("MxBaseCompute"),
42         function(.Object, id) {
43                 .Object@id <- id
44                 .Object
45         })
46
47 setGeneric("getFreeVarGroup",
48         function(.Object) {
49                 return(standardGeneric("getFreeVarGroup"))
50         })
51
52 setMethod("getFreeVarGroup", signature("MxBaseCompute"),
53         function(.Object) {
54                 list()
55         })
56
57 setMethod("updateFromBackend", signature("MxBaseCompute"),
58         function(.Object, computes) {
59                 if (length(computes)) {
60                         mystuff <- which(.Object@id == computes[seq(1,length(computes),2)])
61                         if (length(mystuff)) {
62                                 got <- computes[[2 * mystuff]]
63                                 for (sl in names(got)) {
64                                         slot(.Object, sl) <- got[[sl]]
65                                 }
66                         }
67                 }
68                 .Object
69         })
70
71 #----------------------------------------------------
72
73 setClass(Class = "MxComputeOperation",
74          contains = "MxBaseCompute",
75          representation = representation(
76            free.set = "MxOptionalChar"))
77
78 setMethod("qualifyNames", signature("MxComputeOperation"),
79         function(.Object, modelname, namespace) {
80                 .Object@name <- imxIdentifier(modelname, .Object@name)
81                 .Object@free.set <- imxConvertIdentifier(.Object@free.set, modelname, namespace)
82                 .Object
83         })
84
85 setMethod("getFreeVarGroup", signature("MxComputeOperation"),
86         function(.Object) {
87                 if (length(.Object@free.set)) {
88                         list(.Object@id, .Object@free.set)
89                 } else {
90                         list()
91                 }
92         })
93
94 setMethod("convertForBackend", signature("MxComputeOperation"),
95         function(.Object, flatModel, model) {
96                 name <- .Object@name
97                 .Object
98         })
99
100 #----------------------------------------------------
101
102 setClass(Class = "MxComputeOnce",
103          contains = "MxComputeOperation",
104          representation = representation(
105            what = "MxCharOrNumber",
106            verbose = "integer",
107            context = "character",
108            maxAbsChange = "logical",
109            fit = "logical",
110            gradient = "logical",
111            hessian = "logical",
112              information = "logical",
113              info.method = "MxOptionalChar",
114            ihessian = "logical",
115            hgprod = "logical"))
116
117 setMethod("qualifyNames", signature("MxComputeOnce"),
118         function(.Object, modelname, namespace) {
119                 .Object <- callNextMethod()
120                 for (sl in c('what')) {
121                         slot(.Object, sl) <- imxConvertIdentifier(slot(.Object, sl), modelname, namespace)
122                 }
123                 .Object
124         })
125
126 setMethod("convertForBackend", signature("MxComputeOnce"),
127         function(.Object, flatModel, model) {
128                 .Object <- callNextMethod()
129                 name <- .Object@name
130                 if (any(!is.integer(.Object@what))) {
131                         expNum <- match(.Object@what, names(flatModel@expectations))
132                         algNum <- match(.Object@what, append(names(flatModel@algebras),
133                                                              names(flatModel@fitfunctions)))
134                         if (any(is.na(expNum)) && any(is.na(algNum))) {
135                                 stop(paste("Can only apply MxComputeOnce to MxAlgebra or MxExpectation not",
136                                            deparse(.Object@what)))
137                         }
138                         if (!any(is.na(expNum))) {
139                                         # Usually negative numbers indicate matrices; not here
140                                 .Object@what <- - expNum
141                         } else {
142                                 if (any(algNum > length(flatModel@algebras)) && length(algNum) > 1) {
143                                         stop("MxComputeOnce cannot evaluate more than 1 fit function")
144                                 }
145                                 .Object@what <- algNum - 1L
146                         }
147                 }
148                 if (length(.Object@what) == 0) warning("MxComputeOnce with nothing will have no effect")
149                 if (all(.Object@what >= 0) && !.Object@maxAbsChange && !.Object@fit && !.Object@gradient &&
150                             !.Object@hessian && !.Object@information && !.Object@ihessian && !.Object@hgprod) {
151                         warning("MxComputeOnce with no action")
152                 }
153                 .Object
154         })
155
156 setMethod("initialize", "MxComputeOnce",
157           function(.Object, what, free.set, context, maxAbsChange, fit, gradient,
158                    hessian, information, info.method, ihessian, hgprod, verbose) {
159                   .Object@name <- 'compute'
160                   .Object@what <- what
161                   .Object@verbose = verbose
162                   .Object@free.set <- free.set
163                   .Object@context <- context
164                   .Object@maxAbsChange <- maxAbsChange
165                   .Object@fit <- fit
166                   .Object@gradient <- gradient
167                   .Object@hessian <- hessian
168                   .Object@information <- information
169                   .Object@info.method <- info.method
170                   .Object@ihessian <- ihessian
171                   .Object@hgprod <- hgprod
172                   .Object
173           })
174
175 ##' Compute something once
176 ##'
177 ##' The information matrix is only valid when parameters are at the
178 ##' maximum likelihood estimate. The information matrix is returned in
179 ##' model@output$hessian. You cannot request both the information
180 ##' matrix and the Hessian. The information matrix is invarient to the
181 ##' sign of the log likelihood scale whereas the Hessian is not.
182 ##'
183 ##' Some models are optimized for a sparse Hessian. Therefore, it can
184 ##' be much more efficient to compute the inverse Hessian in
185 ##' comparison to computing the Hessian and then inverting it.
186 ##' 
187 ##' @param what what to compute (a vector of expectation or algebra names)
188 ##' @param free.set names of matrices containing free variables (defaults to all free variables)
189 ##' @param context set the context for expectations
190 ##' @param maxAbsChange compute the maximum absolute change metric
191 ##' @param fit compute the fit statistic (often in log likelihood units)
192 ##' @param gradient compute the analytic gradient
193 ##' @param hessian compute the analytic Hessian
194 ##' @param information compute the analytic information matrix
195 ##' @param info.method which analytic method to use to compute the information
196 ##' matrix (one of "hessian", "sandwich", "bread", and "meat")
197 ##' @param ihessian compute the analytic inverse Hessian
198 ##' @param hgprod not implemented
199 ##' @param verbose the level of debugging output
200 ##' @aliases
201 ##' MxComputeOnce-class
202 ##' @examples
203 ##' data(demoOneFactor)
204 ##' factorModel <- mxModel(name ="One Factor",
205 ##'   mxMatrix(type="Full", nrow=5, ncol=1, free=TRUE, values=0.2, name="A"),
206 ##'     mxMatrix(type="Symm", nrow=1, ncol=1, free=FALSE, values=1, name="L"),
207 ##'     mxMatrix(type="Diag", nrow=5, ncol=5, free=TRUE, values=1, name="U"),
208 ##'     mxAlgebra(expression=A %*% L %*% t(A) + U, name="R"),
209 ##'     mxFitFunctionML(),mxExpectationNormal(covariance="R", dimnames=names(demoOneFactor)),
210 ##'     mxData(observed=cov(demoOneFactor), type="cov", numObs=500),
211 ##'     mxComputeOnce('fitfunction', fit=TRUE))
212 ##' factorModelFit <- mxRun(factorModel)
213 ##' factorModelFit@output$Minus2LogLikelihood  # 972.15
214
215 mxComputeOnce <- function(what, free.set=NULL, context=character(0),
216                           maxAbsChange=FALSE, fit=FALSE, gradient=FALSE,
217                           hessian=FALSE, information=FALSE, info.method=NULL,
218                           ihessian=FALSE, hgprod=FALSE, verbose=0L) {
219         new("MxComputeOnce", what, free.set, context, maxAbsChange, fit, gradient,
220             hessian, information, info.method, ihessian, hgprod, verbose)
221 }
222
223 #----------------------------------------------------
224
225 setClass(Class = "MxComputeGradientDescent",
226          contains = "MxComputeOperation",
227          representation = representation(
228            useGradient = "MxOptionalLogical",
229            fitfunction = "MxCharOrNumber",
230            engine = "character",
231            verbose = "integer"))
232
233 setMethod("qualifyNames", signature("MxComputeGradientDescent"),
234         function(.Object, modelname, namespace) {
235                 .Object <- callNextMethod()
236                 for (sl in c('fitfunction')) {
237                         slot(.Object, sl) <- imxConvertIdentifier(slot(.Object, sl), modelname, namespace)
238                 }
239                 .Object
240         })
241
242 setMethod("convertForBackend", signature("MxComputeGradientDescent"),
243         function(.Object, flatModel, model) {
244                 .Object <- callNextMethod()
245                 name <- .Object@name
246                 if (is.character(.Object@fitfunction)) {
247                         .Object@fitfunction <- imxLocateIndex(flatModel, .Object@fitfunction, name)
248                 }
249                 .Object
250         })
251
252 setMethod("initialize", "MxComputeGradientDescent",
253           function(.Object, free.set, engine, fit, useGradient, verbose) {
254                   .Object@name <- 'compute'
255                   .Object@free.set <- free.set
256                   .Object@fitfunction <- fit
257                   .Object@engine <- engine
258                   .Object@useGradient <- useGradient
259                   .Object@verbose <- verbose
260                   .Object
261           })
262
263 ##' imxHasNPSOL
264 ##'
265 ##' @return
266 ##' Returns TRUE if the NPSOL proprietary optimizer is compiled and
267 ##' linked with OpenMx. Otherwise FALSE.
268 imxHasNPSOL <- function() .Call(hasNPSOL_wrapper)
269
270 ##' Optimize parameters using a gradient descent optimizer
271 ##'
272 ##' This optimizer does not require analytic derivatives of the fit
273 ##' function. The open-source version of OpenMx only offers 1 choice,
274 ##' CSOLNP (based on Ye, 1988).  The proprietary version of OpenMx
275 ##' offers the choice of two optimizers, CSOLNP and NPSOL.
276 ##'
277 ##' @param free.set names of matrices containing free variables (defaults to all free variables)
278 ##' @param useGradient whether to use the analytic gradient (if available)
279 ##' @param engine specific NPSOL or CSOLNP
280 ##' @param fitfunction name of the fitfunction (defaults to 'fitfunction')
281 ##' @param verbose level of debugging output
282 ##' @aliases
283 ##' MxComputeGradientDescent-class
284 ##' @references Ye, Y. (1988). \emph{Interior algorithms for linear,
285 ##' quadratic, and linearly constrained convex programming.}
286 ##' (Unpublished doctoral dissertation.) Stanford University, CA.
287 ##' @examples
288 ##' data(demoOneFactor)
289 ##' factorModel <- mxModel(name ="One Factor",
290 ##'   mxMatrix(type="Full", nrow=5, ncol=1, free=FALSE, values=0.2, name="A"),
291 ##'     mxMatrix(type="Symm", nrow=1, ncol=1, free=FALSE, values=1, name="L"),
292 ##'     mxMatrix(type="Diag", nrow=5, ncol=5, free=TRUE, values=1, name="U"),
293 ##'     mxAlgebra(expression=A %*% L %*% t(A) + U, name="R"),
294 ##'   mxExpectationNormal(covariance="R", dimnames=names(demoOneFactor)),
295 ##'   mxFitFunctionML(),
296 ##'     mxData(observed=cov(demoOneFactor), type="cov", numObs=500),
297 ##'      mxComputeSequence(steps=list(
298 ##'      mxComputeGradientDescent(),
299 ##'      mxComputeEstimatedHessian(),
300 ##'      mxComputeStandardError(),
301 ##'      mxComputeHessianQuality()
302 ##'     )))
303 ##' factorModelFit <- mxRun(factorModel)
304 ##' factorModelFit@output$conditionNumber # 29.5
305
306 mxComputeGradientDescent <- function(free.set=NULL, useGradient=NULL,
307                                      engine=NULL, fitfunction='fitfunction', verbose=0L) {
308
309         if (missing(engine)) {
310                 engine <- options()$mxOptions[["Default optimizer"]]
311         }
312
313         new("MxComputeGradientDescent", free.set, engine, fitfunction, useGradient, verbose)
314 }
315
316 #----------------------------------------------------
317
318 setClass(Class = "MxComputeNewtonRaphson",
319          contains = "MxComputeOperation",
320          representation = representation(
321            fitfunction = "MxCharOrNumber",
322            maxIter = "integer",
323            tolerance = "numeric",
324            verbose = "integer",
325            carefully = "logical",
326              #output
327              iterations = "integer",
328              inform = "integer"))
329
330 setMethod("qualifyNames", signature("MxComputeNewtonRaphson"),
331         function(.Object, modelname, namespace) {
332                 .Object <- callNextMethod()
333                 for (sl in c('fitfunction')) {
334                         slot(.Object, sl) <- imxConvertIdentifier(slot(.Object, sl), modelname, namespace)
335                 }
336                 .Object
337         })
338
339 setMethod("convertForBackend", signature("MxComputeNewtonRaphson"),
340         function(.Object, flatModel, model) {
341                 .Object <- callNextMethod()
342                 name <- .Object@name
343                 if (is.character(.Object@fitfunction)) {
344                         .Object@fitfunction <- imxLocateIndex(flatModel, .Object@fitfunction, name)
345                 }
346                 .Object
347         })
348
349 setMethod("initialize", "MxComputeNewtonRaphson",
350           function(.Object, free.set, fit, maxIter, tolerance, verbose, carefully) {
351                   .Object@name <- 'compute'
352                   .Object@free.set <- free.set
353                   .Object@fitfunction <- fit
354                   .Object@maxIter <- maxIter
355                   .Object@tolerance <- tolerance
356                   .Object@verbose <- verbose
357                   .Object@carefully <- carefully
358                   .Object
359           })
360
361 ##' Optimize parameters using the Newton-Raphson algorithm
362 ##'
363 ##' This optimizer requires analytic 1st and 2nd derivatives of the
364 ##' fit function. Ramsay (1975) is used to speed convergence. Ramsay
365 ##' can be differentially applied to different groups of parameters.
366 ##' Comprehensive diagnostics are available by increasing the verbose
367 ##' level. The carefully option should not be needed if the
368 ##' derivatives are well behaved.
369 ##'
370 ##' @param free.set names of matrices containing free variables (defaults to all free variables)
371 ##' @param fitfunction name of the fitfunction (defaults to 'fitfunction')
372 ##' @param maxIter maximum number of iterations
373 ##' @param tolerance optimization is considered converged when maximum absolute change in parameters is less than tolerance
374 ##' @param verbose level of debugging output
375 ##' @param carefully whether to compute the fit statistic and enforce monotonicity (defaults to FALSE)
376 ##' @aliases
377 ##' MxComputeNewtonRaphson-class
378 ##' @references
379 ##' Luenberger, D. G. & Ye, Y. (2008). \emph{Linear and nonlinear programming.} Springer.
380 ##' 
381 ##' Ramsay, J. O. (1975). Solving implicit equations in psychometric data analysis.
382 ##' \emph{Psychometrika, 40}(3), 337-360.
383
384 mxComputeNewtonRaphson <- function(free.set=NULL, fitfunction='fitfunction', maxIter = 100L,
385                                    tolerance=1e-7, verbose=0L, carefully=FALSE)
386 {
387         new("MxComputeNewtonRaphson", free.set, fitfunction, maxIter, tolerance, verbose, carefully)
388 }
389
390 #----------------------------------------------------
391
392 setClass(Class = "MxComputeSteps",
393          contains = "MxBaseCompute",
394          representation = representation(
395            steps = "list"))
396
397 setMethod("getFreeVarGroup", signature("MxComputeSteps"),
398         function(.Object) {
399                 result <- list()
400                 for (step in .Object@steps) {
401                         got <- getFreeVarGroup(step)
402                         if (length(got)) result <- append(result, got)
403                 }
404                 result
405         })
406
407 setMethod("assignId", signature("MxComputeSteps"),
408         function(.Object, id) {
409                 steps <- .Object@steps
410                 for (sx in 1:length(steps)) {
411                         steps[[sx]] <- assignId(steps[[sx]], id)
412                         id <- steps[[sx]]@id + 1L
413                 }
414                 .Object@steps <- steps
415                 .Object@id <- id
416                 .Object
417         })
418
419 setMethod("qualifyNames", signature("MxComputeSteps"),
420         function(.Object, modelname, namespace) {
421                 .Object@name <- imxIdentifier(modelname, .Object@name)
422                 .Object@steps <- lapply(.Object@steps, function (c) qualifyNames(c, modelname, namespace))
423                 .Object
424         })
425
426 setMethod("convertForBackend", signature("MxComputeSteps"),
427         function(.Object, flatModel, model) {
428                 .Object@steps <- lapply(.Object@steps, function (c) convertForBackend(c, flatModel, model))
429                 .Object
430         })
431
432 setMethod("updateFromBackend", signature("MxComputeSteps"),
433         function(.Object, computes) {
434                 .Object <- callNextMethod()
435                 .Object@steps <- lapply(.Object@steps, function (c) updateFromBackend(c, computes))
436                 .Object
437         })
438
439 #----------------------------------------------------
440
441 setClass(Class = "MxComputeIterate",
442          contains = "MxComputeSteps",
443          representation = representation(
444            maxIter = "integer",
445            tolerance = "numeric",
446            verbose = "integer"))
447
448 setMethod("initialize", "MxComputeIterate",
449           function(.Object, steps, maxIter, tolerance, verbose) {
450                   .Object@name <- 'compute'
451                   .Object@steps <- steps
452                   .Object@maxIter <- maxIter
453                   .Object@tolerance <- tolerance
454                   .Object@verbose <- verbose
455                   .Object
456           })
457
458 ##' Repeatedly invoke a series of compute objects until change is less than tolerance
459 ##'
460 ##' One step (typically the last) must compute the fit or maxAbsChange.
461 ##'
462 ##' @param steps a list of compute objects
463 ##' @param maxIter the maximum number of iterations
464 ##' @param tolerance iterates until change is less than tolerance
465 ##' @param verbose level of debugging output
466 ##' @aliases
467 ##' MxComputeIterate-class
468 mxComputeIterate <- function(steps, maxIter=500L, tolerance=1e-4, verbose=0L) {
469         new("MxComputeIterate", steps=steps, maxIter=maxIter, tolerance=tolerance, verbose)
470 }
471
472 displayMxComputeIterate <- function(opt) {
473         cat(class(opt), omxQuotes(opt@name), '\n')
474         cat("@tolerance :", omxQuotes(opt@tolerance), '\n')
475         cat("@maxIter :", omxQuotes(opt@maxIter), '\n')
476         for (step in 1:length(opt@steps)) {
477                 cat("[[", step, "]] :", class(opt@steps[[step]]), '\n')
478         }
479         invisible(opt)
480 }
481
482 setMethod("print", "MxComputeIterate", function(x, ...) displayMxComputeIterate(x))
483 setMethod("show",  "MxComputeIterate", function(object) displayMxComputeIterate(object))
484
485 #----------------------------------------------------
486
487 setClass(Class = "MxComputeEM",
488          contains = "MxComputeOperation",
489          representation = representation(
490              what = "MxCharOrNumber",
491              mstep.fit = "MxCompute",
492              fit = "MxCompute",
493              maxIter = "integer",
494              tolerance = "numeric",
495              verbose = "integer",
496              ramsay="logical",
497              information="logical",
498              semMethod="numeric",
499              semDebug="logical",
500              noiseTarget="numeric",
501              noiseTolerance="numeric",
502              #output
503              semProbeCount="integer",
504              probeOffset="matrix",
505              semDiff="matrix",
506              paramHistLen="integer"))
507
508 setMethod("assignId", signature("MxComputeEM"),
509         function(.Object, id) {
510                 .Object@mstep.fit <- assignId(.Object@mstep.fit, id)
511                 .Object@fit <- assignId(.Object@fit, id + 1L)
512                 .Object@id <- id + 2L
513                 .Object
514         })
515
516 setMethod("getFreeVarGroup", signature("MxComputeEM"),
517         function(.Object) {
518                 result <- callNextMethod()
519                 for (step in c(.Object@mstep.fit, .Object@fit)) {
520                         got <- getFreeVarGroup(step)
521                         if (length(got)) result <- append(result, got)
522                 }
523                 result
524         })
525
526 setMethod("qualifyNames", signature("MxComputeEM"),
527         function(.Object, modelname, namespace) {
528                 .Object <- callNextMethod()
529                 for (sl in c('what')) {
530                         slot(.Object, sl) <- imxConvertIdentifier(slot(.Object, sl), modelname, namespace)
531                 }
532                 for (sl in c('mstep.fit', 'fit')) {
533                         slot(.Object, sl) <- qualifyNames(slot(.Object, sl), modelname, namespace)
534                 }
535                 .Object
536         })
537
538 setMethod("convertForBackend", signature("MxComputeEM"),
539         function(.Object, flatModel, model) {
540                 .Object <- callNextMethod()
541                 name <- .Object@name
542                 if (any(!is.integer(.Object@what))) {
543                         expNum <- match(.Object@what, names(flatModel@expectations))
544                         if (any(is.na(expNum))) {
545                                 stop(paste("Can only apply MxComputeEM to MxExpectation",
546                                            deparse(.Object@what)))
547                         }
548                         .Object@what <- expNum - 1L
549                 }
550                 if (length(.Object@what) == 0) warning("MxComputeEM with nothing will have no effect")
551                 for (sl in c('mstep.fit', 'fit')) {
552                         slot(.Object, sl) <- convertForBackend(slot(.Object, sl), flatModel, model)
553                 }
554                 .Object
555         })
556
557 setMethod("updateFromBackend", signature("MxComputeEM"),
558         function(.Object, computes) {
559                 .Object <- callNextMethod()
560                 .Object@mstep.fit <- updateFromBackend(.Object@mstep.fit, computes)
561                 .Object@fit <- updateFromBackend(.Object@fit, computes)
562                 .Object
563         })
564
565 setMethod("initialize", "MxComputeEM",
566           function(.Object, what, mstep.fit, fit, maxIter, tolerance,
567                    verbose, ramsay, information, noiseTarget, noiseTolerance, semDebug, semMethod) {
568                   .Object@name <- 'compute'
569                   .Object@what <- what
570                   .Object@mstep.fit <- mstep.fit
571                   .Object@fit <- fit
572                   .Object@maxIter <- maxIter
573                   .Object@tolerance <- tolerance
574                   .Object@verbose <- verbose
575                   .Object@ramsay <- ramsay
576                   .Object@information <- information
577                   .Object@noiseTarget <- noiseTarget
578                   .Object@noiseTolerance <- noiseTolerance
579                   .Object@semDebug <- semDebug
580                   .Object@semMethod <- semMethod
581                   .Object
582           })
583
584 mxComputeEM <- function(expectation, mstep.fit, fit, maxIter=500L, tolerance=1e-4,
585                         verbose=0L, ramsay=TRUE, information=FALSE, noiseTarget=1, noiseTolerance=5,
586                         semDebug=FALSE, semMethod=1.0) {
587         new("MxComputeEM", what=expectation, mstep.fit, fit, maxIter=maxIter,
588             tolerance=tolerance, verbose, ramsay, information, noiseTarget,
589             noiseTolerance, semDebug, semMethod)
590 }
591
592 displayMxComputeEM <- function(opt) {
593         cat(class(opt), omxQuotes(opt@name), '\n')
594         cat("@tolerance :", omxQuotes(opt@tolerance), '\n')
595         cat("@maxIter :", omxQuotes(opt@maxIter), '\n')
596         invisible(opt)
597 }
598
599 setMethod("print", "MxComputeEM", function(x, ...) displayMxComputeEM(x))
600 setMethod("show",  "MxComputeEM", function(object) displayMxComputeEM(object))
601
602 #----------------------------------------------------
603
604 setClass(Class = "MxComputeEstimatedHessian",
605          contains = "MxComputeOperation",
606          representation = representation(
607            fitfunction = "MxCharOrNumber",
608              parallel = "logical",
609              stepSize = "numeric",
610              iterations = "integer"))
611
612 setMethod("qualifyNames", signature("MxComputeEstimatedHessian"),
613         function(.Object, modelname, namespace) {
614                 .Object <- callNextMethod();
615                 .Object@fitfunction <- imxConvertIdentifier(.Object@fitfunction, modelname, namespace)
616                 .Object
617         })
618
619 setMethod("convertForBackend", signature("MxComputeEstimatedHessian"),
620         function(.Object, flatModel, model) {
621                 .Object <- callNextMethod();
622                 name <- .Object@name
623                 if (is.character(.Object@fitfunction)) {
624                         .Object@fitfunction <- imxLocateIndex(flatModel, .Object@fitfunction, name)
625                 }
626                 .Object
627         })
628
629 setMethod("initialize", "MxComputeEstimatedHessian",
630           function(.Object, free.set, fit, parallel, stepSize, iterations) {
631                   .Object@name <- 'compute'
632                   .Object@free.set <- free.set
633                   .Object@fitfunction <- fit
634                   .Object@parallel <- parallel
635                   .Object@stepSize <- stepSize
636                   .Object@iterations <- iterations
637                   .Object
638           })
639
640 ##' Numerically estimate Hessian using Richardson extrapolation
641 ##'
642 ##' For N free parameters, Richardson extrapolation requires
643 ##' (iterations * (N^2 + N)) function evaluations.
644 ##' 
645 ##' The implementation is closely based on the numDeriv R package.
646 ##' 
647 ##' @param free.set names of matrices containing free variables (defaults to all free variables)
648 ##' @param fitfunction name of the fitfunction (defaults to 'fitfunction')
649 ##' @param parallel whether to evaluate the fitfunction in parallel (defaults to TRUE)
650 ##' @param stepSize starting set size (defaults to 0.0001)
651 ##' @param iterations number of Richardson extrapolation iterations (defaults to 4L)
652 ##' @aliases
653 ##' MxComputeEstimatedHessian-class
654 ##' @examples
655 ##' data(demoOneFactor)
656 ##' factorModel <- mxModel(name ="One Factor",
657 ##'   mxMatrix(type="Full", nrow=5, ncol=1, free=FALSE, values=0.2, name="A"),
658 ##'     mxMatrix(type="Symm", nrow=1, ncol=1, free=FALSE, values=1, name="L"),
659 ##'     mxMatrix(type="Diag", nrow=5, ncol=5, free=TRUE, values=1, name="U"),
660 ##'     mxAlgebra(expression=A %*% L %*% t(A) + U, name="R"),
661 ##'   mxExpectationNormal(covariance="R", dimnames=names(demoOneFactor)),
662 ##'   mxFitFunctionML(),
663 ##'     mxData(observed=cov(demoOneFactor), type="cov", numObs=500),
664 ##'     mxComputeEstimatedHessian())
665 ##' factorModelFit <- mxRun(factorModel)
666 ##' factorModelFit@output$hessian
667
668 mxComputeEstimatedHessian <- function(free.set=NULL, fitfunction='fitfunction',
669                                       parallel=TRUE, stepSize=0.0001, iterations=4L) {
670         new("MxComputeEstimatedHessian", free.set, fitfunction, parallel, stepSize, iterations)
671 }
672
673 #----------------------------------------------------
674
675 setClass(Class = "MxComputeStandardError",
676          contains = "MxComputeOperation")
677
678 setMethod("initialize", "MxComputeStandardError",
679           function(.Object, free.set) {
680                   .Object@name <- 'compute'
681                   .Object@free.set <- free.set
682                   .Object
683           })
684
685 ##' Compute standard errors given the Hessian or inverse Hessian
686 ##'
687 ##' @param free.set names of matrices containing free variables (defaults to all free variables)
688 ##' @aliases
689 ##' MxComputeStandardError-class
690
691 mxComputeStandardError <- function(free.set=NULL) {
692         new("MxComputeStandardError", free.set)
693 }
694
695 #----------------------------------------------------
696
697 setClass(Class = "MxComputeHessianQuality",
698          contains = "MxComputeOperation")
699
700 setMethod("initialize", "MxComputeHessianQuality",
701           function(.Object, free.set) {
702                   .Object@name <- 'compute'
703                   .Object@free.set <- free.set
704                   .Object
705           })
706
707 ##' Compute the quality of the Hessian
708 ##'
709 ##' Tests whether the Hessian is positive definite
710 ##' (model@output$infoDefinite) and, if so, computes the condition
711 ##' number (model@output$conditionNumber). See Luenberger & Ye (2008)
712 ##' Second Order Test (p. 190) and Condition Number (p. 239).
713 ##' 
714 ##' @param free.set names of matrices containing free variables (defaults to all free variables)
715 ##' @aliases
716 ##' MxComputeHessianQuality-class
717 ##' @references
718 ##' Luenberger, D. G. & Ye, Y. (2008). Linear and nonlinear programming. Springer.
719
720 mxComputeHessianQuality <- function(free.set=NULL) {
721         new("MxComputeHessianQuality", free.set)
722 }
723
724 #----------------------------------------------------
725
726 setClass(Class = "MxComputeSequence",
727          contains = "MxComputeSteps")
728
729 setMethod("initialize", "MxComputeSequence",
730           function(.Object, steps) {
731                   .Object@name <- 'compute'
732                   .Object@steps <- steps
733                   .Object
734           })
735
736 ##' Invoke a series of compute objects in sequence
737 ##'
738 ##' @param steps a list of compute objects
739 ##' @aliases
740 ##' MxComputeSequence-class
741 mxComputeSequence <- function(steps) {
742         new("MxComputeSequence", steps=steps)
743 }
744
745 displayMxComputeSequence <- function(opt) {
746         cat(class(opt), omxQuotes(opt@name), '\n')
747         for (step in 1:length(opt@steps)) {
748                 cat("[[", step, "]] :", class(opt@steps[[step]]), '\n')
749         }
750         invisible(opt)
751 }
752
753 setMethod("print", "MxComputeSequence", function(x, ...) displayMxComputeSequence(x))
754 setMethod("show",  "MxComputeSequence", function(object) displayMxComputeSequence(object))
755
756 #----------------------------------------------------
757
758 displayMxComputeOperation <- function(opt) {
759         cat(class(opt), omxQuotes(opt@name), '\n')
760         cat("@id :", opt@id, '\n')
761         cat("@free.set :", omxQuotes(opt@free.set), '\n')
762         invisible(opt)
763 }
764
765 setMethod("print", "MxComputeOperation", function(x, ...) displayMxComputeOperation(x))
766 setMethod("show",  "MxComputeOperation", function(object) displayMxComputeOperation(object))
767
768 displayMxComputeGradientDescent <- function(opt) {
769         cat("@type :", omxQuotes(opt@type), '\n')
770         cat("@engine :", omxQuotes(opt@engine), '\n')
771         cat("@fitfunction :", omxQuotes(opt@fitfunction), '\n')
772         invisible(opt)
773 }
774
775 setMethod("print", "MxComputeGradientDescent",
776           function(x, ...) { callNextMethod(); displayMxComputeGradientDescent(x) })
777 setMethod("show",  "MxComputeGradientDescent",
778           function(object) { callNextMethod(); displayMxComputeGradientDescent(object) })
779
780 convertComputes <- function(flatModel, model) {
781         if (is.null(flatModel@compute)) return()
782         convertForBackend(flatModel@compute, flatModel, model)
783 }
784
785 updateModelCompute <- function(model, computes) {
786         if (is.null(model@compute)) return()
787         updateFromBackend(model@compute, computes)
788 }