R/MxCompute.R
In OpenMx: Extended Structural Equation Modelling

Documented in imxHasNPSOL imxSparseInvert mxComputeBootstrap mxComputeCheckpoint mxComputeConfidenceInterval mxComputeDefault mxComputeEM mxComputeGenerateData mxComputeGradientDescent mxComputeHessianQuality mxComputeIterate mxComputeJacobian mxComputeLoadContext mxComputeLoadData mxComputeLoadMatrix mxComputeLoop mxComputeNelderMead mxComputeNewtonRaphson mxComputeNothing mxComputeNumericDeriv mxComputeOnce mxComputePenaltySearch mxComputeReportDeriv mxComputeReportExpectation mxComputeSequence mxComputeSetOriginalStarts mxComputeSimAnnealing mxComputeStandardError mxComputeTryCatch mxComputeTryHard omxDefaultComputePlan omxHasDefaultComputePlan

#
#   Copyright 2013-2021 by the individuals mentioned in the source code history
#
#   Licensed under the Apache License, Version 2.0 (the "License");
#   you may not use this file except in compliance with the License.
#   You may obtain a copy of the License at
#
#        http://www.apache.org/licenses/LICENSE-2.0
#
#   Unless required by applicable law or agreed to in writing, software
#   distributed under the License is distributed on an "AS IS" BASIS,
#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
#   See the License for the specific language governing permissions and
#   limitations under the License.

##' BaseCompute
##'
##' This is an internal class and should not be used directly.
##'
##' @aliases
##' $,BaseCompute-method
##' $<-,BaseCompute-method
##' print,BaseCompute-method
##' show,BaseCompute-method
##' @seealso
##' \link{mxComputeEM}, \link{mxComputeGradientDescent}, \link{mxComputeHessianQuality},
##' \link{mxComputeIterate}, \link{mxComputeNewtonRaphson}, \link{mxComputeNumericDeriv}
##' @rdname BaseCompute-class
setClass(Class = "BaseCompute",
	 representation = representation(
	   id = "integer",
	     freeSet = "MxOptionalChar",
	     output = "list",
	     debug = "list",
	     .persist = "logical",
	   "VIRTUAL"),
	 contains = "MxBaseNamed")

##' @title MxCompute
##' @name MxCompute-class
##'
##' @description
##' This is an internal class and should not be used directly.
##'
##' @aliases
##' MxCompute
##' MxCompute-class
##' @rdname MxCompute-class
setClassUnion("MxCompute", c("NULL", "BaseCompute"))

setGeneric("displayCompute",
	   function(Ob, indent) {
		   return(standardGeneric("displayCompute"))
	   })

setMethod("displayCompute", signature(Ob="BaseCompute", indent="integer"),
	  function(Ob, indent) {
		  sp <- paste(rep('  ', indent), collapse="")
		  cat(sp, class(Ob), omxQuotes(Ob@name), '\n')
#		  cat(sp, "$id :", Ob@id, '\n')   # only of interest to developers and introduces visual noise
		  cat(sp, "$freeSet :", omxQuotes(Ob@freeSet), '\n')
		  if (length(Ob$output)) {
			  for (elem in names(Ob$output)) {
				  stuff <- Ob@output[[elem]]
				  if (is.list(stuff)) {
					  cat(sp, "$output[[", omxQuotes(elem), "]] : ...", '\n')
				  } else {
					  cat(sp, "$output[[", omxQuotes(elem), "]] :", stuff, '\n')
				  }
			  }
		  }
		  if (length(Ob$debug)) {
			  for (elem in names(Ob$debug)) {
				  cat(sp, "$debug[[", omxQuotes(elem), "]] : ...", '\n')
			  }
		  }
		  invisible(Ob)
	  })

setMethod("print", "BaseCompute", function(x, ...) displayCompute(x, 1L))
setMethod("show",  "BaseCompute", function(object) displayCompute(object, 1L))

setGeneric("convertForBackend",
	function(.Object, flatModel, model) {
		return(standardGeneric("convertForBackend"))
	})

setMethod("convertForBackend", signature("BaseCompute"),
       function(.Object, flatModel, model) { .Object })

setGeneric("updateFromBackend",
	function(.Object, computes) {
		return(standardGeneric("updateFromBackend"))
	})

setGeneric("assignId",
	function(.Object, id, defaultFreeSet) {
		return(standardGeneric("assignId"))
	})

setMethod("assignId", signature("BaseCompute"),
	function(.Object, id, defaultFreeSet) {
		.Object@id <- id
		if (length(.Object@freeSet) == 1 && is.na(.Object@freeSet)) .Object@freeSet <- defaultFreeSet
		.Object
	})

setGeneric("getFreeVarGroup",
	function(.Object) {
		return(standardGeneric("getFreeVarGroup"))
	})

setMethod("getFreeVarGroup", signature("BaseCompute"),
	function(.Object) {
		if (length(.Object@freeSet) == 0 || (length(.Object@freeSet) == 1 && .Object@freeSet == '.')) {
			# none or all variables
		} else {
			list(.Object@id, .Object@freeSet)
		}
	})

setMethod("qualifyNames", signature("BaseCompute"),
	function(.Object, modelname, namespace) {
		.Object@name <- imxIdentifier(modelname, .Object@name)
		.Object@freeSet <- imxConvertIdentifier(.Object@freeSet, modelname, namespace)
		.Object
	})

setMethod("updateFromBackend", signature("BaseCompute"),
	function(.Object, computes) {
		if (length(computes)) {
			mystuff <- which(.Object@id == computes[seq(1,length(computes),2)])
			if (length(mystuff)) {
				got <- computes[[2 * mystuff]]
				for (sl in names(got)) {
					slot(.Object, sl) <- got[[sl]]
				}
			}
		}
		.Object
	})

setMethod("$", "BaseCompute", imxExtractSlot)

setReplaceMethod("$", "BaseCompute",
	function(x, name, value) {
		return(imxReplaceSlot(x, name, value, check=TRUE))
	}
)

setMethod("names", "BaseCompute", slotNames)

#----------------------------------------------------

setClass(Class = "MxComputeOnce",
	 contains = "BaseCompute",
	 representation = representation(
	     from = "MxCharOrNumber",
	     what = "MxOptionalChar",
	     how = "MxOptionalChar",
	     verbose = "integer",
	     .is.bestfit="logical"))

setMethod("qualifyNames", signature("MxComputeOnce"),
	function(.Object, modelname, namespace) {
		.Object <- callNextMethod()
		for (sl in c('from')) {
			slot(.Object, sl) <- imxConvertIdentifier(slot(.Object, sl), modelname, namespace)
		}
		.Object
	})

setMethod("convertForBackend", signature("MxComputeOnce"),
	function(.Object, flatModel, model) {
		name <- .Object@name
		if (any(!is.integer(.Object@from))) {
			expNum <- match(.Object@from, names(flatModel@expectations))
			algNum <- match(.Object@from, append(names(flatModel@algebras),
							     names(flatModel@fitfunctions)))
			if (any(is.na(expNum)) && any(is.na(algNum))) {
				stop(paste("Can only apply MxComputeOnce to MxFitFunction or MxExpectation not",
					   deparse(.Object@from)))
			}
			if (!any(is.na(expNum))) {
					# Usually negative numbers indicate matrices; not here
				.Object@from <- - expNum
			} else {
				.Object@from <- algNum - 1L
			}
		}
		.Object
	})

setMethod("initialize", "MxComputeOnce",
	  function(.Object, from, what, how, freeSet, verbose, .is.bestfit) {
		  .Object@name <- 'compute'
		  .Object@.persist <- TRUE
		  .Object@from <- from
		  .Object@what <- what
		  .Object@how <- how
		  .Object@freeSet <- freeSet
		  .Object@verbose = verbose
		  .Object@.is.bestfit <- .is.bestfit
		  .Object
	  })

##' Compute something once
##'
##' Some models are optimized for a sparse Hessian. Therefore, it can
##' be much more efficient to compute the inverse Hessian in
##' comparison to computing the Hessian and then inverting it.
##'
##' The information matrix is only valid when parameters are at the
##' maximum likelihood estimate. The information matrix is returned in
##' model$output$hessian. You cannot request both the information
##' matrix and the Hessian. The information matrix is invariant to the
##' sign of the log likelihood scale whereas the Hessian is not.
##' Use the \code{how} parameter to specify which approximation to use
##' (one of "default", "hessian", "sandwich", "bread", and "meat").
##'
##' @param from the object to perform the computation (a vector of expectation or fit function names)
##' @param what what to compute
##' @param how to compute it (optional)
##' @param ...  Not used.  Forces remaining arguments to be specified by name.
##' @param freeSet names of matrices containing free variables
##' @template args-verbose
##' @param .is.bestfit do not use; for backward compatibility
##' @aliases
##' MxComputeOnce-class
##' @examples
##' data(demoOneFactor)
##' factorModel <- mxModel(name ="One Factor",
##'   mxMatrix(type="Full", nrow=5, ncol=1, free=TRUE, values=0.2, name="A"),
##'     mxMatrix(type="Symm", nrow=1, ncol=1, free=FALSE, values=1, name="L"),
##'     mxMatrix(type="Diag", nrow=5, ncol=5, free=TRUE, values=1, name="U"),
##'     mxAlgebra(expression=A %*% L %*% t(A) + U, name="R"),
##'     mxFitFunctionML(),mxExpectationNormal(covariance="R", dimnames=names(demoOneFactor)),
##'     mxData(observed=cov(demoOneFactor), type="cov", numObs=500),
##'     mxComputeOnce('fitfunction', 'fit'))
##' factorModelFit <- mxRun(factorModel)
##' factorModelFit$output$fit  # 972.15

mxComputeOnce <- function(from, what=NULL, how=NULL, ...,
			  freeSet=NA_character_, verbose=0L, .is.bestfit=FALSE) {
  prohibitDotdotdot(list(...))
	if (length(from) == 0) warning("mxComputeOnce from nothing will have no effect")
	verbose <- as.integer(verbose)
	new("MxComputeOnce", from, what, how, freeSet, verbose, .is.bestfit)
}

setMethod("displayCompute", signature(Ob="MxComputeOnce", indent="integer"),
	  function(Ob, indent) {
		  callNextMethod();
		  sp <- paste(rep('  ', indent), collapse="")
		  for (sl in c("from", "what", "how", "verbose")) {
			  slname <- paste("$", sl, sep="")
			  if (is.null(slot(Ob, sl))) next
			  if (is.character(slot(Ob, sl))) {
				  cat(sp, slname, ":", omxQuotes(slot(Ob, sl)), '\n')
			  } else {
				  cat(sp, slname, ":", slot(Ob, sl), '\n')
			  }
		  }
		  invisible(Ob)
	  })

#----------------------------------------------------

setClass(Class = "MxComputeGradientDescent",
	 contains = "BaseCompute",
	 representation = representation(
	   fitfunction = "MxCharOrNumber",
	   engine = "character",
	     availableEngines = "character",
	     tolerance = "numeric",
	   nudgeZeroStarts = "MxCharOrLogical",
	   verbose = "integer",
	     maxMajorIter = "integer",
	   defaultCImethod = "character",
	     warmStart = "MxOptionalMatrix"))  # rename to 'preconditioner'?

setMethod("qualifyNames", signature("MxComputeGradientDescent"),
	function(.Object, modelname, namespace) {
		.Object <- callNextMethod()
		for (sl in c('fitfunction')) {
			slot(.Object, sl) <- imxConvertIdentifier(slot(.Object, sl), modelname, namespace)
		}
		.Object
	})

setMethod("convertForBackend", signature("MxComputeGradientDescent"),
	function(.Object, flatModel, model) {
		name <- .Object@name
		if (is.character(.Object@fitfunction)) {
			.Object@fitfunction <- imxLocateIndex(flatModel, .Object@fitfunction, .Object)
		}
		.Object
	})

setMethod("initialize", "MxComputeGradientDescent",
	  function(.Object, freeSet, engine, fit, verbose, tolerance, warmStart,
		   nudgeZeroStarts, maxMajorIter) {
		  .Object@name <- 'compute'
		  .Object@.persist <- TRUE
		  .Object@freeSet <- freeSet
		  .Object@fitfunction <- fit
		  .Object@engine <- engine
		  .Object@defaultCImethod <- 'none'
		  if (engine == 'SLSQP') .Object@defaultCImethod <- 'ineq'
		  .Object@verbose <- verbose
		  .Object@tolerance <- tolerance
		  .Object@warmStart <- warmStart
		  .Object@nudgeZeroStarts <- nudgeZeroStarts
		  .Object@maxMajorIter <- maxMajorIter
		  .Object@availableEngines <- c("CSOLNP", "SLSQP")
		  if (imxHasNPSOL()) {
			  .Object@availableEngines <- c(.Object@availableEngines, "NPSOL")
		  }
		  .Object
	  })

##' imxHasNPSOL
##'
##' @return
##' Returns TRUE if the NPSOL proprietary optimizer is compiled and
##' linked with OpenMx. Otherwise FALSE.
imxHasNPSOL <- function() .Call(hasNPSOL_wrapper)

