Better protect BA81 from parallel execution
[openmx:openmx.git] / src / omxFIMLFitFunction.cpp
1 /*
2  *  Copyright 2007-2014 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 "omxDefines.h"
18 #include "omxAlgebraFunctions.h"
19 #include "omxSymbolTable.h"
20 #include "omxData.h"
21 #include "omxFIMLFitFunction.h"
22 #include "omxFIMLSingleIteration.h"
23 #include "omxSadmvnWrapper.h"
24
25 #define max(a,b) \
26    ({ __typeof__ (a) _a = (a); \
27        __typeof__ (b) _b = (b); \
28      _a > _b ? _a : _b; })
29
30 /* FIML Function body */
31 void omxDestroyFIMLFitFunction(omxFitFunction *off) {
32         if(OMX_DEBUG) { mxLog("Destroying FIML fit function object."); }
33         omxFIMLFitFunction *argStruct = (omxFIMLFitFunction*) (off->argStruct);
34
35         omxFreeMatrix(argStruct->smallMeans);
36         omxFreeMatrix(argStruct->ordMeans);
37         omxFreeMatrix(argStruct->contRow);
38         omxFreeMatrix(argStruct->ordRow);
39         omxFreeMatrix(argStruct->ordCov);
40         omxFreeMatrix(argStruct->ordContCov);
41         omxFreeMatrix(argStruct->halfCov);
42         omxFreeMatrix(argStruct->reduceCov);
43
44         omxFreeMatrix(argStruct->smallRow);
45         omxFreeMatrix(argStruct->smallCov);
46         omxFreeMatrix(argStruct->RCX);
47         omxFreeMatrix(argStruct->rowLikelihoods);
48         omxFreeMatrix(argStruct->rowLogLikelihoods);
49 }
50
51 void omxPopulateFIMLAttributes(omxFitFunction *off, SEXP algebra) {
52         if(OMX_DEBUG) { mxLog("Populating FIML Attributes."); }
53         omxFIMLFitFunction *argStruct = ((omxFIMLFitFunction*)off->argStruct);
54         SEXP expCovExt, expMeanExt, rowLikelihoodsExt;
55         omxMatrix *expCovInt, *expMeanInt, *rowLikelihoodsInt;
56         expCovInt = argStruct->cov;
57         expMeanInt = argStruct->means;
58         rowLikelihoodsInt = argStruct->rowLikelihoods;
59
60         Rf_protect(expCovExt = Rf_allocMatrix(REALSXP, expCovInt->rows, expCovInt->cols));
61         for(int row = 0; row < expCovInt->rows; row++)
62                 for(int col = 0; col < expCovInt->cols; col++)
63                         REAL(expCovExt)[col * expCovInt->rows + row] =
64                                 omxMatrixElement(expCovInt, row, col);
65         if (expMeanInt != NULL) {
66                 Rf_protect(expMeanExt = Rf_allocMatrix(REALSXP, expMeanInt->rows, expMeanInt->cols));
67                 for(int row = 0; row < expMeanInt->rows; row++)
68                         for(int col = 0; col < expMeanInt->cols; col++)
69                                 REAL(expMeanExt)[col * expMeanInt->rows + row] =
70                                         omxMatrixElement(expMeanInt, row, col);
71         } else {
72                 Rf_protect(expMeanExt = Rf_allocMatrix(REALSXP, 0, 0));         
73         }
74         Rf_protect(rowLikelihoodsExt = Rf_allocVector(REALSXP, rowLikelihoodsInt->rows));
75         for(int row = 0; row < rowLikelihoodsInt->rows; row++)
76                 REAL(rowLikelihoodsExt)[row] = omxMatrixElement(rowLikelihoodsInt, row, 0);
77
78         Rf_setAttrib(algebra, Rf_install("expCov"), expCovExt);
79         Rf_setAttrib(algebra, Rf_install("expMean"), expMeanExt);
80         Rf_setAttrib(algebra, Rf_install("likelihoods"), rowLikelihoodsExt);
81
82         Rf_unprotect(3); // expCovExp, expCovInt, rowLikelihoodsExt
83 }
84
85 omxRListElement* omxSetFinalReturnsFIMLFitFunction(omxFitFunction *off, int *numReturns)
86 {
87         // DEPRECATED use omxPopulateFIMLAttributes
88         omxFIMLFitFunction* ofiml = (omxFIMLFitFunction *) (off->argStruct);
89
90         if(!ofiml->returnRowLikelihoods) return NULL;
91
92         omxRListElement *retVal = (omxRListElement*) R_alloc(1, sizeof(omxRListElement));
93         omxData* data = ofiml->data;
94         retVal[0].numValues = data->rows;
95         retVal[0].values = (double*) R_alloc(data->rows, sizeof(double));
96
97         return retVal;
98 }
99
100 void markDefVarDependencies(omxState* os, omxDefinitionVar* defVar) {
101
102         int numDeps = defVar->numDeps;
103         int *deps = defVar->deps;
104
105         for (int i = 0; i < numDeps; i++) {
106                 int value = deps[i];
107
108                 if(value < 0) {
109                         omxMarkDirty(os->matrixList[~value]);
110                 } else {
111                         omxMarkDirty(os->algebraList[value]);
112                 }
113         }
114
115 }
116
117 int handleDefinitionVarList(omxData* data, omxState *state, int row, omxDefinitionVar* defVars, double* oldDefs, int numDefs) {
118
119         if(OMX_DEBUG_ROWS(row)) { mxLog("Processing Definition Vars."); }
120         
121         int numVarsFilled = 0;
122
123         /* Fill in Definition Var Estimates */
124         for(int k = 0; k < numDefs; k++) {
125                 if(defVars[k].source != data) {
126                         Rf_error("Internal Rf_error: definition variable population into incorrect data source");
127                 }
128                 double newDefVar = omxDoubleDataElement(data, row, defVars[k].column);
129                 if(ISNA(newDefVar)) {
130                         Rf_error("Error: NA value for a definition variable is Not Yet Implemented.");
131                 }
132                 if(newDefVar == oldDefs[k]) {
133                         continue;       // NOTE: Potential speedup vs accuracy tradeoff here using epsilon comparison
134                 }
135                 oldDefs[k] = newDefVar;
136                 numVarsFilled++;
137
138                 for(int l = 0; l < defVars[k].numLocations; l++) {
139                         if(OMX_DEBUG_ROWS(row)) {
140                                 mxLog("Populating column %d (value %3.2f) into matrix %d.", defVars[k].column, omxDoubleDataElement(defVars[k].source, row, defVars[k].column), defVars[k].matrices[l]);
141                         }
142                         int matrixNumber = defVars[k].matrices[l];
143                         int matrow = defVars[k].rows[l];
144                         int matcol = defVars[k].cols[l];
145                         omxMatrix *matrix = state->matrixList[matrixNumber];
146                         omxSetMatrixElement(matrix, matrow, matcol, newDefVar);
147                 }
148                 markDefVarDependencies(state, &(defVars[k]));
149         }
150         return numVarsFilled;
151 }
152
153 static void omxCallJointFIMLFitFunction(omxFitFunction *off, int want, FitContext *fc)
154 {
155         // TODO: Figure out how to give access to other per-iteration structures.
156         // TODO: Current implementation is slow: update by filtering correlations and thresholds.
157         // TODO: Current implementation does not implement speedups for sorting.
158         // TODO: Current implementation may fail on all-continuous-missing or all-ordinal-missing rows.
159         
160         if (want & (FF_COMPUTE_PREOPTIMIZE)) return;
161
162     if(OMX_DEBUG) { 
163             mxLog("Beginning Joint FIML Evaluation.");
164     }
165         // Requires: Data, means, covariances, thresholds
166
167         int numDefs;
168
169         int returnRowLikelihoods = 0;
170
171         omxMatrix *cov, *means, *dataColumns;
172
173         omxThresholdColumn *thresholdCols;
174         omxData* data;
175         
176         omxExpectation* expectation;
177
178         omxFIMLFitFunction* ofiml = ((omxFIMLFitFunction*)off->argStruct);
179         omxMatrix* fitMatrix  = off->matrix;
180         omxState* parentState = fitMatrix->currentState;
181         int numChildren = parentState==globalState? globalState->childList.size() : 0;
182
183         cov             = ofiml->cov;
184         means           = ofiml->means;
185         data            = ofiml->data;                            //  read-only
186         numDefs         = ofiml->numDefs;                         //  read-only
187         dataColumns     = ofiml->dataColumns;
188         thresholdCols = ofiml->thresholdCols;
189
190         returnRowLikelihoods = ofiml->returnRowLikelihoods;   //  read-only
191         expectation = off->expectation;
192
193
194     if(numDefs == 0) {
195         if(OMX_DEBUG) {mxLog("Precalculating cov and means for all rows.");}
196                 omxExpectationRecompute(expectation);
197                 // MCN Also do the threshold formulae!
198                 
199                 for(int j=0; j < dataColumns->cols; j++) {
200                         int var = omxVectorElement(dataColumns, j);
201                         if(omxDataColumnIsFactor(data, j) && thresholdCols[var].numThresholds > 0) { // j is an ordinal column
202                                 omxMatrix* nextMatrix = thresholdCols[var].matrix;
203                                 omxRecompute(nextMatrix);
204                                 checkIncreasing(nextMatrix, thresholdCols[var].column);
205                                 for(int index = 0; index < numChildren; index++) {
206                                         omxMatrix *target = omxLookupDuplicateElement(parentState->childList[index], nextMatrix);
207                                         omxCopyMatrix(target, nextMatrix);
208                                 }
209             }
210         }
211                 for(int index = 0; index < numChildren; index++) {
212                         omxMatrix *childFit = omxLookupDuplicateElement(parentState->childList[index], fitMatrix);
213                         omxFIMLFitFunction* childOfiml = ((omxFIMLFitFunction*) childFit->fitFunction->argStruct);
214                         omxCopyMatrix(childOfiml->cov, cov);
215                         omxCopyMatrix(childOfiml->means, means);
216                 }
217         if(OMX_DEBUG) { omxPrintMatrix(cov, "Cov"); }
218         if(OMX_DEBUG) { omxPrintMatrix(means, "Means"); }
219     }
220
221         memset(ofiml->rowLogLikelihoods->data, 0, sizeof(double) * data->rows);
222     
223         int parallelism = (numChildren == 0) ? 1 : numChildren;
224
225         if (parallelism > data->rows) {
226                 parallelism = data->rows;
227         }
228
229         bool failed = false;
230         if (parallelism > 1) {
231                 int stride = (data->rows / parallelism);
232
233 #pragma omp parallel for num_threads(parallelism) reduction(||:failed)
234                 for(int i = 0; i < parallelism; i++) {
235                         omxMatrix *childMatrix = omxLookupDuplicateElement(parentState->childList[i], fitMatrix);
236                         omxFitFunction *childFit = childMatrix->fitFunction;
237                         if (i == parallelism - 1) {
238                                 failed |= omxFIMLSingleIterationJoint(fc, childFit, off, stride * i, data->rows - stride * i);
239                         } else {
240                                 failed |= omxFIMLSingleIterationJoint(fc, childFit, off, stride * i, stride);
241                         }
242                 }
243         } else {
244                 failed |= omxFIMLSingleIterationJoint(fc, off, off, 0, data->rows);
245         }
246         if (failed) {
247                 omxSetMatrixElement(off->matrix, 0, 0, NA_REAL);
248                 return;
249         }
250
251         if(!returnRowLikelihoods) {
252                 double val, sum = 0.0;
253                 // floating-point addition is not associative,
254                 // so we serialized the following reduction operation.
255                 for(int i = 0; i < data->rows; i++) {
256                         val = omxVectorElement(ofiml->rowLogLikelihoods, i);
257 //                      mxLog("%d , %f, %llx\n", i, val, *((unsigned long long*) &val));
258                         sum += val;
259                 }       
260                 if(OMX_VERBOSE || OMX_DEBUG) {mxLog("Total Likelihood is %3.3f", sum);}
261                 omxSetMatrixElement(off->matrix, 0, 0, sum);
262         }
263
264 }
265
266 static void omxCallFIMLFitFunction(omxFitFunction *off, int want, FitContext *fc)
267 {
268         if (want & (FF_COMPUTE_PREOPTIMIZE)) return;
269
270         if(OMX_DEBUG) { mxLog("Beginning FIML Evaluation."); }
271         // Requires: Data, means, covariances.
272         // Potential Problem: Definition variables currently are assumed to be at the end of the data matrix.
273
274         int numDefs, returnRowLikelihoods;      
275         omxExpectation* expectation;
276         
277         omxMatrix *cov, *means;//, *oldInverse;
278         omxData *data;
279
280         omxFIMLFitFunction* ofiml = ((omxFIMLFitFunction*) off->argStruct);
281         omxMatrix* objMatrix  = off->matrix;
282         omxState* parentState = objMatrix->currentState;
283         int numChildren = parentState==globalState? globalState->childList.size() : 0;
284
285         // Locals, for readability.  Should compile out.
286         cov             = ofiml->cov;
287         means           = ofiml->means;
288         data            = ofiml->data;                            //  read-only
289         numDefs         = ofiml->numDefs;                         //  read-only
290         returnRowLikelihoods = ofiml->returnRowLikelihoods;   //  read-only
291         expectation = off->expectation;
292
293         if(numDefs == 0 && strcmp(expectation->expType, "MxExpectationStateSpace")) {
294                 if(OMX_DEBUG) {mxLog("Precalculating cov and means for all rows.");}
295                 omxExpectationCompute(expectation, NULL);
296                 
297                 for(int index = 0; index < numChildren; index++) {
298                         omxMatrix *childFit = omxLookupDuplicateElement(parentState->childList[index], objMatrix);
299                         omxFIMLFitFunction* childOfiml = ((omxFIMLFitFunction*) childFit->fitFunction->argStruct);
300                         omxCopyMatrix(childOfiml->cov, cov);
301                         omxCopyMatrix(childOfiml->means, means);
302                 }
303
304                 if(OMX_DEBUG) { omxPrintMatrix(cov, "Cov"); }
305                 if(OMX_DEBUG) { omxPrintMatrix(means, "Means"); }
306         }
307
308         memset(ofiml->rowLogLikelihoods->data, 0, sizeof(double) * data->rows);
309     
310         int parallelism = (numChildren == 0) ? 1 : numChildren;
311
312         if (parallelism > data->rows) {
313                 parallelism = data->rows;
314         }
315
316         bool failed = false;
317         if (parallelism > 1) {
318                 int stride = (data->rows / parallelism);
319
320 #pragma omp parallel for num_threads(parallelism) reduction(||:failed)
321                 for(int i = 0; i < parallelism; i++) {
322                         omxMatrix *childMatrix = omxLookupDuplicateElement(parentState->childList[i], objMatrix);
323                         omxFitFunction *childFit = childMatrix->fitFunction;
324                         if (i == parallelism - 1) {
325                                 failed |= omxFIMLSingleIteration(fc, childFit, off, stride * i, data->rows - stride * i);
326                         } else {
327                                 failed |= omxFIMLSingleIteration(fc, childFit, off, stride * i, stride);
328                         }
329                 }
330         } else {
331                 failed |= omxFIMLSingleIteration(fc, off, off, 0, data->rows);
332         }
333         if (failed) {
334                 omxSetMatrixElement(off->matrix, 0, 0, NA_REAL);
335                 return;
336         }
337
338         if(!returnRowLikelihoods) {
339                 double val, sum = 0.0;
340                 // floating-point addition is not associative,
341                 // so we serialized the following reduction operation.
342                 for(int i = 0; i < data->rows; i++) {
343                         val = omxVectorElement(ofiml->rowLogLikelihoods, i);
344 //                      mxLog("%d , %f, %llx\n", i, val, *((unsigned long long*) &val));
345                         sum += val;
346                 }       
347                 if(OMX_VERBOSE || OMX_DEBUG) {mxLog("Total Likelihood is %3.3f", sum);}
348                 omxSetMatrixElement(off->matrix, 0, 0, sum);
349         }
350 }
351
352 static void omxCallFIMLOrdinalFitFunction(omxFitFunction *off, int want, FitContext *fc)
353 {
354         if (want & (FF_COMPUTE_PREOPTIMIZE)) return;
355
356         /* TODO: Current implementation is slow: update by filtering correlations and thresholds. */
357         if(OMX_DEBUG) { mxLog("Beginning Ordinal FIML Evaluation.");}
358         // Requires: Data, means, covariances, thresholds
359
360         int numDefs;
361         int returnRowLikelihoods = 0;
362
363         omxMatrix *cov, *means, *dataColumns;
364         omxThresholdColumn *thresholdCols;
365         omxData* data;
366         double *corList, *weights;
367         
368         omxExpectation* expectation;
369
370         omxFIMLFitFunction* ofiml = ((omxFIMLFitFunction*)off->argStruct);
371         omxMatrix* objMatrix  = off->matrix;
372         omxState* parentState = objMatrix->currentState;
373         int numChildren = parentState==globalState? globalState->childList.size() : 0;
374
375         // Locals, for readability.  Compiler should cut through this.
376         cov             = ofiml->cov;
377         means           = ofiml->means;
378         data            = ofiml->data;
379         dataColumns     = ofiml->dataColumns;
380         numDefs         = ofiml->numDefs;
381
382         corList         = ofiml->corList;
383         weights         = ofiml->weights;
384         thresholdCols = ofiml->thresholdCols;
385         returnRowLikelihoods = ofiml->returnRowLikelihoods;
386         
387         expectation = off->expectation;
388         
389         if(numDefs == 0) {
390                 if(OMX_DEBUG_ALGEBRA) { mxLog("No Definition Vars: precalculating."); }
391                 omxExpectationCompute(expectation, NULL);
392                 for(int j = 0; j < dataColumns->cols; j++) {
393                         if(thresholdCols[j].numThresholds > 0) { // Actually an ordinal column
394                                 omxMatrix* nextMatrix = thresholdCols[j].matrix;
395                                 omxRecompute(nextMatrix);
396                                 checkIncreasing(nextMatrix, thresholdCols[j].column);
397                                 for(int index = 0; index < numChildren; index++) {
398                                         omxMatrix *target = omxLookupDuplicateElement(parentState->childList[index], nextMatrix);
399                                         omxCopyMatrix(target, nextMatrix);
400                                 }
401                         }
402                 }
403                 omxStandardizeCovMatrix(cov, corList, weights); // Calculate correlation and covariance
404                 for(int index = 0; index < numChildren; index++) {
405                         omxMatrix *childFit = omxLookupDuplicateElement(parentState->childList[index], objMatrix);
406                         omxFIMLFitFunction* childOfiml = ((omxFIMLFitFunction*) childFit->fitFunction->argStruct);
407                         omxCopyMatrix(childOfiml->cov, cov);
408                         omxCopyMatrix(childOfiml->means, means);
409                         memcpy(childOfiml->weights, weights, sizeof(double) * cov->rows);
410                         memcpy(childOfiml->corList, corList, sizeof(double) * (cov->rows * (cov->rows - 1)) / 2);
411                 }
412         }
413
414         memset(ofiml->rowLogLikelihoods->data, 0, sizeof(double) * data->rows);
415
416         int parallelism = (numChildren == 0) ? 1 : numChildren;
417
418         if (parallelism > data->rows) {
419                 parallelism = data->rows;
420         }
421
422         bool failed = false;
423         if (parallelism > 1) {
424                 int stride = (data->rows / parallelism);
425
426 #pragma omp parallel for num_threads(parallelism) reduction(||:failed)
427                 for(int i = 0; i < parallelism; i++) {
428                         omxMatrix *childMatrix = omxLookupDuplicateElement(parentState->childList[i], objMatrix);
429                         omxFitFunction *childFit = childMatrix->fitFunction;
430                         if (i == parallelism - 1) {
431                                 failed |= omxFIMLSingleIterationOrdinal(fc, childFit, off, stride * i, data->rows - stride * i);
432                         } else {
433                                 failed |= omxFIMLSingleIterationOrdinal(fc, childFit, off, stride * i, stride);
434                         }
435                 }
436         } else {
437                 failed |= omxFIMLSingleIterationOrdinal(fc, off, off, 0, data->rows);
438         }
439         if (failed) {
440                 omxSetMatrixElement(off->matrix, 0, 0, NA_REAL);
441                 return;
442         }
443
444         if(!returnRowLikelihoods) {
445                 double val, sum = 0.0;
446                 // floating-point addition is not associative,
447                 // so we serialized the following reduction operation.
448                 for(int i = 0; i < data->rows; i++) {
449                         val = omxVectorElement(ofiml->rowLogLikelihoods, i);
450 //                      mxLog("%d , %f, %llx\n", i, val, *((unsigned long long*) &val));
451                         sum += val;
452                 }       
453                 if(OMX_VERBOSE || OMX_DEBUG) {mxLog("Total Likelihood is %3.3f", sum);}
454                 omxSetMatrixElement(off->matrix, 0, 0, sum);
455         }
456 }
457
458 void omxInitFIMLFitFunction(omxFitFunction* off)
459 {
460         if(OMX_DEBUG) {
461                 mxLog("Initializing FIML fit function function.");
462         }
463         SEXP rObj = off->rObj;
464
465         int numOrdinal = 0, numContinuous = 0;
466         omxMatrix *cov, *means;
467
468         omxFIMLFitFunction *newObj = (omxFIMLFitFunction*) R_alloc(1, sizeof(omxFIMLFitFunction));
469         omxExpectation* expectation = off->expectation;
470         if(expectation == NULL) {
471                 omxRaiseError("FIML cannot fit without model expectations.");
472                 return;
473         }
474
475         cov = omxGetExpectationComponent(expectation, off, "cov");
476         if(cov == NULL) { 
477                 omxRaiseError("No covariance expectation in FIML evaluation.");
478                 return;
479         }
480
481         means = omxGetExpectationComponent(expectation, off, "means");
482         
483         if(means == NULL) { 
484                 omxRaiseError("No means model in FIML evaluation.");
485                 return;
486         }
487
488         if(OMX_DEBUG) {
489                 mxLog("FIML Initialization Completed.");
490         }
491         
492     newObj->cov = cov;
493     newObj->means = means;
494     newObj->smallMeans = NULL;
495     newObj->ordMeans   = NULL;
496     newObj->contRow    = NULL;
497     newObj->ordRow     = NULL;
498     newObj->ordCov     = NULL;
499     newObj->ordContCov = NULL;
500     newObj->halfCov    = NULL;
501     newObj->reduceCov  = NULL;
502     
503     /* Set default FitFunction calls to FIML FitFunction Calls */
504         off->computeFun = omxCallFIMLFitFunction;
505         off->fitType = "imxFitFunctionFIML";
506         off->setFinalReturns = omxSetFinalReturnsFIMLFitFunction;
507         off->destructFun = omxDestroyFIMLFitFunction;
508         off->populateAttrFun = omxPopulateFIMLAttributes;
509
510         if(OMX_DEBUG) {
511                 mxLog("Accessing data source.");
512         }
513         newObj->data = off->expectation->data;
514
515         if(OMX_DEBUG) {
516                 mxLog("Accessing row likelihood option.");
517         }
518         newObj->returnRowLikelihoods = Rf_asInteger(R_do_slot(rObj, Rf_install("vector")));
519         if(newObj->returnRowLikelihoods) {
520                 omxResizeMatrix(off->matrix, newObj->data->rows, 1);
521         }
522     newObj->rowLikelihoods = omxInitMatrix(newObj->data->rows, 1, TRUE, off->matrix->currentState);
523     newObj->rowLogLikelihoods = omxInitMatrix(newObj->data->rows, 1, TRUE, off->matrix->currentState);
524
525         if(OMX_DEBUG) {
526                 mxLog("Accessing variable mapping structure.");
527         }
528         newObj->dataColumns = off->expectation->dataColumns;
529
530         if(OMX_DEBUG) {
531                 mxLog("Accessing Threshold matrix.");
532         }
533         newObj->thresholdCols = off->expectation->thresholds;
534         numOrdinal = off->expectation->numOrdinal;
535         numContinuous = newObj->dataColumns->cols - off->expectation->numOrdinal;
536
537         omxSetContiguousDataColumns(&(newObj->contiguous), newObj->data, newObj->dataColumns);
538         
539         newObj->numDefs = off->expectation->numDefs;
540         newObj->defVars = off->expectation->defVars;
541         newObj->oldDefs = (double *) R_alloc(newObj->numDefs, sizeof(double));          // Storage for Def Vars
542
543         if(OMX_DEBUG) {
544                 mxLog("Accessing definition variables structure.");
545         }
546         newObj->oldDefs = (double *) R_alloc(newObj->numDefs, sizeof(double));          // Storage for Def Vars
547         memset(newObj->oldDefs, NA_REAL, sizeof(double) * newObj->numDefs);                     // Does this work?
548         // for(nextDef = 0; nextDef < newObj->numDefs; nextDef++) {
549         //      newObj->oldDefs[nextDef] = NA_REAL;                                     // Def Vars default to NA
550         // }
551
552     /* Temporary storage for calculation */
553     int covCols = newObj->cov->cols;
554         if(OMX_DEBUG){mxLog("Number of columns found is %d", covCols);}
555     // int ordCols = omxDataNumFactor(newObj->data);        // Unneeded, since we don't use it.
556     // int contCols = omxDataNumNumeric(newObj->data);
557     newObj->smallRow = omxInitMatrix(1, covCols, TRUE, off->matrix->currentState);
558     newObj->smallCov = omxInitMatrix(covCols, covCols, TRUE, off->matrix->currentState);
559     newObj->RCX = omxInitMatrix(1, covCols, TRUE, off->matrix->currentState);
560 //  newObj->zeros = omxInitMatrix(1, newObj->cov->cols, TRUE, off->matrix->currentState);
561
562     omxCopyMatrix(newObj->smallCov, newObj->cov);          // Will keep its aliased state from here on.
563     off->argStruct = (void*)newObj;
564
565     if(numOrdinal > 0 && numContinuous <= 0) {
566         if(OMX_DEBUG) {
567             mxLog("Ordinal Data detected.  Using Ordinal FIML.");
568         }
569         newObj->weights = (double*) R_alloc(covCols, sizeof(double));
570         newObj->smallMeans = omxInitMatrix(covCols, 1, TRUE, off->matrix->currentState);
571         omxCopyMatrix(newObj->smallMeans, newObj->means);
572         newObj->corList = (double*) R_alloc(covCols * (covCols + 1) / 2, sizeof(double));
573         newObj->smallCor = (double*) R_alloc(covCols * (covCols + 1) / 2, sizeof(double));
574         newObj->lThresh = (double*) R_alloc(covCols, sizeof(double));
575         newObj->uThresh = (double*) R_alloc(covCols, sizeof(double));
576         newObj->Infin = (int*) R_alloc(covCols, sizeof(int));
577
578         off->computeFun = omxCallFIMLOrdinalFitFunction;
579     } else if(numOrdinal > 0) {
580         if(OMX_DEBUG) {
581             mxLog("Ordinal and Continuous Data detected.  Using Joint Ordinal/Continuous FIML.");
582         }
583
584         newObj->weights = (double*) R_alloc(covCols, sizeof(double));
585         newObj->smallMeans = omxInitMatrix(covCols, 1, TRUE, off->matrix->currentState);
586         newObj->ordMeans = omxInitMatrix(covCols, 1, TRUE, off->matrix->currentState);
587         newObj->contRow = omxInitMatrix(covCols, 1, TRUE, off->matrix->currentState);
588         newObj->ordRow = omxInitMatrix(covCols, 1, TRUE, off->matrix->currentState);
589         newObj->ordCov = omxInitMatrix(covCols, covCols, TRUE, off->matrix->currentState);
590         newObj->ordContCov = omxInitMatrix(covCols, covCols, TRUE, off->matrix->currentState);
591         newObj->halfCov = omxInitMatrix(covCols, covCols, TRUE, off->matrix->currentState);
592         newObj->reduceCov = omxInitMatrix(covCols, covCols, TRUE, off->matrix->currentState);
593         omxCopyMatrix(newObj->smallMeans, newObj->means);
594         omxCopyMatrix(newObj->ordMeans, newObj->means);
595         omxCopyMatrix(newObj->contRow, newObj->smallRow );
596         omxCopyMatrix(newObj->ordRow, newObj->smallRow );
597         omxCopyMatrix(newObj->ordCov, newObj->cov);
598         omxCopyMatrix(newObj->ordContCov, newObj->cov);
599         omxCopyMatrix(newObj->smallMeans, newObj->means);
600         newObj->corList = (double*) R_alloc(covCols * (covCols + 1) / 2, sizeof(double));
601         newObj->lThresh = (double*) R_alloc(covCols, sizeof(double));
602         newObj->uThresh = (double*) R_alloc(covCols, sizeof(double));
603         newObj->Infin = (int*) R_alloc(covCols, sizeof(int));
604
605         off->computeFun = omxCallJointFIMLFitFunction;
606     }
607 }