You are here

Revision of HOWTO Make an Objective Function from Sun, 07/10/2011 - 15:44

Revisions allow you to track differences between multiple versions of your content, and revert back to older versions.

Making An Objective Function

This howto provides step-by-step instructions on writing your own objective function for OpenMx. You will first want to prototype your objective function using either an MxAlgebraObjective or an MxRObjective. An algebra objective consists of an arbitrary number of matrix algebra statements that ultimately evaluate to a 1 x 1 matrix. An R objective consists of an arbitrary R function that evaluates a model to return a numeric value. Algebra objectives are computationally faster than R objectives, but algebra objectives are more limited in their expressive power. When you are satisfied prototyping your new objective function, we will guide you through the process of writing the front-end for this objective function in R, and writing the back-end for this objective function in C.

Step 1: Prototyping

Using MxAlgebraObjective

Can you express your objective function using the matrix operators and functions provided by OpenMx? If so, then create an MxAlgebraObjective. The following example is from the SimpleAlgebra demo.

A <- mxMatrix("Full", values = c(1, 2), nrow=1, ncol=2, name="A")
B <- mxMatrix("Full", values = c(3, 4), nrow = 2, ncol = 1, free = c(TRUE, FALSE), name="B")
algebra <- mxAlgebra(A %*% B - 3, "algebra")
algSquared <- mxAlgebra(algebra %*% algebra, "algebraSquared")
 
model <- mxModel("themodel")
model <- mxModel(model, A, B)
model <- mxModel(model, algebra, algSquared)
model <- mxModel(model, mxAlgebraObjective("algebraSquared"))
 
# Test 1: Algebra is a multiply and subtract
 
modelOut <- mxRun(model)

Using MxRObjective

Otherwise, express your objective function in the R front-end. Write a function in R that accepts two arguments. The first argument is the mxModel that should be evaluated, and the second argument is some persistent state information that can be stored between one iteration of optimization to the next iteration. It is valid for the function to simply ignore the second argument. The function must return either a single numeric value, or a list of exactly two elements. If the function returns a list, the first argument must be a single numeric value and the second element will be the new persistent state information to be passed into this function at the next iteration. The single numeric value will be used by the optimizer to perform optimization. For more help on the slots in the MxModel class, see the help provided with the OpenMx package (type '?MxModel' in R). The following example is from the SimpleRObjective demo.

A <- mxMatrix(nrow = 2, ncol = 2, values = c(1:4), free = TRUE, name = 'A')
 
squared <- function(x) { x ^ 2 }
 
objFunction <- function(model, state) {
    values <- mxEval(A, model) 
    return(squared(values[1,1] - 4) + squared(values[1,2] - 3) +
        squared(values[2,1] - 2) + squared(values[2,2] - 1))
}
objective <- mxRObjective(objFunction)
 
model <- mxModel('model', A, objective)
 
modelOut <- mxRun(model)

Step 2: Write Your Objective Function

Step 2a: Implement the front-end

In order to write your own objective function, you will need to create your own objective function class in R. For teaching purposes, we are going to use the MxAlgebraObjective class as a guide. First, your class needs to be a sub-class of the MxBaseObjective class. The MxBaseObjective class contains three slots: 'name' which is a character string, 'data' which is either a character string or a number, and 'result' which is a 1x1 matrix that stores the result of optimization. Your class should contain any additional slots that you need to pass along to the optimizer. In the MxAlgebraObjective, we have an additional slot that will inform the optimizer as to which MxAlgebra to use for optimization.

setClass(Class = "MxAlgebraObjective",
    representation = representation(
        algebra = "MxCharOrNumber"),
    contains = "MxBaseObjective")

Next, you will need to overwrite five generic methods for your new objective function class. These six methods are initialize, genericObjModelConvert, genericObjFunNamespace, genericObjFunConvert, genericObjDependencies, and genericObjRename.

initialize

The initialize method is a constructor method. It will never be directly called by the user, but it will be critical for making new instances of your class. When writing your initialize method, it is required that the default value for the data slot is as.integer(NA) and the default value for the name slot is 'objective'. The data slot will only be modified by the genericObjFunConvert method, and the name slot will only be modified by the genericObjFunNamespace method.