##' Optimize parameters using a gradient descent optimizer
##'
##' This optimizer does not require analytic derivatives of the fit
##' function. The fully open-source CRAN version of OpenMx offers 2 choices,
##' CSOLNP and SLSQP (from the NLOPT collection).  The OpenMx Team's version of
##' OpenMx offers the choice of three optimizers: CSOLNP, SLSQP, and NPSOL.
##'
##' All three optimizers can use analytic gradients, and only NPSOL
##' uses \code{warmStart}. To customize more options, see
##' \link{mxOption}.
##'
##' @param freeSet names of matrices containing free parameters.
##' @param ...  Not used.  Forces remaining arguments to be specified by name.
##' @param engine specific 'CSOLNP', 'SLSQP', or 'NPSOL'
##' @param fitfunction name of the fitfunction (defaults to 'fitfunction')
##' @template args-verbose
##' @param tolerance how close to the optimum is close enough (also known as the optimality tolerance)
##' @param useGradient \lifecycle{soft-deprecated}
##' @param warmStart a Cholesky factored Hessian to use as the NPSOL Hessian starting value (preconditioner)
##' @param nudgeZeroStarts whether to nudge any zero starting values prior to optimization (default TRUE)
##' @param maxMajorIter maximum number of major iterations
##' @param gradientAlgo \lifecycle{soft-deprecated}
##' @param gradientIterations \lifecycle{soft-deprecated}
##' @param gradientStepSize \lifecycle{soft-deprecated}
##' @aliases
##' MxComputeGradientDescent-class
##' @references
##' Luenberger, D. G. & Ye, Y. (2008). \emph{Linear and nonlinear programming.} Springer.
##' @examples
##' data(demoOneFactor)
##' factorModel <- mxModel(name ="One Factor",
##'   mxMatrix(type="Full", nrow=5, ncol=1, free=FALSE, values=0.2, name="A"),
##'     mxMatrix(type="Symm", nrow=1, ncol=1, free=FALSE, values=1, name="L"),
##'     mxMatrix(type="Diag", nrow=5, ncol=5, free=TRUE, values=1, name="U"),
##'     mxAlgebra(expression=A %*% L %*% t(A) + U, name="R"),
##'   mxExpectationNormal(covariance="R", dimnames=names(demoOneFactor)),
##'   mxFitFunctionML(),
##'     mxData(observed=cov(demoOneFactor), type="cov", numObs=500),
##'      mxComputeSequence(steps=list(
##'      mxComputeGradientDescent(),
##'      mxComputeNumericDeriv(),
##'      mxComputeStandardError(),
##'      mxComputeHessianQuality()
##'     )))
##' factorModelFit <- mxRun(factorModel)
##' factorModelFit$output$conditionNumber # 29.5

mxComputeGradientDescent <- function(freeSet=NA_character_, ...,
				     engine=NULL, fitfunction='fitfunction', verbose=0L,
				     tolerance=NA_real_, useGradient=deprecated(), warmStart=NULL,
				     nudgeZeroStarts=mxOption(NULL,"Nudge zero starts"), maxMajorIter=NULL,
				     gradientAlgo=deprecated(),
				     gradientIterations=deprecated(),
				     gradientStepSize=deprecated()) {

  prohibitDotdotdot(list(...))
	if (missing(engine)) {
		engine <- options()$mxOptions[["Default optimizer"]]
	}
  if (lifecycle::is_present(useGradient)) {
    deprecate_soft("2.18.2", "OpenMx::mxComputeGradientDescent(useGradient = )",
                   details="For similar functionality, see mxOption(key='Analytic Gradients')")
  }
  if (lifecycle::is_present(gradientAlgo)) {
    deprecate_soft("2.18.2", "OpenMx::mxComputeGradientDescent(gradientAlgo = )",
                   details="For similar functionality, see mxOption(key='Gradient algorithm')")
  }
  if (lifecycle::is_present(gradientIterations)) {
    deprecate_soft("2.18.2", "OpenMx::mxComputeGradientDescent(gradientIterations = )",
                   details="For similar functionality, see mxOption(key='Gradient iterations')")
  }
  if (lifecycle::is_present(gradientStepSize)) {
    deprecate_soft("2.18.2", "OpenMx::mxComputeGradientDescent(gradientStepSize = )",
                   details="For similar functionality, see mxOption(key='Gradient step size')")
  }

	if (!is.null(warmStart) && engine != "NPSOL") {
		stop("Only NPSOL supports warmStart")
	}
	verbose <- as.integer(verbose)
	maxMajorIter <- as.integer(maxMajorIter)

	new("MxComputeGradientDescent", freeSet, engine, fitfunction, verbose,
	    tolerance, warmStart, nudgeZeroStarts, maxMajorIter)
}

setMethod("displayCompute", signature(Ob="MxComputeGradientDescent", indent="integer"),
	  function(Ob, indent) {
		  callNextMethod();
		  sp <- paste(rep('  ', indent), collapse="")
		  for (sl in c("engine", "fitfunction", "verbose", "tolerance",
			       "nudgeZeroStarts", "maxMajorIter")) {
			  val <- slot(Ob, sl)
			  if (length(val)==0 || is.na(val)) next
			  slname <- paste("$", sl, sep="")
			  if (is.character(slot(Ob, sl))) {
				  cat(sp, slname, ":", omxQuotes(slot(Ob, sl)), '\n')
			  } else {
				  cat(sp, slname, ":", slot(Ob, sl), '\n')
			  }
		  }
		  invisible(Ob)
	  })

#----------------------------------------------------

setClass(Class = "MxComputeTryHard",
	 contains = "BaseCompute",
	 representation = representation(
	     plan = "MxCompute",
	     verbose = "integer",
	     location = "numeric",
	     scale = "numeric",
	     maxRetries = "integer"))

setMethod("assignId", signature("MxComputeTryHard"),
	function(.Object, id, defaultFreeSet) {
		.Object <- callNextMethod()
		defaultFreeSet <- .Object@freeSet
		id <- .Object@id
		for (sl in c('plan')) {
			slot(.Object, sl) <- assignId(slot(.Object, sl), id, defaultFreeSet)
			id <- slot(.Object, sl)@id + 1L
		}
		.Object@id <- id
		.Object
	})

setMethod("getFreeVarGroup", signature("MxComputeTryHard"),
	function(.Object) {
		result <- callNextMethod()
		for (step in c(.Object@plan)) {
			got <- getFreeVarGroup(step)
			if (length(got)) result <- append(result, got)
		}
		result
	})

setMethod("qualifyNames", signature("MxComputeTryHard"),
	function(.Object, modelname, namespace) {
		.Object <- callNextMethod()
		for (sl in c('plan')) {
			slot(.Object, sl) <- qualifyNames(slot(.Object, sl), modelname, namespace)
		}
		.Object
	})

setMethod("convertForBackend", signature("MxComputeTryHard"),
	function(.Object, flatModel, model) {
		name <- .Object@name
		for (sl in c('plan')) {
			slot(.Object, sl) <- convertForBackend(slot(.Object, sl), flatModel, model)
		}
		.Object
	})

setMethod("updateFromBackend", signature("MxComputeTryHard"),
	function(.Object, computes) {
		.Object <- callNextMethod()
		for (sl in c('plan')) {
			slot(.Object, sl) <- updateFromBackend(slot(.Object, sl), computes)
		}
		.Object
	})

setMethod("initialize", "MxComputeTryHard",
	  function(.Object, freeSet, plan, verbose, location, scale, maxRetries) {
		  .Object@name <- 'compute'
		  .Object@.persist <- TRUE
		  .Object@freeSet <- freeSet
		  .Object@plan <- plan
		  .Object@verbose <- verbose
		  .Object@location <- location
		  .Object@scale <- scale
		  .Object@maxRetries <- maxRetries
		  .Object
	  })

##' Repeatedly attempt a compute plan until successful
##'
##' The provided compute plan is run until the status code indicates
##' success (0 or 1). It gives up after a small number of retries.
##'
##' Upon failure, start values are randomly perturbed.  Currently only
##' the uniform distribution is implemented.  The distribution is
##' parameterized by arguments \code{location} and \code{scale}.  The
##' location parameter is the distribution's median.  For the uniform
##' distribution, \code{scale} is the absolute difference between its
##' median and extrema (i.e., half the width of the rectangle).  Each
##' start value is multiplied by a random draw and then added to a
##' random draw from a distribution with the same \code{scale} but
##' with a median of zero.
##'
##' @param plan compute plan to optimize the model
##' @param ...  Not used.  Forces remaining arguments to be specified by name.
##' @param freeSet names of matrices containing free variables
##' @template args-verbose
##' @param location location of the perturbation distribution
##' @param scale scale of the perturbation distribution
##' @param maxRetries maximum number of plan evaluations per invocation (including the first evaluation)
##' @seealso
##' \code{\link{mxTryHard}}
##' @aliases
##' MxComputeTryHard-class
##' @references
##' Shanno, D. F. (1985). On Broyden-Fletcher-Goldfarb-Shanno method. \emph{Journal of
##' Optimization Theory and Applications, 46}(1), 87-94.
mxComputeTryHard <- function(plan, ..., freeSet=NA_character_, verbose=0L,
			     location=1.0, scale=0.25, maxRetries=3L)
{
  prohibitDotdotdot(list(...))
	verbose <- as.integer(verbose)
	maxRetries <- as.integer(maxRetries)
	new("MxComputeTryHard", freeSet, plan, verbose, location, scale, maxRetries)
}

setMethod("displayCompute", signature(Ob="MxComputeTryHard", indent="integer"),
	  function(Ob, indent) {
		  callNextMethod()
		  sp <- paste(rep('  ', indent), collapse="")
		  cat(sp, "$plan :", '\n')
		  displayCompute(Ob@plan, indent+1L)
		  for (sl in c("verbose","location","scale",'maxRetries')) {
			  if (is.na(slot(Ob, sl))) next
			  slname <- paste("$", sl, sep="")
			  if (is.character(slot(Ob, sl))) {
				  cat(sp, slname, ":", omxQuotes(slot(Ob, sl)), '\n')
			  } else {
				  cat(sp, slname, ":", slot(Ob, sl), '\n')
			  }
		  }
		  invisible(Ob)
	  })

#----------------------------------------------------

setClass(Class = "MxComputeTryCatch",
	 contains = "BaseCompute",
	 representation = representation(
	     plan = "MxCompute"))

setMethod("assignId", signature("MxComputeTryCatch"),
	function(.Object, id, defaultFreeSet) {
		.Object <- callNextMethod()
		defaultFreeSet <- .Object@freeSet
		id <- .Object@id
		for (sl in c('plan')) {
			slot(.Object, sl) <- assignId(slot(.Object, sl), id, defaultFreeSet)
			id <- slot(.Object, sl)@id + 1L
		}
		.Object@id <- id
		.Object
	})

setMethod("getFreeVarGroup", signature("MxComputeTryCatch"),
	function(.Object) {
		result <- callNextMethod()
		for (step in c(.Object@plan)) {
			got <- getFreeVarGroup(step)
			if (length(got)) result <- append(result, got)
		}
		result
	})

setMethod("qualifyNames", signature("MxComputeTryCatch"),
	function(.Object, modelname, namespace) {
		.Object <- callNextMethod()
		for (sl in c('plan')) {
			slot(.Object, sl) <- qualifyNames(slot(.Object, sl), modelname, namespace)
		}
		.Object
	})

setMethod("convertForBackend", signature("MxComputeTryCatch"),
	function(.Object, flatModel, model) {
		name <- .Object@name
		for (sl in c('plan')) {
			slot(.Object, sl) <- convertForBackend(slot(.Object, sl), flatModel, model)
		}
		.Object
	})

setMethod("updateFromBackend", signature("MxComputeTryCatch"),
	function(.Object, computes) {
		.Object <- callNextMethod()
		for (sl in c('plan')) {
			slot(.Object, sl) <- updateFromBackend(slot(.Object, sl), computes)
		}
		.Object
	})

setMethod("initialize", "MxComputeTryCatch",
	  function(.Object, freeSet, plan) {
		  .Object@name <- 'compute'
		  .Object@.persist <- TRUE
		  .Object@freeSet <- freeSet
		  .Object@plan <- plan
		  .Object
	  })

##' Execute a sub-compute plan, catching errors
##'
##' \lifecycle{experimental}
##' Any error will be recorded in a subsequent checkpoint. After
##' execution, the context will be reset to continue computation as if
##' no errors has occurred.
##'
##' @param plan compute plan to optimize the model
##' @param ...  Not used.  Forces remaining arguments to be specified by name.
##' @param freeSet names of matrices containing free variables
##' @seealso
##' \link{mxComputeCheckpoint}
##' @aliases
##' MxComputeTryCatch-class
mxComputeTryCatch <- function(plan, ..., freeSet=NA_character_)
{
  prohibitDotdotdot(list(...))
	new("MxComputeTryCatch", freeSet, plan)
}

setMethod("displayCompute", signature(Ob="MxComputeTryCatch", indent="integer"),
	  function(Ob, indent) {
		  callNextMethod()
		  sp <- paste(rep('  ', indent), collapse="")
		  cat(sp, "$plan :", '\n')
		  displayCompute(Ob@plan, indent+1L)
		  invisible(Ob)
	  })

#----------------------------------------------------

setClass(Class = "MxComputeConfidenceInterval",
	 contains = "BaseCompute",
	 representation = representation(
	     plan = "MxCompute",
	   fitfunction = "MxCharOrNumber",
	     constraintType = "character",
	     verbose = "integer"))

setMethod("assignId", signature("MxComputeConfidenceInterval"),
	function(.Object, id, defaultFreeSet) {
		.Object <- callNextMethod()
		defaultFreeSet <- .Object@freeSet
		id <- .Object@id
		for (sl in c('plan')) {
			slot(.Object, sl) <- assignId(slot(.Object, sl), id, defaultFreeSet)
			id <- slot(.Object, sl)@id + 1L
		}
		.Object@id <- id
		.Object
	})

setMethod("getFreeVarGroup", signature("MxComputeConfidenceInterval"),
	function(.Object) {
		result <- callNextMethod()
		for (step in c(.Object@plan)) {
			got <- getFreeVarGroup(step)
			if (length(got)) result <- append(result, got)
		}
		result
	})

setMethod("qualifyNames", signature("MxComputeConfidenceInterval"),
	function(.Object, modelname, namespace) {
		.Object <- callNextMethod()
		for (sl in c('plan')) {
			slot(.Object, sl) <- qualifyNames(slot(.Object, sl), modelname, namespace)
		}
		for (sl in c('fitfunction')) {
			slot(.Object, sl) <- imxConvertIdentifier(slot(.Object, sl), modelname, namespace)
		}
		.Object
	})

setMethod("convertForBackend", signature("MxComputeConfidenceInterval"),
	function(.Object, flatModel, model) {
		name <- .Object@name
		for (sl in c('plan')) {
			slot(.Object, sl) <- convertForBackend(slot(.Object, sl), flatModel, model)
		}
		if (is.character(.Object@fitfunction)) {
			.Object@fitfunction <- imxLocateIndex(flatModel, .Object@fitfunction, .Object)
		}
		.Object
	})

setMethod("initialize", "MxComputeConfidenceInterval",
	  function(.Object, freeSet, plan, verbose, fitfunction, constraintType) {
		  .Object@name <- 'compute'
		  .Object@.persist <- TRUE
		  .Object@freeSet <- freeSet
		  .Object@plan <- plan
		  .Object@verbose <- verbose
		  .Object@fitfunction <- fitfunction
		  .Object@constraintType <- constraintType
		  .Object
	  })

