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