Add regression test for RFitFunction and fix related bugs
[openmx:openmx.git] / src / omxImportFrontendState.c
1 /*
2  *  Copyright 2007-2012 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, int n) {
293         SEXP nextVar, nextLoc;
294         if(OMX_VERBOSE) { Rprintf("Processing Free Parameters.\n"); }
295         globalState->freeVarList = (omxFreeVar*) R_alloc (n, sizeof (omxFreeVar));                      // Data for replacement of free vars
296         for(int freeVarIndex = 0; freeVarIndex < n; freeVarIndex++) {
297                 int numDeps;
298                 PROTECT(nextVar = VECTOR_ELT(varList, freeVarIndex));
299                 int numLocs = length(nextVar) - 3;
300                 globalState->freeVarList[freeVarIndex].numLocations = numLocs;
301                 globalState->freeVarList[freeVarIndex].matrices = (int*) R_alloc(numLocs, sizeof(int));
302                 globalState->freeVarList[freeVarIndex].row               = (int*) R_alloc(numLocs, sizeof(int));
303                 globalState->freeVarList[freeVarIndex].col               = (int*) R_alloc(numLocs, sizeof(int));
304                 globalState->freeVarList[freeVarIndex].name = CHAR(STRING_ELT(GET_NAMES(varList), freeVarIndex));
305
306                 /* Lower Bound */
307                 PROTECT(nextLoc = VECTOR_ELT(nextVar, 0));                                                      // Position 0 is lower bound.
308                 globalState->freeVarList[freeVarIndex].lbound = REAL(nextLoc)[0];
309                 if(ISNA(globalState->freeVarList[freeVarIndex].lbound)) globalState->freeVarList[freeVarIndex].lbound = NEG_INF;
310                 if(globalState->freeVarList[freeVarIndex].lbound == 0.0) globalState->freeVarList[freeVarIndex].lbound = 0.0;
311                 UNPROTECT(1); // nextLoc
312                 /* Upper Bound */
313                 PROTECT(nextLoc = VECTOR_ELT(nextVar, 1));                                                      // Position 1 is upper bound.
314                 globalState->freeVarList[freeVarIndex].ubound = REAL(nextLoc)[0];
315                 if(ISNA(globalState->freeVarList[freeVarIndex].ubound)) globalState->freeVarList[freeVarIndex].ubound = INF;
316                 if(globalState->freeVarList[freeVarIndex].ubound == 0.0) globalState->freeVarList[freeVarIndex].ubound = -0.0;
317                 UNPROTECT(1); // nextLoc
318
319                 PROTECT(nextLoc = VECTOR_ELT(nextVar, 2));                                                      // Position 2 is a vector of dependencies.
320                 numDeps = LENGTH(nextLoc);
321                 globalState->freeVarList[freeVarIndex].numDeps = numDeps;
322                 globalState->freeVarList[freeVarIndex].deps = (int*) R_alloc(numDeps, sizeof(int));
323                 for(int i = 0; i < numDeps; i++) {
324                         globalState->freeVarList[freeVarIndex].deps[i] = INTEGER(nextLoc)[i];
325                 }
326                 UNPROTECT(1); // nextLoc
327
328
329                 if(OMX_DEBUG) { 
330                         Rprintf("Free parameter %d bounded (%f, %f): %d locations\n", freeVarIndex, 
331                                 globalState->freeVarList[freeVarIndex].lbound, 
332                                 globalState->freeVarList[freeVarIndex].ubound, numLocs);
333                 }
334                 for(int locIndex = 0; locIndex < globalState->freeVarList[freeVarIndex].numLocations; locIndex++) {
335                         PROTECT(nextLoc = VECTOR_ELT(nextVar, locIndex+3));
336                         int* theVarList = INTEGER(nextLoc);                     // These come through as integers.
337
338                         int theMat = theVarList[0];                     // Matrix is zero-based indexed.
339                         int theRow = theVarList[1];                     // Row is zero-based.
340                         int theCol = theVarList[2];                     // Column is zero-based.
341
342                         globalState->freeVarList[freeVarIndex].matrices[locIndex] = theMat;
343                         globalState->freeVarList[freeVarIndex].row[locIndex] = theRow;
344                         globalState->freeVarList[freeVarIndex].col[locIndex] = theCol;
345                         UNPROTECT(1); // nextLoc
346                 }
347                 UNPROTECT(1); // nextVar
348         }
349
350         int numMats = globalState->numMats;
351
352         for(int freeVarIndex = 0; freeVarIndex < n; freeVarIndex++) {
353                 omxFreeVar* freeVar = globalState->freeVarList + freeVarIndex;
354                 int *deps   = freeVar->deps;
355                 int numDeps = freeVar->numDeps;
356                 for (int index = 0; index < numDeps; index++) {
357                         globalState->markMatrices[deps[index] + numMats] = 1;
358                 }
359         }
360
361 }
362
363 /*
364         intervalList is a list().  Each element refers to one confidence interval request.
365         Each interval request is a length 5 vector of REAL.
366         The first three elements are the matrixPointer, Row, and Column of the element
367         for which bounds are to be calculated, and are cast to ints here for speed.
368         The last two are the upper and lower boundaries for the confidence space (respectively).
369 */
370 void omxProcessConfidenceIntervals(SEXP intervalList)  {
371         SEXP nextVar;
372         if(OMX_VERBOSE) { Rprintf("Processing Confidence Interval Requests.\n");}
373         globalState->numIntervals = length(intervalList);
374         if(OMX_DEBUG) {Rprintf("Found %d requests.\n", globalState->numIntervals); }
375         globalState->intervalList = (omxConfidenceInterval*) R_alloc(globalState->numIntervals, sizeof(omxConfidenceInterval));
376         for(int index = 0; index < globalState->numIntervals; index++) {
377                 omxConfidenceInterval *oCI = &(globalState->intervalList[index]);
378                 PROTECT(nextVar = VECTOR_ELT(intervalList, index));
379                 double* intervalInfo = REAL(nextVar);
380                 oCI->matrix = omxNewMatrixFromMxIndex( nextVar, globalState);   // Expects an R object
381                 oCI->row = (int) intervalInfo[1];               // Cast to int in C to save memory/Protection ops
382                 oCI->col = (int) intervalInfo[2];               // Cast to int in C to save memory/Protection ops
383                 oCI->lbound = intervalInfo[3];
384                 oCI->ubound = intervalInfo[4];
385                 UNPROTECT(1);
386                 oCI->max = R_NaReal;                                    // NAs, in case something goes wrong
387                 oCI->min = R_NaReal;
388         }
389         if(OMX_VERBOSE) { Rprintf("Processed.\n"); }
390         if(OMX_DEBUG) { Rprintf("%d intervals requested.\n", globalState->numIntervals); }
391 }
392
393 void omxProcessConstraints(SEXP constraints)  {
394         int ncnln = 0; 
395         if(OMX_VERBOSE) { Rprintf("Processing Constraints.\n");}
396         omxMatrix *arg1, *arg2;
397         SEXP nextVar, nextLoc;
398         globalState->numConstraints = length(constraints);
399         if(OMX_DEBUG) {Rprintf("Found %d constraints.\n", globalState->numConstraints); }
400         globalState->conList = (omxConstraint*) R_alloc(globalState->numConstraints, sizeof(omxConstraint));
401         ncnln = 0;
402         for(int constraintIndex = 0; constraintIndex < globalState->numConstraints; constraintIndex++) {
403                 PROTECT(nextVar = VECTOR_ELT(constraints, constraintIndex));
404                 PROTECT(nextLoc = VECTOR_ELT(nextVar, 0));
405                 arg1 = omxNewMatrixFromMxIndex(nextLoc, globalState);
406                 PROTECT(nextLoc = VECTOR_ELT(nextVar, 1));
407                 arg2 = omxNewMatrixFromMxIndex(nextLoc, globalState);
408                 PROTECT(nextLoc = AS_INTEGER(VECTOR_ELT(nextVar, 2)));
409                 globalState->conList[constraintIndex].opCode = INTEGER(nextLoc)[0];
410                 UNPROTECT(4);
411                 omxMatrix *args[2] = {arg1, arg2};
412                 globalState->conList[constraintIndex].result = omxNewAlgebraFromOperatorAndArgs(10, args, 2, globalState); // 10 = binary subtract
413                 omxRecompute(globalState->conList[constraintIndex].result);
414                 int nrows = globalState->conList[constraintIndex].result->rows;
415                 int ncols = globalState->conList[constraintIndex].result->cols;
416                 globalState->conList[constraintIndex].size = nrows * ncols;
417                 ncnln += globalState->conList[constraintIndex].size;
418         }
419         if(OMX_VERBOSE) { Rprintf("Processed.\n"); }
420         if(OMX_DEBUG) { Rprintf("%d effective constraints.\n", ncnln); }
421         globalState->ncnln = ncnln;
422         globalState->nclin = 0;
423 }
424
425 void omxSetupBoundsAndConstraints(double * bl, double * bu, int n, int nclin) {
426         /* Set min and max limits */
427         for(int index = 0; index < n; index++) {
428                 bl[index] = globalState->freeVarList[index].lbound;                             // -Infinity'd be -10^20.
429                 bu[index] = globalState->freeVarList[index].ubound;                             // Infinity would be at 10^20.
430         }
431
432         for(int index = n; index < n+nclin; index++) {                                          // At present, nclin == 0.
433                 bl[index] = NEG_INF;                                                    // Linear constraints have no bounds.
434                 bu[index] = INF;                                                                // Because there are no linear constraints.
435         }                                                                                                   // But if there were, they would go here.
436
437         int index = n + nclin;
438         for(int constraintIndex = 0; constraintIndex < globalState->numConstraints; constraintIndex++) {                // Nonlinear constraints:
439                 if(OMX_DEBUG) { Rprintf("Constraint %d: ", constraintIndex);}
440                 switch(globalState->conList[constraintIndex].opCode) {
441                 case 0:                                                                 // Less than: Must be strictly less than 0.
442                         if(OMX_DEBUG) { Rprintf("Bounded at (-INF, 0).\n");}
443                         for(int offset = 0; offset < globalState->conList[constraintIndex].size; offset++) {
444                                 bl[index] = NEG_INF;
445                                 bu[index] = -0.0;
446                                 index++;
447                         }
448                         break;
449                 case 1:                                                                 // Equal: Must be roughly equal to 0.
450                         if(OMX_DEBUG) { Rprintf("Bounded at (-0, 0).\n");}
451                         for(int offset = 0; offset < globalState->conList[constraintIndex].size; offset++) {
452                                 bl[index] = -0.0;
453                                 bu[index] = 0.0;
454                                 index++;
455                         }
456                         break;
457                 case 2:                                                                 // Greater than: Must be strictly greater than 0.
458                         if(OMX_DEBUG) { Rprintf("Bounded at (0, INF).\n");}
459                         for(int offset = 0; offset < globalState->conList[constraintIndex].size; offset++) {
460                                 if(OMX_DEBUG) { Rprintf("\tBounds set for constraint %d.%d.\n", constraintIndex, offset);}
461                                 bl[index] = 0.0;
462                                 bu[index] = INF;
463                                 index++;
464                         }
465                         break;
466                 default:
467                         if(OMX_DEBUG) { Rprintf("Bounded at (-INF, INF).\n");}
468                         for(int offset = 0; offset < globalState->conList[constraintIndex].size; offset++) {
469                                 bl[index] = NEG_INF;
470                                 bu[index] = INF;
471                                 index++;
472                         }
473                         break;
474                 }
475         }
476 }
477
478 /*
479 *  Acknowledgement:
480 *  This function is duplicated from the function of the same name in the R source code.
481 *  The function appears in src/main/sysutils.c
482 *  Thanks to Michael Spiegel for finding it.
483 */
484 static int isDir(const char *path)
485 {
486     struct stat sb;
487     int isdir = 0;
488     if(!path) return 0;
489     if(stat(path, &sb) == 0) {
490         isdir = (sb.st_mode & S_IFDIR) > 0; /* is a directory */
491         /* We want to know if the directory is writable by this user,
492            which mode does not tell us */
493         isdir &= (access(path, W_OK) == 0);
494     }
495     return isdir;
496 }
497