##' Find likelihood-based confidence intervals
##'
##' There are various equivalent ways to pose the optimization
##' problems required to estimate confidence intervals. Most accurate
##' solutions are achieved when the problem is posed using non-linear
##' constraints. However, the available optimizers (CSOLNP, SLSQP, and NPSOL) often have difficulty with non-linear
##' constraints.
##'
##' @param plan compute plan to optimize the model
##' @param ...  Not used.  Forces remaining arguments to be specified by name.
##' @param freeSet names of matrices containing free variables
##' @template args-verbose
##' @param engine \lifecycle{deprecated}
##' @param fitfunction the name of the deviance function
##' @param tolerance \lifecycle{deprecated}
##' @param constraintType one of c('ineq', 'none')
##' @references
##' Neale, M. C. & Miller M. B. (1997). The use of likelihood based
##' confidence intervals in genetic models.  \emph{Behavior Genetics,
##' 27}(2), 113-120.
##'
##' Pek, J. & Wu, H. (2015). Profile likelihood-based confidence intervals and regions for structural equation models.
##' \emph{Psychometrika, 80}(4), 1123-1145.
##'
##' Wu, H. & Neale, M. C. (2012). Adjusted confidence intervals for a
##' bounded parameter. \emph{Behavior genetics, 42}(6), 886-898.
##' @aliases
##' MxComputeConfidenceInterval-class

mxComputeConfidenceInterval <- function(plan, ..., freeSet=NA_character_, verbose=0L,
					engine=NULL, fitfunction='fitfunction',
					tolerance=NA_real_, constraintType='none') {

  prohibitDotdotdot(list(...))
	verbose <- as.integer(verbose)
	new("MxComputeConfidenceInterval", freeSet, plan, verbose, fitfunction, constraintType)
}

setMethod("updateFromBackend", signature("MxComputeConfidenceInterval"),
	function(.Object, computes) {
		.Object <- callNextMethod()
		for (sl in c('plan')) {
			slot(.Object, sl) <- updateFromBackend(slot(.Object, sl), computes)
		}
		.Object
	})

setMethod("displayCompute", signature(Ob="MxComputeConfidenceInterval", indent="integer"),
	  function(Ob, indent) {
		  callNextMethod()
		  sp <- paste(rep('  ', indent), collapse="")
		  cat(sp, "$plan :", '\n')
		  displayCompute(Ob@plan, indent+1L)
		  for (sl in c("verbose")) {
			  if (is.na(slot(Ob, sl))) next
			  slname <- paste("$", sl, sep="")
			  if (is.character(slot(Ob, sl))) {
				  cat(sp, slname, ":", omxQuotes(slot(Ob, sl)), '\n')
			  } else {
				  cat(sp, slname, ":", slot(Ob, sl), '\n')
			  }
		  }
		  invisible(Ob)
	  })

#----------------------------------------------------

setClass(Class = "MxComputePenaltySearch",
	 contains = "BaseCompute",
	 representation = representation(
	     plan = "MxCompute",
       fitfunction = "MxCharOrNumber",
	     verbose = "integer",
       approach = "character",
       ebicGamma = "numeric"))

setMethod("assignId", signature("MxComputePenaltySearch"),
	function(.Object, id, defaultFreeSet) {
		.Object <- callNextMethod()
		defaultFreeSet <- .Object@freeSet
		id <- .Object@id
		for (sl in c('plan')) {
			slot(.Object, sl) <- assignId(slot(.Object, sl), id, defaultFreeSet)
			id <- slot(.Object, sl)@id + 1L
		}
		.Object@id <- id
		.Object
	})

setMethod("getFreeVarGroup", signature("MxComputePenaltySearch"),
	function(.Object) {
		result <- callNextMethod()
		for (step in c(.Object@plan)) {
			got <- getFreeVarGroup(step)
			if (length(got)) result <- append(result, got)
		}
		result
	})

setMethod("qualifyNames", signature("MxComputePenaltySearch"),
	function(.Object, modelname, namespace) {
		.Object <- callNextMethod()
		for (sl in c('plan')) {
			slot(.Object, sl) <- qualifyNames(slot(.Object, sl), modelname, namespace)
		}
		for (sl in c('fitfunction')) {
			slot(.Object, sl) <- imxConvertIdentifier(slot(.Object, sl), modelname, namespace)
		}
		.Object
	})

setMethod("convertForBackend", signature("MxComputePenaltySearch"),
	function(.Object, flatModel, model) {
		name <- .Object@name
		for (sl in c('plan')) {
			slot(.Object, sl) <- convertForBackend(slot(.Object, sl), flatModel, model)
		}
		if (is.character(.Object@fitfunction)) {
			.Object@fitfunction <- imxLocateIndex(flatModel, .Object@fitfunction, .Object)
		}
		.Object
	})

setMethod("initialize", "MxComputePenaltySearch",
	  function(.Object, freeSet, plan, verbose, fitfunction, approach, ebicGamma) {
		  .Object@name <- 'compute'
		  .Object@.persist <- TRUE
		  .Object@freeSet <- freeSet
		  .Object@plan <- plan
		  .Object@verbose <- verbose
		  .Object@fitfunction <- fitfunction
      .Object@approach <- approach
      .Object@ebicGamma <- ebicGamma
		  .Object
	  })

##' Regularize parameter estimates
##'
##' Add a penalty to push some subset of the parameter estimates toward zero.
##'
##' @param plan compute plan to optimize the model
##' @param ...  Not used.  Forces remaining arguments to be specified by name.
##' @param freeSet names of matrices containing free variables
##' @template args-verbose
##' @param fitfunction the name of the deviance function
##' @param approach what fit function to use to compare regularized models? Currently only EBIC is available
##' @param ebicGamma what Gamma value to use for EBIC? Must be between 0 and 1
##' @references
##' Jacobucci, R., Grimm, K. J., & McArdle, J. J. (2016).
##' Regularized structural equation modeling.
##' <i>Structural equation modeling: a multidisciplinary journal, 23</i>(4), 555-566.
##' @aliases
##' MxComputePenaltySearch-class
mxComputePenaltySearch <- function(plan, ..., freeSet=NA_character_, verbose=0L,
                                fitfunction='fitfunction',
                                approach='EBIC', ebicGamma = 0.5) {
  prohibitDotdotdot(list(...))
	verbose <- as.integer(verbose)
	approach <- match.arg(approach)
  if (ebicGamma < 0 || ebicGamma > 1) stop("ebicGamma must be between 0 and 1")
	new("MxComputePenaltySearch", freeSet, plan, verbose, fitfunction, approach,
      ebicGamma)
}

setMethod("updateFromBackend", signature("MxComputePenaltySearch"),
	function(.Object, computes) {
		.Object <- callNextMethod()
		for (sl in c('plan')) {
			slot(.Object, sl) <- updateFromBackend(slot(.Object, sl), computes)
		}
		.Object
	})

setMethod("displayCompute", signature(Ob="MxComputePenaltySearch", indent="integer"),
	  function(Ob, indent) {
		  callNextMethod()
		  sp <- paste(rep('  ', indent), collapse="")
		  cat(sp, "$plan :", '\n')
		  displayCompute(Ob@plan, indent+1L)
		  for (sl in c("verbose")) {
			  if (is.na(slot(Ob, sl))) next
			  slname <- paste("$", sl, sep="")
			  if (is.character(slot(Ob, sl))) {
				  cat(sp, slname, ":", omxQuotes(slot(Ob, sl)), '\n')
			  } else {
				  cat(sp, slname, ":", slot(Ob, sl), '\n')
			  }
		  }
		  invisible(Ob)
	  })

#----------------------------------------------------

setClass(Class = "MxComputeNewtonRaphson",
	 contains = "BaseCompute",
	 representation = representation(
	   fitfunction = "MxCharOrNumber",
	   maxIter = "integer",
	   tolerance = "numeric",
	   verbose = "integer"))

setMethod("qualifyNames", signature("MxComputeNewtonRaphson"),
	function(.Object, modelname, namespace) {
		.Object <- callNextMethod()
		for (sl in c('fitfunction')) {
			slot(.Object, sl) <- imxConvertIdentifier(slot(.Object, sl), modelname, namespace)
		}
		.Object
	})

setMethod("convertForBackend", signature("MxComputeNewtonRaphson"),
	function(.Object, flatModel, model) {
		name <- .Object@name
		if (is.character(.Object@fitfunction)) {
			.Object@fitfunction <- imxLocateIndex(flatModel, .Object@fitfunction, .Object)
		}
		.Object
	})

setMethod("initialize", "MxComputeNewtonRaphson",
	  function(.Object, freeSet, fit, maxIter, tolerance, verbose) {
		  .Object@name <- 'compute'
		  .Object@.persist <- TRUE
		  .Object@freeSet <- freeSet
		  .Object@fitfunction <- fit
		  .Object@maxIter <- maxIter
		  .Object@tolerance <- tolerance
		  .Object@verbose <- verbose
		  .Object
	  })

##' Optimize parameters using the Newton-Raphson algorithm
##'
##' This optimizer requires analytic 1st and 2nd derivatives of the
##' fit function. Box constraints are supported. Parameters can
##' approach box constraints but will not leave the feasible region
##' (even by some small epsilon>0). Non-finite fit values are
##' interpreted as soft feasibility constraints. That is, when a
##' non-finite fit is encountered, line search is continued after the
##' step size is multiplied by 10\%. Comprehensive diagnostics are
##' available by increasing the verbose level.
##'
##' @param freeSet names of matrices containing free variables
##' @param ...  Not used.  Forces remaining arguments to be specified by name.
##' @param fitfunction name of the fitfunction (defaults to 'fitfunction')
##' @param maxIter maximum number of iterations
##' @param tolerance optimization is considered converged when the maximum relative change in fit is less than tolerance
##' @template args-verbose
##' @aliases
##' MxComputeNewtonRaphson-class
##' @references
##' Luenberger, D. G. & Ye, Y. (2008). \emph{Linear and nonlinear programming.} Springer.

mxComputeNewtonRaphson <- function(freeSet=NA_character_, ..., fitfunction='fitfunction', maxIter = 100L,
				   tolerance=1e-12, verbose=0L)
{
  prohibitDotdotdot(list(...))
	verbose <- as.integer(verbose)
	maxIter <- as.integer(maxIter)
	new("MxComputeNewtonRaphson", freeSet, fitfunction, maxIter, tolerance, verbose)
}

setMethod("displayCompute", signature(Ob="MxComputeNewtonRaphson", indent="integer"),
	  function(Ob, indent) {
		  callNextMethod();
		  sp <- paste(rep('  ', indent), collapse="")
		  for (sl in c("fitfunction", "maxIter", "tolerance", "verbose")) {
			  slname <- paste("$", sl, sep="")
			  if (is.character(slot(Ob, sl))) {
				  cat(sp, slname, ":", omxQuotes(slot(Ob, sl)), '\n')
			  } else {
				  cat(sp, slname, ":", slot(Ob, sl), '\n')
			  }
		  }
		  invisible(Ob)
	  })

#----------------------------------------------------

setClass(Class = "MxComputeSimAnnealing",
	 contains = "BaseCompute",
	 representation = representation(
	   fitfunction = "MxCharOrNumber",
	   verbose = "integer",
	   plan = "MxCompute",
	   method = "character",
	   control = "list",
	   defaultGradientStepSize = "numeric",
	   defaultFunctionPrecision = "numeric"))

setMethod("assignId", signature("MxComputeSimAnnealing"),
	function(.Object, id, defaultFreeSet) {
		.Object <- callNextMethod()
		defaultFreeSet <- .Object@freeSet
		id <- .Object@id
		for (sl in c('plan')) {
			slot(.Object, sl) <- assignId(slot(.Object, sl), id, defaultFreeSet)
			id <- slot(.Object, sl)@id + 1L
		}
		.Object@id <- id
		.Object
	})

setMethod("getFreeVarGroup", signature("MxComputeSimAnnealing"),
	function(.Object) {
		result <- callNextMethod()
		for (step in c(.Object@plan)) {
			got <- getFreeVarGroup(step)
			if (length(got)) result <- append(result, got)
		}
		result
	})

setMethod("qualifyNames", signature("MxComputeSimAnnealing"),
	function(.Object, modelname, namespace) {
		.Object <- callNextMethod()
		for (sl in c('fitfunction')) {
			slot(.Object, sl) <- imxConvertIdentifier(slot(.Object, sl), modelname, namespace)
		}
		for (sl in c('plan')) {
			slot(.Object, sl) <- qualifyNames(slot(.Object, sl), modelname, namespace)
		}
		.Object
	})

setMethod("convertForBackend", signature("MxComputeSimAnnealing"),
	function(.Object, flatModel, model) {
		name <- .Object@name
		if (is.character(.Object@fitfunction)) {
			.Object@fitfunction <- imxLocateIndex(flatModel, .Object@fitfunction, .Object)
		}
		for (sl in c('plan')) {
			slot(.Object, sl) <- convertForBackend(slot(.Object, sl), flatModel, model)
		}
		.Object
	})

setMethod("updateFromBackend", signature("MxComputeSimAnnealing"),
	function(.Object, computes) {
		.Object <- callNextMethod()
		for (sl in c('plan')) {
			slot(.Object, sl) <- updateFromBackend(slot(.Object, sl), computes)
		}
		.Object
	})

setMethod("initialize", "MxComputeSimAnnealing",
	function(.Object, freeSet, fit, verbose, plan, method, control,
		 defaultGradientStepSize, defaultFunctionPrecision) {
		  .Object@name <- 'compute'
		  .Object@.persist <- TRUE
		  .Object@freeSet <- freeSet
		  .Object@fitfunction <- fit
		  .Object@verbose <- verbose
		  .Object@plan <- plan
		  .Object@method <- method
		  .Object@control <- control
		  .Object@defaultGradientStepSize <- defaultGradientStepSize
		  .Object@defaultFunctionPrecision <- defaultFunctionPrecision
		  .Object
	  })

