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