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