Store dataList in std::vector
[openmx:openmx.git] / src / omxData.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 /***********************************************************
18  *
19  *  omxData.cc
20  *
21  *  Created: Timothy R. Brick   Date: 2009-07-15
22  *
23  *      Contains code for the omxData class
24  *   omxData objects hold data in whatever form it takes
25  *
26  **********************************************************/
27 #include "omxData.h"
28 #include "npsolWrap.h"
29 #include "omxGlobalState.h"
30
31 omxData* omxInitData(omxState* os) {
32
33         omxData *od = Calloc(1, omxData);
34
35         od->currentState = os;
36
37         if(OMX_DEBUG) {Rprintf("Data's state object is at 0x%x.\n", od->currentState);}
38
39         if (os != globalState) error("Too late to create omxData");
40         os->dataList.push_back(od);
41
42         return od;
43
44 }
45
46 omxData* omxDataLookupFromState(SEXP dataObject, omxState* state) {
47         int dataIdx = INTEGER(dataObject)[0];
48
49         return state->dataList[dataIdx];
50 }
51
52 omxData* omxNewDataFromMxData(SEXP dataObject, omxState* state) {
53         if(OMX_DEBUG) {Rprintf("Initializing data Element.\n");}
54         if(dataObject == NULL) {
55                 error("Null Data Object detected.  This is an internal error, and should be reported on the forums.\n");
56         }
57
58         omxData* od = omxInitData(state);
59
60         SEXP dataLoc, dataVal;
61         int numCols;
62
63         // PARSE MxData Structure
64         if(OMX_DEBUG) {Rprintf("Processing Data Type.\n");}
65         PROTECT(dataLoc = GET_SLOT(dataObject, install("type")));
66         if(dataLoc == NULL) { error("Data has no type.  Sorry.\nThis is an internal error, and should be reported on the forums.\n");}
67         PROTECT(dataVal = STRING_ELT(dataLoc,0));
68         od->_type = CHAR(dataVal);
69         if(OMX_DEBUG) {Rprintf("Element is type %s.\n", od->_type);}
70
71         PROTECT(dataLoc = GET_SLOT(dataObject, install("observed")));
72         if(OMX_DEBUG) {Rprintf("Processing Data Elements.\n");}
73         if(isFrame(dataLoc)) {
74                 if(OMX_DEBUG) {Rprintf("Data is a frame.\n");}
75                 // Process Data Frame into Columns
76                 od->cols = length(dataLoc);
77                 if(OMX_DEBUG) {Rprintf("Data has %d columns.\n", od->cols);}
78                 numCols = od->cols;
79                 od->realData = (double**) R_alloc(numCols, sizeof(double*));
80                 od->intData = (int**) R_alloc(numCols, sizeof(int*));
81                 od->location = (int*) R_alloc(numCols, sizeof(int));
82                 for(int j = 0; j < numCols; j++) {
83                         SEXP rcol;
84                         PROTECT(rcol = VECTOR_ELT(dataLoc, j));
85                         if(isFactor(rcol)) {
86                                 if(OMX_DEBUG) {Rprintf("Column %d is a factor.\n", j);}
87                                 od->intData[j] = INTEGER(rcol);
88                                 od->location[j] = ~j;
89                                 od->numFactor++;
90                         } else if (isInteger(rcol)) {
91                                 error("Internal error: Column %d is in integer format.", j);
92                         } else {
93                                 if(OMX_DEBUG) {Rprintf("Column %d is a numeric.\n", j);}
94                                 od->realData[j] = REAL(rcol);
95                                 od->location[j] = j;
96                                 od->numNumeric++;
97                         }
98                 }
99                 od->rows = length(VECTOR_ELT(dataLoc, 0));
100                 if(OMX_DEBUG) {Rprintf("And %d rows.\n", od->rows);}
101         } else {
102                 if(OMX_DEBUG) {Rprintf("Data contains a matrix.\n");}
103                 od->dataMat = omxNewMatrixFromRPrimitive(dataLoc, od->currentState, 0, 0);
104                 
105                 if (od->dataMat->colMajor && strncmp(od->_type, "raw", 3) == 0) { 
106                         omxToggleRowColumnMajor(od->dataMat);
107                 }
108                 od->cols = od->dataMat->cols;
109                 od->rows = od->dataMat->rows;
110                 od->numNumeric = od->cols;
111         }
112
113         if(OMX_DEBUG) {Rprintf("Processing Means Matrix.\n");}
114         PROTECT(dataLoc = GET_SLOT(dataObject, install("means")));
115         od->meansMat = omxNewMatrixFromRPrimitive(dataLoc, od->currentState, 0, 0);
116         if(od->meansMat->rows == 1 && od->meansMat->cols == 1 && 
117            (!R_finite(omxMatrixElement(od->meansMat, 0, 0)) ||
118             !isfinite(omxMatrixElement(od->meansMat, 0, 0)))) {
119                 omxFreeMatrixData(od->meansMat); // Clear just-allocated memory.
120                 od->meansMat = NULL;  // 1-by-1 matrix of NAs is a null means matrix.
121                 // FIXME: The above check may cause problems for dynamic data if the means
122                 //          originally is a 1x1 that has not yet been calculated.  This should be
123                 //          adjusted.
124         }
125         
126         if(OMX_DEBUG) {
127                 if(od->meansMat == NULL) {Rprintf("No means found.\n");}
128                 else {omxPrint(od->meansMat, "Means Matrix is:");}
129         }
130
131         if(strncmp(od->_type, "raw", 3) != 0) {
132                 if(OMX_DEBUG) {Rprintf("Processing Observation Count.\n");}
133                 PROTECT(dataLoc = GET_SLOT(dataObject, install("numObs")));
134                 od->numObs = REAL(dataLoc)[0];
135         } else {
136                 od->numObs = od->rows;
137                 if(OMX_DEBUG) {Rprintf("Processing presort metadata.\n");}
138                 /* For raw data, process sorting metadata. */
139                 // Process unsorted indices:  // TODO: Generate reverse lookup table
140                 PROTECT(dataLoc = GET_SLOT(dataObject, install("indexVector")));
141                 od->indexVector = INTEGER(dataLoc);
142                 if(od->indexVector[0] == R_NaInt) od->indexVector = NULL;
143                 // Process pre-computed identicality checks
144                 if(OMX_DEBUG) {Rprintf("Processing definition variable identicality.\n");}
145                 PROTECT(dataLoc = GET_SLOT(dataObject, install("identicalDefVars")));
146                 od->identicalDefs = INTEGER(dataLoc);
147                 if(od->identicalDefs[0] == R_NaInt) od->identicalDefs = NULL;
148                 if(OMX_DEBUG) {Rprintf("Processing missingness identicality.\n");}
149                 PROTECT(dataLoc = GET_SLOT(dataObject, install("identicalMissingness")));
150                 od->identicalMissingness = INTEGER(dataLoc);
151                 if(od->identicalMissingness[0] == R_NaInt) od->identicalMissingness = NULL;
152                 if(OMX_DEBUG) {Rprintf("Processing row identicality.\n");}
153                 PROTECT(dataLoc = GET_SLOT(dataObject, install("identicalRows")));
154                 od->identicalRows = INTEGER(dataLoc);
155                 if(od->identicalRows[0] == R_NaInt) od->identicalRows = NULL;
156         }
157
158         return od;
159 }
160
161 void resetDefinitionVariables(double *oldDefs, int numDefs) {
162         int nextDef;
163
164         for(nextDef = 0; nextDef < numDefs; nextDef++) {
165                 oldDefs[nextDef] = NA_REAL;                                     // Def Vars default to NA
166         }
167
168 }
169
170 void omxFreeData(omxData* od) {
171         omxFreeAllMatrixData(od->dataMat);
172         omxFreeAllMatrixData(od->meansMat);
173         Free(od);
174 }
175
176 double omxDoubleDataElement(omxData *od, int row, int col) {
177         if(od->dataMat != NULL) {
178                 return omxMatrixElement(od->dataMat, row, col);
179         }
180         int location = od->location[col];
181         if(location < 0) {
182                 return (double)(od->intData[~location][row]);
183         } else {
184                 return od->realData[location][row];
185         }
186 }
187
188 int omxIntDataElement(omxData *od, int row, int col) {
189         if(od->dataMat != NULL) {
190                 error("Use a data frame for factor data");
191         }
192
193         int location = od->location[col];
194         if(location < 0) {
195                 return (od->intData[~location][row]);
196         } else {
197                 return (int)(od->realData[location][row]);
198         }
199 }
200
201 omxMatrix* omxDataMatrix(omxData *od, omxMatrix* om) {
202         double dataElement;
203
204         if(od->dataMat != NULL) {               // Data was entered as a matrix.
205                 if(om != NULL) {                        // It stays as such
206                         omxCopyMatrix(om, od->dataMat);
207                         return om;
208                 }
209                 return od->dataMat;
210         }
211         // Otherwise, we must construct the matrix.
212         int numRows = od->rows, numCols = od->cols;
213
214         if(om == NULL) {
215                 om = omxInitMatrix(om, numRows, numCols, TRUE, od->currentState);
216         }
217
218         if(om->rows != numRows || om->cols != numCols) {
219                 omxResizeMatrix(om, numRows, numCols, FALSE);
220         }
221
222         for(int j = 0; j < numCols; j++) {
223                 for(int k = 0; k < numRows; k++) {
224                         int location = od->location[j];
225                         if(location < 0) {
226                                 dataElement = (double) od->intData[~location][k];
227                         } else {
228                                 dataElement = od->realData[location][k];
229                         }
230                         omxSetMatrixElement(om, k, j, dataElement);
231                 }
232         }
233         return om;
234 }
235
236 unsigned short int omxDataColumnIsFactor(omxData *od, int col) {
237         if(od->dataMat != NULL) return FALSE;
238         if(col <= od->cols) return (od->location[col] < 0);
239
240         error("Attempted to access column %d of a %d-column data object", col, od->cols);
241         return 0; // not reached
242 }
243
244 omxMatrix* omxDataMeans(omxData *od, omxMatrix* colList, omxMatrix* om) {
245
246         if(od->meansMat == NULL) return NULL;
247
248         if(colList == NULL) {
249                 if(om == NULL) return od->meansMat;
250                 omxCopyMatrix(om, od->meansMat);
251                 return om;
252         }
253
254         int cols = colList->cols;
255
256         if(colList == NULL || cols == 0 || cols > od->cols) {
257                 cols = od->cols;
258                 if(om == NULL) return od->meansMat;
259                 omxCopyMatrix(om, od->meansMat);
260                 return om;
261         }
262
263         if(om == NULL) {
264                 om = omxInitMatrix(om, 1, cols, TRUE, od->currentState);
265         }
266
267         for(int i = 0; i < cols; i++) {
268                 omxSetMatrixElement(om, 1, i, omxVectorElement(od->meansMat, omxVectorElement(colList, i)));
269         }
270
271         return om;
272 }
273
274
275 void omxSetContiguousDataColumns(omxContiguousData* contiguous, omxData* data, omxMatrix* colList) {
276
277         contiguous->isContiguous = FALSE;   // Assume not contiguous
278
279         if (data->dataMat == NULL) return; // Data has no matrix elements, so skip.
280
281         omxMatrix* dataMat = data->dataMat;
282         if (dataMat->colMajor) return;      // If data matrix is column-major, there's no continuity
283         
284         int colListLength = colList->cols;              // # of columns in the cov matrix
285         double start = omxVectorElement(colList, 0);    // Data col of first column of the covariance
286         contiguous->start = (int) start;                // That's our starting point.
287         contiguous->length = colListLength;             // And the length is ncol(cov)
288         
289         for(int i = 1; i < colListLength; i++) {        // Make sure that the col list is 
290                 double next = omxVectorElement(colList, i); // contiguously increasing in column number
291                 if (next != (start + i)) return;            // If it isn't, it's not contiguous data
292         }
293         
294         contiguous->isContiguous = TRUE;    // Passed.  This is contiguous.
295 }
296
297 void omxContiguousDataRow(omxData *od, int row, int start, int length, omxMatrix* om) {
298         // TODO: Might be better to combine this with omxDataRow to make a single accessor omxDataRow with a second signature that accepts an omxContiguousData argument.
299         if(row >= od->rows) error("Invalid row");
300
301         if(om == NULL) error("Must provide an output matrix");
302         
303         int numcols = od->cols;
304         omxMatrix* dataMat = od->dataMat;
305         double *dest = om->data;
306         double *source = dataMat->data + row * numcols + start;
307         memcpy(dest, source, sizeof(double) * length);
308 }
309
310 void omxDataRow(omxData *od, int row, omxMatrix* colList, omxMatrix* om) {
311
312         if(colList == NULL || row >= od->rows) error("Invalid row or colList");
313
314         if(om == NULL) error("Must provide an output matrix");
315
316         int numcols = om->cols;
317         if(od->dataMat != NULL) { // Matrix Object
318                 omxMatrix* dataMat = od->dataMat;
319                 for(int j = 0; j < numcols; j++) {
320                         omxSetMatrixElement(om, 0, j, omxMatrixElement(dataMat, row, 
321                                                                        omxVectorElement(colList, j)));
322                 }
323         } else {                // Data Frame object
324                 double dataElement;
325                 int* locations = od->location;
326                 int** intDataColumns = od->intData;
327                 double **realDataColumns = od->realData;
328                 for(int j = 0; j < numcols; j++) {
329                         int location = locations[(int)omxVectorElement(colList, j)];
330                         if(location < 0) {
331                                 dataElement = (double) intDataColumns[~location][row];
332                         } else {
333                                 dataElement = realDataColumns[location][row];
334                         }
335                         omxSetMatrixElement(om, 0, j, dataElement);
336                 }
337         }
338 }
339
340 int omxDataIndex(omxData *od, int row) {
341         if(od->indexVector != NULL)
342                 return od->indexVector[row];
343         else return row;
344 }
345
346 int omxDataNumIdenticalRows(omxData *od, int row) {
347         if(od->identicalRows != NULL)
348                 return od->identicalRows[row];
349         else return 1;
350 }
351 int omxDataNumIdenticalMissingness(omxData *od, int row) {
352         if(od->identicalMissingness != NULL)
353                 return od->identicalMissingness[row];
354         else return 1;
355 }
356
357 int omxDataNumIdenticalDefs(omxData *od, int row){
358         if(od->identicalDefs != NULL)
359                 return od->identicalDefs[row];
360         else return 1;
361 }
362
363 int omxDataNumIdenticalContinuousRows(omxData *od, int row) {
364         if(od->numNumeric <= 0) {
365                 return od->rows;
366         }
367         return omxDataNumIdenticalRows(od, row);
368 }
369
370 int omxDataNumIdenticalContinuousMissingness(omxData *od, int row) {
371         if(od->numNumeric <= 0) {
372                 return od->rows;
373         }
374         return omxDataNumIdenticalMissingness(od, row);
375 }
376
377 int omxDataNumIdenticalOrdinalRows(omxData *od, int row) {
378         if(od->numFactor <= 0) {
379                 return od->rows;
380         }
381         return omxDataNumIdenticalRows(od, row);
382 }
383
384 int omxDataNumIdenticalOrdinalMissingness(omxData *od, int row) {
385         if(od->numFactor <= 0) {
386                 return od->rows;
387         }
388         return omxDataNumIdenticalMissingness(od, row);
389 }
390
391
392 double omxDataNumObs(omxData *od) {
393         return od->numObs;
394 }
395
396 int omxDataNumFactor(omxData *od) {
397         return od->numFactor;
398 }
399
400 int omxDataNumNumeric(omxData *od) {
401         return od->numNumeric;
402 }
403
404
405 const char *omxDataType(omxData *od) {
406         return od->_type;
407 }
408
409 int elementEqualsDataframe(SEXP column, int offset1, int offset2) {
410         switch (TYPEOF(column)) {
411         case REALSXP:
412                 if(ISNA(REAL(column)[offset1])) return ISNA(REAL(column)[offset2]);
413                 if(ISNA(REAL(column)[offset2])) return ISNA(REAL(column)[offset1]);
414                 return(REAL(column)[offset1] == REAL(column)[offset2]);
415         case LGLSXP:
416         case INTSXP:
417                 return(INTEGER(column)[offset1] == INTEGER(column)[offset2]);           
418         }
419         return(0);
420 }
421
422 int testRowDataframe(SEXP data, int numrow, int numcol, int i, int *row, int base) {
423         SEXP column;
424         int j, equal = TRUE;
425
426         if (i == numrow) {
427                 equal = FALSE;
428         } else {
429                 for(j = 0; j < numcol && equal; j++) {
430                         column = VECTOR_ELT(data, j);
431                         equal = elementEqualsDataframe(column, base, i);
432                 }
433         }
434
435         if (!equal) {
436                 int gap = i - base;
437                 for(j = 0; j < gap; j++) {
438                         row[base + j] = gap - j;
439                 }
440                 base = i;
441         }
442         return(base);
443 }
444
445 int elementEqualsMatrix(SEXP data, int row1, int row2, int numrow, int col) {
446         int coloffset = col * numrow;
447         switch (TYPEOF(data)) {
448         case REALSXP:
449                 if(ISNA(REAL(data)[row1 + coloffset])) return ISNA(REAL(data)[row2 + coloffset]);
450                 if(ISNA(REAL(data)[row2 + coloffset])) return ISNA(REAL(data)[row1 + coloffset]);
451                 return(REAL(data)[row1 + coloffset] == REAL(data)[row2 + coloffset]);
452         case LGLSXP:
453         case INTSXP:
454                 return(INTEGER(data)[row1 + coloffset] == INTEGER(data)[row2 + coloffset]);
455         }
456         return(0);
457 }
458
459 int testRowMatrix(SEXP data, int numrow, int numcol, int i, int *row, int base) {
460         int j, equal = TRUE;
461
462         if (i == numrow) {
463                 equal = FALSE;
464         } else {
465                 for(j = 0; j < numcol && equal; j++) {
466                         equal = elementEqualsMatrix(data, i, base, numrow, j);
467                 }
468         }
469
470         if (!equal) {
471                 int gap = i - base;
472                 for(j = 0; j < gap; j++) {
473                         row[base + j] = gap - j;
474                 }
475                 base = i;
476         }
477         return(base);
478 }
479
480 SEXP findIdenticalMatrix(SEXP data, SEXP missing, SEXP defvars,
481                          SEXP skipMissingExp, SEXP skipDefvarsExp) {
482
483         SEXP retval, identicalRows, identicalMissing, identicalDefvars;
484         int i, numrow, numcol, defvarcol;
485         int *irows, *imissing, *idefvars;
486         int baserows, basemissing, basedefvars;
487         int skipMissing, skipDefvars;
488
489         skipMissing = LOGICAL(skipMissingExp)[0];
490         skipDefvars = LOGICAL(skipDefvarsExp)[0];
491         numrow = nrows(data);
492         numcol = ncols(data);
493         defvarcol = ncols(defvars);
494         PROTECT(retval = allocVector(VECSXP, 3));
495         PROTECT(identicalRows = allocVector(INTSXP, numrow));
496         PROTECT(identicalMissing = allocVector(INTSXP, numrow));
497         PROTECT(identicalDefvars = allocVector(INTSXP, numrow));
498         irows = INTEGER(identicalRows);
499         imissing = INTEGER(identicalMissing);
500         idefvars = INTEGER(identicalDefvars);
501         if (skipMissing) {
502                 for(i = 0; i < numrow; i++) {
503                         imissing[i] = numrow - i;
504                 }
505         }
506         if (skipDefvars) {
507                 for(i = 0; i < numrow; i++) {
508                         idefvars[i] = numrow - i;
509                 }
510         }
511         baserows = 0;
512         basemissing = 0;
513         basedefvars = 0;
514         for(i = 1; i <= numrow; i++) {
515                 baserows = testRowMatrix(data, numrow, numcol, i, irows, baserows); 
516                 if (!skipMissing) {
517                         basemissing = testRowMatrix(missing, numrow, numcol, i, imissing, basemissing); 
518                 }
519                 if (!skipDefvars) {
520                         basedefvars = testRowMatrix(defvars, numrow, defvarcol, i, idefvars, basedefvars);
521                 }
522         }
523         SET_VECTOR_ELT(retval, 0, identicalRows);
524         SET_VECTOR_ELT(retval, 1, identicalMissing);
525         SET_VECTOR_ELT(retval, 2, identicalDefvars);
526         UNPROTECT(4); // retval, identicalRows, identicalMissing, identicalDefvars
527         return retval;
528 }
529
530 SEXP findIdenticalDataFrame(SEXP data, SEXP missing, SEXP defvars,
531                             SEXP skipMissingExp, SEXP skipDefvarsExp) {
532
533         SEXP retval, identicalRows, identicalMissing, identicalDefvars;
534         int i, numrow, numcol, defvarcol;
535         int *irows, *imissing, *idefvars;
536         int baserows, basemissing, basedefvars;
537         int skipMissing, skipDefvars;
538
539         skipMissing = LOGICAL(skipMissingExp)[0];
540         skipDefvars = LOGICAL(skipDefvarsExp)[0];
541         numrow = length(VECTOR_ELT(data, 0));
542         numcol = length(data);
543         defvarcol = length(defvars);
544         PROTECT(retval = allocVector(VECSXP, 3));
545         PROTECT(identicalRows = allocVector(INTSXP, numrow));
546         PROTECT(identicalMissing = allocVector(INTSXP, numrow));
547         PROTECT(identicalDefvars = allocVector(INTSXP, numrow));
548         irows = INTEGER(identicalRows);
549         imissing = INTEGER(identicalMissing);
550         idefvars = INTEGER(identicalDefvars);
551         if (skipMissing) {
552                 for(i = 0; i < numrow; i++) {
553                         imissing[i] = numrow - i;
554                 }
555         }
556         if (skipDefvars) {
557                 for(i = 0; i < numrow; i++) {
558                         idefvars[i] = numrow - i;
559                 }
560         }
561         baserows = 0;
562         basemissing = 0;
563         basedefvars = 0;
564         for(i = 1; i <= numrow; i++) {
565                 baserows = testRowDataframe(data, numrow, numcol, i, irows, baserows); 
566                 if (!skipMissing) {
567                         basemissing = testRowMatrix(missing, numrow, numcol, i, imissing, basemissing);
568                 }
569                 if (!skipDefvars) {
570                         basedefvars = testRowDataframe(defvars, numrow, defvarcol, i, idefvars, basedefvars);
571                 }
572         }
573         SET_VECTOR_ELT(retval, 0, identicalRows);
574         SET_VECTOR_ELT(retval, 1, identicalMissing);
575         SET_VECTOR_ELT(retval, 2, identicalDefvars);
576         UNPROTECT(4); // retval, identicalRows, identicalMissing, identicalDefvars
577         return retval;
578 }
579
580 SEXP findIdenticalRowsData2(SEXP data, SEXP missing, SEXP defvars,
581                            SEXP skipMissingness, SEXP skipDefvars) {
582         if (isMatrix(data)) {
583                 return(findIdenticalMatrix(data, missing, defvars, skipMissingness, skipDefvars));
584         } else {
585                 return(findIdenticalDataFrame(data, missing, defvars, skipMissingness, skipDefvars));
586         }
587 }
588
589 SEXP findIdenticalRowsData(SEXP data, SEXP missing, SEXP defvars,
590                            SEXP skipMissingness, SEXP skipDefvars)
591 {
592         omxManageProtectInsanity protectManager;
593
594         try {
595                 return findIdenticalRowsData2(data, missing, defvars,
596                                               skipMissingness, skipDefvars);
597         } catch( std::exception& __ex__ ) {
598                 exception_to_try_error( __ex__ );
599         } catch(...) {
600                 string_to_try_error( "c++ exception (unknown reason)" );
601         }
602 }
603
604
605 void omxPrintData(omxData *od, const char *header) {
606         if (!header) error("omxPrintData: header is NULL");
607
608         if (!od) {
609                 Rprintf("%s: NULL\n", header);
610                 return;
611         }
612
613         Rprintf("%s(%s): %f observations %d x %d\n", header, od->_type, od->numObs,
614                 od->rows, od->cols);
615         Rprintf("numNumeric %d numFactor %d\n", od->numNumeric, od->numFactor);
616
617         if (od->location) {
618                 for(int j = 0; j < od->cols; j++) {
619                         int loc = od->location[j];
620                         if (loc < 0) {
621                                 Rprintf("Integer[%d]:", j);
622                                 int *val = od->intData[~loc];
623                                 for (int vx=0; vx < od->numObs; vx++) {
624                                         Rprintf(" %d", val[vx]);
625                                 }
626                         } else {
627                                 Rprintf("Numeric[%d]:", j);
628                                 double *val = od->realData[loc];
629                                 for (int vx=0; vx < od->numObs; vx++) {
630                                         Rprintf(" %.3g", val[vx]);
631                                 }
632                         }
633                         Rprintf("\n");
634                 }
635         }
636
637         if (od->identicalRows) {
638                 Rprintf("\trow\tmissing\tdefvars\n");
639                 for(int j = 0; j < od->rows; j++) {
640                         Rprintf("[%d]\t%d\t%d\t%d\n", j, od->identicalRows[j],
641                                 od->identicalMissingness[j], od->identicalDefs[j]);
642                 }
643         }
644
645         if (od->dataMat) omxPrintMatrix(od->dataMat, "dataMat");
646         if (od->meansMat) omxPrintMatrix(od->meansMat, "meansMat");
647 }
648