Reduce use of protect stack
[openmx:openmx.git] / src / omxImportFrontendState.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 <R.h>
18 #include <Rinternals.h>
19 #include <Rdefines.h>
20
21 #include <sys/stat.h>
22 #include <errno.h>
23
24 #include "omxDefines.h"
25 #include "glue.h"
26 #include "omxState.h"
27 #include "omxNPSOLSpecific.h"
28 #include "Compute.h"
29
30 /* Outside R Functions */
31 static int isDir(const char *path);
32
33 int matchCaseInsensitive(const char *source, const char *target) {
34         return strcasecmp(source, target) == 0;
35 }
36
37 void omxProcessMxDataEntities(SEXP data) {
38         SEXP nextLoc;
39         if(OMX_DEBUG) { mxLog("Processing %d data source(s).", length(data));}
40
41         for(int index = 0; index < length(data); index++) {
42                 PROTECT(nextLoc = VECTOR_ELT(data, index));                     // Retrieve the data object
43                 omxNewDataFromMxData(nextLoc, globalState);
44                 if(OMX_DEBUG) {
45                         mxLog("Data initialized at %p = (%d x %d).",
46                                 globalState->dataList[index], globalState->dataList[index]->rows, globalState->dataList[index]->cols);
47                 }
48         }
49 }
50
51 void omxProcessMxMatrixEntities(SEXP matList) {
52         if(OMX_DEBUG) { mxLog("Processing %d matrix(ces).", length(matList));}
53         SEXP nextLoc, nextMat;
54         globalState->matrixList.clear();
55         SEXP matListNames = getAttrib(matList, R_NamesSymbol);
56
57         for(int index = 0; index < length(matList); index++) {
58                 omxManageProtectInsanity protectManager;
59                 PROTECT(nextLoc = VECTOR_ELT(matList, index));          // This is the matrix + populations
60                 PROTECT(nextMat = VECTOR_ELT(nextLoc, 0));              // The first element of the list is the matrix of values
61                 omxMatrix *mat = omxNewMatrixFromRPrimitive(nextMat, globalState, 1, -index - 1);
62                 globalState->matrixList.push_back(mat);
63                 globalState->matrixList[index]->name = CHAR(STRING_ELT(matListNames, index));
64                 if(OMX_DEBUG) {
65                         mxLog("Matrix initialized at %p = (%d x %d).",
66                                 globalState->matrixList[index], globalState->matrixList[index]->rows, globalState->matrixList[index]->cols);
67                 }
68                 if (isErrorRaised(globalState)) return;
69         }
70 }
71
72 void omxProcessMxAlgebraEntities(SEXP algList) {
73         SEXP nextAlgTuple;
74         SEXP algListNames = getAttrib(algList, R_NamesSymbol);
75
76         if(OMX_DEBUG) { mxLog("Processing %d algebras.", length(algList)); }
77
78         for(int index = 0; index < length(algList); index++) {
79                 globalState->algebraList.push_back(omxInitMatrix(NULL, 0, 0, TRUE, globalState));
80         }
81
82         for(int index = 0; index < length(algList); index++) {
83                 omxManageProtectInsanity protectManager;
84                 PROTECT(nextAlgTuple = VECTOR_ELT(algList, index));             // The next algebra or fit function to process
85                 if(IS_S4_OBJECT(nextAlgTuple)) {
86                         // delay until expectations are ready
87                 } else {                                                                // This is an algebra spec.
88                         SEXP initialValue, formula;
89                         PROTECT(initialValue = VECTOR_ELT(nextAlgTuple, 0));
90                         omxFillMatrixFromRPrimitive(globalState->algebraList[index],
91                                 initialValue, globalState, 1, index);
92                         PROTECT(formula = VECTOR_ELT(nextAlgTuple, 1));
93                         omxFillMatrixFromMxAlgebra(globalState->algebraList[index],
94                                 formula, CHAR(STRING_ELT(algListNames, index)));
95                 }
96                 if (isErrorRaised(globalState)) return;
97         }
98 }
99
100 void omxProcessMxFitFunction(SEXP algList)
101 {
102         SEXP nextAlgTuple;
103         SEXP algListNames = getAttrib(algList, R_NamesSymbol);
104
105         for(int index = 0; index < length(algList); index++) {
106                 PROTECT(nextAlgTuple = VECTOR_ELT(algList, index));             // The next algebra or fit function to process
107                 if(IS_S4_OBJECT(nextAlgTuple)) {
108                         SEXP fitFunctionClass;
109                         PROTECT(fitFunctionClass = STRING_ELT(getAttrib(nextAlgTuple, install("class")), 0));
110                         const char *fitType = CHAR(fitFunctionClass);
111                         omxMatrix *fm = globalState->algebraList[index];
112                         omxFillMatrixFromMxFitFunction(fm, fitType, index);
113                         fm->fitFunction->rObj = nextAlgTuple;
114                         fm->name = CHAR(STRING_ELT(algListNames, index));
115                         UNPROTECT(1);   // fitFunctionClass
116                 }
117                 if (isErrorRaised(globalState)) return;
118                 UNPROTECT(1); //nextAlgTuple
119         }
120 }
121
122 void omxCompleteMxFitFunction(SEXP algList)
123 {
124         SEXP nextAlgTuple;
125
126         for(int index = 0; index < length(algList); index++) {
127                 PROTECT(nextAlgTuple = VECTOR_ELT(algList, index));             // The next algebra or fit function to process
128                 if(IS_S4_OBJECT(nextAlgTuple)) {
129                         omxMatrix *fm = globalState->algebraList[index];
130                         setFreeVarGroup(fm->fitFunction, Global->freeGroup[0]);
131                         omxCompleteFitFunction(fm);
132                 }
133                 UNPROTECT(1);
134         }
135 }
136
137 void omxProcessMxExpectationEntities(SEXP expList) {
138         if(OMX_DEBUG) { mxLog("Initializing %d Model Expectation(s).", length(expList));}
139         SEXP nextExp;
140         SEXP eNames = getAttrib(expList, R_NamesSymbol);
141
142         for(int index = 0; index < length(expList); index++) {
143                 PROTECT(nextExp = VECTOR_ELT(expList, index));
144                 omxExpectation *ex = omxNewIncompleteExpectation(nextExp, index, globalState);
145                 ex->name = CHAR(STRING_ELT(eNames, index));
146                 globalState->expectationList.push_back(ex);
147                 if(OMX_DEBUG) {
148                         mxLog("%s incomplete expectation set up at %p.",
149                                 (globalState->expectationList[index]->expType
150                                         == NULL ? "Untyped" : globalState->expectationList[index]->expType),
151                                          globalState->expectationList[index]);
152                 }
153                 if (isErrorRaised(globalState)) return;
154         }
155 }
156
157
158 void omxCompleteMxExpectationEntities() {
159         if(OMX_DEBUG) { mxLog("Completing %lu Model Expectation(s).", globalState->expectationList.size());}
160         
161         for(size_t index = 0; index < globalState->expectationList.size(); index++) {
162                 omxCompleteExpectation(globalState->expectationList[index]);
163                 if(OMX_DEBUG) {
164                         mxLog("%s expectation completed at %p.",
165                                 (globalState->expectationList[index]->expType
166                                         == NULL ? "Untyped" : globalState->expectationList[index]->expType),
167                                          globalState->expectationList[index]);
168                 }
169                 if (isErrorRaised(globalState)) return;
170         }
171 }
172
173 void omxProcessMxComputeEntities(SEXP computeList)
174 {
175         SEXP rObj, s4class;
176
177         for(int index = 0; index < length(computeList); index++) {
178                 PROTECT(rObj = VECTOR_ELT(computeList, index));
179                 PROTECT(s4class = STRING_ELT(getAttrib(rObj, install("class")), 0));
180                 omxCompute *compute = omxNewCompute(globalState, CHAR(s4class));
181                 compute->initFromFrontend(rObj);
182                 Global->computeList.push_back(compute);
183         }
184 }
185
186 void omxInitialMatrixAlgebraCompute() {
187         size_t numMats = globalState->matrixList.size();
188         int numAlgs = globalState->algebraList.size();
189
190         if(OMX_DEBUG) {mxLog("Completed Algebras and Matrices.  Beginning Initial Compute.");}
191         omxStateNextEvaluation(globalState);
192
193         for(size_t index = 0; index < numMats; index++) {
194                 omxRecompute(globalState->matrixList[index]);
195         }
196
197         for(int index = 0; index < numAlgs; index++) {
198                 omxRecompute(globalState->algebraList[index]);
199         }
200 }
201
202 /*
203 checkpointList is a list().  Each element refers to one checkpointing request.
204 Each interval request is a list of length 5.
205 The first element is an integer that specifies type: 0 = file, 1 = socket, 2=R_connection
206 For a file, the next two are the directory(string)  and file name (string).
207 For a socket, they are server (string) and port number (int).
208 For a connection, the next one is the R_connection SEXP object.
209 After that is an integer <type> specifier.  0 means minutes, 1 means iterations.
210 The last element is an integer count, indicating the number of <type>s per checkpoint.
211 */
212 void omxProcessCheckpointOptions(SEXP checkpointList) {
213         if(OMX_VERBOSE) { mxLog("Processing Checkpoint Requests.");}
214         globalState->numCheckpoints = length(checkpointList);
215         if(OMX_DEBUG) {mxLog("Found %d checkpoints.", globalState->numCheckpoints); }
216         globalState->checkpointList = (omxCheckpoint*) R_alloc(globalState->numCheckpoints, sizeof(omxCheckpoint));
217         SEXP nextLoc;
218
219         for(int index = 0; index < globalState->numCheckpoints; index++) {
220                 omxCheckpoint *oC = &(globalState->checkpointList[index]);
221
222                 /* Initialize Checkpoint object */
223                 oC->file = NULL;
224                 oC->connection = NULL;
225                 oC->time = 0;
226                 oC->numIterations = 0;
227                 oC->lastCheckpoint = 0;
228
229                 const char *pathName, *fileName;
230                 const char __attribute__((unused)) *serverName;
231
232                 PROTECT(nextLoc = VECTOR_ELT(checkpointList, index));
233                 int next = 0;
234                 oC->type = (omxCheckpointType) INTEGER(VECTOR_ELT(nextLoc, next++))[0];
235                 switch(oC->type) {
236                 case OMX_FILE_CHECKPOINT:{
237                         pathName = CHAR(STRING_ELT(VECTOR_ELT(nextLoc, next++), 0));                    // FIXME: Might need PROTECT()ion
238                         fileName = CHAR(STRING_ELT(VECTOR_ELT(nextLoc, next++), 0));
239                         char sep ='/';
240
241                         if(!isDir(pathName)) {
242                                 error("Unable to open directory %s for checkpoint storage.\n", pathName);
243                         }
244
245                         char* fullname = Calloc(strlen(pathName) + strlen(fileName) + 5, char);
246                         sprintf(fullname, "%s%c%s", pathName, sep, fileName);
247                         if(OMX_VERBOSE) { mxLog("Opening File: %s", fullname); }
248                         oC->file = fopen(fullname, "w");
249                         if(!oC->file) {
250                                 error("Unable to open file %s for checkpoint storage: %s.\n", fullname, strerror(errno));
251                         }
252                         Free(fullname);
253                         oC->saveHessian = FALSE;        // TODO: Decide if this should be true.
254                         break;}
255
256                 case OMX_CONNECTION_CHECKPOINT:{        // NYI :::DEBUG:::
257                         oC->connection = VECTOR_ELT(nextLoc, next++);
258                         error("Warning NYI: Socket checkpoints Not Yet Implemented.\n");
259                         oC->saveHessian = FALSE;
260                         break;}
261                 }
262
263                 int isCount = INTEGER(VECTOR_ELT(nextLoc, next++))[0];
264                 if(isCount) {
265                         oC->numIterations = INTEGER(AS_INTEGER(VECTOR_ELT(nextLoc, next++)))[0];
266                 } else {
267                         oC->time = REAL(AS_NUMERIC(VECTOR_ELT(nextLoc, next++)))[0] * 60;       // Constrained to seconds.
268                         if(oC->time < 1) oC->time = 1;                                                                          // Constrained to at least one.
269                 }
270         }
271 }
272
273 void omxProcessFreeVarList(SEXP varList, std::vector<double> *startingValues)
274 {
275         if(OMX_VERBOSE) { mxLog("Processing Free Parameters."); }
276
277         {
278                 FreeVarGroup *fvg = new FreeVarGroup;
279                 fvg->id = FREEVARGROUP_ALL;   // all variables
280                 Global->freeGroup.push_back(fvg);
281
282                 fvg = new FreeVarGroup;
283                 fvg->id = FREEVARGROUP_NONE;  // no variables
284                 Global->freeGroup.push_back(fvg);
285         }
286
287         SEXP nextVar, nextLoc;
288         int numVars = length(varList);
289         startingValues->resize(numVars);
290         for (int fx = 0; fx < numVars; fx++) {
291                 omxManageProtectInsanity mpi;
292
293                 omxFreeVar *fv = new omxFreeVar;
294                 // default group has free all variables
295                 Global->freeGroup[0]->vars.push_back(fv);
296
297                 fv->name = CHAR(STRING_ELT(GET_NAMES(varList), fx));
298                 PROTECT(nextVar = VECTOR_ELT(varList, fx));
299
300                 PROTECT(nextLoc = VECTOR_ELT(nextVar, 0));
301                 fv->lbound = REAL(nextLoc)[0];
302                 if (ISNA(fv->lbound)) fv->lbound = NEG_INF;
303                 if (fv->lbound == 0.0) fv->lbound = 0.0;
304
305                 PROTECT(nextLoc = VECTOR_ELT(nextVar, 1));
306                 fv->ubound = REAL(nextLoc)[0];
307                 if (ISNA(fv->ubound)) fv->ubound = INF;
308                 if (fv->ubound == 0.0) fv->ubound = -0.0;
309
310                 PROTECT(nextLoc = VECTOR_ELT(nextVar, 2));
311                 int groupCount = length(nextLoc);
312                 for (int gx=0; gx < groupCount; ++gx) {
313                         int group = INTEGER(nextLoc)[gx];
314                         if (group == 0) continue;
315                         Global->findOrCreateVarGroup(group)->vars.push_back(fv);
316                 }
317
318                 PROTECT(nextLoc = VECTOR_ELT(nextVar, 3));
319                 int numDeps = LENGTH(nextLoc);
320                 fv->numDeps = numDeps;
321                 fv->deps = (int*) R_alloc(numDeps, sizeof(int));
322                 for (int i = 0; i < numDeps; i++) {
323                         fv->deps[i] = INTEGER(nextLoc)[i];
324                 }
325
326                 int numLocs = length(nextVar) - 5;
327                 if(OMX_DEBUG) { 
328                         mxLog("Free parameter %d bounded (%f, %f): %d locations", fx, 
329                               fv->lbound, fv->ubound, numLocs);
330                 }
331                 for(int locIndex = 0; locIndex < numLocs; locIndex++) {
332                         PROTECT(nextLoc = VECTOR_ELT(nextVar, locIndex+4));
333                         int* theVarList = INTEGER(nextLoc);
334
335                         omxFreeVarLocation loc;
336                         loc.matrix = theVarList[0];
337                         loc.row = theVarList[1];
338                         loc.col = theVarList[2];
339
340                         fv->locations.push_back(loc);
341                 }
342                 PROTECT(nextLoc = VECTOR_ELT(nextVar, length(nextVar)-1));
343                 double sv = REAL(nextLoc)[0];
344                 /*if (sv < fv->lbound) {
345                         warning("Moving starting value of parameter '%s' within bounds %f -> %f",
346                                 fv->name, sv, fv->lbound);
347                         sv = fv->lbound;
348                 } else if (sv > fv->ubound) {
349                         warning("Moving starting value of parameter '%s' within bounds %f -> %f",
350                                 fv->name, sv, fv->ubound);
351                         sv = fv->ubound;
352                 }*/
353                 (*startingValues)[fx] = sv;
354         }
355 }
356
357 /*
358         intervalList is a list().  Each element refers to one confidence interval request.
359         Each interval request is a length 5 vector of REAL.
360         The first three elements are the matrixPointer, Row, and Column of the element
361         for which bounds are to be calculated, and are cast to ints here for speed.
362         The last two are the upper and lower boundaries for the confidence space (respectively).
363 */
364 void omxProcessConfidenceIntervals(SEXP intervalList)  {
365         SEXP nextVar;
366         if(OMX_VERBOSE) { mxLog("Processing Confidence Interval Requests.");}
367         Global->numIntervals = length(intervalList);
368         if(OMX_DEBUG) {mxLog("Found %d requests.", Global->numIntervals); }
369         Global->intervalList = (omxConfidenceInterval*) R_alloc(Global->numIntervals, sizeof(omxConfidenceInterval));
370         for(int index = 0; index < Global->numIntervals; index++) {
371                 omxConfidenceInterval *oCI = &(Global->intervalList[index]);
372                 PROTECT(nextVar = VECTOR_ELT(intervalList, index));
373                 double* intervalInfo = REAL(nextVar);
374                 oCI->matrix = omxMatrixLookupFromState1( nextVar, globalState); // Expects an R object
375                 oCI->row = (int) intervalInfo[1];               // Cast to int in C to save memory/Protection ops
376                 oCI->col = (int) intervalInfo[2];               // Cast to int in C to save memory/Protection ops
377                 oCI->lbound = intervalInfo[3];
378                 oCI->ubound = intervalInfo[4];
379                 oCI->max = R_NaReal;                                    // NAs, in case something goes wrong
380                 oCI->min = R_NaReal;
381         }
382         if(OMX_VERBOSE) { mxLog("Processed."); }
383         if(OMX_DEBUG) { mxLog("%d intervals requested.", Global->numIntervals); }
384 }
385
386 void omxProcessConstraints(SEXP constraints)  {
387         int ncnln = 0; 
388         if(OMX_VERBOSE) { mxLog("Processing Constraints.");}
389         omxMatrix *arg1, *arg2;
390         SEXP nextVar, nextLoc;
391         globalState->numConstraints = length(constraints);
392         if(OMX_DEBUG) {mxLog("Found %d constraints.", globalState->numConstraints); }
393         globalState->conList = (omxConstraint*) R_alloc(globalState->numConstraints, sizeof(omxConstraint));
394         ncnln = 0;
395         for(int constraintIndex = 0; constraintIndex < globalState->numConstraints; constraintIndex++) {
396                 PROTECT(nextVar = VECTOR_ELT(constraints, constraintIndex));
397                 PROTECT(nextLoc = VECTOR_ELT(nextVar, 0));
398                 arg1 = omxMatrixLookupFromState1(nextLoc, globalState);
399                 PROTECT(nextLoc = VECTOR_ELT(nextVar, 1));
400                 arg2 = omxMatrixLookupFromState1(nextLoc, globalState);
401                 PROTECT(nextLoc = AS_INTEGER(VECTOR_ELT(nextVar, 2)));
402                 globalState->conList[constraintIndex].opCode = INTEGER(nextLoc)[0];
403                 omxMatrix *args[2] = {arg1, arg2};
404                 globalState->conList[constraintIndex].result = omxNewAlgebraFromOperatorAndArgs(10, args, 2, globalState); // 10 = binary subtract
405                 omxRecompute(globalState->conList[constraintIndex].result);
406                 int nrows = globalState->conList[constraintIndex].result->rows;
407                 int ncols = globalState->conList[constraintIndex].result->cols;
408                 globalState->conList[constraintIndex].size = nrows * ncols;
409                 ncnln += globalState->conList[constraintIndex].size;
410         }
411         if(OMX_VERBOSE) { mxLog("Processed."); }
412         if(OMX_DEBUG) { mxLog("%d effective constraints.", ncnln); }
413         globalState->ncnln = ncnln;
414 }
415
416 /*
417 *  Acknowledgement:
418 *  This function is duplicated from the function of the same name in the R source code.
419 *  The function appears in src/main/sysutils.c
420 *  Thanks to Michael Spiegel for finding it.
421 *  This code is licensed under the terms of the GNU General Public License.
422 */
423 static int isDir(const char *path)
424 {
425     struct stat sb;
426     int isdir = 0;
427     if(!path) return 0;
428     if(stat(path, &sb) == 0) {
429         isdir = (sb.st_mode & S_IFDIR) > 0; /* is a directory */
430         /* We want to know if the directory is writable by this user,
431            which mode does not tell us */
432         isdir &= (access(path, W_OK) == 0);
433     }
434     return isdir;
435 }
436