mxComputeSimAnnealing <- function(freeSet=NA_character_, ..., fitfunction='fitfunction',
			   plan=mxComputeOnce('fitfunction','fit'),
			   verbose=0L, method=c("tsallis1996", "ingber2012"),
			   control=list(),
			   defaultGradientStepSize=imxAutoOptionValue("Gradient step size"),
			   defaultFunctionPrecision=imxAutoOptionValue("Function precision"))
{
  prohibitDotdotdot(list(...))
	method <- match.arg(method)
	verbose <- as.integer(verbose)
	new("MxComputeSimAnnealing", freeSet, fitfunction, verbose, plan, method, control,
		defaultGradientStepSize, defaultFunctionPrecision)
}

#----------------------------------------------------

setClass(Class = "ComputeSteps",
	 contains = "BaseCompute",
	 representation = representation(
	   steps = "list"))

setMethod("getFreeVarGroup", signature("ComputeSteps"),
	function(.Object) {
		result <- callNextMethod()
		for (step in .Object@steps) {
			got <- getFreeVarGroup(step)
			if (length(got)) result <- append(result, got)
		}
		result
	})

setMethod("assignId", signature("ComputeSteps"),
	function(.Object, id, defaultFreeSet) {
		if (length(.Object@freeSet) == 1 && is.na(.Object@freeSet)) .Object@freeSet <- defaultFreeSet
		defaultFreeSet <- .Object@freeSet
		steps <- .Object@steps
		if (length(steps)) for (sx in 1:length(steps)) {
			steps[[sx]] <- assignId(steps[[sx]], id, defaultFreeSet)
			id <- steps[[sx]]@id + 1L
		}
		.Object@steps <- steps
		.Object@id <- id
		.Object
	})

setMethod("qualifyNames", signature("ComputeSteps"),
	function(.Object, modelname, namespace) {
		.Object <- callNextMethod()
		.Object@name <- imxIdentifier(modelname, .Object@name)
		.Object@steps <- lapply(.Object@steps, function (c) qualifyNames(c, modelname, namespace))
		.Object
	})

setMethod("convertForBackend", signature("ComputeSteps"),
	function(.Object, flatModel, model) {
		.Object@steps <- lapply(.Object@steps, function (c) convertForBackend(c, flatModel, model))
		.Object
	})

setMethod("updateFromBackend", signature("ComputeSteps"),
	function(.Object, computes) {
		.Object <- callNextMethod()
		.Object@steps <- lapply(.Object@steps, function (c) updateFromBackend(c, computes))
		.Object
	})

#----------------------------------------------------

setClass(Class = "MxComputeIterate",
	 contains = "ComputeSteps",
	 representation = representation(
	   maxIter = "integer",
	   tolerance = "numeric",
	   verbose = "integer",
	   maxDuration = "numeric"))

setMethod("initialize", "MxComputeIterate",
	  function(.Object, steps, maxIter, tolerance, verbose, freeSet, maxDuration) {
		  .Object@name <- 'compute'
		  .Object@.persist <- TRUE
		  .Object@steps <- steps
		  .Object@maxIter <- maxIter
		  .Object@tolerance <- tolerance
		  .Object@verbose <- verbose
		  .Object@freeSet <- freeSet
		  .Object@maxDuration <- maxDuration
		  .Object
	  })

##' Repeatedly invoke a series of compute objects until change is less than tolerance
##'
##' One step (typically the last) must compute the fit or maxAbsChange.
##'
##' @param steps a list of compute objects
##' @param ...  Not used.  Forces remaining arguments to be specified by name.
##' @param maxIter the maximum number of iterations
##' @param tolerance iterates until maximum relative change is less than tolerance
##' @template args-verbose
##' @param freeSet Names of matrices containing free variables.
##' @param maxDuration the maximum amount of time (in seconds) to iterate
##' @aliases
##' MxComputeIterate-class
mxComputeIterate <- function(steps, ..., maxIter=500L, tolerance=1e-9, verbose=0L, freeSet=NA_character_,
			     maxDuration=as.numeric(NA)) {
  prohibitDotdotdot(list(...))
	verbose <- as.integer(verbose)
	maxIter <- as.integer(maxIter)
	new("MxComputeIterate", steps=steps, maxIter=maxIter, tolerance=tolerance,
	    verbose, freeSet, maxDuration)
}

setMethod("displayCompute", signature(Ob="MxComputeIterate", indent="integer"),
	  function(Ob, indent) {
		  callNextMethod();
		  sp <- paste(rep('  ', indent), collapse="")
		  cat(sp, "maxIter :", Ob@maxIter, '\n')
		  cat(sp, "maxDuration :", Ob@maxDuration, '\n')
		  cat(sp, "tolerance :", Ob@tolerance, '\n')
		  cat(sp, "verbose :", Ob@verbose, '\n')
		  for (step in 1:length(Ob@steps)) {
			  cat(sp, "steps[[", step, "]] :", '\n')
			  displayCompute(Ob@steps[[step]], indent+1L)
		  }
		  invisible(Ob)
	  })

#----------------------------------------------------
setClass(Class = "MxComputeLoop",
	 contains = "ComputeSteps",
	 representation = representation(
		 indices = "integer",
	   maxIter = "integer",
	   maxDuration = "numeric",
	   verbose="integer",
	   startFrom="integer"))

setMethod("initialize", "MxComputeLoop",
	  function(.Object, steps, indices, maxIter, freeSet, maxDuration, verbose, startFrom) {
		  .Object@name <- 'compute'
		  .Object@.persist <- TRUE
		  .Object@steps <- steps
		  .Object@indices <- indices
		  .Object@maxIter <- maxIter
		  .Object@freeSet <- freeSet
		  .Object@maxDuration <- maxDuration
		  .Object@verbose <- verbose
		  .Object@startFrom <- startFrom
		  .Object
	  })

##' Repeatedly invoke a series of compute objects
##'
##' @param steps a list of compute objects
##' @param ...  Not used.  Forces remaining arguments to be specified
##'         by name.
##' @param i the values to iterate over
##' @param maxIter the maximum number of iterations
##' @param freeSet Names of matrices containing free variables.
##' @param maxDuration the maximum amount of time (in seconds) to
##'         iterate
##' @param startFrom When \code{i=NULL}, permits starting from an index greater than 1.
##' @template args-verbose
##' @description When \code{i} is given then these values are iterated
##'         over instead of the sequence 1 to the number of
##'         iterations.
##' @aliases MxComputeLoop-class mxComputeBenchmark
mxComputeLoop <- function(steps, ..., i=NULL, maxIter=as.integer(NA), freeSet=NA_character_,
			     maxDuration=as.numeric(NA), verbose=0L, startFrom=1L) {
  prohibitDotdotdot(list(...))
	if (length(i) && startFrom != 1L) {
	  warning("Argument startFrom is ignored when i is provided")
	}
	maxIter <- as.integer(maxIter)
	new("MxComputeLoop", steps=steps, indices=as.integer(i), maxIter=maxIter,
	    freeSet, maxDuration, as.integer(verbose), as.integer(startFrom))
}

setMethod("displayCompute", signature(Ob="MxComputeLoop", indent="integer"),
	  function(Ob, indent) {
		  callNextMethod();
		  sp <- paste(rep('  ', indent), collapse="")
		  cat(sp, "indices :", Ob@indices, '\n')
		  cat(sp, "maxIter :", Ob@maxIter, '\n')
		  cat(sp, "maxDuration :", Ob@maxDuration, '\n')
		  for (step in 1:length(Ob@steps)) {
			  cat(sp, "steps[[", step, "]] :", '\n')
			  displayCompute(Ob@steps[[step]], indent+1L)
		  }
		  invisible(Ob)
	  })

mxComputeBenchmark <- mxComputeLoop

#----------------------------------------------------

setClass(Class = "MxComputeEM",
	 contains = "BaseCompute",
	 representation = representation(
	     estep = "MxCompute",
	     mstep = "MxCompute",
	     observedFit = "MxCharOrNumber",
	     maxIter = "integer",
	     tolerance = "numeric",
	     verbose = "integer",
	     accel="character",
	     information="character",
	     infoArgs="list"))

setMethod("assignId", signature("MxComputeEM"),
	function(.Object, id, defaultFreeSet) {
		.Object <- callNextMethod()
		defaultFreeSet <- .Object@freeSet
		id <- .Object@id
		for (sl in c('estep', 'mstep')) {
			slot(.Object, sl) <- assignId(slot(.Object, sl), id, defaultFreeSet)
			id <- slot(.Object, sl)@id + 1L
		}
		.Object@id <- id
		.Object
	})

setMethod("getFreeVarGroup", signature("MxComputeEM"),
	function(.Object) {
		result <- callNextMethod()
		for (step in c(.Object@estep, .Object@mstep)) {
			got <- getFreeVarGroup(step)
			if (length(got)) result <- append(result, got)
		}
		result
	})

setMethod("qualifyNames", signature("MxComputeEM"),
	function(.Object, modelname, namespace) {
		.Object <- callNextMethod()
		for (sl in c('observedFit')) {
			slot(.Object, sl) <- imxConvertIdentifier(slot(.Object, sl), modelname, namespace)
		}
		for (sl in c('estep', 'mstep')) {
			slot(.Object, sl) <- qualifyNames(slot(.Object, sl), modelname, namespace)
		}
		.Object@infoArgs$fitfunction <-
		    imxConvertIdentifier(.Object@infoArgs$fitfunction, modelname, namespace)
		.Object
	})

setMethod("convertForBackend", signature("MxComputeEM"),
	function(.Object, flatModel, model) {
		name <- .Object@name
		if (length(.Object@observedFit) != 1) stop("MxComputeEM requires a single observedFit function")
		if (any(!is.integer(.Object@observedFit))) {
			algNum <- match(.Object@observedFit, append(names(flatModel@algebras),
								    names(flatModel@fitfunctions)))
			if (any(is.na(algNum))) {
				stop(paste("MxComputeEM: observedFit fit function", .Object@observedFit, "not found"))
			}
			.Object@observedFit <- algNum - 1L
		}
		for (sl in c('estep', 'mstep')) {
			slot(.Object, sl) <- convertForBackend(slot(.Object, sl), flatModel, model)
		}
		fit <- match(.Object@infoArgs$fitfunction,
			     append(names(flatModel@algebras), names(flatModel@fitfunctions)))
		if (any(is.na(fit))) {
			stop(paste("ComputeEM: cannot find fitfunction",
				   omxQuotes(.Object@infoArgs$fitfunction[is.na(fit)]), "in infoArgs"))
		}
		.Object@infoArgs$fitfunction <-fit - 1L

		.Object
	})

setMethod("updateFromBackend", signature("MxComputeEM"),
	function(.Object, computes) {
		.Object <- callNextMethod()
		for (sl in c('estep', 'mstep')) {
			slot(.Object, sl) <- updateFromBackend(slot(.Object, sl), computes)
		}
		.Object
	})

setMethod("initialize", "MxComputeEM",
	  function(.Object, estep, mstep, observedFit, maxIter, tolerance,
		   verbose, accel, information, freeSet, infoArgs) {
		  .Object@name <- 'compute'
		  .Object@.persist <- TRUE
		  .Object@estep <- estep
		  .Object@mstep <- mstep
		  .Object@observedFit <- observedFit
		  .Object@maxIter <- maxIter
		  .Object@tolerance <- tolerance
		  .Object@verbose <- verbose
		  .Object@accel <- accel
		  .Object@information <- information
		  .Object@freeSet <- freeSet
		  .Object@infoArgs <- infoArgs
		  .Object
	  })

