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