ifa: Use R_GetCCallable to access libirt-rpf
[openmx:openmx.git] / src / omxExpectationBA81.c
1 /*
2   Copyright 2012 Joshua Nathaniel Pritikin and contributors
3
4   This is free software: you can redistribute it and/or modify
5   it under the terms of the GNU General Public License as published by
6   the Free Software Foundation, either version 3 of the License, or
7   (at your option) any later version.
8
9   This program is distributed in the hope that it will be useful,
10   but WITHOUT ANY WARRANTY; without even the implied warranty of
11   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12   GNU General Public License for more details.
13
14   You should have received a copy of the GNU General Public License
15   along with this program.  If not, see <http://www.gnu.org/licenses/>.
16 */
17
18 // Consider replacing log() with log2() in some places? Not worth it?
19
20 #include "omxExpectation.h"
21 #include "omxOpenmpWrap.h"
22 #include "npsolWrap.h"
23 #include "libirt-rpf.h"
24 #include "merge.h"
25
26 static const char *NAME = "ExpectationBA81";
27
28 static const struct rpf *rpf_model = NULL;
29 static int rpf_numModels;
30
31 typedef double *(*rpf_fn_t)(omxExpectation *oo, omxMatrix *itemParam, const int *quad);
32
33 typedef struct {
34
35         // data characteristics
36         omxData *data;
37         int numUnique;
38         double *logNumIdentical;  // length numUnique
39         int *rowMap;              // length numUnique
40
41         // item description related
42         omxMatrix *itemSpec;
43         int maxOutcomes;
44         int maxDims;
45         int maxAbilities;
46         int numSpecific;
47         int *Sgroup;              // item's specific group 0..numSpecific-1
48         omxMatrix *design;        // items * maxDims
49
50         // quadrature related
51         int numQpoints;
52         double *Qpoint;
53         double *Qarea;
54         double *logQarea;
55         long *quadGridSize;       // maxDims
56         long totalPrimaryPoints;  // product of quadGridSize except specific dim
57         long totalQuadPoints;     // product of quadGridSize
58
59         // estimation related
60         omxMatrix *EitemParam;    // E step version
61         omxMatrix *itemParam;     // M step version
62         SEXP rpf;
63         rpf_fn_t computeRPF;
64         omxMatrix *customPrior;
65         int *paramMap;
66         int cacheLXK;             // w/cache,  numUnique * #specific quad points * totalQuadPoints
67         double *lxk;              // wo/cache, numUnique * thread
68         double *allSlxk;          // numUnique * thread
69         double *Slxk;             // numUnique * #specific dimensions * thread
70         double *patternLik;       // length numUnique
71         double ll;                // the most recent finite ll
72
73         int gradientCount;
74         int fitCount;
75 } omxBA81State;
76
77 /*
78 static void
79 pda(const double *ar, int rows, int cols) {
80         for (int rx=0; rx < rows; rx++) {
81                 for (int cx=0; cx < cols; cx++) {
82                         Rprintf("%.6g ", ar[cx * rows + rx]);
83                 }
84                 Rprintf("\n");
85         }
86
87 }
88 */
89
90 static int
91 findFreeVarLocation(omxMatrix *itemParam, const omxFreeVar *fv)
92 {
93         for (int lx=0; lx < fv->numLocations; lx++) {
94                 if (~fv->matrices[lx] == itemParam->matrixNumber) return lx;
95         }
96         return -1;
97 }
98
99 static int
100 compareFV(const int *fv1x, const int *fv2x, omxExpectation* oo)
101 {
102         omxState* currentState = oo->currentState;
103         omxBA81State *state = (omxBA81State *) oo->argStruct;
104         omxMatrix *itemParam = state->itemParam;
105         omxFreeVar *fv1 = currentState->freeVarList + *fv1x;
106         omxFreeVar *fv2 = currentState->freeVarList + *fv2x;
107         int l1 = findFreeVarLocation(itemParam, fv1);
108         int l2 = findFreeVarLocation(itemParam, fv2);
109         if (l1 == -1 && l2 == -1) return 0;
110         if ((l1 == -1) ^ (l2 == -1)) return l1 == -1? 1:-1;  // TODO reversed?
111         // Columns are items. Sort columns together
112         return fv1->col[l1] - fv2->col[l1];
113 }
114
115 static void buildParamMap(omxExpectation* oo)
116 {
117         omxState* currentState = oo->currentState;
118         omxBA81State *state = (omxBA81State *) oo->argStruct;
119         int numFreeParams = currentState->numFreeParams;
120         state->paramMap = Realloc(NULL, numFreeParams, int);
121         for (int px=0; px < numFreeParams; px++) { state->paramMap[px] = px; }
122         freebsd_mergesort(state->paramMap, numFreeParams, sizeof(int),
123                           (mergesort_cmp_t)compareFV, oo);
124 }
125
126 OMXINLINE static void
127 pointToWhere(omxBA81State *state, const int *quad, double *where, int upto)
128 {
129         for (int dx=0; dx < upto; dx++) {
130                 where[dx] = state->Qpoint[quad[dx]];
131         }
132 }
133
134 OMXINLINE static void
135 assignDims(omxMatrix *itemSpec, omxMatrix *design, int dims, int maxDims, int ix,
136            const double *restrict theta, double *restrict ptheta)
137 {
138         for (int dx=0; dx < dims; dx++) {
139                 int ability = (int)omxMatrixElement(design, dx, ix) - 1;
140                 if (ability >= maxDims) ability = maxDims-1;
141                 ptheta[dx] = theta[ability];
142         }
143 }
144
145 /**
146  * This is the main function needed to generate simulated data from
147  * the model. It could be argued that the rest of the estimation
148  * machinery belongs in the fitfunction.
149  *
150  * \param theta Vector of ability parameters, one per ability
151  * \returns A numItems by maxOutcomes colMajor vector of doubles. Caller must Free it.
152  */
153 static double *
154 standardComputeRPF(omxExpectation *oo, omxMatrix *itemParam, const int *quad)
155 {
156         omxBA81State *state = (omxBA81State*) oo->argStruct;
157         omxMatrix *itemSpec = state->itemSpec;
158         int numItems = itemSpec->cols;
159         omxMatrix *design = state->design;
160         int maxDims = state->maxDims;
161
162         double theta[maxDims];
163         pointToWhere(state, quad, theta, maxDims);
164
165         double *outcomeProb = Realloc(NULL, numItems * state->maxOutcomes, double);
166
167         for (int ix=0; ix < numItems; ix++) {
168                 const double *spec = omxMatrixColumn(itemSpec, ix);
169                 double *iparam = omxMatrixColumn(itemParam, ix);
170                 double *out = outcomeProb + ix * state->maxOutcomes;
171                 int id = spec[RPF_ISpecID];
172                 int dims = spec[RPF_ISpecDims];
173                 double ptheta[dims];
174                 assignDims(itemSpec, design, dims, maxDims, ix, theta, ptheta);
175                 (*rpf_model[id].logprob)(spec, iparam, ptheta, out);
176         }
177
178         return outcomeProb;
179 }
180
181 static double *
182 RComputeRPF1(omxExpectation *oo, omxMatrix *itemParam, const int *quad)
183 {
184         omxBA81State *state = (omxBA81State*) oo->argStruct;
185         int maxOutcomes = state->maxOutcomes;
186         omxMatrix *design = state->design;
187         omxMatrix *itemSpec = state->itemSpec;
188         int maxDims = state->maxDims;
189
190         double theta[maxDims];
191         pointToWhere(state, quad, theta, maxDims);
192
193         SEXP invoke;
194         PROTECT(invoke = allocVector(LANGSXP, 4));
195         SETCAR(invoke, state->rpf);
196         SETCADR(invoke, omxExportMatrix(itemParam));
197         SETCADDR(invoke, omxExportMatrix(itemSpec));
198
199         SEXP where;
200         PROTECT(where = allocMatrix(REALSXP, maxDims, itemParam->cols));
201         double *ptheta = REAL(where);
202         for (int ix=0; ix < itemParam->cols; ix++) {
203                 int dims = omxMatrixElement(itemSpec, RPF_ISpecDims, ix);
204                 assignDims(itemSpec, design, dims, maxDims, ix, theta, ptheta + ix*maxDims);
205                 for (int dx=dims; dx < maxDims; dx++) {
206                         ptheta[ix*maxDims + dx] = NA_REAL;
207                 }
208         }
209         SETCADDDR(invoke, where);
210
211         SEXP matrix;
212         PROTECT(matrix = eval(invoke, R_GlobalEnv));
213
214         if (!isMatrix(matrix)) {
215                 omxRaiseError(oo->currentState, -1,
216                               "RPF must return an item by outcome matrix");
217                 return NULL;
218         }
219
220         SEXP matrixDims;
221         PROTECT(matrixDims = getAttrib(matrix, R_DimSymbol));
222         int *dimList = INTEGER(matrixDims);
223         int numItems = state->itemSpec->cols;
224         if (dimList[0] != maxOutcomes || dimList[1] != numItems) {
225                 const int errlen = 200;
226                 char errstr[errlen];
227                 snprintf(errstr, errlen, "RPF must return a %d outcomes by %d items matrix",
228                          maxOutcomes, numItems);
229                 omxRaiseError(oo->currentState, -1, errstr);
230                 return NULL;
231         }
232
233         // Unlikely to be of type INTSXP, but just to be safe
234         PROTECT(matrix = coerceVector(matrix, REALSXP));
235         double *restrict got = REAL(matrix);
236
237         // Need to copy because threads cannot share SEXP
238         double *restrict outcomeProb = Realloc(NULL, numItems * maxOutcomes, double);
239
240         // Double check there aren't NAs in the wrong place
241         for (int ix=0; ix < numItems; ix++) {
242                 int numOutcomes = omxMatrixElement(state->itemSpec, RPF_ISpecOutcomes, ix);
243                 for (int ox=0; ox < numOutcomes; ox++) {
244                         int vx = ix * maxOutcomes + ox;
245                         if (isnan(got[vx])) {
246                                 const int errlen = 200;
247                                 char errstr[errlen];
248                                 snprintf(errstr, errlen, "RPF returned NA in [%d,%d]", ox,ix);
249                                 omxRaiseError(oo->currentState, -1, errstr);
250                         }
251                         outcomeProb[vx] = got[vx];
252                 }
253         }
254
255         return outcomeProb;
256 }
257
258 static double *
259 RComputeRPF(omxExpectation *oo, omxMatrix *itemParam, const int *quad)
260 {
261         omx_omp_set_lock(&GlobalRLock);
262         PROTECT_INDEX pi = omxProtectSave();
263         double *ret = RComputeRPF1(oo, itemParam, quad);
264         omxProtectRestore(pi);
265         omx_omp_unset_lock(&GlobalRLock);  // hope there was no exception!
266         return ret;
267 }
268
269 OMXINLINE static long
270 encodeLocation(const int dims, const long *restrict grid, const int *restrict quad)
271 {
272         long qx = 0;
273         for (int dx=dims-1; dx >= 0; dx--) {
274                 qx = qx * grid[dx];
275                 qx += quad[dx];
276         }
277         return qx;
278 }
279
280 #define CALC_LXK_CACHED(state, numUnique, quad, tqp, specific) \
281         ((state)->lxk + \
282          (numUnique) * encodeLocation((state)->maxDims, (state)->quadGridSize, quad) + \
283          (numUnique) * (tqp) * (specific))
284
285 OMXINLINE static double *
286 ba81Likelihood(omxExpectation *oo, int specific, const int *restrict quad)
287 {
288         omxBA81State *state = (omxBA81State*) oo->argStruct;
289         int numUnique = state->numUnique;
290         int maxOutcomes = state->maxOutcomes;
291         omxData *data = state->data;
292         int numItems = state->itemSpec->cols;
293         rpf_fn_t rpf_fn = state->computeRPF;
294         int *restrict Sgroup = state->Sgroup;
295         double *restrict lxk;
296
297         if (!state->cacheLXK) {
298                 lxk = state->lxk + numUnique * omx_absolute_thread_num();
299         } else {
300                 lxk = CALC_LXK_CACHED(state, numUnique, quad, state->totalQuadPoints, specific);
301         }
302
303         const double *outcomeProb = (*rpf_fn)(oo, state->EitemParam, quad);
304         if (!outcomeProb) {
305                 OMXZERO(lxk, numUnique);
306                 return lxk;
307         }
308
309         const int *rowMap = state->rowMap;
310         for (int px=0; px < numUnique; px++) {
311                 double lxk1 = 0;
312                 for (int ix=0; ix < numItems; ix++) {
313                         if (specific != Sgroup[ix]) continue;
314                         int pick = omxIntDataElementUnsafe(data, rowMap[px], ix);
315                         if (pick == NA_INTEGER) continue;
316                         lxk1 += outcomeProb[ix * maxOutcomes + pick-1];
317                 }
318                 lxk[px] = lxk1;
319         }
320
321         Free(outcomeProb);
322
323         return lxk;
324 }
325
326 OMXINLINE static double *
327 ba81LikelihoodFast(omxExpectation *oo, int specific, const int *restrict quad)
328 {
329         omxBA81State *state = (omxBA81State*) oo->argStruct;
330         if (!state->cacheLXK) {
331                 return ba81LikelihoodFast(oo, specific, quad);
332         } else {
333                 return CALC_LXK_CACHED(state, state->numUnique, quad, state->totalQuadPoints, specific);
334         }
335
336 }
337
338 #define CALC_ALLSLXK(state, numUnique) \
339         (state->allSlxk + omx_absolute_thread_num() * (numUnique))
340
341 #define CALC_SLXK(state, numUnique, numSpecific) \
342         (state->Slxk + omx_absolute_thread_num() * (numUnique) * (numSpecific))
343
344 OMXINLINE static void
345 cai2010(omxExpectation* oo, int recompute, const int *restrict primaryQuad,
346         double *restrict allSlxk, double *restrict Slxk)
347 {
348         omxBA81State *state = (omxBA81State*) oo->argStruct;
349         int numUnique = state->numUnique;
350         int numSpecific = state->numSpecific;
351         int maxDims = state->maxDims;
352         int sDim = maxDims-1;
353
354         int quad[maxDims];
355         memcpy(quad, primaryQuad, sizeof(int)*sDim);
356
357         OMXZERO(Slxk, numUnique * numSpecific);
358         OMXZERO(allSlxk, numUnique);
359
360         for (int sx=0; sx < numSpecific; sx++) {
361                 double *eis = Slxk + numUnique * sx;
362                 int quadGridSize = state->quadGridSize[sDim];
363
364                 for (int qx=0; qx < quadGridSize; qx++) {
365                         quad[sDim] = qx;
366                         double *lxk;
367                         if (recompute) {
368                                 lxk = ba81Likelihood(oo, sx, quad);
369                         } else {
370                                 lxk = CALC_LXK_CACHED(state, numUnique, quad, state->totalQuadPoints, sx);
371                         }
372
373                         for (int ix=0; ix < numUnique; ix++) {
374                                 eis[ix] += exp(lxk[ix] + state->logQarea[qx]);
375                         }
376                 }
377
378                 for (int px=0; px < numUnique; px++) {
379                         eis[px] = log(eis[px]);
380                         allSlxk[px] += eis[px];
381                 }
382         }
383 }
384
385 OMXINLINE static double
386 logAreaProduct(omxBA81State *state, const int *restrict quad, const int upto)
387 {
388         double logArea = 0;
389         for (int dx=0; dx < upto; dx++) {
390                 logArea += state->logQarea[quad[dx]];
391         }
392         return logArea;
393 }
394
395 // The idea of this API is to allow passing in a number larger than 1.
396 OMXINLINE static void
397 areaProduct(omxBA81State *state, const int *restrict quad, const int upto, double *restrict out)
398 {
399         for (int dx=0; dx < upto; dx++) {
400                 *out *= state->Qarea[quad[dx]];
401         }
402 }
403
404 OMXINLINE static void
405 decodeLocation(long qx, const int dims, const long *restrict grid,
406                int *restrict quad)
407 {
408         for (int dx=0; dx < dims; dx++) {
409                 quad[dx] = qx % grid[dx];
410                 qx = qx / grid[dx];
411         }
412 }
413
414 static void
415 ba81Estep(omxExpectation *oo) {
416         if(OMX_DEBUG_MML) {Rprintf("Beginning %s Computation.\n", NAME);}
417
418         omxBA81State *state = (omxBA81State*) oo->argStruct;
419         double *patternLik = state->patternLik;
420         int numUnique = state->numUnique;
421         int numSpecific = state->numSpecific;
422
423         omxCopyMatrix(state->EitemParam, state->itemParam);
424
425         OMXZERO(patternLik, numUnique);
426
427         // E-step, marginalize person ability
428         //
429         // Note: In the notation of Bock & Aitkin (1981) and
430         // Cai~(2010), these loops are reversed.  That is, the inner
431         // loop is over quadrature points and the outer loop is over
432         // all response patterns.
433         //
434         if (numSpecific == 0) {
435 #pragma omp parallel for num_threads(oo->currentState->numThreads)
436                 for (long qx=0; qx < state->totalQuadPoints; qx++) {
437                         int quad[state->maxDims];
438                         decodeLocation(qx, state->maxDims, state->quadGridSize, quad);
439
440                         double *lxk = ba81Likelihood(oo, 0, quad);
441
442                         double logArea = logAreaProduct(state, quad, state->maxDims);
443 #pragma omp critical(EstepUpdate)
444                         for (int px=0; px < numUnique; px++) {
445                                 double tmp = exp(lxk[px] + logArea);
446                                 patternLik[px] += tmp;
447                         }
448                 }
449         } else {
450                 int sDim = state->maxDims-1;
451
452 #pragma omp parallel for num_threads(oo->currentState->numThreads)
453                 for (long qx=0; qx < state->totalPrimaryPoints; qx++) {
454                         int quad[state->maxDims];
455                         decodeLocation(qx, sDim, state->quadGridSize, quad);
456
457                         double *allSlxk = CALC_ALLSLXK(state, numUnique);
458                         double *Slxk = CALC_SLXK(state, numUnique, numSpecific);
459                         cai2010(oo, TRUE, quad, allSlxk, Slxk);
460
461                         double logArea = logAreaProduct(state, quad, sDim);
462 #pragma omp critical(EstepUpdate)
463                         for (int px=0; px < numUnique; px++) {
464                                 double tmp = exp(allSlxk[px] + logArea);
465                                 patternLik[px] += tmp;
466                         }
467                 }
468         }
469
470         for (int px=0; px < numUnique; px++) {
471                 patternLik[px] = log(patternLik[px]);
472         }
473 }
474
475 OMXINLINE static void
476 expectedUpdate(omxData *restrict data, const int *rowMap, const int px, const int item,
477                const double observed, const int outcomes, double *out)
478 {
479         int pick = omxIntDataElementUnsafe(data, rowMap[px], item);
480         if (pick == NA_INTEGER) {
481                 double slice = exp(observed - log(outcomes));
482                 for (int ox=0; ox < outcomes; ox++) {
483                         out[ox] += slice;
484                 }
485         } else {
486                 out[pick-1] += exp(observed);
487         }
488 }
489
490 /** 
491  * \param quad a vector that indexes into a multidimensional quadrature
492  * \param out points to an array numOutcomes wide
493  */
494 OMXINLINE static void
495 ba81Weight(omxExpectation* oo, const int item, const int *quad, int outcomes, double *out)
496 {
497         omxBA81State *state = (omxBA81State*) oo->argStruct;
498         omxData *data = state->data;
499         const int *rowMap = state->rowMap;
500         int specific = state->Sgroup[item];
501         double *patternLik = state->patternLik;
502         double *logNumIdentical = state->logNumIdentical;
503         int numUnique = state->numUnique;
504         int numSpecific = state->numSpecific;
505         int sDim = state->maxDims-1;
506
507         OMXZERO(out, outcomes);
508
509         if (numSpecific == 0) {
510                 double *lxk = ba81LikelihoodFast(oo, specific, quad);
511                 for (int px=0; px < numUnique; px++) {
512                         double observed = logNumIdentical[px] + lxk[px] - patternLik[px];
513                         expectedUpdate(data, rowMap, px, item, observed, outcomes, out);
514                 }
515         } else {
516                 double *allSlxk = CALC_ALLSLXK(state, numUnique);
517                 double *Slxk = CALC_SLXK(state, numUnique, numSpecific);
518                 if (quad[sDim] == 0) {
519                         // allSlxk, Slxk only depend on the ordinate of the primary dimensions
520                         cai2010(oo, !state->cacheLXK, quad, allSlxk, Slxk);
521                 }
522                 double *eis = Slxk + numUnique * specific;
523
524                 // Avoid recalc when cache disabled with modest buffer? TODO
525                 double *lxk = ba81LikelihoodFast(oo, specific, quad);
526
527                 for (int px=0; px < numUnique; px++) {
528                         double observed = logNumIdentical[px] + (allSlxk[px] - eis[px]) +
529                                 (lxk[px] - patternLik[px]);
530                         expectedUpdate(data, rowMap, px, item, observed, outcomes, out);
531                 }
532         }
533 }
534
535 OMXINLINE static double
536 ba81Fit1Ordinate(omxExpectation* oo, const int *quad)
537 {
538         omxBA81State *state = (omxBA81State*) oo->argStruct;
539         omxMatrix *itemParam = state->itemParam;
540         int numItems = itemParam->cols;
541         rpf_fn_t rpf_fn = state->computeRPF;
542         int maxOutcomes = state->maxOutcomes;
543         int maxDims = state->maxDims;
544
545         double *outcomeProb = (*rpf_fn)(oo, itemParam, quad);
546         if (!outcomeProb) return 0;
547
548         double thr_ll = 0;
549         for (int ix=0; ix < numItems; ix++) {
550                 int outcomes = omxMatrixElement(state->itemSpec, RPF_ISpecOutcomes, ix);
551                 double out[outcomes];
552                 ba81Weight(oo, ix, quad, outcomes, out);
553                 for (int ox=0; ox < outcomes; ox++) {
554                         double got = out[ox] * outcomeProb[ix * maxOutcomes + ox];
555                         areaProduct(state, quad, maxDims, &got);
556                         thr_ll += got;
557                 }
558         }
559
560         Free(outcomeProb);
561         return thr_ll;
562 }
563
564 static double
565 ba81ComputeFit1(omxExpectation* oo)
566 {
567         omxBA81State *state = (omxBA81State*) oo->argStruct;
568         ++state->fitCount;
569         omxMatrix *customPrior = state->customPrior;
570         int numSpecific = state->numSpecific;
571         int maxDims = state->maxDims;
572
573         double ll = 0;
574         if (customPrior) {
575                 omxRecompute(customPrior);
576                 ll = customPrior->data[0];
577         } else {
578                 omxMatrix *itemSpec = state->itemSpec;
579                 omxMatrix *itemParam = state->itemParam;
580                 int numItems = itemSpec->cols;
581                 for (int ix=0; ix < numItems; ix++) {
582                         const double *spec = omxMatrixColumn(itemSpec, ix);
583                         int id = spec[RPF_ISpecID];
584                         double *iparam = omxMatrixColumn(itemParam, ix);
585                         ll += (*rpf_model[id].prior)(spec, iparam);
586                 }
587         }
588
589         if (numSpecific == 0) {
590 #pragma omp parallel for num_threads(oo->currentState->numThreads)
591                 for (long qx=0; qx < state->totalQuadPoints; qx++) {
592                         int quad[maxDims];
593                         decodeLocation(qx, maxDims, state->quadGridSize, quad);
594                         double thr_ll = ba81Fit1Ordinate(oo, quad);
595
596 #pragma omp atomic
597                         ll += thr_ll;
598                 }
599         } else {
600                 int sDim = state->maxDims-1;
601                 long *quadGridSize = state->quadGridSize;
602
603 #pragma omp parallel for num_threads(oo->currentState->numThreads)
604                 for (long qx=0; qx < state->totalPrimaryPoints; qx++) {
605                         int quad[maxDims];
606                         decodeLocation(qx, maxDims, quadGridSize, quad);
607
608                         double thr_ll = 0;
609                         long specificPoints = quadGridSize[sDim];
610                         for (long sx=0; sx < specificPoints; sx++) {
611                                 quad[sDim] = sx;
612                                 thr_ll += ba81Fit1Ordinate(oo, quad);
613                         }
614 #pragma omp atomic
615                         ll += thr_ll;
616                 }
617         }
618
619         if (isinf(ll)) {
620                 return 2*state->ll;
621         } else {
622                 ll = -2 * ll;
623                 state->ll = ll;
624                 return ll;
625         }
626 }
627
628 double
629 ba81ComputeFit(omxExpectation* oo)
630 {
631         double got = ba81ComputeFit1(oo);
632         return got;
633 }
634
635 OMXINLINE static void
636 ba81ItemGradientOrdinate(omxExpectation* oo, omxBA81State *state,
637                          int maxDims, int *quad, int item, const double *spec, int id,
638                          int outcomes, double *iparam, int numParam, int *paramMask, double *gq)
639 {
640         double where[maxDims];
641         pointToWhere(state, quad, where, maxDims);
642         double weight[outcomes];
643         ba81Weight(oo, item, quad, outcomes, weight);
644
645         (*rpf_model[id].gradient)(spec, iparam, paramMask, where, weight, gq);
646
647         for (int ox=0; ox < numParam; ox++) {
648                 if (paramMask[ox] == -1) continue;
649                 areaProduct(state, quad, maxDims, gq+ox);
650         }
651 }
652
653 OMXINLINE static void
654 ba81ItemGradient(omxExpectation* oo, omxBA81State *state, const double *spec, omxMatrix *itemParam,
655                  int item, int id, int outcomes, int numParam, int *paramMask, double *out)
656 {
657         int maxDims = state->maxDims;
658         double *iparam = omxMatrixColumn(itemParam, item);
659         double gradient[numParam];
660         OMXZERO(gradient, numParam);
661
662         if (state->numSpecific == 0) {
663 #pragma omp parallel for num_threads(oo->currentState->numThreads)
664                 for (long qx=0; qx < state->totalQuadPoints; qx++) {
665                         int quad[maxDims];
666                         decodeLocation(qx, maxDims, state->quadGridSize, quad);
667                         double gq[numParam];
668                         //                      for (int gx=0; gx < numParam; gx++) gq[gx] = 3.14159; // debugging TODO
669
670                         ba81ItemGradientOrdinate(oo, state, maxDims, quad, item, spec, id,
671                                                  outcomes, iparam, numParam, paramMask, gq);
672
673 #pragma omp critical(GradientUpdate)
674                         for (int ox=0; ox < numParam; ox++) {
675                                 gradient[ox] += gq[ox];
676                         }
677                 }
678         } else {
679                 int sDim = state->maxDims-1;
680                 long *quadGridSize = state->quadGridSize;
681 #pragma omp parallel for num_threads(oo->currentState->numThreads)
682                 for (long qx=0; qx < state->totalPrimaryPoints; qx++) {
683                         int quad[maxDims];
684                         decodeLocation(qx, maxDims, quadGridSize, quad);
685                         double gsubtotal[numParam];
686                         OMXZERO(gsubtotal, numParam);
687
688                         long specificPoints = quadGridSize[sDim];
689                         for (long sx=0; sx < specificPoints; sx++) {
690                                 double gq[numParam];
691                                 quad[sDim] = sx;
692                                 ba81ItemGradientOrdinate(oo, state, maxDims, quad, item, spec, id,
693                                                          outcomes, iparam, numParam, paramMask, gq);
694                                 for (int gx=0; gx < numParam; gx++) {
695                                         gsubtotal[gx] += gq[gx];
696                                 }
697                         }
698 #pragma omp critical(GradientUpdate)
699                         for (int gx=0; gx < numParam; gx++) {
700                                 gradient[gx] += gsubtotal[gx];
701                         }
702                 }
703         }
704
705         (*rpf_model[id].gradient)(spec, iparam, paramMask, NULL, NULL, gradient);
706
707         for (int px=0; px < numParam; px++) {
708                 int loc = paramMask[px];
709                 if (loc == -1) continue;
710                 out[loc] = -2 * gradient[px];
711         }
712 }
713
714 void ba81Gradient(omxExpectation* oo, double *out)
715 {
716         omxState* currentState = oo->currentState;
717         int numFreeParams = currentState->numFreeParams;
718         omxBA81State *state = (omxBA81State *) oo->argStruct;
719         if (!state->paramMap) buildParamMap(oo);
720         ++state->gradientCount;
721         omxMatrix *itemSpec = state->itemSpec;
722         omxMatrix *itemParam = state->itemParam;
723
724         int vx = 0;
725         while (vx < numFreeParams) {
726             omxFreeVar *fv = currentState->freeVarList + state->paramMap[vx];
727             int vloc = findFreeVarLocation(itemParam, fv);
728             if (vloc < 0) {
729                     ++vx;
730                     continue;
731             }
732
733             int item = fv->col[vloc];
734             const double *spec = omxMatrixColumn(itemSpec, item);
735             int id = spec[RPF_ISpecID];
736             int outcomes = spec[RPF_ISpecOutcomes];
737             int numParam = (*rpf_model[id].numParam)(spec);
738
739             int paramMask[numParam];
740             for (int px=0; px < numParam; px++) { paramMask[px] = -1; }
741
742             if (fv->row[vloc] >= numParam) {
743                     warning("Item %d has too many free parameters", item);
744                     continue;
745             }
746             paramMask[fv->row[vloc]] = vx;
747
748             while (++vx < numFreeParams) {
749                     omxFreeVar *fv = currentState->freeVarList + state->paramMap[vx];
750                     int vloc = findFreeVarLocation(itemParam, fv);
751                     if (fv->col[vloc] != item) break;
752                     if (fv->row[vloc] >= numParam) {
753                             warning("Item %d has too many free parameters", item);
754                             continue;
755                     }
756                     paramMask[fv->row[vloc]] = vx;
757             }
758
759             ba81ItemGradient(oo, state, spec, itemParam, item,
760                              id, outcomes, numParam, paramMask, out);
761         }
762 }
763
764 static int
765 getNumThreads(omxExpectation* oo)
766 {
767         int numThreads = oo->currentState->numThreads;
768         if (numThreads < 1) numThreads = 1;
769         return numThreads;
770 }
771
772 static void
773 ba81SetupQuadrature(omxExpectation* oo, int numPoints, double *points, double *area)
774 {
775         omxBA81State *state = (omxBA81State *) oo->argStruct;
776         int numUnique = state->numUnique;
777         int numThreads = getNumThreads(oo);
778
779         state->numQpoints = numPoints;
780
781         Free(state->Qpoint);
782         Free(state->Qarea);
783         state->Qpoint = Realloc(NULL, numPoints, double);
784         state->Qarea = Realloc(NULL, numPoints, double);
785         memcpy(state->Qpoint, points, sizeof(double)*numPoints);
786         memcpy(state->Qarea, area, sizeof(double)*numPoints);
787
788         Free(state->logQarea);
789
790         state->logQarea = Realloc(NULL, state->numQpoints, double);
791         for (int px=0; px < state->numQpoints; px++) {
792                 state->logQarea[px] = log(state->Qarea[px]);
793         }
794
795         state->totalQuadPoints = 1;
796         state->totalPrimaryPoints = 1;
797         state->quadGridSize = (long*) R_alloc(state->maxDims, sizeof(long));
798         for (int dx=0; dx < state->maxDims; dx++) {
799                 state->quadGridSize[dx] = state->numQpoints;
800                 state->totalQuadPoints *= state->quadGridSize[dx];
801                 if (dx < state->maxDims-1) {
802                         state->totalPrimaryPoints *= state->quadGridSize[dx];
803                 }
804         }
805
806         Free(state->lxk);
807
808         if (!state->cacheLXK) {
809                 state->lxk = Realloc(NULL, numUnique * numThreads, double);
810         } else {
811                 int ns = state->numSpecific;
812                 if (ns == 0) ns = 1;
813                 state->lxk = Realloc(NULL, numUnique * state->totalQuadPoints * ns, double);
814         }
815 }
816
817 static void
818 ba81EAP1(omxExpectation *oo, double *workspace, long qx, int maxDims, int numUnique,
819          double *ability, double *cov, double *spstats)
820 {
821         omxBA81State *state = (omxBA81State *) oo->argStruct;
822         double *patternLik = state->patternLik;
823         int quad[maxDims];
824         decodeLocation(qx, maxDims, state->quadGridSize, quad);
825         double where[maxDims];
826         pointToWhere(state, quad, where, maxDims);
827         double logArea = logAreaProduct(state, quad, maxDims);
828         double *lxk = ba81LikelihoodFast(oo, 0, quad);
829         double *myspace = workspace + 2 * maxDims * numUnique * omx_absolute_thread_num();
830
831         for (int px=0; px < numUnique; px++) {
832                 double *piece = myspace + px * 2 * maxDims;
833                 double plik = exp(lxk[px] + logArea - patternLik[px]);
834                 for (int dx=0; dx < maxDims; dx++) {
835                         piece[dx] = where[dx] * plik;
836                 }
837                 /*
838                 for (int d1=0; d1 < maxDims; d1++) {
839                         for (int d2=0; d2 <= d1; d2++) {
840                                 covPiece[d1 * maxDims + d2] += piece[d1] * piece[d2];
841                         }
842                 }
843                 */
844         }
845 #pragma omp critical(EAP1Update)
846         for (int px=0; px < numUnique; px++) {
847                 double *piece = myspace + px * 2 * maxDims;
848                 double *arow = ability + px * 2 * maxDims;
849                 for (int dx=0; dx < maxDims; dx++) {
850                         arow[dx*2] += piece[dx];
851                 }
852                 /*
853                 for (int d1=0; d1 < maxDims; d1++) {
854                         for (int d2=0; d2 <= d1; d2++) {
855                                 int loc = d1 * maxDims + d2;
856                                 cov[loc] += covPiece[loc];
857                         }
858                 }
859                 */
860         }
861 }
862
863 static void
864 ba81EAP2(omxExpectation *oo, double *workspace, long qx, int maxDims, int numUnique,
865          double *ability, double *spstats)
866 {
867         omxBA81State *state = (omxBA81State *) oo->argStruct;
868         double *patternLik = state->patternLik;
869         int quad[maxDims];
870         decodeLocation(qx, maxDims, state->quadGridSize, quad);
871         double where[maxDims];
872         pointToWhere(state, quad, where, maxDims);
873         double logArea = logAreaProduct(state, quad, maxDims);
874         double *lxk = ba81LikelihoodFast(oo, 0, quad);
875
876         for (int px=0; px < numUnique; px++) {
877                 double psd[maxDims];
878                 double *arow = ability + px * 2 * maxDims;
879                 for (int dx=0; dx < maxDims; dx++) {
880                         // is this just sqrt(variance) and redundant with the covariance matrix? TODO
881                         double ldiff = log(fabs(where[dx] - arow[dx*2]));
882                         psd[dx] = exp(2 * ldiff + lxk[px] + logArea - patternLik[px]);
883                 }
884 #pragma omp critical(EAP1Update)
885                 for (int dx=0; dx < maxDims; dx++) {
886                         arow[dx*2+1] += psd[dx];
887                 }
888         }
889 }
890
891 omxRListElement *
892 ba81EAP(omxExpectation *oo, int *numReturns)
893 {
894         omxBA81State *state = (omxBA81State *) oo->argStruct;
895         int maxDims = state->maxDims;
896         //int numSpecific = state->numSpecific;
897
898         *numReturns = 2; // + (maxDims > 1) + (numSpecific > 1);
899         omxRListElement *out = (omxRListElement*) R_alloc(*numReturns, sizeof(omxRListElement));
900
901         out[0].numValues = 1;
902         out[0].values = (double*) R_alloc(1, sizeof(double));
903         strcpy(out[0].label, "Minus2LogLikelihood");
904         out[0].values[0] = state->ll;
905
906         omxData *data = state->data;
907         int numUnique = state->numUnique;
908
909         // TODO Wainer & Thissen. (1987). Estimating ability with the wrong
910         // model. Journal of Educational Statistics, 12, 339-368.
911
912         int numQpoints = state->numQpoints * 2;  // make configurable TODO
913         double Qpoint[numQpoints];
914         double Qarea[numQpoints];
915         const double Qwidth = 4;
916         for (int qx=0; qx < numQpoints; qx++) {
917                 Qpoint[qx] = Qwidth - qx * Qwidth*2 / (numQpoints-1);
918                 Qarea[qx] = 1.0/numQpoints;
919         }
920         ba81SetupQuadrature(oo, numQpoints, Qpoint, Qarea);
921         ba81Estep(oo);   // recalc patternLik with a flat prior
922
923         double *cov = NULL;
924         /*
925         if (maxDims > 1) {
926                 strcpy(out[2].label, "ability.cov");
927                 out[2].numValues = -1;
928                 out[2].rows = maxDims;
929                 out[2].cols = maxDims;
930                 out[2].values = (double*) R_alloc(out[2].rows * out[2].cols, sizeof(double));
931                 cov = out[2].values;
932                 OMXZERO(cov, out[2].rows * out[2].cols);
933         }
934         */
935         double *spstats = NULL;
936         /*
937         if (numSpecific) {
938                 strcpy(out[3].label, "specific");
939                 out[3].numValues = -1;
940                 out[3].rows = numSpecific;
941                 out[3].cols = 2;
942                 out[3].values = (double*) R_alloc(out[3].rows * out[3].cols, sizeof(double));
943                 spstats = out[3].values;
944         }
945         */
946
947         // allocation of workspace could be optional
948         int numThreads = getNumThreads(oo);
949         double *workspace = Realloc(NULL, numUnique * maxDims * 2 * numThreads, double);
950
951         // Need a separate work space because the destination needs
952         // to be in unsorted order with duplicated rows.
953         double *ability = Calloc(numUnique * maxDims * 2, double);
954
955 #pragma omp parallel for num_threads(oo->currentState->numThreads)
956         for (long qx=0; qx < state->totalQuadPoints; qx++) {
957                 ba81EAP1(oo, workspace, qx, maxDims, numUnique, ability, cov, spstats);
958         }
959
960         /*
961         // make symmetric
962         for (int d1=0; d1 < maxDims; d1++) {
963                 for (int d2=0; d2 < d1; d2++) {
964                         cov[d2 * maxDims + d1] = cov[d1 * maxDims + d2];
965                 }
966         }
967         */
968
969 #pragma omp parallel for num_threads(oo->currentState->numThreads)
970         for (long qx=0; qx < state->totalQuadPoints; qx++) {
971                 ba81EAP2(oo, workspace, qx, maxDims, numUnique, ability, spstats);
972         }
973
974         for (int px=0; px < numUnique; px++) {
975                 double *arow = ability + px * 2 * maxDims;
976                 for (int dx=0; dx < maxDims; dx++) {
977                         arow[dx*2+1] = sqrt(arow[dx*2+1]);
978                 }
979         }
980
981         strcpy(out[1].label, "ability");
982         out[1].numValues = -1;
983         out[1].rows = data->rows;
984         out[1].cols = 2 * maxDims;
985         out[1].values = (double*) R_alloc(out[1].rows * out[1].cols, sizeof(double));
986
987         for (int rx=0; rx < numUnique; rx++) {
988                 double *pa = ability + rx * 2 * maxDims;
989
990                 int dups = omxDataNumIdenticalRows(state->data, state->rowMap[rx]);
991                 for (int dup=0; dup < dups; dup++) {
992                         int dest = omxDataIndex(data, state->rowMap[rx]+dup);
993                         int col=0;
994                         for (int dx=0; dx < maxDims; dx++) {
995                                 out[1].values[col * out[1].rows + dest] = pa[col]; ++col;
996                                 out[1].values[col * out[1].rows + dest] = pa[col]; ++col;
997                         }
998                 }
999         }
1000         Free(ability);
1001         Free(workspace);
1002         return out;
1003 }
1004
1005 static void ba81Destroy(omxExpectation *oo) {
1006         if(OMX_DEBUG) {
1007                 Rprintf("Freeing %s function.\n", NAME);
1008         }
1009         omxBA81State *state = (omxBA81State *) oo->argStruct;
1010         Rprintf("fit %d gradient %d\n", state->fitCount, state->gradientCount);
1011         omxFreeAllMatrixData(state->itemSpec);
1012         omxFreeAllMatrixData(state->itemParam);
1013         omxFreeAllMatrixData(state->EitemParam);
1014         omxFreeAllMatrixData(state->design);
1015         omxFreeAllMatrixData(state->customPrior);
1016         Free(state->logNumIdentical);
1017         Free(state->Qpoint);
1018         Free(state->Qarea);
1019         Free(state->logQarea);
1020         Free(state->rowMap);
1021         Free(state->patternLik);
1022         Free(state->lxk);
1023         Free(state->Slxk);
1024         Free(state->allSlxk);
1025         Free(state->Sgroup);
1026         Free(state->paramMap);
1027         Free(state);
1028 }
1029
1030 int ba81ExpectationHasGradients(omxExpectation* oo)
1031 {
1032         omxBA81State *state = (omxBA81State *) oo->argStruct;
1033         return state->computeRPF == standardComputeRPF;
1034 }
1035
1036 void omxInitExpectationBA81(omxExpectation* oo) {
1037         omxState* currentState = oo->currentState;      
1038         SEXP rObj = oo->rObj;
1039         SEXP tmp;
1040         
1041         if(OMX_DEBUG) {
1042                 Rprintf("Initializing %s.\n", NAME);
1043         }
1044         if (!rpf_model) {
1045                 get_librpf_t get_librpf = (get_librpf_t) R_GetCCallable("rpf", "get_librpf_model");
1046                 (*get_librpf)(&rpf_numModels, &rpf_model);
1047         }
1048         
1049         omxBA81State *state = Calloc(1, omxBA81State);
1050         oo->argStruct = (void*) state;
1051
1052         state->ll = 1e20;   // finite but big
1053         
1054         PROTECT(tmp = GET_SLOT(rObj, install("data")));
1055         state->data = omxNewDataFromMxDataPtr(tmp, currentState);
1056         UNPROTECT(1);
1057
1058         if (strcmp(omxDataType(state->data), "raw") != 0) {
1059                 omxRaiseErrorf(currentState, "%s unable to handle data type %s", NAME, omxDataType(state->data));
1060                 return;
1061         }
1062
1063         PROTECT(state->rpf = GET_SLOT(rObj, install("RPF")));
1064         if (state->rpf == R_NilValue) {
1065                 state->computeRPF = standardComputeRPF;
1066         } else {
1067                 state->computeRPF = RComputeRPF;
1068         }
1069
1070         state->itemSpec =
1071                 omxNewMatrixFromIndexSlot(rObj, currentState, "ItemSpec");
1072         state->design =
1073                 omxNewMatrixFromIndexSlot(rObj, currentState, "Design");
1074         state->itemParam =
1075                 omxNewMatrixFromIndexSlot(rObj, currentState, "ItemParam");
1076         state->EitemParam =
1077                 omxInitTemporaryMatrix(NULL, state->itemParam->rows, state->itemParam->cols,
1078                                        TRUE, currentState);
1079         state->customPrior =
1080                 omxNewMatrixFromIndexSlot(rObj, currentState, "CustomPrior");
1081         
1082         oo->computeFun = ba81Estep;
1083         oo->destructFun = ba81Destroy;
1084         
1085         // TODO: Exactly identical rows do not contribute any information.
1086         // The sorting algorithm ought to remove them so we don't waste RAM.
1087         // The following summary stats would be cheaper to calculate too.
1088
1089         int numUnique = 0;
1090         omxData *data = state->data;
1091         if (omxDataNumFactor(data) != data->cols) {
1092                 // verify they are ordered factors TODO
1093                 omxRaiseErrorf(currentState, "%s: all columns must be factors", NAME);
1094                 return;
1095         }
1096
1097         for (int rx=0; rx < data->rows;) {
1098                 rx += omxDataNumIdenticalRows(state->data, rx);
1099                 ++numUnique;
1100         }
1101         state->numUnique = numUnique;
1102
1103         state->rowMap = Realloc(NULL, numUnique, int);
1104         state->logNumIdentical = Realloc(NULL, numUnique, double);
1105
1106         int numItems = state->itemParam->cols;
1107
1108         for (int rx=0, ux=0; rx < data->rows; ux++) {
1109                 if (rx == 0) {
1110                         // all NA rows will sort to the top
1111                         int na=0;
1112                         for (int ix=0; ix < numItems; ix++) {
1113                                 if (omxIntDataElement(data, 0, ix) == NA_INTEGER) { ++na; }
1114                         }
1115                         if (na == numItems) {
1116                                 omxRaiseErrorf(currentState, "Remove rows with all NAs");
1117                                 return;
1118                         }
1119                 }
1120                 int dups = omxDataNumIdenticalRows(state->data, rx);
1121                 state->logNumIdentical[ux] = log(dups);
1122                 state->rowMap[ux] = rx;
1123                 rx += dups;
1124         }
1125
1126         state->patternLik = Realloc(NULL, numUnique, double);
1127
1128         int numThreads = getNumThreads(oo);
1129
1130         int maxSpec = 0;
1131         int maxParam = 0;
1132         state->maxDims = 0;
1133         state->maxOutcomes = 0;
1134
1135         for (int cx = 0; cx < data->cols; cx++) {
1136                 const double *spec = omxMatrixColumn(state->itemSpec, cx);
1137                 int id = spec[RPF_ISpecID];
1138                 if (id < 0 || id >= rpf_numModels) {
1139                         omxRaiseErrorf(currentState, "ItemSpec column %d has unknown item model %d", cx, id);
1140                         return;
1141                 }
1142
1143                 int dims = spec[RPF_ISpecDims];
1144                 if (state->maxDims < dims)
1145                         state->maxDims = dims;
1146
1147                 // TODO verify that item model can have requested number of outcomes
1148                 int no = spec[RPF_ISpecOutcomes];
1149                 if (state->maxOutcomes < no)
1150                         state->maxOutcomes = no;
1151
1152                 int numSpec = (*rpf_model[id].numSpec)(spec);
1153                 if (maxSpec < numSpec)
1154                         maxSpec = numSpec;
1155
1156                 int numParam = (*rpf_model[id].numParam)(spec);
1157                 if (maxParam < numParam)
1158                         maxParam = numParam;
1159         }
1160
1161         if (state->itemSpec->cols != data->cols || state->itemSpec->rows != maxSpec) {
1162                 omxRaiseErrorf(currentState, "ItemSpec must have %d item columns and %d rows",
1163                                data->cols, maxSpec);
1164                 return;
1165         }
1166         if (state->itemParam->rows != maxParam) {
1167                 omxRaiseErrorf(currentState, "ItemParam should have %d rows", maxParam);
1168                 return;
1169         }
1170
1171         if (state->design == NULL) {
1172                 state->maxAbilities = state->maxDims;
1173                 state->design = omxInitTemporaryMatrix(NULL, state->maxDims, numItems,
1174                                        TRUE, currentState);
1175                 for (int ix=0; ix < numItems; ix++) {
1176                         for (int dx=0; dx < state->maxDims; dx++) {
1177                                 omxSetMatrixElement(state->design, dx, ix, (double)dx+1);
1178                         }
1179                 }
1180         } else {
1181                 omxMatrix *design = state->design;
1182                 if (design->cols != numItems ||
1183                     design->rows != state->maxDims) {
1184                         omxRaiseErrorf(currentState, "Design matrix should have %d rows and %d columns",
1185                                        state->maxDims, numItems);
1186                         return;
1187                 }
1188
1189                 state->maxAbilities = 0;
1190                 for (int ix=0; ix < design->rows * design->cols; ix++) {
1191                         double got = design->data[ix];
1192                         if (!R_FINITE(got)) continue;
1193                         if (round(got) != got) error("Design matrix can only contain integers"); // TODO better way?
1194                         if (state->maxAbilities < got)
1195                                 state->maxAbilities = got;
1196                 }
1197         }
1198         if (state->maxAbilities <= state->maxDims) {
1199                 state->Sgroup = Calloc(numItems, int);
1200         } else {
1201                 int Sgroup0 = state->maxDims;
1202                 state->Sgroup = Realloc(NULL, numItems, int);
1203                 for (int ix=0; ix < numItems; ix++) {
1204                         int ss=-1;
1205                         for (int dx=0; dx < state->maxDims; dx++) {
1206                                 int ability = omxMatrixElement(state->design, dx, ix);
1207                                 if (ability >= Sgroup0) {
1208                                         if (ss == -1) {
1209                                                 ss = ability;
1210                                         } else {
1211                                                 omxRaiseErrorf(currentState, "Item %d cannot belong to more than "
1212                                                                "1 specific dimension (both %d and %d)",
1213                                                                ix, ss, ability);
1214                                                 return;
1215                                         }
1216                                 }
1217                         }
1218                         if (ss == -1) ss = 0;
1219                         state->Sgroup[ix] = ss - Sgroup0;
1220                 }
1221                 state->numSpecific = state->maxAbilities - state->maxDims + 1;
1222                 state->allSlxk = Realloc(NULL, numUnique * numThreads, double);
1223                 state->Slxk = Realloc(NULL, numUnique * state->numSpecific * numThreads, double);
1224         }
1225
1226         PROTECT(tmp = GET_SLOT(rObj, install("cache")));
1227         state->cacheLXK = asLogical(tmp);
1228
1229         PROTECT(tmp = GET_SLOT(rObj, install("GHpoints")));
1230         double *qpoints = REAL(tmp);
1231         int numQPoints = length(tmp);
1232
1233         PROTECT(tmp = GET_SLOT(rObj, install("GHarea")));
1234         double *qarea = REAL(tmp);
1235         if (numQPoints != length(tmp)) error("length(GHpoints) != length(GHarea)");
1236
1237         ba81SetupQuadrature(oo, numQPoints, qpoints, qarea);
1238
1239         // verify data bounded between 1 and numOutcomes TODO
1240         // hm, looks like something could be added to omxData for column summary stats?
1241 }