Fix build without NPSOL
[openmx:openmx.git] / src / glue.cpp
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 <stdio.h>
18 #include <sys/types.h>
19 #include <errno.h>
20
21 #include <R.h>
22 #include <Rinternals.h>
23 #include <Rdefines.h>
24 #include <R_ext/Rdynload.h>
25 #include <R_ext/BLAS.h>
26 #include <R_ext/Lapack.h>
27
28 #include "omxDefines.h"
29 #include "types.h"
30 #include "glue.h"
31 #include "omxOpenmpWrap.h"
32 #include "omxState.h"
33 #include "omxMatrix.h"
34 #include "omxAlgebra.h"
35 #include "omxFitFunction.h"
36 #include "omxExpectation.h"
37 #include "omxNPSOLSpecific.h"
38 #include "omxImportFrontendState.h"
39 #include "omxExportBackendState.h"
40 #include "Compute.h"
41 #include "dmvnorm.h"
42 #include "npsolswitch.h"
43
44 static SEXP has_NPSOL()
45 { return ScalarLogical(HAS_NPSOL); }
46
47 static R_CallMethodDef callMethods[] = {
48         {"omxBackend", (DL_FUNC) omxBackend, 11},
49         {"omxCallAlgebra", (DL_FUNC) omxCallAlgebra, 3},
50         {"findIdenticalRowsData", (DL_FUNC) findIdenticalRowsData, 5},
51         {"imxDmvnorm_wrapper", (DL_FUNC) dmvnorm_wrapper, 3},
52         {"imxHasNPSOL", (DL_FUNC) has_NPSOL, 0},
53         {NULL, NULL, 0}
54 };
55
56 #ifdef  __cplusplus
57 extern "C" {
58 #endif
59
60 void R_init_OpenMx(DllInfo *info) {
61         R_registerRoutines(info, NULL, callMethods, NULL, NULL);
62
63         // There is no code that will change behavior whether openmp
64         // is set for nested or not. I'm just keeping this in case it
65         // makes a difference with older versions of openmp. 2012-12-24 JNP
66 #if defined(_OPENMP) && _OPENMP <= 200505
67         omp_set_nested(0);
68 #endif
69 }
70
71 void R_unload_OpenMx(DllInfo *) {
72         // keep this stub in case we need it
73 }
74
75 #ifdef  __cplusplus
76 }
77 #endif
78
79 void string_to_try_error( const std::string& str )
80 {
81         error("%s", str.c_str());
82 }
83
84 void exception_to_try_error( const std::exception& ex )
85 {
86         string_to_try_error(ex.what());
87 }
88
89 SEXP asR(MxRList *out)
90 {
91         // detect duplicate keys TODO
92         SEXP names, ans;
93         int len = out->size();
94         PROTECT(names = allocVector(STRSXP, len));
95         PROTECT(ans = allocVector(VECSXP, len));
96         for (int lx=0; lx < len; ++lx) {
97                 SET_STRING_ELT(names, lx, (*out)[lx].first);
98                 SET_VECTOR_ELT(ans,   lx, (*out)[lx].second);
99         }
100         namesgets(ans, names);
101         return ans;
102 }
103
104 /* Main functions */
105 SEXP omxCallAlgebra2(SEXP matList, SEXP algNum, SEXP) {
106
107         omxManageProtectInsanity protectManager;
108
109         if(OMX_DEBUG) { mxLog("-----------------------------------------------------------------------");}
110         if(OMX_DEBUG) { mxLog("Explicit call to algebra %d.", INTEGER(algNum)[0]);}
111
112         int j,k,l;
113         omxMatrix* algebra;
114         int algebraNum = INTEGER(algNum)[0];
115         SEXP ans, nextMat;
116         char output[MAX_STRING_LEN];
117
118         FitContext::setRFitFunction(NULL);
119         Global = new omxGlobal;
120
121         globalState = new omxState;
122         omxInitState(globalState);
123         if(OMX_DEBUG) { mxLog("Created state object at %p.", globalState);}
124
125         /* Retrieve All Matrices From the MatList */
126
127         if(OMX_DEBUG) { mxLog("Processing %d matrix(ces).", length(matList));}
128
129         omxMatrix *args[length(matList)];
130         for(k = 0; k < length(matList); k++) {
131                 PROTECT(nextMat = VECTOR_ELT(matList, k));      // This is the matrix + populations
132                 args[k] = omxNewMatrixFromRPrimitive(nextMat, globalState, 1, - k - 1);
133                 globalState->matrixList.push_back(args[k]);
134                 if(OMX_DEBUG) {
135                         mxLog("Matrix initialized at %p = (%d x %d).",
136                                 globalState->matrixList[k], globalState->matrixList[k]->rows, globalState->matrixList[k]->cols);
137                 }
138         }
139
140         algebra = omxNewAlgebraFromOperatorAndArgs(algebraNum, args, length(matList), globalState);
141
142         if(algebra==NULL) {
143                 error(globalState->statusMsg);
144         }
145
146         if(OMX_DEBUG) {mxLog("Completed Algebras and Matrices.  Beginning Initial Compute.");}
147         omxStateNextEvaluation(globalState);
148
149         omxRecompute(algebra);
150
151         PROTECT(ans = allocMatrix(REALSXP, algebra->rows, algebra->cols));
152         for(l = 0; l < algebra->rows; l++)
153                 for(j = 0; j < algebra->cols; j++)
154                         REAL(ans)[j * algebra->rows + l] =
155                                 omxMatrixElement(algebra, l, j);
156
157         if(OMX_DEBUG) { mxLog("All Algebras complete."); }
158
159         output[0] = 0;
160         if (isErrorRaised(globalState)) {
161                 strncpy(output, globalState->statusMsg, MAX_STRING_LEN);
162         }
163
164         omxFreeAllMatrixData(algebra);
165         omxFreeState(globalState);
166         delete Global;
167
168         if(output[0]) error(output);
169
170         return ans;
171 }
172
173 SEXP omxCallAlgebra(SEXP matList, SEXP algNum, SEXP options)
174 {
175         try {
176                 return omxCallAlgebra2(matList, algNum, options);
177         } catch( std::exception& __ex__ ) {
178                 exception_to_try_error( __ex__ );
179         } catch(...) {
180                 string_to_try_error( "c++ exception (unknown reason)" );
181         }
182 }
183
184 static void
185 friendlyStringToLogical(const char *key, const char *str, int *out)
186 {
187         int understood = FALSE;
188         int newVal;
189         if (matchCaseInsensitive(str, "Yes")) {
190                 understood = TRUE;
191                 newVal = 1;
192         } else if (matchCaseInsensitive(str, "No")) {
193                 understood = TRUE;
194                 newVal = 0;
195         } else if (isdigit(str[0]) && (atoi(str) == 1 || atoi(str) == 0)) {
196                 understood = TRUE;
197                 newVal = atoi(str);
198         }
199         if (!understood) {
200                 warning("Expecting 'Yes' or 'No' for '%s' but got '%s', ignoring", key, str);
201                 return;
202         }
203         if(OMX_DEBUG) { mxLog("%s=%d", key, newVal); }
204         *out = newVal;
205 }
206
207 static void readOpts(SEXP options, int *ciMaxIterations, int *numThreads,
208                      int *analyticGradients)
209 {
210                 char optionCharArray[250] = "";                 // For setting options
211                 int numOptions = length(options);
212                 SEXP optionNames;
213                 PROTECT(optionNames = GET_NAMES(options));
214                 for(int i = 0; i < numOptions; i++) {
215                         const char *nextOptionName = CHAR(STRING_ELT(optionNames, i));
216                         const char *nextOptionValue = STRING_VALUE(VECTOR_ELT(options, i));
217                         if (matchCaseInsensitive(nextOptionName, "CI Max Iterations")) {
218                                 int newvalue = atoi(nextOptionValue);
219                                 if (newvalue > 0) *ciMaxIterations = newvalue;
220                         } else if(matchCaseInsensitive(nextOptionName, "Analytic Gradients")) {
221                                 friendlyStringToLogical(nextOptionName, nextOptionValue, analyticGradients);
222                         } else if(matchCaseInsensitive(nextOptionName, "Number of Threads")) {
223                                 *numThreads = atoi(nextOptionValue);
224                                 if (*numThreads < 1) {
225                                         warning("Computation will be too slow with %d threads; using 1 thread instead", *numThreads);
226                                         *numThreads = 1;
227                                 }
228                         } else {
229                                 // ignore
230                         }
231                 }
232                 UNPROTECT(1); // optionNames
233 }
234
235 SEXP omxBackend2(SEXP computeIndex, SEXP constraints, SEXP matList,
236                  SEXP varList, SEXP algList, SEXP expectList, SEXP computeList,
237                  SEXP data, SEXP intervalList, SEXP checkpointList, SEXP options)
238 {
239         SEXP nextLoc;
240
241         /* Sanity Check and Parse Inputs */
242         /* TODO: Need to find a way to account for nullness in these.  For now, all checking is done on the front-end. */
243 //      if(!isVector(matList)) error ("matList must be a list");
244 //      if(!isVector(algList)) error ("algList must be a list");
245
246         omxManageProtectInsanity protectManager;
247
248         FitContext::setRFitFunction(NULL);
249         Global = new omxGlobal;
250
251         /* Create new omxState for current state storage and initialize it. */
252         globalState = new omxState;
253         omxInitState(globalState);
254         if(OMX_DEBUG) { mxLog("Created state object at %p.", globalState);}
255
256         Global->ciMaxIterations = 5;
257         Global->numThreads = 1;
258         Global->analyticGradients = 0;
259         Global->numChildren = 0;
260         readOpts(options, &Global->ciMaxIterations, &Global->numThreads, 
261                         &Global->analyticGradients);
262 #if HAS_NPSOL
263         omxSetNPSOLOpts(options);
264 #endif
265
266         omxProcessMxDataEntities(data);
267         if (isErrorRaised(globalState)) error(globalState->statusMsg);
268     
269         omxProcessMxMatrixEntities(matList);
270         if (isErrorRaised(globalState)) error(globalState->statusMsg);
271
272         std::vector<double> startingValues;
273         omxProcessFreeVarList(varList, &startingValues);
274         if (isErrorRaised(globalState)) error(globalState->statusMsg);
275
276         omxProcessMxExpectationEntities(expectList);
277         if (isErrorRaised(globalState)) error(globalState->statusMsg);
278
279         omxProcessMxAlgebraEntities(algList);
280         if (isErrorRaised(globalState)) error(globalState->statusMsg);
281
282         omxProcessMxFitFunction(algList);
283         if (isErrorRaised(globalState)) error(globalState->statusMsg);
284
285         omxProcessMxComputeEntities(computeList);
286         if (isErrorRaised(globalState)) error(globalState->statusMsg);
287
288         omxCompleteMxExpectationEntities();
289         if (isErrorRaised(globalState)) error(globalState->statusMsg);
290
291         omxCompleteMxFitFunction(algList);
292         if (isErrorRaised(globalState)) error(globalState->statusMsg);
293
294         // This is the chance to check for matrix
295         // conformability, etc.  Any errors encountered should
296         // be reported using R's error() function, not
297         // omxRaiseErrorf.
298
299         omxInitialMatrixAlgebraCompute();
300         omxResetStatus(globalState);
301
302         for(size_t index = 0; index < globalState->matrixList.size(); index++) {
303                 omxMarkDirty(globalState->matrixList[index]);
304         }
305         for(size_t index = 0; index < globalState->algebraList.size(); index++) {
306                 omxMarkDirty(globalState->algebraList[index]);
307         }
308
309         // maybe require a Compute object? TODO
310         omxCompute *topCompute = NULL;
311         if (!isNull(computeIndex)) {
312                 int ox = INTEGER(computeIndex)[0];
313                 topCompute = Global->computeList[ox];
314         }
315
316         /* Process Matrix and Algebra Population Function */
317         /*
318           Each matrix is a list containing a matrix and the other matrices/algebras that are
319           populated into it at each iteration.  The first element is already processed, above.
320           The rest of the list will be processed here.
321         */
322         for(int j = 0; j < length(matList); j++) {
323                 PROTECT(nextLoc = VECTOR_ELT(matList, j));              // This is the matrix + populations
324                 omxProcessMatrixPopulationList(globalState->matrixList[j], nextLoc);
325         }
326
327         omxProcessConstraints(constraints);
328
329         omxProcessConfidenceIntervals(intervalList);
330
331         omxProcessCheckpointOptions(checkpointList);
332
333         for (size_t vg=0; vg < Global->freeGroup.size(); ++vg) {
334                 Global->freeGroup[vg]->cacheDependencies();
335         }
336
337         FitContext fc(startingValues);
338
339         if (topCompute && !isErrorRaised(globalState)) {
340                 // switch varGroup, if necessary TODO
341                 topCompute->compute(&fc);
342         }
343
344         SEXP evaluations;
345         PROTECT(evaluations = NEW_NUMERIC(2));
346
347         REAL(evaluations)[0] = globalState->computeCount;
348
349         if (topCompute && !isErrorRaised(globalState) && globalState->stale) {
350                 fc.copyParamToModel(globalState);
351         }
352
353         MxRList result;
354
355         omxExportResults(globalState, &result); 
356
357         REAL(evaluations)[1] = globalState->computeCount;
358
359         double optStatus = NA_REAL;
360         if (topCompute && !isErrorRaised(globalState)) {
361                 topCompute->reportResults(&fc, &result);
362                 optStatus = topCompute->getOptimizerStatus();
363         }
364
365         MxRList backwardCompatStatus;
366         backwardCompatStatus.push_back(std::make_pair(mkChar("code"), ScalarReal(optStatus)));
367         backwardCompatStatus.push_back(std::make_pair(mkChar("status"),
368                                                       ScalarInteger(-isErrorRaised(globalState))));
369
370         if (isErrorRaised(globalState)) {
371                 SEXP msg;
372                 PROTECT(msg = allocVector(STRSXP, 1));
373                 SET_STRING_ELT(msg, 0, mkChar(globalState->statusMsg));
374                 result.push_back(std::make_pair(mkChar("error"), msg));
375                 backwardCompatStatus.push_back(std::make_pair(mkChar("statusMsg"), msg));
376         }
377
378         result.push_back(std::make_pair(mkChar("status"), asR(&backwardCompatStatus)));
379         result.push_back(std::make_pair(mkChar("evaluations"), evaluations));
380
381         omxFreeState(globalState);
382         delete Global;
383
384         return asR(&result);
385 }
386
387 SEXP omxBackend(SEXP computeIndex, SEXP constraints, SEXP matList,
388                 SEXP varList, SEXP algList, SEXP expectList, SEXP computeList,
389                 SEXP data, SEXP intervalList, SEXP checkpointList, SEXP options)
390 {
391         try {
392                 return omxBackend2(computeIndex, constraints, matList,
393                                    varList, algList, expectList, computeList,
394                                    data, intervalList, checkpointList, options);
395         } catch( std::exception& __ex__ ) {
396                 exception_to_try_error( __ex__ );
397         } catch(...) {
398                 string_to_try_error( "c++ exception (unknown reason)" );
399         }
400 }
401