ifa: Fixup
[openmx:openmx.git] / src / omxImportFrontendState.cpp
1 /*
2  *  Copyright 2007-2013 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 #include "Compute.h"
28
29 /* Outside R Functions */
30 static int isDir(const char *path);
31
32 int matchCaseInsensitive(const char *source, const char *target) {
33         return strcasecmp(source, target) == 0;
34 }
35
36 void omxProcessMxDataEntities(SEXP data) {
37         SEXP nextLoc;
38         if(OMX_DEBUG) { mxLog("Processing %d data source(s).", length(data));}
39
40         for(int index = 0; index < length(data); index++) {
41                 PROTECT(nextLoc = VECTOR_ELT(data, index));                     // Retrieve the data object
42                 omxNewDataFromMxData(nextLoc, globalState);
43                 if(OMX_DEBUG) {
44                         mxLog("Data initialized at 0x%0xd = (%d x %d).",
45                                 globalState->dataList[index], globalState->dataList[index]->rows, globalState->dataList[index]->cols);
46                 }
47         }
48 }
49
50 void omxProcessMxMatrixEntities(SEXP matList) {
51         if(OMX_DEBUG) { mxLog("Processing %d matrix(ces).", length(matList));}
52         SEXP nextLoc, nextMat;
53         globalState->matrixList.clear();
54         SEXP matListNames = getAttrib(matList, R_NamesSymbol);
55
56         for(int index = 0; index < length(matList); index++) {
57                 PROTECT(nextLoc = VECTOR_ELT(matList, index));          // This is the matrix + populations
58                 PROTECT(nextMat = VECTOR_ELT(nextLoc, 0));              // The first element of the list is the matrix of values
59                 omxMatrix *mat = omxNewMatrixFromRPrimitive(nextMat, globalState, 1, -index - 1);
60                 globalState->matrixList.push_back(mat);
61                 globalState->matrixList[index]->name = CHAR(STRING_ELT(matListNames, index));
62                 if(OMX_DEBUG) {
63                         mxLog("Matrix initialized at 0x%0xd = (%d x %d).",
64                                 globalState->matrixList[index], globalState->matrixList[index]->rows, globalState->matrixList[index]->cols);
65                 }
66                 if (isErrorRaised(globalState)) return;
67         }
68 }
69
70 void omxProcessMxAlgebraEntities(SEXP algList) {
71         SEXP nextAlgTuple;
72         SEXP algListNames = getAttrib(algList, R_NamesSymbol);
73
74         if(OMX_DEBUG) { mxLog("Processing %d algebras.", length(algList)); }
75
76         for(int index = 0; index < length(algList); index++) {
77                 globalState->algebraList.push_back(omxInitMatrix(NULL, 0, 0, TRUE, globalState));
78         }
79
80         for(int index = 0; index < length(algList); index++) {
81                 PROTECT(nextAlgTuple = VECTOR_ELT(algList, index));             // The next algebra or fit function to process
82                 if(IS_S4_OBJECT(nextAlgTuple)) {
83                         // delay until expectations are ready
84                 } else {                                                                // This is an algebra spec.
85                         SEXP initialValue, formula, dependencies;
86                         PROTECT(initialValue = VECTOR_ELT(nextAlgTuple, 0));
87                         omxFillMatrixFromRPrimitive(globalState->algebraList[index],
88                                 initialValue, globalState, 1, index);
89                         PROTECT(formula = VECTOR_ELT(nextAlgTuple, 1));
90                         omxFillMatrixFromMxAlgebra(globalState->algebraList[index],
91                                 formula, CHAR(STRING_ELT(algListNames, index)));
92                         PROTECT(dependencies = VECTOR_ELT(nextAlgTuple, 2));
93                 }
94                 if (isErrorRaised(globalState)) return;
95         }
96 }
97
98 void omxProcessMxFitFunction(SEXP algList)
99 {
100         SEXP nextAlgTuple;
101
102         for(int index = 0; index < length(algList); index++) {
103                 PROTECT(nextAlgTuple = VECTOR_ELT(algList, index));             // The next algebra or fit function to process
104                 if(IS_S4_OBJECT(nextAlgTuple)) {
105                         SEXP fitFunctionClass;
106                         PROTECT(fitFunctionClass = STRING_ELT(getAttrib(nextAlgTuple, install("class")), 0));
107                         const char *fitType = CHAR(fitFunctionClass);
108                         omxFillMatrixFromMxFitFunction(globalState->algebraList[index], fitType, index);
109                         UNPROTECT(1);   // fitFunctionClass
110                 }
111                 if (isErrorRaised(globalState)) return;
112                 UNPROTECT(1); //nextAlgTuple
113         }
114 }
115
116 void omxCompleteMxFitFunction(SEXP algList)
117 {
118         SEXP nextAlgTuple;
119
120         for(int index = 0; index < length(algList); index++) {
121                 PROTECT(nextAlgTuple = VECTOR_ELT(algList, index));             // The next algebra or fit function to process
122                 if(IS_S4_OBJECT(nextAlgTuple)) {
123                         omxCompleteFitFunction(globalState->algebraList[index], nextAlgTuple);
124                 }
125                 UNPROTECT(1);
126         }
127 }
128
129 void omxProcessMxExpectationEntities(SEXP expList) {
130         if(OMX_DEBUG) { mxLog("Initializing %d Model Expectation(s).", length(expList));}
131         SEXP nextExp;
132
133         for(int index = 0; index < length(expList); index++) {
134                 PROTECT(nextExp = VECTOR_ELT(expList, index));
135                 globalState->expectationList.push_back(omxNewIncompleteExpectation(nextExp, index, globalState));
136                 if(OMX_DEBUG) {
137                         mxLog("%s incomplete expectation set up at 0x%0xd.",
138                                 (globalState->expectationList[index]->expType
139                                         == NULL ? "Untyped" : globalState->expectationList[index]->expType),
140                                          globalState->expectationList[index]);
141                 }
142                 if (isErrorRaised(globalState)) return;
143         }
144 }
145
146
147 void omxCompleteMxExpectationEntities() {
148         if(OMX_DEBUG) { mxLog("Completing %d Model Expectation(s).", globalState->expectationList.size());}
149         
150         for(size_t index = 0; index < globalState->expectationList.size(); index++) {
151                 omxCompleteExpectation(globalState->expectationList[index]);
152                 if(OMX_DEBUG) {
153                         mxLog("%s expectation completed at 0x%0xd.",
154                                 (globalState->expectationList[index]->expType
155                                         == NULL ? "Untyped" : globalState->expectationList[index]->expType),
156                                          globalState->expectationList[index]);
157                 }
158                 if (isErrorRaised(globalState)) return;
159         }
160 }
161
162 void omxProcessMxComputeEntities(SEXP computeList)
163 {
164         SEXP rObj, s4class;
165
166         for(int index = 0; index < length(computeList); index++) {
167                 PROTECT(rObj = VECTOR_ELT(computeList, index));
168                 PROTECT(s4class = STRING_ELT(getAttrib(rObj, install("class")), 0));
169                 omxCompute *compute = omxNewCompute(globalState, CHAR(s4class));
170                 compute->initFromFrontend(rObj);
171                 Global->computeList.push_back(compute);
172         }
173 }
174
175 void omxInitialMatrixAlgebraCompute() {
176         size_t numMats = globalState->matrixList.size();
177         int numAlgs = globalState->algebraList.size();
178
179         if(OMX_DEBUG) {mxLog("Completed Algebras and Matrices.  Beginning Initial Compute.");}
180         omxStateNextEvaluation(globalState);
181
182         for(size_t index = 0; index < numMats; index++) {
183                 omxRecompute(globalState->matrixList[index]);
184         }
185
186         for(int index = 0; index < numAlgs; index++) {
187                 omxRecompute(globalState->algebraList[index]);
188         }
189 }
190
191 /*
192 checkpointList is a list().  Each element refers to one checkpointing request.
193 Each interval request is a list of length 5.
194 The first element is an integer that specifies type: 0 = file, 1 = socket, 2=R_connection
195 For a file, the next two are the directory(string)  and file name (string).
196 For a socket, they are server (string) and port number (int).
197 For a connection, the next one is the R_connection SEXP object.
198 After that is an integer <type> specifier.  0 means minutes, 1 means iterations.
199 The last element is an integer count, indicating the number of <type>s per checkpoint.
200 */
201 void omxProcessCheckpointOptions(SEXP checkpointList) {
202         if(OMX_VERBOSE) { mxLog("Processing Checkpoint Requests.");}
203         globalState->numCheckpoints = length(checkpointList);
204         if(OMX_DEBUG) {mxLog("Found %d checkpoints.", globalState->numCheckpoints); }
205         globalState->checkpointList = (omxCheckpoint*) R_alloc(globalState->numCheckpoints, sizeof(omxCheckpoint));
206         SEXP nextLoc;
207
208         for(int index = 0; index < globalState->numCheckpoints; index++) {
209                 omxCheckpoint *oC = &(globalState->checkpointList[index]);
210
211                 /* Initialize Checkpoint object */
212                 oC->file = NULL;
213                 oC->connection = NULL;
214                 oC->time = 0;
215                 oC->numIterations = 0;
216                 oC->lastCheckpoint = 0;
217
218                 const char *pathName, *fileName;
219                 const char __attribute__((unused)) *serverName;
220
221                 PROTECT(nextLoc = VECTOR_ELT(checkpointList, index));
222                 int next = 0;
223                 oC->type = (omxCheckpointType) INTEGER(VECTOR_ELT(nextLoc, next++))[0];
224                 switch(oC->type) {
225                 case OMX_FILE_CHECKPOINT:{
226                         pathName = CHAR(STRING_ELT(VECTOR_ELT(nextLoc, next++), 0));                    // FIXME: Might need PROTECT()ion
227                         fileName = CHAR(STRING_ELT(VECTOR_ELT(nextLoc, next++), 0));
228                         char sep ='/';
229
230                         if(!isDir(pathName)) {
231                                 error("Unable to open directory %s for checkpoint storage.\n", pathName);
232                         }
233
234                         char* fullname = Calloc(strlen(pathName) + strlen(fileName) + 5, char);
235                         sprintf(fullname, "%s%c%s", pathName, sep, fileName);
236                         if(OMX_VERBOSE) { mxLog("Opening File: %s", fullname); }
237                         oC->file = fopen(fullname, "w");
238                         if(!oC->file) {
239                                 error("Unable to open file %s for checkpoint storage: %s.\n", fullname, strerror(errno));
240                         }
241                         Free(fullname);
242                         oC->saveHessian = FALSE;        // TODO: Decide if this should be true.
243                         break;}
244
245                 case OMX_CONNECTION_CHECKPOINT:{        // NYI :::DEBUG:::
246                         oC->connection = VECTOR_ELT(nextLoc, next++);
247                         error("Warning NYI: Socket checkpoints Not Yet Implemented.\n");
248                         oC->saveHessian = FALSE;
249                         break;}
250                 }
251
252                 int isCount = INTEGER(VECTOR_ELT(nextLoc, next++))[0];
253                 if(isCount) {
254                         oC->numIterations = INTEGER(AS_INTEGER(VECTOR_ELT(nextLoc, next++)))[0];
255                 } else {
256                         oC->time = REAL(AS_NUMERIC(VECTOR_ELT(nextLoc, next++)))[0] * 60;       // Constrained to seconds.
257                         if(oC->time < 1) oC->time = 1;                                                                          // Constrained to at least one.
258                 }
259         }
260 }
261
262 void omxProcessFreeVarList(SEXP fgNames, SEXP varList)
263 {
264         if(OMX_VERBOSE) { mxLog("Processing Free Parameters."); }
265
266         int ng = length(fgNames);
267         for (int gx=0; gx < ng; ++gx) {
268                 FreeVarGroup *fvg = new FreeVarGroup;
269                 fvg->name = CHAR(STRING_ELT(fgNames, gx));
270                 Global->freeGroup.push_back(fvg);
271         }
272
273         SEXP nextVar, nextLoc;
274         int numVars = length(varList);
275         for (int fx = 0; fx < numVars; fx++) {
276                 omxFreeVar *fv = new omxFreeVar;
277                 Global->freeGroup[0]->vars.push_back(fv);
278
279                 fv->name = CHAR(STRING_ELT(GET_NAMES(varList), fx));
280                 PROTECT(nextVar = VECTOR_ELT(varList, fx));
281
282                 PROTECT(nextLoc = VECTOR_ELT(nextVar, 0));
283                 fv->lbound = REAL(nextLoc)[0];
284                 if (ISNA(fv->lbound)) fv->lbound = NEG_INF;
285                 if (fv->lbound == 0.0) fv->lbound = 0.0;
286
287                 PROTECT(nextLoc = VECTOR_ELT(nextVar, 1));
288                 fv->ubound = REAL(nextLoc)[0];
289                 if (ISNA(fv->ubound)) fv->ubound = INF;
290                 if (fv->ubound == 0.0) fv->ubound = -0.0;
291
292                 PROTECT(nextLoc = VECTOR_ELT(nextVar, 2));
293                 int groupCount = length(nextLoc);
294                 if (groupCount == 0) error("Free variable %d not assigned to any group", fx);
295                 for (int gx=0; gx < groupCount; ++gx) {
296                         int group = INTEGER(nextLoc)[gx];
297                         if (group == 0) continue; // all variables are already in group 0 (default)
298                         Global->freeGroup[group]->vars.push_back(fv);
299                 }
300
301                 PROTECT(nextLoc = VECTOR_ELT(nextVar, 3));
302                 int numDeps = LENGTH(nextLoc);
303                 fv->numDeps = numDeps;
304                 fv->deps = (int*) R_alloc(numDeps, sizeof(int));
305                 for (int i = 0; i < numDeps; i++) {
306                         fv->deps[i] = INTEGER(nextLoc)[i];
307                 }
308
309                 int numLocs = length(nextVar) - 5;
310                 if(OMX_DEBUG) { 
311                         mxLog("Free parameter %d bounded (%f, %f): %d locations", fx, 
312                               fv->lbound, fv->ubound, numLocs);
313                 }
314                 for(int locIndex = 0; locIndex < numLocs; locIndex++) {
315                         PROTECT(nextLoc = VECTOR_ELT(nextVar, locIndex+4));
316                         int* theVarList = INTEGER(nextLoc);
317
318                         omxFreeVarLocation loc;
319                         loc.matrix = theVarList[0];
320                         loc.row = theVarList[1];
321                         loc.col = theVarList[2];
322
323                         fv->locations.push_back(loc);
324                 }
325                 PROTECT(nextLoc = VECTOR_ELT(nextVar, length(nextVar)-1));
326                 fv->start = REAL(nextLoc)[0];
327         }
328 }
329
330 /*
331         intervalList is a list().  Each element refers to one confidence interval request.
332         Each interval request is a length 5 vector of REAL.
333         The first three elements are the matrixPointer, Row, and Column of the element
334         for which bounds are to be calculated, and are cast to ints here for speed.
335         The last two are the upper and lower boundaries for the confidence space (respectively).
336 */
337 void omxProcessConfidenceIntervals(SEXP intervalList)  {
338         SEXP nextVar;
339         if(OMX_VERBOSE) { mxLog("Processing Confidence Interval Requests.");}
340         Global->numIntervals = length(intervalList);
341         if(OMX_DEBUG) {mxLog("Found %d requests.", Global->numIntervals); }
342         Global->intervalList = (omxConfidenceInterval*) R_alloc(Global->numIntervals, sizeof(omxConfidenceInterval));
343         for(int index = 0; index < Global->numIntervals; index++) {
344                 omxConfidenceInterval *oCI = &(Global->intervalList[index]);
345                 PROTECT(nextVar = VECTOR_ELT(intervalList, index));
346                 double* intervalInfo = REAL(nextVar);
347                 oCI->matrix = omxMatrixLookupFromState1( nextVar, globalState); // Expects an R object
348                 oCI->row = (int) intervalInfo[1];               // Cast to int in C to save memory/Protection ops
349                 oCI->col = (int) intervalInfo[2];               // Cast to int in C to save memory/Protection ops
350                 oCI->lbound = intervalInfo[3];
351                 oCI->ubound = intervalInfo[4];
352                 oCI->max = R_NaReal;                                    // NAs, in case something goes wrong
353                 oCI->min = R_NaReal;
354         }
355         if(OMX_VERBOSE) { mxLog("Processed."); }
356         if(OMX_DEBUG) { mxLog("%d intervals requested.", Global->numIntervals); }
357 }
358
359 void omxProcessConstraints(SEXP constraints)  {
360         int ncnln = 0; 
361         if(OMX_VERBOSE) { mxLog("Processing Constraints.");}
362         omxMatrix *arg1, *arg2;
363         SEXP nextVar, nextLoc;
364         globalState->numConstraints = length(constraints);
365         if(OMX_DEBUG) {mxLog("Found %d constraints.", globalState->numConstraints); }
366         globalState->conList = (omxConstraint*) R_alloc(globalState->numConstraints, sizeof(omxConstraint));
367         ncnln = 0;
368         for(int constraintIndex = 0; constraintIndex < globalState->numConstraints; constraintIndex++) {
369                 PROTECT(nextVar = VECTOR_ELT(constraints, constraintIndex));
370                 PROTECT(nextLoc = VECTOR_ELT(nextVar, 0));
371                 arg1 = omxMatrixLookupFromState1(nextLoc, globalState);
372                 PROTECT(nextLoc = VECTOR_ELT(nextVar, 1));
373                 arg2 = omxMatrixLookupFromState1(nextLoc, globalState);
374                 PROTECT(nextLoc = AS_INTEGER(VECTOR_ELT(nextVar, 2)));
375                 globalState->conList[constraintIndex].opCode = INTEGER(nextLoc)[0];
376                 omxMatrix *args[2] = {arg1, arg2};
377                 globalState->conList[constraintIndex].result = omxNewAlgebraFromOperatorAndArgs(10, args, 2, globalState); // 10 = binary subtract
378                 omxRecompute(globalState->conList[constraintIndex].result);
379                 int nrows = globalState->conList[constraintIndex].result->rows;
380                 int ncols = globalState->conList[constraintIndex].result->cols;
381                 globalState->conList[constraintIndex].size = nrows * ncols;
382                 ncnln += globalState->conList[constraintIndex].size;
383         }
384         if(OMX_VERBOSE) { mxLog("Processed."); }
385         if(OMX_DEBUG) { mxLog("%d effective constraints.", ncnln); }
386         globalState->ncnln = ncnln;
387 }
388
389 /*
390 *  Acknowledgement:
391 *  This function is duplicated from the function of the same name in the R source code.
392 *  The function appears in src/main/sysutils.c
393 *  Thanks to Michael Spiegel for finding it.
394 *  This code is licensed under the terms of the GNU General Public License.
395 */
396 static int isDir(const char *path)
397 {
398     struct stat sb;
399     int isdir = 0;
400     if(!path) return 0;
401     if(stat(path, &sb) == 0) {
402         isdir = (sb.st_mode & S_IFDIR) > 0; /* is a directory */
403         /* We want to know if the directory is writable by this user,
404            which mode does not tell us */
405         isdir &= (access(path, W_OK) == 0);
406     }
407     return isdir;
408 }
409