##' Fit a model using DLR's (1977) Expectation-Maximization (EM) algorithm
##'
##' The EM algorithm constitutes the following steps: Start with an
##' initial parameter vector. Predict the missing data to form a
##' completed data model. Optimize the completed data model to obtain
##' a new parameter vector. Repeat these steps until convergence
##' criteria are met.
##'
##' The arguments to this function have evolved.  The old style
##' \code{mxComputeEM(e,p,mstep=m)} is equivalent to the new style
##' \code{mxComputeEM(estep=mxComputeOnce(e,p), mstep=m)}. This change
##' allows the API to more closely match the literature on the E-M
##' method.  You might use \code{mxAlgebra(..., recompute='onDemand')} to
##' contain the results of the E-step and then cause this algebra to
##' be recomputed using \code{mxComputeOnce}.
##'
##' This compute plan does not work with any and all expectations. It
##' requires a special kind of expectation that can predict its
##' missing data to create a completed data model.
##'
##' The EM algorithm does not produce a parameter covariance matrix
##' for standard errors. The Oakes (1999) direct method and S-EM, an
##' implementation of Meng & Rubin (1991), are included.
##'
##' Ramsay (1975) was recommended in Bock, Gibbons, & Muraki (1988).
##'
##' @param expectation a vector of expectation names \lifecycle{deprecated}
##' @param predict what to predict from the observed data \lifecycle{deprecated}
##' @param mstep a compute plan to optimize the completed data model
##' @param observedFit the name of the observed data fit function (defaults to "fitfunction")
##' @param ...  Not used.  Forces remaining arguments to be specified by name.
##' @param maxIter maximum number of iterations
##' @param tolerance optimization is considered converged when the maximum relative change in fit is less than tolerance
##' @template args-verbose
##' @param freeSet names of matrices containing free variables
##' @param accel name of acceleration method ("varadhan2008" or "ramsay1975")
##' @param information name of information matrix approximation method
##' @param infoArgs arguments to control the information matrix method
##' @param estep a compute plan to perform the expectation step
##' @seealso
##' \link[=mxAlgebra]{MxAlgebra}, \link{mxComputeOnce}
##' @aliases
##' MxComputeEM-class
##' @references
##'
##' Bock, R. D., Gibbons, R., & Muraki, E. (1988). Full-information
##' item factor analysis. \emph{Applied Psychological Measurement,
##' 6}(4), 431-444.
##'
##' Dempster, A. P., Laird, N. M., & Rubin, D. B. (1977). Maximum likelihood from
##' incomplete data via the EM algorithm. \emph{Journal of the Royal Statistical Society.
##' Series B (Methodological)}, 1-38.
##'
##' Meng, X.-L. & Rubin, D. B. (1991). Using EM to obtain asymptotic variance-covariance
##' matrices: The SEM algorithm. \emph{Journal of the American Statistical Association,
##' 86} (416), 899-909.
##'
##' Oakes, D. (1999). Direct calculation of the information matrix via
##' the EM algorithm.  \emph{Journal of the Royal Statistical Society:
##' Series B (Statistical Methodology), 61}(2), 479-482.
##'
##' Ramsay, J. O. (1975). Solving implicit equations in psychometric data analysis.
##' \emph{Psychometrika, 40} (3), 337-360.
##'
##' Varadhan, R. & Roland, C. (2008). Simple and globally convergent
##' methods for accelerating the convergence of any EM
##' algorithm. \emph{Scandinavian Journal of Statistics, 35}, 335-353.
##' @examples
##' library(OpenMx)
##' set.seed(190127)
##'
##' N <- 200
##' x <- matrix(c(rnorm(N/2,0,1),
##'               rnorm(N/2,3,1)),ncol=1,dimnames=list(NULL,"x"))
##' data4mx <- mxData(observed=x,type="raw")
##'
##' class1 <- mxModel("Class1",
##' 	mxMatrix(type="Full",nrow=1,ncol=1,free=TRUE,values=0,name="Mu"),
##' 	mxMatrix(type="Full",nrow=1,ncol=1,free=TRUE,values=4,name="Sigma"),
##' 	mxExpectationNormal(covariance="Sigma",means="Mu",dimnames="x"),
##' 	mxFitFunctionML(vector=TRUE))
##'
##' class2 <- mxRename(class1, "Class2")
##'
##' mm <- mxModel(
##' 	"Mixture", data4mx, class1, class2,
##' 	mxAlgebra((1-Posteriors) * Class1.fitfunction, name="PL1"),
##' 	mxAlgebra(Posteriors * Class2.fitfunction, name="PL2"),
##' 	mxAlgebra(PL1 + PL2, name="PL"),
##' 	mxAlgebra(PL2 / PL,  recompute='onDemand',
##' 	          initial=matrix(runif(N,.4,.6), nrow=N, ncol = 1), name="Posteriors"),
##' 	mxAlgebra(-2*sum(log(PL)), name="FF"),
##' 	mxFitFunctionAlgebra(algebra="FF"),
##' 	mxComputeEM(
##' 	  estep=mxComputeOnce("Mixture.Posteriors"),
##' 	  mstep=mxComputeGradientDescent(fitfunction="Mixture.fitfunction")))
##'
##' mm <- mxOption(mm, "Max minutes", 1/20)  # remove this line to find optimum
##' mmfit <- mxRun(mm)
##' summary(mmfit)
mxComputeEM <- function(expectation=NULL, predict=NA_character_, mstep, observedFit="fitfunction", ...,
			maxIter=500L, tolerance=1e-9, verbose=0L, freeSet=NA_character_,
			accel="varadhan2008", information=NA_character_, infoArgs=list(), estep=NULL) {
  prohibitDotdotdot(list(...))
	verbose <- as.integer(verbose)
	maxIter <- as.integer(maxIter)
	accel <- as.character(accel)

	# backward compatibility for original API
	if (length(expectation)) {
		if (length(estep)) stop("You cannot provide both 'expectation' and 'estep' arguments")
		estep <- mxComputeOnce(expectation, predict)
		mstep <- mxComputeSequence(list(mstep, mxComputeOnce(expectation)))
	}

	new("MxComputeEM", estep, mstep, observedFit, maxIter=maxIter,
	    tolerance=tolerance, verbose, accel, information, freeSet, infoArgs)
}

setMethod("displayCompute", signature(Ob="MxComputeEM", indent="integer"),
	  function(Ob, indent) {
		  callNextMethod();
		  sp <- paste(rep('  ', indent), collapse="")
		  cat(sp, "$estep :", '\n')
		  displayCompute(Ob@estep, indent+1L)
		  cat(sp, "$mstep :", '\n')
		  displayCompute(Ob@mstep, indent+1L)
		  for (sl in c("observedFit", "maxIter", "tolerance", "verbose", "accel")) {
			  slname <- paste("$", sl, sep="")
			  if (is.character(slot(Ob, sl))) {
				  cat(sp, slname, ":", omxQuotes(slot(Ob, sl)), '\n')
			  } else {
				  cat(sp, slname, ":", slot(Ob, sl), '\n')
			  }
		  }
		  if (!is.na(Ob@information)) {
			  cat(sp, "$information :", Ob@information, '\n')
		  }
		  invisible(Ob)
	  })

#----------------------------------------------------
#Mike Hunter's "better match.arg" function
match.barg <- function (arg, choices, several.ok = FALSE)
{
	if (missing(choices)) {
		formal.args <- formals(sys.function(sys.parent()))
		choices <- eval(formal.args[[deparse(substitute(arg))]])
	}
	if (is.null(arg))
		return(choices[1L])
	else if (!is.character(arg))
		stop("'arg' must be NULL or a character vector")
	if (!several.ok) {
		if (identical(arg, choices))
			return(arg[1L])
		if (length(arg) > 1L)
			stop("'arg' must be of length 1")
	}
	else if (length(arg) == 0L)
		stop("'arg' must be of length >= 1")
	i <- pmatch(arg, choices, nomatch = 0L, duplicates.ok = TRUE)
	if (all(i == 0L))
		stop(gettextf("%s should be one of %s", omxQuotes(arg), omxQuotes(choices)
		), domain = NA)
	i <- i[i > 0L]
	if (!several.ok && length(i) > 1)
		stop("there is more than one match in 'match.arg'")
	choices[i]
}

mxComputeNelderMead <- function(
	freeSet=NA_character_, fitfunction="fitfunction", verbose=0L,
	nudgeZeroStarts=mxOption(NULL,"Nudge zero starts"),
	maxIter=NULL,	...,
	alpha=1, betao=0.5, betai=0.5, gamma=2, sigma=0.5, bignum=1e35,
	iniSimplexType=c("regular","right","smartRight","random"),
	iniSimplexEdge=1, iniSimplexMat=NULL, greedyMinimize=FALSE,
	altContraction=FALSE, degenLimit=0, stagnCtrl=c(-1L,-1L),
	validationRestart=TRUE,
	xTolProx=1e-8, fTolProx=1e-8,
	doPseudoHessian=TRUE,
	ineqConstraintMthd=c("soft","eqMthd"),
	eqConstraintMthd=c("GDsearch","soft","backtrack","l1p"),
	backtrackCtrl=c(0.5,5),
	centerIniSimplex=FALSE
	){
  prohibitDotdotdot(list(...))
	verbose <- as.integer(verbose[1])
	maxIter <- as.integer(maxIter[1])
	if(is.character(nudgeZeroStarts[1])){
		if(substr(nudgeZeroStarts[1],1,1) %in% c("Y","y")){nudgeZeroStarts <- TRUE}
		else if(substr(nudgeZeroStarts[1],1,1) %in% c("N","n")){nudgeZeroStarts <- FALSE}
		else{stop("unrecognized character string provided as argument 'nudgeZeroStarts'")}
	}
	alpha <- as.numeric(alpha[1])
	if(alpha<=0){stop("reflection coefficient 'alpha' must be positive")}
	betao <- as.numeric(betao[1])
	betai <- as.numeric(betai[1])
	if(any(betao<=0, betao>=1, betai<=0, betai>=1)){
		stop("contraction coefficients 'betao' and 'betai' must both be within unit interval (0,1)")
	}
	gamma <- as.numeric(gamma[1])
	#Allow user to provide non-positive gamma to "turn off" expanstion transformations:
	if(gamma>0 && gamma<=alpha){
		stop("if positive, expansion coefficient 'gamma' must be greater than reflection coefficient 'alpha'")
	}
	sigma <- as.numeric(sigma[1])
	#Allow user to provide non-positive sigma to "turn off" shrinks:
	if(sigma>=1){stop("shrink coefficient 'sigma' must be less than 1.0")}
	bignum <- as.numeric(bignum[1])
	iniSimplexType <- as.character(match.barg(iniSimplexType,c("regular","right","smartRight","random")))
	iniSimplexEdge <- as.numeric(iniSimplexEdge[1])
	if(length(iniSimplexMat)){
		iniSimplexMat <- as.matrix(iniSimplexMat)
		iniSimplexColnames <- colnames(iniSimplexMat)
	}
	else{
		iniSimplexColnames <- NULL
		iniSimplexMat <- NULL
	}
	greedyMinimize <- as.logical(greedyMinimize[1])
	altContraction <- as.logical(altContraction[1])
	degenLimit <- as.numeric(degenLimit[1])
	if(degenLimit<0 || degenLimit>pi){
		stop("'degenLimit' must be within interval [0,pi]")
	}
	if(length(stagnCtrl)<2){stop("'stagnCtrl' must be an integer vector of length 2")}
	stagnCtrl <- as.integer(stagnCtrl[1:2])
	validationRestart <- as.logical(validationRestart[1])
	xTolProx <- as.numeric(xTolProx[1])
	fTolProx <- as.numeric(fTolProx[1])
	doPseudoHessian <- as.logical(doPseudoHessian[1])
	ineqConstraintMthd <- as.character(match.barg(ineqConstraintMthd,c("soft","eqMthd")))
	eqConstraintMthd <- as.character(match.barg(eqConstraintMthd,c("GDsearch","soft","backtrack","l1p")))
	if(length(backtrackCtrl)<2){stop("'backtrackCtrl' must be a numeric vector of length 2")}
	backtrackCtrl1 <- as.numeric(backtrackCtrl[1])
	backtrackCtrl2 <- as.integer(backtrackCtrl[2])
	centerIniSimplex <- as.logical(centerIniSimplex[1])
	return(new("MxComputeNelderMead", freeSet, fitfunction, verbose, nudgeZeroStarts, maxIter, alpha,
						 betao, betai, gamma, sigma, bignum, iniSimplexType, iniSimplexEdge, iniSimplexMat,
						 iniSimplexColnames, validationRestart,
						 greedyMinimize, altContraction, degenLimit, stagnCtrl, xTolProx, fTolProx,
						 ineqConstraintMthd, eqConstraintMthd, backtrackCtrl1, backtrackCtrl2, doPseudoHessian,
						 centerIniSimplex))
}

setClass(
	Class="MxComputeNelderMead",
	contains="BaseCompute",
	representation=representation(
		fitfunction="MxCharOrNumber",
		verbose="integer",
		nudgeZeroStarts="MxCharOrLogical",
		maxIter="integer",
		defaultMaxIter="logical",
		alpha="numeric",
		betao="numeric",
		betai="numeric",
		gamma="numeric",
		sigma="numeric",
		bignum="numeric",
		iniSimplexType="character",
		iniSimplexEdge="numeric",
		iniSimplexMat="MxOptionalMatrix",
		.iniSimplexColnames="MxOptionalChar",
		greedyMinimize="logical",
		altContraction="logical",
		degenLimit="numeric",
		stagnCtrl="integer",
		validationRestart="logical",
		xTolProx="numeric",
		fTolProx="numeric",
		doPseudoHessian="logical",
		ineqConstraintMthd="character",
		eqConstraintMthd="character",
		backtrackCtrl1="numeric",
		backtrackCtrl2="integer",
		centerIniSimplex="logical"))

#TODO: a user or developer might someday want to directly use this low-level 'initialize' method instead of the high-level constructor function,
#so typecasting should also occur here:
setMethod(
	"initialize", "MxComputeNelderMead",
	function(.Object, freeSet, fitfunction, verbose, nudgeZeroStarts, maxIter, alpha,
					 betao, betai, gamma, sigma, bignum, iniSimplexType, iniSimplexEdge, iniSimplexMat,
					 iniSimplexColnames, validationRestart,
					 greedyMinimize, altContraction, degenLimit, stagnCtrl, xTolProx, fTolProx,
					 ineqConstraintMthd, eqConstraintMthd, backtrackCtrl1, backtrackCtrl2, doPseudoHessian,
					 centerIniSimplex){
		.Object@name <- 'compute'
		.Object@.persist <- TRUE
		.Object@freeSet <- freeSet
		.Object@fitfunction <- fitfunction
		.Object@verbose <- verbose
		.Object@nudgeZeroStarts <- nudgeZeroStarts
		.Object@defaultMaxIter <- ifelse(length(maxIter),FALSE,TRUE)
		.Object@maxIter <- as.integer(maxIter)

		.Object@alpha <- alpha
		.Object@betao <- betao
		.Object@betai <- betai
		.Object@gamma <- gamma
		.Object@sigma <- sigma
		.Object@bignum <- bignum
		.Object@iniSimplexType <- iniSimplexType
		.Object@iniSimplexEdge <- iniSimplexEdge
		if(!length(iniSimplexMat)){.Object@iniSimplexMat <- NULL}
		else{.Object@iniSimplexMat <- iniSimplexMat}
		.Object@.iniSimplexColnames <- iniSimplexColnames
		.Object@greedyMinimize <- greedyMinimize
		.Object@altContraction <- altContraction
		.Object@degenLimit <- degenLimit
		.Object@stagnCtrl <- stagnCtrl
		.Object@validationRestart <- validationRestart
		.Object@xTolProx <- xTolProx
		.Object@fTolProx <- fTolProx
		.Object@doPseudoHessian <- doPseudoHessian
		.Object@ineqConstraintMthd <- ineqConstraintMthd
		.Object@eqConstraintMthd <- eqConstraintMthd
		.Object@backtrackCtrl1 <- backtrackCtrl1
		.Object@backtrackCtrl2 <- backtrackCtrl2
		.Object@centerIniSimplex <- centerIniSimplex
		.Object
	})

setMethod("qualifyNames", signature("MxComputeNelderMead"),
					function(.Object, modelname, namespace) {
						.Object <- callNextMethod()
						for (sl in c('fitfunction')) {
							slot(.Object, sl) <- imxConvertIdentifier(slot(.Object, sl), modelname, namespace)
						}
						if(length(.Object@iniSimplexMat)){.Object@.iniSimplexColnames <- colnames(.Object@iniSimplexMat)}
						.Object
					})

setMethod("convertForBackend", signature("MxComputeNelderMead"),
					function(.Object, flatModel, model) {
						name <- .Object@name
						if (is.character(.Object@fitfunction)) {
							.Object@fitfunction <- imxLocateIndex(flatModel, .Object@fitfunction, .Object)
						}
						.Object
					})

#----------------------------------------------------

