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