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