setMethod("initialize", "MxAlgebraObjective",
    function(.Object, algebra, 
        data = as.integer(NA), name = 'objective') {
        .Object@name <- name
        .Object@algebra <- algebra
        .Object@data <- data
        return(.Object)
    }
)

genericObjModelConvert

The genericObjModelConvert method is used in the event that you need to perform transformations on the model that are specific to your objective function. One example use-case of this method is to transform RAM objective functions to FIML objective functions, if raw data is provided to the RAM objective function. During that particular transformation, several additional matrices and algebras are added to the model, and the RAM objective function is replaced with a FIML objective function. For now, we will show you the default behavior of this function, which is the most common behavior. You do not need to copy/paste the default behavior into your objective function, it is inherited from "MxBaseObjective".The function genericObjModelConvert accepts five parameters: '.Object' is the current objective function, 'job' is the entire job that was passed to the mxRun() function, 'model' is the individual model associated with the current objective function, 'namespace' is the model namespace, and 'flatJob' is a flattened version of 'job'.  This function MUST return the job. Note: As of OpenMx 1.2, this function must return list(job, flatJob). This change was implemented to improve performance. If you set job@.newobjects, job@.newobjective, and job@.newtree to FALSE (see next paragraph), then you must manually update flatJob with any changes you make to job. Otherwise, you do not need to update flatJob as the job will be re-flattened.
The function is responsible for updating three private slots in the job object. All three of the slots contain boolean values. job@.newobjects should report whether additional matrices or algebras were added anywhere in the model. This triggers a regeneration of the namespace. job@.newobjective should report whether the current objective function was transformed into another type of objective function. This triggers the genericObjModelConvert method on the new objective function. job@.newtree should report whether any new submodels were introduced in the model. This triggers genericObjModelConvert on the new submodels.

setMethod("genericObjModelConvert", "MxBaseObjective",
    function(.Object, job, model, namespace, flatJob) {
        job@.newobjects <- FALSE
        job@.newobjective <- FALSE
        job@.newtree <- FALSE
        return(job) # or return(list(job, flatJob)) in the new API
})

genericObjFunNamespace The genericObjFunNamespace method is used to transform an objective function from its local model namespace, into the flattened model namespace. You must use the function omxConvertIdentifier on any slot that is storing name information. Remember to always invoke omxConvertIdentifier on the 'name' slot and the 'data' slot. Shown below is the implementation for the FIML objective function.

setMethod("genericObjFunNamespace", signature("MxFIMLObjective"), 
    function(.Object, modelname, namespace) {
        .Object@name <- omxIdentifier(modelname, .Object@name)
        .Object@covariance <- omxConvertIdentifier(.Object@covariance, 
            modelname, namespace)
        .Object@means <- omxConvertIdentifier(.Object@means, 
            modelname, namespace)
        .Object@data <- omxConvertIdentifier(.Object@data, 
            modelname, namespace)
        return(.Object)
})

genericObjFunConvert

The genericObjFunConvert method is invoked to transform an objective function from a form that interacts with the R front-end, to a form that is suitable for processing by the C back-end. Any character strings in the objective function are converted into matrix/algebra pointers for the back-end to process. You should read the documentation for the function omxLocateIndex to familiarise yourself with matrix/algebra pointers. Two arguments are passed into this function, a flattened version of the model, and the original model itself. When at all possible, it is recommended that you interact with the flattened version of the model. Below we have included the genericObjFunConvert method for both the algebra objective and the FIML objective. The FIML objective has been included so that you may see how the data slot is processed. NOTE: .Object@algebra is a string on input, and is converted into a matrix/algebra pointer.

setMethod("genericObjFunConvert", signature("MxAlgebraObjective"), 
    function(.Object, flatModel, model) {
        name <- .Object@name
        algebra <- .Object@algebra
        if (is.na(algebra)) {
            modelname <- omxReverseIdentifier(model, .Object@name)[[1]]
            msg <- paste("The algebra name cannot be NA",
            "for the algebra objective of model", omxQuotes(modelname))
            stop(msg, call. = FALSE)
        }
        algebraIndex <- omxLocateIndex(flatModel, algebra, name)
        .Object@algebra <- algebraIndex
        return(.Object)
})

NOTE: .Object@data, .Object@means, and .Object@covariance are all strings at the beginning of this function, and are converted into data pointers or matrix/algebra pointers.

