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