Updating Copyright line (happy new year!)
[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
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].matrices = (int*) R_alloc(numLocs, sizeof(int));
252                 currentState->freeVarList[freeVarIndex].row              = (int*) R_alloc(numLocs, sizeof(int));
253                 currentState->freeVarList[freeVarIndex].col              = (int*) R_alloc(numLocs, sizeof(int));
254                 currentState->freeVarList[freeVarIndex].name = CHAR(STRING_ELT(GET_NAMES(varList), freeVarIndex));
255
256                 /* Lower Bound */
257                 PROTECT(nextLoc = VECTOR_ELT(nextVar, 0));                                                      // Position 0 is lower bound.
258                 currentState->freeVarList[freeVarIndex].lbound = REAL(nextLoc)[0];
259                 if(ISNA(currentState->freeVarList[freeVarIndex].lbound)) currentState->freeVarList[freeVarIndex].lbound = NEG_INF;
260                 if(currentState->freeVarList[freeVarIndex].lbound == 0.0) currentState->freeVarList[freeVarIndex].lbound = 0.0;
261                 UNPROTECT(1); // NextLoc
262                 /* Upper Bound */
263                 PROTECT(nextLoc = VECTOR_ELT(nextVar, 1));                                                      // Position 1 is upper bound.
264                 currentState->freeVarList[freeVarIndex].ubound = REAL(nextLoc)[0];
265                 if(ISNA(currentState->freeVarList[freeVarIndex].ubound)) currentState->freeVarList[freeVarIndex].ubound = INF;
266                 if(currentState->freeVarList[freeVarIndex].ubound == 0.0) currentState->freeVarList[freeVarIndex].ubound = -0.0;
267                 UNPROTECT(1); // NextLoc
268
269                 if(OMX_DEBUG) { 
270                         Rprintf("Free parameter %d bounded (%f, %f): %d locations\n", freeVarIndex, 
271                                 currentState->freeVarList[freeVarIndex].lbound, 
272                                 currentState->freeVarList[freeVarIndex].ubound, numLocs);
273                 }
274                 for(int locIndex = 0; locIndex < currentState->freeVarList[freeVarIndex].numLocations; locIndex++) {
275                         PROTECT(nextLoc = VECTOR_ELT(nextVar, locIndex+2));
276                         int* theVarList = INTEGER(nextLoc);                     // These come through as integers.
277
278                         int theMat = theVarList[0];                     // Matrix is zero-based indexed.
279                         int theRow = theVarList[1];                     // Row is zero-based.
280                         int theCol = theVarList[2];                     // Column is zero-based.
281
282                         currentState->freeVarList[freeVarIndex].matrices[locIndex] = theMat;
283                         currentState->freeVarList[freeVarIndex].row[locIndex] = theRow;
284                         currentState->freeVarList[freeVarIndex].col[locIndex] = theCol;
285                         UNPROTECT(1); // nextLoc
286                 }
287                 UNPROTECT(1); // nextVar
288         }
289 }
290
291 /*
292         intervalList is a list().  Each element refers to one confidence interval request.
293         Each interval request is a length 5 vector of REAL.
294         The first three elements are the matrixPointer, Row, and Column of the element
295         for which bounds are to be calculated, and are cast to ints here for speed.
296         The last two are the upper and lower boundaries for the confidence space (respectively).
297 */
298 void omxProcessConfidenceIntervals(SEXP intervalList)  {
299         SEXP nextVar;
300         if(OMX_VERBOSE) { Rprintf("Processing Confidence Interval Requests.\n");}
301         currentState->numIntervals = length(intervalList);
302         if(OMX_DEBUG) {Rprintf("Found %d requests.\n", currentState->numIntervals); }
303         currentState->intervalList = (omxConfidenceInterval*) R_alloc(currentState->numIntervals, sizeof(omxConfidenceInterval));
304         for(int index = 0; index < currentState->numIntervals; index++) {
305                 omxConfidenceInterval *oCI = &(currentState->intervalList[index]);
306                 PROTECT(nextVar = VECTOR_ELT(intervalList, index));
307                 double* intervalInfo = REAL(nextVar);
308                 oCI->matrix = omxNewMatrixFromMxIndex( nextVar, currentState);  // Expects an R object
309                 oCI->row = (int) intervalInfo[1];               // Cast to int in C to save memory/Protection ops
310                 oCI->col = (int) intervalInfo[2];               // Cast to int in C to save memory/Protection ops
311                 oCI->lbound = intervalInfo[3];
312                 oCI->ubound = intervalInfo[4];
313                 UNPROTECT(1);
314                 oCI->max = R_NaReal;                                    // NAs, in case something goes wrong
315                 oCI->min = R_NaReal;
316         }
317         if(OMX_VERBOSE) { Rprintf("Processed.\n"); }
318         if(OMX_DEBUG) { Rprintf("%d intervals requested.\n", currentState->numIntervals); }
319 }
320
321 int omxProcessConstraints(SEXP constraints)  {
322         int ncnln = 0; 
323         if(OMX_VERBOSE) { Rprintf("Processing Constraints.\n");}
324         omxMatrix *arg1, *arg2;
325         SEXP nextVar, nextLoc;
326         currentState->numConstraints = length(constraints);
327         if(OMX_DEBUG) {Rprintf("Found %d constraints.\n", currentState->numConstraints); }
328         currentState->conList = (omxConstraint*) R_alloc(currentState->numConstraints, sizeof(omxConstraint));
329         ncnln = 0;
330         for(int constraintIndex = 0; constraintIndex < currentState->numConstraints; constraintIndex++) {
331                 PROTECT(nextVar = VECTOR_ELT(constraints, constraintIndex));
332                 PROTECT(nextLoc = VECTOR_ELT(nextVar, 0));
333                 arg1 = omxNewMatrixFromMxIndex(nextLoc, currentState);
334                 PROTECT(nextLoc = VECTOR_ELT(nextVar, 1));
335                 arg2 = omxNewMatrixFromMxIndex(nextLoc, currentState);
336                 PROTECT(nextLoc = AS_INTEGER(VECTOR_ELT(nextVar, 2)));
337                 currentState->conList[constraintIndex].opCode = INTEGER(nextLoc)[0];
338                 UNPROTECT(4);
339                 omxMatrix *args[2] = {arg1, arg2};
340                 currentState->conList[constraintIndex].result = omxNewAlgebraFromOperatorAndArgs(10, args, 2, currentState); // 10 = binary subtract
341                 omxRecompute(currentState->conList[constraintIndex].result);
342                 int nrows = currentState->conList[constraintIndex].result->rows;
343                 int ncols = currentState->conList[constraintIndex].result->cols;
344                 currentState->conList[constraintIndex].size = nrows * ncols;
345                 ncnln += currentState->conList[constraintIndex].size;
346         }
347         if(OMX_VERBOSE) { Rprintf("Processed.\n"); }
348         if(OMX_DEBUG) { Rprintf("%d effective constraints.\n", ncnln); }
349         return(ncnln);
350 }
351
352 void omxSetupBoundsAndConstraints(double * bl, double * bu, int n, int nclin) {
353         /* Set min and max limits */
354         for(int index = 0; index < n; index++) {
355                 bl[index] = currentState->freeVarList[index].lbound;                            // -Infinity'd be -10^20.
356                 bu[index] = currentState->freeVarList[index].ubound;                            // Infinity would be at 10^20.
357         }
358
359         for(int index = n; index < n+nclin; index++) {                                          // At present, nclin == 0.
360                 bl[index] = NEG_INF;                                                    // Linear constraints have no bounds.
361                 bu[index] = INF;                                                                // Because there are no linear constraints.
362         }                                                                                                   // But if there were, they would go here.
363
364         int index = n + nclin;
365         for(int constraintIndex = 0; constraintIndex < currentState->numConstraints; constraintIndex++) {               // Nonlinear constraints:
366                 if(OMX_DEBUG) { Rprintf("Constraint %d: ", constraintIndex);}
367                 switch(currentState->conList[constraintIndex].opCode) {
368                 case 0:                                                                 // Less than: Must be strictly less than 0.
369                         if(OMX_DEBUG) { Rprintf("Bounded at (-INF, 0).\n");}
370                         for(int offset = 0; offset < currentState->conList[constraintIndex].size; offset++) {
371                                 bl[index] = NEG_INF;
372                                 bu[index] = -0.0;
373                                 index++;
374                         }
375                         break;
376                 case 1:                                                                 // Equal: Must be roughly equal to 0.
377                         if(OMX_DEBUG) { Rprintf("Bounded at (-0, 0).\n");}
378                         for(int offset = 0; offset < currentState->conList[constraintIndex].size; offset++) {
379                                 bl[index] = -0.0;
380                                 bu[index] = 0.0;
381                                 index++;
382                         }
383                         break;
384                 case 2:                                                                 // Greater than: Must be strictly greater than 0.
385                         if(OMX_DEBUG) { Rprintf("Bounded at (0, INF).\n");}
386                         for(int offset = 0; offset < currentState->conList[constraintIndex].size; offset++) {
387                                 if(OMX_DEBUG) { Rprintf("\tBounds set for constraint %d.%d.\n", constraintIndex, offset);}
388                                 bl[index] = 0.0;
389                                 bu[index] = INF;
390                                 index++;
391                         }
392                         break;
393                 default:
394                         if(OMX_DEBUG) { Rprintf("Bounded at (-INF, INF).\n");}
395                         for(int offset = 0; offset < currentState->conList[constraintIndex].size; offset++) {
396                                 bl[index] = NEG_INF;
397                                 bu[index] = INF;
398                                 index++;
399                         }
400                         break;
401                 }
402         }
403 }
404
405 /*
406 *  Acknowledgement:
407 *  This function is duplicated from the function of the same name in the R source code.
408 *  The function appears in src/main/sysutils.c
409 *  Thanks to Michael Spiegel for finding it.
410 */
411 static int isDir(const char *path)
412 {
413     struct stat sb;
414     int isdir = 0;
415     if(!path) return 0;
416     if(stat(path, &sb) == 0) {
417         isdir = (sb.st_mode & S_IFDIR) > 0; /* is a directory */
418         /* We want to know if the directory is writable by this user,
419            which mode does not tell us */
420         isdir &= (access(path, W_OK) == 0);
421     }
422     return isdir;
423 }
424