setClass(Class = "MxComputeNumericDeriv",
	 contains = "BaseCompute",
	 representation = representation(
	   fitfunction = "MxCharOrNumber",
	     parallel = "logical",
	     stepSize = "numeric",
	     iterations = "integer",
	     verbose="integer",
	     knownHessian="MxOptionalMatrix",
     checkGradient="logical",
     hessian="logical",
     analytic="logical"))

setMethod("qualifyNames", signature("MxComputeNumericDeriv"),
	function(.Object, modelname, namespace) {
		.Object <- callNextMethod();
		.Object@fitfunction <- imxConvertIdentifier(.Object@fitfunction, modelname, namespace)
		.Object
	})

setMethod("convertForBackend", signature("MxComputeNumericDeriv"),
	function(.Object, flatModel, model) {
		name <- .Object@name
		if (is.character(.Object@fitfunction)) {
			.Object@fitfunction <- imxLocateIndex(flatModel, .Object@fitfunction, .Object)
		}
		.Object
	})

setMethod("initialize", "MxComputeNumericDeriv",
	  function(.Object, freeSet, fit, parallel, stepSize, iterations, verbose, knownHessian,
		   checkGradient, hessian, analytic) {
		  .Object@name <- 'compute'
		  .Object@.persist <- TRUE
		  .Object@freeSet <- freeSet
		  .Object@fitfunction <- fit
		  .Object@parallel <- parallel
		  .Object@stepSize <- stepSize
		  .Object@iterations <- iterations
		  .Object@verbose <- verbose
		  .Object@knownHessian <- knownHessian
		  .Object@checkGradient <- checkGradient
		  .Object@hessian <- hessian
		  .Object@analytic <- analytic
		  .Object
	  })

adjustDefaultNumericDeriv <- function(m, iterations, stepSize) {
	for (nd in 1:length(m@compute@steps)) {
		if (is(m@compute@steps[[nd]], "MxComputeNumericDeriv")) {
			m@compute@steps[[nd]]$iterations <- iterations
			m@compute@steps[[nd]]$stepSize <- stepSize
			break
		}
	}
	m
}

##' Numerically estimate Hessian using Richardson extrapolation
##'
##' For N free parameters, Richardson extrapolation requires
##' (iterations * (N^2 + N)) function evaluations.
##' The implementation is closely based on the numDeriv R package.
##'
##' In addition to an estimate of the Hessian, forward, central, and
##' backward estimates of the gradient are made available in this
##' compute plan's output slot.
##'
##' When \code{checkGradient=TRUE}, the central difference estimate of
##' the gradient is used to determine whether the first order
##' convergence criterion is met. In addition, the forward and
##' backward difference estimates of the gradient are compared for
##' symmetry. When sufficient asymmetry is detected, the standard
##' error is flagged. In the case, profile likelihood confidence
##' intervals should be used for inference instead of standard errors
##' (see \code{mxComputeConfidenceInterval}).
##'
##' If provided, the square matrix \code{knownHessian} should have
##' dimnames set to the names of some subset of the free
##' parameters. Entries of the matrix set to NA will be estimated
##' numerically while entries containing finite values will be copied
##' to the Hessian result.
##'
##' @param freeSet names of matrices containing free variables
##' @param ...  Not used.  Forces remaining arguments to be specified by name.
##' @param fitfunction name of the fitfunction (defaults to 'fitfunction')
##' @param parallel whether to evaluate the fitfunction in parallel (defaults to TRUE)
##' @param stepSize starting set size (defaults to 0.0001)
##' @param iterations number of Richardson extrapolation iterations (defaults to 4L)
##' @template args-verbose
##' @param knownHessian an optional matrix of known Hessian entries
##' @param checkGradient whether to check the first order convergence criterion (gradient is near zero)
##' @param hessian whether to estimate the Hessian. If FALSE then only the gradient is estimated.
##' @param analytic Use the analytic Hessian, if available.
##' @aliases
##' MxComputeNumericDeriv-class
##' @examples
##' library(OpenMx)
##' data(demoOneFactor)
##' factorModel <- mxModel(name ="One Factor",
##' 	mxMatrix(type = "Full", nrow = 5, ncol = 1, free = FALSE, values = .2, name = "A"),
##' 	mxMatrix(type = "Symm", nrow = 1, ncol = 1, free = FALSE, values = 1 , name = "L"),
##' 	mxMatrix(type = "Diag", nrow = 5, ncol = 5, free = TRUE , values = 1 , name = "U"),
##' 	mxAlgebra(A %*% L %*% t(A) + U, name = "R"),
##' 	mxExpectationNormal(covariance = "R", dimnames = names(demoOneFactor)),
##' 	mxFitFunctionML(),
##' 	mxData(cov(demoOneFactor), type = "cov", numObs = 500),
##' 	mxComputeSequence(
##' 		list(mxComputeNumericDeriv(), mxComputeReportDeriv())
##' 	)
##' )
##' factorModelFit <- mxRun(factorModel)
##' factorModelFit$output$hessian

mxComputeNumericDeriv <- function(freeSet=NA_character_, ..., fitfunction='fitfunction',
				  parallel=TRUE,
				  stepSize=imxAutoOptionValue("Gradient step size"),
				  iterations=4L, verbose=0L,
				  knownHessian=NULL, checkGradient=TRUE, hessian=TRUE, analytic=TRUE)
{
  prohibitDotdotdot(list(...))
	verbose <- as.integer(verbose)
	iterations <- as.integer(iterations)

	if (!is.null(knownHessian)) {
		if (nrow(knownHessian) != ncol(knownHessian)) {
			stop("knownHessian must be square")
		}
		if (length(dimnames(knownHessian)) != 2 ||
		    any(rownames(knownHessian) != colnames(knownHessian))) {
			stop("knownHessian must have matching row and column names")
		}
	}

	new("MxComputeNumericDeriv", freeSet, fitfunction, parallel, stepSize, iterations,
	    verbose, knownHessian, checkGradient, hessian, analytic)
}

setMethod("displayCompute", signature(Ob="MxComputeNumericDeriv", indent="integer"),
	  function(Ob, indent) {
		  callNextMethod();
		  sp <- paste(rep('  ', indent), collapse="")
		  for (sl in c("fitfunction", "parallel", "stepSize", "iterations",
			       "verbose", "knownHessian", 'checkGradient', 'hessian')) {
			  slname <- paste("$", sl, sep="")
			  if (is.character(slot(Ob, sl))) {
				  cat(sp, slname, ":", omxQuotes(slot(Ob, sl)), '\n')
			  } else {
				  cat(sp, slname, ":", slot(Ob, sl), '\n')
			  }
		  }
		  invisible(Ob)
	  })

#----------------------------------------------------

setClass(Class = "MxComputeJacobian",
	contains = "BaseCompute",
	representation = representation(
		of = "MxCharOrNumber",
		data = "MxCharOrNumber",
		defvar.row = "integer"))

setMethod("initialize", "MxComputeJacobian",
	  function(.Object, freeSet, of, defvar.row, data) {
		  .Object@name <- 'compute'
		  .Object@.persist <- TRUE
		  .Object@freeSet <- freeSet
		  .Object@of <- of
		  .Object@data <- data
		  .Object@defvar.row <- defvar.row
		  .Object
	  })

setMethod("qualifyNames", signature("MxComputeJacobian"),
	function(.Object, modelname, namespace) {
		.Object <- callNextMethod()
		for (sl in c('of','data')) {
			slot(.Object, sl) <- imxConvertIdentifier(slot(.Object, sl), modelname, namespace)
		}
		.Object
	})

setMethod("convertForBackend", signature("MxComputeJacobian"),
	function(.Object, flatModel, model) {
		name <- .Object@name
		if (any(!is.integer(.Object@of))) {
			expNum <- match(.Object@of, names(flatModel@expectations))
			algNum <- match(.Object@of, names(flatModel@algebras))
			if (any(is.na(expNum)) && any(is.na(algNum))) {
				stop(paste("Can only apply MxComputeJacobian to MxAlgebra *or* MxExpectation not",
					   deparse(.Object@of)))
			}
			if (!any(is.na(expNum))) {
					# Usually negative numbers indicate matrices; not here
				.Object@of <- - expNum
			} else {
				.Object@of <- algNum - 1L
			}
		}
		if (any(!is.integer(.Object@data))) {
			if (is.na(.Object@defvar.row)) {
				.Object@data <- as.integer(NA)
			} else {
				dataNum <- match(.Object@data, names(flatModel@datasets))
				if (any(is.na(dataNum))) {
					stop(paste(class(.Object), omxQuotes(.Object@data),
						"not recognized as MxData"))
				}
				.Object@data <- dataNum - 1L
			}
		}
		.Object
	})

mxComputeJacobian <-
	function(freeSet=NA_character_, ..., of="expectation", defvar.row=as.integer(NA), data='data')
{
	new("MxComputeJacobian", freeSet, of, as.integer(defvar.row), data)
}

#----------------------------------------------------

setClass(Class = "MxComputeStandardError",
	contains = "BaseCompute",
	representation = representation(
		fitfunction = "MxCharOrNumber"
	))

setMethod("qualifyNames", signature("MxComputeStandardError"),
	function(.Object, modelname, namespace) {
		.Object <- callNextMethod();
		.Object@fitfunction <- imxConvertIdentifier(.Object@fitfunction, modelname, namespace)
		.Object
	})

setMethod("convertForBackend", signature("MxComputeStandardError"),
	function(.Object, flatModel, model) {
		name <- .Object@name
		if (is.character(.Object@fitfunction)) {
			.Object@fitfunction <- imxLocateIndex(flatModel, .Object@fitfunction, .Object)
		}
		.Object
	})

setMethod("initialize", "MxComputeStandardError",
	  function(.Object, freeSet, fitfunction) {
		  .Object@name <- 'compute'
		  .Object@.persist <- TRUE
		  .Object@freeSet <- freeSet
		  .Object@fitfunction <- fitfunction
		  .Object
	  })

##' Compute standard errors
##'
##' When the fit is in -2 log likelihood units, the SEs are derived
##' from the diagonal of the Hessian or inverse Hessian. The Hessian
##' (in some form) must already be available.
##'
##' If there are active MxConstraints and the fit is in -2logL units,
##' the SEs are derived from the Hessian and the Jacobian of the
##' constraint functions (see references).
##'
##' @param freeSet names of matrices containing free variables
##' @param fitfunction name of the fitfunction (defaults to 'fitfunction')
##' @aliases
##' MxComputeStandardError-class
##' @references
##' Moore T & Sadler B.  (2006).  \emph{Maximum-Likelihood Estimation and
##'      Scoring Under Parametric Constraints}.  Army Research Laboratory
##'      report ARL-TR-3805.
##' Schoenberg R.  (1997).  Constrained maximum likelihood.
##'      \emph{Computational Economics, 10}, p. 251-266.

mxComputeStandardError <- function(freeSet=NA_character_, fitfunction='fitfunction') {
	new("MxComputeStandardError", freeSet, fitfunction)
}

#----------------------------------------------------

setClass(Class = "MxComputeHessianQuality",
	 contains = "BaseCompute",
	 representation = representation(
	     verbose = "integer"))

setMethod("initialize", "MxComputeHessianQuality",
	  function(.Object, freeSet, verbose) {
		  .Object@name <- 'compute'
		  .Object@.persist <- TRUE
		  .Object@freeSet <- freeSet
		  .Object@verbose <- verbose
		  .Object
	  })

##' Compute the quality of the Hessian
##'
##' Tests whether the Hessian is positive definite
##' (model$output$infoDefinite) and, if so, computes the approximate condition
##' number (model$output$conditionNumber). See Luenberger & Ye (2008)
##' Second Order Test (p. 190) and Condition Number (p. 239).
##'
##' The condition number is approximated by \eqn{\mathrm{norm}(H) *
##' \mathrm{norm}(H^{-1})}{norm(H) * norm(solve(H))} where H is the
##' Hessian. The norm is either the 1- or infinity-norm (both obtain
##' the same result due to symmetry).
##'
##' @param freeSet names of matrices containing free variables
##' @param ...  Not used.  Forces remaining arguments to be specified by name.
##' @template args-verbose
##' @aliases
##' MxComputeHessianQuality-class
##' @references
##' Luenberger, D. G. & Ye, Y. (2008). Linear and nonlinear programming. Springer.

mxComputeHessianQuality <- function(freeSet=NA_character_, ..., verbose=0L) {
  prohibitDotdotdot(list(...))
	new("MxComputeHessianQuality", freeSet, as.integer(verbose))
}

#----------------------------------------------------

setClass(Class = "MxComputeReportDeriv",
	 contains = "BaseCompute")

setMethod("initialize", "MxComputeReportDeriv",
	  function(.Object, freeSet) {
		  .Object@name <- 'compute'
		  .Object@.persist <- TRUE
		  .Object@freeSet <- freeSet
		  .Object
	  })

##' Report derivatives
##'
##' Copy the internal gradient and Hessian back to R.
##'
##' @param freeSet names of matrices containing free variables
##' @aliases
##' MxComputeReportDeriv-class

mxComputeReportDeriv <- function(freeSet=NA_character_) {
	new("MxComputeReportDeriv", freeSet)
}

#----------------------------------------------------

setClass(Class = "MxComputeBootstrap",
	 contains = "BaseCompute",
	 representation = representation(
		 data = "MxCharOrNumber",
	     plan = "MxCompute",
	     replications = "integer",
	     verbose = "integer",
	     parallel = "logical",
	     OK = "ordered",
	     only = "integer"
	 ))

setMethod("initialize", "MxComputeBootstrap",
	  function(.Object, freeSet, data, plan, replications,
		   verbose, parallel, OK, only) {
		  .Object@name <- 'compute'
		  .Object@.persist <- TRUE
		  .Object@freeSet <- freeSet
		  .Object@data <- data
		  .Object@plan <- plan
		  .Object@replications <- replications
		  .Object@verbose <- verbose
		  .Object@parallel <- parallel
		  .Object@OK <- OK
		  .Object@only <- only
		  .Object
	  })

setMethod("getFreeVarGroup", signature("MxComputeBootstrap"),
	function(.Object) {
		result <- callNextMethod()
		for (step in c(.Object@plan)) {
			got <- getFreeVarGroup(step)
			if (length(got)) result <- append(result, got)
		}
		result
	})

