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