Performance improvements to parallel executions
[openmx:openmx.git] / R / MxRun.R
1 #
2 #   Copyright 2007-2010 The OpenMx Project
3 #
4 #   Licensed under the Apache License, Version 2.0 (the "License");
5 #   you may not use this file except in compliance with the License.
6 #   You may obtain a copy of the License at
7
8 #        http://www.apache.org/licenses/LICENSE-2.0
9
10 #   Unless required by applicable law or agreed to in writing, software
11 #   distributed under the License is distributed on an "AS IS" BASIS,
12 #   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13 #   See the License for the specific language governing permissions and
14 #   limitations under the License.
15
16 omxLapply <- function(x, fun, ...) {
17         libraries <- search()
18         if ("package:snowfall" %in% libraries) {
19                 return(sfLapply(x, fun, ...))
20         } else {
21                 return(lapply(x, fun, ...))
22         }
23 }
24
25 omxSapply <- function(x, fun, ...) {
26         libraries <- search()
27         if ("package:snowfall" %in% libraries) {
28                 return(sfSapply(x, fun, ...))
29         } else {
30                 return(sapply(x, fun, ...))
31         }
32 }
33
34 mxRun <- function(model, silent = FALSE) {
35         frontendStart <- Sys.time()
36         if(!silent) cat("Running", model@name, "\n")
37         namespace <- omxGenerateNamespace(model)
38         omxCheckNamespace(model, namespace)
39         omxCheckMatrices(model)
40         omxVerifyModel(model)
41         model <- translateObjectives(model, namespace)
42         # Regenerate the namespace
43         namespace <- omxGenerateNamespace(model)
44         dshare <- shareData(model)
45         independents <- getAllIndependents(dshare)
46         independents <- omxLapply(independents, mxRun, silent)
47         independents <- lapply(independents, omxFreezeModel)
48         model <- omxReplaceModels(model, independents)
49         if(is.null(model@objective) && 
50                         length(model@matrices) == 0 && 
51                         length(model@algebras) == 0 &&
52                         length(omxDependentModels(model)) == 0) {
53                 frontendStop1 <- Sys.time()
54                 backendStop <- frontendStop1
55                 frontendStop2 <- frontendStop1
56                 model@output <- calculateTiming(model@output, frontendStart,
57                         frontendStop1, backendStop, frontendStop2)              
58                 return(model)
59         }
60         model <- convertSquareBracketLabels(model)
61         flatModel <- omxFlattenModel(model, namespace)
62         data <- convertDatasets(flatModel, model)
63         freeFixedValues <- omxCheckVariables(flatModel, namespace)
64         flatModel <- populateDefInitialValues(flatModel)
65         oldFlatModel <- flatModel
66         flatModel <- convertAlgebras(flatModel, list(startvals=freeFixedValues, 
67                 values=namespace$values, parameters=namespace$parameters))
68         cycleDetection(flatModel)
69         flatModel <- checkEvaluation(model, flatModel, oldFlatModel)
70         parameters <- generateParameterList(flatModel)
71         matrices <- generateMatrixList(flatModel)
72         algebras <- generateAlgebraList(flatModel)
73         startVals <- generateValueList(matrices, parameters)
74         objectives <- convertObjectives(flatModel, model)
75         algebras <- append(algebras, objectives)
76         constraints <- convertConstraints(flatModel)
77         state <- c()
78         objective <- getObjectiveIndex(flatModel)
79         options <- generateOptionsList(model@options)
80         frontendStop1 <- Sys.time()
81         output <- .Call("callNPSOL", objective, startVals, 
82                 constraints, matrices, parameters, 
83                 algebras, data, options, state, PACKAGE = "OpenMx")
84         backendStop <- Sys.time()
85         model <- updateModelMatrices(model, flatModel, output$matrices)
86         model <- updateModelAlgebras(model, flatModel, output$algebras)
87         model <- undoSquareBracketLabels(model)
88         model@output <- processOptimizerOutput(flatModel, names(matrices),
89                 names(algebras), names(parameters), output)
90         frontendStop2 <- Sys.time()
91         model@output <- calculateTiming(model@output, frontendStart,
92                 frontendStop1, backendStop, frontendStop2)      
93         return(model)
94 }
95
96 processOptimizerOutput <- function(flatModel, matrixNames, 
97                 algebraNames, parameterNames, output) {
98         output$mxVersion <- mxVersion()
99         if (length(output$estimate) == length(parameterNames)) {
100                 names(output$estimate) <- parameterNames
101         }
102         if (length(output$gradient) == length(parameterNames)) {
103                 names(output$gradient) <- parameterNames
104         }
105         output$hessian <- t(output$hessianCholesky) %*% output$hessianCholesky
106         if (nrow(output$hessian) == length(parameterNames) &&
107                 ncol(output$hessian) == length(parameterNames)) {
108                 dimnames(output$hessian) <- list(parameterNames, parameterNames)
109         }
110         if (length(output$matrices) == length(matrixNames)) {
111                 names(output$matrices) <- matrixNames
112         }
113         if (length(output$algebras) == length(algebraNames)) {
114                 names(output$algebras) <- algebraNames
115         }
116     if (output$status[[1]] > 0) {
117         npsolWarnings(flatModel@name, output$status[[1]])
118     } else if (output$status[[1]] < 0) {
119         stop(paste("The job for model", omxQuotes(flatModel@name),
120             "exited abnormally with the error message:",
121             output$status[[3]]), call. = FALSE)
122     }
123         return(output)
124 }
125
126 calculateTiming <- function(output, frontendStart,
127         frontendStop1, backendStop, frontendStop2) {
128         output$frontendTime <- (frontendStop1 - frontendStart) +
129                 (frontendStop2 - backendStop)
130         output$backendTime <- (backendStop - frontendStop1)
131         return(output)
132 }
133
134 npsolMessages <- list('1' = paste('The final iterate satisfies',
135                 'the optimality conditions to the accuracy requested,',
136                 'but the sequence of iterates has not yet converged.',
137                 'NPSOL was terminated because no further improvement',
138                 'could be made in the merit function (Mx status GREEN).'),
139                 '2' = paste('The linear constraints and bounds could not be satisfied.',
140                 'The problem has no feasible solution.'),
141                 '3' = paste('The nonlinear constraints and bonuds could not be satisfied.',
142                 'The problem may have no feasible solution.'),
143                 '4' = 'The Major iteration limit was reached (Mx status BLUE).',
144                 '6' = paste('model does not satisfy the first-order optimality conditions',
145                 'to the required accuracy, and no improved point for the',
146                 'merit function could be found during the final linesearch (Mx status RED)'),
147                 '7' = paste('The function derivates returned by funcon or funobj',
148                 'appear to be incorrect.'),
149                 '9' = 'An input parameter was invalid')
150
151 npsolWarnings <- function(name, status) {
152         message <- npsolMessages[[as.character(status)]]
153         if(!is.null(message)) {
154                 warning(paste("In model", omxQuotes(name), 
155                         "NPSOL returned a non-zero status code", 
156                         paste(status, '.', sep = ''), message), call. = FALSE)
157         }
158 }