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