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