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