Revert "Refrain from duplicating the model unless required by the fitfunction"
[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 PROTECT_INDEX omxProtectSave()
69 {
70         PROTECT_INDEX initialpix;
71         PROTECT_WITH_INDEX(R_NilValue, &initialpix);
72         UNPROTECT(1);
73         return initialpix;
74 }
75
76 void omxProtectRestore(PROTECT_INDEX initialpix)
77 {
78         PROTECT_INDEX pix;
79         PROTECT_WITH_INDEX(R_NilValue, &pix);
80         PROTECT_INDEX diff = pix - initialpix;
81         if (diff < 0) {
82                 warning("omxProtectRestore diff=%d; you forgot to PROTECT something", diff);
83                 diff = 0;
84         }
85         UNPROTECT(1 + diff);
86 }
87
88 /* Main functions */
89 SEXP omxCallAlgebra(SEXP matList, SEXP algNum, SEXP options) {
90
91         PROTECT_INDEX initialpix = omxProtectSave();
92
93         if(OMX_DEBUG) { Rprintf("-----------------------------------------------------------------------\n");}
94         if(OMX_DEBUG) { Rprintf("Explicit call to algebra %d.\n", INTEGER(algNum));}
95
96         int j,k,l;
97         omxMatrix* algebra;
98         int algebraNum = INTEGER(algNum)[0];
99         SEXP ans, nextMat;
100         char output[250];
101         int errOut = 0;
102
103         /* Create new omxState for current state storage and initialize it. */
104         
105         globalState = (omxState*) R_alloc(1, sizeof(omxState));
106         omxInitState(globalState, NULL, 1);
107         if(OMX_DEBUG) { Rprintf("Created state object at 0x%x.\n", globalState);}
108
109         /* Retrieve All Matrices From the MatList */
110
111         if(OMX_DEBUG) { Rprintf("Processing %d matrix(ces).\n", length(matList));}
112         globalState->numMats = length(matList);
113         globalState->matrixList = (omxMatrix**) R_alloc(length(matList), sizeof(omxMatrix*));
114
115         for(k = 0; k < length(matList); k++) {
116                 PROTECT(nextMat = VECTOR_ELT(matList, k));      // This is the matrix + populations
117                 globalState->matrixList[k] = omxNewMatrixFromRPrimitive(nextMat, globalState, 1, - k - 1);
118                 if(OMX_DEBUG) {
119                         Rprintf("Matrix initialized at 0x%0xd = (%d x %d).\n",
120                                 globalState->matrixList[k], globalState->matrixList[k]->rows, globalState->matrixList[k]->cols);
121                 }
122                 UNPROTECT(1); // nextMat
123         }
124
125         algebra = omxNewAlgebraFromOperatorAndArgs(algebraNum, globalState->matrixList, globalState->numMats, globalState);
126
127         if(algebra==NULL) {
128                 error(globalState->statusMsg);
129         }
130
131         if(OMX_DEBUG) {Rprintf("Completed Algebras and Matrices.  Beginning Initial Compute.\n");}
132         omxStateNextEvaluation(globalState);
133
134         omxRecompute(algebra);
135
136         PROTECT(ans = allocMatrix(REALSXP, algebra->rows, algebra->cols));
137         for(l = 0; l < algebra->rows; l++)
138                 for(j = 0; j < algebra->cols; j++)
139                         REAL(ans)[j * algebra->rows + l] =
140                                 omxMatrixElement(algebra, l, j);
141
142         UNPROTECT(1);   /* algebra */
143
144         if(OMX_DEBUG) { Rprintf("All Algebras complete.\n"); }
145
146         if(globalState->statusCode != 0) {
147                 errOut = globalState->statusCode;
148                 strncpy(output, globalState->statusMsg, 250);
149         }
150
151         omxFreeAllMatrixData(algebra);
152         omxFreeState(globalState);
153
154         if(errOut != 0) {
155                 error(output);
156         }
157
158         omxProtectRestore(initialpix);
159
160         return ans;
161 }
162
163 SEXP omxBackend(SEXP fitfunction, SEXP startVals, SEXP constraints,
164         SEXP matList, SEXP varList, SEXP algList, SEXP expectList,
165         SEXP data, SEXP intervalList, SEXP checkpointList, SEXP options, SEXP state) {
166
167         /* Helpful variables */
168
169         int errOut = 0;                 // Error state: Clear
170
171         SEXP nextLoc;
172
173         int n;
174         int calculateStdErrors = FALSE;
175         int numHessians = 0;
176         int ciMaxIterations = 5;
177         int disableOptimizer = 0;
178         int numThreads = 1;
179         int analyticGradients = 0;
180
181         /* Sanity Check and Parse Inputs */
182         /* TODO: Need to find a way to account for nullness in these.  For now, all checking is done on the front-end. */
183 //      if(!isVector(startVals)) error ("startVals must be a vector");
184 //      if(!isVector(matList)) error ("matList must be a list");
185 //      if(!isVector(algList)) error ("algList must be a list");
186
187         PROTECT_INDEX initialpix = omxProtectSave();
188
189         /*      Set NPSOL options */
190         omxSetNPSOLOpts(options, &numHessians, &calculateStdErrors, 
191                 &ciMaxIterations, &disableOptimizer, &numThreads, 
192                 &analyticGradients, length(startVals));
193
194         /* Create new omxState for current state storage and initialize it. */
195         globalState = (omxState*) R_alloc(1, sizeof(omxState));
196         omxInitState(globalState, NULL, numThreads);
197         globalState->numFreeParams = length(startVals);
198         globalState->analyticGradients = analyticGradients;
199         if(OMX_DEBUG) { Rprintf("Created state object at 0x%x.\n", globalState);}
200
201         /* Retrieve Data Objects */
202         if(!errOut) errOut = omxProcessMxDataEntities(data);
203     
204         /* Retrieve All Matrices From the MatList */
205         if(!errOut) errOut = omxProcessMxMatrixEntities(matList);
206
207         globalState->numAlgs = length(algList);
208         globalState->markMatrices = (int*) R_alloc(globalState->numMats + globalState->numAlgs, sizeof(int));
209         
210         if (length(startVals) != length(varList)) error("varList and startVals must be the same length");
211
212         /* Process Free Var List */
213         omxProcessFreeVarList(varList);
214
215         /* Initialize all Expectations Here */
216         if(!errOut) errOut = omxProcessMxExpectationEntities(expectList);
217
218         if(!errOut) {
219                 omxProcessMxAlgebraEntities(algList);
220                 errOut = globalState->statusMsg[0];
221         }
222
223         /* Complete Expectations */
224         if(!errOut) errOut = omxCompleteMxExpectationEntities();
225
226         if(!errOut) {
227                 // This is the chance to check for matrix
228                 // conformability, etc.  Any errors encountered should
229                 // be reported using R's error() function, not
230                 // omxRaiseErrorf.
231
232                 // disable parallelism until omxDuplicateState() can be invoked
233                 globalState->numChildren = 0;
234
235                 omxInitialMatrixAlgebraCompute();
236
237                 globalState->numChildren = (numThreads > 1) ? numThreads : 0;
238                 omxResetStatus(globalState);
239         }
240
241         if(!errOut) {
242                 omxProcessFitFunction(fitfunction);
243                 errOut = globalState->statusMsg[0];
244         }
245         
246         // TODO: Make calculateHessians an option instead.
247
248         if(errOut) error(globalState->statusMsg);
249
250         int numChildren = globalState->numChildren;
251
252         /* Process Matrix and Algebra Population Function */
253         /*
254           Each matrix is a list containing a matrix and the other matrices/algebras that are
255           populated into it at each iteration.  The first element is already processed, above.
256           The rest of the list will be processed here.
257         */
258         for(int j = 0; j < globalState->numMats; j++) {
259                 PROTECT(nextLoc = VECTOR_ELT(matList, j));              // This is the matrix + populations
260                 omxProcessMatrixPopulationList(globalState->matrixList[j], nextLoc);
261                 UNPROTECT(1);
262         }
263
264         /* Processing Constraints */
265         omxProcessConstraints(constraints);
266
267         /* Process Confidence Interval List */
268         omxProcessConfidenceIntervals(intervalList);
269
270         /* Process Checkpoint List */
271         omxProcessCheckpointOptions(checkpointList);
272
273         for(int i = 0; i < numChildren; i++) {
274                 omxDuplicateState(globalState->childList[i], globalState);
275         }
276
277         n = globalState->numFreeParams;
278
279         SEXP minimum, estimate, gradient, hessian;
280         PROTECT(minimum = NEW_NUMERIC(1));
281         PROTECT(estimate = allocVector(REALSXP, n));
282         PROTECT(gradient = allocVector(REALSXP, n));
283         PROTECT(hessian = allocMatrix(REALSXP, n, n));
284
285         if (n>0) { memcpy(REAL(estimate), REAL(startVals), sizeof(double)*n); }
286         
287         omxInvokeNPSOL(REAL(minimum), REAL(estimate), REAL(gradient), REAL(hessian), disableOptimizer);
288
289         SEXP code, status, statusMsg, iterations;
290         SEXP evaluations, ans=NULL, names=NULL, algebras, matrices, expectations, optimizer;
291         SEXP intervals, NAmat, intervalCodes, calculatedHessian, stdErrors;
292
293         int numReturns = 14;
294
295         PROTECT(code = NEW_NUMERIC(1));
296         PROTECT(status = allocVector(VECSXP, 3));
297         PROTECT(iterations = NEW_NUMERIC(1));
298         PROTECT(evaluations = NEW_NUMERIC(2));
299         PROTECT(matrices = NEW_LIST(globalState->numMats));
300         PROTECT(algebras = NEW_LIST(globalState->numAlgs));
301         PROTECT(expectations = NEW_LIST(globalState->numExpects));
302
303         PROTECT(optimizer = allocVector(VECSXP, 2));
304         PROTECT(calculatedHessian = allocMatrix(REALSXP, n, n));
305         PROTECT(stdErrors = allocMatrix(REALSXP, n, 1)); // for optimizer
306         PROTECT(names = allocVector(STRSXP, 2)); // for optimizer
307         PROTECT(intervals = allocMatrix(REALSXP, globalState->numIntervals, 2)); // for optimizer
308         PROTECT(intervalCodes = allocMatrix(INTSXP, globalState->numIntervals, 2)); // for optimizer
309         PROTECT(NAmat = allocMatrix(REALSXP, 1, 1)); // In case of missingness
310         REAL(NAmat)[0] = R_NaReal;
311
312         if(fitfunction == NULL) {
313                 REAL(minimum)[0] = R_NaReal;
314         }
315
316         omxSaveState(globalState, REAL(estimate), REAL(minimum)[0]);
317
318         /* Fill in details from the optimizer */
319         SET_VECTOR_ELT(optimizer, 0, gradient);
320         SET_VECTOR_ELT(optimizer, 1, hessian);
321
322         SET_STRING_ELT(names, 0, mkChar("minimum"));
323         SET_STRING_ELT(names, 1, mkChar("estimate"));
324         namesgets(optimizer, names);
325
326         REAL(code)[0] = globalState->inform;
327         REAL(iterations)[0] = globalState->iter;
328         REAL(evaluations)[0] = globalState->computeCount;
329
330         /* Fill Status code. */
331         SET_VECTOR_ELT(status, 0, code);
332         PROTECT(code = NEW_NUMERIC(1));
333         REAL(code)[0] = globalState->statusCode;
334         SET_VECTOR_ELT(status, 1, code);
335         PROTECT(statusMsg = allocVector(STRSXP, 1));
336         SET_STRING_ELT(statusMsg, 0, mkChar(globalState->statusMsg));
337         SET_VECTOR_ELT(status, 2, statusMsg);
338
339         if(numHessians && globalState->fitMatrix != NULL && globalState->optimumStatus >= 0) {          // No hessians or standard errors if the optimum is invalid
340                 if(globalState->numConstraints == 0) {
341                         if(OMX_DEBUG) { Rprintf("Calculating Hessian for Fit Function.\n");}
342                         int gotHessians = omxEstimateHessian(numHessians, .0001, 4, globalState);
343                         if(gotHessians) {
344                                 if(calculateStdErrors) {
345                                         for(int j = 0; j < numHessians; j++) {          //TODO: Fix Hessian calculation to allow more if requested
346                                                 if(OMX_DEBUG) { Rprintf("Calculating Standard Errors for Fit Function.\n");}
347                                                 omxFitFunction* oo = globalState->fitMatrix->fitFunction;
348                                                 if(oo->getStandardErrorFun != NULL) {
349                                                         oo->getStandardErrorFun(oo);
350                                                 } else {
351                                                         omxCalculateStdErrorFromHessian(2.0, oo);
352                                                 }
353                                         }
354                                 }
355                         } else {
356                                 numHessians = 0;
357                         }
358                 } else {
359                         numHessians = 0;
360                 }
361         } else {
362                 numHessians = 0;
363         }
364
365         /* Likelihood-based Confidence Interval Calculation */
366         if(globalState->numIntervals) {
367                 omxNPSOLConfidenceIntervals(REAL(minimum), REAL(estimate), REAL(gradient), REAL(hessian), ciMaxIterations);
368         }  
369
370         handleFreeVarList(globalState, globalState->optimalValues, n);  // Restore to optima for final compute
371         if(!errOut) omxFinalAlgebraCalculation(globalState, matrices, algebras, expectations); 
372
373         omxPopulateFitFunction(globalState, numReturns, &ans, &names);
374
375         if(numHessians) {
376                 omxPopulateHessians(numHessians, globalState->fitMatrix, 
377                         calculatedHessian, stdErrors, calculateStdErrors, n);
378         }
379
380         if(globalState->numIntervals) { // Populate CIs
381                 omxPopulateConfidenceIntervals(globalState, intervals, intervalCodes);
382         }
383         
384         REAL(evaluations)[1] = globalState->computeCount;
385
386         int nextEl = 0;
387
388         SET_STRING_ELT(names, nextEl++, mkChar("minimum"));
389         SET_STRING_ELT(names, nextEl++, mkChar("estimate"));
390         SET_STRING_ELT(names, nextEl++, mkChar("gradient"));
391         SET_STRING_ELT(names, nextEl++, mkChar("hessianCholesky"));
392         SET_STRING_ELT(names, nextEl++, mkChar("status"));
393         SET_STRING_ELT(names, nextEl++, mkChar("iterations"));
394         SET_STRING_ELT(names, nextEl++, mkChar("evaluations"));
395         SET_STRING_ELT(names, nextEl++, mkChar("matrices"));
396         SET_STRING_ELT(names, nextEl++, mkChar("algebras"));
397         SET_STRING_ELT(names, nextEl++, mkChar("expectations"));
398         SET_STRING_ELT(names, nextEl++, mkChar("confidenceIntervals"));
399         SET_STRING_ELT(names, nextEl++, mkChar("confidenceIntervalCodes"));
400         SET_STRING_ELT(names, nextEl++, mkChar("calculatedHessian"));
401         SET_STRING_ELT(names, nextEl++, mkChar("standardErrors"));
402
403         nextEl = 0;
404
405         SET_VECTOR_ELT(ans, nextEl++, minimum);
406         SET_VECTOR_ELT(ans, nextEl++, estimate);
407         SET_VECTOR_ELT(ans, nextEl++, gradient);
408         SET_VECTOR_ELT(ans, nextEl++, hessian);
409         SET_VECTOR_ELT(ans, nextEl++, status);
410         SET_VECTOR_ELT(ans, nextEl++, iterations);
411         SET_VECTOR_ELT(ans, nextEl++, evaluations);
412         SET_VECTOR_ELT(ans, nextEl++, matrices);
413         SET_VECTOR_ELT(ans, nextEl++, algebras);
414         SET_VECTOR_ELT(ans, nextEl++, expectations);
415         SET_VECTOR_ELT(ans, nextEl++, intervals);
416         SET_VECTOR_ELT(ans, nextEl++, intervalCodes);
417         if(numHessians == 0) {
418                 SET_VECTOR_ELT(ans, nextEl++, NAmat);
419         } else {
420                 SET_VECTOR_ELT(ans, nextEl++, calculatedHessian);
421         }
422         if(!calculateStdErrors) {
423                 SET_VECTOR_ELT(ans, nextEl++, NAmat);
424         } else {
425                 SET_VECTOR_ELT(ans, nextEl++, stdErrors);
426         }
427         namesgets(ans, names);
428
429         if(OMX_VERBOSE) {
430                 Rprintf("Inform Value: %d\n", globalState->optimumStatus);
431                 Rprintf("--------------------------\n");
432         }
433
434         /* Free data memory */
435         omxFreeState(globalState);
436
437         if(OMX_DEBUG) {Rprintf("All vectors freed.\n");}
438
439         omxProtectRestore(initialpix);
440
441         return(ans);
442
443 }
444