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