Reorder initialization so free parameters are available for expectations
[openmx:openmx.git] / src / npsolWrap.c
1 /*
2  *  Copyright 2007-2013 The OpenMx Project
3  *
4  *  Licensed under the Apache License, Version 2.0 (the "License");
5  *  you may not use this file except in compliance with the License.
6  *  You may obtain a copy of the License at
7  *
8  *       http://www.apache.org/licenses/LICENSE-2.0
9  *
10  *   Unless required by applicable law or agreed to in writing, software
11  *   distributed under the License is distributed on an "AS IS" BASIS,
12  *   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13  *  See the License for the specific language governing permissions and
14  *  limitations under the License.
15  */
16
17 #include <R.h>
18 #include <Rinternals.h>
19 #include <Rdefines.h>
20 #include <R_ext/Rdynload.h>
21 #include <R_ext/BLAS.h>
22 #include <R_ext/Lapack.h>
23 #include "omxDefines.h"
24 #include "npsolWrap.h"
25 #include "omxOpenmpWrap.h"
26
27 #include <stdio.h>
28 #include <sys/types.h>
29 #include <errno.h>
30 #include "omxState.h"
31 #include "omxGlobalState.h"
32 #include "omxMatrix.h"
33 #include "omxAlgebra.h"
34 #include "omxFitFunction.h"
35 #include "omxExpectation.h"
36 #include "omxNPSOLSpecific.h"
37 #include "omxImportFrontendState.h"
38 #include "omxExportBackendState.h"
39 #include "omxHessianCalculation.h"
40 #include "omxOptimizer.h"
41
42 omp_lock_t GlobalRLock;
43
44 static R_CallMethodDef callMethods[] = {
45         {"omxBackend", (DL_FUNC) omxBackend, 12},
46         {"omxCallAlgebra", (DL_FUNC) omxCallAlgebra, 3},
47         {"findIdenticalRowsData", (DL_FUNC) findIdenticalRowsData, 5},
48         {NULL, NULL, 0}
49 };
50
51 void R_init_OpenMx(DllInfo *info) {
52         R_registerRoutines(info, NULL, callMethods, NULL, NULL);
53
54         omx_omp_init_lock(&GlobalRLock);
55
56         // There is no code that will change behavior whether openmp
57         // is set for nested or not. I'm just keeping this in case it
58         // makes a difference with older versions of openmp. 2012-12-24 JNP
59 #if defined(_OPENMP) && _OPENMP <= 200505
60         omp_set_nested(0);
61 #endif
62 }
63
64 void R_unload_OpenMx(DllInfo *info) {
65         omx_omp_destroy_lock(&GlobalRLock);
66 }
67
68 /* Main functions */
69 SEXP omxCallAlgebra(SEXP matList, SEXP algNum, SEXP options) {
70
71         if(OMX_DEBUG) { Rprintf("-----------------------------------------------------------------------\n");}
72         if(OMX_DEBUG) { Rprintf("Explicit call to algebra %d.\n", INTEGER(algNum));}
73
74         int j,k,l;
75         omxMatrix* algebra;
76         int algebraNum = INTEGER(algNum)[0];
77         SEXP ans, nextMat;
78         char output[250];
79         int errOut = 0;
80
81         /* Create new omxState for current state storage and initialize it. */
82         
83         globalState = (omxState*) R_alloc(1, sizeof(omxState));
84         omxInitState(globalState, NULL, 1);
85         if(OMX_DEBUG) { Rprintf("Created state object at 0x%x.\n", globalState);}
86
87         /* Retrieve All Matrices From the MatList */
88
89         if(OMX_DEBUG) { Rprintf("Processing %d matrix(ces).\n", length(matList));}
90         globalState->numMats = length(matList);
91         globalState->matrixList = (omxMatrix**) R_alloc(length(matList), sizeof(omxMatrix*));
92
93         for(k = 0; k < length(matList); k++) {
94                 PROTECT(nextMat = VECTOR_ELT(matList, k));      // This is the matrix + populations
95                 globalState->matrixList[k] = omxNewMatrixFromRPrimitive(nextMat, globalState, 1, - k - 1);
96                 if(OMX_DEBUG) {
97                         Rprintf("Matrix initialized at 0x%0xd = (%d x %d).\n",
98                                 globalState->matrixList[k], globalState->matrixList[k]->rows, globalState->matrixList[k]->cols);
99                 }
100                 UNPROTECT(1); // nextMat
101         }
102
103         algebra = omxNewAlgebraFromOperatorAndArgs(algebraNum, globalState->matrixList, globalState->numMats, globalState);
104
105         if(algebra==NULL) {
106                 error(globalState->statusMsg);
107         }
108
109         if(OMX_DEBUG) {Rprintf("Completed Algebras and Matrices.  Beginning Initial Compute.\n");}
110         omxStateNextEvaluation(globalState);
111
112         omxRecompute(algebra);
113
114         PROTECT(ans = allocMatrix(REALSXP, algebra->rows, algebra->cols));
115         for(l = 0; l < algebra->rows; l++)
116                 for(j = 0; j < algebra->cols; j++)
117                         REAL(ans)[j * algebra->rows + l] =
118                                 omxMatrixElement(algebra, l, j);
119
120         UNPROTECT(1);   /* algebra */
121
122         if(OMX_DEBUG) { Rprintf("All Algebras complete.\n"); }
123
124         if(globalState->statusCode != 0) {
125                 errOut = globalState->statusCode;
126                 strncpy(output, globalState->statusMsg, 250);
127         }
128
129         omxFreeAllMatrixData(algebra);
130         omxFreeState(globalState);
131
132         if(errOut != 0) {
133                 error(output);
134         }
135
136         return ans;
137 }
138
139 SEXP omxBackend(SEXP fitfunction, SEXP startVals, SEXP constraints,
140         SEXP matList, SEXP varList, SEXP algList, SEXP expectList,
141         SEXP data, SEXP intervalList, SEXP checkpointList, SEXP options, SEXP state) {
142
143         /* Helpful variables */
144
145         int errOut = 0;                 // Error state: Clear
146
147         SEXP nextLoc;
148
149         int n;
150         int calculateStdErrors = FALSE;
151         int numHessians = 0;
152         int ciMaxIterations = 5;
153         int disableOptimizer = 0;
154         int numThreads = 1;
155         int analyticGradients = 0;
156
157         /* Sanity Check and Parse Inputs */
158         /* TODO: Need to find a way to account for nullness in these.  For now, all checking is done on the front-end. */
159 //      if(!isVector(startVals)) error ("startVals must be a vector");
160 //      if(!isVector(matList)) error ("matList must be a list");
161 //      if(!isVector(algList)) error ("algList must be a list");
162
163         /*      Set NPSOL options */
164         omxSetNPSOLOpts(options, &numHessians, &calculateStdErrors, 
165                 &ciMaxIterations, &disableOptimizer, &numThreads, 
166                 &analyticGradients, length(startVals));
167
168         /* Create new omxState for current state storage and initialize it. */
169         globalState = (omxState*) R_alloc(1, sizeof(omxState));
170         omxInitState(globalState, NULL, numThreads);
171         globalState->numFreeParams = length(startVals);
172         globalState->analyticGradients = analyticGradients;
173         if(OMX_DEBUG) { Rprintf("Created state object at 0x%x.\n", globalState);}
174
175         /* Retrieve Data Objects */
176         if(!errOut) errOut = omxProcessMxDataEntities(data);
177     
178         /* Retrieve All Matrices From the MatList */
179         if(!errOut) errOut = omxProcessMxMatrixEntities(matList);
180
181         globalState->numAlgs = length(algList);
182         globalState->markMatrices = (int*) R_alloc(globalState->numMats + globalState->numAlgs, sizeof(int));
183         
184         if (length(startVals) != length(varList)) error("varList and startVals must be the same length");
185
186         /* Process Free Var List */
187         omxProcessFreeVarList(varList);
188
189         /* Initialize all Expectations Here */
190         if(!errOut) errOut = omxProcessMxExpectationEntities(expectList);
191
192         if(!errOut) {
193                 omxProcessMxAlgebraEntities(algList);
194                 errOut = globalState->statusMsg[0];
195         }
196
197         /* Complete Expectations */
198         if(!errOut) errOut = omxCompleteMxExpectationEntities();
199
200         if(!errOut) {
201                 // This is the chance to check for matrix
202                 // conformability, etc.  Any errors encountered should
203                 // be reported using R's error() function, not
204                 // omxRaiseErrorf.
205
206                 // disable parallelism until omxDuplicateState() can be invoked
207                 globalState->numChildren = 0;
208
209                 omxInitialMatrixAlgebraCompute();
210
211                 globalState->numChildren = (numThreads > 1) ? numThreads : 0;
212                 omxResetStatus(globalState);
213         }
214
215         if(!errOut && !isNull(fitfunction)) {
216                 if(OMX_DEBUG) { Rprintf("Processing fit function.\n"); }
217                 globalState->fitMatrix = omxMatrixLookupFromState1(fitfunction, globalState);
218                 errOut = globalState->statusMsg[0];
219         }
220         
221         // TODO: Make calculateHessians an option instead.
222
223         if(errOut) error(globalState->statusMsg);
224
225         int numChildren = globalState->numChildren;
226
227         /* Process Matrix and Algebra Population Function */
228         /*
229           Each matrix is a list containing a matrix and the other matrices/algebras that are
230           populated into it at each iteration.  The first element is already processed, above.
231           The rest of the list will be processed here.
232         */
233         for(int j = 0; j < globalState->numMats; j++) {
234                 PROTECT(nextLoc = VECTOR_ELT(matList, j));              // This is the matrix + populations
235                 omxProcessMatrixPopulationList(globalState->matrixList[j], nextLoc);
236                 UNPROTECT(1);
237         }
238
239         /* Processing Constraints */
240         omxProcessConstraints(constraints);
241
242         /* Process Confidence Interval List */
243         omxProcessConfidenceIntervals(intervalList);
244
245         /* Process Checkpoint List */
246         omxProcessCheckpointOptions(checkpointList);
247
248         for(int i = 0; i < numChildren; i++) {
249                 omxDuplicateState(globalState->childList[i], globalState);
250         }
251
252         n = globalState->numFreeParams;
253
254         SEXP minimum, estimate, gradient, hessian;
255         PROTECT(minimum = NEW_NUMERIC(1));
256         PROTECT(estimate = allocVector(REALSXP, n));
257         PROTECT(gradient = allocVector(REALSXP, n));
258         PROTECT(hessian = allocMatrix(REALSXP, n, n));
259
260         if (n>0) { memcpy(REAL(estimate), REAL(startVals), sizeof(double)*n); }
261         
262         omxInvokeNPSOL(REAL(minimum), REAL(estimate), REAL(gradient), REAL(hessian), disableOptimizer);
263
264         SEXP code, status, statusMsg, iterations;
265         SEXP evaluations, ans=NULL, names=NULL, algebras, matrices, expectations, optimizer;
266         SEXP intervals, NAmat, intervalCodes, calculatedHessian, stdErrors;
267
268         int numReturns = 14;
269
270         PROTECT(code = NEW_NUMERIC(1));
271         PROTECT(status = allocVector(VECSXP, 3));
272         PROTECT(iterations = NEW_NUMERIC(1));
273         PROTECT(evaluations = NEW_NUMERIC(2));
274         PROTECT(matrices = NEW_LIST(globalState->numMats));
275         PROTECT(algebras = NEW_LIST(globalState->numAlgs));
276         PROTECT(expectations = NEW_LIST(globalState->numExpects));
277
278         PROTECT(optimizer = allocVector(VECSXP, 2));
279         PROTECT(calculatedHessian = allocMatrix(REALSXP, n, n));
280         PROTECT(stdErrors = allocMatrix(REALSXP, n, 1)); // for optimizer
281         PROTECT(names = allocVector(STRSXP, 2)); // for optimizer
282         PROTECT(intervals = allocMatrix(REALSXP, globalState->numIntervals, 2)); // for optimizer
283         PROTECT(intervalCodes = allocMatrix(INTSXP, globalState->numIntervals, 2)); // for optimizer
284         PROTECT(NAmat = allocMatrix(REALSXP, 1, 1)); // In case of missingness
285         REAL(NAmat)[0] = R_NaReal;
286
287         omxSaveState(globalState, REAL(estimate), REAL(minimum)[0]);
288
289         /* Fill in details from the optimizer */
290         SET_VECTOR_ELT(optimizer, 0, gradient);
291         SET_VECTOR_ELT(optimizer, 1, hessian);
292
293         SET_STRING_ELT(names, 0, mkChar("minimum"));
294         SET_STRING_ELT(names, 1, mkChar("estimate"));
295         namesgets(optimizer, names);
296
297         REAL(code)[0] = globalState->inform;
298         REAL(iterations)[0] = globalState->iter;
299         REAL(evaluations)[0] = globalState->computeCount;
300
301         /* Fill Status code. */
302         SET_VECTOR_ELT(status, 0, code);
303         PROTECT(code = NEW_NUMERIC(1));
304         REAL(code)[0] = globalState->statusCode;
305         SET_VECTOR_ELT(status, 1, code);
306         PROTECT(statusMsg = allocVector(STRSXP, 1));
307         SET_STRING_ELT(statusMsg, 0, mkChar(globalState->statusMsg));
308         SET_VECTOR_ELT(status, 2, statusMsg);
309
310         if(numHessians && globalState->fitMatrix != NULL && globalState->optimumStatus >= 0) {          // No hessians or standard errors if the optimum is invalid
311                 if(globalState->numConstraints == 0) {
312                         if(OMX_DEBUG) { Rprintf("Calculating Hessian for Fit Function.\n");}
313                         int gotHessians = omxEstimateHessian(numHessians, .0001, 4, globalState);
314                         if(gotHessians) {
315                                 if(calculateStdErrors) {
316                                         for(int j = 0; j < numHessians; j++) {          //TODO: Fix Hessian calculation to allow more if requested
317                                                 if(OMX_DEBUG) { Rprintf("Calculating Standard Errors for Fit Function.\n");}
318                                                 omxFitFunction* oo = globalState->fitMatrix->fitFunction;
319                                                 if(oo->getStandardErrorFun != NULL) {
320                                                         oo->getStandardErrorFun(oo);
321                                                 } else {
322                                                         omxCalculateStdErrorFromHessian(2.0, oo);
323                                                 }
324                                         }
325                                 }
326                         } else {
327                                 numHessians = 0;
328                         }
329                 } else {
330                         numHessians = 0;
331                 }
332         } else {
333                 numHessians = 0;
334         }
335
336         /* Likelihood-based Confidence Interval Calculation */
337         if(globalState->numIntervals) {
338                 omxNPSOLConfidenceIntervals(REAL(minimum), REAL(estimate), REAL(gradient), REAL(hessian), ciMaxIterations);
339         }  
340
341         handleFreeVarList(globalState, globalState->optimalValues, n);  // Restore to optima for final compute
342         if(!errOut) omxFinalAlgebraCalculation(globalState, matrices, algebras, expectations); 
343
344         omxPopulateFitFunction(globalState, numReturns, &ans, &names);
345
346         if(numHessians) {
347                 omxPopulateHessians(numHessians, globalState->fitMatrix, 
348                         calculatedHessian, stdErrors, calculateStdErrors, n);
349         }
350
351         if(globalState->numIntervals) { // Populate CIs
352                 omxPopulateConfidenceIntervals(globalState, intervals, intervalCodes);
353         }
354         
355         REAL(evaluations)[1] = globalState->computeCount;
356
357         int nextEl = 0;
358
359         SET_STRING_ELT(names, nextEl++, mkChar("minimum"));
360         SET_STRING_ELT(names, nextEl++, mkChar("estimate"));
361         SET_STRING_ELT(names, nextEl++, mkChar("gradient"));
362         SET_STRING_ELT(names, nextEl++, mkChar("hessianCholesky"));
363         SET_STRING_ELT(names, nextEl++, mkChar("status"));
364         SET_STRING_ELT(names, nextEl++, mkChar("iterations"));
365         SET_STRING_ELT(names, nextEl++, mkChar("evaluations"));
366         SET_STRING_ELT(names, nextEl++, mkChar("matrices"));
367         SET_STRING_ELT(names, nextEl++, mkChar("algebras"));
368         SET_STRING_ELT(names, nextEl++, mkChar("expectations"));
369         SET_STRING_ELT(names, nextEl++, mkChar("confidenceIntervals"));
370         SET_STRING_ELT(names, nextEl++, mkChar("confidenceIntervalCodes"));
371         SET_STRING_ELT(names, nextEl++, mkChar("calculatedHessian"));
372         SET_STRING_ELT(names, nextEl++, mkChar("standardErrors"));
373
374         nextEl = 0;
375
376         SET_VECTOR_ELT(ans, nextEl++, minimum);
377         SET_VECTOR_ELT(ans, nextEl++, estimate);
378         SET_VECTOR_ELT(ans, nextEl++, gradient);
379         SET_VECTOR_ELT(ans, nextEl++, hessian);
380         SET_VECTOR_ELT(ans, nextEl++, status);
381         SET_VECTOR_ELT(ans, nextEl++, iterations);
382         SET_VECTOR_ELT(ans, nextEl++, evaluations);
383         SET_VECTOR_ELT(ans, nextEl++, matrices);
384         SET_VECTOR_ELT(ans, nextEl++, algebras);
385         SET_VECTOR_ELT(ans, nextEl++, expectations);
386         SET_VECTOR_ELT(ans, nextEl++, intervals);
387         SET_VECTOR_ELT(ans, nextEl++, intervalCodes);
388         if(numHessians == 0) {
389                 SET_VECTOR_ELT(ans, nextEl++, NAmat);
390         } else {
391                 SET_VECTOR_ELT(ans, nextEl++, calculatedHessian);
392         }
393         if(!calculateStdErrors) {
394                 SET_VECTOR_ELT(ans, nextEl++, NAmat);
395         } else {
396                 SET_VECTOR_ELT(ans, nextEl++, stdErrors);
397         }
398         namesgets(ans, names);
399
400         if(OMX_VERBOSE) {
401                 Rprintf("Inform Value: %d\n", globalState->optimumStatus);
402                 Rprintf("--------------------------\n");
403         }
404
405         /* Free data memory */
406         omxFreeState(globalState);
407
408         UNPROTECT(numReturns);                                          // Unprotect Output Parameters
409         UNPROTECT(8);                                                           // Unprotect internals
410
411         if(OMX_DEBUG) {Rprintf("All vectors freed.\n");}
412
413         return(ans);
414
415 }
416