setMethod("assignId", signature("MxComputeBootstrap"),
	function(.Object, id, defaultFreeSet) {
		.Object <- callNextMethod()
		defaultFreeSet <- .Object@freeSet
		id <- .Object@id
		for (sl in c('plan')) {
			slot(.Object, sl) <- assignId(slot(.Object, sl), id, defaultFreeSet)
			id <- slot(.Object, sl)@id + 1L
		}
		.Object@id <- id
		.Object
	})

setMethod("qualifyNames", signature("MxComputeBootstrap"),
	function(.Object, modelname, namespace) {
		.Object <- callNextMethod()
		for (sl in c('data')) {
			slot(.Object, sl) <- imxConvertIdentifier(slot(.Object, sl), modelname, namespace)
		}
		for (sl in c('plan')) {
			slot(.Object, sl) <- qualifyNames(slot(.Object, sl), modelname, namespace)
		}
		.Object
	})

setMethod("convertForBackend", signature("MxComputeBootstrap"),
	function(.Object, flatModel, model) {
		name <- .Object@name
		if (any(!is.integer(.Object@data))) {
			dataNum <- match(.Object@data, names(flatModel@datasets))
			if (any(is.na(dataNum))) {
				stop(paste("MxComputeBootstrap:", omxQuotes(.Object@data),
					   "not recognized as MxData"))
			}
			.Object@data <- dataNum - 1L
		}
		for (sl in c('plan')) {
			slot(.Object, sl) <- convertForBackend(slot(.Object, sl), flatModel, model)
		}
		.Object
	})

setMethod("updateFromBackend", signature("MxComputeBootstrap"),
	function(.Object, computes) {
		.Object <- callNextMethod()
		for (sl in c('plan')) {
			slot(.Object, sl) <- updateFromBackend(slot(.Object, sl), computes)
		}
		.Object
	})

mxComputeBootstrap <- function(data, plan, replications=200, ...,
			       verbose=0L, parallel=TRUE, freeSet=NA_character_,
			       OK=c("OK", "OK/green"), only=NA_integer_) {
  prohibitDotdotdot(list(...))
	data <- vapply(data, function(e1) {
		path <- unlist(strsplit(e1, imxSeparatorChar, fixed = TRUE))
		if (length(path) == 1) {
			e1 <- paste(path, "data", sep=imxSeparatorChar)
		}
		e1
	}, "")

	new("MxComputeBootstrap", freeSet, data, plan, as.integer(replications),
	    as.integer(verbose), parallel, as.statusCode(OK), only)
}

setMethod("displayCompute", signature(Ob="MxComputeBootstrap", indent="integer"),
	  function(Ob, indent) {
		  callNextMethod();
		  sp <- paste(rep('  ', indent), collapse="")
		  cat(sp, "$plan :", '\n')
		  displayCompute(Ob@plan, indent+1L)
		  for (sl in c("data", "replications", "OK",
			       "verbose", "parallel")) {
			  slname <- paste("$", sl, sep="")
			  if (is.character(slot(Ob, sl))) {
				  cat(sp, slname, ":", omxQuotes(slot(Ob, sl)), '\n')
			  } else {
				  cat(sp, slname, ":", slot(Ob, sl), '\n')
			  }
		  }
		  invisible(Ob)
	  })

#----------------------------------------------------

setClass(Class = "MxComputeReportExpectation",
	 contains = "BaseCompute")

setMethod("initialize", "MxComputeReportExpectation",
	  function(.Object, freeSet) {
		  .Object@name <- 'compute'
		  .Object@.persist <- TRUE
		  .Object@freeSet <- freeSet
		  .Object
	  })

##' Report expectation
##'
##' Copy the internal model expectations back to R.
##'
##' @param freeSet names of matrices containing free variables
##' @aliases
##' MxComputeReportExpectation-class

mxComputeReportExpectation <- function(freeSet=NA_character_) {
	new("MxComputeReportExpectation", freeSet)
}

#----------------------------------------------------

setClass(Class = "MxComputeSetOriginalStarts",
	 contains = "BaseCompute")

setMethod("initialize", "MxComputeSetOriginalStarts",
	  function(.Object, freeSet) {
		  .Object@name <- 'compute'
		  .Object@.persist <- TRUE
		  .Object@freeSet <- freeSet
		  .Object
	  })

##' Reset parameter starting values
##'
##' Sets the current parameter vector back to the original starting values.
##'
##' @param freeSet names of matrices containing free variables
##' @aliases
##' MxComputeSetOriginalStarts-class

mxComputeSetOriginalStarts <- function(freeSet=NA_character_) {
	new("MxComputeSetOriginalStarts", freeSet)
}

#----------------------------------------------------

setClass(Class = "MxComputeGenerateData",
	 contains = "BaseCompute",
	 representation = representation(
		 expectation = "MxCharOrNumber"
	 ))

setMethod("initialize", "MxComputeGenerateData",
	  function(.Object, expectation) {
		  .Object@name <- 'compute'
		  .Object@.persist <- TRUE
		  .Object@freeSet <- NA_character_
		  .Object@expectation <- expectation
		  .Object
	  })

setMethod("qualifyNames", signature("MxComputeGenerateData"),
	function(.Object, modelname, namespace) {
		.Object <- callNextMethod();
		.Object@expectation <- imxConvertIdentifier(.Object@expectation, modelname, namespace)
		.Object
	})

setMethod("convertForBackend", signature("MxComputeGenerateData"),
	function(.Object, flatModel, model) {
		name <- .Object@name
		if (is.character(.Object@expectation)) {
			.Object@expectation <- imxLocateIndex(flatModel, .Object@expectation, .Object)
		}
		.Object
	})

##' Generate data
##'
##' Generate data specified by the model expectations.
##'
##' @param expectation a character vector of expectations to generate data for
##' @aliases
##' MxComputeGenerateData-class

mxComputeGenerateData <- function(expectation='expectation') {
	new("MxComputeGenerateData", expectation)
}

#----------------------------------------------------

setClass(Class = "MxComputeLoadData",
	 contains = "BaseCompute",
	 representation = representation(
		 dest = "MxCharOrNumber",
		 column = "MxOptionalChar",
		 path = "MxOptionalChar",
		 originalDataIsIndexOne = "logical",
		 byrow = "logical",
		 row.names = "MxOptionalInteger",
		 col.names = "MxOptionalInteger",
		 verbose = "integer",
		 cacheSize = "integer",
		 method = "character",
		 checkpointMetadata = "logical",
		 skip.rows = "integer",
		 skip.cols = "integer",
		 na.strings = "character",
		 observed = "MxOptionalDataFrame",
     rowFilter = "MxOptionalLogical"
	 ))

setMethod("initialize", "MxComputeLoadData",
	function(.Object, dest, column, path, originalDataIsIndexOne,
		 row.names, col.names, skip.rows, skip.cols, byrow, verbose,
		 cacheSize, method, checkpointMetadata, na.strings, observed, rowFilter) {
		  .Object@name <- 'compute'
		  .Object@.persist <- TRUE
		  .Object@freeSet <- NA_character_
		  .Object@dest <- dest
		  .Object@column <- column
		  .Object@path <- path
		  .Object@originalDataIsIndexOne <- originalDataIsIndexOne
		  .Object@byrow <- byrow
		  .Object@row.names <- row.names
		  .Object@col.names <- col.names
		  .Object@verbose <- verbose
		  .Object@cacheSize <- cacheSize
		  .Object@method <- method
		  .Object@checkpointMetadata <- checkpointMetadata
		  .Object@skip.rows <- skip.rows
		  .Object@skip.cols <- skip.cols
		  .Object@na.strings <- na.strings
		  .Object@observed <- observed
      .Object@rowFilter <- rowFilter
		  .Object
	  })

setMethod("convertForBackend", signature("MxComputeLoadData"),
	function(.Object, flatModel, model) {
		name <- .Object@name
    if (length(.Object@column) == 0) .Object@dest <- -1L
		if (is.character(.Object@dest)) {
			full <- grepl('.', .Object@dest, fixed=TRUE)
			for (dx in 1:length(.Object@dest)) {
				if (full[dx]) next
				.Object@dest[dx] <- paste0(.Object@dest[dx], '.data')
			}
			.Object@dest <- imxLocateIndex(flatModel, .Object@dest, .Object)
		}
		.Object
	})

##' Load columns into an MxData object
##'
##' \lifecycle{experimental}
##'
##' The purpose of this compute step is to help quickly perform many
##' similar analyses. For example, if we are given a sample of people
##' with a few million SNPs (single-nucleotide polymorphism) per
##' person then we could fit a separate model for each SNP by iterating
##' over the SNP data.
##'
##' The column names given in the \code{column} parameter must already
##' exist in the model's MxData object. Pre-existing data is assumed to be
##' a placeholder and is not used unless
##' \code{originalDataIsIndexOne} is set to TRUE.
##'
##' For \code{method='csv'}, the highest performance arrangement is
##' \code{byrow=TRUE} because entire columns are stored in single
##' chunks (rows) on the disk and can be easily loaded. For
##' \code{byrow=FALSE}, the data requires transposition. To load a
##' single column of observed data, it is necessary to read through
##' the whole file. This can be slow for large files. To amortize the
##' cost of transposition, \code{cacheSize} columns are loaded on
##' every pass through the file.
##'
##' After \code{mxRun} returns, the \code{dest} mxData object will
##' contain the most recently loaded data. Hence, any single analysis
##' of a series can be reproduced by issuing \code{mxComputeLoadData}
##' with the single index associated with a particular dataset,
##' replacing the compute plan with something like
##' \code{omxDefaultComputePlan}, and then passing the model back
##' through \code{mxRun}. This can be a helpful approach when
##' investigating unexpected results.
##'
##' @param dest the name of the model where the columns will be loaded
##' @param column a character vector. The column names to replace.
##' @param method name of the conduit used to load the columns.
##' @param ...  Not used.  Forces remaining arguments to be specified by name.
##' @param path the path to the file containing the data
##' @param originalDataIsIndexOne logical. Whether to use the initial data for index 1
##' @param byrow logical. Whether the data columns are stored in rows.
##' @param row.names optional integer. Column containing the row names.
##' @param col.names optional integer. Row containing the column names.
##' @param skip.rows integer. Number of rows to skip before reading data.
##' @param skip.cols integer. Number of columns to skip before reading data.
##' @template args-verbose
##' @param cacheSize integer. How many columns to cache per
##' scan through the data. Only used when byrow=FALSE.
##' @param checkpointMetadata logical. Whether to add per record metadata to the checkpoint
##' @param na.strings character vector. A vector of strings that denote a missing value.
##' @param observed data frame. The reservoir of data for \code{method='data.frame'}.
##' @param rowFilter logical vector. Whether to skip the source row.
##' @aliases
##' MxComputeLoadData-class
##' @seealso
##' \link{mxComputeLoadMatrix}, \link{mxComputeCheckpoint}, \link{mxRun}, \link{omxDefaultComputePlan}
mxComputeLoadData <- function(dest, column, method=c('csv', 'data.frame'), ..., path=c(),
			      originalDataIsIndexOne=FALSE, byrow=TRUE,
			      row.names=c(), col.names=c(),
			      skip.rows=0, skip.cols=0,
			      verbose=0L,
			      cacheSize=100L, checkpointMetadata=TRUE, na.strings=c('NA'),
			      observed=NULL, rowFilter=c()) {
  prohibitDotdotdot(list(...))
	if (cacheSize < 1L) stop("cacheSize must be a positive integer")
	new("MxComputeLoadData", dest, column, path, originalDataIsIndexOne,
		as.integer(row.names), as.integer(col.names),
		as.integer(skip.rows), as.integer(skip.cols), byrow,
		as.integer(verbose), as.integer(cacheSize), method,
		as.logical(checkpointMetadata), as.character(na.strings), observed,
    as.logical(rowFilter))
}

#----------------------------------------------------

setClass(Class = "MxComputeLoadContext",
	 contains = "BaseCompute",
	 representation = representation(
		 method = "character",
		 path = "character",
		 column = "integer",
	   sep = "character",
	   verbose = "integer",
	   header = "logical",
	   col.names = "character"
	 ))

setMethod("initialize", "MxComputeLoadContext",
	function(.Object, method, path, column, sep, verbose, header, col.names) {
		  .Object@name <- 'compute'
		  .Object@.persist <- TRUE
		  .Object@freeSet <- NA_character_
		  .Object@method <- method
		  .Object@path <- path
		  .Object@column <- column
		  .Object@sep <- sep
		  .Object@verbose <- verbose
		  .Object@header <- header
		  .Object@col.names <- col.names
		  .Object
	  })

##' Load contextual data to supplement checkpoint
##'
##' \lifecycle{experimental}
##'
##' Currently, this only supports comma separated value format and no
##' row names. If \code{header=TRUE} and \code{col.names} are
##' provided, the \code{col.names} take precedence. If
##' \code{header=FALSE} and no \code{col.names} are provided then
##' the column names consist of the file name and column offset.
##'
##' An \code{originalDataIsIndexOne} option is not offered. You'll need to
##' add an extra line at the start on your file if you wish to make
##' use of \code{originalDataIsIndexOne} in \code{mxComputeLoad*}.
##'
##' @param method name of the conduit used to load the columns.
##' @param path the path to the file containing the data
##' @param column a character vector. The column names to log.
##' @param ...  Not used.  Forces remaining arguments to be specified by name.
##' @param sep the field separator character. Values on each line of the file are separated by this character.
##' @param header logical. Whether the first row contains column headers.
##' @param col.names character vector. Column names
##' @template args-verbose
##' @aliases
##' MxComputeLoadContext-class
##' @seealso
##' \link{mxComputeCheckpoint}, \link{mxComputeLoadData}, \link{mxComputeLoadMatrix}
mxComputeLoadContext <- function(method=c('csv'), path=c(), column, ..., sep=' ',
				 verbose=0L, header=TRUE, col.names=NULL) {
  prohibitDotdotdot(list(...))
	method <- match.arg(method)
	if (length(column) > 1 && any(diff(column) < 0))
	  stop("Columns must be ordered from left to right")
	if (length(col.names) && length(col.names) != length(column)) {
	  stop("If col.names provided, the length must match the number of columns")
	}
	new("MxComputeLoadContext", method, as.character(path), as.integer(column), sep,
	    as.integer(verbose), as.logical(header), as.character(col.names))
}

#----------------------------------------------------

