Initialize freeVarList and numFreeParam together
[openmx:openmx.git] / src / omxImportFrontendState.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
21 #include <sys/stat.h>
22 #include <errno.h>
23
24 #include "omxDefines.h"
25 #include "omxState.h"
26 #include "omxGlobalState.h"
27 #include "omxNPSOLSpecific.h"
28
29 /* Outside R Functions */
30 static int isDir(const char *path);
31
32 int matchCaseInsensitive(const char *source, const char *target) {
33         return strcasecmp(source, target) == 0;
34 }
35
36 int omxProcessMxDataEntities(SEXP data) {
37         int errOut = FALSE;
38         SEXP nextLoc;
39         if(OMX_DEBUG) { Rprintf("Processing %d data source(s).\n", length(data));}
40         globalState->numData = length(data);
41         globalState->dataList = (omxData**) R_alloc(length(data), sizeof(omxData*));
42
43         for(int index = 0; index < length(data); index++) {
44                 PROTECT(nextLoc = VECTOR_ELT(data, index));                     // Retrieve the data object
45                 globalState->dataList[index] = omxNewDataFromMxData(NULL, nextLoc, globalState);
46                 if(OMX_DEBUG) {
47                         Rprintf("Data initialized at 0x%0xd = (%d x %d).\n",
48                                 globalState->dataList[index], globalState->dataList[index]->rows, globalState->dataList[index]->cols);
49                 }
50                 UNPROTECT(1); // nextMat
51                 if(globalState->statusCode < 0) {
52                         errOut = TRUE;
53                         globalState->numData = index+1;
54                         break;
55                 }
56         }
57         return(errOut);
58 }
59
60 int omxProcessMxMatrixEntities(SEXP matList) {
61         if(OMX_DEBUG) { Rprintf("Processing %d matrix(ces).\n", length(matList));}
62         int errOut = FALSE;
63         SEXP nextLoc, nextMat;
64         globalState->numMats = length(matList);
65         globalState->matrixList = (omxMatrix**) R_alloc(length(matList), sizeof(omxMatrix*));
66         SEXP matListNames = getAttrib(matList, R_NamesSymbol);
67
68         for(int index = 0; index < length(matList); index++) {
69                 PROTECT(nextLoc = VECTOR_ELT(matList, index));          // This is the matrix + populations
70                 PROTECT(nextMat = VECTOR_ELT(nextLoc, 0));              // The first element of the list is the matrix of values
71                 globalState->matrixList[index] = omxNewMatrixFromRPrimitive(
72                         nextMat, globalState, 1, -index - 1);
73                 globalState->matrixList[index]->name = CHAR(STRING_ELT(matListNames, index));
74                 if(OMX_DEBUG) {
75                         Rprintf("Matrix initialized at 0x%0xd = (%d x %d).\n",
76                                 globalState->matrixList[index], globalState->matrixList[index]->rows, globalState->matrixList[index]->cols);
77                 }
78                 UNPROTECT(2); // nextLoc and nextMat
79                 if(globalState->statusCode < 0) {
80                         if(OMX_DEBUG) { Rprintf("Initialization Error processing %dth matrix.\n", index+1);}
81                         errOut = TRUE;
82                         globalState->numMats = index+1;
83                         break;
84                 }
85         }
86         return(errOut);
87 }
88
89 void omxProcessMxAlgebraEntities(SEXP algList) {
90         SEXP nextAlgTuple;
91         globalState->numAlgs = length(algList);
92         SEXP algListNames = getAttrib(algList, R_NamesSymbol);
93
94         if(OMX_DEBUG) { Rprintf("Processing %d algebras.\n", globalState->numAlgs); }
95         globalState->algebraList = (omxMatrix**) R_alloc(globalState->numAlgs, sizeof(omxMatrix*));
96
97         for(int index = 0; index < globalState->numAlgs; index++) {
98                 globalState->algebraList[index] = omxInitMatrix(NULL, 0, 0, TRUE, globalState);
99         }
100
101         for(int index = 0; index < globalState->numAlgs; index++) {
102                 PROTECT(nextAlgTuple = VECTOR_ELT(algList, index));             // The next algebra or fit function to process
103                 if(OMX_DEBUG) { Rprintf("Initializing algebra %d at location 0x%0x.\n", index, globalState->algebraList + index); }
104                 if(IS_S4_OBJECT(nextAlgTuple)) {                // This is a fit function object.
105                         omxFillMatrixFromMxFitFunction(globalState->algebraList[index], nextAlgTuple, 1, index);
106                 } else {                                                                // This is an algebra spec.
107                         SEXP initialValue, formula, dependencies;
108                         PROTECT(initialValue = VECTOR_ELT(nextAlgTuple, 0));
109                         omxFillMatrixFromRPrimitive(globalState->algebraList[index],
110                                 initialValue, globalState, 1, index);
111                         UNPROTECT(1);   // initialValue
112                         PROTECT(formula = VECTOR_ELT(nextAlgTuple, 1));
113                         omxFillMatrixFromMxAlgebra(globalState->algebraList[index],
114                                 formula, CHAR(STRING_ELT(algListNames, index)));
115                         UNPROTECT(1);   // formula
116                         PROTECT(dependencies = VECTOR_ELT(nextAlgTuple, 2));
117                         UNPROTECT(1);   // dependencies
118                 }
119                 UNPROTECT(1);   // nextAlgTuple
120                 if (globalState->statusMsg[0]) return;
121         }
122 }
123
124 int omxProcessMxExpectationEntities(SEXP expList) {
125         if(OMX_DEBUG) { Rprintf("Initializing %d Model Expectation(s).\n", length(expList));}
126         int errOut = FALSE;
127         SEXP nextExp;
128         globalState->numExpects = length(expList);
129         globalState->expectationList = (omxExpectation**) R_alloc(length(expList), sizeof(omxExpectation*));
130
131         for(int index = 0; index < length(expList); index++) {
132                 PROTECT(nextExp = VECTOR_ELT(expList, index));
133                 globalState->expectationList[index] = omxNewIncompleteExpectation(nextExp, index, globalState);
134                 if(OMX_DEBUG) {
135                         Rprintf("%s incomplete expectation set up at 0x%0xd.\n",
136                                 (globalState->expectationList[index]->expType
137                                         == NULL ? "Untyped" : globalState->expectationList[index]->expType),
138                                          globalState->expectationList[index]);
139                 }
140                 UNPROTECT(1); // nextExp
141                 if(globalState->statusCode < 0) {
142                         if(OMX_DEBUG) { Rprintf("Initialization Error processing %dth Expectation.\n", index+1);}
143                         errOut = TRUE;
144                         globalState->numMats = index+1;
145                         break;
146                 }
147         }
148         return(errOut);
149 }
150
151
152 int omxCompleteMxExpectationEntities() {
153         if(OMX_DEBUG) { Rprintf("Completing %d Model Expectation(s).\n", globalState->numExpects);}
154         int errOut = FALSE;
155         
156         for(int index = 0; index < globalState->numExpects; index++) {
157                 omxCompleteExpectation(globalState->expectationList[index]);
158                 if(OMX_DEBUG) {
159                         Rprintf("%s expectation completed at 0x%0xd.\n",
160                                 (globalState->expectationList[index]->expType
161                                         == NULL ? "Untyped" : globalState->expectationList[index]->expType),
162                                          globalState->expectationList[index]);
163                 }
164                 if(globalState->statusCode < 0) {
165                         if(OMX_DEBUG) { Rprintf("Initialization Error processing %dth Expectation.\n", index+1);}
166                         errOut = TRUE;
167                         globalState->numMats = index+1;
168                         break;
169                 }
170         }
171         return(errOut);
172 }
173
174
175 void omxInitialMatrixAlgebraCompute() {
176         int numMats = globalState->numMats;
177         int numAlgs = globalState->numAlgs;
178
179         if(OMX_DEBUG) {Rprintf("Completed Algebras and Matrices.  Beginning Initial Compute.\n");}
180         omxStateNextEvaluation(globalState);
181
182         globalState->markMatrices = (int*) R_alloc(numMats + numAlgs, sizeof(int));
183
184         for(int index = 0; index < numMats; index++) {
185                 omxRecompute(globalState->matrixList[index]);
186         }
187
188         for(int index = 0; index < numAlgs; index++) {
189                 omxRecompute(globalState->algebraList[index]);
190         }
191 }
192
193 void omxProcessFitFunction(SEXP fitFunction) {
194         if(!isNull(fitFunction)) {
195                 if(OMX_DEBUG) { Rprintf("Processing fit function.\n"); }
196                 globalState->fitMatrix = omxNewMatrixFromMxIndex(fitFunction, globalState);
197         } else {
198                 globalState->fitMatrix = NULL;
199                 globalState->numFreeParams = 0;
200         }
201 }
202
203 /*
204 checkpointList is a list().  Each element refers to one checkpointing request.
205 Each interval request is a list of length 5.
206 The first element is an integer that specifies type: 0 = file, 1 = socket, 2=R_connection
207 For a file, the next two are the directory(string)  and file name (string).
208 For a socket, they are server (string) and port number (int).
209 For a connection, the next one is the R_connection SEXP object.
210 After that is an integer <type> specifier.  0 means minutes, 1 means iterations.
211 The last element is an integer count, indicating the number of <type>s per checkpoint.
212 */
213 void omxProcessCheckpointOptions(SEXP checkpointList) {
214         if(OMX_VERBOSE) { Rprintf("Processing Checkpoint Requests.\n");}
215         globalState->numCheckpoints = length(checkpointList);
216         if(OMX_DEBUG) {Rprintf("Found %d checkpoints.\n", globalState->numCheckpoints); }
217         globalState->checkpointList = (omxCheckpoint*) R_alloc(globalState->numCheckpoints, sizeof(omxCheckpoint));
218         SEXP nextLoc;
219
220         for(int index = 0; index < globalState->numCheckpoints; index++) {
221                 omxCheckpoint *oC = &(globalState->checkpointList[index]);
222
223                 /* Initialize Checkpoint object */
224                 oC->socket = -1;
225                 oC->file = NULL;
226                 oC->connection = NULL;
227                 oC->time = 0;
228                 oC->numIterations = 0;
229                 oC->lastCheckpoint = 0;
230
231                 const char *pathName, *fileName;
232                 const char __attribute__((unused)) *serverName;
233
234                 PROTECT(nextLoc = VECTOR_ELT(checkpointList, index));
235                 int next = 0;
236                 oC->type = INTEGER(VECTOR_ELT(nextLoc, next++))[0];
237                 switch(oC->type) {
238                 case OMX_FILE_CHECKPOINT:
239                         pathName = CHAR(STRING_ELT(VECTOR_ELT(nextLoc, next++), 0));                    // FIXME: Might need PROTECT()ion
240                         fileName = CHAR(STRING_ELT(VECTOR_ELT(nextLoc, next++), 0));
241                         char sep ='/';
242
243                         if(!isDir(pathName)) {
244                                 error("Unable to open directory %s for checkpoint storage.\n", pathName);
245                         }
246
247                         char* fullname = Calloc(strlen(pathName) + strlen(fileName) + 5, char);
248                         sprintf(fullname, "%s%c%s", pathName, sep, fileName);
249                         if(OMX_VERBOSE) { Rprintf("Opening File: %s\n", fullname); }
250                         oC->file = fopen(fullname, "w");
251                         if(!oC->file) {
252                                 error("Unable to open file %s for checkpoint storage: %s.\n", fullname, strerror(errno));
253                         }
254                         Free(fullname);
255                         oC->saveHessian = FALSE;        // TODO: Decide if this should be true.
256                         break;
257
258                 case OMX_SOCKET_CHECKPOINT:
259                         serverName = CHAR(VECTOR_ELT(nextLoc, next++));
260                         int __attribute__((unused)) portno = INTEGER(AS_INTEGER(VECTOR_ELT(nextLoc, next++)))[0];
261                         Rprintf("Warning NYI: Socket checkpoints Not Yet Implemented.\n");
262                         oC->saveHessian = FALSE;
263                         break;
264
265                 case OMX_CONNECTION_CHECKPOINT: // NYI :::DEBUG:::
266                         oC->connection = VECTOR_ELT(nextLoc, next++);
267                         Rprintf("Warning NYI: Socket checkpoints Not Yet Implemented.\n");
268                         oC->saveHessian = FALSE;
269                         break;
270                 }
271
272                 int isCount = INTEGER(VECTOR_ELT(nextLoc, next++))[0];
273                 if(isCount) {
274                         oC->numIterations = INTEGER(AS_INTEGER(VECTOR_ELT(nextLoc, next++)))[0];
275                 } else {
276                         oC->time = REAL(AS_NUMERIC(VECTOR_ELT(nextLoc, next++)))[0] * 60;       // Constrained to seconds.
277                         if(oC->time < 1) oC->time = 1;                                                                          // Constrained to at least one.
278                 }
279
280                 UNPROTECT(1); /* nextLoc */
281
282         }
283 }
284
285 /*
286 varList is a list().  Each element of this list corresponds to one free parameter.
287 Each free parameter is a list.  The first element of this list is the lower bound.
288 The second element of the list is the upper bound.  The third element of the list
289 is a vector of mxIndices specifying the dependencies of the free parameter. 
290 The remaining elements of the list are 3-tuples.  These 3-tuples are (mxIndex, row, col).
291 */
292 void omxProcessFreeVarList(SEXP varList) {
293         int n = globalState->numFreeParams = length(varList);
294         SEXP nextVar, nextLoc;
295         if(OMX_VERBOSE) { Rprintf("Processing Free Parameters.\n"); }
296         globalState->freeVarList = (omxFreeVar*) R_alloc (n, sizeof (omxFreeVar));                      // Data for replacement of free vars
297         for(int freeVarIndex = 0; freeVarIndex < n; freeVarIndex++) {
298                 int numDeps;
299                 PROTECT(nextVar = VECTOR_ELT(varList, freeVarIndex));
300                 int numLocs = length(nextVar) - 3;
301                 globalState->freeVarList[freeVarIndex].numLocations = numLocs;
302                 globalState->freeVarList[freeVarIndex].matrices = (int*) R_alloc(numLocs, sizeof(int));
303                 globalState->freeVarList[freeVarIndex].row               = (int*) R_alloc(numLocs, sizeof(int));
304                 globalState->freeVarList[freeVarIndex].col               = (int*) R_alloc(numLocs, sizeof(int));
305                 globalState->freeVarList[freeVarIndex].name = CHAR(STRING_ELT(GET_NAMES(varList), freeVarIndex));
306
307                 /* Lower Bound */
308                 PROTECT(nextLoc = VECTOR_ELT(nextVar, 0));                                                      // Position 0 is lower bound.
309                 globalState->freeVarList[freeVarIndex].lbound = REAL(nextLoc)[0];
310                 if(ISNA(globalState->freeVarList[freeVarIndex].lbound)) globalState->freeVarList[freeVarIndex].lbound = NEG_INF;
311                 if(globalState->freeVarList[freeVarIndex].lbound == 0.0) globalState->freeVarList[freeVarIndex].lbound = 0.0;
312                 UNPROTECT(1); // nextLoc
313                 /* Upper Bound */
314                 PROTECT(nextLoc = VECTOR_ELT(nextVar, 1));                                                      // Position 1 is upper bound.
315                 globalState->freeVarList[freeVarIndex].ubound = REAL(nextLoc)[0];
316                 if(ISNA(globalState->freeVarList[freeVarIndex].ubound)) globalState->freeVarList[freeVarIndex].ubound = INF;
317                 if(globalState->freeVarList[freeVarIndex].ubound == 0.0) globalState->freeVarList[freeVarIndex].ubound = -0.0;
318                 UNPROTECT(1); // nextLoc
319
320                 PROTECT(nextLoc = VECTOR_ELT(nextVar, 2));                                                      // Position 2 is a vector of dependencies.
321                 numDeps = LENGTH(nextLoc);
322                 globalState->freeVarList[freeVarIndex].numDeps = numDeps;
323                 globalState->freeVarList[freeVarIndex].deps = (int*) R_alloc(numDeps, sizeof(int));
324                 for(int i = 0; i < numDeps; i++) {
325                         globalState->freeVarList[freeVarIndex].deps[i] = INTEGER(nextLoc)[i];
326                 }
327                 UNPROTECT(1); // nextLoc
328
329
330                 if(OMX_DEBUG) { 
331                         Rprintf("Free parameter %d bounded (%f, %f): %d locations\n", freeVarIndex, 
332                                 globalState->freeVarList[freeVarIndex].lbound, 
333                                 globalState->freeVarList[freeVarIndex].ubound, numLocs);
334                 }
335                 for(int locIndex = 0; locIndex < globalState->freeVarList[freeVarIndex].numLocations; locIndex++) {
336                         PROTECT(nextLoc = VECTOR_ELT(nextVar, locIndex+3));
337                         int* theVarList = INTEGER(nextLoc);                     // These come through as integers.
338
339                         int theMat = theVarList[0];                     // Matrix is zero-based indexed.
340                         int theRow = theVarList[1];                     // Row is zero-based.
341                         int theCol = theVarList[2];                     // Column is zero-based.
342
343                         globalState->freeVarList[freeVarIndex].matrices[locIndex] = theMat;
344                         globalState->freeVarList[freeVarIndex].row[locIndex] = theRow;
345                         globalState->freeVarList[freeVarIndex].col[locIndex] = theCol;
346                         UNPROTECT(1); // nextLoc
347                 }
348                 UNPROTECT(1); // nextVar
349         }
350
351         int numMats = globalState->numMats;
352
353         for(int freeVarIndex = 0; freeVarIndex < n; freeVarIndex++) {
354                 omxFreeVar* freeVar = globalState->freeVarList + freeVarIndex;
355                 int *deps   = freeVar->deps;
356                 int numDeps = freeVar->numDeps;
357                 for (int index = 0; index < numDeps; index++) {
358                         globalState->markMatrices[deps[index] + numMats] = 1;
359                 }
360         }
361
362 }
363
364 /*
365         intervalList is a list().  Each element refers to one confidence interval request.
366         Each interval request is a length 5 vector of REAL.
367         The first three elements are the matrixPointer, Row, and Column of the element
368         for which bounds are to be calculated, and are cast to ints here for speed.
369         The last two are the upper and lower boundaries for the confidence space (respectively).
370 */
371 void omxProcessConfidenceIntervals(SEXP intervalList)  {
372         SEXP nextVar;
373         if(OMX_VERBOSE) { Rprintf("Processing Confidence Interval Requests.\n");}
374         globalState->numIntervals = length(intervalList);
375         if(OMX_DEBUG) {Rprintf("Found %d requests.\n", globalState->numIntervals); }
376         globalState->intervalList = (omxConfidenceInterval*) R_alloc(globalState->numIntervals, sizeof(omxConfidenceInterval));
377         for(int index = 0; index < globalState->numIntervals; index++) {
378                 omxConfidenceInterval *oCI = &(globalState->intervalList[index]);
379                 PROTECT(nextVar = VECTOR_ELT(intervalList, index));
380                 double* intervalInfo = REAL(nextVar);
381                 oCI->matrix = omxNewMatrixFromMxIndex( nextVar, globalState);   // Expects an R object
382                 oCI->row = (int) intervalInfo[1];               // Cast to int in C to save memory/Protection ops
383                 oCI->col = (int) intervalInfo[2];               // Cast to int in C to save memory/Protection ops
384                 oCI->lbound = intervalInfo[3];
385                 oCI->ubound = intervalInfo[4];
386                 UNPROTECT(1);
387                 oCI->max = R_NaReal;                                    // NAs, in case something goes wrong
388                 oCI->min = R_NaReal;
389         }
390         if(OMX_VERBOSE) { Rprintf("Processed.\n"); }
391         if(OMX_DEBUG) { Rprintf("%d intervals requested.\n", globalState->numIntervals); }
392 }
393
394 void omxProcessConstraints(SEXP constraints)  {
395         int ncnln = 0; 
396         if(OMX_VERBOSE) { Rprintf("Processing Constraints.\n");}
397         omxMatrix *arg1, *arg2;
398         SEXP nextVar, nextLoc;
399         globalState->numConstraints = length(constraints);
400         if(OMX_DEBUG) {Rprintf("Found %d constraints.\n", globalState->numConstraints); }
401         globalState->conList = (omxConstraint*) R_alloc(globalState->numConstraints, sizeof(omxConstraint));
402         ncnln = 0;
403         for(int constraintIndex = 0; constraintIndex < globalState->numConstraints; constraintIndex++) {
404                 PROTECT(nextVar = VECTOR_ELT(constraints, constraintIndex));
405                 PROTECT(nextLoc = VECTOR_ELT(nextVar, 0));
406                 arg1 = omxNewMatrixFromMxIndex(nextLoc, globalState);
407                 PROTECT(nextLoc = VECTOR_ELT(nextVar, 1));
408                 arg2 = omxNewMatrixFromMxIndex(nextLoc, globalState);
409                 PROTECT(nextLoc = AS_INTEGER(VECTOR_ELT(nextVar, 2)));
410                 globalState->conList[constraintIndex].opCode = INTEGER(nextLoc)[0];
411                 UNPROTECT(4);
412                 omxMatrix *args[2] = {arg1, arg2};
413                 globalState->conList[constraintIndex].result = omxNewAlgebraFromOperatorAndArgs(10, args, 2, globalState); // 10 = binary subtract
414                 omxRecompute(globalState->conList[constraintIndex].result);
415                 int nrows = globalState->conList[constraintIndex].result->rows;
416                 int ncols = globalState->conList[constraintIndex].result->cols;
417                 globalState->conList[constraintIndex].size = nrows * ncols;
418                 ncnln += globalState->conList[constraintIndex].size;
419         }
420         if(OMX_VERBOSE) { Rprintf("Processed.\n"); }
421         if(OMX_DEBUG) { Rprintf("%d effective constraints.\n", ncnln); }
422         globalState->ncnln = ncnln;
423         globalState->nclin = 0;
424 }
425
426 void omxSetupBoundsAndConstraints(double * bl, double * bu, int n, int nclin) {
427         /* Set min and max limits */
428         for(int index = 0; index < n; index++) {
429                 bl[index] = globalState->freeVarList[index].lbound;                             // -Infinity'd be -10^20.
430                 bu[index] = globalState->freeVarList[index].ubound;                             // Infinity would be at 10^20.
431         }
432
433         for(int index = n; index < n+nclin; index++) {                                          // At present, nclin == 0.
434                 bl[index] = NEG_INF;                                                    // Linear constraints have no bounds.
435                 bu[index] = INF;                                                                // Because there are no linear constraints.
436         }                                                                                                   // But if there were, they would go here.
437
438         int index = n + nclin;
439         for(int constraintIndex = 0; constraintIndex < globalState->numConstraints; constraintIndex++) {                // Nonlinear constraints:
440                 if(OMX_DEBUG) { Rprintf("Constraint %d: ", constraintIndex);}
441                 switch(globalState->conList[constraintIndex].opCode) {
442                 case 0:                                                                 // Less than: Must be strictly less than 0.
443                         if(OMX_DEBUG) { Rprintf("Bounded at (-INF, 0).\n");}
444                         for(int offset = 0; offset < globalState->conList[constraintIndex].size; offset++) {
445                                 bl[index] = NEG_INF;
446                                 bu[index] = -0.0;
447                                 index++;
448                         }
449                         break;
450                 case 1:                                                                 // Equal: Must be roughly equal to 0.
451                         if(OMX_DEBUG) { Rprintf("Bounded at (-0, 0).\n");}
452                         for(int offset = 0; offset < globalState->conList[constraintIndex].size; offset++) {
453                                 bl[index] = -0.0;
454                                 bu[index] = 0.0;
455                                 index++;
456                         }
457                         break;
458                 case 2:                                                                 // Greater than: Must be strictly greater than 0.
459                         if(OMX_DEBUG) { Rprintf("Bounded at (0, INF).\n");}
460                         for(int offset = 0; offset < globalState->conList[constraintIndex].size; offset++) {
461                                 if(OMX_DEBUG) { Rprintf("\tBounds set for constraint %d.%d.\n", constraintIndex, offset);}
462                                 bl[index] = 0.0;
463                                 bu[index] = INF;
464                                 index++;
465                         }
466                         break;
467                 default:
468                         if(OMX_DEBUG) { Rprintf("Bounded at (-INF, INF).\n");}
469                         for(int offset = 0; offset < globalState->conList[constraintIndex].size; offset++) {
470                                 bl[index] = NEG_INF;
471                                 bu[index] = INF;
472                                 index++;
473                         }
474                         break;
475                 }
476         }
477 }
478
479 /*
480 *  Acknowledgement:
481 *  This function is duplicated from the function of the same name in the R source code.
482 *  The function appears in src/main/sysutils.c
483 *  Thanks to Michael Spiegel for finding it.
484 */
485 static int isDir(const char *path)
486 {
487     struct stat sb;
488     int isdir = 0;
489     if(!path) return 0;
490     if(stat(path, &sb) == 0) {
491         isdir = (sb.st_mode & S_IFDIR) > 0; /* is a directory */
492         /* We want to know if the directory is writable by this user,
493            which mode does not tell us */
494         isdir &= (access(path, W_OK) == 0);
495     }
496     return isdir;
497 }
498