Added omxState duplication process.
[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 #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);
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);
106                 } else {                                                                // This is an algebra spec.
107                         PROTECT(nextAlg = VECTOR_ELT(nextAlgTuple, 0));
108                         omxFillMatrixFromRPrimitive(currentState->algebraList[index],
109                                 nextAlg, currentState);
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, *serverName;
184
185                 PROTECT(nextLoc = VECTOR_ELT(checkpointList, index));
186                 int next = 0;
187                 oC->type = INTEGER(VECTOR_ELT(nextLoc, next++))[0];
188                 switch(oC->type) {
189                 case OMX_FILE_CHECKPOINT:
190                         pathName = CHAR(STRING_ELT(VECTOR_ELT(nextLoc, next++), 0));                    // FIXME: Might need PROTECT()ion
191                         fileName = CHAR(STRING_ELT(VECTOR_ELT(nextLoc, next++), 0));
192                         char sep ='/';
193
194                         if(!isDir(pathName)) {
195                                 error("Unable to open directory %s for checkpoint storage.\n", pathName);
196                         }
197
198                         char* fullname = Calloc(strlen(pathName) + strlen(fileName) + 5, char);
199                         sprintf(fullname, "%s%c%s", pathName, sep, fileName);
200                         if(OMX_VERBOSE) { Rprintf("Opening File: %s\n", fullname); }
201                         oC->file = fopen(fullname, "w");
202                         if(!oC->file) {
203                                 error("Unable to open file %s for checkpoint storage: %s.\n", fullname, strerror(errno));
204                         }
205                         Free(fullname);
206                         oC->saveHessian = FALSE;        // TODO: Decide if this should be true.
207                         break;
208
209                 case OMX_SOCKET_CHECKPOINT:
210                         serverName = CHAR(VECTOR_ELT(nextLoc, next++));
211                         // int portno = INTEGER(AS_INTEGER(VECTOR_ELT(nextLoc, next++)))[0]; // Commented to suppress "unused variable" warning.
212                         Rprintf("Warning NYI: Socket checkpoints Not Yet Implemented.\n");
213                         oC->saveHessian = FALSE;
214                         break;
215
216                 case OMX_CONNECTION_CHECKPOINT: // NYI :::DEBUG:::
217                         oC->connection = VECTOR_ELT(nextLoc, next++);
218                         Rprintf("Warning NYI: Socket checkpoints Not Yet Implemented.\n");
219                         oC->saveHessian = FALSE;
220                         break;
221                 }
222
223                 int isCount = INTEGER(VECTOR_ELT(nextLoc, next++))[0];
224                 if(isCount) {
225                         oC->numIterations = INTEGER(AS_INTEGER(VECTOR_ELT(nextLoc, next++)))[0];
226                 } else {
227                         oC->time = REAL(AS_NUMERIC(VECTOR_ELT(nextLoc, next++)))[0] * 60;       // Constrained to seconds.
228                         if(oC->time < 1) oC->time = 1;                                                                          // Constrained to at least one.
229                 }
230
231                 UNPROTECT(1); /* nextLoc */
232
233         }
234 }
235
236 /*
237 varList is a list().  Each element of this list corresponds to one free parameter.
238 Each free parameter is a list.  The first element of this list is the lower bound.
239 The second element of the list is the upper bound.  The remaining elements of this
240 list are 3-tuples.  These 3-tuples are (mxIndex, row, col).
241 */
242 void omxProcessFreeVarList(SEXP varList, int n) {
243         SEXP nextVar, nextLoc;
244         if(OMX_VERBOSE) { Rprintf("Processing Free Parameters.\n"); }
245         currentState->freeVarList = (omxFreeVar*) R_alloc (n, sizeof (omxFreeVar));                     // Data for replacement of free vars
246         for(int freeVarIndex = 0; freeVarIndex < n; freeVarIndex++) {
247                 PROTECT(nextVar = VECTOR_ELT(varList, freeVarIndex));
248                 int numLocs = length(nextVar) - 2;
249                 currentState->freeVarList[freeVarIndex].numLocations = numLocs;
250                 currentState->freeVarList[freeVarIndex].location = (double**) R_alloc(numLocs, sizeof(double*));
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                         double* locationMatrixElement = omxLocationOfMatrixElement(currentState->matrixList[theMat], theRow, theCol);
283                         currentState->freeVarList[freeVarIndex].location[locIndex] = locationMatrixElement;
284                         currentState->freeVarList[freeVarIndex].matrices[locIndex] = theMat;
285                         currentState->freeVarList[freeVarIndex].row[locIndex] = theRow;
286                         currentState->freeVarList[freeVarIndex].col[locIndex] = theCol;
287                         UNPROTECT(1); // nextLoc
288                 }
289                 UNPROTECT(1); // nextVar
290         }
291 }
292
293 /*
294         intervalList is a list().  Each element refers to one confidence interval request.
295         Each interval request is a length 5 vector of REAL.
296         The first three elements are the matrixPointer, Row, and Column of the element
297         for which bounds are to be calculated, and are cast to ints here for speed.
298         The last two are the upper and lower boundaries for the confidence space (respectively).
299 */
300 void omxProcessConfidenceIntervals(SEXP intervalList)  {
301         SEXP nextVar;
302         if(OMX_VERBOSE) { Rprintf("Processing Confidence Interval Requests.\n");}
303         currentState->numIntervals = length(intervalList);
304         if(OMX_DEBUG) {Rprintf("Found %d requests.\n", currentState->numIntervals); }
305         currentState->intervalList = (omxConfidenceInterval*) R_alloc(currentState->numIntervals, sizeof(omxConfidenceInterval));
306         for(int index = 0; index < currentState->numIntervals; index++) {
307                 omxConfidenceInterval *oCI = &(currentState->intervalList[index]);
308                 PROTECT(nextVar = VECTOR_ELT(intervalList, index));
309                 double* intervalInfo = REAL(nextVar);
310                 oCI->matrix = omxNewMatrixFromMxIndex( nextVar, currentState);  // Expects an R object
311                 oCI->row = (int) intervalInfo[1];               // Cast to int in C to save memory/Protection ops
312                 oCI->col = (int) intervalInfo[2];               // Cast to int in C to save memory/Protection ops
313                 oCI->lbound = intervalInfo[3];
314                 oCI->ubound = intervalInfo[4];
315                 UNPROTECT(1);
316                 oCI->max = R_NaReal;                                    // NAs, in case something goes wrong
317                 oCI->min = R_NaReal;
318         }
319         if(OMX_VERBOSE) { Rprintf("Processed.\n"); }
320         if(OMX_DEBUG) { Rprintf("%d intervals requested.\n", currentState->numIntervals); }
321 }
322
323 int omxProcessConstraints(SEXP constraints)  {
324         int ncnln = 0; 
325         if(OMX_VERBOSE) { Rprintf("Processing Constraints.\n");}
326         omxMatrix *arg1, *arg2;
327         SEXP nextVar, nextLoc;
328         currentState->numConstraints = length(constraints);
329         if(OMX_DEBUG) {Rprintf("Found %d constraints.\n", currentState->numConstraints); }
330         currentState->conList = (omxConstraint*) R_alloc(currentState->numConstraints, sizeof(omxConstraint));
331         ncnln = 0;
332         for(int constraintIndex = 0; constraintIndex < currentState->numConstraints; constraintIndex++) {
333                 PROTECT(nextVar = VECTOR_ELT(constraints, constraintIndex));
334                 PROTECT(nextLoc = VECTOR_ELT(nextVar, 0));
335                 arg1 = omxNewMatrixFromMxIndex(nextLoc, currentState);
336                 PROTECT(nextLoc = VECTOR_ELT(nextVar, 1));
337                 arg2 = omxNewMatrixFromMxIndex(nextLoc, currentState);
338                 PROTECT(nextLoc = AS_INTEGER(VECTOR_ELT(nextVar, 2)));
339                 currentState->conList[constraintIndex].opCode = INTEGER(nextLoc)[0];
340                 UNPROTECT(4);
341                 omxMatrix *args[2] = {arg1, arg2};
342                 currentState->conList[constraintIndex].result = omxNewAlgebraFromOperatorAndArgs(10, args, 2, currentState); // 10 = binary subtract
343                 omxRecompute(currentState->conList[constraintIndex].result);
344                 int nrows = currentState->conList[constraintIndex].result->rows;
345                 int ncols = currentState->conList[constraintIndex].result->cols;
346                 currentState->conList[constraintIndex].size = nrows * ncols;
347                 ncnln += currentState->conList[constraintIndex].size;
348         }
349         if(OMX_VERBOSE) { Rprintf("Processed.\n"); }
350         if(OMX_DEBUG) { Rprintf("%d effective constraints.\n", ncnln); }
351         return(ncnln);
352 }
353
354 void omxSetupBoundsAndConstraints(double * bl, double * bu, int n, int nclin) {
355         /* Set min and max limits */
356         for(int index = 0; index < n; index++) {
357                 bl[index] = currentState->freeVarList[index].lbound;                            // -Infinity'd be -10^20.
358                 bu[index] = currentState->freeVarList[index].ubound;                            // Infinity would be at 10^20.
359         }
360
361         for(int index = n; index < n+nclin; index++) {                                          // At present, nclin == 0.
362                 bl[index] = NEG_INF;                                                    // Linear constraints have no bounds.
363                 bu[index] = INF;                                                                // Because there are no linear constraints.
364         }                                                                                                   // But if there were, they would go here.
365
366         int index = n + nclin;
367         for(int constraintIndex = 0; constraintIndex < currentState->numConstraints; constraintIndex++) {               // Nonlinear constraints:
368                 if(OMX_DEBUG) { Rprintf("Constraint %d: ", constraintIndex);}
369                 switch(currentState->conList[constraintIndex].opCode) {
370                 case 0:                                                                 // Less than: Must be strictly less than 0.
371                         if(OMX_DEBUG) { Rprintf("Bounded at (-INF, 0).\n");}
372                         for(int offset = 0; offset < currentState->conList[constraintIndex].size; offset++) {
373                                 bl[index] = NEG_INF;
374                                 bu[index] = -0.0;
375                                 index++;
376                         }
377                         break;
378                 case 1:                                                                 // Equal: Must be roughly equal to 0.
379                         if(OMX_DEBUG) { Rprintf("Bounded at (-0, 0).\n");}
380                         for(int offset = 0; offset < currentState->conList[constraintIndex].size; offset++) {
381                                 bl[index] = -0.0;
382                                 bu[index] = 0.0;
383                                 index++;
384                         }
385                         break;
386                 case 2:                                                                 // Greater than: Must be strictly greater than 0.
387                         if(OMX_DEBUG) { Rprintf("Bounded at (0, INF).\n");}
388                         for(int offset = 0; offset < currentState->conList[constraintIndex].size; offset++) {
389                                 if(OMX_DEBUG) { Rprintf("\tBounds set for constraint %d.%d.\n", constraintIndex, offset);}
390                                 bl[index] = 0.0;
391                                 bu[index] = INF;
392                                 index++;
393                         }
394                         break;
395                 default:
396                         if(OMX_DEBUG) { Rprintf("Bounded at (-INF, INF).\n");}
397                         for(int offset = 0; offset < currentState->conList[constraintIndex].size; offset++) {
398                                 bl[index] = NEG_INF;
399                                 bu[index] = INF;
400                                 index++;
401                         }
402                         break;
403                 }
404         }
405 }
406
407 /*
408 *  Acknowledgement:
409 *  This function is duplicated from the function of the same name in the R source code.
410 *  The function appears in src/main/sysutils.c
411 *  Thanks to Michael Spiegel for finding it.
412 */
413 static int isDir(const char *path)
414 {
415     struct stat sb;
416     int isdir = 0;
417     if(!path) return 0;
418     if(stat(path, &sb) == 0) {
419         isdir = (sb.st_mode & S_IFDIR) > 0; /* is a directory */
420         /* We want to know if the directory is writable by this user,
421            which mode does not tell us */
422         isdir &= (access(path, W_OK) == 0);
423     }
424     return isdir;
425 }
426