setClass(Class = "MxComputeLoadMatrix",
	 contains = "BaseCompute",
	 representation = representation(
		 dest = "MxCharOrNumber",
		 method = "character",
		 path = "MxOptionalChar",
		 originalDataIsIndexOne = "logical",
		 row.names = "logical",
		 col.names = "logical",
		 observed = "MxOptionalDataFrame"
	 ))

setMethod("initialize", "MxComputeLoadMatrix",
	function(.Object, dest, method, path, originalDataIsIndexOne, row.names, col.names,
		 observed) {
		  .Object@name <- 'compute'
		  .Object@.persist <- TRUE
		  .Object@freeSet <- NA_character_
		  .Object@dest <- dest
		  .Object@method <- method
		  .Object@path <- path
		  .Object@originalDataIsIndexOne <- originalDataIsIndexOne
		  .Object@row.names <- row.names
		  .Object@col.names <- col.names
		  .Object@observed <- observed
		  .Object
	  })

setMethod("qualifyNames", signature("MxComputeLoadMatrix"),
	function(.Object, modelname, namespace) {
		.Object <- callNextMethod();
		.Object@dest <- imxConvertIdentifier(.Object@dest, modelname, namespace)
		.Object
	})

setMethod("convertForBackend", signature("MxComputeLoadMatrix"),
	function(.Object, flatModel, model) {
		name <- .Object@name
		if (any(is.character(.Object@dest))) {
			.Object@dest <- sapply(.Object@dest,
				function(dd) imxLocateIndex(flatModel, dd, .Object))
		}
		.Object
	})

mxComputeLoadMatrix <- function(dest, method=c('csv','data.frame'), ..., path=NULL,
				originalDataIsIndexOne=FALSE,
				row.names=FALSE, col.names=FALSE, observed=NULL) {
  prohibitDotdotdot(list(...))
	method <- match.arg(method)
	new("MxComputeLoadMatrix", dest, method, path, originalDataIsIndexOne,
		as.logical(row.names), as.logical(col.names), observed)
}

#----------------------------------------------------

setClass(Class = "MxComputeCheckpoint",
	 contains = "BaseCompute",
	 representation = representation(
		 what = "MxCharOrNumber",
		 toReturn = "logical",
		 path = "MxOptionalChar",
		 append = "logical",
		 header = "logical",
		 log = "MxListOrNull",
		 parameters = "logical",
		 loopIndices = "logical",
		 fit = "logical",
		 counters = "logical",
		 status = "logical",
		 standardErrors = "logical",
		 gradient = "logical",
		 vcov = "logical",
		 vcovWLS = "logical",
     vcovFilter = "MxOptionalChar",
     sampleSize = "logical",
     useVcovFilter = "logical"
	 ))

setMethod("initialize", "MxComputeCheckpoint",
	function(.Object, what, path, append, header, toReturn, parameters,
           loopIndices, fit, counters, status, standardErrors, gradient, vcov,
           vcovFilter, sampleSize, vcovWLS, useVcovFilter) {
		  .Object@name <- 'compute'
		  .Object@.persist <- TRUE
		  .Object@freeSet <- NA_character_
		  .Object@what <- what
		  .Object@path <- path
		  .Object@append <- append
		  .Object@header <- header
		  .Object@toReturn <- toReturn
		  .Object@parameters <- parameters
		  .Object@loopIndices <- loopIndices
		  .Object@fit <- fit
		  .Object@counters <- counters
		  .Object@status <- status
		  .Object@standardErrors <- standardErrors
		  .Object@gradient <- gradient
		  .Object@vcov <- vcov
		  .Object@vcovFilter <- vcovFilter
		  .Object@sampleSize <- sampleSize
		  .Object@vcovWLS <- vcovWLS
		  .Object@useVcovFilter <- useVcovFilter
		  .Object
	  })

setMethod("qualifyNames", signature("MxComputeCheckpoint"),
	function(.Object, modelname, namespace) {
		.Object <- callNextMethod()
		for (sl in c('what')) {
			slot(.Object, sl) <- imxConvertIdentifier(slot(.Object, sl), modelname, namespace)
		}
		.Object
	})

setMethod("convertForBackend", signature("MxComputeCheckpoint"),
	function(.Object, flatModel, model) {
		name <- .Object@name
		if (is.character(.Object@what)) {
			.Object@what <- imxLocateIndex(flatModel, .Object@what, .Object)
		}
		.Object
	})

#' Log parameters and state to disk or memory
#'
#' @param what a character vector of algebra names to include in each checkpoint
#' @template args-dots-barrier
#' @param path a character vector of where to write the checkpoint file
#' @param append if FALSE, truncates the checkpoint file upon open. If TRUE, existing data is preserved and checkpoints are appended.
#' @param header whether to write the header that describes the content of each column
#' @param toReturn logical. Whether to store the checkpoint in memory and return it after the model is run
#' @param parameters logical. Whether to include the parameter vector
#' @param loopIndices logical. Whether to include the loop indices
#' @param fit logical. Whether to include the fit value
#' @param counters logical. Whether to include counters (number of evaluations and iterations)
#' @param status logical. Whether to include the status code
#' @param standardErrors logical. Whether to include the standard errors
#' @param gradient logical. Whether to include the gradients
#' @param vcov logical. Whether to include the vcov in half-vectorized order
#' @param vcovWLS logical. Whether to include the vcov from WLS residualizing regressions in half-vectorized order
#' @param vcovFilter character vector. Vector of parameters indicating
#'   which parameter covariances to include. Only the variance is
#'   included for those parameters not mentioned.
#' @param useVcovFilter logical. Whether to use the vcovFilter (TRUE) or include all entries (FALSE)
#' @param sampleSize logical. Whether to include the sample size of the mxData. \lifecycle{experimental}
#'
#' @description
#' Captures the current state of the backend. When \code{path} is set, the
#' state is written to disk in a single row. When \code{toReturn} is set,
#' the state is recorded in memory and returned after \code{mxRun}.
#'
#' @aliases
#' MxComputeCheckpoint-class
#' @seealso
#' \code{\link{mxComputeLoadData}}, \code{\link{mxComputeLoadMatrix}},
#' \code{\link{mxComputeLoadContext}}, \code{\link{mxComputeLoop}}
#'
#' @family model state
#' @examples
#' library(OpenMx)
#'
#' m1 <- mxModel(
#'   "poly22", # Eqn 22 from Tsallis & Stariolo (1996)
#'   mxMatrix(type='Full', values=runif(4, min=-1e6, max=1e6),
#'            ncol=1, nrow=4, free=TRUE, name='x'),
#'   mxAlgebra(sum((x*x-8)^2) + 5*sum(x) + 57.3276, name="fit"),
#'   mxFitFunctionAlgebra('fit'))
#'
#' plan <- mxComputeLoop(list(
#'   mxComputeSetOriginalStarts(),
#'     mxComputeSimAnnealing(method="tsallis1996",
#'                           control=list(tempEnd=1)),
#'     mxComputeCheckpoint(path = "result.log")),
#'   i=1:4)
#'
#' m1 <- mxRun(mxModel(m1, plan)) # see the file 'result.log'
mxComputeCheckpoint <- function(what=NULL, ..., path=NULL, append=FALSE, header=TRUE, toReturn=FALSE,
				parameters=TRUE, loopIndices=TRUE, fit=TRUE, counters=TRUE,
				status=TRUE, standardErrors=FALSE, gradient=FALSE, vcov=FALSE,
        vcovFilter=c(), sampleSize=FALSE, vcovWLS=FALSE, useVcovFilter=FALSE) {
  prohibitDotdotdot(list(...))
	what <- as.character(what)
	path <- as.character(path)
	new("MxComputeCheckpoint", what, path, as.logical(append), as.logical(header), as.logical(toReturn),
		as.logical(parameters), as.logical(loopIndices), as.logical(fit), as.logical(counters),
		as.logical(status), as.logical(standardErrors), as.logical(gradient), as.logical(vcov),
    vcovFilter, as.logical(sampleSize), as.logical(vcovWLS), as.logical(useVcovFilter))
}

#----------------------------------------------------

setClass(Class = "MxComputeSequence",
	 contains = "ComputeSteps",
	 representation = representation(
	     independent="logical"
	     ))

setMethod("initialize", "MxComputeSequence",
	  function(.Object, steps, freeSet, independent) {
		  .Object@name <- 'compute'
		  .Object@.persist <- TRUE
		  .Object@steps <- steps
		  .Object@freeSet <- freeSet
		  .Object@independent <- independent
		  .Object
	  })

##' Invoke a series of compute objects in sequence
##'
##' @param steps a list of compute objects
##' @param ... Not used; forces argument 'freeSet' to be specified by name.
##' @param freeSet Names of matrices containing free parameters.
##' @param independent Whether the steps could be executed out-of-order.
##' @aliases
##' MxComputeSequence-class
mxComputeSequence <- function(steps=list(), ..., freeSet=NA_character_, independent=FALSE) {
  prohibitDotdotdot(list(...))
	new("MxComputeSequence", steps=steps, freeSet, independent)
}

setClass(Class = "MxComputeDefault",
	 contains = "BaseCompute")

setMethod("initialize", "MxComputeDefault",
	  function(.Object, freeSet) {
		  .Object@name <- 'compute'
		  .Object@.persist <- TRUE
		  .Object@freeSet <- freeSet
		  .Object
	  })

##' Default compute plan
##'
##' This is an empty placeholder for the default compute plan.
##' To create an actual plan, use \link{omxDefaultComputePlan}.
##'
##' @param freeSet names of matrices containing free variables
##' @aliases
##' MxComputeDefault-class
mxComputeDefault <- function(freeSet=NA_character_) {
	new("MxComputeDefault", freeSet)
}

##' Compute nothing
##'
##' Note that this compute plan actually does nothing whereas
##' \code{mxComputeOnce("expectation", "nothing")} may remove the
##' prediction of an expectation.
##'
mxComputeNothing <- function() {
	mxComputeSequence(freeSet=c())
}

setMethod("displayCompute", signature(Ob="MxComputeSequence", indent="integer"),
	  function(Ob, indent) {
		  callNextMethod();
		  sp <- paste(rep('  ', indent), collapse="")
		  cat(sp, "independent :", Ob@independent, '\n')
		  stepName <- paste("'", names(Ob@steps), "'",sep='')
		  if (length(stepName) != length(Ob@steps)) stepName <- 1:length(Ob@steps)
		  if (length(Ob@steps)) for (step in 1:length(Ob@steps)) {
			  cat(sp, "steps[[", stepName[step], "]] :", '\n')
			  displayCompute(Ob@steps[[step]], indent+1L)
		  }
		  invisible(Ob)
	  })

convertComputes <- function(flatModel, model) {
	if (is.null(flatModel@compute)) return()
	convertForBackend(flatModel@compute, flatModel, model)
}

updateModelCompute <- function(model, computes) {
	if (is.null(model@compute)) return()
	updateFromBackend(model@compute, computes)
}

##' Sparse symmetric matrix invert
##'
##' This API is visible to permit testing. Please do not use.
##'
##' @param mat the matrix to invert
imxSparseInvert <- function(mat) .Call(sparseInvert_wrapper, mat)

ProcessCheckHess <- function(model, checkHess) {
  if (omxHasDefaultComputePlan(model)) {
    optList <- options()$mxOption
    if (is.na(checkHess)) checkHess <- FALSE
    checkHessYes <- ifelse(checkHess, 'Yes', 'No')
    optList[['Calculate Hessian']] <- checkHessYes
    optList[['Standard Errors']] <- checkHessYes
    model <- mxModel(model, omxDefaultComputePlan(optionList=optList))
  } else {
    if (!is.na(checkHess)) {
      stop(paste("Model", model$name, "has a custom compute plan.",
                 "Cannot act on checkHess=",checkHess))
    }
  }
  model
}

##' omxHasDefaultComputePlan
##'
##' Determine whether the model has a default complete plan (i.e., not custom).
##'
##' @param model model
omxHasDefaultComputePlan <- function(model) {
	if (is.null(model@compute)) return(TRUE)
	if (!.hasSlot(model@compute, '.persist') || !model@compute@.persist) return(TRUE)
	FALSE
}

omxDefaultComputePlan <- function(modelName=NULL, intervals=FALSE, useOptimizer=TRUE,
				  optionList=options()$mxOption, penaltySearch=FALSE) {
	if(length(modelName) && !is.character(modelName[1])){stop("argument 'modelName' must be a character string")}
	compute <- NULL
	fitNum <- ifelse(length(modelName), paste(modelName, 'fitfunction', sep="."), "fitfunction")
	if (!useOptimizer) {
		compute <- mxComputeSequence(list(CO=mxComputeOnce(from=fitNum, 'fit', .is.bestfit=TRUE),
																			RE=mxComputeReportExpectation()))
  } else{
    if (penaltySearch) {
      steps <- list(PS=mxComputePenaltySearch(plan=mxComputeSequence(list(
        SV=mxComputeSetOriginalStarts(),
        GD=mxComputeGradientDescent(fitfunction=fitNum)))))
    } else {
      steps <- list(GD=mxComputeGradientDescent(fitfunction=fitNum, verbose=0L))
    }
    if (intervals){
      ciOpt <- mxComputeGradientDescent(
        verbose=0L,
        fitfunction=fitNum,
        nudgeZeroStarts=FALSE)
      cType <- ciOpt$defaultCImethod
      if (cType == 'ineq') {
        ciOpt <- mxComputeTryHard(plan=ciOpt, scale=0.05)
      }
      steps <- c(steps, CI=mxComputeConfidenceInterval(
        fitfunction=fitNum,
        constraintType=cType,
        verbose=0L, plan=ciOpt))
    }
    if (optionList[["Calculate Hessian"]] == "Yes") {
      steps <- c(steps, ND=mxComputeNumericDeriv(
        fitfunction=fitNum,
        stepSize=imxAutoOptionValue('Gradient step size',optionList)))
    }
    if (optionList[["Standard Errors"]] == "Yes") {
      steps <- c(steps, SE=mxComputeStandardError(), HQ=mxComputeHessianQuality())
    }
    compute <- mxComputeSequence(c(steps,
                                   RD=mxComputeReportDeriv(),
                                   RE=mxComputeReportExpectation()))
	}
	compute@.persist <- TRUE
	return(compute)
}