setMethod("genericObjFunConvert", signature("MxFIMLObjective"), 
    function(.Object, flatModel, model) {
        modelname <- omxReverseIdentifier(model, .Object@name)[[1]]
        name <- .Object@name
        if(is.na(.Object@data)) {
            msg <- paste("The FIML objective",
                "does not have a dataset associated with it in model",
                omxQuotes(modelname))
            stop(msg, call.=FALSE)
        }
        mxDataObject <- flatModel@datasets[[.Object@data]]
        if (mxDataObject@type != 'raw') {
            msg <- paste("The dataset associated with the FIML objective", 
                "in model", omxQuotes(modelname), "is not raw data.")
            stop(msg, call.=FALSE)
        }
        verifyObservedNames(mxDataObject@observed, mxDataObject@type, flatModel, modelname, "FIML")
        checkNumericData(mxDataObject)
        meansName <- .Object@means
        covName <- .Object@covariance
        dataName <- .Object@data
        threshNames <- .Object@thresholds
        .Object@means <- omxLocateIndex(flatModel, .Object@means, name)
        .Object@covariance <- omxLocateIndex(flatModel, .Object@covariance, name)
        .Object@data <- omxLocateIndex(flatModel, .Object@data, name)
        verifyExpectedNames(covName, meansName, flatModel, modelname, "FIML")
        .Object@definitionVars <- generateDefinitionList(flatModel)
        .Object@dataColumns <- generateDataColumns(flatModel, covName, dataName)
        .Object@thresholdColumns <- generateThresholdColumns(flatModel, threshNames, dataName)
        if (length(mxDataObject@observed) == 0) {
            .Object@data <- as.integer(NA)
        }
        return(.Object)
})

genericObjDependencies

The OpenMx library performs cycle detection to determine if the algebraic relationships in a model contain any cycles. You must specify the dependencies of your objective function type. Below we see how the dependencies for the RAM objective function are specified.

setMethod("genericObjDependencies", signature("MxRAMObjective"), 
    function(.Object, dependencies) {
    sources <- c(.Object@A, .Object@S, .Object@F, .Object@M, .Object@thresholds)
    sources <- sources[!is.na(sources)]
    dependencies <- omxAddDependency(sources, .Object@name, dependencies)
    return(dependencies)
})

genericObjRename

This method is invoked by the mxRename() function whenever an objective function is encountered. Renaming allows an existing model name to be transformed into a new name, and all identifiers using the original model name are transformed into the new name. A typical need for renaming pieces of the objective function concern any definition variables that are named "oldname.dataset.defvar" must be renamed to "newname.dataset.defvar".

setMethod("genericObjRename", signature("MxFIMLObjective"), 
    function(.Object, oldname, newname) {
    .Object@covariance <- renameReference(.Object@covariance, oldname, newname)
    .Object@means <- renameReference(.Object@means, oldname, newname)
    .Object@data <- renameReference(.Object@data, oldname, newname)
    .Object@thresholds <- sapply(.Object@thresholds, renameReference, oldname, newname)
    return(.Object)
})

User-interface

Finally, you should write your own method that creates instances of your new objective function class. This interface will invoke the 'initialize' constructor that you wrote earlier, but should also include some basic sanity-checking before invoking the constructor. Here is the algebra objective function interface. Since this method does not have access to any parts of the model, you will not be able to perform full error-checking at this stage. Full error-checking should be performed in the genericObjFunConvert method.

mxAlgebraObjective <- function(algebra) {
    if (missing(algebra) || typeof(algebra) != "character") {
        stop("Algebra argument is not a string (the name of the algebra)")
    }
    return(new("MxAlgebraObjective", algebra))
}

mxRun() Workflow

The following is the workflow of the mxRun() function, showing where each objective function method occurs in the workflow.

  1. check namespace for errors
  2. translate objectives (genericObjModelConvert)
  3. run independent sub-models
  4. flatten model (genericObjFunNamespace)
  5. check algebras for errors (genericObjDependencies)
  6. convert global variables and constants in algebra expressions
  7. convert model pieces to back-end (genericObjFunConvert)
  8. invoke the back-end
  9. repopulate the model results from the optimizer

Step 2(b): Implement the back-end

TODO: write this section