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