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